2323#include " Common/Core/trackUtilities.h"
2424#include " Common/DataModel/TrackSelectionTables.h"
2525
26- #include " Framework/ASoAHelpers.h"
27- #include " Framework/AnalysisDataModel.h"
28- #include " Framework/AnalysisTask.h"
29- #include " Framework/HistogramRegistry.h"
30- #include " Framework/RunningWorkflowInfo.h"
31- #include " Framework/runDataProcessing.h"
32- #include " ReconstructionDataFormats/Track.h"
3326#include < CCDB/BasicCCDBManager.h>
27+ #include < CommonConstants/MathConstants.h>
28+ #include < CommonConstants/PhysicsConstants.h>
29+ #include < Framework/ASoAHelpers.h>
30+ #include < Framework/AnalysisDataModel.h>
31+ #include < Framework/AnalysisTask.h>
32+ #include < Framework/HistogramRegistry.h>
33+ #include < Framework/RunningWorkflowInfo.h>
34+ #include < Framework/runDataProcessing.h>
35+ #include < ReconstructionDataFormats/Track.h>
3436
35- #include < TF1.h>
3637#include < TPDGCode.h>
3738#include < TProfile.h>
3839#include < TRandom3.h>
4344using namespace o2 ;
4445using namespace o2 ::framework;
4546using namespace o2 ::framework::expressions;
47+ using LorentzVectorM = ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double >>;
4648
4749#define O2_DEFINE_CONFIGURABLE (NAME, TYPE, DEFAULT, HELP ) Configurable<TYPE> NAME{#NAME, DEFAULT, HELP};
4850
@@ -96,10 +98,10 @@ struct FlowMcUpc {
9698 // }
9799
98100 template <typename TTrack>
99- bool trackSelected (TTrack track)
101+ bool trackSelected (TTrack const & track)
100102 {
101- auto momentum = std::array<double , 3 >{track.px (), track.py (), track.pz ()};
102- double pt = RecoDecay:: pt (momentum );
103+ // auto momentum = std::array<double, 3>{track.px(), track.py(), track.pz()};
104+ double pt = track. pt ();
103105 if (pt < cfgPtCutMin || pt > cfgPtCutMax) {
104106 return false ;
105107 }
@@ -124,15 +126,28 @@ struct FlowMcUpc {
124126 histos.fill (HIST (" hImpactParameter" ), imp);
125127
126128 for (auto const & mcParticle : mcParticles) {
127- auto momentum = std::array<double , 3 >{mcParticle.px (), mcParticle.py (), mcParticle.pz ()};
128- double pt = RecoDecay::pt (momentum);
129- double phi = RecoDecay::phi (momentum);
130- double eta = RecoDecay::eta (momentum);
131- // focus on bulk: e, mu, pi, k, p
129+ // auto momentum = std::array<double, 3>{mcParticle.px(), mcParticle.py(), mcParticle.pz()};
130+ LorentzVectorM pMC (mcParticle.px (), mcParticle.py (), mcParticle.pz (), 0 ); // double phi = RecoDecay::phi(momentum); // focus on bulk: e, mu, pi, k, p
132131 int pdgCode = std::abs (mcParticle.pdgCode ());
133- if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != kKPlus && pdgCode != PDG_t::kProton )
132+ if (pdgCode == PDG_t::kElectron )
133+ pMC.SetM (o2::constants::physics::MassElectron);
134+ else if (pdgCode == PDG_t::kMuonMinus )
135+ pMC.SetM (o2::constants::physics::MassMuon);
136+ else if (pdgCode == PDG_t::kPiPlus )
137+ pMC.SetM (o2::constants::physics::MassPionCharged);
138+ else if (pdgCode == PDG_t::kKPlus )
139+ pMC.SetM (o2::constants::physics::MassKaonCharged);
140+ else if (pdgCode == PDG_t::kProton )
141+ pMC.SetM (o2::constants::physics::MassProton);
142+ else
134143 continue ;
135144
145+ double pt = pMC.Pt ();
146+ double eta = pMC.Eta ();
147+
148+ // if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != PDG_t::kKPlus && pdgCode != PDG_t::kProton)
149+ // continue;
150+
136151 if (!mcParticle.isPhysicalPrimary ())
137152 continue ;
138153 // if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance
@@ -164,16 +179,34 @@ struct FlowMcUpc {
164179 float vtxz = collision.posZ ();
165180
166181 for (const auto & track : tracks) {
167- auto momentum = std::array<double , 3 >{track.px (), track.py (), track.pz ()};
168- double pt = RecoDecay::pt (momentum);
169- double phi = RecoDecay::phi (momentum);
170- double eta = RecoDecay::eta (momentum);
182+ LorentzVectorM recoMC (track.px (), track.py (), track.pz (), 0 );
183+ // double phi = RecoDecay::phi(momentum);
184+ // focus on bulk: e, mu, pi, k, p
185+ // auto momentum = std::array<double, 3>{track.px(), track.py(), track.pz()};
186+ // double pt = track.pt();
187+ // double phi = RecoDecay::phi(momentum);
188+ // double eta = track.eta();
171189 if (!trackSelected (track) || (!track.has_udMcParticle ()))
172190 continue ;
173191 auto mcParticle = track.udMcParticle ();
174192 int pdgCode = std::abs (mcParticle.pdgCode ());
175- if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != kKPlus && pdgCode != PDG_t::kProton )
193+ if (pdgCode == PDG_t::kElectron )
194+ recoMC.SetM (o2::constants::physics::MassElectron);
195+ else if (pdgCode == PDG_t::kMuonMinus )
196+ recoMC.SetM (o2::constants::physics::MassMuon);
197+ else if (pdgCode == PDG_t::kPiPlus )
198+ recoMC.SetM (o2::constants::physics::MassPionCharged);
199+ else if (pdgCode == PDG_t::kKPlus )
200+ recoMC.SetM (o2::constants::physics::MassKaonCharged);
201+ else if (pdgCode == PDG_t::kProton )
202+ recoMC.SetM (o2::constants::physics::MassProton);
203+ else
176204 continue ;
205+
206+ double pt = recoMC.Pt ();
207+ double eta = recoMC.Eta ();
208+ // if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != PDG_t::kKPlus && pdgCode != PDG_t::kProton)
209+ // continue;
177210 // if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance
178211 // continue;
179212 if (!mcParticle.isPhysicalPrimary ())
0 commit comments