[u/mrichter/AliRoot.git] / PWG4 / AliGammaMCDataReader.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.2  2007/08/17 12:40:04  schutz
 * New analysis classes by Gustavo Conesa
 *
 * Revision 1.1.2.1  2007/07/26 10:32:09  schutz
 * new analysis classes in the the new analysis framework
 *
 *
 */

//_________________________________________________________________________
// Class for reading data (Kinematics and ESDs) 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 "Riostream.h"
#include <TParticle.h>

//---- ANALYSIS system ----
#include "AliGammaMCDataReader.h" 
#include "AliLog.h"
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliESDCaloCluster.h"
#include "AliStack.h"

ClassImp(AliGammaMCDataReader)

//____________________________________________________________________________
  AliGammaMCDataReader::AliGammaMCDataReader() : 
    AliGammaReader()
{
  //Default Ctor
  
  //Initialize parameters
  fDataType=kMCData;
  
}

//____________________________________________________________________________
AliGammaMCDataReader::AliGammaMCDataReader(const AliGammaMCDataReader & g) :   
  AliGammaReader(g)
{
  // cpy ctor
}

//_________________________________________________________________________
AliGammaMCDataReader & AliGammaMCDataReader::operator = (const AliGammaMCDataReader & source)
{
  // assignment operator
  
  if(&source == this) return *this;
  

  return *this;
  
}

//____________________________________________________________________________
void AliGammaMCDataReader::CreateParticleList(TObject * data, TObject * kine,
					      TClonesArray * plCTS, 
					      TClonesArray * plEMCAL,  
					      TClonesArray * plPHOS,   
					      TClonesArray * plPrimCTS, 
					      TClonesArray * plPrimEMCAL, 
					      TClonesArray * plPrimPHOS){
  
  //Create a list of particles from the ESD. These particles have been measured 
  //by the Central Tracking system (TPC+ITS), PHOS and EMCAL 

  AliESDEvent* esd = (AliESDEvent*) data;
  AliStack* stack = (AliStack*) kine;
  
  Int_t npar  = 0 ;
  Double_t *pid = new Double_t[AliPID::kSPECIESN];  
  AliDebug(3,"Fill particle lists");
  
  //Get vertex for momentum calculation  
  Double_t v[3] ; //vertex ;
  esd->GetVertex()->GetXYZ(v) ; 
  
  //########### PHOS ##############
  
  Int_t begphos = 0;  
  Int_t endphos = esd->GetNumberOfCaloClusters() ;  
  Int_t indexPH = plPHOS->GetEntries() ;
  Int_t begem = 0;
  Int_t endem = endphos;
  
  AliDebug(3,Form("First Calorimeter particle %d, last particle %d", begphos,endphos));
  
  for (npar =  begphos; npar <  endphos; npar++) {//////////////PHOS track loop
    AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd

    //Check if it is PHOS cluster
    if(!clus->IsEMCAL())
      begem=npar+1;
    else
      continue;

    AliDebug(4,Form("PHOS clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
    if(clus->GetTrackMatched()==-1 ){
      //Create a TParticle to fill the particle list
      TLorentzVector momentum ;
      clus->GetMomentum(momentum, v);
      Double_t phi = momentum.Phi();
      if(phi<0) phi+=TMath::TwoPi() ;
      if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fPHOSEtaCut &&
	 phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] ) {
	
	pid=clus->GetPid();	
	Int_t pdg = 22;

    	if(IsPHOSPIDOn()){
	  AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ",
			  momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
			  pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
	  
	  Float_t wPhoton =  fPHOSPhotonWeight;
	  Float_t wPi0 =  fPHOSPi0Weight;
	  
	  if(fPHOSWeightFormula){
	    wPhoton = fPHOSPhotonWeightFormula->Eval(momentum.E()) ;
	    wPi0 =    fPHOSPi0WeightFormula->Eval(momentum.E());
	  }
	  
	  if(pid[AliPID::kPhoton] > wPhoton) 
	    pdg = kPhoton ;
	  else if(pid[AliPID::kPi0] > wPi0) 
	    pdg = kPi0 ; 
	  else if(pid[AliPID::kElectron] > fPHOSElectronWeight)  
	    pdg = kElectron ;
	  else if(pid[AliPID::kEleCon] > fPHOSElectronWeight) 
	    pdg = kEleCon ;
	  else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fPHOSChargeWeight) 
	    pdg = kChargedHadron ;  
	  else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fPHOSNeutralWeight) 
	    pdg = kNeutralHadron ; 
	  
	  else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton]  >  
		  pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron]) 
	    pdg = kChargedUnknown  ; 
	  else 
	    pdg = kNeutralUnknown ; 
	  //neutral cluster, unidentifed.
	}
	
	if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown ){//keep only neutral particles in the array

	  //###############
	  //Check kinematics
	  //###############
	  Int_t label = clus->GetLabel();
	  TParticle * pmother = GetMotherParticle(label,stack, "PHOS",momentum);

	  TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
					       momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
	  
	  AliDebug(4,Form("PHOS added: pdg %d, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
	  //printf("PHOS added: pdg %d, pt %f, phi %f, eta %f\n", pdg, particle->Pt(),particle->Phi(),particle->Eta());
	  new((*plPHOS)[indexPH])   TParticle(*particle) ;
	  new((*plPrimPHOS)[indexPH])   TParticle(*pmother) ;
	  indexPH++;
	  
	}
	else AliDebug(4,Form("PHOS charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f\n", pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));	

      }//pt, eta, phi cut
    }//not charged
  }//cluster loop

  //########### CTS (TPC+ITS) #####################
  Int_t begtpc   = 0 ;  
  Int_t endtpc   = esd->GetNumberOfTracks() ;
  Int_t indexCh  = plCTS->GetEntries() ;
  AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
  
  for (npar =  begtpc; npar <  endtpc; npar++) {////////////// track loop
    AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
    
    //We want tracks fitted in the detectors:
    ULong_t status=AliESDtrack::kTPCrefit;
    status|=AliESDtrack::kITSrefit;
    
    //We want tracks whose PID bit is set:
    //     ULong_t status =AliESDtrack::kITSpid;
    //     status|=AliESDtrack::kTPCpid;
    
    if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
      // Do something with the tracks which were successfully
      // re-fitted 
      Double_t en = 0; //track ->GetTPCsignal() ;
      Double_t mom[3];
      track->GetPxPyPz(mom) ;
      Double_t px = mom[0];
      Double_t py = mom[1];
      Double_t pz = mom[2]; //Check with TPC people if this is correct.
      Int_t pdg = 11; //Give any charged PDG code, in this case electron.
      //I just want to tag the particle as charged

      //Check kinematics
      Int_t label = track->GetLabel();
      TParticle * pmother = new TParticle();
      if(label>=0){
	pmother = stack->Particle(label);
	// 	    AliDebug(4,("Primary Track: Pt %2.2f, mother: label %d, pdg  %d, status %d,  mother 2 %d, grandmother %d",
	// 			TMath::Sqrt(px*px+py*py), label, pmother->GetPdgCode(),pmother->GetStatusCode(), pmother->GetMother(0), stack->GetPrimary(label)));
      }
      else AliDebug(3, Form("Negative Kinematic label of track:  %d",label));
      
      TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
					   px, py, pz, en, v[0], v[1], v[2], 0);
      
      if(particle->Pt() > fChargedPtCut && TMath::Abs(particle->Eta())<fCTSEtaCut){
	new((*plCTS)[indexCh])       TParticle(*particle) ;   
	new((*plPrimCTS)[indexCh])       TParticle(*pmother) ;
	indexCh++;    
	
      }
    }
  }
  
  //################ EMCAL ##############
  
  Int_t indexEM  = plEMCAL->GetEntries() ; 
  
  
  for (npar =  begem; npar <  endem; npar++) {//////////////EMCAL track loop
    AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
    Int_t clustertype= clus->GetClusterType();
    if(clustertype == AliESDCaloCluster::kEMCALClusterv1 && clus->GetTrackMatched()==-1 ){
      
      TLorentzVector momentum ;
      clus->GetMomentum(momentum, v);          
      Double_t phi = momentum.Phi();
      if(phi<0) phi+=TMath::TwoPi() ;
      if(momentum.Pt() > fNeutralPtCut &&  TMath::Abs(momentum.Eta()) < fEMCALEtaCut &&
	 phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] ) {
	
      pid=clus->GetPid();	
      Int_t pdg = 22;

      if(IsEMCALPIDOn()){
	AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ",
			momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
			pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
	
	if(pid[AliPID::kPhoton] > fEMCALPhotonWeight) 
	    pdg = kPhoton ;
	else if(pid[AliPID::kPi0] > fEMCALPi0Weight) 
	  pdg = kPi0 ; 
	else if(pid[AliPID::kElectron] > fEMCALElectronWeight)  
	  pdg = kElectron ;
	else if(pid[AliPID::kEleCon] > fEMCALElectronWeight) 
	  pdg = kEleCon ;
	else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fEMCALChargeWeight) 
	  pdg = kChargedHadron ;  
	else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fEMCALNeutralWeight) 
	  pdg = kNeutralHadron ; 
	else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton]  >  
		pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron]) 
	  pdg = kChargedUnknown ; 
	else 
	  pdg = kNeutralUnknown ;
      }
      
      if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown){//keep only neutral particles in the array
	
	//###############
	//Check kinematics
	//###############
	Int_t label = clus->GetLabel();
	TParticle * pmother = GetMotherParticle(label,stack, "EMCAL",momentum);
	
	TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1, 
					     momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
	AliDebug(4,Form("EMCAL cluster added: pdg %f, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
	
	new((*plEMCAL)[indexEM])   TParticle(*particle) ;
	new((*plPrimEMCAL)[indexEM])   TParticle(*pmother) ;
	indexEM++;
      }
      else AliDebug(4,Form("EMCAL charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f", 
			   pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));
      
      }//pt, phi, eta cut
      else 	AliDebug(4,"Particle not added");
      
    }//not charged, not pseudocluster
  }//Cluster loop
  
  AliDebug(3,Form("Particle lists filled, PHOS clusters %d, EMCAL clusters %d, CTS tracks %d", 
		  indexPH,indexEM,indexCh));
  
}


TParticle * AliGammaMCDataReader:: GetMotherParticle(Int_t label, AliStack *stack, TString calo,  TLorentzVector momentum)
{
  //Gets the primary particle and do some checks:
  //Check if primary is inside calorimeter and look the mother outsie
  //Check if mother is a decay photon, in which case check if decay was overlapped
  
  Float_t minangle = 0;
  Float_t ipdist = 0;
  TParticle * pmother = new TParticle();

  if(calo == "PHOS"){
    ipdist= fPHOSIPDistance;
    minangle = fPHOSMinAngle; 
  }
  else if (calo == "EMCAL"){
    ipdist = fEMCALIPDistance;
    minangle = fEMCALMinAngle;
  }


  if(label>=0){
    pmother = stack->Particle(label);
    Int_t mostatus = pmother->GetStatusCode();
    Int_t mopdg    = TMath::Abs(pmother->GetPdgCode());
    
    //Check if mother is a conversion inside the calorimeter
    //In such case, return the mother before the calorimeter.
    //First approximation.
    Float_t vy = TMath::Abs(pmother->Vy()) ;

    if( mostatus == 0 && vy >= ipdist){

      //cout<<"Calo: "<<calo<<" vy "<<vy<<" E clus "<<momentum.E()<<" Emother "<<pmother->Energy()<<" "
      //  <<pmother->GetName()<<endl;

      while( vy >= ipdist){//inside calorimeter
	AliDebug(4,Form("Conversion inside calorimeter, mother vertex %0.2f, IP distance %0.2f", vy, ipdist));
	pmother =  stack->Particle(pmother->GetMother(0));
	vy = TMath::Abs(pmother->Vy()) ;
	//cout<<" label "<<label<<" Mother: "<<pmother->GetName()<<" E "<<pmother->Energy()<<" Status "<<pmother->GetStatusCode()<<"  and vertex "<<vy<<endl;
	mostatus = pmother->GetStatusCode();
	mopdg    = TMath::Abs(pmother->GetPdgCode());
      }//while vertex is inside calorimeter
      //cout<<"Calo: "<<calo<<" final vy "<<vy<<" E clus "<<momentum.E()<<" Emother "<<pmother->Energy()<<" "
      //  <<pmother->GetName()<<endl;
    }//check status and vertex

    AliDebug(4,Form("%s, ESD E %2.2f, PID %d,  mother: E %2.2f, label %d, status %d,  vertex %3.2f, mother 2 %d, grandmother %d \n",
		    calo.Data(),momentum.E(),pmother->Energy(), label, pmother->GetPdgCode(),
		    pmother->GetStatusCode(), vy, pmother->GetMother(0), stack->GetPrimary(label)));
    
    //Check Decay photons
    if(mopdg == 22){
      
      //his mother was a pi0?
      TParticle * pmotherpi0 =  stack->Particle(pmother->GetMother(0));
      if( pmotherpi0->GetPdgCode() == 111){

	AliDebug(4,Form(" %s: E cluster %2.2f, E gamma %2.2f, Mother Pi0, E %0.2f, status %d, daughters %d\n",
			calo.Data(), momentum.E(),pmother->Energy(),pmotherpi0->Energy(), 
			pmotherpi0->GetStatusCode(), pmotherpi0->GetNDaughters()));
	
	if(pmotherpi0->GetNDaughters() == 1) mostatus = 2 ; //signal that this photon can only be decay, not overlapped.
	else if(pmotherpi0->GetNDaughters() == 2){
	  
	  TParticle * pd1 = stack->Particle(pmotherpi0->GetDaughter(0));
	  TParticle * pd2 = stack->Particle(pmotherpi0->GetDaughter(1));
	  //if(pmotherpi0->Energy()> 10 ){
// 	  cout<<"Two "<<calo<<" daugthers, pi0 label "<<pmother->GetMother(0)<<" E :"<<pmotherpi0->Energy()<<endl;
// 	  cout<<" 1) label "<<pmotherpi0->GetDaughter(0)<<" pdg "<<pd1->GetPdgCode()<<" E  "<<pd1->Energy()
// 	      <<" phi "<<pd1->Phi()*TMath::RadToDeg()<<" eta "<<pd1->Eta()
// 	      <<" mother label "<<pd1->GetMother(0)<<" n daug "<<pd1->GetNDaughters() <<endl;
// 	    cout<<" 2) label "<<pmotherpi0->GetDaughter(1)<<" pdg "<<pd2->GetPdgCode()<<" E  "<<pd2->Energy()
// 		<<" phi "<<pd2->Phi()*TMath::RadToDeg()<<" eta "<<pd2->Eta()<<" mother label "
// 		<<pd2->GetMother(0)<<" n daug "<<pd2->GetNDaughters() <<endl;
	    //}
	  if(pd1->GetPdgCode() == 22 && pd2->GetPdgCode() == 22){
	    TLorentzVector lv1 , lv2 ;
	    pd1->Momentum(lv1);
	    pd2->Momentum(lv2);
	    Double_t angle = lv1.Angle(lv2.Vect());
// 	    if(pmotherpi0->Energy()> 10 )
// 	      cout<<"angle*ipdist "<<angle*ipdist<<" mindist "<< minangle <<endl;
	    if (angle < minangle){
	      //There is overlapping, pass the pi0 as mother
	      
	      AliDebug(4,Form(">>>> %s cluster with E %2.2f is a overlapped pi0, E %2.2f, angle %2.4f < anglemin %2.4f\n",
			      calo.Data(), momentum.E(), pmotherpi0->Energy(), angle, minangle));	    
	     
	      pmother = pmotherpi0 ;
	      
	    }
	    else mostatus = 2 ; // gamma decay not overlapped
	  }// daughters are photons
	  else mostatus = 2; // daughters are e-gamma or e-e, no overlapped, or charged cluster
	}//2 daughters
	else AliDebug(4,Form("pi0 has %d daughters",pmotherpi0->GetNDaughters()));
      }//pi0 decay photon?
    }//photon 

    pmother->SetStatusCode(mostatus); // if status = 2, decay photon.

  }//label >=0
  else AliInfo(Form("Negative Kinematic label of PHOS cluster:  %d",label));

  return pmother ;

}
Commit	Line	Data
bdcfac30	1
	2	/**************************************************************************
	3	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	4	* *
	5	* Author: The ALICE Off-line Project. *
	6	* Contributors are mentioned in the code where appropriate. *
	7	* *
	8	* Permission to use, copy, modify and distribute this software and its *
	9	* documentation strictly for non-commercial purposes is hereby granted *
	10	* without fee, provided that the above copyright notice appears in all *
	11	* copies and that both the copyright notice and this permission notice *
	12	* appear in the supporting documentation. The authors make no claims *
	13	* about the suitability of this software for any purpose. It is *
	14	* provided "as is" without express or implied warranty. *
	15	**************************************************************************/
	16	/* $Id$ */
	17
	18	/* History of cvs commits:
	19	*
	20	* $Log$
4b707925	21	* Revision 1.2 2007/08/17 12:40:04 schutz
	22	* New analysis classes by Gustavo Conesa
	23	*
bdcfac30	24	* Revision 1.1.2.1 2007/07/26 10:32:09 schutz
	25	* new analysis classes in the the new analysis framework
	26	*
	27	*
	28	*/
	29
	30	//_________________________________________________________________________
	31	// Class for reading data (Kinematics and ESDs) in order to do prompt gamma correlations
	32	// Class created from old AliPHOSGammaJet
	33	// (see AliRoot versions previous Release 4-09)
	34	//
	35	//*-- Author: Gustavo Conesa (LNF-INFN)
	36	//////////////////////////////////////////////////////////////////////////////
	37
	38
	39	// --- ROOT system ---
	40	#include "Riostream.h"
	41	#include <TParticle.h>
	42
	43	//---- ANALYSIS system ----
	44	#include "AliGammaMCDataReader.h"
	45	#include "AliLog.h"
	46	#include "AliESDEvent.h"
	47	#include "AliESDVertex.h"
	48	#include "AliESDCaloCluster.h"
	49	#include "AliStack.h"
	50
	51	ClassImp(AliGammaMCDataReader)
	52
	53	//____________________________________________________________________________
	54	AliGammaMCDataReader::AliGammaMCDataReader() :
4b707925	55	AliGammaReader()
bdcfac30	56	{
	57	//Default Ctor
	58
	59	//Initialize parameters
	60	fDataType=kMCData;
bdcfac30	61
	62	}
	63
	64	//____________________________________________________________________________
	65	AliGammaMCDataReader::AliGammaMCDataReader(const AliGammaMCDataReader & g) :
4b707925	66	AliGammaReader(g)
bdcfac30	67	{
	68	// cpy ctor
	69	}
	70
	71	//_________________________________________________________________________
	72	AliGammaMCDataReader & AliGammaMCDataReader::operator = (const AliGammaMCDataReader & source)
	73	{
	74	// assignment operator
	75
	76	if(&source == this) return *this;
	77
4b707925	78
bdcfac30	79	return *this;
	80
	81	}
	82
	83	//____________________________________________________________________________
	84	void AliGammaMCDataReader::CreateParticleList(TObject * data, TObject * kine,
	85	TClonesArray * plCTS,
	86	TClonesArray * plEMCAL,
4b707925	87	TClonesArray * plPHOS,
	88	TClonesArray * plPrimCTS,
	89	TClonesArray * plPrimEMCAL,
	90	TClonesArray * plPrimPHOS){
bdcfac30	91
	92	//Create a list of particles from the ESD. These particles have been measured
	93	//by the Central Tracking system (TPC+ITS), PHOS and EMCAL
	94
	95	AliESDEvent* esd = (AliESDEvent*) data;
	96	AliStack* stack = (AliStack*) kine;
	97
	98	Int_t npar = 0 ;
4b707925	99	Double_t *pid = new Double_t[AliPID::kSPECIESN];
bdcfac30	100	AliDebug(3,"Fill particle lists");
	101
	102	//Get vertex for momentum calculation
	103	Double_t v[3] ; //vertex ;
	104	esd->GetVertex()->GetXYZ(v) ;
	105
	106	//########### PHOS ##############
	107
4b707925	108	Int_t begphos = 0;
4b707925	109	Int_t endphos = esd->GetNumberOfCaloClusters() ;
bdcfac30	110	Int_t indexPH = plPHOS->GetEntries() ;
4b707925	111	Int_t begem = 0;
	112	Int_t endem = endphos;
	113
	114	AliDebug(3,Form("First Calorimeter particle %d, last particle %d", begphos,endphos));
bdcfac30	115
	116	for (npar = begphos; npar < endphos; npar++) {//////////////PHOS track loop
	117	AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
4b707925	118
	119	//Check if it is PHOS cluster
	120	if(!clus->IsEMCAL())
	121	begem=npar+1;
	122	else
	123	continue;
	124
	125	AliDebug(4,Form("PHOS clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
bdcfac30	126	if(clus->GetTrackMatched()==-1 ){
	127	//Create a TParticle to fill the particle list
	128	TLorentzVector momentum ;
	129	clus->GetMomentum(momentum, v);
4b707925	130	Double_t phi = momentum.Phi();
	131	if(phi<0) phi+=TMath::TwoPi() ;
	132	if(momentum.Pt() > fNeutralPtCut && TMath::Abs(momentum.Eta()) < fPHOSEtaCut &&
	133	phi > fPhiPHOSCut[0] && phi < fPhiPHOSCut[1] ) {
	134
	135	pid=clus->GetPid();
	136	Int_t pdg = 22;
bdcfac30	137
4b707925	138	if(IsPHOSPIDOn()){
	139	AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ",
	140	momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
	141	pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
	142
	143	Float_t wPhoton = fPHOSPhotonWeight;
	144	Float_t wPi0 = fPHOSPi0Weight;
	145
	146	if(fPHOSWeightFormula){
	147	wPhoton = fPHOSPhotonWeightFormula->Eval(momentum.E()) ;
	148	wPi0 = fPHOSPi0WeightFormula->Eval(momentum.E());
	149	}
	150
	151	if(pid[AliPID::kPhoton] > wPhoton)
	152	pdg = kPhoton ;
	153	else if(pid[AliPID::kPi0] > wPi0)
	154	pdg = kPi0 ;
	155	else if(pid[AliPID::kElectron] > fPHOSElectronWeight)
	156	pdg = kElectron ;
	157	else if(pid[AliPID::kEleCon] > fPHOSElectronWeight)
	158	pdg = kEleCon ;
	159	else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fPHOSChargeWeight)
	160	pdg = kChargedHadron ;
	161	else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fPHOSNeutralWeight)
	162	pdg = kNeutralHadron ;
	163
	164	else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] >
	165	pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron])
	166	pdg = kChargedUnknown ;
	167	else
	168	pdg = kNeutralUnknown ;
	169	//neutral cluster, unidentifed.
	170	}
	171
	172	if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown ){//keep only neutral particles in the array
bdcfac30	173
4b707925	174	//###############
	175	//Check kinematics
	176	//###############
	177	Int_t label = clus->GetLabel();
	178	TParticle * pmother = GetMotherParticle(label,stack, "PHOS",momentum);
	179
	180	TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
	181	momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
	182
	183	AliDebug(4,Form("PHOS added: pdg %d, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
	184	//printf("PHOS added: pdg %d, pt %f, phi %f, eta %f\n", pdg, particle->Pt(),particle->Phi(),particle->Eta());
	185	new((plPHOS)[indexPH]) TParticle(particle) ;
	186	new((plPrimPHOS)[indexPH]) TParticle(pmother) ;
	187	indexPH++;
	188
	189	}
	190	else AliDebug(4,Form("PHOS charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f\n", pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));
	191
	192	}//pt, eta, phi cut
bdcfac30	193	}//not charged
	194	}//cluster loop
	195
	196	//########### CTS (TPC+ITS) #####################
	197	Int_t begtpc = 0 ;
	198	Int_t endtpc = esd->GetNumberOfTracks() ;
	199	Int_t indexCh = plCTS->GetEntries() ;
	200	AliDebug(3,Form("First CTS particle %d, last particle %d", begtpc,endtpc));
	201
	202	for (npar = begtpc; npar < endtpc; npar++) {////////////// track loop
	203	AliESDtrack * track = esd->GetTrack(npar) ; // retrieve track from esd
	204
	205	//We want tracks fitted in the detectors:
	206	ULong_t status=AliESDtrack::kTPCrefit;
	207	status\|=AliESDtrack::kITSrefit;
	208
	209	//We want tracks whose PID bit is set:
	210	// ULong_t status =AliESDtrack::kITSpid;
	211	// status\|=AliESDtrack::kTPCpid;
	212
	213	if ( (track->GetStatus() & status) == status) {//Check if the bits we want are set
	214	// Do something with the tracks which were successfully
	215	// re-fitted
	216	Double_t en = 0; //track ->GetTPCsignal() ;
	217	Double_t mom[3];
	218	track->GetPxPyPz(mom) ;
	219	Double_t px = mom[0];
	220	Double_t py = mom[1];
	221	Double_t pz = mom[2]; //Check with TPC people if this is correct.
	222	Int_t pdg = 11; //Give any charged PDG code, in this case electron.
	223	//I just want to tag the particle as charged
4b707925	224
	225	//Check kinematics
	226	Int_t label = track->GetLabel();
	227	TParticle * pmother = new TParticle();
	228	if(label>=0){
	229	pmother = stack->Particle(label);
	230	// AliDebug(4,("Primary Track: Pt %2.2f, mother: label %d, pdg %d, status %d, mother 2 %d, grandmother %d",
	231	// TMath::Sqrt(pxpx+pypy), label, pmother->GetPdgCode(),pmother->GetStatusCode(), pmother->GetMother(0), stack->GetPrimary(label)));
	232	}
	233	else AliDebug(3, Form("Negative Kinematic label of track: %d",label));
	234
	235	TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
	236	px, py, pz, en, v[0], v[1], v[2], 0);
	237
	238	if(particle->Pt() > fChargedPtCut && TMath::Abs(particle->Eta())<fCTSEtaCut){
	239	new((plCTS)[indexCh]) TParticle(particle) ;
	240	new((plPrimCTS)[indexCh]) TParticle(pmother) ;
	241	indexCh++;
	242
	243	}
bdcfac30	244	}
	245	}
	246
	247	//################ EMCAL ##############
	248
	249	Int_t indexEM = plEMCAL->GetEntries() ;
bdcfac30	250
bdcfac30	251
	252	for (npar = begem; npar < endem; npar++) {//////////////EMCAL track loop
	253	AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve track from esd
	254	Int_t clustertype= clus->GetClusterType();
4b707925	255	if(clustertype == AliESDCaloCluster::kEMCALClusterv1 && clus->GetTrackMatched()==-1 ){
bdcfac30	256
bdcfac30	257	TLorentzVector momentum ;
4b707925	258	clus->GetMomentum(momentum, v);
	259	Double_t phi = momentum.Phi();
	260	if(phi<0) phi+=TMath::TwoPi() ;
	261	if(momentum.Pt() > fNeutralPtCut && TMath::Abs(momentum.Eta()) < fEMCALEtaCut &&
	262	phi > fPhiEMCALCut[0] && phi < fPhiEMCALCut[1] ) {
	263
bdcfac30	264	pid=clus->GetPid();
	265	Int_t pdg = 22;
	266
4b707925	267	if(IsEMCALPIDOn()){
	268	AliDebug(5,Form("E %1.2f; PID: ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f,pi %0.2f, k %0.2f, p %0.2f, k0 %0.2f, n %0.2f, mu %0.2f ",
	269	momentum.E(),pid[AliPID::kPhoton],pid[AliPID::kPi0],pid[AliPID::kElectron],pid[AliPID::kEleCon],pid[AliPID::kPion],
	270	pid[AliPID::kKaon],pid[AliPID::kProton], pid[AliPID::kKaon0],pid[AliPID::kNeutron], pid[AliPID::kMuon]));
	271
	272	if(pid[AliPID::kPhoton] > fEMCALPhotonWeight)
	273	pdg = kPhoton ;
	274	else if(pid[AliPID::kPi0] > fEMCALPi0Weight)
	275	pdg = kPi0 ;
	276	else if(pid[AliPID::kElectron] > fEMCALElectronWeight)
	277	pdg = kElectron ;
	278	else if(pid[AliPID::kEleCon] > fEMCALElectronWeight)
	279	pdg = kEleCon ;
	280	else if(pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] > fEMCALChargeWeight)
	281	pdg = kChargedHadron ;
	282	else if(pid[AliPID::kKaon0]+pid[AliPID::kNeutron] > fEMCALNeutralWeight)
	283	pdg = kNeutralHadron ;
	284	else if(pid[AliPID::kElectron]+pid[AliPID::kEleCon]+pid[AliPID::kPion]+pid[AliPID::kKaon]+pid[AliPID::kProton] >
	285	pid[AliPID::kPhoton] + pid[AliPID::kPi0]+pid[AliPID::kKaon0]+pid[AliPID::kNeutron])
	286	pdg = kChargedUnknown ;
	287	else
	288	pdg = kNeutralUnknown ;
bdcfac30	289	}
bdcfac30	290
4b707925	291	if(pdg != kElectron && pdg != kEleCon && pdg !=kChargedHadron && pdg !=kChargedUnknown){//keep only neutral particles in the array
	292
	293	//###############
	294	//Check kinematics
	295	//###############
	296	Int_t label = clus->GetLabel();
	297	TParticle * pmother = GetMotherParticle(label,stack, "EMCAL",momentum);
	298
	299	TParticle * particle = new TParticle(pdg, 1, -1, -1, -1, -1,
	300	momentum.Px(), momentum.Py(), momentum.Pz(), momentum.E(), v[0], v[1], v[2], 0);
	301	AliDebug(4,Form("EMCAL cluster added: pdg %f, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
	302
	303	new((plEMCAL)[indexEM]) TParticle(particle) ;
	304	new((plPrimEMCAL)[indexEM]) TParticle(pmother) ;
	305	indexEM++;
bdcfac30	306	}
4b707925	307	else AliDebug(4,Form("EMCAL charged cluster NOT added: pdg %d, pt %f, phi %f, eta %f",
4b707925	308	pdg, momentum.Pt(),momentum.Phi(),momentum.Eta()));
bdcfac30	309
4b707925	310	}//pt, phi, eta cut
4b707925	311	else AliDebug(4,"Particle not added");
bdcfac30	312
	313	}//not charged, not pseudocluster
	314	}//Cluster loop
	315
4b707925	316	AliDebug(3,Form("Particle lists filled, PHOS clusters %d, EMCAL clusters %d, CTS tracks %d",
4b707925	317	indexPH,indexEM,indexCh));
bdcfac30	318
	319	}
	320
4b707925	321
4b707925	322	TParticle * AliGammaMCDataReader:: GetMotherParticle(Int_t label, AliStack *stack, TString calo, TLorentzVector momentum)
bdcfac30	323	{
4b707925	324	//Gets the primary particle and do some checks:
	325	//Check if primary is inside calorimeter and look the mother outsie
	326	//Check if mother is a decay photon, in which case check if decay was overlapped
	327
	328	Float_t minangle = 0;
	329	Float_t ipdist = 0;
	330	TParticle * pmother = new TParticle();
bdcfac30	331
4b707925	332	if(calo == "PHOS"){
	333	ipdist= fPHOSIPDistance;
	334	minangle = fPHOSMinAngle;
	335	}
	336	else if (calo == "EMCAL"){
	337	ipdist = fEMCALIPDistance;
	338	minangle = fEMCALMinAngle;
	339	}
bdcfac30	340
bdcfac30	341
4b707925	342	if(label>=0){
	343	pmother = stack->Particle(label);
	344	Int_t mostatus = pmother->GetStatusCode();
	345	Int_t mopdg = TMath::Abs(pmother->GetPdgCode());
	346
	347	//Check if mother is a conversion inside the calorimeter
	348	//In such case, return the mother before the calorimeter.
	349	//First approximation.
	350	Float_t vy = TMath::Abs(pmother->Vy()) ;
bdcfac30	351
4b707925	352	if( mostatus == 0 && vy >= ipdist){
bdcfac30	353
4b707925	354	//cout<<"Calo: "<<calo<<" vy "<<vy<<" E clus "<<momentum.E()<<" Emother "<<pmother->Energy()<<" "
	355	// <<pmother->GetName()<<endl;
	356
	357	while( vy >= ipdist){//inside calorimeter
	358	AliDebug(4,Form("Conversion inside calorimeter, mother vertex %0.2f, IP distance %0.2f", vy, ipdist));
	359	pmother = stack->Particle(pmother->GetMother(0));
	360	vy = TMath::Abs(pmother->Vy()) ;
	361	//cout<<" label "<<label<<" Mother: "<<pmother->GetName()<<" E "<<pmother->Energy()<<" Status "<<pmother->GetStatusCode()<<" and vertex "<<vy<<endl;
	362	mostatus = pmother->GetStatusCode();
	363	mopdg = TMath::Abs(pmother->GetPdgCode());
	364	}//while vertex is inside calorimeter
	365	//cout<<"Calo: "<<calo<<" final vy "<<vy<<" E clus "<<momentum.E()<<" Emother "<<pmother->Energy()<<" "
	366	// <<pmother->GetName()<<endl;
	367	}//check status and vertex
	368
	369	AliDebug(4,Form("%s, ESD E %2.2f, PID %d, mother: E %2.2f, label %d, status %d, vertex %3.2f, mother 2 %d, grandmother %d \n",
	370	calo.Data(),momentum.E(),pmother->Energy(), label, pmother->GetPdgCode(),
	371	pmother->GetStatusCode(), vy, pmother->GetMother(0), stack->GetPrimary(label)));
	372
	373	//Check Decay photons
	374	if(mopdg == 22){
	375
	376	//his mother was a pi0?
	377	TParticle * pmotherpi0 = stack->Particle(pmother->GetMother(0));
	378	if( pmotherpi0->GetPdgCode() == 111){
	379
	380	AliDebug(4,Form(" %s: E cluster %2.2f, E gamma %2.2f, Mother Pi0, E %0.2f, status %d, daughters %d\n",
	381	calo.Data(), momentum.E(),pmother->Energy(),pmotherpi0->Energy(),
	382	pmotherpi0->GetStatusCode(), pmotherpi0->GetNDaughters()));
	383
	384	if(pmotherpi0->GetNDaughters() == 1) mostatus = 2 ; //signal that this photon can only be decay, not overlapped.
	385	else if(pmotherpi0->GetNDaughters() == 2){
	386
	387	TParticle * pd1 = stack->Particle(pmotherpi0->GetDaughter(0));
	388	TParticle * pd2 = stack->Particle(pmotherpi0->GetDaughter(1));
	389	//if(pmotherpi0->Energy()> 10 ){
	390	// cout<<"Two "<<calo<<" daugthers, pi0 label "<<pmother->GetMother(0)<<" E :"<<pmotherpi0->Energy()<<endl;
	391	// cout<<" 1) label "<<pmotherpi0->GetDaughter(0)<<" pdg "<<pd1->GetPdgCode()<<" E "<<pd1->Energy()
	392	// <<" phi "<<pd1->Phi()*TMath::RadToDeg()<<" eta "<<pd1->Eta()
	393	// <<" mother label "<<pd1->GetMother(0)<<" n daug "<<pd1->GetNDaughters() <<endl;
	394	// cout<<" 2) label "<<pmotherpi0->GetDaughter(1)<<" pdg "<<pd2->GetPdgCode()<<" E "<<pd2->Energy()
	395	// <<" phi "<<pd2->Phi()*TMath::RadToDeg()<<" eta "<<pd2->Eta()<<" mother label "
	396	// <<pd2->GetMother(0)<<" n daug "<<pd2->GetNDaughters() <<endl;
	397	//}
	398	if(pd1->GetPdgCode() == 22 && pd2->GetPdgCode() == 22){
	399	TLorentzVector lv1 , lv2 ;
	400	pd1->Momentum(lv1);
	401	pd2->Momentum(lv2);
	402	Double_t angle = lv1.Angle(lv2.Vect());
	403	// if(pmotherpi0->Energy()> 10 )
	404	// cout<<"angleipdist "<<angleipdist<<" mindist "<< minangle <<endl;
	405	if (angle < minangle){
	406	//There is overlapping, pass the pi0 as mother
	407
	408	AliDebug(4,Form(">>>> %s cluster with E %2.2f is a overlapped pi0, E %2.2f, angle %2.4f < anglemin %2.4f\n",
	409	calo.Data(), momentum.E(), pmotherpi0->Energy(), angle, minangle));
	410
	411	pmother = pmotherpi0 ;
	412
	413	}
	414	else mostatus = 2 ; // gamma decay not overlapped
	415	}// daughters are photons
	416	else mostatus = 2; // daughters are e-gamma or e-e, no overlapped, or charged cluster
	417	}//2 daughters
418	else AliDebug(4,Form("pi0 has %d daughters",pmotherpi0->GetNDaughters()));
419	}//pi0 decay photon?
420	}//photon
421
422	pmother->SetStatusCode(mostatus); // if status = 2, decay photon.
423
424	}//label >=0
425	else AliInfo(Form("Negative Kinematic label of PHOS cluster: %d",label));
426
427	return pmother ;
bdcfac30	428
bdcfac30	429	}