From b8b12c8ddfd67c1ac3904cf3d778dc4577d068e5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Tue, 15 Jul 2025 11:19:57 +0200 Subject: [PATCH 1/9] Update DelphesO2LutWriter.h --- ALICE3/Core/DelphesO2LutWriter.h | 55 +++++++++++++++++++++++--------- 1 file changed, 40 insertions(+), 15 deletions(-) diff --git a/ALICE3/Core/DelphesO2LutWriter.h b/ALICE3/Core/DelphesO2LutWriter.h index 25faf7f382c..d9f00823076 100644 --- a/ALICE3/Core/DelphesO2LutWriter.h +++ b/ALICE3/Core/DelphesO2LutWriter.h @@ -24,8 +24,11 @@ #include "ALICE3/Core/DelphesO2TrackSmearer.h" #include "ALICE3/Core/FastTracker.h" + #include "TGraph.h" +#include + namespace o2::fastsim { class DelphesO2LutWriter @@ -34,34 +37,56 @@ class DelphesO2LutWriter DelphesO2LutWriter() = default; virtual ~DelphesO2LutWriter() = default; - o2::fastsim::FastTracker fat; - void diagonalise(lutEntry_t& lutEntry); - float etaMaxBarrel = 1.75f; - bool usePara = true; // use fwd parameterisation - bool useDipole = false; // use dipole i.e. flat parametrization for efficiency and momentum resolution - bool useFlatDipole = false; // use dipole i.e. flat parametrization outside of the barrel - - int mAtLeastHits = 4; - int mAtLeastCorr = 4; - int mAtLeastFake = 0; + // Setters + void setBinningNch(bool log, int nbins, float min, float max) { mNchBinning = {log, nbins, min, max}; } + void setBinningRadius(bool log, int nbins, float min, float max) { mRadiusBinning = {log, nbins, min, max}; } + void setBinningEta(bool log, int nbins, float min, float max) { mEtaBinning = {log, nbins, min, max}; } + void setBinningPt(bool log, int nbins, float min, float max) { mPtBinning = {log, nbins, min, max}; } + void setEtaMaxBarrel(float eta) { etaMaxBarrel = eta; } void SetAtLeastHits(int n) { mAtLeastHits = n; } void SetAtLeastCorr(int n) { mAtLeastCorr = n; } void SetAtLeastFake(int n) { mAtLeastFake = n; } - - void printLutWriterConfiguration(); bool fatSolve(lutEntry_t& lutEntry, float pt = 0.1, float eta = 0.0, const float mass = 0.13957000, - int itof = 0, - int otof = 0, + size_t itof = 0, + size_t otof = 0, int q = 1, const float nch = 1); + + void Print() const; bool fwdSolve(float* covm, float pt = 0.1, float eta = 0.0, float mass = 0.13957000); bool fwdPara(lutEntry_t& lutEntry, float pt = 0.1, float eta = 0.0, float mass = 0.13957000, float Bfield = 0.5); - void lutWrite(const char* filename = "lutCovm.dat", int pdg = 211, float field = 0.2, int itof = 0, int otof = 0); + void lutWrite(const char* filename = "lutCovm.dat", int pdg = 211, float field = 0.2, size_t itof = 0, size_t otof = 0); TGraph* lutRead(const char* filename, int pdg, int what, int vs, float nch = 0., float radius = 0., float eta = 0., float pt = 0.); + o2::fastsim::FastTracker fat; + + private: + void diagonalise(lutEntry_t& lutEntry); + float etaMaxBarrel = 1.75f; + bool usePara = true; // use fwd parameterisation + bool useDipole = false; // use dipole i.e. flat parametrization for efficiency and momentum resolution + bool useFlatDipole = false; // use dipole i.e. flat parametrization outside of the barrel + + int mAtLeastHits = 4; + int mAtLeastCorr = 4; + int mAtLeastFake = 0; + + // Binning of the LUT to make + struct LutBinning { + bool log; + int nbins; + float min; + float max; + std::string toString() const; + }; + LutBinning mNchBinning = {true, 20, 0.5f, 3.5f}; + LutBinning mRadiusBinning = {false, 1, 0.0f, 100.0f}; + LutBinning mEtaBinning = {false, 80, -4.0f, 4.0f}; + LutBinning mPtBinning = {true, 200, -2.0f, 2.0f}; + ClassDef(DelphesO2LutWriter, 1); }; } // namespace o2::fastsim From 676a1ce6bd10295f97188a2d57fb6e8a31796c1b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Tue, 15 Jul 2025 11:20:44 +0200 Subject: [PATCH 2/9] U --- ALICE3/Core/DelphesO2LutWriter.cxx | 108 ++++++++++++++++++----------- 1 file changed, 68 insertions(+), 40 deletions(-) diff --git a/ALICE3/Core/DelphesO2LutWriter.cxx b/ALICE3/Core/DelphesO2LutWriter.cxx index f3cba264b11..ae0359bf857 100644 --- a/ALICE3/Core/DelphesO2LutWriter.cxx +++ b/ALICE3/Core/DelphesO2LutWriter.cxx @@ -19,19 +19,22 @@ /// @email: preghenella@bo.infn.it /// -#include +#include "ALICE3/Core/DelphesO2LutWriter.h" #include "ALICE3/Core/DelphesO2TrackSmearer.h" -#include "ALICE3/Core/DelphesO2LutWriter.h" -#include -#include "TMatrixD.h" -#include "TVectorD.h" +#include "ALICE3/Core/FastTracker.h" +#include "ALICE3/Core/TrackUtilities.h" + #include "TAxis.h" -#include "TMatrixDSymEigen.h" #include "TDatabasePDG.h" #include "TLorentzVector.h" -#include "ALICE3/Core/FastTracker.h" -#include "ALICE3/Core/TrackUtilities.h" +#include "TMatrixD.h" +#include "TMatrixDSymEigen.h" +#include "TVectorD.h" + +#include +#include +#include // #define USE_FWD_PARAM #ifdef USE_FWD_PARAM @@ -41,21 +44,42 @@ namespace o2::fastsim { -void DelphesO2LutWriter::printLutWriterConfiguration() +void DelphesO2LutWriter::Print() const { - std::cout << " --- Printing configuration of LUT writer --- " << std::endl; - std::cout << " -> etaMaxBarrel = " << etaMaxBarrel << std::endl; - std::cout << " -> usePara = " << usePara << std::endl; - std::cout << " -> useDipole = " << useDipole << std::endl; - std::cout << " -> useFlatDipole = " << useFlatDipole << std::endl; + LOG(info) << " --- Printing configuration of LUT writer --- "; + LOG(info) << " -> etaMaxBarrel = " << etaMaxBarrel; + LOG(info) << " -> usePara = " << usePara; + LOG(info) << " -> useDipole = " << useDipole; + LOG(info) << " -> useFlatDipole = " << useFlatDipole; + LOG(info) << " -> mAtLeastHits = " << mAtLeastHits; + LOG(info) << " -> mAtLeastCorr = " << mAtLeastCorr; + LOG(info) << " -> mAtLeastFake = " << mAtLeastFake; + LOG(info) << " -> Nch Binning: = " << mNchBinning.toString(); + LOG(info) << " -> Radius Binning: = " << mRadiusBinning.toString(); + LOG(info) << " -> Eta Binning: = " << mEtaBinning.toString(); + LOG(info) << " -> Pt Binning: = " << mPtBinning.toString(); + LOG(info) << " --- End of configuration --- "; +} + +std::string DelphesO2LutWriter::LutBinning::toString() const +{ + std::string str = ""; + str.append(log ? "log" : "lin"); + str.append(" nbins: "); + str.append(std::to_string(nbins)); + str.append(" min: "); + str.append(std::to_string(min)); + str.append(" max: "); + str.append(std::to_string(max)); + return str; } bool DelphesO2LutWriter::fatSolve(lutEntry_t& lutEntry, float pt, float eta, const float mass, - int itof, - int otof, + size_t itof, + size_t otof, int q, const float nch) { @@ -65,13 +89,19 @@ bool DelphesO2LutWriter::fatSolve(lutEntry_t& lutEntry, tlv.SetPtEtaPhiM(pt, eta, 0., mass); o2::track::TrackParCov trkIn; o2::upgrade::convertTLorentzVectorToO2Track(q, tlv, {0., 0., 0.}, trkIn); + // tlv.Print(); + // return fmt::format("X:{:+.4e} Alp:{:+.3e} Par: {:+.4e} {:+.4e} {:+.4e} {:+.4e} {:+.4e} |Q|:{:d} {:s}\n", + // getX(), getAlpha(), getY(), getZ(), getSnp(), getTgl(), getQ2Pt(), getAbsCharge(), getPID().getName()); + // trkIn.print(); o2::track::TrackParCov trkOut; const int status = fat.FastTrack(trkIn, trkOut, nch); - if (status <= 0) { + if (status <= mAtLeastHits) { Printf(" --- fatSolve: FastTrack failed --- \n"); - tlv.Print(); + // tlv.Print(); return false; } + Printf(" --- fatSolve: FastTrack succeeded %d --- \n", status); + // trkOut.print(); lutEntry.valid = true; lutEntry.itof = fat.GetGoodHitProb(itof); lutEntry.otof = fat.GetGoodHitProb(otof); @@ -81,13 +111,15 @@ bool DelphesO2LutWriter::fatSolve(lutEntry_t& lutEntry, // define the efficiency auto totfake = 0.; lutEntry.eff = 1.; - for (int i = 1; i < 20; ++i) { + for (size_t i = 1; i < fat.GetNLayers(); ++i) { + if (fat.IsLayerInert(i)) + continue; // skip inert layers auto igoodhit = fat.GetGoodHitProb(i); if (igoodhit <= 0. || i == itof || i == otof) continue; lutEntry.eff *= igoodhit; auto pairfake = 0.; - for (int j = i + 1; j < 20; ++j) { + for (size_t j = i + 1; j < fat.GetNLayers(); ++j) { auto jgoodhit = fat.GetGoodHitProb(j); if (jgoodhit <= 0. || j == itof || j == otof) continue; @@ -180,7 +212,7 @@ bool DelphesO2LutWriter::fwdPara(lutEntry_t& lutEntry, float pt, float eta, floa return true; } -void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, int itof, int otof) +void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, size_t itof, size_t otof) { if (useFlatDipole && useDipole) { @@ -206,26 +238,21 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, in return; } lutHeader.field = field; + auto setMap = [](map_t& map, LutBinning b) { + map.log = b.log; + map.nbins = b.nbins; + map.min = b.min; + map.max = b.max; + }; // nch - lutHeader.nchmap.log = true; - lutHeader.nchmap.nbins = 20; - lutHeader.nchmap.min = 0.5; - lutHeader.nchmap.max = 3.5; + setMap(lutHeader.nchmap, mNchBinning); // radius - lutHeader.radmap.log = false; - lutHeader.radmap.nbins = 1; - lutHeader.radmap.min = 0.; - lutHeader.radmap.max = 100.; + setMap(lutHeader.radmap, mRadiusBinning); // eta - lutHeader.etamap.log = false; - lutHeader.etamap.nbins = 80; - lutHeader.etamap.min = -4.; - lutHeader.etamap.max = 4.; + setMap(lutHeader.etamap, mEtaBinning); // pt - lutHeader.ptmap.log = true; - lutHeader.ptmap.nbins = 200; - lutHeader.ptmap.min = -2; - lutHeader.ptmap.max = 2.; + setMap(lutHeader.ptmap, mPtBinning); + lutFile.write(reinterpret_cast(&lutHeader), sizeof(lutHeader)); // entries @@ -247,11 +274,11 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, in for (int irad = 0; irad < nrad; ++irad) { Printf(" --- writing irad = %d/%d", irad, nrad); for (int ieta = 0; ieta < neta; ++ieta) { - nCalls++; Printf(" --- writing ieta = %d/%d", ieta, neta); auto eta = lutHeader.etamap.eval(ieta); lutEntry.eta = lutHeader.etamap.eval(ieta); for (int ipt = 0; ipt < npt; ++ipt) { + nCalls++; Printf(" --- writing ipt = %d/%d", ipt, npt); lutEntry.pt = lutHeader.ptmap.eval(ipt); lutEntry.valid = true; @@ -325,7 +352,7 @@ void DelphesO2LutWriter::diagonalise(lutEntry_t& lutEntry) } } - m.Print(); + // m.Print(); TMatrixDSymEigen eigen(m); // eigenvalues vector TVectorD eigenVal = eigen.GetEigenValues(); @@ -345,7 +372,7 @@ void DelphesO2LutWriter::diagonalise(lutEntry_t& lutEntry) TGraph* DelphesO2LutWriter::lutRead(const char* filename, int pdg, int what, int vs, float nch, float radius, float eta, float pt) { - Printf(" --- reading LUT file %s", filename); + LOGF(info, " --- reading LUT file %s", filename); // vs static const int kNch = 0; static const int kEta = 1; @@ -363,6 +390,7 @@ TGraph* DelphesO2LutWriter::lutRead(const char* filename, int pdg, int what, int o2::delphes::DelphesO2TrackSmearer smearer; smearer.loadTable(pdg, filename); auto lutHeader = smearer.getLUTHeader(pdg); + lutHeader->print(); map_t lutMap; switch (vs) { case kNch: From 5bf9a68bdbce666af1d6a28dfbf05c6f77a4d335 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Tue, 15 Jul 2025 11:25:39 +0200 Subject: [PATCH 3/9] U --- ALICE3/Core/DelphesO2LutWriter.cxx | 71 +++++++++++++++--------------- ALICE3/Core/DelphesO2LutWriter.h | 8 ++-- 2 files changed, 39 insertions(+), 40 deletions(-) diff --git a/ALICE3/Core/DelphesO2LutWriter.cxx b/ALICE3/Core/DelphesO2LutWriter.cxx index ae0359bf857..f95650e765c 100644 --- a/ALICE3/Core/DelphesO2LutWriter.cxx +++ b/ALICE3/Core/DelphesO2LutWriter.cxx @@ -33,7 +33,6 @@ #include "TVectorD.h" #include -#include #include // #define USE_FWD_PARAM @@ -96,11 +95,11 @@ bool DelphesO2LutWriter::fatSolve(lutEntry_t& lutEntry, o2::track::TrackParCov trkOut; const int status = fat.FastTrack(trkIn, trkOut, nch); if (status <= mAtLeastHits) { - Printf(" --- fatSolve: FastTrack failed --- \n"); + LOGF(info, " --- fatSolve: FastTrack failed --- \n"); // tlv.Print(); return false; } - Printf(" --- fatSolve: FastTrack succeeded %d --- \n", status); + LOGF(info, " --- fatSolve: FastTrack succeeded %d --- \n", status); // trkOut.print(); lutEntry.valid = true; lutEntry.itof = fat.GetGoodHitProb(itof); @@ -176,7 +175,7 @@ bool DelphesO2LutWriter::fwdPara(lutEntry_t& lutEntry, float pt, float eta, floa const double dcaz_ms = sigma_alpha * r0 * std::cosh(eta); const double dcaz2 = dca_pos * dca_pos + dcaz_ms * dcaz_ms; - const float Leta = 2.8 / sinh(eta) - 0.01 * r0; // m + const float Leta = 2.8 / std::sinh(eta) - 0.01 * r0; // m const double relmomres_pos = 10e-6 * pt / 0.3 / Bfield / Leta / Leta * std::sqrt(720. / 15.); const float relmomres_barrel = std::sqrt(covmbarrel[14]) * pt; @@ -186,7 +185,7 @@ bool DelphesO2LutWriter::fwdPara(lutEntry_t& lutEntry, float pt, float eta, floa // interpolate MS contrib (rel resolution 0.4 at eta = 4) const float relmomres_MS_eta4 = 0.4 / beta * 0.5 / Bfield; - const float relmomres_MS = relmomres_MS_eta4 * pow(relmomres_MS_eta4 / relmomres_MS_barrel, (std::fabs(eta) - 4.) / (4. - etaMaxBarrel)); + const float relmomres_MS = relmomres_MS_eta4 * std::pow(relmomres_MS_eta4 / relmomres_MS_barrel, (std::fabs(eta) - 4.) / (4. - etaMaxBarrel)); const float momres_tot = pt * std::sqrt(relmomres_pos * relmomres_pos + relmomres_MS * relmomres_MS); // total absolute mom reso // Fill cov matrix diag @@ -205,7 +204,7 @@ bool DelphesO2LutWriter::fwdPara(lutEntry_t& lutEntry, float pt, float eta, floa // Check that all numbers are numbers for (int i = 0; i < 15; ++i) { if (std::isnan(lutEntry.covm[i])) { - Printf(" --- lutEntry.covm[%d] is NaN", i); + LOGF(info, " --- lutEntry.covm[%d] is NaN", i); return false; } } @@ -216,14 +215,14 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si { if (useFlatDipole && useDipole) { - Printf("Both dipole and dipole flat flags are on, please use only one of them"); + LOGF(info, "Both dipole and dipole flat flags are on, please use only one of them"); return; } // output file std::ofstream lutFile(filename, std::ofstream::binary); if (!lutFile.is_open()) { - Printf("Did not manage to open output file!!"); + LOGF(info, "Did not manage to open output file!!"); return; } @@ -234,7 +233,7 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si lutHeader.mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); const int q = std::abs(TDatabasePDG::Instance()->GetParticle(pdg)->Charge()) / 3; if (q <= 0) { - Printf("Negative or null charge (%f) for pdg code %i. Fix the charge!", TDatabasePDG::Instance()->GetParticle(pdg)->Charge(), pdg); + LOGF(info, "Negative or null charge (%f) for pdg code %i. Fix the charge!", TDatabasePDG::Instance()->GetParticle(pdg)->Charge(), pdg); return; } lutHeader.field = field; @@ -267,23 +266,23 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si int successfullCalls = 0; int failedCalls = 0; for (int inch = 0; inch < nnch; ++inch) { - Printf(" --- writing nch = %d/%d", inch, nnch); + LOGF(info, " --- writing nch = %d/%d", inch, nnch); auto nch = lutHeader.nchmap.eval(inch); lutEntry.nch = nch; fat.SetdNdEtaCent(nch); for (int irad = 0; irad < nrad; ++irad) { - Printf(" --- writing irad = %d/%d", irad, nrad); + LOGF(info, " --- writing irad = %d/%d", irad, nrad); for (int ieta = 0; ieta < neta; ++ieta) { - Printf(" --- writing ieta = %d/%d", ieta, neta); + LOGF(info, " --- writing ieta = %d/%d", ieta, neta); auto eta = lutHeader.etamap.eval(ieta); lutEntry.eta = lutHeader.etamap.eval(ieta); for (int ipt = 0; ipt < npt; ++ipt) { nCalls++; - Printf(" --- writing ipt = %d/%d", ipt, npt); + LOGF(info, " --- writing ipt = %d/%d", ipt, npt); lutEntry.pt = lutHeader.ptmap.eval(ipt); lutEntry.valid = true; if (std::fabs(eta) <= etaMaxBarrel) { // full lever arm ends at etaMaxBarrel - Printf("Solving in the barrel"); + LOGF(info, "Solving in the barrel"); // printf(" --- fatSolve: pt = %f, eta = %f, mass = %f, field=%f \n", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field); successfullCalls++; if (!fatSolve(lutEntry, lutEntry.pt, lutEntry.eta, lutHeader.mass, itof, otof, q)) { @@ -298,7 +297,7 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si failedCalls++; } } else { - Printf("Solving outside the barrel"); + LOGF(info, "Solving outside the barrel"); // printf(" --- fwdSolve: pt = %f, eta = %f, mass = %f, field=%f \n", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field); lutEntry.eff = 1.; lutEntry.eff2 = 1.; @@ -329,16 +328,16 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si failedCalls++; } } - Printf("Diagonalizing"); + LOGF(info, "Diagonalizing"); diagonalise(lutEntry); - Printf("Writing"); + LOGF(info, "Writing"); lutFile.write(reinterpret_cast(&lutEntry), sizeof(lutEntry_t)); } } } } - Printf(" --- finished writing LUT file %s", filename); - Printf(" --- successfull calls: %d/%d, failed calls: %d/%d", successfullCalls, nCalls, failedCalls, nCalls); + LOGF(info, " --- finished writing LUT file %s", filename); + LOGF(info, " --- successfull calls: %d/%d, failed calls: %d/%d", successfullCalls, nCalls, failedCalls, nCalls); lutFile.close(); } @@ -409,52 +408,52 @@ TGraph* DelphesO2LutWriter::lutRead(const char* filename, int pdg, int what, int g->SetTitle(Form("LUT for %s, pdg %d, vs %d, what %d", filename, pdg, vs, what)); switch (vs) { case kNch: - Printf(" --- vs = kNch"); + LOGF(info, " --- vs = kNch"); g->GetXaxis()->SetTitle("Nch"); break; case kEta: - Printf(" --- vs = kEta"); + LOGF(info, " --- vs = kEta"); g->GetXaxis()->SetTitle("#eta"); break; case kPt: - Printf(" --- vs = kPt"); + LOGF(info, " --- vs = kPt"); g->GetXaxis()->SetTitle("p_{T} (GeV/c)"); break; default: - Printf(" --- error: unknown vs %d", vs); + LOGF(info, " --- error: unknown vs %d", vs); return nullptr; } switch (what) { case kEfficiency: - Printf(" --- what = kEfficiency"); + LOGF(info, " --- what = kEfficiency"); g->GetYaxis()->SetTitle("Efficiency (%)"); break; case kEfficiency2: - Printf(" --- what = kEfficiency2"); + LOGF(info, " --- what = kEfficiency2"); g->GetYaxis()->SetTitle("Efficiency2 (%)"); break; case kEfficiencyInnerTOF: - Printf(" --- what = kEfficiencyInnerTOF"); + LOGF(info, " --- what = kEfficiencyInnerTOF"); g->GetYaxis()->SetTitle("Inner TOF Efficiency (%)"); break; case kEfficiencyOuterTOF: - Printf(" --- what = kEfficiencyOuterTOF"); + LOGF(info, " --- what = kEfficiencyOuterTOF"); g->GetYaxis()->SetTitle("Outer TOF Efficiency (%)"); break; case kPtResolution: - Printf(" --- what = kPtResolution"); + LOGF(info, " --- what = kPtResolution"); g->GetYaxis()->SetTitle("p_{T} Resolution (%)"); break; case kRPhiResolution: - Printf(" --- what = kRPhiResolution"); + LOGF(info, " --- what = kRPhiResolution"); g->GetYaxis()->SetTitle("R#phi Resolution (#mum)"); break; case kZResolution: - Printf(" --- what = kZResolution"); + LOGF(info, " --- what = kZResolution"); g->GetYaxis()->SetTitle("Z Resolution (#mum)"); break; default: - Printf(" --- error: unknown what %d", what); + LOGF(info, " --- error: unknown what %d", what); return nullptr; } @@ -475,7 +474,7 @@ TGraph* DelphesO2LutWriter::lutRead(const char* filename, int pdg, int what, int auto lutEntry = smearer.getLUTEntry(pdg, nch, radius, eta, pt, eff); if (!lutEntry->valid || lutEntry->eff == 0.) { if (!canBeInvalid) { - Printf(" --- warning: it cannot be invalid"); + LOGF(info, " --- warning: it cannot be invalid"); } continue; } @@ -508,16 +507,16 @@ TGraph* DelphesO2LutWriter::lutRead(const char* filename, int pdg, int what, int val = lutEntry->otof * 100.; // efficiency (%) break; case kPtResolution: - val = sqrt(lutEntry->covm[14]) * lutEntry->pt * 100.; // pt resolution (%) + val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt * 100.; // pt resolution (%) break; case kRPhiResolution: - val = sqrt(lutEntry->covm[0]) * 1.e4; // rphi resolution (um) + val = std::sqrt(lutEntry->covm[0]) * 1.e4; // rphi resolution (um) break; case kZResolution: - val = sqrt(lutEntry->covm[1]) * 1.e4; // z resolution (um) + val = std::sqrt(lutEntry->covm[1]) * 1.e4; // z resolution (um) break; default: - Printf(" --- error: unknown what %d", what); + LOGF(info, " --- error: unknown what %d", what); break; } g->AddPoint(cen, val); diff --git a/ALICE3/Core/DelphesO2LutWriter.h b/ALICE3/Core/DelphesO2LutWriter.h index d9f00823076..d60ab5e6afa 100644 --- a/ALICE3/Core/DelphesO2LutWriter.h +++ b/ALICE3/Core/DelphesO2LutWriter.h @@ -43,9 +43,9 @@ class DelphesO2LutWriter void setBinningEta(bool log, int nbins, float min, float max) { mEtaBinning = {log, nbins, min, max}; } void setBinningPt(bool log, int nbins, float min, float max) { mPtBinning = {log, nbins, min, max}; } void setEtaMaxBarrel(float eta) { etaMaxBarrel = eta; } - void SetAtLeastHits(int n) { mAtLeastHits = n; } - void SetAtLeastCorr(int n) { mAtLeastCorr = n; } - void SetAtLeastFake(int n) { mAtLeastFake = n; } + void setAtLeastHits(int n) { mAtLeastHits = n; } + void setAtLeastCorr(int n) { mAtLeastCorr = n; } + void setAtLeastFake(int n) { mAtLeastFake = n; } bool fatSolve(lutEntry_t& lutEntry, float pt = 0.1, float eta = 0.0, @@ -55,7 +55,7 @@ class DelphesO2LutWriter int q = 1, const float nch = 1); - void Print() const; + void print() const; bool fwdSolve(float* covm, float pt = 0.1, float eta = 0.0, float mass = 0.13957000); bool fwdPara(lutEntry_t& lutEntry, float pt = 0.1, float eta = 0.0, float mass = 0.13957000, float Bfield = 0.5); void lutWrite(const char* filename = "lutCovm.dat", int pdg = 211, float field = 0.2, size_t itof = 0, size_t otof = 0); From d8ef5392d8c6bd81c01b5f6d048bd57db0d5f58a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Tue, 15 Jul 2025 11:29:24 +0200 Subject: [PATCH 4/9] U --- ALICE3/Core/DelphesO2LutWriter.h | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/ALICE3/Core/DelphesO2LutWriter.h b/ALICE3/Core/DelphesO2LutWriter.h index d60ab5e6afa..a9889b2be73 100644 --- a/ALICE3/Core/DelphesO2LutWriter.h +++ b/ALICE3/Core/DelphesO2LutWriter.h @@ -8,15 +8,14 @@ // In applying this license CERN does not waive the privileges and immunities // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. - /// -/// @file DelphesO2LutWriter.h -/// @brief Porting to O2Physics of DelphesO2 code. +/// \file DelphesO2LutWriter.h +/// \brief Porting to O2Physics of DelphesO2 code. /// Minimal changes have been made to the original code for adaptation purposes, formatting and commented parts have been considered. /// Relevant sources: /// DelphesO2/src/lutWrite.cc https://github.com/AliceO2Group/DelphesO2/blob/master/src/lutWrite.cc -/// @author: Roberto Preghenella -/// @email: preghenella@bo.infn.it +/// \author: Roberto Preghenella +/// \email: preghenella@bo.infn.it /// #ifndef ALICE3_CORE_DELPHESO2LUTWRITER_H_ @@ -25,6 +24,8 @@ #include "ALICE3/Core/DelphesO2TrackSmearer.h" #include "ALICE3/Core/FastTracker.h" +#include "ReconstructionDataFormats/PID.h" + #include "TGraph.h" #include @@ -49,15 +50,15 @@ class DelphesO2LutWriter bool fatSolve(lutEntry_t& lutEntry, float pt = 0.1, float eta = 0.0, - const float mass = 0.13957000, + const float mass = o2::track::pid_constants::sMasses[o2::track::PID::Pion], size_t itof = 0, size_t otof = 0, int q = 1, const float nch = 1); void print() const; - bool fwdSolve(float* covm, float pt = 0.1, float eta = 0.0, float mass = 0.13957000); - bool fwdPara(lutEntry_t& lutEntry, float pt = 0.1, float eta = 0.0, float mass = 0.13957000, float Bfield = 0.5); + bool fwdSolve(float* covm, float pt = 0.1, float eta = 0.0, float mass = o2::track::pid_constants::sMasses[o2::track::PID::Pion]); + bool fwdPara(lutEntry_t& lutEntry, float pt = 0.1, float eta = 0.0, float mass = o2::track::pid_constants::sMasses[o2::track::PID::Pion], float Bfield = 0.5); void lutWrite(const char* filename = "lutCovm.dat", int pdg = 211, float field = 0.2, size_t itof = 0, size_t otof = 0); TGraph* lutRead(const char* filename, int pdg, int what, int vs, float nch = 0., float radius = 0., float eta = 0., float pt = 0.); From 51cd285f619cf731d56c14d08ab65bd679d37750 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Tue, 15 Jul 2025 11:31:00 +0200 Subject: [PATCH 5/9] I --- ALICE3/Core/DelphesO2LutWriter.cxx | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/ALICE3/Core/DelphesO2LutWriter.cxx b/ALICE3/Core/DelphesO2LutWriter.cxx index f95650e765c..f88acefd2b9 100644 --- a/ALICE3/Core/DelphesO2LutWriter.cxx +++ b/ALICE3/Core/DelphesO2LutWriter.cxx @@ -163,30 +163,30 @@ bool DelphesO2LutWriter::fwdPara(lutEntry_t& lutEntry, float pt, float eta, floa // parametrisation at eta = 4 const double beta = 1. / std::sqrt(1 + mass * mass / pt / pt / std::cosh(eta) / std::cosh(eta)); - const float dca_pos = 2.5e-4 / std::sqrt(3); // 2.5 micron/sqrt(3) + const float dcaPos = 2.5e-4 / std::sqrt(3); // 2.5 micron/sqrt(3) const float r0 = 0.5; // layer 0 radius [cm] const float r1 = 1.3; const float r2 = 2.5; const float x0layer = 0.001; // material budget (rad length) per layer - const double sigma_alpha = 0.0136 / beta / pt * std::sqrt(x0layer * std::cosh(eta)) * (1 + 0.038 * std::log(x0layer * std::cosh(eta))); - const double dcaxy_ms = sigma_alpha * r0 * std::sqrt(1 + r1 * r1 / (r2 - r0) / (r2 - r0)); - const double dcaxy2 = dca_pos * dca_pos + dcaxy_ms * dcaxy_ms; + const double sigmaAlpha = 0.0136 / beta / pt * std::sqrt(x0layer * std::cosh(eta)) * (1 + 0.038 * std::log(x0layer * std::cosh(eta))); + const double dcaxyMs = sigmaAlpha * r0 * std::sqrt(1 + r1 * r1 / (r2 - r0) / (r2 - r0)); + const double dcaxy2 = dcaPos * dcaPos + dcaxyMs * dcaxyMs; - const double dcaz_ms = sigma_alpha * r0 * std::cosh(eta); - const double dcaz2 = dca_pos * dca_pos + dcaz_ms * dcaz_ms; + const double dcazMs = sigmaAlpha * r0 * std::cosh(eta); + const double dcaz2 = dcaPos * dcaPos + dcazMs * dcazMs; const float Leta = 2.8 / std::sinh(eta) - 0.01 * r0; // m - const double relmomres_pos = 10e-6 * pt / 0.3 / Bfield / Leta / Leta * std::sqrt(720. / 15.); + const double relmomresPos = 10e-6 * pt / 0.3 / Bfield / Leta / Leta * std::sqrt(720. / 15.); - const float relmomres_barrel = std::sqrt(covmbarrel[14]) * pt; - const float Router = 1; // m - const float relmomres_pos_barrel = 10e-6 * pt / 0.3 / Bfield / Router / Router / std::sqrt(720. / 15.); - const float relmomres_MS_barrel = std::sqrt(relmomres_barrel * relmomres_barrel - relmomres_pos_barrel * relmomres_pos_barrel); + const float relmomresBarrel = std::sqrt(covmbarrel[14]) * pt; + const float rOuter = 1; // m + const float relmomresPosBarrel = 10e-6 * pt / 0.3 / Bfield / rOuter / rOuter / std::sqrt(720. / 15.); + const float relmomresMSBarrel = std::sqrt(relmomresBarrel * relmomresBarrel - relmomresPosBarrel * relmomresPosBarrel); // interpolate MS contrib (rel resolution 0.4 at eta = 4) - const float relmomres_MS_eta4 = 0.4 / beta * 0.5 / Bfield; - const float relmomres_MS = relmomres_MS_eta4 * std::pow(relmomres_MS_eta4 / relmomres_MS_barrel, (std::fabs(eta) - 4.) / (4. - etaMaxBarrel)); - const float momres_tot = pt * std::sqrt(relmomres_pos * relmomres_pos + relmomres_MS * relmomres_MS); // total absolute mom reso + const float relmomresMSEta4 = 0.4 / beta * 0.5 / Bfield; + const float relmomresMS = relmomresMSEta4 * std::pow(relmomresMSEta4 / relmomresMSBarrel, (std::fabs(eta) - 4.) / (4. - etaMaxBarrel)); + const float momresTot = pt * std::sqrt(relmomresPos * relmomresPos + relmomresMS * relmomresMS); // total absolute mom reso // Fill cov matrix diag for (int i = 0; i < 15; ++i) @@ -200,7 +200,7 @@ bool DelphesO2LutWriter::fwdPara(lutEntry_t& lutEntry, float pt, float eta, floa lutEntry.covm[2] = dcaz2; lutEntry.covm[5] = covmbarrel[5]; // sigma^2 sin(phi) lutEntry.covm[9] = covmbarrel[9]; // sigma^2 tanl - lutEntry.covm[14] = momres_tot * momres_tot / pt / pt / pt / pt; // sigma^2 1/pt + lutEntry.covm[14] = momresTot * momresTot / pt / pt / pt / pt; // sigma^2 1/pt // Check that all numbers are numbers for (int i = 0; i < 15; ++i) { if (std::isnan(lutEntry.covm[i])) { From f20d99e3f1cf0e6cda7c206b0632f4b7c3671696 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Tue, 15 Jul 2025 11:34:44 +0200 Subject: [PATCH 6/9] U --- ALICE3/Core/DelphesO2LutWriter.cxx | 29 ++++++++++++++++------------- ALICE3/Core/DelphesO2LutWriter.h | 1 - 2 files changed, 16 insertions(+), 14 deletions(-) diff --git a/ALICE3/Core/DelphesO2LutWriter.cxx b/ALICE3/Core/DelphesO2LutWriter.cxx index f88acefd2b9..a0de645de36 100644 --- a/ALICE3/Core/DelphesO2LutWriter.cxx +++ b/ALICE3/Core/DelphesO2LutWriter.cxx @@ -95,16 +95,17 @@ bool DelphesO2LutWriter::fatSolve(lutEntry_t& lutEntry, o2::track::TrackParCov trkOut; const int status = fat.FastTrack(trkIn, trkOut, nch); if (status <= mAtLeastHits) { - LOGF(info, " --- fatSolve: FastTrack failed --- \n"); + LOGF(info, " --- fatSolve: FastTrack failed ---"); // tlv.Print(); return false; } - LOGF(info, " --- fatSolve: FastTrack succeeded %d --- \n", status); + LOGF(info, " --- fatSolve: FastTrack succeeded %d ---", status); // trkOut.print(); lutEntry.valid = true; lutEntry.itof = fat.GetGoodHitProb(itof); lutEntry.otof = fat.GetGoodHitProb(otof); - for (int i = 0; i < 15; ++i) + static constexpr int nCov = 15; + for (int i = 0; i < nCov; ++i) lutEntry.covm[i] = trkOut.getCov()[i]; // define the efficiency @@ -151,20 +152,22 @@ bool DelphesO2LutWriter::fwdPara(lutEntry_t& lutEntry, float pt, float eta, floa lutEntry.valid = false; // parametrised forward response; interpolates between FAT at eta = 1.75 and a fixed parametrisation at eta = 4; only diagonal elements - if (std::fabs(eta) < etaMaxBarrel || std::fabs(eta) > 4) + static constexpr float etaLimit = 4.0f; + if (std::fabs(eta) < etaMaxBarrel || std::fabs(eta) > etaLimit) return false; if (!fatSolve(lutEntry, pt, etaMaxBarrel, mass)) return false; - float covmbarrel[15] = {0}; - for (int i = 0; i < 15; ++i) { + static constexpr int nCov = 15; + float covmbarrel[nCov] = {0}; + for (int i = 0; i < nCov; ++i) { covmbarrel[i] = lutEntry.covm[i]; } // parametrisation at eta = 4 const double beta = 1. / std::sqrt(1 + mass * mass / pt / pt / std::cosh(eta) / std::cosh(eta)); const float dcaPos = 2.5e-4 / std::sqrt(3); // 2.5 micron/sqrt(3) - const float r0 = 0.5; // layer 0 radius [cm] + const float r0 = 0.5; // layer 0 radius [cm] const float r1 = 1.3; const float r2 = 2.5; const float x0layer = 0.001; // material budget (rad length) per layer @@ -198,8 +201,8 @@ bool DelphesO2LutWriter::fwdPara(lutEntry_t& lutEntry, float pt, float eta, floa lutEntry.covm[2] = covmbarrel[2]; if (dcaz2 > lutEntry.covm[2]) lutEntry.covm[2] = dcaz2; - lutEntry.covm[5] = covmbarrel[5]; // sigma^2 sin(phi) - lutEntry.covm[9] = covmbarrel[9]; // sigma^2 tanl + lutEntry.covm[5] = covmbarrel[5]; // sigma^2 sin(phi) + lutEntry.covm[9] = covmbarrel[9]; // sigma^2 tanl lutEntry.covm[14] = momresTot * momresTot / pt / pt / pt / pt; // sigma^2 1/pt // Check that all numbers are numbers for (int i = 0; i < 15; ++i) { @@ -283,10 +286,10 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si lutEntry.valid = true; if (std::fabs(eta) <= etaMaxBarrel) { // full lever arm ends at etaMaxBarrel LOGF(info, "Solving in the barrel"); - // printf(" --- fatSolve: pt = %f, eta = %f, mass = %f, field=%f \n", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field); + // LOGF(info, " --- fatSolve: pt = %f, eta = %f, mass = %f, field=%f", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field); successfullCalls++; if (!fatSolve(lutEntry, lutEntry.pt, lutEntry.eta, lutHeader.mass, itof, otof, q)) { - // printf(" --- fatSolve: error \n"); + // LOGF(info, " --- fatSolve: error"); lutEntry.valid = false; lutEntry.eff = 0.; lutEntry.eff2 = 0.; @@ -298,7 +301,7 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si } } else { LOGF(info, "Solving outside the barrel"); - // printf(" --- fwdSolve: pt = %f, eta = %f, mass = %f, field=%f \n", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field); + // LOGF(info, " --- fwdSolve: pt = %f, eta = %f, mass = %f, field=%f", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field); lutEntry.eff = 1.; lutEntry.eff2 = 1.; bool retval = true; @@ -319,7 +322,7 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si lutEntry.eff2 = lutEntryBarrel.eff2; } if (!retval) { - printf(" --- fwdSolve: error \n"); + LOGF(info, " --- fwdSolve: error"); lutEntry.valid = false; for (int i = 0; i < 15; ++i) { lutEntry.covm[i] = 0.; diff --git a/ALICE3/Core/DelphesO2LutWriter.h b/ALICE3/Core/DelphesO2LutWriter.h index a9889b2be73..1528ab80bac 100644 --- a/ALICE3/Core/DelphesO2LutWriter.h +++ b/ALICE3/Core/DelphesO2LutWriter.h @@ -8,7 +8,6 @@ // In applying this license CERN does not waive the privileges and immunities // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// /// \file DelphesO2LutWriter.h /// \brief Porting to O2Physics of DelphesO2 code. /// Minimal changes have been made to the original code for adaptation purposes, formatting and commented parts have been considered. From 9174dbd1562e00b556f09c1fa8d25ef222041e3c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Tue, 15 Jul 2025 11:39:37 +0200 Subject: [PATCH 7/9] U --- ALICE3/Core/DelphesO2LutWriter.cxx | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/ALICE3/Core/DelphesO2LutWriter.cxx b/ALICE3/Core/DelphesO2LutWriter.cxx index a0de645de36..faff1c34431 100644 --- a/ALICE3/Core/DelphesO2LutWriter.cxx +++ b/ALICE3/Core/DelphesO2LutWriter.cxx @@ -346,8 +346,9 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si void DelphesO2LutWriter::diagonalise(lutEntry_t& lutEntry) { - TMatrixDSym m(5); - for (int i = 0, k = 0; i < 5; ++i) { + static constexpr int kEig = 5; + TMatrixDSym m(kEig); + for (int i = 0, k = 0; i < kEig; ++i) { for (int j = 0; j < i + 1; ++j, ++k) { m(i, j) = lutEntry.covm[k]; m(j, i) = lutEntry.covm[k]; @@ -358,17 +359,17 @@ void DelphesO2LutWriter::diagonalise(lutEntry_t& lutEntry) TMatrixDSymEigen eigen(m); // eigenvalues vector TVectorD eigenVal = eigen.GetEigenValues(); - for (int i = 0; i < 5; ++i) + for (int i = 0; i < kEig; ++i) lutEntry.eigval[i] = eigenVal[i]; // eigenvectors matrix TMatrixD eigenVec = eigen.GetEigenVectors(); - for (int i = 0; i < 5; ++i) - for (int j = 0; j < 5; ++j) + for (int i = 0; i < kEig; ++i) + for (int j = 0; j < kEig; ++j) lutEntry.eigvec[i][j] = eigenVec[i][j]; // inverse eigenvectors matrix eigenVec.Invert(); - for (int i = 0; i < 5; ++i) - for (int j = 0; j < 5; ++j) + for (int i = 0; i < kEig; ++i) + for (int j = 0; j < kEig; ++j) lutEntry.eiginv[i][j] = eigenVec[i][j]; } From 9b434cdc5fae24d514f7cd0881bedfa19b5c319e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Tue, 15 Jul 2025 11:40:05 +0200 Subject: [PATCH 8/9] U --- ALICE3/Core/DelphesO2LutWriter.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ALICE3/Core/DelphesO2LutWriter.cxx b/ALICE3/Core/DelphesO2LutWriter.cxx index faff1c34431..41339883cb9 100644 --- a/ALICE3/Core/DelphesO2LutWriter.cxx +++ b/ALICE3/Core/DelphesO2LutWriter.cxx @@ -43,7 +43,7 @@ namespace o2::fastsim { -void DelphesO2LutWriter::Print() const +void DelphesO2LutWriter::print() const { LOG(info) << " --- Printing configuration of LUT writer --- "; LOG(info) << " -> etaMaxBarrel = " << etaMaxBarrel; From f8b094dba819f5332ca976dea390780644e8f5cb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Tue, 15 Jul 2025 11:43:24 +0200 Subject: [PATCH 9/9] U --- ALICE3/Core/DelphesO2LutWriter.cxx | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/ALICE3/Core/DelphesO2LutWriter.cxx b/ALICE3/Core/DelphesO2LutWriter.cxx index 41339883cb9..987383092ec 100644 --- a/ALICE3/Core/DelphesO2LutWriter.cxx +++ b/ALICE3/Core/DelphesO2LutWriter.cxx @@ -233,10 +233,15 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si lutHeader_t lutHeader; // pid lutHeader.pdg = pdg; - lutHeader.mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); - const int q = std::abs(TDatabasePDG::Instance()->GetParticle(pdg)->Charge()) / 3; + const TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pdg); + if (!particle) { + LOG(fatal) << "Cannot find particle with PDG code " << pdg; + return; + } + lutHeader.mass = particle->Mass(); + const int q = std::abs(particle->Charge()) / 3; if (q <= 0) { - LOGF(info, "Negative or null charge (%f) for pdg code %i. Fix the charge!", TDatabasePDG::Instance()->GetParticle(pdg)->Charge(), pdg); + LOGF(info, "Negative or null charge (%f) for pdg code %i. Fix the charge!", particle->Charge(), pdg); return; } lutHeader.field = field;