]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/Fluctuations/AliEbyEMultFluctuationTask.cxx
1) Updated AddTask for efficiency with different filterbit notation, 2) Added an...
[u/mrichter/AliRoot.git] / PWGCF / EBYE / Fluctuations / AliEbyEMultFluctuationTask.cxx
CommitLineData
e70360ad 1
2#include "TChain.h"
3#include "TTree.h"
4#include "TNtuple.h"
5#include "TH1F.h"
6#include "TH2D.h"
7#include "TH1D.h"
8#include "TCanvas.h"
9#include "TArrayD.h"
10#include "AliAnalysisTask.h"
11#include "AliAnalysisManager.h"
12#include "AliAODHeader.h"
13#include "AliAODTrack.h"
14#include "AliAODTracklets.h"
15#include "AliAODpidUtil.h"
16#include "AliVTrack.h"
17#include "AliAODEvent.h"
18#include "AliAODInputHandler.h"
19#include "AliCentrality.h"
20#include "AliEbyEMultFluctuationTask.h"
21#include "AliAODVertex.h"
22#include "AliAODVZERO.h"
23#include "AliCentrality.h"
24#include "AliPID.h"
25#include "AliPIDResponse.h"
26#include "AliPIDCombined.h"
27#include "AliAODpidUtil.h"
28#include "TRandom.h"
29#include "AliTPCPIDResponse.h"
30#include "AliTRDPIDResponse.h"
31
32
33ClassImp(AliEbyEMultFluctuationTask)
34
35//________________________________________________________________________
36AliEbyEMultFluctuationTask::AliEbyEMultFluctuationTask(const char *name)
37: AliAnalysisTaskSE(name), fAOD(0), fOutputList(0),fHistNchPt(0),fHistNchEta(0),fHistNchEtaCent(0),fHistNchPhi(0),fHistDCAxy(0),fHistDCAz(0),fHistnclus(0),fHistchi2ndf(0),fDeDx(NULL),fHistMult(0),fAODVertex(0),fPID(0),fEventCounter(0),histcounter(0),fHistVz(0),fHistchi2ndfvscs(0),My_ntuple(0),fHistMultV0A(0),fHistMultV0C(0),fHistMultV0total(0),fCentralityCounter(0),fCentralityEstimator("V0M"),Counter(0)
38
39
40
41{
42 // Constructor
43 for(Int_t ibin=0;ibin<91;ibin++)
44 {
45 fMult[ibin]=NULL;
46 }
47 for(Int_t jbin=0;jbin<46;jbin++)
48 {
49 fMultTwo[jbin]=NULL;
50 }
51 for(Int_t kbin=0;kbin<15;kbin++)
52 {
53 fMultFive[kbin]=NULL;
54 }
55
56
57 DefineInput(0, TChain::Class());
58
59 DefineOutput(1, TList::Class());
60}
61
62//________________________________________________________________________
63Bool_t AliEbyEMultFluctuationTask :: SelectEvent(AliAODVertex* vertex)
64{
65if(vertex)
66
67 if(vertex->GetNContributors() < 0) return kFALSE;
68
69
70 Double_t lvx = vertex->GetX();
71 Double_t lvy = vertex->GetY();
72 Double_t lvz = vertex->GetZ();
73
74 fEventCounter->Fill(3);
75
76
77 if(TMath::Abs(lvx) > 0.3) return kFALSE;
78 if(TMath::Abs(lvy) > 0.3) return kFALSE;
79 if(TMath::Abs(lvz) > 10) return kFALSE;
80 if(vertex->GetType()==AliAODVertex::kPileupSPD) return kFALSE;
81
82 fEventCounter->Fill(5);
83 fHistVz->Fill(lvz);
84return kTRUE;
85}
86//_____________________________________________________________________________
87Int_t AliEbyEMultFluctuationTask :: SelectTrack(AliAODTrack* track)
88{
89
90 Double_t eta = track->Eta();
91 Double_t phi = track->Phi();
92
93 if(TMath :: Abs(eta)>0.8) return 0;
94
95 if(!track->TestFilterBit(128)) return 0;
96
97 Float_t dxy, dz;
98
99 dxy = track->DCA();
100 dz = track->ZAtDCA();
101 if(TMath ::Abs(dxy) >2.4 || TMath ::Abs(dz)>3.0) return 0;
102 // cout<<dxy<<dz<<endl;
103
104
105 Double_t nclus = track->GetTPCClusterInfo(2,1);
106 if(nclus<80) return 0;
107
108//chi2cut
109
110 Double_t chi2ndf = track->Chi2perNDF();
111 if(chi2ndf>4) return 0;
112
113return 1;
114
115
116}
117
118//___________________________________________________________________________________
119void AliEbyEMultFluctuationTask::UserCreateOutputObjects()
120{
121 fOutputList = new TList();
122
123
124
125 My_ntuple = new TNtuple("My_ntuple", "My_ntuple", "Mult:Mult1:Mult2:nCentrality:fV0total:nBin2:nBin5:spdmult0:spdmult1:tracklets:run");
126 fOutputList->Add(My_ntuple);
127
128fOutputList->SetOwner(kTRUE);
129fEventCounter = new TH1D("fEventCounter","EventCounter", 10, 0.5,10.5);
130fOutputList->Add(fEventCounter);
131
132
133histcounter = new TH1D("histcounter","histcounter", 10, 0.5,10.5);
134fOutputList->Add(histcounter);
135
136
137TString BinLebel;
138fCentralityCounter = new TH1D("fEventCounter","EventCounter", 100, 0.5,100.5);
139TAxis *xaxis = fCentralityCounter->GetXaxis();
140 for(Int_t ibin=0;ibin<100;ibin++)
141{
142BinLebel="cent";BinLebel+=ibin;
143xaxis->SetBinLabel(ibin,BinLebel.Data());
144
145}
146fOutputList->Add(fCentralityCounter);
147
148fHistMultV0A = new TH1F("fHistMultV0A","Mult-VOA",22000,0,22000);
149fHistMultV0C = new TH1F("fHistMultV0C","Mult-VOC",22000,0,22000);
150fHistMultV0total = new TH1F("fHistMultV0total","Mult-VOTotal",22000,0,22000);
151
152
153 fHistNchEta = new TH1D("fHistNchEta", "#eta distribution of charged tracks(for hybrid)", 200, -0.8, 0.8);
154fHistNchEta->GetXaxis()->SetTitle("#eta");
155fHistNchEta->GetYaxis()->SetTitle("dNch/d#eta");
156
157
158fHistNchEtaCent = new TH1D("fHistNchEtaCent", "Central #eta distribution of charged tracks(for hybrid)", 200, -0.8, 0.8);
159fHistNchEtaCent->GetXaxis()->SetTitle("#eta");
160fHistNchEtaCent->GetYaxis()->SetTitle("dNch/d#eta(Central)");
161
162
163
164fHistNchPhi = new TH1D("fHistNchPhi", "#phi distribution of charged tracks(for hybrid)", 200, -10, 10);
165fHistNchPhi->GetXaxis()->SetTitle("#phi");
166fHistNchPhi->GetYaxis()->SetTitle("dNch/d#phi");
167
168fHistNchPt = new TH1D("fHistNchPt", "P_{T} distribution(charged tracks for hybrid)", 150, 0.1, 3.1);
169fHistNchPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
170fHistNchPt->GetYaxis()->SetTitle("dNch/dP_{T}(c/GeV)");
171
172
173 fHistchi2ndfvscs = new TH2D("fHistchi2ndfvscs","Chi2/ndf from Gaussian fits of Multiplicity Distribution",95,0,95,60,0,6);
174 fHistchi2ndfvscs->SetXTitle("%cs");
175 fHistchi2ndfvscs->SetYTitle("Chi2/ndf");
176
177
178fHistnclus=new TH1D("no of clusters","no of tracks after TPCcluster cut",200,0,200);
179fHistchi2ndf= new TH1D("chi_square/ndf","no of tracks after chi_square/ndf cut",10,0,10);
180
181 fHistVz = new TH1D("fHistVz","Primary vertex distribution - z coordinate",90,-15,15);
182 fHistVz->SetXTitle("Vz");
183 fHistVz->SetYTitle("Entries");
184
185
186 fOutputList->Add(fHistMultV0A);
187 fOutputList->Add(fHistMultV0C);
188 fOutputList->Add(fHistMultV0total);
189 fOutputList->Add(fHistnclus);
190 fOutputList->Add(fHistchi2ndf);
191 fOutputList->Add(fHistNchPt);
192 fOutputList->Add(fHistNchEta);
193 fOutputList->Add(fHistNchEtaCent);
194
195 fOutputList->Add(fHistNchPhi);
196 fOutputList->Add(fHistVz);
197 fOutputList->Add(fHistchi2ndfvscs);
198
199
200
201TString hname;
202TString htitle;
203
204 for(Int_t i = 0; i < 25; i++) {
205 hname = "fMult"; hname+=i;
206 htitle = "Multiplicity in Cent Bin "; htitle+=i;
207 fMult[i] = new TH1F(hname.Data(),htitle.Data(),5000,0.5,5000.5);
208 fOutputList->Add(fMult[i]);
209 }
210
211 for(Int_t i = 25; i < 50; i++) {
212 hname = "fMult"; hname+=i;
213 htitle = "Multiplicity in Cent Bin "; htitle+=i;
214 fMult[i] = new TH1F(hname.Data(),htitle.Data(),3000,0.5,3000.5);
215 fOutputList->Add(fMult[i]);
216 }
217
218 for(Int_t i = 50; i < 91; i++) {
219 hname = "fMult"; hname+=i;
220 htitle = "Multiplicity in Cent Bin "; htitle+=i;
221 fMult[i] = new TH1F(hname.Data(),htitle.Data(),2000,0.5,2000.5);
222 fOutputList->Add(fMult[i]);
223 }
224
225//****************for 2% centrality bin**************************//
226for(Int_t i = 0; i < 12; i++) {
227 hname = "fMultTwo"; hname+=i;
228 htitle = "Multiplicity in Cent Bin "; htitle+=i;
229 fMultTwo[i] = new TH1F(hname.Data(),htitle.Data(),5000,0.5,5000.5);
230 fOutputList->Add(fMultTwo[i]);
231 }
232
233 for(Int_t i = 12; i < 25; i++) {
234 hname = "fMultTwo"; hname+=i;
235 htitle = "Multiplicity in Cent Bin "; htitle+=i;
236 fMultTwo[i] = new TH1F(hname.Data(),htitle.Data(),3000,0.5,3000.5);
237 fOutputList->Add(fMultTwo[i]);
238 }
239
240 for(Int_t i = 25; i <= 45; i++) {
241 hname = "fMultTwo"; hname+=i;
242 htitle = "Multiplicity in Cent Bin "; htitle+=i;
243 fMultTwo[i] = new TH1F(hname.Data(),htitle.Data(),2000,0.5,2000.5);
244 fOutputList->Add(fMultTwo[i]);
245 }
246
247
248
249//*****************************for 5 % centrality bin********************//
250for(Int_t i = 0; i < 5; i++) {
251 hname = "fMultFive"; hname+=i;
252 htitle = "Multiplicity in Cent Bin "; htitle+=i;
253 fMultFive[i] = new TH1F(hname.Data(),htitle.Data(),5000,0.5,5000.5);
254 fOutputList->Add(fMultFive[i]);
255 }
256
257 for(Int_t i = 5; i < 10; i++) {
258 hname = "fMultFive"; hname+=i;
259 htitle = "Multiplicity in Cent Bin "; htitle+=i;
260 fMultFive[i] = new TH1F(hname.Data(),htitle.Data(),3000,0.5,3000.5);
261 fOutputList->Add(fMultFive[i]);
262 }
263
264 for(Int_t i = 10; i <= 14; i++) {
265 hname = "fMultFive"; hname+=i;
266 htitle = "Multiplicity in Cent Bin "; htitle+=i;
267 fMultFive[i] = new TH1F(hname.Data(),htitle.Data(),2000,0.5,2000.5);
268 fOutputList->Add(fMultFive[i]);
269 }
270
271
272 }
273
274//________________________________________________________________________
275void AliEbyEMultFluctuationTask::UserExec(Option_t *)
276{
277
278
279
280// Post output data.
281 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
282 if (!fAOD) {
283 printf("ERROR: fAOD not available\n");
284 return;
285
286 }
287
288 // Int_t Counter = 0;
289Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
290if(!isSelected) return;
291 fEventCounter->Fill(1);
292 //Counter++;
293
294
295// Get the event vertex.
296 fAODVertex = fAOD->GetPrimaryVertex();
297
298 if(TMath::Abs((fAODVertex->GetZ())-((fAOD->GetPrimaryVertexSPD())->GetZ())) > 0.5) return;
299
300 if (!fAODVertex) {
301
302 return;
303 }
304
305
306 if (!SelectEvent(fAODVertex)) return;
307
308 //***************VZERO-AMPLITUDE*************************//
309
310 AliAODVZERO *aodV0 = fAOD->GetVZEROData();
311 Float_t fV0Amult = aodV0->GetMTotV0A();
312 Float_t fV0Cmult = aodV0->GetMTotV0C();
313 Float_t fV0total = fV0Amult + fV0Cmult;
314
315 //*****************SPD-CLUSTERS*************************//
316
317
318
319 AliAODHeader *fHeader = fAOD->GetHeader();
320 Int_t spdmult0 = fHeader->GetNumberOfITSClusters(0);
321 Int_t spdmult1 = fHeader->GetNumberOfITSClusters(1);
322 Int_t run = fHeader->GetRunNumber();
323
324 AliAODTracklets *fTracklets = fAOD->GetTracklets();
325 Int_t tracklets = fTracklets->GetNumberOfTracklets();
326
327
328 //********************CENTRALITY******************************************//
329
330 AliCentrality *centrality= fAOD->GetCentrality();
331 Double_t cent=centrality->GetCentralityPercentile("V0M");
332 // cout<<"cent= "<<cent<<endl;
333
334 if(cent < 0 || cent >= 91) return;
335
336Int_t nCentrality = (Int_t)cent;
337
338Int_t nBin2= (Int_t)( nCentrality/2.0);
339
340Int_t nBin5 = (Int_t)( nCentrality/5.0);
341 if(nBin5==0){
342
343histcounter->Fill(1);
344 }
345
346
347
348fCentralityCounter->Fill(nCentrality);
349
350fEventCounter->Fill(7);
351
352
353Double_t Mult=0.0;
354 Double_t Mult1=0.0;
355 Double_t Mult2=0.0;
356
357 //printf("There are %d tracks in this event\n", fAOD->GetNumberOfTracks());
358
359
360 for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
361 AliAODTrack* track = fAOD->GetTrack(iTracks);
362 if (!track) {
363 printf("ERROR: Could not receive track %d\n", iTracks);
364 continue;
365 }
366
367
368
369
370// Find the track classification.
371Int_t tracktype = SelectTrack(track);
372
373
374 if (tracktype==0) {
375 continue;
376 }
377
378 if(track->GetType()== AliAODTrack::kPrimary)
379{
380 Short_t Charge = track->Charge();
381 Double_t eta = track->Eta();
382 if(TMath :: Abs(Charge)==1 && nBin5==0){
383 fHistNchEtaCent->Fill(eta);
384
385 }
386
387 Double_t pt = track->Pt();
388 if( pt < 0.3 || pt > 2.0 ) continue;
389 Double_t nclus = track->GetTPCClusterInfo(2,1);
390 Double_t chi2ndf = track->Chi2perNDF();
391 fHistnclus->Fill(nclus);
392 fHistchi2ndf->Fill(chi2ndf);
393 fHistchi2ndfvscs->Fill(cent,chi2ndf);
394
395 //Count the charge particles-----
396 if(TMath :: Abs(Charge)==1) Mult +=1;
397 if(Charge==1) Mult1 +=1;
398 if(Charge==-1) Mult2 +=1;
399
400
401 Double_t phi = track->Phi();
402 if(TMath :: Abs(Charge)==1) fHistNchPt->Fill(pt);
403
404
405 if(TMath :: Abs(Charge)==1) fHistNchEta->Fill(eta);
406 if(TMath :: Abs(Charge)==1) fHistNchPhi->Fill(phi);
407
408
409 }
410
411
412 } //track loop
413
414
415 //****************Filling Tree*************************//
416
417 My_ntuple->Fill(Mult,Mult1,Mult2,nCentrality,fV0total,nBin2,nBin5,spdmult0,spdmult1,tracklets,run);
418
419
420fMult[nCentrality]->Fill(Mult);
421
422if(nBin2<=45) {
423fMultTwo[nBin2]->Fill(Mult);
424 }
425if(nBin5<=14){
426fMultFive[nBin5]->Fill(Mult);
427 }
428
429 fHistMultV0A->Fill(fV0Amult);
430 fHistMultV0C->Fill(fV0Cmult);
431 fHistMultV0total->Fill(fV0total);
432
433 PostData(1, fOutputList);
434}
435//________________________________________________________________________
436void AliEbyEMultFluctuationTask::Terminate(Option_t *)
437{
438
439
440 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
441 if (!fOutputList) {
442 printf("ERROR: Output list not available\n");
443 return;
444 }
445
446}