]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDtracker.cxx
Moving the alignment-related static methods from AliAlignObj to the new geometry...
[u/mrichter/AliRoot.git] / TRD / AliTRDtracker.cxx
index aa8fd853a6f885353611faa6cc0eb507bbe49a44..32a621d5e1954a9d197762926efb21f2e3035fea 100644 (file)
 
 
 #include <Riostream.h>
-#include <TFile.h>
 #include <TBranch.h>
-#include <TTree.h>  
-#include <TObjArray.h> 
-#include <TTreeStream.h>
+#include <TFile.h>
 #include <TGraph.h>
-#include <TLinearFitter.h>
 #include <TH1D.h>
 #include <TH2D.h>
+#include <TLinearFitter.h>
+#include <TObjArray.h> 
+#include <TROOT.h>
+#include <TTree.h>  
+#include <TTreeStream.h>
 
 #include "AliESD.h"
 #include "AliAlignObj.h"
 #include "AliTRDCommonParam.h"
 #include "AliTRDtracker.h"
 #include "AliTRDReconstructor.h"
-#include "AliTRDCalibra.h"
+#include "AliTRDCalibraFillHisto.h"
+
 ClassImp(AliTRDtracker)
 
-const  Float_t  AliTRDtracker::fgkMinClustersInTrack =  0.5;  
-const  Float_t  AliTRDtracker::fgkLabelFraction      =  0.8;  // ??
-const  Double_t AliTRDtracker::fgkMaxChi2            = 12.0; 
+const  Float_t  AliTRDtracker::fgkMinClustersInTrack =  0.5;  //
+const  Float_t  AliTRDtracker::fgkLabelFraction      =  0.8;  //
+const  Double_t AliTRDtracker::fgkMaxChi2            = 12.0;  //
 const  Double_t AliTRDtracker::fgkMaxSnp             =  0.95; // Corresponds to tan = 3
 const  Double_t AliTRDtracker::fgkMaxStep            =  2.0;  // Maximal step size in propagation 
 
 //_____________________________________________________________________________
 AliTRDtracker::AliTRDtracker()
-  :AliTracker(),
-  fGeom(0),
-  fNclusters(0),
-  fClusters(0),
-  fNseeds(0),
-  fSeeds(0),
-  fNtracks(0),
-  fTracks(0),
-  fTimeBinsPerPlane(0),
-  fAddTRDseeds(kFALSE),
-  fNoTilt(kFALSE),
-  fDebugStreamer(0)
+  :AliTracker()
+  ,fHBackfit(0x0)
+  ,fHClSearch(0x0)
+  ,fHRefit(0x0)
+  ,fHX(0x0)
+  ,fHNCl(0x0)
+  ,fHNClTrack(0x0)
+  ,fHMinYPos(0x0)
+  ,fHMinYNeg(0x0)
+  ,fHMinZ(0x0)
+  ,fHMinD(0x0)
+  ,fHDeltaX(0x0)
+  ,fHXCl(0x0)
+  ,fGeom(0)
+  ,fNclusters(0)
+  ,fClusters(0)
+  ,fNseeds(0)
+  ,fSeeds(0)
+  ,fNtracks(0)
+  ,fTracks(0)
+  ,fTimeBinsPerPlane(0)
+  ,fAddTRDseeds(kFALSE)
+  ,fNoTilt(kFALSE)
+  ,fDebugStreamer(0)
 {
   //
   // Default constructor
@@ -91,11 +105,24 @@ AliTRDtracker::AliTRDtracker()
   }
 
   InitLogHists();
+
 } 
 
 //_____________________________________________________________________________
 AliTRDtracker::AliTRDtracker(const AliTRDtracker &t)
   :AliTracker(t)
+  ,fHBackfit(0x0)
+  ,fHClSearch(0x0)
+  ,fHRefit(0x0)
+  ,fHX(0x0)
+  ,fHNCl(0x0)
+  ,fHNClTrack(0x0)
+  ,fHMinYPos(0x0)
+  ,fHMinYNeg(0x0)
+  ,fHMinZ(0x0)
+  ,fHMinD(0x0)
+  ,fHDeltaX(0x0)
+  ,fHXCl(0x0)
   ,fGeom(0)
   ,fNclusters(0)
   ,fClusters(0)
@@ -111,11 +138,24 @@ AliTRDtracker::AliTRDtracker(const AliTRDtracker &t)
   //
   // Copy constructor
   //
+
 }
 
 //_____________________________________________________________________________
 AliTRDtracker::AliTRDtracker(const TFile *geomfile)
   :AliTracker()
+  ,fHBackfit(0x0)
+  ,fHClSearch(0x0)
+  ,fHRefit(0x0)
+  ,fHX(0x0)
+  ,fHNCl(0x0)
+  ,fHNClTrack(0x0)
+  ,fHMinYPos(0x0)
+  ,fHMinYNeg(0x0)
+  ,fHMinZ(0x0)
+  ,fHMinD(0x0)
+  ,fHDeltaX(0x0)
+  ,fHXCl(0x0)
   ,fGeom(0)
   ,fNclusters(0)
   ,fClusters(new TObjArray(2000))
@@ -153,7 +193,7 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile)
   savedir->cd();  
 
   for (Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
-    Int_t trS   = CookSectorIndex(geomS);
+    Int_t trS   = geomS;
     fTrSec[trS] = new AliTRDtrackingSector(fGeom,geomS);
     for (Int_t icham = 0; icham < AliTRDgeometry::kNcham; icham++) {
       fHoles[icham][trS] = fGeom->IsHole(0,icham,geomS);
@@ -166,7 +206,10 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile)
     fNoTilt = kTRUE;
   }
 
-  fTimeBinsPerPlane =  AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
+  if (!AliTRDcalibDB::Instance()) {
+    AliFatal("Could not get calibration object");
+  }
+  fTimeBinsPerPlane = AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
 
   fDebugStreamer    = new TTreeSRedirector("TRDdebug.root");
 
@@ -221,30 +264,30 @@ Int_t  AliTRDtracker::LocalToGlobalID(Int_t lid)
   Int_t  ichamber = fGeom->GetChamber(lid);
   Int_t  iplan    = fGeom->GetPlane(lid);
 
-  AliAlignObj::ELayerID iLayer = AliAlignObj::kTRD1;
+  AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1;
   switch (iplan) {
   case 0:
-    iLayer = AliAlignObj::kTRD1;
+    iLayer = AliGeomManager::kTRD1;
     break;
   case 1:
-    iLayer = AliAlignObj::kTRD2;
+    iLayer = AliGeomManager::kTRD2;
     break;
   case 2:
-    iLayer = AliAlignObj::kTRD3;
+    iLayer = AliGeomManager::kTRD3;
     break;
   case 3:
-    iLayer = AliAlignObj::kTRD4;
+    iLayer = AliGeomManager::kTRD4;
     break;
   case 4:
-    iLayer = AliAlignObj::kTRD5;
+    iLayer = AliGeomManager::kTRD5;
     break;
   case 5:
-    iLayer = AliAlignObj::kTRD6;
+    iLayer = AliGeomManager::kTRD6;
     break;
   };
 
   Int_t    modId = isector * fGeom->Ncham() + ichamber;
-  UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,modId);
+  UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,modId);
 
   return volid;
 
@@ -258,29 +301,29 @@ Int_t  AliTRDtracker::GlobalToLocalID(Int_t gid)
   // 
 
   Int_t modId    = 0;
-  AliAlignObj::ELayerID layerId = AliAlignObj::VolUIDToLayer(gid,modId);
+  AliGeomManager::ELayerID layerId = AliGeomManager::VolUIDToLayer(gid,modId);
 
   Int_t isector  = modId / fGeom->Ncham();
   Int_t ichamber = modId % fGeom->Ncham();
   Int_t iLayer   = -1;
 
   switch (layerId) {
-  case AliAlignObj::kTRD1:
+  case AliGeomManager::kTRD1:
     iLayer = 0;
     break;
-  case AliAlignObj::kTRD2:
+  case AliGeomManager::kTRD2:
     iLayer = 1;
     break;
-  case AliAlignObj::kTRD3:
+  case AliGeomManager::kTRD3:
     iLayer = 2;
     break;
-  case AliAlignObj::kTRD4:
+  case AliGeomManager::kTRD4:
     iLayer = 3;
     break;
-  case AliAlignObj::kTRD5:
+  case AliGeomManager::kTRD5:
     iLayer = 4;
     break;
-  case AliAlignObj::kTRD6:
+  case AliGeomManager::kTRD6:
     iLayer = 5;
     break;
   default:
@@ -323,7 +366,7 @@ Bool_t  AliTRDtracker::Transform(AliTRDcluster *cluster)
   // ExB correction
   //
   Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0);
-  Double_t exB    = AliTRDcalibDB::Instance()->GetOmegaTau(vdrift, -AliTracker::GetBz()*0.1);
+  Double_t exB    = AliTRDcalibDB::Instance()->GetOmegaTau(vdrift,-AliTracker::GetBz()*0.1);
 
   AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();  
   AliTRDpadPlane    *padPlane    = commonParam->GetPadPlane(plane,chamber);
@@ -336,9 +379,9 @@ Bool_t  AliTRDtracker::Transform(AliTRDcluster *cluster)
 
   cluster->SetY(cluster->GetY() - driftX*exB);
   Double_t xplane = (Double_t) AliTRDgeometry::GetTime0(plane); 
-  cluster->SetX(xplane- cluster->GetX());
+  cluster->SetX(xplane - cluster->GetX());
 
-  TGeoHMatrix *matrix =  fGeom->GetCorrectionMatrix(cluster->GetDetector());
+  TGeoHMatrix *matrix = fGeom->GetCorrectionMatrix(cluster->GetDetector());
   if (!matrix) {
     // No matrix found - if somebody used geometry with holes
     AliError("Invalid Geometry - Default Geometry used\n");
@@ -361,7 +404,7 @@ Bool_t  AliTRDtracker::Transform(AliTRDcluster *cluster)
   }
 
   if (plane == 5) {
-     cluster->SetX(localPosTracker[0]+kX0shift5);
+    cluster->SetX(localPosTracker[0]+kX0shift5);
   }
   else {
     cluster->SetX(localPosTracker[0]+kX0shift);
@@ -401,7 +444,7 @@ Bool_t  AliTRDtracker::Transform(AliTRDcluster *cluster)
 //   // ExB correction
 //   //
 //   Double_t vdrift = AliTRDcalibDB::Instance()->GetVdrift(cluster->GetDetector(),0,0);
-//   Double_t exB =   AliTRDcalibDB::Instance()->GetOmegaTau(vdrift);
+//   Double_t exB =   AliTRDcalibDB::Instance()->GetOmegaTau(vdrift,-AliTracker::GetBz()*0.1);
 //   //
 
 //   AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();  
@@ -752,11 +795,11 @@ Int_t AliTRDtracker::PropagateBack(AliESD *event)
     /**/
     // Debug part of tracking
     TTreeSRedirector &cstream = *fDebugStreamer;
-    Int_t eventNr = event->GetEventNumber();
+    Int_t eventNrInFile = event->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
     if (AliTRDReconstructor::StreamLevel() > 0) {
       if (track->GetBackupTrack()) {
        cstream << "Tracks"
-               << "EventNr="  << eventNr
+               << "EventNrInFile="  << eventNrInFile
                << "ESD.="     << seed
                << "trd.="     << track
                << "trdback.=" << track->GetBackupTrack()
@@ -764,7 +807,7 @@ Int_t AliTRDtracker::PropagateBack(AliESD *event)
       }
       else {
        cstream << "Tracks"
-               << "EventNr="  << eventNr
+               << "EventNrInFile="  << eventNrInFile
                << "ESD.="     << seed
                << "trd.="     << track
                << "trdback.=" << track
@@ -1135,9 +1178,10 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack &t)
          }
          maxChi2 = t.GetPredictedChi2(cl,h01);
           
-         if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
-           // ????
-         }
+         if (maxChi2<1e+10)
+           if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
+             // ????
+           }
 
        }                       
 
@@ -1174,7 +1218,7 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t)
   AliTRDtracklet tracklet;
 
   // Calibration fill 2D
-  AliTRDCalibra *calibra = AliTRDCalibra::Instance();
+  AliTRDCalibraFillHisto *calibra = AliTRDCalibraFillHisto::Instance();
   if (!calibra) {
     AliInfo("Could not get Calibra instance\n");
   }
@@ -1295,11 +1339,12 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack &t)
          t.PropagateTo(xcluster,radLength,rho);
          maxChi2 = t.GetPredictedChi2(cl,h01);
 
-         if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
-           if (!t.Update(cl,maxChi2,index,h01)) {
-             // ????
-           }
-          }  
+         if (maxChi2<1e+10)
+           if (!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
+             if (!t.Update(cl,maxChi2,index,h01)) {
+               // ????
+             }
+           }  
 
           if (calibra->GetMITracking()) {
             calibra->UpdateHistograms(cl,&t);
@@ -1416,7 +1461,7 @@ Int_t AliTRDtracker::LoadClusters(TTree *cTree)
     Int_t localTimeBin   = c->GetLocalTimeBin();
     Int_t sector         = fGeom->GetSector(detector);
     Int_t plane          = fGeom->GetPlane(detector);
-    Int_t trackingSector = CookSectorIndex(sector);
+    Int_t trackingSector = sector;
 
     //if (c->GetLabel(0) > 0) {
     if (c->GetQ() > 10) {
@@ -2583,7 +2628,7 @@ void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESD *esd)
          }
        }
 
-       Int_t eventNr = esd->GetEventNumber();
+       Int_t eventNrInFile = esd->GetEventNumberInFile(); // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number.
        TTreeSRedirector &cstream = *fDebugStreamer;
 
        //
@@ -2602,7 +2647,7 @@ void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESD *esd)
          TTreeSRedirector &cstream = *fDebugStreamer;
          if (AliTRDReconstructor::StreamLevel() > 0) {
            cstream << "Tracks"
-                   << "EventNr="  << eventNr
+                   << "EventNrInFile="  << eventNrInFile
                    << "ESD.="     << &esdtrack
                    << "trd.="     << track
                    << "trdback.=" << track
@@ -2632,7 +2677,7 @@ void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/, AliESD *esd)
                  << "Findable="  << findable
                  << "NUsed="     << nused
                  << "sLayer="    << sLayer
-                 << "EventNr="   << eventNr
+                 << "EventNrInFile="   << eventNrInFile
                  << "\n";
        }
 
@@ -2713,29 +2758,29 @@ Bool_t AliTRDtracker::GetTrackPoint(Int_t index, AliTrackPoint &p) const
   Double_t global[3];
   fGeom->RotateBack(idet,local,global);
   p.SetXYZ(global[0],global[1],global[2]);
-  AliAlignObj::ELayerID iLayer = AliAlignObj::kTRD1;
+  AliGeomManager::ELayerID iLayer = AliGeomManager::kTRD1;
   switch (iplan) {
   case 0:
-    iLayer = AliAlignObj::kTRD1;
+    iLayer = AliGeomManager::kTRD1;
     break;
   case 1:
-    iLayer = AliAlignObj::kTRD2;
+    iLayer = AliGeomManager::kTRD2;
     break;
   case 2:
-    iLayer = AliAlignObj::kTRD3;
+    iLayer = AliGeomManager::kTRD3;
     break;
   case 3:
-    iLayer = AliAlignObj::kTRD4;
+    iLayer = AliGeomManager::kTRD4;
     break;
   case 4:
-    iLayer = AliAlignObj::kTRD5;
+    iLayer = AliGeomManager::kTRD5;
     break;
   case 5:
-    iLayer = AliAlignObj::kTRD6;
+    iLayer = AliGeomManager::kTRD6;
     break;
   };
   Int_t    modId = isector * fGeom->Ncham() + ichamber;
-  UShort_t volid = AliAlignObj::LayerToVolUID(iLayer,modId);
+  UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,modId);
   p.SetVolumeID(volid);
 
   return kTRUE;
@@ -3021,7 +3066,7 @@ AliTRDtracker::AliTRDtrackingSector
     }
 
     dx        = AliTRDcalibDB::Instance()->GetVdrift(0,0,0)
-              / AliTRDcalibDB::Instance()->GetSamplingFrequency();
+              / commonParam->GetSamplingFrequency();
     rho       = 0.00295 * 0.85; //????
     radLength = 11.0;  
 
@@ -3463,7 +3508,16 @@ void AliTRDtracker::CookdEdxTimBin(AliTRDtrack &TRDtrack)
     }
     Int_t detector  = pTRDcluster->GetDetector();
     Int_t iPlane    = fGeom->GetPlane(detector);
-    Int_t iSlice    = tb * AliESDtrack::kNSlice / AliTRDtrack::kNtimeBins;
+    if (iPlane >= AliESDtrack::kNPlane) {
+      AliError(Form("Wrong plane %d",iPlane));
+      continue;
+    }
+    Int_t iSlice    = tb * AliESDtrack::kNSlice 
+                         / AliTRDcalibDB::Instance()->GetNumberOfTimeBins();
+    if (iSlice >= AliESDtrack::kNSlice) {
+      AliError(Form("Wrong slice %d",iSlice));
+      continue;
+    }
     clscharge[iPlane][iSlice] = clscharge[iPlane][iSlice] + charge;
     if (charge > maxclscharge[iPlane]) {
       maxclscharge[iPlane] = charge;
@@ -3972,8 +4026,8 @@ Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1
   //if (tchi2s[bestiter]>25.) sigma2*=tchi2s[bestiter]/25.;
   //if (tchi2s[bestiter]>25.) sigma2=1000.;  // dont'accept
 
-  Double_t exB         = AliTRDcalibDB::Instance()->GetOmegaTau(
-                          AliTRDcalibDB::Instance()->GetVdrift(0,0,0), -AliTracker::GetBz()*0.1);
+  Double_t exB         = AliTRDcalibDB::Instance()->GetOmegaTau(AliTRDcalibDB::Instance()->GetVdrift(0,0,0)
+                                                               ,-AliTracker::GetBz()*0.1);
   Double_t expectederr = sigma2*sigma2 + 0.01*0.01;
   if (mpads > 3.5) {
     expectederr += (mpads - 3.5) * 0.04;