]>
Commit | Line | Data |
---|---|---|
a6fd3cfe | 1 | #include "AliAnalysisTaskPIDFluctuation.h" |
2 | #include "AliInputEventHandler.h" | |
3 | #include "AliAnalysisManager.h" | |
4 | #include "AliESDEvent.h" | |
5 | #include "AliESDtrack.h" | |
6 | #include "AliAODEvent.h" | |
7 | #include "AliAODTrack.h" | |
8 | #include "AliCentrality.h" | |
9 | #include "AliPIDResponse.h" | |
10 | #include "AliESDpid.h" | |
11 | #include "AliAODpidUtil.h" | |
12 | #include "TFile.h" | |
13 | #include "TKey.h" | |
14 | #include "TList.h" | |
15 | #include "TH1F.h" | |
16 | #include "TH2F.h" | |
17 | #include "TH3F.h" | |
18 | #include "TCanvas.h" | |
19 | #include "AliESDtrackCuts.h" | |
20 | ||
21 | /* | |
22 | * Event by event PID fluctuation analysis | |
23 | * author: Roberto Preghenella (R+) | |
24 | * email: preghenella@bo.infn.it | |
25 | * | |
26 | */ | |
27 | ||
28 | ClassImp(AliAnalysisTaskPIDFluctuation) | |
29 | ||
30 | //________________________________________________________________________ | |
31 | ||
32 | const Char_t *AliAnalysisTaskPIDFluctuation::fgkEventCounterName[AliAnalysisTaskPIDFluctuation::kNEventCounters] = { | |
33 | "AllEvents", | |
34 | "PhysicsSelection", | |
35 | "PrimayVertex", | |
36 | "PrimayVertexSPD", | |
37 | "VertexAccepted", | |
38 | "GoodCentrality", | |
39 | "AcceptedEvents" | |
40 | }; | |
41 | ||
42 | const Char_t *AliAnalysisTaskPIDFluctuation::fgkEventCounterTitle[AliAnalysisTaskPIDFluctuation::kNEventCounters] = { | |
43 | "all events", | |
44 | "physics selection", | |
45 | "primary vertex", | |
46 | "primary vertex SPD", | |
47 | "vertex accepted", | |
48 | "good centrality", | |
49 | "accepted events" | |
50 | }; | |
51 | ||
52 | const Char_t *AliAnalysisTaskPIDFluctuation::fgkSparseDataName[AliAnalysisTaskPIDFluctuation::kNSparseData] = { | |
53 | "centV0M", | |
54 | "centTRK", | |
55 | "Nch", | |
56 | "Nplus", | |
57 | "Nminus", | |
58 | "Npi", | |
59 | "Npiplus", | |
60 | "Npiminus", | |
61 | "Nka", | |
62 | "Nkaplus", | |
63 | "Nkaminus", | |
64 | "Npr", | |
65 | "Nprplus", | |
66 | "Nprminus" | |
67 | }; | |
68 | ||
69 | const Char_t *AliAnalysisTaskPIDFluctuation::fgkSparseDataTitle[AliAnalysisTaskPIDFluctuation::kNSparseData] = { | |
70 | "centrality percentile (V0M)", | |
71 | "centrality percentile (TRK)", | |
72 | "N_{charged}", | |
73 | "N_{+}", | |
74 | "N_{-}", | |
75 | "N_{#pi}", | |
76 | "N_{#pi^{+}}", | |
77 | "N_{#pi^{-}}", | |
78 | "N_{K}", | |
79 | "N_{K^{+}}", | |
80 | "N_{K^{-}}", | |
81 | "N_{p+#bar{p}}", | |
82 | "N_{p}", | |
83 | "N_{#bar{p}}" | |
84 | }; | |
85 | ||
86 | //________________________________________________________________________ | |
87 | ||
88 | AliAnalysisTaskPIDFluctuation::AliAnalysisTaskPIDFluctuation(const Char_t *name) : | |
89 | AliAnalysisTaskSE(name), | |
90 | fPIDMethod(kTPCTOF), | |
91 | fESDtrackCuts(NULL), | |
92 | fAODfilterBit(AliAODTrack::kTrkGlobal), | |
93 | fEtaMin(-0.8), | |
94 | fEtaMax(0.8), | |
95 | fPtMin(0.3), | |
96 | fPtMax(1.5), | |
97 | fPID(NULL), | |
036e5409 | 98 | fHistoList(NULL), |
99 | fHistoEventCounter(NULL), | |
100 | fHistoAcceptedTracks(NULL), | |
101 | fHistoTOFMatchedTracks(NULL), | |
102 | fHistoTPCdEdx(NULL), | |
103 | fHistoTPCdEdx_inclusive(NULL), | |
104 | fHistoTOFbeta(NULL), | |
105 | fHistoCorrelation(NULL) | |
a6fd3cfe | 106 | { |
107 | ||
108 | /* | |
109 | * default constructor | |
110 | */ | |
111 | ||
112 | DefineOutput(1, TList::Class()); | |
113 | } | |
114 | ||
115 | //________________________________________________________________________ | |
116 | ||
117 | AliAnalysisTaskPIDFluctuation::~AliAnalysisTaskPIDFluctuation() | |
118 | { | |
119 | ||
120 | /* | |
121 | * default destructor | |
122 | */ | |
123 | ||
124 | if (fPID) delete fPID; | |
125 | if (fHistoList) delete fHistoList; | |
126 | ||
127 | } | |
128 | ||
129 | //________________________________________________________________________ | |
130 | ||
131 | void AliAnalysisTaskPIDFluctuation::UserCreateOutputObjects() | |
132 | { | |
133 | ||
134 | /* | |
135 | * user create output objects | |
136 | */ | |
137 | ||
138 | fHistoList = new TList(); | |
139 | fHistoList->SetOwner(kTRUE); | |
140 | ||
141 | fHistoEventCounter = new TH1F("hHistoEventCounter", "", kNEventCounters, 0, kNEventCounters); | |
142 | for (Int_t ievc = 0; ievc < kNEventCounters; ievc++) | |
143 | fHistoEventCounter->GetXaxis()->SetBinLabel(ievc + 1, fgkEventCounterTitle[ievc]); | |
144 | fHistoList->Add(fHistoEventCounter); | |
145 | ||
146 | fHistoAcceptedTracks = new TH2F("hHistoAcceptedTracks", ";p_{T} (GeV/c)", 10, 0., 100., 50, 0., 5.); | |
147 | fHistoList->Add(fHistoAcceptedTracks); | |
148 | ||
149 | fHistoTOFMatchedTracks = new TH2F("hHistoTOFMatchedTracks", ";p_{T} (GeV/c)", 10, 0., 100., 50, 0., 5.); | |
150 | fHistoList->Add(fHistoTOFMatchedTracks); | |
151 | ||
152 | fHistoTPCdEdx = new TH3F("hHistoTPCdEdx", ";centrality percentile;p_{TPCin} (GeV/c);dE/dx (au.)", 10, 0., 100., 50, 0., 5., 500, 0., 500.); | |
153 | fHistoList->Add(fHistoTPCdEdx); | |
154 | ||
155 | fHistoTPCdEdx_inclusive = new TH3F("hHistoTPCdEdx_inclusive", ";centrality percentile;p_{TPCin} (GeV/c);dE/dx (au.)", 10, 0., 100., 50, 0., 5., 500, 0., 500.); | |
156 | fHistoList->Add(fHistoTPCdEdx_inclusive); | |
157 | ||
158 | fHistoTOFbeta = new TH3F(Form("hHistoTOFbeta"), ";centrality percentile;p (GeV/c);v/c", 10, 0., 100., 50, 0., 5., 500, 0.1, 1.1); | |
159 | fHistoList->Add(fHistoTOFbeta); | |
160 | ||
161 | /* loop over species */ | |
162 | for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) { | |
163 | ||
164 | fHistoTPCdEdx_selected[ipart] = new TH3F(Form("hHistoTPCdEdx_selected_%s", AliPID::ParticleName(ipart)), ";centrality percentile;p_{TPCin} (GeV/c);dE/dx (au.)", 10, 0., 100., 50, 0., 5., 500, 0., 500.); | |
165 | fHistoList->Add(fHistoTPCdEdx_selected[ipart]); | |
166 | ||
167 | fHistoTOFbeta_selected[ipart] = new TH3F(Form("hHistoTOFbeta_selected_%s", AliPID::ParticleName(ipart)), ";centrality percentile;p (GeV/c);v/c", 10, 0., 100., 50, 0., 5., 500, 0.1, 1.1); | |
168 | fHistoList->Add(fHistoTOFbeta_selected[ipart]); | |
169 | ||
170 | fHistoNSigmaTPC[ipart] = new TH3F(Form("hHistoNSigmaTPC_%s", AliPID::ParticleName(ipart)), Form(";centrality percentile;p_{T} (GeV/c); N_{#sigma-TPC}^{%s}", AliPID::ParticleLatexName(ipart)), 10, 0., 100., 50, 0., 5., 200, -10., 10.); | |
171 | fHistoList->Add(fHistoNSigmaTPC[ipart]); | |
172 | ||
173 | fHistoNSigmaTOF[ipart] = new TH3F(Form("hHistoNSigmaTOF_%s", AliPID::ParticleName(ipart)), Form(";centrality percentile;p_{T} (GeV/c); N_{#sigma-TOF}^{%s}", AliPID::ParticleLatexName(ipart)), 10, 0., 100., 50, 0., 5., 200, -10., 10.); | |
174 | fHistoList->Add(fHistoNSigmaTOF[ipart]); | |
175 | ||
176 | } /* end of loop over species */ | |
177 | ||
178 | Int_t fgSparseDataBins[kNSparseData] = { | |
179 | 20, | |
180 | 20, | |
181 | 5000, | |
182 | 2500, | |
183 | 2500, | |
184 | 3000, | |
185 | 1500, | |
186 | 1500, | |
187 | 1000, | |
188 | 500, | |
189 | 500, | |
190 | 500, | |
191 | 250, | |
192 | 250 | |
193 | }; | |
194 | Double_t fgSparseDataMin[kNSparseData] = { | |
195 | 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. | |
196 | }; | |
197 | Double_t fgSparseDataMax[kNSparseData] = { | |
198 | 100., 100., 5000., 2500., 2500., 3000., 1500., 1500., 1000., 500., 500., 500., 250., 250. | |
199 | }; | |
200 | ||
201 | ||
202 | fHistoCorrelation = new THnSparseI("hHistoCorrelation", "", kNSparseData, fgSparseDataBins, fgSparseDataMin, fgSparseDataMax); | |
203 | for (Int_t iaxis = 0; iaxis < kNSparseData; iaxis++) | |
204 | fHistoCorrelation->GetAxis(iaxis)->SetTitle(fgkSparseDataTitle[iaxis]); | |
205 | fHistoList->Add(fHistoCorrelation); | |
206 | ||
207 | PostData(1, fHistoList); | |
208 | } | |
209 | ||
210 | //________________________________________________________________________ | |
211 | ||
212 | void AliAnalysisTaskPIDFluctuation::UserExec(Option_t *) | |
213 | { | |
214 | ||
215 | /* | |
216 | * user exec | |
217 | */ | |
218 | ||
219 | /* get ESD event */ | |
220 | AliVEvent *event = InputEvent(); | |
221 | if (!event) return; | |
222 | ||
223 | /* accept event */ | |
224 | if (!AcceptEvent(event)) return; | |
225 | ||
226 | /* get centrality object and centrality */ | |
227 | AliCentrality *centrality = event->GetCentrality(); | |
228 | Float_t cent_v0m = centrality->GetCentralityPercentileUnchecked("V0M"); | |
229 | Float_t cent_trk = centrality->GetCentralityPercentileUnchecked("TRK"); | |
230 | ||
231 | /* init PID */ | |
a40bdd85 | 232 | if (!InitPID(event)) return; |
a6fd3cfe | 233 | |
234 | Bool_t pidFlag[AliPID::kSPECIES]; | |
235 | Int_t icharge; | |
236 | Int_t chargedCounts[2]; | |
237 | Int_t pidCounts[AliPID::kSPECIES][2]; | |
238 | for (Int_t i = 0; i < 2; i++) { | |
239 | chargedCounts[i] = 0; | |
240 | for (Int_t ii = 0; ii < 5; ii++) { | |
241 | pidCounts[ii][i] = 0; | |
242 | } | |
243 | } | |
244 | ||
245 | /* loop over tracks */ | |
246 | for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) { | |
247 | ||
248 | /* get track */ | |
249 | AliVTrack *track = (AliVTrack *)event->GetTrack(itrk); | |
250 | if (!track) continue; | |
251 | ||
252 | /* accept track */ | |
253 | if (!AcceptTrack(track)) continue; | |
254 | ||
255 | /* get charge */ | |
256 | icharge = track->Charge() > 0 ? 0 : 1; | |
257 | chargedCounts[icharge]++; | |
258 | ||
259 | /* make PID */ | |
260 | MakePID(track, pidFlag, cent_v0m); | |
261 | for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) | |
262 | if (pidFlag[ipart]) | |
263 | pidCounts[ipart][icharge]++; | |
264 | ||
265 | } | |
266 | ||
267 | /* fill histogram */ | |
268 | Double_t vsparse[kNSparseData]; | |
269 | vsparse[kCent_V0M] = cent_v0m; | |
270 | vsparse[kCent_TRK] = cent_trk; | |
271 | vsparse[kNch] = chargedCounts[0] + chargedCounts[1]; | |
272 | vsparse[kNch_plus] = chargedCounts[0]; | |
273 | vsparse[kNch_minus] = chargedCounts[1]; | |
274 | vsparse[kNpi] = pidCounts[AliPID::kPion][0] + pidCounts[AliPID::kPion][1]; | |
275 | vsparse[kNpi_plus] = pidCounts[AliPID::kPion][0]; | |
276 | vsparse[kNpi_minus] = pidCounts[AliPID::kPion][1]; | |
277 | vsparse[kNka] = pidCounts[AliPID::kKaon][0] + pidCounts[AliPID::kKaon][1]; | |
278 | vsparse[kNka_plus] = pidCounts[AliPID::kKaon][0]; | |
279 | vsparse[kNka_minus] = pidCounts[AliPID::kKaon][1]; | |
280 | vsparse[kNpr] = pidCounts[AliPID::kProton][0] + pidCounts[AliPID::kProton][1]; | |
281 | vsparse[kNpr_plus] = pidCounts[AliPID::kProton][0]; | |
282 | vsparse[kNpr_minus] = pidCounts[AliPID::kProton][1]; | |
283 | fHistoCorrelation->Fill(vsparse); | |
284 | ||
285 | PostData(1, fHistoList); | |
286 | } | |
287 | ||
288 | //___________________________________________________________ | |
289 | ||
290 | Bool_t | |
291 | AliAnalysisTaskPIDFluctuation::AcceptEvent(AliVEvent *event) const | |
292 | { | |
293 | /* | |
294 | * accept event | |
295 | */ | |
296 | ||
297 | /* fill histo event counter */ | |
298 | fHistoEventCounter->Fill(kAllEvents); | |
299 | ||
300 | /* check physics selection */ | |
301 | if (!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB)) return kFALSE; | |
302 | fHistoEventCounter->Fill(kPhysicsSelection); | |
303 | ||
304 | /* check primary vertex */ | |
305 | const AliVVertex *vertex = event->GetPrimaryVertex(); | |
306 | if (vertex->GetNContributors() < 1) { | |
307 | /* get ESD vertex SPD */ | |
308 | if (event->InheritsFrom("AliESDEvent")) { | |
309 | AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(event); | |
a40bdd85 | 310 | if (!esdevent) return kFALSE; |
a6fd3cfe | 311 | vertex = esdevent->GetPrimaryVertexSPD(); |
312 | } | |
313 | /* get AOD vertex SPD */ | |
314 | else if (event->InheritsFrom("AliAODEvent")) { | |
315 | AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(event); | |
a40bdd85 | 316 | if (!aodevent) return kFALSE; |
a6fd3cfe | 317 | vertex = aodevent->GetPrimaryVertexSPD(); |
318 | } | |
319 | if (vertex->GetNContributors() < 1) return kFALSE; | |
320 | fHistoEventCounter->Fill(kPrimaryVertexSPD); | |
321 | } | |
322 | fHistoEventCounter->Fill(kPrimaryVertex); | |
323 | ||
324 | /* check vertex position */ | |
325 | if (TMath::Abs(vertex->GetZ()) > 10.) return kFALSE; | |
326 | fHistoEventCounter->Fill(kVertexAccepted); | |
327 | ||
328 | /* get centrality object and check quality */ | |
329 | AliCentrality *centrality = event->GetCentrality(); | |
330 | if (centrality->GetQuality() != 0) return kFALSE; | |
331 | fHistoEventCounter->Fill(kGoodCentrality); | |
332 | ||
333 | /* event accepted */ | |
334 | fHistoEventCounter->Fill(kAcceptedEvents); | |
335 | return kTRUE; | |
336 | } | |
337 | ||
338 | //___________________________________________________________ | |
339 | ||
340 | Bool_t | |
341 | AliAnalysisTaskPIDFluctuation::AcceptTrack(AliVTrack *track) const | |
342 | { | |
343 | /* | |
344 | * accept track | |
345 | */ | |
346 | ||
347 | /* check ESD track cuts */ | |
348 | if (track->InheritsFrom("AliESDtrack")) { | |
349 | AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track); | |
a40bdd85 | 350 | if (!esdtrack) return kFALSE; |
a6fd3cfe | 351 | if (!fESDtrackCuts->AcceptTrack(esdtrack)) return kFALSE; |
352 | } | |
353 | /* check AOD filter bit */ | |
354 | else if (track->InheritsFrom("AliAODTrack")) { | |
355 | AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track); | |
a40bdd85 | 356 | if (!aodtrack) return kFALSE; |
a6fd3cfe | 357 | if (!aodtrack->TestFilterBit(fAODfilterBit)) return kFALSE; |
358 | } | |
359 | ||
360 | /* check eta range */ | |
361 | if (track->Eta() < fEtaMin || | |
362 | track->Eta() > fEtaMax) return kFALSE; | |
363 | /* check pt range */ | |
364 | if (track->Pt() < fPtMin || | |
365 | track->Pt() > fPtMax) return kFALSE; | |
366 | ||
367 | /* track accepted */ | |
368 | return kTRUE; | |
369 | } | |
370 | ||
371 | //___________________________________________________________ | |
372 | ||
373 | Bool_t | |
374 | AliAnalysisTaskPIDFluctuation::HasTPCPID(AliVTrack *track) const | |
375 | { | |
376 | /* | |
377 | * has TPC PID | |
378 | */ | |
379 | ||
380 | /* check PID signal */ | |
381 | if (track->GetTPCsignal() <= 0. || | |
382 | track->GetTPCsignalN() == 0) return kFALSE; | |
383 | return kTRUE; | |
384 | ||
385 | } | |
386 | ||
387 | //___________________________________________________________ | |
388 | ||
389 | Bool_t | |
390 | AliAnalysisTaskPIDFluctuation::HasTOFPID(AliVTrack *track) const | |
391 | { | |
392 | /* | |
393 | * has TOF PID | |
394 | */ | |
395 | ||
396 | /* check TOF matched track */ | |
397 | if (!(track->GetStatus() & AliESDtrack::kTOFout)|| | |
398 | !(track->GetStatus() & AliESDtrack::kTIME)) return kFALSE; | |
399 | return kTRUE; | |
400 | ||
401 | } | |
402 | ||
403 | //___________________________________________________________ | |
404 | ||
405 | Double_t | |
406 | AliAnalysisTaskPIDFluctuation::MakeTPCPID(AliVTrack *track, Double_t *nSigma) const | |
407 | { | |
408 | /* | |
409 | * make TPC PID | |
410 | * returns measured dEdx if PID available, otherwise -1. | |
411 | * fills nSigma array with TPC nsigmas for e, mu, pi, K, p | |
412 | */ | |
413 | ||
414 | /* check TPC PID */ | |
415 | if (!HasTPCPID(track)) return -1.; | |
416 | ||
417 | /* get TPC info */ | |
418 | Double_t ptpc = track->GetTPCmomentum(); | |
419 | Double_t dEdx = track->GetTPCsignal(); | |
420 | Int_t dEdxN = track->GetTPCsignalN(); | |
421 | ||
422 | /* loop over particles */ | |
423 | for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) { | |
424 | Double_t bethe = fPID->GetTPCResponse().GetExpectedSignal(ptpc, (AliPID::EParticleType)ipart); | |
425 | Double_t diff = dEdx - bethe; | |
426 | Double_t sigma = fPID->GetTPCResponse().GetExpectedSigma(ptpc, dEdxN, (AliPID::EParticleType)ipart); | |
427 | nSigma[ipart] = diff / sigma; | |
428 | } | |
429 | return dEdx; | |
430 | ||
431 | } | |
432 | ||
433 | //___________________________________________________________ | |
434 | ||
435 | Double_t | |
436 | AliAnalysisTaskPIDFluctuation::MakeTOFPID(AliVTrack *track, Double_t *nSigma) const | |
437 | { | |
438 | /* | |
439 | * make TOF PID | |
440 | * returns measured beta if PID available, otherwise -1. | |
441 | * fills nSigma array with TOF nsigmas for e, mu, pi, K, p | |
442 | */ | |
443 | ||
444 | /* check TOF PID */ | |
445 | if (!HasTOFPID(track)) return -1.; | |
446 | ||
447 | /* get TOF info */ | |
448 | Double_t p = track->P(); | |
449 | Double_t time = track->GetTOFsignal() - fPID->GetTOFResponse().GetStartTime(p); | |
450 | Double_t timei[5]; | |
451 | track->GetIntegratedTimes(timei); | |
452 | ||
453 | /* loop over particles */ | |
454 | for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) { | |
455 | Double_t timez = time - timei[ipart]; | |
456 | Double_t sigma = fPID->GetTOFResponse().GetExpectedSigma(p, timei[ipart], AliPID::ParticleMass(ipart)); | |
457 | nSigma[ipart] = timez / sigma; | |
458 | } | |
459 | ||
460 | return timei[0] / time; | |
461 | ||
462 | } | |
463 | ||
464 | //___________________________________________________________ | |
465 | ||
466 | void | |
467 | AliAnalysisTaskPIDFluctuation::MakePID(AliVTrack *track, Bool_t *pidFlag, Float_t centrality) const | |
468 | { | |
469 | /* | |
470 | * make PID | |
471 | * fills PID QA plots | |
472 | * fills pidFlag array with PID flags for e, mu, pi, K, p | |
473 | */ | |
474 | ||
475 | /* cut definitions | |
476 | (better put them as static variables so they can be changed from outside) */ | |
477 | Double_t fgTPCPIDmomcut[AliPID::kSPECIES] = {0., 0., 0.5, 0.5, 0.7}; | |
478 | Double_t fgTPCPIDsigmacut[AliPID::kSPECIES] = {0., 0., 2., 2., 2.}; | |
479 | Double_t fgTPCPIDcompcut[AliPID::kSPECIES] = {0., 0., 3., 3., 3.}; | |
480 | Double_t fgTOFPIDmomcut[AliPID::kSPECIES] = {0., 0., 1.5, 1.5, 2.0}; | |
481 | Double_t fgTOFPIDsigmacut[AliPID::kSPECIES] = {0., 0., 2., 2., 2.}; | |
482 | ||
483 | /* get momentum information */ | |
484 | Double_t p = track->P(); | |
485 | Double_t pt = track->Pt(); | |
486 | Double_t ptpc = track->GetTPCmomentum(); | |
487 | ||
488 | /* make pid and check if available */ | |
489 | Double_t nsigmaTPC[AliPID::kSPECIES]; | |
490 | Double_t nsigmaTOF[AliPID::kSPECIES]; | |
491 | Double_t dEdx = MakeTPCPID(track, nsigmaTPC); | |
492 | Double_t beta = MakeTOFPID(track, nsigmaTOF); | |
493 | Bool_t hasTPCPID = dEdx > 0.; | |
494 | Bool_t hasTOFPID = beta > 0.; | |
495 | ||
496 | /* check PID method */ | |
497 | if (fPIDMethod == kTPConly) hasTOFPID = kFALSE; // inhibit TOF PID | |
498 | if (fPIDMethod == kTOFonly) hasTPCPID = kFALSE; // inhibit TPC PID | |
499 | ||
500 | /* fill qa histos */ | |
501 | fHistoAcceptedTracks->Fill(centrality, pt); | |
502 | if (hasTPCPID) | |
503 | fHistoTPCdEdx->Fill(centrality, ptpc, dEdx); | |
504 | if (hasTOFPID) { | |
505 | fHistoTOFMatchedTracks->Fill(centrality, pt); | |
506 | fHistoTOFbeta->Fill(centrality, p, beta); | |
507 | } | |
508 | if (hasTPCPID && !hasTOFPID) | |
509 | fHistoTPCdEdx_inclusive->Fill(centrality, ptpc, dEdx); | |
510 | ||
511 | /* loop over species */ | |
512 | for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) { | |
513 | ||
514 | /* reset PID flag */ | |
515 | pidFlag[ipart] = kFALSE; | |
516 | ||
517 | /* fill qa histos */ | |
518 | if (hasTPCPID) | |
519 | fHistoNSigmaTPC[ipart]->Fill(centrality, pt, nsigmaTPC[ipart]); | |
520 | if (hasTOFPID) | |
521 | fHistoNSigmaTOF[ipart]->Fill(centrality, pt, nsigmaTOF[ipart]); | |
522 | ||
523 | /* combined PID cuts */ | |
524 | if (hasTPCPID && hasTOFPID) { | |
525 | if (pt < fgTOFPIDmomcut[ipart] && | |
526 | TMath::Abs(nsigmaTOF[ipart]) < fgTOFPIDsigmacut[ipart] && | |
527 | TMath::Abs(nsigmaTPC[ipart]) < fgTPCPIDcompcut[ipart]) { | |
528 | fHistoTOFbeta_selected[ipart]->Fill(centrality, p, beta); | |
529 | pidFlag[ipart] = kTRUE; | |
530 | } | |
531 | } | |
532 | /* TPC-only PID cuts */ | |
533 | else if (hasTPCPID && !hasTOFPID) { | |
534 | if (pt < fgTPCPIDmomcut[ipart] && | |
535 | TMath::Abs(nsigmaTPC[ipart]) < fgTPCPIDsigmacut[ipart]) { | |
536 | fHistoTPCdEdx_selected[ipart]->Fill(centrality, ptpc, dEdx); | |
537 | pidFlag[ipart] = kTRUE; | |
538 | } | |
539 | } | |
540 | /* TOF-only PID cuts */ | |
541 | else if (!hasTPCPID && hasTOFPID) { | |
542 | if (pt < fgTOFPIDmomcut[ipart] && | |
543 | TMath::Abs(nsigmaTOF[ipart]) < fgTOFPIDsigmacut[ipart]) { | |
544 | fHistoTOFbeta_selected[ipart]->Fill(centrality, p, beta); | |
545 | pidFlag[ipart] = kTRUE; | |
546 | } | |
547 | } | |
548 | ||
549 | } /* end of loop over species */ | |
550 | ||
551 | } | |
552 | ||
553 | //___________________________________________________________ | |
554 | ||
a40bdd85 | 555 | Bool_t |
a6fd3cfe | 556 | AliAnalysisTaskPIDFluctuation::InitPID(AliVEvent *event) |
557 | { | |
558 | /* | |
559 | * init PID | |
560 | */ | |
561 | ||
562 | /* create PID object if not there yet */ | |
563 | if (!fPID) { | |
564 | ||
565 | /* instance object */ | |
566 | Bool_t mcFlag = kFALSE; /*** WARNING: check whether is MC ***/ | |
567 | if (event->InheritsFrom("AliESDEvent")) | |
568 | fPID = new AliESDpid(mcFlag); | |
569 | else if (event->InheritsFrom("AliAODEvent")) | |
570 | fPID = new AliAODpidUtil(mcFlag); | |
a40bdd85 | 571 | else |
572 | return kFALSE; | |
573 | ||
a6fd3cfe | 574 | /* set OADB path */ |
575 | fPID->SetOADBPath("$ALICE_ROOT/OADB"); | |
576 | } | |
577 | ||
578 | /* init ESD PID */ | |
579 | Int_t recoPass = 2; /*** WARNING: need to set the recoPass somehow better ***/ | |
580 | fPID->InitialiseEvent(event, recoPass); /* warning: this apparently sets TOF time | |
581 | * resolution to some hardcoded value, | |
582 | * therefore we have to set correct | |
583 | * resolution value after this call */ | |
584 | ||
585 | /* set TOF response */ | |
586 | Double_t tofReso = 85.; /* ps */ /*** WARNING: need to set tofReso somehow better ***/ | |
587 | fPID->GetTOFResponse().SetTimeResolution(tofReso); | |
588 | fPID->GetTOFResponse().SetTrackParameter(0, 0.007); | |
589 | fPID->GetTOFResponse().SetTrackParameter(1, 0.007); | |
590 | fPID->GetTOFResponse().SetTrackParameter(2, 0.0); | |
591 | fPID->GetTOFResponse().SetTrackParameter(3, 30); | |
592 | ||
a40bdd85 | 593 | return kTRUE; |
a6fd3cfe | 594 | } |
595 | ||
596 | //___________________________________________________________ | |
597 | ||
598 | void | |
599 | AliAnalysisTaskPIDFluctuation::MeasureNuDyn(const Char_t *filename, Int_t i1, Int_t i2, Int_t centralityEstimator) | |
600 | { | |
601 | ||
602 | /* | |
603 | * measure nu | |
604 | */ | |
605 | ||
606 | printf("MeasureNuDyn: running for %s vs. %s\n", fgkSparseDataName[i1], fgkSparseDataName[i2]); | |
607 | ||
608 | /* get data */ | |
609 | TFile *filein = TFile::Open(filename); | |
610 | /* output */ | |
611 | TFile *fileout = TFile::Open(Form("MeasureNuDyn_%s_%s.%s", fgkSparseDataName[i1], fgkSparseDataName[i2], filename), "RECREATE"); | |
612 | ||
613 | /* loop over available containers */ | |
614 | TList *keylist = filein->GetListOfKeys(); | |
615 | for (Int_t ikey = 0; ikey < keylist->GetEntries(); ikey++) { | |
616 | ||
617 | /* get key and check */ | |
618 | TKey *key = (TKey *)keylist->At(ikey); | |
619 | TString contname = key->GetName(); | |
620 | if (!contname.BeginsWith("PIDFluctuation")) continue; | |
621 | ||
622 | /* get data */ | |
623 | TList *list = (TList *)filein->Get(contname.Data()); | |
624 | THnSparse *hsparse = (THnSparse *)list->FindObject("hHistoCorrelation"); | |
625 | ||
626 | /* create output directory and cd there */ | |
627 | fileout->mkdir(contname.Data()); | |
628 | fileout->cd(contname.Data()); | |
629 | ||
630 | /* loop over centralities */ | |
631 | Double_t centBin[kNCentralityBins + 1] = { | |
632 | 0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90. | |
633 | }; | |
634 | TH1D *hNu = new TH1D("hNu", ";centrality percentile;#nu", kNCentralityBins, centBin); | |
635 | TH1D *hNuStat = new TH1D("hNuStat", ";centrality percentile;#nu_{stat}", kNCentralityBins, centBin); | |
636 | TH1D *hNuDyn = new TH1D("hNuDyn", ";centrality percentile;#nu_{dyn}", kNCentralityBins, centBin); | |
637 | for (Int_t icent = 0; icent < kNCentralityBins; icent++) { | |
638 | ||
639 | /* select centrality range */ | |
640 | hsparse->GetAxis(centralityEstimator)->SetRangeUser(centBin[icent] + 0.001, centBin[icent + 1] - 0.001); | |
641 | /* projections */ | |
642 | TH1D *hcent = hsparse->Projection(centralityEstimator); | |
643 | TH1D *h1 = hsparse->Projection(i1); | |
644 | TH1D *h2 = hsparse->Projection(i2); | |
645 | TH2D *hcorr = hsparse->Projection(i2, i1); | |
646 | TH1D *hnu = new TH1D("hnu", "", 2000, -1., 1.); | |
647 | ||
648 | Double_t n, a, b, nev, meanx, rmsx, meany, rmsy, meannu, rmsnu; | |
649 | Double_t meanx_err, meany_err, meannu_err; | |
650 | ||
651 | /* compute mean values */ | |
652 | nev = 0.; meanx = 0.; meany = 0.; | |
653 | for (Int_t ibinx = 0; ibinx < hcorr->GetNbinsX(); ibinx++) | |
654 | for (Int_t ibiny = 0; ibiny < hcorr->GetNbinsY(); ibiny++) { | |
655 | n = hcorr->GetBinContent(ibinx + 1, ibiny + 1); | |
656 | if (n <= 0.) continue; | |
657 | meanx += n * ibinx; | |
658 | meany += n * ibiny; | |
659 | nev += n; | |
660 | } | |
661 | meanx /= nev; | |
662 | meany /= nev; | |
663 | // printf("nev = %f, meanx = %f, meany = %f\n", nev, meanx, meany); | |
664 | ||
665 | /* compute RMS values */ | |
666 | nev = 0.; rmsx = 0.; rmsy = 0.; | |
667 | for (Int_t ibinx = 0; ibinx < hcorr->GetNbinsX(); ibinx++) | |
668 | for (Int_t ibiny = 0; ibiny < hcorr->GetNbinsY(); ibiny++) { | |
669 | n = hcorr->GetBinContent(ibinx + 1, ibiny + 1); | |
670 | if (n <= 0.) continue; | |
671 | a = ibinx - meanx; | |
672 | rmsx += n * a * a; | |
673 | a = ibiny - meany; | |
674 | rmsy += n * a * a; | |
675 | nev += n; | |
676 | } | |
677 | rmsx /= nev; | |
678 | rmsx = TMath::Sqrt(rmsx); | |
679 | rmsy /= nev; | |
680 | rmsy = TMath::Sqrt(rmsy); | |
681 | // printf("nev = %f, rmsx = %f, rmsy = %f\n", nev, rmsx, rmsy); | |
682 | meanx_err = rmsx / TMath::Sqrt(nev); | |
683 | meany_err = rmsy / TMath::Sqrt(nev); | |
684 | // printf("nev = %f, meanx_err = %f, meany_err = %f\n", nev, meanx_err, meany_err); | |
685 | ||
686 | /* compute mean nu */ | |
687 | nev = 0.; meannu = 0.; | |
688 | for (Int_t ibinx = 0; ibinx < hcorr->GetNbinsX(); ibinx++) | |
689 | for (Int_t ibiny = 0; ibiny < hcorr->GetNbinsY(); ibiny++) { | |
690 | n = hcorr->GetBinContent(ibinx + 1, ibiny + 1); | |
691 | if (n <= 0.) continue; | |
692 | a = ibinx / meanx - ibiny / meany; | |
693 | meannu += n * a * a; | |
694 | hnu->Fill(a * a, n); | |
695 | nev += n; | |
696 | } | |
697 | meannu /= nev; | |
698 | // printf("nev = %f, meannu = %f\n", nev, meannu); | |
699 | ||
700 | /* compute RMS nu */ | |
701 | nev = 0.; rmsnu = 0.; | |
702 | for (Int_t ibinx = 0; ibinx < hcorr->GetNbinsX(); ibinx++) | |
703 | for (Int_t ibiny = 0; ibiny < hcorr->GetNbinsY(); ibiny++) { | |
704 | n = hcorr->GetBinContent(ibinx + 1, ibiny + 1); | |
705 | if (n <= 0.) continue; | |
706 | a = ibinx / meanx - ibiny / meany; | |
707 | b = a * a - meannu; | |
708 | rmsnu += n * b * b; | |
709 | nev += n; | |
710 | } | |
711 | rmsnu /= nev; | |
712 | rmsnu = TMath::Sqrt(rmsnu); | |
713 | // printf("nev = %f, rmsnu = %f\n", nev, rmsnu); | |
714 | meannu_err = rmsnu / TMath::Sqrt(nev); | |
715 | // printf("nev = %f, meannu_err = %f\n", nev, meannu_err); | |
716 | ||
717 | /* final calculations */ | |
718 | Double_t nu = meannu; | |
719 | Double_t nu_err = meannu_err; | |
720 | Double_t nu_stat = 1. / meanx + 1. / meany; | |
721 | Double_t meanx4 = meanx * meanx * meanx * meanx; | |
722 | Double_t meanx_err2 = meanx_err * meanx_err; | |
723 | Double_t meany4 = meany * meany * meany * meany; | |
724 | Double_t meany_err2 = meany_err * meany_err; | |
725 | Double_t nu_stat_err = TMath::Sqrt(meanx_err2 / meanx4 + meany_err2 / meany4); | |
726 | Double_t nu_dyn = nu - nu_stat; | |
727 | Double_t nu_dyn_err = TMath::Sqrt(nu_err * nu_err + nu_stat_err * nu_stat_err); | |
728 | ||
729 | /* setup final plots */ | |
730 | hNu->SetBinContent(icent + 1, nu); | |
731 | hNu->SetBinError(icent + 1, nu_err); | |
732 | hNuStat->SetBinContent(icent + 1, nu_stat); | |
733 | hNuStat->SetBinError(icent + 1, nu_stat_err); | |
734 | hNuDyn->SetBinContent(icent + 1, nu_dyn); | |
735 | hNuDyn->SetBinError(icent + 1, nu_dyn_err); | |
736 | ||
737 | hcent->Write(Form("hcent_cent%d", icent)); | |
738 | h1->Write(Form("h1_cent%d", icent)); | |
739 | h2->Write(Form("h2_cent%d", icent)); | |
740 | hcorr->Write(Form("hcorr_cent%d", icent)); | |
741 | hnu->Write(Form("hnu_cent%d", icent)); | |
742 | ||
743 | /* clean-up */ | |
744 | delete hcent; | |
745 | delete h1; | |
746 | delete h2; | |
747 | delete hcorr; | |
748 | delete hnu; | |
749 | } | |
750 | ||
751 | hNu->Write(); | |
752 | hNuStat->Write(); | |
753 | hNuDyn->Write(); | |
754 | } | |
755 | ||
756 | fileout->Close(); | |
757 | filein->Close(); | |
758 | ||
759 | } | |
760 | ||
761 | //___________________________________________________________ |