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

//_________________________________________________________________________
// Class to collect two-photon invariant mass distributions for
// extractin raw pi0 yield.
//
//-- Author: Dmitri Peressounko (RRC "KI") 
//-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
//-- and Gustavo Conesa (INFN-Frascati)
//_________________________________________________________________________


// --- ROOT system ---
#include "TH3.h"
//#include "Riostream.h"

//---- AliRoot system ----
#include "AliAnaPi0.h"
#include "AliCaloTrackReader.h"
#include "AliCaloPID.h"

#ifdef __PHOSGEO__
	#include "AliPHOSGeoUtils.h"
#endif

ClassImp(AliAnaPi0)

//________________________________________________________________________________________________________________________________________________  
AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
fNCentrBin(0),fNZvertBin(0),fNrpBin(0),
fNPID(0),fNmaxMixEv(0), fZvtxCut(0.),fCalorimeter(""),
fEventsList(0x0), fhEtalon(0x0),
fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0),
fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0)
{
//Default Ctor
 InitParameters();
 
}

//________________________________________________________________________________________________________________________________________________
AliAnaPi0::AliAnaPi0(const AliAnaPi0 & ex) : AliAnaPartCorrBaseClass(ex),  
fNCentrBin(ex.fNCentrBin),fNZvertBin(ex.fNZvertBin),fNrpBin(ex.fNrpBin),
fNPID(ex.fNPID),fNmaxMixEv(ex.fNmaxMixEv),fZvtxCut(ex.fZvtxCut), fCalorimeter(ex.fCalorimeter), 
fEventsList(ex.fEventsList), fhEtalon(ex.fhEtalon),
fhRe1(ex.fhRe1),fhMi1(ex.fhMi1),fhRe2(ex.fhRe2),fhMi2(ex.fhMi2),fhRe3(ex.fhRe3),fhMi3(ex.fhMi3),fhEvents(ex.fhEvents),
fhPrimPt(ex.fhPrimPt), fhPrimAccPt(ex.fhPrimAccPt), fhPrimY(ex.fhPrimY), 
fhPrimAccY(ex.fhPrimAccY), fhPrimPhi(ex.fhPrimPhi), fhPrimAccPhi(ex.fhPrimAccPhi)
{
	// cpy ctor
	//Do not need it
}

//________________________________________________________________________________________________________________________________________________
AliAnaPi0 & AliAnaPi0::operator = (const AliAnaPi0 & ex)
{
	// assignment operator
	
	if(this == &ex)return *this;
	((AliAnaPartCorrBaseClass *)this)->operator=(ex);
	
	fNCentrBin = ex.fNCentrBin  ; fNZvertBin = ex.fNZvertBin  ; fNrpBin = ex.fNrpBin  ; 
	fNPID = ex.fNPID  ; fNmaxMixEv = ex.fNmaxMixEv  ; fZvtxCut = ex.fZvtxCut  ;  fCalorimeter = ex.fCalorimeter  ;  
	fEventsList = ex.fEventsList  ;  fhEtalon = ex.fhEtalon  ; 
	fhRe1 = ex.fhRe1  ; fhMi1 = ex.fhMi1  ; fhRe2 = ex.fhRe2  ; fhMi2 = ex.fhMi2  ; 
	fhRe3 = ex.fhRe3  ; fhMi3 = ex.fhMi3  ; fhEvents = ex.fhEvents  ; 
	fhPrimPt = ex.fhPrimPt  ;  fhPrimAccPt = ex.fhPrimAccPt  ;  fhPrimY = ex.fhPrimY  ;  
	fhPrimAccY = ex.fhPrimAccY  ;  fhPrimPhi = ex.fhPrimPhi  ;  fhPrimAccPhi = ex.fhPrimAccPhi ;
	
	return *this;
	
}

//________________________________________________________________________________________________________________________________________________
AliAnaPi0::~AliAnaPi0() {
	// Remove event containers
	if(fEventsList){
		for(Int_t ic=0; ic<fNCentrBin; ic++){
			for(Int_t iz=0; iz<fNZvertBin; iz++){
				for(Int_t irp=0; irp<fNrpBin; irp++){
					fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp]->Delete() ;
					delete fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] ;
				}
			}
		}
		delete[] fEventsList; 
		fEventsList=0 ;
	}
	
#ifdef __PHOSGEO__
    if(fPHOSGeo) delete fPHOSGeo ;
#endif	
}

//________________________________________________________________________________________________________________________________________________
void AliAnaPi0::InitParameters()
{
//Init parameters when first called the analysis
//Set default parameters
	SetInputAODName("photons");
	fNCentrBin = 1;
	fNZvertBin = 1;
	fNrpBin    = 1;
    fNPID      = 9;
	fNmaxMixEv = 10;
	fZvtxCut   = 40;
	fCalorimeter  = "PHOS";
}
//________________________________________________________________________________________________________________________________________________
void AliAnaPi0::Init()
{  
//Init some data members needed in analysis

//Histograms binning and range
	if(!fhEtalon){                                                   //  p_T      alpha   d m_gg    
		fhEtalon = new TH3D("hEtalon","Histo with binning parameters",50,0.,25.,10,0.,1.,200,0.,1.) ; 
		fhEtalon->SetXTitle("P_{T} (GeV)") ;
		fhEtalon->SetYTitle("#alpha") ;
		fhEtalon->SetZTitle("m_{#gamma#gamma} (GeV)") ;
	}


//create event containers
	fEventsList = new TList*[fNCentrBin*fNZvertBin*fNrpBin] ;
	
	for(Int_t ic=0; ic<fNCentrBin; ic++){
		for(Int_t iz=0; iz<fNZvertBin; iz++){
			for(Int_t irp=0; irp<fNrpBin; irp++){
				fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] = new TList() ;
			}
		}
	}
	
		//If Geometry library loaded, do geometry selection during analysis.
#ifdef __PHOSGEO__
	printf("PHOS geometry initialized!\n");
	fPHOSGeo = new AliPHOSGeoUtils("PHOSgeo") ;
#endif	
	
}

//________________________________________________________________________________________________________________________________________________
TList * AliAnaPi0::GetCreateOutputObjects()
{  
	// Create histograms to be saved in output file and 
	// store them in fOutputContainer

	
	TList * outputContainer = new TList() ; 
	outputContainer->SetName(GetName()); 
	
	fhRe1=new TH3D*[fNCentrBin*fNPID] ;
	fhRe2=new TH3D*[fNCentrBin*fNPID] ;
	fhRe3=new TH3D*[fNCentrBin*fNPID] ;
	fhMi1=new TH3D*[fNCentrBin*fNPID] ;
	fhMi2=new TH3D*[fNCentrBin*fNPID] ;
	fhMi3=new TH3D*[fNCentrBin*fNPID] ;
	
	char key[255] ;
	char title[255] ;
	
	for(Int_t ic=0; ic<fNCentrBin; ic++){
		for(Int_t ipid=0; ipid<fNPID; ipid++){
					
			//Distance to bad module 1
			sprintf(key,"hRe_cen%d_pid%d_dist1",ic,ipid) ;
			sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
			
			fhEtalon->Clone(key);
			fhRe1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
			fhRe1[ic*fNPID+ipid]->SetName(key) ;
			fhRe1[ic*fNPID+ipid]->SetTitle(title) ;
			outputContainer->Add(fhRe1[ic*fNPID+ipid]) ;
			
			sprintf(key,"hMi_cen%d_pid%d_dist1",ic,ipid) ;
			sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
			fhMi1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
			fhMi1[ic*fNPID+ipid]->SetName(key) ;
			fhMi1[ic*fNPID+ipid]->SetTitle(title) ;
			outputContainer->Add(fhMi1[ic*fNPID+ipid]) ;
			
			//Distance to bad module 2
			sprintf(key,"hRe_cen%d_pid%d_dist2",ic,ipid) ;
			sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
			fhRe2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
			fhRe2[ic*fNPID+ipid]->SetName(key) ;
			fhRe2[ic*fNPID+ipid]->SetTitle(title) ;
			outputContainer->Add(fhRe2[ic*fNPID+ipid]) ;
			
			sprintf(key,"hMi_cen%d_pid%d_dist2",ic,ipid) ;
			sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
			fhMi2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
			fhMi2[ic*fNPID+ipid]->SetName(key) ;
			fhMi2[ic*fNPID+ipid]->SetTitle(title) ;
			outputContainer->Add(fhMi2[ic*fNPID+ipid]) ;
			
			//Distance to bad module 3
			sprintf(key,"hRe_cen%d_pid%d_dist3",ic,ipid) ;
			sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
			fhRe3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
			fhRe3[ic*fNPID+ipid]->SetName(key) ; 
			fhRe3[ic*fNPID+ipid]->SetTitle(title) ;
			outputContainer->Add(fhRe3[ic*fNPID+ipid]) ;
			
			sprintf(key,"hMi_cen%d_pid%d_dist3",ic,ipid) ;
			sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
			fhMi3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
			fhMi3[ic*fNPID+ipid]->SetName(key) ;
			fhMi3[ic*fNPID+ipid]->SetTitle(title) ;
			outputContainer->Add(fhMi3[ic*fNPID+ipid]) ;
		}
    }
	
	
	fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
					  fNZvertBin,0.,1.*fNZvertBin,fNrpBin,0.,1.*fNrpBin) ;
	outputContainer->Add(fhEvents) ;

	//Histograms filled only if MC data is requested 	
	if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
		if(fhEtalon->GetXaxis()->GetXbins() && fhEtalon->GetXaxis()->GetXbins()->GetSize()){ //Variable bin size
			fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
			fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",fhEtalon->GetXaxis()->GetNbins(),
								   fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
		}
		else{
			fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
			fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",
								   fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
		}
		outputContainer->Add(fhPrimPt) ;
		outputContainer->Add(fhPrimAccPt) ;
		
		fhPrimY = new TH1D("hPrimaryRapidity","Rapidity of primary pi0",100,-5.,5.) ; 
		outputContainer->Add(fhPrimY) ;
		
		fhPrimAccY = new TH1D("hPrimAccRapidity","Rapidity of primary pi0",100,-5.,5.) ; 
		outputContainer->Add(fhPrimAccY) ;
		
		fhPrimPhi = new TH1D("hPrimaryPhi","Azimithal of primary pi0",180,0.,360.) ; 
		outputContainer->Add(fhPrimPhi) ;
		
		fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",180,-0.,360.) ; 
		outputContainer->Add(fhPrimAccPhi) ;
	}
	
	//Save parameters used for analysis
	TString parList ; //this will be list of parameters used for this analysis.
	char onePar[255] ;
	sprintf(onePar,"--- AliAnaPi0 ---\n") ;
	parList+=onePar ;	
	sprintf(onePar,"Number of bins in Centrality:  %d \n",fNCentrBin) ;
	parList+=onePar ;
	sprintf(onePar,"Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
	parList+=onePar ;
	sprintf(onePar,"Number of bins in Reac. Plain: %d \n",fNrpBin) ;
	parList+=onePar ;
	sprintf(onePar,"Depth of event buffer: %d \n",fNmaxMixEv) ;
	parList+=onePar ;
	sprintf(onePar,"Number of different PID used:  %d \n",fNPID) ;
	parList+=onePar ;
	sprintf(onePar,"Cuts: \n") ;
	parList+=onePar ;
	sprintf(onePar,"Z vertex position: -%f < z < %f \n",fZvtxCut,fZvtxCut) ;
	parList+=onePar ;
	sprintf(onePar,"Calorimeter: %s \n",fCalorimeter.Data()) ;
	parList+=onePar ;
	
	TObjString *oString= new TObjString(parList) ;
	outputContainer->Add(oString);
	
	return outputContainer;
}

//_________________________________________________________________________________________________________________________________________________
void AliAnaPi0::Print(const Option_t * /*opt*/) const
{
	//Print some relevant parameters set for the analysis
	AliAnaPartCorrBaseClass::Print(" ");
	printf("Class AliAnaPi0 for gamma-gamma inv.mass construction \n") ;
	printf("Number of bins in Centrality:  %d \n",fNCentrBin) ;
	printf("Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
	printf("Number of bins in Reac. Plain: %d \n",fNrpBin) ;
	printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
	printf("Number of different PID used:  %d \n",fNPID) ;
	printf("Cuts: \n") ;
	printf("Z vertex position: -%2.3f < z < %2.3f \n",fZvtxCut,fZvtxCut) ;
	printf("------------------------------------------------------\n") ;
} 


//____________________________________________________________________________________________________________________________________________________
void AliAnaPi0:: MakeAnalysisFillHistograms() 
{
	//Process one event and extract photons from AOD branch 
	// filled with AliAnaPhoton and fill histos with invariant mass
	
	//Apply some cuts on event: vertex position and centrality range  
	Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
	if(IsBadRun(iRun)) return ;	
	
	Double_t vert[]={0,0,0} ; //vertex ;
	GetReader()->GetVertex(vert);
	if(vert[2]<-fZvtxCut || vert[2]> fZvtxCut) return ; //Event can not be used (vertex, centrality,... cuts not fulfilled)

	//Get Centrality and calculate centrality bin
	//Does not exist in ESD yet???????
	Int_t curCentrBin=0 ;
 
	//Get Reaction Plain position and calculate RP bin
	//does not exist in ESD yet????
	Int_t curRPBin=0 ;	
	
	Int_t curZvertBin=(Int_t)(0.5*fNZvertBin*(vert[2]+fZvtxCut)/fZvtxCut) ;
		
	fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
			
	Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
	if(GetDebug() > 1) printf("AliAnaPi0::FillHistos: photon entries %d\n", nPhot);
	
	for(Int_t i1=0; i1<nPhot-1; i1++){
		AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
		TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
		for(Int_t i2=i1+1; i2<nPhot; i2++){
			AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
			TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
			Double_t m  = (photon1 + photon2).M() ;
			Double_t pt = (photon1 + photon2).Pt();
			Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
			if(GetDebug() > 2)
				printf("AliAnaPi0::FillHistos: Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
						p1->Pt(), p2->Pt(), pt,m,a);
			for(Int_t ipid=0; ipid<fNPID; ipid++)
			{
				if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){ 
					fhRe1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
					if(p1->DistToBad()>0 && p2->DistToBad()>0){
						fhRe2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
						if(p1->DistToBad()>1 && p2->DistToBad()>1){
							fhRe3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
						}
					}
				}
			} 
		}
	}
		
	//Fill mixed
	
	TList * evMixList=fEventsList[curCentrBin*fNZvertBin*fNrpBin+curZvertBin*fNrpBin+curRPBin] ;
	Int_t nMixed = evMixList->GetSize() ;
	for(Int_t ii=0; ii<nMixed; ii++){  
		TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
		Int_t nPhot2=ev2->GetEntriesFast() ;
		Double_t m = -999;
		if(GetDebug() > 1) printf("AliAnaPi0::FillHistos: Mixed event %d photon entries %d\n", ii, nPhot);
		
		for(Int_t i1=0; i1<nPhot; i1++){
			AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
			TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
			for(Int_t i2=0; i2<nPhot2; i2++){
				AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
				
				TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
				m =           (photon1+photon2).M() ; 
				Double_t pt = (photon1 + photon2).Pt();
				Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
				if(GetDebug() > 2)
					printf("AliAnaPi0::FillHistos: Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
							p1->Pt(), p2->Pt(), pt,m,a);			
				for(Int_t ipid=0; ipid<fNPID; ipid++){ 
				if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){ 
						fhMi1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
						if(p1->DistToBad()>0 && p2->DistToBad()>0){
							fhMi2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
							if(p1->DistToBad()>1 && p2->DistToBad()>1){
								fhMi3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
							}
							
						}
					}
				}
			}
		}
	}
	
	TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
	//Add current event to buffer and Remove redandant events 
	if(currentEvent->GetEntriesFast()>0){
		evMixList->AddFirst(currentEvent) ;
		currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
		if(evMixList->GetSize()>=fNmaxMixEv)
		{
			TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
			evMixList->RemoveLast() ;
			delete tmp ;
		}
	} 
	else{ //empty event
		delete currentEvent ;
		currentEvent=0 ; 
	}
		
	//Acceptance
	AliStack * stack = GetMCStack();
	if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
		for(Int_t i=0 ; i<stack->GetNprimary(); i++){
			TParticle * prim = stack->Particle(i) ;
			if(prim->GetPdgCode() == 111){
				Double_t pi0Pt = prim->Pt() ;
				//printf("pi0, pt %2.2f\n",pi0Pt);
				Double_t pi0Y  = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
				Double_t phi   = TMath::RadToDeg()*prim->Phi() ;
				if(TMath::Abs(pi0Y) < 0.5){
					fhPrimPt->Fill(pi0Pt) ;
				}
				fhPrimY  ->Fill(pi0Y) ;
				fhPrimPhi->Fill(phi) ;
				
				//Check if both photons hit Calorimeter
				Int_t iphot1=prim->GetFirstDaughter() ;
				Int_t iphot2=prim->GetLastDaughter() ;
				if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
					TParticle * phot1 = stack->Particle(iphot1) ;
					TParticle * phot2 = stack->Particle(iphot2) ;
					if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
						//printf("2 photons: photon 1: pt %2.2f, phi %3.2f, eta %1.2f; photon 2: pt %2.2f, phi %3.2f, eta %1.2f\n",
						//	phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
						Bool_t inacceptance = kFALSE;
#ifdef __PHOSGEO__
						Int_t mod ;
						Double_t x,z ;
						if(fCalorimeter == "PHOS" && fPHOSGeo->ImpactOnEmc(phot1,mod,z,x) && fPHOSGeo->ImpactOnEmc(phot1,mod,z,x)) 
							inacceptance = kTRUE;
						//printf("In REAL PHOS acceptance? %d\n",inacceptance);
#else
						TLorentzVector lv1, lv2;
						phot1->Momentum(lv1);
						phot2->Momentum(lv2);
						if(GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter) && GetFidutialCut()->IsInFidutialCut(lv2,fCalorimeter)) 
							inacceptance = kTRUE ;
						//printf("In %s fidutial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
#endif							  			 			  
						if(inacceptance){
							fhPrimAccPt->Fill(pi0Pt) ;
							fhPrimAccPhi->Fill(phi) ;
							fhPrimAccY->Fill(pi0Y) ;
						}//Accepted
					}// 2 photons      
				}//Check daughters exist
			}// Primary pi0
		}//loop on primaries	
	}//stack exists and data is MC
	
}
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 is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15	/* $Id: $ */
	16
	17	//_________________________________________________________________________
	18	// Class to collect two-photon invariant mass distributions for
	19	// extractin raw pi0 yield.
	20	//
	21	//-- Author: Dmitri Peressounko (RRC "KI")
	22	//-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
	23	//-- and Gustavo Conesa (INFN-Frascati)
	24	//_________________________________________________________________________
	25
	26
	27	// --- ROOT system ---
	28	#include "TH3.h"
	29	//#include "Riostream.h"
	30
	31	//---- AliRoot system ----
	32	#include "AliAnaPi0.h"
	33	#include "AliCaloTrackReader.h"
	34	#include "AliCaloPID.h"
	35
	36	#ifdef __PHOSGEO__
	37	#include "AliPHOSGeoUtils.h"
	38	#endif
	39
	40	ClassImp(AliAnaPi0)
	41
	42	//________________________________________________________________________________________________________________________________________________
	43	AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
	44	fNCentrBin(0),fNZvertBin(0),fNrpBin(0),
	45	fNPID(0),fNmaxMixEv(0), fZvtxCut(0.),fCalorimeter(""),
	46	fEventsList(0x0), fhEtalon(0x0),
	47	fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0),
	48	fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0)
	49	{
	50	//Default Ctor
	51	InitParameters();
	52
	53	}
	54
	55	//________________________________________________________________________________________________________________________________________________
	56	AliAnaPi0::AliAnaPi0(const AliAnaPi0 & ex) : AliAnaPartCorrBaseClass(ex),
	57	fNCentrBin(ex.fNCentrBin),fNZvertBin(ex.fNZvertBin),fNrpBin(ex.fNrpBin),
	58	fNPID(ex.fNPID),fNmaxMixEv(ex.fNmaxMixEv),fZvtxCut(ex.fZvtxCut), fCalorimeter(ex.fCalorimeter),
	59	fEventsList(ex.fEventsList), fhEtalon(ex.fhEtalon),
	60	fhRe1(ex.fhRe1),fhMi1(ex.fhMi1),fhRe2(ex.fhRe2),fhMi2(ex.fhMi2),fhRe3(ex.fhRe3),fhMi3(ex.fhMi3),fhEvents(ex.fhEvents),
	61	fhPrimPt(ex.fhPrimPt), fhPrimAccPt(ex.fhPrimAccPt), fhPrimY(ex.fhPrimY),
	62	fhPrimAccY(ex.fhPrimAccY), fhPrimPhi(ex.fhPrimPhi), fhPrimAccPhi(ex.fhPrimAccPhi)
	63	{
	64	// cpy ctor
65	//Do not need it
66	}
67
68	//________________________________________________________________________________________________________________________________________________
69	AliAnaPi0 & AliAnaPi0::operator = (const AliAnaPi0 & ex)
70	{
71	// assignment operator
72
73	if(this == &ex)return *this;
74	((AliAnaPartCorrBaseClass *)this)->operator=(ex);
75
76	fNCentrBin = ex.fNCentrBin ; fNZvertBin = ex.fNZvertBin ; fNrpBin = ex.fNrpBin ;
77	fNPID = ex.fNPID ; fNmaxMixEv = ex.fNmaxMixEv ; fZvtxCut = ex.fZvtxCut ; fCalorimeter = ex.fCalorimeter ;
78	fEventsList = ex.fEventsList ; fhEtalon = ex.fhEtalon ;
79	fhRe1 = ex.fhRe1 ; fhMi1 = ex.fhMi1 ; fhRe2 = ex.fhRe2 ; fhMi2 = ex.fhMi2 ;
80	fhRe3 = ex.fhRe3 ; fhMi3 = ex.fhMi3 ; fhEvents = ex.fhEvents ;
81	fhPrimPt = ex.fhPrimPt ; fhPrimAccPt = ex.fhPrimAccPt ; fhPrimY = ex.fhPrimY ;
82	fhPrimAccY = ex.fhPrimAccY ; fhPrimPhi = ex.fhPrimPhi ; fhPrimAccPhi = ex.fhPrimAccPhi ;
83
84	return *this;
85
86	}
87
88	//________________________________________________________________________________________________________________________________________________
89	AliAnaPi0::~AliAnaPi0() {
90	// Remove event containers
91	if(fEventsList){
92	for(Int_t ic=0; ic<fNCentrBin; ic++){
93	for(Int_t iz=0; iz<fNZvertBin; iz++){
94	for(Int_t irp=0; irp<fNrpBin; irp++){
95	fEventsList[icfNZvertBinfNrpBin+iz*fNrpBin+irp]->Delete() ;
96	delete fEventsList[icfNZvertBinfNrpBin+iz*fNrpBin+irp] ;
97	}
98	}
99	}
100	delete[] fEventsList;
101	fEventsList=0 ;
102	}
103
104	#ifdef __PHOSGEO__
105	if(fPHOSGeo) delete fPHOSGeo ;
106	#endif
107	}
108
109	//________________________________________________________________________________________________________________________________________________
110	void AliAnaPi0::InitParameters()
111	{
112	//Init parameters when first called the analysis
113	//Set default parameters
114	SetInputAODName("photons");
115	fNCentrBin = 1;
116	fNZvertBin = 1;
117	fNrpBin = 1;
118	fNPID = 9;
119	fNmaxMixEv = 10;
120	fZvtxCut = 40;
121	fCalorimeter = "PHOS";
122	}
123	//________________________________________________________________________________________________________________________________________________
124	void AliAnaPi0::Init()
125	{
126	//Init some data members needed in analysis
127
128	//Histograms binning and range
129	if(!fhEtalon){ // p_T alpha d m_gg
130	fhEtalon = new TH3D("hEtalon","Histo with binning parameters",50,0.,25.,10,0.,1.,200,0.,1.) ;
131	fhEtalon->SetXTitle("P_{T} (GeV)") ;
132	fhEtalon->SetYTitle("#alpha") ;
133	fhEtalon->SetZTitle("m_{#gamma#gamma} (GeV)") ;
134	}
135
136
137	//create event containers
138	fEventsList = new TList[fNCentrBinfNZvertBin*fNrpBin] ;
139
140	for(Int_t ic=0; ic<fNCentrBin; ic++){
141	for(Int_t iz=0; iz<fNZvertBin; iz++){
142	for(Int_t irp=0; irp<fNrpBin; irp++){
143	fEventsList[icfNZvertBinfNrpBin+iz*fNrpBin+irp] = new TList() ;
144	}
145	}
146	}
147
148	//If Geometry library loaded, do geometry selection during analysis.
149	#ifdef __PHOSGEO__
150	printf("PHOS geometry initialized!\n");
151	fPHOSGeo = new AliPHOSGeoUtils("PHOSgeo") ;
152	#endif
153
154	}
155
156	//________________________________________________________________________________________________________________________________________________
157	TList * AliAnaPi0::GetCreateOutputObjects()
158	{
159	// Create histograms to be saved in output file and
160	// store them in fOutputContainer
161
162
163	TList * outputContainer = new TList() ;
164	outputContainer->SetName(GetName());
165
166	fhRe1=new TH3D[fNCentrBinfNPID] ;
167	fhRe2=new TH3D[fNCentrBinfNPID] ;
168	fhRe3=new TH3D[fNCentrBinfNPID] ;
169	fhMi1=new TH3D[fNCentrBinfNPID] ;
170	fhMi2=new TH3D[fNCentrBinfNPID] ;
171	fhMi3=new TH3D[fNCentrBinfNPID] ;
172
173	char key[255] ;
174	char title[255] ;
175
176	for(Int_t ic=0; ic<fNCentrBin; ic++){
177	for(Int_t ipid=0; ipid<fNPID; ipid++){
178
179	//Distance to bad module 1
180	sprintf(key,"hRe_cen%d_pid%d_dist1",ic,ipid) ;
181	sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
182
183	fhEtalon->Clone(key);
184	fhRe1[icfNPID+ipid]=(TH3D)fhEtalon->Clone(key) ;
185	fhRe1[ic*fNPID+ipid]->SetName(key) ;
186	fhRe1[ic*fNPID+ipid]->SetTitle(title) ;
187	outputContainer->Add(fhRe1[ic*fNPID+ipid]) ;
188
189	sprintf(key,"hMi_cen%d_pid%d_dist1",ic,ipid) ;
190	sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
191	fhMi1[icfNPID+ipid]=(TH3D)fhEtalon->Clone(key) ;
192	fhMi1[ic*fNPID+ipid]->SetName(key) ;
193	fhMi1[ic*fNPID+ipid]->SetTitle(title) ;
194	outputContainer->Add(fhMi1[ic*fNPID+ipid]) ;
195
196	//Distance to bad module 2
197	sprintf(key,"hRe_cen%d_pid%d_dist2",ic,ipid) ;
198	sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
199	fhRe2[icfNPID+ipid]=(TH3D)fhEtalon->Clone(key) ;
200	fhRe2[ic*fNPID+ipid]->SetName(key) ;
201	fhRe2[ic*fNPID+ipid]->SetTitle(title) ;
202	outputContainer->Add(fhRe2[ic*fNPID+ipid]) ;
203
204	sprintf(key,"hMi_cen%d_pid%d_dist2",ic,ipid) ;
205	sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
206	fhMi2[icfNPID+ipid]=(TH3D)fhEtalon->Clone(key) ;
207	fhMi2[ic*fNPID+ipid]->SetName(key) ;
208	fhMi2[ic*fNPID+ipid]->SetTitle(title) ;
209	outputContainer->Add(fhMi2[ic*fNPID+ipid]) ;
210
211	//Distance to bad module 3
212	sprintf(key,"hRe_cen%d_pid%d_dist3",ic,ipid) ;
213	sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
214	fhRe3[icfNPID+ipid]=(TH3D)fhEtalon->Clone(key) ;
215	fhRe3[ic*fNPID+ipid]->SetName(key) ;
216	fhRe3[ic*fNPID+ipid]->SetTitle(title) ;
217	outputContainer->Add(fhRe3[ic*fNPID+ipid]) ;
218
219	sprintf(key,"hMi_cen%d_pid%d_dist3",ic,ipid) ;
220	sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
221	fhMi3[icfNPID+ipid]=(TH3D)fhEtalon->Clone(key) ;
222	fhMi3[ic*fNPID+ipid]->SetName(key) ;
223	fhMi3[ic*fNPID+ipid]->SetTitle(title) ;
224	outputContainer->Add(fhMi3[ic*fNPID+ipid]) ;
225	}
226	}
227
228
229	fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
230	fNZvertBin,0.,1.fNZvertBin,fNrpBin,0.,1.fNrpBin) ;
231	outputContainer->Add(fhEvents) ;
232
233	//Histograms filled only if MC data is requested
234	if(IsDataMC() \|\| (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
235	if(fhEtalon->GetXaxis()->GetXbins() && fhEtalon->GetXaxis()->GetXbins()->GetSize()){ //Variable bin size
236	fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
237	fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",fhEtalon->GetXaxis()->GetNbins(),
238	fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
239	}
240	else{
241	fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
242	fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",
243	fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
244	}
245	outputContainer->Add(fhPrimPt) ;
246	outputContainer->Add(fhPrimAccPt) ;
247
248	fhPrimY = new TH1D("hPrimaryRapidity","Rapidity of primary pi0",100,-5.,5.) ;
249	outputContainer->Add(fhPrimY) ;
250
251	fhPrimAccY = new TH1D("hPrimAccRapidity","Rapidity of primary pi0",100,-5.,5.) ;
252	outputContainer->Add(fhPrimAccY) ;
253
254	fhPrimPhi = new TH1D("hPrimaryPhi","Azimithal of primary pi0",180,0.,360.) ;
255	outputContainer->Add(fhPrimPhi) ;
256
257	fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",180,-0.,360.) ;
258	outputContainer->Add(fhPrimAccPhi) ;
259	}
260
261	//Save parameters used for analysis
262	TString parList ; //this will be list of parameters used for this analysis.
263	char onePar[255] ;
264	sprintf(onePar,"--- AliAnaPi0 ---\n") ;
265	parList+=onePar ;
266	sprintf(onePar,"Number of bins in Centrality: %d \n",fNCentrBin) ;
267	parList+=onePar ;
268	sprintf(onePar,"Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
269	parList+=onePar ;
270	sprintf(onePar,"Number of bins in Reac. Plain: %d \n",fNrpBin) ;
271	parList+=onePar ;
272	sprintf(onePar,"Depth of event buffer: %d \n",fNmaxMixEv) ;
273	parList+=onePar ;
274	sprintf(onePar,"Number of different PID used: %d \n",fNPID) ;
275	parList+=onePar ;
276	sprintf(onePar,"Cuts: \n") ;
277	parList+=onePar ;
278	sprintf(onePar,"Z vertex position: -%f < z < %f \n",fZvtxCut,fZvtxCut) ;
279	parList+=onePar ;
280	sprintf(onePar,"Calorimeter: %s \n",fCalorimeter.Data()) ;
281	parList+=onePar ;
282
283	TObjString *oString= new TObjString(parList) ;
284	outputContainer->Add(oString);
285
286	return outputContainer;
287	}
288
289	//_________________________________________________________________________________________________________________________________________________
290	void AliAnaPi0::Print(const Option_t * /opt/) const
291	{
292	//Print some relevant parameters set for the analysis
293	AliAnaPartCorrBaseClass::Print(" ");
294	printf("Class AliAnaPi0 for gamma-gamma inv.mass construction \n") ;
295	printf("Number of bins in Centrality: %d \n",fNCentrBin) ;
296	printf("Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
297	printf("Number of bins in Reac. Plain: %d \n",fNrpBin) ;
298	printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
299	printf("Number of different PID used: %d \n",fNPID) ;
300	printf("Cuts: \n") ;
301	printf("Z vertex position: -%2.3f < z < %2.3f \n",fZvtxCut,fZvtxCut) ;
302	printf("------------------------------------------------------\n") ;
303	}
304
305
306	//____________________________________________________________________________________________________________________________________________________
307	void AliAnaPi0:: MakeAnalysisFillHistograms()
308	{
309	//Process one event and extract photons from AOD branch
310	// filled with AliAnaPhoton and fill histos with invariant mass
311
312	//Apply some cuts on event: vertex position and centrality range
313	Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
314	if(IsBadRun(iRun)) return ;
315
316	Double_t vert[]={0,0,0} ; //vertex ;
317	GetReader()->GetVertex(vert);
318	if(vert[2]<-fZvtxCut \|\| vert[2]> fZvtxCut) return ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
319
320	//Get Centrality and calculate centrality bin
321	//Does not exist in ESD yet???????
322	Int_t curCentrBin=0 ;
323
324	//Get Reaction Plain position and calculate RP bin
325	//does not exist in ESD yet????
326	Int_t curRPBin=0 ;
327
328	Int_t curZvertBin=(Int_t)(0.5fNZvertBin(vert[2]+fZvtxCut)/fZvtxCut) ;
329
330	fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
331
332	Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
333	if(GetDebug() > 1) printf("AliAnaPi0::FillHistos: photon entries %d\n", nPhot);
334
335	for(Int_t i1=0; i1<nPhot-1; i1++){
336	AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
337	TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
338	for(Int_t i2=i1+1; i2<nPhot; i2++){
339	AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
340	TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
341	Double_t m = (photon1 + photon2).M() ;
342	Double_t pt = (photon1 + photon2).Pt();
343	Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
344	if(GetDebug() > 2)
345	printf("AliAnaPi0::FillHistos: Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
346	p1->Pt(), p2->Pt(), pt,m,a);
347	for(Int_t ipid=0; ipid<fNPID; ipid++)
348	{
349	if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
350	fhRe1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
351	if(p1->DistToBad()>0 && p2->DistToBad()>0){
352	fhRe2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
353	if(p1->DistToBad()>1 && p2->DistToBad()>1){
354	fhRe3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
355	}
356	}
357	}
358	}
359	}
360	}
361
362	//Fill mixed
363
364	TList * evMixList=fEventsList[curCentrBinfNZvertBinfNrpBin+curZvertBin*fNrpBin+curRPBin] ;
365	Int_t nMixed = evMixList->GetSize() ;
366	for(Int_t ii=0; ii<nMixed; ii++){
367	TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
368	Int_t nPhot2=ev2->GetEntriesFast() ;
369	Double_t m = -999;
370	if(GetDebug() > 1) printf("AliAnaPi0::FillHistos: Mixed event %d photon entries %d\n", ii, nPhot);
371
372	for(Int_t i1=0; i1<nPhot; i1++){
373	AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
374	TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
375	for(Int_t i2=0; i2<nPhot2; i2++){
376	AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
377
378	TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
379	m = (photon1+photon2).M() ;
380	Double_t pt = (photon1 + photon2).Pt();
381	Double_t a = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
382	if(GetDebug() > 2)
383	printf("AliAnaPi0::FillHistos: Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
384	p1->Pt(), p2->Pt(), pt,m,a);
385	for(Int_t ipid=0; ipid<fNPID; ipid++){
386	if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
387	fhMi1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
388	if(p1->DistToBad()>0 && p2->DistToBad()>0){
389	fhMi2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
390	if(p1->DistToBad()>1 && p2->DistToBad()>1){
391	fhMi3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
392	}
393
394	}
395	}
396	}
397	}
398	}
399	}
400
401	TClonesArray currentEvent = new TClonesArray(GetInputAODBranch());
402	//Add current event to buffer and Remove redandant events
403	if(currentEvent->GetEntriesFast()>0){
404	evMixList->AddFirst(currentEvent) ;
405	currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
406	if(evMixList->GetSize()>=fNmaxMixEv)
407	{
408	TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
409	evMixList->RemoveLast() ;
410	delete tmp ;
411	}
412	}
413	else{ //empty event
414	delete currentEvent ;
415	currentEvent=0 ;
416	}
417
418	//Acceptance
419	AliStack * stack = GetMCStack();
420	if(stack && (IsDataMC() \|\| (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
421	for(Int_t i=0 ; i<stack->GetNprimary(); i++){
422	TParticle * prim = stack->Particle(i) ;
423	if(prim->GetPdgCode() == 111){
424	Double_t pi0Pt = prim->Pt() ;
425	//printf("pi0, pt %2.2f\n",pi0Pt);
426	Double_t pi0Y = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
427	Double_t phi = TMath::RadToDeg()*prim->Phi() ;
428	if(TMath::Abs(pi0Y) < 0.5){
429	fhPrimPt->Fill(pi0Pt) ;
430	}
431	fhPrimY ->Fill(pi0Y) ;
432	fhPrimPhi->Fill(phi) ;
433
434	//Check if both photons hit Calorimeter
435	Int_t iphot1=prim->GetFirstDaughter() ;
436	Int_t iphot2=prim->GetLastDaughter() ;
437	if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
438	TParticle * phot1 = stack->Particle(iphot1) ;
439	TParticle * phot2 = stack->Particle(iphot2) ;
440	if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
441	//printf("2 photons: photon 1: pt %2.2f, phi %3.2f, eta %1.2f; photon 2: pt %2.2f, phi %3.2f, eta %1.2f\n",
442	// phot1->Pt(), phot1->Phi()180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()180./3.1415, phot2->Eta());
443	Bool_t inacceptance = kFALSE;
444	#ifdef __PHOSGEO__
445	Int_t mod ;
446	Double_t x,z ;
447	if(fCalorimeter == "PHOS" && fPHOSGeo->ImpactOnEmc(phot1,mod,z,x) && fPHOSGeo->ImpactOnEmc(phot1,mod,z,x))
448	inacceptance = kTRUE;
449	//printf("In REAL PHOS acceptance? %d\n",inacceptance);
450	#else
451	TLorentzVector lv1, lv2;
452	phot1->Momentum(lv1);
453	phot2->Momentum(lv2);
454	if(GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter) && GetFidutialCut()->IsInFidutialCut(lv2,fCalorimeter))
455	inacceptance = kTRUE ;
456	//printf("In %s fidutial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
457	#endif
458	if(inacceptance){
459	fhPrimAccPt->Fill(pi0Pt) ;
460	fhPrimAccPhi->Fill(phi) ;
461	fhPrimAccY->Fill(pi0Y) ;
462	}//Accepted
463	}// 2 photons
464	}//Check daughters exist
465	}// Primary pi0
466	}//loop on primaries
467	}//stack exists and data is MC
468
469	}
470
471
472
473
474
475
476
477
478
479
480