Skip to content

Commit e6841fa

Browse files
dajones2dajones2
authored
[PWGJE] Adding task for denominator consistency in effieicncy calcula… (#13255)
1 parent 6899681 commit e6841fa

File tree

1 file changed

+161
-5
lines changed

1 file changed

+161
-5
lines changed

‎PWGJE/Tasks/jetHadronRecoil.cxx‎

Lines changed: 161 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -143,7 +143,7 @@ struct JetHadronRecoil {
143143

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

146-
if (doprocessData || doprocessDataWithRhoSubtraction || doprocessMCD || doprocessMCDWithRhoSubtraction || doprocessMCDWeighted || doprocessMCDWeightedWithRhoSubtraction || doprocessMCP || doprocessMCPWeighted) {
146+
if (doprocessData || doprocessDataWithRhoSubtraction || doprocessMCD || doprocessMCDWithRhoSubtraction || doprocessMCDWeighted || doprocessMCDWeightedWithRhoSubtraction || doprocessMCP || doprocessMCPWeighted || doprocessMCPWeightedWithMatchedTracks) {
147147
registry.add("hNtrig", "number of triggers;trigger type;entries", {HistType::kTH1F, {{2, 0, 2}}}, doSumw);
148148
registry.add("hSignalTriggersPtHard", "Signal triggers vs PtHard", {HistType::kTH1F, {pThatAxis}}, doSumw);
149149
registry.add("hReferenceTriggersPtHard", "Reference triggers vs PtHard", {HistType::kTH1F, {pThatAxis}}, doSumw);
@@ -184,7 +184,7 @@ struct JetHadronRecoil {
184184
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);
185185
}
186186

187-
if (doprocessMCP || doprocessMCPWeighted) {
187+
if (doprocessMCP || doprocessMCPWeighted || doprocessMCPWeightedWithMatchedTracks) {
188188
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);
189189
registry.add("hPtPart", "Particle p_{T};p_{T};entries", {HistType::kTH1F, {{200, 0, 200}}}, doSumw);
190190
registry.add("hEtaPart", "Particle #eta;#eta;entries", {HistType::kTH1F, {{100, -1.0, 1.0}}}, doSumw);
@@ -389,11 +389,139 @@ struct JetHadronRecoil {
389389
phiTTAr.push_back(particle.phi());
390390
ptTTAr.push_back(particle.pt());
391391
nTT++;
392+
registry.fill(HIST("hSignalTriggers"), particle.pt(), weight);
392393
}
393394
if (!isSigCol && particle.pt() < ptTTrefMax && particle.pt() > ptTTrefMin) {
394395
phiTTAr.push_back(particle.phi());
395396
ptTTAr.push_back(particle.pt());
396397
nTT++;
398+
registry.fill(HIST("hReferenceTriggers"), particle.pt(), weight);
399+
}
400+
registry.fill(HIST("hPtPart"), particle.pt(), weight);
401+
registry.fill(HIST("hEtaPart"), particle.eta(), weight);
402+
registry.fill(HIST("hPhiPart"), particle.phi(), weight);
403+
registry.fill(HIST("hPart3D"), particle.pt(), particle.eta(), particle.phi(), weight);
404+
registry.fill(HIST("hPtPartPtHard"), particle.pt(), particle.pt() / pTHat, weight);
405+
}
406+
407+
if (nTT > 0) {
408+
trigNumber = rand->Integer(nTT);
409+
phiTT = phiTTAr[trigNumber];
410+
ptTT = ptTTAr[trigNumber];
411+
if (isSigCol) {
412+
registry.fill(HIST("hNtrig"), 1.5, weight);
413+
registry.fill(HIST("hSigEventTriggers"), nTT, weight);
414+
registry.fill(HIST("hSignalTriggersPtHard"), ptTT / pTHat, weight);
415+
}
416+
if (!isSigCol) {
417+
registry.fill(HIST("hNtrig"), 0.5, weight);
418+
registry.fill(HIST("hRefEventTriggers"), nTT, weight);
419+
registry.fill(HIST("hReferenceTriggersPtHard"), ptTT / pTHat, weight);
420+
}
421+
}
422+
423+
for (const auto& jet : jets) {
424+
if (jet.pt() > leadingJetPt) {
425+
leadingJetPt = jet.pt();
426+
}
427+
if (jet.pt() > pTHatMaxMCP * pTHat) {
428+
if (outlierRejectEvent) {
429+
return;
430+
} else {
431+
continue;
432+
}
433+
}
434+
for (const auto& constituent : jet.template tracks_as<U>()) {
435+
registry.fill(HIST("hConstituents3D"), constituent.pt(), constituent.eta(), constituent.phi());
436+
}
437+
registry.fill(HIST("hJetPt"), jet.pt(), weight);
438+
registry.fill(HIST("hJetEta"), jet.eta(), weight);
439+
registry.fill(HIST("hJetPhi"), jet.phi(), weight);
440+
registry.fill(HIST("hJet3D"), jet.pt(), jet.eta(), jet.phi(), weight);
441+
442+
if (nTT > 0) {
443+
float dphi = RecoDecay::constrainAngle(jet.phi() - phiTT);
444+
double dR = getWTAaxisDifference(jet, particles);
445+
if (isSigCol) {
446+
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {
447+
registry.fill(HIST("hDeltaRpTSignalPart"), jet.pt(), dR, weight);
448+
registry.fill(HIST("hDeltaRSignalPart"), dR, weight);
449+
}
450+
registry.fill(HIST("hDeltaRpTDPhiSignalPart"), jet.pt(), dphi, dR, weight);
451+
registry.fill(HIST("hSignalPtDPhi"), dphi, jet.pt(), weight);
452+
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {
453+
registry.fill(HIST("hSignalPt"), jet.pt(), weight);
454+
registry.fill(HIST("hSignalPtHard"), jet.pt(), ptTT / pTHat, weight);
455+
}
456+
}
457+
if (!isSigCol) {
458+
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {
459+
registry.fill(HIST("hDeltaRpTPartReference"), jet.pt(), dR, weight);
460+
registry.fill(HIST("hDeltaRPartReference"), dR, weight);
461+
}
462+
registry.fill(HIST("hDeltaRpTDPhiReferencePart"), jet.pt(), dphi, dR, weight);
463+
registry.fill(HIST("hReferencePtDPhi"), dphi, jet.pt(), weight);
464+
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {
465+
registry.fill(HIST("hReferencePt"), jet.pt(), weight);
466+
registry.fill(HIST("hReferencePtHard"), jet.pt(), ptTT / pTHat, weight);
467+
}
468+
}
469+
}
470+
}
471+
registry.fill(HIST("hPartvsJets"), leadingPartPt, leadingJetPt, pTHat, weight);
472+
}
473+
474+
template <typename T, typename U, typename V>
475+
void fillMCPHistogramsWithMatchedTracks(T const& jets, U const& particles, V const& tracks, float weight = 1.0, float pTHat = 999.0)
476+
{
477+
bool isSigCol;
478+
std::vector<double> phiTTAr;
479+
std::vector<double> ptTTAr;
480+
double phiTT = 0;
481+
double ptTT = 0;
482+
int trigNumber = 0;
483+
int nTT = 0;
484+
double leadingPartPt = 0;
485+
double leadingJetPt = 0;
486+
float dice = rand->Rndm();
487+
if (dice < fracSig)
488+
isSigCol = true;
489+
else
490+
isSigCol = false;
491+
492+
for (const auto& track : tracks) {
493+
if (!track.has_mcParticle()) {
494+
continue;
495+
}
496+
auto particle = track.template mcParticle_as<U>();
497+
if (particle.pt() > leadingPartPt) {
498+
leadingPartPt = particle.pt();
499+
}
500+
if (particle.pt() > pTHatTrackMaxMCD * pTHat) {
501+
if (outlierRejectEvent) {
502+
return;
503+
} else {
504+
continue;
505+
}
506+
}
507+
auto pdgParticle = pdg->GetParticle(particle.pdgCode());
508+
if (!pdgParticle) {
509+
continue;
510+
}
511+
if ((pdgParticle->Charge() == 0.0) || (!particle.isPhysicalPrimary())) {
512+
continue;
513+
}
514+
if (isSigCol && particle.pt() < ptTTsigMax && particle.pt() > ptTTsigMin && track.pt() < ptTTsigMax && track.pt() > ptTTsigMin) {
515+
phiTTAr.push_back(particle.phi());
516+
ptTTAr.push_back(particle.pt());
517+
nTT++;
518+
registry.fill(HIST("hSignalTriggers"), particle.pt(), weight);
519+
}
520+
if (!isSigCol && particle.pt() < ptTTrefMax && particle.pt() > ptTTrefMin && track.pt() < ptTTrefMax && track.pt() > ptTTrefMin) {
521+
phiTTAr.push_back(particle.phi());
522+
ptTTAr.push_back(particle.pt());
523+
nTT++;
524+
registry.fill(HIST("hReferenceTriggers"), particle.pt(), weight);
397525
}
398526
registry.fill(HIST("hPtPart"), particle.pt(), weight);
399527
registry.fill(HIST("hEtaPart"), particle.eta(), weight);
@@ -537,10 +665,19 @@ struct JetHadronRecoil {
537665
}
538666
}
539667
if (track.pt() < ptTTsigMax && track.pt() > ptTTsigMin) {
540-
phiTTAr.push_back(track.phi());
541-
nTT++;
542668
auto particle = track.template mcParticle_as<Y>();
543-
phiTTArPart.push_back(particle.phi());
669+
auto pdgParticle = pdg->GetParticle(particle.pdgCode());
670+
if (!pdgParticle) {
671+
continue;
672+
}
673+
if ((pdgParticle->Charge() == 0.0) || (!particle.isPhysicalPrimary())) {
674+
continue;
675+
}
676+
if (particle.pt() < ptTTsigMax && particle.pt() > ptTTsigMin) {
677+
nTT++;
678+
phiTTAr.push_back(track.phi());
679+
phiTTArPart.push_back(particle.phi());
680+
}
544681
}
545682
}
546683

@@ -761,6 +898,25 @@ struct JetHadronRecoil {
761898
}
762899
PROCESS_SWITCH(JetHadronRecoil, processMCPWeighted, "process MC particle level with event weights", false);
763900

901+
void processMCPWeightedWithMatchedTracks(aod::JetMcCollision const& collision,
902+
soa::Filtered<soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents>> const& jets,
903+
soa::Filtered<aod::JetParticles> const& particles,
904+
soa::Filtered<aod::JetTracksMCD> const& tracks)
905+
{
906+
if (std::abs(collision.posZ()) > vertexZCut) {
907+
return;
908+
}
909+
if (skipMBGapEvents && collision.subGeneratorId() == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) {
910+
return;
911+
}
912+
if (collision.ptHard() < pTHatMinEvent) {
913+
return;
914+
}
915+
registry.fill(HIST("hZvtxSelected"), collision.posZ(), collision.weight());
916+
fillMCPHistogramsWithMatchedTracks(jets, particles, tracks, collision.weight(), collision.ptHard());
917+
}
918+
PROCESS_SWITCH(JetHadronRecoil, processMCPWeightedWithMatchedTracks, "process MC particle level with event weights - only triggers matched with detector level tracks", false);
919+
764920
void processJetsMCPMCDMatched(soa::Filtered<soa::Join<aod::JetCollisionsMCD, aod::JMcCollisionLbs>>::iterator const& collision,
765921
soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets>> const& mcdjets,
766922
aod::JetTracks const& tracks,

0 commit comments

Comments
 (0)