[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaPhoton.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 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: AliAnaPhoton.cxx 28688 2008-09-11 15:04:07Z gconesab $ */

//_________________________________________________________________________
//
// Class for the photon identification.
// Clusters from calorimeters are identified as photons
// and kept in the AOD. Few histograms produced.
//
// -- Author: Gustavo Conesa (LNF-INFN) 
//////////////////////////////////////////////////////////////////////////////
  
  
// --- ROOT system --- 
#include <TH2F.h>

// --- Analysis system --- 
#include "AliAnaPhoton.h" 
#include "AliCaloTrackReader.h"
#include "AliCaloPID.h"

ClassImp(AliAnaPhoton)
  
//____________________________________________________________________________
  AliAnaPhoton::AliAnaPhoton() : 
    AliAnaPartCorrBaseClass(), fCalorimeter(""), 
	fMinDist(0.),fMinDist2(0.),fMinDist3(0.),
	fhPtPhoton(0),fhPhiPhoton(0),fhEtaPhoton(0),
	//MC
    fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0), 
    fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0), 
    fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0), 
    fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0), 
    fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0), 
    fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0)
{
  //default ctor
  
  //Initialize parameters
  InitParameters();

}

//____________________________________________________________________________
AliAnaPhoton::AliAnaPhoton(const AliAnaPhoton & g) : 
  AliAnaPartCorrBaseClass(g), fCalorimeter(g.fCalorimeter),
   fMinDist(g.fMinDist),fMinDist2(g.fMinDist2), fMinDist3(g.fMinDist3),
   fhPtPhoton(g.fhPtPhoton),fhPhiPhoton(g.fhPhiPhoton),fhEtaPhoton(g.fhEtaPhoton), 
  //MC
  fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt), 
  fhPtFragmentation(g.fhPtFragmentation),fhPhiFragmentation(g.fhPhiFragmentation),fhEtaFragmentation(g.fhEtaFragmentation), 
  fhPtPi0Decay(g.fhPtPi0Decay),fhPhiPi0Decay(g.fhPhiPi0Decay),fhEtaPi0Decay(g.fhEtaPi0Decay), 
  fhPtOtherDecay(g.fhPtOtherDecay),fhPhiOtherDecay(g.fhPhiOtherDecay),fhEtaOtherDecay(g.fhEtaOtherDecay), 
  fhPtConversion(g. fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion), 
  fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown) 
{
  // cpy ctor
  
}

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

  fCalorimeter = g.fCalorimeter ;
  fMinDist  = g.fMinDist;
  fMinDist2 = g.fMinDist2;
  fMinDist3 = g.fMinDist3;
  
  fhPtPhoton = g.fhPtPhoton ; 
  fhPhiPhoton = g.fhPhiPhoton ;
  fhEtaPhoton = g.fhEtaPhoton ;
 
  fhPtPrompt = g.fhPtPrompt;
  fhPhiPrompt = g.fhPhiPrompt;
  fhEtaPrompt = g.fhEtaPrompt; 
  fhPtFragmentation = g.fhPtFragmentation;
  fhPhiFragmentation = g.fhPhiFragmentation;
  fhEtaFragmentation = g.fhEtaFragmentation; 
  fhPtPi0Decay = g.fhPtPi0Decay;
  fhPhiPi0Decay = g.fhPhiPi0Decay;
  fhEtaPi0Decay = g.fhEtaPi0Decay; 
  fhPtOtherDecay = g.fhPtOtherDecay;
  fhPhiOtherDecay = g.fhPhiOtherDecay;
  fhEtaOtherDecay = g.fhEtaOtherDecay; 
  fhPtConversion = g. fhPtConversion;
  fhPhiConversion = g.fhPhiConversion;
  fhEtaConversion = g.fhEtaConversion; 
  fhPtUnknown = g.fhPtUnknown;
  fhPhiUnknown = g.fhPhiUnknown;
  fhEtaUnknown = g.fhEtaUnknown; 

  return *this;
  
}

//____________________________________________________________________________
AliAnaPhoton::~AliAnaPhoton() 
{
  //dtor

}


//________________________________________________________________________
TList *  AliAnaPhoton::GetCreateOutputObjects()
{  
	// Create histograms to be saved in output file and 
	// store them in outputContainer
	TList * outputContainer = new TList() ; 
	outputContainer->SetName("PhotonHistos") ; 
	
	Int_t nptbins  = GetHistoNPtBins();
	Int_t nphibins = GetHistoNPhiBins();
	Int_t netabins = GetHistoNEtaBins();
	Float_t ptmax  = GetHistoPtMax();
	Float_t phimax = GetHistoPhiMax();
	Float_t etamax = GetHistoEtaMax();
	Float_t ptmin  = GetHistoPtMin();
	Float_t phimin = GetHistoPhiMin();
	Float_t etamin = GetHistoEtaMin();	

	//Histograms of highest Photon identified in Event
	fhPtPhoton  = new TH1F("hPtPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
	fhPtPhoton->SetYTitle("N");
	fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
	outputContainer->Add(fhPtPhoton) ; 
	
	fhPhiPhoton  = new TH2F
    ("hPhiPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
	fhPhiPhoton->SetYTitle("#phi");
	fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
	outputContainer->Add(fhPhiPhoton) ; 
	
	fhEtaPhoton  = new TH2F
    ("hEtaPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
	fhEtaPhoton->SetYTitle("#eta");
	fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
	outputContainer->Add(fhEtaPhoton) ;
	
	if(IsDataMC()){
		
		fhPtPrompt  = new TH1F("hPtPrompt","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
		fhPtPrompt->SetYTitle("N");
		fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
		outputContainer->Add(fhPtPrompt) ; 
		
		fhPhiPrompt  = new TH2F
		("hPhiPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
		fhPhiPrompt->SetYTitle("#phi");
		fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
		outputContainer->Add(fhPhiPrompt) ; 
		
		fhEtaPrompt  = new TH2F
		("hEtaPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
		fhEtaPrompt->SetYTitle("#eta");
		fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
		outputContainer->Add(fhEtaPrompt) ;
		
		fhPtFragmentation  = new TH1F("hPtFragmentation","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
		fhPtFragmentation->SetYTitle("N");
		fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
		outputContainer->Add(fhPtFragmentation) ; 
		
		fhPhiFragmentation  = new TH2F
		("hPhiFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
		fhPhiFragmentation->SetYTitle("#phi");
		fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
		outputContainer->Add(fhPhiFragmentation) ; 
		
		fhEtaFragmentation  = new TH2F
		("hEtaFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
		fhEtaFragmentation->SetYTitle("#eta");
		fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
		outputContainer->Add(fhEtaFragmentation) ;
		
		fhPtPi0Decay  = new TH1F("hPtPi0Decay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
		fhPtPi0Decay->SetYTitle("N");
		fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
		outputContainer->Add(fhPtPi0Decay) ; 
		
		fhPhiPi0Decay  = new TH2F
		("hPhiPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
		fhPhiPi0Decay->SetYTitle("#phi");
		fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
		outputContainer->Add(fhPhiPi0Decay) ; 
		
		fhEtaPi0Decay  = new TH2F
		("hEtaPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
		fhEtaPi0Decay->SetYTitle("#eta");
		fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
		outputContainer->Add(fhEtaPi0Decay) ;
		
		fhPtOtherDecay  = new TH1F("hPtOtherDecay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
		fhPtOtherDecay->SetYTitle("N");
		fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
		outputContainer->Add(fhPtOtherDecay) ; 
		
		fhPhiOtherDecay  = new TH2F
		("hPhiOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
		fhPhiOtherDecay->SetYTitle("#phi");
		fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
		outputContainer->Add(fhPhiOtherDecay) ; 
		
		fhEtaOtherDecay  = new TH2F
		("hEtaOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
		fhEtaOtherDecay->SetYTitle("#eta");
		fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
		outputContainer->Add(fhEtaOtherDecay) ;
		
		fhPtConversion  = new TH1F("hPtConversion","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
		fhPtConversion->SetYTitle("N");
		fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
		outputContainer->Add(fhPtConversion) ; 
		
		fhPhiConversion  = new TH2F
		("hPhiConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
		fhPhiConversion->SetYTitle("#phi");
		fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
		outputContainer->Add(fhPhiConversion) ; 
		
		fhEtaConversion  = new TH2F
		("hEtaConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
		fhEtaConversion->SetYTitle("#eta");
		fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
		outputContainer->Add(fhEtaConversion) ;
		
		fhPtUnknown  = new TH1F("hPtUnknown","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); 
		fhPtUnknown->SetYTitle("N");
		fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
		outputContainer->Add(fhPtUnknown) ; 
		
		fhPhiUnknown  = new TH2F
		("hPhiUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); 
		fhPhiUnknown->SetYTitle("#phi");
		fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
		outputContainer->Add(fhPhiUnknown) ; 
		
		fhEtaUnknown  = new TH2F
		("hEtaUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); 
		fhEtaUnknown->SetYTitle("#eta");
		fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
		outputContainer->Add(fhEtaUnknown) ;
	}//Histos with MC
	
 	//Save parameters used for analysis
	TString parList ; //this will be list of parameters used for this analysis.
	char onePar[255] ;
	
	sprintf(onePar,"--- AliAnaPhoton ---\n") ;
	parList+=onePar ;	
	sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;
	parList+=onePar ;
	sprintf(onePar,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
	parList+=onePar ;
	sprintf(onePar,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
	parList+=onePar ;
	sprintf(onePar,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
	parList+=onePar ;
	
 	//Get parameters set in base class.
	parList += GetBaseParametersList() ;
	
	//Get parameters set in PID class.
	parList += GetCaloPID()->GetPIDParametersList() ;
	
 	//Get parameters set in FidutialCut class (not available yet)
	//parlist += GetCaloPID()->GetFidCutParametersList() 
	
 	TObjString *oString= new TObjString(parList) ;
	outputContainer->Add(oString);
	
	return outputContainer ;
	
}

//____________________________________________________________________________
void AliAnaPhoton::InitParameters()
{
  
  //Initialize the parameters of the analysis.
  SetOutputAODClassName("AliAODPWG4Particle");
  SetOutputAODName("photons");
  fCalorimeter = "PHOS" ;
  fMinDist  = 2.;
  fMinDist2 = 4.;
  fMinDist3 = 5.;
  
}

//__________________________________________________________________
void  AliAnaPhoton::MakeAnalysisFillAOD() 
{
	//Do analysis and fill aods
	//Search for photons in fCalorimeter 
	
	TClonesArray * pl = new TClonesArray; 

	//Get vertex for photon momentum calculation
	Double_t vertex[]={0,0,0} ; //vertex ;
	if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);

	//Select the Calorimeter of the photon
	if(fCalorimeter == "PHOS")
		pl = GetAODPHOS();
	else if (fCalorimeter == "EMCAL")
		pl = GetAODEMCAL();
	//Fill AODCaloClusters and AODParticle with PHOS aods
	TLorentzVector mom ;
	for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
		AliAODCaloCluster * calo =  (AliAODCaloCluster*) (pl->At(icalo));	

		//Cluster selection, not charged, with photon id and in fidutial cut
		//Get Momentum vector, 
		calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
		//If too small or big pt, skip it
		if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ; 
		//Check acceptance selection
		if(IsFidutialCutOn()){
			Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter) ;
			if(! in ) continue ;
		}
		
		//Create AOD for analysis
		AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
		aodph.SetLabel(calo->GetLabel(0));
		//printf("Index %d, Id %d\n",icalo, calo->GetID());
		//Set the indeces of the original caloclusters  
		aodph.SetCaloLabel(calo->GetID(),-1);
		aodph.SetDetector(fCalorimeter);
		if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: Min pt cut and fidutial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta());	
			
		//Check Distance to Bad channel, set bit.
		Double_t distBad=calo->GetDistToBadChannel() ; //Distance to bad channel
		if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
		if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
		continue ;
		
		if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: Bad channel cut passed %4.2f\n",distBad);
				
		if(distBad > fMinDist3) aodph.SetDistToBad(2) ;
		else if(distBad > fMinDist2) aodph.SetDistToBad(1) ; 
		else aodph.SetDistToBad(0) ;
		
		//Check PID
		//PID selection or bit setting
		if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
			//Get most probable PID, check PID weights (in MC this option is mandatory)
			aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
			if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: PDG of identified particle %d\n",aodph.GetPdg());
			//If primary is not photon, skip it.
			if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
		}					
		else if(IsCaloPIDOn()){
			//Skip matched clusters with tracks
			if(calo->GetNTracksMatched() > 0) continue ;
		
			//Get most probable PID, 2 options check PID weights 
			//or redo PID, recommended option for EMCal.		
			if(!IsCaloPIDRecalculationOn())
				aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
			else
				aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
				
			if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: PDG of identified particle %d\n",aodph.GetPdg());
			
			//If cluster does not pass pid, not photon, skip it.
			if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;			
				
		}
		else{
			//Set PID bits for later selection (AliAnaPi0 for example)
			//GetPDG already called in SetPIDBits.
			GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph);
			if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: PID Bits set \n");		
		}
		
		if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
		
		//Play with the MC stack if available
		//Check origin of the candidates
		if(IsDataMC()){
			aodph.SetTag(GetCaloPID()->CheckOrigin(calo->GetLabel(0),GetMCStack()));
			if(GetDebug() > 0) printf("AliAnaPhoton::FillAOD: Origin of candidate %d\n",aodph.GetTag());
		}//Work with stack also   
		
		//Add AOD with photon object to aod branch
		AddAODParticle(aodph);
	
	}//loop

if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: End fill AODs \n");  

}

//__________________________________________________________________
void  AliAnaPhoton::MakeAnalysisFillHistograms() 
{
	//Do analysis and fill histograms

	//Get vertex for photon momentum calculation
	Double_t v[]={0,0,0} ; //vertex ;
	//if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) 
	GetReader()->GetVertex(v);

	//Loop on stored AOD photons
	Int_t naod = GetOutputAODBranch()->GetEntriesFast();
	if(GetDebug() > 0) printf("AliAnaPhoton::FillHistos: aod branch entries %d\n", naod);

	for(Int_t iaod = 0; iaod < naod ; iaod++){
		AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
		Int_t pdg = ph->GetPdg();

		if(GetDebug() > 3) printf("AliAnaPhoton::FillHistos: PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
		
		//If PID used, fill histos with photons in Calorimeter fCalorimeter
		if(IsCaloPIDOn())  
			if(pdg != AliCaloPID::kPhoton) continue; 
		if(ph->GetDetector() != fCalorimeter) continue;
		
		if(GetDebug() > 1) printf("AliAnaPhoton::FillHistos: ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
		
		//Fill photon histograms 
		Float_t ptcluster = ph->Pt();
		Float_t phicluster = ph->Phi();
		Float_t etacluster = ph->Eta();
		
		fhPtPhoton   ->Fill(ptcluster);
		fhPhiPhoton ->Fill(ptcluster,phicluster);
		fhEtaPhoton ->Fill(ptcluster,etacluster);
		
		if(IsDataMC()){
			Int_t tag =ph->GetTag();
			if(tag == AliCaloPID::kMCPrompt){
				fhPtPrompt   ->Fill(ptcluster);
				fhPhiPrompt ->Fill(ptcluster,phicluster);
				fhEtaPrompt ->Fill(ptcluster,etacluster);
			}
			else if(tag==AliCaloPID::kMCFragmentation)
			{
				fhPtFragmentation   ->Fill(ptcluster);
				fhPhiFragmentation ->Fill(ptcluster,phicluster);
				fhEtaFragmentation ->Fill(ptcluster,etacluster);
			}
			else if(tag==AliCaloPID::kMCPi0Decay)
			{
				fhPtPi0Decay   ->Fill(ptcluster);
				fhPhiPi0Decay ->Fill(ptcluster,phicluster);
				fhEtaPi0Decay ->Fill(ptcluster,etacluster);
			}
			else if(tag==AliCaloPID::kMCEtaDecay || tag==AliCaloPID::kMCOtherDecay)
			{
				fhPtOtherDecay   ->Fill(ptcluster);
				fhPhiOtherDecay ->Fill(ptcluster,phicluster);
				fhEtaOtherDecay ->Fill(ptcluster,etacluster);
			}
			else if(tag==AliCaloPID::kMCConversion)
			{
				fhPtConversion   ->Fill(ptcluster);
				fhPhiConversion ->Fill(ptcluster,phicluster);
				fhEtaConversion ->Fill(ptcluster,etacluster);
			}
			else{
				
				fhPtUnknown   ->Fill(ptcluster);
				fhPhiUnknown ->Fill(ptcluster,phicluster);
				fhEtaUnknown ->Fill(ptcluster,etacluster);
			}
		}//Histograms with MC
		
	}// aod loop

}


//__________________________________________________________________
void AliAnaPhoton::Print(const Option_t * opt) const
{
	//Print some relevant parameters set for the analysis
	
	if(! opt)
		return;

	printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
	AliAnaPartCorrBaseClass::Print(" ");
    printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
	printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
	printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
	printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
	printf("    \n") ;
	
}
Commit	Line	Data
1c5acb87	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 hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15	/* $Id: AliAnaPhoton.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
	16
	17	//_________________________________________________________________________
	18	//
	19	// Class for the photon identification.
	20	// Clusters from calorimeters are identified as photons
	21	// and kept in the AOD. Few histograms produced.
	22	//
	23	// -- Author: Gustavo Conesa (LNF-INFN)
	24	//////////////////////////////////////////////////////////////////////////////
	25
	26
	27	// --- ROOT system ---
	28	#include <TH2F.h>
	29
	30	// --- Analysis system ---
	31	#include "AliAnaPhoton.h"
	32	#include "AliCaloTrackReader.h"
	33	#include "AliCaloPID.h"
	34
	35	ClassImp(AliAnaPhoton)
	36
	37	//____________________________________________________________________________
	38	AliAnaPhoton::AliAnaPhoton() :
	39	AliAnaPartCorrBaseClass(), fCalorimeter(""),
	40	fMinDist(0.),fMinDist2(0.),fMinDist3(0.),
	41	fhPtPhoton(0),fhPhiPhoton(0),fhEtaPhoton(0),
	42	//MC
	43	fhPtPrompt(0),fhPhiPrompt(0),fhEtaPrompt(0),
	44	fhPtFragmentation(0),fhPhiFragmentation(0),fhEtaFragmentation(0),
	45	fhPtPi0Decay(0),fhPhiPi0Decay(0),fhEtaPi0Decay(0),
	46	fhPtOtherDecay(0),fhPhiOtherDecay(0),fhEtaOtherDecay(0),
	47	fhPtConversion(0),fhPhiConversion(0),fhEtaConversion(0),
	48	fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0)
	49	{
	50	//default ctor
	51
	52	//Initialize parameters
	53	InitParameters();
	54
	55	}
	56
	57	//____________________________________________________________________________
	58	AliAnaPhoton::AliAnaPhoton(const AliAnaPhoton & g) :
	59	AliAnaPartCorrBaseClass(g), fCalorimeter(g.fCalorimeter),
	60	fMinDist(g.fMinDist),fMinDist2(g.fMinDist2), fMinDist3(g.fMinDist3),
	61	fhPtPhoton(g.fhPtPhoton),fhPhiPhoton(g.fhPhiPhoton),fhEtaPhoton(g.fhEtaPhoton),
	62	//MC
	63	fhPtPrompt(g.fhPtPrompt),fhPhiPrompt(g.fhPhiPrompt),fhEtaPrompt(g.fhEtaPrompt),
	64	fhPtFragmentation(g.fhPtFragmentation),fhPhiFragmentation(g.fhPhiFragmentation),fhEtaFragmentation(g.fhEtaFragmentation),
65	fhPtPi0Decay(g.fhPtPi0Decay),fhPhiPi0Decay(g.fhPhiPi0Decay),fhEtaPi0Decay(g.fhEtaPi0Decay),
66	fhPtOtherDecay(g.fhPtOtherDecay),fhPhiOtherDecay(g.fhPhiOtherDecay),fhEtaOtherDecay(g.fhEtaOtherDecay),
67	fhPtConversion(g. fhPtConversion),fhPhiConversion(g.fhPhiConversion),fhEtaConversion(g.fhEtaConversion),
68	fhPtUnknown(g.fhPtUnknown),fhPhiUnknown(g.fhPhiUnknown),fhEtaUnknown(g.fhEtaUnknown)
69	{
70	// cpy ctor
71
72	}
73
74	//_________________________________________________________________________
75	AliAnaPhoton & AliAnaPhoton::operator = (const AliAnaPhoton & g)
76	{
77	// assignment operator
78
79	if(&g == this) return *this;
80
81	fCalorimeter = g.fCalorimeter ;
82	fMinDist = g.fMinDist;
83	fMinDist2 = g.fMinDist2;
84	fMinDist3 = g.fMinDist3;
85
86	fhPtPhoton = g.fhPtPhoton ;
87	fhPhiPhoton = g.fhPhiPhoton ;
88	fhEtaPhoton = g.fhEtaPhoton ;
89
90	fhPtPrompt = g.fhPtPrompt;
91	fhPhiPrompt = g.fhPhiPrompt;
92	fhEtaPrompt = g.fhEtaPrompt;
93	fhPtFragmentation = g.fhPtFragmentation;
94	fhPhiFragmentation = g.fhPhiFragmentation;
95	fhEtaFragmentation = g.fhEtaFragmentation;
96	fhPtPi0Decay = g.fhPtPi0Decay;
97	fhPhiPi0Decay = g.fhPhiPi0Decay;
98	fhEtaPi0Decay = g.fhEtaPi0Decay;
99	fhPtOtherDecay = g.fhPtOtherDecay;
100	fhPhiOtherDecay = g.fhPhiOtherDecay;
101	fhEtaOtherDecay = g.fhEtaOtherDecay;
102	fhPtConversion = g. fhPtConversion;
103	fhPhiConversion = g.fhPhiConversion;
104	fhEtaConversion = g.fhEtaConversion;
105	fhPtUnknown = g.fhPtUnknown;
106	fhPhiUnknown = g.fhPhiUnknown;
107	fhEtaUnknown = g.fhEtaUnknown;
108
109	return *this;
110
111	}
112
113	//____________________________________________________________________________
114	AliAnaPhoton::~AliAnaPhoton()
115	{
116	//dtor
117
118	}
119
120
121	//________________________________________________________________________
122	TList * AliAnaPhoton::GetCreateOutputObjects()
123	{
124	// Create histograms to be saved in output file and
125	// store them in outputContainer
126	TList * outputContainer = new TList() ;
127	outputContainer->SetName("PhotonHistos") ;
128
129	Int_t nptbins = GetHistoNPtBins();
130	Int_t nphibins = GetHistoNPhiBins();
131	Int_t netabins = GetHistoNEtaBins();
132	Float_t ptmax = GetHistoPtMax();
133	Float_t phimax = GetHistoPhiMax();
134	Float_t etamax = GetHistoEtaMax();
135	Float_t ptmin = GetHistoPtMin();
136	Float_t phimin = GetHistoPhiMin();
137	Float_t etamin = GetHistoEtaMin();
138
139	//Histograms of highest Photon identified in Event
140	fhPtPhoton = new TH1F("hPtPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
141	fhPtPhoton->SetYTitle("N");
142	fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
143	outputContainer->Add(fhPtPhoton) ;
144
145	fhPhiPhoton = new TH2F
146	("hPhiPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
147	fhPhiPhoton->SetYTitle("#phi");
148	fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
149	outputContainer->Add(fhPhiPhoton) ;
150
151	fhEtaPhoton = new TH2F
152	("hEtaPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
153	fhEtaPhoton->SetYTitle("#eta");
154	fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
155	outputContainer->Add(fhEtaPhoton) ;
156
157	if(IsDataMC()){
158
159	fhPtPrompt = new TH1F("hPtPrompt","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
160	fhPtPrompt->SetYTitle("N");
161	fhPtPrompt->SetXTitle("p_{T #gamma}(GeV/c)");
162	outputContainer->Add(fhPtPrompt) ;
163
164	fhPhiPrompt = new TH2F
165	("hPhiPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
166	fhPhiPrompt->SetYTitle("#phi");
167	fhPhiPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
168	outputContainer->Add(fhPhiPrompt) ;
169
170	fhEtaPrompt = new TH2F
171	("hEtaPrompt","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
172	fhEtaPrompt->SetYTitle("#eta");
173	fhEtaPrompt->SetXTitle("p_{T #gamma} (GeV/c)");
174	outputContainer->Add(fhEtaPrompt) ;
175
176	fhPtFragmentation = new TH1F("hPtFragmentation","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
177	fhPtFragmentation->SetYTitle("N");
178	fhPtFragmentation->SetXTitle("p_{T #gamma}(GeV/c)");
179	outputContainer->Add(fhPtFragmentation) ;
180
181	fhPhiFragmentation = new TH2F
182	("hPhiFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
183	fhPhiFragmentation->SetYTitle("#phi");
184	fhPhiFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
185	outputContainer->Add(fhPhiFragmentation) ;
186
187	fhEtaFragmentation = new TH2F
188	("hEtaFragmentation","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
189	fhEtaFragmentation->SetYTitle("#eta");
190	fhEtaFragmentation->SetXTitle("p_{T #gamma} (GeV/c)");
191	outputContainer->Add(fhEtaFragmentation) ;
192
193	fhPtPi0Decay = new TH1F("hPtPi0Decay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
194	fhPtPi0Decay->SetYTitle("N");
195	fhPtPi0Decay->SetXTitle("p_{T #gamma}(GeV/c)");
196	outputContainer->Add(fhPtPi0Decay) ;
197
198	fhPhiPi0Decay = new TH2F
199	("hPhiPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
200	fhPhiPi0Decay->SetYTitle("#phi");
201	fhPhiPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
202	outputContainer->Add(fhPhiPi0Decay) ;
203
204	fhEtaPi0Decay = new TH2F
205	("hEtaPi0Decay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
206	fhEtaPi0Decay->SetYTitle("#eta");
207	fhEtaPi0Decay->SetXTitle("p_{T #gamma} (GeV/c)");
208	outputContainer->Add(fhEtaPi0Decay) ;
209
210	fhPtOtherDecay = new TH1F("hPtOtherDecay","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
211	fhPtOtherDecay->SetYTitle("N");
212	fhPtOtherDecay->SetXTitle("p_{T #gamma}(GeV/c)");
213	outputContainer->Add(fhPtOtherDecay) ;
214
215	fhPhiOtherDecay = new TH2F
216	("hPhiOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
217	fhPhiOtherDecay->SetYTitle("#phi");
218	fhPhiOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
219	outputContainer->Add(fhPhiOtherDecay) ;
220
221	fhEtaOtherDecay = new TH2F
222	("hEtaOtherDecay","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
223	fhEtaOtherDecay->SetYTitle("#eta");
224	fhEtaOtherDecay->SetXTitle("p_{T #gamma} (GeV/c)");
225	outputContainer->Add(fhEtaOtherDecay) ;
226
227	fhPtConversion = new TH1F("hPtConversion","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
228	fhPtConversion->SetYTitle("N");
229	fhPtConversion->SetXTitle("p_{T #gamma}(GeV/c)");
230	outputContainer->Add(fhPtConversion) ;
231
232	fhPhiConversion = new TH2F
233	("hPhiConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
234	fhPhiConversion->SetYTitle("#phi");
235	fhPhiConversion->SetXTitle("p_{T #gamma} (GeV/c)");
236	outputContainer->Add(fhPhiConversion) ;
237
238	fhEtaConversion = new TH2F
239	("hEtaConversion","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
240	fhEtaConversion->SetYTitle("#eta");
241	fhEtaConversion->SetXTitle("p_{T #gamma} (GeV/c)");
242	outputContainer->Add(fhEtaConversion) ;
243
244	fhPtUnknown = new TH1F("hPtUnknown","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
245	fhPtUnknown->SetYTitle("N");
246	fhPtUnknown->SetXTitle("p_{T #gamma}(GeV/c)");
247	outputContainer->Add(fhPtUnknown) ;
248
249	fhPhiUnknown = new TH2F
250	("hPhiUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
251	fhPhiUnknown->SetYTitle("#phi");
252	fhPhiUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
253	outputContainer->Add(fhPhiUnknown) ;
254
255	fhEtaUnknown = new TH2F
256	("hEtaUnknown","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
257	fhEtaUnknown->SetYTitle("#eta");
258	fhEtaUnknown->SetXTitle("p_{T #gamma} (GeV/c)");
259	outputContainer->Add(fhEtaUnknown) ;
260	}//Histos with MC
261
262	//Save parameters used for analysis
263	TString parList ; //this will be list of parameters used for this analysis.
264	char onePar[255] ;
265
266	sprintf(onePar,"--- AliAnaPhoton ---\n") ;
267	parList+=onePar ;
268	sprintf(onePar,"Calorimeter: %s\n",fCalorimeter.Data()) ;
269	parList+=onePar ;
270	sprintf(onePar,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
271	parList+=onePar ;
272	sprintf(onePar,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
273	parList+=onePar ;
274	sprintf(onePar,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
275	parList+=onePar ;
276
277	//Get parameters set in base class.
278	parList += GetBaseParametersList() ;
279
280	//Get parameters set in PID class.
281	parList += GetCaloPID()->GetPIDParametersList() ;
282
283	//Get parameters set in FidutialCut class (not available yet)
284	//parlist += GetCaloPID()->GetFidCutParametersList()
285
286	TObjString *oString= new TObjString(parList) ;
287	outputContainer->Add(oString);
288
289	return outputContainer ;
290
291	}
292
293	//____________________________________________________________________________
294	void AliAnaPhoton::InitParameters()
295	{
296
297	//Initialize the parameters of the analysis.
298	SetOutputAODClassName("AliAODPWG4Particle");
299	SetOutputAODName("photons");
300	fCalorimeter = "PHOS" ;
301	fMinDist = 2.;
302	fMinDist2 = 4.;
303	fMinDist3 = 5.;
304
305	}
306
307	//__________________________________________________________________
308	void AliAnaPhoton::MakeAnalysisFillAOD()
309	{
310	//Do analysis and fill aods
311	//Search for photons in fCalorimeter
312
313	TClonesArray * pl = new TClonesArray;
314
315	//Get vertex for photon momentum calculation
316	Double_t vertex[]={0,0,0} ; //vertex ;
317	if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC) GetReader()->GetVertex(vertex);
318
319	//Select the Calorimeter of the photon
320	if(fCalorimeter == "PHOS")
321	pl = GetAODPHOS();
322	else if (fCalorimeter == "EMCAL")
323	pl = GetAODEMCAL();
324	//Fill AODCaloClusters and AODParticle with PHOS aods
325	TLorentzVector mom ;
326	for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
327	AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo));
328
329	//Cluster selection, not charged, with photon id and in fidutial cut
330	//Get Momentum vector,
331	calo->GetMomentum(mom,vertex);//Assume that come from vertex in straight line
332	//If too small or big pt, skip it
333	if(mom.Pt() < GetMinPt() \|\| mom.Pt() > GetMaxPt() ) continue ;
334	//Check acceptance selection
335	if(IsFidutialCutOn()){
336	Bool_t in = GetFidutialCut()->IsInFidutialCut(mom,fCalorimeter) ;
337	if(! in ) continue ;
338	}
339
340	//Create AOD for analysis
341	AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
342	aodph.SetLabel(calo->GetLabel(0));
343	//printf("Index %d, Id %d\n",icalo, calo->GetID());
344	//Set the indeces of the original caloclusters
345	aodph.SetCaloLabel(calo->GetID(),-1);
346	aodph.SetDetector(fCalorimeter);
347	if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: Min pt cut and fidutial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta());
348
349	//Check Distance to Bad channel, set bit.
350	Double_t distBad=calo->GetDistToBadChannel() ; //Distance to bad channel
351	if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
352	if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
353	continue ;
354
355	if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: Bad channel cut passed %4.2f\n",distBad);
356
357	if(distBad > fMinDist3) aodph.SetDistToBad(2) ;
358	else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
359	else aodph.SetDistToBad(0) ;
360
361	//Check PID
362	//PID selection or bit setting
363	if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
364	//Get most probable PID, check PID weights (in MC this option is mandatory)
365	aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
366	if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: PDG of identified particle %d\n",aodph.GetPdg());
367	//If primary is not photon, skip it.
368	if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
369	}
370	else if(IsCaloPIDOn()){
371	//Skip matched clusters with tracks
372	if(calo->GetNTracksMatched() > 0) continue ;
373
374	//Get most probable PID, 2 options check PID weights
375	//or redo PID, recommended option for EMCal.
376	if(!IsCaloPIDRecalculationOn())
377	aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
378	else
379	aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
380
381	if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: PDG of identified particle %d\n",aodph.GetPdg());
382
383	//If cluster does not pass pid, not photon, skip it.
384	if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
385
386	}
387	else{
388	//Set PID bits for later selection (AliAnaPi0 for example)
389	//GetPDG already called in SetPIDBits.
390	GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph);
391	if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: PID Bits set \n");
392	}
393
394	if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
395
396	//Play with the MC stack if available
397	//Check origin of the candidates
398	if(IsDataMC()){
399	aodph.SetTag(GetCaloPID()->CheckOrigin(calo->GetLabel(0),GetMCStack()));
400	if(GetDebug() > 0) printf("AliAnaPhoton::FillAOD: Origin of candidate %d\n",aodph.GetTag());
401	}//Work with stack also
402
403	//Add AOD with photon object to aod branch
404	AddAODParticle(aodph);
405
406	}//loop
407
408	if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: End fill AODs \n");
409
410	}
411
412	//__________________________________________________________________
413	void AliAnaPhoton::MakeAnalysisFillHistograms()
414	{
415	//Do analysis and fill histograms
416
417	//Get vertex for photon momentum calculation
418	Double_t v[]={0,0,0} ; //vertex ;
419	//if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC)
420	GetReader()->GetVertex(v);
421
422	//Loop on stored AOD photons
423	Int_t naod = GetOutputAODBranch()->GetEntriesFast();
424	if(GetDebug() > 0) printf("AliAnaPhoton::FillHistos: aod branch entries %d\n", naod);
425
426	for(Int_t iaod = 0; iaod < naod ; iaod++){
427	AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
428	Int_t pdg = ph->GetPdg();
429
430	if(GetDebug() > 3) printf("AliAnaPhoton::FillHistos: PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
431
432	//If PID used, fill histos with photons in Calorimeter fCalorimeter
433	if(IsCaloPIDOn())
434	if(pdg != AliCaloPID::kPhoton) continue;
435	if(ph->GetDetector() != fCalorimeter) continue;
436
437	if(GetDebug() > 1) printf("AliAnaPhoton::FillHistos: ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
438
439	//Fill photon histograms
440	Float_t ptcluster = ph->Pt();
441	Float_t phicluster = ph->Phi();
442	Float_t etacluster = ph->Eta();
443
444	fhPtPhoton ->Fill(ptcluster);
445	fhPhiPhoton ->Fill(ptcluster,phicluster);
446	fhEtaPhoton ->Fill(ptcluster,etacluster);
447
448	if(IsDataMC()){
449	Int_t tag =ph->GetTag();
450	if(tag == AliCaloPID::kMCPrompt){
451	fhPtPrompt ->Fill(ptcluster);
452	fhPhiPrompt ->Fill(ptcluster,phicluster);
453	fhEtaPrompt ->Fill(ptcluster,etacluster);
454	}
455	else if(tag==AliCaloPID::kMCFragmentation)
456	{
457	fhPtFragmentation ->Fill(ptcluster);
458	fhPhiFragmentation ->Fill(ptcluster,phicluster);
459	fhEtaFragmentation ->Fill(ptcluster,etacluster);
460	}
461	else if(tag==AliCaloPID::kMCPi0Decay)
462	{
463	fhPtPi0Decay ->Fill(ptcluster);
464	fhPhiPi0Decay ->Fill(ptcluster,phicluster);
465	fhEtaPi0Decay ->Fill(ptcluster,etacluster);
466	}
467	else if(tag==AliCaloPID::kMCEtaDecay \|\| tag==AliCaloPID::kMCOtherDecay)
468	{
469	fhPtOtherDecay ->Fill(ptcluster);
470	fhPhiOtherDecay ->Fill(ptcluster,phicluster);
471	fhEtaOtherDecay ->Fill(ptcluster,etacluster);
472	}
473	else if(tag==AliCaloPID::kMCConversion)
474	{
475	fhPtConversion ->Fill(ptcluster);
476	fhPhiConversion ->Fill(ptcluster,phicluster);
477	fhEtaConversion ->Fill(ptcluster,etacluster);
478	}
479	else{
480
481	fhPtUnknown ->Fill(ptcluster);
482	fhPhiUnknown ->Fill(ptcluster,phicluster);
483	fhEtaUnknown ->Fill(ptcluster,etacluster);
484	}
485	}//Histograms with MC
486
487	}// aod loop
488
489	}
490
491
492	//__________________________________________________________________
493	void AliAnaPhoton::Print(const Option_t * opt) const
494	{
495	//Print some relevant parameters set for the analysis
496
497	if(! opt)
498	return;
499
500	printf("** Print %s %s **\n", GetName(), GetTitle() ) ;
501	AliAnaPartCorrBaseClass::Print(" ");
502	printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
503	printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
504	printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
505	printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
506	printf(" \n") ;
507
508	}