2222#include " PWGJE/Core/JetUtilities.h"
2323#include " PWGJE/DataModel/Jet.h"
2424
25- #include " Framework/AnalysisTask.h"
26- #include " Framework/HistogramRegistry.h"
27- #include " Framework/runDataProcessing.h"
25+ #include < Framework/AnalysisTask.h>
26+ #include < Framework/HistogramRegistry.h>
27+ #include < Framework/runDataProcessing.h>
2828
29- #include " Math/GenVector/Boost.h"
30- #include " Math/Vector3D.h"
31- #include " Math/Vector4D.h"
32- #include " TRandom3.h"
33- #include " TVector3.h"
29+ #include < Math/GenVector/Boost.h>
30+ #include < Math/Vector3D.h>
31+ #include < Math/Vector4D.h>
32+ #include < TRandom3.h>
33+ #include < TVector3.h>
3434
3535#include < string>
3636#include < vector>
@@ -62,9 +62,10 @@ enum DecayChannel : uint8_t {
6262struct HfTaskDstarPolarisationInJet {
6363
6464 float massPi{0 .f };
65- float massProton{0 .f };
6665 float massKaon{0 .f };
66+ float massProton{0 .f };
6767 float massDstar{0 .f };
68+
6869 float bkgRotationAngleStep{0 .f };
6970
7071 uint8_t nMassHypos{0u };
@@ -85,11 +86,11 @@ struct HfTaskDstarPolarisationInJet {
8586 Configurable<bool > activateTHnSparseCosThStarJetAxis{" activateTHnSparseCosThStarJetAxis" , true , " Activate the THnSparse with cosThStar w.r.t. production axis" };
8687 Configurable<bool > activateTHnSparseCosThStarProduction{" activateTHnSparseCosThStarProduction" , true , " Activate the THnSparse with cosThStar w.r.t. production axis" };
8788 Configurable<bool > activatePartRecoDstar{" activatePartRecoDstar" , false , " Activate the study of partly reconstructed D*+ -> D0 (-> KPiPi0) Pi decays" };
88- float minInvMass {0 .f };
89- float maxInvMass {1000 .f };
89+ float invMassMin {0 .f };
90+ float invMassMax {1000 .f };
9091
9192 // / Application of rapidity cut for reconstructed candidates
92- Configurable<float > rapidityCut{ " rapidityCut " , 999 .f , " Max. value of reconstructed candidate rapidity (abs. value)" };
93+ Configurable<float > maxAbsRapidityCut{ " maxAbsRapidityCut " , 999 .f , " Max. value of reconstructed candidate rapidity (abs. value)" };
9394
9495 // Tables for MC jet matching
9596 using DstarJets = soa::Join<aod::DstarChargedJets, aod::DstarChargedJetConstituents>;
@@ -125,31 +126,36 @@ struct HfTaskDstarPolarisationInJet {
125126 void init (InitContext&)
126127 {
127128 // initialise event selection:
128- eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits (static_cast <std::string>( eventSelections) );
129+ eventSelectionBits = jetderiveddatautilities::initialiseEventSelectionBits (eventSelections. value );
129130
130131 // / check process functions
131- std::array<int , 4 > processes = {doprocessDstar, doprocessDstarWithMl, doprocessDstarMc, doprocessDstarMcWithMl};
132- const int nProcesses = std::accumulate (processes.begin (), processes.end (), 0 );
132+ const int nProcesses =
133+ static_cast <int >(doprocessDstar) +
134+ static_cast <int >(doprocessDstarWithMl) +
135+ static_cast <int >(doprocessDstarMc) +
136+ static_cast <int >(doprocessDstarMcWithMl);
137+ // std::array<int, 4> processes = {doprocessDstar, doprocessDstarWithMl, doprocessDstarMc, doprocessDstarMcWithMl};
138+ // const int nProcesses = std::accumulate(processes.begin(), processes.end(), 0);
133139 if (nProcesses > 1 ) {
134140 LOGP (fatal, " Only one process function should be enabled at a time, please check your configuration" );
135- } else if (nProcesses == 0 ) {
141+ }
142+ if (nProcesses == 0 ) {
136143 LOGP (fatal, " No process function enabled" );
137144 }
138145
139146 // / check output THnSparses
140147 std::array<int , 3 > sparses = {activateTHnSparseCosThStarHelicity, activateTHnSparseCosThStarJetAxis, activateTHnSparseCosThStarProduction};
141148 if (std::accumulate (sparses.begin (), sparses.end (), 0 ) == 0 ) {
142149 LOGP (fatal, " No output THnSparses enabled" );
143- } else {
144- if (activateTHnSparseCosThStarHelicity) {
145- LOGP (info, " THnSparse with cosThStar w.r.t. helicity axis active." );
146- }
147- if (activateTHnSparseCosThStarJetAxis) {
148- LOGP (info, " THnSparse with cosThStar w.r.t. jet axis active." );
149- }
150- if (activateTHnSparseCosThStarProduction) {
151- LOGP (info, " THnSparse with cosThStar w.r.t. production axis active." );
152- }
150+ }
151+ if (activateTHnSparseCosThStarHelicity) {
152+ LOGP (info, " THnSparse with cosThStar w.r.t. helicity axis active." );
153+ }
154+ if (activateTHnSparseCosThStarJetAxis) {
155+ LOGP (info, " THnSparse with cosThStar w.r.t. jet axis active." );
156+ }
157+ if (activateTHnSparseCosThStarProduction) {
158+ LOGP (info, " THnSparse with cosThStar w.r.t. production axis active." );
153159 }
154160
155161 if (activatePartRecoDstar && !(doprocessDstarMc || doprocessDstarMcWithMl)) {
@@ -188,8 +194,8 @@ struct HfTaskDstarPolarisationInJet {
188194 const AxisSpec thnAxisSelFlag{2 , -0 .5f , 1 .5f , " Sel flag" };
189195
190196 auto invMassBins = thnAxisInvMass.binEdges ;
191- minInvMass = invMassBins.front ();
192- maxInvMass = invMassBins.back ();
197+ invMassMin = invMassBins.front ();
198+ invMassMax = invMassBins.back ();
193199
194200 std::vector<AxisSpec> thnRecDataAxes = {thnAxisInvMass, thnAxisPt, thnAxisY, thnAxisInvMassD0, thnAxisCosThetaStar};
195201 if (activateTrackingSys) {
@@ -591,10 +597,9 @@ struct HfTaskDstarPolarisationInJet {
591597 // / \param ptBhadMother is the pt of the b-hadron mother (only in case of non-prompt)
592598 // / \param areDausInAcc is a flag indicating whether the daughters are in acceptance or not
593599 // / \param isPartRecoDstar is a flag indicating if it is a partly reconstructed Dstar->D0pi->Kpipipi0 meson (MC only)
594- template <charm_polarisation::CosThetaStarType cosThetaStarType>
595- void fillGenHistos (float ptCharmHad, float rapCharmHad, float cosThetaStar, int8_t origin, float ptBhadMother, bool areDausInAcc, int8_t charge, bool isPartRecoDstar, float zParallel, float jetPt)
600+ void fillGenHistos (charm_polarisation::CosThetaStarType cosThetaStarType, float ptCharmHad, float rapCharmHad, float cosThetaStar, int8_t origin, float ptBhadMother, bool areDausInAcc, int8_t charge, bool isPartRecoDstar, float zParallel, float jetPt)
596601 {
597- if constexpr (cosThetaStarType == charm_polarisation::CosThetaStarType::Helicity) { // Helicity
602+ if (cosThetaStarType == charm_polarisation::CosThetaStarType::Helicity) { // Helicity
598603 if (origin == RecoDecay::OriginType::Prompt) { // prompt
599604 if (!isPartRecoDstar) {
600605 registry.fill (HIST (" hGenPromptHelicity" ), ptCharmHad, rapCharmHad, cosThetaStar, areDausInAcc, charge, zParallel, jetPt);
@@ -608,7 +613,7 @@ struct HfTaskDstarPolarisationInJet {
608613 registry.fill (HIST (" hGenPartRecoNonPromptHelicity" ), ptCharmHad, rapCharmHad, cosThetaStar, ptBhadMother, areDausInAcc, charge, zParallel, jetPt);
609614 }
610615 }
611- } else if constexpr (cosThetaStarType == charm_polarisation::CosThetaStarType::Production) { // Production
616+ } else if (cosThetaStarType == charm_polarisation::CosThetaStarType::Production) { // Production
612617 if (origin == RecoDecay::OriginType::Prompt) { // prompt
613618 if (!isPartRecoDstar) {
614619 registry.fill (HIST (" hGenPromptProduction" ), ptCharmHad, rapCharmHad, cosThetaStar, areDausInAcc, charge, zParallel, jetPt);
@@ -622,7 +627,7 @@ struct HfTaskDstarPolarisationInJet {
622627 registry.fill (HIST (" hGenPartRecoNonPromptProduction" ), ptCharmHad, rapCharmHad, cosThetaStar, ptBhadMother, areDausInAcc, charge, zParallel, jetPt);
623628 }
624629 }
625- } else if constexpr (cosThetaStarType == charm_polarisation::CosThetaStarType::JetAxis) { // JetAxis
630+ } else if (cosThetaStarType == charm_polarisation::CosThetaStarType::JetAxis) { // JetAxis
626631 if (origin == RecoDecay::OriginType::Prompt) { // prompt
627632 if (!isPartRecoDstar) {
628633 registry.fill (HIST (" hGenPromptJetAxis" ), ptCharmHad, rapCharmHad, cosThetaStar, areDausInAcc, charge, zParallel, jetPt);
@@ -645,9 +650,9 @@ struct HfTaskDstarPolarisationInJet {
645650 bool isInSignalRegion (float invMass)
646651 {
647652 if constexpr (channel == charm_polarisation::DecayChannel::DstarToDzeroPi) { // D*+
648- float invMassMin = 0 .142f ;
649- float invMassMax = 0 .15f ;
650- if (invMassMin < invMass && invMass < invMassMax ) {
653+ float invMassSigMin = 0 .142f ;
654+ float invMassSigMax = 0 .15f ;
655+ if (invMassSigMin < invMass && invMass < invMassSigMax ) {
651656 return true ;
652657 }
653658 }
@@ -726,43 +731,45 @@ struct HfTaskDstarPolarisationInJet {
726731 pxDau = threeVecSoftPi[0 ];
727732 pyDau = threeVecSoftPi[1 ];
728733 pzDau = threeVecSoftPi[2 ];
729- pxCharmHad = threeVecSoftPi[0 ] + threeVecD0Prong0[0 ] + threeVecD0Prong1[0 ];
730- pyCharmHad = threeVecSoftPi[1 ] + threeVecD0Prong0[1 ] + threeVecD0Prong1[1 ];
731- pzCharmHad = threeVecSoftPi[2 ] + threeVecD0Prong0[2 ] + threeVecD0Prong1[2 ];
734+ std::array<float , 3 > threeVecCand = RecoDecay::pVec (threeVecSoftPi, threeVecD0Prong0, threeVecD0Prong1);
735+ pxCharmHad = threeVecCand[0 ];
736+ pyCharmHad = threeVecCand[1 ];
737+ pzCharmHad = threeVecCand[2 ];
732738 if (candidate.signProng1 () > 0 ) {
733739 invMassCharmHad = RecoDecay::m (std::array{threeVecD0Prong0, threeVecD0Prong1, threeVecSoftPi}, std::array{massPi, massKaon, massPi});
734740 invMassD0 = RecoDecay::m (std::array{threeVecD0Prong0, threeVecD0Prong1}, std::array{massPi, massKaon});
735741 } else {
736742 invMassCharmHad = RecoDecay::m (std::array{threeVecD0Prong0, threeVecD0Prong1, threeVecSoftPi}, std::array{massKaon, massPi, massPi});
737743 invMassD0 = RecoDecay::m (std::array{threeVecD0Prong0, threeVecD0Prong1}, std::array{massKaon, massPi});
738744 }
739- rapidity = RecoDecay::y (std::array{pxCharmHad, pyCharmHad, pzCharmHad} , massDstar);
745+ rapidity = RecoDecay::y (threeVecCand , massDstar);
740746 } else {
741747 isRotatedCandidate = 0 ;
742748 pxDau = candidate.pxProng1 ();
743749 pyDau = candidate.pyProng1 ();
744750 pzDau = candidate.pzProng1 ();
745- pxCharmHad = pxDau + candidate.pxProng0Charm () + candidate.pxProng1Charm ();
746- pyCharmHad = pyDau + candidate.pyProng0Charm () + candidate.pyProng1Charm ();
747- pzCharmHad = pzDau + candidate.pzProng0Charm () + candidate.pzProng1Charm ();
751+ std::array<float , 3 > threeVecCand = RecoDecay::pVec (std::array{candidate.pxProng1 (), candidate.pyProng1 (), candidate.pzProng1 ()},
752+ std::array{candidate.pyProng0Charm (), candidate.pxProng0Charm (), candidate.pzProng0Charm ()},
753+ std::array{candidate.pxProng1Charm (), candidate.pyProng1Charm (), candidate.pzProng1Charm ()});
754+ pxCharmHad = threeVecCand[0 ];
755+ pyCharmHad = threeVecCand[1 ];
756+ pzCharmHad = threeVecCand[2 ];
748757 invMassCharmHad = candidate.m ();
749758 invMassD0 = candidate.invMassCharm ();
750759 rapidity = candidate.y ();
751760 }
752761 invMassCharmHadForSparse = invMassCharmHad - invMassD0;
753762
754763 if constexpr (withMl) {
755- outputMl[0 ] = candidate.mlScores ()[0 ];
756- outputMl[1 ] = candidate.mlScores ()[1 ];
757- outputMl[2 ] = candidate.mlScores ()[2 ];
764+ std::copy_n (candidate.mlScores ().begin (), outputMl.size (), outputMl.begin ());
758765 }
759766 }
760- if (invMassCharmHadForSparse < minInvMass || invMassCharmHadForSparse > maxInvMass ) {
767+ if (invMassCharmHadForSparse < invMassMin || invMassCharmHadForSparse > invMassMax ) {
761768 continue ;
762769 }
763770
764771 // / apply rapidity selection on the reconstructed candidate
765- if (std::abs (rapidity) > rapidityCut ) {
772+ if (std::abs (rapidity) > maxAbsRapidityCut ) {
766773 continue ;
767774 }
768775
0 commit comments