AliAODEvent::GetHeader() returns AliVHeader
[u/mrichter/AliRoot.git] / PWGCF / EBYE / PIDFluctuation / task / AliEbyEParticleRatioFluctuationTask.cxx
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 = dynamic_cast<AliAODHeader*>(event->GetHeader());
147   if(!aodHeader) AliFatal("Not a standard AOD");
148   gCent = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
149   gRefMul = aodHeader->GetRefMultiplicity();
150   if (gCent < 0 || gCent > 100) return;
151   if (isQA) fEventCounter->Fill(2);  
152
153   if (!AcceptEvent(event,gCent)) return;
154   
155   Int_t icharge = -1;
156   Int_t gCharge[2];
157   Int_t gPid[3][2];
158   
159   for (Int_t i = 0; i < 2; i++) {
160     gCharge[i] = 0;
161     for (Int_t ii = 0; ii < 3; ii++) {
162       gPid[ii][i] = 0;
163     }
164   }
165
166
167   if(fAnalysisType == "AOD") {
168     if (isQA) {
169         fEventCounter->Fill(5);
170         fEventCounter->Fill(50+gCent);
171     }
172     
173     for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) {
174       AliAODTrack* track = dynamic_cast<AliAODTrack *>(event->GetTrack(itrk));
175       if (!track) continue;
176       
177       if (!AcceptTrack(track, gCent)) continue;
178             
179       Int_t a = fHelperPID->GetParticleSpecies((AliVTrack*) track,kTRUE);
180
181       if(a < 0 || a > 2) continue;
182       icharge = track->Charge() > 0 ? 0 : 1;
183       gCharge[icharge]++;
184       gPid[a][icharge]++;
185     }
186   }
187   else if(fAnalysisType == "MCAOD") {
188     TClonesArray *arrayMC= 0; 
189     arrayMC = dynamic_cast<TClonesArray*> (event->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
190     if (!arrayMC) {
191       Printf("Error: MC particles branch not found!\n");
192       return;
193     }
194     AliAODMCHeader *mcHdr=0;
195     mcHdr=(AliAODMCHeader*)event->GetList()->FindObject(AliAODMCHeader::StdBranchName());  
196     if(!mcHdr) {
197       Printf("MC header branch not found!\n");
198       return;
199     }
200     
201     if (isQA) {
202       fEventCounter->Fill(5);
203       fEventCounter->Fill(50+gCent);
204     }
205     
206     Int_t nMC = arrayMC->GetEntries();
207     for (Int_t iMC = 0; iMC < nMC; iMC++) {
208       AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
209       if(!AcceptMCTrack(partMC, gCent)) continue;
210       Int_t a = fHelperPID->GetMCParticleSpecie((AliVEvent*) event,(AliVTrack*)partMC,1);
211       if(a < 0 || a > 2) continue;
212       icharge = partMC->Charge() > 0 ? 0 : 1;
213       gCharge[icharge]++;
214       gPid[a][icharge]++;
215     }
216   }
217   else {
218     printf(" No Event Type is Defined ");
219     return;
220   }
221   
222   if( (gCharge[0] + gCharge[1]) != 0 ) {
223     if (isQA) {
224       fEventCounter->Fill(6); 
225       fEventCounter->Fill(160 + gCent);
226     }
227     else { 
228       Double_t vsparse[kNSparseData];
229       vsparse[0]   = gCent;
230       vsparse[1]   = gRefMul;
231       vsparse[2]   = gCharge[0] + gCharge[1];
232       vsparse[3]   = gCharge[0];
233       vsparse[4]   = gCharge[1];
234       vsparse[5]   = gPid[0][0] + gPid[0][1];
235       vsparse[6]   = gPid[0][0];
236       vsparse[7]   = gPid[0][1];
237       vsparse[8]   = gPid[1][0] + gPid[1][0];
238       vsparse[9]   = gPid[1][0];
239       vsparse[10]  = gPid[1][1];
240       vsparse[11]  = gPid[2][0] + gPid[2][1];
241       vsparse[12]  = gPid[2][0];
242       vsparse[13]  = gPid[2][1];
243       fHistoCorrelation->Fill(vsparse);
244     }
245   }
246   
247   if(fDebug && isQA)  Printf(" %6d  %6d %6d  %6d %6d  %6d %6d %6d %6d %6d %6d %6d %6d", 
248                              (Int_t)fEventCounter->GetBinContent(1), 
249                              (Int_t)fEventCounter->GetBinContent(2), 
250                              (Int_t)fEventCounter->GetBinContent(3), 
251                              (Int_t)gCent, (Int_t)gRefMul, 
252                              gCharge[0], gCharge[1], 
253                              gPid[0][0], gPid[0][1],  gPid[1][0], 
254                              gPid[1][1],  gPid[2][0], gPid[2][1]);
255   
256   PostData(1, fThnList);
257 }
258
259 void AliEbyEParticleRatioFluctuationTask::Terminate( Option_t * ){
260   Info("AliEbyEParticleRatioFluctuationTask"," Task Successfully finished");
261 }
262
263 //___________________________________________________________
264 Bool_t AliEbyEParticleRatioFluctuationTask::AcceptEvent(AliAODEvent *event, Int_t cent) const {
265   Bool_t ver = kFALSE;
266   const AliAODVertex *vertex = event->GetPrimaryVertex();
267   if(vertex) {
268     Double32_t fCov[6];
269     vertex->GetCovarianceMatrix(fCov);
270     if(vertex->GetNContributors() > 0) {
271       if(fCov[5] != 0) {
272         
273         if(isQA) {      
274           fEventCounter->Fill(3);
275           fHistQA[3]->Fill(cent,vertex->GetX());
276           fHistQA[4]->Fill(cent,vertex->GetY());
277           fHistQA[5]->Fill(cent,vertex->GetZ());
278         }
279         
280         if(TMath::Abs(vertex->GetX()) < fVxMax) {
281           if(TMath::Abs(vertex->GetY()) < fVyMax) {
282             if(TMath::Abs(vertex->GetZ()) < fVzMax) {
283               ver = kTRUE;
284               if(isQA) {        
285                 fEventCounter->Fill(4);
286                 fHistQA[0]->Fill(cent,vertex->GetX());
287                 fHistQA[1]->Fill(cent,vertex->GetY());
288                 fHistQA[2]->Fill(cent,vertex->GetZ());
289               }
290             }
291           }
292         }
293       }
294     }
295   }
296   
297   AliCentrality *centrality = event->GetCentrality();
298   if (centrality->GetQuality() != 0) ver = kFALSE;
299   return ver;
300 }
301
302
303 //___________________________________________________________
304 Bool_t AliEbyEParticleRatioFluctuationTask::AcceptTrack(AliAODTrack *track, Int_t cent) const {
305   if(!track)                                  return kFALSE;
306   if (track->Charge() == 0 )                  return kFALSE;
307   
308   if(isQA) {
309     fHistQA[6]->Fill(cent,track->DCA());
310     fHistQA[7]->Fill(cent,track->ZAtDCA());
311     fHistQA[8]->Fill(cent,track->Pt());
312     fHistQA[9]->Fill(cent,track->Eta());
313   }
314   
315   if (!track->TestFilterBit(fAODtrackCutBit)) return kFALSE;
316
317   if (track->Eta() < fEtaLowerLimit ||
318       track->Eta() > fEtaHigherLimit)         return kFALSE;
319   if (track->Pt() < fPtLowerLimit ||
320       track->Pt() > fPtHigherLimit)           return kFALSE;  
321   if ( track->DCA() > fDCAxy )                return kFALSE; 
322   if ( track->ZAtDCA() > fDCAz )              return kFALSE;
323    
324   if(isQA) {
325     fHistQA[10]->Fill(cent,track->DCA());
326     fHistQA[11]->Fill(cent,track->ZAtDCA());
327     fHistQA[12]->Fill(cent,track->Pt());
328     fHistQA[13]->Fill(cent,track->Eta());
329   }
330  
331   return kTRUE;
332 }
333
334
335 //___________________________________________________________
336 Bool_t AliEbyEParticleRatioFluctuationTask::AcceptMCTrack(AliAODMCParticle *track, Int_t cent) const {
337   if(!track) return kFALSE;
338   if(!track->IsPhysicalPrimary()) return kFALSE;
339   if (track->Charge() == 0 ) return kFALSE;
340   if(isQA) {
341     fHistQA[8]->Fill(cent,track->Pt());
342     fHistQA[9]->Fill(cent,track->Eta());
343   }
344
345   if (track->Eta() < fEtaLowerLimit ||
346       track->Eta() > fEtaHigherLimit) return kFALSE;
347   if (track->Pt() < fPtLowerLimit ||
348       track->Pt() > fPtHigherLimit) return kFALSE;  
349     
350   if(isQA) {
351     fHistQA[12]->Fill(cent,track->Pt());
352     fHistQA[13]->Fill(cent,track->Eta());
353   }
354  
355   return kTRUE;
356 }