3232#include < Framework/InitContext.h>
3333#include < Framework/runDataProcessing.h>
3434
35- #include < Math/Vector4D.h>
3635#include < TH1.h>
37- #include < TVector3.h>
3836
37+ #include < array>
3938#include < cmath>
4039#include < string>
4140#include < vector>
@@ -44,6 +43,12 @@ using namespace o2;
4443using namespace o2 ::framework;
4544using namespace o2 ::framework::expressions;
4645
46+ namespace
47+ {
48+ constexpr int kOriginMcPrompt = 1 ; // MC origin flag: prompt
49+ constexpr int kOriginMcNonPrompt = 2 ; // MC origin flag: non-prompt
50+ } // namespace
51+
4752namespace o2 ::aod
4853{
4954
@@ -190,25 +195,14 @@ struct JetHFAngularityTask {
190195
191196 // DATA
192197
193- using D0CandidatesData = soa::Join<aod::HfD0Bases,
194- aod::HfD0Pars,
195- aod::HfD0ParEs,
196- aod::HfD0Sels,
197- aod::HfD0Mls,
198- aod::JD0Ids>;
199-
198+ using D0CandidatesData = aod::CandidatesD0Data;
199+
200200 using D0DataJets = soa::Join<aod::D0ChargedJets,
201201 aod::D0ChargedJetConstituents>;
202202
203203 // MC-DETECTOR LEVEL(MCD)
204204
205- using D0CandidatesMCD = soa::Join<aod::HfD0Bases,
206- aod::HfD0Pars,
207- aod::HfD0ParEs,
208- aod::HfD0Sels,
209- aod::HfD0Mls,
210- aod::HfD0Mcs,
211- aod::JD0Ids>;
205+ using D0CandidatesMCD = aod::CandidatesD0MCD;
212206
213207 using D0MCDJets = soa::Join<aod::D0ChargedMCDetectorLevelJets,
214208 aod::D0ChargedMCDetectorLevelJetConstituents>;
@@ -219,8 +213,7 @@ struct JetHFAngularityTask {
219213
220214 // MC PARTICLE LEVEL (MCP)
221215
222- using D0CandidatesMCP = soa::Join<aod::HfD0PBases,
223- aod::JD0PIds>;
216+ using D0CandidatesMCP = aod::CandidatesD0MCP;
224217
225218 using D0MCPJets = soa::Join<aod::D0ChargedMCParticleLevelJets,
226219 aod::D0ChargedMCParticleLevelJetConstituents>;
@@ -393,10 +386,8 @@ struct JetHFAngularityTask {
393386
394387 void init (InitContext const &)
395388 {
396- eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits (
397- static_cast <std::string>(eventSelections));
398- trackSelection = jetderiveddatautilities::initialiseTrackSelection (
399- static_cast <std::string>(trackSelections));
389+ eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits (eventSelections.value );
390+ trackSelection = jetderiveddatautilities::initialiseTrackSelection (trackSelections.value );
400391
401392 massD0MCP = jetcandidateutilities::getTablePDGMass<D0CandidatesMCP>();
402393
@@ -486,57 +477,45 @@ struct JetHFAngularityTask {
486477 }
487478
488479 // Helper: jet invariant mass
489-
480+
490481 template <typename TRACKS , typename CANDIDATES >
491482 float computeJetMass (TRACKS const & tracks,
492483 CANDIDATES const & candidates,
493484 double candMass = -1 .)
494485 {
495- double sumPx = 0 ., sumPy = 0 ., sumPz = 0 ., sumE = 0 .;
486+ std::array<double , 3 > momTotal{0 ., 0 ., 0 .};
487+ double energyTot = 0 .;
496488
497489 for (auto const & trk : tracks) {
498- const double px = trk.pt () * std::cos (trk.phi ());
499- const double py = trk.pt () * std::sin (trk.phi ());
500- const double pz = trk.pt () * std::sinh (trk.eta ());
501- const double p = std::sqrt (px * px + py * py + pz * pz);
502-
503- sumPx += px;
504- sumPy += py;
505- sumPz += pz;
506- sumE += p;
490+ const std::array<double , 3 > mom{trk.px (), trk.py (), trk.pz ()};
491+ momTotal[0 ] += mom[0 ];
492+ momTotal[1 ] += mom[1 ];
493+ momTotal[2 ] += mom[2 ];
494+ energyTot += RecoDecay::e (mom, 0 .); // massless approximation for ordinary tracks
507495 }
508496
509497 for (auto const & cand : candidates) {
510-
511- const double px = cand.px ();
512- const double py = cand.py ();
513- const double pz = cand.pz ();
498+ const std::array<double , 3 > mom{cand.px (), cand.py (), cand.pz ()};
514499
515500 double m;
516501 if (candMass > 0 .) {
517502 m = candMass;
503+ } else if constexpr (requires { cand.m (); }) {
504+ m = cand.m (); // DATA / MCD: reconstructed invariant mass
518505 } else {
519-
520- if constexpr (requires { cand.m (); }) {
521- m = cand.m ();
522- } else {
523- m = candMass;
524- }
506+ m = 0 .;
525507 }
526508
527- const double p = std::sqrt (px * px + py * py + pz * pz);
528- const double e = std::sqrt (p * p + m * m);
529-
530- sumPx += px;
531- sumPy += py;
532- sumPz += pz;
533- sumE += e;
509+ momTotal[0 ] += mom[0 ];
510+ momTotal[1 ] += mom[1 ];
511+ momTotal[2 ] += mom[2 ];
512+ energyTot += RecoDecay::e (mom, m);
534513 }
535514
536- const double m2 = sumE * sumE - (sumPx * sumPx + sumPy * sumPy + sumPz * sumPz);
537-
538- return (m2 > 0 .) ? static_cast <float >(std::sqrt (m2)) : 0 .f ;
515+ const double mass2 = RecoDecay::m2 (momTotal, energyTot);
516+ return (mass2 > 0 .) ? static_cast <float >(std::sqrt (mass2)) : 0 .f ;
539517 }
518+
540519 // Process: collision QA (DATA)
541520
542521 void processCollisions (aod::JetCollision const & collision,
@@ -595,7 +574,7 @@ struct JetHFAngularityTask {
595574 const float girth = computeLambda (jet, jetTracks, jetCandidates, 2 .f , 1 .f ); // λ_2^1
596575 const float mjet = computeJetMass (jetTracks, jetCandidates);
597576
598- TVector3 jetVector ( jet.px (), jet.py (), jet.pz ()) ;
577+ const std::array< double , 3 > jetMom{ jet.px (), jet.py (), jet.pz ()} ;
599578
600579 bool hasD0 = false ;
601580
@@ -604,10 +583,10 @@ struct JetHFAngularityTask {
604583
605584 hasD0 = true ;
606585
607- TVector3 d0Vector ( d0Candidate.px (), d0Candidate.py (), d0Candidate.pz ()) ;
586+ const std::array< double , 3 > d0Mom{ d0Candidate.px (), d0Candidate.py (), d0Candidate.pz ()} ;
608587
609- // Longitudinal momentum fraction
610- const float zParallel = jetVector. Dot (d0Vector ) / jetVector. Dot (jetVector );
588+ // Longitudinal momentum fraction: z_|| = (p_D0 . p_jet) / |p_jet|^2
589+ const float zParallel = RecoDecay::dotProd (d0Mom, jetMom ) / RecoDecay::mag2 (jetMom );
611590 // Angular separation between D0 and jet axis
612591 const float axisDistance = jetutilities::deltaR (jet, d0Candidate);
613592
@@ -700,7 +679,7 @@ struct JetHFAngularityTask {
700679 const float girth = computeLambda (jet, jetTracks, jetCandidates, 2 .f , 1 .f ); // λ_2^1
701680 const float mjet = computeJetMass (jetTracks, jetCandidates);
702681
703- TVector3 jetVector ( jet.px (), jet.py (), jet.pz ()) ;
682+ const std::array< double , 3 > jetMom{ jet.px (), jet.py (), jet.pz ()} ;
704683
705684 bool hasD0 = false ;
706685
@@ -709,9 +688,9 @@ struct JetHFAngularityTask {
709688
710689 hasD0 = true ;
711690
712- TVector3 d0Vector ( d0Candidate.px (), d0Candidate.py (), d0Candidate.pz ()) ;
691+ const std::array< double , 3 > d0Mom{ d0Candidate.px (), d0Candidate.py (), d0Candidate.pz ()} ;
713692
714- const float zParallel = jetVector. Dot (d0Vector ) / jetVector. Dot (jetVector );
693+ const float zParallel = RecoDecay::dotProd (d0Mom, jetMom ) / RecoDecay::mag2 (jetMom );
715694 const float axisDistance = jetutilities::deltaR (jet, d0Candidate);
716695
717696 const int8_t flagMcMatch = d0Candidate.flagMcMatchRec ();
@@ -807,7 +786,7 @@ struct JetHFAngularityTask {
807786 const float girth = computeLambda (jet, jetTracks, jetCandidates, 2 .f , 1 .f ); // λ_2^1
808787 const float mjet = computeJetMass (jetTracks, jetCandidates);
809788
810- TVector3 jetVector ( jet.px (), jet.py (), jet.pz ()) ;
789+ const std::array< double , 3 > jetMom{ jet.px (), jet.py (), jet.pz ()} ;
811790
812791 bool isGeoMatched = false ;
813792 float matchedPt = -1 .f ;
@@ -879,16 +858,16 @@ struct JetHFAngularityTask {
879858 registry.fill (HIST (" h_jet_matching_dr_mcd" ), matchedDR, mcWeight);
880859 }
881860
882- const int8_t geoStatus = static_cast <int8_t >(isGeoMatched ? 1 : 0 );
883- const int8_t candStatus = static_cast <int8_t >(isCandMatched ? 1 : 0 );
884- const int8_t cleanStatus = static_cast <int8_t >(isCleanMatched ? 1 : 0 );
861+ const int8_t geoStatus = static_cast <int8_t >(isGeoMatched);
862+ const int8_t candStatus = static_cast <int8_t >(isCandMatched);
863+ const int8_t cleanStatus = static_cast <int8_t >(isCleanMatched);
885864
886865 // ---- D0 candidate loop
887866 for (const auto & d0Candidate : jetCandidates) {
888867
889- TVector3 d0Vector ( d0Candidate.px (), d0Candidate.py (), d0Candidate.pz ()) ;
868+ const std::array< double , 3 > d0Mom{ d0Candidate.px (), d0Candidate.py (), d0Candidate.pz ()} ;
890869
891- const float zParallel = jetVector. Dot (d0Vector ) / jetVector. Dot (jetVector );
870+ const float zParallel = RecoDecay::dotProd (d0Mom, jetMom ) / RecoDecay::mag2 (jetMom );
892871 const float axisDistance = jetutilities::deltaR (jet, d0Candidate);
893872
894873 const int8_t flagMcMatch = d0Candidate.flagMcMatchRec ();
@@ -967,7 +946,7 @@ struct JetHFAngularityTask {
967946 const float girth = computeLambda (jet, jetParticles, jetCandidates, 2 .f , 1 .f ); // λ_2^1
968947 const float mjet = computeJetMass (jetParticles, jetCandidates, massD0MCP);
969948
970- TVector3 jetVector ( jet.px (), jet.py (), jet.pz ()) ;
949+ const std::array< double , 3 > jetMom{ jet.px (), jet.py (), jet.pz ()} ;
971950
972951 bool hasD0 = false ;
973952
@@ -976,9 +955,9 @@ struct JetHFAngularityTask {
976955
977956 hasD0 = true ;
978957
979- TVector3 d0Vector ( d0Particle.px (), d0Particle.py (), d0Particle.pz ()) ;
958+ const std::array< double , 3 > d0Mom{ d0Particle.px (), d0Particle.py (), d0Particle.pz ()} ;
980959
981- const float zParallel = jetVector. Dot (d0Vector ) / jetVector. Dot (jetVector );
960+ const float zParallel = RecoDecay::dotProd (d0Mom, jetMom ) / RecoDecay::mag2 (jetMom );
982961 const float axisDistance = jetutilities::deltaR (jet, d0Particle);
983962
984963 const int8_t flagMcMatch = d0Particle.flagMcMatchGen ();
@@ -1041,5 +1020,5 @@ struct JetHFAngularityTask {
10411020WorkflowSpec defineDataProcessing (ConfigContext const & cfgc)
10421021{
10431022 return WorkflowSpec{
1044- adaptAnalysisTask<JetHFAngularityTask>(cfgc, TaskName{ " jet-hf-ang-substructure " } )}; // o2-linter: disable=name/o2-task (templated struct)
1023+ adaptAnalysisTask<JetHFAngularityTask>(cfgc)};
10451024}
0 commit comments