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];
111 for (Int_t i = 0; i < nmctracks; i++)
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();
132 if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16)continue;
134 if (pdgPart->Charge() == 0) continue;
146 Double_t numerador = 0;
147 Double_t numerador2 = 0;
148 Double_t phimax = -1;
152 Double_t thrust = 80;
153 Double_t thrustminor = 80;
156 Double_t phiparam = 0;
158 for(Int_t i = 0; i < 360; ++i){
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.
169 pFull=numerador / ptsuma;
170 if(pFull > pMax)//maximization of pFull
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})
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
188 evsh->AddAt(thrust, 0);
189 evsh->AddAt(thrustminor, 1);
190 evsh->AddAt(recoil, 2);
201 Double_t AliEventShape::GetCircularityMC(AliMCEvent* mcEvent, Int_t nstudymin, Double_t ptcutoff, Double_t etacutoff, Bool_t chom)
204 This function returns the circularity value of the event
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).
218 stack = mcEvent->Stack();
227 Double_t circularity = -2;
229 Int_t nPrim = stack->GetNprimary();
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();
245 // skip photons and neutrinos
246 if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16) continue;
249 if (pdgPart->Charge() == 0)continue;
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);
261 if (nmctracks < nstudymin) {
262 Printf("Too few particles, stopping");
273 const TMatrixDSymEigen eigen(s);
274 const TVectorD eigenVal=eigen.GetEigenValues();
275 circularity = 2 * (1 - eigenVal(0));