]>
Commit | Line | Data |
---|---|---|
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 | ||
53 | ClassImp(AliHFEtrdPIDqa) | |
54 | ||
55 | const Double_t AliHFEtrdPIDqa::fgkElectronEff[kNElectronEffs] = {0.7,0.75, 0.8, 0.85, 0.9, 0.95}; | |
56 | //_______________________________________________________________ | |
57 | // Definition of the common binning | |
58 | const Int_t AliHFEtrdPIDqa::fgkNBinsCommon[kQuantitiesCommon] = { | |
59 | AliPID::kSPECIES + 1, // species | |
60 | 40, // p-bins | |
61 | AliESDtrack::kTRDnPlanes + 1 // tracklets including 0 | |
62 | }; | |
63 | const Double_t AliHFEtrdPIDqa::fgkMinBinCommon[kQuantitiesCommon] = { | |
64 | -1, // species | |
65 | 0.1, // p-bins | |
66 | 0 // tracklets including 0 | |
67 | }; | |
68 | ||
69 | const Double_t AliHFEtrdPIDqa::fgkMaxBinCommon[kQuantitiesCommon] = { | |
70 | AliPID::kSPECIES, // species | |
71 | 10., // p-bins | |
72 | AliESDtrack::kTRDnPlanes + 1 // tracklets including 0 | |
73 | }; | |
74 | //_______________________________________________________________ | |
75 | ||
76 | //__________________________________________________________________ | |
77 | AliHFEtrdPIDqa::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 | //__________________________________________________________________ | |
93 | AliHFEtrdPIDqa::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 | //__________________________________________________________________ | |
109 | AliHFEtrdPIDqa::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 | //__________________________________________________________________ | |
126 | AliHFEtrdPIDqa &AliHFEtrdPIDqa::operator=(const AliHFEtrdPIDqa &ref){ | |
127 | // | |
128 | // Assignment operator | |
129 | // | |
130 | if(this != &ref) | |
131 | ref.Copy(*this); | |
132 | return *this; | |
133 | } | |
134 | ||
135 | //__________________________________________________________________ | |
136 | AliHFEtrdPIDqa::~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 | //__________________________________________________________________ | |
148 | void 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 | //__________________________________________________________________ | |
160 | Long64_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 | ||
bf892a6a | 183 | //__________________________________________________________________ |
184 | void 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 | ||
c04c80e6 | 198 | //__________________________________________________________________ |
199 | void 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 | //__________________________________________________________________ | |
215 | void 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 | //__________________________________________________________________ | |
231 | void 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 | //__________________________________________________________________ | |
251 | void 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 | //__________________________________________________________________ | |
275 | void 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 | //__________________________________________________________________ | |
301 | void 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 | //__________________________________________________________________ | |
315 | void 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 | //__________________________________________________________________ | |
327 | void 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 | //__________________________________________________________________ | |
338 | void 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 | //__________________________________________________________________ | |
350 | void 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 | //__________________________________________________________________ | |
373 | void 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 | //__________________________________________________________________ | |
459 | void 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 | //__________________________________________________________________ | |
490 | void 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 | //__________________________________________________________________ | |
504 | void 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 | //__________________________________________________________________ | |
569 | void 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 | //__________________________________________________________________ | |
679 | Int_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 | //__________________________________________________________________ | |
697 | Bool_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 | 713 | void 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 | //__________________________________________________________________ | |
795 | TF1 *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 | //__________________________________________________________________ | |
806 | void 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 | //__________________________________________________________________ | |
828 | Double_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 | //__________________________________________________________________ | |
837 | Double_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 | //__________________________________________________________________ | |
846 | Double_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 |