TRD/Cal/AliTRDCalPIDLQ.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 //////////////////////////////////////////////////////////////////////////
  19 //                                                                      //
  20 // Container for the distributions of dE/dx and the time bin of the     //
  21 // max. cluster for electrons and pions                                 //
  22 //                                                                      //
  23 // Authors:                                                             //
  24 //   Prashant Shukla <shukla@pi0.physi.uni-heidelberg.de>               //
  25 //   Alex Bercuci <a.bercuci@gsi.de>                                    //
  26 //                                                                      //
  27 //////////////////////////////////////////////////////////////////////////
  28
  29 #include <TH1F.h>
  30 #include <TH2F.h>
  31 #include <TFile.h>
  32 #include <TTree.h>
  33 #include <TROOT.h>
  34 #include <TPDGCode.h>
  35 #include <TParticle.h>
  36 #include <TPrincipal.h>
  37
  38 #include "AliLog.h"
  39 #include "AliHeader.h"
  40 #include "AliGenEventHeader.h"
  41 #include "AliRun.h"
  42 #include "AliRunLoader.h"
  43 #include "AliStack.h"
  44 #include "AliPID.h"
  45 #include "AliESD.h"
  46 #include "AliESDtrack.h"
  47
  48 #include "AliTRDCalPIDLQ.h"
  49 #include "AliTRDCalPIDLQRef.h"
  50 #include "AliTRDpidESD.h"
  51 #include "AliTRDcalibDB.h"
  52 #include "AliTRDtrack.h"
  53 #include "AliTRDgeometry.h"
  54
  55 ClassImp(AliTRDCalPIDLQ)
  56
  57 Char_t* AliTRDCalPIDLQ::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
  58 Char_t* AliTRDCalPIDLQ::fpartSymb[AliPID::kSPECIES] = {"EL", "MU", "PI", "KA", "PR"};
  59
  60 //_________________________________________________________________________
  61 AliTRDCalPIDLQ::AliTRDCalPIDLQ()
  62   :TNamed("pid", "PID for TRD")
  63   ,fNMom(0)
  64   ,fNLength(0)
  65   ,fTrackMomentum(0x0)
  66   ,fTrackSegLength(0x0)
  67   ,fNTimeBins(0)
  68   ,fMeanChargeRatio(0)
  69   ,fNbins(0)
  70   ,fBinSize(0)
  71   ,fHistdEdx(0x0)
  72   ,fHistTimeBin(0x0)
  73   ,fEstimator(0x0)
  74 {
  75   //
  76   //  The Default constructor
  77   //
  78
  79         Init();
  80 }
  81
  82 //_________________________________________________________________________
  83 AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title)
  84   :TNamed(name,title)
  85   ,fNMom(0)
  86   ,fNLength(0)
  87   ,fTrackMomentum(0x0)
  88   ,fTrackSegLength(0x0)
  89   ,fNTimeBins(0)
  90   ,fMeanChargeRatio(0)
  91   ,fNbins(0)
  92   ,fBinSize(0)
  93   ,fHistdEdx(0x0)
  94   ,fHistTimeBin(0x0)
  95   ,fEstimator(0x0)
  96 {
  97   //
  98   //  The main constructor
  99   //
 100
 101   Init();
 102
 103 }
 104
 105 //_____________________________________________________________________________
 106 AliTRDCalPIDLQ::AliTRDCalPIDLQ(const AliTRDCalPIDLQ &c)
 107   :TNamed(c)
 108   ,fNMom(c.fNMom)
 109   ,fNLength(c.fNLength)
 110   ,fTrackMomentum(0x0)
 111   ,fTrackSegLength(0x0)
 112   ,fNTimeBins(c.fNTimeBins)
 113   ,fMeanChargeRatio(c.fMeanChargeRatio)
 114   ,fNbins(c.fNbins)
 115   ,fBinSize(c.fBinSize)
 116   ,fHistdEdx(0x0)
 117   ,fHistTimeBin(0x0)
 118   ,fEstimator(0x0)
 119 {
 120   //
 121   // Copy constructor
 122   //
 123
 124   if (this != &c) ((AliTRDCalPIDLQ &) c).Copy(*this);
 125
 126 }
 127
 128 //_________________________________________________________________________
 129 AliTRDCalPIDLQ::~AliTRDCalPIDLQ()
 130 {
 131   //
 132   // Destructor
 133   //
 134
 135   CleanUp();
 136
 137 }
 138
 139 //_________________________________________________________________________
 140 void AliTRDCalPIDLQ::CleanUp()
 141 {
 142   //
 143   // Delets all newly created objects
 144   //
 145
 146   if (fHistdEdx) {
 147     delete fHistdEdx;
 148     fHistdEdx = 0x0;
 149   }
 150
 151   if (fHistTimeBin) {
 152     delete fHistTimeBin;
 153     fHistTimeBin = 0x0;
 154   }
 155
 156   if (fTrackMomentum) {
 157     delete[] fTrackMomentum;
 158     fTrackMomentum = 0x0;
 159   }
 160
 161   if (fTrackSegLength) {
 162     delete[] fTrackSegLength;
 163     fTrackSegLength = 0x0;
 164   }
 165
 166   if (fEstimator) {
 167     delete fEstimator;
 168     fEstimator = 0x0;
 169   }
 170
 171 }
 172
 173 //_____________________________________________________________________________
 174 AliTRDCalPIDLQ &AliTRDCalPIDLQ::operator=(const AliTRDCalPIDLQ &c)
 175 {
 176   //
 177   // Assignment operator
 178   //
 179
 180   if (this != &c) ((AliTRDCalPIDLQ &) c).Copy(*this);
 181   return *this;
 182
 183 }
 184
 185 //_____________________________________________________________________________
 186 void AliTRDCalPIDLQ::Copy(TObject &c) const
 187 {
 188   //
 189   // Copy function
 190   //
 191
 192   AliTRDCalPIDLQ& target = (AliTRDCalPIDLQ &) c;
 193
 194   target.CleanUp();
 195
 196   target.fNbins           = fNbins;
 197   target.fBinSize         = fBinSize;
 198   target.fMeanChargeRatio = fMeanChargeRatio;
 199   target.fNTimeBins       = fNTimeBins;
 200
 201   //target.fNMom            = fNMom;
 202   target.fTrackMomentum = new Double_t[fNMom];
 203   for (Int_t i=0; i<fNMom; ++i) {
 204     target.fTrackMomentum[i] = fTrackMomentum[i];
 205   }
 206
 207   //target.fNLength         = fNLength;
 208   target.fTrackSegLength = new Double_t[fNLength];
 209   for (Int_t i=0; i<fNLength; ++i) {
 210     target.fTrackSegLength[i] = fTrackSegLength[i];
 211   }
 212
 213   if (fHistdEdx) {
 214     target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
 215   }
 216   if (fHistTimeBin) {
 217     target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
 218   }
 219
 220   target.fEstimator = new AliTRDCalPIDLQRef(*fEstimator);
 221
 222   TObject::Copy(c);
 223
 224 }
 225
 226 //_________________________________________________________________________
 227 void AliTRDCalPIDLQ::Init()
 228 {
 229   //
 230   // Initialization
 231   //
 232
 233   fNMom    = 11;
 234   fNLength =  4;
 235
 236   fTrackMomentum = new Double_t[fNMom];
 237   fTrackMomentum[0] = 0.6;
 238   fTrackMomentum[1] = 0.8;
 239   fTrackMomentum[2] = 1.0;
 240   fTrackMomentum[3] = 1.5;
 241   fTrackMomentum[4] = 2.0;
 242   fTrackMomentum[5] = 3.0;
 243   fTrackMomentum[6] = 4.0;
 244   fTrackMomentum[7] = 5.0;
 245   fTrackMomentum[8] = 6.0;
 246   fTrackMomentum[9] = 8.0;
 247   fTrackMomentum[10] = 10.0;
 248
 249   fTrackSegLength = new Double_t[fNLength];
 250   fTrackSegLength[0] = 3.7;
 251   fTrackSegLength[1] = 3.9;
 252   fTrackSegLength[2] = 4.2;
 253   fTrackSegLength[3] = 5.0;
 254
 255   fHistdEdx    = new TObjArray(AliPID::kSPECIES * fNMom/* * fNLength*/);
 256   fHistdEdx->SetOwner();
 257   fHistTimeBin = new TObjArray(2 * fNMom);
 258   fHistTimeBin->SetOwner();
 259         // Initialization of estimator at object instantiation because late
 260         // initialization in function GetProbability() is not working due to
 261         // constantness of this function.
 262         fEstimator = new AliTRDCalPIDLQRef();
 263
 264         // Number of Time bins
 265         AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
 266         if(!calibration){
 267                 AliWarning("No AliTRDcalibDB available. Using 22 time bins.");
 268                 fNTimeBins = 22;
 269         } else fNTimeBins = calibration->GetNumberOfTimeBins();
 270
 271   // ADC Gain normalization
 272   fMeanChargeRatio = 1.0;
 273
 274   // Number of bins and bin size
 275   fNbins   = 0;
 276   fBinSize = 0.0;
 277 }
 278
 279 //_________________________________________________________________________
 280 Bool_t AliTRDCalPIDLQ::ReadReferences(Char_t *responseFile)
 281 {
 282   //
 283   // Read the TRD dEdx histograms.
 284   //
 285
 286   // Read histogram Root file
 287   TFile *histFile = new TFile(responseFile, "READ");
 288   if (!histFile || !histFile->IsOpen()) {
 289     AliError(Form("Opening TRD histgram file %s failed", responseFile));
 290     return kFALSE;
 291   }
 292   gROOT->cd();
 293
 294   // Read histograms
 295   for (Int_t iparticle = 0; iparticle < AliPID::kSPECIES; iparticle++){
 296     for (Int_t imom = 0; imom < fNMom; imom++){
 297                         TH2D* hist = (TH2D*)histFile->Get(Form("h2dEdx%s%02d", fpartSymb[iparticle], imom/*, ilength*/))->Clone();
 298                         hist->Scale(1./hist->Integral());
 299                         fHistdEdx->AddAt(hist, GetHistID(iparticle, imom));
 300
 301                         if (iparticle != AliPID::kElectron && iparticle != AliPID::kPion) continue;
 302
 303                         TH1F* ht = (TH1F*)histFile->Get(Form("h1MaxTB%s%02d", fpartSymb[iparticle], imom))->Clone();
 304                         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));
 305                         ht->Scale(1./ht->Integral());
 306                         fHistTimeBin->AddAt(ht, ((iparticle==AliPID::kElectron)?0:1)*fNMom + imom);
 307                 }
 308   }
 309
 310   histFile->Close();
 311   delete histFile;
 312
 313   // Number of bins and bin size
 314   //TH1F* hist = (TH1F*) fHistdEdx->At(GetHistID(AliPID::kPion, 1));
 315   //fNbins   = hist->GetNbinsX();
 316   //fBinSize = hist->GetBinWidth(1);
 317
 318   return kTRUE;
 319
 320 }
 321
 322 // //_________________________________________________________________________
 323 // Double_t  AliTRDCalPIDLQ::GetMean(Int_t k, Int_t ip) const
 324 // {
 325 //   //
 326 //   // Gets mean of de/dx dist. of e
 327 //   //
 328 //
 329 //   AliInfo(Form("Mean for particle = %s and momentum = %.2f is:\n"
 330 //               ,fpartName[k]
 331 //               ,fTrackMomentum[ip]));
 332 //   if (k < 0 || k > AliPID::kSPECIES) {
 333 //     return 0;
 334 //   }
 335 //
 336 //   return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->GetMean();
 337 //
 338 // }
 339 //
 340 // //_________________________________________________________________________
 341 // Double_t  AliTRDCalPIDLQ::GetNormalization(Int_t k, Int_t ip) const
 342 // {
 343 //   //
 344 //   // Gets Normalization of de/dx dist. of e
 345 //   //
 346 //
 347 //   AliInfo(Form("Normalization for particle = %s and momentum = %.2f is:\n"
 348 //               ,fpartName[k]
 349 //               ,fTrackMomentum[ip]));
 350 //   if (k < 0 || k > AliPID::kSPECIES) {
 351 //     return 0;
 352 //   }
 353 //
 354 //   return ((TH1F*) fHistdEdx->At(GetHistID(k,ip)))->Integral();
 355 //
 356 // }
 357
 358 //_________________________________________________________________________
 359 TH1* AliTRDCalPIDLQ::GetHistogram(Int_t k, Int_t ip/*, Int_t il*/) const
 360 {
 361   //
 362   // Returns one selected dEdx histogram
 363   //
 364
 365   if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
 366         if(ip<0 || ip>= fNMom ) return 0x0;
 367
 368         AliInfo(Form("Retrive dEdx histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));
 369
 370   return (TH1*)fHistdEdx->At(GetHistID(k, ip));
 371
 372 }
 373
 374 //_________________________________________________________________________
 375 TH1* AliTRDCalPIDLQ::GetHistogramT(Int_t k, Int_t ip) const
 376 {
 377   //
 378   // Returns one selected time bin max histogram
 379   //
 380
 381   if (k < 0 || k >= AliPID::kSPECIES) return 0x0;
 382         if(ip<0 || ip>= fNMom ) return 0x0;
 383
 384         AliInfo(Form("Retrive MaxTB histogram for %s of %5.2f GeV/c", fpartName[k], fTrackMomentum[ip]));
 385
 386         return (TH1*)fHistTimeBin->At(((k==AliPID::kElectron)?0:1)*fNMom+ip);
 387 }
 388
 389 //_________________________________________________________________________
 390 Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Double_t mom, Double_t *dedx, Double_t length) const
 391 {
 392   //
 393         // Core function of AliTRDCalPIDLQ class for calculating the
 394         // likelihood for species "spec" (see AliTRDtrack::kNspecie) of a
 395         // given momentum "mom" and a given dE/dx (slice "dedx") yield per
 396         // layer
 397   //
 398
 399         if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
 400
 401         //Double_t dedx   = dedx1/fMeanChargeRatio;
 402
 403         // find the interval in momentum and track segment length which applies for this data
 404         Int_t ilength = 1;
 405   while(ilength<fNLength-1 && length>fTrackSegLength[ilength]){
 406                 ilength++;
 407         }
 408         Int_t imom = 1;
 409   while(imom<fNMom-1 && mom>fTrackMomentum[imom]) imom++;
 410
 411         Int_t nbinsx, nbinsy;
 412         TAxis *ax = 0x0, *ay = 0x0;
 413         Double_t LQ1, LQ2;
 414         Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
 415         TH2 *hist = 0x0;
 416         if(!(hist = (TH2D*)fHistdEdx->At(GetHistID(spec, imom-1/*, ilength*/)))){
 417                 AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
 418                 AliError(Form("EHistogram id %d not found in DB.", GetHistID(spec, imom-1)));
 419                 return 0.;
 420         }
 421         ax = hist->GetXaxis(); nbinsx = ax->GetNbins();
 422         ay = hist->GetYaxis(); nbinsy = ay->GetNbins();
 423         Float_t x = dedx[0]+dedx[1], y = dedx[2];
 424   Bool_t kX = (x < ax->GetBinUpEdge(nbinsx));
 425         Bool_t kY = (y < ay->GetBinUpEdge(nbinsy));
 426         if(kX)
 427                 if(kY) LQ1 = hist->GetBinContent( hist->FindBin(x, y));
 428     //fEstimator->Estimate2D2(hist, x, y);
 429                 else LQ1 = hist->GetBinContent(ax->FindBin(x), nbinsy);
 430         else
 431                 if(kY) LQ1 = hist->GetBinContent(nbinsx, ay->FindBin(y));
 432                 else LQ1 = hist->GetBinContent(nbinsx, nbinsy);
 433
 434
 435         if(!(hist = (TH2D*)fHistdEdx->At(GetHistID(spec, imom/*, ilength*/)))){
 436                 AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
 437                 AliError(Form("EHistogram id %d not found in DB.", GetHistID(spec, imom)));
 438                 return LQ1;
 439         }
 440         if(kX)
 441                 if(kY) LQ2 = hist->GetBinContent( hist->FindBin(x, y));
 442     //fEstimator->Estimate2D2(hist, x, y);
 443                 else LQ2 = hist->GetBinContent(ax->FindBin(x), nbinsy);
 444         else
 445                 if(kY) LQ2 = hist->GetBinContent(nbinsx, ay->FindBin(y));
 446                 else LQ2 = hist->GetBinContent(nbinsx, nbinsy);
 447
 448
 449         // return interpolation over momentum binning
 450   return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
 451
 452 }
 453
 454 //_________________________________________________________________________
 455 Double_t AliTRDCalPIDLQ::GetProbabilityT(Int_t spec, Double_t mom, Int_t timbin) const
 456 {
 457   //
 458   // Gets the Probability of having timbin at a given momentum (mom)
 459   // and particle type (spec) (0 for e) and (2 for pi)
 460   // from the precalculated timbin distributions
 461   //
 462
 463         if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
 464   if (timbin<0 || timbin >= fNTimeBins) return 0.;
 465
 466   Int_t iTBin = timbin+1;
 467
 468   // Everything which is not an electron counts as a pion for time bin max
 469   //if(spec != AliPID::kElectron) spec = AliPID::kPion;
 470
 471
 472
 473         Int_t imom = 1;
 474   while(imom<fNMom-1 && mom>fTrackMomentum[imom]) imom++;
 475
 476         Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
 477         TH1F *hist = 0x0;
 478         if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*fNMom+imom-1))){
 479                 AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
 480                 AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*fNMom+imom-1));
 481                 return 0.;
 482         }
 483         Double_t LQ1 = hist->GetBinContent(iTBin);
 484
 485         if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*fNMom+imom))){
 486                 AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
 487                 AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*fNMom+imom));
 488                 return LQ1;
 489         }
 490         Double_t LQ2 = hist->GetBinContent(iTBin);
 491
 492         // return interpolation over momentum binning
 493   return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
 494 }
 495
 496 //__________________________________________________________________
 497 Bool_t AliTRDCalPIDLQ::WriteReferences(Char_t *File, Char_t *dir)
 498 {
 499         // Build, Fill and write to file the histograms used for PID.
 500         // The simulations are looked in the
 501         // directories with the general form Form("p%3.1f", momentum)
 502         // starting from dir (default .). Here momentum belongs to the list
 503         // of known momentum to PID (see fTrackMomentum).
 504         // The output histograms are
 505         // written to the file "File" in cwd (default
 506         // TRDPIDHistograms.root). In order to build a DB entry
 507         // consider running $ALICE_ROOT/Cal/AliTRDCreateDummyCDB.C
 508         //
 509         // Author:
 510         // Alex Bercuci (A.Bercuci@gsi.de)
 511
 512         const Int_t nDirs = 1;
 513   Int_t partCode[AliPID::kSPECIES] =
 514     {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
 515
 516         // minimal test of simulation location
 517         TFile *f = new TFile(Form("%s/p%3.1f/galice.root", dir, 2.));
 518         if(!f || f->IsZombie()){
 519                 AliError(Form("Could not access file galice in directry \"%s/p%3.1f\". Please check the location and try again.", dir, 2.));
 520                 return kFALSE;
 521         }
 522         f->Close(); delete f;
 523         f = new TFile(Form("%s/p%3.1f/AliESDs.root", dir, 2.));
 524         if(!f || f->IsZombie()){
 525                 AliError(Form("Could not access file AliESDs in directry \"%s/p%3.1f\". Please check the location and try again.", dir, 2.));
 526                 return kFALSE;
 527         }
 528         f->Close(); delete f;
 529
 530
 531         // Init statistics
 532         Int_t nPart[AliPID::kSPECIES], nTotPart;
 533         printf("P[GeV/c] ");
 534         for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %s[%%] ", fpartSymb[ispec]);
 535         printf("\n-----------------------------------------------\n");
 536
 537
 538
 539         // Build PID reference histograms and reference object
 540         const Int_t color[] = {4, 3, 2, 7, 6};
 541         for (Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 542                 if (ispec != AliPID::kElectron && ispec != AliPID::kPion) continue;
 543
 544                 h1MaxTB[(ispec>0)?1:0] = new TH1F(Form("h1%s", fpartSymb[ispec]), "", fNTimeBins, -.5, fNTimeBins-.5);
 545                 h1MaxTB[(ispec>0)?1:0]->SetLineColor(color[ispec]);
 546   }
 547         AliTRDCalPIDLQRef ref;
 548
 549         // momentum loop
 550         AliRunLoader *fRunLoader = 0x0;
 551         TFile *esdFile = 0x0;
 552         TTree *esdTree = 0x0;
 553         AliESD *esd = 0x0;
 554         AliESDtrack *esdTrack = 0x0;
 555         for (Int_t imom = 0; imom < fNMom; imom++) {
 556         //for (Int_t imom = 4; imom < 5; imom++) {
 557                 ref.Reset();
 558
 559                 for(Int_t idir = 0; idir<nDirs; idir++){
 560                         // open run loader and load gAlice, kinematics and header
 561                         fRunLoader = AliRunLoader::Open(Form("%s/p%3.1f/galice.root", dir, fTrackMomentum[imom]));
 562                         if (!fRunLoader) {
 563                                 AliError(Form("Getting run loader for momentum %3.1f failed.", fTrackMomentum[imom]));
 564                                 return kFALSE;
 565                         }
 566                         TString s; s.Form("%s/p%3.1f/", dir, fTrackMomentum[imom]);
 567                         fRunLoader->SetDirName(s);
 568                         fRunLoader->LoadgAlice();
 569                         gAlice = fRunLoader->GetAliRun();
 570                         if (!gAlice) {
 571                                 AliError(Form("galice object not found for momentum %3.1f.", fTrackMomentum[imom]));
 572                                 return kFALSE;
 573                         }
 574                         fRunLoader->LoadKinematics();
 575                         fRunLoader->LoadHeader();
 576
 577                         // open the ESD file
 578                         esdFile = TFile::Open(Form("%s/p%3.1f/AliESDs.root", dir, fTrackMomentum[imom]));
 579                         if (!esdFile || esdFile->IsZombie()) {
 580                                 AliError(Form("Opening ESD file failed for momentum", fTrackMomentum[imom]));
 581                                 return kFALSE;
 582                         }
 583                         esd = new AliESD;
 584                         esdTree = (TTree*)esdFile->Get("esdTree");
 585                         if (!esdTree) {
 586                                 AliError(Form("ESD tree not found for momentum %3.1f.", fTrackMomentum[imom]));
 587                                 return kFALSE;
 588                         }
 589                         esdTree->SetBranchAddress("ESD", &esd);
 590                         nTotPart = 0;
 591                         for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) nPart[ispec] = 0;
 592
 593
 594                         // Event loop
 595                         for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
 596                                 fRunLoader->GetEvent(iEvent);
 597
 598                                 // read stack info
 599                                 AliStack* stack = gAlice->Stack();
 600                                 TArrayF vertex(3);
 601                                 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
 602
 603                                 // Load event summary data
 604                                 esdTree->GetEvent(iEvent);
 605                                 if (!esd) {
 606                                         AliWarning(Form("ESD object not found for event %d. [@ momentum %3.1f]", iEvent, imom));
 607                                         continue;
 608                                 }
 609
 610                                 for(Int_t iTrack=0; iTrack<esd->GetNumberOfTracks(); iTrack++){
 611                                         esdTrack = esd->GetTrack(iTrack);
 612
 613                                         if(!AliTRDpidESD::CheckTrack(esdTrack)) continue;
 614                                         //if((esdTrack->GetStatus() & AliESDtrack::kITSrefit) == 0) continue;
 615                                         //if(esdTrack->GetConstrainedChi2() > 1E9) continue;
 616                                         //if ((esdTrack->GetStatus() & AliESDtrack::kESDpid) == 0) continue;
 617                                         if (esdTrack->GetTRDsignal() == 0.) continue;
 618
 619                                         // read MC info
 620                                         Int_t label = esdTrack->GetLabel();
 621                                         if(label<0) continue;
 622                                         if (label > stack->GetNtrack()) continue;     // background
 623                                         TParticle* particle = stack->Particle(label);
 624                                         if(!particle){
 625                                                 AliWarning(Form("Retriving particle with index %d from AliStack failed. [@ event %d track %d]", label, iEvent, iTrack));
 626                                                 continue;
 627                                         }
 628                                         if(particle->Pt() < 1.E-3) continue;
 629                                         //      if (TMath::Abs(particle->Eta()) > 0.3) continue;
 630                                         TVector3 dVertex(particle->Vx() - vertex[0],
 631                                                                                 particle->Vy() - vertex[1],
 632                                                                                 particle->Vz() - vertex[2]);
 633                                         if (dVertex.Mag() > 1.E-4){
 634                                                 //AliInfo(Form("Particle with index %d generated too far from vertex. Skip from analysis. Details follows. [@ event %d track %d]", label, iEvent, iTrack));
 635                                                 //particle->Print();
 636                                                 continue;
 637                                         }
 638                                         Int_t iGen = -1;
 639                                         for (Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++)
 640                                                 if(TMath::Abs(particle->GetPdgCode()) == partCode[ispec]){
 641                                                         iGen = ispec;
 642                                                         break;
 643                                                 }
 644                                         if(iGen<0) continue;
 645
 646                                         nPart[iGen]++; nTotPart++;
 647
 648                                         Float_t mom, length;
 649                                         Double_t dedx[AliTRDtrack::kNslice], dEdx;
 650                                         Int_t timebin;
 651                                         for (Int_t iPlane=0; iPlane<AliTRDgeometry::kNplan; iPlane++){
 652                                                 // read data for track segment
 653                                                 for(int iSlice=0; iSlice<AliTRDtrack::kNslice; iSlice++)
 654                                                         dedx[iSlice] = esdTrack->GetTRDsignals(iPlane, iSlice);
 655                                                 dEdx    = esdTrack->GetTRDsignals(iPlane, -1);
 656                                                 timebin = esdTrack->GetTRDTimBin(iPlane);
 657
 658                                                 // check data
 659                                                 if ((dEdx <=  0.) || (timebin <= -1.)) continue;
 660
 661                                                 // retrive kinematic info for this track segment
 662                                                 if(!AliTRDpidESD::GetTrackSegmentKine(esdTrack, iPlane, mom, length)) continue;
 663
 664                                                 // find segment length and momentum bin
 665                                                 Int_t jmom = 1, refMom = -1;
 666                                                 while(jmom<fNMom-1 && mom>fTrackMomentum[jmom]) jmom++;
 667                                                 if(TMath::Abs(fTrackMomentum[jmom-1] - mom) < fTrackMomentum[jmom-1] * .2) refMom = jmom-1;
 668                                                 else if(TMath::Abs(fTrackMomentum[jmom] - mom) < fTrackMomentum[jmom] * .2) refMom = jmom;
 669                                                 if(refMom<0){
 670                                                         AliInfo(Form("Momentum at plane %d entrance not in momentum window. [@ event %d track %d]", iPlane, iEvent, iTrack));
 671                                                         continue;
 672                                                 }
 673                                                 /*while(jleng<fNLength-1 && length>fTrackSegLength[jleng]) jleng++;*/
 674
 675                                                 // this track segment has fulfilled all requierments
 676                                                 //nPlanePID++;
 677
 678                                                 if(dedx[0] > 0. && dedx[1] > 0.){
 679                                                         dedx[0] = log(dedx[0]); dedx[1] = log(dedx[1]);
 680                                                         ref.GetPrincipal(iGen)->AddRow(dedx);
 681                                                 }
 682                                                 h1MaxTB[(iGen>0)?1:0]->Fill(timebin);
 683                                         } // end plane loop
 684                                 } // end track loop
 685                         } // end events loop
 686
 687                         delete esd; esd = 0x0;
 688                         esdFile->Close();
 689                         delete esdFile; esdFile = 0x0;
 690
 691                         fRunLoader->UnloadHeader();
 692                         fRunLoader->UnloadKinematics();
 693                         delete fRunLoader; fRunLoader = 0x0;
 694                 } // end directory loop
 695
 696                 // use data to prepare references
 697                 ref.Prepare2DReferences();
 698                 // save this dEdx references
 699                 ref.SaveReferences(imom, File);
 700                 // save MaxTB references
 701                 SaveMaxTimeBin(imom, File);
 702
 703
 704                 // print momentum statistics
 705                 printf("  %3.1f  ", fTrackMomentum[imom]);
 706                 for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %5.2f ", 100.*nPart[ispec]/nTotPart);
 707                 printf("\n");
 708         } // end momentum loop
 709
 710         TFile *fSave = 0x0;
 711         TListIter it((TList*)gROOT->GetListOfFiles());
 712         while((fSave=(TFile*)it.Next()))
 713                 if(strcmp(File, fSave->GetName())==0) break;
 714
 715         fSave->cd();
 716         fSave->Close();
 717         delete fSave;
 718
 719         return kTRUE;
 720 }
 721
 722 //__________________________________________________________________
 723 void    AliTRDCalPIDLQ::SaveMaxTimeBin(const Int_t mom, const char *fn)
 724 {
 725   //
 726   // Save the histograms
 727   //
 728
 729         TFile *fSave = 0x0;
 730         TListIter it((TList*)gROOT->GetListOfFiles());
 731         TDirectory *pwd = gDirectory;
 732         Bool_t kFOUND = kFALSE;
 733         while((fSave=(TFile*)it.Next()))
 734                 if(strcmp(fn, fSave->GetName())==0){
 735                         kFOUND = kTRUE;
 736                         break;
 737                 }
 738         if(!kFOUND) fSave = new TFile(fn, "RECREATE");
 739         fSave->cd();
 740
 741         TH1 *h;
 742         h = (TH1F*)h1MaxTB[0]->Clone(Form("h1MaxTBEL%02d", mom));
 743         h->SetTitle(Form("Maximum Time Bin distribution for electrons @ %4.1f GeV", fTrackMomentum[mom]));
 744         h->GetXaxis()->SetTitle("time [100 ns]");
 745         h->GetYaxis()->SetTitle("Probability");
 746         h->Write();
 747
 748         h = (TH1F*)h1MaxTB[1]->Clone(Form("h1MaxTBPI%02d", mom));
 749         h->SetTitle(Form("Maximum Time Bin distribution for pions @ %4.1f GeV", fTrackMomentum[mom]));
 750         h->GetXaxis()->SetTitle("time [100 ns]");
 751         h->GetYaxis()->SetTitle("Probability");
 752         h->Write();
 753
 754         pwd->cd();
 755 }
 756