AliAODEvent::GetHeader() returns AliVHeader
[u/mrichter/AliRoot.git] / PWGCF / EBYE / Fluctuations / AliEbyEHigherMomentsTask.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 //                                                                         //
19 //           Analysis Task for Net-Charge Higher Moment Analysis           //
20 //              Author: Satyajit Jena || Nirbhay K. Behera                 //
21 //                      sjena@cern.ch || nbehera@cern.ch                   //
22 //           Last Modified: 30/01/2014: only for net-charge part           //
23 //=========================================================================//
24 #include "TChain.h"
25 #include "TList.h"
26 #include "TFile.h"
27 #include "TTree.h"
28 #include "TH1D.h"
29 #include "TH2D.h"
30 #include "TH3D.h"
31 #include "THnSparse.h"
32 #include "TCanvas.h"
33
34 #include "AliAnalysisTaskSE.h"
35 #include "AliAnalysisManager.h"
36 #include "AliAODEvent.h"
37 #include "AliVTrack.h"
38 #include "AliAODTrack.h"
39 #include "AliAODInputHandler.h"
40 #include "AliAODEvent.h"
41 #include "AliStack.h"
42 #include "AliAODMCHeader.h"
43 #include "AliAODMCParticle.h"
44 #include "AliMCEventHandler.h"
45 #include "AliMCEvent.h"
46 #include "AliStack.h"
47 #include "AliGenHijingEventHeader.h"
48 #include "AliGenEventHeader.h"
49 #include "AliPID.h"
50 #include "AliAODPid.h"                
51 #include "AliPIDResponse.h"
52 #include "AliAODpidUtil.h"        
53 #include "AliPIDCombined.h"
54 #include "AliHelperPID.h"
55
56 #include "AliEbyEHigherMomentsTask.h"
57
58 using std::cout;
59 using std::endl;
60
61 ClassImp(AliEbyEHigherMomentsTask)
62
63
64 //-----------------------------------------------------------------------
65 AliEbyEHigherMomentsTask::AliEbyEHigherMomentsTask( const char *name )
66 : AliAnalysisTaskSE( name ),
67   fListOfHistos(0),
68   fArrayMC(0),
69   fAnalysisType(0),
70   fCentralityEstimator("V0M"),
71   fCentrality(0),
72   fVxMax(3.),
73   fVyMax(3.),
74   fVzMax(10.),
75   fPtLowerLimit(0.3),
76   fPtHigherLimit(1.5),
77   fNptBins(0),
78   fBin(0),
79   fEtaLowerLimit(-0.8),
80   fEtaHigherLimit(0.8),
81   fAODtrackCutBit(128),
82   fEventCounter(0),
83   fHistVxVy(0),
84   fTHnCentNplusNminusCh(0),
85   fTHnCentNplusNminusChTruth(0),
86   fPtBinNplusNminusCh(0),
87   fPtBinNplusNminusChTruth(0)
88
89   
90   for ( Int_t i = 0; i < 4; i++) { 
91     fHistQA[i] = NULL;
92   }
93  
94   DefineOutput(1, TList::Class()); // Outpput....
95   //DefineOutput(2, TList::Class()); 
96 }
97
98 AliEbyEHigherMomentsTask::~AliEbyEHigherMomentsTask() {
99
100   if(fListOfHistos) delete fListOfHistos;
101 }
102
103 //---------------------------------------------------------------------------------
104 void AliEbyEHigherMomentsTask::UserCreateOutputObjects() {
105   
106   fListOfHistos = new TList();
107   fListOfHistos->SetOwner(kTRUE);
108   
109   fEventCounter = new TH1D("fEventCounter","EventCounter", 10, 0.5,10.5);
110   fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
111   fEventCounter->GetXaxis()->SetBinLabel(2,"Within 0-80% centrality");
112   fEventCounter->GetXaxis()->SetBinLabel(5,"Have a vertex");
113   fEventCounter->GetXaxis()->SetBinLabel(6,"After vertex Cut");
114   fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
115   fListOfHistos->Add(fEventCounter);
116   
117   //For QA-Histograms
118   fHistQA[0] = new TH1D("fHistQAvz", "Histo Vz After Cut", 500, -25., 25.);  
119   fHistQA[1] = new TH1D("fHistQAPt","p_{T} distribution",600,0,6);
120   fHistQA[2] = new TH1D("fHistQAEta","#eta distribution",240,-1.2,1.2);
121   fHistQA[3] = new TH1D("fHistQAPhi","#phi distribution",340,0,6.8);
122  
123   for(Int_t i = 0; i < 4; i++)
124     {
125       fListOfHistos->Add(fHistQA[i]);
126     }
127   
128   fHistVxVy = new TH2D("fHistVxVy","Vertex-x Vs Vertex-y", 200, -1., 1., 200, -1., 1.);
129   fListOfHistos->Add(fHistVxVy);
130  
131  
132   const Int_t nDim = 3;
133   
134   Int_t fBinsCh[nDim] = {100, 1500, 1500};
135   Double_t fMinCh[nDim] = { -0.5, -0.5, -0.5 };
136   Double_t fMaxCh[nDim] = { 99.5, 1499.5, 1499.5};
137   
138   const Int_t dim = fNptBins*2 + 1;
139   //const Int_t dim = ;
140   Int_t bin[dim];
141   Double_t min[dim];
142   Double_t max[dim];
143   bin[0] = 100; min[0] = -0.5; max[0] = 99.5;
144
145   for(Int_t i = 1; i < dim; i++){
146     
147     bin[i] = 900;
148     min[i] = -0.5;
149     max[i] = 899.5;
150  
151   } 
152   
153   fTHnCentNplusNminusCh = new THnSparseD("fTHnCentNplusNminusCh","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh); 
154   fTHnCentNplusNminusCh->GetAxis(0)->SetTitle("Centrality");
155   fTHnCentNplusNminusCh->GetAxis(1)->SetTitle("Nplus");
156   fTHnCentNplusNminusCh->GetAxis(2)->SetTitle("Nminus");
157   fListOfHistos->Add(fTHnCentNplusNminusCh);
158
159   fPtBinNplusNminusCh = new THnSparseI("fPtBinNplusNminusCh","cent-nplus-nminus", dim, bin, min, max);
160   
161   fListOfHistos->Add(fPtBinNplusNminusCh);
162   
163   if( fAnalysisType == "MCAOD" ){   
164     
165     fTHnCentNplusNminusChTruth = new THnSparseD("fTHnCentNplusNminusChTruth","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh); 
166     fTHnCentNplusNminusChTruth->GetAxis(0)->SetTitle("Centrality");
167     fTHnCentNplusNminusChTruth->GetAxis(1)->SetTitle("Nplus");
168     fTHnCentNplusNminusChTruth->GetAxis(2)->SetTitle("Nminus");
169     fListOfHistos->Add(fTHnCentNplusNminusChTruth);
170     
171     fPtBinNplusNminusChTruth = new THnSparseI("fPtBinNplusNminusChTruth","cent-nplus-nminus", dim, bin, min, max);
172     fListOfHistos->Add(fPtBinNplusNminusChTruth);
173     
174   }//MCAOD---
175   
176  
177   
178   //PostData(1, fListOfHistosQA);
179   PostData(1, fListOfHistos);   
180   
181   
182 }
183
184
185 //----------------------------------------------------------------------------------
186 void AliEbyEHigherMomentsTask::UserExec( Option_t * ){
187   
188   fEventCounter->Fill(1);
189   
190   if(fAnalysisType == "AOD") {
191     
192     doAODEvent();
193     
194   }//AOD--analysis-----
195   
196   else if(fAnalysisType == "MCAOD") {
197     
198     doMCAODEvent();
199     
200   }
201   
202   else return;
203  
204   
205   
206 }
207
208 //--------------------------------------------------------------------------------------
209 void AliEbyEHigherMomentsTask::doAODEvent(){
210   
211   Double_t positiveSum = 0.;
212   Double_t negativeSum = 0.;
213   const Int_t dim = fNptBins*2;
214   Double_t PtCh[dim];
215
216   for(Int_t idx = 0; idx < dim; idx++){
217     PtCh[idx] = 0.;
218   }
219
220
221   AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
222   if (!manager) {
223     cout<<"ERROR: Analysis manager not found."<<endl;
224     return;
225   }
226   //coneect to the inputHandler------------
227   AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
228   if (!inputHandler) {
229     cout<<"ERROR: Input handler not found."<<endl;
230     return;
231   }
232   
233   AliAODEvent *fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
234
235   if (!fAOD)
236     {
237       cout<< "ERROR: AOD not found " <<endl;
238       return;
239     }
240   
241   if(!ProperVertex(fAOD)) return;   
242   
243   AliAODHeader *aodHeader = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
244   if(!aodHeader) AliFatal("Not a standard AOD");
245   
246   fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
247   
248   if(fCentrality < 0 || fCentrality >= 81) return;
249   
250   fEventCounter->Fill(2);
251   
252   Int_t nTracks = fAOD->GetNumberOfTracks();
253   
254   for(Int_t i = 0; i < nTracks; i++) {
255     
256     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i)); 
257     
258     if(!aodTrack) {
259       AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
260       continue;
261     }
262     
263     if(!AcceptTrack(aodTrack) ) continue;
264
265     
266     fBin = GetPtBin(aodTrack->Pt());
267     
268     Short_t gCharge = aodTrack->Charge();
269     if(gCharge > 0){
270       positiveSum += 1.;
271       PtCh[fBin] += 1;
272     }
273     if(gCharge < 0){
274       negativeSum += 1.;
275       PtCh[fNptBins+fBin] += 1;
276     }
277    
278    
279   }//--------- Track Loop to select with filterbit
280   
281   Double_t fContainerCh[3] = { static_cast<Double_t>(fCentrality), positiveSum, negativeSum};
282   
283   fTHnCentNplusNminusCh->Fill(fContainerCh);
284
285   //cout << fCentrality << " " << positiveSum <<" " << negativeSum << endl;
286   
287   Double_t ptContainer[dim+1];
288
289   ptContainer[0] = fCentrality;
290
291   for(Int_t i = 1; i <= dim; i++){
292     ptContainer[i] = PtCh[i-1];
293     //cout << PtCh[i-1] <<" ,";
294   }
295   //cout << endl;
296
297   fPtBinNplusNminusCh->Fill(ptContainer);
298   
299   
300   fEventCounter->Fill(7);
301   PostData(1, fListOfHistos);
302   return;
303   
304 }
305 //--------------------------------------------------------------------------------------
306 //--------------------------------------------------------------------------------------
307 void AliEbyEHigherMomentsTask::doMCAODEvent(){
308   
309   
310   //---------
311   Double_t positiveSumMCRec = 0.;
312   Double_t negativeSumMCRec = 0.;
313  
314   Double_t positiveSumMCTruth = 0.;
315   Double_t negativeSumMCTruth = 0.;
316  
317   const Int_t dim = fNptBins*2;
318   Double_t PtCh[dim];
319   Double_t ptChMC[dim];
320   for(Int_t idx = 0; idx < dim; idx++){
321     PtCh[idx] = 0;
322     ptChMC[idx] = 0;
323   }
324   
325   //Connect to Anlaysis manager------
326   AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
327   if (!manager) {
328     cout<<"ERROR: Analysis manager not found."<<endl;
329     return;
330   }
331   
332   AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
333   if (!inputHandler) {
334     cout<<"ERROR: Input handler not found."<<endl;
335     return;
336   }
337   
338   //AOD event------
339   AliAODEvent *fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
340   
341   if (!fAOD)
342     {
343       cout<< "ERROR: AOD not found " <<endl;
344       return;
345     }
346
347   // -- Get MC header
348   // ------------------------------------------------------------------
349   
350   fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
351   if (!fArrayMC) {
352     AliFatal("Error: MC particles branch not found!\n");
353     return;
354   }
355   
356   AliAODMCHeader *mcHdr=NULL;
357   mcHdr=(AliAODMCHeader*)fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());  
358   if(!mcHdr) {
359     printf("MC header branch not found!\n");
360     return;
361   }
362   
363   if(!ProperVertex(fAOD)) return;
364
365   AliAODHeader *aodHeader = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
366   if(!aodHeader) AliFatal("Not a standard AOD");
367   
368   fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
369
370   if( fCentrality < 0 || fCentrality >= 81) return;
371   
372   fEventCounter->Fill(2);
373  
374   Int_t nTracks = fAOD->GetNumberOfTracks();
375   
376   for(Int_t i = 0; i < nTracks; i++) {
377     
378     AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i)); 
379     
380     if(!aodTrack) {
381       AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
382       continue;
383     }
384     
385     if(!AcceptTrack(aodTrack) ) continue;
386
387     
388     fBin = GetPtBin(aodTrack->Pt());
389     
390     //cout << "Pt Bin " << fBin << endl;
391
392     Short_t gCharge = aodTrack->Charge();
393     if(gCharge > 0){
394       positiveSumMCRec += 1.;
395      PtCh[fBin] += 1;
396     }
397     if(gCharge < 0){
398       negativeSumMCRec += 1.;    
399       PtCh[fNptBins+fBin] += 1;
400     }
401    
402   }//--------- Track Loop to select with filterbit
403   
404   //===============================================================================
405   //---------------------MC Truth--------------------------------------------------
406   //===============================================================================
407   
408   Int_t nMCTrack = fArrayMC->GetEntriesFast();
409   
410   for (Int_t iMC = 0; iMC < nMCTrack; iMC++) {
411     
412     AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
413     
414     if(!partMC){
415       AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
416       continue;
417     }
418     
419     if(partMC->Charge() == 0) continue;
420     if(!partMC->IsPhysicalPrimary()) continue;
421     
422     if (partMC->Eta() < fEtaLowerLimit || partMC->Eta() > fEtaHigherLimit) continue;
423     if (partMC->Pt() < fPtLowerLimit ||  partMC->Pt() > fPtHigherLimit) continue;
424     
425     Short_t gCharge = partMC->Charge();
426     Int_t mcbin = GetPtBin(partMC->Pt());
427     
428     if(gCharge > 0){
429       positiveSumMCTruth += 1.;
430       ptChMC[mcbin] += 1.;
431     }
432     if(gCharge < 0){
433       negativeSumMCTruth += 1.;
434       ptChMC[fNptBins+mcbin] += 1.;
435     }
436    
437    
438   }//MC-Truth Track loop--
439
440   Double_t fContainerCh[3] = { static_cast<Double_t>(fCentrality), positiveSumMCRec, negativeSumMCRec};//Reco. values ch. hadrons
441   Double_t fContainerChTruth[3] = { static_cast<Double_t>(fCentrality), positiveSumMCTruth, negativeSumMCTruth};  
442  
443   fTHnCentNplusNminusCh->Fill(fContainerCh);//Fill the rec. ch. particles---
444   fTHnCentNplusNminusChTruth->Fill(fContainerChTruth);//MC -Truth ch. particles
445   
446   //cout << fCentrality << " " << positiveSumMCRec << " " << negativeSumMCRec <<endl;
447   //cout <<fCentrality<<"   " << positiveSumMCTruth << " " << negativeSumMCTruth << endl;
448   //cout <<"   " << posPidSumMCRec << " " << negPidSumMCRec << endl;
449   //cout <<"   " << posPidSumMCTruth << " " << negPidSumMCTruth << endl;
450   //cout <<"---------------------------------" << endl;
451   
452   Double_t ptContainer[dim+1];
453   Double_t ptContainerMC[dim+1];
454   ptContainer[0] = fCentrality;
455   ptContainerMC[0] = fCentrality;
456
457   for(Int_t i = 1; i <= dim; i++){
458     ptContainer[i] = PtCh[i-1];
459     ptContainerMC[i] = ptChMC[i-1];
460     //cout <<" "<< PtCh[i-1]<<endl;
461     //  cout<< " MC=" << ptChMC[i-1];
462   }
463
464   //cout << endl;
465
466   fPtBinNplusNminusCh->Fill(ptContainer);
467   fPtBinNplusNminusChTruth->Fill(ptContainerMC);
468
469   fEventCounter->Fill(7);
470   PostData(1, fListOfHistos);
471   return;
472   
473 }
474
475 //---------------------------------------------------------------------------------------
476 Bool_t AliEbyEHigherMomentsTask::ProperVertex(AliAODEvent *fAOD) const{
477   
478   Bool_t IsVtx = kFALSE;
479   
480   const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
481   
482   if(vertex) {
483     Double32_t fCov[6];
484     vertex->GetCovarianceMatrix(fCov);
485     if(vertex->GetNContributors() > 0) {
486       if(fCov[5] != 0) {
487         
488         Double_t lvx = vertex->GetX();
489         Double_t lvy = vertex->GetY();
490         Double_t lvz = vertex->GetZ();
491         
492         fEventCounter->Fill(5);
493       
494         if(TMath::Abs(lvx) < fVxMax) {
495           if(TMath::Abs(lvy) < fVyMax) {
496             if(TMath::Abs(lvz) < fVzMax) {
497               
498               fEventCounter->Fill(6);
499               fHistQA[0]->Fill(lvz);
500               fHistVxVy->Fill(lvx,lvy);
501               IsVtx = kTRUE; 
502               
503             }//Z-Vertex cut---
504           }//Y-vertex cut--
505         }//X-vertex cut---
506       }//Covariance------
507     }//Contributors check---
508   }//If vertex-----------
509
510   AliCentrality *centrality = fAOD->GetCentrality();
511   if (centrality->GetQuality() != 0) IsVtx = kFALSE;
512   
513   return IsVtx;
514 }
515 //---------------------------------------------------------------------------------------
516
517 //---------------------------------------------------------------------------------------
518 Bool_t AliEbyEHigherMomentsTask::AcceptTrack(AliAODTrack* track) const{
519   
520   if(!track) return kFALSE;
521   if( track->Charge() == 0 ) return kFALSE;
522
523   Double_t pt = track->Pt();
524   Double_t eta = track->Eta();
525   if(!track->TestFilterBit(fAODtrackCutBit)) return kFALSE;
526   if( pt < fPtLowerLimit || pt > fPtHigherLimit ) return kFALSE;
527   if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) return kFALSE;
528   
529   fHistQA[1]->Fill(pt);
530   fHistQA[2]->Fill(eta);
531   fHistQA[3]->Fill(track->Phi());
532  
533   return kTRUE;
534 }
535 //------------------------------------------------------------------------
536
537 //------------------------------------------------------------------------
538 Int_t AliEbyEHigherMomentsTask::GetPtBin(Double_t pt){
539   
540   Int_t bin = 0;
541
542   Double_t BinSize = (fPtHigherLimit - fPtLowerLimit)/fNptBins;
543   
544   for(Int_t iBin = 0; iBin < fNptBins; iBin++){
545     
546     Double_t xLow = fPtLowerLimit + iBin*BinSize;
547     Double_t xHigh = fPtLowerLimit + (iBin+1)*BinSize;
548     
549     if( pt >= xLow && pt < xHigh){
550       bin = iBin;
551       //cout << "Pt Bin #" << bin <<"pt=" << pt << endl;
552       break;
553     }
554     
555   }//for 
556   
557   
558   return bin;
559 }
560 //------------------------------------------------------------------------
561 //------------------------------------------------------------------------
562 void AliEbyEHigherMomentsTask::Terminate( Option_t * ){
563   
564   Info("AliEbyEHigherMomentTask"," Task Successfully finished");
565   
566 }
567
568 //------------------------------------------------------------------------