Skip to content

Commit 1bd02f7

Browse files
dajones2dajones2
committed
PWGJE: Adding task for denominator consistency for effieicncy calculation
1 parent 5d2ea8a commit 1bd02f7

File tree

1 file changed

+145
-5
lines changed

1 file changed

+145
-5
lines changed

‎PWGJE/Tasks/jetHadronRecoil.cxx‎

Lines changed: 145 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -140,7 +140,7 @@ struct JetHadronRecoil {
140140

141141
registry.add("hZvtxSelected", "Z vertex position;Z_{vtx};entries", {HistType::kTH1F, {{80, -20, 20}}}, doSumw);
142142

143-
if (doprocessData || doprocessDataWithRhoSubtraction || doprocessMCD || doprocessMCDWithRhoSubtraction || doprocessMCDWeighted || doprocessMCDWeightedWithRhoSubtraction || doprocessMCP || doprocessMCPWeighted) {
143+
if (doprocessData || doprocessDataWithRhoSubtraction || doprocessMCD || doprocessMCDWithRhoSubtraction || doprocessMCDWeighted || doprocessMCDWeightedWithRhoSubtraction || doprocessMCP || doprocessMCPWeighted || doprocessMCPWeightedWithMatchedTracks) {
144144
registry.add("hNtrig", "number of triggers;trigger type;entries", {HistType::kTH1F, {{2, 0, 2}}}, doSumw);
145145
registry.add("hSignalTriggersPtHard", "Signal triggers vs PtHard", {HistType::kTH1F, {{20, 0, 5}}}, doSumw);
146146
registry.add("hReferenceTriggersPtHard", "Reference triggers vs PtHard", {HistType::kTH1F, {{20, 0, 5}}}, doSumw);
@@ -181,7 +181,7 @@ struct JetHadronRecoil {
181181
registry.add("hDeltaRpTDPhiReferenceShifts", "testing shifts;p_{T,jet};#Delta#phi;#DeltaR;shifts", {HistType::kTHnSparseD, {{500, -100, 400}, {100, 0, o2::constants::math::TwoPI}, dRAxis, {20, 0.0, 2.0}}}, doSumw);
182182
}
183183

184-
if (doprocessMCP || doprocessMCPWeighted) {
184+
if (doprocessMCP || doprocessMCPWeighted || doprocessMCPWeightedWithMatchedTracks) {
185185
registry.add("hPartvsJets", "comparing leading particles and jets;p_{T,part};p_{T,jet};#hat{p}", {HistType::kTH3F, {{200, 0, 200}, {500, -100, 400}, {195, 5, 200}}}, doSumw);
186186
registry.add("hPtPart", "Particle p_{T};p_{T};entries", {HistType::kTH1F, {{200, 0, 200}}}, doSumw);
187187
registry.add("hEtaPart", "Particle #eta;#eta;entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}, doSumw);
@@ -386,11 +386,13 @@ struct JetHadronRecoil {
386386
phiTTAr.push_back(particle.phi());
387387
ptTTAr.push_back(particle.pt());
388388
nTT++;
389+
registry.fill(HIST("hSignalTriggers"), track.pt(), weight);
389390
}
390391
if (!isSigCol && particle.pt() < ptTTrefMax && particle.pt() > ptTTrefMin) {
391392
phiTTAr.push_back(particle.phi());
392393
ptTTAr.push_back(particle.pt());
393394
nTT++;
395+
registry.fill(HIST("hReferenceTriggers"), track.pt(), weight);
394396
}
395397
registry.fill(HIST("hPtPart"), particle.pt(), weight);
396398
registry.fill(HIST("hEtaPart"), particle.eta(), weight);
@@ -466,6 +468,116 @@ struct JetHadronRecoil {
466468
registry.fill(HIST("hPartvsJets"), leadingPartPt, leadingJetPt, pTHat, weight);
467469
}
468470

471+
template <typename T, typename U, typename V>
472+
void fillMCPHistogramsWithMatchedTracks(T const& jets, U const&, V const& tracks, float weight = 1.0, float pTHat = 999.0)
473+
{
474+
bool isSigCol;
475+
std::vector<double> phiTTAr;
476+
std::vector<double> ptTTAr;
477+
double phiTT = 0;
478+
double ptTT = 0;
479+
int trigNumber = 0;
480+
int nTT = 0;
481+
double leadingPartPt = 0;
482+
double leadingJetPt = 0;
483+
float dice = rand->Rndm();
484+
if (dice < fracSig)
485+
isSigCol = true;
486+
else
487+
isSigCol = false;
488+
489+
for (const auto& track : tracks) {
490+
if (!track.has_mcParticle()) {
491+
continue;
492+
}
493+
auto particle = track.template mcParticle_as<U>();
494+
if (particle.pt() > leadingPartPt) {
495+
leadingPartPt = particle.pt();
496+
}
497+
if (particle.pt() > pTHatTrackMaxMCD * pTHat) {
498+
if (outlierRejectEvent) {
499+
return;
500+
} else {
501+
continue;
502+
}
503+
}
504+
auto pdgParticle = pdg->GetParticle(particle.pdgCode());
505+
if (!pdgParticle) {
506+
continue;
507+
}
508+
if ((pdgParticle->Charge() == 0.0) || (!particle.isPhysicalPrimary())) {
509+
continue;
510+
}
511+
if (isSigCol && particle.pt() < ptTTsigMax && particle.pt() > ptTTsigMin) {
512+
phiTTAr.push_back(particle.phi());
513+
ptTTAr.push_back(particle.pt());
514+
nTT++;
515+
registry.fill(HIST("hSignalTriggers"), track.pt(), weight);
516+
}
517+
if (!isSigCol && particle.pt() < ptTTrefMax && particle.pt() > ptTTrefMin) {
518+
phiTTAr.push_back(particle.phi());
519+
ptTTAr.push_back(particle.pt());
520+
nTT++;
521+
registry.fill(HIST("hReferenceTriggers"), track.pt(), weight);
522+
}
523+
registry.fill(HIST("hPtPart"), particle.pt(), weight);
524+
registry.fill(HIST("hEtaPart"), particle.eta(), weight);
525+
registry.fill(HIST("hPhiPart"), particle.phi(), weight);
526+
registry.fill(HIST("hPart3D"), particle.pt(), particle.eta(), particle.phi(), weight);
527+
registry.fill(HIST("hPtPartPtHard"), particle.pt(), particle.pt() / pTHat, weight);
528+
}
529+
530+
for (const auto& jet : jets) {
531+
if (jet.pt() > leadingJetPt) {
532+
leadingJetPt = jet.pt();
533+
}
534+
if (jet.pt() > pTHatMaxMCP * pTHat) {
535+
if (outlierRejectEvent) {
536+
return;
537+
} else {
538+
continue;
539+
}
540+
}
541+
for (const auto& constituent : jet.template tracks_as<U>()) {
542+
registry.fill(HIST("hConstituents3D"), constituent.pt(), constituent.eta(), constituent.phi());
543+
}
544+
registry.fill(HIST("hJetPt"), jet.pt(), weight);
545+
registry.fill(HIST("hJetEta"), jet.eta(), weight);
546+
registry.fill(HIST("hJetPhi"), jet.phi(), weight);
547+
registry.fill(HIST("hJet3D"), jet.pt(), jet.eta(), jet.phi(), weight);
548+
549+
if (nTT > 0) {
550+
float dphi = RecoDecay::constrainAngle(jet.phi() - phiTT);
551+
double dR = getWTAaxisDifference(jet, particles);
552+
if (isSigCol) {
553+
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {
554+
registry.fill(HIST("hDeltaRpTSignalPart"), jet.pt(), dR, weight);
555+
registry.fill(HIST("hDeltaRSignalPart"), dR, weight);
556+
}
557+
registry.fill(HIST("hDeltaRpTDPhiSignalPart"), jet.pt(), dphi, dR, weight);
558+
registry.fill(HIST("hSignalPtDPhi"), dphi, jet.pt(), weight);
559+
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {
560+
registry.fill(HIST("hSignalPt"), jet.pt(), weight);
561+
registry.fill(HIST("hSignalPtHard"), jet.pt(), ptTT / pTHat, weight);
562+
}
563+
}
564+
if (!isSigCol) {
565+
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {
566+
registry.fill(HIST("hDeltaRpTPartReference"), jet.pt(), dR, weight);
567+
registry.fill(HIST("hDeltaRPartReference"), dR, weight);
568+
}
569+
registry.fill(HIST("hDeltaRpTDPhiReferencePart"), jet.pt(), dphi, dR, weight);
570+
registry.fill(HIST("hReferencePtDPhi"), dphi, jet.pt(), weight);
571+
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {
572+
registry.fill(HIST("hReferencePt"), jet.pt(), weight);
573+
registry.fill(HIST("hReferencePtHard"), jet.pt(), ptTT / pTHat, weight);
574+
}
575+
}
576+
}
577+
}
578+
registry.fill(HIST("hPartvsJets"), leadingPartPt, leadingJetPt, pTHat, weight);
579+
}
580+
469581
template <typename T, typename U, typename X, typename Y>
470582
void fillMatchedHistograms(T const& jetsBase, U const&, X const& tracks, Y const& particles, float weight = 1.0, float rho = 0.0, float pTHat = 999.0)
471583
{
@@ -534,10 +646,19 @@ struct JetHadronRecoil {
534646
}
535647
}
536648
if (track.pt() < ptTTsigMax && track.pt() > ptTTsigMin) {
537-
phiTTAr.push_back(track.phi());
538-
nTT++;
539649
auto particle = track.template mcParticle_as<Y>();
540-
phiTTArPart.push_back(particle.phi());
650+
auto pdgParticle = pdg->GetParticle(particle.pdgCode());
651+
if (!pdgParticle) {
652+
continue;
653+
}
654+
if ((pdgParticle->Charge() == 0.0) || (!particle.isPhysicalPrimary())) {
655+
continue;
656+
}
657+
if (particle.pt() < ptTTsigMax && particle.pt() > ptTTsigMin) {
658+
nTT++;
659+
phiTTAr.push_back(track.phi());
660+
phiTTArPart.push_back(particle.phi());
661+
}
541662
}
542663
}
543664

@@ -758,6 +879,25 @@ struct JetHadronRecoil {
758879
}
759880
PROCESS_SWITCH(JetHadronRecoil, processMCPWeighted, "process MC particle level with event weights", false);
760881

882+
void processMCPWeightedWithMatchedTracks(aod::JetMcCollision const& collision,
883+
soa::Filtered<soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents>> const& jets,
884+
soa::Filtered<aod::JetParticles> const& particles,
885+
soa::Filtered<aod::JetTracks> const& tracks)
886+
{
887+
if (std::abs(collision.posZ()) > vertexZCut) {
888+
return;
889+
}
890+
if (skipMBGapEvents && collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) {
891+
return;
892+
}
893+
if (collision.ptHard() < pTHatMinEvent) {
894+
return;
895+
}
896+
registry.fill(HIST("hZvtxSelected"), collision.posZ(), collision.weight());
897+
fillMCPHistogramsWithMatchedTracks(jets, particles, tracks, collision.weight(), collision.ptHard());
898+
}
899+
PROCESS_SWITCH(JetHadronRecoil, processMCPWeightedWithMatchedTracks, "process MC particle level with event weights - only triggers matched with detector level tracks", false);
900+
761901
void processJetsMCPMCDMatched(soa::Filtered<soa::Join<aod::JetCollisionsMCD, aod::JMcCollisionLbs>>::iterator const& collision,
762902
soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets>> const& mcdjets,
763903
aod::JetTracks const& tracks,

0 commit comments

Comments
 (0)