Coverity fix
[u/mrichter/AliRoot.git] / JETAN / AliEventShape.cxx
CommitLineData
d2f1e153 1/**************************************************************************
2 * Copyright(c) 1998-1999, 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
16/* $Id$ */
17//---------------------------------------------------------------------
18// Event shape utility class
19// Circularity, Thrust, ...
20// Authors: Antonio Ortiz Velasquez <Antonio.Ortiz.Velasquez@cern.ch>
21//
22//---------------------------------------------------------------------
23
24
25#include "AliEventShape.h"
26
27#include "AliStack.h"
28#include "AliLog.h"
29#include "AliMCEventHandler.h"
30#include "AliMCEvent.h"
31
32#include <TMatrixDSym.h>
33#include <TVectorD.h>
34#include <TMatrixDSymEigen.h>
35#include <TParticle.h>
36#include <TParticlePDG.h>
37
38
39
40
41//____________________________________________________________________
42ClassImp(AliEventShape)
43
44//___________________________________________________________________
45TArrayD * AliEventShape::GetThrustParamMC(AliMCEvent* mcEvent, Int_t nstudymin, Double_t ptcutoff, Double_t etacutoff, Bool_t chom)
46{
47 /*
48 This function returns an array of values of thrust. To get these values you have to do:
49 TArrayD* eventshapes = 0;
50 eventshapes = AliShapeEvent::GetThrustParamMC(mcEvent, 3, 1, 1, kFALSE);
51 Double_t thrust=eventshapes->GetAt(0);
52 Double_t thrustmin=eventshapes->GetAt(1);
53 Double_t recoil=eventshapes->GetAt(2);
54 The calculus uses primary particles. The input parameters:
55 1. nstudymin, is the minumum number of particles which you want that participate in the calculus.(default:nstudymin =3)
56 2. ptcutoff, is the cut in pt applied to participants to calculate the variables.(default: ptcutoff=1)
57 3. etacutoff, is the cut in acceptance applied to participants to calculate the variables.(default: etacutoff=1)
58 4. if chom=kTRUE, then the calculus includes neutral particles (rejecting photons and neutrinos).
59 if chom=kFALSE, then the calculus includes only charged particles (rejecting photons and neutrinos).
60 Returned values: thrust->0: 2-jet event, thrust->0.5: isotropic event
61 Recoil is a term which is sensitive to radiation outside from acceptance, 1>=recoil>=0,
62 thrustmin, is a measure of the radiation which is perpendicular to the plane formed by beam axis and thrust axis, 2/TMath::Pi()>thrustmin>0. In the limit of 2 back-to-back jets thrusmin->0, while in the case of a uniformly distributed event thrustmin->2/TMath::Pi();
63 */
64
65 AliStack* stack = 0;
66
67 stack = mcEvent->Stack();
68 Double_t * ptT = 0;
69 Double_t * pxT = 0;
70 Double_t * pyT = 0;
71 Double_t ptsuma = 0;
72 Double_t pxsuma = 0;
73 Double_t pysuma = 0;
74
75 TArrayD* evsh = new TArrayD(3);
76 Int_t nPrim = stack->GetNprimary();
77 Int_t nmctracks = 0;
78 for (Int_t iMCTracks = 0; iMCTracks < nPrim; iMCTracks++) {
79 TParticle* trackmc = stack->Particle(iMCTracks);
80 if (!trackmc) continue;
81 Double_t etamc =trackmc ->Eta();
82 Double_t ptmc=trackmc->Pt();
83 Int_t pdgCode = TMath::Abs(trackmc->GetPdgCode());
84 if (TMath::Abs(etamc) > etacutoff) continue; //only particles in |eta|<=etacutoff
85 if(ptmc < ptcutoff) continue; // PT cut
86 Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks); // Check if particle is charged, and primary
87 if(isprimary == 0) continue; // only primary particles
88 TParticlePDG* pdgPart =trackmc ->GetPDG();
89 if(chom == 1){//include neutral particles
90 // skip photons and neutrinos
91 if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16) continue;
92 nmctracks++;
93 }
94 else{ //only charged particles
95 if (pdgPart->Charge() == 0)continue;
96 nmctracks++;
97 }
98 }
99 // Minimum number of particles used in the analysis
100 if(nmctracks < nstudymin){
101 evsh->AddAt(-2,0);
102 evsh->AddAt(-2,1);
103 evsh->AddAt(-2,2);
104 return evsh;
105 }
106
107 Int_t j=0;
108 pxT = new Double_t[nmctracks];
109 pyT = new Double_t[nmctracks];
110 ptT = new Double_t[nmctracks];
1dfd2819 111 for (Int_t i = 0; i < nmctracks; i++)
112 {
113 pxT[i] = 0;
114 pyT[i] = 0;
115 ptT[i] = 0;
116 }
d2f1e153 117 for (Int_t iMCTracks = 0; iMCTracks < nPrim; ++iMCTracks) {
118 TParticle* trackmc = stack->Particle(iMCTracks);
119 if (!trackmc) continue;
120 Double_t etamc = trackmc ->Eta();
121 Double_t pxmc = trackmc->Px();
122 Double_t pymc = trackmc->Py();
123 Double_t ptmc = trackmc->Pt();
124 Int_t pdgCode = TMath::Abs(trackmc->GetPdgCode());
125 if (TMath::Abs(etamc) > etacutoff) continue;
126 if(ptmc < ptcutoff) continue;
127 Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks);
128 if(isprimary==0) continue;
129 TParticlePDG* pdgPart =trackmc ->GetPDG();
130
131 if(chom==1){
132 if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16)continue;
133 } else {
134 if (pdgPart->Charge() == 0) continue;
135 }
136
137 ptT[j] = ptmc;
138 pxT[j] = pxmc;
139 pyT[j] = pymc;
140 ptsuma += ptmc;
141 pxsuma+=pxmc;
142 pysuma+=pymc;
143 j++;
144 }
145
146 Double_t numerador = 0;
147 Double_t numerador2 = 0;
148 Double_t phimax = -1;
149 Double_t pFull = -1;
150 Double_t pMax = 0;
151 Double_t phi = 0;
152 Double_t thrust = 80;
153 Double_t thrustminor = 80;
154 Double_t nx = 0;
155 Double_t ny = 0;
156 Double_t phiparam = 0;
157 //Getting thrust
158 for(Int_t i = 0; i < 360; ++i){
159 numerador = 0;
160 phiparam = 0;
161 nx = 0;
162 ny = 0;
163 phiparam=((TMath::Pi()) * i) / 180; // parametrization of the angle
164 nx = TMath::Cos(phiparam); // x component of an unitary vector n
165 ny = TMath::Sin(phiparam); // y component of an unitary vector n
166 for(Int_t i1 = 0; i1 < nmctracks; ++i1){
167 numerador += TMath::Abs(nx * pxT[i1] + ny * pyT[i1]);//product between momentum proyection in XY plane and the unitari vector.
168 }
169 pFull=numerador / ptsuma;
170 if(pFull > pMax)//maximization of pFull
171 {
172 pMax = pFull;
173 phi = phiparam;
174 }
175 }
176
177 phimax=(phi * 180) / TMath::Pi();//angular parameter of the unitary vector which maximiza thrust
178 //if n vector and beam axis form a plane, then we can calculate a second unitary vector perpendicular to that plane
179 Double_t nx1 = TMath::Cos(phi);
180 Double_t ny1 = TMath::Sin(phi);
181 for(Int_t i2 =0; i2 < nmctracks; ++i2){
182 numerador2 += TMath::Abs(pxT[i2] * ny1 - nx1 * pyT[i2]);//cross product: P_{i} X n, P_{i}=(px_{i},py_{i})
183 }
184 thrust = 1 - pMax;//this is the value of thrust
185 thrustminor = numerador2 / ptsuma;//this is the value of thrust minor
186 Double_t recoil = TMath::Abs(TMath::Sqrt(pxsuma * pxsuma + pysuma * pysuma)) / (ptsuma);//factor sentsitive to radiation outside from acceptance
187
188 evsh->AddAt(thrust, 0);
189 evsh->AddAt(thrustminor, 1);
190 evsh->AddAt(recoil, 2);
191
192
4ce766eb 193 delete [] ptT;
194 delete [] pxT;
195 delete [] pyT;
d2f1e153 196
197 return evsh;
198}
199
200
201Double_t AliEventShape::GetCircularityMC(AliMCEvent* mcEvent, Int_t nstudymin, Double_t ptcutoff, Double_t etacutoff, Bool_t chom)
202{
203 /*
204 This function returns the circularity value of the event
205
206 The calculus uses primary particles. The input parameters:
207 1. nstudymin, is the minumum number of particles which you want that participate in the calculus.(default:nstudymin =3)
208 2. ptcutoff, is the cut in pt applied to participants to calculate the variables.(default: ptcutoff=1)
209 3. etacutoff, is the cut in acceptance applied to participants to calculate the variables.(default: etacutoff=1)
210 4. if chom=kTRUE, then the calculus includes neutral particles (rejecting photons and neutrinos).
211 if chom=kFALSE, then the calculus includes only charged particles (rejecting photons and neutrinos).
212 1>=circularity>=0
213 */
214
215
216 AliStack* stack = 0;
217
218 stack = mcEvent->Stack();
219
220 TMatrixDSym s(2);
221
222 Double_t s00 = 0;
223 Double_t s01 = 0;
224 Double_t s10 = 0;
225 Double_t s11 = 0;
226 Double_t ptot = 0;
227 Double_t circularity = -2;
228 Int_t nmctracks = 0;
229 Int_t nPrim = stack->GetNprimary();
230
231 for (Int_t iMCTracks = 0; iMCTracks < nPrim; ++iMCTracks) {
232 TParticle* trackmc = stack->Particle(iMCTracks);
233 if (!trackmc) continue;
234 Double_t etamc = trackmc ->Eta();
235 Double_t ptmc = trackmc->Pt();
236 Double_t pxmc = trackmc->Px();
237 Double_t pymc = trackmc->Py();
238 Int_t pdgCode = TMath::Abs(trackmc->GetPdgCode());
239 if (TMath::Abs(etamc) > etacutoff) continue;
240 if (ptmc < ptcutoff) continue;
241 Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks);
242 if (isprimary == 0) continue;
243 TParticlePDG* pdgPart = trackmc ->GetPDG();
244 if(chom == kTRUE){
245 // skip photons and neutrinos
246 if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16) continue;
247 }
248 else{
249 if (pdgPart->Charge() == 0)continue;
250 }
251
252 ptot = ptot + (ptmc * ptmc);
253 s00 = s00 + (pxmc * pxmc);
254 s01 = s01 + (pxmc * pymc);
255 s10 = s10 + (pymc * pxmc);
256 s11 = s11 + (pymc * pymc);
257 nmctracks++;
258 } //track loop
259
260
261 if (nmctracks < nstudymin) {
1240edf5 262 Printf("Too few particles, stopping");
d2f1e153 263 return -2;
264 }
265
266
267
268 if(ptot != 0){
269 s(0,0) = s00 / ptot;
270 s(0,1) = s01 / ptot;
271 s(1,0) = s10 / ptot;
272 s(1,1) = s11 / ptot;
273 const TMatrixDSymEigen eigen(s);
274 const TVectorD eigenVal=eigen.GetEigenValues();
275 circularity = 2 * (1 - eigenVal(0));
276 }
277 return circularity;
278}
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311