Changes to increase HFT hit efficiency (by Hao)

Proposed change 1

divide ladder into 2 halves, which go to 2 Sti layers
Implemented for PXL and IST. The code is on the MAIN branch in CVS. As we agreed I will implement the same for SST too but it is not ready yet. There are issues with StSstFastSimMaker which need to be addressed first.
This feature has been discussed earlier and the general agreement was to implement it. This should prevent tracks from propagating back to Sti volumes at larger radii hence loosing hits
Implemented for PXL and IST. Pending SST
Committed to MAIN branch

Proposed change 2

use central sensor in a ladder to represent the ladder position
Implemented for PXL, IST, and SST
This feature has been implemented for PXL, IST, and SST
Committed to MAIN branch

Proposed change 3

diff --git a/Sti/StiHit.cxx b/Sti/StiHit.cxx
index d0b77be..d8172cc 100644
--- a/Sti/StiHit.cxx
+++ b/Sti/StiHit.cxx
@@ -188,7 +188,7 @@ StiDebug::Break(-46);

   double difX = fabs(mx-posX)-deltaX/2;
   double difY = fabs(my-posY)-deltaY;
-  double difZ = fabs(my-posZ)-deltaZ;
+  double difZ = fabs(mz-posZ)-deltaZ;

   double accX = posX*0.05+1e-3;
   double accY = deltaY*0.1+1e-3;
larger tolerance (>= 1mm) for hits out of volume bundary
is supposed to fix a bug (my -> mz) that is not directly related to the proposed increase in tolerance for hits measured outside of the detector volume?
Victor: It is a bug, for sure, but it was never affected our production. We never see such error messages. Gates in Z are very big.

The change fixes a bug not directly related to HFT hit efficiency
Not committed to MAIN branch

Proposed change 4

diff --git a/Sti/StiHit.cxx b/Sti/StiHit.cxx
index d8172cc..8505c96 100644
--- a/Sti/StiHit.cxx
+++ b/Sti/StiHit.cxx
@@ -190,9 +190,9 @@ StiDebug::Break(-46);
   double difY = fabs(my-posY)-deltaY;
   double difZ = fabs(mz-posZ)-deltaZ;

-  double accX = posX*0.05+1e-3;
-  double accY = deltaY*0.1+1e-3;
-  double accZ = deltaZ*0.1+1e-3;
+  double accX = posX*0.05+0.1;
+  double accY = deltaY*0.1+0.1;
+  double accZ = deltaZ*0.1+0.1;

   if (difX<accX &&  difY<accY && difZ<accZ) return;
   static int counts = 0; counts++; if (counts>100) return;
larger tolerance (>= 1mm) for hits out of volume bundary
I don't see how the change affects the hit efficiency. Looks like only local variables are modified in this routine. Does the increase (0.001 to 0.1 cm) in the minimum allowed distance from the hit to the detector boundary (the accX/Y/Z local variables) affect the decision on whether the considered hit should be picked up by the track or not. To my eye, the change would only affect the number of "hit outside of detector" warnings generated by StiHit::setGlobal(). I don't see how the decision is transferred to either the hit or the track
I think (not very sure after ~month) it's only the inner pxl layer that has problem. poxX~=2 cm. posX*0.05 ~= 1 mm. PXL could have ~mm misalignment. Maybe you are right it's nothing more than get rid of the warnings.
The change does not affect the hit efficiency. Can be implemented to reduce the number of warnings
Committed onto StiHFT_1b branch (Victor)
Not committed to MAIN branch

Proposed change 5
diff --git a/Sti/StiKalmanTrack.cxx b/Sti/StiKalmanTrack.cxx
index 6270843..4e82a3e 100644
--- a/Sti/StiKalmanTrack.cxx
+++ b/Sti/StiKalmanTrack.cxx
@@ -1607,6 +1607,8 @@ static int nCall=0;nCall++;
       if (!targetNode->isValid())   continue;
+    if(targetNode->getHit()) targetNode->nudge(targetNode->getHit());
     status = sTNH.makeFit(0);
     if (status) {restIsWrong = 2005; targetNode->setInvalid();}
@@ -1625,6 +1627,8 @@ static int nCall=0;nCall++;
       if (!targetNode->isValid())   continue;
+    if(targetNode->getHit()) targetNode->nudge(targetNode->getHit());
     status = sTNH.makeFit(1);
     if (status) {restIsWrong = 2005; targetNode->setInvalid();}

nudge to hit local x before fit node clearer logic for active / inactive loop
It is always done in Sti. Tpc hits also sometimes centimetres outside of the volume
Good we all agree
I checked again. This nudge is still not there in dev. Unless we know it's done somewhere else I think it still need to be

Victor: In HFT branch these nudge are allready there. But I am going to remove them. There is a method makeFit is called.
          inside nudge is called.
Hao: I checked this point again by removing this change from my codes and running.


I see a pxl hit not picked up in the test example event:

Eventst_cosmic_adc_15088054_raw_1000006.root has a track with 2 pxl hits (both in layer 2);

Eventst_cosmic_adc_15088054_raw_1000006_noNudge.root has a track with 1 pxl hit. The other pxl hit is cut off by node max chi2 cut in refit(). Larger chi2 is observed without nudge() before makeFit().


By looking at the code, my change use

int StiKalmanTrackNode::nudge(StiHit *hitp)

while in makeFit(), 

int StiTrackNodeHelper::nudge()

is used.

They look like for the same purpose but written in different ways. 

int StiKalmanTrackNode::nudge(StiHit *hitp)     looks more precise.

I didn't dig further. It will be helpful if somebody can test the result of these 2 nudge() function for the same node.

I am not sure whether this is only a matter for overlap (2 hits on one track in one layer as the example here).


I measured efficiency without modification 5 (nudge), all other proposed modification in except the very new "use used hits".
IST: 0.966218 +- 0.00341143
PXL2: 0.900527 +- 0.0136213
PXL1: 0.913769 +- 0.0298183
total: 80%
Comparing with all modifications in except "use used hits":
ist: 0.966231 +- 0.00341087
pxl2: 0.90213 +- 0.013543
pxl1: 1.01406 +- 0.151816
total: 88%
There is a 8% lost, all on PXL layer 1. This is understood, that when nudge on PXL2 is not precise, fit on PXL2 is not precise, and track is not precise enough to find PXL1 hit.
And efficiency loss may even be the less important affect here. The effect on track precision is even more important. When we observe track not fit precisely on layer 2, neither will it be fit precisely on layer 1. This will influence HFT's ability to identify short decay length particles even with tracks that find all the hits.
Some additional comment on this. 
PXL hit can be ~1mm from PXL detector volume plane in Sti. For a typical track, 0.3 GeV/c, R = 2000 mm. Cosmic tracks usually have much higher pT than that. The displacement due to curvature for s=1 mm is 1./2000.*1./2 = 2.5e-4 mm = 0.25 micron. This is much smaller even than PXL hit error, not large enough to explain the efficiency loss.
I think there are 2 possibilities: 1. StiTrackNodeHelper::nudge() not only ignores curvature, but also has some other bug; 2. for some reason node x if far from both the detector plane and the hit before nudge().
But finding this out is not urgent as long as StiKalmanTrackNode::nudge(StiHit *hitp) is always used.


Proposed change 6a, 6b

diff --git a/Sti/StiKalmanTrackNode.cxx b/Sti/StiKalmanTrackNode.cxx
index 30fb9c1..9ed735b 100644
--- a/Sti/StiKalmanTrackNode.cxx
+++ b/Sti/StiKalmanTrackNode.cxx
@@ -484,7 +484,7 @@ using namespace std;
 static const double kMaxEta = 1.25; // 72 degrees for laser tracks
 static const double kMaxSinEta = sin(kMaxEta);
 static const double kMaxCur = 0.2;
-static const double kFarFromBeam = 10.;
+static const double kFarFromBeam = 2.;
 static const Double_t kMaxZ = 250;
 static const Double_t kMaxR = 250;
 StiNodeStat StiKalmanTrackNode::mgP;
@@ -1094,12 +1094,21 @@ StiDebug::StiDebug::Break(nCall);
     mFP._sinCA   = mgP.sinCA2;
     mFP._cosCA   = mgP.cosCA2;
     ians = locate();
-    if (ians<=kEdgeZplus && ians>=0) break;
-  }
-  if (ians>kEdgeZplus || ians<0)                return ians;
-  if (mFP.x()> kFarFromBeam) {
-    if (fabs(mFP.eta())>kMaxEta)           return kEnded;
-    if (mFP.x()*mgP.cosCA2+mFP.y()*mgP.sinCA2<=0)  return kEnded; 
+    if (!ians) break;
+  }  
+  if (ians)                                   return ians;
+  const double tpcInnerRadius = 46.107;
+  if(mFP.x()> kFarFromBeam) {
+      if (fabs(mFP.eta())>kMaxEta)                        
+          {
+              if(mFP.x()<tpcInnerRadius) return -7; 
+              return kEnded;
+          }
+      if (mFP.x()*mgP.cosCA2+mFP.y()*mgP.sinCA2<=0)
+          {
+              if(mFP.x()<tpcInnerRadius) return -8;
+              return kEnded; 
+          }
   mFP.hz()      = getHz();
   if (fabs(mFP.hz()) > 1e-10)       { mFP.curv() = mFP.hz()*mFP.ptin();}

Hao: kFarFromBeam 10 cm -> 2 cm to avoid interfering HFT tracking
        This cut is used in propagate(), the first step for a track to look at any volumes. Basically it's possible that when looking at some pxl volume, the tracker thinks that the track has gone over (0,0) point, but in fact there are other more proper volumes to look at. 

Hao: Avoid terminating the track in HFT region, when angle difference between track and volume is large, or track go through from backside.
Instead continue to look at other possible detetor volumes.

Victor:Termination of track is not related to any detector. It is terminated when
         it is passed DCA point to (0,0). It was made, because there is a danger that
         track near zero could catch another track and will go outside.
         This track is non physical. But there is loopers, which will be cutted.
         It was desided that loopers are not important. In Stv loopers are followed,
         but this problem is not yet solved. We do not know the direction of the looper.
         So I decreased kFarFromBeam 10==> 2. I think it is not dangerous, but I am not
         sure that it will be better

Hao: I didn't think about the problem in the more grand view. I just found how the hits are missed and fixed it. But in view of dealing with the loopers, I think the proper way is the following:
       let the track go from outside-in layer by layer and search for proper volumes, as it is now;
       return a value when propagating has some problem (going from the back side of the volume, large angle ...), but don't kill/end the track;
       ignore the volume that has problem with progagate, don't look for hits in it;
       when all layers of detectors are looped, the track will end naturally.
In this way no wrong hits are included, but also the track is not ended just because it unluckly looked at a wrong volume.


Hao, you misunderstand me. We are tracking from outside towards to beam. We are adding hit by hit to our track. Then our track missed the DCA point to beam and started to move from the beam to outside. Should we pick up hits now?
For loopers yes. But for normal track not. It is not physical track. Thats why tracking is stopped. All other hits do not belong to our track. And not only hits. Prolongation of track after nearest to beam point is not physical

Victor: No track is not ended just because it unluckly looked at a wrong
     volume, it is ended because it started to move from the beam to outside

codes without changes about track end when going to a volume from a wrong direction or with big crossing angle:
ist: 0.341374 +- 0.00893978
pxl2: 0.86833 +- 0.0234455
pxl1: 0.650348 +- 0.0842694
total: 19%

Both these 2 changes increase efficiency.

Victor: OK I understand it now. The logics of stopping tracking is not related to comic tracks.
cosmic track are not coming from the zero.

Proposed change 6c

diff --git a/Sti/StiKalmanTrackNode.cxx b/Sti/StiKalmanTrackNode.cxx
index 9ed735b..147bea6 100644
--- a/Sti/StiKalmanTrackNode.cxx
+++ b/Sti/StiKalmanTrackNode.cxx
@@ -1940,6 +1940,7 @@ int StiKalmanTrackNode::locate()
     yAbsOff = fabs(yOff);
     yAbsOff -=kNStd*sqrt((mFE._cXX+mFE._cYY)/(mFP.x()*mFP.x()+mFP.y()*mFP.y()));
     if (yAbsOff<0) yAbsOff=0;
+    if (yAbsOff > sh->getOpeningAngle()/2) return -1;
     detHW = ((StiCylindricalShape *) sh)->getOpeningAngle()/2.;
     innerY = outerY = detHW;
@@ -1948,6 +1949,7 @@ int StiKalmanTrackNode::locate()
     yOff = mFP.y() - place->getNormalYoffset();
     yAbsOff = fabs(yOff) - kNStd*sqrt(mFE._cYY);
     if (yAbsOff<0) yAbsOff=0;
+    if (yAbsOff > sh->getHalfWidth()) return -1;
     detHW = sh->getHalfWidth();
     innerY = detHW - edge;
     //outerY = innerY + 2*edge;
@@ -1982,6 +1984,12 @@ int StiKalmanTrackNode::locate()
       comment += ::Form(" missed %2d y0/z0 %8.3f/%8.3f dY/dZ %8.3f/%8.3f",
                        position, yOff, zOff, detHW, detHD);
+  const double tpcInnerRadius = 46.107;
+  if(place->getLayerRadius() < tpcInnerRadius)
+      zAbsOff -= kNStd*sqrt(mFE._cZZ);
+  if (zAbsOff < 0) zAbsOff=0;
+  if (zAbsOff > sh->getHalfDepth()) return -1;
+  return 0;
   return position;
Hao: Consider track error when deciding whether a track go through an active volume
Hao: I test efficiency with cosmic for 3 set-ups
codes with all my proposed changes:
ist: 0.966231 +- 0.00341087
pxl2: 0.90213 +- 0.013543
pxl1: 1.01406 +- 0.151816
total: 88%

codes without considering track error when decider whether the track go through an active volume (below):
ist: 0.915096 +- 0.00525854
pxl2: 0.863403 +- 0.015732
pxl1: 1.01618 +- 0.15714
total: 80%
Simulation results wiht pile-up show negative effect on good match efficiency due to this change, shown in Page 2 of
This may be due to increased mismatch. Suggest not to use this change.   


Proposed change 7

diff --git a/Sti/StiTrackNode.cxx b/Sti/StiTrackNode.cxx
index 929e500..e76c435 100644
--- a/Sti/StiTrackNode.cxx
+++ b/Sti/StiTrackNode.cxx
@@ -19,7 +19,7 @@ static const double MAX2ERR[]={MAX1ERR[0]*MAX1ERR[0]

-static const double MIN1ERR[]={1e-4,1e-4,1e-4,1e-4,1e-3,1e-4};
+static const double MIN1ERR[]={1e-4,1e-4,1e-4,1e-5,1e-3,1e-5};
static const double MIN2ERR[]={MIN1ERR[0]*MIN1ERR[0]

static const double recvCORRMAX  = 0.90;
-static const double chekCORRMAX  = 0.99;
+static const double chekCORRMAX  = 0.9999999;
static double MAXPARS[]={250,250,250,1.5,100,100};
Hao: smaller minimal and larger maximum tracking errors due to HFT 
Victor: It is not clear why max errs should be increased for HFT. Max errors were introduced, because error propagation is based in linear approximation. Sometimes it works too incorrectly, and errors became too big. In this case, track accepts wrong hit and in result is failed. Same for too small errors.
You propose to decrease minimal accuracy for azimutal angle from 0.1 milliradian to 0.01 milliradian and the same for tan(lambda) from 1e-4 to 1e-5. It is accuracy 10 microns on the length of 1 meter.
Are you sure that HFT is so fantastically good?

Dmitri:  Is this patch complete? I don't see where you increase the maximum errors

Hao: I copied the changes on chekCORRMAX above. I think this is cut on covariance matrix, not absolute error value. Increasing it allow stronger correlation between 2 fitting parameters. Maybe it's already done for the version of code you compared with.
For the two angle error cut, pixel hit error is 20 microns, distance between inner and outer pxl layer is 6 cm. The error on angle for a track along R direction is 2.e-3/6 = 3.e-4. When a track go through the detector plane with a large angle, error perpendicular to track direction is smaller (because we assign hit errors only along the two directions on the detector plane), angle error is possible to go below 1.e-4. 
In general I got all these proposed changes based on running cosmic data, that tracks are lost due to these cuts. I found and modified them to let tracks go through. I didn't imagine them, data taught me.

Victor      The error on angle for a track along R direction is (20+20)*1e-4/6 =
          6e-4. So it is still bigger than 1e-4. But it is close. So I can
      decrease it to 1e-5.
      "tracks are lost due to these cuts" It is strange, it is not a cut.
      If errors less than minimal, errors increased. Why you lost?
      About chekCORRMAX.You probably look on wrong code.
      StiTrackNode.cxx:static const double chekCORRMAX  = 0.9999;
          StiTrackNode.cxx:static const double recvCORRMAX  = 0.999;
I plot the distribution of values that are cut on.

 is distribution of sqrt(A[idx66[3][3]]), being cut by MINIERR[3]. We can
see 1.e-4 is not small enough.
Same for
 distribution of sqrt(A[idx66[5][5]]), being cut by MINIERR[5]. When 2 pxl
hits with 20 micron error fit together with ~40 TPC hits > 1 m away, such
angle resolution can be achieved.

 are two example distributions of aij*aij/(aii*ajj) that can be very close
to 1, meaning two track parameters are strongly correlated.
 plot 1-aij*aij/(aii*ajj) in log scale to see more clearly. Proposed new
cut placed on 1-1.e-7 is proper.


Proposed change 8

diff --git a/Sti/StiTrackNode.cxx b/Sti/StiTrackNode.cxx
index 929e500..2a6cd44 100644
--- a/Sti/StiTrackNode.cxx
+++ b/Sti/StiTrackNode.cxx
@@ -404,12 +404,12 @@ RETN:
   if (!pri )   return kase;
   switch(kase) {

-    case 1: LOG_DEBUG << Form("StiNodeErrs::check(%s) FAILED: Negative diagonal %g[%d][%d]",pri,aii,i,i)<< endm;  
+    case 1: LOG_WARN << Form("StiNodeErrs::check(%s) FAILED: Negative diagonal %g[%d][%d]",pri,aii,i,i)<< endm;  
-    case 2: LOG_DEBUG << Form("StiNodeErrs::check(%s) FAILED: Correlation too big %g[%d][%d]>%g"
+    case 2: LOG_WARN << Form("StiNodeErrs::check(%s) FAILED: Correlation too big %g[%d][%d]>%g"
-    case 3: LOG_DEBUG << Form("StiNodeErrs::check(%s) FAILED: Non Positive matrix",pri)<<endm;  
+    case 3: LOG_WARN << Form("StiNodeErrs::check(%s) FAILED: Non Positive matrix",pri)<<endm;  
 assert(sign()>0); ///???
   return kase;

Hao: DEBUG -> WARN when error check failed, so that no track is lost silently for this

Victor: This is used for debug only. In real life it happens when hits are wrong and track becames crazy. It happens rather often and output file will be overfull by not important information

Dmitri: I also don't see much benefit in changing these messages to warnings. They clearly will not increase the hit efficiency :) and seem to be for experts only so, should not be written into log files during normal event reconstruction.

         Do not commit


Proposed change 9

diff --git a/Sti/StiKalmanTrackFinder.cxx b/Sti/StiKalmanTrackFinder.cxx
index 9737e2e..44c3caa 100644
--- a/Sti/StiKalmanTrackFinder.cxx
+++ b/Sti/StiKalmanTrackFinder.cxx
@@ -623,9 +623,20 @@ static int nKount = 0; nKount++;
     if (debug() > 2 && nDets==0) cout << "no detector of interest on this layer"<<endl;
     if (!nDets) continue;
     if (nDets>1) sort(detectors.begin(),detectors.end(),CloserAngle(projAngle) );
+//            There is additional loop. 1st loop for active only, second for non active
+    int foundInDetLoop = 0;
+static int activeNonActiveLoop = StiDebug::iFlag("activeNonActiveLoop");
+    activeNonActiveLoop = 1;
+    for (int nowActive=activeNonActiveLoop; nowActive>=0; nowActive--) { //Additional activeNonActive loop
     for (vector<StiDetector*>::const_iterator d=detectors.begin();d!=detectors.end();++d)
       tDet = *d;
+      if ((tDet->isActive() != nowActive) && activeNonActiveLoop ) continue;
       if (debug() > 2) {
        cout << endl<< "target det:"<< *tDet;
        cout << endl<< "lead angle:" << projAngle*radToDeg

Hao: nudge to hit local x before fit node

         clearer logic for active / inactive loop

        getCenterRefAngle() instead of getNormalRefAngle as a better estimate of volume phi position.

       Victor commentted nudge is done with the new makeFit function. If so these nudge() are no longer needed. But need to check with the new codes.
       Victor also state in an email that active-inactive volume loop are in StiHFT_1b branch. Need to make sure it's in the new codes.


Yes, getNormalRefAngle() is not good. But to use getCenterAngle() better but still wrong. Objects are sorted by LayerAngle, provided by user in geometry definition. So correctly to use: getLayerRaduius() & getLayerAngle.
In double search loop you removed test for environment variable. In future it will be removed, but not now.
Also you fix a bug. Yes, index must start from 1. For single loop it is not important but for double one, it is a bug. node->nudge(stiHit);      

 node->nudge(stiHit);  In standard Sti it was node->nudge()  

 And it was correct. Now nudge slightly modified and I have in branch call with a hit as in your case.

Looks like there is an agreement on some of the changes in this patch. A cleaner patch needs to be prepared for the MAIN branch
Not committed to MAIN branch

Proposed change 10

diff --git a/Sti/StiHitContainer.cxx b/Sti/StiHitContainer.cxx
index 9759e4f..6e381c7 100644
--- a/Sti/StiHitContainer.cxx
+++ b/Sti/StiHitContainer.cxx
@@ -22,6 +22,7 @@ using std::stable_partition;
 ostream& operator<<(ostream& os, const StiHit& hit);
 ostream& operator<<(ostream&, const HitMapKey&);

+static const double kRMinTpc =55;
 int VectorAndEnd::fIdCounter = 0;
 VectorAndEnd::VectorAndEnd(): fEffectiveEndValid(false)
@@ -238,7 +239,8 @@ vector<StiHit*> & StiHitContainer::getHits(StiHit& ref, double dY, double dZ, bo
       hit = *cit;
       if (fabs( hit->y() - ref.y() ) < dY)
-        if (fetchAll || (hit->timesUsed()==0 && hit->detector()->isActive()) )
+            //          if (fetchAll || (hit->timesUsed()==0 && hit->detector()->isActive()) )
+            if (fetchAll || ((hit->timesUsed()==0 || hit->x()<kRMinTpc) && hit->detector()->isActive()) )
diff --git a/Sti/StiKalmanTrack.cxx b/Sti/StiKalmanTrack.cxx
index 4e82a3e..18d4fc4 100644
--- a/Sti/StiKalmanTrack.cxx
+++ b/Sti/StiKalmanTrack.cxx
@@ -1520,7 +1520,7 @@ int StiKalmanTrack::refit()
     if (worstNode && worstNode->getChi2()>StiKalmanTrackFitterParameters::instance()->getMaxChi2())
       worstNode->setHit(0); worstNode->setChi2(3e33); continue;}
-    if (rejectByHitSet()) { releaseHits()            ; continue;}
+    //    if (rejectByHitSet()) { releaseHits()            ; continue;}

     if (!fail)                                                         break;

@@ -1558,12 +1558,10 @@ int StiKalmanTrack::refit()
       if (node == vertexNode)                          continue;
       StiHit *hit = node->getHit();
       if(!hit)                                                 continue;
-      hit->setTimesUsed(0);
       if (!node->isValid())                                 continue;
       if (node->getChi2()>10000.)                        continue;
-      hit->setTimesUsed(1);
@@ -1874,7 +1872,8 @@ int StiKalmanTrack::releaseHits(double rMin,double rMax)
     if (hit->x()>rMax)           break;
-    hit->setTimesUsed(0);
+    if(hit->timesUsed()>0)
+        hit->setTimesUsed(hit->timesUsed()-1);
   return sum;
Hao: Use HFT hits that are used by earlier reconstructed tracks in tracking of later tracks. When a hit is within tpc min radius, even it's used, it's considered for the current track.
The idea is to get back the efficiency loss due to hits used by an early track wrongly (mismatch). For HFT, layers are far away from each other, and the hit density is high, mismatch especially for low pt, low nHits tracks are common. Also, all HFT layers are necessary in order to have a good track fitting. It's not like TPC that missing several layers doesn't matter. And since Sti tracks from outside in, no additional tracks will appear due to this change, only tracks have higher propability to pick up a HFT hit.
A test sample of 10 events has number of tracks with hits in all HFT layers increasing from 938 to 1122. This ~20% increase may not be all good matchings. 
Simulation results show that this change increase good match efficiency by near 10%, as shown in Page 3 of