Simple macro for ITS tracks vs BC in ROF

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

@@ -0,0 +1,397 @@
+// Copyright 2019-2026 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#if !defined(__CLING__) || defined(__ROOTCLING__)
+#include <TSystem.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TKey.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TProfile.h>
+#include <TCanvas.h>
+#include <CCDB/BasicCCDBManager.h>
+#include <CommonConstants/LHCConstants.h>
+#include "DataFormatsParameters/GRPLHCIFData.h"
+#include "DataFormatsFIT/Triggers.h"
+#include <DataFormatsFT0/Digit.h>
+#include <Framework/DataTypes.h>
+#include <DataFormatsParameters/AggregatedRunInfo.h>
+#endif
+
+int32_t nBCsPerOrbit = o2::constants::lhc::LHCMaxBunches;
+int32_t nBCsPerITSROF = 198;
+int32_t offsetITSROF = 64; // from current analysis of delay of ITS ROF with respect to other detectors
+constexpr float minCollTimeFT0A = -5.f;
+constexpr float maxCollTimeFT0A = 5.f;
+constexpr float minCollTimeFT0C = -5.f;
+constexpr float maxCollTimeFT0C = 5.f;
+constexpr float minCollTimeFV0A = -5.f;
+constexpr float maxCollTimeFV0A = 5.f;
+constexpr float maxZv = 10.f;
+constexpr float maxTrackEta = 0.8f;
+
+// Simple macro to check ITS tracks in AO2D files as a function of bunch crossing within a ROF
+
+void CheckITSTracksVsROF(int maxFiles = 6, int maxDF = 99999)
+{
+
+  printf("nBCsPerOrbit = %d\n", nBCsPerOrbit);
+  // define histograms
+  TH1F* hTrNoColl = new TH1F("hTrNoColl", "", 2, -0.5, 1.5);
+  hTrNoColl->GetXaxis()->SetBinLabel(1, "No collision assoc.");
+  hTrNoColl->GetXaxis()->SetBinLabel(2, "With collision assoc.");
+  TH1F* hTrTime = new TH1F("hTrTime", ";time;entries", 100, -10000., 10000.);
+  TH1F* hCollTime = new TH1F("hCollTime", ";time;entries", 100, -20., 20.);
+  TH1F* hVz = new TH1F("hVz", ";z_{v} (cm);entries", 100, -20., 20.);
+  TH2F* hTrVsCollTime = new TH2F("hTrVsCollTime", ";coll time;track time", 100, -20., 20., 100, -10000., 10000.);
+  int nBinsROF = 395;
+  float minROF = -197.5;
+  float maxROF = 197.5;
+  TH2F* hCluTrackVsBCinROF = new TH2F("hCluTrackVsBC", ";bcInROF;N_{ITSclusters}", nBinsROF, minROF, maxROF, 8, -0.5, 7.5);
+  TH1F* hNEventsVsBC = new TH1F("hNEventsVsBC", ";bc;N_{Collisions}", 3601, -0.5, 3600.5);
+  TH1F* hNTracksVsBC = new TH1F("hNTracksVsBC", ";bc;N_{Tracks}", 3601, -0.5, 3600.5);
+  TH1F* hNTracksPerEventVsBC = new TH1F("hNTracksPerEventVsBC", ";bc;N_{Tracks} per event", 3601, -0.5, 3600.5);
+  TH1F* hNEventsVsBCinROF = new TH1F("hNEventsVsBCinROF", ";bcInROF;N_{Collisions}", nBinsROF, minROF, maxROF);
+  TH1F* hNTracksVsBCinROF = new TH1F("hNTracksVsBCinROF", ";bcInROF;N_{Tracks}", nBinsROF, minROF, maxROF);
+  TH1F* hNTracksPerEventVsBCinROF = new TH1F("hNTracksPerEventVsBCinROF", ";bcInROF;N_{Tracks} per event", nBinsROF, minROF, maxROF);
+  TH2F* hNPVContribVsBCinROF = new TH2F("hNPVContribVsBC", ";bcInROF;N_{PVcontrib}", nBinsROF, minROF, maxROF, 151, -0.5, 150.5);
+
+  // process all AO2D files in subdirectories of the current directory
+  gSystem->Exec("find ./ -name AO2D.root | sort > filelist.txt");
+  int nFiles = 0;
+  std::ifstream aodList("filelist.txt");
+  std::string line;
+  std::vector<std::string> filNames;
+  while (std::getline(aodList, line)) {
+    if (!line.empty()) {
+      filNames.push_back(line);
+      ++nFiles;
+    }
+  }
+  if (nFiles > maxFiles)
+    nFiles = maxFiles;
+
+  for (int jf = 0; jf < nFiles; ++jf) {
+    TFile* fil = TFile::Open(filNames[jf].c_str());
+    if (!fil || fil->IsZombie()) {
+      printf("Cannot open %s\n", filNames[jf].c_str());
+      continue;
+    }
+    TList* l = (TList*)fil->GetListOfKeys();
+    int nKeys = l->GetEntries();
+    for (int j = 0; j < nKeys; ++j) {
+      if (j >= maxDF)
+        break;
+      TKey* k = (TKey*)l->At(j);
+      TString cname = k->GetClassName();
+      if (cname == "TDirectoryFile") {
+        TDirectoryFile* d = (TDirectoryFile*)fil->Get(k->GetName());
+        printf("%s %s\n", fil->GetName(), d->GetName());
+
+        // Access TTree with BC info
+        TTree* tb = (TTree*)d->Get("O2bc_001");
+        if (!tb) {
+          printf("TTree O2bc_001 missing\n");
+          continue;
+        }
+        printf(" BC Tree entries = %lld\n", tb->GetEntries());
+        ULong64_t gloBC;
+        int nRun;
+        tb->SetBranchAddress("fRunNumber", &nRun);
+        tb->SetBranchAddress("fGlobalBC", &gloBC);
+
+        // Access TTree with collision info
+        TTree* tc = (TTree*)d->Get("O2collision_001");
+        if (!tc) {
+          printf("TTree O2collision_001 missing\n");
+          continue;
+        }
+        printf(" Collision Tree entries = %lld\n", tc->GetEntries());
+        float xv, yv, zv, ctime;
+        int iBC;
+        ushort nPVcontrib;
+        tc->SetBranchAddress("fPosX", &xv);
+        tc->SetBranchAddress("fPosY", &yv);
+        tc->SetBranchAddress("fPosZ", &zv);
+        tc->SetBranchAddress("fCollisionTime", &ctime);
+        tc->SetBranchAddress("fNumContrib", &nPVcontrib);
+        tc->SetBranchAddress("fIndexBCs", &iBC);
+
+        // Access TTree with FT0 info
+        TTree* tft = (TTree*)d->Get("O2ft0");
+        if (!tft) {
+          printf("TTree O2ft0 missing\n");
+          continue;
+        }
+        printf(" FT0 Tree entries = %lld\n", tft->GetEntries());
+        float timeft0a, timeft0c;
+        int iBCft0;
+        UChar_t trigMask;
+        tft->SetBranchAddress("fIndexBCs", &iBCft0);
+        tft->SetBranchAddress("fTimeA", &timeft0a);
+        tft->SetBranchAddress("fTimeC", &timeft0c);
+        tft->SetBranchAddress("fTriggerMask", &trigMask);
+
+        // Access TTree with FV0 info
+        TTree* tfv = (TTree*)d->Get("O2fv0a");
+        if (!tfv) {
+          printf("TTree O2fv0a missing\n");
+          continue;
+        }
+        printf(" FV0A Tree entries = %lld\n", tfv->GetEntries());
+        float timefv0a;
+        int iBCfv0;
+        tfv->SetBranchAddress("fIndexBCs", &iBCfv0);
+        tfv->SetBranchAddress("fTime", &timefv0a);
+
+        // Access TTree with track parameters at their innermost update
+        TTree* tt = (TTree*)d->Get("O2track_iu");
+        if (!tt) {
+          printf("TTree O2track_iu missing\n");
+          continue;
+        }
+        printf(" Track Tree entries = %lld\n", tt->GetEntries());
+        // 5 parameters defining the track helix: y, z, sin(phi), tan(lambda), q/pt
+        // x is a coordinate along the track
+        // alpha: angle between track reference system and global reference system
+        float x, alp, y, z, snp, tgl, qpt;
+        int iColl;
+        tt->SetBranchAddress("fIndexCollisions", &iColl);
+        tt->SetBranchAddress("fX", &x);
+        tt->SetBranchAddress("fAlpha", &alp);
+        tt->SetBranchAddress("fY", &y);
+        tt->SetBranchAddress("fZ", &z);
+        tt->SetBranchAddress("fSnp", &snp);
+        tt->SetBranchAddress("fTgl", &tgl);
+        tt->SetBranchAddress("fSigned1Pt", &qpt);
+
+        // Access TTree with track extra information
+        TTree* te = (TTree*)d->Get("O2trackextra_002");
+        if (!te)
+          te = (TTree*)d->Get("O2trackextra_001");
+        if (!te) {
+          printf("TTree O2trackextra_002 and O2trackextra_001 both missing\n");
+          continue;
+        }
+        uint itsclusiz;
+        uint8_t nTPCclu;
+        uint trflag;
+        float itschi2, tpcchi2, trtime;
+        te->SetBranchAddress("fITSClusterSizes", &itsclusiz);
+        te->SetBranchAddress("fITSChi2NCl", &itschi2);
+        te->SetBranchAddress("fTPCNClsFindable", &nTPCclu);
+        te->SetBranchAddress("fTPCChi2NCl", &tpcchi2);
+        te->SetBranchAddress("fTrackTime", &trtime);
+        te->SetBranchAddress("fFlags", &trflag);
+
+        // get run information
+        if (tb->GetEntries() == 0) {
+          printf("Empty BC tree\n");
+          continue;
+        }
+        tb->GetEntry(0);
+        auto runInfo = o2::parameters::AggregatedRunInfo::buildAggregatedRunInfo(o2::ccdb::BasicCCDBManager::instance(), nRun);
+        //        uint64_t sorTimestamp = runInfo.sor;
+        //        uint64_t eorTimestamp = runInfo.eor;
+        int64_t bcSOR = runInfo.orbitSOR * nBCsPerOrbit;
+        int64_t nBCsPerTF = runInfo.orbitsPerTF * nBCsPerOrbit;
+        auto calcBcInROF = [&](ULong64_t bcGlobal) {
+          int bc = (bcGlobal + nBCsPerOrbit) % nBCsPerITSROF;
+          bc -= offsetITSROF;
+          return bc;
+        };
+
+        // fill maps for BC and FT0
+        int nFT0 = tft->GetEntries();
+        std::unordered_map<int, int> bcToFT0;
+        bcToFT0.reserve(nFT0);
+        for (int jft = 0; jft < nFT0; ++jft) {
+          tft->GetEntry(jft);
+          bcToFT0[iBCft0] = jft;
+        }
+        // fill maps for BC and FV0
+        int nFV0 = tfv->GetEntries();
+        std::unordered_map<int, int> bcToFV0;
+        bcToFV0.reserve(nFV0);
+        for (int jfv = 0; jfv < nFV0; ++jfv) {
+          tfv->GetEntry(jfv);
+          bcToFV0[iBCfv0] = jfv;
+        }
+        // tag good collisions
+        int nColl = tc->GetEntries();
+        std::vector<int> ft0indices(nColl, -1);
+        std::vector<int> fv0indices(nColl, -1);
+        std::vector<bool> isGoodColl(nColl, false);
+        for (int i = 0; i < nColl; ++i) {
+          ft0indices[i] = -1;
+          tc->GetEntry(i);
+          if (iBC < 0 || iBC >= tb->GetEntries()) {
+            printf("ERROR: iBC out of range\n");
+            continue;
+          }
+          tb->GetEntry(iBC);
+          auto it = bcToFT0.find(iBC);
+          if (it != bcToFT0.end())
+            ft0indices[i] = it->second;
+          auto iv = bcToFV0.find(iBC);
+          if (iv != bcToFV0.end())
+            fv0indices[i] = iv->second;
+          if (std::abs(zv) > maxZv)
+            continue;
+          int jft0 = ft0indices[i];
+          int jfv0 = fv0indices[i];
+          if (jft0 >= 0 && jfv0 >= 0) {
+            tft->GetEntry(jft0);
+            tfv->GetEntry(jfv0);
+            bool isTVX = TESTBIT(trigMask, o2::ft0::Triggers::bitVertex);
+            int64_t bcInTF = (static_cast<int64_t>(gloBC) - bcSOR) % nBCsPerTF;
+            bool isTFBorderOK = (bcInTF > 300 && bcInTF < nBCsPerTF - 4000);
+            if (isTFBorderOK && isTVX &&
+                timeft0a > minCollTimeFT0A && timeft0a < maxCollTimeFT0A &&
+                timeft0c > minCollTimeFT0C && timeft0c < maxCollTimeFT0C &&
+                timefv0a > minCollTimeFV0A && timefv0a < maxCollTimeFV0A) {
+              isGoodColl[i] = true;
+              int bcInITSROF = calcBcInROF(gloBC);
+              hNEventsVsBC->Fill(gloBC % nBCsPerOrbit);
+              hNEventsVsBCinROF->Fill(bcInITSROF);
+              hCollTime->Fill(ctime);
+              hVz->Fill(zv);
+            }
+          }
+        }
+
+        // loop over tracks
+        int nTracks = tt->GetEntries();
+        for (int i = 0; i < nTracks; ++i) {
+          // O2track_iu and O2trackextra_xxx entries are aligned 1-to-1
+          tt->GetEntry(i);
+          te->GetEntry(i);
+          // compute the ITS cluster map, which can be used in track selecion
+          int nITSclu = 0;
+          uint itsCluMap = 0;
+          for (int jLay = 0; jLay < 7; ++jLay) {
+            if ((itsclusiz >> (jLay * 4)) & 0xf) {
+              nITSclu++;
+              itsCluMap |= (1 << jLay);
+            }
+          }
+          // track selections
+          if (!(trflag & o2::aod::track::PVContributor))
+            continue;
+          if (nITSclu == 0 || itschi2 < 0. || tpcchi2 < 0.)
+            continue;
+          // compute eta from track dip angle lambda via tgl = tan(lambda)
+          float cl = 1. / std::sqrt(1. + tgl * tgl);
+          float sl = cl * tgl;
+          float eta = 0.5 * std::log((1 + sl) / (1 - sl));
+          if (std::abs(eta) > maxTrackEta)
+            continue;
+          // collision selections
+          if (iColl >= 0 && iColl < nColl && isGoodColl[iColl]) {
+            tc->GetEntry(iColl);
+            hTrTime->Fill(trtime);
+            hTrVsCollTime->Fill(ctime, trtime);
+            tb->GetEntry(iBC);
+            int bcInITSROF = calcBcInROF(gloBC);
+            hNTracksVsBC->Fill(gloBC % nBCsPerOrbit);
+            hCluTrackVsBCinROF->Fill(bcInITSROF, nITSclu);
+            hNTracksVsBCinROF->Fill(bcInITSROF);
+            hNPVContribVsBCinROF->Fill(bcInITSROF, nPVcontrib);
+            hTrNoColl->Fill(1);
+          } else {
+            hTrNoColl->Fill(0);
+          }
+        }
+      }
+    }
+  }
+
+  // plotting
+  TProfile* pctr = hNPVContribVsBCinROF->ProfileX();
+  TProfile* pclu = hCluTrackVsBCinROF->ProfileX();
+  TH1F* hFrac7 = new TH1F("hFrac7", ";bcInROF;Fraction of tracks with 7 clusters", hCluTrackVsBCinROF->GetXaxis()->GetNbins(), hCluTrackVsBCinROF->GetXaxis()->GetXmin(), hCluTrackVsBCinROF->GetXaxis()->GetXmax());
+  for (int jx = 1; jx <= hCluTrackVsBCinROF->GetXaxis()->GetNbins(); jx++) {
+    double tot = 0.;
+    for (int jy = 1; jy <= hCluTrackVsBCinROF->GetYaxis()->GetNbins(); jy++) {
+      double c = hCluTrackVsBCinROF->GetBinContent(jx, jy);
+      tot += c;
+    }
+    double seven = hCluTrackVsBCinROF->GetBinContent(jx, 8);
+    double f = 0.;
+    double ef = 0.;
+    if (tot > 0) {
+      f = seven / tot;
+      ef = std::sqrt(f * (1 - f) / tot);
+    }
+    hFrac7->SetBinContent(jx, f);
+    hFrac7->SetBinError(jx, ef);
+  }
+  for (int jx = 1; jx <= hNTracksVsBC->GetXaxis()->GetNbins(); jx++) {
+    double nt = hNTracksVsBC->GetBinContent(jx);
+    double ne = hNEventsVsBC->GetBinContent(jx);
+    if (ne > 0)
+      hNTracksPerEventVsBC->SetBinContent(jx, nt / ne);
+  }
+  for (int jx = 1; jx <= hNTracksVsBCinROF->GetXaxis()->GetNbins(); jx++) {
+    double nt = hNTracksVsBCinROF->GetBinContent(jx);
+    double ne = hNEventsVsBCinROF->GetBinContent(jx);
+    if (ne > 0)
+      hNTracksPerEventVsBCinROF->SetBinContent(jx, nt / ne);
+  }
+
+  hCluTrackVsBCinROF->SetStats(0);
+  hFrac7->SetStats(0);
+  hNPVContribVsBCinROF->SetStats(0);
+
+  TCanvas* ctr = new TCanvas("ctr", "", 1400, 500);
+  ctr->Divide(3, 1);
+  ctr->cd(1);
+  hTrNoColl->Draw();
+  ctr->cd(2);
+  hCollTime->Draw();
+  ctr->cd(3);
+  hVz->Draw();
+
+  TCanvas* cbc = new TCanvas("cbc", "", 1400, 800);
+  cbc->Divide(1, 3);
+  cbc->cd(1);
+  hNEventsVsBC->Draw();
+  cbc->cd(2);
+  hNTracksVsBC->Draw();
+  cbc->cd(3);
+  hNTracksPerEventVsBC->Draw();
+
+  TCanvas* crof = new TCanvas("crof", "", 1400, 800);
+  crof->Divide(3, 2);
+  crof->cd(1);
+  hNEventsVsBCinROF->Draw();
+  crof->cd(2);
+  hNTracksVsBCinROF->Draw();
+  crof->cd(3);
+  hNTracksPerEventVsBCinROF->Draw();
+  crof->cd(4);
+  gPad->SetLogz();
+  hCluTrackVsBCinROF->Draw("colz");
+  pclu->SetLineWidth(2);
+  pclu->Draw("same");
+  crof->cd(5);
+  hFrac7->SetMinimum(0.);
+  hFrac7->SetMaximum(1.);
+  hFrac7->SetLineWidth(2);
+  hFrac7->Draw();
+  crof->cd(6);
+  gPad->SetLogz();
+  hNPVContribVsBCinROF->Draw("colz");
+  pctr->SetLineWidth(2);
+  pctr->Draw("same");
+  crof->SaveAs("ITSTracksVsBCinROF.png");
+}