Coverity fix: 8443 - sjena
[u/mrichter/AliRoot.git] / PWGCF / EBYE / PIDFluctuation / task / AliEbyEPidRatioTaskOnFly.cxx
CommitLineData
5a289b9a 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#include "TProfile.h"
31
32#include "AliAnalysisTask.h"
33#include "AliAnalysisManager.h"
34
35#include "AliVEvent.h"
36#include "AliESDEvent.h"
37#include "AliMCEvent.h"
38#include "AliAODEvent.h"
39
40#include "AliStack.h"
41#include "AliGenEventHeader.h"
42#include "AliGenPythiaEventHeader.h"
43#include "AliGenHijingEventHeader.h"
44#include "AliGenDPMjetEventHeader.h"
45#include "TDatabasePDG.h"
46
47#include "AliEbyEPidRatioTaskOnFly.h"
48
49ClassImp(AliEbyEPidRatioTaskOnFly)
50
51//-----------------------------------------------------------------------
52AliEbyEPidRatioTaskOnFly::AliEbyEPidRatioTaskOnFly( const char *name )
53: AliAnalysisTaskSE( name ),
54 fThnList(0),
55 fPtLowerLimit(0.2),
56 fPtHigherLimit(5.),
57 fEtaLowerLimit(-1.),
58 fEtaHigherLimit(1.),
baacd542 59 fCentrality(-1),
5a289b9a 60 fOrder(8),
61 fRedFactp(NULL)
62
63
64{
65
66 DefineOutput(1, TList::Class());
67}
68
69AliEbyEPidRatioTaskOnFly::~AliEbyEPidRatioTaskOnFly()
70{
71 if (fThnList) delete fThnList;
72
73 for (Int_t ii = 0; ii <= fOrder; ++ii)
74 if (fRedFactp[ii]) delete[] fRedFactp[ii];
75 if (fRedFactp) delete[] fRedFactp;
76
77}
78
79
80
81
82//---------------------------------------------------------------------------------
83void AliEbyEPidRatioTaskOnFly::UserCreateOutputObjects() {
84 fThnList = new TList();
85 fThnList->SetOwner(kTRUE);
86
87 const Char_t *name = "Mc";
88 const Char_t *title = Form(" #eta [%2.1f-%2.1f]",fEtaLowerLimit,fEtaHigherLimit);
89
90
91 TString sName(name);
92 TString sTitle(title);
93
94
95 fRedFactp = new Double_t*[fOrder+1];
96 for (Int_t ii = 0 ; ii <= fOrder; ++ii)
97 fRedFactp[ii] = new Double_t[2];
98
99 //TList *list[4];
100 fThnList->Add(new TList);
101 TList *list = static_cast<TList*>(fThnList->Last());
102 list->SetName(Form("f%s",name));
103 list->SetOwner(kTRUE);
104
105 for (Int_t iPid = 0; iPid < 4; ++iPid) {
106 TString sNetTitle(Form("%s - %s",fgkPidLatex[iPid][1],fgkPidLatex[iPid][0]));
107
108 list->Add(new TProfile(Form("fProfTot%sPlus%s",fgkPidName[iPid],name),
109 Form("(%s) : %s;Centrality(100);(%s)",fgkPidName[iPid], sTitle.Data(), sNetTitle.Data()),
110 20,-0.5,19.5));
111
112 list->Add(new TProfile(Form("fProfTot%sMinus%s",fgkPidName[iPid],name),
113 Form("(%s) : %s;Centrality(100);(%s)",fgkPidName[iPid], sTitle.Data(), sNetTitle.Data()),
114 20,-0.5,19.5));
115
116 for (Int_t idx = 1; idx <= fOrder; ++idx) {
117 list->Add(new TProfile(Form("fProf%s%sNet%dM",fgkPidName[iPid],name, idx),
118 Form("(%s)^{%d} : %s;Centrality(100);(%s)^{%d}",sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx),
119 20,-0.5,19.5));
120 }
121
122 for (Int_t ii = 0; ii <= fOrder; ++ii) {
123 for (Int_t kk = 0; kk <= fOrder; ++kk) {
124 list->Add(new TProfile(Form("fProf%s%sNetF%02d%02d",fgkPidName[iPid], name, ii, kk),
125 Form("f_{%02d%02d} : %s;Centrality(100);f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk),
126 20,-0.5,19.5));
127 }
128 }
129
130 }
131
132
133 for (Int_t iPhy = 0; iPhy < 46; ++iPhy) {
134 list->Add(new TProfile(Form("fProf%sNu%02d",name,iPhy),Form("Physics Variable for index %d | %s ; Centrality;",iPhy,name),20,-0.5,19.5));
135 }
136
137
138 PostData(1, fThnList);
139}
140
141
142//----------------------------------------------------------------------------------
143void AliEbyEPidRatioTaskOnFly::UserExec( Option_t * ){
144 AliMCEvent* mcEvent = MCEvent();
145 if (!mcEvent) {
146 Printf("ERROR: Could not retrieve MC event");
147 return;
148 }
149
150 for (Int_t ii = 0 ; ii < 4; ++ii)
151 for (Int_t kk = 0 ; kk < 2; ++kk)
152 fNp[ii][kk] = 0;
153
154
155 const AliVVertex *vtxMC = mcEvent->GetPrimaryVertex();
156 if (vtxMC->GetZ() > 10.) return;
157
158 AliStack *stack = mcEvent->Stack();
159
160 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
161 if(!genHeader){
162 printf(" Event generator header not available!!!\n");
163 return;
164 }
165
166 Double_t imp = 21;
167
168 if(genHeader->InheritsFrom(AliGenHijingEventHeader::Class())){
169 imp = ((AliGenHijingEventHeader*) genHeader)->ImpactParameter();
170 }
171 if (imp < 20) fCentrality = Int_t(imp);
172 else return;
173
174 Int_t nParticles = stack->GetNtrack();
175 for (Int_t iParticle = 0; iParticle < nParticles; iParticle++) {
176 TParticle* part = stack->Particle(iParticle);
177 if (!part) {
178 Printf(" No Particle Available ");
179 continue;
180 }
181
182 Float_t pt = part->Pt();
183 if(pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
184 Float_t gEta = part->Eta();
185 if(gEta < fEtaLowerLimit || gEta > fEtaHigherLimit ) continue;
186
187 // Float_t gCharge = (pdgPart ? pdgPart->Charge() : 0);
188
189 Int_t pid = part->GetPdgCode();
190 if(pid == 211) { fNp[1][1]++; fNp[0][1]++; }
191 else if(pid == -211) { fNp[1][0]++; fNp[0][0]++; }
192 else if(pid == 321) { fNp[2][1]++; fNp[0][1]++; }
193 else if(pid == -321) { fNp[2][0]++; fNp[0][0]++; }
194 else if(pid == 2212) { fNp[3][1]++; fNp[0][1]++; }
195 else if(pid == -2212) { fNp[3][0]++; fNp[0][0]++; }
196
197 }
198
199 /*Printf("%6d %6d %6d %6d %6d %6d %6d %6d %6d", fCentrality,
200 fNp[0][0], fNp[0][1],
201 fNp[1][0], fNp[1][1],
202 fNp[2][0], fNp[2][1],
203 fNp[3][0], fNp[3][1]);*/
204
205
206 FillHistSetCent();
207 PostData(1, fThnList);
208}
209
210void AliEbyEPidRatioTaskOnFly::Terminate( Option_t * ){
211 Info("AliEbyEPidRatioTaskOnFly"," Task Successfully finished");
212}
213//________________________________________________________________________
214const Char_t* AliEbyEPidRatioTaskOnFly::fgkPidName[4] = {"Nch","Npi","Nka","Npr"};
215//________________________________________________________________________
216const Char_t* AliEbyEPidRatioTaskOnFly::fgkPidLatex[4][2] = {{"N_{-}","N_{+}"}, {"N_{#pi^{-}}","N_{#pi^{+}}"},{"N_{K^{-}}","N_{K^{+}}"}, {"N_{#bar{p}}","N_{p}"}};
217//________________________________________________________________________
218const Char_t* AliEbyEPidRatioTaskOnFly::fgkPidTitles[4][2] = {{"Negative","Positive"},{"Anti-Pions","Pions"},{"Anti-Kaons","Kaons"}, {"Anti-Protons","Protons"}};
219//________________________________________________________________________
220void AliEbyEPidRatioTaskOnFly::FillHistSetCent() {
221
222 const Char_t *name = "Mc";
223
224 TList *list = static_cast<TList*>(fThnList->FindObject(Form("f%s",name)));
225
226 for (Int_t iPid = 0; iPid < 4; ++iPid) {
227 Int_t deltaNp = fNp[iPid][1]-fNp[iPid][0];
228 Double_t delta = 1.;
229 (static_cast<TProfile*>(list->FindObject(Form("fProfTot%sPlus%s", fgkPidName[iPid], name))))->Fill(fCentrality, fNp[iPid][1]);
230 (static_cast<TProfile*>(list->FindObject(Form("fProfTot%sMinus%s", fgkPidName[iPid], name))))->Fill(fCentrality, fNp[iPid][0]);
231 for (Int_t idxOrder = 1; idxOrder <= fOrder; ++idxOrder) {
232 delta *= deltaNp;
233 (static_cast<TProfile*>(list->FindObject(Form("fProf%s%sNet%dM", fgkPidName[iPid], name, idxOrder))))->Fill(fCentrality, delta);
234 }
235
236 for (Int_t idxOrder = 0; idxOrder <= fOrder; ++ idxOrder) {
237 fRedFactp[idxOrder][0] = 1.;
238 fRedFactp[idxOrder][1] = 1.;
239 }
240
241 for (Int_t idxOrder = 1; idxOrder <= fOrder; ++ idxOrder) {
242 fRedFactp[idxOrder][0] = fRedFactp[idxOrder-1][0] * Double_t(fNp[iPid][0]-(idxOrder-1));
243 fRedFactp[idxOrder][1] = fRedFactp[idxOrder-1][1] * Double_t(fNp[iPid][1]-(idxOrder-1));
244 }
245
246 for (Int_t ii = 0; ii <= fOrder; ++ii) {
247 for (Int_t kk = 0; kk <= fOrder; ++kk) {
248 Double_t fik = fRedFactp[ii][1] * fRedFactp[kk][0];
249 (static_cast<TProfile*>(list->FindObject(Form("fProf%s%sNetF%02d%02d", fgkPidName[iPid], name, ii, kk))))->Fill(fCentrality, fik);
250 }
251 }
252 }
253
254 //Printf("%6d %20s %6.2f %6d %6d %6d %6d %6d %6d %6d %6d", idx, name, centralityBin,
255 // fNp[0][1], fNp[0][0],
256 // fNp[1][1], fNp[1][0],
257 /// fNp[2][1], fNp[2][0],
258 // fNp[3][1], fNp[3][0]);
259 //
260
261 Int_t a[6][4]; Int_t b[22];
262 for (Int_t iPid = 0; iPid < 4; ++iPid) {
263 a[0][iPid] = fNp[iPid][1]+fNp[iPid][0]; // 0 n+ + n-
264 a[1][iPid] = fNp[iPid][1]; // 1 n+
265 a[2][iPid] = fNp[iPid][0]; // 2 n-
266 a[3][iPid] = fNp[iPid][1]*fNp[iPid][0]; // 3 n+ . n-
267 a[4][iPid] = fNp[iPid][1]*(fNp[iPid][1]-1); // 4 n+ (n+ - 1)
268 a[5][iPid] = fNp[iPid][0]*(fNp[iPid][0]-1); // 5 n- (n- - 1)
269
270 // Printf("%6d %20s %6.2f %6d %6d %6d ", idx, name, centralityBin,
271 // a[0][iPid], a[1][iPid], a[2][iPid]);
272
273 }
274
275 b[0] = a[0][0]*a[0][2]; // 24 N K
276 b[1] = a[0][1]*a[0][2]; // 25 Pi K
277 b[2] = a[1][1]*a[1][2]; // 26 pi+ k+
278 b[3] = a[1][1]*a[2][2]; // 27 pi+ k-
279 b[4] = a[2][1]*a[1][2]; // 28 pi- k+
280 b[5] = a[2][1]*a[2][2]; // 29 pi- k-
281
282 b[6] = a[0][0]*a[0][3]; // 30 N P
283 b[7] = a[0][2]*a[0][3]; // 31 K P
284 b[8] = a[1][2]*a[1][3]; // 32 k+ p+
285 b[9] = a[1][2]*a[2][3]; // 33 k+ p-
286 b[10] = a[2][2]*a[1][3]; // 34 k- p+
287 b[11] = a[2][2]*a[2][3]; // 35 k- p-
288
289 b[12] = a[0][0]*a[0][1]; // 36 N Pi
290 b[13] = a[0][3]*a[0][1]; // 37 P Pi
291 b[14] = a[1][3]*a[1][1]; // 38 p+ pi+
292 b[15] = a[1][3]*a[2][1]; // 39 p+ pi-
293 b[16] = a[2][3]*a[1][1]; // 40 p- pi+
294 b[17] = a[2][3]*a[2][1]; // 41 p- pi-
295
296 b[18] = a[0][0]*(a[0][0] - 1); // 42 N ( N - 1 )
297 b[19] = a[0][1]*(a[0][1] - 1); // 43 Pi( Pi- 1 )
298 b[20] = a[0][2]*(a[0][1] - 1); // 44 K ( K - 1 )
299 b[21] = a[0][3]*(a[0][3] - 1); // 45 P ( P - 1 )
300 // TList *list_nu = static_cast<TList*>(fOutList->FindObject(Form("f%s_nu",name)));
301 Int_t k = 0;
302 for (Int_t j = 0; j < 4; j++) {
303 for (Int_t i = 0; i < 6; i++) {
304 (static_cast<TProfile*>(list->FindObject(Form("fProf%sNu%02d", name,k))))->Fill(fCentrality,a[i][j]);
305 k++;
306 }
307 }
308
309 for (Int_t j = 0; j < 22; j++) {
310 (static_cast<TProfile*>(list->FindObject(Form("fProf%sNu%02d", name,j+24))))->Fill(fCentrality,b[j]);
311 }
312
313 /*
314 Printf("%6d %6d %6d %6d %6d %6d %6d %6d %6d", fCentrality,
315 fNp[0][0], fNp[0][1],
316 fNp[1][0], fNp[1][1],
317 fNp[2][0], fNp[2][1],
318 fNp[3][0], fNp[3][1]);
319 */
320
321 return;
322}
323