@@ -73,7 +73,7 @@ struct OnTheFlyDecayer {
7373
7474 o2::upgrade::Decayer decayer;
7575 Service<o2::framework::O2DatabasePDG> pdgDB;
76- std::map<int64_t , std::vector<o2::upgrade::OTFParticle>> mDecayDaughters ;
76+ std::map<int , std::vector<o2::upgrade::OTFParticle>> mDecayDaughters ;
7777
7878 Configurable<int > seed{" seed" 0 , " Set seed for particle decayer"
7979 Configurable<float > magneticField{" magneticField" 20 ., " Magnetic field (kG)"
@@ -84,9 +84,62 @@ struct OnTheFlyDecayer {
8484 HistogramRegistry histos{" histos"
8585 ConfigurableAxis axisPt{" axisPt" 0 .0f , 0 .1f , 0 .2f , 0 .3f , 0 .4f , 0 .5f , 0 .6f , 0 .7f , 0 .8f , 0 .9f , 1 .0f , 1 .1f , 1 .2f , 1 .3f , 1 .4f , 1 .5f , 1 .6f , 1 .7f , 1 .8f , 1 .9f , 2 .0f , 2 .2f , 2 .4f , 2 .6f , 2 .8f , 3 .0f , 3 .2f , 3 .4f , 3 .6f , 3 .8f , 4 .0f , 4 .4f , 4 .8f , 5 .2f , 5 .6f , 6 .0f , 6 .5f , 7 .0f , 7 .5f , 8 .0f , 9 .0f , 10 .0f , 11 .0f , 12 .0f , 13 .0f , 14 .0f , 15 .0f , 17 .0f , 19 .0f , 21 .0f , 23 .0f , 25 .0f , 30 .0f , 35 .0f , 40 .0f , 50 .0f }, " pt axis for QA histograms"
8686
87+ struct McParticleAlice3 {
88+ McParticleAlice3 () = default ;
89+ ~McParticleAlice3 () = default ;
90+ McParticleAlice3 (const McParticleAlice3& src) = default ;
91+ McParticleAlice3 (int collisionId,
92+ int pdgCode,
93+ int statusCode,
94+ int flags,
95+ int mother0,
96+ int mother1,
97+ int daughter0,
98+ int daughter1,
99+ float weight,
100+ float px, float py, float pz, float e,
101+ float vx, float vy, float vz, float vt,
102+ float phi, float eta, float pt, float p, float y,
103+ bool isAlive, bool isPrimary) : collisionId(collisionId),
104+ pdgCode (pdgCode),
105+ statusCode(statusCode),
106+ flags(flags),
107+ mothersIds{mother0, mother1},
108+ daughtersIdSlice{daughter0, daughter1},
109+ weight (weight),
110+ px (px),
111+ py (py),
112+ pz (pz),
113+ e (e),
114+ vx (vx),
115+ vy (vy),
116+ vz (vz),
117+ vt (vt),
118+ phi (phi),
119+ eta (eta),
120+ pt (pt),
121+ p (p),
122+ y (y),
123+ isAlive (isAlive),
124+ isPrimary (isPrimary) {}
125+ int collisionId;
126+ int pdgCode;
127+ int statusCode;
128+ int flags;
129+ int mothersIds[2 ];
130+ int daughtersIdSlice[2 ];
131+ float weight;
132+ float px, py, pz, e;
133+ float vx, vy, vz, vt;
134+ float phi, eta, pt, p, y;
135+ bool isAlive;
136+ bool isPrimary;
137+ };
138+
87139 std::vector<int > mEnabledDecays ;
88140 void init (o2::framework::InitContext&)
89141 {
142+ LOG (info) << " Initializing on-the-fly-decayer."
90143 decayer.setSeed (seed);
91144 decayer.setBField (magneticField);
92145 for (int i = 0 ; i < NumDecays; ++i) {
@@ -115,11 +168,13 @@ struct OnTheFlyDecayer {
115168 return std::find (mEnabledDecays .begin (), mEnabledDecays .end (), pdgCode) != mEnabledDecays .end ();
116169 }
117170
171+ std::vector<McParticleAlice3> mcParticlesAlice3;
118172 void process (aod::McCollision const &, aod::McParticles const & mcParticles)
119173 {
120174 mDecayDaughters .clear ();
175+ mcParticlesAlice3.clear ();
121176 u_int64_t nStoredDaughters = 0 ;
122- for (int64_t index{0 }; index < mcParticles.size (); ++index) {
177+ for (int index{0 }; index < static_cast < int >( mcParticles.size () ); ++index) {
123178 const auto & particle = mcParticles.iteratorAt (index);
124179 std::vector<o2::upgrade::OTFParticle> decayDaughters;
125180 static constexpr int MaxNestedDecays = 10 ;
@@ -129,13 +184,14 @@ struct OnTheFlyDecayer {
129184 o2::upgrade::convertMCParticleToO2Track (particle, o2track, pdgDB);
130185 decayDaughters = decayer.decayParticle (pdgDB, o2track, particle.pdgCode ());
131186 for (size_t idau{0 }; idau < decayDaughters.size (); ++idau) {
132- o2::upgrade::OTFParticle& dau = decayDaughters[idau];
187+ o2::upgrade::OTFParticle dau = decayDaughters[idau];
133188 o2::track::TrackParCov dauTrack;
134189 o2::upgrade::convertOTFParticleToO2Track (dau, dauTrack, pdgDB);
135190 if (canDecay (dau.pdgCode ())) {
136191 dau.setIsAlive (false );
137192 std::vector<o2::upgrade::OTFParticle> cascadingDaughers = decayer.decayParticle (pdgDB, dauTrack, dau.pdgCode ());
138- for (const auto & daudau : cascadingDaughers) {
193+ for (size_t idaudau{0 }; idaudau < cascadingDaughers.size (); ++idaudau) {
194+ o2::upgrade::OTFParticle daudau = cascadingDaughers[idaudau];
139195 decayDaughters.push_back (daudau);
140196 if (MaxNestedDecays < ++nDecays) {
141197 LOG (error) << " Seemingly stuck trying to perpetually decay products from pdg: " pdgCode ();
@@ -194,14 +250,13 @@ struct OnTheFlyDecayer {
194250 daughtersIdSlice[1 ] = static_cast <int >(particle.daughtersIds ()[1 ]);
195251 }
196252
197- std::span<const int > motherSpan (particle.mothersIds ().data (), particle.mothersIds ().size ());
198253 mDecayDaughters .emplace (index, decayDaughters);
199254 nStoredDaughters += decayDaughters.size ();
200255
201- float phi = o2::constants::math::PI + std::atan2 (-1 .0f * particle.py (), -1 .0f * particle.px ());
256+ const float phi = o2::constants::math::PI + std::atan2 (-1 .0f * particle.py (), -1 .0f * particle.px ());
202257 float eta; // As https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1922
203- float pt = std::sqrt (particle.px () * particle.px () + particle.py () * particle.py ());
204- float p = std::sqrt (particle.px () * particle.px () + particle.py () * particle.py () + particle.pz () * particle.pz ());
258+ const float pt = std::sqrt (particle.px () * particle.px () + particle.py () * particle.py ());
259+ const float p = std::sqrt (particle.px () * particle.px () + particle.py () * particle.py () + particle.pz () * particle.pz ());
205260 float y; // As https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1943
206261
207262 if ((p - particle.pz ()) < Tolerance) {
@@ -218,29 +273,29 @@ struct OnTheFlyDecayer {
218273
219274 // TODO: Particle status code
220275 // TODO: Expression columns
221- tableMcParticlesWithDau (particle.mcCollisionId (), particle.pdgCode (), particle.statusCode (),
222- particle.flags (), motherSpan, daughtersIdSlice, particle.weight (),
223- particle.px (), particle.py (), particle.pz (), particle.e (),
224- particle.vx (), particle.vy (), particle.vz (), particle.vt (),
225- phi, eta, pt, p, y, !canDecay (particle.pdgCode ()), true );
276+ auto mothers = particle.mothersIds ();
277+ int mother0 = mothers.size () > 0 ? mothers[0 ] : -1 ;
278+ int mother1 = mothers.size () > 1 ? mothers[1 ] : mother0;
279+ mcParticlesAlice3.push_back (McParticleAlice3{particle.mcCollisionId (), particle.pdgCode (), particle.statusCode (),
280+ particle.flags (), mother0, mother1,
281+ daughtersIdSlice[0 ], daughtersIdSlice[1 ], particle.weight (),
282+ particle.px (), particle.py (), particle.pz (), particle.e (),
283+ particle.vx (), particle.vy (), particle.vz (), particle.vt (),
284+ phi, eta, pt, p, y, !canDecay (particle.pdgCode ()), true });
226285 }
227286
228287 int daughtersIdSlice[2 ] = {-1 , -1 };
229288 for (const auto & [index, decayDaughters] : mDecayDaughters ) {
230289 for (const auto & dau : decayDaughters) {
231290 if (index >= mcParticles.size ()) {
232- LOG (warn ) << " --- Index " " out of bounds for mcParticles table of size " size ();
291+ LOG (error ) << " --- Index " " out of bounds for mcParticles table of size " size () << std::endl ;
233292 continue ;
234293 }
235294
236- auto mother = mcParticles.iteratorAt (index);
237- std::vector<int > motherIds = {static_cast <int >(index)};
238- std::span<const int > motherSpan (motherIds.data (), motherIds.size ());
239-
240- float phi = o2::constants::math::PI + std::atan2 (-1 .0f * dau.py (), -1 .0f * dau.px ());
295+ const float phi = o2::constants::math::PI + std::atan2 (-1 .0f * dau.py (), -1 .0f * dau.px ());
241296 float eta; // Conditional as https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1922
242- float pt = std::sqrt (dau.px () * dau.px () + dau.py () * dau.py ());
243- float p = std::sqrt (dau.px () * dau.px () + dau.py () * dau.py () + dau.pz () * dau.pz ());
297+ const float pt = std::sqrt (dau.px () * dau.px () + dau.py () * dau.py ());
298+ const float p = std::sqrt (dau.px () * dau.px () + dau.py () * dau.py () + dau.pz () * dau.pz ());
244299 float y; // Conditional as https://github.com/AliceO2Group/AliceO2/blob/dev/Framework/Core/include/Framework/AnalysisDataModel.h#L1943
245300
246301 if ((p - dau.pz ()) < Tolerance) {
@@ -283,13 +338,25 @@ struct OnTheFlyDecayer {
283338 // TODO: Particle status code
284339 // TODO: Expression columns
285340 // TODO: vt
286- tableMcParticlesWithDau (mother.mcCollisionId (), dau.pdgCode (), 1 ,
287- mother.flags (), motherSpan, daughtersIdSlice, mother.weight (),
288- dau.px (), dau.py (), dau.pz (), dau.e (),
289- dau.vx (), dau.vy (), dau.vz (), mother.vt (),
290- phi, eta, pt, p, y, dau.isAlive (), false );
341+ auto mother = mcParticles.iteratorAt (index);
342+ mcParticlesAlice3.push_back (McParticleAlice3{mother.mcCollisionId (), dau.pdgCode (), 1 ,
343+ -1 , index, index, daughtersIdSlice[0 ], daughtersIdSlice[1 ], mother.weight (),
344+ dau.px (), dau.py (), dau.pz (), dau.e (),
345+ dau.vx (), dau.vy (), dau.vz (), mother.vt (),
346+ phi, eta, pt, p, y, dau.isAlive (), false });
291347 }
292348 }
349+
350+ for (const auto & particle : mcParticlesAlice3) {
351+ std::span<const int > motherSpan (particle.mothersIds , 2 );
352+
353+ tableMcParticlesWithDau (particle.collisionId , particle.pdgCode , particle.statusCode ,
354+ particle.flags , motherSpan, particle.daughtersIdSlice , particle.weight ,
355+ particle.px , particle.py , particle.pz , particle.e ,
356+ particle.vx , particle.vy , particle.vz , particle.vt ,
357+ particle.phi , particle.eta , particle.pt , particle.p , particle.y ,
358+ particle.isAlive , particle.isPrimary );
359+ }
293360 }
294361};
295362
0 commit comments