Skip to content

Commit 1a41563

Browse files
jackal1-66jackal1-66
committed
Multiple improvements
- Set number of possible loopers from short int to unsigned int - If time constrained are applied now values are always replaced inside the loopers - Added time histogram plot for loopers saved during generator destruction - Added option to set limits through collisiocontext - Added falt distribution option using time ratios from collisioncontext
1 parent 46dbb49 commit 1a41563

1 file changed

Lines changed: 126 additions & 42 deletions

File tree

‎MC/config/common/external/generator/TPCLoopers.C‎

Lines changed: 126 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -218,11 +218,11 @@ class GenTPCLoopers : public Generator
218218
// Set number of loopers if poissonian params are available
219219
if (mPoissonSet)
220220
{
221-
mNLoopersPairs = static_cast<short int>(std::round(mMultiplier[0] * PoissonPairs()));
221+
mNLoopersPairs = static_cast<unsigned int>(std::round(mMultiplier[0] * PoissonPairs()));
222222
}
223223
if (mGaussSet)
224224
{
225-
mNLoopersCompton = static_cast<short int>(std::round(mMultiplier[1] * GaussianElectrons()));
225+
mNLoopersCompton = static_cast<unsigned int>(std::round(mMultiplier[1] * GaussianElectrons()));
226226
}
227227
// Generate pairs
228228
for (int i = 0; i < mNLoopersPairs; ++i)
@@ -252,11 +252,11 @@ class GenTPCLoopers : public Generator
252252
// Set number of loopers if poissonian params are available
253253
if (mPoissonSet)
254254
{
255-
mNLoopersPairs = static_cast<short int>(std::round(mMultiplier[0] * PoissonPairs()));
255+
mNLoopersPairs = static_cast<unsigned int>(std::round(mMultiplier[0] * PoissonPairs()));
256256
}
257257
if (mGaussSet)
258258
{
259-
mNLoopersCompton = static_cast<short int>(std::round(mMultiplier[1] * GaussianElectrons()));
259+
mNLoopersCompton = static_cast<unsigned int>(std::round(mMultiplier[1] * GaussianElectrons()));
260260
}
261261
// Find the time histogram edge
262262
double time_constraint = t_hist->GetXaxis()->GetXmax();
@@ -268,10 +268,10 @@ class GenTPCLoopers : public Generator
268268
// Apply the inverse transformation using the scaler
269269
std::vector<double> transformed_pair = mScaler_pair->inverse_transform(pair);
270270
// Check if the time constraint is satisfied
271-
while ((transformed_pair[9] / 1e9) > time_constraint)
272-
{
271+
//while ((transformed_pair[9] / 1e9) > time_constraint)
272+
//{
273273
transformed_pair[9] = t_hist->GetRandom() * 1e9; // Regenerate time if it exceeds the constraint, multiplied by 1e9 to account for ratio in importParticles
274-
}
274+
//}
275275
LOG(info) << "Transformed pair time: " << transformed_pair[9] / 1e9 << " ns";
276276
mGenPairs.push_back(transformed_pair);
277277
}
@@ -281,10 +281,10 @@ class GenTPCLoopers : public Generator
281281
std::vector<double> electron = mONNX_compton->generate_sample();
282282
// Apply the inverse transformation using the scaler
283283
std::vector<double> transformed_electron = mScaler_compton->inverse_transform(electron);
284-
while ((transformed_electron[6] / 1e9) > time_constraint)
285-
{
284+
//while ((transformed_electron[6] / 1e9) > time_constraint)
285+
//{
286286
transformed_electron[6] = t_hist->GetRandom() * 1e9; // Regenerate time if it exceeds the constraint, multiplied by 1e9 to account for ratio in importParticles
287-
}
287+
//}
288288
LOG(info) << "Transformed electron time: " << transformed_electron[6] / 1e9 << " ns";
289289
mGenElectrons.push_back(transformed_electron);
290290
}
@@ -301,11 +301,11 @@ class GenTPCLoopers : public Generator
301301
// Set number of loopers if poissonian params are available
302302
if (mPoissonSet)
303303
{
304-
mNLoopersPairs = static_cast<short int>(std::round(mMultiplier[0] * PoissonPairs()));
304+
mNLoopersPairs = static_cast<unsigned int>(std::round(mMultiplier[0] * PoissonPairs()));
305305
}
306306
if (mGaussSet)
307307
{
308-
mNLoopersCompton = static_cast<short int>(std::round(mMultiplier[1] * GaussianElectrons()));
308+
mNLoopersCompton = static_cast<unsigned int>(std::round(mMultiplier[1] * GaussianElectrons()));
309309
}
310310
// Find the time histogram edge
311311
double time_constraint = time_limit;
@@ -317,11 +317,11 @@ class GenTPCLoopers : public Generator
317317
// Apply the inverse transformation using the scaler
318318
std::vector<double> transformed_pair = mScaler_pair->inverse_transform(pair);
319319
// Check if the time constraint is satisfied
320-
while ((transformed_pair[9]) > time_constraint)
321-
{
320+
//while ((transformed_pair[9]) > time_constraint)
321+
//{
322322
// Regenerate time using random number generator with maximum time_constraint
323323
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-
}
324+
//}
325325
LOG(info) << "Transformed pair time: " << transformed_pair[9] << " ns";
326326
mGenPairs.push_back(transformed_pair);
327327
}
@@ -331,10 +331,10 @@ class GenTPCLoopers : public Generator
331331
std::vector<double> electron = mONNX_compton->generate_sample();
332332
// Apply the inverse transformation using the scaler
333333
std::vector<double> transformed_electron = mScaler_compton->inverse_transform(electron);
334-
while ((transformed_electron[6] / 1e9) > time_constraint)
335-
{
334+
//while ((transformed_electron[6] / 1e9) > time_constraint)
335+
//{
336336
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-
}
337+
//}
338338
LOG(info) << "Transformed electron time: " << transformed_electron[6] << " ns";
339339
mGenElectrons.push_back(transformed_electron);
340340
}
@@ -399,9 +399,9 @@ class GenTPCLoopers : public Generator
399399
return true;
400400
}
401401

402-
short int PoissonPairs()
402+
unsigned int PoissonPairs()
403403
{
404-
short int poissonValue;
404+
unsigned int poissonValue;
405405
do
406406
{
407407
// Generate a Poisson-distributed random number with mean mPoisson[0]
@@ -411,9 +411,9 @@ class GenTPCLoopers : public Generator
411411
return poissonValue;
412412
}
413413

414-
short int GaussianElectrons()
414+
unsigned int GaussianElectrons()
415415
{
416-
short int gaussValue;
416+
unsigned int gaussValue;
417417
do
418418
{
419419
// Generate a Normal-distributed random number with mean mGass[0] and stddev mGauss[1]
@@ -423,7 +423,7 @@ class GenTPCLoopers : public Generator
423423
return gaussValue;
424424
}
425425

426-
void SetNLoopers(short int &nsig_pair, short int &nsig_compton)
426+
void SetNLoopers(unsigned int &nsig_pair, unsigned int &nsig_compton)
427427
{
428428
if(mPoissonSet) {
429429
LOG(info) << "Poissonian parameters correctly loaded.";
@@ -461,8 +461,8 @@ class GenTPCLoopers : public Generator
461461
double mGauss[4] = {0.0, 0.0, 0.0, 0.0}; // Mean, Std, Min, Max
462462
std::vector<std::vector<double>> mGenPairs;
463463
std::vector<std::vector<double>> mGenElectrons;
464-
short int mNLoopersPairs = -1;
465-
short int mNLoopersCompton = -1;
464+
unsigned int mNLoopersPairs = -1;
465+
unsigned int mNLoopersCompton = -1;
466466
std::array<float, 2> mMultiplier = {1., 1.};
467467
bool mPoissonSet = false;
468468
bool mGaussSet = false;
@@ -495,7 +495,32 @@ class GenLoopersInjector : public Generator
495495
exit(1);
496496
} else {
497497
LOG(info) << "Interaction Time records has " << mInteractionTimeRecords.size() << " entries.";
498+
mCollisionContext->printCollisionSummary();
498499
}
500+
mTimeHist = new TH1D("time_hist", "Time Histogram", 10000, 0., 0.0007); // Histogram to store absolute time values of loopers with respect to collision time in ns
501+
mStartTime = mInteractionTimeRecords[0].bc2ns();
502+
}
503+
504+
~GenLoopersInjector()
505+
{
506+
TFile *outputFile = TFile::Open("loopers_output.root", "RECREATE");
507+
//if (mTimeHist) {
508+
// mTimeHist->Write();
509+
// delete mTimeHist;
510+
//}
511+
TCanvas *canvas = new TCanvas("c1", "Looper Time Histogram", 800, 600);
512+
mTimeHist->Draw();
513+
// Create lines for each interaction time record
514+
for (const auto &record : mInteractionTimeRecords)
515+
{
516+
TLine *line = new TLine((record.bc2ns() - mStartTime)/1e9, 0, (record.bc2ns() - mStartTime)/1e9, mTimeHist->GetMaximum());
517+
line->SetLineColor(kRed);
518+
line->SetLineStyle(2); // Dashed line
519+
line->SetLineWidth(2);
520+
line->Draw("same");
521+
}
522+
canvas->Write();
523+
outputFile->Close();
499524
}
500525

501526
void setAdaptiveLoopers(Bool_t &adaptive)
@@ -522,6 +547,29 @@ class GenLoopersInjector : public Generator
522547
}
523548
}
524549

550+
void setFlatGas(Bool_t &flat, const Int_t &number = -1)
551+
{
552+
mFlatGas = flat;
553+
if (mFlatGas) {
554+
if (number < 0) {
555+
LOG(warn) << "Warning: Number of loopers per event must be non-negative! Switching option off.";
556+
mFlatGas = false;
557+
mFlatGasNumber = -1;
558+
} else {
559+
mFlatGasNumber = number;
560+
// Calculate the number of loopers to inject adaptively
561+
unsigned int nLoopers = static_cast<unsigned int>(number);
562+
unsigned int nLoopersPairs = static_cast<unsigned int>(std::round(nLoopers * mLoopsFractionPairs));
563+
unsigned int nLoopersCompton = nLoopers - nLoopersPairs;
564+
mGenTPCLoopers->SetNLoopers(nLoopersPairs, nLoopersCompton);
565+
mStartingTimeLoopersRatio = number / (mInteractionTimeRecords[1].bc2ns() - mInteractionTimeRecords[0].bc2ns()); // Ratio of loopers to time in seconds
566+
}
567+
} else {
568+
mFlatGasNumber = -1;
569+
}
570+
LOG(info) << "Flat gas loopers: " << (mFlatGas ? "ON" : "OFF") << ", Number of loopers per event: " << mFlatGasNumber;
571+
}
572+
525573
void setLoopsFractions(float &fraction, float &fractionPairs)
526574
{
527575
if (fraction < 0 || fraction >= 1)
@@ -544,7 +592,9 @@ class GenLoopersInjector : public Generator
544592
// Trivial, real work in importParticles
545593
LOG(info) << "mCurrentEvent is " << mCurrentEvent;
546594
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;
595+
LOG(info) << "Current time offset wrt BC: " << mInteractionTimeRecords[mCurrentEvent].getTimeOffsetWrtBC() << " ns";
596+
mTimeLimit = (mCurrentEvent < mInteractionTimeRecords.size() - 1) ? mInteractionTimeRecords[mCurrentEvent + 1].bc2ns() - mInteractionTimeRecords[mCurrentEvent].bc2ns() : mInteractionTimeRecords[1].bc2ns() - mInteractionTimeRecords[0].bc2ns();
597+
mTimeLoopersRatio = mFlatGas ? (mFlatGasNumber / mTimeLimit) : 1.0; // Ratio of loopers to time in nanoseconds, if flat gas is used
548598
mCurrentEvent++;
549599
return true;
550600
}
@@ -632,25 +682,43 @@ class GenLoopersInjector : public Generator
632682
}
633683
// Check size of mParticles stack and set loopers accordingly if mAdaptiveLoopers is true
634684
if (mAdaptiveLoopers) {
635-
int nParticles = mParticles.size();
636-
if (nParticles > 0)
637-
{
638-
// Calculate the number of loopers to inject adaptively
639-
short int nLoopers = static_cast<short int>(std::round((nParticles * mLoopsFraction) / (1 - mLoopsFraction)));
640-
short int nLoopersPairs = static_cast<short int>(std::round(nLoopers * mLoopsFractionPairs));
641-
short int nLoopersCompton = nLoopers - nLoopersPairs;
642-
mGenTPCLoopers->SetNLoopers(nLoopersPairs, nLoopersCompton);
643-
LOG(info) << "Adaptive loopers: " << nLoopers << " (pairs: " << nLoopersPairs << ", compton: " << nLoopersCompton << ")";
644-
if (mDecreasingLoopers) {
645-
mLoopsFraction = mLoopsFraction - mSlopeStep;
685+
unsigned int nLoopers, nLoopersPairs, nLoopersCompton;
686+
if (!mFlatGas){
687+
int nParticles = mParticles.size();
688+
if (nParticles > 0)
689+
{
690+
// Calculate the number of loopers to inject adaptively
691+
nLoopers = static_cast<unsigned int>(std::round((nParticles * mLoopsFraction) / (1 - mLoopsFraction)));
692+
nLoopersPairs = static_cast<unsigned int>(std::round(nLoopers * mLoopsFractionPairs));
693+
nLoopersCompton = nLoopers - nLoopersPairs;
694+
mGenTPCLoopers->SetNLoopers(nLoopersPairs, nLoopersCompton);
695+
LOG(info) << "Adaptive loopers: " << nLoopers << " (pairs: " << nLoopersPairs << ", compton: " << nLoopersCompton << ")";
696+
if (mDecreasingLoopers)
697+
{
698+
mLoopsFraction = mLoopsFraction - mSlopeStep;
699+
}
700+
}
701+
else
702+
{
703+
LOG(info) << "No particles found in O2 Kinematics, no loopers will be generated";
704+
return false;
646705
}
647706
} else {
648-
LOG(info) << "No particles found in O2 Kinematics, no loopers will be generated";
649-
return false;
707+
// If flat gas loopers are used, set the number of loopers to mFlatGasNumber
708+
LOG(info) << "Flat gas loopers: " << mFlatGasNumber << " per event " << "with startingRatio " << mStartingTimeLoopersRatio << " with mTimeLoopersRatio: " << mTimeLoopersRatio;
709+
nLoopers = mFlatGasNumber*(mStartingTimeLoopersRatio / mTimeLoopersRatio);
710+
nLoopersPairs = static_cast<unsigned int>(std::round(nLoopers * mLoopsFractionPairs));
711+
nLoopersCompton = nLoopers - nLoopersPairs;
712+
mGenTPCLoopers->SetNLoopers(nLoopersPairs, nLoopersCompton);
713+
LOG(info) << "Flat gas loopers: " << nLoopers << " (pairs: " << nLoopersPairs << ", compton: " << nLoopersCompton << ")";
650714
}
651715
}
652716
// Generate loopers using GenTPCLoopers
653717
// this is valid also when number of loopers is fixed
718+
if (mTimeLimit == 0.0){
719+
mTimeLimit = time_hist->GetXaxis()->GetXmax();
720+
LOG(info) << "Time limit for last event set to: " << mTimeLimit << " ns";
721+
}
654722
if (mTimeConstraint && mTimeLimit == 0.0) {
655723
mGenTPCLoopers->generateEvent(time_hist.get());
656724
} else if (mTimeLimit > 0.) {
@@ -667,6 +735,14 @@ class GenLoopersInjector : public Generator
667735
auto loopers = genProcessor(mGenTPCLoopers.get());
668736
LOG(info) << "Size of loopers: " << loopers.size();
669737
mParticles.insert(mParticles.end(), loopers.begin(), loopers.end());
738+
// Fill the time histogram with loopers time values
739+
LOG(info) << "Filling time histogram with loopers time values starting from " << (mInteractionTimeRecords[mCurrentEvent - 1].bc2ns() - mStartTime) << " ns";
740+
for (const auto &p : loopers) {
741+
LOG(info) << "Looper time: " << p.T() << " s";
742+
// Fill the histogram with absolute time values of loopers with respect to collision time
743+
// Normalise for loopers.size() the y axis
744+
mTimeHist->Fill(((mInteractionTimeRecords[mCurrentEvent - 1].bc2ns() - mStartTime)/1e9) + p.T());
745+
}
670746
LOG(info) << "Size mParticles at loopers stage " << mParticles.size();
671747
} else {
672748
LOG(error) << "Failed to import particles from O2 Kinematics";
@@ -691,6 +767,12 @@ class GenLoopersInjector : public Generator
691767
std::vector<o2::InteractionTimeRecord> mInteractionTimeRecords; // Interaction time records from collision context
692768
int mCurrentEvent = 0; // Current event number, used for decreasing loopers
693769
double mTimeLimit = 0.0; // Time limit for the current event, used for decreasing loopers
770+
TH1D *mTimeHist = nullptr; // Histogram to store absolute time values of loopers with respect to collision time
771+
double mStartTime = 0.0; // Start time of the event, used for decreasing loopers
772+
Bool_t mFlatGas = false; // Flag to indicate if flat gas loopers are used
773+
Int_t mFlatGasNumber = -1; // Number of flat gas loopers per event
774+
double mStartingTimeLoopersRatio = 1.0; // Ratio of number of loopers to the time interval between collision events
775+
double mTimeLoopersRatio = 1.0; // Ratio of number of loopers to the time interval between collision events
694776
};
695777

696778
} // namespace eventgen
@@ -704,8 +786,8 @@ class GenLoopersInjector : public Generator
704786
FairGenerator *
705787
Generator_TPCLoopers(std::string model_pairs = "tpcloopmodel.onnx", std::string model_compton = "tpcloopmodelcompton.onnx",
706788
std::string poisson = "poisson.csv", std::string gauss = "gauss.csv", std::string scaler_pair = "scaler_pair.json",
707-
std::string scaler_compton = "scaler_compton.json", std::array<float, 2> mult = {1., 1.}, short int nloopers_pairs = 1,
708-
short int nloopers_compton = 1)
789+
std::string scaler_compton = "scaler_compton.json", std::array<float, 2> mult = {1., 1.}, unsigned int nloopers_pairs = 1,
790+
unsigned int nloopers_compton = 1)
709791
{
710792
// Expand all environment paths
711793
model_pairs = gSystem->ExpandPathName(model_pairs.c_str());
@@ -773,7 +855,8 @@ FairGenerator *
773855
// 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
774856
FairGenerator *
775857
GeneratorLoopersInjector(std::string kineFileName = "genevents_Kine.root", std::string model_pairs = "tpcloopmodel.onnx", std::string model_compton = "tpcloopmodelcompton.onnx",
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)
858+
std::string scaler_pair = "scaler_pair.json", std::string scaler_compton = "scaler_compton.json", bool time_constraint = true, bool decreasing_loopers = false,
859+
bool flat_gas = true, const int loops_num = 2000, float loopers_fraction = 0.05, float fraction_pairs = 0.08)
777860
{
778861
// Expand all environment paths
779862
model_pairs = gSystem->ExpandPathName(model_pairs.c_str());
@@ -871,5 +954,6 @@ GeneratorLoopersInjector(std::string kineFileName = "genevents_Kine.root", std::
871954
generator->setTimeConstraint(time_constraint);
872955
// set decreasing loopers distribution
873956
generator->setDecreasingLoopers(decreasing_loopers);
957+
generator->setFlatGas(flat_gas, loops_num);
874958
return generator;
875959
}

0 commit comments

Comments
 (0)