Fixes from Marian:
authorshahoian <ruben.shahoyan@cern.ch>
Sat, 3 May 2014 15:52:16 +0000 (17:52 +0200)
committerhristov <Peter.Hristov@cern.ch>
Thu, 29 May 2014 14:36:01 +0000 (16:36 +0200)
AliTRDseedV1.cxx    - adding correction for pad crossing 0 with effective factor for incoplete model
AliTRDtrackerV1.cxx  - adding chi2 to the streamer
AliTRDrecoParam.h AliTRDrecoParam.cxx AliTRDtrackerV1.cxx  - parameterised chi2 cut /tracklet to track
AliTRDseedV1.cxx - increase error estimate for crossing tracklets  -taking incount unisochronity, and crosscorelation of errors

TRD/AliTRDrecoParam.cxx
TRD/AliTRDrecoParam.h
TRD/AliTRDseedV1.cxx
TRD/AliTRDtrackerV1.cxx

index 0716dc7..6d9f598 100644 (file)
@@ -52,6 +52,7 @@ AliTRDrecoParam::AliTRDrecoParam()
   ,fkChi2Y(.25)
   ,fkChi2YSlope(7.73)
   ,fkChi2ZSlope(0.069)
+  ,fChi2Cut(25)
   ,fkChi2YCut(0.5)
   ,fkPhiSlope(10.6)
   ,fkNMeanClusters(20.)
@@ -121,6 +122,7 @@ AliTRDrecoParam::AliTRDrecoParam(const AliTRDrecoParam &ref)
   ,fkChi2Y(ref.fkChi2Y)
   ,fkChi2YSlope(ref.fkChi2YSlope)
   ,fkChi2ZSlope(ref.fkChi2ZSlope)
+  ,fChi2Cut(ref.fChi2Cut)
   ,fkChi2YCut(ref.fkChi2YCut)
   ,fkPhiSlope(ref.fkPhiSlope)
   ,fkNMeanClusters(ref.fkNMeanClusters)
@@ -175,6 +177,7 @@ AliTRDrecoParam& AliTRDrecoParam::operator=(const AliTRDrecoParam &ref)
   fkChi2Y               = ref.fkChi2Y;
   fkChi2YSlope          = ref.fkChi2YSlope;
   fkChi2ZSlope          = ref.fkChi2ZSlope;
+  fChi2Cut            = ref.fChi2Cut;
   fkChi2YCut            = ref.fkChi2YCut;
   fkPhiSlope            = ref.fkPhiSlope;
   fkNMeanClusters       = ref.fkNMeanClusters;
@@ -269,6 +272,7 @@ AliTRDrecoParam *AliTRDrecoParam::GetCosmicTestParam()
   par->fSysCovMatrix[1] = 2.; // z direction (1 cm)
   par->fkChi2YSlope     = 0.11853;
   par->fkChi2ZSlope     = 0.04527;
+  par->fkChi2YCut       = 25.;
   par->fkChi2YCut       = 1.;
   par->fkPhiSlope       = 10.; //3.17954;
   par->fkMaxTheta       = 2.1445;
index 2c6b587..78cbf45 100644 (file)
@@ -56,6 +56,7 @@ public:
   Double_t GetChi2Z() const                 { return fkChi2Z;    }
   Double_t GetChi2YSlope() const            { return fkChi2YSlope; }
   Double_t GetChi2ZSlope() const            { return fkChi2ZSlope; }
+       Double_t GetChi2Cut() const              { return fChi2Cut; }
        Double_t GetChi2YCut() const              { return fkChi2YCut; }
   Double_t GetPhiSlope() const              { return fkPhiSlope;   }
   Float_t  GetNClusters() const;
@@ -128,6 +129,7 @@ public:
   void     SetChi2Z(Double_t chi2)                            {fkChi2Z = chi2;}
   void     SetChi2YSlope(Double_t chi2YSlope)                 {fkChi2YSlope = chi2YSlope;}
   void     SetChi2ZSlope(Double_t chi2ZSlope)                 {fkChi2ZSlope = chi2ZSlope;}
+       void       SetChi2Cut(Double_t chi2Cut)                      {fChi2Cut = chi2Cut; }
        void       SetChi2YCut(Double_t chi2Cut)                      {fkChi2YCut = chi2Cut; }
   void     SetPhiSlope(Double_t phiSlope)                     {fkPhiSlope = phiSlope;}
   void     SetNMeanClusters(Double_t meanNclusters)           {fkNMeanClusters = meanNclusters;}
@@ -170,6 +172,7 @@ private:
   Double_t  fkChi2Y;                 // Max chi2 on the y direction for seeding clusters Rieman fit
   Double_t  fkChi2YSlope;            // Slope of the chi2-distribution in y-direction
   Double_t  fkChi2ZSlope;            // Slope of the chi2-distribution in z-direction
+  Double_t  fChi2Cut;               // Cut on the Chi2 track/tracklet 0 used to diecide if the kalman track should be updated
   Double_t  fkChi2YCut;                                                         // Cut on the Chi2 in y-direction in the likelihood filter
   Double_t  fkPhiSlope;              // Slope of the distribution of the deviation between track angle and tracklet angle
   Double_t  fkNMeanClusters;         // Mean number of clusters per tracklet
index dfa0acd..01002a5 100644 (file)
@@ -638,6 +638,14 @@ Double_t AliTRDseedV1::EstimatedCrossPoint(AliTRDpadPlane *pp)
 //  fS2Z     = 0.05+0.4*TMath::Abs(fZfit[1]); 
   Float_t s2dzdx(0.0245-0.0014*fZfit[1]+0.0557*fZfit[1]*fZfit[1]); s2dzdx*=s2dzdx;
   fS2Z  = fZfit[1]*fZfit[1]*sx*sx+fS2Y*fS2Y*s2dzdx; 
+  //
+  // MI modification 02.05.2014 -  add wire step + unisochronity contribution + diffusion contribution
+  // Normalization factor for fS2Z+=(factor*wirepitch/sqrt(12))^2
+  // MI question - fS2Y definetion looks strange. Is it  used somewhere at all???
+  //
+  const Double_t kWirePitch=0.125;
+  const Double_t kWirePitchFactor=12;
+  fS2Z+=kWirePitch*kWirePitch/kWirePitchFactor;
 
   return fS2Y;
 }
@@ -922,16 +930,22 @@ void AliTRDseedV1::GetCovAt(Double_t x, Double_t *cov) const
   }
 
   // rotate covariance matrix if no RC
+  Double_t t2 = GetTilt()*GetTilt();
   if(!IsRowCross()){
-    Double_t t2 = GetTilt()*GetTilt();
     Double_t correction = 1./(1. + t2);
     cov[0] = (sy2+t2*sz2)*correction;
     cov[1] = GetTilt()*(sz2 - sy2)*correction;
     cov[2] = (t2*sy2+sz2)*correction;
    } else {
      cov[0] = sy2; cov[1] = 0.; cov[2] = sz2;
+    //
+    // MI change (03.05.2014)  
+    // current implemenatation of the linear fit does not take into account step for padrow crossing
+    // bug to be fixed - step should be taken properly into account in the fit, rather than increase of the error estimate
+     cov[0]+=t2*GetPadLength()*GetPadLength()/24.; //this line will disappear after fix of linear FIT 
    }
 
+
   AliDebug(4, Form("C(%6.1f %+6.3f %6.1f)  RC[%c]", 1.e4*TMath::Sqrt(cov[0]), cov[1], 1.e4*TMath::Sqrt(cov[2]), IsRowCross()?'y':'n'));
 }
 
@@ -1368,7 +1382,7 @@ Bool_t  AliTRDseedV1::AttachClusters(AliTRDtrackingChamber *const chamber, Bool_
   // initialize debug streamer
   TTreeSRedirector *pstreamer(NULL);
   if((recoParam->GetStreamLevel(AliTRDrecoParam::kTracker) > 3 && fkReconstructor->IsDebugStreaming())||
-     AliTRDReconstructor::GetStreamLevel()>3) pstreamer = fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
+     AliTRDReconstructor::GetStreamLevel()>30) pstreamer = fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
   if(pstreamer){
     // save config. for calibration
     TVectorD vdy[2], vdx[2], vs2[2];
@@ -2266,9 +2280,25 @@ Bool_t AliTRDseedV1::FitRobust(AliTRDpadPlane *pp, Int_t opt)
   fYfit[1] = -par[1];
   //printf(" yfit: %f [%f] x[%e] dydx[%f]\n", fYfit[0], par[0], fX, par[1]);
   // store covariance
-  fCov[0] = kScalePulls*cov[0]; // variance of y0
-  fCov[1] = kScalePulls*cov[2]; // covariance of y0, dydx
-  fCov[2] = kScalePulls*cov[1]; // variance of dydx
+
+   //
+  // MI comment modification (03.05.2014)
+  // In currently used code variance is used to estimate the tracklet error 
+  //     see AliTRDtrackerV1::AliTRDLeastSquare::Eval()
+  // 2 problems:
+  //    a.) Assumption about the intdpendent error not valid beacause of time response function 
+  //    b.) Constuction of cluster errors used in cl->GetSigmaY2() and AliTRDtrackerV1::AliTRDLeastSquare::Eval
+  //        to be parameterized
+  //    c.) For crossing pad-rows the linear fit AliTRDtrackerV1::AliTRDLeastSquare::Eval do not take into account step
+  //        Therefore y possition is biased, bias is y dependent
+  // Effective correction factor used to rescale error estimates
+  // In future additional angular dependent error should be investigated ( as in TPC and in previous versions of the TRD tracking) 
+  // 
+  const Double_t kCorrelError=1.5;
+
+  fCov[0] = kCorrelError*kScalePulls*cov[0]; // variance of y0
+  fCov[1] = kCorrelError*kScalePulls*cov[2]; // covariance of y0, dydx
+  fCov[2] = kCorrelError*kScalePulls*cov[1]; // variance of dydx
   // check radial position
   Float_t xs=fX+.5*AliTRDgeometry::CamHght();
   if(xs < 0. || xs > AliTRDgeometry::CamHght()+AliTRDgeometry::CdrHght()){
index 1dc6365..0a1bdca 100644 (file)
@@ -624,7 +624,34 @@ Int_t AliTRDtrackerV1::FollowProlongation(AliTRDtrackV1 &t)
     Double_t cov[3]; tracklet->GetCovAt(x, cov);
     Double_t p[2] = { tracklet->GetY(), tracklet->GetZ()};
     Double_t chi2 = ((AliExternalTrackParam)t).GetPredictedChi2(p, cov);
-    if (chi2 < 1e+10 && ((AliExternalTrackParam&)t).Update(p, cov)){ 
+
+    if(fkReconstructor->IsDebugStreaming()){
+      Int_t eventNumber = AliTRDtrackerDebug::GetEventNumber();
+      TTreeSRedirector &cstreamer = *fkReconstructor->GetDebugStream(AliTRDrecoParam::kTracker);
+      AliExternalTrackParam param0(t);
+      AliExternalTrackParam param1(t);
+      param1.Update(p, cov);
+      TVectorD vcov(3,cov);
+      TVectorD vpar(3,p);
+      cstreamer << "FollowProlongationInfo"
+               << "EventNumber="       << eventNumber
+               << "iplane="<<iplane
+               << "vcov.="<<&vcov
+               << "vpar.="<<&vpar
+               << "tracklet.="      << tracklet
+               << "chi2="<< chi2 
+               << "param0.="           << &param0
+               << "param1.="           << &param1
+               << "\n";
+    }
+    /*
+    AliInfo(Form("Pl:%d X:%+e : %+e P: %+e %+e Cov:%+e %+e %+e -> dXY: %+e %+e | chi2:%.2f pT:%.2f alp:%.3f",
+                iplane,x,t.GetX(),p[0],p[1],cov[0],cov[1],cov[2],
+                p[0]-t.GetY(),p[1]-t.GetZ(),
+                chi2,t.Pt()*t.Charge(),t.GetAlpha()));
+    */
+    if (chi2 < fkRecoParam->GetChi2Cut() && ((AliExternalTrackParam&)t).Update(p, cov)){  // MI parameterizad chi2 cut 03.05.2014
+      //    if (chi2 < 1e+10 && ((AliExternalTrackParam&)t).Update(p, cov)){ 
       // Register info to track
       t.SetNumberOfClusters();
       t.UpdateChi2(chi2);
@@ -985,11 +1012,12 @@ Int_t AliTRDtrackerV1::FollowBackProlongation(AliTRDtrackV1 &t)
                << "param0.="           << &param0
                << "param1.="           << &param1
                << "\n";
-    }
-
-    // update Kalman with the TRD measurement
-    if(chi2>10){ // RS
-      //    if(chi2>1e+10){ // TODO
+     }
+     
+     // update Kalman with the TRD measurement
+     if (chi2> fkRecoParam->GetChi2Cut()){ // MI parameterizad chi2 cut 03.05.2014
+       //       if(chi2>10){ // RS
+       //    if(chi2>1e+10){ // TODO
       t.SetErrStat(AliTRDtrackV1::kChi2, ily);
       if(debugLevel > 2){
         UChar_t status(t.GetStatusTRD());