diff --git a/Detectors/ITSMFT/ITS/macros/test/CheckITSTracksVsROF.C b/Detectors/ITSMFT/ITS/macros/test/CheckITSTracksVsROF.C new file mode 100644 index 0000000000000..5128b4a373757 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "DataFormatsParameters/GRPLHCIFData.h" +#include "DataFormatsFIT/Triggers.h" +#include +#include +#include +#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 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 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 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 ft0indices(nColl, -1); + std::vector fv0indices(nColl, -1); + std::vector 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(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"); +}