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