Event Shape Utility Class (Antonio Ortiz Velasquez <Antonio.Ortiz.Velasquez@cern...
[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
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
197 Double_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