[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 i = 0; i < nmctracks; i++)
    {
      pxT[i] = 0;
      pyT[i] = 0;
      ptT[i] = 0;
    }
  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, stopping");
      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
	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	for (Int_t i = 0; i < nmctracks; i++)
	112	{
	113	pxT[i] = 0;
	114	pyT[i] = 0;
	115	ptT[i] = 0;
	116	}
	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();
	130
	131	if(chom==1){
	132	if (pdgCode == 22 \|\| pdgCode == 12 \|\| pdgCode == 14 \|\| pdgCode == 16)continue;
	133	} else {
	134	if (pdgPart->Charge() == 0) continue;
	135	}
	136
	137	ptT[j] = ptmc;
	138	pxT[j] = pxmc;
	139	pyT[j] = pymc;
	140	ptsuma += ptmc;
	141	pxsuma+=pxmc;
	142	pysuma+=pymc;
	143	j++;
	144	}
	145
	146	Double_t numerador = 0;
	147	Double_t numerador2 = 0;
	148	Double_t phimax = -1;
	149	Double_t pFull = -1;
	150	Double_t pMax = 0;
	151	Double_t phi = 0;
	152	Double_t thrust = 80;
	153	Double_t thrustminor = 80;
	154	Double_t nx = 0;
	155	Double_t ny = 0;
	156	Double_t phiparam = 0;
	157	//Getting thrust
	158	for(Int_t i = 0; i < 360; ++i){
	159	numerador = 0;
	160	phiparam = 0;
	161	nx = 0;
	162	ny = 0;
	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.
	168	}
	169	pFull=numerador / ptsuma;
	170	if(pFull > pMax)//maximization of pFull
	171	{
	172	pMax = pFull;
	173	phi = phiparam;
	174	}
	175	}
	176
	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})
	183	}
	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
	187
	188	evsh->AddAt(thrust, 0);
	189	evsh->AddAt(thrustminor, 1);
	190	evsh->AddAt(recoil, 2);
	191
	192
	193	delete [] ptT;
	194	delete [] pxT;
	195	delete [] pyT;
	196
	197	return evsh;
	198	}
	199
	200
	201	Double_t AliEventShape::GetCircularityMC(AliMCEvent* mcEvent, Int_t nstudymin, Double_t ptcutoff, Double_t etacutoff, Bool_t chom)
	202	{
	203	/*
	204	This function returns the circularity value of the event
	205
	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).
	212	1>=circularity>=0
	213	*/
	214
	215
	216	AliStack* stack = 0;
	217
	218	stack = mcEvent->Stack();
	219
	220	TMatrixDSym s(2);
	221
	222	Double_t s00 = 0;
	223	Double_t s01 = 0;
	224	Double_t s10 = 0;
	225	Double_t s11 = 0;
	226	Double_t ptot = 0;
	227	Double_t circularity = -2;
	228	Int_t nmctracks = 0;
	229	Int_t nPrim = stack->GetNprimary();
	230
	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();
	244	if(chom == kTRUE){
	245	// skip photons and neutrinos
	246	if (pdgCode == 22 \|\| pdgCode == 12 \|\| pdgCode == 14 \|\| pdgCode == 16) continue;
	247	}
	248	else{
	249	if (pdgPart->Charge() == 0)continue;
	250	}
	251
	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);
	257	nmctracks++;
	258	} //track loop
	259
	260
	261	if (nmctracks < nstudymin) {
	262	Printf("Too few particles, stopping");
	263	return -2;
	264	}
	265
	266
	267
	268	if(ptot != 0){
	269	s(0,0) = s00 / ptot;
	270	s(0,1) = s01 / ptot;
	271	s(1,0) = s10 / ptot;
	272	s(1,1) = s11 / ptot;
	273	const TMatrixDSymEigen eigen(s);
	274	const TVectorD eigenVal=eigen.GetEigenValues();
	275	circularity = 2 * (1 - eigenVal(0));
	276	}
	277	return circularity;
	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
	308
	309
	310
	311