[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / 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.                  *
 **************************************************************************/

//_________________________________________________________________________
//
// 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 "AliAODEvent.h"


ClassImp(AliAnaShowerParameter)

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