[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 += GetFidCut()->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(""),
7175a03a	40	fMinDist(0.),fMinDist2(0.),fMinDist3(0.),
	41	fhPtPhoton(0),fhPhiPhoton(0),fhEtaPhoton(0),
	42	//MC
1c5acb87	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
7175a03a	146	("hPhiPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1c5acb87	147	fhPhiPhoton->SetYTitle("#phi");
	148	fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
	149	outputContainer->Add(fhPhiPhoton) ;
	150
	151	fhEtaPhoton = new TH2F
7175a03a	152	("hEtaPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1c5acb87	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)
7175a03a	284	//parlist += GetFidCut()->GetFidCutParametersList()
1c5acb87	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);
7175a03a	347	if(GetDebug() > 1)
7175a03a	348	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());
1c5acb87	349
	350	//Check Distance to Bad channel, set bit.
	351	Double_t distBad=calo->GetDistToBadChannel() ; //Distance to bad channel
	352	if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
	353	if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
	354	continue ;
	355
	356	if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: Bad channel cut passed %4.2f\n",distBad);
	357
	358	if(distBad > fMinDist3) aodph.SetDistToBad(2) ;
	359	else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
	360	else aodph.SetDistToBad(0) ;
	361
	362	//Check PID
	363	//PID selection or bit setting
	364	if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
	365	//Get most probable PID, check PID weights (in MC this option is mandatory)
	366	aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
7175a03a	367	if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: PDG of identified particle %d\n",aodph.GetPdg());
1c5acb87	368	//If primary is not photon, skip it.
	369	if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
	370	}
	371	else if(IsCaloPIDOn()){
	372	//Skip matched clusters with tracks
	373	if(calo->GetNTracksMatched() > 0) continue ;
	374
	375	//Get most probable PID, 2 options check PID weights
	376	//or redo PID, recommended option for EMCal.
	377	if(!IsCaloPIDRecalculationOn())
	378	aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->PID(),mom.E()));//PID with weights
	379	else
	380	aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
	381
	382	if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: PDG of identified particle %d\n",aodph.GetPdg());
	383
	384	//If cluster does not pass pid, not photon, skip it.
	385	if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
	386
	387	}
	388	else{
	389	//Set PID bits for later selection (AliAnaPi0 for example)
	390	//GetPDG already called in SetPIDBits.
	391	GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph);
	392	if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: PID Bits set \n");
	393	}
	394
	395	if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
	396
	397	//Play with the MC stack if available
	398	//Check origin of the candidates
	399	if(IsDataMC()){
	400	aodph.SetTag(GetCaloPID()->CheckOrigin(calo->GetLabel(0),GetMCStack()));
	401	if(GetDebug() > 0) printf("AliAnaPhoton::FillAOD: Origin of candidate %d\n",aodph.GetTag());
	402	}//Work with stack also
	403
	404	//Add AOD with photon object to aod branch
	405	AddAODParticle(aodph);
	406
	407	}//loop
	408
	409	if(GetDebug() > 1) printf("AliAnaPhoton::FillAOD: End fill AODs \n");
	410
	411	}
	412
	413	//__________________________________________________________________
	414	void AliAnaPhoton::MakeAnalysisFillHistograms()
	415	{
	416	//Do analysis and fill histograms
	417
	418	//Get vertex for photon momentum calculation
7175a03a	419	//Double_t v[]={0,0,0} ; //vertex ;
1c5acb87	420	//if(!GetReader()->GetDataType()== AliCaloTrackReader::kMC)
7175a03a	421	//GetReader()->GetVertex(v);
1c5acb87	422
	423	//Loop on stored AOD photons
	424	Int_t naod = GetOutputAODBranch()->GetEntriesFast();
	425	if(GetDebug() > 0) printf("AliAnaPhoton::FillHistos: aod branch entries %d\n", naod);
	426
	427	for(Int_t iaod = 0; iaod < naod ; iaod++){
	428	AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
	429	Int_t pdg = ph->GetPdg();
	430
7175a03a	431	if(GetDebug() > 3)
7175a03a	432	printf("AliAnaPhoton::FillHistos: PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
1c5acb87	433
	434	//If PID used, fill histos with photons in Calorimeter fCalorimeter
	435	if(IsCaloPIDOn())
	436	if(pdg != AliCaloPID::kPhoton) continue;
	437	if(ph->GetDetector() != fCalorimeter) continue;
	438
7175a03a	439	if(GetDebug() > 1)
7175a03a	440	printf("AliAnaPhoton::FillHistos: ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1c5acb87	441
	442	//Fill photon histograms
	443	Float_t ptcluster = ph->Pt();
	444	Float_t phicluster = ph->Phi();
	445	Float_t etacluster = ph->Eta();
	446
	447	fhPtPhoton ->Fill(ptcluster);
	448	fhPhiPhoton ->Fill(ptcluster,phicluster);
	449	fhEtaPhoton ->Fill(ptcluster,etacluster);
	450
	451	if(IsDataMC()){
	452	Int_t tag =ph->GetTag();
	453	if(tag == AliCaloPID::kMCPrompt){
	454	fhPtPrompt ->Fill(ptcluster);
	455	fhPhiPrompt ->Fill(ptcluster,phicluster);
	456	fhEtaPrompt ->Fill(ptcluster,etacluster);
	457	}
	458	else if(tag==AliCaloPID::kMCFragmentation)
	459	{
	460	fhPtFragmentation ->Fill(ptcluster);
	461	fhPhiFragmentation ->Fill(ptcluster,phicluster);
	462	fhEtaFragmentation ->Fill(ptcluster,etacluster);
	463	}
	464	else if(tag==AliCaloPID::kMCPi0Decay)
	465	{
	466	fhPtPi0Decay ->Fill(ptcluster);
	467	fhPhiPi0Decay ->Fill(ptcluster,phicluster);
	468	fhEtaPi0Decay ->Fill(ptcluster,etacluster);
	469	}
	470	else if(tag==AliCaloPID::kMCEtaDecay \|\| tag==AliCaloPID::kMCOtherDecay)
	471	{
	472	fhPtOtherDecay ->Fill(ptcluster);
	473	fhPhiOtherDecay ->Fill(ptcluster,phicluster);
	474	fhEtaOtherDecay ->Fill(ptcluster,etacluster);
	475	}
	476	else if(tag==AliCaloPID::kMCConversion)
	477	{
	478	fhPtConversion ->Fill(ptcluster);
	479	fhPhiConversion ->Fill(ptcluster,phicluster);
	480	fhEtaConversion ->Fill(ptcluster,etacluster);
	481	}
	482	else{
1c5acb87	483	fhPtUnknown ->Fill(ptcluster);
	484	fhPhiUnknown ->Fill(ptcluster,phicluster);
	485	fhEtaUnknown ->Fill(ptcluster,etacluster);
	486	}
	487	}//Histograms with MC
	488
	489	}// aod loop
	490
	491	}
	492
	493
	494	//__________________________________________________________________
	495	void AliAnaPhoton::Print(const Option_t * opt) const
	496	{
	497	//Print some relevant parameters set for the analysis
	498
	499	if(! opt)
	500	return;
	501
	502	printf("** Print %s %s **\n", GetName(), GetTitle() ) ;
	503	AliAnaPartCorrBaseClass::Print(" ");
	504	printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
	505	printf("Min Distance to Bad Channel = %2.1f\n",fMinDist);
	506	printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
	507	printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
	508	printf(" \n") ;
	509
	510	}