TRD: Introduce strict matching mode

martenole · shahor02 · commit 356097fa0de2 · 2021-08-10T20:15:37.000+02:00
diff --git a/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx b/Detectors/TRD/workflow/src/TRDGlobalTrackingSpec.cxx
@@ -85,6 +85,9 @@ void TRDGlobalTracking::init(InitContext& ic)
 
   mTracker = new GPUTRDTracker();
   mTracker->SetNCandidates(mRec->GetProcessingSettings().trdNCandidates); // must be set before initialization
+  if (mStrict && mRec->GetProcessingSettings().trdNCandidates == 1) {
+    LOG(ERROR) << "Strict matching mode requested, but tracks with another close hypothesis will not be rejected. Please set trdNCandidates to at least 3.";
+  }
   mTracker->SetProcessPerTimeFrame(true);
   mTracker->SetGenerateSpacePoints(false); // set to true to force space point calculation by the TRD tracker itself
 
@@ -228,6 +231,10 @@ void TRDGlobalTracking::run(ProcessingContext& pc)
       // skip tracks without TRD tracklets (the collision ID for the TRD tracks is initialized to -1 and only changed if a tracklet is attached to the track)
       continue;
     }
+    if (mStrict && (trdTrack.getIsAmbiguous() || trdTrack.getReducedChi2() > mTracker->Param().rec.trd.chi2StrictCut)) {
+      // skip tracks which have another hypothesis close to the best one or which do are above strict chi2 threshold
+      continue;
+    }
     nTrackletsAttached += trdTrack.getNtracklets();
     auto trackGID = trdTrack.getRefGlobalTrackId();
     if (trackGID.includesDet(GTrackID::Source::ITS)) {
diff --git a/GPU/GPUTracking/DataTypes/GPUTRDTrack.cxx b/GPU/GPUTracking/DataTypes/GPUTRDTrack.cxx
@@ -33,7 +33,7 @@ GPUd() void GPUTRDTrack_t<T>::initialize()
   mChi2 = 0.f;
   mRefGlobalTrackId = 0;
   mCollisionId = -1;
-  mIsFindable = 0;
+  mFlags = 0;
   for (int i = 0; i < kNLayers; ++i) {
     mAttachedTracklets[i] = -1;
   }
@@ -80,7 +80,7 @@ GPUd() void GPUTRDTrack_t<T>::ConvertFrom(const GPUTRDTrackDataRecord& t)
   T::set(t.fX, t.mAlpha, &(t.fY), t.fC);
   setRefGlobalTrackIdRaw(t.fTPCTrackID);
   mChi2 = 0.f;
-  mIsFindable = 0;
+  mFlags = 0;
   mCollisionId = -1;
   for (int iLayer = 0; iLayer < kNLayers; iLayer++) {
     mAttachedTracklets[iLayer] = t.fAttachedTracklets[iLayer];
@@ -109,7 +109,7 @@ GPUd() GPUTRDTrack_t<T>::GPUTRDTrack_t(const o2::tpc::TrackTPC& t) : T(t)
 
 template <typename T>
 GPUd() GPUTRDTrack_t<T>::GPUTRDTrack_t(const GPUTRDTrack_t<T>& t)
-  : T(t), mChi2(t.mChi2), mRefGlobalTrackId(t.mRefGlobalTrackId), mCollisionId(t.mCollisionId), mIsFindable(t.mIsFindable)
+  : T(t), mChi2(t.mChi2), mRefGlobalTrackId(t.mRefGlobalTrackId), mCollisionId(t.mCollisionId), mFlags(t.mFlags)
 {
   // copy constructor
   for (int i = 0; i < kNLayers; ++i) {
@@ -135,7 +135,7 @@ GPUd() GPUTRDTrack_t<T>& GPUTRDTrack_t<T>::operator=(const GPUTRDTrack_t<T>& t)
   mChi2 = t.mChi2;
   mRefGlobalTrackId = t.mRefGlobalTrackId;
   mCollisionId = t.mCollisionId;
-  mIsFindable = t.mIsFindable;
+  mFlags = t.mFlags;
   for (int i = 0; i < kNLayers; ++i) {
     mAttachedTracklets[i] = t.mAttachedTracklets[i];
   }
@@ -148,7 +148,7 @@ GPUd() int GPUTRDTrack_t<T>::getNlayersFindable() const
   // returns number of layers in which the track is in active area of TRD
   int retVal = 0;
   for (int iLy = 0; iLy < kNLayers; iLy++) {
-    if ((mIsFindable >> iLy) & 0x1) {
+    if ((mFlags >> iLy) & 0x1) {
       ++retVal;
     }
   }
diff --git a/GPU/GPUTracking/DataTypes/GPUTRDTrack.h b/GPU/GPUTracking/DataTypes/GPUTRDTrack.h
@@ -62,6 +62,7 @@ class GPUTRDTrack_t : public T
  public:
   enum EGPUTRDTrack {
     kNLayers = 6,
+    kAmbiguousFlag = 6,
     kStopFlag = 7
   };
 
@@ -76,7 +77,7 @@ class GPUTRDTrack_t : public T
   GPUd() GPUTRDTrack_t(const T& t);
   GPUd() GPUTRDTrack_t& operator=(const GPUTRDTrack_t& t);
 
-  // attach a tracklet to this track; this overwrites the mIsFindable flag to true for this layer
+  // attach a tracklet to this track; this overwrites the mFlags flag to true for this layer
   GPUd() void addTracklet(int iLayer, int idx) { mAttachedTracklets[iLayer] = idx; }
 
   // getters
@@ -89,8 +90,9 @@ class GPUTRDTrack_t : public T
   GPUd() int getNtracklets() const;
   GPUd() float getChi2() const { return mChi2; }
   GPUd() float getReducedChi2() const { return getNlayersFindable() == 0 ? mChi2 : mChi2 / getNlayersFindable(); }
-  GPUd() bool getIsStopped() const { return (mIsFindable >> kStopFlag) & 0x1; }
-  GPUd() bool getIsFindable(int iLayer) const { return (mIsFindable >> iLayer) & 0x1; }
+  GPUd() bool getIsStopped() const { return (mFlags >> kStopFlag) & 0x1; }
+  GPUd() bool getIsAmbiguous() const { return (mFlags >> kAmbiguousFlag) & 0x1; }
+  GPUd() bool getIsFindable(int iLayer) const { return (mFlags >> iLayer) & 0x1; }
   GPUd() int getNmissingConsecLayers(int iLayer) const;
   // for AliRoot compatibility. To be removed once HLT/global/AliHLTGlobalEsdConverterComponent.cxx does not require them anymore
   GPUd() int GetTPCtrackId() const { return mRefGlobalTrackId; }
@@ -102,8 +104,9 @@ class GPUTRDTrack_t : public T
   // This method is only defined in TrackTRD.h and is intended to be used only with that TRD track type
   GPUd() void setRefGlobalTrackId(o2::dataformats::GlobalTrackID id);
   GPUd() void setCollisionId(short id) { mCollisionId = id; }
-  GPUd() void setIsFindable(int iLayer) { mIsFindable |= (1U << iLayer); }
-  GPUd() void setIsStopped() { mIsFindable |= (1U << kStopFlag); }
+  GPUd() void setIsFindable(int iLayer) { mFlags |= (1U << iLayer); }
+  GPUd() void setIsStopped() { mFlags |= (1U << kStopFlag); }
+  GPUd() void setIsAmbiguous() { mFlags |= (1U << kAmbiguousFlag); }
   GPUd() void setChi2(float chi2) { mChi2 = chi2; }
 
   // conversion to / from HLT track structure (only for AliRoot)
@@ -115,12 +118,12 @@ class GPUTRDTrack_t : public T
   unsigned int mRefGlobalTrackId;   // raw GlobalTrackID of the seeding track (either ITS-TPC or TPC)
   int mAttachedTracklets[kNLayers]; // indices of the tracklets attached to this track; -1 means no tracklet in that layer
   short mCollisionId;               // the collision ID of the tracklets attached to this track; is used to retrieve the BC information for this track after the tracking is done
-  unsigned char mIsFindable;        // bitfield; LSB indicates whether track is findable in layer 0; MSB flags whether the track is stopped in the TRD; one bit is currently not used
+  unsigned char mFlags;             // bits 0 to 5 indicate whether track is findable in layer 0 to 5, bit 6 indicates an ambiguous track and bit 8 flags if the track is stopped in the TRD
 
  private:
   GPUd() void initialize();
 #if !defined(GPUCA_STANDALONE) && !defined(GPUCA_ALIROOT_LIB)
-  ClassDefNV(GPUTRDTrack_t, 1);
+  ClassDefNV(GPUTRDTrack_t, 2);
 #endif
 };
 
diff --git a/GPU/GPUTracking/Definitions/GPUSettingsList.h b/GPU/GPUTracking/Definitions/GPUSettingsList.h
@@ -83,6 +83,8 @@ BeginSubConfig(GPUSettingsRecTRD, trd, configStandalone.rec, "RECTRD", 0, "Recon
 AddOptionRTC(minTrackPt, float, .5f, "", 0, "Min Pt for tracks to be propagated through the TRD")
 AddOptionRTC(maxChi2, float, 15.f, "", 0, "Max chi2 for TRD tracklets to be matched to a track")
 AddOptionRTC(penaltyChi2, float, 12.f, "", 0, "Chi2 penalty for no available TRD tracklet (effective chi2 cut value)")
+AddOptionRTC(chi2StrictCut, float, 10.f, "", 0, "Chi2 cut for strict matching mode")
+AddOptionRTC(chi2SeparationCut, float, 3.f, "", 0, "Minimum difference between chi2 of winner match and chi2 of second best match")
 AddOptionRTC(nSigmaTerrITSTPC, float, 4.f, "", 0, "Number of sigmas for ITS-TPC track time error estimate")
 AddOptionRTC(stopTrkAfterNMissLy, unsigned char, 6, "", 0, "Abandon track following after N layers without a TRD match")
 AddHelp("help", 'h')
diff --git a/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx b/GPU/GPUTracking/TRDTracking/GPUTRDTracker.cxx
@@ -454,7 +454,10 @@ GPUd() void GPUTRDTracker_t<TRDTRK, PROP>::DoTrackingThread(int iTrk, int thread
     auto trkCopy = trkStart;
     prop.setTrack(&trkCopy);
     prop.setFitInProjections(true);
-    FollowProlongation(&prop, &trkCopy, iTrk, threadId, collisionIds[iColl]);
+    if (!FollowProlongation(&prop, &trkCopy, iTrk, threadId, collisionIds[iColl])) {
+      // track following failed
+      continue;
+    }
     if (trkCopy.getReducedChi2() < mTracks[iTrk].getReducedChi2()) {
       mTracks[iTrk] = trkCopy; // copy back the resulting track
     }
@@ -585,7 +588,8 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
     mCandidates[candidateIdxOffset] = *t;
   }
 
-  int nCandidates = 1;
+  int nCandidates = 1;     // we always start with one candidate
+  int nCurrHypothesis = 0; // the number of track hypothesis in given iLayer
 
   // search window
   float roadY = 0.f;
@@ -595,7 +599,7 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
   mDebug->SetGeneralInfo(mNEvents, mNTracks, iTrk, t->getPt());
 
   for (int iLayer = 0; iLayer < kNLayers; ++iLayer) {
-    int nCurrHypothesis = 0;
+    nCurrHypothesis = 0;
     bool isOK = false; // if at least one candidate could be propagated or the track was stopped this becomes true
     int currIdx = candidateIdxOffset + iLayer % 2;
     int nextIdx = candidateIdxOffset + (iLayer + 1) % 2;
@@ -745,12 +749,9 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
       if (mHypothesis[iUpdate + hypothesisIdxOffset].mCandidateId == -1) {
         // no more candidates
         if (iUpdate == 0) {
-          if (ENABLE_WARNING) {
-            GPUWarning("No valid candidates for track %i in layer %i", iTrk, iLayer);
-          }
-          nCandidates = 0;
+          return false; // no valid candidates for this track (probably propagation failed)
         }
-        break;
+        break; // go to next layer
       }
       nCandidates = iUpdate + 1;
       if (mNCandidates > 1) {
@@ -876,6 +877,11 @@ GPUd() bool GPUTRDTracker_t<TRDTRK, PROP>::FollowProlongation(PROP* prop, TRDTRK
   if (ENABLE_INFO) {
     GPUInfo("Ended track following for track %i at x=%f with pt=%f. Attached %i tracklets", t->getRefGlobalTrackIdRaw(), t->getX(), t->getPt(), t->getNtracklets());
   }
+  if (nCurrHypothesis > 1) {
+    if (CAMath::Abs(mHypothesis[hypothesisIdxOffset + 1].GetReducedChi2() - mHypothesis[hypothesisIdxOffset].GetReducedChi2()) < Param().rec.trd.chi2SeparationCut) {
+      t->setIsAmbiguous();
+    }
+  }
   return true;
 }
 
