Merge cbf50d0 into sapling-pr-archive-ktf

Author: ktf

Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TrackMCStudyTypes.h

@@ -86,6 +86,28 @@ struct RecTrack {
   ClassDefNV(RecTrack, 1);
 };
 
+struct TrackPairInfo {
+  RecTrack tr0;
+  RecTrack tr1;
+  uint8_t nshTPC = 0;
+  uint8_t nshTPCRow = 0;
+
+  int getComb() const { return tr0.track.getSign() != tr1.track.getSign() ? 0 : (tr0.track.getSign() > 0 ? 1 : 2); }
+  float getDPhi() const
+  {
+    float dphi = tr0.track.getPhi() - tr1.track.getPhi();
+    if (dphi < -o2::constants::math::PI) {
+      dphi += o2::constants::math::TwoPI;
+    } else if (dphi > o2::constants::math::PI) {
+      dphi -= o2::constants::math::TwoPI;
+    }
+    return dphi;
+  }
+  float getDTgl() const { return tr0.track.getTgl() - tr1.track.getTgl(); }
+
+  ClassDefNV(TrackPairInfo, 1)
+};
+
 struct TrackFamily { // set of tracks related to the same MC label
   MCTrackInfo mcTrackInfo{};
   std::vector<RecTrack> recTracks{};

Detectors/GlobalTrackingWorkflow/study/src/GlobalTrackingStudyLinkDef.h

@@ -36,5 +36,7 @@
 #pragma link C++ class o2::trackstudy::ClResTPC + ;
 #pragma link C++ class o2::trackstudy::ClResTPCCont + ;
 #pragma link C++ class std::vector < o2::trackstudy::ClResTPCCont> + ;
+#pragma link C++ class o2::trackstudy::TrackPairInfo + ;
+#pragma link C++ class std::vector < o2::trackstudy::TrackPairInfo> + ;
 
 #endif

Detectors/GlobalTrackingWorkflow/study/src/TrackingStudy.cxx

@@ -33,6 +33,7 @@
 #include "DetectorsBase/GRPGeomHelper.h"
 #include "GlobalTrackingStudy/TrackingStudy.h"
 #include "GlobalTrackingStudy/TrackInfoExt.h"
+#include "GlobalTrackingStudy/TrackMCStudyTypes.h"
 #include "TPCBase/ParameterElectronics.h"
 #include "ReconstructionDataFormats/PrimaryVertex.h"
 #include "ReconstructionDataFormats/PrimaryVertexExt.h"
@@ -108,6 +109,7 @@ class TrackingStudySpec : public Task
   int mNHBPerTF = 0;
   float mNTPCOccBinLengthInv;
   bool mStoreWithITSOnly = false;
+  bool mDoPairsCorr = false;
   std::string mDCAYFormula = "0.0105 + 0.0350 / pow(x, 1.1)";
   std::string mDCAZFormula = "0.0105 + 0.0350 / pow(x, 1.1)";
   GTrackID::mask_t mTracksSrc{};
@@ -136,6 +138,7 @@ void TrackingStudySpec::init(InitContext& ic)
   mMinTPCClusters = ic.options().get<int>("min-tpc-clusters");
   mDCAYFormula = ic.options().get<std::string>("dcay-vs-pt");
   mDCAZFormula = ic.options().get<std::string>("dcaz-vs-pt");
+  mDoPairsCorr = ic.options().get<bool>("pair-correlations");
 }
 
 void TrackingStudySpec::run(ProcessingContext& pc)
@@ -245,6 +248,7 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
   float tBiasITS = alpParams.roFrameBiasInBC * o2::constants::lhc::LHCBunchSpacingMUS;
   const o2::ft0::InteractionTag& ft0Params = o2::ft0::InteractionTag::Instance();
   std::vector<o2::dataformats::TrackInfoExt> trcExtVec;
+  std::vector<o2::trackstudy::TrackPairInfo> trcPairsVec;
   auto vdrit = mTPCVDriftHelper.getVDriftObject().getVDrift();
   bool tpcTrackOK = recoData.isTrackSourceLoaded(GTrackID::TPC);
 
@@ -278,6 +282,82 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
     }
   };
 
+  auto getTPCPairSharing = [&recoData, this](const o2::tpc::TrackTPC& trc0, const o2::tpc::TrackTPC& trc1) {
+    const auto clRefs = recoData.getTPCTracksClusterRefs();
+    const auto shMap = recoData.clusterShMapTPC.data();
+    uint8_t nsh = 0, nshRows = 0, lastSharedRow = -1;
+    if (recoData.inputsTPCclusters) {
+      uint8_t clSect0 = 0, clRow0 = 0, clSect1 = 0, clRow1 = 0;
+      uint32_t clIdx0 = 0, clIdx1 = 0;
+      int ic1Start = 0;
+      for (int ic0 = 0; ic0 < trc0.getNClusterReferences(); ic0++) { // outside -> inside
+        trc0.getClusterReference(clRefs, ic0, clSect0, clRow0, clIdx0);
+        for (int ic1 = ic1Start; ic1 < trc1.getNClusterReferences(); ic1++) { // outside -> inside
+          trc1.getClusterReference(clRefs, ic1, clSect1, clRow1, clIdx1);
+          if (clRow1 > clRow0) {
+            ic1Start = ic1 + 1;
+            continue; // catch up ic0
+          }
+          if (clRow1 == clRow0) {
+            if (clSect0 == clSect1 && clIdx0 == clIdx1) {
+              nsh++;
+              if (lastSharedRow != clRow0) {
+                lastSharedRow = clRow0;
+                nshRows++;
+              }
+              ic1Start = ic1 + 1;
+              break; // check next ic0
+            }
+          }
+        }
+      }
+    }
+    return std::make_pair(nsh, nshRows);
+  };
+
+  auto assignRecTrack = [&recoData, this](const o2::dataformats::TrackInfoExt& src, o2::trackstudy::RecTrack& dst) {
+    dst.track = src.track;
+    dst.gid = src.gid;
+    dst.ts.setTimeStamp(src.ttime);
+    dst.ts.setTimeStampError(src.ttimeE);
+    dst.nClITS = src.nClITS;
+    dst.nClTPC = src.nClTPC;
+    dst.pattITS = src.pattITS;
+    if (src.q2ptITS == 0. && dst.nClITS > 0) {
+      dst.pattITS |= 0x1 << 7;
+    }
+    dst.lowestPadRow = src.rowMinTPC;
+    if (this->mUseMC) {
+      auto gidSet = recoData.getSingleDetectorRefs(src.gid);
+      if (recoData.getTrackMCLabel(src.gid).isFake()) {
+        dst.flags |= RecTrack::FakeGLO;
+      }
+      auto msk = src.gid.getSourceDetectorsMask();
+      if (msk[DetID::ITS]) {
+        if (gidSet[GTrackID::ITS].isSourceSet()) { // has ITS track rather than AB tracklet
+          auto lblITS = recoData.getTrackMCLabel(gidSet[GTrackID::ITS]);
+          if (lblITS.isFake()) {
+            dst.flags |= RecTrack::FakeITS;
+          }
+        } else { // AB ITS tracklet
+          if (recoData.getTrackMCLabel(gidSet[GTrackID::ITSAB]).isFake()) {
+            dst.flags |= RecTrack::FakeITS;
+          }
+        }
+        if (msk[DetID::TPC]) { // has both ITS and TPC contribution
+          if (recoData.getTrackMCLabel(gidSet[GTrackID::ITSTPC]).isFake()) {
+            dst.flags |= RecTrack::FakeITSTPC;
+          }
+        }
+      }
+      if (msk[DetID::TPC]) {
+        if (recoData.getTrackMCLabel(gidSet[GTrackID::TPC]).isFake()) {
+          dst.flags |= RecTrack::FakeTPC;
+        }
+      }
+    }
+  };
+
   for (int iv = 0; iv < nv; iv++) {
     LOGP(debug, "processing PV {} of {}", iv, nv);
     const auto& vtref = vtxRefs[iv];
@@ -309,6 +389,7 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
       pve.VtxID = iv;
     }
     trcExtVec.clear();
+    trcPairsVec.clear();
     float q2ptITS, q2ptTPC, q2ptITSTPC, q2ptITSTPCTRD;
     for (int is = 0; is < GTrackID::NSources; is++) {
       DetID::mask_t dm = GTrackID::getSourceDetectorsMask(is);
@@ -444,6 +525,42 @@ void TrackingStudySpec::process(o2::globaltracking::RecoContainer& recoData)
                << "orbit=" << recoData.startIR.orbit << "tfID=" << TFCount
                << "tpcOccBef=" << tpcOccBef << "tpcOccAft=" << tpcOccAft
                << "pve=" << pveVec[iv] << "trc=" << trcExtVec << "\n";
+
+    if (mDoPairsCorr) {
+      for (int it0 = 0; it0 < (int)trcExtVec.size(); it0++) {
+        const auto& tr0 = trcExtVec[it0];
+        if (tr0.nClTPC < 1) {
+          continue;
+        }
+        for (int it1 = it0 + 1; it1 < (int)trcExtVec.size(); it1++) {
+          const auto& tr1 = trcExtVec[it1];
+          if (tr1.nClTPC < 1) {
+            continue;
+          }
+
+          if (std::abs(tr0.track.getTgl() - tr1.track.getTgl()) > 0.25) {
+            continue;
+          }
+          auto dphi = tr0.track.getPhi() - tr1.track.getPhi();
+          if (dphi < -o2::constants::math::PI) {
+            dphi += o2::constants::math::TwoPI;
+          } else if (dphi > o2::constants::math::PI) {
+            dphi -= o2::constants::math::TwoPI;
+          }
+          if (std::abs(dphi) > 0.25) {
+            continue;
+          }
+          auto& pr = trcPairsVec.emplace_back();
+          assignRecTrack(tr0, pr.tr0);
+          assignRecTrack(tr1, pr.tr1);
+          auto shinfo = getTPCPairSharing(recoData.getTPCTrack(recoData.getTPCContributorGID(tr0.gid)), recoData.getTPCTrack(recoData.getTPCContributorGID(tr1.gid)));
+          pr.nshTPC = shinfo.first;
+          pr.nshTPCRow = shinfo.second;
+        }
+      }
+    }
+    (*mDBGOut) << "pairs"
+               << "pr=" << trcPairsVec << "\n";
   }
 
   int nvtot = mMaxNeighbours < 0 ? -1 : (int)pveVec.size();
@@ -600,6 +717,7 @@ DataProcessorSpec getTrackingStudySpec(GTrackID::mask_t srcTracks, GTrackID::mas
     {"max-eta", VariantType::Float, 1.0f, {"Cut on track eta"}},
     {"min-pt", VariantType::Float, 0.1f, {"Cut on track pT"}},
     {"with-its-only", VariantType::Bool, false, {"Store tracks with ITS only"}},
+    {"pair-correlations", VariantType::Bool, false, {"Do pairs correlation"}},
     {"min-x-prop", VariantType::Float, 100.f, {"track should be propagated to this X at least"}},
   };
   o2::tpc::VDriftHelper::requestCCDBInputs(dataRequest->inputs);

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

@@ -87,6 +87,7 @@ void CheckTracksCA(bool doFakeClStud = false,
   TTree* mcTree = (TTree*)gFile->Get("o2sim");
   mcTree->SetBranchStatus("*", 0); // disable all branches
   mcTree->SetBranchStatus("MCTrack*", 1);
+  mcTree->SetBranchStatus("MCEventHeader*", 1);
 
   std::vector<o2::MCTrack>* mcArr = nullptr;
   mcTree->SetBranchAddress("MCTrack", &mcArr);
@@ -115,10 +116,13 @@ void CheckTracksCA(bool doFakeClStud = false,
   std::cout << "** Filling particle table ... " << std::flush;
   int lastEventIDcl = -1, cf = 0;
   int nev = mcTree->GetEntriesFast();
-  std::vector<std::vector<ParticleInfo>> info(nev);
+  std::vector<std::vector<ParticleInfo>> info;
+  info.resize(nev);
+  TH1D* hZvertex = new TH1D("hZvertex", "Z vertex", 100, -20, 20);
   for (int n = 0; n < nev; n++) { // loop over MC events
     mcTree->GetEvent(n);
     info[n].resize(mcArr->size());
+    hZvertex->Fill(mcEvent->GetZ());
     for (unsigned int mcI{0}; mcI < mcArr->size(); ++mcI) {
       auto part = mcArr->at(mcI);
       info[n][mcI].event = n;
@@ -196,7 +200,6 @@ void CheckTracksCA(bool doFakeClStud = false,
         info[evID][trackID].track.getImpactParams(info[evID][trackID].pvx, info[evID][trackID].pvy, info[evID][trackID].pvz, bz, ip);
         info[evID][trackID].dcaxy = ip[0];
         info[evID][trackID].dcaz = ip[1];
-        Info("", "dcaxy=%f dcaz=%f bz=%f", ip[0], ip[1], bz);
       }
 
       fakes += fake;
@@ -286,6 +289,10 @@ void CheckTracksCA(bool doFakeClStud = false,
   clone->Divide(clone, den, 1, 1, "b");
   clone->SetLineColor(3);
   clone->Draw("histesame");
+  TCanvas* c2 = new TCanvas;
+  c2->SetGridx();
+  c2->SetGridy();
+  hZvertex->DrawClone();
 
   std::cout << "** Streaming output TTree to file ... " << std::flush;
   TFile file("CheckTracksCA.root", "recreate");

Detectors/ITSMFT/ITS/tracking/include/ITStracking/Configuration.h

@@ -94,12 +94,14 @@ struct TrackingParameters {
   unsigned long MaxMemory = 12000000000UL;
   float MaxChi2ClusterAttachment = 60.f;
   float MaxChi2NDF = 30.f;
-  float MinPt = 0.f;
+  std::vector<float> MinPt = {0.f, 0.f, 0.f, 0.f};
   unsigned char StartLayerMask = 0x7F;
   bool FindShortTracks = false;
   bool PerPrimaryVertexProcessing = false;
   bool SaveTimeBenchmarks = false;
   bool DoUPCIteration = false;
+  bool FataliseUponFailure = true;
+  bool DropTFUponFailure = false;
   /// Cluster attachment
   bool UseTrackFollower = false;
   bool UseTrackFollowerTop = false;

Detectors/ITSMFT/ITS/tracking/include/ITStracking/TimeFrame.h

@@ -259,7 +259,7 @@ class TimeFrame
   void printCellLUTonLayer(int i);
   void printTrackletLUTs();
   void printCellLUTs();
-  void printROFInfo(const int rofId);
+  void printSliceInfo(const int, const int);
 
   IndexTableUtils mIndexTableUtils;
 
@@ -297,6 +297,13 @@ class TimeFrame
   std::vector<uint8_t> mMultiplicityCutMask;
 
   const o2::base::PropagatorImpl<float>* mPropagatorDevice = nullptr; // Needed only for GPU
+  void dropTracks()
+  {
+    for (auto& v : mTracks) {
+      deepVectorClear(v);
+    }
+  }
+
  protected:
   template <typename T>
   void deepVectorClear(std::vector<T>& vec)

Detectors/ITSMFT/ITS/tracking/include/ITStracking/TrackingConfigParam.h

@@ -22,13 +22,12 @@ namespace its
 
 struct VertexerParamConfig : public o2::conf::ConfigurableParamHelper<VertexerParamConfig> {
 
-  int nIterations = 1;         // Number of vertexing passes to perform
-  int vertPerRofThreshold = 0; // Maximum number of vertices per ROF to trigger second a round
-  bool allowSingleContribClusters = false;
-  // Number of ROFs to be considered for the vertexing
-  int deltaRof = 0;
+  int nIterations = 1;                     // Number of vertexing passes to perform.
+  int vertPerRofThreshold = 0;             // Maximum number of vertices per ROF to trigger second a iteration.
+  bool allowSingleContribClusters = false; // attempt to find vertices in case of a single tracklet found.
+  int deltaRof = 0;                        // Number of ROFs to be considered for the vertexing.
 
-  // geometrical cuts
+  // geometrical cuts for tracklet selection
   float zCut = 0.002f;
   float phiCut = 0.005f;
   float pairCut = 0.04f;
@@ -42,12 +41,12 @@ struct VertexerParamConfig : public o2::conf::ConfigurableParamHelper<VertexerPa
   float maxZPositionAllowed = 25.f; // 4x sZ of the beam
 
   // Artefacts selections
-  int clusterContributorsCut = 16;
+  int clusterContributorsCut = 16; // minimum number of contributors for the second vertex found in the same ROF (pileup cut)
   int maxTrackletsPerCluster = 1e2;
   int phiSpan = -1;
   int zSpan = -1;
-  int ZBins = 1;
-  int PhiBins = 128;
+  int ZBins = 1;     // z-phi index table configutation: number of z bins
+  int PhiBins = 128; // z-phi index table configutation: number of phi bins
 
   int nThreads = 1;
 
@@ -56,11 +55,11 @@ struct VertexerParamConfig : public o2::conf::ConfigurableParamHelper<VertexerPa
 
 struct TrackerParamConfig : public o2::conf::ConfigurableParamHelper<TrackerParamConfig> {
   // Use TGeo for mat. budget
-  bool useMatCorrTGeo = false;
-  bool useFastMaterial = false;
-  int deltaRof = 0;
-  float sysErrY2[7] = {0}; // systematic error^2 in Y per layer
-  float sysErrZ2[7] = {0}; // systematic error^2 in Z per layer
+  bool useMatCorrTGeo = false;  // use full geometry to corect for material budget accounting in the fits. Default is to use the material budget LUT.
+  bool useFastMaterial = false; // use faster material approximation for material budget accounting in the fits.
+  int deltaRof = 0;             // configure the width of the window in ROFs to be considered for the tracking.
+  float sysErrY2[7] = {0};      // systematic error^2 in Y per layer
+  float sysErrZ2[7] = {0};      // systematic error^2 in Z per layer
   float maxChi2ClusterAttachment = -1.f;
   float maxChi2NDF = -1.f;
   float nSigmaCut = -1.f;
@@ -69,23 +68,25 @@ struct TrackerParamConfig : public o2::conf::ConfigurableParamHelper<TrackerPara
   float pvRes = -1.f;
   int LUTbinsPhi = -1;
   int LUTbinsZ = -1;
-  float diamondPos[3] = {0.f, 0.f, 0.f};
-  bool useDiamond = false;
-  unsigned long maxMemory = 0;
-  int useTrackFollower = -1;          // bit 0: allow mixing implies bits 1&2; bit 1: topwards; bit2: downwards; => 0 off
-  float trackFollowerNSigmaZ = 1.f;   // sigma in z-cut for track-following search rectangle
-  float trackFollowerNSigmaPhi = 1.f; // sigma in phi-cut for track-following search rectangle
+  float diamondPos[3] = {0.f, 0.f, 0.f}; // override the position of the vertex
+  bool useDiamond = false;               // enable overriding the vertex position
+  unsigned long maxMemory = 0;           // override default protections on the maximum memory to be used by the tracking
+  int useTrackFollower = -1;             // bit 0: allow mixing implies bits 1&2; bit 1: topwards; bit2: downwards; => 0 off
+  float trackFollowerNSigmaZ = 1.f;      // sigma in z-cut for track-following search rectangle
+  float trackFollowerNSigmaPhi = 1.f;    // sigma in phi-cut for track-following search rectangle
   float cellsPerClusterLimit = -1.f;
   float trackletsPerClusterLimit = -1.f;
   int findShortTracks = -1;
-  int nThreads = 1;
-  int nOrbitsPerIterations = 0;
-  int nROFsPerIterations = 0;
-  bool perPrimaryVertexProcessing = false;
-  bool saveTimeBenchmarks = false;
-  bool overrideBeamEstimation = false; // used by gpuwf only
-  int trackingMode = -1;               // -1: unset, 0=sync, 1=async, 2=cosmics used by gpuwf only
-  bool doUPCIteration = false;
+  int nThreads = 1;                        // number of threads to perform the operations in parallel.
+  int nROFsPerIterations = 0;              // size of the slice of ROFs to be processed at a time, preferably integer divisors of nROFs per TF, to balance the iterations.
+  int nOrbitsPerIterations = 0;            // not implemented: size of the slice of ROFs to be processed at a time, computed using the number of ROFs per orbit.
+  bool perPrimaryVertexProcessing = false; // perform the full tracking considering the vertex hypotheses one at the time.
+  bool saveTimeBenchmarks = false;         // dump metrics on file
+  bool overrideBeamEstimation = false;     // use beam position from meanVertex CCDB object
+  int trackingMode = -1;                   // -1: unset, 0=sync, 1=async, 2=cosmics used by gpuwf only
+  bool doUPCIteration = false;             // Perform an additional iteration for UPC events on tagged vertices. You want to combine this config with VertexerParamConfig.nIterations=2
+  bool fataliseUponFailure = true;         // granular management of the fatalisation in async mode
+  bool dropTFUponFailure = false;
 
   O2ParamDef(TrackerParamConfig, "ITSCATrackerParam");
 };