]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/Cal/AliTRDCalPIDLQ.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
coding conventions
[u/mrichter/AliRoot.git] / TRD / Cal / AliTRDCalPIDLQ.cxx
CommitLineData
7754cd1f 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
9edefa04 16/* $Id$ */
2745a409 17
dab811d0 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//////////////////////////////////////////////////////////////////////////
7754cd1f 28
7754cd1f 29#include <TH1F.h>
dab811d0 30#include <TH2F.h>
7754cd1f 31#include <TFile.h>
dab811d0 32#include <TTree.h>
9edefa04 33#include <TROOT.h>
dab811d0 34#include <TPDGCode.h>
35#include <TParticle.h>
36#include <TPrincipal.h>
7754cd1f 37
b92cdef5 38#include "AliLog.h"
dab811d0 39#include "AliHeader.h"
40#include "AliGenEventHeader.h"
41#include "AliRun.h"
42#include "AliRunLoader.h"
43#include "AliStack.h"
2745a409 44#include "AliPID.h"
dab811d0 45#include "AliESD.h"
46#include "AliESDtrack.h"
b92cdef5 47
48#include "AliTRDCalPIDLQ.h"
dab811d0 49#include "AliTRDCalPIDLQRef.h"
50#include "AliTRDpidESD.h"
51#include "AliTRDcalibDB.h"
52#include "AliTRDtrack.h"
53#include "AliTRDgeometry.h"
b92cdef5 54
7754cd1f 55ClassImp(AliTRDCalPIDLQ)
56
57Char_t* AliTRDCalPIDLQ::fpartName[AliPID::kSPECIES] = {"electron", "muon", "pion", "kaon", "proton"};
dab811d0 58Char_t* AliTRDCalPIDLQ::fpartSymb[AliPID::kSPECIES] = {"EL", "MU", "PI", "KA", "PR"};
7754cd1f 59
60//_________________________________________________________________________
2745a409 61AliTRDCalPIDLQ::AliTRDCalPIDLQ()
dab811d0 62 :TNamed("pid", "PID for TRD")
b2bc6e56 63 ,fNMom(0)
dab811d0 64 ,fTrackMomentum(0x0)
b2bc6e56 65 ,fNLength(0)
dab811d0 66 ,fTrackSegLength(0x0)
67 ,fNTimeBins(0)
2745a409 68 ,fMeanChargeRatio(0)
69 ,fNbins(0)
70 ,fBinSize(0)
dab811d0 71 ,fHistdEdx(0x0)
72 ,fHistTimeBin(0x0)
73 ,fEstimator(0x0)
7754cd1f 74{
75 //
76 // The Default constructor
77 //
dab811d0 78
79 Init();
7754cd1f 80}
81
82//_________________________________________________________________________
2745a409 83AliTRDCalPIDLQ::AliTRDCalPIDLQ(const Text_t *name, const Text_t *title)
84 :TNamed(name,title)
b2bc6e56 85 ,fNMom(0)
dab811d0 86 ,fTrackMomentum(0x0)
b2bc6e56 87 ,fNLength(0)
dab811d0 88 ,fTrackSegLength(0x0)
89 ,fNTimeBins(0)
2745a409 90 ,fMeanChargeRatio(0)
91 ,fNbins(0)
92 ,fBinSize(0)
dab811d0 93 ,fHistdEdx(0x0)
94 ,fHistTimeBin(0x0)
95 ,fEstimator(0x0)
7754cd1f 96{
97 //
98 // The main constructor
99 //
100
101 Init();
b92cdef5 102
7754cd1f 103}
104
105//_____________________________________________________________________________
2745a409 106AliTRDCalPIDLQ::AliTRDCalPIDLQ(const AliTRDCalPIDLQ &c)
107 :TNamed(c)
b2bc6e56 108 ,fNMom(c.fNMom)
dab811d0 109 ,fTrackMomentum(0x0)
b2bc6e56 110 ,fNLength(c.fNLength)
dab811d0 111 ,fTrackSegLength(0x0)
112 ,fNTimeBins(c.fNTimeBins)
2745a409 113 ,fMeanChargeRatio(c.fMeanChargeRatio)
114 ,fNbins(c.fNbins)
115 ,fBinSize(c.fBinSize)
dab811d0 116 ,fHistdEdx(0x0)
117 ,fHistTimeBin(0x0)
118 ,fEstimator(0x0)
7754cd1f 119{
120 //
b92cdef5 121 // Copy constructor
7754cd1f 122 //
123
dab811d0 124 if (this != &c) ((AliTRDCalPIDLQ &) c).Copy(*this);
7754cd1f 125
7754cd1f 126}
127
128//_________________________________________________________________________
129AliTRDCalPIDLQ::~AliTRDCalPIDLQ()
130{
131 //
132 // Destructor
133 //
134
135 CleanUp();
b92cdef5 136
7754cd1f 137}
138
139//_________________________________________________________________________
140void AliTRDCalPIDLQ::CleanUp()
141{
b92cdef5 142 //
143 // Delets all newly created objects
144 //
145
2745a409 146 if (fHistdEdx) {
7754cd1f 147 delete fHistdEdx;
dab811d0 148 fHistdEdx = 0x0;
7754cd1f 149 }
150
2745a409 151 if (fHistTimeBin) {
7754cd1f 152 delete fHistTimeBin;
dab811d0 153 fHistTimeBin = 0x0;
7754cd1f 154 }
155
2745a409 156 if (fTrackMomentum) {
7754cd1f 157 delete[] fTrackMomentum;
dab811d0 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;
7754cd1f 169 }
b92cdef5 170
7754cd1f 171}
172
173//_____________________________________________________________________________
174AliTRDCalPIDLQ &AliTRDCalPIDLQ::operator=(const AliTRDCalPIDLQ &c)
175{
176 //
177 // Assignment operator
178 //
179
180 if (this != &c) ((AliTRDCalPIDLQ &) c).Copy(*this);
181 return *this;
b92cdef5 182
7754cd1f 183}
184
185//_____________________________________________________________________________
186void AliTRDCalPIDLQ::Copy(TObject &c) const
187{
188 //
189 // Copy function
190 //
191
192 AliTRDCalPIDLQ& target = (AliTRDCalPIDLQ &) c;
193
194 target.CleanUp();
195
2745a409 196 target.fNbins = fNbins;
197 target.fBinSize = fBinSize;
198 target.fMeanChargeRatio = fMeanChargeRatio;
dab811d0 199 target.fNTimeBins = fNTimeBins;
200
201 //target.fNMom = fNMom;
7754cd1f 202 target.fTrackMomentum = new Double_t[fNMom];
2745a409 203 for (Int_t i=0; i<fNMom; ++i) {
7754cd1f 204 target.fTrackMomentum[i] = fTrackMomentum[i];
2745a409 205 }
7754cd1f 206
dab811d0 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
2745a409 213 if (fHistdEdx) {
7754cd1f 214 target.fHistdEdx = (TObjArray*) fHistdEdx->Clone();
2745a409 215 }
216 if (fHistTimeBin) {
7754cd1f 217 target.fHistTimeBin = (TObjArray*) fHistTimeBin->Clone();
2745a409 218 }
219
dab811d0 220 target.fEstimator = new AliTRDCalPIDLQRef(*fEstimator);
221
7754cd1f 222 TObject::Copy(c);
b92cdef5 223
7754cd1f 224}
225
226//_________________________________________________________________________
227void AliTRDCalPIDLQ::Init()
228{
b92cdef5 229 //
230 // Initialization
231 //
232
dab811d0 233 fNMom = 11;
234 fNLength = 4;
b92cdef5 235
7754cd1f 236 fTrackMomentum = new Double_t[fNMom];
dab811d0 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*/);
7754cd1f 256 fHistdEdx->SetOwner();
dab811d0 257 fHistTimeBin = new TObjArray(2 * fNMom);
7754cd1f 258 fHistTimeBin->SetOwner();
dab811d0 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
7754cd1f 271 // ADC Gain normalization
b92cdef5 272 fMeanChargeRatio = 1.0;
7754cd1f 273
274 // Number of bins and bin size
b92cdef5 275 fNbins = 0;
276 fBinSize = 0.0;
7754cd1f 277}
278
279//_________________________________________________________________________
dab811d0 280Bool_t AliTRDCalPIDLQ::ReadReferences(Char_t *responseFile)
7754cd1f 281{
282 //
283 // Read the TRD dEdx histograms.
284 //
b92cdef5 285
7754cd1f 286 // Read histogram Root file
287 TFile *histFile = new TFile(responseFile, "READ");
288 if (!histFile || !histFile->IsOpen()) {
2745a409 289 AliError(Form("Opening TRD histgram file %s failed", responseFile));
7754cd1f 290 return kFALSE;
291 }
292 gROOT->cd();
293
294 // Read histograms
dab811d0 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 }
7754cd1f 308 }
309
310 histFile->Close();
311 delete histFile;
312
313 // Number of bins and bin size
dab811d0 314 //TH1F* hist = (TH1F*) fHistdEdx->At(GetHistID(AliPID::kPion, 1));
315 //fNbins = hist->GetNbinsX();
316 //fBinSize = hist->GetBinWidth(1);
7754cd1f 317
318 return kTRUE;
b92cdef5 319
7754cd1f 320}
321
dab811d0 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// }
7754cd1f 357
358//_________________________________________________________________________
dab811d0 359TH1* AliTRDCalPIDLQ::GetHistogram(Int_t k, Int_t ip/*, Int_t il*/) const
7754cd1f 360{
361 //
b92cdef5 362 // Returns one selected dEdx histogram
7754cd1f 363 //
b92cdef5 364
dab811d0 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]));
7754cd1f 369
dab811d0 370 return (TH1*)fHistdEdx->At(GetHistID(k, ip));
b92cdef5 371
372}
373
374//_________________________________________________________________________
dab811d0 375TH1* AliTRDCalPIDLQ::GetHistogramT(Int_t k, Int_t ip) const
b92cdef5 376{
377 //
378 // Returns one selected time bin max histogram
379 //
380
dab811d0 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]));
b92cdef5 385
dab811d0 386 return (TH1*)fHistTimeBin->At(((k==AliPID::kElectron)?0:1)*fNMom+ip);
7754cd1f 387}
388
389//_________________________________________________________________________
dab811d0 390Double_t AliTRDCalPIDLQ::GetProbability(Int_t spec, Double_t mom, Double_t *dedx, Double_t length) const
7754cd1f 391{
392 //
dab811d0 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
b92cdef5 397 //
7754cd1f 398
dab811d0 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);
b92cdef5 448
7754cd1f 449}
450
451//_________________________________________________________________________
dab811d0 452Double_t AliTRDCalPIDLQ::GetProbabilityT(Int_t spec, Double_t mom, Int_t timbin) const
7754cd1f 453{
454 //
455 // Gets the Probability of having timbin at a given momentum (mom)
dab811d0 456 // and particle type (spec) (0 for e) and (2 for pi)
7754cd1f 457 // from the precalculated timbin distributions
b92cdef5 458 //
7754cd1f 459
dab811d0 460 if (spec < 0 || spec >= AliPID::kSPECIES) return 0.;
461 if (timbin<0 || timbin >= fNTimeBins) return 0.;
2745a409 462
463 Int_t iTBin = timbin+1;
7754cd1f 464
2745a409 465 // Everything which is not an electron counts as a pion for time bin max
dab811d0 466 //if(spec != AliPID::kElectron) spec = AliPID::kPion;
467
7754cd1f 468
dab811d0 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}
2745a409 492
dab811d0 493//__________________________________________________________________
494Bool_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]);
7754cd1f 543 }
dab811d0 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}
b92cdef5 718
dab811d0 719//__________________________________________________________________
720void 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();
7754cd1f 752}
dab811d0 753