PWG4/PartCorrDep/AliAnaShowerParameter.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes hereby granted      *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id: AliAnaPhoton.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
  16
  17 //_________________________________________________________________________
  18 //
  19 // Class cloned from AliAnaPhoton, main aim is shower shape studies
  20 //
  21 //
  22 //
  23 //-- Author: Jocelyn Mlynarz (WSU) and Gustavo Conesa (LPSC-CNRS)
  24 //////////////////////////////////////////////////////////////////////////////
  25
  26
  27 // --- ROOT system ---
  28 #include <TH3F.h>
  29 #include <TH2F.h>
  30 #include <TClonesArray.h>
  31 //#include <TObjString.h>
  32 #include <Riostream.h>
  33 #include "TParticle.h"
  34 //#include <fstream>
  35
  36 // --- Analysis system ---
  37 #include "AliAnaShowerParameter.h"
  38 #include "AliCaloTrackReader.h"
  39 #include "AliStack.h"
  40 #include "AliCaloPID.h"
  41 #include "AliMCAnalysisUtils.h"
  42 #include "AliFiducialCut.h"
  43 #include "AliAODCaloCluster.h"
  44 #include "AliAODMCParticle.h"
  45 #include "AliAnaShowerParameter.h"
  46 #include "AliEMCALGeoUtils.h"
  47 #include "AliAODCaloCells.h"
  48 #include "AliAODEvent.h"
  49
  50
  51 ClassImp(AliAnaShowerParameter)
  52
  53 //____________________________________________________________________________
  54 AliAnaShowerParameter::AliAnaShowerParameter() :
  55 AliAnaPartCorrBaseClass(), fCalorimeter(""),
  56 fMinDist(0.),fMinDist2(0.),fMinDist3(0.),fRejectTrackMatch(0),
  57 fCheckConversion(kFALSE),fAddConvertedPairsToAOD(kFALSE), fMassCut(0), fNCellsCut(0),
  58 fTimeCutMin(-1), fTimeCutMax(9999999),fEMCALGeo(0),fNumClusters(0),
  59 fhNClusters(0),fhNCellCluster(0),fhPtCluster(0),fhPhiCluster(0),fhEtaCluster(0),fhDeltaPhiClusters(0),fhDeltaEtaClusters(0),
  60 fhLambdaCluster(0),fhDispersionCluster(0),fhELambdaCluster(0),fhELambdaCellCluster(0),
  61 fhNCellPhoton(0),fhLambdaPhoton(0),fhDispersionPhoton(0),fhELambdaPhoton(0),fhELambdaCellPhoton(0),
  62 fhNCellPi0(0),fhLambdaPi0(0),fhDispersionPi0(0),fhELambdaPi0(0),fhELambdaCellPi0(0),
  63 fhNCellChargedHadron(0),fhLambdaChargedHadron(0),fhDispersionChargedHadron(0),fhELambdaChargedHadron(0),fhELambdaCellChargedHadron(0),
  64 //MC
  65 fhDeltaE(0), fhDeltaPt(0),fhRatioE(0), fhRatioPt(0),fh2E(0),fh2Pt(0),fhMCPdg(0),
  66 fhEMCPhoton(0),fhPhiMCPhoton(0),fhEtaMCPhoton(0),fhLambdaMCPhoton(0),fhPhotTrueE(0),fhRatioEPhoton(0),
  67 fhEMCPi0(0),fhPhiMCPi0(0),fhEtaMCPi0(0),fhLambdaMCPi0(0),fhPi0TrueE(0),
  68 fhProductionDistance(0),fhProductionRadius(0),fhDecayAngle(0),fhRatioEPi0(0),
  69 fhEMCPion(0),fhPhiMCPion(0),fhEtaMCPion(0),fhLambdaMCPion(0),fhPionTrueE(0),fhRatioEPion(0),
  70 fhEMCProton(0),fhPhiMCProton(0),fhEtaMCProton(0),fhLambdaMCProton(0),fhProtonTrueE(0),fhRatioEProton(0),
  71 fhEMCAntiProton(0),fhPhiMCAntiProton(0),fhEtaMCAntiProton(0),fhLambdaMCAntiProton(0),fhAntiProtonTrueE(0),fhRatioEAntiProton(0),
  72 fhEMCNeutron(0),fhPhiMCNeutron(0),fhEtaMCNeutron(0),fhLambdaMCNeutron(0),fhNeutronTrueE(0),fhRatioENeutron(0),
  73 fhEMCEta(0),fhPhiMCEta(0),fhEtaMCEta(0),fhLambdaMCEta(0),
  74 fhPrimPt(0x0)
  75 {
  76   //default ctor
  77
  78   //Initialize parameters
  79   InitParameters();
  80
  81 }
  82
  83 //____________________________________________________________________________
  84 AliAnaShowerParameter::~AliAnaShowerParameter()
  85 {
  86   //dtor
  87
  88 }
  89
  90 //________________________________________________________________________
  91 TObjString *  AliAnaShowerParameter::GetAnalysisCuts()
  92 {
  93   //Save parameters used for analysis
  94   TString parList ; //this will be list of parameters used for this analysis.
  95   const Int_t buffersize = 255;
  96   char onePar[buffersize] ;
  97
  98   snprintf(onePar,buffersize,"--- AliAnaShowerParameter ---\n") ;
  99   parList+=onePar ;
 100   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
 101   parList+=onePar ;
 102   snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
 103   parList+=onePar ;
 104   snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
 105   parList+=onePar ;
 106   snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
 107   parList+=onePar ;
 108   snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
 109   parList+=onePar ;
 110
 111   //Get parameters set in base class.
 112   parList += GetBaseParametersList() ;
 113
 114   //Get parameters set in PID class.
 115   parList += GetCaloPID()->GetPIDParametersList() ;
 116
 117   //Get parameters set in FiducialCut class (not available yet)
 118   //parlist += GetFidCut()->GetFidCutParametersList()
 119
 120   return new TObjString(parList) ;
 121 }
 122
 123
 124 //________________________________________________________________________
 125 TList *  AliAnaShowerParameter::GetCreateOutputObjects()
 126 {
 127   // Create histograms to be saved in output file and
 128   // store them in outputContainer
 129   TList * outputContainer = new TList() ;
 130   outputContainer->SetName("PhotonHistos") ;
 131
 132   Int_t nptbins  = GetHistoPtBins();
 133   Int_t nphibins = GetHistoPhiBins();
 134   Int_t netabins = GetHistoEtaBins();
 135   Float_t ptmax  = GetHistoPtMax();
 136   Float_t phimax = GetHistoPhiMax();
 137   Float_t etamax = GetHistoEtaMax();
 138   Float_t ptmin  = GetHistoPtMin();
 139   Float_t phimin = GetHistoPhiMin();
 140   Float_t etamin = GetHistoEtaMin();
 141
 142   //General non-MC Cluster histograms
 143   fhNClusters = new TH1F ("hNClusters","NClusters",21,-0.5,20.5);
 144   fhNClusters->SetXTitle("N_{Clusters}");
 145   outputContainer->Add(fhNClusters) ;
 146
 147   fhNCellCluster = new TH2F ("hNCellCluster","Number of cell per cluster",nptbins,ptmin,ptmax,21,-0.5,20.5);
 148   fhNCellCluster->SetXTitle("p_{T}");
 149   fhNCellCluster->SetYTitle("N_{Cells}");
 150   outputContainer->Add(fhNCellCluster) ;
 151
 152   fhPtCluster  = new TH1F("hPtCluster","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 153   fhPtCluster->SetXTitle("p_{Cluster}(GeV/c)");
 154   outputContainer->Add(fhPtCluster) ;
 155
 156   fhPhiCluster  = new TH2F
 157   ("hPhiCluster","#phi_{Cluster}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 158   fhPhiCluster->SetYTitle("#phi");
 159   fhPhiCluster->SetXTitle("E_{Cluster} (GeV/c)");
 160   outputContainer->Add(fhPhiCluster) ;
 161
 162   fhEtaCluster  = new TH2F
 163   ("hEtaCluster","#phi_{Cluster}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 164   fhEtaCluster->SetYTitle("#eta");
 165   fhEtaCluster->SetXTitle("E_{Cluster} (GeV/c)");
 166   outputContainer->Add(fhEtaCluster) ;
 167
 168   fhDeltaEtaClusters  = new TH1F("hDeltaEtaClusters","#Delta #eta of clusters over calorimeter",netabins,etamin,etamax);
 169   fhDeltaEtaClusters->SetXTitle("#Delta #eta_{Clusters}");
 170   outputContainer->Add(fhDeltaEtaClusters) ;
 171
 172   fhDeltaPhiClusters  = new TH1F("hDeltaPhiClusters","#Delta #phi of clusters over calorimeter",nphibins,phimin,phimax);
 173   fhDeltaPhiClusters->SetXTitle("#Delta #phi_{Clusters}");
 174   outputContainer->Add(fhDeltaPhiClusters) ;
 175
 176   fhLambdaCluster  = new TH3F
 177   ("hLambdaCluster","#lambda_{Cluster}",nptbins,ptmin,ptmax,200,0,2,200,0,2);
 178   fhLambdaCluster->SetZTitle("#lambda_{1}^{2}");
 179   fhLambdaCluster->SetYTitle("#lambda_{0}^{2}");
 180   fhLambdaCluster->SetXTitle("p_{T Cluster} (GeV/c)");
 181   outputContainer->Add(fhLambdaCluster) ;
 182
 183   fhELambdaCluster  = new TH3F
 184   ("hELambdaCluster","#lambda_{Cluster}",nptbins,ptmin,ptmax,200,0,2,200,0,2);
 185   fhELambdaCluster->SetZTitle("#lambda_{1}^{2}");
 186   fhELambdaCluster->SetYTitle("#lambda_{0}^{2}");
 187   fhELambdaCluster->SetXTitle("p_{T Cluster} (GeV)");
 188   outputContainer->Add(fhELambdaCluster) ;
 189
 190   fhDispersionCluster  = new TH2F
 191   ("hDispersionCluster","Dispersion_{Cluster}",nptbins,ptmin,ptmax,200,0,2);
 192   fhDispersionCluster->SetYTitle("Dispersion");
 193   fhDispersionCluster->SetXTitle("p_{T Cluster} (GeV/c)");
 194   outputContainer->Add(fhDispersionCluster) ;
 195
 196   fhELambdaCellCluster  = new TH3F
 197   ("hELambdaCellCluster","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5);
 198   fhELambdaCellCluster->SetZTitle("N Cells");
 199   fhELambdaCellCluster->SetYTitle("#lambda_{0}^{2}");
 200   fhELambdaCellCluster->SetXTitle("E_{Cluster} (GeV)");
 201   outputContainer->Add(fhELambdaCellCluster) ;
 202
 203   //Identified photon histograms
 204   fhNCellPhoton = new TH2F ("hNCellPhoton","Number of cell per cluster",nptbins,ptmin,ptmax,21,-0.5,20.5);
 205   fhNCellPhoton->SetXTitle("E");
 206   fhNCellPhoton->SetYTitle("N_{Cells}");
 207   outputContainer->Add(fhNCellPhoton) ;
 208
 209   fhLambdaPhoton  = new TH3F
 210   ("hLambdaPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2);
 211   fhLambdaPhoton->SetZTitle("#lambda_{1}^{2}");
 212   fhLambdaPhoton->SetYTitle("#lambda_{0}^{2}");
 213   fhLambdaPhoton->SetXTitle("p_{T Photon} (GeV/c)");
 214   outputContainer->Add(fhLambdaPhoton) ;
 215
 216   fhELambdaPhoton  = new TH3F
 217   ("hELambdaPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2);
 218   fhELambdaPhoton->SetZTitle("#lambda_{1}^{2}");
 219   fhELambdaPhoton->SetYTitle("#lambda_{0}^{2}");
 220   fhELambdaPhoton->SetXTitle("E_{Photon} (GeV)");
 221   outputContainer->Add(fhELambdaPhoton) ;
 222
 223   fhDispersionPhoton  = new TH2F
 224   ("hLambdaDispersion","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2);
 225   fhDispersionPhoton->SetYTitle("Dispersion");
 226   fhDispersionPhoton->SetXTitle("p_{T Photon} (GeV/c)");
 227   outputContainer->Add(fhDispersionPhoton) ;
 228
 229   fhELambdaCellPhoton  = new TH3F
 230   ("hELambdaCellPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5);
 231   fhELambdaCellPhoton->SetZTitle("N Cells");
 232   fhELambdaCellPhoton->SetYTitle("#lambda_{0}^{2}");
 233   fhELambdaCellPhoton->SetXTitle("E_{Photon} (GeV)");
 234   outputContainer->Add(fhELambdaCellPhoton) ;
 235
 236   //Identified Pi0 histograms
 237   fhNCellPi0 = new TH2F ("hNCellPi0","Number of cell per cluster",nptbins,ptmin,ptmax,21,-0.5,20.5);
 238   fhNCellPi0->SetXTitle("E");
 239   fhNCellPi0->SetYTitle("N_{Cells}");
 240   outputContainer->Add(fhNCellPi0) ;
 241
 242   fhLambdaPi0  = new TH3F
 243   ("hLambdaPi0","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2);
 244   fhLambdaPi0->SetZTitle("#lambda_{1}^{2}");
 245   fhLambdaPi0->SetYTitle("#lambda_{0}^{2}");
 246   fhLambdaPi0->SetXTitle("p_{T Pi0} (GeV/c)");
 247   outputContainer->Add(fhLambdaPi0) ;
 248
 249   fhELambdaPi0  = new TH3F
 250   ("hELambdaPi0","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2);
 251   fhELambdaPi0->SetZTitle("#lambda_{1}^{2}");
 252   fhELambdaPi0->SetYTitle("#lambda_{0}^{2}");
 253   fhELambdaPi0->SetXTitle("E_{Pi0} (GeV)");
 254   outputContainer->Add(fhELambdaPi0) ;
 255
 256   fhDispersionPi0  = new TH2F
 257   ("hLambdaDispersion","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2);
 258   fhDispersionPi0->SetYTitle("Dispersion");
 259   fhDispersionPi0->SetXTitle("p_{T Pi0} (GeV/c)");
 260   outputContainer->Add(fhDispersionPi0) ;
 261
 262   fhELambdaCellPi0  = new TH3F
 263   ("hELambdaCellPi0","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5);
 264   fhELambdaCellPi0->SetZTitle("N Cells");
 265   fhELambdaCellPi0->SetYTitle("#lambda_{0}^{2}");
 266   fhELambdaCellPi0->SetXTitle("E_{Pi0} (GeV)");
 267   outputContainer->Add(fhELambdaCellPi0) ;
 268
 269   //Identified charged hadron histograms
 270   fhNCellChargedHadron = new TH2F ("hNCellChargedHadron","Number of cell per cluster",nptbins,ptmin,ptmax,21,-0.5,20.5);
 271   fhNCellChargedHadron->SetXTitle("p_{T}");
 272   fhNCellChargedHadron->SetYTitle("N_{Cells}");
 273   outputContainer->Add(fhNCellChargedHadron) ;
 274
 275   fhLambdaChargedHadron  = new TH3F
 276   ("hLambdaChargedHadron","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2);
 277   fhLambdaChargedHadron->SetZTitle("#lambda_{1}^{2}");
 278   fhLambdaChargedHadron->SetYTitle("#lambda_{0}^{2}");
 279   fhLambdaChargedHadron->SetXTitle("p_{T ChargedHadron} (GeV/c)");
 280   outputContainer->Add(fhLambdaChargedHadron) ;
 281
 282   fhELambdaChargedHadron  = new TH3F
 283   ("hELambdaChargedHadron","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,200,0,2);
 284   fhELambdaChargedHadron->SetZTitle("#lambda_{1}^{2}");
 285   fhELambdaChargedHadron->SetYTitle("#lambda_{0}^{2}");
 286   fhELambdaChargedHadron->SetXTitle("E_{ChargedHadron} (GeV)");
 287   outputContainer->Add(fhELambdaChargedHadron) ;
 288
 289   fhDispersionChargedHadron  = new TH2F
 290   ("hLambdaDispersion","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2);
 291   fhDispersionChargedHadron->SetYTitle("Dispersion");
 292   fhDispersionChargedHadron->SetXTitle("p_{T ChargedHadron} (GeV/c)");
 293   outputContainer->Add(fhDispersionChargedHadron) ;
 294
 295   fhELambdaCellChargedHadron  = new TH3F
 296   ("hELambdaCellChargedHadron","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5);
 297   fhELambdaCellChargedHadron->SetZTitle("N Cells");
 298   fhELambdaCellChargedHadron->SetYTitle("#lambda_{0}^{2}");
 299   fhELambdaCellChargedHadron->SetXTitle("E_{ChargedHadron} (GeV)");
 300   outputContainer->Add(fhELambdaCellChargedHadron) ;
 301
 302   if(IsDataMC()){
 303     fhDeltaE  = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50);
 304     fhDeltaE->SetXTitle("#Delta E (GeV)");
 305     outputContainer->Add(fhDeltaE);
 306
 307     fhDeltaPt  = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50);
 308     fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
 309     outputContainer->Add(fhDeltaPt);
 310
 311     fhRatioE  = new TH1F ("hRatioE","Reco/MC E ", 1000,0,2);
 312     fhRatioE->SetXTitle("E_{reco}/E_{gen}");
 313     outputContainer->Add(fhRatioE);
 314
 315     fhRatioPt  = new TH1F ("hRatioPt","Reco/MC p_{T} ", 1000,0,2);
 316     fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
 317     outputContainer->Add(fhRatioPt);
 318
 319     fh2E  = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 320     fh2E->SetXTitle("E_{rec} (GeV)");
 321     fh2E->SetYTitle("E_{gen} (GeV)");
 322     outputContainer->Add(fh2E);
 323
 324     fh2Pt  = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 325     fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
 326     fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
 327     outputContainer->Add(fh2Pt);
 328
 329     fhMCPdg  = new TH2F ("hMCPdg","PDG Code distribution",nptbins,ptmin,ptmax,6001,-3000.5,3000.5);
 330     fhMCPdg->SetXTitle("E_{Reco Cluster} (GeV)");
 331     fhMCPdg->SetYTitle("PDG Code");
 332     outputContainer->Add(fhMCPdg);
 333
 334     fhEMCPhoton  = new TH1F("hEMCPhoton","Number of #gamma over calorimeter",nptbins,ptmin,ptmax);
 335     //fhEMCPhoton->SetYTitle("N");
 336     fhEMCPhoton->SetXTitle("E_{#gamma}(GeV)");
 337     outputContainer->Add(fhEMCPhoton) ;
 338
 339     fhPhiMCPhoton  = new TH2F
 340     ("hPhiMCPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 341     fhPhiMCPhoton->SetYTitle("#phi");
 342     fhPhiMCPhoton->SetXTitle("E_{ #gamma} (GeV)");
 343     outputContainer->Add(fhPhiMCPhoton) ;
 344
 345     fhEtaMCPhoton  = new TH2F
 346     ("hEtaMCPhoton","#phi_{#gamma}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 347     fhEtaMCPhoton->SetYTitle("#eta");
 348     fhEtaMCPhoton->SetXTitle("E_{ #gamma} (GeV)");
 349     outputContainer->Add(fhEtaMCPhoton) ;
 350
 351     fhLambdaMCPhoton  = new TH3F
 352     ("hLambdaMCPhoton","#lambda_{#gamma}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5);
 353     fhLambdaMCPhoton->SetZTitle("N_{Cells}");
 354     fhLambdaMCPhoton->SetYTitle("#lambda_{0}^{2}");
 355     fhLambdaMCPhoton->SetXTitle("E_{#gamma} (GeV)");
 356     outputContainer->Add(fhLambdaMCPhoton) ;
 357
 358     fhPhotTrueE  = new TH3F
 359     ("hPhotTrueE","#lambda_{#gamma}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2);
 360     fhPhotTrueE->SetZTitle("#lambda_{0}^{2}");
 361     fhPhotTrueE->SetYTitle("Recons. E_{#gamma}");
 362     fhPhotTrueE->SetXTitle("MC Truth E_{#gamma} (GeV)");
 363     outputContainer->Add(fhPhotTrueE) ;
 364
 365     fhRatioEPhoton  = new TH1F ("hRatioEPhoton","E_{Reco #gamma}/E_{MC truth #gamma}", 1000,0,2);
 366     fhRatioEPhoton->SetXTitle("E_{reco #gamma}/E_{gen #gamma}");
 367     outputContainer->Add(fhRatioEPhoton);
 368
 369     fhEMCPi0  = new TH1F("hEMCPi0","Number of #pi^0 over calorimeter",nptbins,ptmin,ptmax);
 370     fhEMCPi0->SetYTitle("N");
 371     fhEMCPi0->SetXTitle("E_{ #pi^0}(GeV)");
 372     outputContainer->Add(fhEMCPi0) ;
 373
 374     fhPhiMCPi0  = new TH2F
 375     ("hPhiMCPi0","#phi_{#pi^0}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 376     fhPhiMCPi0->SetYTitle("#phi");
 377     fhPhiMCPi0->SetXTitle("E_{ #pi^0} (GeV)");
 378     outputContainer->Add(fhPhiMCPi0) ;
 379
 380     fhEtaMCPi0  = new TH2F
 381     ("hEtaMCPi0","#phi_{#pi^0}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 382     fhEtaMCPi0->SetYTitle("#eta");
 383     fhEtaMCPi0->SetXTitle("E_{ #pi^0} (GeV)");
 384     outputContainer->Add(fhEtaMCPi0) ;
 385
 386     fhLambdaMCPi0  = new TH3F
 387     ("hLambdaMCPi0","#lambda_{#pi^0}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5);
 388     fhLambdaMCPi0->SetZTitle("N_{Cells}");
 389     fhLambdaMCPi0->SetYTitle("#lambda_{0}^{2}");
 390     fhLambdaMCPi0->SetXTitle("E_{ #pi^0} (GeV)");
 391     outputContainer->Add(fhLambdaMCPi0) ;
 392
 393     fhProductionDistance  = new TH2F
 394     ("hProductionDistance","#lambda_{#pi^0}",nptbins,ptmin,ptmax,160,0,800);
 395     fhProductionDistance->SetYTitle("Distance of production of Pi0 from vertex (cm)");
 396     fhProductionDistance->SetXTitle("E_{ #pi^0} (GeV)");
 397     outputContainer->Add(fhProductionDistance) ;
 398
 399     fhProductionRadius  = new TH2F
 400     ("hProductionRadius","#lambda_{#pi^0}",nptbins,ptmin,ptmax,160,0,800);
 401     fhProductionRadius->SetYTitle("Radius of production of Pi0 from vertex (cm)");
 402     fhProductionRadius->SetXTitle("E_{ #pi^0} (GeV)");
 403     outputContainer->Add(fhProductionRadius) ;
 404
 405     fhDecayAngle = new TH2F
 406     ("hDecayAngle","#lambda_{#pi^0}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 407     fhDecayAngle->SetYTitle("Opening angle of the #pi^{0} decay (radians)");
 408     fhDecayAngle->SetXTitle("E_{ #pi^0} (GeV)");
 409     outputContainer->Add(fhDecayAngle) ;
 410
 411     fhPi0TrueE  = new TH3F
 412     ("hPi0TrueE","#lambda_{#pi^0}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2);
 413     fhPi0TrueE->SetZTitle("#lambda_{0}^{2}");
 414     fhPi0TrueE->SetYTitle("Recons. E_{#pi^0}");
 415     fhPi0TrueE->SetXTitle("MC Truth E_{#pi^0} (GeV)");
 416     outputContainer->Add(fhPi0TrueE) ;
 417
 418     fhRatioEPi0  = new TH1F ("hRatioEPi0","E_{Reco #pi^{0}}/E_{MC truth #pi^{0}}", 1000,0,2);
 419     fhRatioEPi0->SetXTitle("E_{reco #pi^{0}}/E_{gen #pi^{0}}");
 420     outputContainer->Add(fhRatioEPi0);
 421
 422     fhEMCPion  = new TH1F("hEMCPion","Number of #pi over calorimeter",nptbins,ptmin,ptmax);
 423     fhEMCPion->SetYTitle("N");
 424     fhEMCPion->SetXTitle("E_{ #pi}(GeV)");
 425     outputContainer->Add(fhEMCPion) ;
 426
 427     fhPhiMCPion  = new TH2F
 428     ("hPhiMCPion","#phi_{#pi}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 429     fhPhiMCPion->SetYTitle("#phi");
 430     fhPhiMCPion->SetXTitle("E_{ #pi} (GeV)");
 431     outputContainer->Add(fhPhiMCPion) ;
 432
 433     fhEtaMCPion  = new TH2F
 434     ("hEtaMCPion","#phi_{#pi}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 435     fhEtaMCPion->SetYTitle("#eta");
 436     fhEtaMCPion->SetXTitle("E_{ #pi} (GeV)");
 437     outputContainer->Add(fhEtaMCPion) ;
 438
 439     fhLambdaMCPion  = new TH3F
 440     ("hLambdaMCPion","#lambda_{#pi}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5);
 441     fhLambdaMCPion->SetZTitle("N_{Cells}");
 442     fhLambdaMCPion->SetYTitle("#lambda_{0}^{2}");
 443     fhLambdaMCPion->SetXTitle("E_{ #pi} (GeV)");
 444     outputContainer->Add(fhLambdaMCPion) ;
 445
 446     fhPionTrueE  = new TH3F
 447     ("hPionTrueE","#lambda_{#pi}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2);
 448     fhPionTrueE->SetZTitle("#lambda_{0}^{2}");
 449     fhPionTrueE->SetYTitle("Recons. E_{#pi}");
 450     fhPionTrueE->SetXTitle("MC Truth E_{#pi} (GeV)");
 451     outputContainer->Add(fhPionTrueE) ;
 452
 453     fhRatioEPion  = new TH1F ("hRatioEPion","E_{Reco #pi^{#pm}}/E_{MC truth #pi^{#pm}}", 1000,0,2);
 454     fhRatioEPion->SetXTitle("E_{reco #pi^{#pm}}/E_{gen #pi^{#pm}}");
 455     outputContainer->Add(fhRatioEPion);
 456
 457     fhEMCProton  = new TH1F("hEMCProton","Number of p over calorimeter",nptbins,ptmin,ptmax);
 458     fhEMCProton->SetYTitle("N");
 459     fhEMCProton->SetXTitle("E_{ p}(GeV)");
 460     outputContainer->Add(fhEMCProton) ;
 461
 462     fhPhiMCProton  = new TH2F
 463     ("hPhiMCProton","#phi_{p}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 464     fhPhiMCProton->SetYTitle("#phi");
 465     fhPhiMCProton->SetXTitle("E_{ p} (GeV)");
 466     outputContainer->Add(fhPhiMCProton) ;
 467
 468     fhEtaMCProton  = new TH2F
 469     ("hEtaMCProton","#phi_{p}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 470     fhEtaMCProton->SetYTitle("#eta");
 471     fhEtaMCProton->SetXTitle("E_{ p} (GeV)");
 472     outputContainer->Add(fhEtaMCProton) ;
 473
 474     fhLambdaMCProton  = new TH3F
 475     ("hLambdaMCProton","#lambda_{p}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5);
 476     fhLambdaMCProton->SetZTitle("N_{Cells}");
 477     fhLambdaMCProton->SetYTitle("#lambda_{0}^{2}");
 478     fhLambdaMCProton->SetXTitle("E_{ p} (GeV)");
 479     outputContainer->Add(fhLambdaMCProton) ;
 480
 481     fhProtonTrueE  = new TH3F
 482     ("hProtonTrueE","#lambda_{p}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2);
 483     fhProtonTrueE->SetZTitle("#lambda_{0}^{2}");
 484     fhProtonTrueE->SetYTitle("Recons. E_{p}");
 485     fhProtonTrueE->SetXTitle("MC Truth E_{p} (GeV)");
 486     outputContainer->Add(fhProtonTrueE) ;
 487
 488     fhRatioEProton  = new TH1F ("hRatioEProton","E_{Reco p}/E_{MC truth p}", 1000,0,2);
 489     fhRatioEProton->SetXTitle("E_{reco p}/E_{gen p}");
 490     outputContainer->Add(fhRatioEProton);
 491
 492     fhEMCAntiProton  = new TH1F("hEMCAntiProton","Number of #bar{p} over calorimeter",nptbins,ptmin,ptmax);
 493     fhEMCAntiProton->SetYTitle("N");
 494     fhEMCAntiProton->SetXTitle("E_{#bar{p}}(GeV)");
 495     outputContainer->Add(fhEMCAntiProton) ;
 496
 497     fhPhiMCAntiProton  = new TH2F
 498     ("hPhiMCAntiProton","#phi_{#bar{p}}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 499     fhPhiMCAntiProton->SetYTitle("#phi");
 500     fhPhiMCAntiProton->SetXTitle("E_{#bar{p}} (GeV)");
 501     outputContainer->Add(fhPhiMCAntiProton) ;
 502
 503     fhEtaMCAntiProton  = new TH2F
 504     ("hEtaMCAntiProton","#phi_{#bar{p}}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 505     fhEtaMCAntiProton->SetYTitle("#eta");
 506     fhEtaMCAntiProton->SetXTitle("E_{#bar{p}} (GeV)");
 507     outputContainer->Add(fhEtaMCAntiProton) ;
 508
 509     fhLambdaMCAntiProton  = new TH3F
 510     ("hLambdaMCAntiProton","#lambda_{#bar{p}}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5);
 511     fhLambdaMCAntiProton->SetZTitle("N_{Cells}");
 512     fhLambdaMCAntiProton->SetYTitle("#lambda_{0}^{2}");
 513     fhLambdaMCAntiProton->SetXTitle("E_{#bar{p}} (GeV)");
 514     outputContainer->Add(fhLambdaMCAntiProton) ;
 515
 516     fhAntiProtonTrueE  = new TH3F
 517     ("hAntiProtonTrueE","#lambda_{#bar{p}}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2);
 518     fhAntiProtonTrueE->SetZTitle("#lambda_{0}^{2}");
 519     fhAntiProtonTrueE->SetYTitle("Recons. E_{#bar{p}}");
 520     fhAntiProtonTrueE->SetXTitle("MC Truth E_{#bar{p}} (GeV)");
 521     outputContainer->Add(fhAntiProtonTrueE) ;
 522
 523     fhRatioEAntiProton  = new TH1F ("hRatioEAntiProton","E_{Reco #bar{p}}/E_{MC truth #bar{p}}", 1000,0,2);
 524     fhRatioEAntiProton->SetXTitle("E_{reco #bar{p}}/E_{gen #bar{p}}");
 525     outputContainer->Add(fhRatioEAntiProton);
 526
 527     fhEMCNeutron  = new TH1F("hEMCNeutron","Number of p over calorimeter",nptbins,ptmin,ptmax);
 528     fhEMCNeutron->SetYTitle("N");
 529     fhEMCNeutron->SetXTitle("E_{ p}(GeV)");
 530     outputContainer->Add(fhEMCNeutron) ;
 531
 532     fhPhiMCNeutron  = new TH2F
 533     ("hPhiMCNeutron","#phi_{p}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 534     fhPhiMCNeutron->SetYTitle("#phi");
 535     fhPhiMCNeutron->SetXTitle("E_{ p} (GeV)");
 536     outputContainer->Add(fhPhiMCNeutron) ;
 537
 538     fhEtaMCNeutron  = new TH2F
 539     ("hEtaMCNeutron","#phi_{p}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 540     fhEtaMCNeutron->SetYTitle("#eta");
 541     fhEtaMCNeutron->SetXTitle("E_{ p} (GeV)");
 542     outputContainer->Add(fhEtaMCNeutron) ;
 543
 544     fhLambdaMCNeutron  = new TH3F
 545     ("hLambdaMCNeutron","#lambda_{p}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5);
 546     fhLambdaMCNeutron->SetZTitle("N_{Cells}");
 547     fhLambdaMCNeutron->SetYTitle("#lambda_{0}^{2}");
 548     fhLambdaMCNeutron->SetXTitle("E_{ p} (GeV)");
 549     outputContainer->Add(fhLambdaMCNeutron) ;
 550
 551     fhNeutronTrueE  = new TH3F
 552     ("hNeutronTrueE","#lambda_{n}",nptbins,ptmin,ptmax,nptbins,ptmin,ptmax,200,0,2);
 553     fhNeutronTrueE->SetZTitle("#lambda_{0}^{2}");
 554     fhNeutronTrueE->SetYTitle("Recons. E_{n}");
 555     fhNeutronTrueE->SetXTitle("MC Truth E_{n} (GeV)");
 556     outputContainer->Add(fhNeutronTrueE) ;
 557
 558     fhRatioENeutron  = new TH1F ("hRatioENeutron","E_{Reco n}/E_{MC truth n}", 1000,0,2);
 559     fhRatioENeutron->SetXTitle("E_{reco n}/E_{gen n}");
 560     outputContainer->Add(fhRatioENeutron);
 561
 562     fhEMCEta  = new TH1F("hEMCEta","Number of #eta over calorimeter",nptbins,ptmin,ptmax);
 563     fhEMCEta->SetYTitle("N");
 564     fhEMCEta->SetXTitle("E_{ #eta}(GeV)");
 565     outputContainer->Add(fhEMCEta) ;
 566
 567     fhPhiMCEta  = new TH2F
 568     ("hPhiMCEta","#phi_{#eta}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 569     fhPhiMCEta->SetYTitle("#phi");
 570     fhPhiMCEta->SetXTitle("E_{ #eta} (GeV)");
 571     outputContainer->Add(fhPhiMCEta) ;
 572
 573     fhEtaMCEta  = new TH2F
 574     ("hEtaMCEta","#phi_{#eta}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 575     fhEtaMCEta->SetYTitle("#eta");
 576     fhEtaMCEta->SetXTitle("E_{ #eta} (GeV)");
 577     outputContainer->Add(fhEtaMCEta) ;
 578
 579     fhLambdaMCEta  = new TH3F
 580     ("hLambdaMCEta","#lambda_{#eta}",nptbins,ptmin,ptmax,200,0,2,21,-0.5,20.5);
 581     fhLambdaMCEta->SetZTitle("N_{Cells}");
 582     fhLambdaMCEta->SetYTitle("#lambda_{0}^{2}");
 583     fhLambdaMCEta->SetXTitle("E_{ #eta} (GeV)");
 584     outputContainer->Add(fhLambdaMCEta) ;
 585
 586     fhPrimPt     = new TH1D("hPrimPt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
 587     outputContainer->Add(fhPrimPt) ;
 588
 589   }//Histos with MC
 590
 591   return outputContainer ;
 592
 593 }
 594
 595 //____________________________________________________________________________
 596 void AliAnaShowerParameter::Init()
 597 {
 598
 599   //Init
 600   //Do some checks
 601   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
 602     printf("AliAnaShowerParameter::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
 603     abort();
 604   }
 605   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
 606     printf("AliAnaShowerParameter::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
 607     abort();
 608   }
 609
 610 }
 611
 612 //____________________________________________________________________________
 613 void AliAnaShowerParameter::InitParameters()
 614 {
 615
 616   //Initialize the parameters of the analysis.
 617   AddToHistogramsName("AnaPhoton_");
 618
 619   fCalorimeter = "PHOS" ;
 620   fMinDist  = 2.;
 621   fMinDist2 = 4.;
 622   fMinDist3 = 5.;
 623   fMassCut  = 0.03; //30 MeV
 624
 625   fTimeCutMin  = -1;
 626   fTimeCutMax  = 9999999;
 627   fNCellsCut = 0;
 628
 629   fRejectTrackMatch       = kTRUE ;
 630   fCheckConversion        = kFALSE;
 631   fAddConvertedPairsToAOD = kFALSE;
 632
 633   fNumClusters = -1; // By Default, select all events.
 634
 635 }
 636
 637 //__________________________________________________________________
 638 void  AliAnaShowerParameter::MakeAnalysisFillAOD()
 639 {
 640   //Do analysis and fill aods
 641   //Search for photons in fCalorimeter
 642
 643   //Get vertex for photon momentum calculation
 644
 645   for (Int_t iev = 0; iev < GetNMixedEvent(); iev++) {
 646     if (!GetMixedEvent())
 647       GetReader()->GetVertex(GetVertex(iev));
 648     else
 649       GetMixedEvent()->GetVertexOfEvent(iev)->GetXYZ(GetVertex(iev));
 650   }
 651
 652   //Select the Calorimeter of the photon
 653   TObjArray * pl = 0x0;
 654   if(fCalorimeter == "PHOS")
 655     pl = GetAODPHOS();
 656   else if (fCalorimeter == "EMCAL")
 657     pl = GetAODEMCAL();
 658
 659   if(!pl){
 660     printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Careful cluster array NULL!!\n");
 661     return;
 662   }
 663
 664   //Fill AODCaloClusters and AODParticle with PHOS/EMCAL aods
 665   TLorentzVector mom, mom2 ;
 666   Int_t nCaloClusters = pl->GetEntriesFast();
 667   //Cut on the number of clusters in the event.
 668   if ((fNumClusters !=-1) && (nCaloClusters != fNumClusters)) return;
 669   Bool_t * indexConverted = new Bool_t[nCaloClusters];
 670   for (Int_t i = 0; i < nCaloClusters; i++)
 671     indexConverted[i] = kFALSE;
 672
 673   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
 674
 675     AliAODCaloCluster * calo =  (AliAODCaloCluster*) (pl->At(icalo));
 676     Int_t evtIndex = 0 ;
 677     if (GetMixedEvent()) {
 678       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
 679     }
 680     //Cluster selection, not charged, with photon id and in fiducial cut
 681
 682     //Input from second AOD?
 683     Int_t input = 0;
 684
 685     //Get Momentum vector,
 686     if (input == 0)
 687       calo->GetMomentum(mom,GetVertex(evtIndex)) ;//Assume that come from vertex in straight line
 688
 689     //Skip the cluster if it doesn't fit inside the cuts.
 690     if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
 691     Double_t tof = calo->GetTOF()*1e9;
 692     if(tof < fTimeCutMin || tof > fTimeCutMax) continue;
 693     if(calo->GetNCells() <= fNCellsCut) continue;
 694
 695     //Check acceptance selection
 696     if(IsFiducialCutOn()){
 697       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
 698       if(! in ) continue ;
 699     }
 700
 701     //Create AOD for analysis
 702     AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
 703     Int_t label = calo->GetLabel();
 704     aodph.SetLabel(label);
 705     aodph.SetInputFileIndex(input);
 706
 707     //Set the indices of the original caloclusters
 708     aodph.SetCaloLabel(calo->GetID(),-1);
 709     aodph.SetDetector(fCalorimeter);
 710     if(GetDebug() > 1)
 711       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());
 712
 713     //Check Distance to Bad channel, set bit.
 714     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
 715     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
 716     if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
 717       continue ;
 718
 719     if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Bad channel cut passed %4.2f\n",distBad);
 720
 721     if     (distBad > fMinDist3) aodph.SetDistToBad(2) ;
 722     else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
 723     else                         aodph.SetDistToBad(0) ;
 724
 725     //Skip matched clusters with tracks
 726     if(fRejectTrackMatch && IsTrackMatched(calo)) continue ;
 727     if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - TrackMatching cut passed \n");
 728
 729     //Set PID bits for later selection (AliAnaPi0 for example)
 730     //GetPDG already called in SetPIDBits.
 731     GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils());
 732     if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - PID Bits set \n");
 733
 734     //Play with the MC stack if available
 735     //Check origin of the candidates
 736     if(IsDataMC()){
 737       aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
 738       if(GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
 739     }
 740
 741     // Check if cluster comes from a conversion in the material in front of the calorimeter
 742     // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
 743
 744     if(fCheckConversion && nCaloClusters > 1){
 745       Bool_t bConverted = kFALSE;
 746       Int_t id2 = -1;
 747
 748       //Check if set previously as converted couple, if so skip its use.
 749       if (indexConverted[icalo]) continue;
 750
 751       for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
 752         //Check if set previously as converted couple, if so skip its use.
 753         if (indexConverted[jcalo]) continue;
 754         //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
 755         AliAODCaloCluster * calo2 =  (AliAODCaloCluster*) (pl->At(jcalo));              //Get cluster kinematics
 756         Int_t evtIndex2 = 0 ;
 757         if (GetMixedEvent()) {
 758           evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ;
 759         }
 760         calo2->GetMomentum(mom2,GetVertex(evtIndex2));
 761         //Check only certain regions
 762         Bool_t in2 = kTRUE;
 763         if(IsFiducialCutOn()) in2 =  GetFiducialCut()->IsInFiducialCut(mom2,fCalorimeter) ;
 764         if(!in2) continue;
 765
 766         //Get mass of pair, if small, take this pair.
 767         //printf("\t both in calo, mass %f, cut %f\n",(mom+mom2).M(),fMassCut);
 768         if((mom+mom2).M() < fMassCut){
 769           bConverted = kTRUE;
 770           id2 = calo2->GetID();
 771           indexConverted[jcalo]=kTRUE;
 772           break;
 773         }
 774
 775       }//Mass loop
 776
 777       if(bConverted){
 778         if(fAddConvertedPairsToAOD){
 779           //Create AOD of pair analysis
 780           TLorentzVector mpair = mom+mom2;
 781           AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
 782           aodpair.SetLabel(aodph.GetLabel());
 783           aodpair.SetInputFileIndex(input);
 784
 785           //printf("Index %d, Id %d\n",icalo, calo->GetID());
 786           //Set the indeces of the original caloclusters
 787           aodpair.SetCaloLabel(calo->GetID(),id2);
 788           aodpair.SetDetector(fCalorimeter);
 789           aodpair.SetPdg(aodph.GetPdg());
 790           aodpair.SetTag(aodph.GetTag());
 791
 792           //Add AOD with pair object to aod branch
 793           AddAODParticle(aodpair);
 794           //printf("\t \t both added pair\n");
 795         }
 796
 797         //Do not add the current calocluster
 798         continue;
 799       }//converted pair
 800     }//check conversion
 801
 802     //Add AOD with photon object to aod branch
 803     AddAODParticle(aodph);
 804
 805   }//loop;
 806   delete [] indexConverted;
 807
 808   if(GetDebug() > 1) printf("AliAnaShowerParameter::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
 809
 810 }
 811
 812 //__________________________________________________________________
 813 void  AliAnaShowerParameter::MakeAnalysisFillHistograms()
 814 {
 815
 816   //Do analysis and fill histograms
 817
 818   // Access MC information in stack if requested, check that it exists.
 819   AliStack * stack = 0x0;
 820   TParticle * primary = 0x0;
 821   TClonesArray * mcparticles0 = 0x0;
 822   //TClonesArray * mcparticles1 = 0x0;
 823   AliAODMCParticle * aodprimary = 0x0;
 824   TObjArray * pl = 0x0;
 825   Int_t iNumCell=0;
 826
 827   //Check if the stack is available when analysing MC data.
 828   if(IsDataMC()){
 829
 830     if(GetReader()->ReadStack()){
 831       stack =  GetMCStack() ;
 832       if(!stack) {
 833         printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
 834                                 abort();
 835       }
 836
 837     }
 838     else if(GetReader()->ReadAODMCParticles()){
 839
 840       //Get the list of MC particles
 841       mcparticles0 = GetReader()->GetAODMCParticles(0);
 842       if(!mcparticles0 && GetDebug() > 0)       {
 843         printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
 844       }
 845     }// is data and MC
 846   }
 847   //Loop on stored AOD photons
 848   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 849   fhNClusters->Fill(naod);
 850   if(GetDebug() > 0) printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
 851
 852   for(Int_t iaod = 0; iaod < naod ; iaod++){
 853     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 854
 855     if(ph->GetDetector() != fCalorimeter) continue;
 856
 857     if(GetDebug() > 2)
 858       printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
 859
 860     //Fill Cluster histograms
 861     Float_t ptcluster  = ph->Pt();
 862     Float_t phicluster = ph->Phi();
 863     Float_t etacluster = ph->Eta();
 864     Float_t ecluster   = ph->E();
 865     Float_t lambdaMainCluster = -1; //Can't set their values here, but need a default -1 value in case something is wrong
 866     Float_t lambdaSecondCluster = -1;
 867     Float_t dispcluster = -1;
 868
 869     //Get the list of clusters from the AOD and loop over them to find the onces corresponding to the reconstructed 'photons'.
 870     if(fCalorimeter == "PHOS")
 871       pl = GetAODPHOS();
 872     else if (fCalorimeter == "EMCAL")
 873       pl = GetAODEMCAL();
 874     if(pl){
 875       //Some values are stored in AliAODCaloCluster objects only; we need to fetch them.
 876       for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++){
 877         AliAODCaloCluster * calo =  (AliAODCaloCluster*) (pl->At(icalo));
 878         if (calo->GetLabel()==ph->GetLabel()) {  //The Cluster is the right one for this particle
 879           lambdaMainCluster = calo->GetM02();   //lambda_0
 880           lambdaSecondCluster = calo->GetM20();     //lambda_1
 881           dispcluster = calo->GetDispersion();    //Dispersion
 882           iNumCell = calo->GetNCells();
 883           if(GetDebug() > 2)
 884             printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Cluster Lambda0: %3.2f, Lambda1: %3.2f, Dispersion: %3.2f, NCells: %d \n",lambdaMainCluster,lambdaSecondCluster,dispcluster,iNumCell) ;
 885         }
 886       }
 887     }
 888
 889     //Fill the basic non-MC information about the cluster.
 890     fhPtCluster  ->Fill(ecluster);
 891     fhPhiCluster ->Fill(ecluster,phicluster);
 892     fhEtaCluster ->Fill(ecluster,etacluster);
 893     fhDispersionCluster->Fill(ptcluster,dispcluster);
 894     fhLambdaCluster ->Fill(ptcluster,lambdaMainCluster,lambdaSecondCluster);
 895     fhELambdaCluster ->Fill(ecluster,lambdaMainCluster,lambdaSecondCluster);
 896     fhELambdaCellCluster ->Fill(ecluster,lambdaMainCluster,iNumCell);
 897     fhNCellCluster->Fill(ecluster,iNumCell);
 898
 899     //In the case of an event with several clusters, calculate DeltaEta and DeltaPhi for the cluster pairs.
 900     for (Int_t jaod=iaod+1;jaod<naod;jaod++){
 901       AliAODPWG4Particle* phSecond = (AliAODPWG4Particle*) (GetOutputAODBranch()->At(jaod));
 902       if(phSecond->GetDetector() != fCalorimeter) continue;
 903       fhDeltaPhiClusters->Fill(phicluster-(phSecond->Phi()));
 904       fhDeltaEtaClusters->Fill(etacluster-(phSecond->Eta()));
 905     }
 906
 907     //Fill the lambda distribution for identified particles
 908     if(ph->GetPdg() ==  AliCaloPID::kPhoton){
 909       fhDispersionPhoton->Fill(ptcluster,dispcluster);
 910       fhLambdaPhoton ->Fill(ptcluster,lambdaMainCluster,lambdaSecondCluster);
 911       fhELambdaPhoton ->Fill(ecluster,lambdaMainCluster,lambdaSecondCluster);
 912       fhELambdaCellPhoton ->Fill(ecluster,lambdaMainCluster,iNumCell);
 913       fhNCellPhoton->Fill(ecluster,iNumCell);
 914     }
 915     if(ph->GetPdg() ==  AliCaloPID::kPi0){
 916       fhDispersionPi0->Fill(ptcluster,dispcluster);
 917       fhLambdaPi0 ->Fill(ptcluster,lambdaMainCluster,lambdaSecondCluster);
 918       fhELambdaPi0 ->Fill(ecluster,lambdaMainCluster,lambdaSecondCluster);
 919       fhELambdaCellPi0 ->Fill(ecluster,lambdaMainCluster,iNumCell);
 920       fhNCellPi0->Fill(ecluster,iNumCell);
 921     }
 922     if(ph->GetPdg() ==  AliCaloPID::kChargedHadron){
 923       fhDispersionChargedHadron->Fill(ptcluster,dispcluster);
 924       fhLambdaChargedHadron ->Fill(ptcluster,lambdaMainCluster,lambdaSecondCluster);
 925       fhELambdaChargedHadron ->Fill(ecluster,lambdaMainCluster,lambdaSecondCluster);
 926       fhELambdaCellChargedHadron ->Fill(ecluster,lambdaMainCluster,iNumCell);
 927       fhNCellChargedHadron->Fill(ecluster,iNumCell);
 928     }
 929
 930     //Play with the MC data if available
 931     if(IsDataMC()){
 932       if(GetReader()->ReadStack() && !stack) return;
 933       if(GetReader()->ReadAODMCParticles() && !mcparticles0) return;
 934
 935       //Get the tag from AliMCAnalysisUtils for PID
 936       Int_t tag =ph->GetTag();
 937       if ( ph->GetLabel() < 0){
 938         if(GetDebug() > 0)
 939           printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - Label is -1; problem in the MC ESD? ");
 940         continue;
 941       }
 942
 943
 944       //Check if the tag matches on of the different particle types and fill the corresponding histograms
 945       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).
 946       {
 947         fhEMCPhoton  ->Fill(ecluster);
 948         fhPhiMCPhoton ->Fill(ecluster,phicluster);
 949         fhEtaMCPhoton ->Fill(ecluster,etacluster);
 950         fhLambdaMCPhoton ->Fill(ecluster,lambdaMainCluster,iNumCell);
 951         Double_t ePhot;
 952         //Get the true energy of the photon
 953         Int_t iCurrent = ph->GetLabel();
 954         TParticle* pCurrent = stack->Particle(iCurrent);
 955         ePhot = pCurrent->Energy();
 956         fhPhotTrueE->Fill(ePhot,ecluster,lambdaMainCluster);
 957         fhRatioEPhoton->Fill(ecluster/ePhot);
 958         fhMCPdg->Fill(ecluster,22);
 959         if(GetDebug() > 3)
 960           printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCPhoton: True E: %3.2f, Reco E: %3.2f, Lambda0:  %3.2f\n",ePhot,ecluster,lambdaMainCluster);
 961       }//kMCPhoton
 962
 963       if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)) //kMCPi0 : single cluster created by 2 photons from the same Pi0 decay
 964       {
 965
 966         //Fill the basic information about the Pi0 cluster
 967         fhEMCPi0  ->Fill(ecluster);
 968         fhPhiMCPi0 ->Fill(ecluster,phicluster);
 969         fhEtaMCPi0 ->Fill(ecluster,etacluster);
 970         fhLambdaMCPi0 ->Fill(ecluster,lambdaMainCluster,iNumCell);
 971
 972         //Check the position of the Pi0: does it come from the vertex or was it created by hadron-material collision?
 973         Int_t iCurrent = ph->GetLabel();
 974         Double_t ePi0 = 0;
 975         TParticle* pCurrent = stack->Particle(iCurrent);
 976         while (pCurrent->GetPdgCode()!=111)
 977         {
 978           iCurrent = pCurrent->GetFirstMother();
 979           pCurrent = stack->Particle(iCurrent);
 980         }
 981         Float_t fDistance = pCurrent->Rho();
 982         Float_t fRadius = pCurrent->R();
 983         ePi0 = pCurrent->Energy();
 984         fhProductionDistance->Fill(ecluster,fDistance);
 985         fhProductionRadius->Fill(ecluster,fRadius);
 986         //Compare the cluster energy and true energy
 987         fhPi0TrueE->Fill(ePi0,ecluster,lambdaMainCluster);
 988         fhRatioEPi0->Fill(ecluster/ePi0);
 989         fhMCPdg->Fill(ecluster,111);
 990         if(GetDebug() > 3)
 991           printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCPi0: True E: %3.2f, Reco E: %3.2f, Lambda0:  %3.2f\n",ePi0,ecluster,lambdaMainCluster);
 992       }
 993       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion))//kMCPion
 994       {
 995         fhEMCPion  ->Fill(ecluster);
 996         fhPhiMCPion ->Fill(ecluster,phicluster);
 997         fhEtaMCPion ->Fill(ecluster,etacluster);
 998         fhLambdaMCPion ->Fill(ecluster,lambdaMainCluster,iNumCell);
 999
1000         Int_t iCurrent = ph->GetLabel();
1001         Double_t ePion = 0;
1002         TParticle* pCurrent = stack->Particle(iCurrent);
1003         while ((TMath::Abs(pCurrent->GetPdgCode())!=211)&&(iCurrent>=0))
1004         {
1005           if (iCurrent!=pCurrent->GetFirstMother()) iCurrent = pCurrent->GetFirstMother(); //Avoiding infinite loops
1006           if (iCurrent>=0) pCurrent = stack->Particle(iCurrent);
1007         }
1008         if ((TMath::Abs(pCurrent->GetPdgCode())==211)&&(iCurrent>=0))
1009         {
1010           ePion = pCurrent->Energy();
1011           fhPionTrueE->Fill(ePion,ecluster,lambdaMainCluster);
1012           fhRatioEPion->Fill(ecluster/ePion);
1013         }
1014         fhMCPdg->Fill(ecluster,211);
1015         if(GetDebug() > 3)
1016           printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCPion: True E: %3.2f, Reco E: %3.2f, Lambda0:  %3.2f\n",ePion,ecluster,lambdaMainCluster);
1017       }
1018       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton))//kMCProton
1019       {
1020         fhEMCProton  ->Fill(ecluster);
1021         fhPhiMCProton ->Fill(ecluster,phicluster);
1022         fhEtaMCProton ->Fill(ecluster,etacluster);
1023         fhLambdaMCProton ->Fill(ecluster,lambdaMainCluster,iNumCell);
1024
1025         Int_t iCurrent = ph->GetLabel();
1026         Double_t eProton = 0;
1027         TParticle* pCurrent = stack->Particle(iCurrent);
1028         while ((TMath::Abs(pCurrent->GetPdgCode())!=2212)&&(iCurrent>=0))
1029         {
1030           if (iCurrent!=pCurrent->GetFirstMother()) iCurrent = pCurrent->GetFirstMother(); //Avoiding infinite loops
1031           if (iCurrent>=0) pCurrent = stack->Particle(iCurrent);
1032         }
1033         if ((TMath::Abs(pCurrent->GetPdgCode())==2212)&&(iCurrent>=0))
1034         {
1035           eProton = pCurrent->Energy();
1036           fhProtonTrueE->Fill(eProton,ecluster,lambdaMainCluster);
1037           fhRatioEProton->Fill(ecluster/eProton);
1038         }
1039         fhMCPdg->Fill(ecluster,2212);
1040         if(GetDebug() > 3)
1041           printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCProton: True E: %3.2f, Reco E: %3.2f, Lambda0:  %3.2f\n",eProton,ecluster,lambdaMainCluster);
1042       }
1043       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))//kMCAntiProton
1044       {
1045         fhEMCAntiProton  ->Fill(ecluster);
1046         fhPhiMCAntiProton ->Fill(ecluster,phicluster);
1047         fhEtaMCAntiProton ->Fill(ecluster,etacluster);
1048         fhLambdaMCAntiProton ->Fill(ecluster,lambdaMainCluster,iNumCell);
1049
1050         Int_t iCurrent = ph->GetLabel();
1051         Double_t eAntiProton = 0;
1052         TParticle* pCurrent = stack->Particle(iCurrent);
1053         while ((TMath::Abs(pCurrent->GetPdgCode())!=-2212)&&(iCurrent>=0))
1054         {
1055           if (iCurrent!=pCurrent->GetFirstMother()) iCurrent = pCurrent->GetFirstMother(); //Avoiding infinite loops
1056           if (iCurrent>=0) pCurrent = stack->Particle(iCurrent);
1057         }
1058         if ((TMath::Abs(pCurrent->GetPdgCode())==-2212)&&(iCurrent>=0))
1059         {
1060           eAntiProton = pCurrent->Energy();
1061           fhAntiProtonTrueE->Fill(eAntiProton,ecluster,lambdaMainCluster);
1062           fhRatioEAntiProton->Fill(ecluster/eAntiProton);
1063         }
1064         fhMCPdg->Fill(ecluster,-2212);
1065         if(GetDebug() > 3)
1066           printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCAntiProton: True E: %3.2f, Reco E: %3.2f, Lambda0:  %3.2f\n",eAntiProton,ecluster,lambdaMainCluster);
1067       }
1068       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCNeutron))//kMCNeutron
1069       {
1070         fhEMCNeutron  ->Fill(ecluster);
1071         fhPhiMCNeutron ->Fill(ecluster,phicluster);
1072         fhEtaMCNeutron ->Fill(ecluster,etacluster);
1073         fhLambdaMCNeutron ->Fill(ecluster,lambdaMainCluster,iNumCell);
1074
1075         Int_t iCurrent = ph->GetLabel();
1076         Double_t eNeutron = 0;
1077         TParticle* pCurrent = stack->Particle(iCurrent);
1078         while ((TMath::Abs(pCurrent->GetPdgCode())!=2112)&&(iCurrent>=0))
1079         {
1080           if (iCurrent!=pCurrent->GetFirstMother()) iCurrent = pCurrent->GetFirstMother(); //Avoiding infinite loops
1081           if (iCurrent>=0) pCurrent = stack->Particle(iCurrent);
1082         }
1083         if ((TMath::Abs(pCurrent->GetPdgCode())==2112)&&(iCurrent>=0))
1084         {
1085           eNeutron = pCurrent->Energy();
1086           fhNeutronTrueE->Fill(eNeutron,ecluster,lambdaMainCluster);
1087           fhRatioENeutron->Fill(ecluster/eNeutron);
1088         }
1089         fhMCPdg->Fill(ecluster,2112);
1090                   if(GetDebug() > 3)
1091           printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() - MC Cluster identified as kMCNeutron: True E: %3.2f, Reco E: %3.2f, Lambda0:  %3.2f\n",eNeutron,ecluster,lambdaMainCluster);
1092       }
1093       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))//kMCEta
1094       {
1095         fhEMCEta  ->Fill(ecluster);
1096         fhPhiMCEta ->Fill(ecluster,phicluster);
1097         fhEtaMCEta ->Fill(ecluster,etacluster);
1098         fhLambdaMCEta ->Fill(ecluster,lambdaMainCluster,iNumCell);
1099       }
1100
1101       // Access MC information in stack if requested, check that it exists.
1102       Int_t label =ph->GetLabel();
1103       if(label < 0) {
1104         printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
1105         continue;
1106       }
1107
1108       //Calculate Pi0 decay angles
1109       if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
1110         for(Int_t i=0 ; i<stack->GetNprimary(); i++){
1111           TParticle * prim = stack->Particle(i) ;
1112           if(prim->GetPdgCode() == 111){
1113             TLorentzVector mom1 ;
1114             (stack->Particle((prim->GetFirstDaughter())))->Momentum(mom1);
1115             TLorentzVector mom2 ;
1116             (stack->Particle((prim->GetLastDaughter())))->Momentum(mom2);
1117             Float_t fDecayAngle = mom1.Angle(mom2.Vect());
1118             fhDecayAngle->Fill(ecluster,fDecayAngle);
1119           }
1120         }
1121       }
1122
1123
1124       Float_t eprim   = 0;
1125       Float_t ptprim  = 0;
1126       if(GetReader()->ReadStack()){
1127
1128         if(label >=  stack->GetNtrack()) {
1129           if(GetDebug() > 2)  printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
1130           continue ;
1131         }
1132
1133         primary = stack->Particle(label);
1134         if(!primary){
1135           printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
1136           continue;
1137         }
1138         eprim   = primary->Energy();
1139         ptprim  = primary->Pt();
1140       }
1141       else if(GetReader()->ReadAODMCParticles()){
1142         //Check which is the input
1143         if(ph->GetInputFileIndex() == 0){
1144           if(!mcparticles0) continue;
1145           if(label >=  mcparticles0->GetEntriesFast()) {
1146             if(GetDebug() > 2)  printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",label, mcparticles0->GetEntriesFast());
1147             continue ;
1148           }
1149           //Get the particle
1150           aodprimary = (AliAODMCParticle*) mcparticles0->At(label);
1151
1152         }
1153         //        else {//Second input
1154         //          if(!mcparticles1) continue;
1155         //          if(label >=  mcparticles1->GetEntriesFast()) {
1156         //            if(GetDebug() > 2)  printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",label, mcparticles1->GetEntriesFast());
1157         //            continue ;
1158         //          }
1159         //          //Get the particle
1160         //          aodprimary = (AliAODMCParticle*) mcparticles1->At(label);
1161         //        }//second input
1162
1163         if(!aodprimary){
1164           printf("AliAnaShowerParameter::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
1165           continue;
1166         }
1167
1168         eprim   = aodprimary->E();
1169         ptprim  = aodprimary->Pt();
1170       }
1171       //Write down MC truth information concerning all clusters
1172       fh2E     ->Fill(ecluster, eprim);
1173       fh2Pt    ->Fill(ptcluster, ptprim);
1174       fhDeltaE ->Fill(eprim-ecluster);
1175       fhDeltaPt->Fill(ptprim-ptcluster);
1176       if(eprim > 0)  fhRatioE  ->Fill(ecluster/eprim);
1177       if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
1178     }//Histograms with MC
1179   }// aod loop
1180 }
1181
1182
1183 //__________________________________________________________________
1184 void AliAnaShowerParameter::Print(const Option_t * opt) const
1185 {
1186   //Print some relevant parameters set for the analysis
1187
1188   if(! opt)
1189     return;
1190
1191   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1192   AliAnaPartCorrBaseClass::Print(" ");
1193
1194   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
1195   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
1196   printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1197   printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1198   printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
1199   printf("Check Pair Conversion                = %d\n",fCheckConversion);
1200   printf("Add conversion pair to AOD           = %d\n",fAddConvertedPairsToAOD);
1201   printf("Conversion pair mass cut             = %f\n",fMassCut);
1202   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
1203   printf("Number of cells in cluster is        > %f \n", fNCellsCut);
1204   printf("Number of clusters in envent is        > %d \n", fNumClusters);
1205   printf("    \n") ;
1206
1207 }