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   ,fTrackMomentum(0x0)
  65   ,fNLength(0)
  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   ,fTrackMomentum(0x0)
  87   ,fNLength(0)
  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   ,fTrackMomentum(0x0)
 110   ,fNLength(c.fNLength)
 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         Bool_t kX = (dedx[0] < ax->GetBinUpEdge(nbinsx));
 424         Bool_t kY = (dedx[1] < ay->GetBinUpEdge(nbinsy));
 425         if(kX)
 426                 if(kY) LQ1 = hist->GetBinContent( hist->FindBin(dedx[0], dedx[1])); //fEstimator->Estimate2D2(hist, (Float_t&)dedx[0], (Float_t&)dedx[1]);
 427                 else LQ1 = hist->GetBinContent(ax->FindBin(dedx[0]), nbinsy);
 428         else
 429                 if(kY) LQ1 = hist->GetBinContent(nbinsx, ay->FindBin(dedx[1]));
 430                 else LQ1 = hist->GetBinContent(nbinsx, nbinsy);
 431
 432
 433         if(!(hist = (TH2D*)fHistdEdx->At(GetHistID(spec, imom/*, ilength*/)))){
 434                 AliInfo(Form("Looking for spec(%d) mom(%f) Ex(%f) Ey(%f) length(%f)", spec, mom, dedx[0], dedx[1], length));
 435                 AliError(Form("EHistogram id %d not found in DB.", GetHistID(spec, imom)));
 436                 return LQ1;
 437         }
 438         if(kX)
 439                 if(kY) LQ2 = hist->GetBinContent( hist->FindBin(dedx[0], dedx[1])); //fEstimator->Estimate2D2(hist, (Float_t&)dedx[0], (Float_t&)dedx[1]);
 440                 else LQ2 = hist->GetBinContent(ax->FindBin(dedx[0]), nbinsy);
 441         else
 442                 if(kY) LQ2 = hist->GetBinContent(nbinsx, ay->FindBin(dedx[1]));
 443                 else LQ2 = hist->GetBinContent(nbinsx, nbinsy);
 444
 445
 446         // return interpolation over momentum binning
 447   return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
 448
 449 }
 450
 451 //_________________________________________________________________________
 452 Double_t AliTRDCalPIDLQ::GetProbabilityT(Int_t spec, Double_t mom, Int_t timbin) const
 453 {
 454   //
 455   // Gets the Probability of having timbin at a given momentum (mom)
 456   // and particle type (spec) (0 for e) and (2 for pi)
 457   // from the precalculated timbin distributions
 458   //
 459
 460         if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
 461   if (timbin<0 || timbin >= fNTimeBins) return 0.;
 462
 463   Int_t iTBin = timbin+1;
 464
 465   // Everything which is not an electron counts as a pion for time bin max
 466   //if(spec != AliPID::kElectron) spec = AliPID::kPion;
 467
 468
 469
 470         Int_t imom = 1;
 471   while(imom<fNMom-1 && mom>fTrackMomentum[imom]) imom++;
 472
 473         Double_t mom1 = fTrackMomentum[imom-1], mom2 = fTrackMomentum[imom];
 474         TH1F *hist = 0x0;
 475         if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*fNMom+imom-1))){
 476                 AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
 477                 AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*fNMom+imom-1));
 478                 return 0.;
 479         }
 480         Double_t LQ1 = hist->GetBinContent(iTBin);
 481
 482         if(!(hist = (TH1F*) fHistTimeBin->At(((spec==AliPID::kElectron)?0:1)*fNMom+imom))){
 483                 AliInfo(Form("Looking for spec(%d) mom(%f) timbin(%d)", spec, mom, timbin));
 484                 AliError(Form("THistogram id %d not found in DB.", ((spec==AliPID::kElectron)?0:1)*fNMom+imom));
 485                 return LQ1;
 486         }
 487         Double_t LQ2 = hist->GetBinContent(iTBin);
 488
 489         // return interpolation over momentum binning
 490   return LQ1 + (LQ2 - LQ1)*(mom - mom1)/(mom2 - mom1);
 491 }
 492
 493 //__________________________________________________________________
 494 Bool_t AliTRDCalPIDLQ::WriteReferences(Char_t *File, Char_t *dir)
 495 {
 496         // Build, Fill and write to file the histograms used for PID.
 497         // The simulations are looked in the
 498         // directories with the general form Form("p%3.1f", momentum)
 499         // starting from dir (default .). Here momentum belongs to the list
 500         // of known momentum to PID (see fTrackMomentum).
 501         // The output histograms are
 502         // written to the file "File" in cwd (default
 503         // TRDPIDHistograms.root). In order to build a DB entry
 504         // consider running $ALICE_ROOT/Cal/AliTRDCreateDummyCDB.C
 505         //
 506         // Author:
 507         // Alex Bercuci (A.Bercuci@gsi.de)
 508
 509         const Int_t nDirs = 1;
 510   Int_t partCode[AliPID::kSPECIES] =
 511     {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
 512
 513         // minimal test of simulation location
 514         TFile *f = new TFile(Form("%s/p%3.1f/galice.root", dir, 2.));
 515         if(!f || f->IsZombie()){
 516                 AliError(Form("Could not access file galice in directry \"%s/p%3.1f\". Please check the location and try again.", dir, 2.));
 517                 return kFALSE;
 518         }
 519         f->Close(); delete f;
 520         f = new TFile(Form("%s/p%3.1f/AliESDs.root", dir, 2.));
 521         if(!f || f->IsZombie()){
 522                 AliError(Form("Could not access file AliESDs in directry \"%s/p%3.1f\". Please check the location and try again.", dir, 2.));
 523                 return kFALSE;
 524         }
 525         f->Close(); delete f;
 526
 527
 528         // Init statistics
 529         Int_t nPart[AliPID::kSPECIES], nTotPart;
 530         printf("P[GeV/c] ");
 531         for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %s[%%] ", fpartSymb[ispec]);
 532         printf("\n-----------------------------------------------\n");
 533
 534
 535
 536         // Build PID reference histograms and reference object
 537         const Int_t color[] = {4, 3, 2, 7, 6};
 538         for (Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
 539                 if (ispec != AliPID::kElectron && ispec != AliPID::kPion) continue;
 540
 541                 h1MaxTB[(ispec>0)?1:0] = new TH1F(Form("h1%s", fpartSymb[ispec]), "", fNTimeBins, -.5, fNTimeBins-.5);
 542                 h1MaxTB[(ispec>0)?1:0]->SetLineColor(color[ispec]);
 543   }
 544         AliTRDCalPIDLQRef ref;
 545
 546         // momentum loop
 547         AliRunLoader *fRunLoader = 0x0;
 548         TFile *esdFile = 0x0;
 549         TTree *esdTree = 0x0;
 550         AliESD *esd = 0x0;
 551         AliESDtrack *esdTrack = 0x0;
 552         for (Int_t imom = 0; imom < fNMom; imom++) {
 553         //for (Int_t imom = 4; imom < 5; imom++) {
 554                 ref.Reset();
 555
 556                 for(Int_t idir = 0; idir<nDirs; idir++){
 557                         // open run loader and load gAlice, kinematics and header
 558                         fRunLoader = AliRunLoader::Open(Form("%s/p%3.1f/galice.root", dir, fTrackMomentum[imom]));
 559                         if (!fRunLoader) {
 560                                 AliError(Form("Getting run loader for momentum %3.1f failed.", fTrackMomentum[imom]));
 561                                 return kFALSE;
 562                         }
 563                         TString s; s.Form("%s/p%3.1f/", dir, fTrackMomentum[imom]);
 564                         fRunLoader->SetDirName(s);
 565                         fRunLoader->LoadgAlice();
 566                         gAlice = fRunLoader->GetAliRun();
 567                         if (!gAlice) {
 568                                 AliError(Form("galice object not found for momentum %3.1f.", fTrackMomentum[imom]));
 569                                 return kFALSE;
 570                         }
 571                         fRunLoader->LoadKinematics();
 572                         fRunLoader->LoadHeader();
 573
 574                         // open the ESD file
 575                         esdFile = TFile::Open(Form("%s/p%3.1f/AliESDs.root", dir, fTrackMomentum[imom]));
 576                         if (!esdFile || esdFile->IsZombie()) {
 577                                 AliError(Form("Opening ESD file failed for momentum", fTrackMomentum[imom]));
 578                                 return kFALSE;
 579                         }
 580                         esd = new AliESD;
 581                         esdTree = (TTree*)esdFile->Get("esdTree");
 582                         if (!esdTree) {
 583                                 AliError(Form("ESD tree not found for momentum %3.1f.", fTrackMomentum[imom]));
 584                                 return kFALSE;
 585                         }
 586                         esdTree->SetBranchAddress("ESD", &esd);
 587                         nTotPart = 0;
 588                         for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) nPart[ispec] = 0;
 589
 590
 591                         // Event loop
 592                         for (Int_t iEvent = 0; iEvent < fRunLoader->GetNumberOfEvents(); iEvent++) {
 593                                 fRunLoader->GetEvent(iEvent);
 594
 595                                 // read stack info
 596                                 AliStack* stack = gAlice->Stack();
 597                                 TArrayF vertex(3);
 598                                 fRunLoader->GetHeader()->GenEventHeader()->PrimaryVertex(vertex);
 599
 600                                 // Load event summary data
 601                                 esdTree->GetEvent(iEvent);
 602                                 if (!esd) {
 603                                         AliWarning(Form("ESD object not found for event %d. [@ momentum %3.1f]", iEvent, imom));
 604                                         continue;
 605                                 }
 606
 607                                 for(Int_t iTrack=0; iTrack<esd->GetNumberOfTracks(); iTrack++){
 608                                         esdTrack = esd->GetTrack(iTrack);
 609
 610                                         if(!AliTRDpidESD::CheckTrack(esdTrack)) continue;
 611                                         //if((esdTrack->GetStatus() & AliESDtrack::kITSrefit) == 0) continue;
 612                                         //if(esdTrack->GetConstrainedChi2() > 1E9) continue;
 613                                         //if ((esdTrack->GetStatus() & AliESDtrack::kESDpid) == 0) continue;
 614                                         if (esdTrack->GetTRDsignal() == 0.) continue;
 615
 616                                         // read MC info
 617                                         Int_t label = esdTrack->GetLabel();
 618                                         if(label<0) continue;
 619                                         if (label > stack->GetNtrack()) continue;     // background
 620                                         TParticle* particle = stack->Particle(label);
 621                                         if(!particle){
 622                                                 AliWarning(Form("Retriving particle with index %d from AliStack failed. [@ event %d track %d]", label, iEvent, iTrack));
 623                                                 continue;
 624                                         }
 625                                         if(particle->Pt() < 1.E-3) continue;
 626                                         //      if (TMath::Abs(particle->Eta()) > 0.3) continue;
 627                                         TVector3 dVertex(particle->Vx() - vertex[0],
 628                                                                                 particle->Vy() - vertex[1],
 629                                                                                 particle->Vz() - vertex[2]);
 630                                         if (dVertex.Mag() > 1.E-4){
 631                                                 //AliInfo(Form("Particle with index %d generated too far from vertex. Skip from analysis. Details follows. [@ event %d track %d]", label, iEvent, iTrack));
 632                                                 //particle->Print();
 633                                                 continue;
 634                                         }
 635                                         Int_t iGen = -1;
 636                                         for (Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++)
 637                                                 if(TMath::Abs(particle->GetPdgCode()) == partCode[ispec]){
 638                                                         iGen = ispec;
 639                                                         break;
 640                                                 }
 641                                         if(iGen<0) continue;
 642
 643                                         nPart[iGen]++; nTotPart++;
 644
 645                                         Float_t mom, length;
 646                                         Double_t dedx[AliTRDtrack::kNslice], dEdx;
 647                                         Int_t timebin;
 648                                         for (Int_t iPlane=0; iPlane<AliTRDgeometry::kNplan; iPlane++){
 649                                                 // read data for track segment
 650                                                 for(int iSlice=0; iSlice<AliTRDtrack::kNslice; iSlice++)
 651                                                         dedx[iSlice] = esdTrack->GetTRDsignals(iPlane, iSlice);
 652                                                 dEdx    = esdTrack->GetTRDsignals(iPlane, -1);
 653                                                 timebin = esdTrack->GetTRDTimBin(iPlane);
 654
 655                                                 // check data
 656                                                 if ((dEdx <=  0.) || (timebin <= -1.)) continue;
 657
 658                                                 // retrive kinematic info for this track segment
 659                                                 if(!AliTRDpidESD::GetTrackSegmentKine(esdTrack, iPlane, mom, length)) continue;
 660
 661                                                 // find segment length and momentum bin
 662                                                 Int_t jmom = 1, refMom = -1;
 663                                                 while(jmom<fNMom-1 && mom>fTrackMomentum[jmom]) jmom++;
 664                                                 if(TMath::Abs(fTrackMomentum[jmom-1] - mom) < fTrackMomentum[jmom-1] * .2) refMom = jmom-1;
 665                                                 else if(TMath::Abs(fTrackMomentum[jmom] - mom) < fTrackMomentum[jmom] * .2) refMom = jmom;
 666                                                 if(refMom<0){
 667                                                         AliInfo(Form("Momentum at plane %d entrance not in momentum window. [@ event %d track %d]", iPlane, iEvent, iTrack));
 668                                                         continue;
 669                                                 }
 670                                                 /*while(jleng<fNLength-1 && length>fTrackSegLength[jleng]) jleng++;*/
 671
 672                                                 // this track segment has fulfilled all requierments
 673                                                 //nPlanePID++;
 674
 675                                                 if(dedx[0] > 0. && dedx[1] > 0.){
 676                                                         dedx[0] = log(dedx[0]); dedx[1] = log(dedx[1]);
 677                                                         ref.GetPrincipal(iGen)->AddRow(dedx);
 678                                                 }
 679                                                 h1MaxTB[(iGen>0)?1:0]->Fill(timebin);
 680                                         } // end plane loop
 681                                 } // end track loop
 682                         } // end events loop
 683
 684                         delete esd; esd = 0x0;
 685                         esdFile->Close();
 686                         delete esdFile; esdFile = 0x0;
 687
 688                         fRunLoader->UnloadHeader();
 689                         fRunLoader->UnloadKinematics();
 690                         delete fRunLoader; fRunLoader = 0x0;
 691                 } // end directory loop
 692
 693                 // use data to prepare references
 694                 ref.Prepare2DReferences();
 695                 // save this dEdx references
 696                 ref.SaveReferences(imom, File);
 697                 // save MaxTB references
 698                 SaveMaxTimeBin(imom, File);
 699
 700
 701                 // print momentum statistics
 702                 printf("  %3.1f  ", fTrackMomentum[imom]);
 703                 for(Int_t ispec=0; ispec<AliPID::kSPECIES; ispec++) printf(" %5.2f ", 100.*nPart[ispec]/nTotPart);
 704                 printf("\n");
 705         } // end momentum loop
 706
 707         TFile *fSave = 0x0;
 708         TListIter it((TList*)gROOT->GetListOfFiles());
 709         while((fSave=(TFile*)it.Next()))
 710                 if(strcmp(File, fSave->GetName())==0) break;
 711
 712         fSave->cd();
 713         fSave->Close();
 714         delete fSave;
 715
 716         return kTRUE;
 717 }
 718
 719 //__________________________________________________________________
 720 void    AliTRDCalPIDLQ::SaveMaxTimeBin(const Int_t mom, const char *fn)
 721 {
 722   //
 723   // Save the histograms
 724   //
 725
 726         TFile *fSave = 0x0;
 727         TListIter it((TList*)gROOT->GetListOfFiles());
 728         TDirectory *pwd = gDirectory;
 729         Bool_t kFOUND = kFALSE;
 730         while((fSave=(TFile*)it.Next()))
 731                 if(strcmp(fn, fSave->GetName())==0){
 732                         kFOUND = kTRUE;
 733                         break;
 734                 }
 735         if(!kFOUND) fSave = new TFile(fn, "RECREATE");
 736         fSave->cd();
 737
 738         TH1 *h;
 739         h = (TH1F*)h1MaxTB[0]->Clone(Form("h1MaxTBEL%02d", mom));
 740         h->SetTitle(Form("Maximum Time Bin distribution for electrons @ %4.1f GeV", fTrackMomentum[mom]));
 741         h->GetXaxis()->SetTitle("time [100 ns]");
 742         h->GetYaxis()->SetTitle("Probability");
 743         h->Write();
 744
 745         h = (TH1F*)h1MaxTB[1]->Clone(Form("h1MaxTBPI%02d", mom));
 746         h->SetTitle(Form("Maximum Time Bin distribution for pions @ %4.1f GeV", fTrackMomentum[mom]));
 747         h->GetXaxis()->SetTitle("time [100 ns]");
 748         h->GetYaxis()->SetTitle("Probability");
 749         h->Write();
 750
 751         pwd->cd();
 752 }
 753