Skip to content

Commit a4cbe1b

Browse files
dajones2dajones2
committed
adding separate mcd function
1 parent 7ebc2af commit a4cbe1b

File tree

1 file changed

+148
-14
lines changed

1 file changed

+148
-14
lines changed

‎PWGJE/Tasks/jetHadronRecoil.cxx‎

Lines changed: 148 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -263,13 +263,147 @@ struct JetHadronRecoil {
263263
registry.fill(HIST("hPhiTrack"), track.phi(), weight);
264264
registry.fill(HIST("hTrack3D"), track.pt(), track.eta(), track.phi(), weight);
265265
registry.fill(HIST("hPtTrackPtHard"), track.pt() / pTHat, track.pt(), weight);
266-
if (isMC) {
267-
if (track.has_mcParticle()) {
268-
registry.fill(HIST("hPtTrackMatched"), track.pt(), weight);
269-
auto particle = track.template mcParticle_as<U>();
270-
if (track.collisionId() == particle.collisionId()) {
271-
registry.fill(HIST("hPtTrackMatchedToCollisions"), track.pt(), weight);
266+
}
267+
if (nTT > 0) {
268+
trigNumber = rand->Integer(nTT);
269+
phiTT = phiTTAr[trigNumber];
270+
ptTT = ptTTAr[trigNumber];
271+
if (isSigCol) {
272+
registry.fill(HIST("hNtrig"), 1.5, weight);
273+
registry.fill(HIST("hSigEventTriggers"), nTT, weight);
274+
registry.fill(HIST("hRhoSignal"), rho, weight);
275+
registry.fill(HIST("hSignalTriggersPtHard"), ptTT / pTHat, weight);
276+
}
277+
if (!isSigCol) {
278+
registry.fill(HIST("hNtrig"), 0.5, weight);
279+
registry.fill(HIST("hRefEventTriggers"), nTT, weight);
280+
registry.fill(HIST("hRhoReference"), rhoReference, weight);
281+
for (double shift = 0.0; shift <= 2.0; shift += 0.1) {
282+
registry.fill(HIST("hRhoReferenceShift"), rho + shift, shift, weight);
283+
}
284+
registry.fill(HIST("hReferenceTriggersPtHard"), ptTT / pTHat, weight);
285+
}
286+
}
287+
for (const auto& jet : jets) {
288+
if (jet.pt() > leadingJetPt) {
289+
leadingJetPt = jet.pt();
290+
}
291+
if (jet.pt() > pTHatMaxMCD * pTHat) {
292+
if (outlierRejectEvent) {
293+
return;
294+
} else {
295+
continue;
296+
}
297+
}
298+
for (const auto& constituent : jet.template tracks_as<U>()) {
299+
if (constituent.pt() > leadingPT) {
300+
leadingPT = constituent.pt();
301+
}
302+
registry.fill(HIST("hConstituents3D"), constituent.pt(), constituent.eta(), constituent.phi());
303+
}
304+
if (leadingPT > maxLeadingTrackPt) {
305+
continue;
306+
}
307+
registry.fill(HIST("hJetPt"), jet.pt() - (rho * jet.area()), weight);
308+
registry.fill(HIST("hJetEta"), jet.eta(), weight);
309+
registry.fill(HIST("hJetPhi"), jet.phi(), weight);
310+
registry.fill(HIST("hJet3D"), jet.pt() - (rho * jet.area()), jet.eta(), jet.phi(), weight);
311+
312+
if (nTT > 0) {
313+
float dphi = RecoDecay::constrainAngle(jet.phi() - phiTT);
314+
double dR = getWTAaxisDifference(jet, tracks);
315+
if (isSigCol) {
316+
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {
317+
registry.fill(HIST("hDeltaRpTSignal"), jet.pt() - (rho * jet.area()), dR, weight);
318+
registry.fill(HIST("hDeltaRSignal"), dR, weight);
319+
}
320+
registry.fill(HIST("hDeltaRpTDPhiSignal"), jet.pt() - (rho * jet.area()), dphi, dR, weight);
321+
registry.fill(HIST("hSignalPtDPhi"), dphi, jet.pt() - (rho * jet.area()), weight);
322+
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {
323+
registry.fill(HIST("hSignalPt"), jet.pt() - (rho * jet.area()), weight);
324+
registry.fill(HIST("hSignalPtHard"), jet.pt() - (rho * jet.area()), ptTT / pTHat, weight);
325+
}
326+
}
327+
if (!isSigCol) {
328+
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {
329+
registry.fill(HIST("hDeltaRpTReference"), jet.pt() - (rhoReference * jet.area()), dR, weight);
330+
registry.fill(HIST("hDeltaRReference"), dR, weight);
331+
}
332+
registry.fill(HIST("hDeltaRpTDPhiReference"), jet.pt() - (rhoReference * jet.area()), dphi, dR, weight);
333+
for (double shift = 0.0; shift <= 2.0; shift += 0.1) {
334+
registry.fill(HIST("hDeltaRpTDPhiReferenceShifts"), jet.pt() - ((rho + shift) * jet.area()), dphi, dR, shift, weight);
272335
}
336+
registry.fill(HIST("hReferencePtDPhi"), dphi, jet.pt() - (rhoReference * jet.area()), weight);
337+
for (double shift = 0.0; shift <= 2.0; shift += 0.1) {
338+
registry.fill(HIST("hReferencePtDPhiShifts"), dphi, jet.pt() - ((rho + shift) * jet.area()), shift, weight);
339+
}
340+
if (std::abs(dphi - o2::constants::math::PI) < 0.6) {
341+
registry.fill(HIST("hReferencePt"), jet.pt() - (rhoReference * jet.area()), weight);
342+
registry.fill(HIST("hReferencePtHard"), jet.pt() - (rhoReference * jet.area()), ptTT / pTHat, weight);
343+
}
344+
}
345+
}
346+
}
347+
registry.fill(HIST("hTracksvsJets"), leadingTrackPt, leadingJetPt, pTHat, weight);
348+
}
349+
350+
template <typename T, typename U>
351+
void fillHistogramsMCD(T const& jets, U const& tracks, float weight = 1.0, float rho = 0.0, float pTHat = 999.0)
352+
{
353+
bool isSigCol;
354+
std::vector<double> phiTTAr;
355+
std::vector<double> ptTTAr;
356+
double phiTT = 0;
357+
double ptTT = 0;
358+
int trigNumber = 0;
359+
int nTT = 0;
360+
double leadingPT = 0;
361+
double leadingTrackPt = 0;
362+
double leadingJetPt = 0;
363+
float rhoReference = rho + rhoReferenceShift;
364+
365+
float dice = rand->Rndm();
366+
if (dice < fracSig)
367+
isSigCol = true;
368+
else
369+
isSigCol = false;
370+
371+
for (const auto& track : tracks) {
372+
if (!jetderiveddatautilities::selectTrack(track, trackSelection)) {
373+
continue;
374+
}
375+
if (track.pt() > leadingTrackPt) {
376+
leadingTrackPt = track.pt();
377+
}
378+
if (track.pt() > pTHatTrackMaxMCD * pTHat) {
379+
if (outlierRejectEvent) {
380+
return;
381+
} else {
382+
continue;
383+
}
384+
}
385+
if (isSigCol && track.pt() < ptTTsigMax && track.pt() > ptTTsigMin) {
386+
phiTTAr.push_back(track.phi());
387+
ptTTAr.push_back(track.pt());
388+
registry.fill(HIST("hSignalTriggers"), track.pt(), weight);
389+
nTT++;
390+
}
391+
if (!isSigCol && track.pt() < ptTTrefMax && track.pt() > ptTTrefMin) {
392+
phiTTAr.push_back(track.phi());
393+
ptTTAr.push_back(track.pt());
394+
registry.fill(HIST("hReferenceTriggers"), track.pt(), weight);
395+
nTT++;
396+
}
397+
registry.fill(HIST("hPtTrack"), track.pt(), weight);
398+
registry.fill(HIST("hEtaTrack"), track.eta(), weight);
399+
registry.fill(HIST("hPhiTrack"), track.phi(), weight);
400+
registry.fill(HIST("hTrack3D"), track.pt(), track.eta(), track.phi(), weight);
401+
registry.fill(HIST("hPtTrackPtHard"), track.pt() / pTHat, track.pt(), weight);
402+
if (track.has_mcParticle()) {
403+
registry.fill(HIST("hPtTrackMatched"), track.pt(), weight);
404+
auto particle = track.template mcParticle_as<U>();
405+
if (track.collisionId() == particle.collisionId()) {
406+
registry.fill(HIST("hPtTrackMatchedToCollisions"), track.pt(), weight);
273407
}
274408
}
275409
}
@@ -772,7 +906,7 @@ struct JetHadronRecoil {
772906
void processMCD(soa::Filtered<soa::Join<aod::JetCollisions, aod::JMcCollisionLbs>>::iterator const& collision,
773907
aod::JMcCollisions const&,
774908
soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents>> const& jets,
775-
soa::Filtered<aod::JetTracks> const& tracks)
909+
soa::Filtered<aod::JetTracksMCD> const& tracks)
776910
{
777911
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
778912
return;
@@ -790,14 +924,14 @@ struct JetHadronRecoil {
790924
return;
791925
}
792926
registry.fill(HIST("hZvtxSelected"), collision.posZ());
793-
fillHistograms(jets, tracks, 1.0, 0.0, collision.mcCollision().ptHard());
927+
fillHistogramsMCD(jets, tracks, 1.0, 0.0, collision.mcCollision().ptHard());
794928
}
795929
PROCESS_SWITCH(JetHadronRecoil, processMCD, "process MC detector level", false);
796930

797931
void processMCDWithRhoSubtraction(soa::Filtered<soa::Join<aod::JetCollisions, aod::BkgChargedRhos, aod::JMcCollisionLbs>>::iterator const& collision,
798932
aod::JMcCollisions const&,
799933
soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents>> const& jets,
800-
soa::Filtered<aod::JetTracks> const& tracks)
934+
soa::Filtered<aod::JetTracksMCD> const& tracks)
801935
{
802936
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
803937
return;
@@ -815,14 +949,14 @@ struct JetHadronRecoil {
815949
return;
816950
}
817951
registry.fill(HIST("hZvtxSelected"), collision.posZ());
818-
fillHistograms(jets, tracks, 1.0, collision.rho(), collision.mcCollision().ptHard());
952+
fillHistogramsMCD(jets, tracks, 1.0, collision.rho(), collision.mcCollision().ptHard());
819953
}
820954
PROCESS_SWITCH(JetHadronRecoil, processMCDWithRhoSubtraction, "process MC detector level with rho subtraction", false);
821955

822956
void processMCDWeighted(soa::Filtered<soa::Join<aod::JetCollisions, aod::JMcCollisionLbs>>::iterator const& collision,
823957
aod::JMcCollisions const&,
824958
soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents>> const& jets,
825-
soa::Filtered<aod::JetTracks> const& tracks)
959+
soa::Filtered<aod::JetTracksMCD> const& tracks)
826960
{
827961
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
828962
return;
@@ -840,14 +974,14 @@ struct JetHadronRecoil {
840974
return;
841975
}
842976
registry.fill(HIST("hZvtxSelected"), collision.posZ(), collision.mcCollision().weight());
843-
fillHistograms(jets, tracks, collision.mcCollision().weight(), 0.0, collision.mcCollision().ptHard());
977+
fillHistogramsMCD(jets, tracks, collision.mcCollision().weight(), 0.0, collision.mcCollision().ptHard());
844978
}
845979
PROCESS_SWITCH(JetHadronRecoil, processMCDWeighted, "process MC detector level with event weights", false);
846980

847981
void processMCDWeightedWithRhoSubtraction(soa::Filtered<soa::Join<aod::JetCollisions, aod::JMcCollisionLbs, aod::BkgChargedRhos>>::iterator const& collision,
848982
aod::JMcCollisions const&,
849983
soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents>> const& jets,
850-
soa::Filtered<aod::JetTracks> const& tracks)
984+
soa::Filtered<aod::JetTracksMCD> const& tracks)
851985
{
852986
if (!jetderiveddatautilities::selectCollision(collision, eventSelectionBits)) {
853987
return;
@@ -865,7 +999,7 @@ struct JetHadronRecoil {
865999
return;
8661000
}
8671001
registry.fill(HIST("hZvtxSelected"), collision.posZ(), collision.mcCollision().weight());
868-
fillHistograms(jets, tracks, collision.mcCollision().weight(), collision.rho(), collision.mcCollision().ptHard());
1002+
fillHistogramsMCD(jets, tracks, collision.mcCollision().weight(), collision.rho(), collision.mcCollision().ptHard());
8691003
}
8701004
PROCESS_SWITCH(JetHadronRecoil, processMCDWeightedWithRhoSubtraction, "process MC detector level with event weights and rho subtraction", false);
8711005

0 commit comments

Comments
 (0)