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