[ALICE3] A3: update protection of DetLayer

ALICE3/Core/DelphesO2LutWriter.cxx

@@ -26,9 +26,10 @@
 #include "iostream"
 #include "TMatrixD.h"
 #include "TVectorD.h"
+#include "TAxis.h"
 #include "TMatrixDSymEigen.h"
 #include "TDatabasePDG.h"
 #include "TLorentzVector.h"
 #include "ALICE3/Core/FastTracker.h"
 #include "ALICE3/Core/TrackUtilities.h"
 
@@ -55,17 +56,18 @@
                                   const float mass,
                                   int itof,
                                   int otof,
-                                  int q)
+                                  int q,
+                                  const float nch)
 {
   lutEntry.valid = false;
 
   static TLorentzVector tlv;
   tlv.SetPtEtaPhiM(pt, eta, 0., mass);
   o2::track::TrackParCov trkIn;
   o2::upgrade::convertTLorentzVectorToO2Track(q, tlv, {0., 0., 0.}, trkIn);
-
   o2::track::TrackParCov trkOut;
-  if (fat.FastTrack(trkIn, trkOut, 1) < 0) {
+  const int status = fat.FastTrack(trkIn, trkOut, nch);
+  if (status <= 0) {
     Printf(" --- fatSolve: FastTrack failed --- \n");
     tlv.Print();
     return false;
@@ -142,7 +144,7 @@
   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 double relmomres_pos = 10e-6 * pt / 0.3 / Bfield / Leta / Leta * std::sqrt(720. / 15.);
 
   const float relmomres_barrel = std::sqrt(covmbarrel[14]) * pt;
@@ -152,7 +154,7 @@
 
   // 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 momres_tot = pt * std::sqrt(relmomres_pos * relmomres_pos + relmomres_MS * relmomres_MS); // total absolute mom reso
 
   // Fill cov matrix diag
@@ -197,7 +199,7 @@
   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;
   if (q <= 0) {
     Printf("Negative or null charge (%f) for pdg code %i. Fix the charge!", TDatabasePDG::Instance()->GetParticle(pdg)->Charge(), pdg);
@@ -234,6 +236,9 @@
   lutEntry_t lutEntry;
 
   // write entries
+  int nCalls = 0;
+  int successfullCalls = 0;
+  int failedCalls = 0;
   for (int inch = 0; inch < nnch; ++inch) {
     Printf(" --- writing nch = %d/%d", inch, nnch);
     auto nch = lutHeader.nchmap.eval(inch);
@@ -242,6 +247,7 @@
     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);
@@ -252,6 +258,7 @@
           if (std::fabs(eta) <= etaMaxBarrel) { // full lever arm ends at etaMaxBarrel
             Printf("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)) {
               // printf(" --- fatSolve: error \n");
               lutEntry.valid = false;
@@ -260,13 +267,16 @@
               for (int i = 0; i < 15; ++i) {
                 lutEntry.covm[i] = 0.;
               }
+              successfullCalls--;
+              failedCalls++;
             }
           } else {
             Printf("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.;
             bool retval = true;
+            successfullCalls++;
             if (useFlatDipole) { // Using the parametrization at the border of the barrel
               retval = fatSolve(lutEntry, lutEntry.pt, etaMaxBarrel, lutHeader.mass, itof, otof, q);
             } else if (usePara) {
@@ -288,6 +298,8 @@
               for (int i = 0; i < 15; ++i) {
                 lutEntry.covm[i] = 0.;
               }
+              successfullCalls--;
+              failedCalls++;
             }
           }
           Printf("Diagonalizing");
@@ -298,6 +310,8 @@
       }
     }
   }
+  Printf(" --- finished writing LUT file %s", filename);
+  Printf(" --- successfull calls: %d/%d, failed calls: %d/%d", successfullCalls, nCalls, failedCalls, nCalls);
   lutFile.close();
 }
 
@@ -331,10 +345,13 @@
 
 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);
+  // vs
   static const int kNch = 0;
   static const int kEta = 1;
   static const int kPt = 2;
 
+  // what
   static const int kEfficiency = 0;
   static const int kEfficiency2 = 1;
   static const int kEfficiencyInnerTOF = 2;
@@ -360,6 +377,58 @@
   }
   auto nbins = lutMap.nbins;
   auto g = new TGraph();
+  g->SetName(Form("lut_%s_%d_vs_%d_what_%d", filename, pdg, vs, what));
+  g->SetTitle(Form("LUT for %s, pdg %d, vs %d, what %d", filename, pdg, vs, what));
+  switch (vs) {
+    case kNch:
+      Printf(" --- vs = kNch");
+      g->GetXaxis()->SetTitle("Nch");
+      break;
+    case kEta:
+      Printf(" --- vs = kEta");
+      g->GetXaxis()->SetTitle("#eta");
+      break;
+    case kPt:
+      Printf(" --- vs = kPt");
+      g->GetXaxis()->SetTitle("p_{T} (GeV/c)");
+      break;
+    default:
+      Printf(" --- error: unknown vs %d", vs);
+      return nullptr;
+  }
+  switch (what) {
+    case kEfficiency:
+      Printf(" --- what = kEfficiency");
+      g->GetYaxis()->SetTitle("Efficiency (%)");
+      break;
+    case kEfficiency2:
+      Printf(" --- what = kEfficiency2");
+      g->GetYaxis()->SetTitle("Efficiency2 (%)");
+      break;
+    case kEfficiencyInnerTOF:
+      Printf(" --- what = kEfficiencyInnerTOF");
+      g->GetYaxis()->SetTitle("Inner TOF Efficiency (%)");
+      break;
+    case kEfficiencyOuterTOF:
+      Printf(" --- what = kEfficiencyOuterTOF");
+      g->GetYaxis()->SetTitle("Outer TOF Efficiency (%)");
+      break;
+    case kPtResolution:
+      Printf(" --- what = kPtResolution");
+      g->GetYaxis()->SetTitle("p_{T} Resolution (%)");
+      break;
+    case kRPhiResolution:
+      Printf(" --- what = kRPhiResolution");
+      g->GetYaxis()->SetTitle("R#phi Resolution (#mum)");
+      break;
+    case kZResolution:
+      Printf(" --- what = kZResolution");
+      g->GetYaxis()->SetTitle("Z Resolution (#mum)");
+      break;
+    default:
+      Printf(" --- error: unknown what %d", what);
+      return nullptr;
+  }
 
   bool canBeInvalid = true;
   for (int i = 0; i < nbins; ++i) {
@@ -411,13 +480,13 @@
         val = lutEntry->otof * 100.; // efficiency (%)
         break;
       case kPtResolution:
         val = sqrt(lutEntry->covm[14]) * lutEntry->pt * 100.; // pt resolution (%)
         break;
       case kRPhiResolution:
         val = sqrt(lutEntry->covm[0]) * 1.e4; // rphi resolution (um)
         break;
       case kZResolution:
         val = sqrt(lutEntry->covm[1]) * 1.e4; // z resolution (um)
         break;
       default:
         Printf(" --- error: unknown what %d", what);

ALICE3/Core/DelphesO2LutWriter.h

@@ -55,7 +55,8 @@ class DelphesO2LutWriter
                 const float mass = 0.13957000,
                 int itof = 0,
                 int otof = 0,
-                int q = 1);
+                int q = 1,
+                const float nch = 1);
   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);

ALICE3/Core/DetLayer.cxx

@@ -16,4 +16,76 @@
 /// \brief  Basic struct to hold information regarding a detector layer to be used in fast simulation
 ///
 
+#include <vector>
+#include <string>
+
 #include "DetLayer.h"
+
+namespace o2::fastsim
+{
+
+// Parametric constructor
+DetLayer::DetLayer(const TString& name_,
+                   float r_,
+                   float z_,
+                   float x0_,
+                   float xrho_,
+                   float resRPhi_,
+                   float resZ_,
+                   float eff_,
+                   int type_)
+  : name(name_),
+    r(r_),
+    z(z_),
+    x0(x0_),
+    xrho(xrho_),
+    resRPhi(resRPhi_),
+    resZ(resZ_),
+    eff(eff_),
+    type(type_)
+{
+}
+
+DetLayer::DetLayer(const DetLayer& other)
+  : name(other.name), r(other.r), z(other.z), x0(other.x0), xrho(other.xrho), resRPhi(other.resRPhi), resZ(other.resZ), eff(other.eff), type(other.type)
+{
+}
+
+std::string DetLayer::toString() const
+{
+  std::string out = "";
+  out.append("DetLayer: ");
+  out.append(name.Data());
+  out.append(" | r: ");
+  out.append(std::to_string(r));
+  out.append(" cm | z: ");
+  out.append(std::to_string(z));
+  out.append(" cm | x0: ");
+  out.append(std::to_string(x0));
+  out.append(" cm | xrho: ");
+  out.append(std::to_string(xrho));
+  out.append(" g/cm^3 | resRPhi: ");
+  out.append(std::to_string(resRPhi));
+  out.append(" cm | resZ: ");
+  out.append(std::to_string(resZ));
+  out.append(" cm | eff: ");
+  out.append(std::to_string(eff));
+  out.append(" | type: ");
+  switch (type) {
+    case layerInert:
+      out.append("Inert");
+      break;
+    case layerSilicon:
+      out.append("Silicon");
+      break;
+    case layerGas:
+      out.append("Gas/TPC");
+      break;
+    default:
+      out.append("Unknown");
+      break;
+  }
+  return out;
+}
+
+} // namespace o2::fastsim

ALICE3/Core/DetLayer.h

@@ -19,12 +19,59 @@
 #ifndef ALICE3_CORE_DETLAYER_H_
 #define ALICE3_CORE_DETLAYER_H_
 
+#include <string>
+
 #include "TString.h"
 
 namespace o2::fastsim
 {
 
 struct DetLayer {
+ public:
+  // Default constructor
+  DetLayer() = default;
+  // Parametric constructor
+  DetLayer(const TString& name_, float r_, float z_, float x0_, float xrho_,
+           float resRPhi_ = 0.0f, float resZ_ = 0.0f, float eff_ = 0.0f, int type_ = layerInert);
+  // Copy constructor
+  DetLayer(const DetLayer& other);
+
+  // Setters
+  void setName(const TString& name_) { name = name_; }
+  void setRadius(float r_) { r = r_; }
+  void setZ(float z_) { z = z_; }
+  void setRadiationLength(float x0_) { x0 = x0_; }
+  void setDensity(float xrho_) { xrho = xrho_; }
+  void setResolutionRPhi(float resRPhi_) { resRPhi = resRPhi_; }
+  void setResolutionZ(float resZ_) { resZ = resZ_; }
+  void setEfficiency(float eff_) { eff = eff_; }
+  void setType(int type_) { type = type_; }
+
+  // Getters
+  float getRadius() const { return r; }
+  float getZ() const { return z; }
+  float getRadiationLength() const { return x0; }
+  float getDensity() const { return xrho; }
+  float getResolutionRPhi() const { return resRPhi; }
+  float getResolutionZ() const { return resZ; }
+  float getEfficiency() const { return eff; }
+  int getType() const { return type; }
+  const TString& getName() const { return name; }
+
+  // Check layer type
+  bool isInert() const { return type == layerInert; }
+  bool isSilicon() const { return type == layerSilicon; }
+  bool isGas() const { return type == layerGas; }
+
+  // Utilities
+  std::string toString() const;
+  friend std::ostream& operator<<(std::ostream& os, const DetLayer& layer)
+  {
+    os << layer.toString();
+    return os;
+  }
+
+ private:
   // TString for holding name
   TString name;
 
@@ -44,7 +91,10 @@ struct DetLayer {
   float eff; // detection efficiency
 
   // layer type
-  int type; // 0: undefined/inert, 1: silicon, 2: gas/tpc
+  int type;                              // 0: undefined/inert, 1: silicon, 2: gas/tpc
+  static constexpr int layerInert = 0;   // inert/undefined layer
+  static constexpr int layerSilicon = 1; // silicon layer
+  static constexpr int layerGas = 2;     // gas/tpc layer
 };
 
 } // namespace o2::fastsim