[u/mrichter/AliRoot.git] / PWG4 / AliAnaGammaDirect.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
/* $Id$ */

/* History of cvs commits:
 *
 * $Log$
 * Revision 1.4.4.4  2007/07/26 10:32:09  schutz
 * new analysis classes in the the new analysis framework
 *
 *
 */

//_________________________________________________________________________
// Class for the prompt gamma analysis, isolation cut
//
//  Class created from old AliPHOSGammaJet 
//  (see AliRoot versions previous Release 4-09)
//
//*-- Author: Gustavo Conesa (LNF-INFN) 
//////////////////////////////////////////////////////////////////////////////
  
  
// --- ROOT system --- 
#include <TParticle.h>
#include <TH2.h>
#include <TList.h>
#include "AliAnaGammaDirect.h" 
#include "Riostream.h"
#include "AliLog.h"
  
ClassImp(AliAnaGammaDirect)
  
//____________________________________________________________________________
  AliAnaGammaDirect::AliAnaGammaDirect() : 
    TObject(),
    fMinGammaPt(0.),
    fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0), 
    fICMethod(0),fhNGamma(0),fhPhiGamma(0),fhEtaGamma(0),  
    //kSeveralIC
    fNCones(0),fNPtThres(0), fConeSizes(),  fPtThresholds(), 
    fhPtThresIsolated(), fhPtSumIsolated()
{
  //default ctor
  
  //Initialize parameters
  InitParameters();

}

//____________________________________________________________________________
AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) : 
  TObject(g),
  fMinGammaPt(g.fMinGammaPt), 
  fConeSize(g.fConeSize),
  fPtThreshold(g.fPtThreshold),
  fPtSumThreshold(g.fPtSumThreshold), 
  fICMethod(g.fICMethod),
  fhNGamma(g.fhNGamma),fhPhiGamma(g.fhPhiGamma),fhEtaGamma(g.fhEtaGamma),  
  //kSeveralIC
  fNCones(g.fNCones),fNPtThres(g.fNPtThres), fConeSizes(),fPtThresholds(), 
  fhPtThresIsolated(), fhPtSumIsolated()
{
  // cpy ctor
  
  //kSeveralIC
  for(Int_t i = 0; i < fNCones ; i++){
    fConeSizes[i] =  g.fConeSizes[i];
    fhPtSumIsolated[i] = g.fhPtSumIsolated[i]; 
    for(Int_t j = 0; j < fNPtThres ; j++)
      fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j]; 
  }
  
  for(Int_t i = 0; i < fNPtThres ; i++)
    fPtThresholds[i]=   g.fPtThresholds[i];
}

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

  fMinGammaPt = source.fMinGammaPt ;   
  fConeSize = source.fConeSize ;
  fPtThreshold = source.fPtThreshold ;
  fPtSumThreshold = source.fPtSumThreshold ; 
  fICMethod = source.fICMethod ;
  fhNGamma = source.fhNGamma ; 
  fhPhiGamma = source.fhPhiGamma ;
  fhEtaGamma = source.fhEtaGamma ;
  
  //kSeveralIC
  fNCones = source.fNCones ;
  fNPtThres = source.fNPtThres ; 
   
  for(Int_t i = 0; i < fNCones ; i++){
    fConeSizes[i] =  source.fConeSizes[i];
    fhPtSumIsolated[i] = source.fhPtSumIsolated[i] ;
    for(Int_t j = 0; j < fNPtThres ; j++)
      fhPtThresIsolated[i][j] = source.fhPtThresIsolated[i][j] ;
  }
  
  for(Int_t i = 0; i < fNPtThres ; i++)
    fPtThresholds[i]=   source.fPtThresholds[i];
  
  return *this;
  
}

//____________________________________________________________________________
AliAnaGammaDirect::~AliAnaGammaDirect() 
{
  // Remove all pointers
  
  delete fhNGamma    ;  
  delete fhPhiGamma  ; 
  delete fhEtaGamma   ;  
  
  //kSeveralIC
  delete [] fhPtThresIsolated ;
  delete [] fhPtSumIsolated ;
  
}

//________________________________________________________________________
TList *  AliAnaGammaDirect::GetCreateOutputObjects()
{  

  // Create histograms to be saved in output file and 
  // store them in outputContainer
  TList * outputContainer = new TList() ; 
  outputContainer->SetName("DirectGammaHistos") ; 
  
  //Histograms of highest gamma identified in Event
  fhNGamma  = new TH1F("NGamma","Number of #gamma over calorimeter",240,0,120); 
  fhNGamma->SetYTitle("N");
  fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
  outputContainer->Add(fhNGamma) ; 
  
  fhPhiGamma  = new TH2F
    ("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7); 
  fhPhiGamma->SetYTitle("#phi");
  fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
  outputContainer->Add(fhPhiGamma) ; 
  
  fhEtaGamma  = new TH2F
    ("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8); 
  fhEtaGamma->SetYTitle("#eta");
  fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
  outputContainer->Add(fhEtaGamma) ;

  if(fICMethod == kSeveralIC){
    char name[128];
    char title[128];
    for(Int_t icone = 0; icone<fNCones; icone++){
      sprintf(name,"PtSumIsolated_Cone_%d",icone);
      sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
      fhPtSumIsolated[icone]  = new TH2F(name, title,240,0,120,120,0,10);
      fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
      fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
      outputContainer->Add(fhPtSumIsolated[icone]) ; 
      
      for(Int_t ipt = 0; ipt<fNPtThres;ipt++){ 
	sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
	sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
	fhPtThresIsolated[icone][ipt]  = new TH1F(name, title,240,0,120);
	fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
	outputContainer->Add(fhPtThresIsolated[icone][ipt]) ; 
      }//icone loop
    }//ipt loop
  }

  return outputContainer ;

}

//____________________________________________________________________________
void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS, TParticle *pGamma, Bool_t &found) const 
{
  //Search for the prompt photon in Calorimeter with pt > fMinGammaPt

  Double_t pt = 0;
  Int_t index = -1; 
  for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
    TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;

    if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt)){
      index = ipr ;
      pt = particle->Pt();
      pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
      found  = kTRUE;
    }
  }

  //Do Isolation?
  if( ( fICMethod == kPtIC  ||  fICMethod == kSumPtIC )  && found)
    {
      Float_t coneptsum = 0 ;
      Bool_t  icPtThres   = kFALSE;
      Bool_t  icPtSum     = kFALSE;
      MakeIsolationCut(plCTS,pl, pGamma, index, 
		       icPtThres, icPtSum,coneptsum);
      if(fICMethod == kPtIC) //Pt thres method
	found = icPtThres ;
      if(fICMethod == kSumPtIC) //Pt cone sum method
	found = icPtSum ;
    }
  
  if(found){
    AliDebug(1,Form("Cluster with p_{T} larger than %f found in calorimeter ", fMinGammaPt)) ;
    AliDebug(1,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
    //Fill prompt gamma histograms 
    fhNGamma->Fill(pGamma->Pt());
    fhPhiGamma->Fill( pGamma->Pt(),pGamma->Phi());
    fhEtaGamma->Fill(pGamma->Pt(),pGamma->Eta());
  }
  else
    AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
}

  //____________________________________________________________________________
void AliAnaGammaDirect::InitParameters()
{
 
  //Initialize the parameters of the analysis.
  fMinGammaPt  = 5. ;

  //Fill particle lists when PID is ok
  fConeSize             = 0.2 ; 
  fPtThreshold         = 2.0; 
  fPtSumThreshold  = 1.; 

  fICMethod = kNoIC; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method

 //-----------kSeveralIC-----------------
  fNCones           = 4 ; 
  fNPtThres         = 4 ; 
  fConeSizes[0] = 0.1; fConeSizes[1] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4;
  fPtThresholds[0]=1.; fPtThresholds[1]=2.; fPtThresholds[2]=3.; fPtThresholds[3]=4.;

}

//__________________________________________________________________
void  AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS, 
					   TClonesArray * plNe, 
					   TParticle * pCandidate, 
					   Int_t index, 
					   Bool_t  &icmpt,  Bool_t  &icms, 
					   Float_t &coneptsum) const 
{  
  //Search in cone around a candidate particle if it is isolated 
  Float_t phiC  = pCandidate->Phi() ;
  Float_t etaC = pCandidate->Eta() ;
  Float_t pt     = -100. ;
  Float_t eta   = -100.  ;
  Float_t phi    = -100.  ;
  Float_t rad   = -100 ;
  Int_t    n        = 0 ;
  TParticle * particle  = new TParticle;

  coneptsum = 0; 
  icmpt = kFALSE;
  icms  = kFALSE;

  //Check charged particles in cone.
  for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
    particle = dynamic_cast<TParticle *>(plCTS->At(ipr)) ;
    pt    = particle->Pt();
    eta  = particle->Eta();
    phi  = particle->Phi() ;
    
    //Check if there is any particle inside cone with pt larger than  fPtThreshold
    rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
		      TMath::Power((phi-phiC),2));
    if(rad<fConeSize){
      AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
      coneptsum+=pt;
      if(pt > fPtThreshold ) n++;
    }
  }// charged particle loop
  
  //Check neutral particles in cone.
  for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
    if(ipr != index){//Do not count the candidate
      particle = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
      pt    = particle->Pt();
      eta  = particle->Eta();
      phi  = particle->Phi() ;
      
      //Check if there is any particle inside cone with pt larger than  fPtThreshold
      rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
			TMath::Power((phi-phiC),2));
      if(rad<fConeSize){
	AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
	coneptsum+=pt;
	if(pt > fPtThreshold ) n++;
      }
    }
  }// neutral particle loop
  
  if(n == 0) 
    icmpt =  kTRUE ;
  if(coneptsum < fPtSumThreshold)
    icms  =  kTRUE ;

}

//__________________________________________________________________
void  AliAnaGammaDirect::MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArray * plCTS) 
{
  //Isolation Cut Analysis for both methods and different pt cuts and cones

  if (fICMethod != kSeveralIC)
    AliFatal("Remember to set in config file: directGamma->SetICMethod(kSeveralIC)");
  
  for(Int_t ipr = 0; ipr < plCalo->GetEntries() ; ipr ++ ){
    TParticle * pCandidate = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
    
    if(pCandidate->Pt() > fMinGammaPt){
      
      Bool_t  icPtThres   = kFALSE;
      Bool_t  icPtSum     = kFALSE;
      
      Float_t ptC      = pCandidate->Pt() ;
   
      fhNGamma->Fill(ptC);
      fhPhiGamma->Fill( ptC,pCandidate->Phi());
      fhEtaGamma->Fill(ptC,pCandidate->Eta());
    
      for(Int_t icone = 0; icone<fNCones; icone++){
	fConeSize=fConeSizes[icone] ;
	Float_t coneptsum = 0 ;
	for(Int_t ipt = 0; ipt<fNPtThres;ipt++){ 
	  fPtThreshold=fPtThresholds[ipt] ;
	  MakeIsolationCut(plCTS,plCalo, pCandidate, ipr, icPtThres, icPtSum,coneptsum);
	  AliDebug(1,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
			  pCandidate->Pt(), coneptsum, icPtThres, icPtSum));

	  fhPtThresIsolated[icone][ipt]->Fill(ptC); 
	}//pt thresh loop
	fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
      }//cone size loop
    }//min pt candidate
  }//candidate loop
}

void AliAnaGammaDirect::Print(const Option_t * opt) const
{
  
  //Print some relevant parameters set for the analysis
  if(! opt)
    return;
  
  Info("Print", "%s %s", GetName(), GetTitle() ) ;
  
  printf("Min Gamma pT      =     %f\n",  fMinGammaPt) ;
  printf("IC method               =     %d\n", fICMethod) ; 
  printf("Cone Size               =     %f\n", fConeSize) ; 
   if(fICMethod == kPtIC) printf("pT threshold           =     %f\n", fPtThreshold) ;
   if(fICMethod == kSumPtIC) printf("pT sum threshold   =     %f\n", fPtSumThreshold) ;
   
  if(fICMethod == kSeveralIC){
    printf("N Cone Sizes               =     %d\n", fNCones) ; 
    printf("N pT thresholds           =     %d\n", fNPtThres) ;
    printf("Cone Sizes                  =    \n") ;
    for(Int_t i = 0; i < fNCones; i++)
      printf("   %f;",  fConeSizes[i]) ;
    printf("    \n") ;
    for(Int_t i = 0; i < fNPtThres; i++)
      printf("   %f;",  fPtThresholds[i]) ;
  }

  printf("    \n") ;
  
}
Commit	Line	Data
f9cea31c	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	/* History of cvs commits:
	18	*
	19	* $Log$
bdcfac30	20	* Revision 1.4.4.4 2007/07/26 10:32:09 schutz
bdcfac30	21	* new analysis classes in the the new analysis framework
2a1d8a29	22	*
f9cea31c	23	*
	24	*/
	25
	26	//_________________________________________________________________________
bdcfac30	27	// Class for the prompt gamma analysis, isolation cut
f9cea31c	28	//
	29	// Class created from old AliPHOSGammaJet
	30	// (see AliRoot versions previous Release 4-09)
	31	//
	32	//*-- Author: Gustavo Conesa (LNF-INFN)
	33	//////////////////////////////////////////////////////////////////////////////
bdcfac30	34
	35
	36	// --- ROOT system ---
f9cea31c	37	#include <TParticle.h>
f9cea31c	38	#include <TH2.h>
bdcfac30	39	#include <TList.h>
f9cea31c	40	#include "AliAnaGammaDirect.h"
f9cea31c	41	#include "Riostream.h"
f9cea31c	42	#include "AliLog.h"
bdcfac30	43
f9cea31c	44	ClassImp(AliAnaGammaDirect)
bdcfac30	45
f9cea31c	46	//____________________________________________________________________________
bdcfac30	47	AliAnaGammaDirect::AliAnaGammaDirect() :
	48	TObject(),
	49	fMinGammaPt(0.),
f9cea31c	50	fConeSize(0.),fPtThreshold(0.),fPtSumThreshold(0),
bdcfac30	51	fICMethod(0),fhNGamma(0),fhPhiGamma(0),fhEtaGamma(0),
	52	//kSeveralIC
	53	fNCones(0),fNPtThres(0), fConeSizes(), fPtThresholds(),
	54	fhPtThresIsolated(), fhPtSumIsolated()
f9cea31c	55	{
bdcfac30	56	//default ctor
bdcfac30	57
463ee300	58	//Initialize parameters
	59	InitParameters();
	60
f9cea31c	61	}
f9cea31c	62
f9cea31c	63	//____________________________________________________________________________
f9cea31c	64	AliAnaGammaDirect::AliAnaGammaDirect(const AliAnaGammaDirect & g) :
bdcfac30	65	TObject(g),
bdcfac30	66	fMinGammaPt(g.fMinGammaPt),
f9cea31c	67	fConeSize(g.fConeSize),
	68	fPtThreshold(g.fPtThreshold),
	69	fPtSumThreshold(g.fPtSumThreshold),
bdcfac30	70	fICMethod(g.fICMethod),
	71	fhNGamma(g.fhNGamma),fhPhiGamma(g.fhPhiGamma),fhEtaGamma(g.fhEtaGamma),
	72	//kSeveralIC
	73	fNCones(g.fNCones),fNPtThres(g.fNPtThres), fConeSizes(),fPtThresholds(),
	74	fhPtThresIsolated(), fhPtSumIsolated()
f9cea31c	75	{
f9cea31c	76	// cpy ctor
bdcfac30	77
	78	//kSeveralIC
	79	for(Int_t i = 0; i < fNCones ; i++){
	80	fConeSizes[i] = g.fConeSizes[i];
	81	fhPtSumIsolated[i] = g.fhPtSumIsolated[i];
	82	for(Int_t j = 0; j < fNPtThres ; j++)
	83	fhPtThresIsolated[i][j] = g.fhPtThresIsolated[i][j];
	84	}
	85
	86	for(Int_t i = 0; i < fNPtThres ; i++)
	87	fPtThresholds[i]= g.fPtThresholds[i];
f9cea31c	88	}
f9cea31c	89
463ee300	90	//_________________________________________________________________________
	91	AliAnaGammaDirect & AliAnaGammaDirect::operator = (const AliAnaGammaDirect & source)
	92	{
	93	// assignment operator
bdcfac30	94
463ee300	95	if(&source == this) return *this;
463ee300	96
463ee300	97	fMinGammaPt = source.fMinGammaPt ;
463ee300	98	fConeSize = source.fConeSize ;
	99	fPtThreshold = source.fPtThreshold ;
	100	fPtSumThreshold = source.fPtSumThreshold ;
bdcfac30	101	fICMethod = source.fICMethod ;
463ee300	102	fhNGamma = source.fhNGamma ;
	103	fhPhiGamma = source.fhPhiGamma ;
	104	fhEtaGamma = source.fhEtaGamma ;
bdcfac30	105
	106	//kSeveralIC
	107	fNCones = source.fNCones ;
	108	fNPtThres = source.fNPtThres ;
	109
	110	for(Int_t i = 0; i < fNCones ; i++){
	111	fConeSizes[i] = source.fConeSizes[i];
	112	fhPtSumIsolated[i] = source.fhPtSumIsolated[i] ;
	113	for(Int_t j = 0; j < fNPtThres ; j++)
	114	fhPtThresIsolated[i][j] = source.fhPtThresIsolated[i][j] ;
	115	}
	116
	117	for(Int_t i = 0; i < fNPtThres ; i++)
	118	fPtThresholds[i]= source.fPtThresholds[i];
	119
463ee300	120	return *this;
bdcfac30	121
463ee300	122	}
463ee300	123
f9cea31c	124	//____________________________________________________________________________
	125	AliAnaGammaDirect::~AliAnaGammaDirect()
	126	{
	127	// Remove all pointers
bdcfac30	128
f9cea31c	129	delete fhNGamma ;
	130	delete fhPhiGamma ;
	131	delete fhEtaGamma ;
2a1d8a29	132
bdcfac30	133	//kSeveralIC
	134	delete [] fhPtThresIsolated ;
	135	delete [] fhPtSumIsolated ;
2a1d8a29	136
2a1d8a29	137	}
	138
	139	//________________________________________________________________________
bdcfac30	140	TList * AliAnaGammaDirect::GetCreateOutputObjects()
2a1d8a29	141	{
2a1d8a29	142
463ee300	143	// Create histograms to be saved in output file and
bdcfac30	144	// store them in outputContainer
	145	TList * outputContainer = new TList() ;
	146	outputContainer->SetName("DirectGammaHistos") ;
	147
2a1d8a29	148	//Histograms of highest gamma identified in Event
bdcfac30	149	fhNGamma = new TH1F("NGamma","Number of #gamma over calorimeter",240,0,120);
2a1d8a29	150	fhNGamma->SetYTitle("N");
2a1d8a29	151	fhNGamma->SetXTitle("p_{T #gamma}(GeV/c)");
bdcfac30	152	outputContainer->Add(fhNGamma) ;
2a1d8a29	153
	154	fhPhiGamma = new TH2F
	155	("PhiGamma","#phi_{#gamma}",200,0,120,200,0,7);
	156	fhPhiGamma->SetYTitle("#phi");
	157	fhPhiGamma->SetXTitle("p_{T #gamma} (GeV/c)");
bdcfac30	158	outputContainer->Add(fhPhiGamma) ;
2a1d8a29	159
	160	fhEtaGamma = new TH2F
	161	("EtaGamma","#phi_{#gamma}",200,0,120,200,-0.8,0.8);
	162	fhEtaGamma->SetYTitle("#eta");
	163	fhEtaGamma->SetXTitle("p_{T #gamma} (GeV/c)");
bdcfac30	164	outputContainer->Add(fhEtaGamma) ;
	165
	166	if(fICMethod == kSeveralIC){
	167	char name[128];
	168	char title[128];
	169	for(Int_t icone = 0; icone<fNCones; icone++){
	170	sprintf(name,"PtSumIsolated_Cone_%d",icone);
	171	sprintf(title,"Candidate cone sum p_{T} for cone size %d vs candidate p_{T}",icone);
	172	fhPtSumIsolated[icone] = new TH2F(name, title,240,0,120,120,0,10);
	173	fhPtSumIsolated[icone]->SetYTitle("#Sigma p_{T} (GeV/c)");
	174	fhPtSumIsolated[icone]->SetXTitle("p_{T} (GeV/c)");
	175	outputContainer->Add(fhPtSumIsolated[icone]) ;
	176
	177	for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
	178	sprintf(name,"PtThresIsol_Cone_%d_Pt%d",icone,ipt);
	179	sprintf(title,"Isolated candidate p_{T} distribution for cone size %d and p_{T}^{th} %d",icone,ipt);
	180	fhPtThresIsolated[icone][ipt] = new TH1F(name, title,240,0,120);
	181	fhPtThresIsolated[icone][ipt]->SetXTitle("p_{T} (GeV/c)");
	182	outputContainer->Add(fhPtThresIsolated[icone][ipt]) ;
	183	}//icone loop
	184	}//ipt loop
f9cea31c	185	}
f9cea31c	186
bdcfac30	187	return outputContainer ;
463ee300	188
f9cea31c	189	}
f9cea31c	190
f9cea31c	191	//____________________________________________________________________________
bdcfac30	192	void AliAnaGammaDirect::GetPromptGamma(TClonesArray * pl, TClonesArray * plCTS, TParticle *pGamma, Bool_t &found) const
f9cea31c	193	{
	194	//Search for the prompt photon in Calorimeter with pt > fMinGammaPt
	195
	196	Double_t pt = 0;
	197	Int_t index = -1;
	198	for(Int_t ipr = 0;ipr < pl->GetEntries() ; ipr ++ ){
	199	TParticle * particle = dynamic_cast<TParticle *>(pl->At(ipr)) ;
	200
	201	if((particle->Pt() > fMinGammaPt) && (particle->Pt() > pt)){
	202	index = ipr ;
	203	pt = particle->Pt();
	204	pGamma->SetMomentum(particle->Px(),particle->Py(),particle->Pz(),particle->Energy());
bdcfac30	205	found = kTRUE;
f9cea31c	206	}
	207	}
	208
	209	//Do Isolation?
bdcfac30	210	if( ( fICMethod == kPtIC \|\| fICMethod == kSumPtIC ) && found)
f9cea31c	211	{
	212	Float_t coneptsum = 0 ;
	213	Bool_t icPtThres = kFALSE;
	214	Bool_t icPtSum = kFALSE;
	215	MakeIsolationCut(plCTS,pl, pGamma, index,
	216	icPtThres, icPtSum,coneptsum);
bdcfac30	217	if(fICMethod == kPtIC) //Pt thres method
	218	found = icPtThres ;
	219	if(fICMethod == kSumPtIC) //Pt cone sum method
	220	found = icPtSum ;
f9cea31c	221	}
f9cea31c	222
bdcfac30	223	if(found){
	224	AliDebug(1,Form("Cluster with p_{T} larger than %f found in calorimeter ", fMinGammaPt)) ;
	225	AliDebug(1,Form("Gamma: pt %f, phi %f, eta %f", pGamma->Pt(),pGamma->Phi(),pGamma->Eta())) ;
	226	//Fill prompt gamma histograms
	227	fhNGamma->Fill(pGamma->Pt());
f9cea31c	228	fhPhiGamma->Fill( pGamma->Pt(),pGamma->Phi());
	229	fhEtaGamma->Fill(pGamma->Pt(),pGamma->Eta());
	230	}
	231	else
	232	AliDebug(1,Form("NO Cluster with pT larger than %f found in calorimeter ", fMinGammaPt)) ;
	233	}
	234
	235	//____________________________________________________________________________
2a1d8a29	236	void AliAnaGammaDirect::InitParameters()
f9cea31c	237	{
f9cea31c	238
2a1d8a29	239	//Initialize the parameters of the analysis.
2a1d8a29	240	fMinGammaPt = 5. ;
f9cea31c	241
f9cea31c	242	//Fill particle lists when PID is ok
f9cea31c	243	fConeSize = 0.2 ;
	244	fPtThreshold = 2.0;
	245	fPtSumThreshold = 1.;
	246
bdcfac30	247	fICMethod = kNoIC; // 0 don't isolate, 1 pt thresh method, 2 cone pt sum method
	248
	249	//-----------kSeveralIC-----------------
	250	fNCones = 4 ;
	251	fNPtThres = 4 ;
	252	fConeSizes[0] = 0.1; fConeSizes[1] = 0.2; fConeSizes[2] = 0.3; fConeSizes[3] = 0.4;
	253	fPtThresholds[0]=1.; fPtThresholds[1]=2.; fPtThresholds[2]=3.; fPtThresholds[3]=4.;
	254
f9cea31c	255	}
	256
	257	//__________________________________________________________________
	258	void AliAnaGammaDirect::MakeIsolationCut(TClonesArray * plCTS,
	259	TClonesArray * plNe,
	260	TParticle * pCandidate,
	261	Int_t index,
	262	Bool_t &icmpt, Bool_t &icms,
	263	Float_t &coneptsum) const
	264	{
	265	//Search in cone around a candidate particle if it is isolated
	266	Float_t phiC = pCandidate->Phi() ;
	267	Float_t etaC = pCandidate->Eta() ;
	268	Float_t pt = -100. ;
	269	Float_t eta = -100. ;
	270	Float_t phi = -100. ;
	271	Float_t rad = -100 ;
	272	Int_t n = 0 ;
bdcfac30	273	TParticle * particle = new TParticle;
f9cea31c	274
	275	coneptsum = 0;
	276	icmpt = kFALSE;
	277	icms = kFALSE;
	278
	279	//Check charged particles in cone.
	280	for(Int_t ipr = 0;ipr < plCTS->GetEntries() ; ipr ++ ){
	281	particle = dynamic_cast<TParticle *>(plCTS->At(ipr)) ;
	282	pt = particle->Pt();
	283	eta = particle->Eta();
	284	phi = particle->Phi() ;
	285
	286	//Check if there is any particle inside cone with pt larger than fPtThreshold
	287	rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
	288	TMath::Power((phi-phiC),2));
	289	if(rad<fConeSize){
	290	AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
	291	coneptsum+=pt;
	292	if(pt > fPtThreshold ) n++;
	293	}
	294	}// charged particle loop
	295
	296	//Check neutral particles in cone.
	297	for(Int_t ipr = 0;ipr < plNe->GetEntries() ; ipr ++ ){
	298	if(ipr != index){//Do not count the candidate
	299	particle = dynamic_cast<TParticle *>(plNe->At(ipr)) ;
	300	pt = particle->Pt();
	301	eta = particle->Eta();
	302	phi = particle->Phi() ;
	303
	304	//Check if there is any particle inside cone with pt larger than fPtThreshold
	305	rad = TMath::Sqrt(TMath::Power((eta-etaC),2) +
	306	TMath::Power((phi-phiC),2));
	307	if(rad<fConeSize){
	308	AliDebug(3,Form("charged in cone pt %f, phi %f, eta %f, R %f ",pt,phi,eta,rad));
	309	coneptsum+=pt;
	310	if(pt > fPtThreshold ) n++;
	311	}
	312	}
	313	}// neutral particle loop
	314
	315	if(n == 0)
	316	icmpt = kTRUE ;
	317	if(coneptsum < fPtSumThreshold)
	318	icms = kTRUE ;
	319
	320	}
	321
bdcfac30	322	//__________________________________________________________________
bdcfac30	323	void AliAnaGammaDirect::MakeSeveralICAnalysis(TClonesArray * plCalo, TClonesArray * plCTS)
f9cea31c	324	{
bdcfac30	325	//Isolation Cut Analysis for both methods and different pt cuts and cones
f9cea31c	326
bdcfac30	327	if (fICMethod != kSeveralIC)
	328	AliFatal("Remember to set in config file: directGamma->SetICMethod(kSeveralIC)");
	329
	330	for(Int_t ipr = 0; ipr < plCalo->GetEntries() ; ipr ++ ){
	331	TParticle * pCandidate = dynamic_cast<TParticle *>(plCalo->At(ipr)) ;
	332
	333	if(pCandidate->Pt() > fMinGammaPt){
	334
	335	Bool_t icPtThres = kFALSE;
	336	Bool_t icPtSum = kFALSE;
	337
	338	Float_t ptC = pCandidate->Pt() ;
	339
	340	fhNGamma->Fill(ptC);
	341	fhPhiGamma->Fill( ptC,pCandidate->Phi());
	342	fhEtaGamma->Fill(ptC,pCandidate->Eta());
	343
	344	for(Int_t icone = 0; icone<fNCones; icone++){
	345	fConeSize=fConeSizes[icone] ;
	346	Float_t coneptsum = 0 ;
	347	for(Int_t ipt = 0; ipt<fNPtThres;ipt++){
	348	fPtThreshold=fPtThresholds[ipt] ;
	349	MakeIsolationCut(plCTS,plCalo, pCandidate, ipr, icPtThres, icPtSum,coneptsum);
	350	AliDebug(1,Form("Candidate pt %f, pt in cone %f, Isolated? ICPt %d, ICSum %d",
	351	pCandidate->Pt(), coneptsum, icPtThres, icPtSum));
	352
	353	fhPtThresIsolated[icone][ipt]->Fill(ptC);
	354	}//pt thresh loop
	355	fhPtSumIsolated[icone]->Fill(ptC,coneptsum) ;
	356	}//cone size loop
	357	}//min pt candidate
	358	}//candidate loop
	359	}
	360
	361	void AliAnaGammaDirect::Print(const Option_t * opt) const
	362	{
	363
f9cea31c	364	//Print some relevant parameters set for the analysis
	365	if(! opt)
	366	return;
bdcfac30	367
f9cea31c	368	Info("Print", "%s %s", GetName(), GetTitle() ) ;
bdcfac30	369
	370	printf("Min Gamma pT = %f\n", fMinGammaPt) ;
	371	printf("IC method = %d\n", fICMethod) ;
f9cea31c	372	printf("Cone Size = %f\n", fConeSize) ;
bdcfac30	373	if(fICMethod == kPtIC) printf("pT threshold = %f\n", fPtThreshold) ;
	374	if(fICMethod == kSumPtIC) printf("pT sum threshold = %f\n", fPtSumThreshold) ;
	375
	376	if(fICMethod == kSeveralIC){
	377	printf("N Cone Sizes = %d\n", fNCones) ;
	378	printf("N pT thresholds = %d\n", fNPtThres) ;
	379	printf("Cone Sizes = \n") ;
	380	for(Int_t i = 0; i < fNCones; i++)
	381	printf(" %f;", fConeSizes[i]) ;
	382	printf(" \n") ;
	383	for(Int_t i = 0; i < fNPtThres; i++)
	384	printf(" %f;", fPtThresholds[i]) ;
	385	}
f9cea31c	386
bdcfac30	387	printf(" \n") ;
	388
	389	}