[u/mrichter/AliRoot.git] / JETAN / AliEventShape.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
   
/* $Id$ */
//---------------------------------------------------------------------
// Event shape utility class
// Circularity, Thrust, ... 
// Authors: Antonio Ortiz Velasquez <Antonio.Ortiz.Velasquez@cern.ch>
//          
//---------------------------------------------------------------------


#include "AliEventShape.h"

#include "AliStack.h"
#include "AliLog.h"
#include "AliMCEventHandler.h"
#include "AliMCEvent.h"

#include <TMatrixDSym.h>
#include <TVectorD.h>
#include <TMatrixDSymEigen.h>
#include <TParticle.h>
#include <TParticlePDG.h>


//____________________________________________________________________
ClassImp(AliEventShape)

//___________________________________________________________________
TArrayD * AliEventShape::GetThrustParamMC(AliMCEvent* mcEvent, Int_t  nstudymin, Double_t ptcutoff, Double_t etacutoff, Bool_t chom)
{
  /*
    This function returns an array of values of thrust. To get these values you have to do:
    TArrayD* eventshapes = 0;
    eventshapes = AliShapeEvent::GetThrustParamMC(mcEvent, 3, 1, 1, kFALSE);
    Double_t thrust=eventshapes->GetAt(0);
    Double_t thrustmin=eventshapes->GetAt(1);
    Double_t recoil=eventshapes->GetAt(2);  
    The calculus uses  primary  particles. The input parameters:
    1. nstudymin, is the minumum number of particles which you want that participate in the calculus.(default:nstudymin =3)
    2. ptcutoff, is the cut in pt applied to participants to calculate the variables.(default: ptcutoff=1)
    3. etacutoff, is the cut in acceptance applied to participants to calculate the variables.(default: etacutoff=1)
    4. if chom=kTRUE, then the calculus includes neutral particles (rejecting photons and neutrinos).
       if chom=kFALSE, then the calculus includes only charged particles (rejecting photons and neutrinos).
       Returned values: thrust->0: 2-jet event, thrust->0.5: isotropic event
       Recoil is a term which is sensitive to radiation outside from acceptance, 1>=recoil>=0,
       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();
  */
  
  AliStack* stack = 0;

  stack = mcEvent->Stack();
  Double_t * ptT = 0;
  Double_t * pxT = 0;
  Double_t * pyT = 0;
  Double_t ptsuma = 0;
  Double_t pxsuma = 0;
  Double_t pysuma = 0;

  TArrayD* evsh = new TArrayD(3);  
  Int_t nPrim  = stack->GetNprimary();
  Int_t  nmctracks = 0;  
  for (Int_t iMCTracks = 0; iMCTracks < nPrim; iMCTracks++) {    
      TParticle* trackmc = stack->Particle(iMCTracks);
      if (!trackmc) continue;
      Double_t etamc =trackmc ->Eta();
      Double_t ptmc=trackmc->Pt();
      Int_t pdgCode = TMath::Abs(trackmc->GetPdgCode());
      if (TMath::Abs(etamc) > etacutoff) continue; //only particles in |eta|<=etacutoff
      if(ptmc < ptcutoff) continue;  // PT cut
      Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks);    // Check if particle is charged, and primary
      if(isprimary == 0) continue;  // only primary particles
      TParticlePDG* pdgPart =trackmc ->GetPDG();
      if(chom == 1){//include neutral particles
	  // skip photons and neutrinos
	  if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16) continue;
	  nmctracks++;  
      }
      else{ //only charged particles
	  if (pdgPart->Charge() == 0)continue;
	  nmctracks++;     
      }
  }
  // Minimum number of particles used in the analysis
  if(nmctracks < nstudymin){
      evsh->AddAt(-2,0);
      evsh->AddAt(-2,1);
      evsh->AddAt(-2,2);
      return evsh;
  }
  
  Int_t j=0;
  pxT = new Double_t[nmctracks];
  pyT = new Double_t[nmctracks];
  ptT = new Double_t[nmctracks];

  for (Int_t iMCTracks = 0; iMCTracks < nPrim; ++iMCTracks) {    
      TParticle* trackmc = stack->Particle(iMCTracks);
      if (!trackmc) continue;
      Double_t etamc = trackmc ->Eta();
      Double_t pxmc  = trackmc->Px();
      Double_t pymc  = trackmc->Py();
      Double_t ptmc  = trackmc->Pt();
      Int_t pdgCode  = TMath::Abs(trackmc->GetPdgCode());
      if (TMath::Abs(etamc) > etacutoff) continue;
      if(ptmc < ptcutoff) continue;
      Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks); 
      if(isprimary==0) continue;
      TParticlePDG* pdgPart =trackmc ->GetPDG();

      if(chom==1){
	  if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16)continue;
      } else {
	  if (pdgPart->Charge() == 0) continue;
      }
      
      ptT[j] = ptmc;
      pxT[j] = pxmc;
      pyT[j] = pymc;
      ptsuma += ptmc;
      pxsuma+=pxmc;
      pysuma+=pymc;
      j++;    
  }

  Double_t numerador = 0;
  Double_t numerador2 = 0;
  Double_t phimax = -1;  
  Double_t pFull = -1;
  Double_t pMax = 0;
  Double_t phi = 0;
  Double_t thrust = 80;
  Double_t thrustminor = 80;
  Double_t nx = 0;
  Double_t ny = 0;
  Double_t phiparam = 0;
  //Getting thrust
  for(Int_t i = 0; i < 360; ++i){
      numerador = 0;
      phiparam  = 0;
      nx = 0;
      ny = 0;
      phiparam=((TMath::Pi()) * i) / 180; // parametrization of the angle
      nx = TMath::Cos(phiparam);            // x component of an unitary vector n
      ny = TMath::Sin(phiparam);            // y component of an unitary vector n
      for(Int_t i1 = 0; i1 < nmctracks; ++i1){
	  numerador += TMath::Abs(nx * pxT[i1] + ny * pyT[i1]);//product between momentum proyection in XY plane and the unitari vector.
      }
      pFull=numerador / ptsuma;
      if(pFull > pMax)//maximization of pFull
      {
	  pMax = pFull;
	  phi = phiparam;
      }
  }

  phimax=(phi * 180) / TMath::Pi();//angular parameter of the unitary vector which maximiza thrust
  //if n vector and beam axis form a plane, then we can calculate a second unitary vector perpendicular to that plane
  Double_t nx1 = TMath::Cos(phi);
  Double_t ny1 = TMath::Sin(phi);
  for(Int_t i2 =0; i2 < nmctracks; ++i2){
      numerador2 += TMath::Abs(pxT[i2] * ny1 - nx1 * pyT[i2]);//cross product: P_{i} X n, P_{i}=(px_{i},py_{i})
  }
  thrust = 1 - pMax;//this is the value of thrust
  thrustminor = numerador2 / ptsuma;//this is the value of thrust minor
  Double_t recoil = TMath::Abs(TMath::Sqrt(pxsuma * pxsuma + pysuma * pysuma)) / (ptsuma);//factor sentsitive to radiation outside from acceptance 

  evsh->AddAt(thrust, 0);
  evsh->AddAt(thrustminor, 1);
  evsh->AddAt(recoil, 2);


  delete evsh;
  delete ptT;
  delete pxT;
  delete pyT;

  return evsh;  
}


Double_t AliEventShape::GetCircularityMC(AliMCEvent* mcEvent, Int_t  nstudymin, Double_t ptcutoff, Double_t etacutoff, Bool_t chom)
{
  /*
    This function returns the circularity value of the event 

    The calculus uses  primary  particles. The input parameters:
    1. nstudymin, is the minumum number of particles which you want that participate in the calculus.(default:nstudymin =3)
    2. ptcutoff, is the cut in pt applied to participants to calculate the variables.(default: ptcutoff=1)
    3. etacutoff, is the cut in acceptance applied to participants to calculate the variables.(default: etacutoff=1)
    4. if chom=kTRUE, then the calculus includes neutral particles (rejecting photons and neutrinos).
       if chom=kFALSE, then the calculus includes only charged particles (rejecting photons and neutrinos).
       1>=circularity>=0
  */


  AliStack* stack = 0;

  stack = mcEvent->Stack();

  TMatrixDSym s(2);

  Double_t s00 =  0;
  Double_t s01 =  0;
  Double_t s10 =  0;
  Double_t s11 =  0;
  Double_t ptot = 0;
  Double_t circularity = -2;
  Int_t  nmctracks = 0;  
  Int_t nPrim  = stack->GetNprimary();

  for (Int_t iMCTracks = 0; iMCTracks < nPrim; ++iMCTracks) {
    TParticle* trackmc = stack->Particle(iMCTracks);
    if (!trackmc) continue;
    Double_t etamc = trackmc ->Eta();
    Double_t ptmc  = trackmc->Pt();
    Double_t pxmc  = trackmc->Px();
    Double_t pymc  = trackmc->Py();
    Int_t pdgCode = TMath::Abs(trackmc->GetPdgCode());
    if (TMath::Abs(etamc) > etacutoff) continue;
    if (ptmc < ptcutoff) continue;
    Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks);
    if (isprimary == 0) continue;
    TParticlePDG* pdgPart = trackmc ->GetPDG();
    if(chom == kTRUE){
      // skip photons and neutrinos
      if (pdgCode == 22 || pdgCode == 12 || pdgCode == 14 || pdgCode == 16) continue; 
    }
    else{
      if (pdgPart->Charge() == 0)continue;
    }

    ptot = ptot + (ptmc * ptmc);
    s00 = s00 + (pxmc * pxmc);
    s01 = s01 + (pxmc * pymc);
    s10 = s10 + (pymc * pxmc);
    s11 = s11 + (pymc * pymc); 
    nmctracks++;
  } //track loop 


  if (nmctracks < nstudymin) {
      Printf("Too few particles, abort");
      return -2;
  }


  if(ptot != 0){
    s(0,0) = s00 / ptot;
    s(0,1) = s01 / ptot;
    s(1,0) = s10 / ptot;
    s(1,1) = s11 / ptot;
    const TMatrixDSymEigen eigen(s);
    const TVectorD eigenVal=eigen.GetEigenValues();
    circularity = 2 * (1 - eigenVal(0));
  }
  return circularity;
}
Commit	Line	Data
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	//____________________________________________________________________
	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