[u/mrichter/AliRoot.git] / PWG4 / PartCorrDep / AliAnaShowerParameter.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 cloned from AliAnaPhoton, main aim is shower shape studies
// 
// 
//
//-- Author: Jocelyn Mlynarz (WSU) and Gustavo Conesa (LPSC-CNRS)
//////////////////////////////////////////////////////////////////////////////


// --- ROOT system --- 
#include <TH3F.h>
#include <TH2F.h>
#include <TClonesArray.h>
//#include <TObjString.h>
#include <Riostream.h>
#include "TParticle.h"
//#include <fstream>

// --- Analysis system --- 
#include "AliAnaShowerParameter.h" 
#include "AliCaloTrackReader.h"
#include "AliStack.h"
#include "AliCaloPID.h"
#include "AliMCAnalysisUtils.h"
#include "AliFiducialCut.h"
#include "AliAODCaloCluster.h"
#include "AliAODMCParticle.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisTaskParticleCorrelation.h"
#include "AliAODEvent.h"


ClassImp(AliAnaShowerParameter)

//____________________________________________________________________________
AliAnaShowerParameter::AliAnaShowerParameter() : 
AliAnaPartCorrBaseClass(), fCalorimeter(""), fNCellsCutMin(0),
fNCellsCutMax(0), fLambdaCut(0), fTimeCutMin(-1), fTimeCutMax(9999999),
fhNClusters(0), fhNCellCluster(0), fhEtaPhiPtCluster(0), 
fhLambdaPtCluster(0), 

//MC

fhLambdaPtPhoton(0), fhLambdaPtPi0(0), fhLambdaPtPion(0), fhPtTruthPi0(0)

{
  //default ctor
  
  //Initialize parameters
  InitParameters();
  
}

//____________________________________________________________________________
AliAnaShowerParameter::~AliAnaShowerParameter() 
{
  //dtor
  
}

//________________________________________________________________________
TObjString *  AliAnaShowerParameter::GetAnalysisCuts()
{  	
  //Save parameters used for analysis
  TString parList ; //this will be list of parameters used for this analysis.
  const Int_t buffersize = 255;
  char onePar[buffersize] ;
  
  snprintf(onePar,buffersize,"--- AliAnaShowerParameter ---\n") ;
  parList+=onePar ;	
  snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
  parList+=onePar ;
  snprintf(onePar,buffersize,"fNCellsCutMin: (Cut in the minimum number of cells) %e\n",fNCellsCutMin) ;
  parList+=onePar ;  
  snprintf(onePar,buffersize,"fNCellsCutMax: (Cut in the maximum number of cells) %e\n",fNCellsCutMax) ;
  parList+=onePar ;
  snprintf(onePar,buffersize,"fLambdaCut: (Cut in the minimum lambda) %e\n",fLambdaCut) ;
  parList+=onePar ;
  
  //Get parameters set in base class.
  parList += GetBaseParametersList() ;
  
  //Get parameters set in PID class.
  parList += GetCaloPID()->GetPIDParametersList() ;
  
  //Get parameters set in FiducialCut class (not available yet)
  //parlist += GetFidCut()->GetFidCutParametersList() 
  
  return new TObjString(parList) ;
}


//________________________________________________________________________
TList *  AliAnaShowerParameter::GetCreateOutputObjects()
{  
  // Create histograms to be saved in output file and 
  // store them in outputContainer
  TList * outputContainer = new TList() ; 
  outputContainer->SetName("PhotonHistos") ; 
	
  Int_t nptbins  = GetHistoPtBins();
  Int_t nphibins = GetHistoPhiBins();
  Int_t netabins = GetHistoEtaBins();
  Float_t ptmax  = GetHistoPtMax();
  Float_t phimax = GetHistoPhiMax();
  Float_t etamax = GetHistoEtaMax();
  Float_t ptmin  = GetHistoPtMin();
  Float_t phimin = GetHistoPhiMin();
  Float_t etamin = GetHistoEtaMin();	
  
  //General non-MC Cluster histograms
  fhNClusters = new TH1F ("hNClusters","NClusters",21,-0.5,20.5); 
  fhNClusters->SetXTitle("N_{Clusters}");
  outputContainer->Add(fhNClusters) ;
  
  fhNCellCluster = new TH2F ("hNCellCluster","Number of cell per cluster",nptbins,ptmin,ptmax,21,-0.5,20.5); 
  fhNCellCluster->SetYTitle("N_{Cells}");  
  fhNCellCluster->SetXTitle("p_{T}");
  outputContainer->Add(fhNCellCluster) ;
  
  fhEtaPhiPtCluster = new TH3F
  ("hEtaPhiPtCluster","#phi_{Cluster} and #eta_{Cluster}",nptbins,ptmin,ptmax,nphibins,phimin,phimax,netabins,etamin,etamax); 
  fhEtaPhiPtCluster->SetZTitle("#eta");  
  fhEtaPhiPtCluster->SetYTitle("#phi");
  fhEtaPhiPtCluster->SetXTitle("pT_{Cluster} (GeV/c)");
  outputContainer->Add(fhEtaPhiPtCluster) ; 
  
  fhLambdaPtCluster  = new TH2F
  ("hLambdaCluster","#lambda_{Cluster}",nptbins,ptmin,ptmax,300,0,3); 
  fhLambdaPtCluster->SetYTitle("#lambda_{0}^{2}");
  fhLambdaPtCluster->SetXTitle("p_{T Cluster} (GeV/c)");
  outputContainer->Add(fhLambdaPtCluster) ;
  
  if(IsDataMC()){
    
    fhLambdaPtPhoton  = new TH2F
    ("hLambdaPtPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2); 
    fhLambdaPtPhoton->SetYTitle("#lambda_{0}^{2}");
    fhLambdaPtPhoton->SetXTitle("pT_{#gamma, Reco} (GeV)");
    outputContainer->Add(fhLambdaPtPhoton) ;
    
    fhLambdaPtPi0  = new TH2F
    ("hLambdaPtPi0","#lambda_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,2); 
    fhLambdaPtPi0->SetYTitle("#lambda_{0}^{2}");
    fhLambdaPtPi0->SetXTitle("pT_{#pi^{0}, Reco} (GeV)");
    outputContainer->Add(fhLambdaPtPi0) ;
    
    fhLambdaPtPion  = new TH2F
    ("hLambdaPtPion","#lambda_{#pi^{+}}",nptbins,ptmin,ptmax,200,0,2); 
    fhLambdaPtPion->SetYTitle("#lambda_{0}^{2}");
    fhLambdaPtPion->SetXTitle("pT_{#pi^{+}, Reco} (GeV)");
    outputContainer->Add(fhLambdaPtPion) ;
   
    fhPtTruthPi0  = new TH1D("hPtTruthPi0","#pi^{0} MC truth pT",nptbins,ptmin,ptmax) ;
    fhPtTruthPi0->SetXTitle("pT_{#pi^{0}, Reco} (GeV)");
    outputContainer->Add(fhPtTruthPi0) ;
    
  }//Histos with MC
  
  return outputContainer ;
  
}

//____________________________________________________________________________
void AliAnaShowerParameter::Init()
{
  
  //Init
  //Do some checks
  if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
    printf("AliAnaShowerParameter::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
    abort();
  }
  else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
    printf("AliAnaShowerParameter::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
    abort();
  }
  
}

//____________________________________________________________________________
void AliAnaShowerParameter::InitParameters()
{
  
  //Initialize the parameters of the analysis.
  AddToHistogramsName("AnaLambda_");
  
  fCalorimeter = "EMCAL" ;
	
  fTimeCutMin  = -1;
  fTimeCutMax  = 9999999;
  fNCellsCutMin = 0 ;
  fNCellsCutMax = 0 ;
  fLambdaCut = 0.01 ;
	
}

//__________________________________________________________________
void  AliAnaShowerParameter::MakeAnalysisFillHistograms() 
{
  
  //Do analysis and fill histograms
  if ( GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Starting analysis.\n");
  
  // Access MC information in stack if requested, check that it exists.	
  AliStack * stack = 0x0;
  TClonesArray * mcparticles0 = 0x0;
  Int_t NClusters = 0 ;
  TLorentzVector momCluster ;
  
  //Check if the stack is available when analysing MC data.
  if(IsDataMC()){
    
    if(GetReader()->ReadStack()){
      stack =  GetMCStack() ;
      if(!stack) {
        printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
				abort();
      }
      
    }
    else if(GetReader()->ReadAODMCParticles()){
      
      //Get the list of MC particles
      mcparticles0 = GetReader()->GetAODMCParticles(0);
      if(!mcparticles0 && GetDebug() > 0) 	{
        printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() -  Standard MCParticles not available !\n"); 
      }
    }// is data and MC
  }	
  //Loop on stored AOD photons
  
  TClonesArray *  clustArray = 0x0 ; 
  clustArray = GetAODCaloClusters() ;
  NClusters = clustArray->GetEntriesFast() ;    
  fhNClusters->Fill(NClusters) ;
  
  if(GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - aod branch entries %d\n", NClusters);
  
  for(Int_t iCluster = 0 ; iCluster < NClusters ; iCluster++){    
    Int_t input = 0 ;
    
    AliAODCaloCluster * caloCluster =  (AliAODCaloCluster*) (clustArray->At(iCluster)) ;
    caloCluster->GetMomentum(momCluster,GetVertex(0)) ;
    AliAODPWG4Particle * aodCluster = new AliAODPWG4Particle(momCluster) ;
    Int_t labelCluster = caloCluster->GetLabel() ;
    aodCluster->SetLabel(labelCluster) ;;
    aodCluster->SetInputFileIndex(input) ;
    aodCluster->SetCaloLabel(caloCluster->GetID(),-1) ;
    aodCluster->SetDetector(fCalorimeter) ;
    aodCluster->SetTag(GetMCAnalysisUtils()->CheckOrigin(caloCluster->GetLabels(),caloCluster->GetNLabels(),GetReader(), aodCluster->GetInputFileIndex())) ;
    
    if (GetDebug() > 3) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Cluster particle label = %d\n",labelCluster) ;
    
    //Fill Cluster histograms 
    Float_t ptcluster  = aodCluster->Pt() ;
    Float_t phicluster = aodCluster->Phi() ;
    Float_t etacluster = aodCluster->Eta() ;
    Float_t lambdaMainCluster = caloCluster->GetM02() ;
    Int_t iNumCell = caloCluster->GetNCells() ;	
    
    if (iNumCell>=fNCellsCutMin&&iNumCell<=fNCellsCutMax&&lambdaMainCluster>=fLambdaCut){
      
      //Fill the basic non-MC information about the cluster.
      fhNCellCluster->Fill(ptcluster,iNumCell) ;
      fhEtaPhiPtCluster->Fill(ptcluster,phicluster,etacluster) ;
      fhLambdaPtCluster->Fill(ptcluster,lambdaMainCluster) ;
      
      //Play with the MC data if available
      if(IsDataMC()){
        
        //Get the tag from AliMCAnalysisUtils for PID
        Int_t tag = aodCluster->GetTag();
        if (GetDebug() > 3) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Cluster particle tag= %d\n",tag) ;
        if ( aodCluster->GetLabel() < 0){
          if(GetDebug() > 0) 
            printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Label is -1; problem in the MC ESD? ");
          continue;
        }
        
        //Check if the tag matches one of the different particle types and fill the corresponding histograms
        if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)&&!(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))) {
          fhLambdaPtPhoton->Fill(ptcluster,lambdaMainCluster) ;
        }//kMCPhoton
        
        if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)) {     
          fhLambdaPtPi0->Fill(ptcluster,lambdaMainCluster) ;        
        }
        if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion)) {
          fhLambdaPtPion->Fill(ptcluster,lambdaMainCluster) ;
        }
        
        // Access MC information in stack if requested, check that it exists.
        Int_t label = aodCluster->GetLabel();
        if(label < 0) {
          printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
          continue;
        }
      }
    }
  }
  
  Float_t fPtGener = 0 ;
  if(IsDataMC() && stack){
    TParticle* pGener = stack->Particle(0) ;
    fPtGener = pGener->Pt() ;
  }
  
  
}


//__________________________________________________________________
void AliAnaShowerParameter::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 number of cells in cluster is        > %f \n", fNCellsCutMin);
  printf("Max number of cells in cluster is        > %f \n", fNCellsCutMax);
  printf("Min lambda of cluster is        > %f \n", fLambdaCut);
  printf("    \n") ;
  
}
Commit	Line	Data
2ef5608f	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 cloned from AliAnaPhoton, main aim is shower shape studies
	20	//
	21	//
	22	//
	23	//-- Author: Jocelyn Mlynarz (WSU) and Gustavo Conesa (LPSC-CNRS)
	24	//////////////////////////////////////////////////////////////////////////////
	25
	26
	27	// --- ROOT system ---
	28	#include <TH3F.h>
	29	#include <TH2F.h>
	30	#include <TClonesArray.h>
	31	//#include <TObjString.h>
	32	#include <Riostream.h>
	33	#include "TParticle.h"
	34	//#include <fstream>
	35
	36	// --- Analysis system ---
	37	#include "AliAnaShowerParameter.h"
	38	#include "AliCaloTrackReader.h"
	39	#include "AliStack.h"
	40	#include "AliCaloPID.h"
	41	#include "AliMCAnalysisUtils.h"
	42	#include "AliFiducialCut.h"
	43	#include "AliAODCaloCluster.h"
	44	#include "AliAODMCParticle.h"
	45	#include "AliAnalysisManager.h"
	46	#include "AliAnalysisTaskParticleCorrelation.h"
	47	#include "AliAODEvent.h"
	48
	49
	50	ClassImp(AliAnaShowerParameter)
	51
	52	//____________________________________________________________________________
	53	AliAnaShowerParameter::AliAnaShowerParameter() :
	54	AliAnaPartCorrBaseClass(), fCalorimeter(""), fNCellsCutMin(0),
	55	fNCellsCutMax(0), fLambdaCut(0), fTimeCutMin(-1), fTimeCutMax(9999999),
	56	fhNClusters(0), fhNCellCluster(0), fhEtaPhiPtCluster(0),
	57	fhLambdaPtCluster(0),
	58
	59	//MC
	60
	61	fhLambdaPtPhoton(0), fhLambdaPtPi0(0), fhLambdaPtPion(0), fhPtTruthPi0(0)
	62
	63	{
	64	//default ctor
65
66	//Initialize parameters
67	InitParameters();
68
69	}
70
71	//____________________________________________________________________________
72	AliAnaShowerParameter::~AliAnaShowerParameter()
73	{
74	//dtor
75
76	}
77
78	//________________________________________________________________________
79	TObjString * AliAnaShowerParameter::GetAnalysisCuts()
80	{
81	//Save parameters used for analysis
82	TString parList ; //this will be list of parameters used for this analysis.
83	const Int_t buffersize = 255;
84	char onePar[buffersize] ;
85
86	snprintf(onePar,buffersize,"--- AliAnaShowerParameter ---\n") ;
87	parList+=onePar ;
88	snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
89	parList+=onePar ;
90	snprintf(onePar,buffersize,"fNCellsCutMin: (Cut in the minimum number of cells) %e\n",fNCellsCutMin) ;
91	parList+=onePar ;
92	snprintf(onePar,buffersize,"fNCellsCutMax: (Cut in the maximum number of cells) %e\n",fNCellsCutMax) ;
93	parList+=onePar ;
94	snprintf(onePar,buffersize,"fLambdaCut: (Cut in the minimum lambda) %e\n",fLambdaCut) ;
95	parList+=onePar ;
96
97	//Get parameters set in base class.
98	parList += GetBaseParametersList() ;
99
100	//Get parameters set in PID class.
101	parList += GetCaloPID()->GetPIDParametersList() ;
102
103	//Get parameters set in FiducialCut class (not available yet)
104	//parlist += GetFidCut()->GetFidCutParametersList()
105
106	return new TObjString(parList) ;
107	}
108
109
110	//________________________________________________________________________
111	TList * AliAnaShowerParameter::GetCreateOutputObjects()
112	{
113	// Create histograms to be saved in output file and
114	// store them in outputContainer
115	TList * outputContainer = new TList() ;
116	outputContainer->SetName("PhotonHistos") ;
117
118	Int_t nptbins = GetHistoPtBins();
119	Int_t nphibins = GetHistoPhiBins();
120	Int_t netabins = GetHistoEtaBins();
121	Float_t ptmax = GetHistoPtMax();
122	Float_t phimax = GetHistoPhiMax();
123	Float_t etamax = GetHistoEtaMax();
124	Float_t ptmin = GetHistoPtMin();
125	Float_t phimin = GetHistoPhiMin();
126	Float_t etamin = GetHistoEtaMin();
127
128	//General non-MC Cluster histograms
129	fhNClusters = new TH1F ("hNClusters","NClusters",21,-0.5,20.5);
130	fhNClusters->SetXTitle("N_{Clusters}");
131	outputContainer->Add(fhNClusters) ;
132
133	fhNCellCluster = new TH2F ("hNCellCluster","Number of cell per cluster",nptbins,ptmin,ptmax,21,-0.5,20.5);
134	fhNCellCluster->SetYTitle("N_{Cells}");
135	fhNCellCluster->SetXTitle("p_{T}");
136	outputContainer->Add(fhNCellCluster) ;
137
138	fhEtaPhiPtCluster = new TH3F
139	("hEtaPhiPtCluster","#phi_{Cluster} and #eta_{Cluster}",nptbins,ptmin,ptmax,nphibins,phimin,phimax,netabins,etamin,etamax);
140	fhEtaPhiPtCluster->SetZTitle("#eta");
141	fhEtaPhiPtCluster->SetYTitle("#phi");
142	fhEtaPhiPtCluster->SetXTitle("pT_{Cluster} (GeV/c)");
143	outputContainer->Add(fhEtaPhiPtCluster) ;
144
145	fhLambdaPtCluster = new TH2F
146	("hLambdaCluster","#lambda_{Cluster}",nptbins,ptmin,ptmax,300,0,3);
147	fhLambdaPtCluster->SetYTitle("#lambda_{0}^{2}");
148	fhLambdaPtCluster->SetXTitle("p_{T Cluster} (GeV/c)");
149	outputContainer->Add(fhLambdaPtCluster) ;
150
151	if(IsDataMC()){
152
153	fhLambdaPtPhoton = new TH2F
154	("hLambdaPtPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2);
155	fhLambdaPtPhoton->SetYTitle("#lambda_{0}^{2}");
156	fhLambdaPtPhoton->SetXTitle("pT_{#gamma, Reco} (GeV)");
157	outputContainer->Add(fhLambdaPtPhoton) ;
158
159	fhLambdaPtPi0 = new TH2F
160	("hLambdaPtPi0","#lambda_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,2);
161	fhLambdaPtPi0->SetYTitle("#lambda_{0}^{2}");
162	fhLambdaPtPi0->SetXTitle("pT_{#pi^{0}, Reco} (GeV)");
163	outputContainer->Add(fhLambdaPtPi0) ;
164
165	fhLambdaPtPion = new TH2F
166	("hLambdaPtPion","#lambda_{#pi^{+}}",nptbins,ptmin,ptmax,200,0,2);
167	fhLambdaPtPion->SetYTitle("#lambda_{0}^{2}");
168	fhLambdaPtPion->SetXTitle("pT_{#pi^{+}, Reco} (GeV)");
169	outputContainer->Add(fhLambdaPtPion) ;
170
171	fhPtTruthPi0 = new TH1D("hPtTruthPi0","#pi^{0} MC truth pT",nptbins,ptmin,ptmax) ;
172	fhPtTruthPi0->SetXTitle("pT_{#pi^{0}, Reco} (GeV)");
173	outputContainer->Add(fhPtTruthPi0) ;
174
175	}//Histos with MC
176
177	return outputContainer ;
178
179	}
180
181	//____________________________________________________________________________
182	void AliAnaShowerParameter::Init()
183	{
184
185	//Init
186	//Do some checks
187	if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
188	printf("AliAnaShowerParameter::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
189	abort();
190	}
191	else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
192	printf("AliAnaShowerParameter::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
193	abort();
194	}
195
196	}
197
198	//____________________________________________________________________________
199	void AliAnaShowerParameter::InitParameters()
200	{
201
202	//Initialize the parameters of the analysis.
203	AddToHistogramsName("AnaLambda_");
204
205	fCalorimeter = "EMCAL" ;
206
207	fTimeCutMin = -1;
208	fTimeCutMax = 9999999;
209	fNCellsCutMin = 0 ;
210	fNCellsCutMax = 0 ;
211	fLambdaCut = 0.01 ;
212
213	}
214
2ef5608f	215	//__________________________________________________________________
	216	void AliAnaShowerParameter::MakeAnalysisFillHistograms()
	217	{
	218
	219	//Do analysis and fill histograms
	220	if ( GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Starting analysis.\n");
	221
	222	// Access MC information in stack if requested, check that it exists.
	223	AliStack * stack = 0x0;
2ef5608f	224	TClonesArray * mcparticles0 = 0x0;
2ef5608f	225	Int_t NClusters = 0 ;
	226	TLorentzVector momCluster ;
	227
2ef5608f	228	//Check if the stack is available when analysing MC data.
	229	if(IsDataMC()){
	230
	231	if(GetReader()->ReadStack()){
	232	stack = GetMCStack() ;
	233	if(!stack) {
	234	printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
	235	abort();
	236	}
	237
	238	}
	239	else if(GetReader()->ReadAODMCParticles()){
	240
	241	//Get the list of MC particles
	242	mcparticles0 = GetReader()->GetAODMCParticles(0);
	243	if(!mcparticles0 && GetDebug() > 0) {
	244	printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Standard MCParticles not available !\n");
	245	}
	246	}// is data and MC
	247	}
	248	//Loop on stored AOD photons
	249
	250	TClonesArray * clustArray = 0x0 ;
	251	clustArray = GetAODCaloClusters() ;
	252	NClusters = clustArray->GetEntriesFast() ;
	253	fhNClusters->Fill(NClusters) ;
	254
	255	if(GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - aod branch entries %d\n", NClusters);
	256
	257	for(Int_t iCluster = 0 ; iCluster < NClusters ; iCluster++){
	258	Int_t input = 0 ;
	259
	260	AliAODCaloCluster * caloCluster = (AliAODCaloCluster*) (clustArray->At(iCluster)) ;
	261	caloCluster->GetMomentum(momCluster,GetVertex(0)) ;
	262	AliAODPWG4Particle * aodCluster = new AliAODPWG4Particle(momCluster) ;
	263	Int_t labelCluster = caloCluster->GetLabel() ;
	264	aodCluster->SetLabel(labelCluster) ;;
	265	aodCluster->SetInputFileIndex(input) ;
	266	aodCluster->SetCaloLabel(caloCluster->GetID(),-1) ;
	267	aodCluster->SetDetector(fCalorimeter) ;
	268	aodCluster->SetTag(GetMCAnalysisUtils()->CheckOrigin(caloCluster->GetLabels(),caloCluster->GetNLabels(),GetReader(), aodCluster->GetInputFileIndex())) ;
	269
	270	if (GetDebug() > 3) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Cluster particle label = %d\n",labelCluster) ;
	271
	272	//Fill Cluster histograms
	273	Float_t ptcluster = aodCluster->Pt() ;
	274	Float_t phicluster = aodCluster->Phi() ;
	275	Float_t etacluster = aodCluster->Eta() ;
	276	Float_t lambdaMainCluster = caloCluster->GetM02() ;
	277	Int_t iNumCell = caloCluster->GetNCells() ;
	278
	279	if (iNumCell>=fNCellsCutMin&&iNumCell<=fNCellsCutMax&&lambdaMainCluster>=fLambdaCut){
	280
	281	//Fill the basic non-MC information about the cluster.
	282	fhNCellCluster->Fill(ptcluster,iNumCell) ;
	283	fhEtaPhiPtCluster->Fill(ptcluster,phicluster,etacluster) ;
	284	fhLambdaPtCluster->Fill(ptcluster,lambdaMainCluster) ;
	285
	286	//Play with the MC data if available
	287	if(IsDataMC()){
	288
2ef5608f	289	//Get the tag from AliMCAnalysisUtils for PID
	290	Int_t tag = aodCluster->GetTag();
	291	if (GetDebug() > 3) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Cluster particle tag= %d\n",tag) ;
	292	if ( aodCluster->GetLabel() < 0){
	293	if(GetDebug() > 0)
	294	printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Label is -1; problem in the MC ESD? ");
	295	continue;
	296	}
	297
	298	//Check if the tag matches one of the different particle types and fill the corresponding histograms
	299	if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)&&!(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))) {
	300	fhLambdaPtPhoton->Fill(ptcluster,lambdaMainCluster) ;
	301	}//kMCPhoton
	302
	303	if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)) {
	304	fhLambdaPtPi0->Fill(ptcluster,lambdaMainCluster) ;
	305	}
	306	if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion)) {
	307	fhLambdaPtPion->Fill(ptcluster,lambdaMainCluster) ;
	308	}
	309
	310	// Access MC information in stack if requested, check that it exists.
	311	Int_t label = aodCluster->GetLabel();
	312	if(label < 0) {
	313	printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() * bad label *: label %d \n", label);
	314	continue;
	315	}
	316	}
	317	}
	318	}
	319
	320	Float_t fPtGener = 0 ;
	321	if(IsDataMC() && stack){
	322	TParticle* pGener = stack->Particle(0) ;
	323	fPtGener = pGener->Pt() ;
	324	}
	325
	326
	327	}
	328
	329
	330	//__________________________________________________________________
	331	void AliAnaShowerParameter::Print(const Option_t * opt) const
	332	{
	333	//Print some relevant parameters set for the analysis
	334
	335	if(! opt)
	336	return;
	337
	338	printf("** Print %s %s **\n", GetName(), GetTitle() ) ;
	339	AliAnaPartCorrBaseClass::Print(" ");
	340	printf("Calorimeter = %s\n", fCalorimeter.Data()) ;
	341	printf("Min number of cells in cluster is > %f \n", fNCellsCutMin);
	342	printf("Max number of cells in cluster is > %f \n", fNCellsCutMax);
	343	printf("Min lambda of cluster is > %f \n", fLambdaCut);
	344	printf(" \n") ;
	345
	346	}