Major update of the HFE package (comments inside the code
[u/mrichter/AliRoot.git] / PWG3 / 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 <TClass.h>
26 #include <TCanvas.h>
27 #include <TF1.h>
28 #include <TFile.h>
29 #include <TGraph.h>
30 #include <TGraphErrors.h>
31 #include <THnSparse.h>
32 #include <TH1.h>
33 #include <TH2.h>
34 #include <TIterator.h>
35 #include <TLegend.h>
36 #include <TList.h>
37 #include <TMath.h>
38 #include <TObjArray.h>
39 #include <TString.h>
40
41 #include "AliAODTrack.h"
42 #include "AliAODPid.h"
43 #include "AliESDtrack.h"
44 #include "AliHFEtrdPIDqa.h"
45 #include "AliHFEtools.h"
46 #include "AliHFEpidTRD.h"
47 #include "AliLog.h"
48
49 ClassImp(AliHFEtrdPIDqa)
50
51 const Double_t AliHFEtrdPIDqa::fgkElectronEff[kNElectronEffs] = {0.7,0.75, 0.8, 0.85, 0.9, 0.95};
52 //_______________________________________________________________
53 // Definition of the common binning 
54 const Int_t AliHFEtrdPIDqa::fgkNBinsCommon[kQuantitiesCommon] = {
55   AliPID::kSPECIES + 1,         // species
56   40,                           // p-bins
57   AliESDtrack::kTRDnPlanes + 1 // tracklets including 0
58 };
59 const Double_t AliHFEtrdPIDqa::fgkMinBinCommon[kQuantitiesCommon] = {
60   -1,      // species
61   0.1,     // p-bins
62   0        // tracklets including 0
63 };
64
65 const Double_t AliHFEtrdPIDqa::fgkMaxBinCommon[kQuantitiesCommon] = {
66   AliPID::kSPECIES,               // species
67   10.,                            // p-bins
68   AliESDtrack::kTRDnPlanes + 1    // tracklets including 0
69 };
70 //_______________________________________________________________
71
72 //__________________________________________________________________
73 AliHFEtrdPIDqa::AliHFEtrdPIDqa():
74   TNamed("trdPIDqa", ""),
75   fTRDpid(NULL),
76   fLikeTRD(NULL),
77   fQAtrack(NULL),
78   fQAdEdx(NULL),
79   fTRDtruncMean(NULL),
80   fPionEfficiencies(NULL),
81   fProtonEfficiencies(NULL),
82   fKaonEfficiencies(NULL),
83   fThresholds(NULL)
84 {
85   //
86   // Default Constructor
87   //
88 }
89
90 //__________________________________________________________________
91 AliHFEtrdPIDqa::AliHFEtrdPIDqa(const Char_t *name):
92   TNamed(name, ""),
93   fTRDpid(NULL),
94   fLikeTRD(NULL),
95   fQAtrack(NULL),
96   fQAdEdx(NULL),
97   fTRDtruncMean(NULL),
98   fPionEfficiencies(NULL),
99   fProtonEfficiencies(NULL),
100   fKaonEfficiencies(NULL),
101   fThresholds(NULL)
102 {
103   //
104   // Main Constructor
105   //
106 }
107
108 //__________________________________________________________________
109 AliHFEtrdPIDqa::AliHFEtrdPIDqa(const AliHFEtrdPIDqa &ref):
110   TNamed(ref),
111   fTRDpid(NULL),
112   fLikeTRD(NULL),
113   fQAtrack(NULL),
114   fQAdEdx(NULL),
115   fTRDtruncMean(NULL),
116   fPionEfficiencies(NULL),
117   fProtonEfficiencies(NULL),
118   fKaonEfficiencies(NULL),
119   fThresholds(NULL)
120 {
121   //
122   // Copy constructor
123   //
124   ref.Copy(*this);
125 }
126
127 //__________________________________________________________________
128 AliHFEtrdPIDqa &AliHFEtrdPIDqa::operator=(const AliHFEtrdPIDqa &ref){
129   //
130   // Assignment operator
131   //
132   if(this != &ref)
133     ref.Copy(*this);
134   return *this;
135 }
136
137 //__________________________________________________________________
138 AliHFEtrdPIDqa::~AliHFEtrdPIDqa(){
139   //
140   // Destructor
141   //
142   if(fTRDpid) delete fTRDpid;
143   if(fLikeTRD) delete fLikeTRD;
144   if(fQAtrack) delete fQAtrack;
145   if(fQAdEdx) delete fQAdEdx;
146   if(fTRDtruncMean) delete fTRDtruncMean;
147   
148   if(fPionEfficiencies) delete fPionEfficiencies;
149   if(fProtonEfficiencies) delete fProtonEfficiencies;
150   if(fKaonEfficiencies) delete fKaonEfficiencies;
151 }
152
153 //__________________________________________________________________
154 void AliHFEtrdPIDqa::Copy(TObject &ref) const{
155   //
156   // Copies content of this object into object ref
157   //
158   TNamed::Copy(ref);
159
160   AliHFEtrdPIDqa &target = dynamic_cast<AliHFEtrdPIDqa &>(ref);
161   target.fTRDpid = fTRDpid;
162   target.fLikeTRD = dynamic_cast<THnSparseF *>(fLikeTRD->Clone());
163   target.fQAtrack = dynamic_cast<THnSparseF *>(fQAtrack->Clone());
164   target.fQAdEdx = dynamic_cast<THnSparseF *>(fQAdEdx->Clone());
165   target.fTRDtruncMean = dynamic_cast<THnSparseF *>(fTRDtruncMean->Clone());
166 }
167
168 //__________________________________________________________________
169 Long64_t AliHFEtrdPIDqa::Merge(TCollection *coll){
170   //
171   // Merge objects
172   //
173   if(!coll) return 0;
174   if(coll->IsEmpty()) return 1;
175   
176   TIterator *it = coll->MakeIterator();
177   TObject *o = NULL;
178   Long64_t count = 0;
179   while((o = it->Next())){
180     AliHFEtrdPIDqa *refQA = dynamic_cast<AliHFEtrdPIDqa *>(o);
181     if(!refQA) continue;
182
183     if(fLikeTRD) fLikeTRD->Add(refQA->fLikeTRD);
184     if(fQAtrack) fQAtrack->Add(refQA->fQAtrack);
185     if(fQAdEdx) fQAdEdx->Add(refQA->fQAdEdx);
186     if(fTRDtruncMean) fTRDtruncMean->Add(refQA->fTRDtruncMean);
187     count++;
188   }
189   delete it;
190   return count+1;
191 }
192
193 //__________________________________________________________________
194 void AliHFEtrdPIDqa::Init(){
195   //
196   // Initialize Object
197   //
198
199   CreateLikelihoodHistogram();
200   CreateQAHistogram();
201   CreatedEdxHistogram();
202   CreateHistoTruncatedMean();
203
204   fTRDpid = new AliHFEpidTRD("QAtrdPID");
205 }
206
207 //__________________________________________________________________
208 void AliHFEtrdPIDqa::CreateLikelihoodHistogram(){
209   //
210   // Create Histogram for TRD Likelihood Studies
211   //
212   Int_t nbins[kQuantitiesLike]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
213   nbins[kElectronLike] = 100;
214   Double_t binMin[kQuantitiesLike]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
215   Double_t binMax[kQuantitiesLike]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
216   binMin[kElectronLike] = 0.;      
217   binMax[kElectronLike] = 1.;
218
219   fLikeTRD = new THnSparseF("fLikeTRD","TRD Likelihood Studies", kQuantitiesLike, nbins, binMin, binMax);
220   Double_t *binLog = AliHFEtools::MakeLogarithmicBinning(nbins[kP], binMin[kP], binMax[kP]);
221   fLikeTRD->GetAxis(kP)->Set(nbins[kP], binLog);
222   delete[] binLog;
223 }
224
225 //__________________________________________________________________
226 void AliHFEtrdPIDqa::CreateQAHistogram(){
227   //
228   // Create Histogram for Basic TRD PID QA
229   //
230   AliDebug(1, "Called");
231   Int_t nbins[kQuantitiesQA]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
232   nbins[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1;
233   nbins[kNClusters] = 200;
234   Double_t binMin[kQuantitiesQA]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
235   binMin[kNonZeroTrackletCharge] = 0.;      
236   binMin[kNClusters] = 0.;      
237   Double_t binMax[kQuantitiesQA]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
238   binMax[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1.;
239   binMax[kNClusters] = 200.;
240
241   fQAtrack = new THnSparseF("fQAtrack","TRD QA Histogram", kQuantitiesQA, nbins, binMin, binMax);
242   Double_t *binLog = AliHFEtools::MakeLogarithmicBinning(nbins[kP], binMin[kP], binMax[kP]);
243   fQAtrack->GetAxis(kP)->Set(nbins[kP], binLog);
244   delete[] binLog;
245 }
246
247 //__________________________________________________________________
248 void AliHFEtrdPIDqa::CreatedEdxHistogram(){
249   //
250   // Create QA histogram for dEdx investigations
251   //
252   AliDebug(1, "Called");
253   Int_t nbins[kQuantitiesdEdx]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
254   nbins[kdEdx] = 100;
255   nbins[kNclusters] = 261;
256   nbins[kNonZeroSlices] = 9; 
257   Double_t binMin[kQuantitiesdEdx]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
258   binMin[kdEdx] = 0.;     
259   binMin[kNclusters] = 0;
260   Double_t binMax[kQuantitiesdEdx]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
261   binMax[kdEdx] = 100000.;
262   binMax[kNclusters] = 260.;
263   binMax[kNonZeroSlices] = 8.;
264
265   fQAdEdx = new THnSparseF("fQAdEdx","TRD summed dEdx", kQuantitiesdEdx, nbins, binMin, binMax);
266   Double_t *binLog = AliHFEtools::MakeLogarithmicBinning(nbins[kP], binMin[kP], binMax[kP]);
267   fQAdEdx->GetAxis(kP)->Set(nbins[kP], binLog);
268   delete[] binLog;
269 }
270
271 //__________________________________________________________________
272 void AliHFEtrdPIDqa::CreateHistoTruncatedMean(){
273   //
274   // Create Histogram for Basic TRD PID QA
275   //
276   AliDebug(1, "Called");
277   Int_t nbins[kQuantitiesTruncMean]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
278   nbins[kTPCdEdx] = 600;
279   nbins[kTRDdEdxMethod1] = 1000;
280   nbins[kTRDdEdxMethod2] = 1000;
281   Double_t binMin[kQuantitiesTruncMean]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
282   binMin[kTPCdEdx] = 0.;      
283   binMin[kTRDdEdxMethod1] = 0.;      
284   binMin[kTRDdEdxMethod2] = 0.;      
285   Double_t binMax[kQuantitiesTruncMean]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
286   binMax[kTPCdEdx] = 600;
287   binMax[kTRDdEdxMethod1] = 20000.;
288   binMax[kTRDdEdxMethod2] = 20000.;
289
290   fTRDtruncMean = new THnSparseF("fTRDtruncMean","TRD TruncatedMean studies", kQuantitiesTruncMean, nbins, binMin, binMax);
291   Double_t *binLog = AliHFEtools::MakeLogarithmicBinning(nbins[kP], binMin[kP], binMax[kP]);
292   fTRDtruncMean->GetAxis(kP)->Set(nbins[kP], binLog);
293   delete[] binLog;
294 }
295
296
297 //__________________________________________________________________
298 void AliHFEtrdPIDqa::ProcessTracks(TObjArray * const tracks, Int_t species){
299   //
300   // Process Electron Tracks
301   //
302   if(species < -1 || species >= AliPID::kSPECIES) return;
303   TIterator *it = tracks->MakeIterator();
304   AliVTrack *track = NULL;
305   while((track = dynamic_cast<AliVTrack *>(it->Next()))){
306     if(track) ProcessTrack(track, species);
307   }
308   delete it;
309 }
310
311 //__________________________________________________________________
312 void AliHFEtrdPIDqa::ProcessTrack(AliVTrack *track, Int_t species){
313   //
314   // Process Single Track
315   //
316   if(TString(track->IsA()->GetName()).CompareTo("AliESDtrack") == 0)
317     ProcessTrackESD(dynamic_cast<AliESDtrack *>(track), species);
318   else if(TString(track->IsA()->GetName()).CompareTo("AliAODTrack") == 0)
319     ProcessTrackAOD(dynamic_cast<AliAODTrack *>(track), species);
320 }
321
322
323 //__________________________________________________________________
324 void AliHFEtrdPIDqa::ProcessTrackESD(AliESDtrack *track, Int_t species){
325   //
326   // Process single ESD track
327   //
328   if((track->GetStatus() & AliESDtrack::kTRDout) == 0) return;  // require TRD track
329   FillTRDLikelihoods(track, species);
330   FillTRDQAplots(track, species);
331 }
332
333 //__________________________________________________________________
334 void AliHFEtrdPIDqa::ProcessTrackAOD(AliAODTrack * const track, Int_t /*species*/){
335   //
336   // Process single AOD track
337   // AOD PID object required
338   //
339   AliAODPid *trackPID = track->GetDetPid();
340   if(!trackPID) return;
341
342 }
343
344 //__________________________________________________________________
345 void AliHFEtrdPIDqa::FillTRDLikelihoods(AliESDtrack *track, Int_t species){
346   //
347   // Fill TRD Likelihood Histogram
348   //
349   Double_t trdLike[AliPID::kSPECIES];
350   track->GetTRDpid(trdLike);
351   const AliExternalTrackParam *outerPars = track->GetOuterParam();
352
353   Double_t quantities[kQuantitiesLike]; memset(quantities, 0, sizeof(Double_t) * kQuantitiesLike);
354   // we store:
355   // species
356   // p
357   // ntracklets
358   // Electron Likelihood
359   quantities[kSpecies] = species;
360   quantities[kP] = outerPars ? outerPars->P() : track->P();
361   quantities[kNTracklets] = track->GetTRDntrackletsPID();
362   quantities[kElectronLike] = trdLike[AliPID::kElectron];
363
364   fLikeTRD->Fill(quantities);
365 }
366
367 //__________________________________________________________________
368 void AliHFEtrdPIDqa::FillTRDQAplots(AliESDtrack *track, Int_t species){
369   //
370   // Fill QA Plots containing further information
371   //
372   const AliExternalTrackParam *outerPars = track->GetOuterParam();
373
374   Double_t quantitiesQA[kQuantitiesQA], quantitiesdEdx[kQuantitiesdEdx], quantitiesTruncMean[kQuantitiesTruncMean];
375   // we store:
376   // 1. QA
377   // species
378   // p
379   // ntracklets
380   // Non-zero tracklet charges
381   // Number of clusters / full track
382   // 2. dEdx
383   // species
384   // p
385   // ntracklets
386   // dEdx
387   // 3. Truncated Mean
388   // ...
389   // TPC dEdx
390   // TRD dEdx Method 1
391   // TRD dEdx Method 2
392   quantitiesQA[kSpecies]  = quantitiesdEdx[kSpecies] 
393                           = quantitiesTruncMean[kSpecies] 
394                           = species;
395   quantitiesQA[kP]  = quantitiesdEdx[kP] 
396                     = quantitiesTruncMean[kP]
397                     = outerPars ? outerPars->P() : track->P();
398   quantitiesQA[kNTracklets] = quantitiesdEdx[kNTracklets] 
399                             = quantitiesTruncMean[kNTracklets]
400                             = track->GetTRDntrackletsPID();
401   quantitiesQA[kNClusters] = quantitiesdEdx[kNclusters] = track->GetTRDncls();
402   
403
404   Double_t dEdxSum = 0., qSlice = 0.;
405   // remove the last slice from the histogram
406   Int_t ntrackletsNonZero = 0, nSlices = track->GetNumberOfTRDslices() - 1, nSlicesNonZero = 0;
407   for(Int_t iplane = 0; iplane < AliESDtrack::kTRDnPlanes; iplane++){
408     dEdxSum = 0.;
409     for(Int_t islice = 0; islice < nSlices; islice++){
410       qSlice = track->GetTRDslice(iplane, islice);
411       if(qSlice > 1e-1){
412         // cut out 0 slices
413         nSlicesNonZero++;
414         dEdxSum += qSlice;
415       }
416     }
417     quantitiesdEdx[kNonZeroSlices] = nSlicesNonZero;
418     quantitiesdEdx[kdEdx] = dEdxSum;
419     if(dEdxSum) ntrackletsNonZero++;
420     // Fill dEdx histogram
421     if(dEdxSum > 1e-1) fQAdEdx->Fill(quantitiesdEdx); // Cut out 0 entries
422   }
423   quantitiesQA[kNonZeroTrackletCharge] = ntrackletsNonZero;
424   fQAtrack->Fill(quantitiesQA);
425
426   quantitiesTruncMean[kTPCdEdx] = track->GetTPCsignal();
427   quantitiesTruncMean[kTRDdEdxMethod1] = fTRDpid->GetTRDSignalV1(track);
428   quantitiesTruncMean[kTRDdEdxMethod2] = fTRDpid->GetTRDSignalV2(track);
429   fTRDtruncMean->Fill(quantitiesTruncMean);
430 }
431
432 /////////////////////////////////////////////////////////
433 //
434 // Code for Post Processing
435 //
436 // //////////////////////////////////////////////////////
437
438 //__________________________________________________________________
439 void AliHFEtrdPIDqa::FinishAnalysis(){
440   //
441   // Finish Analysis:
442   // Calculate Electron Efficiency for ntracklets = 4...6
443   // Calculate thresholds for ntracklets = 4...6
444   //
445   
446   if(!fPionEfficiencies){
447     fPionEfficiencies = new TList;
448     fPionEfficiencies->SetName("pionEfficiencies");
449   }
450   if(!fProtonEfficiencies){
451     fProtonEfficiencies = new TList;
452     fProtonEfficiencies->SetName("protonEfficiencies");
453   }
454   if(!fThresholds){
455     fThresholds = new TList;
456     fThresholds->SetName("thresholds");
457   }
458
459   for(Int_t itr = 4; itr <= 6; itr++){
460     printf("========================================\n");
461     printf("Analysing %d trackltes\n", itr);
462     printf("========================================\n");
463     AnalyseNTracklets(itr);
464   }
465 }
466
467 //__________________________________________________________________
468 void AliHFEtrdPIDqa::StoreResults(const Char_t *filename){
469   //
470   // Store histos into a root file
471   //
472   TFile *outfile = new TFile(filename, "RECREATE");
473   outfile->cd();
474   fPionEfficiencies->Clone()->Write(fPionEfficiencies->GetName(), kSingleKey);
475   fProtonEfficiencies->Clone()->Write(fProtonEfficiencies->GetName(), kSingleKey);
476   fThresholds->Clone()->Write(fThresholds->GetName(), kSingleKey);
477   outfile->Close();
478   delete outfile;
479 }
480
481 //__________________________________________________________________
482 void AliHFEtrdPIDqa::SaveThresholdParameters(const Char_t *name){
483   //
484   // Fit the threshold histograms with the given parametrisation
485   // and store the TF1 in the file
486   //
487
488   if(!fThresholds){
489     AliError("Threshold histograms have to be created first");
490     return;
491   }
492
493   printf("========================================\n");
494   printf("Calculating threshold parameters\n");
495   printf("========================================\n");
496
497   TList *outlist = new TList;
498   outlist->SetName("thresholdTRD");
499   
500   TGraph *threshhist = NULL;
501   TF1 *threshparam = NULL;
502   TList *lHistos = NULL, *lFormulas = NULL;
503   for(Int_t itracklet = 4; itracklet <= 6; itracklet++){
504   
505     printf("-------------------------------\n");
506     printf("Processing %d tracklets\n", itracklet);
507     printf("-------------------------------\n");
508
509     lHistos = dynamic_cast<TList *>(fThresholds->FindObject(Form("%dTracklets", itracklet)));
510     if(!lHistos){
511       AliError(Form("Threshold histograms for the case %d tracklets not found", itracklet));
512       continue;
513     }
514     lFormulas = new TList;
515     lFormulas->SetName(Form("%dTracklets", itracklet));
516     outlist->Add(lFormulas);
517     
518     for(Int_t ieff = 0; ieff <  kNElectronEffs; ieff++){
519       
520       printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n");
521       printf("Processing Electron Efficiency %f\n", fgkElectronEff[ieff]);
522       printf(">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n");
523
524       threshhist = dynamic_cast<TGraph *>(lHistos->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
525       threshparam = MakeThresholds(threshhist);
526       threshparam->SetName(Form("thresh_%d_%d", itracklet, static_cast<Int_t>(fgkElectronEff[ieff] * 100)));
527     }
528   }
529
530   // store the output
531   TFile *outfile = new TFile(name, "RECREATE");
532   outfile->cd();
533   outlist->Write(outlist->GetName(), kSingleKey);
534   outfile->Close();
535   delete outfile;
536 }
537
538 //__________________________________________________________________
539 void AliHFEtrdPIDqa::AnalyseNTracklets(Int_t nTracklets){
540   //
541   // Calculate Pion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
542   // elPion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
543   // electron efficiencies
544   //
545   Int_t binTracklets = fLikeTRD->GetAxis(kNTracklets)->FindBin(nTracklets);
546   fLikeTRD->GetAxis(kNTracklets)->SetRange(binTracklets, binTracklets);
547   
548   Int_t binElectrons = fLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kElectron); 
549   AliDebug(1, Form("BinElectrons %d", binElectrons));
550   Int_t binPions = fLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kPion);
551   AliDebug(1, Form("BinPions %d", binPions));
552   Int_t binProtons =  fLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kProton);
553   AliDebug(1, Form("BinProtons %d", binProtons));
554   fLikeTRD->GetAxis(kSpecies)->SetRange(binElectrons, binElectrons);
555   TH2 *likeElectron = fLikeTRD->Projection(kElectronLike, kP);
556   likeElectron->SetName("likeElectron");
557   fLikeTRD->GetAxis(kSpecies)->SetRange(binPions, binPions);
558   TH2 *likePion = fLikeTRD->Projection(kElectronLike, kP);
559   likePion->SetName("likePion");
560   fLikeTRD->GetAxis(kSpecies)->SetRange(binProtons, binProtons);
561   TH2 *likeProton = fLikeTRD->Projection(kElectronLike, kP);
562   likeProton->SetName("likeProton");
563   
564   // Undo ranges
565   fLikeTRD->GetAxis(kSpecies)->SetRange(0, fLikeTRD->GetAxis(kSpecies)->GetNbins());
566   fLikeTRD->GetAxis(kNTracklets)->SetRange(0, fLikeTRD->GetAxis(kNTracklets)->GetNbins());
567
568   // Prepare List for output
569   TList *listPions = new TList; listPions->SetName(Form("%dTracklets", nTracklets));
570   TList *listProtons = new TList; listProtons->SetName(Form("%dTracklets", nTracklets));
571   TList *listThresholds = new TList; listThresholds->SetName(Form("%dTracklets", nTracklets));
572   fPionEfficiencies->Add(listPions);
573   fProtonEfficiencies->Add(listProtons);
574   fThresholds->Add(listThresholds);
575
576   TH1 *probsEl = NULL, *probsPi = NULL, *probsPr = NULL;
577   TGraphErrors *effPi = NULL, *effPr = NULL; TGraph *thresholds = NULL;
578   Double_t p = 0, dp = 0;
579   Int_t threshbin = 0;
580   Double_t noElEff[2]; // value and error
581   for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
582     printf("-----------------------------------------\n");
583     printf("Doing Electron Efficiency %f\n", fgkElectronEff[ieff]);
584     printf("-----------------------------------------\n");
585     effPi = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
586     effPi->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
587     effPr = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
588     effPr->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
589     thresholds = new TGraph(likeElectron->GetXaxis()->GetNbins());
590     thresholds->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
591   
592     // Add to lists
593     listPions->Add(effPi);
594     listProtons->Add(effPr);
595     listThresholds->Add(thresholds);
596
597     for(Int_t imom = 1; imom <= likeElectron->GetXaxis()->GetLast(); imom++){
598       p = likeElectron->GetXaxis()->GetBinCenter(imom);
599       dp = likeElectron->GetXaxis()->GetBinWidth(imom)/2;
600
601       probsEl = likeElectron->ProjectionY("el", imom);
602       if(!probsEl->GetEntries()) continue;
603       probsEl->Scale(1./probsEl->Integral());
604       probsPi = likePion->ProjectionY("pi", imom);
605       if(!probsPi->GetEntries()) continue;
606       probsPi->Scale(1./probsPi->Integral());
607       probsPr = likeProton->ProjectionY("pr", imom);
608       if(!probsPr->GetEntries()) continue;
609       probsPr->Scale(1./probsPr->Integral());
610       AliDebug(1, Form("Calculating Values for p = %f", p));
611
612       // Calculare threshold we need to achive the x% electron Efficiency
613       threshbin = GetThresholdBin(probsEl, fgkElectronEff[ieff]);
614       thresholds->SetPoint(imom - 1, p, probsEl->GetXaxis()->GetBinCenter(threshbin));
615       AliDebug(1, Form("threshold %d|%f", threshbin, probsEl->GetXaxis()->GetBinCenter(threshbin)));
616
617       // Calculate non-electronEfficiency and error
618       CalculateEfficiency(probsPi, threshbin, noElEff);
619       AliDebug(1, Form("Pion Efficiency %f", noElEff[0]));
620       effPi->SetPoint(imom - 1, p, noElEff[0]);
621       effPi->SetPointError(imom - 1, dp, noElEff[1]);
622       CalculateEfficiency(probsPr, threshbin, noElEff);
623       effPr->SetPoint(imom - 1, p, noElEff[0]);
624       effPr->SetPointError(imom - 1, dp, noElEff[1]);
625       AliDebug(1, Form("Proton Efficiency %f", noElEff[0]));
626  
627       // cleanup
628       delete probsEl;
629       delete probsPi;
630       delete probsPr;
631     }
632   }
633
634   // remove temporary histograms
635   delete likeElectron;
636   delete likePion;
637   delete likeProton;
638 }
639
640 //__________________________________________________________________
641 Int_t AliHFEtrdPIDqa::GetThresholdBin(TH1 * const input, Double_t eff){
642   //
643   // Calculate the threshold bin  
644   //
645   Double_t integralEff = 0.;
646   Int_t currentBin = 0;
647   for(Int_t ibin = input->GetXaxis()->GetLast(); ibin >= input->GetXaxis()->GetFirst(); ibin--){
648     currentBin = ibin;
649     integralEff += input->GetBinContent(ibin);
650     if(integralEff >= eff){
651       // we found the matching bin, break the loop
652       break;
653     }
654   }
655   return currentBin;
656 }
657
658 //__________________________________________________________________
659 Bool_t AliHFEtrdPIDqa::CalculateEfficiency(TH1 * const input, Int_t threshbin, Double_t *par){
660   // 
661   // Calculate non-electron efficiency
662   //
663   Double_t integralEff = 0; 
664   for(Int_t ibin = threshbin; ibin <= input->GetXaxis()->GetLast(); ibin++) 
665     integralEff += input->GetBinContent(ibin);
666   par[0] = integralEff;
667
668   // @TODO: Error calculation
669   par[1] = 0;
670
671   return kTRUE;
672 }
673
674 //__________________________________________________________________
675 void AliHFEtrdPIDqa::DrawTracklet(Int_t itracklet){
676   //
677   // Draw efficiencies and threshold as function of p
678   //
679   if(!(fPionEfficiencies && fProtonEfficiencies && fThresholds)){
680     AliError("No graphs to draw available");  
681     return;
682   }
683
684   TList *lpions = dynamic_cast<TList *>(fPionEfficiencies->FindObject(Form("%dTracklets", itracklet))); 
685   TList *lprotons = dynamic_cast<TList *>(fProtonEfficiencies->FindObject(Form("%dTracklets", itracklet))); 
686   
687   TList *lthresholds = dynamic_cast<TList *>(fThresholds->FindObject(Form("%dTracklets", itracklet))); 
688
689   TGraphErrors *pi, *pr;
690   TGraph *tr;
691   TLegend *leg;
692   TCanvas *c1 = new TCanvas(Form("tracklet%d", itracklet), Form("Tracklet %d", itracklet), 1024, 768);
693   c1->Divide(3,2);
694   for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
695     c1->cd(ieff + 1);
696     leg = new TLegend(0.6, 0.7, 0.89, 0.89);
697     leg->SetBorderSize(0);
698     leg->SetFillStyle(0);
699     pi = dynamic_cast<TGraphErrors *>(lpions->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
700     pr = dynamic_cast<TGraphErrors *>(lprotons->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
701     tr = dynamic_cast<TGraph *>(lthresholds->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
702
703     // Define nice plot, draw
704     // Axis Title
705     pi->GetXaxis()->SetTitle("p / GeV/c");
706     pi->GetYaxis()->SetTitle("Efficiency");
707     pr->GetXaxis()->SetTitle("p / GeV/c");
708     pr->GetYaxis()->SetTitle("Efficiency");
709     tr->GetXaxis()->SetTitle("p / GeV/c");
710     tr->GetYaxis()->SetTitle("Efficiency");
711     // Axis Range
712     pi->GetYaxis()->SetRangeUser(0., 1.);
713     pr->GetYaxis()->SetRangeUser(0., 1.);
714     tr->GetYaxis()->SetRangeUser(0., 1.);
715     // Marker
716     pi->SetMarkerColor(kRed);
717     pi->SetMarkerStyle(20);
718     pr->SetMarkerColor(kBlue);
719     pr->SetMarkerStyle(21);
720     tr->SetMarkerColor(kBlack);
721     tr->SetMarkerStyle(22);
722     // Title
723     pi->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
724     pr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
725     tr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
726     // Draw
727     pi->Draw("ape"); pr->Draw("pesame"); tr->Draw("psame");
728
729     // Add entries to legend
730     leg->AddEntry(pi, "Pion Efficiency", "lp");
731     leg->AddEntry(pr, "Proton Efficiency", "lp");
732     leg->AddEntry(tr, "Thresholds", "lp");
733     leg->Draw();
734     c1->Update();
735   }
736 }
737
738 //__________________________________________________________________
739 TF1 *AliHFEtrdPIDqa::MakeThresholds(TGraph *threshist){
740   //
741   // Create TF1 containing the threshold parametrisation
742   //
743
744   TF1 *threshparam = new TF1("thresh", "1-[0]-[1]*x-[2]*TMath::Exp(-[3]*x)", 0.1, 10);
745   threshist->Fit(threshparam, "NE", "", 0, 10);
746   return threshparam;
747 }
748
749 //__________________________________________________________________
750 void AliHFEtrdPIDqa::ClearLists(){
751   //
752   // Clear lists for particle efficiencies and thresholds
753   //
754   if(fPionEfficiencies){
755     fPionEfficiencies->Delete();
756     delete fPionEfficiencies;
757     fPionEfficiencies = NULL;
758   }
759   if(fProtonEfficiencies){
760     fProtonEfficiencies->Delete();
761     delete fProtonEfficiencies;
762     fProtonEfficiencies = NULL;
763   }
764   if(fThresholds){
765     fThresholds->Delete();
766     delete fThresholds;
767     fThresholds = NULL;
768   }
769 }