Centrality estimation using SPD tracklets (#282)

victor-gonzalez · web-flow · commit 9be98dc802e2 · 2021-11-18T15:32:20.000+01:00
* Centrality estimation using SPD tracklets

* Implementation of SPD Clusters, CL0 and CL1 estimators
diff --git a/Common/DataModel/Centrality.h b/Common/DataModel/Centrality.h
@@ -19,15 +19,18 @@ namespace cent
 {
 DECLARE_SOA_COLUMN(CentEstV0M, centV0M, float);                           //! Centrality percentile estimated from V0C+V0A multiplicities
 DECLARE_SOA_COLUMN(CentEstRun2SPDTracklets, centRun2SPDTracklets, float); //! Run2 centrality percentile estimated from SPD tracklets multiplicity
+DECLARE_SOA_COLUMN(CentEstRun2SPDClusters, centRun2SPDClusters, float);   //! Run2 centrality percentile estimated from SPD clusters multiplicity
 DECLARE_SOA_COLUMN(CentEstRun2CL0, centRun2CL0, float);                   //! Run2 centrality percentile estimated from CL0 multiplicity
 DECLARE_SOA_COLUMN(CentEstRun2CL1, centRun2CL1, float);                   //! Run2 centrality percentile estimated from CL1 multiplicity
 } // namespace cent
-DECLARE_SOA_TABLE(CentV0Ms, "AOD", "CENTV0M", cent::CentEstV0M);                      //! V0M estimated centrality table
-DECLARE_SOA_TABLE(CentRun2SPDs, "AOD", "CENTRUN2SPD", cent::CentEstRun2SPDTracklets); //! Run2 SPD tracklets estimated centrality table
-DECLARE_SOA_TABLE(CentRun2CL0s, "AOD", "CENTRUN2CL0", cent::CentEstRun2CL0);          //! Run2 CL0 estimated centrality table
-DECLARE_SOA_TABLE(CentRun2CL1s, "AOD", "CENTRUN2CL1", cent::CentEstRun2CL1);          //! Run2 CL1 estimated centrality table
+DECLARE_SOA_TABLE(CentV0Ms, "AOD", "CENTV0M", cent::CentEstV0M);                            //! V0M estimated centrality table
+DECLARE_SOA_TABLE(CentRun2SPDTrks, "AOD", "CENTRUN2SPDTRK", cent::CentEstRun2SPDTracklets); //! Run2 SPD tracklets estimated centrality table
+DECLARE_SOA_TABLE(CentRun2SPDClss, "AOD", "CENTRUN2SPDCLS", cent::CentEstRun2SPDClusters);  //! Run2 SPD clusters estimated centrality table
+DECLARE_SOA_TABLE(CentRun2CL0s, "AOD", "CENTRUN2CL0", cent::CentEstRun2CL0);                //! Run2 CL0 estimated centrality table
+DECLARE_SOA_TABLE(CentRun2CL1s, "AOD", "CENTRUN2CL1", cent::CentEstRun2CL1);                //! Run2 CL1 estimated centrality table
 using CentV0M = CentV0Ms::iterator;
-using CentRun2SPD = CentRun2SPDs::iterator;
+using CentRun2SPDTrk = CentRun2SPDTrks::iterator;
+using CentRun2SPDCls = CentRun2SPDClss::iterator;
 using CentRun2CL0 = CentRun2CL0s::iterator;
 using CentRun2CL1 = CentRun2CL1s::iterator;
 } // namespace o2::aod
diff --git a/Common/TableProducer/centralityTable.cxx b/Common/TableProducer/centralityTable.cxx
@@ -26,23 +26,49 @@ using namespace o2::framework;
 
 struct CentralityTable {
   Produces<aod::CentV0Ms> centVOM;
-  Produces<aod::CentRun2SPDs> centRun2SPD;
-  Produces<aod::CentRun2CL0s> centRun2CL1;
-  Produces<aod::CentRun2CL1s> centRun2CL2;
+  Produces<aod::CentRun2SPDTrks> centRun2SPDTracklets;
+  Produces<aod::CentRun2SPDClss> centRun2SPDClusters;
+  Produces<aod::CentRun2CL0s> centRun2CL0;
+  Produces<aod::CentRun2CL1s> centRun2CL1;
   Service<o2::ccdb::BasicCCDBManager> ccdb;
 
   Configurable<int> estV0M{"estV0M", -1, {"Produces centrality percentiles using V0 multiplicity. -1: auto, 0: don't, 1: yes. Default: auto (-1)"}};
-  Configurable<int> estRun2SPD{"estRun2SPD", -1, {"Produces Run2 centrality percentiles using SPD tracklets multiplicity. -1: auto, 0: don't, 1: yes. Default: auto (-1)"}};
+  Configurable<int> estRun2SPDTrklets{"estRun2SPDtks", -1, {"Produces Run2 centrality percentiles using SPD tracklets multiplicity. -1: auto, 0: don't, 1: yes. Default: auto (-1)"}};
+  Configurable<int> estRun2SPDClusters{"estRun2SPDcls", -1, {"Produces Run2 centrality percentiles using SPD clusters multiplicity. -1: auto, 0: don't, 1: yes. Default: auto (-1)"}};
   Configurable<int> estRun2CL0{"estRun2CL0", -1, {"Produces Run2 centrality percentiles using CL0 multiplicity. -1: auto, 0: don't, 1: yes. Default: auto (-1)"}};
   Configurable<int> estRun2CL1{"estRun2CL1", -1, {"Produces Run2 centrality percentiles using CL1 multiplicity. -1: auto, 0: don't, 1: yes. Default: auto (-1)"}};
 
   int mRunNumber;
-  bool mV0MCalibrationStored;
-  TH1* mhVtxAmpCorrV0A;
-  TH1* mhVtxAmpCorrV0C;
-  TH1* mhMultSelCalibV0M;
+  struct tagV0MCalibration {
+    bool mCalibrationStored = false;
+    TH1* mhVtxAmpCorrV0A = nullptr;
+    TH1* mhVtxAmpCorrV0C = nullptr;
+    TH1* mhMultSelCalib = nullptr;
+  } V0MInfo;
+  struct tagSPDTrackletsCalibration {
+    bool mCalibrationStored = false;
+    TH1* mhVtxAmpCorr = nullptr;
+    TH1* mhMultSelCalib = nullptr;
+  } SPDTksInfo;
+  struct tagSPDClustersCalibration {
+    bool mCalibrationStored = false;
+    TH1* mhVtxAmpCorrCL0 = nullptr;
+    TH1* mhVtxAmpCorrCL1 = nullptr;
+    TH1* mhMultSelCalib = nullptr;
+  } SPDClsInfo;
+  struct tagCL0Calibration {
+    bool mCalibrationStored = false;
+    TH1* mhVtxAmpCorr = nullptr;
+    TH1* mhMultSelCalib = nullptr;
+  } CL0Info;
+  struct tagCL1Calibration {
+    bool mCalibrationStored = false;
+    TH1* mhVtxAmpCorr = nullptr;
+    TH1* mhMultSelCalib = nullptr;
+  } CL1Info;
 
-  void init(InitContext& context)
+  void
+    init(InitContext& context)
   {
     /* Checking the tables which are requested in the workflow and enabling them */
     auto& workflows = context.services().get<RunningWorkflowInfo const>();
@@ -63,7 +89,8 @@ struct CentralityTable {
           }
         };
         enable("V0M", estV0M);
-        enable("Run2SPD", estRun2SPD);
+        enable("Run2SPDTrk", estRun2SPDTrklets);
+        enable("Run2SPDCls", estRun2SPDClusters);
         enable("Run2CL0", estRun2CL0);
         enable("Run2CL1", estRun2CL1);
       }
@@ -72,64 +99,134 @@ struct CentralityTable {
     ccdb->setCaching(true);
     ccdb->setLocalObjectValidityChecking();
     mRunNumber = 0;
-    mV0MCalibrationStored = false;
-    mhVtxAmpCorrV0A = nullptr;
-    mhVtxAmpCorrV0C = nullptr;
-    mhMultSelCalibV0M = nullptr;
   }
 
-  void process(soa::Join<aod::Collisions, aod::Mults>::iterator const& collision, aod::BCsWithTimestamps const&, aod::Tracks const& tracks)
+  using BCsWithTimestampsAndRun2Infos = soa::Join<aod::BCs, aod::Run2BCInfos, aod::Timestamps>;
+
+  void process(soa::Join<aod::Collisions, aod::Mults>::iterator const& collision, BCsWithTimestampsAndRun2Infos const&)
   {
     /* check the previous run number */
-    auto bc = collision.bc_as<aod::BCsWithTimestamps>();
+    auto bc = collision.bc_as<BCsWithTimestampsAndRun2Infos>();
     if (bc.runNumber() != mRunNumber) {
       LOGF(debug, "timestamp=%llu", bc.timestamp());
       TList* callst = ccdb->getForTimeStamp<TList>("Centrality/Estimators", bc.timestamp());
 
+      V0MInfo.mCalibrationStored = false;
+      SPDTksInfo.mCalibrationStored = false;
+      SPDClsInfo.mCalibrationStored = false;
+      CL0Info.mCalibrationStored = false;
+      CL1Info.mCalibrationStored = false;
       if (callst != nullptr) {
         auto getccdb = [callst](const char* ccdbhname) {
           TH1* h = (TH1*)callst->FindObject(ccdbhname);
           return h;
         };
         if (estV0M == 1) {
           LOGF(debug, "Getting new histograms with %d run number for %d run number", mRunNumber, bc.runNumber());
-          mhVtxAmpCorrV0A = getccdb("hVtx_fAmplitude_V0A_Normalized");
-          mhVtxAmpCorrV0C = getccdb("hVtx_fAmplitude_V0C_Normalized");
-          mhMultSelCalibV0M = getccdb("hMultSelCalib_V0M");
-          if ((mhVtxAmpCorrV0A != nullptr) and (mhVtxAmpCorrV0C != nullptr) and (mhMultSelCalibV0M != nullptr)) {
-            mV0MCalibrationStored = true;
+          V0MInfo.mhVtxAmpCorrV0A = getccdb("hVtx_fAmplitude_V0A_Normalized");
+          V0MInfo.mhVtxAmpCorrV0C = getccdb("hVtx_fAmplitude_V0C_Normalized");
+          V0MInfo.mhMultSelCalib = getccdb("hMultSelCalib_V0M");
+          if ((V0MInfo.mhVtxAmpCorrV0A != nullptr) and (V0MInfo.mhVtxAmpCorrV0C != nullptr) and (V0MInfo.mhMultSelCalib != nullptr)) {
+            V0MInfo.mCalibrationStored = true;
           } else {
-            LOGF(fatal, "Calibration information from V0M for run %d corrupted");
+            LOGF(fatal, "Calibration information from V0M for run %d corrupted", bc.runNumber());
           }
         }
-        if (estRun2SPD == 1) {
-          LOGF(fatal, "Run2 calibration information estimated from SPD tracklets still not available");
+        if (estRun2SPDTrklets == 1) {
+          LOGF(debug, "Getting new histograms with %d run number for %d run number", mRunNumber, bc.runNumber());
+          SPDTksInfo.mhVtxAmpCorr = getccdb("hVtx_fnTracklets_Normalized");
+          SPDTksInfo.mhMultSelCalib = getccdb("hMultSelCalib_SPDTracklets");
+          if ((SPDTksInfo.mhVtxAmpCorr != nullptr) and (SPDTksInfo.mhMultSelCalib != nullptr)) {
+            SPDTksInfo.mCalibrationStored = true;
+          } else {
+            LOGF(fatal, "Calibration information from SPD tracklets for run %d corrupted", bc.runNumber());
+          }
+        }
+        if (estRun2SPDClusters == 1) {
+          LOGF(debug, "Getting new histograms with %d run number for %d run number", mRunNumber, bc.runNumber());
+          SPDClsInfo.mhVtxAmpCorrCL0 = getccdb("hVtx_fnSPDClusters0_Normalized");
+          SPDClsInfo.mhVtxAmpCorrCL1 = getccdb("hVtx_fnSPDClusters1_Normalized");
+          SPDClsInfo.mhMultSelCalib = getccdb("hMultSelCalib_SPDClusters");
+          if ((SPDClsInfo.mhVtxAmpCorrCL0 != nullptr) and (SPDClsInfo.mhVtxAmpCorrCL1 != nullptr) and (SPDClsInfo.mhMultSelCalib != nullptr)) {
+            SPDClsInfo.mCalibrationStored = true;
+          } else {
+            LOGF(fatal, "Calibration information from SPD clusters for run %d corrupted", bc.runNumber());
+          }
         }
         if (estRun2CL0 == 1) {
-          LOGF(fatal, "Run2 calibration information estimated from CL0 still not available");
+          LOGF(debug, "Getting new histograms with %d run number for %d run number", mRunNumber, bc.runNumber());
+          CL0Info.mhVtxAmpCorr = getccdb("hVtx_fnSPDClusters0_Normalized");
+          CL0Info.mhMultSelCalib = getccdb("hMultSelCalib_CL0");
+          if ((CL0Info.mhVtxAmpCorr != nullptr) and (CL0Info.mhMultSelCalib != nullptr)) {
+            CL0Info.mCalibrationStored = true;
+          } else {
+            LOGF(fatal, "Calibration information from CL0 multiplicity for run %d corrupted", bc.runNumber());
+          }
         }
         if (estRun2CL1 == 1) {
-          LOGF(fatal, "Run2 calibration information estimated from CL1 still not available");
+          LOGF(debug, "Getting new histograms with %d run number for %d run number", mRunNumber, bc.runNumber());
+          CL1Info.mhVtxAmpCorr = getccdb("hVtx_fnSPDClusters1_Normalized");
+          CL1Info.mhMultSelCalib = getccdb("hMultSelCalib_CL1");
+          if ((CL1Info.mhVtxAmpCorr != nullptr) and (CL1Info.mhMultSelCalib != nullptr)) {
+            CL1Info.mCalibrationStored = true;
+          } else {
+            LOGF(fatal, "Calibration information from CL1 multiplicity for run %d corrupted", bc.runNumber());
+          }
         }
-        if (mV0MCalibrationStored) {
+        if (V0MInfo.mCalibrationStored or SPDTksInfo.mCalibrationStored or SPDClsInfo.mCalibrationStored or CL0Info.mCalibrationStored or CL1Info.mCalibrationStored) {
           mRunNumber = bc.runNumber();
         }
       } else {
-        /* we dont change the run number to keep trying */
-        mV0MCalibrationStored = false;
         LOGF(fatal, "Centrality calibration is not available in CCDB for run=%d at timestamp=%llu", bc.runNumber(), bc.timestamp());
       }
     }
     if (estV0M == 1) {
-      float centV0M = 105.0f;
-      if (mV0MCalibrationStored) {
-        float v0m = collision.multV0A() * mhVtxAmpCorrV0A->GetBinContent(mhVtxAmpCorrV0A->FindFixBin(collision.posZ())) +
-                    collision.multV0C() * mhVtxAmpCorrV0C->GetBinContent(mhVtxAmpCorrV0C->FindFixBin(collision.posZ()));
-        centV0M = mhMultSelCalibV0M->GetBinContent(mhMultSelCalibV0M->FindFixBin(v0m));
+      float cV0M = 105.0f;
+      if (V0MInfo.mCalibrationStored) {
+        float v0m = collision.multV0A() * V0MInfo.mhVtxAmpCorrV0A->GetBinContent(V0MInfo.mhVtxAmpCorrV0A->FindFixBin(collision.posZ())) +
+                    collision.multV0C() * V0MInfo.mhVtxAmpCorrV0C->GetBinContent(V0MInfo.mhVtxAmpCorrV0C->FindFixBin(collision.posZ()));
+        cV0M = V0MInfo.mhMultSelCalib->GetBinContent(V0MInfo.mhMultSelCalib->FindFixBin(v0m));
       }
-      LOGF(debug, "centV0M=%.0f", centV0M);
+      LOGF(debug, "centV0M=%.0f", cV0M);
       // fill centrality columns
-      centVOM(centV0M);
+      centVOM(cV0M);
+    }
+    if (estRun2SPDTrklets == 1) {
+      float cSPD = 105.0f;
+      if (SPDTksInfo.mCalibrationStored) {
+        float spdm = collision.multTracklets() * SPDTksInfo.mhVtxAmpCorr->GetBinContent(SPDTksInfo.mhVtxAmpCorr->FindFixBin(collision.posZ()));
+        cSPD = SPDTksInfo.mhMultSelCalib->GetBinContent(SPDTksInfo.mhMultSelCalib->FindFixBin(spdm));
+      }
+      LOGF(debug, "centSPDTracklets=%.0f", cSPD);
+      centRun2SPDTracklets(cSPD);
+    }
+    if (estRun2SPDClusters == 1) {
+      float cSPD = 105.0f;
+      if (SPDClsInfo.mCalibrationStored) {
+        float spdm = bc.spdClustersL0() * SPDClsInfo.mhVtxAmpCorrCL0->GetBinContent(SPDClsInfo.mhVtxAmpCorrCL0->FindFixBin(collision.posZ())) +
+                     bc.spdClustersL1() * SPDClsInfo.mhVtxAmpCorrCL1->GetBinContent(SPDClsInfo.mhVtxAmpCorrCL1->FindFixBin(collision.posZ()));
+        cSPD = SPDClsInfo.mhMultSelCalib->GetBinContent(SPDClsInfo.mhMultSelCalib->FindFixBin(spdm));
+      }
+      LOGF(debug, "centSPDClusters=%.0f", cSPD);
+      centRun2SPDClusters(cSPD);
+    }
+    if (estRun2CL0 == 1) {
+      float cCL0 = 105.0f;
+      if (CL0Info.mCalibrationStored) {
+        float cl0m = bc.spdClustersL0() * CL0Info.mhVtxAmpCorr->GetBinContent(CL0Info.mhVtxAmpCorr->FindFixBin(collision.posZ()));
+        cCL0 = CL0Info.mhMultSelCalib->GetBinContent(CL0Info.mhMultSelCalib->FindFixBin(cl0m));
+      }
+      LOGF(debug, "centCL0=%.0f", cCL0);
+      centRun2CL0(cCL0);
+    }
+    if (estRun2CL1 == 1) {
+      float cCL1 = 105.0f;
+      if (CL1Info.mCalibrationStored) {
+        float cl1m = bc.spdClustersL1() * CL1Info.mhVtxAmpCorr->GetBinContent(CL1Info.mhVtxAmpCorr->FindFixBin(collision.posZ()));
+        cCL1 = CL1Info.mhMultSelCalib->GetBinContent(CL1Info.mhMultSelCalib->FindFixBin(cl1m));
+      }
+      LOGF(debug, "centCL1=%.0f", cCL1);
+      centRun2CL1(cCL1);
     }
   }
 };
diff --git a/Common/Tasks/centralityQa.cxx b/Common/Tasks/centralityQa.cxx
@@ -19,8 +19,12 @@ using namespace o2;
 using namespace o2::framework;
 
 struct CentralityQa {
-  OutputObj<TH1F> hCentV0M{TH1F("hCentV0M", "", 21, 0, 105.)};
-  void process(soa::Join<aod::Collisions, aod::EvSels, aod::CentV0Ms>::iterator const& col)
+  OutputObj<TH1F> hCentV0M{TH1F("hCentV0M", "V0M", 21, 0, 105.)};
+  OutputObj<TH1F> hCentSPDTks{TH1F("hCentSPDTks", "SPD Tracklets", 21, 0, 105.)};
+  OutputObj<TH1F> hCentSPDCls{TH1F("hCentSPDCls", "SPD Clusters", 21, 0, 105.)};
+  OutputObj<TH1F> hCentCL0{TH1F("hCentCL0", "CL0", 21, 0, 105.)};
+  OutputObj<TH1F> hCentCL1{TH1F("hCentCL1", "CL1", 21, 0, 105.)};
+  void processPP(soa::Join<aod::Collisions, aod::EvSels, aod::CentV0Ms, aod::CentRun2SPDTrks, aod::CentRun2SPDClss>::iterator const& col)
   {
     if (!col.alias()[kINT7]) {
       return;
@@ -30,9 +34,35 @@ struct CentralityQa {
     }
 
     LOGF(debug, "centV0M=%.0f", col.centV0M());
+    LOGF(debug, "centSPDTracklets=%.0f", col.centRun2SPDTracklets());
+    LOGF(debug, "centSPDClusters=%.0f", col.centRun2SPDClusters());
     // fill centrality histos
     hCentV0M->Fill(col.centV0M());
+    hCentSPDTks->Fill(col.centRun2SPDTracklets());
+    hCentSPDCls->Fill(col.centRun2SPDClusters());
   }
+  PROCESS_SWITCH(CentralityQa, processPP, "Process with SPD clusters centrality/multiplicity estimation", false);
+
+  void processPbPb(soa::Join<aod::Collisions, aod::EvSels, aod::CentV0Ms, aod::CentRun2SPDTrks, aod::CentRun2CL0s, aod::CentRun2CL1s>::iterator const& col)
+  {
+    if (!col.alias()[kINT7]) {
+      return;
+    }
+    if (!col.sel7()) {
+      return;
+    }
+
+    LOGF(debug, "centV0M=%.0f", col.centV0M());
+    LOGF(debug, "centSPDTracklets=%.0f", col.centRun2SPDTracklets());
+    LOGF(debug, "centCL0=%.0f", col.centRun2CL0());
+    LOGF(debug, "centCL1=%.0f", col.centRun2CL1());
+    // fill centrality histos
+    hCentV0M->Fill(col.centV0M());
+    hCentSPDTks->Fill(col.centRun2SPDTracklets());
+    hCentCL0->Fill(col.centRun2CL0());
+    hCentCL1->Fill(col.centRun2CL1());
+  }
+  PROCESS_SWITCH(CentralityQa, processPbPb, "Process with CL0 and CL1 multiplicities centrality/multiplicity  estimation", true);
 };
 
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
