From: rvernet Date: Fri, 25 Mar 2011 18:43:05 +0000 (+0000) Subject: adapted to part corr fw X-Git-Url: http://git.uio.no/git/?a=commitdiff_plain;h=ab07843c222d524809ac2a5859f4e260b2d38b97;p=u%2Fmrichter%2FAliRoot.git adapted to part corr fw --- diff --git a/PWG4/PartCorrDep/AliAnaShowerParameter.cxx b/PWG4/PartCorrDep/AliAnaShowerParameter.cxx index 486ed4fe274..98dba8f19cf 100644 --- a/PWG4/PartCorrDep/AliAnaShowerParameter.cxx +++ b/PWG4/PartCorrDep/AliAnaShowerParameter.cxx @@ -1,1207 +1,529 @@ -/************************************************************************** - * 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 -#include -#include -//#include -#include -#include "TParticle.h" -//#include - -// --- 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 "AliAnaShowerParameter.h" -#include "AliEMCALGeoUtils.h" -#include "AliAODCaloCells.h" -#include "AliAODEvent.h" - - -ClassImp(AliAnaShowerParameter) - -//____________________________________________________________________________ -AliAnaShowerParameter::AliAnaShowerParameter() : -AliAnaPartCorrBaseClass(), fCalorimeter(""), -fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0), -fCheckConversion(kFALSE),fAddConvertedPairsToAOD(kFALSE), fMassCut(0), fNCellsCut(0), -fTimeCutMin(-1), fTimeCutMax(9999999),fEMCALGeo(0),fNumClusters(0), -fhNClusters(0),fhNCellCluster(0),fhPtCluster(0),fhPhiCluster(0),fhEtaCluster(0),fhDeltaPhiClusters(0),fhDeltaEtaClusters(0), -fhLambdaCluster(0),fhDispersionCluster(0),fhELambdaCluster(0),fhELambdaCellCluster(0), -fhNCellPhoton(0),fhLambdaPhoton(0),fhDispersionPhoton(0),fhELambdaPhoton(0),fhELambdaCellPhoton(0), -fhNCellPi0(0),fhLambdaPi0(0),fhDispersionPi0(0),fhELambdaPi0(0),fhELambdaCellPi0(0), -fhNCellChargedHadron(0),fhLambdaChargedHadron(0),fhDispersionChargedHadron(0),fhELambdaChargedHadron(0),fhELambdaCellChargedHadron(0), -//MC -fhDeltaE(0), fhDeltaPt(0),fhRatioE(0), fhRatioPt(0),fh2E(0),fh2Pt(0),fhMCPdg(0), -fhEMCPhoton(0),fhPhiMCPhoton(0),fhEtaMCPhoton(0),fhLambdaMCPhoton(0),fhPhotTrueE(0),fhRatioEPhoton(0), -fhEMCPi0(0),fhPhiMCPi0(0),fhEtaMCPi0(0),fhLambdaMCPi0(0),fhPi0TrueE(0), -fhProductionDistance(0),fhProductionRadius(0),fhDecayAngle(0),fhRatioEPi0(0), -fhEMCPion(0),fhPhiMCPion(0),fhEtaMCPion(0),fhLambdaMCPion(0),fhPionTrueE(0),fhRatioEPion(0), -fhEMCProton(0),fhPhiMCProton(0),fhEtaMCProton(0),fhLambdaMCProton(0),fhProtonTrueE(0),fhRatioEProton(0), -fhEMCAntiProton(0),fhPhiMCAntiProton(0),fhEtaMCAntiProton(0),fhLambdaMCAntiProton(0),fhAntiProtonTrueE(0),fhRatioEAntiProton(0), -fhEMCNeutron(0),fhPhiMCNeutron(0),fhEtaMCNeutron(0),fhLambdaMCNeutron(0),fhNeutronTrueE(0),fhRatioENeutron(0), -fhEMCEta(0),fhPhiMCEta(0),fhEtaMCEta(0),fhLambdaMCEta(0), -fhPrimPt(0x0) -{ - //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,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ; - parList+=onePar ; - snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ; - parList+=onePar ; - snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ; - parList+=onePar ; - snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ; - 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->SetXTitle("p_{T}"); - fhNCellCluster->SetYTitle("N_{Cells}"); - outputContainer->Add(fhNCellCluster) ; - - fhPtCluster = new TH1F("hPtCluster","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); - fhPtCluster->SetXTitle("p_{Cluster}(GeV/c)"); - outputContainer->Add(fhPtCluster) ; - - fhPhiCluster = new TH2F - ("hPhiCluster","#phi_{Cluster}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); - fhPhiCluster->SetYTitle("#phi"); - fhPhiCluster->SetXTitle("E_{Cluster} (GeV/c)"); - outputContainer->Add(fhPhiCluster) ; - - fhEtaCluster = new TH2F - ("hEtaCluster","#phi_{Cluster}",nptbins,ptmin,ptmax,netabins,etamin,etamax); - fhEtaCluster->SetYTitle("#eta"); - fhEtaCluster->SetXTitle("E_{Cluster} (GeV/c)"); - outputContainer->Add(fhEtaCluster) ; - - fhDeltaEtaClusters = new TH1F("hDeltaEtaClusters","#Delta #eta of clusters over calorimeter",netabins,etamin,etamax); - fhDeltaEtaClusters->SetXTitle("#Delta #eta_{Clusters}"); - outputContainer->Add(fhDeltaEtaClusters) ; - - fhDeltaPhiClusters = new TH1F("hDeltaPhiClusters","#Delta #phi of clusters over calorimeter",nphibins,phimin,phimax); - fhDeltaPhiClusters->SetXTitle("#Delta #phi_{Clusters}"); - outputContainer->Add(fhDeltaPhiClusters) ; - - fhLambdaCluster = new TH3F - ("hLambdaCluster","#lambda_{Cluster}",nptbins,ptmin,ptmax,200,0,2,200,0,2); - fhLambdaCluster->SetZTitle("#lambda_{1}^{2}"); - fhLambdaCluster->SetYTitle("#lambda_{0}^{2}"); - fhLambdaCluster->SetXTitle("p_{T Cluster} (GeV/c)"); - outputContainer->Add(fhLambdaCluster) ; - - fhELambdaCluster = new TH3F - ("hELambdaCluster","#lambda_{Cluster}",nptbins,ptmin,ptmax,200,0,2,200,0,2); - fhELambdaCluster->SetZTitle("#lambda_{1}^{2}"); - fhELambdaCluster->SetYTitle("#lambda_{0}^{2}"); - fhELambdaCluster->SetXTitle("p_{T Cluster} (GeV)"); - outputContainer->Add(fhELambdaCluster) ; - - fhDispersionCluster = new TH2F - ("hDispersionCluster","Dispersion_{Cluster}",nptbins,ptmin,ptmax,200,0,2); - fhDispersionCluster->SetYTitle("Dispersion"); - fhDispersionCluster->SetXTitle("p_{T Cluster} (GeV/c)"); - outputContainer->Add(fhDispersionCluster) ; - - fhELambdaCellCluster = new TH3F - ("hELambdaCellCluster","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); - fhELambdaCellCluster->SetZTitle("N Cells"); - fhELambdaCellCluster->SetYTitle("#lambda_{0}^{2}"); - fhELambdaCellCluster->SetXTitle("E_{Cluster} (GeV)"); - outputContainer->Add(fhELambdaCellCluster) ; - - //Identified photon histograms - fhNCellPhoton = new TH2F ("hNCellPhoton","Number of cell per cluster",nptbins,ptmin,ptmax,21,-0.5,20.5); - fhNCellPhoton->SetXTitle("E"); - fhNCellPhoton->SetYTitle("N_{Cells}"); - outputContainer->Add(fhNCellPhoton) ; - - fhLambdaPhoton = new TH3F - ("hLambdaPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2); - fhLambdaPhoton->SetZTitle("#lambda_{1}^{2}"); - fhLambdaPhoton->SetYTitle("#lambda_{0}^{2}"); - fhLambdaPhoton->SetXTitle("p_{T Photon} (GeV/c)"); - outputContainer->Add(fhLambdaPhoton) ; - - fhELambdaPhoton = new TH3F - ("hELambdaPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2); - fhELambdaPhoton->SetZTitle("#lambda_{1}^{2}"); - fhELambdaPhoton->SetYTitle("#lambda_{0}^{2}"); - fhELambdaPhoton->SetXTitle("E_{Photon} (GeV)"); - outputContainer->Add(fhELambdaPhoton) ; - - fhDispersionPhoton = new TH2F - ("hLambdaDispersion","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2); - fhDispersionPhoton->SetYTitle("Dispersion"); - fhDispersionPhoton->SetXTitle("p_{T Photon} (GeV/c)"); - outputContainer->Add(fhDispersionPhoton) ; - - fhELambdaCellPhoton = new TH3F - ("hELambdaCellPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); - fhELambdaCellPhoton->SetZTitle("N Cells"); - fhELambdaCellPhoton->SetYTitle("#lambda_{0}^{2}"); - fhELambdaCellPhoton->SetXTitle("E_{Photon} (GeV)"); - outputContainer->Add(fhELambdaCellPhoton) ; - - //Identified Pi0 histograms - fhNCellPi0 = new TH2F ("hNCellPi0","Number of cell per cluster",nptbins,ptmin,ptmax,21,-0.5,20.5); - fhNCellPi0->SetXTitle("E"); - fhNCellPi0->SetYTitle("N_{Cells}"); - outputContainer->Add(fhNCellPi0) ; - - fhLambdaPi0 = new TH3F - ("hLambdaPi0","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2); - fhLambdaPi0->SetZTitle("#lambda_{1}^{2}"); - fhLambdaPi0->SetYTitle("#lambda_{0}^{2}"); - fhLambdaPi0->SetXTitle("p_{T Pi0} (GeV/c)"); - outputContainer->Add(fhLambdaPi0) ; - - fhELambdaPi0 = new TH3F - ("hELambdaPi0","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2); - fhELambdaPi0->SetZTitle("#lambda_{1}^{2}"); - fhELambdaPi0->SetYTitle("#lambda_{0}^{2}"); - fhELambdaPi0->SetXTitle("E_{Pi0} (GeV)"); - outputContainer->Add(fhELambdaPi0) ; - - fhDispersionPi0 = new TH2F - ("hLambdaDispersion","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2); - fhDispersionPi0->SetYTitle("Dispersion"); - fhDispersionPi0->SetXTitle("p_{T Pi0} (GeV/c)"); - outputContainer->Add(fhDispersionPi0) ; - - fhELambdaCellPi0 = new TH3F - ("hELambdaCellPi0","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); - fhELambdaCellPi0->SetZTitle("N Cells"); - fhELambdaCellPi0->SetYTitle("#lambda_{0}^{2}"); - fhELambdaCellPi0->SetXTitle("E_{Pi0} (GeV)"); - outputContainer->Add(fhELambdaCellPi0) ; - - //Identified charged hadron histograms - fhNCellChargedHadron = new TH2F ("hNCellChargedHadron","Number of cell per cluster",nptbins,ptmin,ptmax,21,-0.5,20.5); - fhNCellChargedHadron->SetXTitle("p_{T}"); - fhNCellChargedHadron->SetYTitle("N_{Cells}"); - outputContainer->Add(fhNCellChargedHadron) ; - - fhLambdaChargedHadron = new TH3F - ("hLambdaChargedHadron","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2); - fhLambdaChargedHadron->SetZTitle("#lambda_{1}^{2}"); - fhLambdaChargedHadron->SetYTitle("#lambda_{0}^{2}"); - fhLambdaChargedHadron->SetXTitle("p_{T ChargedHadron} (GeV/c)"); - outputContainer->Add(fhLambdaChargedHadron) ; - - fhELambdaChargedHadron = new TH3F - ("hELambdaChargedHadron","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2); - fhELambdaChargedHadron->SetZTitle("#lambda_{1}^{2}"); - fhELambdaChargedHadron->SetYTitle("#lambda_{0}^{2}"); - fhELambdaChargedHadron->SetXTitle("E_{ChargedHadron} (GeV)"); - outputContainer->Add(fhELambdaChargedHadron) ; - - fhDispersionChargedHadron = new TH2F - ("hLambdaDispersion","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2); - fhDispersionChargedHadron->SetYTitle("Dispersion"); - fhDispersionChargedHadron->SetXTitle("p_{T ChargedHadron} (GeV/c)"); - outputContainer->Add(fhDispersionChargedHadron) ; - - fhELambdaCellChargedHadron = new TH3F - ("hELambdaCellChargedHadron","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); - fhELambdaCellChargedHadron->SetZTitle("N Cells"); - fhELambdaCellChargedHadron->SetYTitle("#lambda_{0}^{2}"); - fhELambdaCellChargedHadron->SetXTitle("E_{ChargedHadron} (GeV)"); - outputContainer->Add(fhELambdaCellChargedHadron) ; - - if(IsDataMC()){ - fhDeltaE = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50); - fhDeltaE->SetXTitle("#Delta E (GeV)"); - outputContainer->Add(fhDeltaE); - - fhDeltaPt = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50); - fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)"); - outputContainer->Add(fhDeltaPt); - - fhRatioE = new TH1F ("hRatioE","Reco/MC E ", 1000,0,2); - fhRatioE->SetXTitle("E_{reco}/E_{gen}"); - outputContainer->Add(fhRatioE); - - fhRatioPt = new TH1F ("hRatioPt","Reco/MC p_{T} ", 1000,0,2); - fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}"); - outputContainer->Add(fhRatioPt); - - fh2E = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); - fh2E->SetXTitle("E_{rec} (GeV)"); - fh2E->SetYTitle("E_{gen} (GeV)"); - outputContainer->Add(fh2E); - - fh2Pt = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax); - fh2Pt->SetXTitle("p_{T,rec} (GeV/c)"); - fh2Pt->SetYTitle("p_{T,gen} (GeV/c)"); - outputContainer->Add(fh2Pt); - - fhMCPdg = new TH2F ("hMCPdg","PDG Code distribution",nptbins,ptmin,ptmax,6001,-3000.5,3000.5); - fhMCPdg->SetXTitle("E_{Reco Cluster} (GeV)"); - fhMCPdg->SetYTitle("PDG Code"); - outputContainer->Add(fhMCPdg); - - fhEMCPhoton = new TH1F("hEMCPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax); - //fhEMCPhoton->SetYTitle("N"); - fhEMCPhoton->SetXTitle("E_{#gamma}(GeV)"); - outputContainer->Add(fhEMCPhoton) ; - - fhPhiMCPhoton = new TH2F - ("hPhiMCPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); - fhPhiMCPhoton->SetYTitle("#phi"); - fhPhiMCPhoton->SetXTitle("E_{ #gamma} (GeV)"); - outputContainer->Add(fhPhiMCPhoton) ; - - fhEtaMCPhoton = new TH2F - ("hEtaMCPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax); - fhEtaMCPhoton->SetYTitle("#eta"); - fhEtaMCPhoton->SetXTitle("E_{ #gamma} (GeV)"); - outputContainer->Add(fhEtaMCPhoton) ; - - fhLambdaMCPhoton = new TH3F - ("hLambdaMCPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); - fhLambdaMCPhoton->SetZTitle("N_{Cells}"); - fhLambdaMCPhoton->SetYTitle("#lambda_{0}^{2}"); - fhLambdaMCPhoton->SetXTitle("E_{#gamma} (GeV)"); - outputContainer->Add(fhLambdaMCPhoton) ; - - fhPhotTrueE = new TH3F - ("hPhotTrueE","#lambda_{#gamma}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2); - fhPhotTrueE->SetZTitle("#lambda_{0}^{2}"); - fhPhotTrueE->SetYTitle("Recons. E_{#gamma}"); - fhPhotTrueE->SetXTitle("MC Truth E_{#gamma} (GeV)"); - outputContainer->Add(fhPhotTrueE) ; - - fhRatioEPhoton = new TH1F ("hRatioEPhoton","E_{Reco #gamma}/E_{MC truth #gamma}", 1000,0,2); - fhRatioEPhoton->SetXTitle("E_{reco #gamma}/E_{gen #gamma}"); - outputContainer->Add(fhRatioEPhoton); - - fhEMCPi0 = new TH1F("hEMCPi0","Number of #pi^0 over calorimeter",nptbins,ptmin,ptmax); - fhEMCPi0->SetYTitle("N"); - fhEMCPi0->SetXTitle("E_{ #pi^0}(GeV)"); - outputContainer->Add(fhEMCPi0) ; - - fhPhiMCPi0 = new TH2F - ("hPhiMCPi0","#phi_{#pi^0}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); - fhPhiMCPi0->SetYTitle("#phi"); - fhPhiMCPi0->SetXTitle("E_{ #pi^0} (GeV)"); - outputContainer->Add(fhPhiMCPi0) ; - - fhEtaMCPi0 = new TH2F - ("hEtaMCPi0","#phi_{#pi^0}",nptbins,ptmin,ptmax,netabins,etamin,etamax); - fhEtaMCPi0->SetYTitle("#eta"); - fhEtaMCPi0->SetXTitle("E_{ #pi^0} (GeV)"); - outputContainer->Add(fhEtaMCPi0) ; - - fhLambdaMCPi0 = new TH3F - ("hLambdaMCPi0","#lambda_{#pi^0}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); - fhLambdaMCPi0->SetZTitle("N_{Cells}"); - fhLambdaMCPi0->SetYTitle("#lambda_{0}^{2}"); - fhLambdaMCPi0->SetXTitle("E_{ #pi^0} (GeV)"); - outputContainer->Add(fhLambdaMCPi0) ; - - fhProductionDistance = new TH2F - ("hProductionDistance","#lambda_{#pi^0}",nptbins,ptmin,ptmax,160,0,800); - fhProductionDistance->SetYTitle("Distance of production of Pi0 from vertex (cm)"); - fhProductionDistance->SetXTitle("E_{ #pi^0} (GeV)"); - outputContainer->Add(fhProductionDistance) ; - - fhProductionRadius = new TH2F - ("hProductionRadius","#lambda_{#pi^0}",nptbins,ptmin,ptmax,160,0,800); - fhProductionRadius->SetYTitle("Radius of production of Pi0 from vertex (cm)"); - fhProductionRadius->SetXTitle("E_{ #pi^0} (GeV)"); - outputContainer->Add(fhProductionRadius) ; - - fhDecayAngle = new TH2F - ("hDecayAngle","#lambda_{#pi^0}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); - fhDecayAngle->SetYTitle("Opening angle of the #pi^{0} decay (radians)"); - fhDecayAngle->SetXTitle("E_{ #pi^0} (GeV)"); - outputContainer->Add(fhDecayAngle) ; - - fhPi0TrueE = new TH3F - ("hPi0TrueE","#lambda_{#pi^0}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2); - fhPi0TrueE->SetZTitle("#lambda_{0}^{2}"); - fhPi0TrueE->SetYTitle("Recons. E_{#pi^0}"); - fhPi0TrueE->SetXTitle("MC Truth E_{#pi^0} (GeV)"); - outputContainer->Add(fhPi0TrueE) ; - - fhRatioEPi0 = new TH1F ("hRatioEPi0","E_{Reco #pi^{0}}/E_{MC truth #pi^{0}}", 1000,0,2); - fhRatioEPi0->SetXTitle("E_{reco #pi^{0}}/E_{gen #pi^{0}}"); - outputContainer->Add(fhRatioEPi0); - - fhEMCPion = new TH1F("hEMCPion","Number of #pi over calorimeter",nptbins,ptmin,ptmax); - fhEMCPion->SetYTitle("N"); - fhEMCPion->SetXTitle("E_{ #pi}(GeV)"); - outputContainer->Add(fhEMCPion) ; - - fhPhiMCPion = new TH2F - ("hPhiMCPion","#phi_{#pi}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); - fhPhiMCPion->SetYTitle("#phi"); - fhPhiMCPion->SetXTitle("E_{ #pi} (GeV)"); - outputContainer->Add(fhPhiMCPion) ; - - fhEtaMCPion = new TH2F - ("hEtaMCPion","#phi_{#pi}",nptbins,ptmin,ptmax,netabins,etamin,etamax); - fhEtaMCPion->SetYTitle("#eta"); - fhEtaMCPion->SetXTitle("E_{ #pi} (GeV)"); - outputContainer->Add(fhEtaMCPion) ; - - fhLambdaMCPion = new TH3F - ("hLambdaMCPion","#lambda_{#pi}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); - fhLambdaMCPion->SetZTitle("N_{Cells}"); - fhLambdaMCPion->SetYTitle("#lambda_{0}^{2}"); - fhLambdaMCPion->SetXTitle("E_{ #pi} (GeV)"); - outputContainer->Add(fhLambdaMCPion) ; - - fhPionTrueE = new TH3F - ("hPionTrueE","#lambda_{#pi}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2); - fhPionTrueE->SetZTitle("#lambda_{0}^{2}"); - fhPionTrueE->SetYTitle("Recons. E_{#pi}"); - fhPionTrueE->SetXTitle("MC Truth E_{#pi} (GeV)"); - outputContainer->Add(fhPionTrueE) ; - - fhRatioEPion = new TH1F ("hRatioEPion","E_{Reco #pi^{#pm}}/E_{MC truth #pi^{#pm}}", 1000,0,2); - fhRatioEPion->SetXTitle("E_{reco #pi^{#pm}}/E_{gen #pi^{#pm}}"); - outputContainer->Add(fhRatioEPion); - - fhEMCProton = new TH1F("hEMCProton","Number of p over calorimeter",nptbins,ptmin,ptmax); - fhEMCProton->SetYTitle("N"); - fhEMCProton->SetXTitle("E_{ p}(GeV)"); - outputContainer->Add(fhEMCProton) ; - - fhPhiMCProton = new TH2F - ("hPhiMCProton","#phi_{p}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); - fhPhiMCProton->SetYTitle("#phi"); - fhPhiMCProton->SetXTitle("E_{ p} (GeV)"); - outputContainer->Add(fhPhiMCProton) ; - - fhEtaMCProton = new TH2F - ("hEtaMCProton","#phi_{p}",nptbins,ptmin,ptmax,netabins,etamin,etamax); - fhEtaMCProton->SetYTitle("#eta"); - fhEtaMCProton->SetXTitle("E_{ p} (GeV)"); - outputContainer->Add(fhEtaMCProton) ; - - fhLambdaMCProton = new TH3F - ("hLambdaMCProton","#lambda_{p}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); - fhLambdaMCProton->SetZTitle("N_{Cells}"); - fhLambdaMCProton->SetYTitle("#lambda_{0}^{2}"); - fhLambdaMCProton->SetXTitle("E_{ p} (GeV)"); - outputContainer->Add(fhLambdaMCProton) ; - - fhProtonTrueE = new TH3F - ("hProtonTrueE","#lambda_{p}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2); - fhProtonTrueE->SetZTitle("#lambda_{0}^{2}"); - fhProtonTrueE->SetYTitle("Recons. E_{p}"); - fhProtonTrueE->SetXTitle("MC Truth E_{p} (GeV)"); - outputContainer->Add(fhProtonTrueE) ; - - fhRatioEProton = new TH1F ("hRatioEProton","E_{Reco p}/E_{MC truth p}", 1000,0,2); - fhRatioEProton->SetXTitle("E_{reco p}/E_{gen p}"); - outputContainer->Add(fhRatioEProton); - - fhEMCAntiProton = new TH1F("hEMCAntiProton","Number of #bar{p} over calorimeter",nptbins,ptmin,ptmax); - fhEMCAntiProton->SetYTitle("N"); - fhEMCAntiProton->SetXTitle("E_{#bar{p}}(GeV)"); - outputContainer->Add(fhEMCAntiProton) ; - - fhPhiMCAntiProton = new TH2F - ("hPhiMCAntiProton","#phi_{#bar{p}}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); - fhPhiMCAntiProton->SetYTitle("#phi"); - fhPhiMCAntiProton->SetXTitle("E_{#bar{p}} (GeV)"); - outputContainer->Add(fhPhiMCAntiProton) ; - - fhEtaMCAntiProton = new TH2F - ("hEtaMCAntiProton","#phi_{#bar{p}}",nptbins,ptmin,ptmax,netabins,etamin,etamax); - fhEtaMCAntiProton->SetYTitle("#eta"); - fhEtaMCAntiProton->SetXTitle("E_{#bar{p}} (GeV)"); - outputContainer->Add(fhEtaMCAntiProton) ; - - fhLambdaMCAntiProton = new TH3F - ("hLambdaMCAntiProton","#lambda_{#bar{p}}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); - fhLambdaMCAntiProton->SetZTitle("N_{Cells}"); - fhLambdaMCAntiProton->SetYTitle("#lambda_{0}^{2}"); - fhLambdaMCAntiProton->SetXTitle("E_{#bar{p}} (GeV)"); - outputContainer->Add(fhLambdaMCAntiProton) ; - - fhAntiProtonTrueE = new TH3F - ("hAntiProtonTrueE","#lambda_{#bar{p}}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2); - fhAntiProtonTrueE->SetZTitle("#lambda_{0}^{2}"); - fhAntiProtonTrueE->SetYTitle("Recons. E_{#bar{p}}"); - fhAntiProtonTrueE->SetXTitle("MC Truth E_{#bar{p}} (GeV)"); - outputContainer->Add(fhAntiProtonTrueE) ; - - fhRatioEAntiProton = new TH1F ("hRatioEAntiProton","E_{Reco #bar{p}}/E_{MC truth #bar{p}}", 1000,0,2); - fhRatioEAntiProton->SetXTitle("E_{reco #bar{p}}/E_{gen #bar{p}}"); - outputContainer->Add(fhRatioEAntiProton); - - fhEMCNeutron = new TH1F("hEMCNeutron","Number of p over calorimeter",nptbins,ptmin,ptmax); - fhEMCNeutron->SetYTitle("N"); - fhEMCNeutron->SetXTitle("E_{ p}(GeV)"); - outputContainer->Add(fhEMCNeutron) ; - - fhPhiMCNeutron = new TH2F - ("hPhiMCNeutron","#phi_{p}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); - fhPhiMCNeutron->SetYTitle("#phi"); - fhPhiMCNeutron->SetXTitle("E_{ p} (GeV)"); - outputContainer->Add(fhPhiMCNeutron) ; - - fhEtaMCNeutron = new TH2F - ("hEtaMCNeutron","#phi_{p}",nptbins,ptmin,ptmax,netabins,etamin,etamax); - fhEtaMCNeutron->SetYTitle("#eta"); - fhEtaMCNeutron->SetXTitle("E_{ p} (GeV)"); - outputContainer->Add(fhEtaMCNeutron) ; - - fhLambdaMCNeutron = new TH3F - ("hLambdaMCNeutron","#lambda_{p}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); - fhLambdaMCNeutron->SetZTitle("N_{Cells}"); - fhLambdaMCNeutron->SetYTitle("#lambda_{0}^{2}"); - fhLambdaMCNeutron->SetXTitle("E_{ p} (GeV)"); - outputContainer->Add(fhLambdaMCNeutron) ; - - fhNeutronTrueE = new TH3F - ("hNeutronTrueE","#lambda_{n}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2); - fhNeutronTrueE->SetZTitle("#lambda_{0}^{2}"); - fhNeutronTrueE->SetYTitle("Recons. E_{n}"); - fhNeutronTrueE->SetXTitle("MC Truth E_{n} (GeV)"); - outputContainer->Add(fhNeutronTrueE) ; - - fhRatioENeutron = new TH1F ("hRatioENeutron","E_{Reco n}/E_{MC truth n}", 1000,0,2); - fhRatioENeutron->SetXTitle("E_{reco n}/E_{gen n}"); - outputContainer->Add(fhRatioENeutron); - - fhEMCEta = new TH1F("hEMCEta","Number of #eta over calorimeter",nptbins,ptmin,ptmax); - fhEMCEta->SetYTitle("N"); - fhEMCEta->SetXTitle("E_{ #eta}(GeV)"); - outputContainer->Add(fhEMCEta) ; - - fhPhiMCEta = new TH2F - ("hPhiMCEta","#phi_{#eta}",nptbins,ptmin,ptmax,nphibins,phimin,phimax); - fhPhiMCEta->SetYTitle("#phi"); - fhPhiMCEta->SetXTitle("E_{ #eta} (GeV)"); - outputContainer->Add(fhPhiMCEta) ; - - fhEtaMCEta = new TH2F - ("hEtaMCEta","#phi_{#eta}",nptbins,ptmin,ptmax,netabins,etamin,etamax); - fhEtaMCEta->SetYTitle("#eta"); - fhEtaMCEta->SetXTitle("E_{ #eta} (GeV)"); - outputContainer->Add(fhEtaMCEta) ; - - fhLambdaMCEta = new TH3F - ("hLambdaMCEta","#lambda_{#eta}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5); - fhLambdaMCEta->SetZTitle("N_{Cells}"); - fhLambdaMCEta->SetYTitle("#lambda_{0}^{2}"); - fhLambdaMCEta->SetXTitle("E_{ #eta} (GeV)"); - outputContainer->Add(fhLambdaMCEta) ; - - fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",nptbins,ptmin,ptmax) ; - outputContainer->Add(fhPrimPt) ; - - }//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("AnaPhoton_"); - - fCalorimeter = "PHOS" ; - fMinDist = 2.; - fMinDist2 = 4.; - fMinDist3 = 5.; - fMassCut = 0.03; //30 MeV - - fTimeCutMin = -1; - fTimeCutMax = 9999999; - fNCellsCut = 0; - - fRejectTrackMatch = kTRUE ; - fCheckConversion = kFALSE; - fAddConvertedPairsToAOD = kFALSE; - - fNumClusters = -1; // By Default, select all events. - -} - -//__________________________________________________________________ -void AliAnaShowerParameter::MakeAnalysisFillAOD() -{ - //Do analysis and fill aods - //Search for photons in fCalorimeter - - //Get vertex for photon momentum calculation - - for (Int_t iev = 0; iev < GetNMixedEvent(); iev++) { - if (!GetMixedEvent()) - GetReader()->GetVertex(GetVertex(iev)); - else - GetMixedEvent()->GetVertexOfEvent(iev)->GetXYZ(GetVertex(iev)); - } - - //Select the Calorimeter of the photon - TObjArray * pl = 0x0; - if(fCalorimeter == "PHOS") - pl = GetPHOSClusters(); - else if (fCalorimeter == "EMCAL") - pl = GetEMCALClusters(); - - if(!pl){ - printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Careful cluster array NULL!!\n"); - return; - } - - //Fill AODCaloClusters and AODParticle with PHOS/EMCAL aods - TLorentzVector mom, mom2 ; - Int_t nCaloClusters = pl->GetEntriesFast(); - //Cut on the number of clusters in the event. - if ((fNumClusters !=-1) && (nCaloClusters != fNumClusters)) return; - Bool_t * indexConverted = new Bool_t[nCaloClusters]; - for (Int_t i = 0; i < nCaloClusters; i++) - indexConverted[i] = kFALSE; - - for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){ - - AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo)); - Int_t evtIndex = 0 ; - if (GetMixedEvent()) { - evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; - } - //Cluster selection, not charged, with photon id and in fiducial cut - - //Input from second AOD? - Int_t input = 0; - - //Get Momentum vector, - if (input == 0) - calo->GetMomentum(mom,GetVertex(evtIndex)) ;//Assume that come from vertex in straight line - - //Skip the cluster if it doesn't fit inside the cuts. - if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ; - Double_t tof = calo->GetTOF()*1e9; - if(tof < fTimeCutMin || tof > fTimeCutMax) continue; - if(calo->GetNCells() <= fNCellsCut) continue; - - //Check acceptance selection - if(IsFiducialCutOn()){ - Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ; - if(! in ) continue ; - } - - //Create AOD for analysis - AliAODPWG4Particle aodph = AliAODPWG4Particle(mom); - Int_t label = calo->GetLabel(); - aodph.SetLabel(label); - aodph.SetInputFileIndex(input); - - //Set the indices of the original caloclusters - aodph.SetCaloLabel(calo->GetID(),-1); - aodph.SetDetector(fCalorimeter); - if(GetDebug() > 1) - printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta()); - - //Check Distance to Bad channel, set bit. - Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel - if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ; - if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm) - continue ; - - if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad); - - if (distBad > fMinDist3) aodph.SetDistToBad(2) ; - else if(distBad > fMinDist2) aodph.SetDistToBad(1) ; - else aodph.SetDistToBad(0) ; - - //Skip matched clusters with tracks - if(fRejectTrackMatch && IsTrackMatched(calo)) continue ; - if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - TrackMatching cut passed \n"); - - //Set PID bits for later selection (AliAnaPi0 for example) - //GetPDG already called in SetPIDBits. - GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils()); - if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - PID Bits set \n"); - - //Play with the MC stack if available - //Check origin of the candidates - if(IsDataMC()){ - aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex())); - if(GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag()); - } - - // Check if cluster comes from a conversion in the material in front of the calorimeter - // Do invariant mass of all pairs, if mass is close to 0, then it is conversion. - - if(fCheckConversion && nCaloClusters > 1){ - Bool_t bConverted = kFALSE; - Int_t id2 = -1; - - //Check if set previously as converted couple, if so skip its use. - if (indexConverted[icalo]) continue; - - for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) { - //Check if set previously as converted couple, if so skip its use. - if (indexConverted[jcalo]) continue; - //printf("Check Conversion indeces %d and %d\n",icalo,jcalo); - AliAODCaloCluster * calo2 = (AliAODCaloCluster*) (pl->At(jcalo)); //Get cluster kinematics - Int_t evtIndex2 = 0 ; - if (GetMixedEvent()) { - evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ; - } - calo2->GetMomentum(mom2,GetVertex(evtIndex2)); - //Check only certain regions - Bool_t in2 = kTRUE; - if(IsFiducialCutOn()) in2 = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ; - if(!in2) continue; - - //Get mass of pair, if small, take this pair. - //printf("\t both in calo, mass %f, cut %f\n",(mom+mom2).M(),fMassCut); - if((mom+mom2).M() < fMassCut){ - bConverted = kTRUE; - id2 = calo2->GetID(); - indexConverted[jcalo]=kTRUE; - break; - } - - }//Mass loop - - if(bConverted){ - if(fAddConvertedPairsToAOD){ - //Create AOD of pair analysis - TLorentzVector mpair = mom+mom2; - AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair); - aodpair.SetLabel(aodph.GetLabel()); - aodpair.SetInputFileIndex(input); - - //printf("Index %d, Id %d\n",icalo, calo->GetID()); - //Set the indeces of the original caloclusters - aodpair.SetCaloLabel(calo->GetID(),id2); - aodpair.SetDetector(fCalorimeter); - aodpair.SetPdg(aodph.GetPdg()); - aodpair.SetTag(aodph.GetTag()); - - //Add AOD with pair object to aod branch - AddAODParticle(aodpair); - //printf("\t \t both added pair\n"); - } - - //Do not add the current calocluster - continue; - }//converted pair - }//check conversion - - //Add AOD with photon object to aod branch - AddAODParticle(aodph); - - }//loop; - delete [] indexConverted; - - if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast()); - -} - -//__________________________________________________________________ -void AliAnaShowerParameter::MakeAnalysisFillHistograms() -{ - - //Do analysis and fill histograms - - // Access MC information in stack if requested, check that it exists. - AliStack * stack = 0x0; - TParticle * primary = 0x0; - TClonesArray * mcparticles0 = 0x0; - //TClonesArray * mcparticles1 = 0x0; - AliAODMCParticle * aodprimary = 0x0; - TObjArray * pl = 0x0; - Int_t iNumCell=0; - - //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 - Int_t naod = GetOutputAODBranch()->GetEntriesFast(); - fhNClusters->Fill(naod); - if(GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod); - - for(Int_t iaod = 0; iaod < naod ; iaod++){ - AliAODPWG4Particle* ph = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod)); - - if(ph->GetDetector() != fCalorimeter) continue; - - if(GetDebug() > 2) - printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ; - - //Fill Cluster histograms - Float_t ptcluster = ph->Pt(); - Float_t phicluster = ph->Phi(); - Float_t etacluster = ph->Eta(); - Float_t ecluster = ph->E(); - Float_t lambdaMainCluster = -1; //Can't set their values here, but need a default -1 value in case something is wrong - Float_t lambdaSecondCluster = -1; - Float_t dispcluster = -1; - - //Get the list of clusters from the AOD and loop over them to find the onces corresponding to the reconstructed 'photons'. - if(fCalorimeter == "PHOS") - pl = GetPHOSClusters(); - else if (fCalorimeter == "EMCAL") - pl = GetEMCALClusters(); - if(pl){ - //Some values are stored in AliAODCaloCluster objects only; we need to fetch them. - for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){ - AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo)); - if (calo->GetLabel()==ph->GetLabel()) { //The Cluster is the right one for this particle - lambdaMainCluster = calo->GetM02(); //lambda_0 - lambdaSecondCluster = calo->GetM20(); //lambda_1 - dispcluster = calo->GetDispersion(); //Dispersion - iNumCell = calo->GetNCells(); - if(GetDebug() > 2) - printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Cluster Lambda0: %3.2f, Lambda1: %3.2f, Dispersion: %3.2f, NCells: %d \n",lambdaMainCluster,lambdaSecondCluster,dispcluster,iNumCell) ; - } - } - } - - //Fill the basic non-MC information about the cluster. - fhPtCluster ->Fill(ecluster); - fhPhiCluster ->Fill(ecluster,phicluster); - fhEtaCluster ->Fill(ecluster,etacluster); - fhDispersionCluster->Fill(ptcluster,dispcluster); - fhLambdaCluster ->Fill(ptcluster,lambdaMainCluster,lambdaSecondCluster); - fhELambdaCluster ->Fill(ecluster,lambdaMainCluster,lambdaSecondCluster); - fhELambdaCellCluster ->Fill(ecluster,lambdaMainCluster,iNumCell); - fhNCellCluster->Fill(ecluster,iNumCell); - - //In the case of an event with several clusters, calculate DeltaEta and DeltaPhi for the cluster pairs. - for (Int_t jaod=iaod+1;jaodAt(jaod)); - if(phSecond->GetDetector() != fCalorimeter) continue; - fhDeltaPhiClusters->Fill(phicluster-(phSecond->Phi())); - fhDeltaEtaClusters->Fill(etacluster-(phSecond->Eta())); - } - - //Fill the lambda distribution for identified particles - if(ph->GetPdg() == AliCaloPID::kPhoton){ - fhDispersionPhoton->Fill(ptcluster,dispcluster); - fhLambdaPhoton ->Fill(ptcluster,lambdaMainCluster,lambdaSecondCluster); - fhELambdaPhoton ->Fill(ecluster,lambdaMainCluster,lambdaSecondCluster); - fhELambdaCellPhoton ->Fill(ecluster,lambdaMainCluster,iNumCell); - fhNCellPhoton->Fill(ecluster,iNumCell); - } - if(ph->GetPdg() == AliCaloPID::kPi0){ - fhDispersionPi0->Fill(ptcluster,dispcluster); - fhLambdaPi0 ->Fill(ptcluster,lambdaMainCluster,lambdaSecondCluster); - fhELambdaPi0 ->Fill(ecluster,lambdaMainCluster,lambdaSecondCluster); - fhELambdaCellPi0 ->Fill(ecluster,lambdaMainCluster,iNumCell); - fhNCellPi0->Fill(ecluster,iNumCell); - } - if(ph->GetPdg() == AliCaloPID::kChargedHadron){ - fhDispersionChargedHadron->Fill(ptcluster,dispcluster); - fhLambdaChargedHadron ->Fill(ptcluster,lambdaMainCluster,lambdaSecondCluster); - fhELambdaChargedHadron ->Fill(ecluster,lambdaMainCluster,lambdaSecondCluster); - fhELambdaCellChargedHadron ->Fill(ecluster,lambdaMainCluster,iNumCell); - fhNCellChargedHadron->Fill(ecluster,iNumCell); - } - - //Play with the MC data if available - if(IsDataMC()){ - if(GetReader()->ReadStack() && !stack) return; - if(GetReader()->ReadAODMCParticles() && !mcparticles0) return; - - //Get the tag from AliMCAnalysisUtils for PID - Int_t tag =ph->GetTag(); - if ( ph->GetLabel() < 0){ - if(GetDebug() > 0) - printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Label is -1; problem in the MC ESD? "); - continue; - } - - - //Check if the tag matches on of the different particle types and fill the corresponding histograms - if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)&&!(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))) //kMCPhoton; making sure that this is not a Pi0 cluster (it should not happen). - { - fhEMCPhoton ->Fill(ecluster); - fhPhiMCPhoton ->Fill(ecluster,phicluster); - fhEtaMCPhoton ->Fill(ecluster,etacluster); - fhLambdaMCPhoton ->Fill(ecluster,lambdaMainCluster,iNumCell); - Double_t ePhot; - //Get the true energy of the photon - Int_t iCurrent = ph->GetLabel(); - TParticle* pCurrent = stack->Particle(iCurrent); - ePhot = pCurrent->Energy(); - fhPhotTrueE->Fill(ePhot,ecluster,lambdaMainCluster); - fhRatioEPhoton->Fill(ecluster/ePhot); - fhMCPdg->Fill(ecluster,22); - if(GetDebug() > 3) - printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCPhoton: True E: %3.2f, Reco E: %3.2f, Lambda0: %3.2f\n",ePhot,ecluster,lambdaMainCluster); - }//kMCPhoton - - if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)) //kMCPi0 : single cluster created by 2 photons from the same Pi0 decay - { - - //Fill the basic information about the Pi0 cluster - fhEMCPi0 ->Fill(ecluster); - fhPhiMCPi0 ->Fill(ecluster,phicluster); - fhEtaMCPi0 ->Fill(ecluster,etacluster); - fhLambdaMCPi0 ->Fill(ecluster,lambdaMainCluster,iNumCell); - - //Check the position of the Pi0: does it come from the vertex or was it created by hadron-material collision? - Int_t iCurrent = ph->GetLabel(); - Double_t ePi0 = 0; - TParticle* pCurrent = stack->Particle(iCurrent); - while (pCurrent->GetPdgCode()!=111) - { - iCurrent = pCurrent->GetFirstMother(); - pCurrent = stack->Particle(iCurrent); - } - Float_t fDistance = pCurrent->Rho(); - Float_t fRadius = pCurrent->R(); - ePi0 = pCurrent->Energy(); - fhProductionDistance->Fill(ecluster,fDistance); - fhProductionRadius->Fill(ecluster,fRadius); - //Compare the cluster energy and true energy - fhPi0TrueE->Fill(ePi0,ecluster,lambdaMainCluster); - fhRatioEPi0->Fill(ecluster/ePi0); - fhMCPdg->Fill(ecluster,111); - if(GetDebug() > 3) - printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCPi0: True E: %3.2f, Reco E: %3.2f, Lambda0: %3.2f\n",ePi0,ecluster,lambdaMainCluster); - } - if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion))//kMCPion - { - fhEMCPion ->Fill(ecluster); - fhPhiMCPion ->Fill(ecluster,phicluster); - fhEtaMCPion ->Fill(ecluster,etacluster); - fhLambdaMCPion ->Fill(ecluster,lambdaMainCluster,iNumCell); - - Int_t iCurrent = ph->GetLabel(); - Double_t ePion = 0; - TParticle* pCurrent = stack->Particle(iCurrent); - while ((TMath::Abs(pCurrent->GetPdgCode())!=211)&&(iCurrent>=0)) - { - if (iCurrent!=pCurrent->GetFirstMother()) iCurrent = pCurrent->GetFirstMother(); //Avoiding infinite loops - if (iCurrent>=0) pCurrent = stack->Particle(iCurrent); - } - if ((TMath::Abs(pCurrent->GetPdgCode())==211)&&(iCurrent>=0)) - { - ePion = pCurrent->Energy(); - fhPionTrueE->Fill(ePion,ecluster,lambdaMainCluster); - fhRatioEPion->Fill(ecluster/ePion); - } - fhMCPdg->Fill(ecluster,211); - if(GetDebug() > 3) - printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCPion: True E: %3.2f, Reco E: %3.2f, Lambda0: %3.2f\n",ePion,ecluster,lambdaMainCluster); - } - if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton))//kMCProton - { - fhEMCProton ->Fill(ecluster); - fhPhiMCProton ->Fill(ecluster,phicluster); - fhEtaMCProton ->Fill(ecluster,etacluster); - fhLambdaMCProton ->Fill(ecluster,lambdaMainCluster,iNumCell); - - Int_t iCurrent = ph->GetLabel(); - Double_t eProton = 0; - TParticle* pCurrent = stack->Particle(iCurrent); - while ((TMath::Abs(pCurrent->GetPdgCode())!=2212)&&(iCurrent>=0)) - { - if (iCurrent!=pCurrent->GetFirstMother()) iCurrent = pCurrent->GetFirstMother(); //Avoiding infinite loops - if (iCurrent>=0) pCurrent = stack->Particle(iCurrent); - } - if ((TMath::Abs(pCurrent->GetPdgCode())==2212)&&(iCurrent>=0)) - { - eProton = pCurrent->Energy(); - fhProtonTrueE->Fill(eProton,ecluster,lambdaMainCluster); - fhRatioEProton->Fill(ecluster/eProton); - } - fhMCPdg->Fill(ecluster,2212); - if(GetDebug() > 3) - printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCProton: True E: %3.2f, Reco E: %3.2f, Lambda0: %3.2f\n",eProton,ecluster,lambdaMainCluster); - } - if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))//kMCAntiProton - { - fhEMCAntiProton ->Fill(ecluster); - fhPhiMCAntiProton ->Fill(ecluster,phicluster); - fhEtaMCAntiProton ->Fill(ecluster,etacluster); - fhLambdaMCAntiProton ->Fill(ecluster,lambdaMainCluster,iNumCell); - - Int_t iCurrent = ph->GetLabel(); - Double_t eAntiProton = 0; - TParticle* pCurrent = stack->Particle(iCurrent); - while ((TMath::Abs(pCurrent->GetPdgCode())!=-2212)&&(iCurrent>=0)) - { - if (iCurrent!=pCurrent->GetFirstMother()) iCurrent = pCurrent->GetFirstMother(); //Avoiding infinite loops - if (iCurrent>=0) pCurrent = stack->Particle(iCurrent); - } - if ((TMath::Abs(pCurrent->GetPdgCode())==-2212)&&(iCurrent>=0)) - { - eAntiProton = pCurrent->Energy(); - fhAntiProtonTrueE->Fill(eAntiProton,ecluster,lambdaMainCluster); - fhRatioEAntiProton->Fill(ecluster/eAntiProton); - } - fhMCPdg->Fill(ecluster,-2212); - if(GetDebug() > 3) - printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCAntiProton: True E: %3.2f, Reco E: %3.2f, Lambda0: %3.2f\n",eAntiProton,ecluster,lambdaMainCluster); - } - if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCNeutron))//kMCNeutron - { - fhEMCNeutron ->Fill(ecluster); - fhPhiMCNeutron ->Fill(ecluster,phicluster); - fhEtaMCNeutron ->Fill(ecluster,etacluster); - fhLambdaMCNeutron ->Fill(ecluster,lambdaMainCluster,iNumCell); - - Int_t iCurrent = ph->GetLabel(); - Double_t eNeutron = 0; - TParticle* pCurrent = stack->Particle(iCurrent); - while ((TMath::Abs(pCurrent->GetPdgCode())!=2112)&&(iCurrent>=0)) - { - if (iCurrent!=pCurrent->GetFirstMother()) iCurrent = pCurrent->GetFirstMother(); //Avoiding infinite loops - if (iCurrent>=0) pCurrent = stack->Particle(iCurrent); - } - if ((TMath::Abs(pCurrent->GetPdgCode())==2112)&&(iCurrent>=0)) - { - eNeutron = pCurrent->Energy(); - fhNeutronTrueE->Fill(eNeutron,ecluster,lambdaMainCluster); - fhRatioENeutron->Fill(ecluster/eNeutron); - } - fhMCPdg->Fill(ecluster,2112); - if(GetDebug() > 3) - printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCNeutron: True E: %3.2f, Reco E: %3.2f, Lambda0: %3.2f\n",eNeutron,ecluster,lambdaMainCluster); - } - if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))//kMCEta - { - fhEMCEta ->Fill(ecluster); - fhPhiMCEta ->Fill(ecluster,phicluster); - fhEtaMCEta ->Fill(ecluster,etacluster); - fhLambdaMCEta ->Fill(ecluster,lambdaMainCluster,iNumCell); - } - - // Access MC information in stack if requested, check that it exists. - Int_t label =ph->GetLabel(); - if(label < 0) { - printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label); - continue; - } - - //Calculate Pi0 decay angles - if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){ - for(Int_t i=0 ; iGetNprimary(); i++){ - TParticle * prim = stack->Particle(i) ; - if(prim->GetPdgCode() == 111){ - TLorentzVector mom1 ; - (stack->Particle((prim->GetFirstDaughter())))->Momentum(mom1); - TLorentzVector mom2 ; - (stack->Particle((prim->GetLastDaughter())))->Momentum(mom2); - Float_t fDecayAngle = mom1.Angle(mom2.Vect()); - fhDecayAngle->Fill(ecluster,fDecayAngle); - } - } - } - - - Float_t eprim = 0; - Float_t ptprim = 0; - if(GetReader()->ReadStack()){ - - if(label >= stack->GetNtrack()) { - if(GetDebug() > 2) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack()); - continue ; - } - - primary = stack->Particle(label); - if(!primary){ - printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label); - continue; - } - eprim = primary->Energy(); - ptprim = primary->Pt(); - } - else if(GetReader()->ReadAODMCParticles()){ - //Check which is the input - if(ph->GetInputFileIndex() == 0){ - if(!mcparticles0) continue; - if(label >= mcparticles0->GetEntriesFast()) { - if(GetDebug() > 2) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",label, mcparticles0->GetEntriesFast()); - continue ; - } - //Get the particle - aodprimary = (AliAODMCParticle*) mcparticles0->At(label); - - } - // else {//Second input - // if(!mcparticles1) continue; - // if(label >= mcparticles1->GetEntriesFast()) { - // if(GetDebug() > 2) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** large label ***: label %d, n tracks %d \n",label, mcparticles1->GetEntriesFast()); - // continue ; - // } - // //Get the particle - // aodprimary = (AliAODMCParticle*) mcparticles1->At(label); - // }//second input - - if(!aodprimary){ - printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** no primary ***: label %d \n", label); - continue; - } - - eprim = aodprimary->E(); - ptprim = aodprimary->Pt(); - } - //Write down MC truth information concerning all clusters - fh2E ->Fill(ecluster, eprim); - fh2Pt ->Fill(ptcluster, ptprim); - fhDeltaE ->Fill(eprim-ecluster); - fhDeltaPt->Fill(ptprim-ptcluster); - if(eprim > 0) fhRatioE ->Fill(ecluster/eprim); - if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim); - }//Histograms with MC - }// aod loop -} - - -//__________________________________________________________________ -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 Distance to Bad Channel = %2.1f\n",fMinDist); - printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2); - printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3); - printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch); - printf("Check Pair Conversion = %d\n",fCheckConversion); - printf("Add conversion pair to AOD = %d\n",fAddConvertedPairsToAOD); - printf("Conversion pair mass cut = %f\n",fMassCut); - printf("Time Cut: %3.1f < TOF < %3.1f\n", fTimeCutMin, fTimeCutMax); - printf("Number of cells in cluster is > %f \n", fNCellsCut); - printf("Number of clusters in envent is > %d \n", fNumClusters); - printf(" \n") ; - -} +/************************************************************************** + * 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 +#include +#include +//#include +#include +#include "TParticle.h" +//#include + +// --- 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 "AliEMCALGeoUtils.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::MakeAnalysisFillAOD() +{ + /*//Do analysis and fill aods + //Search for photons in fCalorimeter + + //Get vertex for photon momentum calculation + + for (Int_t iev = 0; iev < GetNMixedEvent(); iev++) { + if (!GetMixedEvent()) + GetReader()->GetVertex(GetVertex(iev)); + else + GetMixedEvent()->GetVertexOfEvent(iev)->GetXYZ(GetVertex(iev)); + } + + //Select the Calorimeter of the photon + TObjArray * pl = 0x0; + if(fCalorimeter == "PHOS") + pl = GetPHOSClusters(); + else if (fCalorimeter == "EMCAL") + pl = GetEMCALClusters(); + + if(!pl){ + printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Careful cluster array NULL!!\n"); + return; + } + + //Fill AODCaloClusters and AODParticle with PHOS/EMCAL aods + TLorentzVector mom, mom2 ; + Int_t nCaloClusters = pl->GetEntriesFast(); + //Cut on the number of clusters in the event. + if ((fNumClusters !=-1) && (nCaloClusters != fNumClusters)) return; + Bool_t * indexConverted = new Bool_t[nCaloClusters]; + for (Int_t i = 0; i < nCaloClusters; i++) + indexConverted[i] = kFALSE; + + for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){ + + AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo)); + Int_t evtIndex = 0 ; + if (GetMixedEvent()) { + evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; + } + //Cluster selection, not charged, with photon id and in fiducial cut + + //Input from second AOD? + Int_t input = 0; + + //Get Momentum vector, + if (input == 0) + calo->GetMomentum(mom,GetVertex(evtIndex)) ;//Assume that come from vertex in straight line + + //Skip the cluster if it doesn't fit inside the cuts. + if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ; + Double_t tof = calo->GetTOF()*1e9; + if(tof < fTimeCutMin || tof > fTimeCutMax) continue; + if(calo->GetNCells() <= fNCellsCut) continue; + + //Check acceptance selection + if(IsFiducialCutOn()){ + Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ; + if(! in ) continue ; + } + + //Create AOD for analysis + AliAODPWG4Particle aodph = AliAODPWG4Particle(mom); + Int_t label = calo->GetLabel(); + aodph.SetLabel(label); + aodph.SetInputFileIndex(input); + + //Set the indices of the original caloclusters + aodph.SetCaloLabel(calo->GetID(),-1); + aodph.SetDetector(fCalorimeter); + if(GetDebug() > 1) + printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",aodph.Pt(),aodph.Phi(),aodph.Eta()); + + //Check Distance to Bad channel, set bit. + Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel + if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ; + if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm) + continue ; + + if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad); + + if (distBad > fMinDist3) aodph.SetDistToBad(2) ; + else if(distBad > fMinDist2) aodph.SetDistToBad(1) ; + else aodph.SetDistToBad(0) ; + + //Skip matched clusters with tracks + if(fRejectTrackMatch && IsTrackMatched(calo)) continue ; + if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - TrackMatching cut passed \n"); + + //Set PID bits for later selection (AliAnaPi0 for example) + //GetPDG already called in SetPIDBits. + GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils()); + if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - PID Bits set \n"); + + //Play with the MC stack if available + //Check origin of the candidates + if(IsDataMC()){ + aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex())); + if(GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag()); + } + + // Check if cluster comes from a conversion in the material in front of the calorimeter + // Do invariant mass of all pairs, if mass is close to 0, then it is conversion. + + if(fCheckConversion && nCaloClusters > 1){ + Bool_t bConverted = kFALSE; + Int_t id2 = -1; + + //Check if set previously as converted couple, if so skip its use. + if (indexConverted[icalo]) continue; + + for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) { + //Check if set previously as converted couple, if so skip its use. + if (indexConverted[jcalo]) continue; + //printf("Check Conversion indeces %d and %d\n",icalo,jcalo); + AliAODCaloCluster * calo2 = (AliAODCaloCluster*) (pl->At(jcalo)); //Get cluster kinematics + Int_t evtIndex2 = 0 ; + if (GetMixedEvent()) { + evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ; + } + calo2->GetMomentum(mom2,GetVertex(evtIndex2)); + //Check only certain regions + Bool_t in2 = kTRUE; + if(IsFiducialCutOn()) in2 = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ; + if(!in2) continue; + + //Get mass of pair, if small, take this pair. + //printf("\t both in calo, mass %f, cut %f\n",(mom+mom2).M(),fMassCut); + if((mom+mom2).M() < fMassCut){ + bConverted = kTRUE; + id2 = calo2->GetID(); + indexConverted[jcalo]=kTRUE; + break; + } + + }//Mass loop + + if(bConverted){ + if(fAddConvertedPairsToAOD){ + //Create AOD of pair analysis + TLorentzVector mpair = mom+mom2; + AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair); + aodpair.SetLabel(aodph.GetLabel()); + aodpair.SetInputFileIndex(input); + + //printf("Index %d, Id %d\n",icalo, calo->GetID()); + //Set the indeces of the original caloclusters + aodpair.SetCaloLabel(calo->GetID(),id2); + aodpair.SetDetector(fCalorimeter); + aodpair.SetPdg(aodph.GetPdg()); + aodpair.SetTag(aodph.GetTag()); + + //Add AOD with pair object to aod branch + AddAODParticle(aodpair); + //printf("\t \t both added pair\n"); + } + + //Do not add the current calocluster + continue; + }//converted pair + }//check conversion + + //Add AOD with photon object to aod branch + AddAODParticle(aodph); + + }//loop; + delete [] indexConverted; + + if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast()); + +*/ +} +//__________________________________________________________________ +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; + //TParticle * primary = 0x0; + TClonesArray * mcparticles0 = 0x0; + //AliAODMCParticle * aodprimary = 0x0; + //TObjArray * pl = 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()){ + + if(GetReader()->ReadStack() && !stack) return; + if(GetReader()->ReadAODMCParticles() && !mcparticles0) return; + + //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()){ + 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") ; + +} diff --git a/PWG4/PartCorrDep/AliAnaShowerParameter.h b/PWG4/PartCorrDep/AliAnaShowerParameter.h index 7a2604d86b7..ecb15da2942 100644 --- a/PWG4/PartCorrDep/AliAnaShowerParameter.h +++ b/PWG4/PartCorrDep/AliAnaShowerParameter.h @@ -1,198 +1,93 @@ -#ifndef ALIANASHOWERPARAMETER_H -#define ALIANASHOWERPARAMETER_H -/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * - * See cxx source for full Copyright notice */ -/* $Id: AliAnaShowerParameter.h 27413 2008-07-18 13:28:12Z gconesab $ */ - -//_________________________________________________________________________ -// -// Class cloned from AliAnaPhoton, main aim is shower shape studies -// -// -// -//-- Author: Jocelyn Mlynarz (WSU) and Gustavo Conesa (LPSC) - -// --- ROOT system --- -class TH3F; -class TH2F ; -class TH1F; -class TString ; - -// --- ANALYSIS system --- -#include "AliAnaPartCorrBaseClass.h" -//#include "AliStack.h" -//#include "TParticle.h" -class AliStack; -class TParticle; - -class TList ; -class AliEMCALGeoUtils; -class AliAnaShowerParameter : public AliAnaPartCorrBaseClass { - -public: - - AliAnaShowerParameter() ; // default ctor - AliAnaShowerParameter(const AliAnaShowerParameter & g) ; // cpy ctor - AliAnaShowerParameter & operator = (const AliAnaShowerParameter & g) ;//cpy assignment - virtual ~AliAnaShowerParameter() ; //virtual dtor - - TList * GetCreateOutputObjects(); - - void Init(); - - TObjString* GetAnalysisCuts(); - - void MakeAnalysisFillAOD() ; - - void MakeAnalysisFillHistograms() ; - - void Print(const Option_t * opt)const; - - TString GetCalorimeter() const {return fCalorimeter ; } - void SetCalorimeter(TString det) {fCalorimeter = det ; } - - Bool_t IsTrackMatchRejectionOn() const {return fRejectTrackMatch ; } - void SwitchOnTrackMatchRejection() {fRejectTrackMatch = kTRUE ; } - void SwitchOffTrackMatchRejection() {fRejectTrackMatch = kFALSE ; } - - Bool_t IsCheckConversionOn() const {return fCheckConversion ; } - void SwitchOnConversionChecker() {fCheckConversion = kTRUE ; } - void SwitchOffConversionChecker() {fCheckConversion = kFALSE ; } - - Bool_t AreConvertedPairsInAOD() const {return fAddConvertedPairsToAOD ; } - void SwitchOnAdditionConvertedPairsToAOD() {fAddConvertedPairsToAOD = kTRUE ; } - void SwitchOffAdditionConvertedPairsToAOD() {fAddConvertedPairsToAOD = kFALSE ; } - - void InitParameters(); - - void SetMinDistanceToBadChannel(Float_t m1, Float_t m2, Float_t m3) { - fMinDist = m1; - fMinDist2 = m2; - fMinDist3 = m3; - } - - Float_t GetMassCut() const {return fMassCut ; } - void SetMassCut(Float_t m) {fMassCut = m ; } - void SetNClusterCut(Int_t nc) {fNumClusters = nc;} - Int_t GetNClusterCut() {return fNumClusters;} - void SetTimeCut(Double_t min, Double_t max) {fTimeCutMin = min; fTimeCutMax = max;} - Double_t GetTimeCutMin() const {return fTimeCutMin;} - Double_t GetTimeCutMax() const {return fTimeCutMax;} - - private: - - TString fCalorimeter ; // Calorimeter where the gamma is searched; - Float_t fMinDist ; // Minimal distance to bad channel to accept cluster - Float_t fMinDist2; // Cuts on Minimal distance to study acceptance evaluation - Float_t fMinDist3; // One more cut on distance used for acceptance-efficiency study - Bool_t fRejectTrackMatch ; //If PID on, reject clusters which have an associated TPC track - Bool_t fCheckConversion; // Combine pairs of clusters with mass close to 0 - Bool_t fAddConvertedPairsToAOD; // Put Converted pairs in AOD - Float_t fMassCut; // Mass cut for the conversion pairs selection - Float_t fNCellsCut; - Double_t fTimeCutMin ; // Remove clusters/cells with time smaller than this value, in ns - Double_t fTimeCutMax ; // Remove clusters/cells with time larger than this value, in ns - AliEMCALGeoUtils * fEMCALGeo; - Int_t fNumClusters; //Cut that selects events with a specific number of clusters - - //Histograms - TH1F * fhNClusters ; //! Number of clusters per event. - TH2F * fhNCellCluster;//! Number of cells per cluster - TH1F * fhPtCluster ; //! Number of clusters vs transerse momentum - TH2F * fhPhiCluster ; //! Azimuthal angle of clusters vs transerse momentum - TH2F * fhEtaCluster ; //! Pseudorapidity of clusters vs transerse momentum - TH1F * fhDeltaPhiClusters ; //! Delta phi of cluster pairs - TH1F * fhDeltaEtaClusters ; //! Delta eta of cluster pairs - TH3F * fhLambdaCluster ; //! Shower parameters of clusters vs transerse momentum - TH2F * fhDispersionCluster ; //! Dispersion of the clusters - TH3F * fhELambdaCluster ; //! Shower parameters of clusters vs Energy - TH3F * fhELambdaCellCluster ; //! Shower parameters of clusters vs Energy - - TH2F * fhNCellPhoton;//! Number of cells per photon - TH3F * fhLambdaPhoton ; //! Shower parameters of photons vs transerse momentum - TH2F * fhDispersionPhoton ; //! Dispersion of the photons - TH3F * fhELambdaPhoton ; //! Shower parameters of photons vs Energy - TH3F * fhELambdaCellPhoton ; //! Shower parameters of photons vs Energy - - TH2F * fhNCellPi0;//! Number of cells per neutral pion - TH3F * fhLambdaPi0 ; //! Shower parameters of neutral pions vs transerse momentum - TH2F * fhDispersionPi0 ; //! Dispersion of the neutral pions - TH3F * fhELambdaPi0 ; //! Shower parameters of neutral pions vs Energy - TH3F * fhELambdaCellPi0 ; //! Shower parameters of neutral pions vs Energy - - TH2F * fhNCellChargedHadron;//! Number of cells per charged hadron - TH3F * fhLambdaChargedHadron ; //! Shower parameters of charged hadrons vs transerse momentum - TH2F * fhDispersionChargedHadron ; //! Dispersion of the charged hadrons - TH3F * fhELambdaChargedHadron ; //! Shower parameters of charged hadrons vs Energy - TH3F * fhELambdaCellChargedHadron ; //! Shower parameters of charged hadrons vs Energy - - //MC - TH1F * fhDeltaE ; //! MC-Reco E distribution - TH1F * fhDeltaPt ; //! MC-Reco pT distribution - TH1F * fhRatioE ; //! Reco/MC E distribution - TH1F * fhRatioPt ; //! Reco/MC pT distribution - TH2F * fh2E ; //! E distribution, Reco vs MC - TH2F * fh2Pt ; //! pT distribution, Reco vs MC - TH2F * fhMCPdg; //! Complete list of PDG Codes. - - TH1F * fhEMCPhoton; //! Number of identified gamma - TH2F * fhPhiMCPhoton; //! Phi of identified gamma - TH2F * fhEtaMCPhoton; //! eta of identified gamma - TH3F * fhLambdaMCPhoton ; //! Shower parameters of MC photons vs transerse momentum - TH3F * fhPhotTrueE; // MC truth E vs Recons E vs. Lambda of the cluster for MC photons - TH1F * fhRatioEPhoton ; //! Reco/MC E distribution for photons - - TH1F * fhEMCPi0; //! Number of identified (single shower) Pi0 - TH2F * fhPhiMCPi0; //! Phi of identified Pi0 - TH2F * fhEtaMCPi0; //! eta of identified Pi0 - TH3F * fhLambdaMCPi0 ; //! Shower parameters of MC Pi0 vs transerse momentum - TH3F * fhPi0TrueE; // MC truth E vs Recons E vs. Lambda of the cluster for MC Pi0s - TH2F * fhProductionDistance; //! Distance from beam to production of the Pi0 - TH2F * fhProductionRadius; //! R from beam to Pi0 - TH2F * fhDecayAngle;//! Decay angle of the Pi0 - TH1F * fhRatioEPi0 ; //! Reco/MC E distribution for Pi0s - - TH1F * fhEMCPion; //! Number of identified pions - TH2F * fhPhiMCPion; //! Phi of identified pions - TH2F * fhEtaMCPion; //! eta of identified pions - TH3F * fhLambdaMCPion; //! Shower parameters of MC pions vs transerse momentum - TH3F * fhPionTrueE; // MC truth E vs Recons E vs. Lambda of the cluster for MC pions - TH1F * fhRatioEPion ; //! Reco/MC E distribution for charged pions - - TH1F * fhEMCProton; //! Number of identified protons - TH2F * fhPhiMCProton; //! Phi of identified protons - TH2F * fhEtaMCProton; //! eta of identified protons - TH3F * fhLambdaMCProton ; //! Shower parameters of MC protons vs transerse momentum - TH3F * fhProtonTrueE; // MC truth E vs Recons E vs. Lambda of the cluster for MC protons - TH1F * fhRatioEProton ; //! Reco/MC E distribution for protons - - TH1F * fhEMCAntiProton; //! Number of identified antiprotons - TH2F * fhPhiMCAntiProton; //! Phi of identified antiprotons - TH2F * fhEtaMCAntiProton; //! eta of identified antiprotons - TH3F * fhLambdaMCAntiProton ; //! Shower parameters of MC antiprotons vs transerse momentum - TH3F * fhAntiProtonTrueE; // MC truth E vs Recons E vs. Lambda of the cluster for MC antiprotons - TH1F * fhRatioEAntiProton ; //! Reco/MC E distribution for antiprotons - - TH1F * fhEMCNeutron; //! Number of identified neutrons - TH2F * fhPhiMCNeutron; //! Phi of identified neutrons - TH2F * fhEtaMCNeutron; //! eta of identified neutrons - TH3F * fhLambdaMCNeutron ; //! Shower parameters of MC neutrons vs transerse momentum - TH3F * fhNeutronTrueE; // MC truth E vs Recons E vs. Lambda of the cluster for MC Neutrons - TH1F * fhRatioENeutron ; //! Reco/MC E distribution for photons - - TH1F * fhEMCEta; //! Number of identified etas - TH2F * fhPhiMCEta; //! Phi of identified etas - TH2F * fhEtaMCEta; //! eta of identified etas - TH3F * fhLambdaMCEta ; //! Shower parameters of MC etas vs transerse momentum - - TH1D * fhPrimPt; //! Pi0 pT spectrum truth. - - ClassDef(AliAnaShowerParameter,1) - -} ; - - -#endif//AliAnaShowerParameter_H - - - +#ifndef ALIANASHOWERPARAMETER_H +#define ALIANASHOWERPARAMETER_H +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ +/* $Id: AliAnaShowerParameter.h 27413 2008-07-18 13:28:12Z gconesab $ */ + +//_________________________________________________________________________ +// +// Class cloned from AliAnaPhoton, main aim is shower shape studies +// +// +// +//-- Author: Jocelyn Mlynarz (WSU) and Gustavo Conesa (LPSC) + +// --- ROOT system --- +class TH3F; +class TH2F ; +class TH1F; +class TString ; + +// --- ANALYSIS system --- +#include "AliAnaPartCorrBaseClass.h" +//#include "AliStack.h" +//#include "TParticle.h" +class AliStack; +class TParticle; + +class TList ; +class AliEMCALGeoUtils; +class AliAnaShowerParameter : public AliAnaPartCorrBaseClass { + +public: + + AliAnaShowerParameter() ; // default ctor + AliAnaShowerParameter(const AliAnaShowerParameter & g) ; // cpy ctor + AliAnaShowerParameter & operator = (const AliAnaShowerParameter & g) ;//cpy assignment + virtual ~AliAnaShowerParameter() ; //virtual dtor + + TList * GetCreateOutputObjects(); + + void Init(); + + TObjString* GetAnalysisCuts(); + + void MakeAnalysisFillAOD() ; + + void MakeAnalysisFillHistograms() ; + + void Print(const Option_t * opt)const; + + TString GetCalorimeter() const {return fCalorimeter ; } + void SetCalorimeter(TString det) {fCalorimeter = det ; } + + void InitParameters(); + + void SetTimeCut(Double_t min, Double_t max) {fTimeCutMin = min; fTimeCutMax = max;} + Double_t GetTimeCutMin() const {return fTimeCutMin;} + Double_t GetTimeCutMax() const {return fTimeCutMax;} + + void SetNCellsCut(Double_t min, Double_t max) {fNCellsCutMin = min; fNCellsCutMax = max;} + Double_t GetNCellsCutMin() const {return fNCellsCutMin;} + Double_t GetNCellsCutMax() const {return fNCellsCutMax;} + + private: + + TString fCalorimeter ; // Calorimeter where the gamma is searched; + Float_t fNCellsCutMin ; + Float_t fNCellsCutMax ; + Float_t fLambdaCut ; + Double_t fTimeCutMin ; // Remove clusters/cells with time smaller than this value, in ns + Double_t fTimeCutMax ; // Remove clusters/cells with time larger than this value, in ns + + //Histograms + TH1F * fhNClusters ; + TH2F * fhNCellCluster; + TH3F * fhEtaPhiPtCluster ; + TH2F * fhLambdaPtCluster ; + + //MC + TH2F * fhLambdaPtPhoton ; + TH2F * fhLambdaPtPi0 ; + TH2F * fhLambdaPtPion ; + TH1D * fhPtTruthPi0 ; + + ClassDef(AliAnaShowerParameter,1) + +} ; + + +#endif//AliAnaShowerParameter_H + + +