]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/hfe/AliHFEtrdPIDqa.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Add missing classes (first try)
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEtrdPIDqa.cxx
CommitLineData
c04c80e6 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
49ClassImp(AliHFEtrdPIDqa)
50
51const Double_t AliHFEtrdPIDqa::fgkElectronEff[kNElectronEffs] = {0.7,0.75, 0.8, 0.85, 0.9, 0.95};
52//_______________________________________________________________
53// Definition of the common binning
54const Int_t AliHFEtrdPIDqa::fgkNBinsCommon[kQuantitiesCommon] = {
55 AliPID::kSPECIES + 1, // species
56 40, // p-bins
57 AliESDtrack::kTRDnPlanes + 1 // tracklets including 0
58};
59const Double_t AliHFEtrdPIDqa::fgkMinBinCommon[kQuantitiesCommon] = {
60 -1, // species
61 0.1, // p-bins
62 0 // tracklets including 0
63};
64
65const Double_t AliHFEtrdPIDqa::fgkMaxBinCommon[kQuantitiesCommon] = {
66 AliPID::kSPECIES, // species
67 10., // p-bins
68 AliESDtrack::kTRDnPlanes + 1 // tracklets including 0
69};
70//_______________________________________________________________
71
72//__________________________________________________________________
73AliHFEtrdPIDqa::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//__________________________________________________________________
91AliHFEtrdPIDqa::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//__________________________________________________________________
109AliHFEtrdPIDqa::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//__________________________________________________________________
128AliHFEtrdPIDqa &AliHFEtrdPIDqa::operator=(const AliHFEtrdPIDqa &ref){
129 //
130 // Assignment operator
131 //
132 if(this != &ref)
133 ref.Copy(*this);
134 return *this;
135}
136
137//__________________________________________________________________
138AliHFEtrdPIDqa::~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//__________________________________________________________________
154void 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//__________________________________________________________________
169Long64_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//__________________________________________________________________
194void 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//__________________________________________________________________
208void 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//__________________________________________________________________
226void 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//__________________________________________________________________
248void 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 Double_t binMin[kQuantitiesdEdx]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
256 binMin[kdEdx] = 0.;
257 Double_t binMax[kQuantitiesdEdx]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
258 binMax[kdEdx] = 100000.;
259
260 fQAdEdx = new THnSparseF("fQAdEdx","TRD summed dEdx", kQuantitiesdEdx, nbins, binMin, binMax);
261 Double_t *binLog = AliHFEtools::MakeLogarithmicBinning(nbins[kP], binMin[kP], binMax[kP]);
262 fQAdEdx->GetAxis(kP)->Set(nbins[kP], binLog);
263 delete[] binLog;
264}
265
266//__________________________________________________________________
267void AliHFEtrdPIDqa::CreateHistoTruncatedMean(){
268 //
269 // Create Histogram for Basic TRD PID QA
270 //
271 AliDebug(1, "Called");
272 Int_t nbins[kQuantitiesTruncMean]; memcpy(nbins, fgkNBinsCommon, sizeof(Int_t) * kQuantitiesCommon);
273 nbins[kTPCdEdx] = 600;
274 nbins[kTRDdEdxMethod1] = 1000;
275 nbins[kTRDdEdxMethod2] = 1000;
276 Double_t binMin[kQuantitiesTruncMean]; memcpy(binMin, fgkMinBinCommon, sizeof(Double_t) * kQuantitiesCommon);
277 binMin[kTPCdEdx] = 0.;
278 binMin[kTRDdEdxMethod1] = 0.;
279 binMin[kTRDdEdxMethod2] = 0.;
280 Double_t binMax[kQuantitiesTruncMean]; memcpy(binMax, fgkMaxBinCommon, sizeof(Double_t) * kQuantitiesCommon);
281 binMax[kTPCdEdx] = 600;
282 binMax[kTRDdEdxMethod1] = 20000.;
283 binMax[kTRDdEdxMethod2] = 20000.;
284
285 fTRDtruncMean = new THnSparseF("fTRDtruncMean","TRD TruncatedMean studies", kQuantitiesTruncMean, nbins, binMin, binMax);
286 Double_t *binLog = AliHFEtools::MakeLogarithmicBinning(nbins[kP], binMin[kP], binMax[kP]);
287 fTRDtruncMean->GetAxis(kP)->Set(nbins[kP], binLog);
288 delete[] binLog;
289}
290
291
292//__________________________________________________________________
293void AliHFEtrdPIDqa::ProcessTracks(TObjArray * const tracks, Int_t species){
294 //
295 // Process Electron Tracks
296 //
297 if(species < -1 || species >= AliPID::kSPECIES) return;
298 TIterator *it = tracks->MakeIterator();
299 AliVTrack *track = NULL;
300 while((track = dynamic_cast<AliVTrack *>(it->Next()))){
301 if(track) ProcessTrack(track, species);
302 }
303 delete it;
304}
305
306//__________________________________________________________________
307void AliHFEtrdPIDqa::ProcessTrack(AliVTrack *track, Int_t species){
308 //
309 // Process Single Track
310 //
311 if(TString(track->IsA()->GetName()).CompareTo("AliESDtrack") == 0)
312 ProcessTrackESD(dynamic_cast<AliESDtrack *>(track), species);
313 else if(TString(track->IsA()->GetName()).CompareTo("AliAODTrack") == 0)
314 ProcessTrackAOD(dynamic_cast<AliAODTrack *>(track), species);
315}
316
317
318//__________________________________________________________________
319void AliHFEtrdPIDqa::ProcessTrackESD(AliESDtrack *track, Int_t species){
320 //
321 // Process single ESD track
322 //
323 if((track->GetStatus() & AliESDtrack::kTRDout) == 0) return; // require TRD track
324 FillTRDLikelihoods(track, species);
325 FillTRDQAplots(track, species);
326}
327
328//__________________________________________________________________
329void AliHFEtrdPIDqa::ProcessTrackAOD(AliAODTrack * const track, Int_t /*species*/){
330 //
331 // Process single AOD track
332 // AOD PID object required
333 //
334 AliAODPid *trackPID = track->GetDetPid();
335 if(!trackPID) return;
336
337}
338
339//__________________________________________________________________
340void AliHFEtrdPIDqa::FillTRDLikelihoods(AliESDtrack *track, Int_t species){
341 //
342 // Fill TRD Likelihood Histogram
343 //
344 Double_t trdLike[AliPID::kSPECIES];
345 track->GetTRDpid(trdLike);
346 const AliExternalTrackParam *outerPars = track->GetOuterParam();
347
348 Double_t quantities[kQuantitiesLike]; memset(quantities, 0, sizeof(Double_t) * kQuantitiesLike);
349 // we store:
350 // species
351 // p
352 // ntracklets
353 // Electron Likelihood
354 quantities[kSpecies] = species;
355 quantities[kP] = outerPars ? outerPars->P() : track->P();
356 quantities[kNTracklets] = track->GetTRDntrackletsPID();
357 quantities[kElectronLike] = trdLike[AliPID::kElectron];
358
359 fLikeTRD->Fill(quantities);
360}
361
362//__________________________________________________________________
363void AliHFEtrdPIDqa::FillTRDQAplots(AliESDtrack *track, Int_t species){
364 //
365 // Fill QA Plots containing further information
366 //
367 const AliExternalTrackParam *outerPars = track->GetOuterParam();
368
369 Double_t quantitiesQA[kQuantitiesQA], quantitiesdEdx[kQuantitiesdEdx], quantitiesTruncMean[kQuantitiesTruncMean];
370 // we store:
371 // 1. QA
372 // species
373 // p
374 // ntracklets
375 // Non-zero tracklet charges
376 // Number of clusters / full track
377 // 2. dEdx
378 // species
379 // p
380 // ntracklets
381 // dEdx
382 // 3. Truncated Mean
383 // ...
384 // TPC dEdx
385 // TRD dEdx Method 1
386 // TRD dEdx Method 2
387 quantitiesQA[kSpecies] = quantitiesdEdx[kSpecies]
388 = quantitiesTruncMean[kSpecies]
389 = species;
390 quantitiesQA[kP] = quantitiesdEdx[kP]
391 = quantitiesTruncMean[kP]
392 = outerPars ? outerPars->P() : track->P();
393 quantitiesQA[kNTracklets] = quantitiesdEdx[kNTracklets]
394 = quantitiesTruncMean[kNTracklets]
395 = track->GetTRDntrackletsPID();
396 quantitiesQA[kNClusters] = track->GetTRDncls();
397
398 Double_t dEdxSum;
399 Int_t ntrackletsNonZero = 0, nSlices = track->GetNumberOfTRDslices();
400 for(Int_t iplane = 0; iplane < AliESDtrack::kTRDnPlanes; iplane++){
401 dEdxSum = 0.;
402 for(Int_t islice = 0; islice < nSlices; islice++)
403 dEdxSum += track->GetTRDslice(iplane, islice);
404 quantitiesdEdx[kdEdx] = dEdxSum;
405 if(dEdxSum) ntrackletsNonZero++;
406 // Fill dEdx histogram
407 fQAdEdx->Fill(quantitiesdEdx);
408 }
409 quantitiesQA[kNonZeroTrackletCharge] = ntrackletsNonZero;
410 fQAtrack->Fill(quantitiesQA);
411
412 quantitiesTruncMean[kTPCdEdx] = track->GetTPCsignal();
413 quantitiesTruncMean[kTRDdEdxMethod1] = fTRDpid->GetTRDSignalV1(track, -1);
414 quantitiesTruncMean[kTRDdEdxMethod2] = fTRDpid->GetTRDSignalV2(track, -1);
415 fTRDtruncMean->Fill(quantitiesTruncMean);
416}
417
418/////////////////////////////////////////////////////////
419//
420// Code for Post Processing
421//
422// //////////////////////////////////////////////////////
423
424//__________________________________________________________________
425void AliHFEtrdPIDqa::FinishAnalysis(){
426 //
427 // Finish Analysis:
428 // Calculate Electron Efficiency for ntracklets = 4...6
429 // Calculate thresholds for ntracklets = 4...6
430 //
431
432 if(!fPionEfficiencies){
433 fPionEfficiencies = new TList;
434 fPionEfficiencies->SetName("pionEfficiencies");
435 }
436 if(!fProtonEfficiencies){
437 fProtonEfficiencies = new TList;
438 fProtonEfficiencies->SetName("protonEfficiencies");
439 }
440 if(!fThresholds){
441 fThresholds = new TList;
442 fThresholds->SetName("thresholds");
443 }
444
445 for(Int_t itr = 4; itr <= 6; itr++){
446 printf("========================================\n");
447 printf("Analysing %d trackltes\n", itr);
448 printf("========================================\n");
449 AnalyseNTracklets(itr);
450 }
451}
452
453//__________________________________________________________________
454void AliHFEtrdPIDqa::StoreResults(const Char_t *filename){
455 //
456 // Store histos into a root file
457 //
458 TFile *outfile = new TFile(filename, "RECREATE");
459 outfile->cd();
460 fPionEfficiencies->Clone()->Write(fPionEfficiencies->GetName(), kSingleKey);
461 fProtonEfficiencies->Clone()->Write(fProtonEfficiencies->GetName(), kSingleKey);
462 fThresholds->Clone()->Write(fThresholds->GetName(), kSingleKey);
463 outfile->Close();
464 delete outfile;
465}
466
467//__________________________________________________________________
468void AliHFEtrdPIDqa::SaveThresholdParameters(const Char_t *name){
469 //
470 // Fit the threshold histograms with the given parametrisation
471 // and store the TF1 in the file
472 //
473
474 if(!fThresholds){
475 AliError("Threshold histograms have to be created first");
476 return;
477 }
478
479 printf("========================================\n");
480 printf("Calculating threshold parameters\n");
481 printf("========================================\n");
482
483 TList *outlist = new TList;
484 outlist->SetName("thresholdTRD");
485
486 TGraph *threshhist = NULL;
487 TF1 *threshparam = NULL;
488 TList *lHistos = NULL, *lFormulas = NULL;
489 for(Int_t itracklet = 4; itracklet <= 6; itracklet++){
490
491 printf("-------------------------------\n");
492 printf("Processing %d tracklets\n", itracklet);
493 printf("-------------------------------\n");
494
495 lHistos = dynamic_cast<TList *>(fThresholds->FindObject(Form("%dTracklets", itracklet)));
496 if(!lHistos){
497 AliError(Form("Threshold histograms for the case %d tracklets not found"));
498 continue;
499 }
500 lFormulas = new TList;
501 lFormulas->SetName(Form("%dTracklets", itracklet));
502 outlist->Add(lFormulas);
503
504 for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
505
506 printf("<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<\n");
507 printf("Processing Electron Efficiency %f\n", fgkElectronEff[ieff]);
508 printf(">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>\n");
509
510 threshhist = dynamic_cast<TGraph *>(lHistos->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
511 threshparam = MakeThresholds(threshhist);
512 threshparam->SetName(Form("thresh_%d_%d", itracklet, static_cast<Int_t>(fgkElectronEff[ieff] * 100)));
513 }
514 }
515
516 // store the output
517 TFile *outfile = new TFile(name, "RECREATE");
518 outfile->cd();
519 outlist->Write(outlist->GetName(), kSingleKey);
520 outfile->Close();
521 delete outfile;
522}
523
524//__________________________________________________________________
525void AliHFEtrdPIDqa::AnalyseNTracklets(Int_t nTracklets){
526 //
527 // Calculate Pion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
528 // elPion Efficiency, Proton Efficiency and Kaon Efficiency at discrete
529 // electron efficiencies
530 //
531 Int_t binTracklets = fLikeTRD->GetAxis(kNTracklets)->FindBin(nTracklets);
532 fLikeTRD->GetAxis(kNTracklets)->SetRange(binTracklets, binTracklets);
533
534 Int_t binElectrons = fLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kElectron);
535 AliDebug(1, Form("BinElectrons %d", binElectrons));
536 Int_t binPions = fLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kPion);
537 AliDebug(1, Form("BinPions %d", binPions));
538 Int_t binProtons = fLikeTRD->GetAxis(kSpecies)->FindBin(AliPID::kProton);
539 AliDebug(1, Form("BinProtons %d", binProtons));
540 fLikeTRD->GetAxis(kSpecies)->SetRange(binElectrons, binElectrons);
541 TH2 *likeElectron = fLikeTRD->Projection(kElectronLike, kP);
542 likeElectron->SetName("likeElectron");
543 fLikeTRD->GetAxis(kSpecies)->SetRange(binPions, binPions);
544 TH2 *likePion = fLikeTRD->Projection(kElectronLike, kP);
545 likePion->SetName("likePion");
546 fLikeTRD->GetAxis(kSpecies)->SetRange(binProtons, binProtons);
547 TH2 *likeProton = fLikeTRD->Projection(kElectronLike, kP);
548 likeProton->SetName("likeProton");
549
550 // Undo ranges
551 fLikeTRD->GetAxis(kSpecies)->SetRange(0, fLikeTRD->GetAxis(kSpecies)->GetNbins());
552 fLikeTRD->GetAxis(kNTracklets)->SetRange(0, fLikeTRD->GetAxis(kNTracklets)->GetNbins());
553
554 // Prepare List for output
555 TList *listPions = new TList; listPions->SetName(Form("%dTracklets", nTracklets));
556 TList *listProtons = new TList; listProtons->SetName(Form("%dTracklets", nTracklets));
557 TList *listThresholds = new TList; listThresholds->SetName(Form("%dTracklets", nTracklets));
558 fPionEfficiencies->Add(listPions);
559 fProtonEfficiencies->Add(listProtons);
560 fThresholds->Add(listThresholds);
561
562 TH1 *probsEl = NULL, *probsPi = NULL, *probsPr = NULL;
563 TGraphErrors *effPi = NULL, *effPr = NULL; TGraph *thresholds = NULL;
564 Double_t p = 0, dp = 0;
565 Int_t threshbin = 0;
566 Double_t noElEff[2]; // value and error
567 for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
568 printf("-----------------------------------------\n");
569 printf("Doing Electron Efficiency %f\n", fgkElectronEff[ieff]);
570 printf("-----------------------------------------\n");
571 effPi = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
572 effPi->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
573 effPr = new TGraphErrors(likeElectron->GetXaxis()->GetNbins());
574 effPr->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
575 thresholds = new TGraph(likeElectron->GetXaxis()->GetNbins());
576 thresholds->SetName(Form("eff%d", static_cast<Int_t >(fgkElectronEff[ieff] * 100)));
577
578 // Add to lists
579 listPions->Add(effPi);
580 listProtons->Add(effPr);
581 listThresholds->Add(thresholds);
582
583 for(Int_t imom = 1; imom <= likeElectron->GetXaxis()->GetLast(); imom++){
584 p = likeElectron->GetXaxis()->GetBinCenter(imom);
585 dp = likeElectron->GetXaxis()->GetBinWidth(imom)/2;
586
587 probsEl = likeElectron->ProjectionY("el", imom);
588 if(!probsEl->GetEntries()) continue;
589 probsEl->Scale(1./probsEl->Integral());
590 probsPi = likePion->ProjectionY("pi", imom);
591 if(!probsPi->GetEntries()) continue;
592 probsPi->Scale(1./probsPi->Integral());
593 probsPr = likeProton->ProjectionY("pr", imom);
594 if(!probsPr->GetEntries()) continue;
595 probsPr->Scale(1./probsPr->Integral());
596 AliDebug(1, Form("Calculating Values for p = %f", p));
597
598 // Calculare threshold we need to achive the x% electron Efficiency
599 threshbin = GetThresholdBin(probsEl, fgkElectronEff[ieff]);
600 thresholds->SetPoint(imom - 1, p, probsEl->GetXaxis()->GetBinCenter(threshbin));
601 AliDebug(1, Form("threshold %d|%f", threshbin, probsEl->GetXaxis()->GetBinCenter(threshbin)));
602
603 // Calculate non-electronEfficiency and error
604 CalculateEfficiency(probsPi, threshbin, noElEff);
605 AliDebug(1, Form("Pion Efficiency %f", noElEff[0]));
606 effPi->SetPoint(imom - 1, p, noElEff[0]);
607 effPi->SetPointError(imom - 1, dp, noElEff[1]);
608 CalculateEfficiency(probsPr, threshbin, noElEff);
609 effPr->SetPoint(imom - 1, p, noElEff[0]);
610 effPr->SetPointError(imom - 1, dp, noElEff[1]);
611 AliDebug(1, Form("Proton Efficiency %f", noElEff[0]));
612
613 // cleanup
614 delete probsEl;
615 delete probsPi;
616 delete probsPr;
617 }
618 }
619
620 // remove temporary histograms
621 delete likeElectron;
622 delete likePion;
623 delete likeProton;
624}
625
626//__________________________________________________________________
627Int_t AliHFEtrdPIDqa::GetThresholdBin(TH1 * const input, Double_t eff){
628 //
629 // Calculate the threshold bin
630 //
631 Double_t integralEff = 0.;
632 Int_t currentBin = 0;
633 for(Int_t ibin = input->GetXaxis()->GetLast(); ibin >= input->GetXaxis()->GetFirst(); ibin--){
634 currentBin = ibin;
635 integralEff += input->GetBinContent(ibin);
636 if(integralEff >= eff){
637 // we found the matching bin, break the loop
638 break;
639 }
640 }
641 return currentBin;
642}
643
644//__________________________________________________________________
645Bool_t AliHFEtrdPIDqa::CalculateEfficiency(TH1 * const input, Int_t threshbin, Double_t *par){
646 //
647 // Calculate non-electron efficiency
648 //
649 Double_t integralEff = 0;
650 for(Int_t ibin = threshbin; ibin <= input->GetXaxis()->GetLast(); ibin++)
651 integralEff += input->GetBinContent(ibin);
652 par[0] = integralEff;
653
654 // @TODO: Error calculation
655 par[1] = 0;
656
657 return kTRUE;
658}
659
660//__________________________________________________________________
661void AliHFEtrdPIDqa::DrawTracklet(Int_t itracklet){
662 //
663 // Draw efficiencies and threshold as function of p
664 //
665 if(!(fPionEfficiencies && fProtonEfficiencies && fThresholds)){
666 AliError("No graphs to draw available");
667 return;
668 }
669
670 TList *lpions = dynamic_cast<TList *>(fPionEfficiencies->FindObject(Form("%dTracklets", itracklet)));
671 TList *lprotons = dynamic_cast<TList *>(fProtonEfficiencies->FindObject(Form("%dTracklets", itracklet)));
672
673 TList *lthresholds = dynamic_cast<TList *>(fThresholds->FindObject(Form("%dTracklets", itracklet)));
674
675 TGraphErrors *pi, *pr;
676 TGraph *tr;
677 TLegend *leg;
678 TCanvas *c1 = new TCanvas(Form("tracklet%d", itracklet), Form("Tracklet %d", itracklet));
679 c1->Divide(3,2);
680 for(Int_t ieff = 0; ieff < kNElectronEffs; ieff++){
681 c1->cd(ieff + 1);
682 leg = new TLegend(0.6, 0.7, 0.89, 0.89);
683 leg->SetBorderSize(0);
684 leg->SetFillStyle(0);
685 pi = dynamic_cast<TGraphErrors *>(lpions->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
686 pr = dynamic_cast<TGraphErrors *>(lprotons->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
687 tr = dynamic_cast<TGraph *>(lthresholds->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
688
689 // Define nice plot, draw
690 // Axis Title
691 pi->GetXaxis()->SetTitle("p / GeV/c");
692 pi->GetYaxis()->SetTitle("Efficiency / %");
693 pr->GetXaxis()->SetTitle("p / GeV/c");
694 pr->GetYaxis()->SetTitle("Efficiency / %");
695 tr->GetXaxis()->SetTitle("p / GeV/c");
696 tr->GetYaxis()->SetTitle("Efficiency / %");
697 // Axis Range
698 pi->GetYaxis()->SetRangeUser(0., 1.);
699 pr->GetYaxis()->SetRangeUser(0., 1.);
700 tr->GetYaxis()->SetRangeUser(0., 1.);
701 // Marker
702 pi->SetMarkerColor(kRed);
703 pi->SetMarkerStyle(20);
704 pr->SetMarkerColor(kBlue);
705 pr->SetMarkerStyle(21);
706 tr->SetMarkerColor(kBlack);
707 tr->SetMarkerStyle(22);
708 // Title
709 pi->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
710 pr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
711 tr->SetTitle(Form ("%.2f Electron Efficiency", fgkElectronEff[ieff]));
712 // Draw
713 pi->Draw("ape"); pr->Draw("pesame"); tr->Draw("psame");
714
715 // Add entries to legend
716 leg->AddEntry(pi, "Pion Efficiency", "lp");
717 leg->AddEntry(pr, "Proton Efficiency", "lp");
718 leg->AddEntry(tr, "Thresholds", "lp");
719 leg->Draw();
720 c1->Update();
721 }
722}
723
724//__________________________________________________________________
725TF1 *AliHFEtrdPIDqa::MakeThresholds(TGraph *threshist){
726 //
727 // Create TF1 containing the threshold parametrisation
728 //
729
730 TF1 *threshparam = new TF1("thresh", "1-[0]-[1]*x-[2]*TMath::Exp(-[3]*x)", 0.1, 10);
731 threshist->Fit(threshparam, "NE", "", 0, 10);
732 return threshparam;
733}
734
735//__________________________________________________________________
736void AliHFEtrdPIDqa::ClearLists(){
737 //
738 // Clear lists for particle efficiencies and thresholds
739 //
740 if(fPionEfficiencies){
741 fPionEfficiencies->Delete();
742 delete fPionEfficiencies;
743 fPionEfficiencies = NULL;
744 }
745 if(fProtonEfficiencies){
746 fProtonEfficiencies->Delete();
747 delete fProtonEfficiencies;
748 fProtonEfficiencies = NULL;
749 }
750 if(fThresholds){
751 fThresholds->Delete();
752 delete fThresholds;
753 fThresholds = NULL;
754 }
755}