[u/mrichter/AliRoot.git] / PWG4 / AliGammaMCReader.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$ */

/* History of cvs commits:
 *
 * $Log$
 * Revision 1.4  2007/10/29 13:48:42  gustavo
 * Corrected coding violations
 *
 *
 */

//_________________________________________________________________________
// Class for reading data (Kinematics) in order to do prompt gamma correlations
//  Class created from old AliPHOSGammaJet 
//  (see AliRoot versions previous Release 4-09)
//
//*-- Author: Gustavo Conesa (LNF-INFN) 
//////////////////////////////////////////////////////////////////////////////


// --- ROOT system ---

#include <TFile.h>
#include <TParticle.h>
#include <TH2.h>
#include <TChain.h>
#include <TRandom.h>
#include <TArrayI.h>
#include <TClonesArray.h>

//---- ANALYSIS system ----
#include "AliGammaMCReader.h" 
#include "Riostream.h"
#include "AliLog.h"
#include "AliStack.h"

ClassImp(AliGammaMCReader)

//____________________________________________________________________________
AliGammaMCReader::AliGammaMCReader() : 
  AliGammaReader(), fDecayPi0(0), fCheckOverlapping(kFALSE),  fNeutralParticlesArray(0x0)
{
  //Ctor
  
  //Initialize parameters
  InitParameters();
  fDataType = kMC;  
  
}

//____________________________________________________________________________
AliGammaMCReader::AliGammaMCReader(const AliGammaMCReader & g) :   
  AliGammaReader(g), fDecayPi0(g.fDecayPi0), fCheckOverlapping(g.fCheckOverlapping),
  fNeutralParticlesArray(g.fNeutralParticlesArray?new TArrayI(*g.fNeutralParticlesArray):0x0)
{
  // cpy ctor
}

//_________________________________________________________________________
AliGammaMCReader & AliGammaMCReader::operator = (const AliGammaMCReader & source)
{
  // assignment operator

  if(&source == this) return *this;

  fDecayPi0 = source.fDecayPi0; 
  fCheckOverlapping = source.fCheckOverlapping;
  delete fNeutralParticlesArray;
  fNeutralParticlesArray = source.fNeutralParticlesArray?new TArrayI(*source.fNeutralParticlesArray):0x0;
  return *this;

}

//_________________________________
AliGammaMCReader::~AliGammaMCReader() {
  //Dtor

  delete fNeutralParticlesArray ;

}

//___________________________________________________________________________
void AliGammaMCReader::CaseDecayGamma(Int_t index, TParticle * particle, AliStack * sta,
				  TClonesArray * plEMCAL, Int_t &indexEMCAL,
				  TClonesArray * plPHOS, Int_t &indexPHOS){
  //In case pi0 are decayed by pythia, check if mother is pi0 and in such case look if 
  //there is overlapping. Send particle=pi0 if there is overlapping. 

  TParticle * pmother =sta->Particle(particle->GetFirstMother());
  if(pmother->GetPdgCode() == 111 && pmother->GetNDaughters() == 2) {//Do not consider 3 particle decay case
    Int_t idaug0 = pmother->GetDaughter(0);
    Int_t idaug1 = pmother->GetDaughter(1);
    TParticle * pdaug0 = sta -> Particle(idaug0);
    TParticle * pdaug1 = sta -> Particle(idaug1);

    if((index ==  idaug0 &&  pdaug0->Pt() > fNeutralPtCut) ||
       (index ==  idaug1 &&  pdaug0->Pt() <= fNeutralPtCut)){//Check decay when first daughter arrives, do nothing with second.

      FillListWithDecayGammaOrPi0(pmother, pdaug0, pdaug1, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
    
    }
  }//mother is a pi0 with 2 daughters
  else{
    
    if(IsInPHOS(particle->Phi(),particle->Eta())) {
      TParticle * pmother =sta->Particle(particle->GetFirstMother());	   
      SetPhotonStatus(particle,pmother);     
      new((*plPHOS)[indexPHOS++])  TParticle(*particle) ;
    }
    else if(IsInEMCAL(particle->Phi(),particle->Eta())){
      TParticle * pmother =sta->Particle(particle->GetFirstMother());	   
      SetPhotonStatus(particle,pmother);
      new((*plEMCAL)[indexEMCAL++])  TParticle(*particle) ;
    }
  } 
}

//___________________________________________________________________________
 void AliGammaMCReader::CaseGeantDecay(TParticle * particle, AliStack * stack,
				   TClonesArray * plEMCAL, Int_t &indexEMCAL,
				   TClonesArray * plPHOS, Int_t &indexPHOS){
   //Find decay gamma from pi0, decayed by GEANT and put them in the list.
   
   Int_t ndaug = particle->GetNDaughters() ;
   TParticle * d1 = new TParticle();
   TParticle * d2 = new TParticle();
   
   if(ndaug > 0){//At least 1 daugther
     d1 = stack->Particle(particle->GetDaughter(0));
     if (ndaug > 1 ) 
       d2 = stack->Particle(particle->GetDaughter(1));
     
     FillListWithDecayGammaOrPi0(particle, d1, d2, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
     
   }// ndaugh > 0
 }

//___________________________________________________________________________
void AliGammaMCReader::CasePi0Decay(TParticle * pPi0, 
				    TClonesArray * plEMCAL, Int_t &indexEMCAL,
				    TClonesArray * plPHOS, Int_t &indexPHOS){
  
  //Decays pi0, see if aperture angle is small and then add the pi0 or the 2 gamma
  
  TLorentzVector lvPi0, lvGamma1, lvGamma2 ;
  //Double_t angle = 0;
  
  lvPi0.SetPxPyPzE(pPi0->Px(),pPi0->Py(),pPi0->Pz(),pPi0->Energy());

  //Decay
  MakePi0Decay(lvPi0,lvGamma1,lvGamma2);//,angle);
  
  //Check if Pi0 is in the acceptance of the calorimeters, if aperture angle is small, keep it
  TParticle * pPhoton1 = new TParticle(22,1,0,0,0,0,lvGamma1.Px(),lvGamma1.Py(),
				      lvGamma1.Pz(),lvGamma1.E(),0,0,0,0);   
  TParticle * pPhoton2 = new TParticle(22,1,0,0,0,0,lvGamma2.Px(),lvGamma2.Py(),
				      lvGamma2.Pz(),lvGamma2.E(),0,0,0,0);

  FillListWithDecayGammaOrPi0(pPi0, pPhoton1, pPhoton2, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
 
}

//_______________________________________________________________
void AliGammaMCReader::InitParameters()
{
  
  //Initialize the parameters of the analysis.

  fDecayPi0 = kGeantDecay;
  fCheckOverlapping = kTRUE ;

  Int_t pdgarray[]={12,14,16};// skip neutrinos
  fNeutralParticlesArray = new TArrayI(3, pdgarray);

}

//____________________________________________________________________________
void AliGammaMCReader::CreateParticleList(TObject * data, TObject *, 
					  TClonesArray * plCh, 
					  TClonesArray * plEMCAL,  
					  TClonesArray * plPHOS,
					  TClonesArray *plParton,TClonesArray *,TClonesArray *)
{
  //Create list of particles from EMCAL, PHOS and CTS. 
  AliStack * stack = (AliStack *) data ;
  Int_t indexCh     = plCh->GetEntries() ;
  Int_t indexEMCAL = plEMCAL->GetEntries() ;
  Int_t indexPHOS = plPHOS->GetEntries() ;
  Int_t indexParton = plParton->GetEntries() ;
  Int_t iParticle = 0 ;
  Double_t charge = 0.;
    
  for (iParticle=0 ; iParticle <  stack->GetNprimary() ; iParticle++) {
    TParticle * particle = stack->Particle(iParticle); 

    //Keep partons
    if(particle->GetStatusCode() == 21 && iParticle>=2){//All partons, not nucleus
      new((*plParton)[indexParton++])  TParticle(*particle) ;
    }

    //Keep Stable particles 
    if((particle->GetStatusCode() == 1) && (particle->Pt() > 0)){
      
      charge = TDatabasePDG::Instance()->GetParticle(particle->GetPdgCode())->Charge();
      
      //---------- Charged particles ----------------------
      if( fSwitchOnCTS && (charge != 0) && (particle->Pt() > fChargedPtCut)){
	//Particles in CTS acceptance
	if(TMath::Abs(particle->Eta())<fCTSEtaCut){  
	  //Fill lists
	  new((*plCh)[indexCh++])       TParticle(*particle) ;
	}
      }

      //-------------Neutral particles ----------------------
      else if((charge == 0) && particle->Pt() > fNeutralPtCut ){ //&&  
	//TMath::Abs(particle->GetPdgCode())>16){//Avoid neutrinos
	
	if(particle->GetPdgCode()!=111){
	  //In case that we force PYTHIA to decay pi0, and we want to check the overlapping of 
	  // the decay gamma.
	  if(particle->GetPdgCode() == 22 && fDecayPi0==kDecayGamma){
	    CaseDecayGamma(iParticle,particle, stack,plEMCAL, indexEMCAL, plPHOS, indexPHOS); //If pythia decays pi0
	  }
	  else{
	    //Take out some particles like neutrinos
	    if(!SkipNeutralParticles(particle->GetPdgCode()) && !particle->IsPrimary()){ // protection added (MG)
	      TParticle * pmother =stack->Particle(particle->GetFirstMother());	   	       
	      if(IsInPHOS(particle->Phi(),particle->Eta())){
		if(particle->GetPdgCode()==22)SetPhotonStatus(particle,pmother);     
		new((*plPHOS)[indexPHOS++])  TParticle(*particle) ;
	      }
	      else if(IsInEMCAL(particle->Phi(),particle->Eta())){
		if(particle->GetPdgCode()==22) SetPhotonStatus(particle,pmother); 
		new((*plEMCAL)[indexEMCAL++])  TParticle(*particle) ;
	      }
	    }//skip neutral particles
	  }//Case kDecayGamma
	}//no pi0
	else{
	  if(fDecayPi0 == kNoDecay){//keep the pi0 do not decay
	    if(IsInPHOS(particle->Phi(),particle->Eta()))
	      new((*plPHOS)[indexPHOS++])  TParticle(*particle) ;
	    else if(IsInEMCAL(particle->Phi(),particle->Eta()))
	      new((*plEMCAL)[indexEMCAL++])  TParticle(*particle) ;
	  }//nodecay
	  else if(fDecayPi0 == kDecay)
	    CasePi0Decay(particle,plEMCAL,indexEMCAL,plPHOS, indexPHOS);
	  else if(fDecayPi0 == kGeantDecay)
	    CaseGeantDecay(particle, stack,plEMCAL, indexEMCAL, plPHOS, indexPHOS);
	}//pi0  
      }//neutral particle
    }//stable particle
  }//particle loop
}


//___________________________________________________________________________
void AliGammaMCReader::FillListWithDecayGammaOrPi0(TParticle * pPi0, TParticle * pdaug0, TParticle * pdaug1,
				   TClonesArray * plEMCAL, Int_t &indexEMCAL,
				   TClonesArray * plPHOS, Int_t &indexPHOS){

  //Check if decay gamma overlapp in calorimeter, in such case keep the pi0, if not keep both photons.
  Bool_t  overlap = kFALSE ;
  TLorentzVector lv1 , lv2 ;
  pdaug0->Momentum(lv1);
  pdaug1->Momentum(lv2);
  Double_t angle = lv1.Angle(lv2.Vect());
 
  if(IsInEMCAL(pPi0->Phi(), pPi0->Eta()))
    {
      if (angle < fEMCALMinAngle){
	new((*plEMCAL)[indexEMCAL++])       TParticle(*pPi0) ;
	overlap = kTRUE;
      }
    }

  else if(IsInPHOS(pPi0->Phi(), pPi0->Eta())){
    if (angle < fPHOSMinAngle){
      new((*plPHOS)[indexPHOS++])       TParticle(*pPi0) ;
      overlap = kTRUE;
    }
  }//fCheckOverlapping

  //Fill with gammas if not overlapp
  if(!overlap){
    if(pdaug0->GetPdgCode() == 22 || TMath::Abs(pdaug0->GetPdgCode() ) == 11 ){
      pdaug0->SetStatusCode(kPi0DecayPhoton);
      if(IsInEMCAL(pdaug0->Phi(),pdaug0->Eta()) &&  pdaug0->Pt() > fNeutralPtCut)
	new((*plEMCAL)[indexEMCAL++])       TParticle(*pdaug0) ;
      else if(IsInPHOS(pdaug0->Phi(),pdaug0->Eta()) &&  pdaug0->Pt() > fNeutralPtCut)
	new((*plPHOS)[indexPHOS++])       TParticle(*pdaug0) ;
    }
    
    if(pdaug1->GetPdgCode() == 22 || TMath::Abs(pdaug1->GetPdgCode() ) == 11 ){
      pdaug1->SetStatusCode(kPi0DecayPhoton);
      if(IsInEMCAL(pdaug1->Phi(),pdaug1->Eta()) &&  pdaug1->Pt() > fNeutralPtCut)
	new((*plEMCAL)[indexEMCAL++])       TParticle(*pdaug1) ;
      else if(IsInPHOS(pdaug1->Phi(),pdaug1->Eta()) &&  pdaug1->Pt() > fNeutralPtCut)
	new((*plPHOS)[indexPHOS++])       TParticle(*pdaug1) ;
    }
  }// overlap?

 }


//___________________________________________________________________________
Bool_t  AliGammaMCReader::IsInEMCAL(Double_t phi, Double_t eta) const {
  //Check if particle is in EMCAL acceptance
 
  if(fSwitchOnEMCAL){
    if(phi<0)
      phi+=TMath::TwoPi();
    if( phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] && 
	TMath::Abs(eta)<fEMCALEtaCut) return kTRUE ;
    else  return kFALSE;     
  }
  return kFALSE ;
}

//___________________________________________________________________________
Bool_t  AliGammaMCReader::IsInPHOS(Double_t phi, Double_t eta) const {
  //Check if particle is in PHOS acceptance
  
  if(fSwitchOnPHOS){
    if(phi<0)
      phi+=TMath::TwoPi();
    if( phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] && 
	TMath::Abs(eta)<fPHOSEtaCut) return kTRUE ;
    else  return kFALSE;
  }
  
  return kFALSE ;
}

//____________________________________________________________________________
void AliGammaMCReader::MakePi0Decay(TLorentzVector &p0, TLorentzVector &p1, 
				TLorentzVector &p2){//, Double_t &angle) {
  // Perform isotropic decay pi0 -> 2 photons
  // p0 is pi0 4-momentum (inut)
  // p1 and p2 are photon 4-momenta (output)
  //  cout<<"Boost vector"<<endl;
  Double_t mPi0 = TDatabasePDG::Instance()->GetParticle(111)->Mass();
  TVector3 b = p0.BoostVector();
  //cout<<"Parameters"<<endl;
  //Double_t mPi0   = p0.M();
  Double_t phi    = TMath::TwoPi() * gRandom->Rndm();
  Double_t cosThe = 2 * gRandom->Rndm() - 1;
  Double_t cosPhi = TMath::Cos(phi);
  Double_t sinPhi = TMath::Sin(phi);
  Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
  Double_t ePi0   = mPi0/2.;
  //cout<<"ePi0 "<<ePi0<<endl;
  //cout<<"Components"<<endl;
  p1.SetPx(+ePi0*cosPhi*sinThe);
  p1.SetPy(+ePi0*sinPhi*sinThe);
  p1.SetPz(+ePi0*cosThe);
  p1.SetE(ePi0);
  //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
  //cout<<"p1 Mass: "<<p1.Px()*p1.Px()+p1.Py()*p1.Py()+p1.Pz()*p1.Pz()-p1.E()*p1.E()<<endl;
  p2.SetPx(-ePi0*cosPhi*sinThe);
  p2.SetPy(-ePi0*sinPhi*sinThe);
  p2.SetPz(-ePi0*cosThe);
  p2.SetE(ePi0);
  //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
  //cout<<"p2 Mass: "<<p2.Px()*p2.Px()+p2.Py()*p2.Py()+p2.Pz()*p2.Pz()-p2.E()*p2.E()<<endl;
  //cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
  p1.Boost(b);
  //cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
  p2.Boost(b);
  //cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
  //cout<<"angle"<<endl;
  //angle = p1.Angle(p2.Vect());
  //cout<<angle<<endl;
}

//________________________________________________________________
void AliGammaMCReader::Print(const Option_t * opt) const
{
  
  //Print some relevant parameters set for the analysis
  if(! opt)
    return;
  
  Info("Print", "%s %s", GetName(), GetTitle() ) ;
  
  printf("Decay Pi0?          : %d\n", fDecayPi0) ;
  printf("Check Overlapping?          : %d\n", fCheckOverlapping) ;

}


//________________________________________________________________
void AliGammaMCReader::SetPhotonStatus(TParticle* pphoton, TParticle* pmother)
{

  //Check the origin of the photon and tag it  as decay from pi0, from eta, from other mesons, or prompt or fragmentation.
  
  if(pmother->GetStatusCode() != 21 ){
    if(pmother->GetStatusCode() ==11){
      if(pmother->GetPdgCode()==111) pphoton->SetStatusCode(kPi0DecayPhoton);//decay gamma from pi0
      if(pmother->GetPdgCode()==221) pphoton->SetStatusCode(kEtaDecayPhoton);//decay gamma from eta
      else pphoton->SetStatusCode(kOtherDecayPhoton);// decay gamma from other pphotons	
    }
    else	 pphoton->SetStatusCode(kUnknown);
  }
  else if(pmother->GetPdgCode()==22)	pphoton->SetStatusCode(kPromptPhoton);//Prompt photon
  else pphoton->SetStatusCode(kFragmentPhoton); //Fragmentation photon
  
}

//________________________________________________________________
Bool_t AliGammaMCReader::SkipNeutralParticles(Int_t pdg) const {
  //Check if pdg is equal to one of the neutral particles list
  //These particles will be skipped from analysis.

  for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++)
    if(TMath::Abs(pdg) ==  fNeutralParticlesArray->At(i)) return kTRUE ;
  
  return kFALSE; 
  
}
Commit	Line	Data
bdcfac30	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 *
4b707925	12	* about the suitability of this software for any purpose. It is *
bdcfac30	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15	/* $Id$ */
	16
	17	/* History of cvs commits:
	18	*
	19	* $Log$
2021085e	20	* Revision 1.4 2007/10/29 13:48:42 gustavo
	21	* Corrected coding violations
	22	*
bdcfac30	23	*
	24	*/
	25
	26	//_________________________________________________________________________
	27	// Class for reading data (Kinematics) in order to do prompt gamma correlations
	28	// Class created from old AliPHOSGammaJet
	29	// (see AliRoot versions previous Release 4-09)
	30	//
	31	//*-- Author: Gustavo Conesa (LNF-INFN)
	32	//////////////////////////////////////////////////////////////////////////////
	33
	34
	35	// --- ROOT system ---
	36
	37	#include <TFile.h>
	38	#include <TParticle.h>
	39	#include <TH2.h>
	40	#include <TChain.h>
	41	#include <TRandom.h>
3bb2c538	42	#include <TArrayI.h>
3bb2c538	43	#include <TClonesArray.h>
bdcfac30	44
	45	//---- ANALYSIS system ----
	46	#include "AliGammaMCReader.h"
	47	#include "Riostream.h"
	48	#include "AliLog.h"
3bb2c538	49	#include "AliStack.h"
bdcfac30	50
	51	ClassImp(AliGammaMCReader)
	52
	53	//____________________________________________________________________________
	54	AliGammaMCReader::AliGammaMCReader() :
3bb2c538	55	AliGammaReader(), fDecayPi0(0), fCheckOverlapping(kFALSE), fNeutralParticlesArray(0x0)
bdcfac30	56	{
	57	//Ctor
	58
	59	//Initialize parameters
bdcfac30	60	InitParameters();
4b707925	61	fDataType = kMC;
bdcfac30	62
	63	}
	64
	65	//____________________________________________________________________________
	66	AliGammaMCReader::AliGammaMCReader(const AliGammaMCReader & g) :
3bb2c538	67	AliGammaReader(g), fDecayPi0(g.fDecayPi0), fCheckOverlapping(g.fCheckOverlapping),
3bb2c538	68	fNeutralParticlesArray(g.fNeutralParticlesArray?new TArrayI(*g.fNeutralParticlesArray):0x0)
bdcfac30	69	{
	70	// cpy ctor
	71	}
	72
	73	//_________________________________________________________________________
	74	AliGammaMCReader & AliGammaMCReader::operator = (const AliGammaMCReader & source)
	75	{
	76	// assignment operator
	77
	78	if(&source == this) return *this;
	79
4b707925	80	fDecayPi0 = source.fDecayPi0;
4b707925	81	fCheckOverlapping = source.fCheckOverlapping;
3bb2c538	82	delete fNeutralParticlesArray;
3bb2c538	83	fNeutralParticlesArray = source.fNeutralParticlesArray?new TArrayI(*source.fNeutralParticlesArray):0x0;
bdcfac30	84	return *this;
	85
	86	}
	87
3bb2c538	88	//_________________________________
	89	AliGammaMCReader::~AliGammaMCReader() {
	90	//Dtor
	91
	92	delete fNeutralParticlesArray ;
	93
	94	}
	95
4b707925	96	//___________________________________________________________________________
	97	void AliGammaMCReader::CaseDecayGamma(Int_t index, TParticle * particle, AliStack * sta,
	98	TClonesArray * plEMCAL, Int_t &indexEMCAL,
	99	TClonesArray * plPHOS, Int_t &indexPHOS){
	100	//In case pi0 are decayed by pythia, check if mother is pi0 and in such case look if
	101	//there is overlapping. Send particle=pi0 if there is overlapping.
3bb2c538	102
4b707925	103	TParticle * pmother =sta->Particle(particle->GetFirstMother());
	104	if(pmother->GetPdgCode() == 111 && pmother->GetNDaughters() == 2) {//Do not consider 3 particle decay case
	105	Int_t idaug0 = pmother->GetDaughter(0);
	106	Int_t idaug1 = pmother->GetDaughter(1);
	107	TParticle * pdaug0 = sta -> Particle(idaug0);
	108	TParticle * pdaug1 = sta -> Particle(idaug1);
3bb2c538	109
4b707925	110	if((index == idaug0 && pdaug0->Pt() > fNeutralPtCut) \|\|
3bb2c538	111	(index == idaug1 && pdaug0->Pt() <= fNeutralPtCut)){//Check decay when first daughter arrives, do nothing with second.
3bb2c538	112
4b707925	113	FillListWithDecayGammaOrPi0(pmother, pdaug0, pdaug1, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
4b707925	114
3bb2c538	115	}
4b707925	116	}//mother is a pi0 with 2 daughters
4b707925	117	else{
3bb2c538	118
	119	if(IsInPHOS(particle->Phi(),particle->Eta())) {
	120	TParticle * pmother =sta->Particle(particle->GetFirstMother());
	121	SetPhotonStatus(particle,pmother);
4b707925	122	new((plPHOS)[indexPHOS++]) TParticle(particle) ;
3bb2c538	123	}
	124	else if(IsInEMCAL(particle->Phi(),particle->Eta())){
	125	TParticle * pmother =sta->Particle(particle->GetFirstMother());
	126	SetPhotonStatus(particle,pmother);
4b707925	127	new((plEMCAL)[indexEMCAL++]) TParticle(particle) ;
3bb2c538	128	}
3bb2c538	129	}
4b707925	130	}
	131
	132	//___________________________________________________________________________
	133	void AliGammaMCReader::CaseGeantDecay(TParticle * particle, AliStack * stack,
	134	TClonesArray * plEMCAL, Int_t &indexEMCAL,
	135	TClonesArray * plPHOS, Int_t &indexPHOS){
	136	//Find decay gamma from pi0, decayed by GEANT and put them in the list.
	137
	138	Int_t ndaug = particle->GetNDaughters() ;
	139	TParticle * d1 = new TParticle();
	140	TParticle * d2 = new TParticle();
	141
	142	if(ndaug > 0){//At least 1 daugther
	143	d1 = stack->Particle(particle->GetDaughter(0));
	144	if (ndaug > 1 )
	145	d2 = stack->Particle(particle->GetDaughter(1));
	146
	147	FillListWithDecayGammaOrPi0(particle, d1, d2, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
	148
	149	}// ndaugh > 0
	150	}
	151
	152	//___________________________________________________________________________
	153	void AliGammaMCReader::CasePi0Decay(TParticle * pPi0,
	154	TClonesArray * plEMCAL, Int_t &indexEMCAL,
	155	TClonesArray * plPHOS, Int_t &indexPHOS){
	156
	157	//Decays pi0, see if aperture angle is small and then add the pi0 or the 2 gamma
	158
	159	TLorentzVector lvPi0, lvGamma1, lvGamma2 ;
	160	//Double_t angle = 0;
	161
	162	lvPi0.SetPxPyPzE(pPi0->Px(),pPi0->Py(),pPi0->Pz(),pPi0->Energy());
	163
	164	//Decay
	165	MakePi0Decay(lvPi0,lvGamma1,lvGamma2);//,angle);
	166
	167	//Check if Pi0 is in the acceptance of the calorimeters, if aperture angle is small, keep it
	168	TParticle * pPhoton1 = new TParticle(22,1,0,0,0,0,lvGamma1.Px(),lvGamma1.Py(),
	169	lvGamma1.Pz(),lvGamma1.E(),0,0,0,0);
	170	TParticle * pPhoton2 = new TParticle(22,1,0,0,0,0,lvGamma2.Px(),lvGamma2.Py(),
	171	lvGamma2.Pz(),lvGamma2.E(),0,0,0,0);
	172
	173	FillListWithDecayGammaOrPi0(pPi0, pPhoton1, pPhoton2, plEMCAL, indexEMCAL, plPHOS, indexPHOS);
	174
	175	}
	176
	177	//_______________________________________________________________
	178	void AliGammaMCReader::InitParameters()
	179	{
	180
	181	//Initialize the parameters of the analysis.
	182
	183	fDecayPi0 = kGeantDecay;
	184	fCheckOverlapping = kTRUE ;
3bb2c538	185
	186	Int_t pdgarray[]={12,14,16};// skip neutrinos
	187	fNeutralParticlesArray = new TArrayI(3, pdgarray);
	188
4b707925	189	}
bdcfac30	190
	191	//____________________________________________________________________________
	192	void AliGammaMCReader::CreateParticleList(TObject * data, TObject *,
4b707925	193	TClonesArray * plCh,
4b707925	194	TClonesArray * plEMCAL,
bdcfac30	195	TClonesArray * plPHOS,
4b707925	196	TClonesArray plParton,TClonesArray ,TClonesArray *)
bdcfac30	197	{
	198	//Create list of particles from EMCAL, PHOS and CTS.
	199	AliStack * stack = (AliStack *) data ;
	200	Int_t indexCh = plCh->GetEntries() ;
	201	Int_t indexEMCAL = plEMCAL->GetEntries() ;
	202	Int_t indexPHOS = plPHOS->GetEntries() ;
	203	Int_t indexParton = plParton->GetEntries() ;
	204	Int_t iParticle = 0 ;
	205	Double_t charge = 0.;
	206
	207	for (iParticle=0 ; iParticle < stack->GetNprimary() ; iParticle++) {
	208	TParticle * particle = stack->Particle(iParticle);
4b707925	209
bdcfac30	210	//Keep partons
	211	if(particle->GetStatusCode() == 21 && iParticle>=2){//All partons, not nucleus
	212	new((plParton)[indexParton++]) TParticle(particle) ;
	213	}
	214
	215	//Keep Stable particles
4b707925	216	if((particle->GetStatusCode() == 1) && (particle->Pt() > 0)){
bdcfac30	217
	218	charge = TDatabasePDG::Instance()->GetParticle(particle->GetPdgCode())->Charge();
	219
	220	//---------- Charged particles ----------------------
3bb2c538	221	if( fSwitchOnCTS && (charge != 0) && (particle->Pt() > fChargedPtCut)){
bdcfac30	222	//Particles in CTS acceptance
	223	if(TMath::Abs(particle->Eta())<fCTSEtaCut){
	224	//Fill lists
	225	new((plCh)[indexCh++]) TParticle(particle) ;
	226	}
	227	}
3bb2c538	228
bdcfac30	229	//-------------Neutral particles ----------------------
3bb2c538	230	else if((charge == 0) && particle->Pt() > fNeutralPtCut ){ //&&
3bb2c538	231	//TMath::Abs(particle->GetPdgCode())>16){//Avoid neutrinos
bdcfac30	232
bdcfac30	233	if(particle->GetPdgCode()!=111){
4b707925	234	//In case that we force PYTHIA to decay pi0, and we want to check the overlapping of
	235	// the decay gamma.
	236	if(particle->GetPdgCode() == 22 && fDecayPi0==kDecayGamma){
	237	CaseDecayGamma(iParticle,particle, stack,plEMCAL, indexEMCAL, plPHOS, indexPHOS); //If pythia decays pi0
	238	}
	239	else{
3bb2c538	240	//Take out some particles like neutrinos
2021085e	241	if(!SkipNeutralParticles(particle->GetPdgCode()) && !particle->IsPrimary()){ // protection added (MG)
3bb2c538	242	TParticle * pmother =stack->Particle(particle->GetFirstMother());
	243	if(IsInPHOS(particle->Phi(),particle->Eta())){
	244	if(particle->GetPdgCode()==22)SetPhotonStatus(particle,pmother);
	245	new((plPHOS)[indexPHOS++]) TParticle(particle) ;
	246	}
	247	else if(IsInEMCAL(particle->Phi(),particle->Eta())){
	248	if(particle->GetPdgCode()==22) SetPhotonStatus(particle,pmother);
	249	new((plEMCAL)[indexEMCAL++]) TParticle(particle) ;
	250	}
	251	}//skip neutral particles
	252	}//Case kDecayGamma
bdcfac30	253	}//no pi0
	254	else{
	255	if(fDecayPi0 == kNoDecay){//keep the pi0 do not decay
	256	if(IsInPHOS(particle->Phi(),particle->Eta()))
	257	new((plPHOS)[indexPHOS++]) TParticle(particle) ;
	258	else if(IsInEMCAL(particle->Phi(),particle->Eta()))
	259	new((plEMCAL)[indexEMCAL++]) TParticle(particle) ;
3bb2c538	260	}//nodecay
bdcfac30	261	else if(fDecayPi0 == kDecay)
4b707925	262	CasePi0Decay(particle,plEMCAL,indexEMCAL,plPHOS, indexPHOS);
bdcfac30	263	else if(fDecayPi0 == kGeantDecay)
4b707925	264	CaseGeantDecay(particle, stack,plEMCAL, indexEMCAL, plPHOS, indexPHOS);
bdcfac30	265	}//pi0
	266	}//neutral particle
	267	}//stable particle
	268	}//particle loop
	269	}
	270
4b707925	271
	272	//___________________________________________________________________________
	273	void AliGammaMCReader::FillListWithDecayGammaOrPi0(TParticle * pPi0, TParticle * pdaug0, TParticle * pdaug1,
	274	TClonesArray * plEMCAL, Int_t &indexEMCAL,
	275	TClonesArray * plPHOS, Int_t &indexPHOS){
	276
	277	//Check if decay gamma overlapp in calorimeter, in such case keep the pi0, if not keep both photons.
4b707925	278	Bool_t overlap = kFALSE ;
	279	TLorentzVector lv1 , lv2 ;
	280	pdaug0->Momentum(lv1);
	281	pdaug1->Momentum(lv2);
	282	Double_t angle = lv1.Angle(lv2.Vect());
3bb2c538	283
	284	if(IsInEMCAL(pPi0->Phi(), pPi0->Eta()))
	285	{
4b707925	286	if (angle < fEMCALMinAngle){
	287	new((plEMCAL)[indexEMCAL++]) TParticle(pPi0) ;
	288	overlap = kTRUE;
	289	}
	290	}
3bb2c538	291
	292	else if(IsInPHOS(pPi0->Phi(), pPi0->Eta())){
	293	if (angle < fPHOSMinAngle){
	294	new((plPHOS)[indexPHOS++]) TParticle(pPi0) ;
	295	overlap = kTRUE;
4b707925	296	}
4b707925	297	}//fCheckOverlapping
3bb2c538	298
4b707925	299	//Fill with gammas if not overlapp
	300	if(!overlap){
	301	if(pdaug0->GetPdgCode() == 22 \|\| TMath::Abs(pdaug0->GetPdgCode() ) == 11 ){
3bb2c538	302	pdaug0->SetStatusCode(kPi0DecayPhoton);
4b707925	303	if(IsInEMCAL(pdaug0->Phi(),pdaug0->Eta()) && pdaug0->Pt() > fNeutralPtCut)
	304	new((plEMCAL)[indexEMCAL++]) TParticle(pdaug0) ;
	305	else if(IsInPHOS(pdaug0->Phi(),pdaug0->Eta()) && pdaug0->Pt() > fNeutralPtCut)
	306	new((plPHOS)[indexPHOS++]) TParticle(pdaug0) ;
	307	}
	308
	309	if(pdaug1->GetPdgCode() == 22 \|\| TMath::Abs(pdaug1->GetPdgCode() ) == 11 ){
3bb2c538	310	pdaug1->SetStatusCode(kPi0DecayPhoton);
4b707925	311	if(IsInEMCAL(pdaug1->Phi(),pdaug1->Eta()) && pdaug1->Pt() > fNeutralPtCut)
	312	new((plEMCAL)[indexEMCAL++]) TParticle(pdaug1) ;
	313	else if(IsInPHOS(pdaug1->Phi(),pdaug1->Eta()) && pdaug1->Pt() > fNeutralPtCut)
	314	new((plPHOS)[indexPHOS++]) TParticle(pdaug1) ;
	315	}
	316	}// overlap?
3bb2c538	317
4b707925	318	}
	319
	320
	321
bdcfac30	322	//___________________________________________________________________________
3bb2c538	323	Bool_t AliGammaMCReader::IsInEMCAL(Double_t phi, Double_t eta) const {
bdcfac30	324	//Check if particle is in EMCAL acceptance
3bb2c538	325
	326	if(fSwitchOnEMCAL){
	327	if(phi<0)
	328	phi+=TMath::TwoPi();
	329	if( phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] &&
	330	TMath::Abs(eta)<fEMCALEtaCut) return kTRUE ;
	331	else return kFALSE;
	332	}
bdcfac30	333	return kFALSE ;
	334	}
	335
	336	//___________________________________________________________________________
3bb2c538	337	Bool_t AliGammaMCReader::IsInPHOS(Double_t phi, Double_t eta) const {
	338	//Check if particle is in PHOS acceptance
	339
	340	if(fSwitchOnPHOS){
	341	if(phi<0)
	342	phi+=TMath::TwoPi();
	343	if( phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] &&
	344	TMath::Abs(eta)<fPHOSEtaCut) return kTRUE ;
	345	else return kFALSE;
	346	}
bdcfac30	347
	348	return kFALSE ;
	349	}
	350
bdcfac30	351	//____________________________________________________________________________
4b707925	352	void AliGammaMCReader::MakePi0Decay(TLorentzVector &p0, TLorentzVector &p1,
4b707925	353	TLorentzVector &p2){//, Double_t &angle) {
bdcfac30	354	// Perform isotropic decay pi0 -> 2 photons
	355	// p0 is pi0 4-momentum (inut)
	356	// p1 and p2 are photon 4-momenta (output)
	357	// cout<<"Boost vector"<<endl;
	358	Double_t mPi0 = TDatabasePDG::Instance()->GetParticle(111)->Mass();
	359	TVector3 b = p0.BoostVector();
	360	//cout<<"Parameters"<<endl;
	361	//Double_t mPi0 = p0.M();
	362	Double_t phi = TMath::TwoPi() * gRandom->Rndm();
	363	Double_t cosThe = 2 * gRandom->Rndm() - 1;
	364	Double_t cosPhi = TMath::Cos(phi);
	365	Double_t sinPhi = TMath::Sin(phi);
	366	Double_t sinThe = TMath::Sqrt(1-cosThe*cosThe);
	367	Double_t ePi0 = mPi0/2.;
	368	//cout<<"ePi0 "<<ePi0<<endl;
	369	//cout<<"Components"<<endl;
	370	p1.SetPx(+ePi0cosPhisinThe);
	371	p1.SetPy(+ePi0sinPhisinThe);
	372	p1.SetPz(+ePi0*cosThe);
	373	p1.SetE(ePi0);
	374	//cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
	375	//cout<<"p1 Mass: "<<p1.Px()p1.Px()+p1.Py()p1.Py()+p1.Pz()p1.Pz()-p1.E()p1.E()<<endl;
	376	p2.SetPx(-ePi0cosPhisinThe);
	377	p2.SetPy(-ePi0sinPhisinThe);
	378	p2.SetPz(-ePi0*cosThe);
	379	p2.SetE(ePi0);
	380	//cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
	381	//cout<<"p2 Mass: "<<p2.Px()p2.Px()+p2.Py()p2.Py()+p2.Pz()p2.Pz()-p2.E()p2.E()<<endl;
	382	//cout<<"Boost "<<b.X()<<" "<<b.Y()<<" "<<b.Z()<<endl;
	383	p1.Boost(b);
	384	//cout<<"p1: "<<p1.Px()<<" "<<p1.Py()<<" "<<p1.Pz()<<" "<<p1.E()<<endl;
	385	p2.Boost(b);
	386	//cout<<"p2: "<<p2.Px()<<" "<<p2.Py()<<" "<<p2.Pz()<<" "<<p2.E()<<endl;
	387	//cout<<"angle"<<endl;
4b707925	388	//angle = p1.Angle(p2.Vect());
bdcfac30	389	//cout<<angle<<endl;
	390	}
	391
	392	//________________________________________________________________
	393	void AliGammaMCReader::Print(const Option_t * opt) const
	394	{
	395
	396	//Print some relevant parameters set for the analysis
	397	if(! opt)
	398	return;
	399
	400	Info("Print", "%s %s", GetName(), GetTitle() ) ;
	401
bdcfac30	402	printf("Decay Pi0? : %d\n", fDecayPi0) ;
4b707925	403	printf("Check Overlapping? : %d\n", fCheckOverlapping) ;
4b707925	404
bdcfac30	405	}
4b707925	406
	407
	408
3bb2c538	409	//________________________________________________________________
	410	void AliGammaMCReader::SetPhotonStatus(TParticle* pphoton, TParticle* pmother)
	411	{
4b707925	412
3bb2c538	413	//Check the origin of the photon and tag it as decay from pi0, from eta, from other mesons, or prompt or fragmentation.
	414
	415	if(pmother->GetStatusCode() != 21 ){
	416	if(pmother->GetStatusCode() ==11){
	417	if(pmother->GetPdgCode()==111) pphoton->SetStatusCode(kPi0DecayPhoton);//decay gamma from pi0
	418	if(pmother->GetPdgCode()==221) pphoton->SetStatusCode(kEtaDecayPhoton);//decay gamma from eta
	419	else pphoton->SetStatusCode(kOtherDecayPhoton);// decay gamma from other pphotons
	420	}
	421	else pphoton->SetStatusCode(kUnknown);
	422	}
	423	else if(pmother->GetPdgCode()==22) pphoton->SetStatusCode(kPromptPhoton);//Prompt photon
	424	else pphoton->SetStatusCode(kFragmentPhoton); //Fragmentation photon
	425
	426	}
	427
	428	//________________________________________________________________
	429	Bool_t AliGammaMCReader::SkipNeutralParticles(Int_t pdg) const {
	430	//Check if pdg is equal to one of the neutral particles list
	431	//These particles will be skipped from analysis.
	432
	433	for(Int_t i= 0; i < fNeutralParticlesArray->GetSize(); i++)
	434	if(TMath::Abs(pdg) == fNeutralParticlesArray->At(i)) return kTRUE ;
	435
	436	return kFALSE;
	437
	438	}