-/**************************************************************************
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
- * *
- * Author: The ALICE Off-line Project. *
- * Contributors are mentioned in the code where appropriate. *
- * *
- * Permission to use, copy, modify and distribute this software and its *
- * documentation strictly for non-commercial purposes hereby granted *
- * without fee, provided that the above copyright notice appears in all *
- * copies and that both the copyright notice and this permission notice *
- * appear in the supporting documentation. The authors make no claims *
- * about the suitability of this software for any purpose. It is *
- * provided "as is" without express or implied warranty. *
- **************************************************************************/
-/* $Id: AliAnaPhoton.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
-
-//_________________________________________________________________________
-//
-// Class cloned from AliAnaPhoton, main aim is shower shape studies
-//
-//
-//
-//-- Author: Jocelyn Mlynarz (WSU) and Gustavo Conesa (LPSC-CNRS)
-//////////////////////////////////////////////////////////////////////////////
-
-
-// --- ROOT system ---
-#include <TH3F.h>
-#include <TH2F.h>
-#include <TClonesArray.h>
-//#include <TObjString.h>
-#include <Riostream.h>
-#include "TParticle.h"
-//#include <fstream>
-
-// --- Analysis system ---
-#include "AliAnaShowerParameter.h"
-#include "AliCaloTrackReader.h"
-#include "AliStack.h"
-#include "AliCaloPID.h"
-#include "AliMCAnalysisUtils.h"
-#include "AliFiducialCut.h"
-#include "AliAODCaloCluster.h"
-#include "AliAODMCParticle.h"
-#include "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;jaod<naod;jaod++){
- AliAODPWG4Particle* phSecond = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(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 ; i<stack->GetNprimary(); 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") ;
-
-}
+/**************************************************************************\r
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
+ * *\r
+ * Author: The ALICE Off-line Project. *\r
+ * Contributors are mentioned in the code where appropriate. *\r
+ * *\r
+ * Permission to use, copy, modify and distribute this software and its *\r
+ * documentation strictly for non-commercial purposes hereby granted *\r
+ * without fee, provided that the above copyright notice appears in all *\r
+ * copies and that both the copyright notice and this permission notice *\r
+ * appear in the supporting documentation. The authors make no claims *\r
+ * about the suitability of this software for any purpose. It is *\r
+ * provided "as is" without express or implied warranty. *\r
+ **************************************************************************/\r
+/* $Id: AliAnaPhoton.cxx 28688 2008-09-11 15:04:07Z gconesab $ */\r
+\r
+//_________________________________________________________________________\r
+//\r
+// Class cloned from AliAnaPhoton, main aim is shower shape studies\r
+// \r
+// \r
+//\r
+//-- Author: Jocelyn Mlynarz (WSU) and Gustavo Conesa (LPSC-CNRS)\r
+//////////////////////////////////////////////////////////////////////////////\r
+\r
+\r
+// --- ROOT system --- \r
+#include <TH3F.h>\r
+#include <TH2F.h>\r
+#include <TClonesArray.h>\r
+//#include <TObjString.h>\r
+#include <Riostream.h>\r
+#include "TParticle.h"\r
+//#include <fstream>\r
+\r
+// --- Analysis system --- \r
+#include "AliAnaShowerParameter.h" \r
+#include "AliCaloTrackReader.h"\r
+#include "AliStack.h"\r
+#include "AliCaloPID.h"\r
+#include "AliMCAnalysisUtils.h"\r
+#include "AliFiducialCut.h"\r
+#include "AliAODCaloCluster.h"\r
+#include "AliAODMCParticle.h"\r
+#include "AliAnalysisManager.h"\r
+#include "AliAnalysisTaskParticleCorrelation.h"\r
+#include "AliEMCALGeoUtils.h"\r
+#include "AliAODEvent.h"\r
+\r
+\r
+ClassImp(AliAnaShowerParameter)\r
+\r
+//____________________________________________________________________________\r
+AliAnaShowerParameter::AliAnaShowerParameter() : \r
+AliAnaPartCorrBaseClass(), fCalorimeter(""), fNCellsCutMin(0),\r
+fNCellsCutMax(0), fLambdaCut(0), fTimeCutMin(-1), fTimeCutMax(9999999),\r
+fhNClusters(0), fhNCellCluster(0), fhEtaPhiPtCluster(0), \r
+fhLambdaPtCluster(0), \r
+\r
+//MC\r
+\r
+fhLambdaPtPhoton(0), fhLambdaPtPi0(0), fhLambdaPtPion(0), fhPtTruthPi0(0)\r
+\r
+{\r
+ //default ctor\r
+ \r
+ //Initialize parameters\r
+ InitParameters();\r
+ \r
+}\r
+\r
+//____________________________________________________________________________\r
+AliAnaShowerParameter::~AliAnaShowerParameter() \r
+{\r
+ //dtor\r
+ \r
+}\r
+\r
+//________________________________________________________________________\r
+TObjString * AliAnaShowerParameter::GetAnalysisCuts()\r
+{ \r
+ //Save parameters used for analysis\r
+ TString parList ; //this will be list of parameters used for this analysis.\r
+ const Int_t buffersize = 255;\r
+ char onePar[buffersize] ;\r
+ \r
+ snprintf(onePar,buffersize,"--- AliAnaShowerParameter ---\n") ;\r
+ parList+=onePar ; \r
+ snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"fNCellsCutMin: (Cut in the minimum number of cells) %e\n",fNCellsCutMin) ;\r
+ parList+=onePar ; \r
+ snprintf(onePar,buffersize,"fNCellsCutMax: (Cut in the maximum number of cells) %e\n",fNCellsCutMax) ;\r
+ parList+=onePar ;\r
+ snprintf(onePar,buffersize,"fLambdaCut: (Cut in the minimum lambda) %e\n",fLambdaCut) ;\r
+ parList+=onePar ;\r
+ \r
+ //Get parameters set in base class.\r
+ parList += GetBaseParametersList() ;\r
+ \r
+ //Get parameters set in PID class.\r
+ parList += GetCaloPID()->GetPIDParametersList() ;\r
+ \r
+ //Get parameters set in FiducialCut class (not available yet)\r
+ //parlist += GetFidCut()->GetFidCutParametersList() \r
+ \r
+ return new TObjString(parList) ;\r
+}\r
+\r
+\r
+//________________________________________________________________________\r
+TList * AliAnaShowerParameter::GetCreateOutputObjects()\r
+{ \r
+ // Create histograms to be saved in output file and \r
+ // store them in outputContainer\r
+ TList * outputContainer = new TList() ; \r
+ outputContainer->SetName("PhotonHistos") ; \r
+ \r
+ Int_t nptbins = GetHistoPtBins();\r
+ Int_t nphibins = GetHistoPhiBins();\r
+ Int_t netabins = GetHistoEtaBins();\r
+ Float_t ptmax = GetHistoPtMax();\r
+ Float_t phimax = GetHistoPhiMax();\r
+ Float_t etamax = GetHistoEtaMax();\r
+ Float_t ptmin = GetHistoPtMin();\r
+ Float_t phimin = GetHistoPhiMin();\r
+ Float_t etamin = GetHistoEtaMin(); \r
+ \r
+ //General non-MC Cluster histograms\r
+ fhNClusters = new TH1F ("hNClusters","NClusters",21,-0.5,20.5); \r
+ fhNClusters->SetXTitle("N_{Clusters}");\r
+ outputContainer->Add(fhNClusters) ;\r
+ \r
+ fhNCellCluster = new TH2F ("hNCellCluster","Number of cell per cluster",nptbins,ptmin,ptmax,21,-0.5,20.5); \r
+ fhNCellCluster->SetYTitle("N_{Cells}"); \r
+ fhNCellCluster->SetXTitle("p_{T}");\r
+ outputContainer->Add(fhNCellCluster) ;\r
+ \r
+ fhEtaPhiPtCluster = new TH3F\r
+ ("hEtaPhiPtCluster","#phi_{Cluster} and #eta_{Cluster}",nptbins,ptmin,ptmax,nphibins,phimin,phimax,netabins,etamin,etamax); \r
+ fhEtaPhiPtCluster->SetZTitle("#eta"); \r
+ fhEtaPhiPtCluster->SetYTitle("#phi");\r
+ fhEtaPhiPtCluster->SetXTitle("pT_{Cluster} (GeV/c)");\r
+ outputContainer->Add(fhEtaPhiPtCluster) ; \r
+ \r
+ fhLambdaPtCluster = new TH2F\r
+ ("hLambdaCluster","#lambda_{Cluster}",nptbins,ptmin,ptmax,300,0,3); \r
+ fhLambdaPtCluster->SetYTitle("#lambda_{0}^{2}");\r
+ fhLambdaPtCluster->SetXTitle("p_{T Cluster} (GeV/c)");\r
+ outputContainer->Add(fhLambdaPtCluster) ;\r
+ \r
+ if(IsDataMC()){\r
+ \r
+ fhLambdaPtPhoton = new TH2F\r
+ ("hLambdaPtPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2); \r
+ fhLambdaPtPhoton->SetYTitle("#lambda_{0}^{2}");\r
+ fhLambdaPtPhoton->SetXTitle("pT_{#gamma, Reco} (GeV)");\r
+ outputContainer->Add(fhLambdaPtPhoton) ;\r
+ \r
+ fhLambdaPtPi0 = new TH2F\r
+ ("hLambdaPtPi0","#lambda_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,2); \r
+ fhLambdaPtPi0->SetYTitle("#lambda_{0}^{2}");\r
+ fhLambdaPtPi0->SetXTitle("pT_{#pi^{0}, Reco} (GeV)");\r
+ outputContainer->Add(fhLambdaPtPi0) ;\r
+ \r
+ fhLambdaPtPion = new TH2F\r
+ ("hLambdaPtPion","#lambda_{#pi^{+}}",nptbins,ptmin,ptmax,200,0,2); \r
+ fhLambdaPtPion->SetYTitle("#lambda_{0}^{2}");\r
+ fhLambdaPtPion->SetXTitle("pT_{#pi^{+}, Reco} (GeV)");\r
+ outputContainer->Add(fhLambdaPtPion) ;\r
+ \r
+ fhPtTruthPi0 = new TH1D("hPtTruthPi0","#pi^{0} MC truth pT",nptbins,ptmin,ptmax) ;\r
+ fhPtTruthPi0->SetXTitle("pT_{#pi^{0}, Reco} (GeV)");\r
+ outputContainer->Add(fhPtTruthPi0) ;\r
+ \r
+ }//Histos with MC\r
+ \r
+ return outputContainer ;\r
+ \r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliAnaShowerParameter::Init()\r
+{\r
+ \r
+ //Init\r
+ //Do some checks\r
+ if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){\r
+ printf("AliAnaShowerParameter::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");\r
+ abort();\r
+ }\r
+ else if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){\r
+ printf("AliAnaShowerParameter::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");\r
+ abort();\r
+ }\r
+ \r
+}\r
+\r
+//____________________________________________________________________________\r
+void AliAnaShowerParameter::InitParameters()\r
+{\r
+ \r
+ //Initialize the parameters of the analysis.\r
+ AddToHistogramsName("AnaLambda_");\r
+ \r
+ fCalorimeter = "EMCAL" ;\r
+ \r
+ fTimeCutMin = -1;\r
+ fTimeCutMax = 9999999;\r
+ fNCellsCutMin = 0 ;\r
+ fNCellsCutMax = 0 ;\r
+ fLambdaCut = 0.01 ;\r
+ \r
+}\r
+\r
+//__________________________________________________________________\r
+void AliAnaShowerParameter::MakeAnalysisFillAOD() \r
+{\r
+ /*//Do analysis and fill aods\r
+ //Search for photons in fCalorimeter \r
+ \r
+ //Get vertex for photon momentum calculation\r
+ \r
+ for (Int_t iev = 0; iev < GetNMixedEvent(); iev++) {\r
+ if (!GetMixedEvent()) \r
+ GetReader()->GetVertex(GetVertex(iev));\r
+ else \r
+ GetMixedEvent()->GetVertexOfEvent(iev)->GetXYZ(GetVertex(iev)); \r
+ } \r
+ \r
+ //Select the Calorimeter of the photon\r
+ TObjArray * pl = 0x0; \r
+ if(fCalorimeter == "PHOS")\r
+ pl = GetPHOSClusters();\r
+ else if (fCalorimeter == "EMCAL")\r
+ pl = GetEMCALClusters();\r
+ \r
+ if(!pl){\r
+ printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Careful cluster array NULL!!\n");\r
+ return;\r
+ }\r
+ \r
+ //Fill AODCaloClusters and AODParticle with PHOS/EMCAL aods\r
+ TLorentzVector mom, mom2 ; \r
+ Int_t nCaloClusters = pl->GetEntriesFast(); \r
+ //Cut on the number of clusters in the event.\r
+ if ((fNumClusters !=-1) && (nCaloClusters != fNumClusters)) return;\r
+ Bool_t * indexConverted = new Bool_t[nCaloClusters];\r
+ for (Int_t i = 0; i < nCaloClusters; i++) \r
+ indexConverted[i] = kFALSE;\r
+ \r
+ for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){ \r
+ \r
+ AliAODCaloCluster * calo = (AliAODCaloCluster*) (pl->At(icalo)); \r
+ Int_t evtIndex = 0 ; \r
+ if (GetMixedEvent()) {\r
+ evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ; \r
+ }\r
+ //Cluster selection, not charged, with photon id and in fiducial cut\r
+ \r
+ //Input from second AOD?\r
+ Int_t input = 0;\r
+ \r
+ //Get Momentum vector, \r
+ if (input == 0) \r
+ calo->GetMomentum(mom,GetVertex(evtIndex)) ;//Assume that come from vertex in straight line\r
+ \r
+ //Skip the cluster if it doesn't fit inside the cuts.\r
+ if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ; \r
+ Double_t tof = calo->GetTOF()*1e9; \r
+ if(tof < fTimeCutMin || tof > fTimeCutMax) continue; \r
+ if(calo->GetNCells() <= fNCellsCut) continue;\r
+ \r
+ //Check acceptance selection\r
+ if(IsFiducialCutOn()){\r
+ Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;\r
+ if(! in ) continue ;\r
+ }\r
+ \r
+ //Create AOD for analysis\r
+ AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);\r
+ Int_t label = calo->GetLabel();\r
+ aodph.SetLabel(label);\r
+ aodph.SetInputFileIndex(input);\r
+ \r
+ //Set the indices of the original caloclusters \r
+ aodph.SetCaloLabel(calo->GetID(),-1);\r
+ aodph.SetDetector(fCalorimeter);\r
+ if(GetDebug() > 1) \r
+ 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()); \r
+ \r
+ //Check Distance to Bad channel, set bit.\r
+ Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel\r
+ if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;\r
+ if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)\r
+ continue ;\r
+ \r
+ if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad);\r
+ \r
+ if (distBad > fMinDist3) aodph.SetDistToBad(2) ;\r
+ else if(distBad > fMinDist2) aodph.SetDistToBad(1) ; \r
+ else aodph.SetDistToBad(0) ;\r
+ \r
+ //Skip matched clusters with tracks\r
+ if(fRejectTrackMatch && IsTrackMatched(calo)) continue ;\r
+ if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - TrackMatching cut passed \n");\r
+ \r
+ //Set PID bits for later selection (AliAnaPi0 for example)\r
+ //GetPDG already called in SetPIDBits.\r
+ GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils());\r
+ if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - PID Bits set \n"); \r
+ \r
+ //Play with the MC stack if available\r
+ //Check origin of the candidates\r
+ if(IsDataMC()){\r
+ aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));\r
+ if(GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());\r
+ }\r
+ \r
+ // Check if cluster comes from a conversion in the material in front of the calorimeter\r
+ // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.\r
+ \r
+ if(fCheckConversion && nCaloClusters > 1){\r
+ Bool_t bConverted = kFALSE;\r
+ Int_t id2 = -1;\r
+ \r
+ //Check if set previously as converted couple, if so skip its use.\r
+ if (indexConverted[icalo]) continue;\r
+ \r
+ for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {\r
+ //Check if set previously as converted couple, if so skip its use.\r
+ if (indexConverted[jcalo]) continue;\r
+ //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);\r
+ AliAODCaloCluster * calo2 = (AliAODCaloCluster*) (pl->At(jcalo)); //Get cluster kinematics\r
+ Int_t evtIndex2 = 0 ; \r
+ if (GetMixedEvent()) {\r
+ evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ; \r
+ } \r
+ calo2->GetMomentum(mom2,GetVertex(evtIndex2));\r
+ //Check only certain regions\r
+ Bool_t in2 = kTRUE;\r
+ if(IsFiducialCutOn()) in2 = GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;\r
+ if(!in2) continue; \r
+ \r
+ //Get mass of pair, if small, take this pair.\r
+ //printf("\t both in calo, mass %f, cut %f\n",(mom+mom2).M(),fMassCut);\r
+ if((mom+mom2).M() < fMassCut){ \r
+ bConverted = kTRUE;\r
+ id2 = calo2->GetID();\r
+ indexConverted[jcalo]=kTRUE;\r
+ break;\r
+ }\r
+ \r
+ }//Mass loop\r
+ \r
+ if(bConverted){ \r
+ if(fAddConvertedPairsToAOD){\r
+ //Create AOD of pair analysis\r
+ TLorentzVector mpair = mom+mom2;\r
+ AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);\r
+ aodpair.SetLabel(aodph.GetLabel());\r
+ aodpair.SetInputFileIndex(input);\r
+ \r
+ //printf("Index %d, Id %d\n",icalo, calo->GetID());\r
+ //Set the indeces of the original caloclusters \r
+ aodpair.SetCaloLabel(calo->GetID(),id2);\r
+ aodpair.SetDetector(fCalorimeter);\r
+ aodpair.SetPdg(aodph.GetPdg());\r
+ aodpair.SetTag(aodph.GetTag());\r
+ \r
+ //Add AOD with pair object to aod branch\r
+ AddAODParticle(aodpair);\r
+ //printf("\t \t both added pair\n");\r
+ }\r
+ \r
+ //Do not add the current calocluster\r
+ continue;\r
+ }//converted pair\r
+ }//check conversion\r
+ \r
+ //Add AOD with photon object to aod branch\r
+ AddAODParticle(aodph);\r
+ \r
+ }//loop;\r
+ delete [] indexConverted;\r
+ \r
+ if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast()); \r
+ \r
+*/\r
+}\r
+//__________________________________________________________________\r
+void AliAnaShowerParameter::MakeAnalysisFillHistograms() \r
+{\r
+ \r
+ //Do analysis and fill histograms\r
+ if ( GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Starting analysis.\n");\r
+ \r
+ // Access MC information in stack if requested, check that it exists. \r
+ AliStack * stack = 0x0;\r
+ //TParticle * primary = 0x0; \r
+ TClonesArray * mcparticles0 = 0x0;\r
+ //AliAODMCParticle * aodprimary = 0x0; \r
+ //TObjArray * pl = 0x0;\r
+ Int_t NClusters = 0 ;\r
+ TLorentzVector momCluster ;\r
+\r
+ \r
+ //Check if the stack is available when analysing MC data.\r
+ if(IsDataMC()){\r
+ \r
+ if(GetReader()->ReadStack()){\r
+ stack = GetMCStack() ;\r
+ if(!stack) {\r
+ printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");\r
+ abort();\r
+ }\r
+ \r
+ }\r
+ else if(GetReader()->ReadAODMCParticles()){\r
+ \r
+ //Get the list of MC particles\r
+ mcparticles0 = GetReader()->GetAODMCParticles(0);\r
+ if(!mcparticles0 && GetDebug() > 0) {\r
+ printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Standard MCParticles not available !\n"); \r
+ }\r
+ }// is data and MC\r
+ } \r
+ //Loop on stored AOD photons\r
+\r
+ TClonesArray * clustArray = 0x0 ; \r
+ clustArray = GetAODCaloClusters() ;\r
+ NClusters = clustArray->GetEntriesFast() ; \r
+ fhNClusters->Fill(NClusters) ;\r
+\r
+ if(GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - aod branch entries %d\n", NClusters);\r
+ \r
+ for(Int_t iCluster = 0 ; iCluster < NClusters ; iCluster++){ \r
+ Int_t input = 0 ;\r
+\r
+ AliAODCaloCluster * caloCluster = (AliAODCaloCluster*) (clustArray->At(iCluster)) ;\r
+ caloCluster->GetMomentum(momCluster,GetVertex(0)) ;\r
+ AliAODPWG4Particle * aodCluster = new AliAODPWG4Particle(momCluster) ;\r
+ Int_t labelCluster = caloCluster->GetLabel() ;\r
+ aodCluster->SetLabel(labelCluster) ;;\r
+ aodCluster->SetInputFileIndex(input) ;\r
+ aodCluster->SetCaloLabel(caloCluster->GetID(),-1) ;\r
+ aodCluster->SetDetector(fCalorimeter) ;\r
+ aodCluster->SetTag(GetMCAnalysisUtils()->CheckOrigin(caloCluster->GetLabels(),caloCluster->GetNLabels(),GetReader(), aodCluster->GetInputFileIndex())) ;\r
+\r
+ if (GetDebug() > 3) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Cluster particle label = %d\n",labelCluster) ;\r
+ \r
+ //Fill Cluster histograms \r
+ Float_t ptcluster = aodCluster->Pt() ;\r
+ Float_t phicluster = aodCluster->Phi() ;\r
+ Float_t etacluster = aodCluster->Eta() ;\r
+ Float_t lambdaMainCluster = caloCluster->GetM02() ;\r
+ Int_t iNumCell = caloCluster->GetNCells() ; \r
+ \r
+ if (iNumCell>=fNCellsCutMin&&iNumCell<=fNCellsCutMax&&lambdaMainCluster>=fLambdaCut){\r
+ \r
+ //Fill the basic non-MC information about the cluster.\r
+ fhNCellCluster->Fill(ptcluster,iNumCell) ;\r
+ fhEtaPhiPtCluster->Fill(ptcluster,phicluster,etacluster) ;\r
+ fhLambdaPtCluster->Fill(ptcluster,lambdaMainCluster) ;\r
+ \r
+ //Play with the MC data if available\r
+ if(IsDataMC()){\r
+ \r
+ if(GetReader()->ReadStack() && !stack) return;\r
+ if(GetReader()->ReadAODMCParticles() && !mcparticles0) return;\r
+ \r
+ //Get the tag from AliMCAnalysisUtils for PID\r
+ Int_t tag = aodCluster->GetTag();\r
+ if (GetDebug() > 3) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Cluster particle tag= %d\n",tag) ;\r
+ if ( aodCluster->GetLabel() < 0){\r
+ if(GetDebug() > 0) \r
+ printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Label is -1; problem in the MC ESD? ");\r
+ continue;\r
+ }\r
+ \r
+ //Check if the tag matches one of the different particle types and fill the corresponding histograms\r
+ if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)&&!(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))) {\r
+ fhLambdaPtPhoton->Fill(ptcluster,lambdaMainCluster) ;\r
+ }//kMCPhoton\r
+ \r
+ if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)) { \r
+ fhLambdaPtPi0->Fill(ptcluster,lambdaMainCluster) ; \r
+ }\r
+ if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion)) {\r
+ fhLambdaPtPion->Fill(ptcluster,lambdaMainCluster) ;\r
+ }\r
+ \r
+ // Access MC information in stack if requested, check that it exists.\r
+ Int_t label = aodCluster->GetLabel();\r
+ if(label < 0) {\r
+ printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** bad label ***: label %d \n", label);\r
+ continue;\r
+ }\r
+ }\r
+ }\r
+ }\r
+ \r
+ Float_t fPtGener = 0 ;\r
+ if(IsDataMC()){\r
+ TParticle* pGener = stack->Particle(0) ;\r
+ fPtGener = pGener->Pt() ;\r
+ }\r
+\r
+ \r
+}\r
+\r
+\r
+//__________________________________________________________________\r
+void AliAnaShowerParameter::Print(const Option_t * opt) const\r
+{\r
+ //Print some relevant parameters set for the analysis\r
+ \r
+ if(! opt)\r
+ return;\r
+ \r
+ printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;\r
+ AliAnaPartCorrBaseClass::Print(" ");\r
+ printf("Calorimeter = %s\n", fCalorimeter.Data()) ;\r
+ printf("Min number of cells in cluster is > %f \n", fNCellsCutMin);\r
+ printf("Max number of cells in cluster is > %f \n", fNCellsCutMax);\r
+ printf("Min lambda of cluster is > %f \n", fLambdaCut);\r
+ printf(" \n") ;\r
+ \r
+} \r
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff