@@ -291,6 +291,56 @@ class GenTPCLoopers : public Generator
291291 return true;
292292 }
293293
294+ Bool_t generateEvent (double & time_limit )
295+ {
296+ // Time_limit is in nanoseconds
297+ // Clear the vector of pairs
298+ mGenPairs .clear ();
299+ // Clear the vector of compton electrons
300+ mGenElectrons .clear ();
301+ // Set number of loopers if poissonian params are available
302+ if (mPoissonSet )
303+ {
304+ mNLoopersPairs = static_cast < short int > (std ::round (mMultiplier [0 ] * PoissonPairs ()));
305+ }
306+ if (mGaussSet )
307+ {
308+ mNLoopersCompton = static_cast < short int > (std ::round (mMultiplier [1 ] * GaussianElectrons ()));
309+ }
310+ // Find the time histogram edge
311+ double time_constraint = time_limit ;
312+ LOG (info ) << "Time constraint for loopers: " << time_constraint << " ns" ;
313+ // Generate pairs
314+ for (int i = 0 ; i < mNLoopersPairs ; ++ i )
315+ {
316+ std ::vector < double > pair = mONNX_pair -> generate_sample ();
317+ // Apply the inverse transformation using the scaler
318+ std ::vector < double > transformed_pair = mScaler_pair -> inverse_transform (pair );
319+ // Check if the time constraint is satisfied
320+ while ((transformed_pair [9 ]) > time_constraint )
321+ {
322+ // Regenerate time using random number generator with maximum time_constraint
323+ transformed_pair [9 ] = gRandom -> Uniform (0. , time_constraint ); // Regenerate time if it exceeds the constraint, scaling is not needed because time_limit is already in nanoseconds
324+ }
325+ LOG (info ) << "Transformed pair time: " << transformed_pair [9 ] << " ns" ;
326+ mGenPairs .push_back (transformed_pair );
327+ }
328+ // Generate compton electrons
329+ for (int i = 0 ; i < mNLoopersCompton ; ++ i )
330+ {
331+ std ::vector < double > electron = mONNX_compton -> generate_sample ();
332+ // Apply the inverse transformation using the scaler
333+ std ::vector < double > transformed_electron = mScaler_compton -> inverse_transform (electron );
334+ while ((transformed_electron [6 ] / 1e9 ) > time_constraint )
335+ {
336+ transformed_electron [6 ] = gRandom -> Uniform (0. , time_constraint ); // Regenerate time if it exceeds the constraint, scaling is not needed because time_limit is already in nanoseconds
337+ }
338+ LOG (info ) << "Transformed electron time: " << transformed_electron [6 ] << " ns" ;
339+ mGenElectrons .push_back (transformed_electron );
340+ }
341+ return true;
342+ }
343+
294344 Bool_t importParticles () override
295345 {
296346 // Get looper pairs from the event
@@ -436,6 +486,16 @@ class GenLoopersInjector : public Generator
436486 Generator ::setPositionUnit (1.0 );
437487 Generator ::setMomentumUnit (1.0 );
438488 Generator ::setEnergyUnit (1.0 );
489+ mContextFile = std ::filesystem ::exists ("collisioncontext.root" ) ? TFile ::Open ("collisioncontext.root" ) : nullptr ;
490+ mCollisionContext = mContextFile ? (o2 ::steer ::DigitizationContext * )mContextFile -> Get ("DigitizationContext" ) : nullptr ;
491+ mInteractionTimeRecords = mCollisionContext ? mCollisionContext -> getEventRecords () : std ::vector < o2 ::InteractionTimeRecord > {};
492+ if (mInteractionTimeRecords .empty ())
493+ {
494+ LOG (warn ) << "Error: No interaction time records found in the collision context!" ;
495+ exit (1 );
496+ } else {
497+ LOG (info ) << "Interaction Time records has " << mInteractionTimeRecords .size () << " entries." ;
498+ }
439499 }
440500
441501 void setAdaptiveLoopers (Bool_t & adaptive )
@@ -482,6 +542,10 @@ class GenLoopersInjector : public Generator
482542 Bool_t generateEvent () override
483543 {
484544 // Trivial, real work in importParticles
545+ LOG (info ) << "mCurrentEvent is " << mCurrentEvent ;
546+ LOG (info ) << "Current event time: " << ((mCurrentEvent < mInteractionTimeRecords .size () - 1 ) ? std ::to_string (mInteractionTimeRecords [mCurrentEvent + 1 ].bc2ns () - mInteractionTimeRecords [mCurrentEvent ].bc2ns ()) : "Final Event till the end" ) << " ns" ;
547+ mTimeLimit = (mCurrentEvent < mInteractionTimeRecords .size () - 1 ) ? mInteractionTimeRecords [mCurrentEvent + 1 ].bc2ns () - mInteractionTimeRecords [mCurrentEvent ].bc2ns () : 0.0 ;
548+ mCurrentEvent ++ ;
485549 return true;
486550 }
487551
@@ -587,8 +651,10 @@ class GenLoopersInjector : public Generator
587651 }
588652 // Generate loopers using GenTPCLoopers
589653 // this is valid also when number of loopers is fixed
590- if (mTimeConstraint ) {
654+ if (mTimeConstraint && mTimeLimit == 0.0 ) {
591655 mGenTPCLoopers -> generateEvent (time_hist .get ());
656+ } else if (mTimeLimit > 0. ) {
657+ mGenTPCLoopers -> generateEvent (mTimeLimit );
592658 } else {
593659 mGenTPCLoopers -> generateEvent ();
594660 }
@@ -620,6 +686,11 @@ class GenLoopersInjector : public Generator
620686 Bool_t mDecreasingLoopers = false; // Flag to indicate if decreasing loopers are used
621687 float mSlopeMinimum = 0.02 ; // Fraction of loopers to be injected adaptively
622688 float mSlopeStep = 0.01 ; // Step of decreasing slope
689+ TFile * mContextFile = nullptr ; // Input collision context file
690+ o2 ::steer ::DigitizationContext * mCollisionContext = nullptr ; // Pointer to the digitization context
691+ std ::vector < o2 ::InteractionTimeRecord > mInteractionTimeRecords ; // Interaction time records from collision context
692+ int mCurrentEvent = 0 ; // Current event number, used for decreasing loopers
693+ double mTimeLimit = 0.0 ; // Time limit for the current event, used for decreasing loopers
623694};
624695
625696} // namespace eventgen
@@ -702,7 +773,7 @@ FairGenerator *
702773// Loopers are considered adaptive by default, meaning that the number of loopers is determined by the number of particles in the kinematics file per event
703774FairGenerator *
704775GeneratorLoopersInjector (std ::string kineFileName = "genevents_Kine.root" , std ::string model_pairs = "tpcloopmodel.onnx" , std ::string model_compton = "tpcloopmodelcompton.onnx" ,
705- std ::string scaler_pair = "scaler_pair.json" , std ::string scaler_compton = "scaler_compton.json" , bool time_constraint = false , bool decreasing_loopers = true , float loopers_fraction = 0.05 , float fraction_pairs = 0.08 )
776+ std ::string scaler_pair = "scaler_pair.json" , std ::string scaler_compton = "scaler_compton.json" , bool time_constraint = true , bool decreasing_loopers = false , float loopers_fraction = 0.05 , float fraction_pairs = 0.08 )
706777{
707778 // Expand all environment paths
708779 model_pairs = gSystem -> ExpandPathName (model_pairs .c_str ()) ;
0 commit comments