Skip to content

[PWGJE] Add a new histogram for the DCAxy template fit #14609

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
ynishida-style wants to merge 1 commit into
base: master
Choose a base branch
from
+101 −32
Open

[PWGJE] Add a new histogram for the DCAxy template fit #14609

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
133 changes: 101 additions & 32 deletions PWGJE/Tasks/jetShape.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,8 @@ struct JetShapeTask {
Configurable<int> nBinsPForBeta{"nBinsPForBeta", 500, "Number of pT bins"};
Configurable<int> nBinsTpcDedx{"nBinsTpcDedx", 500, "Number of DEdx bins"};
Configurable<int> nBinsTofBeta{"nBinsTofBeta", 350, "Number of Beta bins"};
Configurable<float> dcaxyMin{"dcaxyMin", -1.0f, "Min value of dcaXY"};
Configurable<float> dcaxyMax{"dcaxyMax", 1.0f, "Max value of dcaXY"};
Configurable<float> pMax{"pMax", 8.0f, "Max value of p"};
Configurable<float> ptMax{"ptMax", 6.0f, "Max value of pT"};
Configurable<float> jetPtMinForCut{"jetPtMinForCut", 0.0f, "Minimum value of jet pT cut"};
Expand All @@ -58,8 +60,11 @@ struct JetShapeTask {
Configurable<float> centralityMaxForCut{"centralityMaxForCut", 100.0f, "Maximum value of the jet pT cut"};
Configurable<float> jetShapeFuncMax{"jetShapeFuncMax", 300, "Maximum value of JetShapeFunction"};
Configurable<int> nBinsJetShapeFunc{"nBinsJetShapeFunc", 900, "Number of JetShapeFunction bins"};
Configurable<int> nBinsP{"nBinsP", 80, "Number of p bins"};
Configurable<int> nBinsPt{"nBinsPt", 60, "Number of pT bins"};
Configurable<int> nBinsDcaxyForData{"nBinsDcaxyForData", 400, "Number of DcaXY bins for data"};
Configurable<int> nBinsDcaxyForMc{"nBinsDcaxyForMc", 400, "Number of DcaXY bins for mc data"};
Configurable<int> nBinsP{"nBinsP", 40, "Number of p bins"};
Configurable<int> nBinsPt{"nBinsPt", 30, "Number of pT bins"};
Configurable<int> nBinsPtForDca{"nBinsPtForDca", 15, "Number of pT bins for dcaXY"};
Configurable<int> nBinsJetPt{"nBinsJetPt", 10, "Number of jet pT bins"};
Configurable<int> nBinsPForCut{"nBinsPForCut", 30, "Number of p track bins"};
Configurable<int> nBinsCentrality{"nBinsCentrality", 10, "Number of centrality bins"};
Expand Down Expand Up @@ -106,6 +111,10 @@ struct JetShapeTask {
{"jetpVsPtForPi", "jetpVsPtPi", {HistType::kTHnSparseD, {{nBinsP, 0, pMax}, {nBinsPt, 0, ptMax}, {nBinsDistance, 0, distanceMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"pVsPtForPrOutOfJet", "pVsPtForPrOutOfJet", {HistType::kTHnSparseD, {{nBinsP, 0, pMax}, {nBinsPt, 0, ptMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"pVsPtForPiOutOfJet", "pVsPtPionOutOfJet", {HistType::kTHnSparseD, {{nBinsP, 0, pMax}, {nBinsPt, 0, ptMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"jetDcaPr", "jetDcaPr", {HistType::kTHnSparseD, {{nBinsPtForDca, 0, ptMax}, {nBinsDcaxyForData, dcaxyMin, dcaxyMax}, {nBinsDistance, 0, distanceMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"jetDcaPi", "jetDcaPi", {HistType::kTHnSparseD, {{nBinsPtForDca, 0, ptMax}, {nBinsDcaxyForData, dcaxyMin, dcaxyMax}, {nBinsDistance, 0, distanceMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"dcaPrOutOfJet", "dcaPrOutOfJet", {HistType::kTHnSparseD, {{nBinsPtForDca, 0, ptMax}, {nBinsDcaxyForData, dcaxyMin, dcaxyMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"dcaPiOutOfJet", "dcaPiOutOfJet", {HistType::kTHnSparseD, {{nBinsPtForDca, 0, ptMax}, {nBinsDcaxyForData, dcaxyMin, dcaxyMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"jetPt", "jet pT;#it{p}_{T,jet} (GeV/#it{c});entries", {HistType::kTH1F, {{200, 0., 200.}}}},
{"jetEta", "jet #eta;#eta_{jet};entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}},
{"jetPhi", "jet #phi;#phi_{jet};entries", {HistType::kTH1F, {{80, -1.0, 7.}}}},
Expand All @@ -118,7 +127,10 @@ struct JetShapeTask {
{"ptSumBg1", "ptSumBg1", {HistType::kTHnSparseD, {{14, 0, 0.7}, {nBinsJetShapeFunc, 0, jetShapeFuncMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"ptSumBg2", "ptSumBg2", {HistType::kTHnSparseD, {{14, 0, 0.7}, {nBinsJetShapeFunc, 0, jetShapeFuncMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"event/vertexz", ";Vtx_{z} (cm);Entries", {HistType::kTH1F, {{100, -20, 20}}}},
{"eventCounter", "eventCounter", {HistType::kTH1F, {{1, 0, +1, ""}}}},
{"eventCounterJetShape", "eventCounterJetShape", {HistType::kTH1F, {{1, 0, +1, ""}}}},
{"eventCounterJet", "eventCounterJet", {HistType::kTH1F, {{1, 0, +1, ""}}}},
{"eventCounterInc", "eventCounterInc", {HistType::kTH1F, {{1, 0, +1, ""}}}},
{"eventCounterMc", "eventCounterMc", {HistType::kTH1F, {{1, 0, +1, ""}}}},
{"ptVsCentrality", "ptvscentrality", {HistType::kTH2F, {{100, 0, 100}, {300, 0, 300}}}},
{"ptResolution", "ptResolution", {HistType::kTH2F, {{nBinsPt, 0, ptMax}, {100, -1.0, +1.0}}}},
{"mcCentralityReco", "mcCentralityReco", {HistType::kTH1F, {{100, 0, 100}}}},
Expand All @@ -129,6 +141,10 @@ struct JetShapeTask {
{"ptHistogramPionTof", "ptHistogramPionTof", {HistType::kTHnSparseD, {{nBinsPt, 0, ptMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"ptHistogramKaonTof", "ptHistogramKaonTof", {HistType::kTHnSparseD, {{nBinsPt, 0, ptMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"ptHistogramProtonTof", "ptHistogramProtonTof", {HistType::kTHnSparseD, {{nBinsPt, 0, ptMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"dcaDecayPion", "dcaDecayPion", {HistType::kTHnSparseD, {{nBinsPt, 0, ptMax}, {nBinsDcaxyForMc, dcaxyMin, dcaxyMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"dcaDecayProton", "dcaDecayProton", {HistType::kTHnSparseD, {{nBinsPt, 0, ptMax}, {nBinsDcaxyForMc, dcaxyMin, dcaxyMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"dcaMaterialPion", "dcaMaterialPion", {HistType::kTHnSparseD, {{nBinsPt, 0, ptMax}, {nBinsDcaxyForMc, dcaxyMin, dcaxyMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"dcaMaterialProton", "dcaMaterialProton", {HistType::kTHnSparseD, {{nBinsPt, 0, ptMax}, {nBinsDcaxyForMc, dcaxyMin, dcaxyMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"ptGeneratedPion", "ptGeneratedPion", {HistType::kTHnSparseD, {{nBinsPt, 0, ptMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"ptGeneratedKaon", "ptGeneratedKaon", {HistType::kTHnSparseD, {{nBinsPt, 0, ptMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}},
{"ptGeneratedProton", "ptGeneratedProton", {HistType::kTHnSparseD, {{nBinsPt, 0, ptMax}, {nBinsJetPt, jetPtMinForCut, jetPtMaxForCut}, {nBinsCentrality, centralityMinForCut, centralityMaxForCut}}}}}};
Expand All @@ -151,7 +167,7 @@ struct JetShapeTask {

// for ppi production
Configurable<float> etaTrUp{"etaTrUp", 0.7f, "maximum track eta"};
Configurable<float> dcaxyMax{"dcaxyMax", 2.0f, "maximum DCA xy"};
Configurable<float> dcaxyCutMax{"dcaxyCutMax", 2.0f, "maximum DCA xy"};
Configurable<float> chi2ItsMax{"chi2ItsMax", 15.0f, "its chi2 cut"};
Configurable<float> chi2TpcMax{"chi2TpcMax", 4.0f, "tpc chi2 cut"};
Configurable<float> nclItsMin{"nclItsMin", 2.0f, "its # of cluster cut"};
Expand All @@ -160,6 +176,7 @@ struct JetShapeTask {
Configurable<float> mcRapidityMax{"mcRapidityMax", 0.5f, "maximum mctrack y"};
Configurable<double> epsilon{"epsilon", 1e-6, "standard for aboid division of zero"};
Configurable<float> maxDeltaEtaSafe{"maxDeltaEtaSafe", 0.9f, "maximum track eta for cut"};
Configurable<float> nSigmaMaxForDcaxy{"nSigmaMaxForDcaxy", 4.0f, "maximum nSigma for DCAxy"};

Configurable<std::string> triggerMasks{"triggerMasks", "", "possible JE Trigger masks: fJetChLowPt,fJetChHighPt,fTrackLowPt,fTrackHighPt,fJetD0ChLowPt,fJetD0ChHighPt,fJetLcChLowPt,fJetLcChHighPt,fEMCALReadout,fJetFullHighPt,fJetFullLowPt,fJetNeutralHighPt,fJetNeutralLowPt,fGammaVeryHighPtEMCAL,fGammaVeryHighPtDCAL,fGammaHighPtEMCAL,fGammaHighPtDCAL,fGammaLowPtEMCAL,fGammaLowPtDCAL,fGammaVeryLowPtEMCAL,fGammaVeryLowPtDCAL"};

Expand Down Expand Up @@ -233,6 +250,8 @@ struct JetShapeTask {
return;
}

registry.fill(HIST("eventCounterJetShape"), 0.5);

size_t nBins = distanceCategory->size() - 1;

float maxDistance = distanceCategory->at(nBins);
Expand Down Expand Up @@ -378,7 +397,7 @@ struct JetShapeTask {
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
return;
}

registry.fill(HIST("eventCounterJet"), 0.5);
registry.fill(HIST("event/vertexz"), collision.posZ());

float rho = collision.rho();
Expand Down Expand Up @@ -419,7 +438,7 @@ struct JetShapeTask {
continue;
if (track.tpcNClsCrossedRows() < nclcrossTpcMin)
continue;
if (std::abs(track.dcaXY()) > dcaxyMax)
if (std::abs(track.dcaXY()) > dcaxyCutMax)
continue;
if (track.itsChi2NCl() > chi2ItsMax)
continue;
Expand Down Expand Up @@ -453,6 +472,9 @@ struct JetShapeTask {
bool isTpcPiRange = (tpcPi > tpcNSigmaPiMin && tpcPi < tpcNSigmaPiMax);
bool isTpcPrRange = (tpcPr > tpcNSigmaPrMin && tpcPr < tpcNSigmaPrMax);

float nSigmaSqPr = tpcPr * tpcPr + tofPr * tofPr;
float nSigmaSqPi = tpcPi * tpcPi + tofPi * tofPi;

for (const auto& jet : cachedJets) {

float dEta = trkEta - jet.eta;
Expand All @@ -472,16 +494,27 @@ struct JetShapeTask {
if (distBg1 < distanceMax || distBg2 < distanceMax) {
registry.fill(HIST("tpcDedxOutOfJet"), trkP, tpcSig);

// dcaXY
if (track.hasTOF()) {
if (nSigmaSqPr < nSigmaMaxForDcaxy) {
registry.fill(HIST("dcaPrOutOfJet"), trkPt, track.dcaXY(), jet.ptCorr, centrality);
}

if (nSigmaSqPi < nSigmaMaxForDcaxy) {
registry.fill(HIST("dcaPiOutOfJet"), trkPt, track.dcaXY(), jet.ptCorr, centrality);
}
}

if (hasTofPi) {
registry.fill(HIST("tpcTofPiOutOfJet"), trkP, tpcPi, jet.ptCorr, centrality);
if (isTpcPiRange) {
registry.fill(HIST("pVsPtForPiOutOfJet"), trkP, trkPt, jet.ptCorr, centrality);
registry.fill(HIST("pVsPtForPiOutOfJet"), trkP, trkPt, jet.ptCorr, centrality, track.dcaXY());
}
}
if (hasTofPr) {
registry.fill(HIST("tpcTofPrOutOfJet"), trkP, tpcPr, jet.ptCorr, centrality);
if (isTpcPrRange) {
registry.fill(HIST("pVsPtForPrOutOfJet"), trkP, trkPt, jet.ptCorr, centrality);
registry.fill(HIST("pVsPtForPrOutOfJet"), trkP, trkPt, jet.ptCorr, centrality, track.dcaXY());
}
}
}
Expand All @@ -491,6 +524,17 @@ struct JetShapeTask {
registry.fill(HIST("jetTpcDedx"), trkP, tpcSig, distance);
registry.fill(HIST("jetTofBeta"), trkP, beta);

// dcaXY
if (track.hasTOF()) {
if (nSigmaSqPr < nSigmaMaxForDcaxy) {
registry.fill(HIST("jetDcaPr"), trkPt, track.dcaXY(), distance, jet.ptCorr, centrality);
}

if (nSigmaSqPi < nSigmaMaxForDcaxy) {
registry.fill(HIST("jetDcaPi"), trkPt, track.dcaXY(), distance, jet.ptCorr, centrality);
}
}

if (hasTofPr) {
registry.fill(HIST("jetTpcTofPr"), trkP, tpcPr, distance, jet.ptCorr, centrality);
if (isTpcPrRange) {
Expand All @@ -514,7 +558,7 @@ struct JetShapeTask {
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
return;
}

registry.fill(HIST("eventCounterJet"), 0.5);
// tracks conditions
for (const auto& jetTrack : tracks) {

Expand All @@ -537,7 +581,7 @@ struct JetShapeTask {
continue;
if (track.tpcNClsCrossedRows() < nclcrossTpcMin)
continue;
if (std::abs(track.dcaXY()) > dcaxyMax)
if (std::abs(track.dcaXY()) > dcaxyCutMax)
continue;
if (track.itsChi2NCl() > chi2ItsMax)
continue;
Expand Down Expand Up @@ -582,11 +626,10 @@ struct JetShapeTask {

(void)mcParticles;

registry.fill(HIST("eventCounter"), 0.5);

float centrality = collision.centFT0M();
float rho = collision.rho();

registry.fill(HIST("eventCounterMc"), 0.5);
registry.fill(HIST("mcCentralityReco"), centrality);

struct CachedJet {
Expand All @@ -599,12 +642,12 @@ struct JetShapeTask {
cachedJets.reserve(jets.size());

for (const auto& jet : jets) {
registry.fill(HIST("jetPt"), jet.pt());

float mcdPtCorr = jet.pt() - rho * jet.area();
cachedJets.push_back({jet.pt(), jet.eta(), jet.phi(), mcdPtCorr});
registry.fill(HIST("jetPt"), jet.pt());
}

// reco track loop
for (const auto& track : tracks) {

if (!jetderiveddatautilities::selectTrack(track, trackSelection)) {
Expand All @@ -618,7 +661,7 @@ struct JetShapeTask {
continue;
if (track.tpcNClsCrossedRows() < nclcrossTpcMin)
continue;
if (std::abs(track.dcaXY()) > dcaxyMax)
if (std::abs(track.dcaXY()) > dcaxyCutMax)
continue;
if (track.itsChi2NCl() > chi2ItsMax)
continue;
Expand All @@ -632,9 +675,15 @@ struct JetShapeTask {
auto mcParticle = track.mcParticle();
registry.fill(HIST("ptResolution"), track.pt(), track.pt() - mcParticle.pt());

if (!mcParticle.isPhysicalPrimary() || std::fabs(mcParticle.y()) >= mcRapidityMax)
if (std::fabs(mcParticle.y()) >= mcRapidityMax)
continue;

const int producedByDecay = 4;

bool isPrimary = mcParticle.isPhysicalPrimary();
bool isSecondDecay = !isPrimary && (mcParticle.getProcess() == producedByDecay);
bool isSecondMaterial = !isPrimary && !isSecondDecay;

int pdg = std::abs(mcParticle.pdgCode());
bool isPion = (pdg == PDG_t::kPiPlus);
bool isKaon = (pdg == PDG_t::kKPlus);
Expand Down Expand Up @@ -663,38 +712,58 @@ struct JetShapeTask {
if (deltaR > distanceMax)
continue;

// TPC (All matched)
if (isPion)
registry.fill(HIST("ptHistogramPion"), mcParticle.pt(), jet.ptCorr, centrality);
else if (isKaon)
registry.fill(HIST("ptHistogramKaon"), mcParticle.pt(), jet.ptCorr, centrality);
else if (isProton)
registry.fill(HIST("ptHistogramProton"), mcParticle.pt(), jet.ptCorr, centrality);

// TOF (Required)
if (hasTof) {
if (isPrimary) {
// Tracking
if (isPion)
registry.fill(HIST("ptHistogramPionTof"), mcParticle.pt(), jet.ptCorr, centrality);
registry.fill(HIST("ptHistogramPion"), mcParticle.pt(), jet.ptCorr, centrality);
else if (isKaon)
registry.fill(HIST("ptHistogramKaonTof"), mcParticle.pt(), jet.ptCorr, centrality);
registry.fill(HIST("ptHistogramKaon"), mcParticle.pt(), jet.ptCorr, centrality);
else if (isProton)
registry.fill(HIST("ptHistogramProtonTof"), mcParticle.pt(), jet.ptCorr, centrality);
registry.fill(HIST("ptHistogramProton"), mcParticle.pt(), jet.ptCorr, centrality);

// TOF matched
if (hasTof) {
if (isPion)
registry.fill(HIST("ptHistogramPionTof"), mcParticle.pt(), jet.ptCorr, centrality);
else if (isKaon)
registry.fill(HIST("ptHistogramKaonTof"), mcParticle.pt(), jet.ptCorr, centrality);
else if (isProton)
registry.fill(HIST("ptHistogramProtonTof"), mcParticle.pt(), jet.ptCorr, centrality);
}
} else { // Secondary
if (isSecondDecay) {
// from Decay
if (isPion)
registry.fill(HIST("dcaDecayPion"), mcParticle.pt(), track.dcaXY(), jet.ptCorr, centrality);
else if (isProton)
registry.fill(HIST("dcaDecayProton"), mcParticle.pt(), track.dcaXY(), jet.ptCorr, centrality);
} else if (isSecondMaterial) {
// from Material
if (isPion)
registry.fill(HIST("dcaMaterialPion"), mcParticle.pt(), track.dcaXY(), jet.ptCorr, centrality);
else if (isProton)
registry.fill(HIST("dcaMaterialProton"), mcParticle.pt(), track.dcaXY(), jet.ptCorr, centrality);
}
}
}
}
}
PROCESS_SWITCH(JetShapeTask, processReco, "process reconstructed simulation information", true);

void processSim(aod::JetMcCollisions::iterator const& mcCollision, aod::ChargedMCParticleLevelJets const& mcpjets, aod::JetParticles const& mcParticles)
void processSim(aod::JetMcCollisions::iterator const& mcCollision, soa::SmallGroups<aod::JetCollisionsMCD> const& collisions, aod::ChargedMCParticleLevelJets const& mcpjets, aod::JetParticles const& mcParticles)
{
if (std::abs(mcCollision.posZ()) > vertexZCut) {
return;
}

if (collisions.size() == 0) {
return;
}

// --- centrality ---
float centrality = mcCollision.centFT0M();
float centrality = collisions.begin().centFT0M();
// float centrality = mcCollision.centFT0M();
registry.fill(HIST("mcCentralitySim"), centrality);

const float maxR2 = distanceMax * distanceMax;

// --- loop over MC particles only once ---
Expand Down
Loading