fix for TPC-TOF matching (thanks to Francesco)

Author: noferini

DataFormats/Reconstruction/include/ReconstructionDataFormats/MatchInfoTOF.h

@@ -28,7 +28,7 @@ class MatchInfoTOF
   using evIdx = o2::dataformats::EvIndex<int, int>;
 
  public:
-  MatchInfoTOF(evIdx evIdxTOFCl, float chi2, o2::track::TrackLTIntegral trkIntLT, evGIdx evIdxTrack) : mEvIdxTOFCl(evIdxTOFCl), mChi2(chi2), mIntLT(trkIntLT), mEvIdxTrack(evIdxTrack){};
+  MatchInfoTOF(evIdx evIdxTOFCl, float chi2, o2::track::TrackLTIntegral trkIntLT, evGIdx evIdxTrack, float dt = 0, float z = 0) : mEvIdxTOFCl(evIdxTOFCl), mChi2(chi2), mIntLT(trkIntLT), mEvIdxTrack(evIdxTrack), mDeltaT(dt), mZatTOF(z){};
   MatchInfoTOF() = default;
   void setEvIdxTOFCl(evIdx index) { mEvIdxTOFCl = index; }
   void setEvIdxTrack(evGIdx index) { mEvIdxTrack = index; }
@@ -46,13 +46,20 @@ class MatchInfoTOF
   const o2::track::TrackLTIntegral& getLTIntegralOut() const { return mIntLT; }
   void print() const;
 
+  void setDeltaT(float val) { mDeltaT = val; }
+  float getDeltaT() const { return mDeltaT; }
+  void setZatTOF(float val) { mZatTOF = val; }
+  float getZatTOF() const { return mZatTOF; }
+
  private:
   float mChi2;                       // chi2 of the pair track-TOFcluster
   o2::track::TrackLTIntegral mIntLT; ///< L,TOF integral calculated during the propagation
   evIdx mEvIdxTOFCl;                 ///< EvIdx for TOF cluster (first: ev index; second: cluster index)
   evGIdx mEvIdxTrack;                ///< EvIdx for track (first: ev index; second: track global index)
+  float mZatTOF = 0.0;               ///< Z position at  TOF
+  float mDeltaT = 0.0;               ///< tTOF - TPC (microsec)
 
-  ClassDefNV(MatchInfoTOF, 1);
+  ClassDefNV(MatchInfoTOF, 2);
 };
 } // namespace dataformats
 } // namespace o2

Detectors/GlobalTracking/include/GlobalTracking/MatchTOF.h

@@ -90,6 +90,13 @@ class MatchTOF
   void init();
   void initTPConly();
 
+  void setCosmics()
+  {
+    mIsCosmics = true;
+    mSpaceTolerance = 150;
+    mTimeTolerance = 50e3;
+  }
+
   ///< attach DPL data and run
   void run(const gsl::span<const o2::dataformats::TrackTPCITS>& trackArray, const gsl::span<const Cluster>& clusterArray, const o2::dataformats::MCTruthContainer<o2::MCCompLabel>& toflab, const gsl::span<const o2::MCCompLabel>& itstpclab);
   void run(const gsl::span<const o2::tpc::TrackTPC>& trackArray, const gsl::span<const Cluster>& clusterArray, const o2::dataformats::MCTruthContainer<o2::MCCompLabel>& toflab, const gsl::span<const o2::MCCompLabel>& tpclab);
@@ -273,8 +280,9 @@ class MatchTOF
   float mTPCTBinMUSInv = 0.; ///< inverse TPC time bin duration in microseconds
   float mTPCBin2Z = 0.;      ///< conversion coeff from TPC time-bin to Z
 
-  float mTimeTolerance = 1e3; ///<tolerance in ns for track-TOF time bracket matching
-  float mSpaceTolerance = 10; ///<tolerance in cm for track-TOF time bracket matching
+  bool mIsCosmics = false;    ///< switch on to reconstruct cosmics and match with TPC
+  float mTimeTolerance = 1e3; ///< tolerance in ns for track-TOF time bracket matching
+  float mSpaceTolerance = 10; ///< tolerance in cm for track-TOF time bracket matching
   int mSigmaTimeCut = 30.;    ///< number of sigmas to cut on time when matching the track to the TOF cluster
 
   TTree* mInputTreeTracks = nullptr; ///< input tree for tracks

Detectors/GlobalTracking/src/MatchTOF.cxx

@@ -505,13 +505,12 @@ bool MatchTOF::prepareTracks()
   std::array<float, 3> globalPos;
   int itmp = 0;
   for (int sec = o2::constants::math::NSectors; sec--;) {
-    Printf("sector %d", sec);
     auto& cacheTrk = mTracksSectIndexCache[sec];   // array of cached tracks indices for this sector; reminder: they are ordered in time!
     for (int itrk = 0; itrk < cacheTrk.size(); itrk++){
       itmp++; 
       auto& trc = mTracksWork[cacheTrk[itrk]];
       trc.getXYZGlo(globalPos);
-      printf("Track %d: Global coordinates After propagating to 371 cm: globalPos[0] = %f, globalPos[1] = %f, globalPos[2] = %f\n", itrk, globalPos[0], globalPos[1], globalPos[2]);
+      //printf("Track %d: Global coordinates After propagating to 371 cm: globalPos[0] = %f, globalPos[1] = %f, globalPos[2] = %f\n", itrk, globalPos[0], globalPos[1], globalPos[2]);
       //      Printf("The phi angle is %f", TMath::ATan2(globalPos[1], globalPos[0]));
     }
   }
@@ -563,10 +562,14 @@ bool MatchTOF::prepareTPCTracks()
     // create working copy of track param
     timeEst timeInfo;
     // set
+    float extraErr = 0;
+    if (mIsCosmics) {
+      extraErr = 100;
+    }
     timeInfo.setTimeStamp(trcOrig.getTime0() * mTPCTBinMUS);
-    timeInfo.setTimeStampError((trcOrig.getDeltaTBwd() + 5) * mTPCTBinMUS);
+    timeInfo.setTimeStampError((trcOrig.getDeltaTBwd() + 5) * mTPCTBinMUS + extraErr);
     mSideTPC.push_back(trcOrig.hasASideClustersOnly() ? 1 : (trcOrig.hasCSideClustersOnly() ? -1 : 0));
-    mExtraTPCFwdTime.push_back((trcOrig.getDeltaTFwd() + 5) * mTPCTBinMUS);
+    mExtraTPCFwdTime.push_back((trcOrig.getDeltaTFwd() + 5) * mTPCTBinMUS + extraErr);
 
     o2::track::TrackLTIntegral intLT0; //mTPCTracksWork.back().getLTIntegralOut(); // we get the integrated length from TPC-ITC outward propagation
     // make a copy of the TPC track that we have to propagate
@@ -641,6 +644,7 @@ bool MatchTOF::prepareTPCTracks()
       auto& trcB = mTracksWork[b].second;
       return ((trcA.getTimeStamp() - trcA.getTimeStampError()) - (trcB.getTimeStamp() - trcB.getTimeStampError()) < 0.);
     });
+
   } // loop over tracks of single sector
 
   // Uncomment for local debug
@@ -653,14 +657,15 @@ bool MatchTOF::prepareTPCTracks()
     auto& cacheTrk = mTracksSectIndexCache[sec];   // array of cached tracks indices for this sector; reminder: they are ordered in time!
     for (int itrk = 0; itrk < cacheTrk.size(); itrk++){
       itmp++; 
-      auto& trc = mTracksWork[cacheTrk[itrk]];
+      auto& trc = mTracksWork[cacheTrk[itrk]].first;
+      auto& trcAttr = mTracksWork[cacheTrk[itrk]].second;
       trc.getXYZGlo(globalPos);
-      printf("Track %d: Global coordinates After propagating to 371 cm: globalPos[0] = %f, globalPos[1] = %f, globalPos[2] = %f\n", itrk, globalPos[0], globalPos[1], globalPos[2]);
+      printf("Track %d: Global coordinates After propagating to 371 cm: globalPos[0] = %f, globalPos[1] = %f, globalPos[2] = %f -- timestamp = %f +/- %f\n", itrk, globalPos[0], globalPos[1], globalPos[2],trcAttr.getTimeStamp(),trcAttr.getTimeStampError()); 
       //      Printf("The phi angle is %f", TMath::ATan2(globalPos[1], globalPos[0]));
     }
   }
   Printf("we have %d tracks",itmp);      
-  */
+*/
 
   return true;
 }
@@ -1116,6 +1121,7 @@ void MatchTOF::doMatchingForTPC(int sec)
   float vdriftInBC = Geo::BC_TIME_INPS * 1E-6 * vdrift;
 
   int bc_grouping = 40;
+  int bc_grouping_tolerance = bc_grouping + mTimeTolerance / 25;
   int bc_grouping_half = (bc_grouping + 1) / 2;
   double BCgranularity = Geo::BC_TIME_INPS * bc_grouping;
 
@@ -1162,17 +1168,20 @@ void MatchTOF::doMatchingForTPC(int sec)
     minTrkTime = int(minTrkTime / BCgranularity) * BCgranularity;                                        // align min to a BC
     double maxTrkTime = (trackWork.second.getTimeStamp() + mExtraTPCFwdTime[cacheTrk[itrk]]) * 1.E6;     // maximum time in ps
 
-    /*
+    //   printf("trk time %f - %f (max shift +/- %f cm)\n",minTrkTime,maxTrkTime,trackWork.second.getTimeStampError()*vdrift );
+
     for (double tBC = minTrkTime; tBC < maxTrkTime; tBC += BCgranularity) {
       unsigned long ibc = (unsigned long)(tBC * Geo::BC_TIME_INPS_INV);
       BCcand.emplace_back(ibc);
       nStripsCrossedInPropagation.emplace_back(0);
     }
-*/
 
+    /*
     for (auto itof = itof0; itof < nTOFCls; itof++) {
       auto& trefTOF = mTOFClusWork[cacheTOF[itof]];
 
+//     printf("clus time = %f\n",trefTOF.getTime());
+
       if (trefTOF.getTime() < minTrkTime) { // this cluster has a time that is too small for the current track, we will get to the next one
         itof0 = itof + 1;
         continue;
@@ -1203,7 +1212,7 @@ void MatchTOF::doMatchingForTPC(int sec)
         nStripsCrossedInPropagation.emplace_back(0);
       }
     }
-
+*/
     //    printf("BC = %ld\n",BCcand.size());
 
     detId.clear();
@@ -1320,8 +1329,8 @@ void MatchTOF::doMatchingForTPC(int sec)
       }
     }
     for (int ibc = 0; ibc < BCcand.size(); ibc++) {
-      float minTime = (BCcand[ibc] - bc_grouping) * Geo::BC_TIME_INPS;
-      float maxTime = (BCcand[ibc] + bc_grouping) * Geo::BC_TIME_INPS;
+      float minTime = (BCcand[ibc] - bc_grouping_tolerance) * Geo::BC_TIME_INPS;
+      float maxTime = (BCcand[ibc] + bc_grouping_tolerance) * Geo::BC_TIME_INPS;
       for (Int_t imatch = 0; imatch < nStripsCrossedInPropagation[ibc]; imatch++) {
         // we take as residual the average of the residuals along the propagation in the same strip
         deltaPos[ibc][imatch][0] /= nStepsInsideSameStrip[ibc][imatch];
@@ -1340,17 +1349,14 @@ void MatchTOF::doMatchingForTPC(int sec)
         //      printf("itof = %d\n", itof);
         auto& trefTOF = mTOFClusWork[cacheTOF[itof]];
         // compare the times of the track and the TOF clusters - remember that they both are ordered in time!
-        //Printf("trefTOF.getTime() = %f, maxTrkTime = %f, minTrkTime = %f", trefTOF.getTime(), maxTrkTime, minTrkTime);
 
         if (trefTOF.getTime() < minTime) { // this cluster has a time that is too small for the current track, we will get to the next one
-          //Printf("In trefTOF.getTime() < minTrkTime");
           itof0 = itof + 1; // but for the next track that we will check, we will ignore this cluster (the time is anyway too small)
           continue;
         }
         if (trefTOF.getTime() > maxTime) { // no more TOF clusters can be matched to this track
           break;
         }
-
         unsigned long bcClus = trefTOF.getTime() * Geo::BC_TIME_INPS_INV;
 
         int mainChannel = trefTOF.getMainContributingChannel();
@@ -1402,8 +1408,13 @@ void MatchTOF::doMatchingForTPC(int sec)
           LOG(DEBUG) << "Propagated Track [" << itrk << ", " << cacheTrk[itrk] << "]: detId[" << iPropagation << "]  = " << detId[ibc][iPropagation][0] << ", " << detId[ibc][iPropagation][1] << ", " << detId[ibc][iPropagation][2] << ", " << detId[ibc][iPropagation][3] << ", " << detId[ibc][iPropagation][4];
           float resX = deltaPos[ibc][iPropagation][0] - (indices[4] - detId[ibc][iPropagation][4]) * Geo::XPAD + posCorr[0]; // readjusting the residuals due to the fact that the propagation fell in a pad that was not exactly the one of the cluster
           float resZ = deltaPos[ibc][iPropagation][2] - (indices[3] - detId[ibc][iPropagation][3]) * Geo::ZPAD + posCorr[2]; // readjusting the residuals due to the fact that the propagation fell in a pad that was not exactly the one of the cluster
-          resZ += (BCcand[ibc] - bcClus) * vdriftInBC * side;                                                                // add bc correction
+          if (BCcand[ibc] > bcClus) {
+            resZ += (BCcand[ibc] - bcClus) * vdriftInBC * side; // add bc correction
+          } else {
+            resZ -= (bcClus - BCcand[ibc]) * vdriftInBC * side;
+          }
           float res = TMath::Sqrt(resX * resX + resZ * resZ);
+
           if (indices[0] != detId[ibc][iPropagation][0]) {
             continue;
           }
@@ -1413,16 +1424,17 @@ void MatchTOF::doMatchingForTPC(int sec)
           if (indices[2] != detId[ibc][iPropagation][2]) {
             continue;
           }
+
           LOG(DEBUG) << "resX = " << resX << ", resZ = " << resZ << ", res = " << res;
-          float chi2 = res; // TODO: take into account also the time!
+          float chi2 = mIsCosmics ? resX : res; // TODO: take into account also the time!
 
           if (res < mSpaceTolerance) { // matching ok!
             LOG(DEBUG) << "MATCHING FOUND: We have a match! between track " << mTracksSectIndexCache[indices[0]][itrk] << " and TOF cluster " << mTOFClusSectIndexCache[indices[0]][itof];
             foundCluster = true;
             // set event indexes (to be checked)
             evIdx eventIndexTOFCluster(trefTOF.getEntryInTree(), mTOFClusSectIndexCache[indices[0]][itof]);
             evGIdx eventIndexTracks(mCurrTracksTreeEntry, {uint32_t(mTracksSectIndexCache[indices[0]][itrk]), o2::dataformats::GlobalTrackID::TPC});
-            mMatchedTracksPairs.emplace_back(eventIndexTOFCluster, chi2, trkLTInt[ibc][iPropagation], eventIndexTracks); // TODO: check if this is correct!
+            mMatchedTracksPairs.emplace_back(eventIndexTOFCluster, chi2, trkLTInt[ibc][iPropagation], eventIndexTracks, resZ / vdrift * side, trefTOF.getZ()); // TODO: check if this is correct!
           }
         }
       }

Detectors/GlobalTrackingWorkflow/tofworkflow/include/TOFWorkflow/RecoWorkflowWithTPCSpec.h

@@ -19,7 +19,7 @@ namespace o2
 namespace tof
 {
 
-o2::framework::DataProcessorSpec getTOFRecoWorkflowWithTPCSpec(bool useMC, bool useFIT, bool doTPCRefit);
+o2::framework::DataProcessorSpec getTOFRecoWorkflowWithTPCSpec(bool useMC, bool useFIT, bool doTPCRefit, bool iscosmics);
 
 } // end namespace tof
 } // end namespace o2

Detectors/GlobalTrackingWorkflow/tofworkflow/src/RecoWorkflowWithTPCSpec.cxx

@@ -49,9 +49,10 @@ class TOFDPLRecoWorkflowWithTPCTask
   bool mUseMC = true;
   bool mUseFIT = false;
   bool mDoTPCRefit = false;
+  bool mIsCosmics = false;
 
  public:
-  explicit TOFDPLRecoWorkflowWithTPCTask(bool useMC, bool useFIT, bool doTPCRefit) : mUseMC(useMC), mUseFIT(useFIT), mDoTPCRefit(doTPCRefit) {}
+  explicit TOFDPLRecoWorkflowWithTPCTask(bool useMC, bool useFIT, bool doTPCRefit, bool iscosmics) : mUseMC(useMC), mUseFIT(useFIT), mDoTPCRefit(doTPCRefit), mIsCosmics(iscosmics) {}
 
   void init(framework::InitContext& ic)
   {
@@ -74,6 +75,11 @@ class TOFDPLRecoWorkflowWithTPCTask
 
   void run(framework::ProcessingContext& pc)
   {
+    if (mIsCosmics) {
+      mMatcher.setCosmics();
+    }
+    mMatcher.print();
+
     mTimer.Start(false);
     //>>>---------- attach input data --------------->>>
     const auto clustersRO = pc.inputs().get<gsl::span<o2::tof::Cluster>>("tofcluster");
@@ -142,7 +148,7 @@ class TOFDPLRecoWorkflowWithTPCTask
   TStopwatch mTimer;
 };
 
-o2::framework::DataProcessorSpec getTOFRecoWorkflowWithTPCSpec(bool useMC, bool useFIT, bool doTPCRefit)
+o2::framework::DataProcessorSpec getTOFRecoWorkflowWithTPCSpec(bool useMC, bool useFIT, bool doTPCRefit, bool iscosmics)
 {
   std::vector<InputSpec> inputs;
   std::vector<OutputSpec> outputs;
@@ -174,7 +180,7 @@ o2::framework::DataProcessorSpec getTOFRecoWorkflowWithTPCSpec(bool useMC, bool
     "TOFRecoWorkflowWithTPC",
     inputs,
     outputs,
-    AlgorithmSpec{adaptFromTask<TOFDPLRecoWorkflowWithTPCTask>(useMC, useFIT, doTPCRefit)},
+    AlgorithmSpec{adaptFromTask<TOFDPLRecoWorkflowWithTPCTask>(useMC, useFIT, doTPCRefit, iscosmics)},
     Options{
       {"material-lut-path", VariantType::String, "", {"Path of the material LUT file"}}}};
 }

Detectors/GlobalTrackingWorkflow/tofworkflow/src/tof-matcher-tpc.cxx

@@ -52,6 +52,7 @@ void customize(std::vector<o2::framework::ConfigParamSpec>& workflowOptions)
   workflowOptions.push_back(ConfigParamSpec{"tpc-refit", o2::framework::VariantType::Bool, false, {"refit matched TPC tracks"}});
   workflowOptions.push_back(ConfigParamSpec{"disable-root-input", o2::framework::VariantType::Bool, false, {"disable root-files input readers"}});
   workflowOptions.push_back(ConfigParamSpec{"disable-root-output", o2::framework::VariantType::Bool, false, {"disable root-files output writers"}});
+  workflowOptions.push_back(ConfigParamSpec{"cosmics", o2::framework::VariantType::Bool, false, {"reco for cosmics"}});
   workflowOptions.push_back(ConfigParamSpec{"configKeyValues", o2::framework::VariantType::String, "", {"Semicolon separated key=value strings ..."}});
 }
 
@@ -96,6 +97,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
   auto doTPCRefit = cfgc.options().get<bool>("tpc-refit");
   bool disableRootInput = cfgc.options().get<bool>("disable-root-input");
   bool disableRootOutput = cfgc.options().get<bool>("disable-root-output");
+  auto isCosmics = cfgc.options().get<bool>("cosmics");
 
   LOG(INFO) << "TOF RECO WORKFLOW configuration";
   LOG(INFO) << "TOF output = " << cfgc.options().get<std::string>("output-type");
@@ -146,7 +148,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
   }
 
   LOG(INFO) << "Insert TOF Matching";
-  specs.emplace_back(o2::tof::getTOFRecoWorkflowWithTPCSpec(useMC, useFIT, doTPCRefit));
+  specs.emplace_back(o2::tof::getTOFRecoWorkflowWithTPCSpec(useMC, useFIT, doTPCRefit, isCosmics));
 
   if (!disableRootOutput) {
     if (writematching) {

Detectors/TOF/reconstruction/src/Clusterer.cxx

@@ -54,7 +54,7 @@ void Clusterer::calibrateStrip()
     double calib = mCalibApi->getTimeCalibration(dig->getChannel(), dig->getTOT() * Geo::TOTBIN_NS);
     //printf("channel %d) isProblematic = %d, fractionUnderPeak = %f\n",dig->getChannel(),mCalibApi->isProblematic(dig->getChannel()),mCalibApi->getFractionUnderPeak(dig->getChannel())); // toberem
     dig->setIsProblematic(mCalibApi->isProblematic(dig->getChannel()));
-    dig->setCalibratedTime(dig->getTDC() * Geo::TDCBIN + dig->getBC() * o2::constants::lhc::LHCBunchSpacingNS * 1E3 - calib); //TODO:  to be checked that "-" is correct, and we did not need "+" instead :-)
+    dig->setCalibratedTime(dig->getTDC() * Geo::TDCBIN + dig->getBC() * o2::constants::lhc::LHCBunchSpacingNS * 1E3 - Geo::LATENCYWINDOW * 1E3 - calib); //TODO:  to be checked that "-" is correct, and we did not need "+" instead :-)
     //printf("calibration correction = %f\n",calib); // toberem
   }
 }

Detectors/TOF/simulation/src/Digitizer.cxx

@@ -105,7 +105,7 @@ int Digitizer::process(const std::vector<HitType>* hits, std::vector<Digit>* dig
   for (auto& hit : *hits) {
     //TODO: put readout window counting/selection
 
-    processHit(hit, mEventTime.getTimeOffsetWrtBC());
+    processHit(hit, mEventTime.getTimeOffsetWrtBC() + Geo::LATENCYWINDOW);
   } // end loop over hits
 
   if (!mContinuous) { // fill output container per event