ITS: dROF study in CheckTracksCA

Signed-off-by: Felix Schlepper <felix.schlepper@cern.ch>

Author: f3sch

Detectors/ITSMFT/ITS/macros/test/CheckTracksCA.C

@@ -35,6 +35,8 @@
 #include "SimulationDataFormat/MCCompLabel.h"
 #include "SimulationDataFormat/MCTruthContainer.h"
 #include "DataFormatsITS/TrackITS.h"
+#include "DataFormatsITSMFT/ROFRecord.h"
+#include "SimulationDataFormat/DigitizationContext.h"
 
 #endif
 #include "DataFormatsITSMFT/CompCluster.h"
@@ -57,19 +59,29 @@ struct ParticleInfo {
   unsigned short clusters = 0u;
   unsigned char isReco = 0u;
   unsigned char isFake = 0u;
-  bool isPrimary = 0u;
+  bool isPrimary = false;
+  int bcInROF{-1};
+  int rofId{-1};
   unsigned char storedStatus = 2; /// not stored = 2, fake = 1, good = 0
   o2::its::TrackITS track;
+
+  void print() const
+  {
+    LOGP(info, "event={} pdg={} pt={} eta={} phi={} mother={} clusters={:7b} isReco={} isFake={} isPrimary={} bcInROF={} rofId={} | {}", event, pdg, pt, eta, phi, mother, clusters, isReco, isFake, isPrimary, bcInROF, rofId, track.asString());
+  }
 };
 
 #pragma link C++ class ParticleInfo + ;
 
 void CheckTracksCA(bool doFakeClStud = false,
+                   bool doDRofStud = true,
                    std::string tracfile = "o2trac_its.root",
                    std::string magfile = "o2sim_grp.root",
                    std::string clusfile = "o2clus_its.root",
                    std::string kinefile = "o2sim_Kine.root")
 {
+  constexpr int64_t roFrameLengthInBC = 198; // for pp=198
+  constexpr int64_t roFrameBiasInBC = 64;    // ITS delay accounted for in digitization
 
   using namespace o2::itsmft;
   using namespace o2::its;
@@ -94,6 +106,10 @@ void CheckTracksCA(bool doFakeClStud = false,
   o2::dataformats::MCEventHeader* mcEvent = nullptr;
   mcTree->SetBranchAddress("MCEventHeader.", &mcEvent);
 
+  auto* dc = o2::steer::DigitizationContext::loadFromFile("collisioncontext.root");
+  const auto& irs = dc->getEventRecords();
+  dc->printCollisionSummary(false, 20);
+
   // Clusters
   TFile::Open(clusfile.data());
   TTree* clusTree = (TTree*)gFile->Get("o2sim");
@@ -112,6 +128,8 @@ void CheckTracksCA(bool doFakeClStud = false,
   // Track MC labels
   std::vector<o2::MCCompLabel>* trkLabArr = nullptr;
   recTree->SetBranchAddress("ITSTrackMCTruth", &trkLabArr);
+  std::vector<o2::itsmft::ROFRecord> rofRecVec, *rofRecVecP = &rofRecVec;
+  recTree->SetBranchAddress("ITSTracksROF", &rofRecVecP);
 
   std::cout << "** Filling particle table ... " << std::flush;
   int lastEventIDcl = -1, cf = 0;
@@ -123,6 +141,7 @@ void CheckTracksCA(bool doFakeClStud = false,
     mcTree->GetEvent(n);
     info[n].resize(mcArr->size());
     hZvertex->Fill(mcEvent->GetZ());
+    const auto& ir = irs[mcEvent->GetEventID() - 1]; // event id start from 1
     for (unsigned int mcI{0}; mcI < mcArr->size(); ++mcI) {
       auto part = mcArr->at(mcI);
       info[n][mcI].event = n;
@@ -134,6 +153,10 @@ void CheckTracksCA(bool doFakeClStud = false,
       info[n][mcI].phi = part.GetPhi();
       info[n][mcI].eta = part.GetEta();
       info[n][mcI].isPrimary = part.isPrimary();
+      if (!ir.isDummy()) {
+        info[n][mcI].bcInROF = (ir.toLong() - roFrameBiasInBC) % roFrameLengthInBC;
+        info[n][mcI].rofId = (ir.toLong() - roFrameBiasInBC) / roFrameLengthInBC;
+      }
     }
   }
   std::cout << "done." << std::endl;
@@ -168,7 +191,7 @@ void CheckTracksCA(bool doFakeClStud = false,
   std::cout << "done." << std::endl;
 
   std::cout << "** Analysing tracks ... " << std::flush;
-  int unaccounted{0}, good{0}, fakes{0}, total{0};
+  int unaccounted{0}, good{0}, fakes{0}, total{0}, length{0};
   for (int frame = 0; frame < recTree->GetEntriesFast(); frame++) { // Cluster frames
     if (!recTree->GetEvent(frame))
       continue;
@@ -204,6 +227,11 @@ void CheckTracksCA(bool doFakeClStud = false,
 
       fakes += fake;
       good += !fake;
+      if (!fake) {
+        for (unsigned int bit{0}; bit < 7; ++bit) {
+          length += bool(info[evID][trackID].clusters & (1 << bit));
+        }
+      }
     }
   }
   std::cout << "done." << std::endl;
@@ -213,6 +241,7 @@ void CheckTracksCA(bool doFakeClStud = false,
   std::cout << "\t- Total number of tracks not corresponding to particles: " << unaccounted << " (" << unaccounted * 100. / total << "%)" << std::endl;
   std::cout << "\t- Total number of fakes: " << fakes << " (" << fakes * 100. / total << "%)" << std::endl;
   std::cout << "\t- Total number of good: " << good << " (" << good * 100. / total << "%)" << std::endl;
+  std::cout << "\t- Average length of good tracks: " << (double)length / (double)good << std::endl;
 
   int nb = 100;
   double xbins[nb + 1], ptcutl = 0.01, ptcuth = 10.;
@@ -271,13 +300,13 @@ void CheckTracksCA(bool doFakeClStud = false,
 
   TCanvas* c1 = new TCanvas;
   c1->SetLogx();
-  c1->SetGridx();
-  c1->SetGridy();
+  c1->SetGrid();
+  gPad->DrawFrame(ptcutl, 0.001, ptcuth, 1.03);
   TH1* sum = (TH1*)num->Clone("sum");
   sum->Add(fak);
   sum->Divide(sum, den, 1, 1);
   sum->SetLineColor(kBlack);
-  sum->Draw("hist");
+  sum->Draw("histesame");
   num->Divide(num, den, 1, 1, "b");
   num->Draw("histesame");
   fak->Divide(fak, den, 1, 1, "b");
@@ -410,4 +439,58 @@ void CheckTracksCA(bool doFakeClStud = false,
     }
     canvas->SaveAs("fakeClusters.png", "recreate");
   }
+
+  if (doDRofStud) {
+    constexpr float rofBins{roFrameLengthInBC + 2.f};
+    TH1F* hBC = new TH1F("hBC", "Distance in BC;bcInROF;counts.", rofBins, -0.5, rofBins - 0.5);
+    TH1F* hBCTracksDen = new TH1F("hBCTracksDen;bcInROF;eff.", "BC Den Tracks", rofBins, -0.5, rofBins - 0.5);
+    TH1F* hBCTracksNum = new TH1F("hBCTracksNum;bcInROF;eff.", "BC Num Tracks", rofBins, -0.5, rofBins - 0.5);
+    TH1F* hBCTracksFake = new TH1F("hBCTracksFake;bcInROF;eff.", "BC Fake Tracks", rofBins, -0.5, rofBins - 0.5);
+    TH1F* hBCTracksSum = new TH1F("hBCTracksSum;bcInROF;eff.", "BC Sum Tracks", rofBins, -0.5, rofBins - 0.5);
+
+    for (auto& evInfo : info) {
+      for (auto& part : evInfo) {
+        int nCl{0};
+        for (unsigned int bit{0}; bit < sizeof(pInfo.clusters) * 8; ++bit) {
+          nCl += bool(part.clusters & (1 << bit));
+        }
+        if (nCl < 5 || part.bcInROF < 0) {
+          continue;
+        }
+        hBC->Fill(part.bcInROF);
+        hBCTracksDen->Fill(part.bcInROF);
+        if (part.isReco) {
+          hBCTracksNum->Fill(part.bcInROF);
+        }
+        if (part.isFake) {
+          hBCTracksFake->Fill(part.bcInROF);
+        }
+      }
+    }
+
+    hBCTracksSum->Add(hBCTracksNum);
+    hBCTracksSum->Add(hBCTracksFake);
+    hBCTracksSum->Divide(hBCTracksSum, hBCTracksDen, 1., 1., "b");
+    hBCTracksNum->Divide(hBCTracksNum, hBCTracksDen, 1., 1., "b");
+    hBCTracksFake->Divide(hBCTracksFake, hBCTracksDen, 1., 1., "b");
+
+    hBCTracksSum->SetLineColor(kBlack);
+    hBCTracksFake->SetLineColor(2);
+
+    auto c = new TCanvas;
+    c->Divide(2, 1);
+    c->cd(1);
+    hBC->Draw();
+    c->cd(2);
+    gPad->DrawFrame(-0.5, 1e-3, rofBins - 0.5, 1.1, "Tracking >4 ITS cls;bcInROF;eff.");
+    gPad->SetGrid();
+    hBCTracksSum->Draw("histe;same");
+    hBCTracksNum->Draw("histe;same");
+    hBCTracksFake->Draw("histe;same");
+    auto leg = new TLegend;
+    leg->AddEntry(hBCTracksSum, "Sum");
+    leg->AddEntry(hBCTracksNum, "Good");
+    leg->AddEntry(hBCTracksFake, "Fake");
+    leg->Draw();
+  }
 }