Place the config and root file at the right place
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEtrdPIDqa.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 // QA class for TRD PID
17 // Plot Pion Efficiency at x electron efficiency
18 // Calculate the threshold parametrisation and save
19 // them in a root file
20 //
21 // Author:
22 //   Markus Fasel <M.Fasel@gsi.de>
23 //
24 #include <TAxis.h>
25 #include <TBrowser.h>
26 #include <TClass.h>
27 #include <TCanvas.h>
28 #include <TF1.h>
29 #include <TFile.h>
30 #include <TGraph.h>
31 #include <TGraphErrors.h>
32 #include <THnSparse.h>
33 #include <TH1.h>
34 #include <TH2.h>
35 #include <TIterator.h>
36 #include <TLegend.h>
37 #include <TList.h>
38 #include <TMath.h>
39 #include <TObjArray.h>
40 #include <TString.h>
41
42 #include "AliAODTrack.h"
43 #include "AliAODPid.h"
44 #include "AliESDtrack.h"
45 #include "AliHFEtrdPIDqa.h"
46 #include "AliHFEtools.h"
47 #include "AliHFEpidTRD.h"
48 #include "AliLog.h"
49
50 ClassImp(AliHFEtrdPIDqa)
51
52 const Double_t AliHFEtrdPIDqa::fgkElectronEff[kNElectronEffs] = {0.7,0.75, 0.8, 0.85, 0.9, 0.95};
53 //_______________________________________________________________
54 // Definition of the common binning 
55 const Int_t AliHFEtrdPIDqa::fgkNBinsCommon[kQuantitiesCommon] = {
56   AliPID::kSPECIES + 1,         // species
57   40,                           // p-bins
58   AliESDtrack::kTRDnPlanes + 1 // tracklets including 0
59 };
60 const Double_t AliHFEtrdPIDqa::fgkMinBinCommon[kQuantitiesCommon] = {
61   -1,      // species
62   0.1,     // p-bins:
63   0        // tracklets including 0
64 };
65
66 const Double_t AliHFEtrdPIDqa::fgkMaxBinCommon[kQuantitiesCommon] = {
67   AliPID::kSPECIES,               // species
68   10.,                            // p-bins
69   AliESDtrack::kTRDnPlanes + 1    // tracklets including 0
70 };
71 //_______________________________________________________________
72
73 //__________________________________________________________________
74 AliHFEtrdPIDqa::AliHFEtrdPIDqa():
75   TNamed("trdPIDqa", ""),
76   fTRDpid(NULL),
77   fHistos(NULL),
78   fPionEfficiencies(NULL),
79   fProtonEfficiencies(NULL),
80   fKaonEfficiencies(NULL),
81   fThresholds(NULL),
82   fShowMessage(kFALSE),
83   fTotalChargeInSlice0(kFALSE),
84   fCentralityBin(-1)
85 {
86   //
87   // Default Constructor
88   //
89 }
90
91 //__________________________________________________________________
92 AliHFEtrdPIDqa::AliHFEtrdPIDqa(const Char_t *name):
93   TNamed(name, ""),
94   fTRDpid(NULL),
95   fHistos(NULL),
96   fPionEfficiencies(NULL),
97   fProtonEfficiencies(NULL),
98   fKaonEfficiencies(NULL),
99   fThresholds(NULL),
100   fShowMessage(kFALSE),
101   fTotalChargeInSlice0(kFALSE),
102   fCentralityBin(-1)
103 {
104   //
105   // Main Constructor
106   //
107 }
108
109 //__________________________________________________________________
110 AliHFEtrdPIDqa::AliHFEtrdPIDqa(const AliHFEtrdPIDqa &ref):
111   TNamed(ref),
112   fTRDpid(NULL),
113   fHistos(NULL),
114   fPionEfficiencies(NULL),
115   fProtonEfficiencies(NULL),
116   fKaonEfficiencies(NULL),
117   fThresholds(NULL),
118   fShowMessage(kFALSE),
119   fTotalChargeInSlice0(ref.fTotalChargeInSlice0),
120   fCentralityBin(ref.fCentralityBin)
121 {
122   //
123   // Copy constructor
124   //
125   ref.Copy(*this);
126 }
127
128 //__________________________________________________________________
129 AliHFEtrdPIDqa &AliHFEtrdPIDqa::operator=(const AliHFEtrdPIDqa &ref){
130   //
131   // Assignment operator
132   //
133   if(this != &ref)
134     ref.Copy(*this);
135   return *this;
136 }
137
138 //__________________________________________________________________
139 AliHFEtrdPIDqa::~AliHFEtrdPIDqa(){
140   //
141   // Destructor
142   //
143   if(fTRDpid) delete fTRDpid;
144   if(fHistos) delete fHistos;
145   if(fPionEfficiencies) delete fPionEfficiencies;
146   if(fProtonEfficiencies) delete fProtonEfficiencies;
147   if(fKaonEfficiencies) delete fKaonEfficiencies;
148 }
149
150 //__________________________________________________________________
151 void AliHFEtrdPIDqa::Copy(TObject &ref) const{
152   //
153   // Copies content of this object into object ref
154   //
155   TNamed::Copy(ref);
156
157   AliHFEtrdPIDqa &target = dynamic_cast<AliHFEtrdPIDqa &>(ref);
158   target.fTRDpid = fTRDpid;
159   target.fHistos = dynamic_cast<AliHFEcollection *>(fHistos->Clone());  
160   target.fTotalChargeInSlice0 = fTotalChargeInSlice0;
161   target.fCentralityBin = fCentralityBin;
162 }
163
164 //__________________________________________________________________
165 Long64_t AliHFEtrdPIDqa::Merge(TCollection *coll){
166   //
167   // Merge objects
168   //
169   if(!coll) return 0;
170   if(coll->IsEmpty()) return 1;
171   
172   AliHFEtrdPIDqa *refQA = NULL;
173   TIter it(coll);
174   TObject *o = NULL;
175   Long64_t count = 0;
176   TList listHistos;
177   while((o = it())){
178     refQA = dynamic_cast<AliHFEtrdPIDqa *>(o);
179     if(!refQA) continue;
180
181     listHistos.Add(refQA->fHistos);
182     count++; 
183   }
184   fHistos->Merge(&listHistos);
185   return count+1;
186 }
187
188 //__________________________________________________________________
189 void AliHFEtrdPIDqa::Browse(TBrowser *b){
190   //
191   // Enable Browser functionality
192   //
193   if(b){
194     // Add objects to the browser
195     if(fHistos) b->Add(fHistos, fHistos->GetName());
196     if(fPionEfficiencies) b->Add(fPionEfficiencies, "Pion Efficiencies");
197     if(fProtonEfficiencies) b->Add(fProtonEfficiencies, "Proton Efficiencies");  
198     if(fKaonEfficiencies) b->Add(fKaonEfficiencies, "Kaon Efficiencies");
199     if(fThresholds) b->Add(fThresholds, "Thresholds");
200   }
201 }
202
203 //__________________________________________________________________
204 void AliHFEtrdPIDqa::Init(){
205   //
206   // Initialize Object
207   //
208   
209   fHistos = new AliHFEcollection("TRDqa", "Histos for TRD QA");
210
211   CreateLikelihoodHistogram();
212   CreateQAHistogram();
213   CreatedEdxHistogram();
214   CreateHistoTruncatedMean();
215
216   fTRDpid = new AliHFEpidTRD("QAtrdPID");
217   if(fTotalChargeInSlice0) fTRDpid->SetTotalChargeInSlice0();
218 }
219
220 //__________________________________________________________________
221 void AliHFEtrdPIDqa::CreateLikelihoodHistogram(){
222   //
223   // Create Histogram for TRD Likelihood Studies
224   //
225   Int_t nbins[kQuantitiesLike]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
226   nbins[kElectronLike] = 100;
227   nbins[kNClustersLike] = 200;
228   nbins[kCentralityBinLike] = 12;
229   Double_t binMin[kQuantitiesLike]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
230   Double_t binMax[kQuantitiesLike]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
231   binMin[kElectronLike] = 0.;      
232   binMin[kNClustersLike] = 0.;
233   binMin[kCentralityBinLike] = -1.;
234   binMax[kElectronLike] = 1.;
235   binMax[kNClustersLike] = 200.;
236   binMax[kCentralityBinLike] = 11.;
237
238   fHistos->CreateTHnSparse("fLikeTRD","TRD Likelihood Studies", kQuantitiesLike, nbins, binMin, binMax);
239   fHistos->BinLogAxis("fLikeTRD", kP);
240 }
241
242 //__________________________________________________________________
243 void AliHFEtrdPIDqa::CreateQAHistogram(){
244   //
245   // Create Histogram for Basic TRD PID QA
246   //
247   AliDebug(1, "Called");
248   Int_t nbins[kQuantitiesQA]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
249   nbins[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1;
250   nbins[kNClusters] = 200;
251   nbins[kCentralityBinQA] = 12;
252   Double_t binMin[kQuantitiesQA]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
253   binMin[kNonZeroTrackletCharge] = 0.;      
254   binMin[kNClusters] = 0.;
255   binMin[kCentralityBinQA] = -1.;
256   Double_t binMax[kQuantitiesQA]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
257   binMax[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1.;
258   binMax[kNClusters] = 200.;
259   binMax[kCentralityBinQA] = 11.;
260   fHistos->CreateTHnSparse("fQAtrack","TRD QA Histogram", kQuantitiesQA, nbins, binMin, binMax);
261   fHistos->BinLogAxis("fQAtrack", kP);
262 }
263
264 //__________________________________________________________________
265 void AliHFEtrdPIDqa::CreatedEdxHistogram(){
266   //
267   // Create QA histogram for dEdx investigations
268   //
269   AliDebug(1, "Called");
270   Int_t nbins[kQuantitiesdEdx]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
271   nbins[kdEdx] = 100;
272   nbins[kNclusters] = 261;
273   nbins[kNonZeroSlices] = 9;
274   nbins[kCentralityBindEdx] = 12;
275   Double_t binMin[kQuantitiesdEdx]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
276   binMin[kdEdx] = 0.;     
277   binMin[kNclusters] = 0;
278   binMin[kNonZeroSlices] = 0.;
279   binMin[kCentralityBindEdx] = -1.;
280   Double_t binMax[kQuantitiesdEdx]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
281   binMax[kdEdx] = 10000.;
282   binMax[kNclusters] = 260.;
283   binMax[kNonZeroSlices] = 8.;
284   binMax[kCentralityBindEdx] = 11.;
285
286   fHistos->CreateTHnSparse("fQAdEdx","TRD summed dEdx", kQuantitiesdEdx, nbins, binMin, binMax);
287   fHistos->BinLogAxis("fQAdEdx", kP);
288   fHistos->Sumw2("fQAdEdx");
289 }
290
291 //__________________________________________________________________
292 void AliHFEtrdPIDqa::CreateHistoTruncatedMean(){
293   //
294   // Create Histogram for Basic TRD PID QA:
295   //
296   AliDebug(1, "Called");
297   Int_t nbins[kQuantitiesTruncMean]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
298   nbins[kTPCdEdx] = 600;
299   nbins[kTRDdEdxMethod1] = 1000;
300   nbins[kTRDdEdxMethod2] = 1000;
301   nbins[kCentralityBinTruncMean] = 12;
302   Double_t binMin[kQuantitiesTruncMean]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
303   binMin[kTPCdEdx] = 0.;      
304   binMin[kTRDdEdxMethod1] = 0.;      
305   binMin[kTRDdEdxMethod2] = 0.;
306   binMin[kCentralityBinTruncMean] = -1.;
307   Double_t binMax[kQuantitiesTruncMean]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
308   binMax[kTPCdEdx] = 600;
309   binMax[kTRDdEdxMethod1] = 20000.;
310   binMax[kTRDdEdxMethod2] = 20000.;
311   binMax[kCentralityBinTruncMean] = 11.;
312
313   fHistos->CreateTHnSparse("fTRDtruncMean","TRD TruncatedMean studies", kQuantitiesTruncMean, nbins, binMin, binMax);
314   fHistos->BinLogAxis("fTRDtruncMean", kP);
315   fHistos->CreateTH2F("fTRDslicesPions","TRD dEdx per slice for Pions", 8, 0, 8, 2000, 0, 8000);
316   fHistos->CreateTH2F("fTRDslicesElectrons","TRD dEdx per slice for Electrons", 8, 0, 8, 2000, 0, 8000);
317 }
318
319
320 //__________________________________________________________________
321 void AliHFEtrdPIDqa::ProcessTracks(const TObjArray * const tracks, Int_t species){
322   //
323   // Process Electron Tracks
324   //
325   if(species < -1 || species >= AliPID::kSPECIES) return;
326   TIter it(tracks);
327   const AliVTrack *track = NULL;
328   while((track = dynamic_cast<const AliVTrack *>(it()))){
329     if(track) ProcessTrack(track, species);
330   }
331 }
332
333 //__________________________________________________________________
334 void AliHFEtrdPIDqa::ProcessTrack(const AliVTrack * const track, Int_t species){
335   //
336   // Process Single Track
337   //
338   if(TString(track->IsA()->GetName()).CompareTo("AliESDtrack") == 0)
339     ProcessTrackESD(dynamic_cast<const AliESDtrack *>(track), species);
340   else if(TString(track->IsA()->GetName()).CompareTo("AliAODTrack") == 0)
341     ProcessTrackAOD(dynamic_cast<const AliAODTrack *>(track), species);
342 }
343
344
345 //__________________________________________________________________
346 void AliHFEtrdPIDqa::ProcessTrackESD(const AliESDtrack *track, Int_t species){
347   //
348   // Process single ESD track
349   //
350   if(!track) return;
351   if((track->GetStatus() & AliESDtrack::kTRDout) == 0) return;  // require TRD track
352   FillTRDLikelihoods(track, species);
353   FillTRDQAplots(track, species);
354 }
355
356 //__________________________________________________________________
357 void AliHFEtrdPIDqa::ProcessTrackAOD(const AliAODTrack * const track, Int_t /*species*/){
358   //
359   // Process single AOD track
360   // AOD PID object required
361   //
362   if(!track) return;
363   AliAODPid *trackPID = track->GetDetPid();
364   if(!trackPID) return;
365
366 }
367
368 //__________________________________________________________________
369 void AliHFEtrdPIDqa::FillTRDLikelihoods(const AliESDtrack * const track, Int_t species){
370   //
371   // Fill TRD Likelihood Histogram
372   //
373   Double_t trdLike[AliPID::kSPECIES];
374   track->GetTRDpid(trdLike);
375   // Renormalize 
376   Double_t norm =trdLike[AliPID::kElectron]+trdLike[AliPID::kPion];
377   Double_t likeEle = norm == 0. ? 0. : trdLike[AliPID::kElectron]/norm;
378   const AliExternalTrackParam *outerPars = track->GetOuterParam();
379
380   //Int_t kQuantitiesLike;
381   Double_t quantities[kQuantitiesLike]; memset(quantities, 0, sizeof(Double_t) * kQuantitiesLike);
382   // we store:
383   // species
384   // p
385   // ntracklets
386   // Electron Likelihood
387   // Centrality
388   quantities[kSpecies] = species;
389   quantities[kP] = outerPars ? outerPars->P() : track->P();
390   quantities[kNTracklets] = track->GetTRDntrackletsPID();
391   quantities[kElectronLike] = likeEle;
392   quantities[kNClustersLike] =  track->GetTRDncls();
393   quantities[kCentralityBinLike] = fCentralityBin;
394   fHistos->Fill("fLikeTRD", quantities);
395
396 }
397
398 //__________________________________________________________________
399 void AliHFEtrdPIDqa::FillTRDQAplots(const AliESDtrack * const track, Int_t species){
400   //
401   // Fill QA Plots containing further information
402   //
403   const AliExternalTrackParam *outerPars = track->GetOuterParam();
404
405   Double_t quantitiesQA[kQuantitiesQA], quantitiesdEdx[kQuantitiesdEdx], quantitiesTruncMean[kQuantitiesTruncMean];
406   // we store:
407   // 1. QA
408   // species
409   // p
410   // ntracklets
411   // Non-zero tracklet charges
412   // Number of clusters / full track
413   // 2. dEdx
414   // species
415   // p
416   // ntracklets
417   // dEdx
418   // 3. Truncated Mean
419   // ...
420   // TPC dEdx
421   // TRD dEdx Method 1
422   // TRD dEdx Method 2
423   quantitiesQA[kSpecies]  = quantitiesdEdx[kSpecies] 
424                           = quantitiesTruncMean[kSpecies] 
425                           = species;
426   quantitiesQA[kP]  = quantitiesTruncMean[kP]
427                     = outerPars ? outerPars->P() : track->P();
428   quantitiesQA[kNTracklets] = quantitiesdEdx[kNTracklets] 
429                             = quantitiesTruncMean[kNTracklets]
430                             = track->GetTRDntrackletsPID();
431   quantitiesQA[kNClusters] = quantitiesdEdx[kNclusters] = track->GetTRDncls();
432   
433
434   Double_t dEdxSum = 0., qSlice = 0.;
435   // remove the last slice from the histogram
436   Int_t ntrackletsNonZero = 0, nSlices = track->GetNumberOfTRDslices(), nSlicesNonZero = 0;
437   TString speciesname = "pions";
438   Bool_t selectedForSlicemon = kFALSE;
439   
440   switch(species){
441     case AliPID::kElectron: speciesname = "Electrons"; selectedForSlicemon = kTRUE; break;
442     case AliPID::kPion: speciesname = "Pions"; selectedForSlicemon = kTRUE; break;
443     default: speciesname = "undefined"; selectedForSlicemon = kFALSE; break;
444   };
445   AliDebug(1, Form("species %d, speciesname %s, momentum %f, selected %s", species, speciesname.Data(), track->P(), selectedForSlicemon ? "yes" : "no"));
446   for(Int_t iplane = 0; iplane < AliESDtrack::kTRDnPlanes; iplane++){
447     quantitiesdEdx[kP] = track->GetTRDmomentum(iplane);
448     dEdxSum = 0.;
449     for(Int_t islice = 0; islice < nSlices; islice++){
450       if(fTotalChargeInSlice0 && islice >= 7) break;
451       qSlice = track->GetTRDslice(iplane, fTotalChargeInSlice0 ? islice + 1 : islice);  // hack by mfasel: For data with the new reconstruction, slice 0 is used to store the total charge, the total number of slices is 7 instead of 8
452       if(qSlice > 1e-1){
453         // cut out 0 slices
454         nSlicesNonZero++;
455         dEdxSum += qSlice;
456         // Reweighting of the slices for the truncated mean: select all pion tracks above
457         // 1.5 GeV and monitor the dEdx as function of slice
458         if(selectedForSlicemon && track->P() > 1.5){
459           AliDebug(2, Form("Filling Histogram fTRDslices%s", speciesname.Data()));
460           fHistos->Fill(Form("fTRDslices%s", speciesname.Data()), static_cast<Double_t>(islice), qSlice);
461         }
462       }
463     }
464     quantitiesdEdx[kNonZeroSlices] = nSlicesNonZero;
465     quantitiesdEdx[kdEdx] = fTotalChargeInSlice0 ? track->GetTRDslice(iplane, 0) : dEdxSum; // hack by mfasel: In the new reconstruction, the total charge is stored in the first slice, in the old reconstruction it has to be calculated from the slice charges.     
466     if(dEdxSum) ntrackletsNonZero++;
467     // Fill dEdx histogram
468     quantitiesdEdx[kCentralityBindEdx] = fCentralityBin;
469     if(dEdxSum > 1e-1) fHistos->Fill("fQAdEdx", quantitiesdEdx); // Cut out 0 entries
470   }
471   quantitiesQA[kNonZeroTrackletCharge] = ntrackletsNonZero;
472   quantitiesQA[kCentralityBinQA] = fCentralityBin;
473   fHistos->Fill("fQAtrack", quantitiesQA);
474
475   quantitiesTruncMean[kTPCdEdx] = track->GetTPCsignal();
476   quantitiesTruncMean[kTRDdEdxMethod1] = fTRDpid->GetTRDSignalV1(track, 0.6);
477   quantitiesTruncMean[kTRDdEdxMethod2] = fTRDpid->GetTRDSignalV2(track, 0.6);
478   quantitiesTruncMean[kCentralityBinTruncMean] = fCentralityBin;
479   fHistos->Fill("fTRDtruncMean", quantitiesTruncMean);
480 }
481
482 /////////////////////////////////////////////////////////
483 //
484 // Code for Post Processing
485 //
486 // //////////////////////////////////////////////////////
487
488 //__________________________________________________________________
489 void AliHFEtrdPIDqa::FinishAnalysis(Int_t fCentrality, Bool_t isGreaterEqual){
490   //
491   // Finish Analysis:
492   // Calculate Electron Efficiency for ntracklets = 4...6
493   // Calculate thresholds for ntracklets = 4...6
494   //
495   
496   if(!fPionEfficiencies){ 
497     fPionEfficiencies = new TList;
498     fPionEfficiencies->SetName("pionEfficiencies");
499     fPionEfficiencies->SetOwner();
500   }
501   if(!fProtonEfficiencies){
502     fProtonEfficiencies = new TList;
503     fProtonEfficiencies->SetName("protonEfficiencies");
504     fProtonEfficiencies->SetOwner();
505   }
506   if(!fThresholds){
507     fThresholds = new TList;
508     fThresholds->SetName("thresholds");
509     fThresholds->SetOwner();
510   }
511
512   for(Int_t itr = 4; itr <= 6; itr++){
513     if(fShowMessage){
514       printf("========================================\n");
515       printf("Analysing %d trackltes centrality %i \n", itr, fCentrality);
516       printf("========================================\n");
517     }
518     AnalyseNTracklets(itr, fCentrality, isGreaterEqual);
519   }
520 }
521
522 //__________________________________________________________________
523 void AliHFEtrdPIDqa::StoreResults(const Char_t *filename){
524   //
525   // Store histos into a root file
526   //
527   TFile *outfile = new TFile(filename, "RECREATE");
528   outfile->cd();
529   fPionEfficiencies->Clone()->Write(fPionEfficiencies->GetName(), kSingleKey);
530   fProtonEfficiencies->Clone()->Write(fProtonEfficiencies->GetName(), kSingleKey);
531   fThresholds->Clone()->Write(fThresholds->GetName(), kSingleKey);
532   outfile->Close();
533   delete outfile;
534 }
535
536 //__________________________________________________________________
537 void AliHFEtrdPIDqa::SaveThresholdParameters(const Char_t *name, Double_t lowerLimit, Double_t upperLimit, Int_t icentrality){
538   //
539   // Fit the threshold histograms with the given parametrisation
540   // and store the TF1 in the file
541   //
542
543   if(!fThresholds){
544     AliError("Threshold histograms have to be created first");
545     return;
546   }
547
548     if(fShowMessage){
549     printf("========================================\n");
550     printf("Calculating threshold parameters\n");
551     printf("========================================\n");
552   }
553
554   TList *outlist = new TList;
555   outlist->SetName("thresholdTRD");
556   
557   TGraph *threshhist = NULL;
558   TF1 *threshparam = NULL;
559   TList *lHistos = NULL, *lFormulas = NULL;
560   for(Int_t itracklet = 4; itracklet <= 6; itracklet++){
561   
562     if(fShowMessage){
563       printf("-------------------------------\n");
564       printf("Processing %d tracklets\n", itracklet);
565       printf("-------------------------------\n");
566     }
567
568     Char_t *listname=Form("%dTracklets", itracklet);
569     if(icentrality!=-1) listname=Form("%dTracklets%dCentrality", itracklet, icentrality);
570     lHistos = dynamic_cast<TList *>(fThresholds->FindObject(listname));
571     if(!lHistos){
572       AliError(Form("Threshold histograms for the case %s not found", listname));
573       continue;
574     }
575     lFormulas = new TList;
576     lFormulas->SetName(listname);
577     outlist->Add(lFormulas);
578     
579     for(Int_t ieff = 0; ieff <  kNElectronEffs; ieff++){
580       
581       if(fShowMessage){
582         printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n");
583         printf("Processing Electron Efficiency %f\n", fgkElectronEff[ieff]);
584         printf(">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n");
585       }
586
587       threshhist = dynamic_cast<TGraph *>(lHistos->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
588       if(!threshhist) continue;
589       threshparam = MakeThresholds(threshhist, lowerLimit, upperLimit);
590       threshparam->SetName(Form("thresh_%d_%d", itracklet, static_cast<Int_t>(fgkElectronEff[ieff] * 100)));
591       if(icentrality!=-1) threshparam->SetName(Form("thresh_%d_%d_%d", itracklet, icentrality, static_cast<Int_t>(fgkElectronEff[ieff] * 100)));
592       lFormulas->Add(threshparam);
593     }
594   }
595
596   // store the output
597   TFile *outfile = new TFile(name, "RECREATE");
598   outfile->cd();
599   outlist->Write(outlist->GetName(), kSingleKey);
600   outfile->Close();
601   delete outfile;
602 }
603
604 //__________________________________________________________________
605 void AliHFEtrdPIDqa::AnalyseNTracklets(Int_t nTracklets, Int_t nCentrality, Bool_t isGreaterEqual){
606   //
607   // Calculate Pion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
608   // elPion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
609   // electron efficiencies
610   //
611   THnSparse *hLikeTRD = dynamic_cast<THnSparseF *>(fHistos->Get("fLikeTRD"));
612   if(!hLikeTRD){
613     AliError("Likelihood Histogram not available");
614     return;
615   }
616
617   Bool_t isPbPb = kFALSE;
618   if(nCentrality==-1)  isPbPb = kFALSE;
619   if(nCentrality!=-1)  isPbPb = kTRUE;
620
621   Int_t binTracklets = hLikeTRD->GetAxis(kNTracklets)->FindBin(nTracklets);
622   hLikeTRD->GetAxis(kNTracklets)->SetRange(binTracklets, isGreaterEqual ? 7 : binTracklets);
623
624   if(isPbPb){
625       Int_t binCentrality = hLikeTRD->GetAxis(kCentralityBinLike)->FindBin(nCentrality);
626       hLikeTRD->GetAxis(kCentralityBinLike)->SetRange(binCentrality, isGreaterEqual ? 11 : binCentrality);
627       /*
628        new TCanvas;
629        TH2 *test = hLikeTRD->Projection(kCentralityBin, kP);
630        test->Draw("colz");
631        */
632   }
633
634   Int_t binElectrons = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kElectron); 
635   AliDebug(1, Form("BinElectrons %d", binElectrons));
636   Int_t binPions = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kPion);
637   AliDebug(1, Form("BinPions %d", binPions));
638   Int_t binProtons =  hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kProton);
639   AliDebug(1, Form("BinProtons %d", binProtons));
640   hLikeTRD->GetAxis(kSpecies)->SetRange(binElectrons, binElectrons);
641   TH2 *likeElectron = hLikeTRD->Projection(kElectronLike, kP);
642   likeElectron->SetName("likeElectron");
643   hLikeTRD->GetAxis(kSpecies)->SetRange(binPions, binPions);
644   TH2 *likePion = hLikeTRD->Projection(kElectronLike, kP);
645   likePion->SetName("likePion");
646   hLikeTRD->GetAxis(kSpecies)->SetRange(binProtons, binProtons);
647   TH2 *likeProton = hLikeTRD->Projection(kElectronLike, kP);
648   likeProton->SetName("likeProton");
649   
650   // Undo ranges
651   hLikeTRD->GetAxis(kSpecies)->SetRange(0, hLikeTRD->GetAxis(kSpecies)->GetNbins());
652   hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kNTracklets)->GetNbins());
653   hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kCentralityBinLike)->GetNbins());
654
655   // Prepare List for output
656   Char_t *listname=Form("%dTracklets", nTracklets);
657   if(isPbPb) listname=Form("%dTracklets%dCentrality", nTracklets, nCentrality);
658
659
660   TList *listPions = new TList; listPions->SetName(listname); listPions->SetOwner();
661   TList *listProtons = new TList; listProtons->SetName(listname); listProtons->SetOwner();
662   TList *listThresholds = new TList; listThresholds->SetName(listname); listThresholds->SetOwner();
663   fPionEfficiencies->Add(listPions);
664   fProtonEfficiencies->Add(listProtons);
665   fThresholds->Add(listThresholds);
666
667   TH1 *probsEl = NULL, *probsPi = NULL, *probsPr = NULL;
668   TGraphErrors *effPi = NULL, *effPr = NULL, *thresholds = NULL;
669   Double_t p = 0, dp = 0;
670   Int_t threshbin = 0;
671   Double_t eff, error; // value and error
672   for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
673     
674     if(fShowMessage){
675       printf("-----------------------------------------\n");
676       printf("Doing Electron Efficiency %f\n", fgkElectronEff[ieff]);
677       printf("-----------------------------------------\n");
678     }
679     effPi = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
680     effPi->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
681     effPr = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
682     effPr->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
683     thresholds = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
684     thresholds->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
685   
686     // Add to lists
687     listPions->Add(effPi);
688     listProtons->Add(effPr);
689     listThresholds->Add(thresholds);
690
691     for(Int_t imom = 1; imom <= likeElectron->GetXaxis()->GetLast(); imom++){
692       p = likeElectron->GetXaxis()->GetBinCenter(imom);
693       dp = likeElectron->GetXaxis()->GetBinWidth(imom)/2;
694
695       probsEl = likeElectron->ProjectionY("el", imom, imom);
696       if(!probsEl->GetEntries()) continue;
697       probsPi = likePion->ProjectionY("pi", imom, imom);
698       if(!probsPi->GetEntries()) continue;
699       probsPr = likeProton->ProjectionY("pr", imom, imom);
700       if(!probsPr->GetEntries()) continue;
701       AliDebug(1, Form("Calculating Values for p = %f", p));
702
703       // Calculate non-electronEfficiency and error
704       eff = CalculateHadronEfficiency(probsPi, probsEl, fgkElectronEff[ieff], threshbin, error);
705       thresholds->SetPoint(imom - 1, p, probsEl->GetXaxis()->GetBinCenter(threshbin));
706       thresholds->SetPointError(imom - 1, dp, EstimateThresholdError(probsEl, threshbin));
707       AliDebug(1, Form("threshold %d|%f", threshbin, probsEl->GetXaxis()->GetBinCenter(threshbin)));
708       AliDebug(1, Form("Pion Efficiency %f +- %f", eff, error));
709       effPi->SetPoint(imom - 1, p, eff);
710       effPi->SetPointError(imom - 1, dp, error);
711       eff = CalculateHadronEfficiency(probsPr, probsEl, fgkElectronEff[ieff] , threshbin, error);
712       AliDebug(1, Form("Proton Efficiency %f", eff));
713       effPr->SetPoint(imom - 1, p, eff);
714       effPr->SetPointError(imom - 1, dp, error);
715  
716       // cleanup
717       delete probsEl;
718       delete probsPi;
719       delete probsPr;
720     }
721   }
722
723   // remove temporary histograms
724   delete likeElectron;
725   delete likePion;
726   delete likeProton;
727 }
728
729 //__________________________________________________________________
730 Double_t AliHFEtrdPIDqa::CalculateHadronEfficiency(const TH1 * const hadron, const TH1 *const electron, Double_t eff, Int_t &threshbin, Double_t &error){
731   // 
732   // Calculate non-electron efficiency
733   // optionally returns sums as second parameter
734   //
735
736   TH1D eleWorking(*((const TH1D *)electron)), hadronWorking(*((const TH1D *)hadron)); // Leave the original histograms untouched and do calculation including scale on copy
737   eleWorking.Scale(1./eleWorking.Integral());
738   hadronWorking.Scale(1./hadronWorking.Integral());
739
740   TArrayD sumsEl(eleWorking.GetNbinsX()), sumsHd(eleWorking.GetNbinsX());
741
742   // calculate threshold and estimated electron efficiency the threshold was taken
743   Double_t elEff = 0.;  // estimated electron efficiency at the end
744   Int_t currentBin = 0, nbins = 0;
745   for(Int_t ibin = eleWorking.GetXaxis()->GetLast(); ibin >= eleWorking.GetXaxis()->GetFirst(); ibin--){
746     currentBin = ibin;
747     nbins++;
748     elEff += eleWorking.GetBinContent(ibin);
749     sumsEl[eleWorking.GetXaxis()->GetLast() - ibin] = elEff;
750     if(elEff >= eff){
751       // we found the matching bin, break the loop
752       break;
753     }
754   }
755   threshbin = currentBin;
756
757   Double_t hdEff = 0; 
758   for(Int_t ibin = hadronWorking.GetXaxis()->GetLast(); ibin >= threshbin; ibin--) {
759     hdEff += hadronWorking.GetBinContent(ibin);
760     sumsHd[hadronWorking.GetXaxis()->GetLast() - ibin] = hdEff;
761   }
762
763   // search sums of electron efficiency for double counts, eliminate in electron and hadron array
764   TArrayD newsumsEl(100), newsumsHd(100);
765   Int_t nusable = 0;
766   for(Int_t ien = 0; ien < nbins; ien++){
767     if(ien==0){
768       newsumsEl[0] = sumsEl[0];
769       nusable++;
770       continue;
771     }
772     Int_t index = TMath::BinarySearch(nusable, newsumsEl.GetArray(), sumsEl[ien]);
773     if(TMath::Abs(sumsEl[ien] - newsumsEl[index]) < 1e-13){
774       // element already counted, don't add to the new arrays
775       continue; 
776     }
777     newsumsEl[nusable] = sumsEl[ien];
778     newsumsHd[nusable] = sumsHd[ien];
779     nusable++;
780   }
781
782   //printf("New array\n");
783   //for(Int_t ib = 0; ib < nusable; ib++){
784   //  printf("Electron Efficiency %f, Pion Efficiency %f\n", newsumsEl[ib], newsumsHd[ib]);
785   //}
786   //printf("Do Fit\n");
787
788   // Calculate error
789   error = 0;
790   if(hadronWorking.GetEntries() > 0 && eleWorking.GetEntries() > 0 && nusable > 2){
791     // Do error calculation in case the bins have enough statistics
792     TGraph gevh(nusable, newsumsEl.GetArray(), newsumsHd.GetArray()); 
793     TF1 evh("evh","pol2", eff-.05, eff+.05);
794     gevh.Fit(&evh, "Q", "", eff-.05, eff+.05);
795   
796     // return the error of the pion efficiency
797     if(((1.-hdEff) < 0) || ((1.- elEff) < 0)){
798       AliError(" ElEffi or HdEffi > 1. Error can not be calculated. Please increase statistics!");
799     }   else {
800       error = TMath::Sqrt(hdEff*(1-hdEff)/hadronWorking.GetEntries()+TMath::Power(evh.Derivative(eff), 2)*elEff*(1-elEff)/eleWorking.GetEntries());
801     }
802     AliDebug(2, Form("Pion Effi at [%f] : [%f +/- %f], Threshold[%f]", elEff, hdEff, error, eleWorking.GetBinCenter(threshbin)));
803     AliDebug(2, Form("Derivative at %4.2f : %f\n", eff, evh.Derivative(eff)));
804   }
805
806   return hdEff;
807 }
808
809 //__________________________________________________________________
810 Double_t AliHFEtrdPIDqa::EstimateThresholdError(const TH1 * const electron, Int_t threshbin){
811   //
812   // Estimate threshold error as sqrt(N_int)/N_ent
813   // where N_int are the counts integrated under the electron selection band
814   // and N_ent are all entries in the histogram
815   //
816   Int_t nInt = 0;
817   for(Int_t ibin = threshbin; ibin <= electron->GetXaxis()->GetNbins(); ibin++){
818     nInt += (Int_t) electron->GetBinContent(ibin);
819   }
820   Int_t nEnt = (Int_t) electron->Integral();
821   if(nEnt == 0) return 1.;
822   return TMath::Sqrt(static_cast<Double_t>(nInt))/static_cast<Double_t>(nEnt);
823 }
824
825 //__________________________________________________________________
826 Double_t AliHFEtrdPIDqa::CalculateIntegratedPionEfficiency(UInt_t nTracklets, Double_t electronEff, Double_t pmin, Double_t pmax, Int_t icentrality, Double_t *error){
827   //
828   // Calculate Pion Efficiency for a given electron efficiency in the specified momentum range
829   //
830   if(nTracklets < 4 || nTracklets > 6){
831     AliError("Pion Efficiency calculation only available for 4, 5, and 6 tracklets");
832     return 0.;
833   }
834   if(electronEff < 0.6 || electronEff > 1.){
835     AliError("Pion Efficiency calculation only available in the electron efficiency range 0.6 to 1");
836     return 0.;
837   }
838   if(pmin < 0.1 || pmin > 10 || pmax < 0.1 || pmax > 10.){
839     AliError("Pion Efficiency calculation only available in the momentum range 0.1 to 10 GeV/c");
840     return 0.;
841   }
842   if(pmax < pmin){
843     AliError("pmin is expected to be >= pmax");
844     return 0.;
845   }
846
847   // prerequierements fullfiled
848   // prepare histos
849   THnSparse *hLikeTRD = dynamic_cast<THnSparseF *>(fHistos->Get("fLikeTRD"));
850   if(!hLikeTRD){
851     AliError("Likelihood Histogram not available");
852     return 0;
853   }
854   Int_t binTracklets = hLikeTRD->GetAxis(kNTracklets)->FindBin(nTracklets);
855   hLikeTRD->GetAxis(kNTracklets)->SetRange(binTracklets, binTracklets);
856
857   if(icentrality!=1){
858       Int_t binCentrality = hLikeTRD->GetAxis(kCentralityBinLike)->FindBin(icentrality);
859       hLikeTRD->GetAxis(kCentralityBinLike)->SetRange(binCentrality, binCentrality);
860   }
861
862   Int_t pbinMin = hLikeTRD->GetAxis(kP)->FindBin(pmax),
863         pbinMax = hLikeTRD->GetAxis(kP)->FindBin(pmax);
864   hLikeTRD->GetAxis(kP)->SetRange(pbinMin, pbinMax);
865   Int_t binElectrons = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kElectron); 
866   Int_t binPions = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kPion);
867   hLikeTRD->GetAxis(kSpecies)->SetRange(binElectrons, binElectrons);
868   TH1 *likeElectron = hLikeTRD->Projection(kElectronLike);
869   likeElectron->Scale(1./likeElectron->Integral());
870   likeElectron->SetName("likeElectron");
871   hLikeTRD->GetAxis(kSpecies)->SetRange(binPions, binPions);
872   TH1 *likePion = hLikeTRD->Projection(kElectronLike);
873   likePion->Scale(1./likePion->Integral());
874   likePion->SetName("likePion");
875   
876   // Undo ranges
877   hLikeTRD->GetAxis(kSpecies)->SetRange(0, hLikeTRD->GetAxis(kSpecies)->GetNbins());
878   hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kNTracklets)->GetNbins());
879   hLikeTRD->GetAxis(kP)->SetRange(0, hLikeTRD->GetAxis(kP)->GetNbins());
880
881   // Do Calculation
882   Int_t thresh; Double_t err;
883   Double_t effpi = CalculateHadronEfficiency(likePion, likeElectron, electronEff, thresh, err);
884   delete likePion; delete likeElectron;
885   if(error) *error = err;
886   return effpi;
887 }
888
889 //__________________________________________________________________
890 void AliHFEtrdPIDqa::DrawTracklet(Int_t itracklet, Double_t pmin, Double_t pmax, Int_t icentrality, Bool_t doFit){
891   //
892   // Draw efficiencies and threshold as function of p
893   //
894   if(!(fPionEfficiencies && fProtonEfficiencies && fThresholds)){
895     AliError("No graphs to draw available");  
896     return;
897   }
898
899   Char_t *listname=Form("%dTracklets", itracklet);
900   if(icentrality!=-1) listname=Form("%dTracklets%dCentrality", itracklet, icentrality);
901
902   TList *lpions = dynamic_cast<TList *>(fPionEfficiencies->FindObject(listname));
903   TList *lprotons = dynamic_cast<TList *>(fProtonEfficiencies->FindObject(listname));
904   
905   TList *lthresholds = dynamic_cast<TList *>(fThresholds->FindObject(listname));
906   if(!(lpions && lprotons && lthresholds)){
907     AliDebug(1, "Relevant lists not found. Did you forget to run FinishAnalysis()?");
908     return;
909   }
910
911   TGraphErrors *pi, *pr;
912   TGraph *tr;
913   TLegend *leg;
914   Char_t *canvasname=Form("tracklet%d", itracklet);
915   if(icentrality!=-1) canvasname=Form("tracklet%dcentrality%d", itracklet, icentrality);
916   TCanvas *c1 = new TCanvas(canvasname, canvasname, 1024, 768);
917   c1->Divide(3,2);
918   TF1 *threshfit = NULL;
919   for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
920     c1->cd(ieff + 1);
921     gPad->SetGrid(0,0);
922     gPad->SetLeftMargin(0.12);
923     gPad->SetRightMargin(0.08);
924     pi = dynamic_cast<TGraphErrors *>(lpions->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
925     pr = dynamic_cast<TGraphErrors *>(lprotons->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
926     tr = dynamic_cast<TGraph *>(lthresholds->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
927     if(!(pi && pr && tr)) continue;
928
929     // Define nice plot, draw
930     // Axis Title
931     pi->GetXaxis()->SetTitle("p / GeV/c");
932     pi->GetYaxis()->SetTitle("Efficiency");
933     pr->GetXaxis()->SetTitle("p / GeV/c");
934     pr->GetYaxis()->SetTitle("Efficiency");
935     tr->GetXaxis()->SetTitle("p / GeV/c");
936     tr->GetYaxis()->SetTitle("Efficiency");
937     pi->GetYaxis()->SetTitleOffset(1.2);
938     pr->GetYaxis()->SetTitleOffset(1.2);
939     tr->GetYaxis()->SetTitleOffset(1.2);
940     pi->GetXaxis()->SetTitleSize(0.045);
941     pi->GetYaxis()->SetTitleSize(0.045);
942     pr->GetXaxis()->SetTitleSize(0.045);
943     pr->GetYaxis()->SetTitleSize(0.045);
944     tr->GetXaxis()->SetTitleSize(0.045);
945     tr->GetYaxis()->SetTitleSize(0.045);
946     // Axis Range
947     pi->GetYaxis()->SetRangeUser(1e-3, 1.);
948     pr->GetYaxis()->SetRangeUser(1e-3, 1.);
949     tr->GetYaxis()->SetRangeUser(1e-3, 1.);
950     if(pmin >= 0 && pmax >= 0.){
951       pi->GetXaxis()->SetRangeUser(pmin, pmax);
952       pr->GetXaxis()->SetRangeUser(pmin, pmax);
953       tr->GetXaxis()->SetRangeUser(pmin, pmax);
954     }
955     // Marker
956     pi->SetMarkerColor(kRed);
957     pi->SetMarkerStyle(20);
958     pr->SetMarkerColor(kBlue);
959     pr->SetMarkerStyle(21);
960     tr->SetMarkerColor(kBlack);
961     tr->SetMarkerStyle(22);
962     // Title
963     pi->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
964     pr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
965     tr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
966     // Draw
967     pi->Draw("ape"); pr->Draw("pesame"); tr->Draw("pesame");
968
969     // Optionally do Fit
970     if(doFit){
971       threshfit = MakeThresholds(tr, pmin, pmax);
972       threshfit->SetLineColor(kBlack);
973       threshfit->Draw("same");
974     }
975
976     // Add entries to legend
977     if(ieff==0){
978       leg = new TLegend(0.5, 0.65, 0.89, 0.85);
979       leg->SetBorderSize(0);
980       leg->SetFillStyle(0);
981       leg->AddEntry(pi, "Pion Efficiency", "lp");
982       leg->AddEntry(pr, "Proton Efficiency", "lp");
983       leg->AddEntry(tr, "Thresholds", "lp");
984       leg->Draw();
985       gPad->Update();
986     }
987   }
988   c1->cd();
989 }
990
991 //__________________________________________________________________
992 TF1 *AliHFEtrdPIDqa::MakeThresholds(TGraph *threshist, Double_t lowerLimit, Double_t upperLimit){
993   //
994   // Create TF1 containing the threshold parametrisation
995   //
996
997   TF1 *threshparam = new TF1("thresh", "1-[0]-[1]*x-[2]*TMath::Exp(-[3]*x)", 0.1, 10);
998   threshist->Fit(threshparam, "NE", "", lowerLimit, upperLimit);
999   return threshparam;
1000 }
1001
1002 //__________________________________________________________________
1003 void AliHFEtrdPIDqa::ClearLists(){
1004   //
1005   // Clear lists for particle efficiencies and thresholds
1006   //
1007   if(fPionEfficiencies){
1008     fPionEfficiencies->Delete();
1009     delete fPionEfficiencies;
1010     fPionEfficiencies = NULL;
1011   }
1012   if(fProtonEfficiencies){
1013     fProtonEfficiencies->Delete();
1014     delete fProtonEfficiencies;
1015     fProtonEfficiencies = NULL;
1016   }
1017   if(fThresholds){
1018     fThresholds->Delete();
1019     delete fThresholds;
1020     fThresholds = NULL;
1021   }
1022 }
1023
1024 //__________________________________________________________________
1025 Double_t AliHFEtrdPIDqa::EvalPionEfficiency(Int_t ntls, Int_t eEff, Double_t p, Int_t icentrality){
1026   //
1027   // calculate pion efficiency
1028   // Arguments:
1029   //   Number of tracklets
1030   //   Electron Efficiency
1031   //   Momentum
1032   //
1033   TGraphErrors *measurement = GetPionEfficiency(ntls, eEff, icentrality);
1034   if(!measurement) return -1.;
1035   return measurement->Eval(p);
1036 }
1037
1038 //__________________________________________________________________
1039 Double_t AliHFEtrdPIDqa::EvalProtonEfficiency(Int_t ntls, Int_t eEff, Double_t p, Int_t icentrality){
1040   //
1041   // calculate proton efficiency
1042   // Arguments:
1043   //   Number of tracklets
1044   //   Electron Efficiency
1045   //   Momentum
1046   //
1047   TGraphErrors *measurement = GetProtonEfficiency(ntls, eEff, icentrality);
1048   if(!measurement) return -1.;
1049   return measurement->Eval(p);
1050 }
1051
1052 //__________________________________________________________________
1053 Double_t AliHFEtrdPIDqa::EvalThreshold(Int_t ntls, Int_t eEff, Double_t p, Int_t icentrality){
1054   //
1055   // Get the threshold to reach a certain electron efficiency
1056   // Arguments:
1057   //   Number of tracklets
1058   //   Electron Efficiency
1059   //   Momentum
1060   //
1061   TGraph *measurement = GetThreshold(ntls, eEff, icentrality);
1062   if(!measurement) return -1.;
1063   return measurement->Eval(p);
1064 }
1065
1066 //__________________________________________________________________
1067 TGraphErrors *AliHFEtrdPIDqa::GetPionEfficiency(Int_t ntracklets, Int_t eleffpercent, Int_t icentrality){
1068   //
1069   // Get Graph with pion efficiencies
1070   //
1071   Char_t *listname=Form("%dTracklets", ntracklets);
1072   if(icentrality!=-1) listname=Form("%dTracklets%dCentrality", ntracklets, icentrality);
1073   TList *graphs = dynamic_cast<TList *>(fPionEfficiencies->FindObject(listname));
1074   if(!graphs) return NULL;
1075   return dynamic_cast<TGraphErrors *>(graphs->FindObject(Form("eff%d", eleffpercent)));
1076 }
1077
1078 //__________________________________________________________________
1079 TGraphErrors *AliHFEtrdPIDqa::GetProtonEfficiency(Int_t ntracklets, Int_t eleffpercent, Int_t icentrality){
1080   // 
1081   // Get Graph with proton efficiencies
1082   //
1083   Char_t *listname=Form("%dTracklets", ntracklets);
1084   if(icentrality!=-1) listname=Form("%dTracklets%dCentrality", ntracklets, icentrality);
1085   TList *graphs = dynamic_cast<TList *>(fProtonEfficiencies->FindObject(listname));
1086   if(!graphs) return NULL;
1087   return dynamic_cast<TGraphErrors *>(graphs->FindObject(Form("eff%d", eleffpercent)));
1088 }
1089
1090 //__________________________________________________________________
1091 TGraph *AliHFEtrdPIDqa::GetThreshold(Int_t ntracklets, Int_t eleffpercent, Int_t icentrality){
1092   //
1093   // Get Graph with threshols
1094   //
1095   Char_t *listname=Form("%dTracklets", ntracklets);
1096   if(icentrality!=-1) listname=Form("%dTracklets%dCentrality", ntracklets, icentrality);
1097   TList *graphs = dynamic_cast<TList *>(fThresholds->FindObject(listname));
1098   if(!graphs) return NULL;
1099   return dynamic_cast<TGraph *>(graphs->FindObject(Form("eff%d", eleffpercent)));
1100 }
1101