Update PR task: drathee
[u/mrichter/AliRoot.git] / PWGCF / EBYE / PIDFluctuation / task / AliEbyEPidTaskFastGen.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 Fluctuaion Analysis for PID                     //
19 //                   Deepika Jena | drathee@cern.ch                        //
20 //                   Satyajit Jena | sjena@cern.ch                         //
21 //=========================================================================//
22
23 #include "TChain.h"
24 #include "TList.h"
25 #include "TFile.h"
26 #include "TTree.h"
27 #include "TH1D.h"
28 #include "TH2F.h"
29 #include "TH3F.h"
30
31 #include "THnSparse.h"
32
33 #include "AliAnalysisTask.h"
34 #include "AliAnalysisManager.h"
35
36 #include "AliVEvent.h"
37 #include "AliESDEvent.h"
38 #include "AliMCEvent.h"
39 #include "AliAODEvent.h"
40
41 #include "AliStack.h"
42 #include "AliGenEventHeader.h"
43 #include "AliGenPythiaEventHeader.h"
44 #include "AliGenHijingEventHeader.h"
45 #include "AliGenDPMjetEventHeader.h"
46 #include "TDatabasePDG.h"
47
48 #include "AliEbyEPidTaskFastGen.h"
49
50 ClassImp(AliEbyEPidTaskFastGen)
51
52 //-----------------------------------------------------------------------
53 AliEbyEPidTaskFastGen::AliEbyEPidTaskFastGen( const char *name )
54 : AliAnalysisTaskSE( name ),
55   fThnList(0),
56   fVxMax(3.),
57   fVyMax(3.),
58   fVzMax(15.),
59   fPtLowerLimit(0.2),
60   fPtHigherLimit(5.),
61   fEtaLowerLimit(-1.),
62   fEtaHigherLimit(1.),
63   fHistoCorrelationMC(NULL)
64
65  
66   DefineOutput(1, TList::Class()); 
67 }
68
69 AliEbyEPidTaskFastGen::~AliEbyEPidTaskFastGen()
70 {
71   if (fThnList) delete fThnList;
72 }
73
74
75 //---------------------------------------------------------------------------------
76 void AliEbyEPidTaskFastGen::UserCreateOutputObjects() {
77   fThnList = new TList();
78   fThnList->SetOwner(kTRUE);
79   Int_t fgSparseDataBins[kNSparseData]   = {50000,   500,     440,  30000,   15000,  15000,   10000,  5000,   5000,  5000,  2500,  2500,  2500,  1250, 1250};
80   Double_t fgSparseDataMin[kNSparseData] = {0.,       0.,      0.,     0.,      0.,     0.,      0.,    0.,     0.,    0.,    0.,    0.,    0.,   0.,    0.};
81   Double_t fgSparseDataMax[kNSparseData] = {50000.,  100.,    440., 30000.,  15000., 15000.,  10000., 5000.,  5000., 5000., 2500., 2500., 2500., 1250., 1250.};
82   
83   const Char_t *fgkSparseDataTitle[] = {
84     "Nnumber of Partiles",
85     "Impact Parameters",
86     "N_{p}",
87     "N_{ch}",
88     "N_{+}",
89     "N_{-}",
90     "N_{#pi}",
91     "N_{#pi^{+}}",
92     "N_{#pi^{-}}",
93     "N_{K}",
94     "N_{K^{+}}",
95     "N_{K^{-}}",
96     "N_{pr}",
97     "N_{p}",
98     "N_{#bar{p}}"
99   };
100   
101   fHistoCorrelationMC = new THnSparseI("hHistoPR", "", kNSparseData, fgSparseDataBins, fgSparseDataMin, fgSparseDataMax);
102   for (Int_t iaxis = 0; iaxis < kNSparseData; iaxis++)
103     fHistoCorrelationMC->GetAxis(iaxis)->SetTitle(fgkSparseDataTitle[iaxis]);
104   fThnList->Add(fHistoCorrelationMC);
105   
106   PostData(1, fThnList);
107 }
108
109
110 //----------------------------------------------------------------------------------
111 void AliEbyEPidTaskFastGen::UserExec( Option_t * ){
112   AliMCEvent* mcEvent = MCEvent();
113   if (!mcEvent) {
114     Printf("ERROR: Could not retrieve MC event");
115     return;
116   }
117   
118   AliStack *stack = mcEvent->Stack();
119  
120   AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
121   if(!genHeader){
122     printf("  Event generator header not available!!!\n");
123     return;
124   }
125   
126   Double_t imp = -1;
127   Double_t np = -1;
128   Double_t nt = -1;
129   
130   if(genHeader->InheritsFrom(AliGenHijingEventHeader::Class())){
131     imp = ((AliGenHijingEventHeader*) genHeader)->ImpactParameter();
132     nt  = ((AliGenHijingEventHeader*) genHeader)->TargetParticipants();
133     np  = ((AliGenHijingEventHeader*) genHeader)->ProjectileParticipants();
134   }  
135   
136
137   Int_t count[4][2];
138   for (Int_t i = 0; i < 4; i++) {
139     for (Int_t j = 0; j < 2; j++) count[i][j] = 0; 
140   }
141   
142   Int_t nParticles = stack->GetNtrack();
143   for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
144     TParticle* part = stack->Particle(iParticle);
145     if (!part) {
146       Printf(" No Particle Available ");
147       continue;
148     }
149     
150     //  TDatabasePDG* pdgDB   = TDatabasePDG::Instance();
151     //  TParticlePDG* pdgPart = pdgDB->GetParticle(part->GetPdgCode());
152     
153     Float_t pt = part->Pt();
154     if(pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
155     Float_t gEta = part->Eta();
156     if(gEta < fEtaLowerLimit || gEta > fEtaHigherLimit ) continue;
157     
158     
159     //  Float_t gCharge = (pdgPart ? pdgPart->Charge() : 0);
160     
161     Int_t pid = part->GetPdgCode();
162     if(pid == 211)         { count[1][0]++;  count[0][0]++; } 
163     else if(pid ==  -211)  { count[1][1]++;  count[0][1]++; } 
164     else if(pid ==   321)  { count[2][0]++;  count[0][0]++; } 
165     else if(pid ==  -321)  { count[2][1]++;  count[0][1]++; } 
166     else if(pid ==  2212)  { count[3][0]++;  count[0][0]++; } 
167     else if(pid == -2212)  { count[3][1]++;  count[0][1]++; } 
168   }
169   
170   Double_t vsparseMC[kNSparseData];
171   vsparseMC[0]  = nParticles;
172   vsparseMC[1]  = imp;
173   vsparseMC[2]  = np + nt;
174   vsparseMC[3]  = count[0][0] + count[0][1];
175   vsparseMC[4]  = count[0][0];
176   vsparseMC[5]  = count[0][1];
177
178   vsparseMC[6]  = count[1][0] + count[1][1];
179   vsparseMC[7]  = count[1][0];
180   vsparseMC[8]  = count[1][1];
181
182   vsparseMC[9]  = count[2][0] + count[2][1];
183   vsparseMC[10]  = count[2][0];
184   vsparseMC[11]  = count[2][1];
185
186   vsparseMC[12] = count[3][0] + count[3][1];
187   vsparseMC[13] = count[3][0];
188   vsparseMC[14] = count[3][1];
189   
190
191    /* Printf(" %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f %10.2f",
192          vsparseMC[0], vsparseMC[1], vsparseMC[2], 
193          vsparseMC[3], vsparseMC[4], vsparseMC[5], 
194          vsparseMC[6], vsparseMC[7], vsparseMC[8], 
195          vsparseMC[9], vsparseMC[10], vsparseMC[11],
196          vsparseMC[12], vsparseMC[13], vsparseMC[13]);
197    */
198   fHistoCorrelationMC->Fill(vsparseMC);
199   PostData(1, fThnList);
200   
201 }
202
203 void AliEbyEPidTaskFastGen::Terminate( Option_t * ){
204
205   Info("AliEbyEPiDTask"," Task Successfully finished");
206   
207 }
208
209