]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDseed.cxx
Make some calculations optional for HLT
[u/mrichter/AliRoot.git] / TRD / AliTRDseed.cxx
index 7d1129f1814592d1c3932893531b5096946926a4..6a77664081e268f8f3d1e31026c5350604859af9 100644 (file)
@@ -30,6 +30,7 @@
 #include "AliTRDcalibDB.h"
 #include "AliTRDcluster.h"
 #include "AliTRDtracker.h"
+#include "AliTRDtrackerV1.h"
 
 ClassImp(AliTRDseed)
 
@@ -126,14 +127,13 @@ AliTRDseed::AliTRDseed(const AliTRDseed &s)
 }
 
 //_____________________________________________________________________________
-void AliTRDseed::Copy(TObject &o) const 
+void AliTRDseed::Copy(TObject &o) const
 {
-  //
-  // Copy function
-  //
+       //printf("AliTRDseed::Copy()\n");
 
-  AliTRDseed &seed = (AliTRDseed &)o;
-  seed.fTilt = fTilt;
+       AliTRDseed &seed = (AliTRDseed &)o;
+  
+       seed.fTilt = fTilt;
   seed.fPadLength = fPadLength;
   seed.fX0 = fX0;
   seed.fSigmaY = fSigmaY;
@@ -150,8 +150,7 @@ void AliTRDseed::Copy(TObject &o) const
   seed.fCC = fCC;
   seed.fChi2 = fChi2;
   seed.fChi2Z = fChi2Z;
-
-  for (Int_t i = 0; i < knTimebins; i++) {
+       for (Int_t i = 0; i < knTimebins; i++) {
     seed.fX[i]        = fX[i];
     seed.fY[i]        = fY[i]; 
     seed.fZ[i]        = fZ[i]; 
@@ -218,14 +217,11 @@ void AliTRDseed::CookLabels()
   // Cook 2 labels for seed
   //
 
-  AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
-  Int_t nTimeBins = cal->GetNumberOfTimeBins();
-
   Int_t labels[200];
   Int_t out[200];
   Int_t nlab = 0;
 
-  for (Int_t i = 0; i < nTimeBins+1; i++) {
+  for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
     if (!fClusters[i]) continue;
     for (Int_t ilab = 0; ilab < 3; ilab++) {
       if (fClusters[i]->GetLabel(ilab) >= 0) {
@@ -251,10 +247,7 @@ void AliTRDseed::UseClusters()
   // Use clusters
   //
 
-  AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
-  Int_t nTimeBins = cal->GetNumberOfTimeBins();
-
-  for (Int_t i = 0; i < nTimeBins+1; i++) {
+  for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
     if (!fClusters[i]) continue;
     if (!(fClusters[i]->IsUsed())) fClusters[i]->Use();
   }
@@ -268,13 +261,25 @@ void AliTRDseed::Update()
   // Update the seed.
   //
 
+
+
+       // linear fit on the y direction with respect to the reference direction. 
+       // The residuals for each x (x = xc - x0) are deduced from:
+       // dy = y - yt             (1)
+       // the tilting correction is written :
+       // y = yc + h*(zc-zt)      (2)
+       // yt = y0+dy/dx*x         (3)
+       // zt = z0+dz/dx*x         (4)
+       // from (1),(2),(3) and (4)
+       // dy = yc - y0 - (dy/dx + h*dz/dx)*x + h*(zc-z0)
+       // the last term introduces the correction on y direction due to tilting pads. There are 2 ways to account for this:
+       // 1. use tilting correction for calculating the y
+       // 2. neglect tilting correction here and account for it in the error parametrization of the tracklet.
+
   const Float_t kRatio  = 0.8;
   const Int_t   kClmin  = 5;
   const Float_t kmaxtan = 2;
 
-  AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
-  Int_t nTimeBins = cal->GetNumberOfTimeBins();
-
   if (TMath::Abs(fYref[1]) > kmaxtan){
                //printf("Exit: Abs(fYref[1]) = %3.3f, kmaxtan = %3.3f\n", TMath::Abs(fYref[1]), kmaxtan);
                return;              // Track inclined too much
@@ -298,15 +303,18 @@ void AliTRDseed::Update()
   Int_t    zouts[2*knTimebins];       
   Float_t  allowedz[knTimebins];         // Allowed z for given time bin
   Float_t  yres[knTimebins];             // Residuals from reference
-  Float_t  anglecor = fTilt * fZref[1];  // Correction to the angle
+  //Float_t  anglecor = fTilt * fZref[1];  // Correction to the angle
   
   
   fN  = 0; 
   fN2 = 0;
-  for (Int_t i = 0; i < nTimeBins; i++) {
+  for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
     yres[i] = 10000.0;
     if (!fClusters[i]) continue;
-    yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i];   // Residual y
+    if(!fClusters[i]->IsInChamber()) continue;
+    // Residual y
+    //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] + fTilt*(fZ[i] - fZref[0]);
+    yres[i] = fY[i] - fTilt*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
     zints[fN] = Int_t(fZ[i]);
     fN++;    
   }
@@ -315,7 +323,7 @@ void AliTRDseed::Update()
                //printf("Exit fN < kClmin: fN = %d\n", fN);
                return; 
        }
-  Int_t nz = AliTRDtracker::Freq(fN,zints,zouts,kFALSE);
+  Int_t nz = AliTRDtracker::Freq(fN, zints, zouts, kFALSE);
   fZProb   = zouts[0];
   if (nz <= 1) zouts[3] = 0;
   if (zouts[1] + zouts[3] < kClmin) {
@@ -324,9 +332,7 @@ void AliTRDseed::Update()
        }
   
   // Z distance bigger than pad - length
-  if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) {
-    zouts[3]=0;           
-  }
+  if (TMath::Abs(zouts[0]-zouts[2]) > 12.0) zouts[3] = 0;
   
   Int_t  breaktime = -1;
   Bool_t mbefore   = kFALSE;
@@ -337,25 +343,25 @@ void AliTRDseed::Update()
 
     //
     // Find the break time allowing one chage on pad-rows
-    // with maximal numebr of accepted clusters
+    // with maximal number of accepted clusters
     //
     fNChange = 1;
-    for (Int_t i = 0; i < nTimeBins; i++) {
+    for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
       cumul[i][0] = counts[0];
       cumul[i][1] = counts[1];
       if (TMath::Abs(fZ[i]-zouts[0]) < 2) counts[0]++;
       if (TMath::Abs(fZ[i]-zouts[2]) < 2) counts[1]++;
     }
     Int_t  maxcount = 0;
-    for (Int_t i = 0; i < nTimeBins; i++) {
-      Int_t after  = cumul[nTimeBins][0] - cumul[i][0];
+    for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
+      Int_t after  = cumul[AliTRDtrackerV1::GetNTimeBins()][0] - cumul[i][0];
       Int_t before = cumul[i][1];
       if (after + before > maxcount) { 
        maxcount  = after + before; 
        breaktime = i;
        mbefore   = kFALSE;
       }
-      after  = cumul[nTimeBins-1][1] - cumul[i][1];
+      after  = cumul[AliTRDtrackerV1::GetNTimeBins()-1][1] - cumul[i][1];
       before = cumul[i][0];
       if (after + before > maxcount) { 
        maxcount  = after + before; 
@@ -368,18 +374,18 @@ void AliTRDseed::Update()
 
   }
 
-  for (Int_t i = 0; i < nTimeBins+1; i++) {
+  for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
     if (i >  breaktime) allowedz[i] =   mbefore  ? zouts[2] : zouts[0];
     if (i <= breaktime) allowedz[i] = (!mbefore) ? zouts[2] : zouts[0];
   }  
 
-  if (((allowedz[0] > allowedz[nTimeBins]) && (fZref[1] < 0)) || 
-      ((allowedz[0] < allowedz[nTimeBins]) && (fZref[1] > 0))) {
+  if (((allowedz[0] > allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] < 0)) ||
+      ((allowedz[0] < allowedz[AliTRDtrackerV1::GetNTimeBins()]) && (fZref[1] > 0))) {
     //
     // Tracklet z-direction not in correspondance with track z direction 
     //
     fNChange = 0;
-    for (Int_t i = 0; i < nTimeBins+1; i++) {
+    for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
       allowedz[i] = zouts[0];  // Only longest taken
     } 
   }
@@ -388,22 +394,26 @@ void AliTRDseed::Update()
     //
     // Cross pad -row tracklet  - take the step change into account
     //
-    for (Int_t i = 0; i < nTimeBins+1; i++) {
+    for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
       if (!fClusters[i]) continue; 
+      if(!fClusters[i]->IsInChamber()) continue;
       if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
-      yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i];   // Residual y
-      if (TMath::Abs(fZ[i] - fZProb) > 2) {
-       if (fZ[i] > fZProb) yres[i] += fTilt * fPadLength;
-       if (fZ[i] < fZProb) yres[i] -= fTilt * fPadLength;
-      }
+      // Residual y
+      //yres[i] = fY[i] - fYref[0] - (fYref[1] + anglecor) * fX[i] /*+ fTilt*(fZ[i] - fZref[0])*/;   
+      yres[i] = fY[i] - fTilt*(fZ[i] - (fZref[0] - fX[i]*fZref[1]));
+/*      if (TMath::Abs(fZ[i] - fZProb) > 2) {
+        if (fZ[i] > fZProb) yres[i] += fTilt * fPadLength;
+        if (fZ[i] < fZProb) yres[i] -= fTilt * fPadLength;
+      }*/
     }
   }
   
   Double_t yres2[knTimebins];
   Double_t mean;
   Double_t sigma;
-  for (Int_t i = 0; i < nTimeBins+1; i++) {
+  for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
     if (!fClusters[i]) continue;
+    if(!fClusters[i]->IsInChamber()) continue;
     if (TMath::Abs(fZ[i] - allowedz[i]) > 2) continue;
     yres2[fN2] = yres[i];
     fN2++;
@@ -432,12 +442,13 @@ void AliTRDseed::Update()
   fMeanz = 0;
   fMPads = 0;
 
-  for (Int_t i = 0; i < nTimeBins+1; i++) {
+  for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
 
     fUsable[i] = kFALSE;
     if (!fClusters[i]) continue;
-    if (TMath::Abs(fZ[i] - allowedz[i]) > 2)  continue;
-    if (TMath::Abs(yres[i] - mean) > 4.0 * sigma) continue;
+    if (!fClusters[i]->IsInChamber()) continue;
+    if (TMath::Abs(fZ[i] - allowedz[i]) > 2){fClusters[i] = 0x0; continue;}
+    if (TMath::Abs(yres[i] - mean) > 4.0 * sigma){fClusters[i] = 0x0;  continue;}
     fUsable[i] = kTRUE;
     fN2++;
     fMPads += fClusters[i]->GetNPads();
@@ -445,7 +456,10 @@ void AliTRDseed::Update()
     if (fClusters[i]->GetNPads() > 4) weight = 0.5;
     if (fClusters[i]->GetNPads() > 5) weight = 0.2;
    
+       
     Double_t x = fX[i];
+    //printf("x = %7.3f dy = %7.3f fit %7.3f\n", x, yres[i], fY[i]-yres[i]);
+    
     sumw   += weight; 
     sumwx  += x * weight; 
     sumwx2 += x*x * weight;
@@ -474,22 +488,26 @@ void AliTRDseed::Update()
   fYfitR[1]    = (sumw   * sumwxy - sumwx * sumwy)  / det;
   
   fSigmaY2 = 0;
-  for (Int_t i = 0; i < nTimeBins+1; i++) {    
+  for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins()+1; i++) {
     if (!fUsable[i]) continue;
     Float_t delta = yres[i] - fYfitR[0] - fYfitR[1] * fX[i];
     fSigmaY2 += delta*delta;
   }
   fSigmaY2 = TMath::Sqrt(fSigmaY2 / Float_t(fN2-2));
+       // TEMPORARY UNTIL covariance properly calculated
+       fSigmaY2 = TMath::Max(fSigmaY2, Float_t(.1));
   
   fZfitR[0]  = (sumwx2 * sumwz  - sumwx * sumwxz) / det;
   fZfitR[1]  = (sumw   * sumwxz - sumwx * sumwz)  / det;
   fZfit[0]   = (sumwx2 * sumwz  - sumwx * sumwxz) / det;
   fZfit[1]   = (sumw   * sumwxz - sumwx * sumwz)  / det;
-  fYfitR[0] += fYref[0] + correction;
-  fYfitR[1] += fYref[1];
+//   fYfitR[0] += fYref[0] + correction;
+//   fYfitR[1] += fYref[1];
   fYfit[0]   = fYfitR[0];
-  fYfit[1]   = fYfitR[1];
-    
+  fYfit[1]   = -fYfitR[1];
+
+       //printf("y0 = %7.3f tgy = %7.3f z0 = %7.3f tgz = %7.3f \n", fYfitR[0], fYfitR[1], fZfitR[0], fZfitR[1]);
+
   UpdateUsed();
 
 }
@@ -501,17 +519,11 @@ void AliTRDseed::UpdateUsed()
   // Update used seed
   //
 
-  AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
-  Int_t nTimeBins = cal->GetNumberOfTimeBins();
-
   fNUsed = 0;
-  for (Int_t i = 0; i < nTimeBins; i++) {
-    if (!fClusters[i]) {
-      continue;
-    }
-    if ((fClusters[i]->IsUsed())) {
-      fNUsed++;
-    }
+  for (Int_t i = 0; i < AliTRDtrackerV1::GetNTimeBins(); i++) {
+    if (!fClusters[i]) continue;
+               if(!fUsable[i]) continue;   
+    if ((fClusters[i]->IsUsed())) fNUsed++;
   }
 
 }
@@ -527,9 +539,6 @@ Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
   TLinearFitter fitterT2(4,"hyp4");  
   fitterT2.StoreData(kTRUE);
        
-  AliTRDcalibDB *cal = AliTRDcalibDB::Instance();
-  Int_t nTimeBins = cal->GetNumberOfTimeBins();
-
   Float_t xref2 = (cseed[2].fX0 + cseed[3].fX0) * 0.5; // Reference x0 for z
   
   Int_t npointsT = 0;
@@ -540,7 +549,7 @@ Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
     if (!cseed[iLayer].IsOK()) continue;
     Double_t tilt = cseed[iLayer].fTilt;
 
-    for (Int_t itime = 0; itime < nTimeBins+1; itime++) {
+    for (Int_t itime = 0; itime < AliTRDtrackerV1::GetNTimeBins()+1; itime++) {
 
       if (!cseed[iLayer].fUsable[itime]) continue;
       // x relative to the midle chamber
@@ -629,7 +638,7 @@ Float_t AliTRDseed::FitRiemanTilt(AliTRDseed * cseed, Bool_t terror)
     if (-params[2]*params[0] + params[1]*params[1] + 1 > 0) {
       Double_t rm1 = params[0] / TMath::Sqrt(-params[2]*params[0] + params[1]*params[1] + 1); 
       if (1.0/(rm1*rm1) - (x-x0) * (x-x0) > 0.0) {
-                               Double_t res = (x - x0) / TMath::Sqrt(1.0 / (rm1*rm1) - (x-x0)*(x-x0));
+                               Double_t res = (x - x0) / TMath::Sqrt((1./rm1-(x-x0))*(1./rm1+(x-x0)));
                                if (params[0] < 0) res *= -1.0;
                                dy = res;
       }