]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCtrackerMI.cxx
correcting baryon mass subtraction for visible energy in MC
[u/mrichter/AliRoot.git] / TPC / AliTPCtrackerMI.cxx
index ec9807fea4cb87b917b2cc02333c038db389b00e..6edb2eb50e2060d93a34f5c196444b6586231d2d 100644 (file)
@@ -297,7 +297,7 @@ Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluste
   Double_t sdistancey2 = sy2+seed->GetSigmaY2();
   Double_t sdistancez2 = sz2+seed->GetSigmaZ2();
   Double_t dy=seed->GetCurrentCluster()->GetY()-yt;
-  Double_t dz=seed->GetCurrentCluster()->GetY()-zt;
+  Double_t dz=seed->GetCurrentCluster()->GetZ()-zt;
   Double_t rdistancey2 = (seed->GetCurrentCluster()->GetY()-yt)*
     (seed->GetCurrentCluster()->GetY()-yt)/sdistancey2;
   Double_t rdistancez2 = (seed->GetCurrentCluster()->GetZ()-zt)*
@@ -339,6 +339,10 @@ Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluste
       "rmsz2p30="<<rmsz2p30<<  
       "rmsy2p30R="<<rmsy2p30R<<
       "rmsz2p30R="<<rmsz2p30R<<        
+      // normalize distance - 
+      "rdisty="<<rdistancey2<<
+      "rdistz="<<rdistancez2<<
+      "rdist="<<rdistance2<< //       
       "\n";
     }
   }
@@ -606,7 +610,7 @@ void AliTPCtrackerMI::FillESD(const TObjArray* arr)
     }
     // << account for suppressed tracks in the kink indices (RS)
   }
-  printf("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks());
+  AliInfo(Form("Number of filled ESDs-\t%d\n",fEvent->GetNumberOfTracks()));
 }
 
 
@@ -980,7 +984,7 @@ Double_t AliTPCtrackerMI::F1(Double_t x1,Double_t y1,
   y2 -=y1;
   //  
   Double_t det = x3*y2-x2*y3;
-  if (TMath::AreEqualRel(det,0,1e-16)){
+  if (TMath::Abs(det)<1e-10){
     return 100;
   }
   //
@@ -1006,7 +1010,7 @@ Double_t AliTPCtrackerMI::F2(Double_t x1,Double_t y1,
   y2 -=y1;
   //  
   Double_t det = x3*y2-x2*y3;
-  if (TMath::AreEqualRel(det,0,1e-16)) {
+  if (TMath::Abs(det)<1e-10) {
     return 100;
   }
   //
@@ -1123,7 +1127,7 @@ Int_t  AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
 {
   //
   // load clusters to the memory
-  AliTPCClustersRow *clrow = 0x0; 
+  AliTPCClustersRow *clrow = new AliTPCClustersRow("AliTPCclusterMI");
   Int_t lower   = arr->LowerBound();
   Int_t entries = arr->GetEntriesFast();
   
@@ -1163,6 +1167,7 @@ Int_t  AliTPCtrackerMI::LoadClusters(const TObjArray *arr)
       for (Int_t j=0;j<tpcrow->GetN2();++j) 
        tpcrow->SetCluster2(j,*(AliTPCclusterMI*)(clrow->GetArray()->At(j)));
     }
+    clrow->GetArray()->Clear("C");
   }
   //
   delete clrow;
@@ -1233,10 +1238,7 @@ Int_t  AliTPCtrackerMI::LoadClusters()
 {
   //
   // load clusters to the memory
-  AliTPCClustersRow *clrow= new AliTPCClustersRow;
-  clrow->SetClass("AliTPCclusterMI");
-  clrow->SetArray(0);
-  clrow->GetArray()->ExpandCreateFast(10000);
+  AliTPCClustersRow *clrow= new AliTPCClustersRow("AliTPCclusterMI");
   //
   //  TTree * tree = fClustersArray.GetTree();
 
@@ -2310,6 +2312,88 @@ void AliTPCtrackerMI::RemoveUsed2(TObjArray * arr, Float_t factor1,  Float_t fac
   delete []indexes;
 }
 
+void AliTPCtrackerMI::DumpClusters(Int_t iter, TObjArray *trackArray) 
+{
+  //
+  // Dump clusters after reco
+  // signed and unsigned cluster can be visualized   
+  // 1. Unsign all cluster
+  // 2. Sign all used clusters
+  // 3. Dump clusters
+  UnsignClusters();
+  Int_t nseed = trackArray->GetEntries();
+  for (Int_t i=0; i<nseed; i++){
+    AliTPCseed *pt=(AliTPCseed*)trackArray->UncheckedAt(i);    
+    if (!pt) {
+      continue;
+    }    
+    Bool_t isKink=pt->GetKinkIndex(0)!=0;
+    for (Int_t j=0; j<160; ++j) {
+      Int_t index=pt->GetClusterIndex2(j);
+      if (index<0) continue;
+      AliTPCclusterMI *c= pt->GetClusterPointer(j);
+      if (!c) continue;
+      if (isKink) c->Use(100);   // kink
+      c->Use(10);                // by default usage 10
+    }
+  }
+  //
+
+  for (Int_t sec=0;sec<fkNIS;sec++){
+    for (Int_t row=0;row<fInnerSec->GetNRows();row++){
+      AliTPCclusterMI *cl = fInnerSec[sec][row].GetClusters1();
+      for (Int_t icl =0;icl< fInnerSec[sec][row].GetN1();icl++){    
+       Float_t gx[3];  cl[icl].GetGlobalXYZ(gx);
+       (*fDebugStreamer)<<"clDump"<< 
+         "iter="<<iter<<
+         "cl.="<<&cl[icl]<<      
+         "gx0="<<gx[0]<<
+         "gx1="<<gx[1]<<
+         "gx2="<<gx[2]<<
+         "\n";
+      }
+      cl = fInnerSec[sec][row].GetClusters2();
+      for (Int_t icl =0;icl< fInnerSec[sec][row].GetN2();icl++){
+       Float_t gx[3];  cl[icl].GetGlobalXYZ(gx);
+       (*fDebugStreamer)<<"clDump"<< 
+         "iter="<<iter<<
+         "cl.="<<&cl[icl]<<
+         "gx0="<<gx[0]<<
+         "gx1="<<gx[1]<<
+         "gx2="<<gx[2]<<
+         "\n";
+      }
+    }
+  }
+  
+  for (Int_t sec=0;sec<fkNOS;sec++){
+    for (Int_t row=0;row<fOuterSec->GetNRows();row++){
+      AliTPCclusterMI *cl = fOuterSec[sec][row].GetClusters1();
+      for (Int_t icl =0;icl< fOuterSec[sec][row].GetN1();icl++){
+       Float_t gx[3];  cl[icl].GetGlobalXYZ(gx);
+       (*fDebugStreamer)<<"clDump"<< 
+         "iter="<<iter<<
+         "cl.="<<&cl[icl]<<
+         "gx0="<<gx[0]<<
+         "gx1="<<gx[1]<<
+         "gx2="<<gx[2]<<
+         "\n";      
+      }
+      cl = fOuterSec[sec][row].GetClusters2();
+      for (Int_t icl =0;icl< fOuterSec[sec][row].GetN2();icl++){
+       Float_t gx[3];  cl[icl].GetGlobalXYZ(gx);
+       (*fDebugStreamer)<<"clDump"<< 
+         "iter="<<iter<<
+         "cl.="<<&cl[icl]<<
+         "gx0="<<gx[0]<<
+         "gx1="<<gx[1]<<
+         "gx2="<<gx[2]<<
+         "\n";      
+      }
+    }
+  }
+  
+}
 void AliTPCtrackerMI::UnsignClusters() 
 {
   //
@@ -2656,6 +2740,7 @@ Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
     }
   }
   //FindKinks(fSeeds,event);
+  if (AliTPCReconstructor::StreamLevel()>3)  DumpClusters(2,fSeeds);
   Info("RefitInward","Number of refitted tracks %d",ntracks);
   return 0;
 }
@@ -2726,6 +2811,7 @@ Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
       }
     }
   }
+  if (AliTPCReconstructor::StreamLevel()>3)  DumpClusters(1,fSeeds);
   //FindKinks(fSeeds,event);
   Info("PropagateBack","Number of back propagated tracks %d",ntracks);
   fEvent =0;
@@ -3423,7 +3509,7 @@ void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,
   }
   
   if (fDebug>3){
-    Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
+    Info("MakeSeeds5","\nSeeding statiistic:\t%d\t%d\t%d\t%d\t%d\t%d\t%d",nin0,nin1,nin2,nin,nout1,nout2,nout3);
   }
   delete seed;
 }
@@ -4295,6 +4381,7 @@ void  AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int
   //0.  Sort tracks according quality
   //1.  Propagate the ext. param to reference radius
   Int_t nseed = array->GetEntriesFast();  
+  if (nseed<=0) return;
   Float_t * quality = new Float_t[nseed];
   Int_t   * indexes = new Int_t[nseed];
   for (Int_t i=0; i<nseed; i++) {
@@ -4410,6 +4497,9 @@ void  AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int
   // 4. Delete temporary array
   //
   delete [] params; 
+  delete [] quality;
+  delete [] indexes;
+
 }
 
 
@@ -5327,7 +5417,7 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
   kinks->Delete();
   delete kinks;
 
-  printf("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall);
+  AliInfo(Form("Ncandidates=\t%d\t%d\t%d\t%d\n",esd->GetNumberOfKinks(),ncandidates,ntracks,nall));
   timer.Print();
 }
 
@@ -6090,6 +6180,7 @@ Int_t AliTPCtrackerMI::Clusters2Tracks (AliESDEvent *const esd)
   Clusters2Tracks();
   if (!fSeeds) return 1;
   FillESD(fSeeds);
+  if (AliTPCReconstructor::StreamLevel()>3)  DumpClusters(0,fSeeds);
   return 0;
   //
 }
@@ -6305,6 +6396,9 @@ TObjArray * AliTPCtrackerMI::Tracking()
 
   TObjArray * seeds = new TObjArray;
   TObjArray * arr=0;
+  Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
+  Int_t gapPrim = AliTPCReconstructor::GetRecoParam()->GetSeedGapPrim();
+  Int_t gapSec = AliTPCReconstructor::GetRecoParam()->GetSeedGapSec();
   
   Int_t gap =20;
   Float_t cuts[4];
@@ -6318,7 +6412,7 @@ TObjArray * AliTPCtrackerMI::Tracking()
   //  
   //find primaries  
   cuts[0]=0.0066;
-  for (Int_t delta = 0; delta<18; delta+=6){
+  for (Int_t delta = 0; delta<18; delta+=gapPrim){
     //
     cuts[0]=0.0070;
     cuts[1] = 1.5;
@@ -6343,7 +6437,7 @@ TObjArray * AliTPCtrackerMI::Tracking()
 
   //find primaries  
   cuts[0]=0.0077;
-  for (Int_t delta = 20; delta<120; delta+=10){
+  for (Int_t delta = 20; delta<120; delta+=gapPrim){
     //
     // seed high pt tracks
     cuts[0]=0.0060;
@@ -6396,9 +6490,22 @@ TObjArray * AliTPCtrackerMI::Tracking()
   SumTracks(seeds,arr);   
   SignClusters(seeds,fnumber,fdensity);   
   //
+  arr = Tracking(4,nup-5,nup-5-gap,cuts,-1);
+  SumTracks(seeds,arr);   
+  SignClusters(seeds,fnumber,fdensity);   
+  //
+  arr = Tracking(4,nup-7,nup-7-gap,cuts,-1);
+  SumTracks(seeds,arr);   
+  SignClusters(seeds,fnumber,fdensity);   
+  //
+  //
+  arr = Tracking(4,nup-9,nup-9-gap,cuts,-1);
+  SumTracks(seeds,arr);   
+  SignClusters(seeds,fnumber,fdensity);   
+  //
 
 
-  for (Int_t delta = 3; delta<30; delta+=5){
+  for (Int_t delta = 9; delta<30; delta+=gapSec){
     //
     cuts[0] = 0.3;
     cuts[1] = 1.5;
@@ -6421,10 +6528,9 @@ TObjArray * AliTPCtrackerMI::Tracking()
   fdensity = 2.;
   cuts[0]=0.0080;
 
-  Int_t fLastSeedRowSec=AliTPCReconstructor::GetRecoParam()->GetLastSeedRowSec();
 
   // find secondaries
-  for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=10){
+  for (Int_t delta = 30; delta<fLastSeedRowSec; delta+=gapSec){
     //
     cuts[0] = 0.3;
     cuts[1] = 3.5;
@@ -6855,8 +6961,8 @@ void AliTPCtrackerMI::CookLabel(AliKalmanTrack *tk, Float_t wrong) const {
     if (TMath::Abs(c->GetLabel(1)) == lab ||
         TMath::Abs(c->GetLabel(2)) == lab ) max++;
   }
-
-  if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
+  if (noc<=0) { lab=-1; return;}
+  if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
 
   else {
      Int_t tail=Int_t(0.10*noc);
@@ -6939,8 +7045,8 @@ Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *const t, Float_t wrong,Int_t first,
     if (TMath::Abs(c->GetLabel(1)) == lab ||
         TMath::Abs(c->GetLabel(2)) == lab ) max++;
   }
-
-  if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
+  if (noc<=0) { lab=-1; return -1;}
+  if ((1.- Float_t(max)/(noc)) > wrong) lab=-lab;
 
   else {
      Int_t tail=Int_t(0.10*noc);
@@ -7030,6 +7136,7 @@ void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
   // !!!! the systematic error for element 4 is in 1/cm not in pt 
 
   const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
+  Double_t *covarIn= (Double_t*)seed->GetCovariance();
   Double_t covar[15];
   for (Int_t i=0;i<15;i++) covar[i]=0;
   // 0
@@ -7043,5 +7150,20 @@ void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
   covar[9] = param[3]*param[3];
   Double_t facC =  AliTracker::GetBz()*kB2C;
   covar[14]= param[4]*param[4]*facC*facC;
+  //
+  covar[1]=TMath::Sqrt((covar[0]*covar[2]))*covarIn[1]/TMath::Sqrt((covarIn[0]*covarIn[2]));
+  //
+  covar[3]=TMath::Sqrt((covar[0]*covar[5]))*covarIn[3]/TMath::Sqrt((covarIn[0]*covarIn[5]));
+  covar[4]=TMath::Sqrt((covar[2]*covar[5]))*covarIn[4]/TMath::Sqrt((covarIn[2]*covarIn[5]));
+  //
+  covar[6]=TMath::Sqrt((covar[0]*covar[9]))*covarIn[6]/TMath::Sqrt((covarIn[0]*covarIn[9]));
+  covar[7]=TMath::Sqrt((covar[2]*covar[9]))*covarIn[7]/TMath::Sqrt((covarIn[2]*covarIn[9]));
+  covar[8]=TMath::Sqrt((covar[5]*covar[9]))*covarIn[8]/TMath::Sqrt((covarIn[5]*covarIn[9]));
+  //
+  covar[10]=TMath::Sqrt((covar[0]*covar[14]))*covarIn[10]/TMath::Sqrt((covarIn[0]*covarIn[14]));
+  covar[11]=TMath::Sqrt((covar[2]*covar[14]))*covarIn[11]/TMath::Sqrt((covarIn[2]*covarIn[14]));
+  covar[12]=TMath::Sqrt((covar[5]*covar[14]))*covarIn[12]/TMath::Sqrt((covarIn[5]*covarIn[14]));
+  covar[13]=TMath::Sqrt((covar[9]*covar[14]))*covarIn[13]/TMath::Sqrt((covarIn[9]*covarIn[14]));
+  //
   seed->AddCovariance(covar);
 }