[ALICE3] make the LUT collection writeable to file

ALICE3/Core/DelphesO2LutWriter.cxx

@@ -73,7 +73,7 @@ std::string DelphesO2LutWriter::LutBinning::toString() const
   return str;
 }
 
-bool DelphesO2LutWriter::fatSolve(lutEntry_t& lutEntry,
+bool DelphesO2LutWriter::fatSolve(o2::delphes::DelphesO2TrackSmearer::lutEntry_t& lutEntry,
                                   float pt,
                                   float eta,
                                   const float mass,
@@ -147,7 +147,7 @@ bool DelphesO2LutWriter::fwdSolve(float*, float, float, float)
 }
 #endif
 
-bool DelphesO2LutWriter::fwdPara(lutEntry_t& lutEntry, float pt, float eta, float mass, float Bfield)
+bool DelphesO2LutWriter::fwdPara(o2::delphes::DelphesO2TrackSmearer::lutEntry_t& lutEntry, float pt, float eta, float mass, float Bfield)
 {
   lutEntry.valid = false;
 
@@ -230,7 +230,7 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si
   }
 
   // write header
-  lutHeader_t lutHeader;
+  o2::delphes::DelphesO2TrackSmearer::lutHeader_t lutHeader;
   // pid
   lutHeader.pdg = pdg;
   const TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pdg);
@@ -245,7 +245,7 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si
     return;
   }
   lutHeader.field = field;
-  auto setMap = [](map_t& map, LutBinning b) {
+  auto setMap = [](o2::delphes::DelphesO2TrackSmearer::map_t& map, LutBinning b) {
     map.log = b.log;
     map.nbins = b.nbins;
     map.min = b.min;
@@ -267,7 +267,7 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si
   const int nrad = lutHeader.radmap.nbins;
   const int neta = lutHeader.etamap.nbins;
   const int npt = lutHeader.ptmap.nbins;
-  lutEntry_t lutEntry;
+  o2::delphes::DelphesO2TrackSmearer::lutEntry_t lutEntry;
 
   // write entries
   int nCalls = 0;
@@ -319,7 +319,7 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si
               retval = fwdSolve(lutEntry.covm, lutEntry.pt, lutEntry.eta, lutHeader.mass);
             }
             if (useDipole) { // Using the parametrization at the border of the barrel only for efficiency and momentum resolution
-              lutEntry_t lutEntryBarrel;
+              o2::delphes::DelphesO2TrackSmearer::lutEntry_t lutEntryBarrel;
               retval = fatSolve(lutEntryBarrel, lutEntry.pt, etaMaxBarrel, lutHeader.mass, itof, otof, q);
               lutEntry.valid = lutEntryBarrel.valid;
               lutEntry.covm[14] = lutEntryBarrel.covm[14];
@@ -339,7 +339,7 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si
           LOGF(info, "Diagonalizing");
           diagonalise(lutEntry);
           LOGF(info, "Writing");
-          lutFile.write(reinterpret_cast<char*>(&lutEntry), sizeof(lutEntry_t));
+          lutFile.write(reinterpret_cast<char*>(&lutEntry), sizeof(o2::delphes::DelphesO2TrackSmearer::lutEntry_t));
         }
       }
     }
@@ -349,7 +349,7 @@ void DelphesO2LutWriter::lutWrite(const char* filename, int pdg, float field, si
   lutFile.close();
 }
 
-void DelphesO2LutWriter::diagonalise(lutEntry_t& lutEntry)
+void DelphesO2LutWriter::diagonalise(o2::delphes::DelphesO2TrackSmearer::lutEntry_t& lutEntry)
 {
   static constexpr int kEig = 5;
   TMatrixDSym m(kEig);
@@ -399,7 +399,7 @@ TGraph* DelphesO2LutWriter::lutRead(const char* filename, int pdg, int what, int
   smearer.loadTable(pdg, filename);
   auto lutHeader = smearer.getLUTHeader(pdg);
   lutHeader->print();
-  map_t lutMap;
+  o2::delphes::DelphesO2TrackSmearer::map_t lutMap;
   switch (vs) {
     case kNch:
       lutMap = lutHeader->nchmap;

ALICE3/Core/DelphesO2LutWriter.h

@@ -46,7 +46,7 @@ class DelphesO2LutWriter
   void setAtLeastHits(int n) { mAtLeastHits = n; }
   void setAtLeastCorr(int n) { mAtLeastCorr = n; }
   void setAtLeastFake(int n) { mAtLeastFake = n; }
-  bool fatSolve(lutEntry_t& lutEntry,
+  bool fatSolve(o2::delphes::DelphesO2TrackSmearer::lutEntry_t& lutEntry,
                 float pt = 0.1,
                 float eta = 0.0,
                 const float mass = o2::track::pid_constants::sMasses[o2::track::PID::Pion],
@@ -57,14 +57,14 @@ class DelphesO2LutWriter
 
   void print() const;
   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);
+  bool fwdPara(o2::delphes::DelphesO2TrackSmearer::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.);
 
   o2::fastsim::FastTracker fat;
 
  private:
-  void diagonalise(lutEntry_t& lutEntry);
+  void diagonalise(o2::delphes::DelphesO2TrackSmearer::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

ALICE3/Core/DelphesO2TrackSmearer.cxx

@@ -13,27 +13,13 @@
 /// \file DelphesO2TrackSmearer.cxx
 /// \author Roberto Preghenella
 /// \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/lutCovm.hh https://github.com/AliceO2Group/DelphesO2/blob/master/src/lutCovm.hh
 ///                 DelphesO2/src/TrackSmearer.cc https://github.com/AliceO2Group/DelphesO2/blob/master/src/TrackSmearer.cc
 ///                 DelphesO2/src/TrackSmearer.hh https://github.com/AliceO2Group/DelphesO2/blob/master/src/TrackSmearer.hh
 /// @email: preghenella@bo.infn.it
 ///
 
-//////////////////////////////
-// DelphesO2/src/lutCovm.cc //
-//////////////////////////////
-
-/// @author: Roberto Preghenella
-/// @email: preghenella@bo.infn.it
-
-// #include "TrackSmearer.hh"
-// #include "TrackUtils.hh"
-// #include "TRandom.h"
-// #include <iostream>
-// #include <fstream>
-
 #include "ALICE3/Core/DelphesO2TrackSmearer.h"
 
 #include "ALICE3/Core/GeometryContainer.h"
@@ -51,7 +37,59 @@ namespace delphes
 
 /*****************************************************************/
 
-bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload)
+int DelphesO2TrackSmearer::getIndexPDG(const int pdg)
+{
+  switch (abs(pdg)) {
+    case 11:
+      return 0; // Electron
+    case 13:
+      return 1; // Muon
+    case 211:
+      return 2; // Pion
+    case 321:
+      return 3; // Kaon
+    case 2212:
+      return 4; // Proton
+    case 1000010020:
+      return 5; // Deuteron
+    case 1000010030:
+      return 6; // Triton
+    case 1000020030:
+      return 7; // Helium3
+    case 1000020040:
+      return 8; // Alphas
+    default:
+      return 2; // Default: pion
+  }
+}
+
+const char* DelphesO2TrackSmearer::getParticleName(int pdg)
+{
+  switch (abs(pdg)) {
+    case 11:
+      return "electron";
+    case 13:
+      return "muon";
+    case 211:
+      return "pion";
+    case 321:
+      return "kaon";
+    case 2212:
+      return "proton";
+    case 1000010020:
+      return "deuteron";
+    case 1000010030:
+      return "triton";
+    case 1000020030:
+      return "helium3";
+    case 1000020040:
+      return "alpha";
+    default:
+      return "pion"; // Default: pion
+  }
+}
+
+bool DelphesO2TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload)
 {
   if (!filename || filename[0] == '\0') {
     LOG(info) << " --- No LUT file provided for PDG " << pdg << ". Skipping load.";
@@ -109,17 +147,17 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload)
   const int nrad = mLUTHeader[ipdg]->radmap.nbins;
   const int neta = mLUTHeader[ipdg]->etamap.nbins;
   const int npt = mLUTHeader[ipdg]->ptmap.nbins;
-  mLUTEntry[ipdg] = new lutEntry_t****[nnch];
+  mLUTEntry[ipdg] = new DelphesO2TrackSmearer::lutEntry_t****[nnch];
   for (int inch = 0; inch < nnch; ++inch) {
-    mLUTEntry[ipdg][inch] = new lutEntry_t***[nrad];
+    mLUTEntry[ipdg][inch] = new DelphesO2TrackSmearer::lutEntry_t***[nrad];
     for (int irad = 0; irad < nrad; ++irad) {
-      mLUTEntry[ipdg][inch][irad] = new lutEntry_t**[neta];
+      mLUTEntry[ipdg][inch][irad] = new DelphesO2TrackSmearer::lutEntry_t**[neta];
       for (int ieta = 0; ieta < neta; ++ieta) {
-        mLUTEntry[ipdg][inch][irad][ieta] = new lutEntry_t*[npt];
+        mLUTEntry[ipdg][inch][irad][ieta] = new DelphesO2TrackSmearer::lutEntry_t*[npt];
         for (int ipt = 0; ipt < npt; ++ipt) {
-          mLUTEntry[ipdg][inch][irad][ieta][ipt] = new lutEntry_t;
-          lutFile.read(reinterpret_cast<char*>(mLUTEntry[ipdg][inch][irad][ieta][ipt]), sizeof(lutEntry_t));
-          if (lutFile.gcount() != sizeof(lutEntry_t)) {
+          mLUTEntry[ipdg][inch][irad][ieta][ipt] = new DelphesO2TrackSmearer::lutEntry_t;
+          lutFile.read(reinterpret_cast<char*>(mLUTEntry[ipdg][inch][irad][ieta][ipt]), sizeof(DelphesO2TrackSmearer::lutEntry_t));
+          if (lutFile.gcount() != sizeof(DelphesO2TrackSmearer::lutEntry_t)) {
             LOG(info) << " --- troubles reading covariance matrix entry for PDG " << pdg << ": " << localFilename << std::endl;
             LOG(info) << " --- expected/detected " << sizeof(lutHeader_t) << "/" << lutFile.gcount() << std::endl;
             return false;
@@ -137,7 +175,7 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload)
 
 /*****************************************************************/
 
-lutEntry_t* TrackSmearer::getLUTEntry(const int pdg, const float nch, const float radius, const float eta, const float pt, float& interpolatedEff)
+DelphesO2TrackSmearer::lutEntry_t* DelphesO2TrackSmearer::getLUTEntry(const int pdg, const float nch, const float radius, const float eta, const float pt, float& interpolatedEff)
 {
   const int ipdg = getIndexPDG(pdg);
   if (!mLUTHeader[ipdg]) {
@@ -210,7 +248,7 @@ lutEntry_t* TrackSmearer::getLUTEntry(const int pdg, const float nch, const floa
 
 /*****************************************************************/
 
-bool TrackSmearer::smearTrack(O2Track& o2track, lutEntry_t* lutEntry, float interpolatedEff)
+bool DelphesO2TrackSmearer::smearTrack(o2::track::TrackParCov& o2track, DelphesO2TrackSmearer::lutEntry_t* lutEntry, float interpolatedEff)
 {
   bool isReconstructed = true;
   // generate efficiency
@@ -263,7 +301,7 @@ bool TrackSmearer::smearTrack(O2Track& o2track, lutEntry_t* lutEntry, float inte
 
 /*****************************************************************/
 
-bool TrackSmearer::smearTrack(O2Track& o2track, int pdg, float nch)
+bool DelphesO2TrackSmearer::smearTrack(o2::track::TrackParCov& o2track, int pdg, float nch)
 {
   auto pt = o2track.getPt();
   switch (pdg) {
@@ -274,67 +312,67 @@ bool TrackSmearer::smearTrack(O2Track& o2track, int pdg, float nch)
   }
   auto eta = o2track.getEta();
   float interpolatedEff = 0.0f;
-  lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, interpolatedEff);
+  DelphesO2TrackSmearer::lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, interpolatedEff);
   if (!lutEntry || !lutEntry->valid)
     return false;
   return smearTrack(o2track, lutEntry, interpolatedEff);
 }
 
 /*****************************************************************/
 // relative uncertainty on pt
-double TrackSmearer::getPtRes(const int pdg, const float nch, const float eta, const float pt)
+double DelphesO2TrackSmearer::getPtRes(const int pdg, const float nch, const float eta, const float pt)
 {
   float dummy = 0.0f;
-  lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
+  DelphesO2TrackSmearer::lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
   auto val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt;
   return val;
 }
 
 /*****************************************************************/
 // relative uncertainty on eta
-double TrackSmearer::getEtaRes(const int pdg, const float nch, const float eta, const float pt)
+double DelphesO2TrackSmearer::getEtaRes(const int pdg, const float nch, const float eta, const float pt)
 {
   float dummy = 0.0f;
-  lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
+  DelphesO2TrackSmearer::lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
   auto sigmatgl = std::sqrt(lutEntry->covm[9]);                                  // sigmatgl2
   auto etaRes = std::fabs(std::sin(2.0 * std::atan(std::exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty
   etaRes /= lutEntry->eta;                                                       // relative uncertainty
   return etaRes;
 }
 /*****************************************************************/
 // absolute uncertainty on pt
-double TrackSmearer::getAbsPtRes(const int pdg, const float nch, const float eta, const float pt)
+double DelphesO2TrackSmearer::getAbsPtRes(const int pdg, const float nch, const float eta, const float pt)
 {
   float dummy = 0.0f;
-  lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
+  DelphesO2TrackSmearer::lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
   auto val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt * lutEntry->pt;
   return val;
 }
 
 /*****************************************************************/
 // absolute uncertainty on eta
-double TrackSmearer::getAbsEtaRes(const int pdg, const float nch, const float eta, const float pt)
+double DelphesO2TrackSmearer::getAbsEtaRes(const int pdg, const float nch, const float eta, const float pt)
 {
   float dummy = 0.0f;
-  lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
+  DelphesO2TrackSmearer::lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy);
   auto sigmatgl = std::sqrt(lutEntry->covm[9]);                                  // sigmatgl2
   auto etaRes = std::fabs(std::sin(2.0 * std::atan(std::exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty
   return etaRes;
 }
 /*****************************************************************/
 // efficiency
-double TrackSmearer::getEfficiency(const int pdg, const float nch, const float eta, const float pt)
+double DelphesO2TrackSmearer::getEfficiency(const int pdg, const float nch, const float eta, const float pt)
 {
   float efficiency = 0.0f;
   getLUTEntry(pdg, nch, 0., eta, pt, efficiency);
   return efficiency;
 }
 /*****************************************************************/
 // Only in DelphesO2
-// bool TrackSmearer::smearTrack(Track& track, bool atDCA)
+// bool DelphesO2TrackSmearer::smearTrack(Track& track, bool atDCA)
 // {
 
-//   O2Track o2track;
+//   o2::track::TrackParCov o2track;
 //   TrackUtils::convertTrackToO2Track(track, o2track, atDCA);
 //   int pdg = track.PID;
 //   float nch = mdNdEta; // use locally stored dNch/deta for the time being
@@ -344,11 +382,11 @@ double TrackSmearer::getEfficiency(const int pdg, const float nch, const float e
 //   return true;
 
 // #if 0
-//   lutEntry_t* lutEntry = getLUTEntry(track.PID, 0., 0., track.Eta, track.PT);
+//   DelphesO2TrackSmearer::lutEntry_t* lutEntry = getLUTEntry(track.PID, 0., 0., track.Eta, track.PT);
 //   if (!lutEntry)
 //     return;
 
-//   O2Track o2track;
+//   o2::track::TrackParCov o2track;
 //   TrackUtils::convertTrackToO2Track(track, o2track, atDCA);
 //   smearTrack(o2track, lutEntry);
 //   TrackUtils::convertO2TrackToTrack(o2track, track, atDCA);