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