]>
Commit | Line | Data |
---|---|---|
70da6c5a | 1 | /************************************************************************** |
faee3b18 | 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 | **************************************************************************/ | |
70da6c5a | 15 | // |
16 | // Class for PID QA | |
17 | // Several studies done on clean samples of electrons, pions and kaons | |
18 | // coming from V0 PID | |
19 | // Compatible with both ESDs and AODs | |
20 | // | |
21 | // Autors: | |
22 | // Matus Kalisky <matus.kalisky@cern.ch> | |
23 | // Markus Heide <mheide@uni-muenster.de> | |
24 | // Markus Fasel <M.Fasel@gsi.de> | |
25 | // | |
faee3b18 | 26 | |
27 | ||
70da6c5a | 28 | #include <TMath.h> |
29 | #include <TObjArray.h> | |
70da6c5a | 30 | #include <TPDGCode.h> |
31 | #include <TString.h> | |
faee3b18 | 32 | #include <TMultiLayerPerceptron.h> |
33 | #include <TFile.h> | |
70da6c5a | 34 | |
35 | #include "AliAODMCParticle.h" | |
faee3b18 | 36 | #include "AliAODEvent.h" |
70da6c5a | 37 | #include "AliAODTrack.h" |
38 | #include "AliESDtrack.h" | |
faee3b18 | 39 | #include "AliESDEvent.h" |
70da6c5a | 40 | #include "AliMCEvent.h" |
41 | #include "AliMCParticle.h" | |
e17c1f86 | 42 | #include "AliMultiplicity.h" |
70da6c5a | 43 | #include "AliPID.h" |
8c1c76e9 | 44 | #include "AliPIDResponse.h" |
70da6c5a | 45 | #include "AliVParticle.h" |
70da6c5a | 46 | |
47 | ||
48 | #include "AliHFEcollection.h" | |
49 | #include "AliHFEpidQA.h" | |
faee3b18 | 50 | #include "AliHFEV0info.h" |
70da6c5a | 51 | #include "AliHFEV0pid.h" |
52 | #include "AliHFEV0pidMC.h" | |
3a72645a | 53 | #include "AliHFEV0cuts.h" |
faee3b18 | 54 | #include "AliHFEtrdPIDqa.h" |
55 | ||
70da6c5a | 56 | ClassImp(AliHFEpidQA) |
57 | ||
faee3b18 | 58 | //__________________________________________ |
59 | AliHFEpidQA::AliHFEpidQA(): | |
60 | fEvent(NULL), | |
61 | fMC(NULL), | |
62 | fV0pid(NULL), | |
63 | fV0pidMC(NULL), | |
64 | fTRDpidQA(NULL), | |
65 | fOutput(NULL), | |
66 | fESDpid(NULL), | |
e156c3bb | 67 | fNNref(NULL), |
e17c1f86 | 68 | fTRDTotalChargeInSlice0(kFALSE), |
69 | fIsPbPb(kFALSE), | |
70 | fIsppMultiBin(kFALSE) | |
70da6c5a | 71 | { |
72 | // | |
73 | // Default constructor | |
74 | // | |
faee3b18 | 75 | for(Int_t mom = 0; mom < 11; mom++){ |
76 | fNet[mom] = NULL; | |
77 | } | |
78 | } | |
79 | ||
80 | //__________________________________________ | |
81 | AliHFEpidQA::AliHFEpidQA(const AliHFEpidQA &ref): | |
82 | TObject(ref), | |
83 | fEvent(NULL), | |
84 | fMC(NULL), | |
85 | fV0pid(NULL), | |
86 | fV0pidMC(NULL), | |
87 | fTRDpidQA(NULL), | |
88 | fOutput(NULL), | |
89 | fESDpid(NULL), | |
e156c3bb | 90 | fNNref(NULL), |
e17c1f86 | 91 | fTRDTotalChargeInSlice0(ref.fTRDTotalChargeInSlice0), |
92 | fIsPbPb(kFALSE), | |
93 | fIsppMultiBin(kFALSE) | |
faee3b18 | 94 | { |
95 | // | |
96 | // Copy constructor | |
97 | // | |
98 | for(Int_t mom = 0; mom < 11; mom++){ | |
99 | fNet[mom] = NULL; | |
100 | } | |
101 | ref.Copy(*this); | |
102 | } | |
103 | ||
104 | //__________________________________________ | |
105 | AliHFEpidQA &AliHFEpidQA::operator=(const AliHFEpidQA &ref){ | |
106 | // | |
107 | // Assignment operator | |
108 | // | |
109 | if(this != &ref) | |
110 | ref.Copy(*this); | |
111 | return *this; | |
70da6c5a | 112 | } |
113 | ||
114 | //__________________________________________ | |
115 | AliHFEpidQA::~AliHFEpidQA(){ | |
116 | // | |
117 | // Destructor | |
118 | // | |
3e03bd67 | 119 | |
120 | // the pointers bellow are not dfeleted to prevend double-deleting of some of the content | |
121 | // these pointers are defined only once during the program call and should not cause a problem, | |
122 | // but cleaner solution is necessary. | |
123 | //if(fV0pid) delete fV0pid; | |
124 | //if(fV0pidMC) delete fV0pidMC; | |
125 | //if(fOutput) delete fOutput; | |
126 | ||
faee3b18 | 127 | // if(fTRDpidResponse) delete fTRDpidResponse; |
128 | } | |
129 | ||
130 | //__________________________________________ | |
131 | void AliHFEpidQA::Copy(TObject &o) const { | |
132 | // | |
133 | // Copy function | |
134 | // | |
135 | ||
136 | TObject::Copy(o); | |
137 | ||
138 | AliHFEpidQA &target = dynamic_cast<AliHFEpidQA &>(o); | |
139 | target.fMC = fMC; | |
e156c3bb | 140 | target.fTRDTotalChargeInSlice0 = fTRDTotalChargeInSlice0; |
faee3b18 | 141 | |
142 | if(target.fESDpid) delete target.fESDpid; | |
8c1c76e9 | 143 | target.fESDpid = fESDpid; |
faee3b18 | 144 | if(target.fV0pid) delete target.fV0pid; |
145 | if(fV0pid) | |
146 | target.fV0pid = dynamic_cast<AliHFEV0pid *>(fV0pid->Clone()); | |
147 | else | |
148 | target.fV0pid = NULL; | |
149 | if(target.fV0pidMC) delete target.fV0pidMC; | |
150 | if(fV0pidMC) | |
151 | target.fV0pidMC = dynamic_cast<AliHFEV0pidMC *>(fV0pidMC->Clone()); | |
152 | else | |
153 | target.fV0pidMC = NULL; | |
154 | if(target.fTRDpidQA) delete target.fTRDpidQA; | |
155 | if(fTRDpidQA) | |
156 | target.fTRDpidQA = dynamic_cast<AliHFEtrdPIDqa *>(fTRDpidQA->Clone()); | |
157 | else | |
158 | target.fTRDpidQA = NULL; | |
159 | if(target.fOutput) delete target.fOutput; | |
160 | if(fOutput) | |
161 | target.fOutput = dynamic_cast<AliHFEcollection *>(fOutput->Clone()); | |
70da6c5a | 162 | } |
163 | ||
164 | //__________________________________________ | |
165 | void AliHFEpidQA::Init(){ | |
166 | // | |
167 | // Prepare task output | |
168 | // | |
169 | ||
3a72645a | 170 | // load networks |
faee3b18 | 171 | if(fNNref){ |
3a72645a | 172 | for(Int_t mom = 0; mom < 11; mom++){ |
faee3b18 | 173 | fNet[mom] = (TMultiLayerPerceptron*) fNNref->Get(Form("NN_Mom%d", mom)); |
174 | if(!fNet[mom]){ | |
175 | AliError(Form("No reference network for momentum bin %d!", mom)); | |
176 | } | |
177 | } | |
178 | } | |
179 | ||
63bffcf1 | 180 | fV0pid = new AliHFEV0pid("fV0pid"); |
70da6c5a | 181 | if(HasV0pidQA()) fV0pid->InitQA(); |
182 | fV0pidMC = new AliHFEV0pidMC(); | |
183 | fV0pidMC->Init(); | |
184 | ||
faee3b18 | 185 | fTRDpidQA = new AliHFEtrdPIDqa; |
e156c3bb | 186 | if(fTRDTotalChargeInSlice0) fTRDpidQA->SetTotalChargeInSlice0(); |
faee3b18 | 187 | fTRDpidQA->Init(); |
188 | ||
70da6c5a | 189 | fOutput = new AliHFEcollection("pidQA", "PID QA output"); |
190 | ||
faee3b18 | 191 | const char *name[4] = {"Electron", "PionK0", "PionL", "Proton"}; |
192 | const char *title[4] = {"Electron", "K0 Pion", "Lambda Pion", "Proton"}; | |
193 | const char *det[4] = {"ITS", "TPC", "TRD", "TOF"}; | |
8c1c76e9 | 194 | const int effs[6] = {70, 75, 80, 85, 90, 95}; |
faee3b18 | 195 | for(Int_t i = 0; i < 4; i++){ |
196 | fOutput->CreateTH2F(Form("purity%s", name[i]), Form("%s Purity", title[i]), 2, -0.5, 1.5, 20, 0.1, 10, 1); | |
197 | ||
198 | for(Int_t idet = 0; idet < 4; idet++){ | |
199 | // create all the histograms which all the detectors have in common | |
200 | if(idet != 2){ // no nSigma histogram for TRD | |
201 | fOutput->CreateTH2F(Form("h%s_nSigma_%s", det[idet], name[i]), Form("%s number of sigmas for %ss; p (GeV/c); number of sigmas", det[idet], title[i]), 20, 0.1, 10, 100, -7, 7, 0); | |
202 | } | |
203 | fOutput->CreateTH2F(Form("h%s_PID_p_%s", det[idet], name[i]), Form("%s PID for %ss; p (GeV/c); ITS PID", det[idet], title[i]), 100, 0.1, 10, 5, -0.5, 4.5, 0); | |
204 | fOutput->CreateTH2F(Form("h%s_El_like_%s", det[idet], name[i]), Form("%s Electron Likelihoods for %ss; p (GeV/c); likelihood", det[idet], title[i]), 25, 0.1, 10, 1000, 0., 1., 0); | |
205 | fOutput->CreateTH2F(Form("h%s_El_like_MC_%s", det[idet], name[i]), Form("%s Electron Likelihoods for MC %ss; p (GeV/c); likelihood", det[idet], title[i]), 25, 0.1, 10, 1000, 0., 1., 0); | |
206 | } | |
70da6c5a | 207 | |
faee3b18 | 208 | // |
209 | // ITS pid response | |
210 | // | |
211 | fOutput->CreateTH2F(Form("hITS_Signal_%s", name[i]), Form("ITS Signal org. for %ss", title[i]), 40, 0.1, 10, 400, 0, 1000, 0); | |
212 | fOutput->CreateTH2Fvector1(5, Form("hITS_dEdx_%s", name[i]), Form("ITS dEdx for %ss; p (GeV/c); dEdx (a.u.)", title[i]), 40, 0.1, 10, 400, 0, 1000, 0); | |
213 | ||
214 | // | |
215 | // TPC pid response | |
216 | // | |
217 | fOutput->CreateTH2F(Form("hTPC_dEdx_%s", name[i]), Form("TPC dEdx for %ss; p (GeV/c); dEdx (a.u.)", title[i]), 20, 0.1, 10, 200, 0, 200, 0); | |
218 | ||
219 | // | |
220 | // TRD pid response | |
221 | // | |
222 | fOutput->CreateTH2F(Form("hTRD_trk_%s", name[i]), Form("%s PID tracklets; p (GeV/c); N TRD tracklets", title[i]), 100, 0.1, 10, 7, -0.5, 6.5, 0); | |
223 | // number of the non 0 slices per tracklet | |
224 | fOutput->CreateTH2F(Form("hTRD_Nslc_%s", name[i]), Form("%s PID slices > 0; p (GeV/c); N slc > 0", title[i]), 100, 0.1, 10, 9, -0.5, 8.5, 0); | |
225 | // location of the slices > 0 - where are the emtpy slices located ? | |
226 | fOutput->CreateTH2F(Form("hTRD_slc_%s", name[i]), Form("%s PID slices > 0 position; p (GeV/c); slice", title[i]), 100, 0.1, 10, 9, -0.5, 8.5, 0); | |
227 | fOutput->CreateTH2F(Form("hTRD_cls_%s", name[i]), Form("TRD clusters for %s candidates; p (GeV/c); N cls", title[i]), 25, 0.1, 10, 1000, 0, 1000); | |
228 | fOutput->CreateTH2F(Form("hTRD_dEdx_%s", name[i]), Form("TRD dEdx (trk) for %s candidates; p (GeV/c); tracklet dEdx (a.u.)", title[i]), 25, 0.1, 10, 1000, 0, 100000, 0); | |
8c1c76e9 | 229 | |
230 | for(Int_t itl = 4; itl < 6; itl++){ | |
231 | for(Int_t ieff = 0; ieff < 6; ieff++){ | |
232 | fOutput->CreateTH2F(Form("hTRD_likeSel_%s_%dtls_%deff", name[i], itl, effs[ieff]), Form(" %s selected as electrons with %d tracklets and %f electron efficiency", title[i], itl, static_cast<Double_t>(effs[ieff])/100.), 44, 0.1, 20., 100, 0., 1.); | |
233 | fOutput->CreateTH2F(Form("hTRD_likeRej_%s_%dtls_%deff", name[i], itl, effs[ieff]), Form(" %s rejected as electrons with %d tracklets and %f electron efficiency", title[i], itl, static_cast<Double_t>(effs[ieff])/100.), 44, 0.1, 20., 100, 0., 1.); | |
234 | } | |
235 | } | |
faee3b18 | 236 | // |
237 | // TOF pid response | |
238 | // | |
239 | fOutput->CreateTH2F(Form("hTOF_beta_%s", name[i]), Form("TOF beta for %s; p (GeV/c); #beta", title[i]), 50, 0.1, 5, 350, 0.4, 1.1, 0); | |
240 | }//.. loop over identified particle species | |
241 | ||
242 | // Global histograms | |
243 | fOutput->CreateTH1F("hITS_dEdx_nSamples", "ITS - number of non 0 dEdx samples", 4, 0.5, 4.5); | |
70da6c5a | 244 | fOutput->CreateTH2F("hTPC_PID", "TPC pid all tracks; tpc pid probability; species",100, 0, 1, 5, -0.5, 4.5 ); |
70da6c5a | 245 | fOutput->CreateTH2F("hTOF_PID", "TOF pid all tracks; tof pid probability; species",100, 0, 1,5, -0.5, 4.5 ); |
faee3b18 | 246 | fOutput->CreateTH2F("hTOF_beta_all", "TOF beta for all nice single tracks; p (GeV/c); #beta", 100, 0.1, 10, 480, 0, 1.2, 0); |
247 | fOutput->CreateTH2F("hTOF_qa_TmT0mT", "TOF (t - t0 - t[pion]) qa verus run", 10000, 114000, 124000, 200, -200, 200); | |
248 | fOutput->CreateTH2F("hTOF_qa_length", "TOF track length verus run", 10000, 114000, 112400, 200, 0, 1000); | |
249 | ||
250 | // | |
251 | // debug histograms | |
252 | // | |
253 | fOutput->CreateTH2F("hITS_kFlags", "ITS flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5); | |
254 | fOutput->CreateTH2F("hTPC_kFlags", "TPC flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5); | |
255 | fOutput->CreateTH2F("hTRD_kFlags", "TRD flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5); | |
256 | fOutput->CreateTH2F("hTOF_kFlags", "TOF flags; flag; V0 candidates", 5, 0.5, 5.5, 5, -0.5, 4.5); | |
70da6c5a | 257 | |
258 | ||
faee3b18 | 259 | // |
260 | // event based histograms | |
261 | // | |
e3fc062d | 262 | Int_t nT0[2] = {10000, 100}; |
263 | Double_t minT0[2] = {114500, -1000}; | |
264 | Double_t maxT0[2] = {124500, 1000}; | |
265 | fOutput->CreateTHnSparse("hEvent_T0", "T0 as a function of run number; run number; T0 (ns)", 2, nT0, minT0, maxT0); | |
faee3b18 | 266 | |
267 | // | |
268 | // test the tender V0 supply | |
269 | // | |
270 | fOutput->CreateTH2F("h_tender_check_01", "tender -vs- HFEpidQA; HFEpidQA V0 candiadates; tender V0 candidates", 4, -1.5, 2.5, 4, -1.5, 2.5); | |
271 | fOutput->CreateTH2Fvector1(3, "h_tender_check_02", "tender -vs- HFEpidQA per type; AliHFEpidQA V0 ; tender V0", 4, -1.5, 2.5, 4, -1.5, 2.5); | |
70da6c5a | 272 | |
70da6c5a | 273 | |
70da6c5a | 274 | // |
faee3b18 | 275 | // THnSpasre objects |
70da6c5a | 276 | // |
277 | ||
faee3b18 | 278 | // Create Illumination Plot |
279 | // bins: species, pt, eta, phi, TPC status, TRD status | |
280 | { | |
281 | Int_t nbins[6] = {4, 40, 16, 72, 2, 2}; | |
282 | Double_t min[6] = { 0, 0.1, -0.8, 0., 0., 0.}; | |
283 | Double_t max[6] = { 4., 20., 0.8, 2*TMath::Pi(), 2., 2.}; | |
284 | fOutput->CreateTHnSparse("hIllumination", "Illumination", 6, nbins, min, max); | |
285 | fOutput->BinLogAxis("hIllumination", 1); | |
286 | } | |
287 | ||
288 | // TPC THnSparse | |
289 | // bins: species, pt, n TPC clusters., TPC electron likelihood, TPC n sigmas, TPC signal | |
290 | { | |
291 | Int_t nbins[6] = { 5, 40, 20, 100, 100, 200}; | |
292 | Double_t min[6] = { -0.5, 0.1, 0., 0., -5., 0.}; | |
293 | Double_t max[6] = { 4.5, 20., 200, 1., 5., 120.}; | |
294 | TString htitle = "TPC info sparse; VO identified species; p (GeV/c); n TPC clusters; TPC N sigma; TPC signal"; | |
295 | fOutput->CreateTHnSparse("hTPCclusters",htitle,6, nbins, min, max); | |
296 | fOutput->BinLogAxis("hTPCclusters", 1); | |
297 | } | |
298 | // TRD THnSparse - entries per tracklet | |
299 | // species, p, tracklet position, number of PID tracklets, number of slices (non 0), number of clusters, electron likelihood, | |
3a72645a | 300 | |
faee3b18 | 301 | { |
302 | Int_t nbins[7] = {5, 20, 6, 7, 9, 45, 100}; | |
303 | Double_t min[7] = {-0.5, 0.5, -0.5, -0.5, -0.5, -0.5, 0.}; | |
304 | Double_t max[7] = {4.5, 10, 5.5, 6.5, 8.5, 179.5, 1.}; | |
305 | TString htitle = "TRD tracklets; VO identified species; p (GeV/c); tracklet position; No. PID tacklets; No. slices; No. clusters; El. likelihood"; | |
306 | fOutput->CreateTHnSparse("hTRDtracklets",htitle,7, nbins, min, max); | |
307 | fOutput->BinLogAxis("hTRDtracklets", 1); | |
70da6c5a | 308 | } |
faee3b18 | 309 | |
310 | ||
311 | } | |
312 | //__________________________________________ | |
313 | void AliHFEpidQA::Process(){ | |
314 | // | |
315 | // Run PID QA | |
316 | // | |
70da6c5a | 317 | |
318 | if(!fV0pid){ | |
319 | AliError("V0pid not available! Forgotten to initialize?"); | |
320 | return; | |
321 | } | |
faee3b18 | 322 | if(!fESDpid){ |
323 | AliError("fESDpid not initialized, I am leaving this!"); | |
324 | return; | |
325 | } | |
326 | ||
327 | // to be udpated to AOD save mdoe | |
328 | if(!fEvent){ | |
329 | AliError("AliVEvent not available, returning"); | |
330 | } | |
70da6c5a | 331 | |
3a72645a | 332 | |
70da6c5a | 333 | if(fMC) fV0pidMC->SetMCEvent(fMC); |
3a72645a | 334 | if(fMC) fV0pid->SetMCEvent(fMC); |
70da6c5a | 335 | |
faee3b18 | 336 | fV0pid->Process(fEvent); |
337 | TObjArray *hfeelectrons = fV0pid->GetListOfElectrons(); | |
338 | TObjArray *hfepionsK0 = fV0pid->GetListOfPionsK0(); | |
339 | TObjArray *hfepionsL = fV0pid->GetListOfPionsL(); | |
340 | TObjArray *hfeprotons = fV0pid->GetListOfProtons(); | |
341 | ||
342 | // Get Track list for normal purpose | |
343 | TObjArray *electrons = MakeTrackList(hfeelectrons); | |
344 | TObjArray *pionsK0 = MakeTrackList(hfepionsK0); | |
345 | TObjArray *pionsL = MakeTrackList(hfepionsL); | |
346 | TObjArray *protons = MakeTrackList(hfeprotons); | |
347 | TObjArray *cleanElectrons = MakeCleanListElectrons(hfeelectrons); | |
70da6c5a | 348 | |
349 | if(fMC){ | |
3a72645a | 350 | fV0pidMC->Process(electrons, AliHFEV0cuts::kRecoElectron); |
351 | fV0pidMC->Process(pionsK0, AliHFEV0cuts::kRecoPionK0); | |
352 | fV0pidMC->Process(pionsL, AliHFEV0cuts::kRecoPionL); | |
353 | fV0pidMC->Process(protons, AliHFEV0cuts::kRecoProton); | |
70da6c5a | 354 | } |
355 | ||
356 | AliDebug(2, Form("Number of Electrons : %d", electrons->GetEntries())); | |
357 | AliDebug(2, Form("Number of K0 Pions : %d", pionsK0->GetEntries())); | |
358 | AliDebug(2, Form("Number of Lambda Pions : %d", pionsL->GetEntries())); | |
359 | AliDebug(2, Form("Number of Protons : %d", protons->GetEntries())); | |
360 | if(fMC){ | |
361 | AliDebug(2, "MC Information available. Doing Purity checks..."); | |
362 | // Calculate the purity of the clean samples using MC | |
3a72645a | 363 | MakePurity(electrons, AliHFEV0cuts::kRecoElectron); |
364 | MakePurity(pionsK0, AliHFEV0cuts::kRecoPionK0); | |
365 | MakePurity(pionsL, AliHFEV0cuts::kRecoPionL); | |
366 | MakePurity(protons, AliHFEV0cuts::kRecoProton); | |
70da6c5a | 367 | } |
faee3b18 | 368 | |
369 | // some event wise checks | |
370 | CheckEvent(); | |
371 | ||
372 | // Make Illumination Plot | |
3a72645a | 373 | FillIllumination(electrons, AliHFEV0cuts::kRecoElectron); |
374 | FillIllumination(pionsK0, AliHFEV0cuts::kRecoPionK0); | |
375 | FillIllumination(pionsL, AliHFEV0cuts::kRecoPionL); | |
376 | FillIllumination(protons, AliHFEV0cuts::kRecoProton); | |
faee3b18 | 377 | |
70da6c5a | 378 | // Now we can do studies on the PID itself |
faee3b18 | 379 | // For TRD use the TRD PID QA object |
e17c1f86 | 380 | Int_t centrality=-1; |
381 | if(fIsPbPb) centrality=GetCentrality(fEvent); | |
382 | if(fIsppMultiBin) centrality=GetMultiplicityITS(fEvent); | |
383 | fTRDpidQA->SetCentralityBin(centrality); | |
faee3b18 | 384 | fTRDpidQA->ProcessTracks(cleanElectrons, AliPID::kElectron); |
385 | fTRDpidQA->ProcessTracks(pionsK0, AliPID::kPion); | |
386 | fTRDpidQA->ProcessTracks(pionsL, AliPID::kPion); | |
387 | fTRDpidQA->ProcessTracks(protons, AliPID::kProton); | |
388 | ||
8c1c76e9 | 389 | // Monitor TRD PID Response |
390 | /*TestTRDResponse(cleanElectrons, AliHFEV0cuts::kRecoElectron); | |
391 | TestTRDResponse(pionsK0, AliHFEV0cuts::kRecoPionK0); | |
392 | TestTRDResponse(pionsL, AliHFEV0cuts::kRecoPionL); | |
393 | TestTRDResponse(protons, AliHFEV0cuts::kRecoProton);*/ | |
394 | ||
3a72645a | 395 | FillElectronLikelihoods(electrons, AliHFEV0cuts::kRecoElectron); |
396 | FillElectronLikelihoods(pionsK0, AliHFEV0cuts::kRecoPionK0); | |
397 | FillElectronLikelihoods(pionsL, AliHFEV0cuts::kRecoPionL); | |
398 | FillElectronLikelihoods(protons, AliHFEV0cuts::kRecoProton); | |
70da6c5a | 399 | |
3a72645a | 400 | FillPIDresponse(electrons, AliHFEV0cuts::kRecoElectron); |
401 | FillPIDresponse(pionsK0, AliHFEV0cuts::kRecoPionK0); | |
402 | FillPIDresponse(pionsL, AliHFEV0cuts::kRecoPionL); | |
403 | FillPIDresponse(protons, AliHFEV0cuts::kRecoProton); | |
70da6c5a | 404 | |
faee3b18 | 405 | // check the tender V0s |
3a72645a | 406 | CheckTenderV0pid(electrons, AliHFEV0cuts::kRecoElectron); |
407 | CheckTenderV0pid(pionsK0, AliHFEV0cuts::kRecoPionK0); | |
408 | CheckTenderV0pid(pionsL, AliHFEV0cuts::kRecoPionL); | |
409 | CheckTenderV0pid(protons, AliHFEV0cuts::kRecoProton); | |
faee3b18 | 410 | |
70da6c5a | 411 | // Analysis done, flush the containers |
412 | fV0pid->Flush(); | |
faee3b18 | 413 | |
414 | delete electrons; | |
415 | delete pionsL; | |
416 | delete pionsK0; | |
417 | delete protons; | |
418 | delete cleanElectrons; | |
419 | } | |
420 | ||
421 | //__________________________________________ | |
3a72645a | 422 | void AliHFEpidQA::FillIllumination(const TObjArray * const tracks, Int_t species){ |
faee3b18 | 423 | // |
424 | // Fill Illumination Plot | |
425 | // | |
426 | THnSparseF *hIllumination = dynamic_cast<THnSparseF *>(fOutput->Get("hIllumination")); | |
427 | if(!hIllumination) return; | |
428 | ||
429 | Double_t quantities[6]; memset(quantities, 0, sizeof(Double_t) *6); | |
430 | TIter trackIter(tracks); | |
431 | ||
432 | quantities[0] = species; | |
433 | TObject *o = NULL; AliESDtrack *esdtrack = NULL; | |
434 | while((o = trackIter())){ | |
435 | if(!TString(o->IsA()->GetName()).CompareTo("AliESDtrack")){ | |
436 | // work on local copy in order to not spoil others | |
e3ae862b | 437 | esdtrack = new AliESDtrack(*(static_cast<AliESDtrack *>(o))); |
bf892a6a | 438 | if(!esdtrack) continue; |
faee3b18 | 439 | } else if(!TString(o->IsA()->GetName()).CompareTo("AliAODrack")){ |
440 | // Bad hack: Fill ESD track with AOD information | |
215ffe88 | 441 | esdtrack = new AliESDtrack(static_cast<AliAODTrack *>(o)); |
bf892a6a | 442 | if(!esdtrack) continue; |
faee3b18 | 443 | } else { |
444 | // Non usable | |
445 | continue; | |
446 | } | |
447 | ||
448 | // Propagate to the entrance of the TRD | |
449 | esdtrack->PropagateTo(300, fEvent->GetMagneticField()); | |
450 | quantities[1] = esdtrack->Pt(); | |
451 | quantities[2] = esdtrack->Eta(); | |
452 | quantities[3] = esdtrack->Phi(); | |
453 | quantities[4] = (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) ? 1 : 0; | |
454 | quantities[5] = (esdtrack->GetStatus() & AliESDtrack::kTRDout) ? 1. : 0.; | |
455 | hIllumination->Fill(quantities); | |
456 | ||
457 | delete esdtrack; | |
458 | } | |
459 | } | |
460 | //__________________________________________ | |
461 | void AliHFEpidQA::FillTPCinfo(AliESDtrack *const esdtrack, Int_t species){ | |
462 | // | |
463 | // Fill TPC Cluster Plots | |
464 | // | |
465 | THnSparseF *hTPCclusters = dynamic_cast<THnSparseF *>(fOutput->Get("hTPCclusters")); | |
466 | if(!hTPCclusters) return; | |
467 | ||
468 | Double_t quantities[6]; memset(quantities, 0, sizeof(Double_t) *6); | |
469 | ||
470 | Double_t pidProbs[5]; | |
471 | const Int_t typePID[5] = {0, 2, 2, 3, 4}; | |
472 | ||
473 | ||
474 | quantities[0] = species; | |
475 | ||
476 | ||
477 | esdtrack->GetTPCpid(pidProbs); | |
478 | ||
479 | quantities[1] = esdtrack->P(); | |
480 | quantities[2] = esdtrack->GetTPCNcls(); | |
481 | quantities[3] = pidProbs[0]; | |
482 | quantities[4] = fESDpid->NumberOfSigmasTPC(esdtrack,(AliPID::EParticleType)typePID[species]); | |
483 | quantities[5] = esdtrack->GetTPCsignal(); | |
484 | hTPCclusters->Fill(quantities); | |
485 | ||
70da6c5a | 486 | } |
487 | ||
8c1c76e9 | 488 | //__________________________________________ |
489 | void AliHFEpidQA::TestTRDResponse(const TObjArray * const tracks, Int_t species){ | |
490 | // | |
491 | // Test PID Response function of the TRD | |
492 | // | |
493 | Int_t effInt[6] = {70, 75, 80, 85, 90, 95}; | |
494 | const char *sname[5] = {"Electron", "PionK0", "PionL", "Kaon", "Proton"}; | |
495 | AliESDtrack *track = NULL; | |
496 | TIter trackIter(tracks); | |
497 | Float_t momenta[6], p; | |
498 | Int_t nmomenta; | |
499 | Double_t probs[5], likeEl; | |
500 | while((track = static_cast<AliESDtrack *>(trackIter()))){ | |
501 | if(track->GetTRDntrackletsPID() < 4) continue; | |
502 | ||
503 | // calculate momentum at TRD level | |
504 | memset(momenta, 0, sizeof(Float_t) * 6); nmomenta = 0; | |
505 | for(Int_t ily = 0; ily < 6; ily++){ | |
506 | if(track->GetTRDmomentum(ily) > 0.01) momenta[nmomenta++] = track->GetTRDmomentum(ily); | |
507 | } | |
508 | p = TMath::Mean(nmomenta, momenta); | |
509 | ||
510 | // Get Electron likelihood | |
511 | track->GetTRDpid(probs); | |
512 | likeEl = probs[0]/(probs[0] + probs[2]); | |
513 | ||
514 | for(Int_t ieff = 0; ieff < 6; ieff++){ | |
515 | if(fESDpid->IdentifiedAsElectronTRD(track, static_cast<Double_t>(effInt[ieff])/100.)) fOutput->Fill(Form("hTRD_likeSel_%s_%dtls_%deff", sname[species], static_cast<Int_t>(track->GetTRDntrackletsPID()), effInt[ieff]), p, likeEl); | |
516 | else fOutput->Fill(Form("hTRD_likeRej_%s_%dtls_%deff", sname[species], static_cast<Int_t>(track->GetTRDntrackletsPID()), effInt[ieff]), p, likeEl); | |
517 | } | |
518 | } | |
519 | } | |
520 | ||
70da6c5a | 521 | //__________________________________________ |
3a72645a | 522 | void AliHFEpidQA::MakePurity(const TObjArray *tracks, Int_t species){ |
70da6c5a | 523 | // |
524 | // Fill the QA histos for a given species | |
525 | // | |
526 | if(!fMC) return; | |
527 | AliDebug(3, Form("Doing Purity checks for species %d", species)); | |
faee3b18 | 528 | Int_t pdg = GetPDG(species); |
e3ae862b | 529 | TString hname; |
faee3b18 | 530 | |
70da6c5a | 531 | switch(species){ |
3a72645a | 532 | case AliHFEV0cuts::kRecoElectron: |
e3ae862b | 533 | hname = "purityElectron"; |
faee3b18 | 534 | break; |
3a72645a | 535 | case AliHFEV0cuts::kRecoPionK0: |
e3ae862b | 536 | hname = "purityPionK0"; |
faee3b18 | 537 | break; |
3a72645a | 538 | case AliHFEV0cuts::kRecoPionL: |
e3ae862b | 539 | hname = "purityPionL"; |
faee3b18 | 540 | break; |
3a72645a | 541 | case AliHFEV0cuts::kRecoProton: |
e3ae862b | 542 | hname = "purityProton"; |
faee3b18 | 543 | break; |
544 | default: // non investigated species | |
545 | AliDebug(3, "Species not investigated"); | |
546 | return; | |
547 | } | |
548 | ||
70da6c5a | 549 | AliDebug(3, Form("Number of tracks: %d", tracks->GetEntries())); |
faee3b18 | 550 | TIter trackIter(tracks); |
70da6c5a | 551 | AliVParticle *recTrack = NULL, *mcTrack = NULL; |
faee3b18 | 552 | while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){ |
70da6c5a | 553 | Int_t label = recTrack->GetLabel(); |
554 | AliDebug(4, Form("MC Label %d", label)); | |
555 | mcTrack =fMC->GetTrack(TMath::Abs(label)); | |
556 | if(!mcTrack){ | |
557 | AliDebug(4, "MC track not available"); | |
558 | continue; // we don't know | |
559 | } | |
560 | ||
561 | // Get the pdg code | |
562 | Int_t trackPdg = 0; | |
563 | if(!TString(mcTrack->IsA()->GetName()).CompareTo("AliMCParticle")){ | |
564 | // case ESD | |
565 | AliMCParticle *mcp = dynamic_cast<AliMCParticle *>(mcTrack); | |
bf892a6a | 566 | if(!mcp) continue; |
70da6c5a | 567 | trackPdg = TMath::Abs(mcp->Particle()->GetPdgCode()); |
568 | } else { | |
569 | // case AOD | |
570 | AliAODMCParticle *aodmcp = dynamic_cast<AliAODMCParticle *>(mcTrack); | |
bf892a6a | 571 | if(!aodmcp) continue; |
70da6c5a | 572 | trackPdg = TMath::Abs(aodmcp->GetPdgCode()); |
573 | } | |
574 | if(trackPdg == pdg) // Correct identification | |
faee3b18 | 575 | { |
576 | fOutput->Fill(hname, 0., recTrack->Pt()); | |
577 | } | |
70da6c5a | 578 | else // Wrong identification |
579 | fOutput->Fill(hname, 1., recTrack->Pt()); | |
580 | } | |
70da6c5a | 581 | } |
582 | ||
583 | //__________________________________________ | |
3a72645a | 584 | void AliHFEpidQA::FillElectronLikelihoods(const TObjArray * const particles, Int_t species){ |
70da6c5a | 585 | // |
faee3b18 | 586 | // Fill electron Likelihoods for the ITS, TPC and TOF |
70da6c5a | 587 | // Required for the calculation of the electron efficiency, |
588 | // pion and proton efficiency and the thresholds | |
589 | // | |
faee3b18 | 590 | Long_t status = 0; |
e3ae862b | 591 | const TString detname[4] = {"ITS", "TPC", "TRD", "TOF"}; |
592 | TString specname; | |
faee3b18 | 593 | |
70da6c5a | 594 | switch(species){ |
3a72645a | 595 | case AliHFEV0cuts::kRecoElectron: |
e3ae862b | 596 | specname = "Electron"; |
faee3b18 | 597 | break; |
3a72645a | 598 | case AliHFEV0cuts::kRecoPionK0: |
e3ae862b | 599 | specname = "PionK0"; |
faee3b18 | 600 | break; |
3a72645a | 601 | case AliHFEV0cuts::kRecoPionL: |
e3ae862b | 602 | specname = "PionL"; |
faee3b18 | 603 | break; |
3a72645a | 604 | case AliHFEV0cuts::kRecoProton: |
e3ae862b | 605 | specname = "Proton"; |
faee3b18 | 606 | break; |
607 | default: | |
608 | AliDebug(2, Form("Species %d not investigated", species)); | |
609 | return; | |
70da6c5a | 610 | }; |
611 | AliVParticle *recTrack = NULL; | |
faee3b18 | 612 | // mcTrack =fMC->GetTrack(TMath::Abs(label)); |
613 | // if(!mcTrack){ | |
614 | // AliDebug(4, "MC track not available"); | |
615 | // continue; // we don't know | |
616 | // } | |
617 | ||
618 | TIter trackIter(particles); | |
619 | ||
620 | Double_t quantities[2]; | |
621 | Double_t pidProbs[5]; | |
622 | ||
623 | while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){ | |
70da6c5a | 624 | if(!TString(recTrack->IsA()->GetName()).CompareTo("AliESDtrack")){ |
625 | // case ESD | |
626 | AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(recTrack); | |
bf892a6a | 627 | if(!esdTrack) continue; |
faee3b18 | 628 | status = esdTrack->GetStatus(); |
629 | ||
630 | //TPC momentum and likelihoods | |
631 | Double_t pTPC = 0.; | |
632 | pTPC = esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P(); | |
633 | Bool_t mcFound = kFALSE; | |
634 | if(fMC){ | |
635 | Int_t label = esdTrack->GetLabel(); | |
636 | AliMCParticle *mcTrack = dynamic_cast<AliMCParticle *>(fMC->GetTrack(label)); | |
637 | Int_t pdg = GetPDG(species); | |
638 | Int_t trackPdg = 0; | |
639 | if(mcTrack){ | |
640 | trackPdg = TMath::Abs(mcTrack->Particle()->GetPdgCode()); | |
641 | } | |
642 | if(pdg == trackPdg) mcFound = kTRUE; | |
70da6c5a | 643 | } |
faee3b18 | 644 | quantities[0] = pTPC; |
645 | Bool_t detFlagSet = kFALSE; | |
646 | for(Int_t idet = 0; idet < 4; idet++){ | |
e3ae862b | 647 | TString histname, histnameMC; |
e97c2edf | 648 | histname = "h" + detname[idet] + "_El_like_" + specname; |
649 | histnameMC = "h" + detname[idet] + "_El_like_MC_" + specname; | |
e3ae862b | 650 | |
faee3b18 | 651 | switch(idet){ |
652 | case kITS: esdTrack->GetITSpid(pidProbs); | |
653 | detFlagSet = status & AliESDtrack::kITSpid; | |
654 | break; | |
655 | case kTPC: esdTrack->GetTPCpid(pidProbs); | |
656 | detFlagSet = status & AliESDtrack::kTPCpid; | |
657 | break; | |
658 | case kTRD: esdTrack->GetTRDpid(pidProbs); | |
659 | detFlagSet = status & AliESDtrack::kTRDpid; | |
660 | break; | |
661 | case kTOF: esdTrack->GetTOFpid(pidProbs); | |
662 | detFlagSet = status & AliESDtrack::kTOFpid; | |
663 | break; | |
664 | }; | |
665 | quantities[1] = pidProbs[AliPID::kElectron]; | |
666 | // in case of TRD require 6 PID tracklets | |
667 | if(kTRD == idet && esdTrack->GetTRDntrackletsPID() != 6) continue; | |
668 | if(detFlagSet){ | |
669 | fOutput->Fill(histname, quantities[0], quantities[1]); | |
670 | if(mcFound) | |
671 | fOutput->Fill(histnameMC, quantities[0], quantities[1]); | |
672 | } | |
70da6c5a | 673 | } |
faee3b18 | 674 | }//.. ESD |
675 | else{ | |
676 | //AOD | |
677 | }//.. aod | |
678 | }//.. while tracks | |
70da6c5a | 679 | } |
680 | //__________________________________________ | |
3a72645a | 681 | void AliHFEpidQA::FillPIDresponse(const TObjArray * const particles, Int_t species){ |
70da6c5a | 682 | // |
683 | // Fill the PID response of different detectors to V0 daughter particles | |
684 | // | |
e3ae862b | 685 | TString hname, hname2, hname3; |
faee3b18 | 686 | |
e3ae862b | 687 | const TString typeName[5] = {"Electron", "PionK0", "PionL", "Kaon", "Proton"}; |
70da6c5a | 688 | const Int_t typePID[5] = {0, 2, 2, 3, 4}; |
faee3b18 | 689 | |
690 | // PID THnSparse | |
691 | // axes: | |
692 | // 0) species, 1) momentum, 2) DCA xy, 3) DCA z | |
693 | // 4) ITS signal | |
694 | // 5) TPC Ncls 6) TPC signal 7) TPC nSigma, | |
695 | // 8) TRD Ntrk, 9) TRD Ncls, 10) TRD dEdx, | |
696 | ||
8c1c76e9 | 697 | //Double_t data[12]; |
bf892a6a | 698 | |
faee3b18 | 699 | |
700 | Int_t run = fEvent->GetRunNumber(); | |
701 | ||
70da6c5a | 702 | AliVParticle *recTrack = NULL; |
faee3b18 | 703 | TIter trackIter(particles); |
704 | while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){ | |
8c1c76e9 | 705 | //for(Int_t i=0; i<12; ++i) data[i] = -99.; |
70da6c5a | 706 | // ESD |
707 | if(!TString(recTrack->IsA()->GetName()).CompareTo("AliESDtrack")){ | |
708 | // case ESD | |
709 | AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(recTrack); | |
faee3b18 | 710 | if(!esdTrack) continue; |
70da6c5a | 711 | |
faee3b18 | 712 | // for the PID THnSparse |
8c1c76e9 | 713 | //data[0] = species; |
714 | //data[1] = esdTrack->P(); | |
faee3b18 | 715 | Float_t impactR = -1.; |
716 | Float_t impactZ = -1.; | |
717 | esdTrack->GetImpactParameters(impactR, impactZ); | |
8c1c76e9 | 718 | //data[2] = impactR; |
719 | //data[3] = impactZ; | |
720 | //data[11] = 0; // initialize the TOF pid cut on elecgrons to false | |
faee3b18 | 721 | // use ONLY tracks with PID flag TRUE |
722 | ULong_t status = 0; | |
723 | status = esdTrack->GetStatus(); | |
724 | ||
725 | // | |
726 | // DEBUG | |
727 | // | |
728 | ||
729 | fOutput->Fill("hITS_kFlags", 5., species); | |
730 | if(status & AliESDtrack::kITSin) fOutput->Fill("hITS_kFlags", 1., species); | |
731 | if(status & AliESDtrack::kITSout) fOutput->Fill("hITS_kFlags", 2., species); | |
732 | if(status & AliESDtrack::kITSrefit) fOutput->Fill("hITS_kFlags", 3., species); | |
733 | if(status & AliESDtrack::kITSpid) fOutput->Fill("hITS_kFlags", 4., species); | |
734 | ||
735 | fOutput->Fill("hTPC_kFlags", 5., species); | |
736 | if(status & AliESDtrack::kTPCin) fOutput->Fill("hTPC_kFlags", 1., species); | |
737 | if(status & AliESDtrack::kTPCout) fOutput->Fill("hTPC_kFlags", 2., species); | |
738 | if(status & AliESDtrack::kTPCrefit) fOutput->Fill("hTPC_kFlags", 3., species); | |
739 | if(status & AliESDtrack::kTPCpid) fOutput->Fill("hTPC_kFlags", 4., species); | |
740 | ||
741 | fOutput->Fill("hTRD_kFlags", 5., species); | |
742 | if(status & AliESDtrack::kTRDin) fOutput->Fill("hTRD_kFlags", 1., species); | |
743 | if(status & AliESDtrack::kTRDout) fOutput->Fill("hTRD_kFlags", 2., species); | |
744 | if(status & AliESDtrack::kTRDrefit) fOutput->Fill("hTRD_kFlags", 3., species); | |
745 | if(status & AliESDtrack::kTRDpid) fOutput->Fill("hTRD_kFlags", 4., species); | |
746 | ||
747 | fOutput->Fill("hTOF_kFlags", 5., species); | |
748 | if(status & AliESDtrack::kTOFin) fOutput->Fill("hTOF_kFlags", 1., species); | |
749 | if(status & AliESDtrack::kTOFout) fOutput->Fill("hTOF_kFlags", 2., species); | |
750 | if(status & AliESDtrack::kTOFrefit) fOutput->Fill("hTOF_kFlags", 3., species); | |
751 | if(status & AliESDtrack::kTOFpid) fOutput->Fill("hTOF_kFlags", 4., species); | |
752 | ||
753 | ||
754 | // | |
755 | // ITS - | |
756 | // | |
757 | if(status & AliESDtrack::kITSpid){ | |
758 | Double_t p = esdTrack->P(); | |
759 | ||
760 | // ITS signal | |
761 | //Double_t itsSignal = esdTrack->GetITSsignal(); | |
762 | ||
763 | // ITS dEdx | |
764 | Double_t dEdxSamples[4]; | |
765 | esdTrack->GetITSdEdxSamples(dEdxSamples); | |
766 | Int_t nSamples = 0; | |
767 | Double_t dEdxSum = 0.; | |
e3ae862b | 768 | hname = "hITS_dEdx_" + typeName[species]; |
faee3b18 | 769 | for(Int_t i=0; i<4; ++i){ |
770 | if(dEdxSamples[i] > 0){ | |
771 | nSamples++; | |
772 | fOutput->Fill(hname, i+1, p, dEdxSamples[i]); | |
773 | dEdxSum += dEdxSamples[i]; | |
774 | } | |
70da6c5a | 775 | } |
faee3b18 | 776 | if(4 == nSamples)fOutput->Fill(hname, 0, p, dEdxSum); |
777 | fOutput->Fill("hITS_dEdx_nSamples", nSamples); | |
778 | ||
779 | Double_t signal = esdTrack->GetITSsignal(); | |
e3ae862b | 780 | hname = "hITS_Signal_" + typeName[species]; |
faee3b18 | 781 | fOutput->Fill(hname, p, signal); |
8c1c76e9 | 782 | //data[4] = signal; |
faee3b18 | 783 | |
784 | // ITS number of signas | |
785 | Double_t nsigma = fESDpid->NumberOfSigmasITS(esdTrack,(AliPID::EParticleType)typePID[species]); | |
e3ae862b | 786 | hname = "hITS_nSigma_" + typeName[species]; |
faee3b18 | 787 | fOutput->Fill(hname, p, nsigma); |
faee3b18 | 788 | // ITS PID response |
789 | Double_t itsPID[5] = {-1, -1, -1, -1, -1}; | |
790 | esdTrack->GetITSpid(itsPID); | |
791 | Int_t ix = GetMaxPID(itsPID); | |
e97c2edf | 792 | hname = "hITS_PID_p_" + typeName[species]; |
91c7e1ec | 793 | fOutput->Fill(hname, p, ix); |
faee3b18 | 794 | }//.. kITSpid |
795 | ||
796 | // | |
797 | // TPC | |
798 | // | |
799 | if(status & AliESDtrack::kTPCpid){ | |
800 | // Make TPC clusters Plot | |
8c1c76e9 | 801 | //data[5] = esdTrack->GetTPCNcls(); |
faee3b18 | 802 | FillTPCinfo(esdTrack, species); |
70da6c5a | 803 | |
faee3b18 | 804 | Double_t p = esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P(); |
805 | // TPC dEdx | |
806 | Double_t dEdx = esdTrack->GetTPCsignal(); | |
e3ae862b | 807 | hname = "hTPC_dEdx_" + typeName[species]; |
faee3b18 | 808 | fOutput->Fill(hname, p, dEdx); |
8c1c76e9 | 809 | //data[6] = dEdx; |
faee3b18 | 810 | |
811 | //TPC number of sigmas | |
812 | Double_t nsigma = fESDpid->NumberOfSigmasTPC(esdTrack,(AliPID::EParticleType)typePID[species]); | |
e3ae862b | 813 | hname = "hTPC_nSigma_" + typeName[species]; |
faee3b18 | 814 | fOutput->Fill(hname, p, nsigma); |
8c1c76e9 | 815 | //data[7] = nsigma; |
faee3b18 | 816 | |
817 | // TPC PID response | |
e3ae862b | 818 | hname = "hTPC_PID_p_" + typeName[species]; |
faee3b18 | 819 | Double_t tpcPID[5] = {-1, -1, -1, -1, -1}; |
820 | esdTrack->GetTPCpid(tpcPID); | |
821 | Int_t ix = GetMaxPID(tpcPID); | |
822 | fOutput->Fill(hname, p, ix); | |
823 | }//.. kTPCpid | |
824 | ||
825 | // | |
826 | // TRD | |
827 | // | |
828 | ||
829 | if(status & AliESDtrack::kTRDpid){ | |
830 | Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P(); | |
831 | ||
832 | // TRD number of tracklets | |
833 | Int_t ntrk = esdTrack->GetTRDntrackletsPID(); | |
e3ae862b | 834 | hname = "hTRD_trk_" + typeName[species]; |
faee3b18 | 835 | fOutput->Fill(hname, p, ntrk); |
8c1c76e9 | 836 | //data[8] = ntrk; |
faee3b18 | 837 | |
838 | // TRD PID response | |
e3ae862b | 839 | hname = "hTRD_PID_p_" + typeName[species]; |
faee3b18 | 840 | Double_t trdPID[5] = {-1., -1., -1., -1., -1}; |
841 | esdTrack->GetTRDpid(trdPID); | |
842 | Int_t ix = GetMaxPID(trdPID); | |
91c7e1ec | 843 | fOutput->Fill(hname, p, ix); |
faee3b18 | 844 | // TRD n clusters |
845 | Int_t ncls = esdTrack->GetTRDncls(); | |
e3ae862b | 846 | hname = "hTRD_cls_" + typeName[species]; |
faee3b18 | 847 | fOutput->Fill(hname, p, ncls); |
8c1c76e9 | 848 | //data[9] = ncls; |
faee3b18 | 849 | |
850 | // TRD - per tracklet - dEdx per, likelihood | |
e3ae862b | 851 | hname = "hTRD_Nslc_" + typeName[species]; |
852 | hname2 = "hTRD_slc_" + typeName[species]; | |
853 | hname3 = "hTRD_dEdx_" + typeName[species]; | |
faee3b18 | 854 | Int_t nSlices = esdTrack->GetNumberOfTRDslices(); |
855 | Double_t sumTot = 0.; | |
856 | Int_t not0Tot = 0; | |
857 | for(Int_t l=0; l< 6; ++l){ | |
858 | Double_t trkData[7] = {-1.,-1, -1, -1, -1, -1, -1}; | |
859 | trkData[0] = species; | |
860 | trkData[1] = p; | |
861 | trkData[2] = l; | |
862 | trkData[3] = ntrk; | |
863 | trkData[5] = ncls; | |
864 | Double_t sum = 0.; | |
865 | Int_t not0 = 0; | |
866 | for(Int_t s=0; s<nSlices; ++s){ | |
867 | Double_t slice = esdTrack->GetTRDslice(l, s); | |
868 | sum += slice; | |
869 | if(slice > 0){ | |
870 | not0 += 1; | |
871 | fOutput->Fill(hname2, p, s); | |
872 | } | |
873 | }//..slices | |
874 | ||
875 | trkData[4] = not0; | |
876 | fOutput->Fill(hname, p, not0); | |
877 | fOutput->Fill(hname3, p, sum); | |
878 | if(sum > 0){ | |
879 | sumTot += sum; | |
880 | not0Tot += 1; | |
881 | } | |
882 | // lkelihoods per layer | |
883 | if(not0Tot > 0 && fNNref){ | |
884 | Double_t likelihoods[5] = {-1., -1., -1., -1., -1}; | |
885 | TRDlikeTracklet(l, esdTrack, likelihoods); | |
886 | trkData[6] = likelihoods[0]; | |
887 | //printf(" -D: species: %i, P; %f : %f, s: %f\n", species, p, likelihoods[0], s); | |
888 | } | |
889 | if(not0Tot) fOutput->Fill("hTRDtracklets", trkData); | |
890 | }//..layers | |
891 | // average dEx per number of tracklets | |
8c1c76e9 | 892 | //if(0 < not0Tot) |
893 | // data[10] = sumTot / not0Tot; | |
faee3b18 | 894 | }//.. kTRDpid |
70da6c5a | 895 | |
faee3b18 | 896 | |
897 | // | |
898 | // TOF | |
899 | // | |
900 | if(status & AliESDtrack::kTOFpid){ | |
901 | Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P(); | |
bf892a6a | 902 | Double_t t0 = fESDpid->GetTOFResponse().GetStartTime(esdTrack->P()); |
faee3b18 | 903 | |
904 | //TOF beta | |
e3ae862b | 905 | hname = "hTOF_beta_" + typeName[species]; |
faee3b18 | 906 | Float_t beta = TOFbeta(esdTrack); |
907 | fOutput->Fill(hname, p, beta); | |
908 | fOutput->Fill("hTOF_beta_all", p, beta); | |
909 | // TOF nSigma | |
8c1c76e9 | 910 | Double_t nsigma = fESDpid->NumberOfSigmasTOF(esdTrack,(AliPID::EParticleType)typePID[species]); |
e3ae862b | 911 | hname = "hTOF_nSigma_" + typeName[species]; |
faee3b18 | 912 | fOutput->Fill(hname, p, nsigma); |
8c1c76e9 | 913 | //if(beta > 0.97 && beta < 1.03){ |
914 | // data[11] = 1; | |
915 | //} | |
faee3b18 | 916 | |
917 | // TOF PID response | |
e3ae862b | 918 | hname = "hTOF_PID_p_" + typeName[species]; |
faee3b18 | 919 | Double_t tofPID[5] = {-1., -1., -1., -1., -1}; |
920 | esdTrack->GetTOFpid(tofPID); | |
921 | Int_t ix = GetMaxPID(tofPID); | |
922 | fOutput->Fill(hname, p, ix); | |
923 | ||
924 | // time of flight QA | |
925 | // - distribution of (time - t0 - pion_time) | |
926 | Double_t times[5]; | |
927 | esdTrack->GetIntegratedTimes(times); | |
928 | Double_t tItrackL = esdTrack->GetIntegratedLength(); | |
929 | Double_t tTOFsignal = esdTrack->GetTOFsignal(); | |
930 | Double_t dT = tTOFsignal - t0 - times[2]; | |
931 | fOutput->Fill("hTOF_qa_TmT0mT", run*1.0, dT); | |
932 | fOutput->Fill("hTOF_qa_length", run*1.0, tItrackL); | |
70da6c5a | 933 | |
faee3b18 | 934 | |
935 | }//.. kTOFpid | |
936 | // temporary - the PIDsparse needs rebuilding | |
937 | //fOutput->Fill("PIDsparse", data); | |
70da6c5a | 938 | } |
939 | // AOD - comming soon | |
940 | else{ | |
941 | continue; | |
942 | } | |
943 | }// .. tracks in TObjArray | |
944 | ||
70da6c5a | 945 | |
946 | } | |
faee3b18 | 947 | //__________________________________________ |
948 | void AliHFEpidQA:: CheckEvent(){ | |
949 | // | |
950 | // check some event variables | |
951 | // | |
70da6c5a | 952 | |
faee3b18 | 953 | // check the T0 as a function of run number (less than one bin per number |
954 | Double_t t0 = fESDpid->GetTOFResponse().GetTimeZero(); | |
955 | Int_t run = fEvent->GetRunNumber(); | |
956 | Double_t data[2] = {run*1.0, t0*1000.}; | |
957 | fOutput->Fill("hEvent_T0", data); | |
958 | ||
959 | ||
960 | } | |
70da6c5a | 961 | //__________________________________________ |
962 | TList *AliHFEpidQA::GetOutput(){ | |
963 | // | |
964 | // Getter for Output histograms | |
965 | // | |
966 | return fOutput->GetList(); | |
967 | } | |
968 | ||
969 | //__________________________________________ | |
970 | TList *AliHFEpidQA::GetV0pidQA(){ | |
971 | // | |
972 | // Getter for V0 PID QA histograms | |
973 | // | |
974 | return fV0pid->GetListOfQAhistograms(); | |
975 | } | |
976 | ||
977 | //__________________________________________ | |
978 | TList *AliHFEpidQA::GetV0pidMC(){ | |
979 | // | |
980 | // Getter for V0 PID QA histograms | |
981 | // | |
982 | if(fV0pidMC) | |
983 | return fV0pidMC->GetListOfQAhistograms(); | |
984 | return NULL; | |
985 | } | |
986 | ||
987 | //__________________________________________ | |
988 | void AliHFEpidQA::RecalculateTRDpid(AliESDtrack * /*track*/, Double_t * /*pidProbs*/) const{ | |
faee3b18 | 989 | // fTRDpidResponse->MakePID(track); |
990 | // track->GetTRDpid(pidProbs); | |
70da6c5a | 991 | } |
992 | ||
993 | //__________________________________________ | |
994 | void AliHFEpidQA::RecalculateTRDpid(AliAODTrack * /*track*/, Double_t * /*pidProbs*/) const{ | |
faee3b18 | 995 | // fTRDpidResponse->MakePID(track, pidProbs); |
70da6c5a | 996 | } |
997 | //___________________________________________________________________ | |
e156c3bb | 998 | Float_t AliHFEpidQA::TOFbeta(const AliESDtrack * const track) const { |
70da6c5a | 999 | // computes the TOF beta |
1000 | Double_t l = track->GetIntegratedLength(); // cm | |
1001 | Double_t t = track->GetTOFsignal(); | |
faee3b18 | 1002 | Double_t t0 = fESDpid->GetTOFResponse().GetTimeZero(); // ps |
1003 | ||
1004 | //printf("-D: l: %f, t: %f, t0: %f\n", l, t, t0); | |
1005 | ||
70da6c5a | 1006 | if(l < 360. || l > 800.) return 0.; |
1007 | if(t <= 0.) return 0.; | |
faee3b18 | 1008 | if(t0 >999990.0) return 0.; |
1009 | ||
70da6c5a | 1010 | |
1011 | t -= t0; // subtract the T0 | |
1012 | ||
1013 | l *= 0.01; // cm ->m | |
1014 | t *= 1e-12; //ps -> s | |
1015 | ||
1016 | ||
1017 | Double_t v = l / t; | |
1018 | Float_t beta = v / TMath::C(); | |
1019 | ||
1020 | return beta; | |
1021 | } | |
faee3b18 | 1022 | //____________________________________________ |
3a72645a | 1023 | Int_t AliHFEpidQA::GetMaxPID(const Double_t *pidProbs) const { |
faee3b18 | 1024 | // |
1025 | // return the index of maximal PID probability | |
1026 | // | |
1027 | Int_t ix = -1; | |
1028 | Double_t tmp = 0.2; | |
1029 | for(Int_t i=0; i<5; ++i){ | |
1030 | if(pidProbs[i] > tmp){ | |
1031 | ix = i; | |
1032 | tmp = pidProbs[i]; | |
1033 | } | |
1034 | } | |
1035 | return ix; | |
1036 | } | |
1037 | //_____________________________________________ | |
1038 | Int_t AliHFEpidQA::GetPDG(Int_t species){ | |
1039 | // | |
1040 | // return the PDG particle code | |
1041 | // | |
1042 | ||
1043 | Int_t pdg = 0; | |
1044 | ||
1045 | switch(species){ | |
3a72645a | 1046 | case AliHFEV0cuts::kRecoElectron: |
faee3b18 | 1047 | pdg = TMath::Abs(kElectron); |
1048 | break; | |
3a72645a | 1049 | case AliHFEV0cuts::kRecoPionK0: |
faee3b18 | 1050 | pdg = TMath::Abs(kPiPlus); |
1051 | break; | |
3a72645a | 1052 | case AliHFEV0cuts::kRecoPionL: |
faee3b18 | 1053 | pdg = TMath::Abs(kPiPlus); |
1054 | break; | |
3a72645a | 1055 | case AliHFEV0cuts::kRecoProton: |
faee3b18 | 1056 | pdg = TMath::Abs(kProton); |
1057 | break; | |
1058 | default: // non investigated species | |
1059 | AliDebug(3, "Species not recognised"); | |
1060 | return 0; | |
1061 | } | |
1062 | ||
1063 | return pdg; | |
1064 | ||
1065 | } | |
1066 | ||
1067 | //_____________________________________________ | |
3a72645a | 1068 | TObjArray * AliHFEpidQA::MakeTrackList(const TObjArray *tracks) const { |
faee3b18 | 1069 | // |
1070 | // convert list of AliHFEV0Info into a list of AliVParticle | |
1071 | // | |
1072 | TObjArray *output = new TObjArray; | |
1073 | TIter trackInfos(tracks); | |
1074 | AliHFEV0info *trackInfo = NULL; | |
1075 | while((trackInfo = dynamic_cast<AliHFEV0info *>(trackInfos()))) | |
1076 | output->Add(trackInfo->GetTrack()); | |
1077 | ||
1078 | return output; | |
1079 | } | |
1080 | ||
1081 | //_____________________________________________ | |
3a72645a | 1082 | TObjArray * AliHFEpidQA::MakeCleanListElectrons(const TObjArray *electrons) const { |
faee3b18 | 1083 | // |
1084 | // Cleanup electron sample using TPC PID | |
1085 | // PID requirement will allways be implemented to the pair | |
1086 | // Strategy | |
1087 | // | |
1088 | TObjArray *tracks = new TObjArray; | |
1089 | TIter candidates(electrons); | |
fd282aa0 | 1090 | AliESDEvent *esd; |
1091 | //AliAODEvent *aod; | |
faee3b18 | 1092 | AliHFEV0info *hfetrack; |
e3fc062d | 1093 | // const Double_t kSigmaTight = 1.; |
1094 | // const Double_t kSigmaLoose = 4.; | |
3a72645a | 1095 | const Double_t kSigmaTight = 2.; |
1096 | const Double_t kSigmaLoose = 2.; | |
1097 | const Double_t shift = 0.57; | |
faee3b18 | 1098 | if((esd = dynamic_cast<AliESDEvent *>(fEvent))){ |
1099 | AliESDtrack *track = NULL, *partnerTrack = NULL; | |
1100 | while((hfetrack = dynamic_cast<AliHFEV0info *>(candidates()))){ | |
1101 | track = dynamic_cast<AliESDtrack *>(hfetrack->GetTrack()); | |
bf892a6a | 1102 | if(!track) continue; |
faee3b18 | 1103 | partnerTrack = esd->GetTrack(hfetrack->GetPartnerID()); |
e3fc062d | 1104 | Double_t nSigmaTrack = TMath::Abs(fESDpid->NumberOfSigmasTPC(track, AliPID::kElectron) - shift); |
1105 | Double_t nSigmaPartner = TMath::Abs(fESDpid->NumberOfSigmasTPC(partnerTrack, AliPID::kElectron) - shift); | |
1106 | if((nSigmaTrack < kSigmaTight && nSigmaPartner < kSigmaLoose) || (nSigmaTrack < kSigmaLoose && nSigmaPartner < kSigmaTight)) | |
faee3b18 | 1107 | tracks->Add(track); |
1108 | } | |
3e03bd67 | 1109 | } |
1110 | /*else { | |
faee3b18 | 1111 | aod = dynamic_cast<AliAODEvent *>(fEvent); |
bf892a6a | 1112 | if(!aod) return NULL; |
1113 | //AliAODTrack *track = NULL, *partnerTrack = NULL; | |
faee3b18 | 1114 | while((hfetrack = dynamic_cast<AliHFEV0info *>(candidates()))){ |
bf892a6a | 1115 | if(!hfetrack) continue; |
1116 | //track = dynamic_cast<AliAODTrack *>(hfetrack->GetTrack()); | |
1117 | //partnerTrack = aod->GetTrack(hfetrack->GetPartnerID()); | |
faee3b18 | 1118 | // will be coming soon |
1119 | } | |
3e03bd67 | 1120 | }*/ |
faee3b18 | 1121 | return tracks; |
1122 | } | |
1123 | //___________________________________________________________ | |
3a72645a | 1124 | void AliHFEpidQA::CheckTenderV0pid(const TObjArray * const particles, Int_t species){ |
faee3b18 | 1125 | |
1126 | // | |
1127 | // retrieve the status bits from the TObject used to flag | |
1128 | // the V0 daughter tracks and return the found PID | |
1129 | // 0 - electron, 1 - pion, 2 - proton | |
1130 | // | |
1131 | ||
1132 | const Int_t id[5] = {0, 1, 1, -1, 2}; // convert AliHFEpid to simple 0, 1, 2 | |
1133 | ||
1134 | AliVParticle *recTrack = NULL; | |
1135 | TIter trackIter(particles); | |
1136 | while((recTrack = dynamic_cast<AliVParticle *>(trackIter()))){ | |
1137 | if(!TString(recTrack->IsA()->GetName()).CompareTo("AliESDtrack")){ | |
1138 | // case ESD | |
1139 | AliESDtrack *track = dynamic_cast<AliESDtrack *>(recTrack); | |
1140 | if(!track) continue; | |
1141 | Int_t tPID = GetTenderV0pid(track); | |
1142 | fOutput->Fill("h_tender_check_01", id[species]*1.0, tPID); | |
1143 | fOutput->Fill("h_tender_check_02", id[species], id[species]*1.0, tPID); | |
1144 | ||
1145 | } //.. case ESD | |
1146 | ||
1147 | }//.. iterate of tracks | |
1148 | } | |
1149 | //___________________________________________________________ | |
1150 | Int_t AliHFEpidQA::GetTenderV0pid(AliESDtrack * const track){ | |
1151 | // | |
1152 | // retrieve the PID nformation stored in the status flags by the train tender | |
1153 | // | |
1154 | ||
1155 | ||
1156 | Int_t pid = -1; | |
1157 | if(!track){ | |
1158 | return pid; | |
1159 | } | |
1160 | ||
1161 | Int_t nTimes = 0; | |
1162 | ||
1163 | if(track->TestBit(2<<14)){ | |
1164 | pid = 0; | |
1165 | nTimes++; | |
1166 | } | |
1167 | if(track->TestBit(2<<15)){ | |
1168 | pid = 1; | |
1169 | nTimes++; | |
1170 | } | |
1171 | if(track->TestBit(2<<16)){ | |
1172 | pid = 2; | |
1173 | nTimes++; | |
1174 | } | |
1175 | ||
1176 | if(nTimes > 1){ | |
1177 | AliWarning("V0 track labeled multiple times by the V0 tender"); | |
1178 | pid = -1; | |
1179 | } | |
1180 | ||
1181 | return pid; | |
1182 | ||
1183 | } | |
1184 | //___________________________________________________________ | |
1185 | Double_t AliHFEpidQA::TRDlikeTracklet(Int_t layer, AliESDtrack * const track, Double_t *likelihood){ | |
1186 | // | |
1187 | // compute the TRD electron likelihoods for 1 tracklet | |
1188 | // based on teh AliTRDpidRecalculator in train/until/tender | |
1189 | // returns sum of the likelihoods (which should be 1) | |
1190 | // | |
1191 | ||
1192 | const Double_t cScaleGain = 1./ 16000.; | |
1193 | const Float_t pBins[11] ={0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0}; // momentum bins | |
1194 | ||
1195 | if(!track) return kFALSE; | |
1196 | Float_t p = track->GetTRDmomentum(layer); // momentum for a tracklet in the ESDtrack | |
1197 | if(p < 0) return kFALSE; | |
1198 | ||
1199 | Int_t mombin = TRDmomBin(p); // momentum bin | |
bf892a6a | 1200 | if(mombin < 0) return -1.; |
faee3b18 | 1201 | Float_t dEdxTRDsum = 0; // dEdxTRDsum for checking if tracklet is available |
1202 | Float_t dEdxTRD[8]; // dEdx for a tracklet in the ESD slices | |
1203 | Double_t ddEdxTRD[8]; // dEdx as Double_t for TMultiLayerPerceptron::Evaluate() | |
1204 | ||
1205 | Double_t prob[AliPID::kSPECIES]; // probabilities for all species in all layers | |
1206 | ||
1207 | for(Int_t is = 0; is < AliPID::kSPECIES; is++){ | |
1208 | likelihood[is] = 0.2; // init probabilities | |
1209 | prob[is] = 0.2; | |
1210 | } | |
1211 | ||
1212 | Double_t sum = 0.; | |
1213 | ||
1214 | for(Int_t islice = 0; islice<8; islice++){ | |
1215 | dEdxTRD[islice]=0.; // init dE/dx | |
1216 | ddEdxTRD[islice]=0.; // init dE/dx | |
1217 | dEdxTRD[islice]=track->GetTRDslice(layer,islice); // get the dE/dx values | |
1218 | dEdxTRDsum += dEdxTRD[islice]; | |
1219 | ddEdxTRD[islice]=(Double_t)dEdxTRD[islice]*cScaleGain; // rescale dE/dx | |
1220 | ||
1221 | } | |
1222 | for(Int_t is = 0; is < AliPID::kSPECIES; is++){ | |
1223 | Double_t probloc1, probloc2; | |
1224 | if(mombin == 0 && mombin < pBins[0]){ // calculate PID for p > 0.6 GeV/c | |
1225 | prob[is] = fNet[mombin]->Evaluate(is, ddEdxTRD); | |
1226 | } | |
1227 | else if(mombin == 10 && mombin >= pBins[10]){ // calculate PID for p >= 10.0 GeV/c | |
1228 | prob[is] = fNet[mombin]->Evaluate(is, ddEdxTRD); | |
1229 | } | |
1230 | else{ // standard calculation | |
1231 | Int_t mombin1 = 0, mombin2 = 0; // lower and upper momentum bin | |
1232 | if(p < pBins[mombin]) {mombin1 = mombin -1; mombin2 = mombin;} | |
1233 | if(p >= pBins[mombin]) {mombin1 = mombin; mombin2 = mombin+1;} | |
1234 | probloc1 = fNet[mombin1]->Evaluate(is, ddEdxTRD); | |
1235 | probloc2 = fNet[mombin2]->Evaluate(is, ddEdxTRD); | |
1236 | // weighting of the probability with regard to the track momentum | |
1237 | prob[is] = probloc1 + (probloc2-probloc1)*(p-pBins[mombin1])/(pBins[mombin2]-pBins[mombin1]); | |
1238 | } | |
1239 | likelihood[is] = prob[is]; | |
1240 | sum += likelihood[is]; | |
1241 | } | |
1242 | ||
1243 | return sum; | |
1244 | } | |
1245 | //__________________________________________________________________________ | |
6555e2ad | 1246 | Int_t AliHFEpidQA::TRDmomBin(Double_t p) const { |
faee3b18 | 1247 | // |
1248 | // compute the momentum bin position | |
1249 | // | |
1250 | ||
1251 | const Float_t pBins[11] ={0.6, 0.8, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 8.0, 10.0}; // momentum bins | |
1252 | ||
1253 | Int_t pBin1 = -1; // check bin1 | |
1254 | Int_t pBin2 = -1; // check bin2 | |
1255 | ||
1256 | if(p < 0) return -1; // return -1 if momentum < 0 | |
1257 | if(p < pBins[0]) return 0; // smallest momentum bin | |
1258 | if(p >= pBins[10]) return 10; // largest momentum bin | |
1259 | ||
1260 | ||
1261 | // calculate momentum bin for non extremal momenta | |
1262 | for(Int_t iMomBin = 1; iMomBin < 11; iMomBin++){ | |
1263 | if(p < pBins[iMomBin]){ | |
1264 | pBin1 = iMomBin - 1; | |
1265 | pBin2 = iMomBin; | |
1266 | } | |
1267 | else | |
1268 | continue; | |
1269 | ||
1270 | if(p - pBins[pBin1] >= pBins[pBin2] - p){ | |
1271 | return pBin2; | |
1272 | } | |
1273 | else{ | |
1274 | return pBin1; | |
1275 | } | |
1276 | } | |
1277 | ||
1278 | return -1; | |
1279 | ||
1280 | ||
1281 | } | |
faee3b18 | 1282 | |
e17c1f86 | 1283 | //___________________________________________________ |
1284 | Int_t AliHFEpidQA::GetCentrality(AliVEvent* const fInputEvent){ | |
1285 | // | |
1286 | // Recover the centrality of the event from ESD or AOD | |
1287 | // | |
1288 | Int_t bin = -1; | |
1289 | ||
1290 | if(fIsPbPb) | |
1291 | { | |
1292 | // Centrality | |
1293 | AliCentrality *centrality = fInputEvent->GetCentrality(); | |
1294 | Float_t fCentralityFtemp = -1; | |
1295 | fCentralityFtemp = centrality->GetCentralityPercentile("V0M"); | |
1296 | Float_t centralityLimits[12] = {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.}; | |
1297 | for(Int_t ibin = 0; ibin < 11; ibin++){ | |
1298 | if(fCentralityFtemp >= centralityLimits[ibin] && fCentralityFtemp < centralityLimits[ibin+1]){ | |
1299 | bin = ibin; | |
1300 | break; | |
1301 | } | |
1302 | } | |
1303 | } | |
1304 | AliDebug(2, Form("Centrality class %d\n", bin)); | |
1305 | ||
1306 | ||
1307 | return bin; | |
1308 | ||
1309 | } | |
1310 | ||
1311 | //___________________________________________________ | |
a8ef1999 | 1312 | Int_t AliHFEpidQA::GetMultiplicityITS(AliVEvent* const fInputEvent) const |
e17c1f86 | 1313 | { |
1314 | // | |
1315 | // Definition of the Multiplicity according to the JPSI group (F. Kramer) | |
1316 | // | |
1317 | Int_t nTracklets = 0; | |
1318 | Int_t nAcc = 0; | |
1319 | Double_t etaRange = 1.6; | |
1320 | Int_t bin = -1; | |
1321 | ||
1322 | nTracklets = ((AliESDEvent*) fInputEvent)->GetMultiplicity()->GetNumberOfTracklets(); | |
1323 | for (Int_t nn = 0; nn < nTracklets; nn++) { | |
1324 | Double_t eta = ((AliESDEvent*)fInputEvent)->GetMultiplicity()->GetEta(nn); | |
1325 | if (TMath::Abs(eta) < etaRange) nAcc++; | |
1326 | } | |
1327 | ||
1328 | Int_t itsMultiplicity = nAcc; | |
1329 | ||
1330 | Int_t multiplicityLimits[8] = {0, 1, 9, 17, 25, 36, 60, 500}; | |
1331 | for(Int_t ibin = 0; ibin < 7; ibin++){ | |
1332 | if(itsMultiplicity >= multiplicityLimits[ibin] && itsMultiplicity < multiplicityLimits[ibin + 1]){ | |
1333 | bin = ibin; | |
1334 | break; | |
1335 | } | |
1336 | } | |
1337 | ||
1338 | return bin; | |
1339 | ||
1340 | } | |
faee3b18 | 1341 |