[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.5  2007/12/14 10:54:15  gustavo
 * Protection against cluster labels larger than kinematics entries
 *
 * Revision 1.4  2007/10/29 13:48:42  gustavo
 * Corrected coding violations
 *
 * 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 <TFormula.h>
#include <TParticle.h>
 
//---- ANALYSIS system ----
#include "AliGammaMCDataReader.h" 
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliESDCaloCluster.h"
#include "AliStack.h"
#include "AliLog.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 
  //Also create particle list with mothers.

  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) ; 
  
  //########### CALORIMETERS ##############  
  Int_t nCaloCluster = esd->GetNumberOfCaloClusters() ;  
  Int_t indexPH = plPHOS->GetEntries() ;
  Int_t indexEM = plEMCAL->GetEntries() ;
  
  for (npar =  0; npar <  nCaloCluster; npar++) {//////////////CaloCluster loop
    AliESDCaloCluster * clus = esd->GetCaloCluster(npar) ; // retrieve cluster from esd
    Int_t type = clus->GetClusterType();
    
    //########### PHOS ##############
    if(fSwitchOnPHOS && type ==  AliESDCaloCluster::kPHOSCluster){
      AliDebug(4,Form("PHOS clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
      
      if(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()) < 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

	    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);	    
	    new((*plPHOS)[indexPH])   TParticle(*particle) ;
	    AliDebug(4,Form("PHOS added: pdg %d, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
 
	    //###############
	    //Check kinematics
	    //###############
	    TParticle * pmother = new TParticle();
	    Int_t label = clus->GetLabel();
	    if(label < stack->GetNtrack())
	      pmother = GetMotherParticle(label,stack, "PHOS",momentum);
	    else 
	      AliInfo(Form("PHOS: Bad label %d, too large, NPrimaries %d",label,stack->GetNtrack()));
	    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
	else 	AliDebug(4,"Particle not added");
      }//track-match?
    }//PHOS cluster

    //################ EMCAL ##############
    else if(fSwitchOnEMCAL &&  type ==  AliESDCaloCluster::kEMCALClusterv1){
      AliDebug(4,Form("EMCAL clusters: E %f, match %d", clus->E(),clus->GetTrackMatched()));
      
      if(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

	    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);
	    new((*plEMCAL)[indexEM])   TParticle(*particle) ;
	    AliDebug(4,Form("EMCAL cluster added: pdg %f, pt %f, phi %f, eta %f", pdg, particle->Pt(),particle->Phi(),particle->Eta()));
	    
	    //###############
	    //Check kinematics
	    //###############
	    TParticle * pmother = new TParticle();
	    Int_t label = clus->GetLabel();
	    if(label < stack->GetNtrack())
	      pmother = GetMotherParticle(label,stack, "EMCAL",momentum);
	    else 
	      AliInfo(Form("EMCAL: Bad label %d, too large, NPrimaries %d",label,stack->GetNtrack()));	    
	    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");
      }//track-matched
    }//EMCAL cluster

  }//cluster loop

 
  //########### CTS (TPC+ITS) #####################
  Int_t nTracks   = esd->GetNumberOfTracks() ;
  Int_t indexCh  = plCTS->GetEntries() ;

  if(fSwitchOnCTS){
    AliDebug(3,Form("Number of tracks %d",nTracks));
  
    for (npar =  0; npar <  nTracks; 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 = TMath::Abs(track->GetLabel());
	TParticle * pmother = new TParticle();
	pmother = stack->Particle(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++;    
	  
	}//kinematic selection
      }//select track from refit
    }//track loop    
  }//CTS

  AliDebug(3,Form("Particle lists filled, tracks  %d , clusters: EMCAL %d, PHOS %d", indexCh,indexEM,indexPH));

}

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