@@ -224,7 +224,8 @@ class GenTPCLoopers : public Generator
224224 LOG (debug ) << "Current time offset wrt BC: " << mInteractionTimeRecords [mCurrentEvent ].getTimeOffsetWrtBC () << " ns" ;
225225 mTimeLimit = (mCurrentEvent < mInteractionTimeRecords .size () - 1 ) ? mInteractionTimeRecords [mCurrentEvent + 1 ].bc2ns () - mInteractionTimeRecords [mCurrentEvent ].bc2ns () : mTimeEnd - mInteractionTimeRecords [mCurrentEvent ].bc2ns ();
226226 // With flat gas the number of loopers are adapted based on time interval widths
227- nLoopers = mFlatGasNumber * (mTimeLimit / mIntTimeRecMean );
227+ // The denominator is either the LHC orbit (if mFlatGasOrbit is true) or the mean interaction time record interval
228+ nLoopers = mFlatGasOrbit ? (mFlatGasNumber * (mTimeLimit / o2 ::constants ::lhc ::LHCOrbitNS )) : (mFlatGasNumber * (mTimeLimit / mIntTimeRecMean ));
228229 nLoopersPairs = static_cast < unsigned int > (std ::round (nLoopers * mLoopsFractionPairs ));
229230 nLoopersCompton = nLoopers - nLoopersPairs ;
230231 SetNLoopers (nLoopersPairs , nLoopersCompton );
@@ -397,18 +398,28 @@ class GenTPCLoopers : public Generator
397398 }
398399 }
399400
400- void setFlatGas (Bool_t & flat , const Int_t & number = -1 )
401+ void setFlatGas (Bool_t & flat , const Int_t & number = -1 , const Int_t & nloopers_orbit = -1 )
401402 {
402403 mFlatGas = flat ;
403404 if (mFlatGas )
404405 {
405- if ( number < 0 )
406+ if ( nloopers_orbit > 0 )
406407 {
407- LOG ( warn ) << "Warning: Number of loopers per event must be non-negative! Switching option off." ;
408- mFlatGas = false ;
409- mFlatGasNumber = -1 ;
408+ mFlatGasOrbit = true ;
409+ mFlatGasNumber = nloopers_orbit ;
410+ LOG ( info ) << "Flat gas loopers will be generated using orbit reference." ;
410411 } else {
411- mFlatGasNumber = number ;
412+ mFlatGasOrbit = false;
413+ if (number < 0 )
414+ {
415+ LOG (warn ) << "Warning: Number of loopers per event must be non-negative! Switching option off." ;
416+ mFlatGas = false;
417+ mFlatGasNumber = -1 ;
418+ } else {
419+ mFlatGasNumber = number ;
420+ }
421+ }
422+ if (mFlatGas ){
412423 mContextFile = std ::filesystem ::exists ("collisioncontext.root" ) ? TFile ::Open ("collisioncontext.root" ) : nullptr ;
413424 mCollisionContext = mContextFile ? (o2 ::steer ::DigitizationContext * )mContextFile -> Get ("DigitizationContext" ) : nullptr ;
414425 mInteractionTimeRecords = mCollisionContext ? mCollisionContext -> getEventRecords () : std ::vector < o2 ::InteractionTimeRecord > {};
@@ -425,17 +436,17 @@ class GenTPCLoopers : public Generator
425436 mIntTimeRecMean += mInteractionTimeRecords [c + 1 ].bc2ns () - mInteractionTimeRecords [c ].bc2ns ();
426437 }
427438 mIntTimeRecMean /= (mInteractionTimeRecords .size () - 1 ); // Average interaction time record used as reference
428- const auto& hbfUtils = o2 ::raw ::HBFUtils ::Instance ();
439+ const auto & hbfUtils = o2 ::raw ::HBFUtils ::Instance ();
429440 // Get the start time of the second orbit after the last interaction record
430- const auto& lastIR = mInteractionTimeRecords .back ();
441+ const auto & lastIR = mInteractionTimeRecords .back ();
431442 o2 ::InteractionRecord finalOrbitIR (0 , lastIR .orbit + 2 ); // Final orbit, BC = 0
432443 mTimeEnd = finalOrbitIR .bc2ns ();
433444 LOG (debug ) << "Final orbit start time: " << mTimeEnd << " ns while last interaction record time is " << mInteractionTimeRecords .back ().bc2ns () << " ns" ;
434445 }
435446 } else {
436447 mFlatGasNumber = -1 ;
437448 }
438- LOG (info ) << "Flat gas loopers: " << (mFlatGas ? "ON" : "OFF" ) << ", Reference loopers number per event : " << mFlatGasNumber ;
449+ LOG (info ) << "Flat gas loopers: " << (mFlatGas ? "ON" : "OFF" ) << ", Reference loopers number per " << ( mFlatGasOrbit ? "orbit " : "event " ) << mFlatGasNumber ;
439450 }
440451
441452 void setFractionPairs (float & fractionPairs )
@@ -449,6 +460,45 @@ class GenTPCLoopers : public Generator
449460 LOG (info ) << "Pairs fraction set to: " << mLoopsFractionPairs ;
450461 }
451462
463+ void SetRate (const std ::string & rateFile , const bool isPbPb = true, const int & intRate = 50000 )
464+ {
465+ // Checking if the rate file exists and is not empty
466+ TFile rate_file (rateFile .c_str (), "READ" );
467+ if (!rate_file .IsOpen () || rate_file .IsZombie ())
468+ {
469+ LOG (fatal ) << "Error: Rate file is empty or does not exist!" ;
470+ exit (1 );
471+ }
472+ const char * fitName = isPbPb ? "fitPbPb" : "fitpp" ;
473+ auto fit = (TF1 * )rate_file .Get (fitName );
474+ if (!fit )
475+ {
476+ LOG (fatal ) << "Error: Could not find fit function '" << fitName << "' in rate file!" ;
477+ exit (1 );
478+ }
479+ auto ref = static_cast < int > (std ::floor (fit -> Eval (intRate / 1000. ))); // fit expects rate in kHz
480+ rate_file .Close ();
481+ if (ref <= 0 )
482+ {
483+ LOG (fatal ) << "Computed flat gas number reference per orbit is <=0" ;
484+ exit (1 );
485+ } else {
486+ LOG (info ) << "Set flat gas number to " << ref << " loopers per orbit using " << fitName << " from " << intRate << " Hz interaction rate." ;
487+ auto flat = true;
488+ setFlatGas (flat , -1 , ref );
489+ }
490+ }
491+
492+ void SetAdjust (const float & adjust = 0.f )
493+ {
494+ if (mFlatGas && mFlatGasOrbit && adjust >= -1.f && adjust != 0.f )
495+ {
496+ LOG (info ) << "Adjusting flat gas number per orbit by " << adjust * 100.f << "%" ;
497+ mFlatGasNumber = static_cast < int > (std ::round (mFlatGasNumber * (1.f + adjust )));
498+ LOG (info ) << "New flat gas number per orbit: " << mFlatGasNumber ;
499+ }
500+ }
501+
452502 private :
453503 std ::unique_ptr < ONNXGenerator > mONNX_pair = nullptr ;
454504 std ::unique_ptr < ONNXGenerator > mONNX_compton = nullptr ;
@@ -474,11 +524,13 @@ class GenTPCLoopers : public Generator
474524 o2 ::steer ::DigitizationContext * mCollisionContext = nullptr ; // Pointer to the digitization context
475525 std ::vector < o2 ::InteractionTimeRecord > mInteractionTimeRecords ; // Interaction time records from collision context
476526 Bool_t mFlatGas = false; // Flag to indicate if flat gas loopers are used
527+ Bool_t mFlatGasOrbit = false; // Flag to indicate if flat gas loopers are per orbit
477528 Int_t mFlatGasNumber = -1 ; // Number of flat gas loopers per event
478529 double mIntTimeRecMean = 1.0 ; // Average interaction time record used for the reference
479530 double mTimeLimit = 0.0 ; // Time limit for the current event
480531 double mTimeEnd = 0.0 ; // Time limit for the last event
481532 float mLoopsFractionPairs = 0.08 ; // Fraction of loopers from Pairs
533+ std ::string mRateFile = "" ; // File with clusters/rate information per orbit
482534};
483535
484536} // namespace eventgen
@@ -562,7 +614,7 @@ FairGenerator *
562614FairGenerator *
563615Generator_TPCLoopersFlat (std ::string model_pairs = "tpcloopmodel.onnx" , std ::string model_compton = "tpcloopmodelcompton.onnx" ,
564616 std ::string scaler_pair = "scaler_pair.json" , std ::string scaler_compton = "scaler_compton.json" ,
565- bool flat_gas = true, const int loops_num = 500 , float fraction_pairs = 0.08 )
617+ bool flat_gas = true, const int loops_num = 500 , float fraction_pairs = 0.08 , const int nloopers_orbit = -1 )
566618{
567619 // Expand all environment paths
568620 model_pairs = gSystem -> ExpandPathName (model_pairs .c_str ())
@@ -624,6 +676,79 @@ Generator_TPCLoopersFlat(std::string model_pairs = "tpcloopmodel.onnx", std::str
624676 model_compton = isAlien [1 ] || isCCDB [1 ] ? local_names [1 ] : model_compton ;
625677 auto generator = new o2 ::eventgen ::GenTPCLoopers (model_pairs , model_compton , "" , "" , scaler_pair , scaler_compton );
626678 generator -> setFractionPairs (fraction_pairs );
627- generator -> setFlatGas (flat_gas , loops_num );
679+ generator -> setFlatGas (flat_gas , loops_num , nloopers_orbit );
628680 return generator ;
629681}
682+
683+ // Generator with flat gas loopers. Reference number of loopers is provided per orbit via external file
684+ FairGenerator *
685+ Generator_TPCLoopersOrbitRef (std ::string model_pairs = "tpcloopmodel.onnx" , std ::string model_compton = "tpcloopmodelcompton.onnx" ,
686+ std ::string scaler_pair = "scaler_pair.json" , std ::string scaler_compton = "scaler_compton.json" ,
687+ std ::string nclxrate = "nclxrate.root" , bool isPbPb = true, const int intrate = 38000 , const float adjust = 0.f )
688+ {
689+ // Expand all environment paths
690+ model_pairs = gSystem -> ExpandPathName (model_pairs .c_str ())
691+ model_compton = gSystem -> ExpandPathName (model_compton .c_str ());
692+ scaler_pair = gSystem -> ExpandPathName (scaler_pair .c_str ());
693+ scaler_compton = gSystem -> ExpandPathName (scaler_compton .c_str ());
694+ nclxrate = gSystem -> ExpandPathName (nclxrate .c_str ());
695+ const std ::array < std ::string , 3 > models = {model_pairs , model_compton , nclxrate };
696+ const std ::array < std ::string , 3 > local_names = {"WGANpair.onnx" , "WGANcompton.onnx" , "nclxrate.root" };
697+ const std ::array < bool , 3 > isAlien = {models [0 ].starts_with ("alien://" ), models [1 ].starts_with ("alien://" ), nclxrate .starts_with ("alien://" )};
698+ const std ::array < bool , 3 > isCCDB = {models [0 ].starts_with ("ccdb://" ), models [1 ].starts_with ("ccdb://" ), nclxrate .starts_with ("ccdb://" )};
699+ if (std ::any_of (isAlien .begin (), isAlien .end (), [](bool v )
700+ { return v ; }))
701+ {
702+ if (!gGrid )
703+ {
704+ TGrid ::Connect ("alien://" );
705+ if (!gGrid )
706+ {
707+ LOG (fatal ) << "AliEn connection failed, check token." ;
708+ exit (1 );
709+ }
710+ }
711+ for (size_t i = 0 ; i < models .size (); ++ i )
712+ {
713+ if (isAlien [i ] && !TFile ::Cp (models [i ].c_str (), local_names [i ].c_str ()))
714+ {
715+ LOG (fatal ) << "Error: Model file " << models [i ] << " does not exist!" ;
716+ exit (1 );
717+ }
718+ }
719+ }
720+ if (std ::any_of (isCCDB .begin (), isCCDB .end (), [](bool v )
721+ { return v ; }))
722+ {
723+ o2 ::ccdb ::CcdbApi ccdb_api ;
724+ ccdb_api .init ("http://alice-ccdb.cern.ch" );
725+ for (size_t i = 0 ; i < models .size (); ++ i )
726+ {
727+ if (isCCDB [i ])
728+ {
729+ auto model_path = models [i ].substr (7 ); // Remove "ccdb://"
730+ // Treat filename if provided in the CCDB path
731+ auto extension = model_path .find (".onnx" );
732+ if (extension != std ::string ::npos )
733+ {
734+ auto last_slash = model_path .find_last_of ('/' );
735+ model_path = model_path .substr (0 , last_slash );
736+ }
737+ std ::map < std ::string , std ::string > filter ;
738+ if (!ccdb_api .retrieveBlob (model_path , "./" , filter , o2 ::ccdb ::getCurrentTimestamp (), false, local_names [i ].c_str ()))
739+ {
740+ LOG (fatal ) << "Error: issues in retrieving " << model_path << " from CCDB!" ;
741+ exit (1 );
742+ }
743+ }
744+ }
745+ }
746+ model_pairs = isAlien [0 ] || isCCDB [0 ] ? local_names [0 ] : model_pairs ;
747+ model_compton = isAlien [1 ] || isCCDB [1 ] ? local_names [1 ] : model_compton ;
748+ nclxrate = isAlien [2 ] || isCCDB [2 ] ? local_names [2 ] : nclxrate ;
749+ auto generator = new o2 ::eventgen ::GenTPCLoopers (model_pairs , model_compton , "" , "" , scaler_pair , scaler_compton );
750+ generator -> SetRate (nclxrate , isPbPb , intrate );
751+ // Adjust can be negative (-1 maximum) or positive to decrease or increase the number of loopers per orbit
752+ generator -> SetAdjust (adjust );
753+ return generator ;
754+ }
0 commit comments