Support for dN/deta dependence of the LUTs (#91)

Author: preghenella

src/TrackSmearer.cc

@@ -30,6 +30,10 @@ TrackSmearer::loadTable(int pdg, const char *filename)
     std::cout << " --- troubles reading covariance matrix header for PDG " << pdg << ": " << filename << std::endl;
     return false;
   }
+  if (mLUTHeader[ipdg]->version != LUTCOVM_VERSION) {
+    std::cout << " --- LUT header version mismatch: expected/detected = " << LUTCOVM_VERSION << "/" << mLUTHeader[ipdg]->version << std::endl;
+    return false;
+  }
   const int nnch = mLUTHeader[ipdg]->nchmap.nbins;
   const int nrad = mLUTHeader[ipdg]->radmap.nbins;
   const int neta = mLUTHeader[ipdg]->etamap.nbins;
@@ -106,12 +110,12 @@ TrackSmearer::smearTrack(O2Track &o2track, lutEntry_t *lutEntry)
 /*****************************************************************/
 
 bool
-TrackSmearer::smearTrack(O2Track &o2track, int pid)
+TrackSmearer::smearTrack(O2Track &o2track, int pid, float nch)
 {
 
   auto pt = o2track.getPt();
   auto eta = o2track.getEta();
-  auto lutEntry = getLUTEntry(pid, 0., 0., eta, pt);
+  auto lutEntry = getLUTEntry(pid, nch, 0., eta, pt);
   if (!lutEntry || !lutEntry->valid) return false;
   return smearTrack(o2track, lutEntry);
 }
@@ -124,7 +128,9 @@ TrackSmearer::smearTrack(Track &track, bool atDCA)
 
   O2Track o2track;
   TrackUtils::convertTrackToO2Track(track, o2track, atDCA);
-  if (!smearTrack(o2track, track.PID)) return false;
+  int pdg = track.PID;
+  float nch = mdNdEta; // use locally stored dNch/deta for the time being
+  if (!smearTrack(o2track, pdg, nch)) return false;
   TrackUtils::convertO2TrackToTrack(o2track, track, atDCA);
   return true;
 

src/TrackSmearer.hh

@@ -25,10 +25,11 @@ public:
   /** LUT methods **/
   bool loadTable(int pdg, const char *filename);
   void useEfficiency(bool val) { mUseEfficiency = val; };
+  lutHeader_t *getLUTHeader(int pdg) { return mLUTHeader[getIndexPDG(pdg)]; };
   lutEntry_t *getLUTEntry(int pdg, float nch, float radius, float eta, float pt);
 
   bool smearTrack(O2Track &o2track, lutEntry_t *lutEntry);
-  bool smearTrack(O2Track &o2track, int pid);
+  bool smearTrack(O2Track &o2track, int pid, float nch);
   bool smearTrack(Track &track, bool atDCA = true);
 
   int getIndexPDG(int pdg) {
@@ -41,12 +42,15 @@ public:
     default: return 2;
     };
   };
+
+  void setdNdEta(float val) { mdNdEta = val; };
 
 protected:
 
   lutHeader_t *mLUTHeader[5] = {nullptr};
   lutEntry_t *****mLUTEntry[5] = {nullptr};
   bool mUseEfficiency = true;
+  float mdNdEta =  1600.;
 
 };
 

src/lutCovm.hh

@@ -2,6 +2,7 @@
 /// @email: [email protected]
 
 #pragma once
+#define LUTCOVM_VERSION 20210611
 
 struct map_t {
   int nbins = 1;
@@ -27,24 +28,30 @@ struct map_t {
 };
 
 struct lutHeader_t {
+  int   version = LUTCOVM_VERSION;
   int   pdg = 0;
   float mass = 0.;
   float field = 0.;
   map_t nchmap;
   map_t radmap;
   map_t etamap;
   map_t ptmap;
+  bool check_version() {
+    return (version == LUTCOVM_VERSION);
+  };
   void print() {
-    printf("    pdg: %d \n", pdg);
-    printf("  field: %f \n", field);
-    printf(" nchmap: "); nchmap.print();
-    printf(" radmap: "); radmap.print();
-    printf(" etamap: "); etamap.print();
-    printf("  ptmap: "); ptmap.print();
+    printf(" version: %d \n", version);
+    printf("     pdg: %d \n", pdg);
+    printf("   field: %f \n", field);
+    printf("  nchmap: "); nchmap.print();
+    printf("  radmap: "); radmap.print();
+    printf("  etamap: "); etamap.print();
+    printf("   ptmap: "); ptmap.print();
   };
 };
 
 struct lutEntry_t {
+  float nch = 0.;
   float eta = 0.;
   float pt = 0.;
   bool valid = false;

src/lutWrite.cc

@@ -68,19 +68,19 @@ lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2,
   lutHeader.field = field;
   // nch
   lutHeader.nchmap.log   = true;
-  lutHeader.nchmap.nbins = 1;
-  lutHeader.nchmap.min   = 0.;
-  lutHeader.nchmap.max   = 4.;
+  lutHeader.nchmap.nbins = 20;
+  lutHeader.nchmap.min   = 0.5;
+  lutHeader.nchmap.max   = 3.5;
   // radius
   lutHeader.radmap.log   = false;
   lutHeader.radmap.nbins = 1;
   lutHeader.radmap.min   = 0.;
   lutHeader.radmap.max   = 100.;
   // eta
   lutHeader.etamap.log   = false;
-  lutHeader.etamap.nbins = 80;
-  lutHeader.etamap.min   = -4.;
-  lutHeader.etamap.max   =  4.;
+  lutHeader.etamap.nbins = 40;
+  lutHeader.etamap.min   = -2.;
+  lutHeader.etamap.max   =  2.;
   // pt
   lutHeader.ptmap.log    = true;
   lutHeader.ptmap.nbins  = 200;
@@ -96,8 +96,12 @@ lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2,
   lutEntry_t lutEntry;
 
   // write entries
-  for (int inch = 0; inch < nnch; ++inch)
-    for (int irad = 0; irad < nrad; ++irad)
+  for (int inch = 0; inch < nnch; ++inch) {
+    auto nch = lutHeader.nchmap.eval(inch);
+    lutEntry.nch = nch;
+    fat.SetdNdEtaCent(nch);
+    std::cout << " --- setting FAT dN/deta: " << nch << std::endl;
+    for (int irad = 0; irad < nrad; ++irad) {
       for (int ieta = 0; ieta < neta; ++ieta) {
 	auto eta = lutHeader.etamap.eval(ieta);
 	lutEntry.eta = lutHeader.etamap.eval(ieta);
@@ -125,7 +129,7 @@ lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2,
 	  }
 	  diagonalise(lutEntry);
 	  lutFile.write(reinterpret_cast<char *>(&lutEntry), sizeof(lutEntry_t));
-	}}
+	}}}}
 
 
   lutFile.close();

src/lutWrite.geometry_v1.cc

@@ -9,7 +9,6 @@ void
 fatInit_geometry_v1(float field = 0.5, float rmin = 100.)
 {
   fat.SetBField(field);
-  fat.SetdNdEtaCent(1600.);
   // new ideal Pixel properties?
   Double_t x0IB     = 0.001;
   Double_t x0OB     = 0.01;

src/lutWrite.werner.cc

@@ -9,7 +9,6 @@ void
 fatInit_werner(float field = 0.5, float rmin = 100.)
 {
   fat.SetBField(field);
-  fat.SetdNdEtaCent(1600.);
   // new ideal Pixel properties?
   Double_t x0IB     = 0.001;
   Double_t x0OB     = 0.01;