1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 //---------------------------------------------------------------------
18 // Event shape utility class
19 // Circularity, Thrust, ...
20 // Authors: Antonio Ortiz Velasquez <Antonio.Ortiz.Velasquez@cern.ch>
22 //---------------------------------------------------------------------
25 #include "AliEventShape.h"
29 #include "AliMCEventHandler.h"
30 #include "AliMCEvent.h"
32 #include <TMatrixDSym.h>
34 #include <TMatrixDSymEigen.h>
35 #include <TParticle.h>
36 #include <TParticlePDG.h>
41 //____________________________________________________________________
42 ClassImp(AliEventShape)
44 //___________________________________________________________________
45 TArrayD * AliEventShape::GetThrustParamMC(AliMCEvent* mcEvent, Int_t nstudymin, Double_t ptcutoff, Double_t etacutoff, Bool_t chom)
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();
67 stack = mcEvent->Stack();
75 TArrayD* evsh = new TArrayD(3);
76 Int_t nPrim = stack->GetNprimary();
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;
94 else{ //only charged particles
95 if (pdgPart->Charge() == 0)continue;
99 // Minimum number of particles used in the analysis
100 if(nmctracks < nstudymin){
108 pxT = new Double_t[nmctracks];
109 pyT = new Double_t[nmctracks];
110 ptT = new Double_t[nmctracks];
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();
127 if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16)continue;
129 if (pdgPart->Charge() == 0) continue;
141 Double_t numerador = 0;
142 Double_t numerador2 = 0;
143 Double_t phimax = -1;
147 Double_t thrust = 80;
148 Double_t thrustminor = 80;
151 Double_t phiparam = 0;
153 for(Int_t i = 0; i < 360; ++i){
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.
164 pFull=numerador / ptsuma;
165 if(pFull > pMax)//maximization of pFull
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})
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
183 evsh->AddAt(thrust, 0);
184 evsh->AddAt(thrustminor, 1);
185 evsh->AddAt(recoil, 2);
197 Double_t AliEventShape::GetCircularityMC(AliMCEvent* mcEvent, Int_t nstudymin, Double_t ptcutoff, Double_t etacutoff, Bool_t chom)
200 This function returns the circularity value of the event
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).
214 stack = mcEvent->Stack();
223 Double_t circularity = -2;
225 Int_t nPrim = stack->GetNprimary();
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();
241 // skip photons and neutrinos
242 if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16) continue;
245 if (pdgPart->Charge() == 0)continue;
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);
257 if (nmctracks < nstudymin) {
258 Printf("Too few particles, abort");
269 const TMatrixDSymEigen eigen(s);
270 const TVectorD eigenVal=eigen.GetEigenValues();
271 circularity = 2 * (1 - eigenVal(0));