3b24486341d5ca94435fe2ff1d521264941c9886
[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 <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 {
84   //
85   // Default Constructor
86   //
87 }
88
89 //__________________________________________________________________
90 AliHFEtrdPIDqa::AliHFEtrdPIDqa(const Char_t *name):
91   TNamed(name, ""),
92   fTRDpid(NULL),
93   fHistos(NULL),
94   fPionEfficiencies(NULL),
95   fProtonEfficiencies(NULL),
96   fKaonEfficiencies(NULL),
97   fThresholds(NULL),
98   fShowMessage(kFALSE)
99 {
100   //
101   // Main Constructor
102   //
103 }
104
105 //__________________________________________________________________
106 AliHFEtrdPIDqa::AliHFEtrdPIDqa(const AliHFEtrdPIDqa &ref):
107   TNamed(ref),
108   fTRDpid(NULL),
109   fHistos(NULL),
110   fPionEfficiencies(NULL),
111   fProtonEfficiencies(NULL),
112   fKaonEfficiencies(NULL),
113   fThresholds(NULL),
114   fShowMessage(kFALSE)
115 {
116   //
117   // Copy constructor
118   //
119   ref.Copy(*this);
120 }
121
122 //__________________________________________________________________
123 AliHFEtrdPIDqa &AliHFEtrdPIDqa::operator=(const AliHFEtrdPIDqa &ref){
124   //
125   // Assignment operator
126   //
127   if(this != &ref)
128     ref.Copy(*this);
129   return *this;
130 }
131
132 //__________________________________________________________________
133 AliHFEtrdPIDqa::~AliHFEtrdPIDqa(){
134   //
135   // Destructor
136   //
137   if(fTRDpid) delete fTRDpid;
138   if(fHistos) delete fHistos;
139   if(fPionEfficiencies) delete fPionEfficiencies;
140   if(fProtonEfficiencies) delete fProtonEfficiencies;
141   if(fKaonEfficiencies) delete fKaonEfficiencies;
142 }
143
144 //__________________________________________________________________
145 void AliHFEtrdPIDqa::Copy(TObject &ref) const{
146   //
147   // Copies content of this object into object ref
148   //
149   TNamed::Copy(ref);
150
151   AliHFEtrdPIDqa &target = dynamic_cast<AliHFEtrdPIDqa &>(ref);
152   target.fTRDpid = fTRDpid;
153   target.fHistos = dynamic_cast<AliHFEcollection *>(fHistos->Clone());
154 }
155
156 //__________________________________________________________________
157 Long64_t AliHFEtrdPIDqa::Merge(TCollection *coll){
158   //
159   // Merge objects
160   //
161   if(!coll) return 0;
162   if(coll->IsEmpty()) return 1;
163   
164   AliHFEtrdPIDqa *refQA = NULL;
165   TIter it(coll);
166   TObject *o = NULL;
167   Long64_t count = 0;
168   TList listHistos;
169   while((o = it())){
170     refQA = dynamic_cast<AliHFEtrdPIDqa *>(o);
171     if(!refQA) continue;
172
173     listHistos.Add(refQA->fHistos);
174     count++; 
175   }
176   fHistos->Merge(&listHistos);
177   return count+1;
178 }
179
180 //__________________________________________________________________
181 void AliHFEtrdPIDqa::Browse(TBrowser *b){
182   //
183   // Enable Browser functionality
184   //
185   if(b){
186     // Add objects to the browser
187     if(fHistos) b->Add(fHistos, fHistos->GetName());
188     if(fPionEfficiencies) b->Add(fPionEfficiencies, "Pion Efficiencies");
189     if(fProtonEfficiencies) b->Add(fProtonEfficiencies, "Proton Efficiencies");  
190     if(fKaonEfficiencies) b->Add(fKaonEfficiencies, "Kaon Efficiencies");
191     if(fThresholds) b->Add(fThresholds, "Thresholds");
192   }
193 }
194
195 //__________________________________________________________________
196 void AliHFEtrdPIDqa::Init(){
197   //
198   // Initialize Object
199   //
200   
201   fHistos = new AliHFEcollection("TRDqa", "Histos for TRD QA");
202
203   CreateLikelihoodHistogram();
204   CreateQAHistogram();
205   CreatedEdxHistogram();
206   CreateHistoTruncatedMean();
207
208   fTRDpid = new AliHFEpidTRD("QAtrdPID");
209 }
210
211 //__________________________________________________________________
212 void AliHFEtrdPIDqa::CreateLikelihoodHistogram(){
213   //
214   // Create Histogram for TRD Likelihood Studies
215   //
216   Int_t nbins[kQuantitiesLike]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
217   nbins[kElectronLike] = 100;
218   Double_t binMin[kQuantitiesLike]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
219   Double_t binMax[kQuantitiesLike]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
220   binMin[kElectronLike] = 0.;      
221   binMax[kElectronLike] = 1.;
222
223   fHistos->CreateTHnSparse("fLikeTRD","TRD Likelihood Studies", kQuantitiesLike, nbins, binMin, binMax);
224   fHistos->BinLogAxis("fLikeTRD", kP);
225 }
226
227 //__________________________________________________________________
228 void AliHFEtrdPIDqa::CreateQAHistogram(){
229   //
230   // Create Histogram for Basic TRD PID QA
231   //
232   AliDebug(1, "Called");
233   Int_t nbins[kQuantitiesQA]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
234   nbins[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1;
235   nbins[kNClusters] = 200;
236   Double_t binMin[kQuantitiesQA]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
237   binMin[kNonZeroTrackletCharge] = 0.;      
238   binMin[kNClusters] = 0.;      
239   Double_t binMax[kQuantitiesQA]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
240   binMax[kNonZeroTrackletCharge] = AliESDtrack::kTRDnPlanes + 1.;
241   binMax[kNClusters] = 200.;
242
243   fHistos->CreateTHnSparse("fQAtrack","TRD QA Histogram", kQuantitiesQA, nbins, binMin, binMax);
244   fHistos->BinLogAxis("fQAtrack", kP);
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   binMin[kNonZeroSlices] = 0.;
261   Double_t binMax[kQuantitiesdEdx]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
262   binMax[kdEdx] = 100000.;
263   binMax[kNclusters] = 260.;
264   binMax[kNonZeroSlices] = 8.;
265
266   fHistos->CreateTHnSparse("fQAdEdx","TRD summed dEdx", kQuantitiesdEdx, nbins, binMin, binMax);
267   fHistos->BinLogAxis("fQAdEdx", kP);
268   fHistos->Sumw2("fQAdEdx");
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   fHistos->CreateTHnSparse("fTRDtruncMean","TRD TruncatedMean studies", kQuantitiesTruncMean, nbins, binMin, binMax);
291   fHistos->BinLogAxis("fTRDtruncMean", kP);
292   fHistos->CreateTH2F("fTRDslicesPions","TRD dEdx per slice for Pions", 8, 0, 8, 500, 0, 2000);
293   fHistos->CreateTH2F("fTRDslicesElectrons","TRD dEdx per slice for Electrons", 8, 0, 8, 500, 0, 2000);
294 }
295
296
297 //__________________________________________________________________
298 void AliHFEtrdPIDqa::ProcessTracks(const TObjArray * const tracks, Int_t species){
299   //
300   // Process Electron Tracks
301   //
302   if(species < -1 || species >= AliPID::kSPECIES) return;
303   TIter it(tracks);
304   const AliVTrack *track = NULL;
305   while((track = dynamic_cast<const AliVTrack *>(it()))){
306     if(track) ProcessTrack(track, species);
307   }
308 }
309
310 //__________________________________________________________________
311 void AliHFEtrdPIDqa::ProcessTrack(const AliVTrack * const track, Int_t species){
312   //
313   // Process Single Track
314   //
315   if(TString(track->IsA()->GetName()).CompareTo("AliESDtrack") == 0)
316     ProcessTrackESD(dynamic_cast<const AliESDtrack *>(track), species);
317   else if(TString(track->IsA()->GetName()).CompareTo("AliAODTrack") == 0)
318     ProcessTrackAOD(dynamic_cast<const AliAODTrack *>(track), species);
319 }
320
321
322 //__________________________________________________________________
323 void AliHFEtrdPIDqa::ProcessTrackESD(const AliESDtrack *track, Int_t species){
324   //
325   // Process single ESD track
326   //
327   if(!track) return;
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(const AliAODTrack * const track, Int_t /*species*/){
335   //
336   // Process single AOD track
337   // AOD PID object required
338   //
339   if(!track) return;
340   AliAODPid *trackPID = track->GetDetPid();
341   if(!trackPID) return;
342
343 }
344
345 //__________________________________________________________________
346 void AliHFEtrdPIDqa::FillTRDLikelihoods(const AliESDtrack * const track, Int_t species){
347   //
348   // Fill TRD Likelihood Histogram
349   //
350   Double_t trdLike[AliPID::kSPECIES];
351   track->GetTRDpid(trdLike);
352   // Renormalize 
353   Double_t norm =trdLike[AliPID::kElectron]+trdLike[AliPID::kPion];
354   Double_t likeEle = norm == 0. ? 0. : trdLike[AliPID::kElectron]/norm;
355   const AliExternalTrackParam *outerPars = track->GetOuterParam();
356
357   Double_t quantities[kQuantitiesLike]; memset(quantities, 0, sizeof(Double_t) * kQuantitiesLike);
358   // we store:
359   // species
360   // p
361   // ntracklets
362   // Electron Likelihood
363   quantities[kSpecies] = species;
364   quantities[kP] = outerPars ? outerPars->P() : track->P();
365   quantities[kNTracklets] = track->GetTRDntrackletsPID();
366   quantities[kElectronLike] = likeEle;
367
368   fHistos->Fill("fLikeTRD", quantities);
369 }
370
371 //__________________________________________________________________
372 void AliHFEtrdPIDqa::FillTRDQAplots(const AliESDtrack * const track, Int_t species){
373   //
374   // Fill QA Plots containing further information
375   //
376   const AliExternalTrackParam *outerPars = track->GetOuterParam();
377
378   Double_t quantitiesQA[kQuantitiesQA], quantitiesdEdx[kQuantitiesdEdx], quantitiesTruncMean[kQuantitiesTruncMean];
379   // we store:
380   // 1. QA
381   // species
382   // p
383   // ntracklets
384   // Non-zero tracklet charges
385   // Number of clusters / full track
386   // 2. dEdx
387   // species
388   // p
389   // ntracklets
390   // dEdx
391   // 3. Truncated Mean
392   // ...
393   // TPC dEdx
394   // TRD dEdx Method 1
395   // TRD dEdx Method 2
396   quantitiesQA[kSpecies]  = quantitiesdEdx[kSpecies] 
397                           = quantitiesTruncMean[kSpecies] 
398                           = species;
399   quantitiesQA[kP]  = quantitiesTruncMean[kP]
400                     = outerPars ? outerPars->P() : track->P();
401   quantitiesQA[kNTracklets] = quantitiesdEdx[kNTracklets] 
402                             = quantitiesTruncMean[kNTracklets]
403                             = track->GetTRDntrackletsPID();
404   quantitiesQA[kNClusters] = quantitiesdEdx[kNclusters] = track->GetTRDncls();
405   
406
407   Double_t dEdxSum = 0., qSlice = 0.;
408   // remove the last slice from the histogram
409   Int_t ntrackletsNonZero = 0, nSlices = track->GetNumberOfTRDslices(), nSlicesNonZero = 0;
410   TString speciesname = "pions";
411   Bool_t selectedForSlicemon = kFALSE;
412   
413   switch(species){
414     case AliPID::kElectron: speciesname = "Electrons"; selectedForSlicemon = kTRUE; break;
415     case AliPID::kPion: speciesname = "Pions"; selectedForSlicemon = kTRUE; break;
416     default: speciesname = "undefined"; selectedForSlicemon = kFALSE; break;
417   };
418   AliDebug(1, Form("species %d, speciesname %s, momentum %f, selected %s", species, speciesname.Data(), track->P(), selectedForSlicemon ? "yes" : "no"));
419   for(Int_t iplane = 0; iplane < AliESDtrack::kTRDnPlanes; iplane++){
420     quantitiesdEdx[kP] = track->GetTRDmomentum(iplane);
421     dEdxSum = 0.;
422     for(Int_t islice = 0; islice < nSlices; islice++){
423       qSlice = track->GetTRDslice(iplane, islice);
424       if(qSlice > 1e-1){
425         // cut out 0 slices
426         nSlicesNonZero++;
427         dEdxSum += qSlice;
428         // Reweighting of the slices for the truncated mean: select all pion tracks above
429         // 1.5 GeV and monitor the dEdx as function of slice
430         if(selectedForSlicemon && track->P() > 1.5){
431           AliDebug(2, Form("Filling Histogram fTRDslices%s", speciesname.Data()));
432           fHistos->Fill(Form("fTRDslices%s", speciesname.Data()), static_cast<Double_t>(islice), qSlice);
433         }
434       }
435     }
436     quantitiesdEdx[kNonZeroSlices] = nSlicesNonZero;
437     quantitiesdEdx[kdEdx] = dEdxSum;
438     if(dEdxSum) ntrackletsNonZero++;
439     // Fill dEdx histogram
440     if(dEdxSum > 1e-1) fHistos->Fill("fQAdEdx", quantitiesdEdx); // Cut out 0 entries
441   }
442   quantitiesQA[kNonZeroTrackletCharge] = ntrackletsNonZero;
443   fHistos->Fill("fQAtrack", quantitiesQA);
444
445   quantitiesTruncMean[kTPCdEdx] = track->GetTPCsignal();
446   quantitiesTruncMean[kTRDdEdxMethod1] = fTRDpid->GetTRDSignalV1(track, 0.6);
447   quantitiesTruncMean[kTRDdEdxMethod2] = fTRDpid->GetTRDSignalV2(track, 0.6);
448   fHistos->Fill("fTRDtruncMean", quantitiesTruncMean);
449 }
450
451 /////////////////////////////////////////////////////////
452 //
453 // Code for Post Processing
454 //
455 // //////////////////////////////////////////////////////
456
457 //__________________________________________________________________
458 void AliHFEtrdPIDqa::FinishAnalysis(){
459   //
460   // Finish Analysis:
461   // Calculate Electron Efficiency for ntracklets = 4...6
462   // Calculate thresholds for ntracklets = 4...6
463   //
464   
465   if(!fPionEfficiencies){
466     fPionEfficiencies = new TList;
467     fPionEfficiencies->SetName("pionEfficiencies");
468   }
469   if(!fProtonEfficiencies){
470     fProtonEfficiencies = new TList;
471     fProtonEfficiencies->SetName("protonEfficiencies");
472   }
473   if(!fThresholds){
474     fThresholds = new TList;
475     fThresholds->SetName("thresholds");
476   }
477
478   for(Int_t itr = 4; itr <= 6; itr++){
479     if(fShowMessage){
480       printf("========================================\n");
481       printf("Analysing %d trackltes\n", itr);
482       printf("========================================\n");
483     }
484     AnalyseNTracklets(itr);
485   }
486 }
487
488 //__________________________________________________________________
489 void AliHFEtrdPIDqa::StoreResults(const Char_t *filename){
490   //
491   // Store histos into a root file
492   //
493   TFile *outfile = new TFile(filename, "RECREATE");
494   outfile->cd();
495   fPionEfficiencies->Clone()->Write(fPionEfficiencies->GetName(), kSingleKey);
496   fProtonEfficiencies->Clone()->Write(fProtonEfficiencies->GetName(), kSingleKey);
497   fThresholds->Clone()->Write(fThresholds->GetName(), kSingleKey);
498   outfile->Close();
499   delete outfile;
500 }
501
502 //__________________________________________________________________
503 void AliHFEtrdPIDqa::SaveThresholdParameters(const Char_t *name){
504   //
505   // Fit the threshold histograms with the given parametrisation
506   // and store the TF1 in the file
507   //
508
509   if(!fThresholds){
510     AliError("Threshold histograms have to be created first");
511     return;
512   }
513
514     if(fShowMessage){
515     printf("========================================\n");
516     printf("Calculating threshold parameters\n");
517     printf("========================================\n");
518   }
519
520   TList *outlist = new TList;
521   outlist->SetName("thresholdTRD");
522   
523   TGraph *threshhist = NULL;
524   TF1 *threshparam = NULL;
525   TList *lHistos = NULL, *lFormulas = NULL;
526   for(Int_t itracklet = 4; itracklet <= 6; itracklet++){
527   
528     if(fShowMessage){
529       printf("-------------------------------\n");
530       printf("Processing %d tracklets\n", itracklet);
531       printf("-------------------------------\n");
532     }
533
534     lHistos = dynamic_cast<TList *>(fThresholds->FindObject(Form("%dTracklets", itracklet)));
535     if(!lHistos){
536       AliError(Form("Threshold histograms for the case %d tracklets not found", itracklet));
537       continue;
538     }
539     lFormulas = new TList;
540     lFormulas->SetName(Form("%dTracklets", itracklet));
541     outlist->Add(lFormulas);
542     
543     for(Int_t ieff = 0; ieff <  kNElectronEffs; ieff++){
544       
545       if(fShowMessage){
546         printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n");
547         printf("Processing Electron Efficiency %f\n", fgkElectronEff[ieff]);
548         printf(">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n");
549       }
550
551       threshhist = dynamic_cast<TGraph *>(lHistos->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
552       if(!threshhist) continue;
553       threshparam = MakeThresholds(threshhist);
554       threshparam->SetName(Form("thresh_%d_%d", itracklet, static_cast<Int_t>(fgkElectronEff[ieff] * 100)));
555       lFormulas->Add(threshparam);
556     }
557   }
558
559   // store the output
560   TFile *outfile = new TFile(name, "RECREATE");
561   outfile->cd();
562   outlist->Write(outlist->GetName(), kSingleKey);
563   outfile->Close();
564   delete outfile;
565 }
566
567 //__________________________________________________________________
568 void AliHFEtrdPIDqa::AnalyseNTracklets(Int_t nTracklets){
569   //
570   // Calculate Pion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
571   // elPion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
572   // electron efficiencies
573   //
574   THnSparse *hLikeTRD = dynamic_cast<THnSparseF *>(fHistos->Get("fLikeTRD"));
575   if(!hLikeTRD){
576     AliError("Likelihood Histogram not available");
577     return;
578   }
579   Int_t binTracklets = hLikeTRD->GetAxis(kNTracklets)->FindBin(nTracklets);
580   hLikeTRD->GetAxis(kNTracklets)->SetRange(binTracklets, binTracklets);
581   
582   Int_t binElectrons = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kElectron); 
583   AliDebug(1, Form("BinElectrons %d", binElectrons));
584   Int_t binPions = hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kPion);
585   AliDebug(1, Form("BinPions %d", binPions));
586   Int_t binProtons =  hLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kProton);
587   AliDebug(1, Form("BinProtons %d", binProtons));
588   hLikeTRD->GetAxis(kSpecies)->SetRange(binElectrons, binElectrons);
589   TH2 *likeElectron = hLikeTRD->Projection(kElectronLike, kP);
590   likeElectron->SetName("likeElectron");
591   hLikeTRD->GetAxis(kSpecies)->SetRange(binPions, binPions);
592   TH2 *likePion = hLikeTRD->Projection(kElectronLike, kP);
593   likePion->SetName("likePion");
594   hLikeTRD->GetAxis(kSpecies)->SetRange(binProtons, binProtons);
595   TH2 *likeProton = hLikeTRD->Projection(kElectronLike, kP);
596   likeProton->SetName("likeProton");
597   
598   // Undo ranges
599   hLikeTRD->GetAxis(kSpecies)->SetRange(0, hLikeTRD->GetAxis(kSpecies)->GetNbins());
600   hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kNTracklets)->GetNbins());
601
602   // Prepare List for output
603   TList *listPions = new TList; listPions->SetName(Form("%dTracklets", nTracklets));
604   TList *listProtons = new TList; listProtons->SetName(Form("%dTracklets", nTracklets));
605   TList *listThresholds = new TList; listThresholds->SetName(Form("%dTracklets", nTracklets));
606   fPionEfficiencies->Add(listPions);
607   fProtonEfficiencies->Add(listProtons);
608   fThresholds->Add(listThresholds);
609
610   TH1 *probsEl = NULL, *probsPi = NULL, *probsPr = NULL;
611   TGraphErrors *effPi = NULL, *effPr = NULL; TGraph *thresholds = NULL;
612   Double_t p = 0, dp = 0;
613   Int_t threshbin = 0;
614   Double_t noElEff[2]; // value and error
615   for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
616     
617     if(fShowMessage){
618       printf("-----------------------------------------\n");
619       printf("Doing Electron Efficiency %f\n", fgkElectronEff[ieff]);
620       printf("-----------------------------------------\n");
621     }
622     effPi = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
623     effPi->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
624     effPr = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
625     effPr->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
626     thresholds = new TGraph(likeElectron->GetXaxis()->GetNbins());
627     thresholds->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
628   
629     // Add to lists
630     listPions->Add(effPi);
631     listProtons->Add(effPr);
632     listThresholds->Add(thresholds);
633
634     for(Int_t imom = 1; imom <= likeElectron->GetXaxis()->GetLast(); imom++){
635       p = likeElectron->GetXaxis()->GetBinCenter(imom);
636       dp = likeElectron->GetXaxis()->GetBinWidth(imom)/2;
637
638       probsEl = likeElectron->ProjectionY("el", imom);
639       if(!probsEl->GetEntries()) continue;
640       probsEl->Scale(1./probsEl->Integral());
641       probsPi = likePion->ProjectionY("pi", imom);
642       if(!probsPi->GetEntries()) continue;
643       probsPi->Scale(1./probsPi->Integral());
644       probsPr = likeProton->ProjectionY("pr", imom);
645       if(!probsPr->GetEntries()) continue;
646       probsPr->Scale(1./probsPr->Integral());
647       AliDebug(1, Form("Calculating Values for p = %f", p));
648
649       // Calculare threshold we need to achive the x% electron Efficiency
650       threshbin = GetThresholdBin(probsEl, fgkElectronEff[ieff]);
651       thresholds->SetPoint(imom - 1, p, probsEl->GetXaxis()->GetBinCenter(threshbin));
652       AliDebug(1, Form("threshold %d|%f", threshbin, probsEl->GetXaxis()->GetBinCenter(threshbin)));
653
654       // Calculate non-electronEfficiency and error
655       CalculateEfficiency(probsPi, threshbin, noElEff);
656       AliDebug(1, Form("Pion Efficiency %f", noElEff[0]));
657       effPi->SetPoint(imom - 1, p, noElEff[0]);
658       effPi->SetPointError(imom - 1, dp, noElEff[1]);
659       CalculateEfficiency(probsPr, threshbin, noElEff);
660       effPr->SetPoint(imom - 1, p, noElEff[0]);
661       effPr->SetPointError(imom - 1, dp, noElEff[1]);
662       AliDebug(1, Form("Proton Efficiency %f", noElEff[0]));
663  
664       // cleanup
665       delete probsEl;
666       delete probsPi;
667       delete probsPr;
668     }
669   }
670
671   // remove temporary histograms
672   delete likeElectron;
673   delete likePion;
674   delete likeProton;
675 }
676
677 //__________________________________________________________________
678 Int_t AliHFEtrdPIDqa::GetThresholdBin(const TH1 * const input, Double_t eff){
679   //
680   // Calculate the threshold bin  
681   //
682   Double_t integralEff = 0.;
683   Int_t currentBin = 0;
684   for(Int_t ibin = input->GetXaxis()->GetLast(); ibin >= input->GetXaxis()->GetFirst(); ibin--){
685     currentBin = ibin;
686     integralEff += input->GetBinContent(ibin);
687     if(integralEff >= eff){
688       // we found the matching bin, break the loop
689       break;
690     }
691   }
692   return currentBin;
693 }
694
695 //__________________________________________________________________
696 Bool_t AliHFEtrdPIDqa::CalculateEfficiency(const TH1 * const input, Int_t threshbin, Double_t * const par){
697   // 
698   // Calculate non-electron efficiency
699   //
700   Double_t integralEff = 0; 
701   for(Int_t ibin = threshbin; ibin <= input->GetXaxis()->GetLast(); ibin++) 
702     integralEff += input->GetBinContent(ibin);
703   par[0] = integralEff;
704
705   // @TODO: Error calculation
706   par[1] = 0;
707
708   return kTRUE;
709 }
710
711 //__________________________________________________________________
712 void AliHFEtrdPIDqa::DrawTracklet(Int_t itracklet, Double_t pmin, Double_t pmax, Bool_t doFit){
713   //
714   // Draw efficiencies and threshold as function of p
715   //
716   if(!(fPionEfficiencies && fProtonEfficiencies && fThresholds)){
717     AliError("No graphs to draw available");  
718     return;
719   }
720
721   TList *lpions = dynamic_cast<TList *>(fPionEfficiencies->FindObject(Form("%dTracklets", itracklet))); 
722   TList *lprotons = dynamic_cast<TList *>(fProtonEfficiencies->FindObject(Form("%dTracklets", itracklet))); 
723   
724   TList *lthresholds = dynamic_cast<TList *>(fThresholds->FindObject(Form("%dTracklets", itracklet))); 
725   if(!(lpions && lprotons && lthresholds)){
726     AliDebug(1, "Relevant lists not found. Did you forget to run FinishAnalysis()?");
727     return;
728   }
729
730   TGraphErrors *pi, *pr;
731   TGraph *tr;
732   TLegend *leg;
733   TCanvas *c1 = new TCanvas(Form("tracklet%d", itracklet), Form("Tracklet %d", itracklet), 1024, 768);
734   c1->Divide(3,2);
735   TF1 *threshfit = NULL;
736   for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
737     c1->cd(ieff + 1);
738     leg = new TLegend(0.6, 0.7, 0.89, 0.89);
739     leg->SetBorderSize(0);
740     leg->SetFillStyle(0);
741     pi = dynamic_cast<TGraphErrors *>(lpions->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
742     pr = dynamic_cast<TGraphErrors *>(lprotons->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
743     tr = dynamic_cast<TGraph *>(lthresholds->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
744     if(!(pi && pr && tr)) continue;
745
746     // Define nice plot, draw
747     // Axis Title
748     pi->GetXaxis()->SetTitle("p / GeV/c");
749     pi->GetYaxis()->SetTitle("Efficiency");
750     pr->GetXaxis()->SetTitle("p / GeV/c");
751     pr->GetYaxis()->SetTitle("Efficiency");
752     tr->GetXaxis()->SetTitle("p / GeV/c");
753     tr->GetYaxis()->SetTitle("Efficiency");
754     // Axis Range
755     pi->GetYaxis()->SetRangeUser(1e-3, 1.);
756     pr->GetYaxis()->SetRangeUser(1e-3, 1.);
757     tr->GetYaxis()->SetRangeUser(1e-3, 1.);
758     if(pmin > 0 && pmax > 0.){
759       pi->GetXaxis()->SetRangeUser(pmin, pmax);
760       pr->GetXaxis()->SetRangeUser(pmin, pmax);
761       tr->GetXaxis()->SetRangeUser(pmin, pmax);
762     }
763     // Marker
764     pi->SetMarkerColor(kRed);
765     pi->SetMarkerStyle(20);
766     pr->SetMarkerColor(kBlue);
767     pr->SetMarkerStyle(21);
768     tr->SetMarkerColor(kBlack);
769     tr->SetMarkerStyle(22);
770     // Title
771     pi->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
772     pr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
773     tr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
774     // Draw
775     pi->Draw("ape"); pr->Draw("pesame"); tr->Draw("psame");
776
777     // Optionally do Fit
778     if(doFit){
779       threshfit = MakeThresholds(tr);
780       threshfit->SetLineColor(kBlack);
781       threshfit->Draw("same");
782     }
783
784     // Add entries to legend
785     leg->AddEntry(pi, "Pion Efficiency", "lp");
786     leg->AddEntry(pr, "Proton Efficiency", "lp");
787     leg->AddEntry(tr, "Thresholds", "lp");
788     leg->Draw();
789     c1->Update();
790   }
791 }
792
793 //__________________________________________________________________
794 TF1 *AliHFEtrdPIDqa::MakeThresholds(TGraph *threshist){
795   //
796   // Create TF1 containing the threshold parametrisation
797   //
798
799   TF1 *threshparam = new TF1("thresh", "1-[0]-[1]*x-[2]*TMath::Exp(-[3]*x)", 0.1, 10);
800   threshist->Fit(threshparam, "NE", "", 0.1, 3.5);
801   return threshparam;
802 }
803
804 //__________________________________________________________________
805 void AliHFEtrdPIDqa::ClearLists(){
806   //
807   // Clear lists for particle efficiencies and thresholds
808   //
809   if(fPionEfficiencies){
810     fPionEfficiencies->Delete();
811     delete fPionEfficiencies;
812     fPionEfficiencies = NULL;
813   }
814   if(fProtonEfficiencies){
815     fProtonEfficiencies->Delete();
816     delete fProtonEfficiencies;
817     fProtonEfficiencies = NULL;
818   }
819   if(fThresholds){
820     fThresholds->Delete();
821     delete fThresholds;
822     fThresholds = NULL;
823   }
824 }
825
826 //__________________________________________________________________
827 Double_t AliHFEtrdPIDqa::EvalPionEfficiency(Int_t ntls, Int_t eEff, Double_t p){
828   //
829   // calculate pion efficiency
830   // Arguments:
831   //   Number of tracklets
832   //   Electron Efficiency
833   //   Momentum
834   //
835   TList *graphs = dynamic_cast<TList *>(fPionEfficiencies->FindObject(Form("%dTracklets", ntls)));
836   if(!graphs) return -1.;
837   TGraph *measurement = dynamic_cast<TGraph *>(graphs->FindObject(Form("eff%d", eEff)));
838   if(!measurement) return -1.;
839   return measurement->Eval(p);
840 }
841
842 //__________________________________________________________________
843 Double_t AliHFEtrdPIDqa::EvalProtonEfficiency(Int_t ntls, Int_t eEff, Double_t p){
844   //
845   // calculate proton efficiency
846   // Arguments:
847   //   Number of tracklets
848   //   Electron Efficiency
849   //   Momentum
850   //
851   TList *graphs = dynamic_cast<TList *>(fProtonEfficiencies->FindObject(Form("%dTracklets", ntls)));
852   if(!graphs) return -1.;
853   TGraph *measurement = dynamic_cast<TGraph *>(graphs->FindObject(Form("eff%d", eEff)));
854   if(!measurement) return -1.;
855   return measurement->Eval(p);
856 }
857
858 //__________________________________________________________________
859 Double_t AliHFEtrdPIDqa::EvalThreshold(Int_t ntls, Int_t eEff, Double_t p){
860   //
861   // Get the threshold to reach a certain electron efficiency
862   // Arguments:
863   //   Number of tracklets
864   //   Electron Efficiency
865   //   Momentum
866   //
867   TList *graphs = dynamic_cast<TList *>(fThresholds->FindObject(Form("%dTracklets", ntls)));
868   if(!graphs) return -1.;
869   TGraph *measurement = dynamic_cast<TGraph *>(graphs->FindObject(Form("eff%d", eEff)));
870   if(!measurement) return -1.;
871   return measurement->Eval(p);
872 }
873