Implement luminosity calculation in SG candidate producer

Author: nburmaso

PWGUD/TableProducer/SGCandProducer.cxx

@@ -21,6 +21,8 @@
 #include "DataFormatsFIT/Triggers.h"
 #include "DataFormatsParameters/GRPMagField.h"
 #include "DataFormatsParameters/GRPObject.h"
+#include "DataFormatsParameters/GRPLHCIFData.h"
+#include "DataFormatsParameters/AggregatedRunInfo.h"
 
 #include "Framework/AnalysisTask.h"
 #include "Framework/ASoAHelpers.h"
@@ -74,6 +76,7 @@ struct SGCandProducer {
   // Configurables to decide which tables are filled
   Configurable<bool> fillTrackTables{"fillTrackTables", true, "Fill track tables"};
   Configurable<bool> fillFwdTrackTables{"fillFwdTrackTables", true, "Fill forward track tables"};
+
   //  SG selector
   SGSelector sgSelector;
   ctpRateFetcher mRateFetcher;
@@ -103,6 +106,8 @@ struct SGCandProducer {
     {}};
   std::map<std::string, HistPtr> histPointers;
 
+  int runNumber = -1;
+
   // function to update UDFwdTracks, UDFwdTracksExtra
   template <typename TFwdTrack>
   void updateUDFwdTrackTables(TFwdTrack const& fwdtrack, uint64_t const& bcnum)
@@ -180,6 +185,78 @@ struct SGCandProducer {
     outputTracksLabel(track.globalIndex());
   }
 
+  // function to process trigger counters, accounting for BC selection bits
+  void processCountersTrg(BCs const& bcs, aod::FT0s const&, aod::Zdcs const&)
+  {
+    const auto& firstBc = bcs.iteratorAt(0);
+    if (runNumber != firstBc.runNumber())
+      runNumber = firstBc.runNumber();
+
+    auto hCountersTrg = getHist(TH1, "reco/hCountersTrg");
+    auto hCountersTrgBcSel = getHist(TH1, "reco/hCountersTrgBcSel");
+    auto hLumi = getHist(TH1, "reco/hLumi");
+    auto hLumiBcSel = getHist(TH1, "reco/hLumiBcSel");
+
+    // Cross sections in ub. Using dummy -1 if lumi estimator is not reliable
+    float csTCE = 10.36e6;
+    float csZEM = 415.2e6; // see AN: https://alice-notes.web.cern.ch/node/1515
+    float csZNC = 214.5e6; // see AN: https://alice-notes.web.cern.ch/node/1515
+    if (runNumber > 543437 && runNumber < 543514) {
+      csTCE = 8.3e6;
+    }
+    if (runNumber >= 543514) {
+      csTCE = 4.10e6; // see AN: https://alice-notes.web.cern.ch/node/1515
+    }
+
+    for (const auto& bc : bcs) {
+      bool hasFT0 = bc.has_foundFT0();
+      bool hasZDC = bc.has_foundZDC();
+      if (!hasFT0 && !hasZDC)
+        continue;
+      bool isSelectedBc = true;
+      if (rejectAtTFBoundary && !bc.selection_bit(aod::evsel::kNoTimeFrameBorder))
+        isSelectedBc = false;
+      if (noITSROFrameBorder && !bc.selection_bit(aod::evsel::kNoITSROFrameBorder))
+        isSelectedBc = false;
+      if (hasFT0) {
+        auto ft0TrgMask = bc.ft0().triggerMask();
+        if (TESTBIT(ft0TrgMask, o2::fit::Triggers::bitVertex)) {
+          hCountersTrg->Fill("TVX", 1);
+          if (isSelectedBc)
+            hCountersTrgBcSel->Fill("TVX", 1);
+        }
+        if (TESTBIT(ft0TrgMask, o2::fit::Triggers::bitVertex) && TESTBIT(ft0TrgMask, o2::fit::Triggers::bitCen)) {
+          hCountersTrg->Fill("TCE", 1);
+          hLumi->Fill("TCE", 1. / csTCE);
+          if (isSelectedBc) {
+            hCountersTrgBcSel->Fill("TCE", 1);
+            hLumiBcSel->Fill("TCE", 1. / csTCE);
+          }
+        }
+      }
+      if (hasZDC) {
+        if (bc.selection_bit(aod::evsel::kIsBBZNA) || bc.selection_bit(aod::evsel::kIsBBZNC)) {
+          hCountersTrg->Fill("ZEM", 1);
+          hLumi->Fill("ZEM", 1. / csZEM);
+          if (isSelectedBc) {
+            hCountersTrgBcSel->Fill("ZEM", 1);
+            hLumiBcSel->Fill("ZEM", 1. / csZEM);
+          }
+        }
+        if (bc.selection_bit(aod::evsel::kIsBBZNC)) {
+          hCountersTrg->Fill("ZNC", 1);
+          hLumi->Fill("ZNC", 1. / csZNC);
+          if (isSelectedBc) {
+            hCountersTrgBcSel->Fill("ZNC", 1);
+            hLumiBcSel->Fill("ZNC", 1. / csZNC);
+          }
+        }
+      }
+    }
+  }
+
+  PROCESS_SWITCH(SGCandProducer, processCountersTrg, "Produce trigger counters and luminosity histograms", true);
+
   // function to process reconstructed data
   template <typename TCol>
   void processReco(std::string histdir, TCol const& collision, BCs const& bcs,
@@ -324,6 +401,22 @@ struct SGCandProducer {
     histPointers.clear();
     if (context.mOptions.get<bool>("processData")) {
       histPointers.insert({"reco/Stat", registry.add("reco/Stat", "Cut statistics; Selection criterion; Collisions", {HistType::kTH1F, {{14, -0.5, 13.5}}})});
+
+      const AxisSpec axisCountersTrg{10, 0.5, 10.5, ""};
+      histPointers.insert({"reco/hCountersTrg", registry.add("reco/hCountersTrg", "Trigger counts before selections; Trigger; Counts", {HistType::kTH1F, {axisCountersTrg}})});
+      histPointers.insert({"reco/hCountersTrgBcSel", registry.add("reco/hCountersTrgSel", "Trigger counts after BC selections; Trigger; Counts", {HistType::kTH1F, {axisCountersTrg}})});
+      histPointers.insert({"reco/hLumi", registry.add("reco/hLumi", "Integrated luminosity before selections; Trigger; Luminosity, 1/#mub", {HistType::kTH1F, {axisCountersTrg}})});
+      histPointers.insert({"reco/hLumiBcSel", registry.add("reco/hLumiBcSel", "Integrated luminosity before selections; Trigger; Luminosity, 1/#mub", {HistType::kTH1F, {axisCountersTrg}})});
+      auto hCountersTrg = getHist(TH1, "reco/hCountersTrg");
+      auto hCountersTrgBcSel = getHist(TH1, "reco/hCountersTrgBcSel");
+      auto hLumi = getHist(TH1, "reco/hLumi");
+      auto hLumiBcSel = getHist(TH1, "reco/hLumiBcSel");
+      for (auto h : {hCountersTrg, hCountersTrgBcSel, hLumi, hLumiBcSel}) {
+        h->GetXaxis()->SetBinLabel(1, "TVX");
+        h->GetXaxis()->SetBinLabel(2, "TCE");
+        h->GetXaxis()->SetBinLabel(3, "ZEM");
+        h->GetXaxis()->SetBinLabel(4, "ZNC");
+      }
     }
     if (context.mOptions.get<bool>("processMcData")) {
       histPointers.insert({"MCreco/Stat", registry.add("MCreco/Stat", "Cut statistics; Selection criterion; Collisions", {HistType::kTH1F, {{14, -0.5, 13.5}}})});