]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCtrackerMI.cxx
Use symetric Mag field
[u/mrichter/AliRoot.git] / TPC / AliTPCtrackerMI.cxx
index a0f1c44f9ef321f1c6a2af6f01854c46595465f1..0b0f2adcc401d8377eb24b47c83aab609fcabf77 100644 (file)
 // The debug level -  different procedure produce tree for numerical debugging
 //                    To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
 //
+
+//
+// Adding systematic errors to the covariance:
+// 
+// The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
+// of the tracks (not to the clusters as they are dependent):
+// The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
+// The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/cm)
+// The default values are 0. 
+//
+// The sytematic errors are added to the covariance matrix in following places:
+//
+// 1. During fisrt itteration - AliTPCtrackerMI::FillESD
+// 2. Second iteration - 
+//      2.a ITS->TPC   - AliTPCtrackerMI::ReadSeeds 
+//      2.b TPC->TRD   - AliTPCtrackerMI::PropagateBack
+// 3. Third iteration  -
+//      3.a TRD->TPC   - AliTPCtrackerMI::ReadSeeds
+//      3.b TPC->ITS   - AliTPCtrackerMI::RefitInward
+//
 // There are several places in the code which can be numerically debuged
 // This code is keeped in order to enable code development and to check the calibration implementtion
 //
@@ -283,7 +303,7 @@ Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluste
   Double_t rdistance2  = rdistancey2+rdistancez2;
   //Int_t  accept =0;
   
-  if (AliTPCReconstructor::StreamLevel()>5) {
+  if (AliTPCReconstructor::StreamLevel()>5 && seed->GetNumberOfClusters()>20) {
     Float_t rmsy2 = seed->GetCurrentSigmaY2();
     Float_t rmsz2 = seed->GetCurrentSigmaZ2();
     Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
@@ -297,6 +317,8 @@ Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluste
     param.GetXYZ(gcl.GetMatrixArray());
     cluster->GetGlobalXYZ(gclf);
     gcl[0]=gclf[0];    gcl[1]=gclf[1];    gcl[2]=gclf[2];
+    
+    if (AliTPCReconstructor::StreamLevel()>0) {
     (*fDebugStreamer)<<"ErrParam"<<
       "Cl.="<<cluster<<
       "T.="<<&param<<
@@ -311,6 +333,7 @@ Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluste
       "rmsy2p30R="<<rmsy2p30R<<
       "rmsz2p30R="<<rmsz2p30R<<        
       "\n";
+    }
   }
   
   if (rdistance2>16) return 3;
@@ -384,7 +407,9 @@ AliTracker(),
     fYMax[i+nrowlow]      = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
   }
 
-  fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
+  if (AliTPCReconstructor::StreamLevel()>0) {
+    fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
+  }
 }
 //________________________________________________________________________
 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
@@ -448,6 +473,12 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr)
       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
       if (!pt) continue; 
       pt->UpdatePoints();
+      AddCovariance(pt);
+      if (AliTPCReconstructor::StreamLevel()>1) {
+       (*fDebugStreamer)<<"Track0"<<
+         "Tr0.="<<pt<<
+         "\n";       
+      }
       //      pt->PropagateTo(fParam->GetInnerRadiusLow());
       if (pt->GetKinkIndex(0)<=0){  //don't propagate daughter tracks 
        pt->PropagateTo(fParam->GetInnerRadiusLow());
@@ -1114,6 +1145,63 @@ Int_t  AliTPCtrackerMI::LoadClusters(TObjArray *arr)
   return 0;
 }
 
+Int_t  AliTPCtrackerMI::LoadClusters(TClonesArray *arr)
+{
+  //
+  // load clusters to the memory from one 
+  // TClonesArray
+  //
+  AliTPCclusterMI *clust=0;
+  Int_t count[72][96] = { {0} , {0} }; 
+
+  // loop over clusters
+  for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
+    clust = (AliTPCclusterMI*)arr->At(icl);
+    if(!clust) continue;
+    //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
+
+    // transform clusters
+    Transform(clust);
+
+    // count clusters per pad row
+    count[clust->GetDetector()][clust->GetRow()]++;
+  }
+
+  // insert clusters to sectors
+  for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
+    clust = (AliTPCclusterMI*)arr->At(icl);
+    if(!clust) continue;
+
+    Int_t sec = clust->GetDetector();
+    Int_t row = clust->GetRow();
+
+    // filter overlapping pad rows needed by HLT
+    if(sec<fkNIS*2) { //IROCs
+     if(row == 30) continue;
+    }
+    else { // OROCs
+      if(row == 27 || row == 76) continue;
+    }
+
+    Int_t left=0;
+    if (sec<fkNIS*2){
+      left = sec/fkNIS;
+      fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fParam);    
+    }
+    else{
+      left = (sec-fkNIS*2)/fkNOS;
+      fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fParam);
+    }
+  }
+
+  // Load functions must be called behind LoadCluster(TClonesArray*)
+  // needed by HLT
+  //LoadOuterSectors();
+  //LoadInnerSectors();
+
+  return 0;
+}
+
 
 Int_t  AliTPCtrackerMI::LoadClusters()
 {
@@ -1238,11 +1326,9 @@ void   AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
   Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
   Int_t i[1]={cluster->GetDetector()};
   transform->Transform(x,i,0,1);  
-  if (!AliTPCReconstructor::GetRecoParam()->GetBYMirror()){
-    if (cluster->GetDetector()%36>17){
-      x[1]*=-1;
-    }
-  }
+  //  if (cluster->GetDetector()%36>17){
+  //  x[1]*=-1;
+  //}
 
   //
   // in debug mode  check the transformation
@@ -2556,6 +2642,7 @@ Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
 
     seed->PropagateTo(fParam->GetInnerRadiusLow());
     seed->UpdatePoints();
+    AddCovariance(seed);
     MakeBitmaps(seed);
     AliESDtrack *esd=event->GetTrack(i);
     seed->CookdEdx(0.02,0.6);
@@ -2621,6 +2708,7 @@ Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
     if (!seed) continue;
     if (seed->GetKinkIndex(0)<0)  UpdateKinkQualityM(seed);  // update quality informations for kinks
     seed->UpdatePoints();
+    AddCovariance(seed);
     AliESDtrack *esd=event->GetTrack(i);
     seed->CookdEdx(0.02,0.6);
     CookLabel(seed,0.1); //For comparison only
@@ -2635,11 +2723,12 @@ Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
       ntracks++;
       Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
       // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number      
-      if (AliTPCReconstructor::StreamLevel()>1) {
+      if (AliTPCReconstructor::StreamLevel()>1 && esd) {
        (*fDebugStreamer)<<"Cback"<<
          "Tr0.="<<seed<<
+         "esd.="<<esd<<
          "EventNrInFile="<<eventnumber<<
-         "\n"; // patch 28 fev 06         
+         "\n";       
       }
     }
   }
@@ -2691,6 +2780,7 @@ void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
     //    AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
     AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
     seed->SetUniqueID(esd->GetID());
+    AddCovariance(seed);   //add systematic ucertainty
     for (Int_t ikink=0;ikink<3;ikink++) {
       Int_t index = esd->GetKinkIndex(ikink);
       seed->GetKinkIndexes()[ikink] = index;
@@ -4071,7 +4161,6 @@ void  AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_
     }
     if (ncl>0) xm[i]/=Float_t(ncl);
   }  
-  TTreeSRedirector &cstream = *fDebugStreamer;
   //
   for (Int_t i0=0;i0<nentries;i0++){
     AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
@@ -4132,6 +4221,8 @@ void  AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_
        }
       }
       //
+      if (AliTPCReconstructor::StreamLevel()>0) {
+      TTreeSRedirector &cstream = *fDebugStreamer;
       cstream<<"Multi"<<
        "iter="<<iter<<
        "lab0="<<lab0<<
@@ -4167,6 +4258,7 @@ void  AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_
        "r1="<<r1<<
        "rc1="<<rc1<<
        "\n";
+       }
     }
   }    
   delete [] helixes;
@@ -4253,7 +4345,6 @@ void  AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int
     }
     if (ncl>0) xm[i]/=Float_t(ncl);
   }  
-  TTreeSRedirector &cstream = *fDebugStreamer;
   //
   for (Int_t is0=0;is0<nentries;is0++){
     Int_t i0 = indexes[is0];
@@ -4335,39 +4426,43 @@ void  AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int
       //
       Int_t lab0=track0->GetLabel();
       Int_t lab1=track1->GetLabel();
-      cstream<<"Splitted"<<
-       "iter="<<iter<<
-       "lab0="<<lab0<<
-       "lab1="<<lab1<<   
-       "Tr0.="<<track0<<       // seed0
-       "Tr1.="<<track1<<       // seed1
-       "h0.="<<&helixes[i0]<<
-       "h1.="<<&helixes[i1]<<
-       //
-       "sum="<<sum<<           //the sum of rows with cl in both
-       "sum0="<<sum0<<           //the sum of rows with cl in 0 track
-       "sum1="<<sum1<<           //the sum of rows with cl in 1 track
-       "sums="<<sums<<         //the sum of shared clusters
-       "xm0="<<xm[i0]<<        // the center of track
-       "xm1="<<xm[i1]<<        // the x center of track
-       // General cut variables                   
-       "dfi="<<dfi<<           // distance in fi angle
-       "dtheta="<<dtheta<<     // distance int theta angle
-       //
-       //
-       "dist0="<<dist[4][0]<<     //distance x
-       "dist1="<<dist[4][1]<<     //distance y
-       "dist2="<<dist[4][2]<<     //distance z
-       "mdist0="<<mdist[0]<<   //distance x
-       "mdist1="<<mdist[1]<<   //distance y
-       "mdist2="<<mdist[2]<<   //distance z
-
-       "\n";
+      if( AliTPCReconstructor::StreamLevel()>5){
+      TTreeSRedirector &cstream = *fDebugStreamer;
+       cstream<<"Splitted"<<
+         "iter="<<iter<<
+         "lab0="<<lab0<<
+         "lab1="<<lab1<<   
+         "Tr0.="<<track0<<       // seed0
+         "Tr1.="<<track1<<       // seed1
+         "h0.="<<&helixes[i0]<<
+         "h1.="<<&helixes[i1]<<
+         //
+         "sum="<<sum<<           //the sum of rows with cl in both
+         "sum0="<<sum0<<           //the sum of rows with cl in 0 track
+         "sum1="<<sum1<<           //the sum of rows with cl in 1 track
+         "sums="<<sums<<         //the sum of shared clusters
+         "xm0="<<xm[i0]<<        // the center of track
+         "xm1="<<xm[i1]<<        // the x center of track
+         // General cut variables                   
+         "dfi="<<dfi<<           // distance in fi angle
+         "dtheta="<<dtheta<<     // distance int theta angle
+         //
+         //
+         "dist0="<<dist[4][0]<<     //distance x
+         "dist1="<<dist[4][1]<<     //distance y
+         "dist2="<<dist[4][2]<<     //distance z
+         "mdist0="<<mdist[0]<<   //distance x
+         "mdist1="<<mdist[1]<<   //distance y
+         "mdist2="<<mdist[2]<<   //distance z    
+         "\n";
+      }
       delete array->RemoveAt(i1);
     }
   }    
   delete [] helixes;
   delete [] xm;
+  delete [] quality;
+  delete [] indexes;
   AliInfo("Time for splitted tracks removal");
   timer.Print();
 }
@@ -4439,7 +4534,6 @@ void  AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_
   //
   // Find tracks
   //
-  TTreeSRedirector &cstream = *fDebugStreamer;
   //
   for (Int_t i0=0;i0<nentries;i0++){
     AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
@@ -4537,6 +4631,7 @@ void  AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_
          //debug stream to tune "fine" cuts      
          Int_t lab0=track0->GetLabel();
          Int_t lab1=track1->GetLabel();
+          TTreeSRedirector &cstream = *fDebugStreamer;
          cstream<<"Curling2"<<
            "iter="<<iter<<
            "lab0="<<lab0<<
@@ -4651,7 +4746,6 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
 
   //
   // Find circling track
-  TTreeSRedirector &cstream = *fDebugStreamer;
   //
   for (Int_t i0=0;i0<nentries;i0++){
     AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
@@ -4738,6 +4832,7 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
          //debug stream          
          Int_t lab0=track0->GetLabel();
          Int_t lab1=track1->GetLabel();
+          TTreeSRedirector &cstream = *fDebugStreamer;
          cstream<<"Curling"<<
            "lab0="<<lab0<<
            "lab1="<<lab1<<   
@@ -5342,7 +5437,6 @@ void  AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
   // 
   //
   // //  
-  TTreeSRedirector &cstream = *fDebugStreamer; 
   Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
   AliV0 vertex; 
   Double_t cradius0 = 10*10;
@@ -5357,6 +5451,7 @@ void  AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
     AliTPCseed * track0 = (AliTPCseed*)array->At(i);
     if (!track0) continue;
     if (AliTPCReconstructor::StreamLevel()>1){
+      TTreeSRedirector &cstream = *fDebugStreamer;
       cstream<<"Tracks"<<
        "Tr0.="<<track0<<
        "dca="<<dca[i]<<
@@ -5470,6 +5565,7 @@ void  AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
          Int_t lab1=track1->GetLabel();
          Char_t c0=track0->GetCircular();
          Char_t c1=track1->GetCircular();
+          TTreeSRedirector &cstream = *fDebugStreamer;
          cstream<<"V0"<<
          "Event="<<eventNr<<
          "vertex.="<<&vertex<<
@@ -5580,6 +5676,7 @@ void  AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
       if (AliTPCReconstructor::StreamLevel()>1) {
        Int_t lab0=track0->GetLabel();
        Int_t lab1=track1->GetLabel();
+        TTreeSRedirector &cstream = *fDebugStreamer;
        cstream<<"V02"<<
        "Event="<<eventNr<<
        "vertex.="<<v0<<        
@@ -6910,3 +7007,27 @@ void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
     }
   }
 }
+
+
+
+void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
+  //
+  // Adding systematic error
+  // !!!! the systematic error for element 4 is in 1/cm not in pt 
+
+  const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
+  Double_t covar[15];
+  for (Int_t i=0;i<15;i++) covar[i]=0;
+  // 0
+  // 1    2
+  // 3    4    5
+  // 6    7    8    9 
+  // 10   11   12   13   14
+  covar[0] = param[0]*param[0];
+  covar[2] = param[1]*param[1];
+  covar[5] = param[2]*param[2];
+  covar[9] = param[3]*param[3];
+  Double_t facC =  AliTracker::GetBz()*kB2C;
+  covar[14]= param[4]*param[4]*facC*facC;
+  seed->AddCovariance(covar);
+}