Keep track of missing DCS points in DDL maps (flagged by 'x')
[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];
111
112 for (Int_t iMCTracks = 0; iMCTracks < nPrim; ++iMCTracks) {
113 TParticle* trackmc = stack->Particle(iMCTracks);
114 if (!trackmc) continue;
115 Double_t etamc = trackmc ->Eta();
116 Double_t pxmc = trackmc->Px();
117 Double_t pymc = trackmc->Py();
118 Double_t ptmc = trackmc->Pt();
119 Int_t pdgCode = TMath::Abs(trackmc->GetPdgCode());
120 if (TMath::Abs(etamc) > etacutoff) continue;
121 if(ptmc < ptcutoff) continue;
122 Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks);
123 if(isprimary==0) continue;
124 TParticlePDG* pdgPart =trackmc ->GetPDG();
125
126 if(chom==1){
127 if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16)continue;
128 } else {
129 if (pdgPart->Charge() == 0) continue;
130 }
131
132 ptT[j] = ptmc;
133 pxT[j] = pxmc;
134 pyT[j] = pymc;
135 ptsuma += ptmc;
136 pxsuma+=pxmc;
137 pysuma+=pymc;
138 j++;
139 }
140
141 Double_t numerador = 0;
142 Double_t numerador2 = 0;
143 Double_t phimax = -1;
144 Double_t pFull = -1;
145 Double_t pMax = 0;
146 Double_t phi = 0;
147 Double_t thrust = 80;
148 Double_t thrustminor = 80;
149 Double_t nx = 0;
150 Double_t ny = 0;
151 Double_t phiparam = 0;
152 //Getting thrust
153 for(Int_t i = 0; i < 360; ++i){
154 numerador = 0;
155 phiparam = 0;
156 nx = 0;
157 ny = 0;
158 phiparam=((TMath::Pi()) * i) / 180; // parametrization of the angle
159 nx = TMath::Cos(phiparam); // x component of an unitary vector n
160 ny = TMath::Sin(phiparam); // y component of an unitary vector n
161 for(Int_t i1 = 0; i1 < nmctracks; ++i1){
162 numerador += TMath::Abs(nx * pxT[i1] + ny * pyT[i1]);//product between momentum proyection in XY plane and the unitari vector.
163 }
164 pFull=numerador / ptsuma;
165 if(pFull > pMax)//maximization of pFull
166 {
167 pMax = pFull;
168 phi = phiparam;
169 }
170 }
171
172 phimax=(phi * 180) / TMath::Pi();//angular parameter of the unitary vector which maximiza thrust
173 //if n vector and beam axis form a plane, then we can calculate a second unitary vector perpendicular to that plane
174 Double_t nx1 = TMath::Cos(phi);
175 Double_t ny1 = TMath::Sin(phi);
176 for(Int_t i2 =0; i2 < nmctracks; ++i2){
177 numerador2 += TMath::Abs(pxT[i2] * ny1 - nx1 * pyT[i2]);//cross product: P_{i} X n, P_{i}=(px_{i},py_{i})
178 }
179 thrust = 1 - pMax;//this is the value of thrust
180 thrustminor = numerador2 / ptsuma;//this is the value of thrust minor
181 Double_t recoil = TMath::Abs(TMath::Sqrt(pxsuma * pxsuma + pysuma * pysuma)) / (ptsuma);//factor sentsitive to radiation outside from acceptance
182
183 evsh->AddAt(thrust, 0);
184 evsh->AddAt(thrustminor, 1);
185 evsh->AddAt(recoil, 2);
186
187
188 delete evsh;
189 delete ptT;
190 delete pxT;
191 delete pyT;
192
193 return evsh;
194}
195
196
197Double_t AliEventShape::GetCircularityMC(AliMCEvent* mcEvent, Int_t nstudymin, Double_t ptcutoff, Double_t etacutoff, Bool_t chom)
198{
199 /*
200 This function returns the circularity value of the event
201
202 The calculus uses primary particles. The input parameters:
203 1. nstudymin, is the minumum number of particles which you want that participate in the calculus.(default:nstudymin =3)
204 2. ptcutoff, is the cut in pt applied to participants to calculate the variables.(default: ptcutoff=1)
205 3. etacutoff, is the cut in acceptance applied to participants to calculate the variables.(default: etacutoff=1)
206 4. if chom=kTRUE, then the calculus includes neutral particles (rejecting photons and neutrinos).
207 if chom=kFALSE, then the calculus includes only charged particles (rejecting photons and neutrinos).
208 1>=circularity>=0
209 */
210
211
212 AliStack* stack = 0;
213
214 stack = mcEvent->Stack();
215
216 TMatrixDSym s(2);
217
218 Double_t s00 = 0;
219 Double_t s01 = 0;
220 Double_t s10 = 0;
221 Double_t s11 = 0;
222 Double_t ptot = 0;
223 Double_t circularity = -2;
224 Int_t nmctracks = 0;
225 Int_t nPrim = stack->GetNprimary();
226
227 for (Int_t iMCTracks = 0; iMCTracks < nPrim; ++iMCTracks) {
228 TParticle* trackmc = stack->Particle(iMCTracks);
229 if (!trackmc) continue;
230 Double_t etamc = trackmc ->Eta();
231 Double_t ptmc = trackmc->Pt();
232 Double_t pxmc = trackmc->Px();
233 Double_t pymc = trackmc->Py();
234 Int_t pdgCode = TMath::Abs(trackmc->GetPdgCode());
235 if (TMath::Abs(etamc) > etacutoff) continue;
236 if (ptmc < ptcutoff) continue;
237 Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks);
238 if (isprimary == 0) continue;
239 TParticlePDG* pdgPart = trackmc ->GetPDG();
240 if(chom == kTRUE){
241 // skip photons and neutrinos
242 if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16) continue;
243 }
244 else{
245 if (pdgPart->Charge() == 0)continue;
246 }
247
248 ptot = ptot + (ptmc * ptmc);
249 s00 = s00 + (pxmc * pxmc);
250 s01 = s01 + (pxmc * pymc);
251 s10 = s10 + (pymc * pxmc);
252 s11 = s11 + (pymc * pymc);
253 nmctracks++;
254 } //track loop
255
256
257 if (nmctracks < nstudymin) {
258 Printf("Too few particles, abort");
259 return -2;
260 }
261
262
263
264 if(ptot != 0){
265 s(0,0) = s00 / ptot;
266 s(0,1) = s01 / ptot;
267 s(1,0) = s10 / ptot;
268 s(1,1) = s11 / ptot;
269 const TMatrixDSymEigen eigen(s);
270 const TVectorD eigenVal=eigen.GetEigenValues();
271 circularity = 2 * (1 - eigenVal(0));
272 }
273 return circularity;
274}
275
276
277
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