coverity fix
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / MultEvShape / AliAnalysisCentralCutEvtMC.cxx
CommitLineData
bde49c8a 1/*************************************************************************
2* Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3* *
4* Author: The ALICE Off-line Project. *
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
410ff8cc 16// ---------------------------------------------------
17// MC event level cuts for azimuthal isotropic
18// expansion in highly central collisions analysis
19// author: Cristian ANDREI
20// acristian@niham.nipne.ro
21// ----------------------------------------------------
22
bde49c8a 23
24#include <TParticle.h>
25
26#include "AliMCEvent.h"
27#include "AliStack.h"
28
29#include "AliAnalysisCentralCutEvtMC.h"
410ff8cc 30
31
32class TObject;
33
bde49c8a 34//____________________________________________________________________
35ClassImp(AliAnalysisCentralCutEvtMC)
36
37//____________________________________________________________________
38AliAnalysisCentralCutEvtMC::AliAnalysisCentralCutEvtMC(const Char_t* name, const Char_t* title)
39 :AliAnalysisCuts(name,title)
40 ,fReqMult(kFALSE)
41 ,fReqDir(kFALSE)
42 ,fReqDirUnit(kFALSE)
43 ,fMultMin(0)
44 ,fMultMax(0)
45 ,fDirMin(0)
46 ,fDirMax(0)
47 ,fDirUMin(0)
48 ,fDirUMax(0)
49
50{
51 //constructor
0fb1f0cf 52
bde49c8a 53}
54
55AliAnalysisCentralCutEvtMC::~AliAnalysisCentralCutEvtMC(){
56//destructor
0fb1f0cf 57
bde49c8a 58}
59
60//___________________________________________________________________________
9eeae5d5 61Bool_t AliAnalysisCentralCutEvtMC::IsSelected(TObject *obj){
bde49c8a 62// Checks whether the event passes the cuts
63
64 if (!obj){
65 printf("Pointer to MC Event is null! \n");
66 return kFALSE;
67 }
68
69 AliMCEvent *mcEvent = dynamic_cast<AliMCEvent*> (obj);
70 if(!mcEvent){
71 printf("AliAnalysisCentralCutEvtMC -> Can't get MCEvt!\n");
72 return kFALSE;
73 }
74
75
76 if(fReqMult){
77 Int_t mult = CalcMult(mcEvent);
78 if((mult<fMultMin)||(mult>fMultMax)){
79 return kFALSE;
80 }
81 }
82
83 if(fReqDir){
84 Double_t dir = CalcDir(mcEvent);
85 if((dir<fDirMin)||(dir>fDirMax)){
86 return kFALSE;
87 }
88 }
89
90 if(fReqDirUnit){
91 Double_t dirU = CalcDirUnit(mcEvent);
92 if((dirU<fDirUMin)||(dirU>fDirUMax)){
93 return kFALSE;
94 }
95 }
96
97 return kTRUE;
98}
99
100//___________________________________________________________________________
101Int_t AliAnalysisCentralCutEvtMC::CalcMult(AliMCEvent* const mcEvent) {
102// Computes the multiplicity of the event
103
104 AliStack *stack=mcEvent->Stack();
105 if(!stack) return -1;
106
107 Int_t nPrim = stack->GetNprimary();
108
109 Int_t charged = 0;
110 Double_t eta;
111
112 for(Int_t i = 0; i < nPrim; i++){//track loop -> compute multiplicity
113
114 TParticle *particle=stack->Particle(i);
115 if(!particle){
116 continue;
117 }
118
119 eta = particle->Eta();
120
121 if((eta > 0.5)||(eta < -0.5)){ //------------mid-rapidity
122 continue;
123 }
124
125 TParticlePDG *particlePDG = particle->GetPDG(0);
126 if(!particlePDG){
127 continue;
128 }
129
130 if(TMath::Abs(particlePDG->Charge()) < 3){
131 continue;
132 }
133
134 charged++;
135
136 }//end track loop
137
138 return charged;
139
140}
141
142//____________________________________________________________________________
143Double_t AliAnalysisCentralCutEvtMC::CalcDir(AliMCEvent* const mcEvent) {
144// computes the directivity of the event
145
146 AliStack *stack=mcEvent->Stack();
147 if(!stack) return -1;
148
149 Int_t nPrim = stack->GetNprimary();
150
151 Int_t goodtrack = 0;
152
153 Double_t eta;
154 Double_t pt;
155
156 Double_t dir;
157 Double_t px,py;
158 Double_t sumaPt = 0;
159 Double_t sumaPx = 0;
160 Double_t sumaPy = 0;
161
162 for(Int_t i = 0; i < nPrim; i++){//track loop -> compute directivity
163
164 TParticle *particle=stack->Particle(i);
165 if(!particle){
166 continue;
167 }
168
169 eta = particle->Eta();
170
171 if((eta > 1.9)||(eta < 0.0)){ // half of the SPD coverage -> directivity
172 continue;
173 }
174
175 TParticlePDG *particlePDG = particle->GetPDG(0);
176 if(!particlePDG){
177 continue;
178 }
179
180 if(TMath::Abs(particlePDG->Charge()) < 3){
181 continue;
182 }
183
184 px = particle->Px();
185 py = particle->Py();
186
187 pt = particle->Pt();
188
189 sumaPx = sumaPx + px;
190 sumaPy = sumaPy + py;
191
192 sumaPt = sumaPt + pt;
193
194 goodtrack++;
195
196 }//end track loop
197
198
199 if(sumaPt < 0.0000001){
200 return -1;
201 }
202
203 dir = (sqrt(pow(sumaPx,2)+pow(sumaPy,2)))/sumaPt;
204
205 return dir;
206}
207
208//__________________________________________________________________
209Double_t AliAnalysisCentralCutEvtMC::CalcDirUnit(AliMCEvent* const mcEvent) {
210// computes the directivity using only the SPD unity vectors
211
212 AliStack *stack=mcEvent->Stack();
213 if(!stack){
214 return -1;
215 }
216
217 Int_t nPrim = stack->GetNprimary();
218
219 Int_t goodtrack = 0;
220
221 Double_t eta;
222 Double_t pt;
223
224 Double_t dirU;
225 Double_t px,py,pxU,pyU;
226 Double_t sumaPxU = 0;
227 Double_t sumaPyU = 0;
228
229
230 for(Int_t i = 0; i < nPrim; i++){//track loop -> compute directivity
231
232 TParticle *particle=stack->Particle(i);
233 if(!particle){
234 continue;
235 }
236
237 eta = particle->Eta();
238
239 if((eta > 1.9)||(eta < 0.0)){ // half of the SPD coverage -> directivity
240 continue;
241 }
242
243 TParticlePDG *particlePDG = particle->GetPDG(0);
244 if(!particlePDG){
245 continue;
246 }
247
248 if(TMath::Abs(particlePDG->Charge()) < 3){
249 continue;
250 }
251
252 px = particle->Px();
253 py = particle->Py();
254 pt = particle->Pt();
255
256 pxU = px/pt;
257 pyU = py/pt;
258
259 sumaPxU = sumaPxU + pxU;
260 sumaPyU = sumaPyU + pyU;
261
262 goodtrack++;
263
264 }//end track loop
265
266 if(goodtrack == 0){
267 return -1;
268 }
269
270 dirU = (sqrt(pow(sumaPxU,2)+pow(sumaPyU,2)))/goodtrack;
271
272
273 return dirU;
274}