Compatibility with ROOT trunk
[u/mrichter/AliRoot.git] / JETAN / AliEventShape.cxx
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 //____________________________________________________________________
42 ClassImp(AliEventShape)
43
44 //___________________________________________________________________
45 TArrayD * 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   for (Int_t i = 0; i < nmctracks; i++)
112     {
113       pxT[i] = 0;
114       pyT[i] = 0;
115       ptT[i] = 0;
116     }
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
193   delete [] ptT;
194   delete [] pxT;
195   delete [] pyT;
196
197   return evsh;  
198 }
199
200
201 Double_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) {
262       Printf("Too few particles, stopping");
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