@@ -349,7 +349,7 @@ struct JetHadronRecoil {
349349 }
350350
351351 template <typename T, typename U>
352- void fillHistogramsMCD (T const & jets, U const & tracks, float weight = 1.0 , float rho = 0.0 , float pTHat = 999.0 , auto collisionID)
352+ void fillHistogramsMCD (T const & jets, U const & tracks, float weight = 1.0 , float rho = 0.0 , float pTHat = 999.0 , auto collisionID = 0 )
353353 {
354354 bool isSigCol;
355355 std::vector<double > phiTTAr;
@@ -403,7 +403,7 @@ struct JetHadronRecoil {
403403 if (track.has_mcParticle ()) {
404404 registry.fill (HIST (" hPtTrackMatched" ), track.pt (), weight);
405405 auto mcParticle = track.mcParticle ();
406- if (mcParticle.mcCollisionId () == collision. mcCollisionId () ) {
406+ if (mcParticle.mcCollisionId () == collisionID ) {
407407 registry.fill (HIST (" hPtTrackMatchedToCollisions" ), track.pt (), weight);
408408 }
409409 }
@@ -613,219 +613,42 @@ struct JetHadronRecoil {
613613 registry.fill (HIST (" hPartvsJets" ), leadingPartPt, leadingJetPt, pTHat, weight);
614614 }
615615
616- template <typename T, typename U, typename V >
617- void fillMCPHistogramsWithMatchedTracks (T const & jets , U const & particles, V const & tracks, float weight = 1.0 , float pTHat = 999.0 )
616+ template <typename T, typename U, typename X, typename Y >
617+ void fillRecoilJetMatchedHistograms (T const & jetsBase , U const &, X const & tracks, Y const & particles, float weight = 1.0 , float rho = 0 .0 , float pTHat = 999.0 )
618618 {
619- bool isSigCol;
620619 std::vector<double > phiTTAr;
621- std::vector<double > phiTTArTrack;
622- std::vector<double > ptTTAr;
620+ std::vector<double > phiTTArPart;
623621 double phiTT = 0 ;
624- double phiTTTrack = 0 ;
625- double ptTT = 0 ;
622+ double phiTTPart = 0 ;
626623 int trigNumber = 0 ;
627624 int nTT = 0 ;
628- double leadingPartPt = 0 ;
629- double leadingJetPt = 0 ;
630- float dice = rand->Rndm ();
631- if (dice < fracSig)
632- isSigCol = true ;
633- else
634- isSigCol = false ;
635625
636626 for (const auto & track : tracks) {
637627 if (!track.has_mcParticle ()) {
638628 continue ;
639629 }
640- auto particle = track.template mcParticle_as <U>();
641- if (particle.pt () > leadingPartPt) {
642- leadingPartPt = particle.pt ();
630+ if (!jetderiveddatautilities::selectTrack (track, trackSelection)) {
631+ continue ;
643632 }
644- if (particle .pt () > pTHatTrackMaxMCD * pTHat) {
633+ if (track .pt () > pTHatTrackMaxMCD * pTHat) {
645634 if (outlierRejectEvent) {
646635 return ;
647636 } else {
648637 continue ;
649638 }
650639 }
640+ auto particle = track.template mcParticle_as <Y>();
651641 auto pdgParticle = pdg->GetParticle (particle.pdgCode ());
652642 if (!pdgParticle) {
653643 continue ;
654644 }
655645 if ((pdgParticle->Charge () == 0.0 ) || (!particle.isPhysicalPrimary ())) {
656646 continue ;
657647 }
658- if (isSigCol && particle.pt () < ptTTsigMax && particle.pt () > ptTTsigMin && track.pt () < ptTTsigMax && track.pt () > ptTTsigMin) {
659- phiTTAr.push_back (particle.phi ());
660- phiTTArTrack.push_back (track.phi ());
661- ptTTAr.push_back (particle.pt ());
648+ if (particle.pt () < ptTTsigMax && particle.pt () > ptTTsigMin) {
662649 nTT++;
663- registry.fill (HIST (" hSignalTriggers" ), particle.pt (), weight);
664- }
665- if (!isSigCol && particle.pt () < ptTTrefMax && particle.pt () > ptTTrefMin && track.pt () < ptTTrefMax && track.pt () > ptTTrefMin) {
666- phiTTAr.push_back (particle.phi ());
667- ptTTAr.push_back (particle.pt ());
668- phiTTArTrack.push_back (track.phi ());
669- nTT++;
670- registry.fill (HIST (" hReferenceTriggers" ), particle.pt (), weight);
671- }
672- registry.fill (HIST (" hPtPart" ), particle.pt (), weight);
673- registry.fill (HIST (" hEtaPart" ), particle.eta (), weight);
674- registry.fill (HIST (" hPhiPart" ), particle.phi (), weight);
675- registry.fill (HIST (" hPart3D" ), particle.pt (), particle.eta (), particle.phi (), weight);
676- registry.fill (HIST (" hPtPartPtHard" ), particle.pt (), particle.pt () / pTHat, weight);
677- }
678-
679- if (nTT > 0 ) {
680- trigNumber = rand->Integer (nTT);
681- phiTT = phiTTAr[trigNumber];
682- phiTTTrack = phiTTArTrack[trigNumber];
683- ptTT = ptTTAr[trigNumber];
684- if (isSigCol) {
685- registry.fill (HIST (" hNtrig" ), 1.5 , weight);
686- registry.fill (HIST (" hSigEventTriggers" ), nTT, weight);
687- registry.fill (HIST (" hSignalTriggersPtHard" ), ptTT / pTHat, weight);
688- }
689- if (!isSigCol) {
690- registry.fill (HIST (" hNtrig" ), 0.5 , weight);
691- registry.fill (HIST (" hRefEventTriggers" ), nTT, weight);
692- registry.fill (HIST (" hReferenceTriggersPtHard" ), ptTT / pTHat, weight);
693- }
694- }
695-
696- for (const auto & jet : jets) {
697- if (jet.pt () > leadingJetPt) {
698- leadingJetPt = jet.pt ();
699- }
700- if (jet.pt () > pTHatMaxMCP * pTHat) {
701- if (outlierRejectEvent) {
702- return ;
703- } else {
704- continue ;
705- }
706- }
707- for (const auto & constituent : jet.template tracks_as <U>()) {
708- registry.fill (HIST (" hConstituents3D" ), constituent.pt (), constituent.eta (), constituent.phi ());
709- }
710- registry.fill (HIST (" hJetPt" ), jet.pt (), weight);
711- registry.fill (HIST (" hJetEta" ), jet.eta (), weight);
712- registry.fill (HIST (" hJetPhi" ), jet.phi (), weight);
713- registry.fill (HIST (" hJet3D" ), jet.pt (), jet.eta (), jet.phi (), weight);
714-
715- if (nTT > 0 ) {
716- float dphi = RecoDecay::constrainAngle (jet.phi () - phiTT);
717- float dphiTrack = RecoDecay::constrainAngle (jet.phi () - phiTTTrack);
718- double dR = getWTAaxisDifference (jet, particles);
719- if (isSigCol) {
720- if (std::abs (dphi - o2::constants::math::PI) < 0.6 ) {
721- registry.fill (HIST (" hDeltaRpTSignalPart" ), jet.pt (), dR, weight);
722- registry.fill (HIST (" hDeltaRSignalPart" ), dR, weight);
723- }
724- registry.fill (HIST (" hDeltaRpTDPhiSignalPart" ), jet.pt (), dphi, dR, weight);
725- registry.fill (HIST (" hSignalPtDPhi" ), dphi, jet.pt (), weight);
726- if (std::abs (dphi - o2::constants::math::PI) < 0.6 && std::abs (dphiTrack - o2::constants::math::PI) < 0.6 ) {
727- registry.fill (HIST (" hSignalPt" ), jet.pt (), weight);
728- registry.fill (HIST (" hSignalPtHard" ), jet.pt (), ptTT / pTHat, weight);
729- }
730- }
731- if (!isSigCol) {
732- if (std::abs (dphi - o2::constants::math::PI) < 0.6 ) {
733- registry.fill (HIST (" hDeltaRpTPartReference" ), jet.pt (), dR, weight);
734- registry.fill (HIST (" hDeltaRPartReference" ), dR, weight);
735- }
736- registry.fill (HIST (" hDeltaRpTDPhiReferencePart" ), jet.pt (), dphi, dR, weight);
737- registry.fill (HIST (" hReferencePtDPhi" ), dphi, jet.pt (), weight);
738- if (std::abs (dphi - o2::constants::math::PI) < 0.6 ) {
739- registry.fill (HIST (" hReferencePt" ), jet.pt (), weight);
740- registry.fill (HIST (" hReferencePtHard" ), jet.pt (), ptTT / pTHat, weight);
741- }
742- }
743- }
744- }
745- registry.fill (HIST (" hPartvsJets" ), leadingPartPt, leadingJetPt, pTHat, weight);
746- }
747-
748- template <typename T, typename U, typename X, typename Y>
749- 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 )
750- {
751- for (const auto & jetBase : jetsBase) {
752- double dR = 0 ;
753- double dRp = 0 ;
754-
755- if (jetBase.pt () > pTHatMaxMCD * pTHat) {
756- if (outlierRejectEvent) {
757- return ;
758- } else {
759- continue ;
760- }
761- }
762-
763- dR = getWTAaxisDifference (jetBase, tracks);
764-
765- if (jetBase.has_matchedJetGeo ()) {
766- for (const auto & jetTag : jetBase.template matchedJetGeo_as <std::decay_t <U>>()) {
767- if (jetTag.pt () > pTHatMaxMCP * pTHat) {
768- if (outlierRejectEvent) {
769- return ;
770- } else {
771- continue ;
772- }
773- }
774-
775- dRp = getWTAaxisDifference (jetTag, particles);
776-
777- registry.fill (HIST (" hPtMatched" ), jetBase.pt () - (rho * jetBase.area ()), jetTag.pt (), weight);
778- registry.fill (HIST (" hPhiMatched" ), jetBase.phi (), jetTag.phi (), weight);
779- registry.fill (HIST (" hPtResolution" ), jetTag.pt (), (jetTag.pt () - (jetBase.pt () - (rho * jetBase.area ()))) / jetTag.pt (), weight);
780- registry.fill (HIST (" hPhiResolution" ), jetTag.pt (), jetTag.phi () - jetBase.phi (), weight);
781- registry.fill (HIST (" hDeltaRMatched" ), dR, dRp, weight);
782- registry.fill (HIST (" hDeltaRResolution" ), jetTag.pt (), dRp - dR, weight);
783- registry.fill (HIST (" hFullMatching" ), jetBase.pt () - (rho * jetBase.area ()), jetTag.pt (), jetBase.phi (), jetTag.phi (), dR, dRp, weight);
784- registry.fill (HIST (" hPtMatched1d" ), jetTag.pt (), weight);
785- registry.fill (HIST (" hDeltaRMatched1d" ), dRp, weight);
786- }
787- }
788- }
789- }
790-
791- template <typename T, typename U, typename X, typename Y>
792- void fillRecoilJetMatchedHistograms (T const & jetsBase, U const &, X const & tracks, Y const & particles, float weight = 1.0 , float rho = 0.0 , float pTHat = 999.0 )
793- {
794- std::vector<double > phiTTAr;
795- std::vector<double > phiTTArPart;
796- double phiTT = 0 ;
797- double phiTTPart = 0 ;
798- int trigNumber = 0 ;
799- int nTT = 0 ;
800-
801- for (const auto & track : tracks) {
802- if (!track.has_mcParticle ()) {
803- continue ;
804- }
805- if (!jetderiveddatautilities::selectTrack (track, trackSelection)) {
806- continue ;
807- }
808- if (track.pt () > pTHatTrackMaxMCD * pTHat) {
809- if (outlierRejectEvent) {
810- return ;
811- } else {
812- continue ;
813- }
814- }
815- if (track.pt () < ptTTsigMax && track.pt () > ptTTsigMin) {
816- auto particle = track.template mcParticle_as <Y>();
817- auto pdgParticle = pdg->GetParticle (particle.pdgCode ());
818- if (!pdgParticle) {
819- continue ;
820- }
821- if ((pdgParticle->Charge () == 0.0 ) || (!particle.isPhysicalPrimary ())) {
822- continue ;
823- }
824- if (particle.pt () < ptTTsigMax && particle.pt () > ptTTsigMin) {
825- nTT++;
826- phiTTAr.push_back (track.phi ());
827- phiTTArPart.push_back (particle.phi ());
828- }
650+ phiTTAr.push_back (track.phi ());
651+ phiTTArPart.push_back (particle.phi ());
829652 }
830653 }
831654
@@ -868,7 +691,7 @@ struct JetHadronRecoil {
868691 registry.fill (HIST (" hPhiMatched2d" ), jetTag.phi (), jetTag.pt (), weight);
869692 registry.fill (HIST (" hPhiResolution" ), jetTag.pt (), dphip - dphi, weight);
870693 registry.fill (HIST (" hFullMatching" ), jetBase.pt () - (rho * jetBase.area ()), jetTag.pt (), dphi, dphip, dR, dRp, weight);
871- if ((std::abs (dphi - o2::constants::math::PI) < 0.6 ) || ( std::abs ( dphip - o2::constants::math::PI) < 0.6 )) {
694+ if ((std::abs (dphip - o2::constants::math::PI) < 0.6 )) {
872695 registry.fill (HIST (" hPtMatched1d" ), jetTag.pt (), weight);
873696 registry.fill (HIST (" hDeltaRMatched1d" ), dRp, weight);
874697 registry.fill (HIST (" hPtMatched" ), jetBase.pt () - (rho * jetBase.area ()), jetTag.pt (), weight);
@@ -1047,25 +870,6 @@ struct JetHadronRecoil {
1047870 }
1048871 PROCESS_SWITCH (JetHadronRecoil, processMCPWeighted, " process MC particle level with event weights" , false );
1049872
1050- void processMCPWeightedWithMatchedTracks (aod::JetMcCollision const & collision,
1051- soa::Filtered<soa::Join<aod::ChargedMCParticleLevelJets, aod::ChargedMCParticleLevelJetConstituents>> const & jets,
1052- soa::Filtered<aod::JetParticles> const & particles,
1053- soa::Filtered<aod::JetTracksMCD> const & tracks)
1054- {
1055- if (std::abs (collision.posZ ()) > vertexZCut) {
1056- return ;
1057- }
1058- if (skipMBGapEvents && collision.subGeneratorId () == jetderiveddatautilities::JCollisionSubGeneratorId::mbGap) {
1059- return ;
1060- }
1061- if (collision.ptHard () < pTHatMinEvent) {
1062- return ;
1063- }
1064- registry.fill (HIST (" hZvtxSelected" ), collision.posZ (), collision.weight ());
1065- fillMCPHistogramsWithMatchedTracks (jets, particles, tracks, collision.weight (), collision.ptHard ());
1066- }
1067- PROCESS_SWITCH (JetHadronRecoil, processMCPWeightedWithMatchedTracks, " process MC particle level with event weights - only triggers matched with detector level tracks" , false );
1068-
1069873 void processJetsMCPMCDMatched (soa::Filtered<soa::Join<aod::JetCollisionsMCD, aod::JMcCollisionLbs>>::iterator const & collision,
1070874 soa::Filtered<soa::Join<aod::ChargedMCDetectorLevelJets, aod::ChargedMCDetectorLevelJetConstituents, aod::ChargedMCDetectorLevelJetsMatchedToChargedMCParticleLevelJets>> const & mcdjets,
1071875 aod::JetTracks const & tracks,
0 commit comments