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