]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/Cal/AliTRDCalPIDLQ.cxx
Rename AliTRDCalPIDLQRef
[u/mrichter/AliRoot.git] / TRD / Cal / AliTRDCalPIDLQ.cxx
index 909339bdc9b6834df092ca41fd8c79570c859246..0fe96d6ed991bdd8227ebf092a74c4fc31835554 100644 (file)
 #include <TH1F.h>
 #include <TH2F.h>
 #include <TFile.h>
-#include <TTree.h>
 #include <TROOT.h>
-#include <TPDGCode.h>
-#include <TParticle.h>
-#include <TPrincipal.h>
 
 #include "AliLog.h"
-#include "AliHeader.h"
-#include "AliGenEventHeader.h"
-#include "AliRun.h"
-#include "AliRunLoader.h"
-#include "AliStack.h"
 #include "AliPID.h"
 #include "AliESD.h"
 #include "AliESDtrack.h"
 
 #include "AliTRDCalPIDLQ.h"
-#include "AliTRDCalPIDLQRef.h"
-#include "AliTRDpidESD.h"
 #include "AliTRDcalibDB.h"
-#include "AliTRDtrack.h"
-#include "AliTRDgeometry.h"
 
 ClassImp(AliTRDCalPIDLQ)
 
 Char_t* AliTRDCalPIDLQ::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
+
 Char_t* AliTRDCalPIDLQ::fpartSymb[AliPID::kSPECIES] = {"EL", "MU", "PI", "KA", "PR"};
+
+Float_t AliTRDCalPIDLQ::fTrackMomentum[kNMom] = {0.6,  0.8,  1.0,  1.5,  2.0,  3.0,  4.0,  5.0,  6.0,  8.0,  10.0};
+  
+Float_t AliTRDCalPIDLQ::fTrackSegLength[kNLength] = {3.7, 3.9, 4.2, 5.0};
+
     
 //_________________________________________________________________________
 AliTRDCalPIDLQ::AliTRDCalPIDLQ()
   :TNamed("pid", "PID for TRD")
-  ,fNMom(0)
-  ,fNLength(0)
-  ,fTrackMomentum(0x0)
-  ,fTrackSegLength(0x0)
-  ,fNTimeBins(0)
   ,fMeanChargeRatio(0)
-  ,fNbins(0)
-  ,fBinSize(0)
   ,fHistdEdx(0x0)
   ,fHistTimeBin(0x0)
-  ,fEstimator(0x0)
 {
   //
   //  The Default constructor
@@ -83,17 +68,9 @@ AliTRDCalPIDLQ::AliTRDCalPIDLQ()
 //_________________________________________________________________________
 AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title) 
   :TNamed(name,title)
-  ,fNMom(0)
-  ,fNLength(0)
-  ,fTrackMomentum(0x0)
-  ,fTrackSegLength(0x0)
-  ,fNTimeBins(0)
   ,fMeanChargeRatio(0)
-  ,fNbins(0)
-  ,fBinSize(0)
   ,fHistdEdx(0x0)
   ,fHistTimeBin(0x0)
-  ,fEstimator(0x0)
 {
   //
   //  The main constructor
@@ -106,17 +83,9 @@ AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title)
 //_____________________________________________________________________________
 AliTRDCalPIDLQ::AliTRDCalPIDLQ(const AliTRDCalPIDLQ &c) 
   :TNamed(c)
-  ,fNMom(c.fNMom)
-  ,fNLength(c.fNLength)
-  ,fTrackMomentum(0x0)
-  ,fTrackSegLength(0x0)
-  ,fNTimeBins(c.fNTimeBins)
   ,fMeanChargeRatio(c.fMeanChargeRatio)
-  ,fNbins(c.fNbins)
-  ,fBinSize(c.fBinSize)
   ,fHistdEdx(0x0)
   ,fHistTimeBin(0x0)
-  ,fEstimator(0x0)
 {
   //
   // Copy constructor
@@ -153,22 +122,6 @@ void AliTRDCalPIDLQ::CleanUp()
     delete fHistTimeBin;
     fHistTimeBin = 0x0;
   }
-
-  if (fTrackMomentum) {
-    delete[] fTrackMomentum;
-    fTrackMomentum = 0x0;
-  }
-
-  if (fTrackSegLength) {
-    delete[] fTrackSegLength;
-    fTrackSegLength = 0x0;
-  }
-
-  if (fEstimator) {
-    delete fEstimator;
-    fEstimator = 0x0;
-  }
-
 }
 
 //_____________________________________________________________________________
@@ -194,22 +147,7 @@ void AliTRDCalPIDLQ::Copy(TObject &c) const
   
   target.CleanUp();
   
-  target.fNbins           = fNbins;
-  target.fBinSize         = fBinSize;
   target.fMeanChargeRatio = fMeanChargeRatio;
-  target.fNTimeBins       = fNTimeBins;
-
-  //target.fNMom            = fNMom;
-  target.fTrackMomentum = new Double_t[fNMom];
-  for (Int_t i=0; i<fNMom; ++i) {
-    target.fTrackMomentum[i] = fTrackMomentum[i];
-  }
-
-  //target.fNLength         = fNLength;
-  target.fTrackSegLength = new Double_t[fNLength];
-  for (Int_t i=0; i<fNLength; ++i) {
-    target.fTrackSegLength[i] = fTrackSegLength[i];
-  }
 
   if (fHistdEdx) {
     target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
@@ -218,8 +156,6 @@ void AliTRDCalPIDLQ::Copy(TObject &c) const
     target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
   }
 
-  target.fEstimator = new AliTRDCalPIDLQRef(*fEstimator);
-
   TObject::Copy(c);
 
 }
@@ -231,79 +167,49 @@ void AliTRDCalPIDLQ::Init()
   // Initialization
   //
 
-  fNMom    = 11;
-  fNLength =  4;
-
-  fTrackMomentum = new Double_t[fNMom];
-  fTrackMomentum[0] = 0.6;
-  fTrackMomentum[1] = 0.8;
-  fTrackMomentum[2] = 1.0;
-  fTrackMomentum[3] = 1.5;
-  fTrackMomentum[4] = 2.0;
-  fTrackMomentum[5] = 3.0;
-  fTrackMomentum[6] = 4.0;
-  fTrackMomentum[7] = 5.0;
-  fTrackMomentum[8] = 6.0;
-  fTrackMomentum[9] = 8.0;
-  fTrackMomentum[10] = 10.0;
-  
-  fTrackSegLength = new Double_t[fNLength];
-  fTrackSegLength[0] = 3.7;
-  fTrackSegLength[1] = 3.9;
-  fTrackSegLength[2] = 4.2;
-  fTrackSegLength[3] = 5.0;
-
-  fHistdEdx    = new TObjArray(AliPID::kSPECIES * fNMom/* * fNLength*/);
+  fHistdEdx    = new TObjArray(AliPID::kSPECIES * kNMom/* * kNLength*/);
   fHistdEdx->SetOwner();
-  fHistTimeBin = new TObjArray(2 * fNMom);
+  fHistTimeBin = new TObjArray(2 * kNMom);
   fHistTimeBin->SetOwner();  
 
        // Initialization of estimator at object instantiation because late
        // initialization in function GetProbability() is not working due to
        // constantness of this function. 
-       fEstimator = new AliTRDCalPIDLQRef();
-
-       // Number of Time bins
-       AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
-       if(!calibration){
-         AliWarning("No AliTRDcalibDB available. Using 22 time bins.");
-         fNTimeBins = 22;
-       } else {
-          if (calibration->GetRun() > -1) {
-            fNTimeBins = calibration->GetNumberOfTimeBins();
-         }
-          else {
-           AliWarning("No run number set. Using 22 time bins.");
-           fNTimeBins = 22;
-         }
-       }
+       // fEstimator = new AliTRDCalPIDRefMaker();
        
   // ADC Gain normalization
   fMeanChargeRatio = 1.0;
-  
-  // Number of bins and bin size
-  fNbins   = 0;
-  fBinSize = 0.0;
 }
 
 //_________________________________________________________________________
-Bool_t AliTRDCalPIDLQ::ReadReferences(Char_t *responseFile)
+Bool_t AliTRDCalPIDLQ::LoadLQReferences(Char_t *refFile)
 {
   //
   // Read the TRD dEdx histograms.
   //
 
+       Int_t nTimeBins = 22;
+       // Get number of time bins from CDB
+       AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
+       if(!calibration){
+               AliWarning(Form("No AliTRDcalibDB available. Using %d time bins.", nTimeBins));
+       }else{
+               if(calibration->GetRun() > -1) nTimeBins = calibration->GetNumberOfTimeBins();
+               else AliWarning(Form("Run number not set. Using %d time bins.", nTimeBins));
+       }
+
+       
   // Read histogram Root file  
-  TFile *histFile = new TFile(responseFile, "READ");
+  TFile *histFile = TFile::Open(refFile, "READ");
   if (!histFile || !histFile->IsOpen()) {
-    AliError(Form("Opening TRD histgram file %s failed", responseFile));    
+    AliError(Form("Opening TRD histgram file %s failed", refFile));
     return kFALSE;
   }
   gROOT->cd();
 
   // Read histograms
   for (Int_t iparticle = 0; iparticle < AliPID::kSPECIES; iparticle++){
-    for (Int_t imom = 0; imom < fNMom; imom++){
+    for (Int_t imom = 0; imom < kNMom; imom++){
                        TH2D* hist = (TH2D*)histFile->Get(Form("h2dEdx%s%02d", fpartSymb[iparticle], imom/*, ilength*/))->Clone();
                        hist->Scale(1./hist->Integral());
                        fHistdEdx->AddAt(hist, GetHistID(iparticle, imom));
@@ -311,9 +217,9 @@ Bool_t AliTRDCalPIDLQ::ReadReferences(Char_t *responseFile)
                        if (iparticle != AliPID::kElectron && iparticle != AliPID::kPion) continue;
 
                        TH1F* ht = (TH1F*)histFile->Get(Form("h1MaxTB%s%02d", fpartSymb[iparticle], imom))->Clone();
-                       if(ht->GetNbinsX() != fNTimeBins) AliWarning(Form("The number of time bins %d defined in h1MaxTB%s%02d differs from calibration value of %d. This may lead to erroneous results.", ht->GetNbinsX(), fpartSymb[iparticle], imom, fNTimeBins));
+                       if(ht->GetNbinsX() != nTimeBins) AliWarning(Form("The number of time bins %d defined in h1MaxTB%s%02d differs from calibration value of %d. This may lead to erroneous results.", ht->GetNbinsX(), fpartSymb[iparticle], imom, nTimeBins));
                        ht->Scale(1./ht->Integral());
-                       fHistTimeBin->AddAt(ht, ((iparticle==AliPID::kElectron)?0:1)*fNMom + imom);
+                       fHistTimeBin->AddAt(ht, ((iparticle==AliPID::kElectron)?0:1)*kNMom + imom);
                }
   }
   
@@ -373,7 +279,7 @@ TH1* AliTRDCalPIDLQ::GetHistogram(Int_t k, Int_t ip/*, Int_t il*/) const
   //
 
   if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
-       if(ip<0 || ip>= fNMom ) return 0x0;
+       if(ip<0 || ip>= kNMom ) return 0x0;
 
        AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));
   
@@ -389,15 +295,17 @@ TH1* AliTRDCalPIDLQ::GetHistogramT(Int_t k, Int_t ip) const
   //
 
   if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
-       if(ip<0 || ip>= fNMom ) return 0x0;
+       if(ip<0 || ip>= kNMom ) return 0x0;
          
        AliInfo(Form("Retrive MaxTB histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));
 
-       return (TH1*)fHistTimeBin->At(((k==AliPID::kElectron)?0:1)*fNMom+ip);
+       return (TH1*)fHistTimeBin->At(((k==AliPID::kElectron)?0:1)*kNMom+ip);
 }
 
+
+
 //_________________________________________________________________________
-Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Double_t mom, Double_t *dedx, Double_t length) const 
+Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Float_t mom, Float_t *dedx, Float_t length) const
 {
   //
        // Core function of AliTRDCalPIDLQ class for calculating the
@@ -412,11 +320,11 @@ Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Double_t mom, Double_t *dedx
        
        // find the interval in momentum and track segment length which applies for this data
        Int_t ilength = 1;
-  while(ilength<fNLength-1 && length>fTrackSegLength[ilength]){
+  while(ilength<kNLength-1 && length>fTrackSegLength[ilength]){
                ilength++;
        }
        Int_t imom = 1;
-  while(imom<fNMom-1 && mom>fTrackMomentum[imom]) imom++;
+  while(imom<kNMom-1 && mom>fTrackMomentum[imom]) imom++;
        
        Int_t nbinsx, nbinsy;
        TAxis *ax = 0x0, *ay = 0x0;
@@ -471,7 +379,6 @@ Double_t AliTRDCalPIDLQ::GetProbabilityT(Int_t spec, Double_t mom, Int_t timbin)
   //
   
        if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
-  if (timbin<0 || timbin >= fNTimeBins) return 0.;
 
   Int_t iTBin = timbin+1;
   
@@ -481,20 +388,20 @@ Double_t AliTRDCalPIDLQ::GetProbabilityT(Int_t spec, Double_t mom, Int_t timbin)
 
   
        Int_t imom = 1;
-  while(imom<fNMom-1 && mom>fTrackMomentum[imom]) imom++;
+  while(imom<kNMom-1 && mom>fTrackMomentum[imom]) imom++;
 
        Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
        TH1F *hist = 0x0;
-       if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*fNMom+imom-1))){
+       if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*kNMom+imom-1))){
                AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
-               AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*fNMom+imom-1));
+               AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*kNMom+imom-1));
                return 0.;
        }
        Double_t LQ1 = hist->GetBinContent(iTBin);
 
-       if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*fNMom+imom))){
+       if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*kNMom+imom))){
                AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
-               AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*fNMom+imom));
+               AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*kNMom+imom));
                return LQ1;
        }
        Double_t LQ2 = hist->GetBinContent(iTBin);
@@ -503,264 +410,3 @@ Double_t AliTRDCalPIDLQ::GetProbabilityT(Int_t spec, Double_t mom, Int_t timbin)
   return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
 }
 
-//__________________________________________________________________
-Bool_t AliTRDCalPIDLQ::WriteReferences(Char_t *File, Char_t *dir)
-{
-       // Build, Fill and write to file the histograms used for PID.
-       // The simulations are looked in the
-       // directories with the general form Form("p%3.1f", momentum)
-       // starting from dir (default .). Here momentum belongs to the list
-       // of known momentum to PID (see fTrackMomentum).
-       // The output histograms are
-       // written to the file "File" in cwd (default
-       // TRDPIDHistograms.root). In order to build a DB entry
-       // consider running $ALICE_ROOT/Cal/AliTRDCreateDummyCDB.C
-       // 
-       // Author:
-       // Alex Bercuci (A.Bercuci@gsi.de)
-
-       const Int_t nDirs = 1;
-  Int_t partCode[AliPID::kSPECIES] =
-    {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
-       
-       // minimal test of simulation location
-       TFile *f = new TFile(Form("%s/p%3.1f/galice.root", dir, 2.));
-       if(!f || f->IsZombie()){
-               AliError(Form("Could not access file galice in directry \"%s/p%3.1f\". Please check the location and try again.", dir, 2.));
-               return kFALSE;
-       }
-       f->Close(); delete f;
-       f = new TFile(Form("%s/p%3.1f/AliESDs.root", dir, 2.));
-       if(!f || f->IsZombie()){
-               AliError(Form("Could not access file AliESDs in directry \"%s/p%3.1f\". Please check the location and try again.", dir, 2.));
-               return kFALSE;
-       }
-       f->Close(); delete f;
-       
-
-       // Init statistics
-       Int_t nPart[AliPID::kSPECIES], nTotPart;
-       printf("P[GeV/c] ");
-       for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %s[%%] ", fpartSymb[ispec]);
-       printf("\n-----------------------------------------------\n");
-       
-       
-
-       // Build PID reference histograms and reference object
-       const Int_t color[] = {4, 3, 2, 7, 6};
-       for (Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
-               if (ispec != AliPID::kElectron && ispec != AliPID::kPion) continue;
-       
-               h1MaxTB[(ispec>0)?1:0] = new TH1F(Form("h1%s", fpartSymb[ispec]), "", fNTimeBins, -.5, fNTimeBins-.5);
-               h1MaxTB[(ispec>0)?1:0]->SetLineColor(color[ispec]);
-  }
-       AliTRDCalPIDLQRef ref;
-       
-       // momentum loop
-       AliRunLoader *fRunLoader = 0x0;
-       TFile *esdFile = 0x0;
-       TTree *esdTree = 0x0;
-       AliESD *esd = 0x0;
-       AliESDtrack *esdTrack = 0x0;
-       for (Int_t imom = 0; imom < fNMom; imom++) {
-       //for (Int_t imom = 4; imom < 5; imom++) {
-               ref.Reset();
-               
-               for(Int_t idir = 0; idir<nDirs; idir++){
-                       // open run loader and load gAlice, kinematics and header
-                       fRunLoader = AliRunLoader::Open(Form("%s/p%3.1f/galice.root", dir, fTrackMomentum[imom]));
-                       if (!fRunLoader) {
-                               AliError(Form("Getting run loader for momentum %3.1f failed.", fTrackMomentum[imom]));
-                               return kFALSE;
-                       }
-                       TString s; s.Form("%s/p%3.1f/", dir, fTrackMomentum[imom]);
-                       fRunLoader->SetDirName(s);
-                       fRunLoader->LoadgAlice();
-                       gAlice = fRunLoader->GetAliRun();
-                       if (!gAlice) {
-                               AliError(Form("galice object not found for momentum %3.1f.", fTrackMomentum[imom]));
-                               return kFALSE;
-                       }
-                       fRunLoader->LoadKinematics();
-                       fRunLoader->LoadHeader();
-       
-                       // open the ESD file
-                       esdFile = TFile::Open(Form("%s/p%3.1f/AliESDs.root", dir, fTrackMomentum[imom]));
-                       if (!esdFile || esdFile->IsZombie()) {
-                               AliError(Form("Opening ESD file failed for momentum", fTrackMomentum[imom]));
-                               return kFALSE;
-                       }
-                       esd = new AliESD;
-                       esdTree = (TTree*)esdFile->Get("esdTree");
-                       if (!esdTree) {
-                               AliError(Form("ESD tree not found for momentum %3.1f.", fTrackMomentum[imom]));
-                               return kFALSE;
-                       }
-                       esdTree->SetBranchAddress("ESD", &esd);
-                       nTotPart = 0;
-                       for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) nPart[ispec] = 0;
-       
-                       
-                       // Event loop
-                       for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
-                               fRunLoader->GetEvent(iEvent);
-       
-                               // read stack info
-                               AliStack* stack = gAlice->Stack();
-                               TArrayF vertex(3);
-                               fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
-                                                       
-                               // Load event summary data
-                               esdTree->GetEvent(iEvent);
-                               if (!esd) {
-                                       AliWarning(Form("ESD object not found for event %d. [@ momentum %3.1f]", iEvent, imom));
-                                       continue;
-                               }
-       
-                               for(Int_t iTrack=0; iTrack<esd->GetNumberOfTracks(); iTrack++){
-                                       esdTrack = esd->GetTrack(iTrack);
-       
-                                       if(!AliTRDpidESD::CheckTrack(esdTrack)) continue;
-                                       //if((esdTrack->GetStatus() & AliESDtrack::kITSrefit) == 0) continue;
-                                       //if(esdTrack->GetConstrainedChi2() > 1E9) continue;
-                                       //if ((esdTrack->GetStatus() & AliESDtrack::kESDpid) == 0) continue;
-                                       if (esdTrack->GetTRDsignal() == 0.) continue;
-       
-                                       // read MC info
-                                       Int_t label = esdTrack->GetLabel();
-                                       if(label<0) continue;
-                                       if (label > stack->GetNtrack()) continue;     // background
-                                       TParticle* particle = stack->Particle(label);
-                                       if(!particle){
-                                               AliWarning(Form("Retriving particle with index %d from AliStack failed. [@ event %d track %d]", label, iEvent, iTrack));
-                                               continue;
-                                       }
-                                       if(particle->Pt() < 1.E-3) continue;
-                                       //      if (TMath::Abs(particle->Eta()) > 0.3) continue;
-                                       TVector3 dVertex(particle->Vx() - vertex[0],
-                                                                               particle->Vy() - vertex[1],
-                                                                               particle->Vz() - vertex[2]);
-                                       if (dVertex.Mag() > 1.E-4){
-                                               //AliInfo(Form("Particle with index %d generated too far from vertex. Skip from analysis. Details follows. [@ event %d track %d]", label, iEvent, iTrack));
-                                               //particle->Print();
-                                               continue;
-                                       }
-                                       Int_t iGen = -1;
-                                       for (Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++)
-                                               if(TMath::Abs(particle->GetPdgCode()) == partCode[ispec]){
-                                                       iGen = ispec;
-                                                       break;
-                                               }
-                                       if(iGen<0) continue;
-       
-                                       nPart[iGen]++; nTotPart++;
-                                       
-                                       Float_t mom, length;
-                                       Double_t dedx[AliTRDtrack::kNslice], dEdx;
-                                       Int_t timebin;
-                                       for (Int_t iPlane=0; iPlane<AliTRDgeometry::kNplan; iPlane++){
-                                               // read data for track segment
-                                               for(int iSlice=0; iSlice<AliTRDtrack::kNslice; iSlice++)
-                                                       dedx[iSlice] = esdTrack->GetTRDsignals(iPlane, iSlice);
-                                               dEdx    = esdTrack->GetTRDsignals(iPlane, -1);
-                                               timebin = esdTrack->GetTRDTimBin(iPlane);
-                       
-                                               // check data
-                                               if ((dEdx <=  0.) || (timebin <= -1.)) continue;
-                       
-                                               // retrive kinematic info for this track segment
-                                               if(!AliTRDpidESD::GetTrackSegmentKine(esdTrack, iPlane, mom, length)) continue;
-                                               
-                                               // find segment length and momentum bin
-                                               Int_t jmom = 1, refMom = -1;
-                                               while(jmom<fNMom-1 && mom>fTrackMomentum[jmom]) jmom++;
-                                               if(TMath::Abs(fTrackMomentum[jmom-1] - mom) < fTrackMomentum[jmom-1] * .2) refMom = jmom-1;
-                                               else if(TMath::Abs(fTrackMomentum[jmom] - mom) < fTrackMomentum[jmom] * .2) refMom = jmom;
-                                               if(refMom<0){
-                                                       AliInfo(Form("Momentum at plane %d entrance not in momentum window. [@ event %d track %d]", iPlane, iEvent, iTrack));
-                                                       continue;
-                                               }
-                                               /*while(jleng<fNLength-1 && length>fTrackSegLength[jleng]) jleng++;*/
-                                               
-                                               // this track segment has fulfilled all requierments
-                                               //nPlanePID++;
-
-                                               if(dedx[0] > 0. && dedx[1] > 0.){
-                                                       dedx[0] = log(dedx[0]); dedx[1] = log(dedx[1]);
-                                                       ref.GetPrincipal(iGen)->AddRow(dedx);
-                                               }
-                                               h1MaxTB[(iGen>0)?1:0]->Fill(timebin);
-                                       } // end plane loop
-                               } // end track loop
-                       } // end events loop
-                       
-                       delete esd; esd = 0x0;
-                       esdFile->Close();
-                       delete esdFile; esdFile = 0x0;
-       
-                       fRunLoader->UnloadHeader();
-                       fRunLoader->UnloadKinematics();
-                       delete fRunLoader; fRunLoader = 0x0;
-               } // end directory loop
-               
-               // use data to prepare references
-               ref.Prepare2DReferences();
-               // save this dEdx references
-               ref.SaveReferences(imom, File);
-               // save MaxTB references
-               SaveMaxTimeBin(imom, File);
-
-                       
-               // print momentum statistics
-               printf("  %3.1f  ", fTrackMomentum[imom]);
-               for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %5.2f ", 100.*nPart[ispec]/nTotPart);
-               printf("\n");
-       } // end momentum loop
-       
-       TFile *fSave = 0x0;
-       TListIter it((TList*)gROOT->GetListOfFiles());
-       while((fSave=(TFile*)it.Next()))
-               if(strcmp(File, fSave->GetName())==0) break;
-
-       fSave->cd();
-       fSave->Close();
-       delete fSave;
-
-       return kTRUE;
-}
-
-//__________________________________________________________________
-void   AliTRDCalPIDLQ::SaveMaxTimeBin(const Int_t mom, const char *fn)
-{
-  //
-  // Save the histograms
-  //
-
-       TFile *fSave = 0x0;
-       TListIter it((TList*)gROOT->GetListOfFiles());
-       TDirectory *pwd = gDirectory;
-       Bool_t kFOUND = kFALSE;
-       while((fSave=(TFile*)it.Next()))
-               if(strcmp(fn, fSave->GetName())==0){
-                       kFOUND = kTRUE;
-                       break;
-               }
-       if(!kFOUND) fSave = new TFile(fn, "RECREATE");
-       fSave->cd();
-
-       TH1 *h;
-       h = (TH1F*)h1MaxTB[0]->Clone(Form("h1MaxTBEL%02d", mom));
-       h->SetTitle(Form("Maximum Time Bin distribution for electrons @ %4.1f GeV", fTrackMomentum[mom]));
-       h->GetXaxis()->SetTitle("time [100 ns]");
-       h->GetYaxis()->SetTitle("Probability");
-       h->Write();
-
-       h = (TH1F*)h1MaxTB[1]->Clone(Form("h1MaxTBPI%02d", mom));
-       h->SetTitle(Form("Maximum Time Bin distribution for pions @ %4.1f GeV", fTrackMomentum[mom]));
-       h->GetXaxis()->SetTitle("time [100 ns]");
-       h->GetYaxis()->SetTitle("Probability");
-       h->Write();
-       
-       pwd->cd();
-}
-