[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 [] 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
4ce766eb	188	delete [] ptT;
	189	delete [] pxT;
	190	delete [] pyT;
d2f1e153	191
	192	return evsh;
	193	}
	194
	195
	196	Double_t AliEventShape::GetCircularityMC(AliMCEvent* mcEvent, Int_t nstudymin, Double_t ptcutoff, Double_t etacutoff, Bool_t chom)
	197	{
	198	/*
	199	This function returns the circularity value of the event
	200
	201	The calculus uses primary particles. The input parameters:
	202	1. nstudymin, is the minumum number of particles which you want that participate in the calculus.(default:nstudymin =3)
	203	2. ptcutoff, is the cut in pt applied to participants to calculate the variables.(default: ptcutoff=1)
	204	3. etacutoff, is the cut in acceptance applied to participants to calculate the variables.(default: etacutoff=1)
	205	4. if chom=kTRUE, then the calculus includes neutral particles (rejecting photons and neutrinos).
	206	if chom=kFALSE, then the calculus includes only charged particles (rejecting photons and neutrinos).
	207	1>=circularity>=0
	208	*/
	209
	210
	211	AliStack* stack = 0;
	212
	213	stack = mcEvent->Stack();
	214
	215	TMatrixDSym s(2);
	216
	217	Double_t s00 = 0;
	218	Double_t s01 = 0;
	219	Double_t s10 = 0;
	220	Double_t s11 = 0;
	221	Double_t ptot = 0;
	222	Double_t circularity = -2;
	223	Int_t nmctracks = 0;
	224	Int_t nPrim = stack->GetNprimary();
	225
	226	for (Int_t iMCTracks = 0; iMCTracks < nPrim; ++iMCTracks) {
	227	TParticle* trackmc = stack->Particle(iMCTracks);
	228	if (!trackmc) continue;
	229	Double_t etamc = trackmc ->Eta();
	230	Double_t ptmc = trackmc->Pt();
	231	Double_t pxmc = trackmc->Px();
	232	Double_t pymc = trackmc->Py();
	233	Int_t pdgCode = TMath::Abs(trackmc->GetPdgCode());
	234	if (TMath::Abs(etamc) > etacutoff) continue;
	235	if (ptmc < ptcutoff) continue;
	236	Bool_t isprimary = stack->IsPhysicalPrimary(iMCTracks);
	237	if (isprimary == 0) continue;
	238	TParticlePDG* pdgPart = trackmc ->GetPDG();
	239	if(chom == kTRUE){
	240	// skip photons and neutrinos
	241	if (pdgCode == 22 \|\| pdgCode == 12 \|\| pdgCode == 14 \|\| pdgCode == 16) continue;
	242	}
	243	else{
	244	if (pdgPart->Charge() == 0)continue;
	245	}
	246
	247	ptot = ptot + (ptmc * ptmc);
	248	s00 = s00 + (pxmc * pxmc);
	249	s01 = s01 + (pxmc * pymc);
	250	s10 = s10 + (pymc * pxmc);
	251	s11 = s11 + (pymc * pymc);
	252	nmctracks++;
	253	} //track loop
	254
255
256	if (nmctracks < nstudymin) {
257	Printf("Too few particles, abort");
258	return -2;
259	}
260
261
262
263	if(ptot != 0){
264	s(0,0) = s00 / ptot;
265	s(0,1) = s01 / ptot;
266	s(1,0) = s10 / ptot;
267	s(1,1) = s11 / ptot;
268	const TMatrixDSymEigen eigen(s);
269	const TVectorD eigenVal=eigen.GetEigenValues();
270	circularity = 2 * (1 - eigenVal(0));
271	}
272	return circularity;
273	}
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