Update track selection for track interpolation

martenole · martenole · commit 0b5d4db5e6dc · 2023-01-18T16:18:35.000+01:00
diff --git a/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx b/Detectors/GlobalTrackingWorkflow/tpcinterpolationworkflow/src/TPCInterpolationSpec.cxx
@@ -155,7 +155,13 @@ void TPCInterpolationDPL::run(ProcessingContext& pc)
         gids.push_back(_origID);
         trkCounters[_origID.getSource()] += 1;
       }
-      return true;
+      if (gidTable[GTrackID::TRD].isIndexSet() && gidTable[GTrackID::TOF].isIndexSet()) {
+        // for ITS-TPC-TRD-TOF tracks we are interested also in the ITS-TPC-TRD part
+        // we want to have both available, so return false for the full barrell tracks
+        return false;
+      } else {
+        return true;
+      }
     } else {
       return false;
     }
diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/SpacePointsCalibConfParam.h
@@ -35,14 +35,21 @@ struct SpacePointsCalibConfParam : public o2::conf::ConfigurableParamHelper<Spac
   int minITSNClsNoOuterPoint = 6;  ///< min number of ITS clusters if no hit in TRD or TOF exists
   int minTRDNTrklts = 3;           ///< min number of TRD space points
   float maxITSChi2 = 20.f;         ///< cut on ITS reduced chi2
+  float maxTRDChi2 = 10.f;         ///< cut on TRD reduced chi2
 
   // other settings for track interpolation
   float sigYZ2TOF{.75f}; ///< for now assume cluster error for TOF equal for all clusters in both Y and Z
   float maxSnp{.85f};    ///< max snp when propagating tracks
   float maxStep{2.f};    ///< maximum step for propagation
 
   // steering of map creation after the residuals have already been written to file
+  bool writeBinnedResiduals{false}; ///< when creating the map from unbinned residuals store the binned residuals together with the voxel results
   bool useTrackData{false}; ///< if we have the track data available, we can redefine the above cuts for the map creation, e.g. minTPCNCls etc
+  bool timeFilter{false};   ///< consider only residuals as input from TFs with a specific time range specified via startTimeMS and endTimeMS
+  long startTimeMS{0L};     ///< the start of the time range in MS
+  long endTimeMS{1999999999999L}; ///< the end of the time range in MS
+  bool cutOnDCA{true};            ///< when creating the map from unbinned residuals cut on DCA estimated from ITS outer parameter
+  float maxDCA = 10.f;            ///< DCA cut value in cm
 
   // parameters for outlier rejection
   bool writeUnfiltered{false};           ///< if set, all residuals and track parameters will be aggregated and dumped additionally without outlier rejection
diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackInterpolation.h
@@ -100,7 +100,8 @@ struct TrackDataCompact {
 /// Structure filled for each track with track quality information and a vector with TPCClusterResiduals
 struct TrackData {
   o2::dataformats::GlobalTrackID gid{}; ///< global track ID for seeding track
-  // the track parameters are taken from the ITS track
+  o2::track::TrackPar par{};            ///< ITS track at inner TPC radius
+  // FIXME remove x, alpha and p once we have a new production where all that information is stored in par
   float x{};                                  ///< track X position
   float alpha{};                              ///< track alpha angle
   std::array<float, o2::track::kNParams> p{}; ///< track parameters
@@ -112,7 +113,7 @@ struct TrackData {
   unsigned short nTrkltsTRD{}; ///< number of attached TRD tracklets
   unsigned short clAvailTOF{}; ///< whether or not track seed has a matched TOF cluster
   o2::dataformats::RangeReference<> clIdx{}; ///< index of first cluster residual and total number of cluster residuals of this track
-  ClassDefNV(TrackData, 3);
+  ClassDefNV(TrackData, 4);
 };
 
 /// \class TrackInterpolation
diff --git a/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackResiduals.h b/Detectors/TPC/calibration/SpacePoints/include/SpacePoints/TrackResiduals.h
@@ -420,6 +420,9 @@ class TrackResiduals
   /// Closes the file with the debug output.
   void closeOutputFile();
 
+  /// Allow to access the output file from outside
+  TFile* getOutputFilePtr() { return mFileOut.get(); }
+
   /// Set the voxel statistics directly from outside
   void setStats(const std::vector<TrackResiduals::VoxStats>& statsIn, int iSec);
 
diff --git a/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx b/Detectors/TPC/calibration/SpacePoints/src/TrackInterpolation.cxx
@@ -82,6 +82,17 @@ void TrackInterpolation::process(const o2::globaltracking::RecoContainer& inp, c
   mTrackData.reserve(nSeeds);
   mClRes.reserve(nSeeds * param::NPadRows);
 
+  // to obtain ITS-TPC-TRD-TOF track from ITS-TPC-TRD track we fill the trkMap
+  std::unordered_map<GTrackID, int> trkMap;
+  for (int iTrk = 0; iTrk < nSeeds; ++iTrk) {
+    const auto& gidTable = (*mGIDtables)[iTrk];
+    if (!(gidTable[GTrackID::TRD].isIndexSet() && gidTable[GTrackID::TOF].isIndexSet())) {
+      // no more ITS-TPC-TRD-TOF tracks
+      break;
+    }
+    trkMap.emplace(std::make_pair(gidTable[GTrackID::ITSTPCTRD], iTrk));
+  }
+
   // In case we have more input tracks available than are required per TF
   // we want to sample them. But we still prefer global ITS-TPC-TRD-TOF tracks
   // over ITS-TPC-TRD tracks and so on. So we have to shuffle the indices
@@ -92,6 +103,7 @@ void TrackInterpolation::process(const o2::globaltracking::RecoContainer& inp, c
   std::mt19937 g(rd());
   std::vector<int> trackIndices(nSeeds);
   std::iota(trackIndices.begin(), trackIndices.end(), 0);
+
   std::shuffle(trackIndices.begin(), trackIndices.begin() + trkCounters.at(GTrackID::Source::ITSTPCTRDTOF), g);
   int nTracks = trkCounters.at(GTrackID::Source::ITSTPCTRDTOF);
   std::shuffle(trackIndices.begin() + nTracks, trackIndices.begin() + nTracks + trkCounters.at(GTrackID::Source::ITSTPCTRD), g);
@@ -100,17 +112,30 @@ void TrackInterpolation::process(const o2::globaltracking::RecoContainer& inp, c
   nTracks += trkCounters.at(GTrackID::Source::ITSTPCTOF);
   std::shuffle(trackIndices.begin() + nTracks, trackIndices.begin() + nTracks + trkCounters.at(GTrackID::Source::ITSTPC), g);
 
-  for (int iSeed = 0; iSeed < nSeeds; ++iSeed) {
+  std::vector<int> globalTracksToCheck;
+  for (int iSeed = trkCounters.at(GTrackID::Source::ITSTPCTRDTOF); iSeed < nSeeds; ++iSeed) {
     if (mMaxTracksPerTF >= 0 && mTrackDataCompact.size() >= mMaxTracksPerTF) {
       LOG(info) << "Maximum number of tracks per TF reached. Skipping the remaining " << nSeeds - iSeed << " tracks.";
       break;
     }
+    if (auto search = trkMap.find((*mGIDs)[iSeed]); search != trkMap.end()) {
+      globalTracksToCheck.push_back(search->second);
+    }
     if (gids[trackIndices[iSeed]].includesDet(DetID::TRD) || gids[trackIndices[iSeed]].includesDet(DetID::TOF)) {
       interpolateTrack(trackIndices[iSeed]);
     } else {
       extrapolateTrack(trackIndices[iSeed]);
     }
   }
+  // irrespective of the number of tracks already processed, interpolate the ITS-TPC-TRD-TOF tracks
+  // which belong to the ITS-TPC-TRD tracks that were already processed, to allow their analysis
+  // offline
+  // Q: Or should we for these tracks just skip the ITS-TPC-TRD part and do only ITS-TPC-TRD-TOF?
+  // might still be interesting to compare the two?
+  LOGP(info, "Processing {} ITS-TPC-TRD-TOF tracks for which the ITS-TPC-TRD track was already done", globalTracksToCheck.size());
+  for (auto iTrk : globalTracksToCheck) {
+    interpolateTrack(iTrk);
+  }
 
   LOG(info) << "Could process " << mTrackData.size() << " tracks successfully";
 }
@@ -284,11 +309,7 @@ void TrackInterpolation::interpolateTrack(int iSeed)
   }
 
   trackData.gid = (*mGIDs)[iSeed];
-  trackData.x = (*mSeeds)[iSeed].getX();
-  trackData.alpha = (*mSeeds)[iSeed].getAlpha();
-  for (int i = 0; i < o2::track::kNParams; ++i) {
-    trackData.p[i] = (*mSeeds)[iSeed].getParam(i);
-  }
+  trackData.par = (*mSeeds)[iSeed];
   trackData.chi2TRD = gidTable[GTrackID::TRD].isIndexSet() ? mRecoCont->getITSTPCTRDTrack<o2::trd::TrackTRD>(gidTable[GTrackID::ITSTPCTRD]).getChi2() : 0;
   trackData.chi2TPC = trkTPC.getChi2();
   trackData.chi2ITS = trkITS.getChi2();
@@ -365,11 +386,7 @@ void TrackInterpolation::extrapolateTrack(int iSeed)
     ++nMeasurements;
   }
   trackData.gid = (*mGIDs)[iSeed];
-  trackData.x = (*mSeeds)[iSeed].getX();
-  trackData.alpha = (*mSeeds)[iSeed].getAlpha();
-  for (int i = 0; i < o2::track::kNParams; ++i) {
-    trackData.p[i] = (*mSeeds)[iSeed].getParam(i);
-  }
+  trackData.par = (*mSeeds)[iSeed];
   trackData.chi2TPC = trkTPC.getChi2();
   trackData.chi2ITS = trkITS.getChi2();
   trackData.nClsTPC = trkTPC.getNClusterReferences();
@@ -437,7 +454,7 @@ bool TrackInterpolation::compareToHelix(const TrackData& trk, TrackParams& param
   std::array<float, param::NPadRows> yLab;
   std::array<float, param::NPadRows> sPath;
 
-  float curvature = fabsf(trk.p[o2::track::ParLabels::kQ2Pt] * mBz * o2::constants::physics::LightSpeedCm2S * 1e-14f);
+  float curvature = fabsf(trk.par.getQ2Pt() * mBz * o2::constants::physics::LightSpeedCm2S * 1e-14f);
   int secFirst = clsRes[0].sec;
   float phiSect = (secFirst + .5f) * o2::constants::math::SectorSpanRad;
   float snPhi = sin(phiSect);
diff --git a/Detectors/TPC/calibration/SpacePoints/src/TrackResiduals.cxx b/Detectors/TPC/calibration/SpacePoints/src/TrackResiduals.cxx
@@ -334,6 +334,7 @@ void TrackResiduals::fillStats(int iSec)
 //______________________________________________________________________________
 void TrackResiduals::processSectorResiduals(int iSec)
 {
+  LOGP(info, "Processing {} voxel residuals for sector {}", mLocalResidualsIn.size(), iSec);
   initResultsContainer(iSec);
   std::vector<unsigned short> binData;
   for (const auto& res : mLocalResidualsIn) {
@@ -451,7 +452,7 @@ void TrackResiduals::processVoxelResiduals(std::vector<float>& dy, std::vector<f
     LOG(info) << "voxel " << getGlbVoxBin(resVox.bvox) << " is skipped due to too few entries (" << nPoints << " < " << mParams->minEntriesPerVoxel << ")";
     return;
   } else {
-    LOGF(info, "Processing voxel %i with %i entries", getGlbVoxBin(resVox.bvox), nPoints);
+    LOGF(debug, "Processing voxel %i with %i entries", getGlbVoxBin(resVox.bvox), nPoints);
   }
   std::array<float, 7> zResults;
   resVox.flags = 0;
