11// External generator requested in https://its.cern.ch/jira/browse/O2-6235
22// for multidimensional performance studies
3+ // Example usage:
4+ // o2-sim -j 8 -o test -n 100 --seed 612 -g hybrid --configKeyValues "GeneratorHybrid.configFile=${O2DPG_MC_CONFIG_ROOT}/MC/config/common/external/generator/perfConf.json"
35namespace o2
46{
57 namespace eventgen
@@ -8,8 +10,7 @@ namespace o2
810 class GenPerf : public Generator
911 {
1012 public :
11- GenPerf (float fraction = 0.03f , unsigned short int nsig = 100 , unsigned short int tag = 1 , bool setDecayZ0 = false)
12- : mDecayZ0 (setDecayZ0 )
13+ GenPerf (float fraction = 0.03f , unsigned short int nsig = 100 , unsigned short int tag = 1 )
1314 {
1415 if (fraction == -1 ) {
1516 LOG (info ) << nsig << " Signal particles will be generated in each event" ;
@@ -36,15 +37,13 @@ namespace o2
3637 mTag = tag ;
3738 LOG (info ) << "Generator with tag " << mTag << " is selected" ;
3839 }
39- if (mDecayZ0 ){
40- LOG (info ) << "Z0 decays will be handled with Pythia8" ;
41- mPythia = std ::make_unique < Pythia8 ::Pythia > ();
42- // Configure Pythia8 with minimal beam setup for standalone decays
43- mPythia -> readString ("WeakSingleBoson:ffbar2gmZ = on" );
44- mPythia -> init (); // Initialize
45- } else {
46- LOG (info ) << "Z0 decays will be created but not transported (incompatible with Geant4 physics list)" ;
47- }
40+ LOG (info ) << "Z0 decays are handled with Pythia8" ;
41+ mPythia = std ::make_unique < Pythia8 ::Pythia > ();
42+ // Turn off all event generation - we only want to decay our Z0
43+ mPythia -> readString ("ProcessLevel:all = off" );
44+ // Disable standard event checks since we're manually building the event
45+ mPythia -> readString ("Check:event = off" );
46+ mPythia -> init (); // Initialize
4847 Generator ::setTimeUnit (1.0 );
4948 Generator ::setPositionUnit (1.0 );
5049 }
@@ -74,14 +73,10 @@ namespace o2
7473 for (int k = 0 ; k < nSig ; k ++ ){
7574 auto part = genMap [mTag ]();
7675 if (part .GetPdgCode () == 23 ) {
77- if (!mDecayZ0 ) {
78- mParticles .push_back (part );
79- } else {
80- auto daughters = decayZ0 (part );
81- for (auto & dau : daughters )
82- {
83- mParticles .push_back (dau );
84- }
76+ auto daughters = decayZ0 (part );
77+ for (auto & dau : daughters )
78+ {
79+ mParticles .push_back (dau );
8580 }
8681 } else {
8782 mParticles .push_back (part );
@@ -95,7 +90,6 @@ namespace o2
9590 unsigned short int mNSig = 0 ; // Number of particles to generate
9691 unsigned int mNUE = 0 ; // Number of tracks in the Underlying event
9792 unsigned short int mTag = 1 ; // Tag to select the generation function
98- bool mDecayZ0 = false; // Whether to decay Z0 with Pythia8
9993 std ::unique_ptr < Pythia8 ::Pythia > mPythia ; // Pythia8 instance for particle decays not present in the physics list of Geant4 (like Z0)
10094 const std ::vector < std ::shared_ptr < o2 ::eventgen ::Generator >>* mGenList = nullptr ; // Cached generators list
10195 std ::map < unsigned short int , std ::function < TParticle ()>> genMap ;
@@ -317,10 +311,22 @@ namespace o2
317311 std ::vector < TParticle > decayZ0 (TParticle & z0 )
318312 {
319313 std ::vector < TParticle > subparts ;
320- // Start a Pythia8 event with Z0
321- mPythia -> next ();
322- // Find the Z0 and its final products
323314 auto & event = mPythia -> event ;
315+ // Reset event record for new decay
316+ event .reset ();
317+ // Add the Z0 particle to the event record
318+ // Arguments: id, status, mother1, mother2, daughter1, daughter2,
319+ // col, acol, px, py, pz, e, m, scale, pol
320+ // Status code: 91 = incoming particles (needed for proper decay handling)
321+ int iZ0 = event .append (23 , 91 , 0 , 0 , 0 , 0 , 0 , 0 ,
322+ z0 .Px (), z0 .Py (), z0 .Pz (), z0 .Energy (), z0 .GetMass ());
323+ // Set production vertex
324+ event [iZ0 ].vProd (z0 .Vx (), z0 .Vy (), z0 .Vz (), 0.0 );
325+ // Forcing decay by calling hadron level function
326+ if (!mPythia -> forceHadronLevel ())
327+ {
328+ cout << "Warning: Z0 decay failed!" << endl ;
329+ }
324330 for (int j = 0 ; j < event .size (); ++ j )
325331 {
326332 const Pythia8 ::Particle & p = event [j ];
@@ -343,39 +349,34 @@ namespace o2
343349 iZ0 = event [iZ0 ].daughter1 ();
344350 }
345351 // Recursively collect all final-state descendants
346- std ::function < void (int )> collectFinalState = [& ](int idx )
352+ std ::function < void (int )> collectAllDescendants = [& ](int idx )
347353 {
348354 const Pythia8 ::Particle & particle = event [idx ];
349-
350- if (particle .isFinal ())
355+ subparts .push_back (TParticle (particle .id (), particle .status (),
356+ -1 , -1 , -1 , -1 ,
357+ particle .px (), particle .py (),
358+ particle .pz (), particle .e (),
359+ p .xProd (), p .yProd (), p .zProd (), p .tProd ()));
360+ subparts .back ().SetStatusCode (o2 ::mcgenstatus ::MCGenStatusEncoding (particle .status (), 0 ).fullEncoding );
361+ subparts .back ().SetUniqueID (mGenID + 1 );
362+ subparts .back ().SetBit (ParticleStatus ::kToBeDone ,
363+ o2 ::mcgenstatus ::getHepMCStatusCode (subparts .back ().GetStatusCode ()) == 1 );
364+ // Not final-state, recurse through daughters
365+ if (!particle .isFinal ())
351366 {
352- subparts .push_back (TParticle (particle .id (), particle .status (),
353- -1 , -1 , -1 , -1 ,
354- particle .px (), particle .py (),
355- particle .pz (), particle .e (),
356- p .xProd () + z0 .Vx (), p .yProd () + z0 .Vy (), p .zProd () + z0 .Vz (), 0.0 ));
357- subparts .back ().SetStatusCode (o2 ::mcgenstatus ::MCGenStatusEncoding (particle .status (), 0 ).fullEncoding );
358- subparts .back ().SetUniqueID (mGenID + 1 );
359- subparts .back ().SetBit (ParticleStatus ::kToBeDone ,
360- o2 ::mcgenstatus ::getHepMCStatusCode (subparts .back ().GetStatusCode ()) == 1 );
361- }
362- else
363- {
364- // Not final-state, recurse through daughters
365367 int d1 = particle .daughter1 ();
366368 int d2 = particle .daughter2 ();
367369 if (d1 > 0 )
368370 {
369371 for (int k = d1 ; k <= d2 ; ++ k )
370372 {
371- collectFinalState (k );
373+ collectAllDescendants (k );
372374 }
373375 }
374376 }
375377 };
376-
377378 // Start collecting from the final Z0
378- collectFinalState (iZ0 );
379+ collectAllDescendants (iZ0 );
379380 break ; // Found and processed the Z0
380381 }
381382 }
@@ -390,8 +391,8 @@ namespace o2
390391// fraction == -1 enables the fixed number of signal particles per event (nsig)
391392// tag selects the generator type to be used
392393FairGenerator *
393- Generator_Performance (const float fraction = 0.03f , const unsigned short int nsig = 100 , unsigned short int tag = 1 , bool setDecayZ0 = false )
394+ Generator_Performance (const float fraction = 0.03f , const unsigned short int nsig = 100 , unsigned short int tag = 1 )
394395{
395- auto generator = new o2 ::eventgen ::GenPerf (fraction , nsig , tag , setDecayZ0 );
396+ auto generator = new o2 ::eventgen ::GenPerf (fraction , nsig , tag );
396397 return generator ;
397398}
0 commit comments