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