]>
Commit | Line | Data |
---|---|---|
c683985a | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: Satyajit Jena. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | ||
17 | //=========================================================================// | |
18 | // AliEbyE Analysis for Particle Ratio Fluctuation // | |
19 | // Deepika Rathee | Satyajit Jena // | |
20 | // drathee@cern.ch | sjena@cern.ch // | |
21 | // | |
22 | // V0.1 2013/03/25 Using THnSparse | |
23 | // V0.2 2013/04/03 Cleanup | |
24 | // V1.0 2013/04/10 Cleanup Bins for Less Memory | |
25 | // V1.1 2013/04/15 Bins Added | |
26 | // Todo: pp and pA, Mix Events | |
27 | //=========================================================================// | |
28 | ||
29 | #include "TChain.h" | |
30 | #include "TList.h" | |
31 | #include "TFile.h" | |
32 | #include "TTree.h" | |
33 | #include "TH1D.h" | |
34 | #include "TH2F.h" | |
35 | #include "TH3F.h" | |
36 | #include "TCanvas.h" | |
37 | #include "AliAnalysisTask.h" | |
38 | #include "AliAnalysisManager.h" | |
39 | #include "AliVEvent.h" | |
40 | #include "AliESD.h" | |
41 | #include "AliESDEvent.h" | |
42 | #include "AliAODEvent.h" | |
43 | #include "AliAODMCParticle.h" | |
44 | #include "AliAODMCHeader.h" | |
45 | #include "AliPIDResponse.h" | |
46 | #include "AliAODHeader.h" | |
47 | #include "AliAODpidUtil.h" | |
48 | #include "AliHelperPID.h" | |
49 | using std::endl; | |
50 | using std::cout; | |
51 | #include "AliEbyEParticleRatioFluctuationTask.h" | |
52 | ||
53 | ClassImp(AliEbyEParticleRatioFluctuationTask) | |
54 | ||
55 | //----------------------------------------------------------------------- | |
56 | AliEbyEParticleRatioFluctuationTask::AliEbyEParticleRatioFluctuationTask( const char *name ) : AliAnalysisTaskSE( name ), | |
57 | fThnList(NULL), | |
58 | fAnalysisType("AOD"), | |
59 | fAnalysisData("PbPb"), | |
60 | fCentralityEstimator("V0M"), | |
61 | fVxMax(3.), | |
62 | fVyMax(3.), | |
63 | fVzMax(10.), | |
64 | fDCAxy(2.4), | |
65 | fDCAz(3.2), | |
66 | fPtLowerLimit(0.2), | |
67 | fPtHigherLimit(5.), | |
68 | fEtaLowerLimit(-1.), | |
69 | fEtaHigherLimit(1.), | |
70 | fTPCNClus(80), | |
71 | fAODtrackCutBit(128), | |
72 | isQA(kFALSE), | |
73 | fDebug(kFALSE), | |
74 | fHelperPID(0x0), | |
75 | fEventCounter(NULL), | |
76 | fHistoCorrelation(NULL) { | |
77 | for(Int_t i = 0; i < 14; i++ ) fHistQA[i] = NULL; | |
78 | DefineOutput(1, TList::Class()); //! Connect Outpput.... | |
79 | } | |
80 | ||
81 | AliEbyEParticleRatioFluctuationTask::~AliEbyEParticleRatioFluctuationTask() { | |
82 | //! Cleaning up | |
83 | if (fThnList) delete fThnList; | |
84 | if (fHelperPID) delete fHelperPID; | |
85 | } | |
86 | ||
87 | //--------------------------------------------------------------------------------- | |
88 | void AliEbyEParticleRatioFluctuationTask::UserCreateOutputObjects() { | |
89 | fThnList = new TList(); | |
90 | fThnList->SetOwner(kTRUE); | |
91 | ||
92 | fEventCounter = new TH1D("fEventCounter","EventCounter", 300, 0.5,300.5); | |
93 | if (isQA) fThnList->Add(fEventCounter); | |
94 | ||
95 | fHistQA[0] = new TH2F("fHistQAvx", "Histo Vx Selected;Centrality;Vx", 100,0,100, 5000, -5., 5.); | |
96 | fHistQA[1] = new TH2F("fHistQAvy", "Histo Vy Selected;Centrality;Vy", 100,0,100, 5000, -5., 5.); | |
97 | fHistQA[2] = new TH2F("fHistQAvz", "Histo Vz Selected;Centrality;Vz", 100,0,100, 5000, -25., 25.); | |
98 | fHistQA[3] = new TH2F("fHistQAvxA", "Histo Vx;Centrality;Vx", 100,0,100, 5000, -5., 5.); | |
99 | fHistQA[4] = new TH2F("fHistQAvyA", "Histo Vy;Centrality;Vy", 100,0,100, 5000, -5., 5.); | |
100 | fHistQA[5] = new TH2F("fHistQAvzA", "Histo Vz;Centrality;Vz", 100,0,100, 5000, -25., 25.); | |
101 | ||
102 | fHistQA[6] = new TH2F("fHistQADcaXyA", "Histo DCAxy;Centrality;DCAxy",100,0,100, 600, -15., 15.); | |
103 | fHistQA[7] = new TH2F("fHistQADcaZA", "Histo DCAz;Centrality;DCAz ",100,0,100, 600, -15., 15.); | |
104 | fHistQA[8] = new TH2F("fHistQAPtA","p_{T} distribution;Centrality;p_{T}",100,0,100,1000,0,10); | |
105 | fHistQA[9] = new TH2F("fHistQAEtaA","#eta distribution;Centrality;#eta",100,0,100,240,-1.2,1.2); | |
106 | ||
107 | fHistQA[10] = new TH2F("fHistQADcaXy", "Histo DCAxy after Selected;Centrality;DCAxy", 100,0,100,600, -15., 15.); | |
108 | fHistQA[11] = new TH2F("fHistQADcaZ", "Histo DCAz Selected;Centrality;DCAz", 100,0,100,600, -15., 15.); | |
109 | fHistQA[12] = new TH2F("fHistQAPt","p_{T} distribution Selected;Centrality;p_{T}",100,0,100,1000,0,10); | |
110 | fHistQA[13] = new TH2F("fHistQAEta","#eta distribution Selected;Centrality;#eta",100,0,100, 240,-1.2,1.2); | |
111 | ||
112 | if (isQA) for(Int_t i = 0; i < 14; i++) fThnList->Add(fHistQA[i]); | |
113 | ||
114 | Int_t fgSparseDataBins[kNSparseData] = {100, 5000, 5000, 2500, 2500, 3000, 1500, 1500, 1000, 500, 500, 500, 250, 250}; | |
115 | Double_t fgSparseDataMin[kNSparseData] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}; | |
116 | Double_t fgSparseDataMax[kNSparseData] = {100.,5000.,5000.,2500.,2500.,3000.,1500.,1500.,1000.,500.,500.,500.,250.,250.}; | |
117 | ||
118 | const Char_t *fgkSparseDataTitle[] = {"centrality","RefMult","N_{ch}", "N_{+}","N_{-}","N_{#pi}", "N_{#pi^{+}}","N_{#pi^{-}}","N_{K}","N_{K^{+}}", "N_{K^{-}}","N_{p}","N_{p}","N_{#bar{p}}"}; | |
119 | ||
120 | fHistoCorrelation = new THnSparseI("fThnCorr", "", kNSparseData, fgSparseDataBins, fgSparseDataMin, fgSparseDataMax); | |
121 | for (Int_t iaxis = 0; iaxis < kNSparseData; iaxis++) | |
122 | fHistoCorrelation->GetAxis(iaxis)->SetTitle(fgkSparseDataTitle[iaxis]); | |
123 | ||
124 | if(!isQA) fThnList->Add(fHistoCorrelation); | |
125 | if(isQA) | |
126 | if (fHelperPID) | |
127 | fThnList->Add(fHelperPID); | |
128 | ||
129 | PostData(1, fThnList); | |
130 | } | |
131 | ||
132 | //---------------------------------------------------------------------------------- | |
133 | void AliEbyEParticleRatioFluctuationTask::UserExec( Option_t * ){ | |
134 | ||
135 | if (isQA) fEventCounter->Fill(1); | |
136 | ||
137 | AliAODEvent* event = dynamic_cast<AliAODEvent*>(InputEvent()); | |
138 | if (!event) { | |
139 | Printf("ERROR 01: AOD not found "); | |
140 | return; | |
141 | } | |
142 | ||
143 | Int_t gCent = -1; | |
144 | Float_t gRefMul = -1; | |
145 | ||
146 | AliAODHeader *aodHeader = event->GetHeader(); | |
147 | gCent = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data()); | |
148 | gRefMul = aodHeader->GetRefMultiplicity(); | |
149 | if (gCent < 0 || gCent > 100) return; | |
150 | if (isQA) fEventCounter->Fill(2); | |
151 | ||
152 | if (!AcceptEvent(event,gCent)) return; | |
153 | ||
154 | Int_t icharge = -1; | |
155 | Int_t gCharge[2]; | |
156 | Int_t gPid[3][2]; | |
157 | ||
158 | for (Int_t i = 0; i < 2; i++) { | |
159 | gCharge[i] = 0; | |
160 | for (Int_t ii = 0; ii < 3; ii++) { | |
161 | gPid[ii][i] = 0; | |
162 | } | |
163 | } | |
164 | ||
165 | ||
166 | if(fAnalysisType == "AOD") { | |
167 | if (isQA) { | |
168 | fEventCounter->Fill(5); | |
169 | fEventCounter->Fill(50+gCent); | |
170 | } | |
171 | ||
172 | for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) { | |
173 | AliAODTrack* track = dynamic_cast<AliAODTrack *>(event->GetTrack(itrk)); | |
174 | if (!track) continue; | |
175 | ||
176 | if (!AcceptTrack(track, gCent)) continue; | |
177 | ||
178 | Int_t a = fHelperPID->GetParticleSpecies((AliVTrack*) track,kTRUE); | |
179 | ||
180 | if(a < 0 || a > 2) continue; | |
181 | icharge = track->Charge() > 0 ? 0 : 1; | |
182 | gCharge[icharge]++; | |
183 | gPid[a][icharge]++; | |
184 | } | |
185 | } | |
186 | else if(fAnalysisType == "MCAOD") { | |
187 | TClonesArray *arrayMC= 0; | |
188 | arrayMC = dynamic_cast<TClonesArray*> (event->GetList()->FindObject(AliAODMCParticle::StdBranchName())); | |
189 | if (!arrayMC) { | |
190 | Printf("Error: MC particles branch not found!\n"); | |
191 | return; | |
192 | } | |
193 | AliAODMCHeader *mcHdr=0; | |
194 | mcHdr=(AliAODMCHeader*)event->GetList()->FindObject(AliAODMCHeader::StdBranchName()); | |
195 | if(!mcHdr) { | |
196 | Printf("MC header branch not found!\n"); | |
197 | return; | |
198 | } | |
199 | ||
200 | if (isQA) { | |
201 | fEventCounter->Fill(5); | |
202 | fEventCounter->Fill(50+gCent); | |
203 | } | |
204 | ||
205 | Int_t nMC = arrayMC->GetEntries(); | |
206 | for (Int_t iMC = 0; iMC < nMC; iMC++) { | |
207 | AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC); | |
208 | if(!AcceptMCTrack(partMC, gCent)) continue; | |
209 | Int_t a = fHelperPID->GetMCParticleSpecie((AliVEvent*) event,(AliVTrack*)partMC,1); | |
210 | if(a < 0 || a > 2) continue; | |
211 | icharge = partMC->Charge() > 0 ? 0 : 1; | |
212 | gCharge[icharge]++; | |
213 | gPid[a][icharge]++; | |
214 | } | |
215 | } | |
216 | else { | |
217 | printf(" No Event Type is Defined "); | |
218 | return; | |
219 | } | |
220 | ||
221 | if( (gCharge[0] + gCharge[1]) != 0 ) { | |
222 | if (isQA) { | |
223 | fEventCounter->Fill(6); | |
224 | fEventCounter->Fill(160 + gCent); | |
225 | } | |
226 | else { | |
227 | Double_t vsparse[kNSparseData]; | |
228 | vsparse[0] = gCent; | |
229 | vsparse[1] = gRefMul; | |
230 | vsparse[2] = gCharge[0] + gCharge[1]; | |
231 | vsparse[3] = gCharge[0]; | |
232 | vsparse[4] = gCharge[1]; | |
233 | vsparse[5] = gPid[0][0] + gPid[0][1]; | |
234 | vsparse[6] = gPid[0][0]; | |
235 | vsparse[7] = gPid[0][1]; | |
236 | vsparse[8] = gPid[1][0] + gPid[1][0]; | |
237 | vsparse[9] = gPid[1][0]; | |
238 | vsparse[10] = gPid[1][1]; | |
239 | vsparse[11] = gPid[2][0] + gPid[2][1]; | |
240 | vsparse[12] = gPid[2][0]; | |
241 | vsparse[13] = gPid[2][1]; | |
242 | fHistoCorrelation->Fill(vsparse); | |
243 | } | |
244 | } | |
245 | ||
246 | if(fDebug && isQA) Printf(" %6d %6d %6d %6d %6d %6d %6d %6d %6d %6d %6d %6d %6d", | |
247 | (Int_t)fEventCounter->GetBinContent(1), | |
248 | (Int_t)fEventCounter->GetBinContent(2), | |
249 | (Int_t)fEventCounter->GetBinContent(3), | |
250 | (Int_t)gCent, (Int_t)gRefMul, | |
251 | gCharge[0], gCharge[1], | |
252 | gPid[0][0], gPid[0][1], gPid[1][0], | |
253 | gPid[1][1], gPid[2][0], gPid[2][1]); | |
254 | ||
255 | PostData(1, fThnList); | |
256 | } | |
257 | ||
258 | void AliEbyEParticleRatioFluctuationTask::Terminate( Option_t * ){ | |
259 | Info("AliEbyEParticleRatioFluctuationTask"," Task Successfully finished"); | |
260 | } | |
261 | ||
262 | //___________________________________________________________ | |
263 | Bool_t AliEbyEParticleRatioFluctuationTask::AcceptEvent(AliAODEvent *event, Int_t cent) const { | |
264 | Bool_t ver = kFALSE; | |
265 | const AliAODVertex *vertex = event->GetPrimaryVertex(); | |
266 | if(vertex) { | |
267 | Double32_t fCov[6]; | |
268 | vertex->GetCovarianceMatrix(fCov); | |
269 | if(vertex->GetNContributors() > 0) { | |
270 | if(fCov[5] != 0) { | |
271 | ||
272 | if(isQA) { | |
273 | fEventCounter->Fill(3); | |
274 | fHistQA[3]->Fill(cent,vertex->GetX()); | |
275 | fHistQA[4]->Fill(cent,vertex->GetY()); | |
276 | fHistQA[5]->Fill(cent,vertex->GetZ()); | |
277 | } | |
278 | ||
279 | if(TMath::Abs(vertex->GetX()) < fVxMax) { | |
280 | if(TMath::Abs(vertex->GetY()) < fVyMax) { | |
281 | if(TMath::Abs(vertex->GetZ()) < fVzMax) { | |
282 | ver = kTRUE; | |
283 | if(isQA) { | |
284 | fEventCounter->Fill(4); | |
285 | fHistQA[0]->Fill(cent,vertex->GetX()); | |
286 | fHistQA[1]->Fill(cent,vertex->GetY()); | |
287 | fHistQA[2]->Fill(cent,vertex->GetZ()); | |
288 | } | |
289 | } | |
290 | } | |
291 | } | |
292 | } | |
293 | } | |
294 | } | |
295 | ||
296 | AliCentrality *centrality = event->GetCentrality(); | |
297 | if (centrality->GetQuality() != 0) ver = kFALSE; | |
298 | return ver; | |
299 | } | |
300 | ||
301 | ||
302 | //___________________________________________________________ | |
303 | Bool_t AliEbyEParticleRatioFluctuationTask::AcceptTrack(AliAODTrack *track, Int_t cent) const { | |
304 | if(!track) return kFALSE; | |
305 | if (track->Charge() == 0 ) return kFALSE; | |
306 | ||
307 | if(isQA) { | |
308 | fHistQA[6]->Fill(cent,track->DCA()); | |
309 | fHistQA[7]->Fill(cent,track->ZAtDCA()); | |
310 | fHistQA[8]->Fill(cent,track->Pt()); | |
311 | fHistQA[9]->Fill(cent,track->Eta()); | |
312 | } | |
313 | ||
314 | if (!track->TestFilterBit(fAODtrackCutBit)) return kFALSE; | |
315 | ||
316 | if (track->Eta() < fEtaLowerLimit || | |
317 | track->Eta() > fEtaHigherLimit) return kFALSE; | |
318 | if (track->Pt() < fPtLowerLimit || | |
319 | track->Pt() > fPtHigherLimit) return kFALSE; | |
320 | if ( track->DCA() > fDCAxy ) return kFALSE; | |
321 | if ( track->ZAtDCA() > fDCAz ) return kFALSE; | |
322 | ||
323 | if(isQA) { | |
324 | fHistQA[10]->Fill(cent,track->DCA()); | |
325 | fHistQA[11]->Fill(cent,track->ZAtDCA()); | |
326 | fHistQA[12]->Fill(cent,track->Pt()); | |
327 | fHistQA[13]->Fill(cent,track->Eta()); | |
328 | } | |
329 | ||
330 | return kTRUE; | |
331 | } | |
332 | ||
333 | ||
334 | //___________________________________________________________ | |
335 | Bool_t AliEbyEParticleRatioFluctuationTask::AcceptMCTrack(AliAODMCParticle *track, Int_t cent) const { | |
336 | if(!track) return kFALSE; | |
337 | if(!track->IsPhysicalPrimary()) return kFALSE; | |
338 | if (track->Charge() == 0 ) return kFALSE; | |
339 | if(isQA) { | |
340 | fHistQA[8]->Fill(cent,track->Pt()); | |
341 | fHistQA[9]->Fill(cent,track->Eta()); | |
342 | } | |
343 | ||
344 | if (track->Eta() < fEtaLowerLimit || | |
345 | track->Eta() > fEtaHigherLimit) return kFALSE; | |
346 | if (track->Pt() < fPtLowerLimit || | |
347 | track->Pt() > fPtHigherLimit) return kFALSE; | |
348 | ||
349 | if(isQA) { | |
350 | fHistQA[12]->Fill(cent,track->Pt()); | |
351 | fHistQA[13]->Fill(cent,track->Eta()); | |
352 | } | |
353 | ||
354 | return kTRUE; | |
355 | } |