Skip to content

Commit 7580453

Browse files
jackal1-66jackal1-66
committed
Implemented loopers dependence on IR and collision system
1 parent 142f769 commit 7580453

File tree

2 files changed

+145
-12
lines changed

2 files changed

+145
-12
lines changed

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

Lines changed: 137 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -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)
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 *
562614
FairGenerator *
563615
Generator_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+
}

‎MC/config/common/ini/GeneratorLoopersFlatFile.ini‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
# TPC loopers injector using fit to calculate reference loopers number per orbit. File with fit function is pulled from the CCDB.
2+
# Three additional parameters are available in the function: (isPbPb = true, intRate = 38000, adjust = 0.)
3+
# isPbPb and intRate must be set in case the collision system is not PbPb at 38 kHz, while adjust can be used to decrease/increase
4+
# the number of loopers per orbit obtained from the reference (e.g. -0.1 reduces the loopers by 10%)
5+
#---> GeneratorTPCloopers
6+
[GeneratorExternal]
7+
fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/common/external/generator/TPCLoopers.C
8+
funcName = Generator_TPCLoopersOrbitRef("ccdb://Users/m/mgiacalo/WGAN_ExtGenPair", "ccdb://Users/m/mgiacalo/WGAN_ExtGenCompton", "${O2DPG_MC_CONFIG_ROOT}/MC/config/common/TPCloopers/ScalerPairParams.json", "${O2DPG_MC_CONFIG_ROOT}/MC/config/common/TPCloopers/ScalerComptonParams.json","ccdb://Users/m/mgiacalo/ClustersTrackRatio")

0 commit comments

Comments
 (0)