Skip to content

Commit b2c6158

Browse files
preghenellasawenzel
preghenella
authored and
sawenzel
committed
Added functions to compute Ncoll, Npart and Nremn in Pythia8
1 parent d821aa3 commit b2c6158

2 files changed

Lines changed: 168 additions & 0 deletions

File tree

‎Generators/include/Generators/GeneratorPythia8.h‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -59,6 +59,24 @@ class GeneratorPythia8 : public Generator
5959
bool readString(std::string val) { return mPythia.readString(val, true); };
6060
bool readFile(std::string val) { return mPythia.readFile(val, true); };
6161

62+
/** utilities **/
63+
void getNcoll(int& nColl)
64+
{
65+
getNcoll(mPythia.info, nColl);
66+
};
67+
void getNpart(int& nPart)
68+
{
69+
getNpart(mPythia.info, nPart);
70+
};
71+
void getNpart(int& nProtonProj, int& nNeutronProj, int& nProtonTarg, int& nNeutronTarg)
72+
{
73+
getNpart(mPythia.info, nProtonProj, nNeutronProj, nProtonTarg, nNeutronTarg);
74+
};
75+
void getNremn(int& nProtonProj, int& nNeutronProj, int& nProtonTarg, int& nNeutronTarg)
76+
{
77+
getNremn(mPythia.event, nProtonProj, nNeutronProj, nProtonTarg, nNeutronTarg);
78+
};
79+
6280
protected:
6381
/** copy constructor **/
6482
GeneratorPythia8(const GeneratorPythia8&);
@@ -73,6 +91,10 @@ class GeneratorPythia8 : public Generator
7391

7492
/** utilities **/
7593
void selectFromAncestor(int ancestor, Pythia8::Event& inputEvent, Pythia8::Event& outputEvent);
94+
void getNcoll(const Pythia8::Info& info, int& nColl);
95+
void getNpart(const Pythia8::Info& info, int& nPart);
96+
void getNpart(const Pythia8::Info& info, int& nProtonProj, int& nNeutronProj, int& nProtonTarg, int& nNeutronTarg);
97+
void getNremn(const Pythia8::Event& event, int& nProtonProj, int& nNeutronProj, int& nProtonTarg, int& nNeutronTarg);
7698

7799
/** Pythia8 **/
78100
Pythia8::Pythia mPythia; //!

‎Generators/src/GeneratorPythia8.cxx‎

Lines changed: 146 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -250,6 +250,152 @@ void GeneratorPythia8::selectFromAncestor(int ancestor, Pythia8::Event& inputEve
250250
}
251251
}
252252

253+
/*****************************************************************/
254+
255+
void GeneratorPythia8::getNcoll(const Pythia8::Info& info, int& nColl)
256+
{
257+
258+
/** compute number of collisions from sub-collision information **/
259+
260+
#if PYTHIA_VERSION_INTEGER < 8300
261+
auto hiinfo = info.hiinfo;
262+
#else
263+
auto hiinfo = info.hiInfo;
264+
#endif
265+
266+
nColl = 0;
267+
268+
if (!hiinfo) {
269+
return;
270+
}
271+
272+
// loop over sub-collisions
273+
auto scptr = hiinfo->subCollisionsPtr();
274+
for (auto sc : *scptr) {
275+
276+
// wounded nucleon flag in projectile/target
277+
auto pW = sc.proj->status() == Pythia8::Nucleon::ABS; // according to C.Bierlich this should be == Nucleon::ABS
278+
auto tW = sc.targ->status() == Pythia8::Nucleon::ABS;
279+
280+
// increase number of collisions if both are wounded
281+
if (pW && tW) {
282+
nColl++;
283+
}
284+
}
285+
}
286+
287+
/*****************************************************************/
288+
289+
void GeneratorPythia8::getNpart(const Pythia8::Info& info, int& nPart)
290+
{
291+
292+
/** compute number of participants as the sum of all participants nucleons **/
293+
294+
int nProtonProj, nNeutronProj, nProtonTarg, nNeutronTarg;
295+
getNpart(info, nProtonProj, nNeutronProj, nProtonTarg, nNeutronTarg);
296+
nPart = nProtonProj + nNeutronProj + nProtonTarg + nNeutronTarg;
297+
}
298+
299+
/*****************************************************************/
300+
301+
void GeneratorPythia8::getNpart(const Pythia8::Info& info, int& nProtonProj, int& nNeutronProj, int& nProtonTarg, int& nNeutronTarg)
302+
{
303+
304+
/** compute number of participants from sub-collision information **/
305+
306+
#if PYTHIA_VERSION_INTEGER < 8300
307+
auto hiinfo = info.hiinfo;
308+
#else
309+
auto hiinfo = info.hiInfo;
310+
#endif
311+
312+
nProtonProj = nNeutronProj = nProtonTarg = nNeutronTarg = 0;
313+
if (!hiinfo) {
314+
return;
315+
}
316+
317+
// keep track of wounded nucleons
318+
std::vector<Pythia8::Nucleon*> projW;
319+
std::vector<Pythia8::Nucleon*> targW;
320+
321+
// loop over sub-collisions
322+
auto scptr = hiinfo->subCollisionsPtr();
323+
for (auto sc : *scptr) {
324+
325+
// wounded nucleon flag in projectile/target
326+
auto pW = sc.proj->status() == Pythia8::Nucleon::ABS || sc.proj->status() == Pythia8::Nucleon::DIFF; // according to C.Bierlich this should be == Nucleon::ABS || Nucleon::DIFF
327+
auto tW = sc.targ->status() == Pythia8::Nucleon::ABS || sc.targ->status() == Pythia8::Nucleon::DIFF;
328+
329+
// increase number of wounded projectile nucleons if not yet in the wounded vector
330+
if (pW && std::find(projW.begin(), projW.end(), sc.proj) == projW.end()) {
331+
projW.push_back(sc.proj);
332+
if (sc.proj->id() == 2212) {
333+
nProtonProj++;
334+
} else if (sc.proj->id() == 2112) {
335+
nNeutronProj++;
336+
}
337+
}
338+
339+
// increase number of wounded target nucleons if not yet in the wounded vector
340+
if (tW && std::find(targW.begin(), targW.end(), sc.targ) == targW.end()) {
341+
targW.push_back(sc.targ);
342+
if (sc.targ->id() == 2212) {
343+
nProtonTarg++;
344+
} else if (sc.targ->id() == 2112) {
345+
nNeutronTarg++;
346+
}
347+
}
348+
}
349+
}
350+
351+
/*****************************************************************/
352+
353+
void GeneratorPythia8::getNremn(const Pythia8::Event& event, int& nProtonProj, int& nNeutronProj, int& nProtonTarg, int& nNeutronTarg)
354+
{
355+
356+
/** compute number of spectators from the nuclear remnant of the beams **/
357+
358+
// reset
359+
nProtonProj = nNeutronProj = nProtonTarg = nNeutronTarg = 0;
360+
auto nNucRem = 0;
361+
362+
// particle loop
363+
auto nparticles = event.size();
364+
for (int ipa = 0; ipa < nparticles; ++ipa) {
365+
const auto particle = event[ipa];
366+
auto pdg = particle.id();
367+
368+
// nuclear remnants have pdg code = ±10LZZZAAA9
369+
if (pdg < 1000000000)
370+
continue; // must be nucleus
371+
if (pdg % 10 != 9)
372+
continue; // first digit must be 9
373+
nNucRem++;
374+
375+
// extract A, Z and L from pdg code
376+
pdg /= 10;
377+
auto A = pdg % 1000;
378+
pdg /= 1000;
379+
auto Z = pdg % 1000;
380+
pdg /= 1000;
381+
auto L = pdg % 10;
382+
383+
if (particle.pz() > 0.) {
384+
nProtonProj = Z;
385+
nNeutronProj = A - Z;
386+
}
387+
if (particle.pz() < 0.) {
388+
nProtonTarg = Z;
389+
nNeutronTarg = A - Z;
390+
}
391+
392+
} // end of particle loop
393+
394+
if (nNucRem > 2) {
395+
LOG(WARNING) << " GeneratorPythia8: found more than two nuclear remnants (weird)";
396+
}
397+
}
398+
253399
/*****************************************************************/
254400
/*****************************************************************/
255401

0 commit comments

Comments
 (0)