]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDtracker.cxx
Code for unified TPC/TRD tracking (S.Radomski)
[u/mrichter/AliRoot.git] / TRD / AliTRDtracker.cxx
index 0287f0678d558def0e396ffd6547e8137fa37ed4..4b666449257e65a443282c19ace600f98c9b2e88 100644 (file)
@@ -15,6 +15,9 @@
                                                       
 /*
 $Log$
+Revision 1.25  2003/03/19 17:14:11  hristov
+Load/UnloadClusters added to the base class and the derived classes changed correspondingly. Possibility to give 2 input files for ITS and TPC tracks in PropagateBack. TRD tracker uses fEventN from the base class (T.Kuhr)
+
 Revision 1.24  2003/02/19 09:02:28  hristov
 Track time measurement (S.Radomski)
 
@@ -126,9 +129,12 @@ ClassImp(AliTRDtracker)
 
   const  Double_t    AliTRDtracker::fMaxChi2            = 12.; 
 
+const Int_t AliTRDtracker::kFirstPlane = 5;
+const Int_t AliTRDtracker::kLastPlane = 17;
+
 
 //____________________________________________________________________
-AliTRDtracker::AliTRDtracker(const TFile *geomfile)
+AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker()
 {
   // 
   //  Main constructor
@@ -213,6 +219,15 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile)
 
   fVocal = kFALSE;
 
+
+  // Barrel Tracks [SR, 03.04.2003]
+
+  fBarrelFile = 0;
+  fBarrelTree = 0;
+  fBarrelArray = 0;
+  fBarrelTrack = 0;
+
+  savedir->cd();
 }   
 
 //___________________________________________________________________
@@ -229,6 +244,106 @@ AliTRDtracker::~AliTRDtracker()
   }
 }   
 
+//_____________________________________________________________________
+
+void AliTRDtracker::SetBarrelTree(const char *mode) {
+  //
+  //
+  //
+
+  if (!IsStoringBarrel()) return;
+
+  TDirectory *sav = gDirectory;
+  if (!fBarrelFile) fBarrelFile = new TFile("AliBarrelTracks.root", "UPDATE");
+
+  char buff[40];
+  sprintf(buff,  "BarrelTRD_%d_%s", GetEventNumber(), mode);
+
+  fBarrelFile->cd();
+  fBarrelTree = new TTree(buff, "Barrel TPC tracks");
+  
+  Int_t nRefs = kLastPlane - kFirstPlane + 1;
+
+  if (!fBarrelArray) fBarrelArray = new TClonesArray("AliBarrelTrack", nRefs);
+  for(Int_t i=0; i<nRefs; i++) new((*fBarrelArray)[i]) AliBarrelTrack();
+  
+  fBarrelTree->Branch("tracks", &fBarrelArray);
+  sav->cd();
+}
+  
+//_____________________________________________________________________
+
+void AliTRDtracker::StoreBarrelTrack(AliTRDtrack *ps, Int_t refPlane, Int_t isIn) {
+  //
+  //
+  //
+  
+  if (!IsStoringBarrel()) return;
+  
+  static Int_t nClusters;
+  static Int_t nWrong;
+  static Double_t chi2;
+  static Int_t index;
+  static Bool_t wasLast = kTRUE;
+  
+  Int_t newClusters, newWrong;
+  Double_t newChi2;
+  
+  if (wasLast) {   
+    fBarrelArray->Clear();
+    nClusters = nWrong = 0;
+    chi2 = 0.0;
+    index = 0;
+    wasLast = kFALSE;
+  }
+  
+  fBarrelTrack = (AliBarrelTrack*)(*fBarrelArray)[index++];
+  ps->GetBarrelTrack(fBarrelTrack);
+  
+  newClusters = ps->GetNumberOfClusters() - nClusters; 
+  newWrong = ps->GetNWrong() - nWrong;
+  newChi2 = ps->GetChi2() - chi2;
+  
+  nClusters =  ps->GetNumberOfClusters();
+  nWrong = ps->GetNWrong();
+  chi2 = ps->GetChi2();  
+
+  if (refPlane != kLastPlane) {
+    fBarrelTrack->SetNClusters(newClusters, newChi2);
+    fBarrelTrack->SetNWrongClusters(newWrong);
+  } else {
+    wasLast = kTRUE;
+  } 
+
+  fBarrelTrack->SetRefPlane(refPlane, isIn);
+}
+
+//_____________________________________________________________________
+
+Bool_t AliTRDtracker::AdjustSector(AliTRDtrack *track) {
+  //
+  // Rotates the track when necessary
+  //
+
+  Double_t alpha = AliTRDgeometry::GetAlpha(); 
+  Double_t y = track->GetY();
+  Double_t ymax = track->GetX()*TMath::Tan(0.5*alpha);
+
+  Int_t ns = AliTRDgeometry::kNsect;
+  //Int_t s=Int_t(track->GetAlpha()/alpha)%ns; 
+
+  if (y > ymax) {
+    //s = (s+1) % ns;
+    if (!track->Rotate(alpha)) return kFALSE;
+  } else if (y <-ymax) {
+    //s = (s-1+ns) % ns;                           
+    if (!track->Rotate(-alpha)) return kFALSE;   
+  } 
+
+  return kTRUE;
+}
+
 //_____________________________________________________________________
 inline Double_t f1trd(Double_t x1,Double_t y1,
                       Double_t x2,Double_t y2,
@@ -530,12 +645,14 @@ Int_t AliTRDtracker::PropagateBack(const TFile *inp, TFile *out) {
   AliTRDtrack *otrack_trd=0;
   trdTree.Branch("tracks","AliTRDtrack",&otrack_trd,32000,0);   
      
-  Int_t found=0;  
+  if (IsStoringBarrel()) SetBarrelTree("back");
+  out->cd();
 
+  Int_t found=0;  
   Int_t nseed=fSeeds->GetEntriesFast();
 
   //  Float_t foundMin = fMinClustersInTrack * fTimeBinsPerPlane * fGeom->Nplan(); 
-  Float_t foundMin = 0;
+  Float_t foundMin = 40;
 
   Int_t outermost_tb  = fTrSec[0]->GetOuterTimeBin();
 
@@ -543,6 +660,12 @@ Int_t AliTRDtracker::PropagateBack(const TFile *inp, TFile *out) {
 
     AliTRDtrack *ps=(AliTRDtrack*)fSeeds->UncheckedAt(i), &s=*ps;
     Int_t expectedClr = FollowBackProlongation(s);
+
+    if (IsStoringBarrel()) {
+      StoreBarrelTrack(ps, kLastPlane, kTrackBack);
+      fBarrelTree->Fill();        
+    }
+
     Int_t foundClr = s.GetNumberOfClusters();
     Int_t last_tb = fTrSec[0]->GetLayerNumber(s.GetX());
 
@@ -607,12 +730,20 @@ Int_t AliTRDtracker::PropagateBack(const TFile *inp, TFile *out) {
     }
   }
   
+  
+  out->cd();
   tofTree1.Write(); 
   tofTree2.Write(); 
   phosTree.Write(); 
   richTree.Write(); 
   trdTree.Write(); 
 
+  if (IsStoringBarrel()) { // [SR, 03.04.2003]
+    fBarrelFile->cd();
+    fBarrelTree->Write();
+    fBarrelFile->Flush();
+  }
+
   savedir->cd();  
   cerr<<"Number of seeds: "<<nseed<<endl;  
   cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
@@ -646,12 +777,9 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf)
   Int_t try_again=fMaxGap;
 
   Double_t alpha=t.GetAlpha();
-
-  if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
-  if (alpha < 0.            ) alpha += 2.*TMath::Pi();
+  TVector2::Phi_0_2pi(alpha);
 
   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;  
-
   Double_t rad_length, rho, x, dx, y, ymax, z;
 
   Int_t expectedNumberOfClusters = 0;
@@ -869,22 +997,15 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
   Float_t  wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
 
   Int_t trackIndex = t.GetLabel();  
-
-  Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);     
-
   Int_t try_again=fMaxGap;
 
   Double_t alpha=t.GetAlpha();
+  TVector2::Phi_0_2pi(alpha);
 
-  if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
-  if (alpha < 0.            ) alpha += 2.*TMath::Pi();
-
-  Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;  
+  Int_t s;
 
   Int_t outerTB = fTrSec[0]->GetOuterTimeBin();
-
-  Double_t rad_length, rho, x, dx, y, ymax, z;
-
+  Double_t rad_length, rho, x, dx, y, ymax = 0, z;
   Bool_t lookForCluster;
 
   Int_t expectedNumberOfClusters = 0;
@@ -892,55 +1013,63 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
 
   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
 
+  Int_t nRefPlane = kFirstPlane;
+  Bool_t isNewLayer = kFALSE; 
 
-  for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr<outerTB; nr++) { 
+  Double_t chi2;
+  Double_t minDY;
 
-    y = t.GetY(); z = t.GetZ();
+  for (Int_t nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr<outerTB+1; nr++) { 
+    
+    y = t.GetY(); 
+    z = t.GetZ();
 
     // first propagate to the outer surface of the current time bin 
 
+    s = t.GetSector();
     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
-    x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; y = t.GetY(); z = t.GetZ();
+    x = fTrSec[s]->GetLayer(nr)->GetX()+dx/2; 
+    y = t.GetY(); 
+    z = t.GetZ();
 
     if(!t.PropagateTo(x,rad_length,rho)) break;
+    if (!AdjustSector(&t)) break;
+    s = t.GetSector();
+    if (!t.PropagateTo(x,rad_length,rho)) break;
 
     y = t.GetY();
-    ymax = x*TMath::Tan(0.5*alpha);
+    z = t.GetZ();
 
-    if (y > ymax) {
-      s = (s+1) % ns;
-      if (!t.Rotate(alpha)) break;
-      if(!t.PropagateTo(x,rad_length,rho)) break;
-    } else if (y <-ymax) {
-      s = (s-1+ns) % ns;                           
-      if (!t.Rotate(-alpha)) break;   
-      if(!t.PropagateTo(x,rad_length,rho)) break;
-    } 
-    y = t.GetY(); z = t.GetZ();
+    // Barrel Tracks [SR, 04.04.2003]
 
-    // now propagate to the middle plane of the next time bin 
-    fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
+    s = t.GetSector();
+    if (fTrSec[s]->GetLayer(nr)->IsSensitive() != 
+        fTrSec[s]->GetLayer(nr+1)->IsSensitive() ) {
 
-    x = fTrSec[s]->GetLayer(nr+1)->GetX(); y = t.GetY(); z = t.GetZ();
+      if (IsStoringBarrel()) StoreBarrelTrack(&t, nRefPlane++, kTrackBack);
+    }
 
-    if(!t.PropagateTo(x,rad_length,rho)) break;
+    if (fTrSec[s]->GetLayer(nr-1)->IsSensitive() && 
+          ! fTrSec[s]->GetLayer(nr)->IsSensitive()) {
+      isNewLayer = kTRUE;
+    } else {isNewLayer = kFALSE;}
 
     y = t.GetY();
+    z = t.GetZ();
 
-    ymax = x*TMath::Tan(0.5*alpha);
-
-    if(fVocal) printf("nr+1=%d, x %f, z %f, y %f, ymax %f\n",nr+1,x,z,y,ymax);
+    // now propagate to the middle plane of the next time bin 
+    fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
 
-    if (y > ymax) {
-      s = (s+1) % ns;
-      if (!t.Rotate(alpha)) break;
+    x = fTrSec[s]->GetLayer(nr+1)->GetX(); 
       if(!t.PropagateTo(x,rad_length,rho)) break;
-    } else if (y <-ymax) {
-      s = (s-1+ns) % ns;              
-      if (!t.Rotate(-alpha)) break;   
+    if (!AdjustSector(&t)) break;
+    s = t.GetSector();
       if(!t.PropagateTo(x,rad_length,rho)) break;
-    } 
 
+    y = t.GetY();
+    z = t.GetZ();
+
+    if(fVocal) printf("nr+1=%d, x %f, z %f, y %f, ymax %f\n",nr+1,x,z,y,ymax);
     //    printf("label %d, pl %d, lookForCluster %d \n",
     //     trackIndex, nr+1, lookForCluster);
 
@@ -983,6 +1112,17 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
 
       Double_t max_chi2=fMaxChi2;
 
+      if (isNewLayer) { 
+        road = 3 * road;
+        //sz2 = 3 * sz2;
+        max_chi2 = 10 * fMaxChi2;
+      }
+      
+      if (nRefPlane == kFirstPlane) max_chi2 = 20 * fMaxChi2; 
+      if (nRefPlane == kFirstPlane+2) max_chi2 = 15 * fMaxChi2;
+      if (t.GetNRotate() > 0) max_chi2 = 3 * max_chi2;
+      
+
       wYclosest = 12345678;
       wYcorrect = 12345678;
       wZclosest = 12345678;
@@ -991,14 +1131,15 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
 
       // Find the closest correct cluster for debugging purposes
       if (time_bin) {
-        Float_t minDY = 1000000;
+        minDY = 1000000;
         for (Int_t i=0; i<time_bin; i++) {
           AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
           if((c->GetLabel(0) != trackIndex) &&
              (c->GetLabel(1) != trackIndex) &&
              (c->GetLabel(2) != trackIndex)) continue;
           if(TMath::Abs(c->GetY() - y) > minDY) continue;
-          minDY = TMath::Abs(c->GetY() - y);
+          //minDY = TMath::Abs(c->GetY() - y);
+          minDY = c->GetY() - y;
           wYcorrect = c->GetY();
           wZcorrect = c->GetZ();
 
@@ -1014,16 +1155,21 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
         for (Int_t i=time_bin.Find(y-road); i<time_bin; i++) {
           AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
           if (c->GetY() > y+road) break;
-         //          if (c->IsUsed() > 0) continue;
+          if (c->IsUsed() > 0) continue;
           if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
 
           Double_t h01 = GetTiltFactor(c);
-          Double_t chi2=t.GetPredictedChi2(c,h01);
+          chi2=t.GetPredictedChi2(c,h01);
           
           if (chi2 > max_chi2) continue;
           max_chi2=chi2;
           cl=c;
           index=time_bin.GetIndex(i);
+
+          //check is correct
+          if((c->GetLabel(0) != trackIndex) &&
+             (c->GetLabel(1) != trackIndex) &&
+             (c->GetLabel(2) != trackIndex)) t.AddNWrong();
         }               
        
         if(!cl) {
@@ -1032,11 +1178,11 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
             AliTRDcluster* c=(AliTRDcluster*)(time_bin[i]);
             
             if (c->GetY() > y+road) break;
-           //            if (c->IsUsed() > 0) continue;
+            if (c->IsUsed() > 0) continue;
             if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue;
             
             Double_t h01 = GetTiltFactor(c);
-            Double_t chi2=t.GetPredictedChi2(c,h01);
+            chi2=t.GetPredictedChi2(c,h01);
             
             if (chi2 > max_chi2) continue;
             max_chi2=chi2;
@@ -1059,8 +1205,15 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
         else {
           if (try_again==0) break; 
           try_again--;
+          
+          //if (minDY < 1000000 && isNewLayer) 
+            //cout << "\t" << nRefPlane << "\t" << "\t" << t.GetNRotate() <<  "\t" << 
+            //  road << "\t" << minDY << "\t" << chi2 << "\t" << wChi2 << "\t" << max_chi2 << endl;
+                                                                     
         }
 
+        isNewLayer = kFALSE;
+
         /*
         if((((Int_t) wTB)%15 == 0) || (((Int_t) wTB)%15 == 14)) {
           
@@ -1090,6 +1243,8 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
     }  
   }
   return expectedNumberOfClusters;
+
+
 }         
 
 //___________________________________________________________________
@@ -1169,9 +1324,7 @@ Int_t AliTRDtracker::PropagateToTPC(AliTRDtrack& t)
   Int_t ns=Int_t(2*TMath::Pi()/AliTRDgeometry::GetAlpha()+0.5);     
 
   Double_t alpha=t.GetAlpha();
-
-  if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
-  if (alpha < 0.            ) alpha += 2.*TMath::Pi();
+  TVector2::Phi_0_2pi(alpha);
 
   Int_t s=Int_t(alpha/AliTRDgeometry::GetAlpha())%AliTRDgeometry::kNsect;  
 
@@ -1181,49 +1334,35 @@ Int_t AliTRDtracker::PropagateToTPC(AliTRDtrack& t)
   x = t.GetX();
 
   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
-
   Int_t plTPC = fTrSec[0]->GetLayerNumber(246.055);
 
   for (Int_t nr=fTrSec[0]->GetLayerNumber(x); nr>plTPC; nr--) { 
 
-    y = t.GetY(); z = t.GetZ();
+    y = t.GetY(); 
+    z = t.GetZ();
 
     // first propagate to the outer surface of the current time bin 
     fTrSec[s]->GetLayer(nr)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
-    x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; y = t.GetY(); z = t.GetZ();
+    x = fTrSec[s]->GetLayer(nr)->GetX()-dx/2; 
+    
     if(!t.PropagateTo(x,rad_length,rho)) return 0;
-    y = t.GetY();
-    ymax = x*TMath::Tan(0.5*alpha);
-    if (y > ymax) {
-      s = (s+1) % ns;
-      if (!t.Rotate(alpha)) return 0;
-    } else if (y <-ymax) {
-      s = (s-1+ns) % ns;                           
-      if (!t.Rotate(-alpha)) return 0;   
-    } 
+    AdjustSector(&t);
     if(!t.PropagateTo(x,rad_length,rho)) return 0;
 
-    y = t.GetY(); z = t.GetZ();
+    y = t.GetY(); 
+    z = t.GetZ();
 
     // now propagate to the middle plane of the next time bin 
     fTrSec[s]->GetLayer(nr-1)->GetPropagationParameters(y,z,dx,rho,rad_length,lookForCluster);
-    x = fTrSec[s]->GetLayer(nr-1)->GetX(); y = t.GetY(); z = t.GetZ();
+    x = fTrSec[s]->GetLayer(nr-1)->GetX(); 
+    
     if(!t.PropagateTo(x,rad_length,rho)) return 0;
-    y = t.GetY();
-    ymax = x*TMath::Tan(0.5*alpha);
-    if (y > ymax) {
-      s = (s+1) % ns;
-      if (!t.Rotate(alpha)) return 0;
-    } else if (y <-ymax) {
-      s = (s-1+ns) % ns;                           
-      if (!t.Rotate(-alpha)) return 0;   
-    } 
+    AdjustSector(&t);
     if(!t.PropagateTo(x,rad_length,rho)) return 0;
   }
   return 1;
 }         
 
-
 //_____________________________________________________________________________
 void AliTRDtracker::LoadEvent()
 {