PWGGA/PHOSTasks/PHOS_Tagging/AliAnalysisTaskTaggedPhotons.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 is 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
  16 /* $Id: */
  17
  18 //_________________________________________________________________________
  19 // Analysis for Tagged Photons
  20 // Prepares all necessary histograms for calculation of
  21 // the yield of pi0 decay photon in calorimeter:
  22 // Marks photons which makes pi0 with some other and
  23 // fill invariant mass distributions for estimate background below pi0
  24 // peak so that estimate proportion of fake pairs.
  25 // Fills as well controll MC histograms with clasification of the photon origin
  26 // and check of the proportion of truly tagged photons.
  27 //
  28 //
  29 //*-- Dmitry Blau, Dmitri Peresunko
  30 //////////////////////////////////////////////////////////////////////////////
  31
  32 #include <TH1.h>
  33 #include <TH2.h>
  34 #include <TH3.h>
  35 #include <THashList.h>
  36 #include <TFile.h>
  37 #include <TROOT.h>
  38
  39 #include "AliAnalysisTaskTaggedPhotons.h"
  40 #include "AliAnalysisManager.h"
  41 #include "AliAODEvent.h"
  42 #include "AliAODEvent.h"
  43 #include "AliVCluster.h"
  44 #include "AliCaloPhoton.h"
  45 #include "AliAODMCParticle.h"
  46 #include "AliAnalysisManager.h"
  47 #include "AliLog.h"
  48 #include "TGeoManager.h"
  49 #include "AliMCEventHandler.h"
  50 #include "AliMCEvent.h"
  51 #include "AliStack.h"
  52 #include "AliPHOSGeometry.h"
  53 #include "AliTriggerAnalysis.h"
  54 #include "AliPHOSTriggerUtils.h"
  55 #include "AliEMCALGeometry.h"
  56 #include "AliAnalysisUtils.h"
  57 #include "AliOADBContainer.h"
  58 #include "AliAODMCHeader.h"
  59
  60 ClassImp(AliAnalysisTaskTaggedPhotons)
  61
  62 //______________________________________________________________________________
  63 AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons() :
  64   AliAnalysisTaskSE(),
  65   fPHOSgeom(0x0),
  66   fOutputContainer(0x0),
  67   fStack(0x0),
  68   fTrackEvent(0x0),
  69   fPHOSEvent(0x0),
  70   fCurrentMixedList(0x0),
  71   fTriggerAnalysis(0x0),
  72   fUtils(0x0),
  73   fPHOSTrigUtils(0x0),
  74   fZmax(0.),
  75   fZmin(0.),
  76   fPhimax(0.),
  77   fPhimin(0.),
  78   fMinBCDistance(0.),
  79   fCentrality(0.),
  80   fCentBin(0),
  81   fIsMB(0),
  82   fIsMC(0),
  83   fIsFastMC(0)
  84 {
  85   //Deafult constructor
  86   //no memory allocations
  87   for(Int_t i=0;i<10;i++)
  88     for(Int_t j=0;j<2;j++)
  89       fPHOSEvents[i][j]=0x0 ;    //Container for PHOS photons
  90   for(Int_t i=0;i<6;i++)
  91     fPHOSBadMap[i]=0x0;
  92 }
  93 //______________________________________________________________________________
  94 AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const char *name) :
  95   AliAnalysisTaskSE(name),
  96   fPHOSgeom(0x0),
  97   fOutputContainer(0x0),
  98   fStack(0x0),
  99   fTrackEvent(0x0),
 100   fPHOSEvent(0x0),
 101   fCurrentMixedList(0x0),
 102   fTriggerAnalysis(new AliTriggerAnalysis),
 103   fUtils(0x0),
 104   fPHOSTrigUtils(0x0),
 105   fZmax(-60.),
 106   fZmin(60.),
 107   fPhimax(250.),
 108   fPhimin(320.),
 109   fMinBCDistance(0.),
 110   fCentrality(0.),
 111   fCentBin(0),
 112   fIsMB(0),
 113   fIsMC(0),
 114   fIsFastMC(0)
 115 {
 116   // Constructor.
 117
 118   // Output slots #0 write into a TH1 container
 119   DefineOutput(1,THashList::Class());
 120   // Set bad channel map (empty so far)
 121   for(Int_t i=0;i<1;i++)
 122     for(Int_t j=0;j<5;j++)
 123       fPHOSEvents[i][j]=0x0 ;    //Container for PHOS photons
 124   for(Int_t i=0;i<6;i++)
 125     fPHOSBadMap[i]=0x0;
 126
 127 }
 128
 129 //____________________________________________________________________________
 130 AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const AliAnalysisTaskTaggedPhotons& ap) :
 131   AliAnalysisTaskSE(ap.GetName()),
 132   fPHOSgeom(0x0),
 133   fOutputContainer(0x0),
 134   fStack(0x0),
 135   fTrackEvent(0x0),
 136   fPHOSEvent(0x0),
 137   fCurrentMixedList(0x0),
 138   fTriggerAnalysis(new AliTriggerAnalysis),
 139   fUtils(0x0),
 140   fPHOSTrigUtils(0x0),
 141   fZmax(-60.),
 142   fZmin(60.),
 143   fPhimax(250.),
 144   fPhimin(320.),
 145   fMinBCDistance(0.),
 146   fCentrality(0.),
 147   fCentBin(0),
 148   fIsMB(0),
 149   fIsMC(0),
 150   fIsFastMC(0)
 151 {
 152   // cpy ctor
 153   fZmax=ap.fZmax ;
 154   fZmin=ap.fZmin ;
 155   fPhimax=ap.fPhimax ;
 156   fPhimin=ap.fPhimin ;
 157   for(Int_t i=0;i<1;i++)
 158     for(Int_t j=0;j<5;j++)
 159       fPHOSEvents[i][j]=0x0 ;    //Container for PHOS photons
 160   for(Int_t i=0;i<6;i++)
 161     fPHOSBadMap[i]=0x0;
 162
 163 }
 164
 165 //_____________________________________________________________________________
 166 AliAnalysisTaskTaggedPhotons& AliAnalysisTaskTaggedPhotons::operator = (const AliAnalysisTaskTaggedPhotons& ap)
 167 {
 168 // assignment operator
 169
 170   this->~AliAnalysisTaskTaggedPhotons();
 171   new(this) AliAnalysisTaskTaggedPhotons(ap);
 172   return *this;
 173 }
 174
 175 //______________________________________________________________________________
 176 AliAnalysisTaskTaggedPhotons::~AliAnalysisTaskTaggedPhotons()
 177 {
 178   // dtor
 179   if(fOutputContainer && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 180     fOutputContainer->Clear() ;
 181     delete fOutputContainer ;
 182   }
 183   for(Int_t i=0;i<10;i++)
 184     for(Int_t j=0;j<2;j++)
 185       if(fPHOSEvents[i][j]){
 186         delete fPHOSEvents[i][j] ;
 187         fPHOSEvents[i][j]=0x0 ;
 188       }
 189 }
 190 //________________________________________________________________________
 191 void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()
 192 {
 193
 194
 195   // Create the outputs containers
 196   if(fOutputContainer != NULL){
 197     delete fOutputContainer;
 198   }
 199   fOutputContainer = new THashList();
 200   fOutputContainer->SetOwner(kTRUE);
 201   fOutputContainer->SetName(GetName()) ;
 202
 203   //QA histograms
 204   fOutputContainer->Add(new TH1F("hSelEvents","Event selection", 10,0.,10.)) ;
 205
 206   //vertex distribution
 207   fOutputContainer->Add(new TH1F("hNvertexTracks","N of primary tracks from the primary vertex",150,0.,150.));
 208   fOutputContainer->Add(new TH1F("hZvertex","Z vertex",200,-50.,+50.));
 209   fOutputContainer->Add(new TH2F("hTrackMult","Charged track multiplicity",100,0.,100.,250,0.,500.));
 210   fOutputContainer->Add(new TH2F("hTrackEtaPhi","Charged track eta vs phi distribution",200,-2.,2.,200,0.,TMath::TwoPi()));
 211   fOutputContainer->Add(new TH2F("hTrackEtaPt","Charged track eta vs pt distribution",200,-2.,2.,200,0.,50.));
 212
 213   //centrality
 214   fOutputContainer->Add(new TH1F("hCentrality","Ccentrality",100,0.,100.));
 215   fOutputContainer->Add(new TH2F("hPHOSCentrality","PHOS vs centrality",100,0.,100.,100,0.,100.));
 216   fOutputContainer->Add(new TH2F("hTOF","cluster TOF",200,0.,20.,300,-3.e-6,6.e-6));
 217
 218   fOutputContainer->Add(new TH2F("hCluNXZM1","Clu (X,Z), M1"   ,64,0.5,64.5, 56,0.5,56.5));
 219   fOutputContainer->Add(new TH2F("hCluNXZM2","Clu (X,Z), M2"   ,64,0.5,64.5, 56,0.5,56.5));
 220   fOutputContainer->Add(new TH2F("hCluNXZM3","Clu (X,Z), M3"   ,64,0.5,64.5, 56,0.5,56.5));
 221   fOutputContainer->Add(new TH2F("hCluEXZM1","Clu E(X,Z), M1"  ,64,0.5,64.5, 56,0.5,56.5));
 222   fOutputContainer->Add(new TH2F("hCluEXZM2","Clu E(X,Z), M2"  ,64,0.5,64.5, 56,0.5,56.5));
 223   fOutputContainer->Add(new TH2F("hCluEXZM3","Clu E(X,Z), M3"  ,64,0.5,64.5, 56,0.5,56.5));
 224
 225   fOutputContainer->Add(new TH2F("hCluArea2M1","Clu (X,Z), M1"   ,64,0.5,64.5, 56,0.5,56.5));
 226   fOutputContainer->Add(new TH2F("hCluArea2M2","Clu (X,Z), M1"   ,64,0.5,64.5, 56,0.5,56.5));
 227   fOutputContainer->Add(new TH2F("hCluArea2M3","Clu (X,Z), M1"   ,64,0.5,64.5, 56,0.5,56.5));
 228   fOutputContainer->Add(new TH2F("hCluArea3M1","Clu (X,Z), M1"   ,64,0.5,64.5, 56,0.5,56.5));
 229   fOutputContainer->Add(new TH2F("hCluArea3M2","Clu (X,Z), M1"   ,64,0.5,64.5, 56,0.5,56.5));
 230   fOutputContainer->Add(new TH2F("hCluArea3M3","Clu (X,Z), M1"   ,64,0.5,64.5, 56,0.5,56.5));
 231
 232   fOutputContainer->Add(new TH2F("hTofM1","TOF in mod1",200,-1.e-6,1.e-6,200,0.,20.)) ;
 233   fOutputContainer->Add(new TH2F("hTofM2","TOF in mod2",200,-1.e-6,1.e-6,200,0.,20.)) ;
 234   fOutputContainer->Add(new TH2F("hTofM3","TOF in mod3",200,-1.e-6,1.e-6,200,0.,20.)) ;
 235
 236   char cPID[4][5] ;
 237   snprintf(cPID[0],5,"All") ;
 238   snprintf(cPID[1],5,"Disp");
 239   snprintf(cPID[2],5,"CPV") ;
 240   snprintf(cPID[3],5,"Both");
 241
 242
 243   const Int_t nPt=500 ;
 244   const Double_t ptMax=50. ;
 245   const Int_t nM=400 ;
 246   const Double_t mMax=1. ;
 247
 248   //QA histograms
 249    for(Int_t iPID=0; iPID<4; iPID++){
 250     fOutputContainer->Add(new TH1F(Form("hPhotM1_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M1",nPt,0.,ptMax)) ;
 251     fOutputContainer->Add(new TH1F(Form("hPhotM2_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M2",nPt,0.,ptMax)) ;
 252     fOutputContainer->Add(new TH1F(Form("hPhotM3_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M3",nPt,0.,ptMax)) ;
 253     fOutputContainer->Add(new TH1F(Form("hPhotM1A2_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M1",nPt,0.,ptMax)) ;
 254     fOutputContainer->Add(new TH1F(Form("hPhotM2A2_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M2",nPt,0.,ptMax)) ;
 255     fOutputContainer->Add(new TH1F(Form("hPhotM3A2_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M3",nPt,0.,ptMax)) ;
 256     fOutputContainer->Add(new TH1F(Form("hPhotM1A3_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M1",nPt,0.,ptMax)) ;
 257     fOutputContainer->Add(new TH1F(Form("hPhotM2A3_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M2",nPt,0.,ptMax)) ;
 258     fOutputContainer->Add(new TH1F(Form("hPhotM3A3_%s",cPID[iPID]),"Spectrum of all reconstructed particles, M3",nPt,0.,ptMax)) ;
 259    }
 260
 261
 262   const Int_t nCenBin=5;
 263   for(Int_t cen=0; cen<nCenBin; cen++){
 264
 265   for(Int_t iPID=0; iPID<4; iPID++){
 266
 267     //Inclusive spectra
 268     fOutputContainer->Add(new TH1F(Form("hPhot_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
 269     fOutputContainer->Add(new TH1F(Form("hPhot_Dist2_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
 270     fOutputContainer->Add(new TH1F(Form("hPhot_Dist3_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
 271
 272     fOutputContainer->Add(new TH1F(Form("hPhot_Area1_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
 273     fOutputContainer->Add(new TH1F(Form("hPhot_Area2_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
 274     fOutputContainer->Add(new TH1F(Form("hPhot_Area3_%s_cent%d",cPID[iPID],cen),"Spectrum of all reconstructed particles",nPt,0.,ptMax)) ;
 275
 276     for(Int_t itag=0; itag<18; itag++){
 277       fOutputContainer->Add(new TH1F(Form("hPhot_TaggedMult%d_Area1_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of multiply tagged photons",nPt,0.,ptMax)) ;
 278       fOutputContainer->Add(new TH1F(Form("hPhot_TaggedMult%d_Area2_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of multiply tagged photons",nPt,0.,ptMax)) ;
 279       fOutputContainer->Add(new TH1F(Form("hPhot_TaggedMult%d_Area3_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of multiply tagged photons",nPt,0.,ptMax)) ;
 280       fOutputContainer->Add(new TH1F(Form("hPhot_Tagged%d_Area1_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
 281       fOutputContainer->Add(new TH1F(Form("hPhot_Tagged%d_Area2_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
 282       fOutputContainer->Add(new TH1F(Form("hPhot_Tagged%d_Area3_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
 283
 284       fOutputContainer->Add(new TH1F(Form("hPhot_TrueTagged%d_%s_cent%d",itag,cPID[iPID],cen),"Spectrum of all tagged particles",nPt,0.,ptMax)) ;
 285     }
 286     for(Int_t kind=1; kind<33; kind*=2){
 287       fOutputContainer->Add(new TH1F(Form("hPhot_Isolation%d_%s_cent%d",kind,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
 288       fOutputContainer->Add(new TH1F(Form("hPhot_Isolation%d_Area1_%s_cent%d",kind,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
 289       fOutputContainer->Add(new TH1F(Form("hPhot_nTagged_Isolation%d_Area1_%s_cent%d",kind,cPID[iPID],cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
 290     }
 291   }
 292   for(Int_t kind=1; kind<33; kind*=2){
 293     fOutputContainer->Add(new TH1F(Form("hPi_Isolation%d_cent%d",kind,cen),"Spectrum of all reconstructed particles, no PID",nPt,0.,ptMax)) ;
 294   }
 295
 296
 297
 298   //Invariant mass distributions for fake corrections
 299
 300   for(Int_t iPID=0; iPID<4; iPID++){
 301     fOutputContainer->Add(new TH2F(Form("hInvM_Re_Emin1_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 302     fOutputContainer->Add(new TH2F(Form("hInvM_Re_Emin2_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 303     fOutputContainer->Add(new TH2F(Form("hInvM_Re_Emin3_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 304     fOutputContainer->Add(new TH2F(Form("hInvM_Mi_Emin1_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 305     fOutputContainer->Add(new TH2F(Form("hInvM_Mi_Emin2_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 306     fOutputContainer->Add(new TH2F(Form("hInvM_Mi_Emin3_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 307
 308     fOutputContainer->Add(new TH2F(Form("hSingleInvM_Re_Emin1_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 309     fOutputContainer->Add(new TH2F(Form("hSingleInvM_Re_Emin2_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 310     fOutputContainer->Add(new TH2F(Form("hSingleInvM_Re_Emin3_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 311     fOutputContainer->Add(new TH2F(Form("hSingleInvM_Re_Emin1_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 312     fOutputContainer->Add(new TH2F(Form("hSingleInvM_Re_Emin2_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 313     fOutputContainer->Add(new TH2F(Form("hSingleInvM_Re_Emin3_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 314
 315     fOutputContainer->Add(new TH2F(Form("hSingleInvM_Mi_Emin1_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 316     fOutputContainer->Add(new TH2F(Form("hSingleInvM_Mi_Emin2_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 317     fOutputContainer->Add(new TH2F(Form("hSingleInvM_Mi_Emin3_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 318     fOutputContainer->Add(new TH2F(Form("hSingleInvM_Mi_Emin1_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 319     fOutputContainer->Add(new TH2F(Form("hSingleInvM_Mi_Emin2_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 320     fOutputContainer->Add(new TH2F(Form("hSingleInvM_Mi_Emin3_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 321
 322
 323
 324   }
 325
 326
 327   fOutputContainer->Add(new TH2F(Form("QA_Cone1_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
 328   fOutputContainer->Add(new TH2F(Form("QA_Cone2_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
 329   fOutputContainer->Add(new TH2F(Form("QA_Cone3_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
 330   fOutputContainer->Add(new TH2F(Form("QA_PCone1_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
 331   fOutputContainer->Add(new TH2F(Form("QA_PCone2_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
 332   fOutputContainer->Add(new TH2F(Form("QA_PCone3_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
 333   fOutputContainer->Add(new TH2F(Form("QA_Pi0Cone1_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
 334   fOutputContainer->Add(new TH2F(Form("QA_Pi0Cone2_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
 335   fOutputContainer->Add(new TH2F(Form("QA_Pi0Cone3_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
 336   fOutputContainer->Add(new TH2F(Form("QA_Pi0PCone1_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
 337   fOutputContainer->Add(new TH2F(Form("QA_Pi0PCone2_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
 338   fOutputContainer->Add(new TH2F(Form("QA_Pi0PCone3_Tracks_cent%d",cen),"Two-photon inv. mass vs first photon pt",50,0.,50.,200,0.,100.)) ;
 339   }//centrality
 340
 341   //MC
 342   char partName[15][10] ;
 343   snprintf(partName[0],10,"gamma") ;
 344   snprintf(partName[1],10,"pi0");
 345   snprintf(partName[2],10,"eta") ;
 346   snprintf(partName[3],10,"omega");
 347   snprintf(partName[4],10,"K0s");
 348   snprintf(partName[5],10,"Kpm");
 349   snprintf(partName[6],10,"pipm");
 350   snprintf(partName[7],10,"n");
 351   snprintf(partName[8],10,"nbar");
 352   snprintf(partName[9],10,"p");
 353   snprintf(partName[10],10,"pbar");
 354
 355
 356   if(fIsMC){
 357
 358       fOutputContainer->Add(new TH1F("hMCConversionRadius","Clusters without label",600,0.,600.)) ;
 359       fOutputContainer->Add(new TH2F("hMCRecPi0Vtx","Secondary pi0s",100,0.,10.,600,0.,600.)) ;
 360       fOutputContainer->Add(new TH2F("hMCRecEtaVtx","Secondary etas",100,0.,10.,600,0.,600.)) ;
 361       fOutputContainer->Add(new TH2F("hMCRecOmegaVtx","Secondary etas",100,0.,10.,600,0.,600.)) ;
 362       fOutputContainer->Add(new TH2F("hMCRecEtaprVtx","Secondary etas",100,0.,10.,600,0.,600.)) ;
 363       fOutputContainer->Add(new TH2F("hMCRecK0sVtx","Secondary K0s",100,0.,10.,600,0.,600.)) ;
 364       fOutputContainer->Add(new TH2F("hMCRecK0lVtx","Secondary K0l",100,0.,10.,600,0.,600.)) ;
 365       fOutputContainer->Add(new TH2F("hMCGammaPi0MisConvR","Converted photons",400,0.,40.,600,0.,600.)) ;
 366
 367   for(Int_t ipart=0; ipart<11; ipart++){
 368     fOutputContainer->Add(new TH2F(Form("hMC%s_ptrap",partName[ipart]),"Spectrum of primary photons",100,0.,10.,200,-1.,1.)) ;
 369     fOutputContainer->Add(new TH2F(Form("hMC%s_ptphi",partName[ipart]),"Spectrum of primary photons",100,0.,10.,100,0.,TMath::TwoPi())) ;
 370     fOutputContainer->Add(new TH2F(Form("hMC_%s_vertex",partName[ipart]),"vertex",nPt,0.,ptMax,150,0.,600.)) ;
 371     for(Int_t cen=0; cen<nCenBin; cen++){
 372        fOutputContainer->Add(new TH1F(Form("hMC_all_%s_cent%d",partName[ipart],cen),"Spectum (full rapifity)",nPt,0.,ptMax)) ;
 373        fOutputContainer->Add(new TH1F(Form("hMC_unitEta_%s_cent%d",partName[ipart],cen),"Spectum, |y|<0.15",nPt,0.,ptMax)) ;
 374        fOutputContainer->Add(new TH1F(Form("hMC_rap_%s_cent%d",partName[ipart],cen),"Rapidity",100,-5.,5.)) ;
 375        fOutputContainer->Add(new TH1F(Form("hMC_phi_%s_cent%d",partName[ipart],cen),"Azimuthal angle",100,0.,TMath::TwoPi())) ;
 376     }
 377   }
 378   for(Int_t cen=0; cen<nCenBin; cen++){
 379        fOutputContainer->Add(new TH1F(Form("hMC_unitEta_GammaPi0_cent%d",cen),"Spectum, |y|<0.15",nPt,0.,ptMax)) ;
 380   }
 381
 382   for(Int_t cen=0; cen<nCenBin; cen++){
 383
 384     fOutputContainer->Add(new TH2F(Form("LabelsNPrim_cent%d",cen),"vertex",nPt,0.,ptMax,20,0.,20.)) ;
 385     fOutputContainer->Add(new TH1F(Form("LabelsGamma_cent%d",cen),"vertex",nPt,0.,ptMax)) ;
 386     fOutputContainer->Add(new TH2F(Form("LabelsGammaE_cent%d",cen),"vertex",nPt,0.,ptMax,100,0.,2.)) ;
 387
 388
 389        //Sort registered particles spectra according MC information
 390        for(Int_t iPID=0; iPID<4; iPID++){
 391          fOutputContainer->Add(new TH1F(Form("hMCRecPhoton_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
 392          fOutputContainer->Add(new TH1F(Form("hMCRecPhoton_Area1_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
 393          fOutputContainer->Add(new TH1F(Form("hMCRecPhoton_Area2_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
 394          fOutputContainer->Add(new TH1F(Form("hMCRecPhoton_Area3_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
 395          fOutputContainer->Add(new TH1F(Form("hMCRecE_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
 396          fOutputContainer->Add(new TH1F(Form("hMCRecPbar_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
 397          fOutputContainer->Add(new TH1F(Form("hMCRecNbar_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
 398          fOutputContainer->Add(new TH1F(Form("hMCRecP_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
 399          fOutputContainer->Add(new TH1F(Form("hMCRecPipm_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
 400          fOutputContainer->Add(new TH1F(Form("hMCRecKpm_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
 401          fOutputContainer->Add(new TH1F(Form("hMCRecN_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
 402          fOutputContainer->Add(new TH1F(Form("hMCRecCharg_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
 403          fOutputContainer->Add(new TH1F(Form("hMCRecNeutral_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
 404          fOutputContainer->Add(new TH1F(Form("hMCRecK0s_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
 405          fOutputContainer->Add(new TH1F(Form("hMCRecNoPRim_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. electrons", nPt,0.,ptMax )) ;
 406
 407          //Decay photons
 408          fOutputContainer->Add(new TH1F(Form("hMCRecPhotPi0_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
 409          fOutputContainer->Add(new TH1F(Form("hMCRecPhotEta_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
 410          fOutputContainer->Add(new TH1F(Form("hMCRecPhotOmega_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
 411          fOutputContainer->Add(new TH1F(Form("hMCRecPhotOther_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
 412          fOutputContainer->Add(new TH1F(Form("hMCRecPhotNoPrim_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax )) ;
 413
 414
 415
 416          //MC tagging: reasons of partner loss etc.
 417          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnStack_%s_cent%d",cPID[iPID],cen),"Decay photons with partner not in Stack", nPt,0.,ptMax )) ;
 418          for(Int_t iType=0; iType<9; iType++)
 419            fOutputContainer->Add(new TH1F(Form("hMCDecWithFoundPartnType%d_%s_cent%d",iType,cPID[iPID],cen),"Decay photon with found partner", nPt,0.,ptMax )) ;
 420          fOutputContainer->Add(new TH1F(Form("hMCDecWithWrongMass_%s_cent%d",cPID[iPID],cen),"Decay photon with wrong mass", nPt,0.,ptMax )) ;
 421          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAccept_%s_cent%d",cPID[iPID],cen),"Decay photon with parttner not in PHOS acc", nPt,0.,ptMax )) ;
 422          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA1_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due geometry Fid. area. 1", nPt,0.,ptMax )) ;
 423          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA2_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due geometry Fid. area. 2", nPt,0.,ptMax)) ;
 424          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA3_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due geometry Fid. area. 3", nPt,0.,ptMax )) ;
 425          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnConv_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due to conversion", nPt,0.,ptMax )) ;
 426          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnEmin_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due to low energy", nPt,0.,ptMax )) ;
 427          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnOther_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due unknown reason", nPt,0.,ptMax )) ;
 428          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAll_%s_cent%d",cPID[iPID],cen),"Decay photons with partner missed due to any reason", nPt,0.,ptMax )) ;
 429          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnNPhot_%s_cent%d",cPID[iPID],cen),"pi0 decay photon with non-photon partner", nPt,0.,ptMax )) ;
 430
 431          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutEmin_%s_cent%d",cPID[iPID],cen),"Decay photons with rec. partner but failed Emin cut", nPt,0.,ptMax )) ;
 432          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutNcell_%s_cent%d",cPID[iPID],cen),"Decay photons with rec. partner but failed Ncell cut", nPt,0.,ptMax )) ;
 433          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutEcross_%s_cent%d",cPID[iPID],cen),"Decay photons with rec. partner but failed Ecross cut", nPt,0.,ptMax )) ;
 434          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutM02_%s_cent%d",cPID[iPID],cen),"Decay photons with rec. partner but failed M02 cut", nPt,0.,ptMax )) ;
 435          fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnDefCuts_%s_cent%d",cPID[iPID],cen),"Decay photons with rec. partner but failed default cuts", nPt,0.,ptMax )) ;
 436          fOutputContainer->Add(new TH1F(Form("hMCDecWRecPartn_%s_cent%d",cPID[iPID],cen),"Decay photons with rec partner", nPt,0.,ptMax )) ;
 437          fOutputContainer->Add(new TH1F(Form("hMCDecWRecUniqPartn_%s_cent%d",cPID[iPID],cen),"Decay photons with rec partner", nPt,0.,ptMax )) ;
 438
 439          fOutputContainer->Add(new TH1F(Form("hMCDecMerged_%s_cent%d",cPID[iPID],cen),"Decay photons with rec partner", nPt,0.,ptMax )) ;
 440          fOutputContainer->Add(new TH1F(Form("hMCDecUnfolded_%s_cent%d",cPID[iPID],cen),"Decay photons with rec partner", nPt,0.,ptMax )) ;
 441
 442
 443          fOutputContainer->Add(new TH2F(Form("hMC_InvM_Re_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 444          fOutputContainer->Add(new TH2F(Form("hMC_InvM_Re_Strict_%s_cent%d",cPID[iPID],cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 445        }
 446        fOutputContainer->Add(new TH2F(Form("hMCmass_cent%d",cen),"Mass with reconstructed decay partner",nM,0.,mMax,nPt,0.,ptMax )) ;
 447     }
 448
 449 /*
 450     fOutputContainer->Add(new TH1F(Form("hMCRecNoLabel_cent%d",cen),"Clusters without label",nPt,0.,ptMax)) ;
 451     fOutputContainer->Add(new TH1F(Form("hMCConversionRadius_cent%d",cen),"Clusters without label",600,0.,600.)) ;
 452     fOutputContainer->Add(new TH2F(Form("hMCRecPi0Vtx_cent%d",cen),"Secondary pi0s",100,0.,10.,600,0.,600.)) ;
 453     fOutputContainer->Add(new TH2F(Form("hMCRecEtaVtx_cent%d",cen),"Secondary etas",100,0.,10.,600,0.,600.)) ;
 454     fOutputContainer->Add(new TH2F(Form("hMCRecOmegaVtx_cent%d",cen),"Secondary etas",100,0.,10.,600,0.,600.)) ;
 455     fOutputContainer->Add(new TH2F(Form("hMCRecEtaprVtx_cent%d",cen),"Secondary etas",100,0.,10.,600,0.,600.)) ;
 456     fOutputContainer->Add(new TH2F(Form("hMCRecK0sVtx_cent%d",cen),"Secondary K0s",100,0.,10.,600,0.,600.)) ;
 457     fOutputContainer->Add(new TH2F(Form("hMCRecK0lVtx_cent%d",cen),"Secondary K0l",100,0.,10.,600,0.,600.)) ;
 458     fOutputContainer->Add(new TH2F(Form("hMCGammaPi0MisConvR_cent%d",cen),"Converted photons",400,0.,40.,600,0.,600.)) ;
 459
 460
 461     fOutputContainer->Add(new TH2F(Form("hMCGammaPi0PrimMgg_cent%d",cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 462     fOutputContainer->Add(new TH2F(Form("hMCGammaPi0RecMgg_cent%d",cen),"Two-photon inv. mass vs first photon pt",nM,0.,mMax,nPt,0.,ptMax)) ;
 463
 464     for(Int_t iPID=0; iPID<4; iPID++){
 465       fOutputContainer->Add(new TH1F(Form("hMCRecNoPrim_%s_cent%d",cPID[iPID],cen),"Reconstructed particles withour primary",nPt,0.,ptMax)) ;
 466       fOutputContainer->Add(new TH1F(Form("hMCRecGamma_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: Gamma",nPt,0.,ptMax)) ;
 467       fOutputContainer->Add(new TH1F(Form("hMCRecPhotConv_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: e+-",nPt,0.,ptMax)) ;
 468       fOutputContainer->Add(new TH1F(Form("hMCRecPi0_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: pi0",nPt,0.,ptMax)) ;
 469       fOutputContainer->Add(new TH1F(Form("hMCRecEta_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: eta",nPt,0.,ptMax)) ;
 470       fOutputContainer->Add(new TH1F(Form("hMCRecOmega_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: Gamma",nPt,0.,ptMax)) ;
 471       fOutputContainer->Add(new TH1F(Form("hMCRecEtapr_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: eta prime",nPt,0.,ptMax)) ;
 472       fOutputContainer->Add(new TH1F(Form("hMCRecPbar_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: bar(p)",nPt,0.,ptMax)) ;
 473       fOutputContainer->Add(new TH1F(Form("hMCRecNbar_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: bar(n)",nPt,0.,ptMax)) ;
 474       fOutputContainer->Add(new TH1F(Form("hMCRecPipm_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: pipm",nPt,0.,ptMax)) ;
 475       fOutputContainer->Add(new TH1F(Form("hMCRecN_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: n",nPt,0.,ptMax)) ;
 476       fOutputContainer->Add(new TH1F(Form("hMCRecP_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: p",nPt,0.,ptMax)) ;
 477       fOutputContainer->Add(new TH1F(Form("hMCRecKpm_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: K+-",nPt,0.,ptMax)) ;
 478       fOutputContainer->Add(new TH1F(Form("hMCRecK0s_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: K0s",nPt,0.,ptMax)) ;
 479       fOutputContainer->Add(new TH1F(Form("hMCRecK0l_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: K0l",nPt,0.,ptMax)) ;
 480       fOutputContainer->Add(new TH1F(Form("hMCRecUnknownCh_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: K0l",nPt,0.,ptMax)) ;
 481       fOutputContainer->Add(new TH1F(Form("hMCRecUnknownNeu_%s_cent%d",cPID[iPID],cen),"Reconstructed particles with primary: K0l",nPt,0.,ptMax)) ;
 482
 483       //Decay photons
 484       fOutputContainer->Add(new TH1F(Form("hMCRecGammaDir_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, no primary",nPt,0.,ptMax)) ;
 485       fOutputContainer->Add(new TH1F(Form("hMCRecGammaPi0_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 486       fOutputContainer->Add(new TH1F(Form("hMCRecGammaEta_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from eta",nPt,0.,ptMax)) ;
 487       fOutputContainer->Add(new TH1F(Form("hMCRecGammaOmega_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from omega",nPt,0.,ptMax)) ;
 488       fOutputContainer->Add(new TH1F(Form("hMCRecGammaOther_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from other",nPt,0.,ptMax)) ;
 489       fOutputContainer->Add(new TH1F(Form("hMCRecPhotNoPrim_%s_cent%d",cPID[iPID],cen),"Spectrum of rec. photons", nPt,0.,ptMax)) ;
 490
 491       //Pi0 decay photons
 492
 493       //MC tagging: reasons of partner loss etc.
 494        fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnStack_cent%d",cen),"Decay photons with partner not in Stack", nPt,0.,ptMax)) ;
 495        fOutputContainer->Add(new TH1F(Form("hMCDecWithFoundPartn_cent%d",cen),"Decay photon with found partner", nPt,0.,ptMax)) ;
 496        fOutputContainer->Add(new TH1F(Form("hMCDecWithWrongMass_cent%d",cen),"Decay photon with wrong mass", nPt,0.,ptMax)) ;
 497        fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAccept_cent%d",cen),"Decay photon with parttner not in PHOS acc", nPt,0.,ptMax)) ;
 498        fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA1_cent%d",cen),"Decay photons with partner missed due geometry Fid. area. 1", nPt,0.,ptMax)) ;
 499        fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA2_cent%d",cen),"Decay photons with partner missed due geometry Fid. area. 2", nPt,0.,ptMax)) ;
 500        fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAcceptFA3_cent%d",cen),"Decay photons with partner missed due geometry Fid. area. 3", nPt,0.,ptMax)) ;
 501        fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnConv_cent%d",cen),"Decay photons with partner missed due to conversion", nPt,0.,ptMax)) ;
 502        fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnEmin_cent%d",cen),"Decay photons with partner missed due to low energy", nPt,0.,ptMax)) ;
 503        fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnOther_cent%d",cen),"Decay photons with partner missed due unknown reason", nPt,0.,ptMax)) ;
 504        fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnAll_cent%d",cen),"Decay photons with partner missed due to any reason", nPt,0.,ptMax)) ;
 505        fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnNPhot_cent%d",cen),"pi0 decay photon with non-photon partner", nPt,0.,ptMax)) ;
 506
 507
 508        fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutEmin_cent%d",cen),"Decay photons with rec. partner but failed Emin cut", nPt,0.,ptMax)) ;
 509        fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutNcell_cent%d",cen),"Decay photons with rec. partner but failed Ncell cut", nPt,0.,ptMax)) ;
 510        fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutEcross_cent%d",cen),"Decay photons with rec. partner but failed Ecross cut", nPt,0.,ptMax)) ;
 511        fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnCutM02_cent%d",cen),"Decay photons with rec. partner but failed M02 cut", nPt,0.,ptMax)) ;
 512        fOutputContainer->Add(new TH1F(Form("hMCDecWMisPartnDefCuts_cent%d",cen),"Decay photons with rec. partner but failed default cuts", nPt,0.,ptMax)) ;
 513        fOutputContainer->Add(new TH1F(Form("hMCDecWRecPartn_cent%d",cen),"Decay photons with rec partner", nPt,0.,ptMax)) ;
 514
 515        fOutputContainer->Add(new TH2F(Form("hMCmass_cent%d",cen),"Mass with reconstructed decay partner",nM,0.,mMax,nPt,0.,ptMax)) ;
 516
 517
 518       fOutputContainer->Add(new TH1F(Form("hMCRecGammaPi0Dalitz_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 519       fOutputContainer->Add(new TH1F(Form("hMCRecGammaPi0NoStack_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 520       fOutputContainer->Add(new TH1F(Form("hMCGammaPi02Gamma_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 521       fOutputContainer->Add(new TH1F(Form("hMCGammaPi0Rec_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 522       fOutputContainer->Add(new TH1F(Form("hMCGammaPi0RecSoft_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 523       fOutputContainer->Add(new TH1F(Form("hMCGammaPi0RecCells_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 524
 525     fOutputContainer->Add(new TH1F(Form("hMCGammaPi0Tagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 526     fOutputContainer->Add(new TH1F(Form("hMCGammaPi0FakeTagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 527     fOutputContainer->Add(new TH1F(Form("hMCGammaPi0TrueTagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 528     fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MultyTagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 529     fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisGeo_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 530     fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisGeoFA1_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 531     fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisGeoFA2_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 532     fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisGeoFA3_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 533     fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisConv_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 534     fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisEmin_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 535     fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisOther_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 536     fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisFakePrim_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 537     fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisFoundPrim_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 538     fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisNPhot_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 539
 540     //all clusters fake tagged
 541     fOutputContainer->Add(new TH1F(Form("hMCAllFakeTagged_%s_cent%d",cPID[iPID],cen),"Reconstructed gammas, from pi0",nPt,0.,ptMax)) ;
 542
 543
 544   }
 545   fOutputContainer->Add(new TH1F(Form("hMCGammaPi0MisPartner_cent%d",cen),"Spectrum of missed partners",nPt,0.,ptMax)) ;
 546   fOutputContainer->Add(new TH2F(Form("hMCGammaPi0MisPartnerEtaPhi_cent%d",cen),"Spectrum of missed partners",100,-0.2,0.2,100,4.5,5.6)) ;
 547     }
 548 */
 549   }
 550
 551   //If we work with MC, need to set Sumw2 - we will use weights
 552   if(fIsMC){
 553       for(Int_t i=0; i<fOutputContainer->GetSize();i++){
 554         ((TH1*)fOutputContainer->At(i))->Sumw2() ;
 555       }
 556   }
 557
 558
 559
 560   for(Int_t i=0;i<10;i++)
 561     for(Int_t j=0;j<5;j++)
 562       fPHOSEvents[i][j]=0x0 ;    //Container for PHOS photons
 563
 564
 565   //Prepare PHOS trigger utils if necessary
 566   if(!fIsMB)
 567    fPHOSTrigUtils = new AliPHOSTriggerUtils("PHOSTrig") ;
 568
 569
 570   PostData(1, fOutputContainer);
 571
 572
 573 }
 574
 575 //______________________________________________________________________________
 576 void AliAnalysisTaskTaggedPhotons::UserExec(Option_t *)
 577 {
 578   //Select events
 579   //Select photons
 580   //Fill QA histograms
 581   //Fill Tagging histogsms
 582
 583   const Double_t kEcrossCut=0.98 ;
 584   const Double_t kTOFMaxCut= 100.e-9 ;
 585   const Double_t kTOFMinCut=-100.e-9 ;
 586
 587   // Event selection flags
 588   //  FillHistogram("hSelEvents",0) ;
 589
 590   AliVEvent* event = (AliVEvent*)InputEvent();
 591   if(!event){
 592     AliDebug(1,"No event") ;
 593     PostData(1, fOutputContainer);
 594     return;
 595   }
 596   FillHistogram("hSelEvents",1) ;
 597
 598   //MC stack init
 599   fStack = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
 600
 601   //read geometry if not read yet
 602   if(fPHOSgeom==0){
 603     InitGeometry() ;
 604   }
 605
 606   if(!fIsFastMC){
 607     if(!fUtils)
 608       fUtils = new AliAnalysisUtils();
 609
 610     Bool_t isMB = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kINT7)  ;
 611     Bool_t isPHI7 = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kPHI7);
 612
 613     if((fIsMB && !isMB) || (!fIsMB && !isPHI7)){
 614       PostData(1, fOutputContainer);
 615       return;
 616     }
 617   }
 618
 619   //If we work with PHOS trigger, prepapre TriggerUtils
 620   if(!fIsMB)
 621    fPHOSTrigUtils->SetEvent(event) ;
 622
 623   FillHistogram("hSelEvents",2) ;
 624
 625   // Checks if we have a primary vertex
 626   // Get primary vertices form AOD
 627
 628   Double_t vtx5[3];
 629   vtx5[0] = event->GetPrimaryVertex()->GetX();
 630   vtx5[1] = event->GetPrimaryVertex()->GetY();
 631   vtx5[2] = event->GetPrimaryVertex()->GetZ();
 632
 633   FillHistogram("hNvertexTracks",event->GetPrimaryVertex()->GetNContributors());
 634   FillHistogram("hZvertex"      ,vtx5[2]);
 635   if(fIsFastMC){ //vertex from header
 636     AliAODMCHeader *cHeaderAOD = dynamic_cast<AliAODMCHeader*>(event->FindListObject(AliAODMCHeader::StdBranchName()));
 637     if(!cHeaderAOD){
 638       PostData(1, fOutputContainer);
 639       return ;
 640     }
 641     cHeaderAOD->GetVertex(vtx5);
 642   }
 643   if (TMath::Abs(vtx5[2]) > 10. ){
 644     PostData(1, fOutputContainer);
 645     return ;
 646   }
 647   FillHistogram("hSelEvents",3) ;
 648   //Vtx class z-bin
 649   Int_t zvtx = TMath::Min(9,Int_t((vtx5[2]+10.)/2.)) ;
 650
 651
 652   if(!fIsFastMC){
 653 //    if (event->IsPileupFromSPD()){
 654 //      PostData(1, fOutputContainer);
 655 //      return ;
 656 //    }
 657
 658     if(!fUtils->IsVertexSelected2013pA(event)){
 659       PostData(1, fOutputContainer);
 660       return ;
 661     }
 662
 663     FillHistogram("hSelEvents",4) ;
 664
 665     if(fUtils->IsPileUpEvent(event)){
 666       PostData(1, fOutputContainer);
 667       return ;
 668     }
 669     FillHistogram("hSelEvents",5) ;
 670   }
 671
 672   //centrality
 673   if(!fIsFastMC){
 674     AliCentrality *centrality = event->GetCentrality();
 675     if( centrality )
 676       fCentrality=centrality->GetCentralityPercentile("V0M");
 677     else {
 678       AliError("Event has 0x0 centrality");
 679       fCentrality = -1.;
 680     }
 681     FillHistogram("hCentrality",fCentrality) ;
 682
 683     if(fCentrality<0. || fCentrality>=100.){
 684       PostData(1, fOutputContainer);
 685       return ;
 686     }
 687   }
 688   else{
 689     fCentrality=1.;
 690   }
 691   fCentBin = (Int_t)(fCentrality/20.) ;
 692
 693   FillHistogram("hSelEvents",6) ;
 694
 695
 696   //Calculate charged multiplicity
 697   Int_t trackMult = 0;
 698   if(fTrackEvent)
 699     fTrackEvent->Clear() ;
 700   else
 701     fTrackEvent = new TClonesArray("AliCaloPhoton",event->GetNumberOfTracks()) ;
 702
 703   for (Int_t i=0;i<event->GetNumberOfTracks();++i) {
 704     AliAODTrack *track = (AliAODTrack*)event->GetTrack(i) ;
 705     if(!track->IsHybridGlobalConstrainedGlobal())
 706       continue ;
 707     if(TMath::Abs(track->Eta())< 0.9){
 708       if(trackMult>=fTrackEvent->GetSize())
 709         fTrackEvent->Expand(2*trackMult) ;
 710       new ((*fTrackEvent)[trackMult]) AliCaloPhoton(track->Px(),track->Py(),track->Pz(),track->P());
 711       trackMult++;
 712      FillHistogram("hTrackEtaPhi",track->Eta(),track->Phi()) ;
 713      FillHistogram("hTrackEtaPt",track->Eta(),track->Pt()) ;
 714     }
 715   }
 716   FillHistogram("hTrackMult",fCentrality,trackMult+0.5) ;
 717
 718   if(!fPHOSEvents[zvtx][fCentBin])
 719     fPHOSEvents[zvtx][fCentBin]=new TList() ;
 720   fCurrentMixedList = fPHOSEvents[zvtx][fCentBin] ;
 721
 722    const Double_t rcut=1. ; //cut on vertex to consider particle as "primary"
 723
 724   //---------Select photons-------------------
 725   Int_t multClust = event->GetNumberOfCaloClusters();
 726   if(!fPHOSEvent)
 727     fPHOSEvent   = new TClonesArray("AliCaloPhoton",multClust);
 728   else
 729     fPHOSEvent->Clear() ;
 730   Int_t inList = 0; //counter of caloClusters
 731
 732   for (Int_t i=0; i<multClust; i++) {
 733     AliVCluster * clu = event->GetCaloCluster(i);
 734
 735     if(!clu->IsPHOS())
 736       continue ;
 737
 738     if(clu->E()<0.1)
 739       continue;
 740
 741     if(clu->GetNCells()<3)
 742       continue ;
 743
 744     if(clu->GetM02()<0.2)
 745       continue ;
 746
 747     if(clu->GetMCEnergyFraction()>kEcrossCut) //Ecross cut, should be filled with Tender
 748      continue ;
 749
 750     if(clu->GetDistanceToBadChannel()<fMinBCDistance)
 751       continue ;
 752
 753     Float_t pos[3] ;
 754     clu->GetPosition(pos) ;
 755     Int_t fidArea=GetFiducialArea(pos) ;
 756 //    if(fidArea==0) //Bad cell
 757 //      continue;
 758
 759     TVector3 global1(pos) ;
 760     Int_t relId[4] ;
 761     fPHOSgeom->GlobalPos2RelId(global1,relId) ;
 762     Int_t mod  = relId[0] ;
 763     Int_t cellX = relId[2];
 764     Int_t cellZ = relId[3] ;
 765
 766     FillHistogram("hTOF",clu->E(),clu->GetTOF()) ;
 767     FillHistogram(Form("hTofM%d",mod),clu->GetTOF(),clu->E()) ;
 768     if((!fIsMC) && (clu->GetTOF() < kTOFMinCut || clu->GetTOF() > kTOFMaxCut))
 769       continue ;
 770
 771
 772
 773     TLorentzVector momentum ;
 774     clu->GetMomentum(momentum, vtx5);
 775
 776     AliCaloPhoton *p = new ((*fPHOSEvent)[inList]) AliCaloPhoton(momentum.Px(),momentum.Py(),momentum.Pz(),clu->E() );
 777 //    momentum*= clu->GetCoreEnergy()/momentum.E() ;
 778 //    AliCaloPhoton *p = new ((*fPHOSEvent)[inList]) AliCaloPhoton(momentum.Px(),momentum.Py(),momentum.Pz(),clu->GetCoreEnergy() );
 779     inList++;
 780
 781     Int_t isolation = EvalIsolation(&momentum,kTRUE) ;
 782     p->SetIsolationTag(isolation) ;
 783
 784     p->SetDistToBad((Int_t)(1.+clu->GetDistanceToBadChannel()/2.2));
 785     p->SetBC(i) ; //reference to CaloCluster
 786     p->SetTagInfo(0); //No pi0 partners found so far
 787     p->SetTagged(kFALSE);   //Reconstructed pairs found
 788
 789     p->SetFiducialArea(fidArea) ;
 790
 791     //Mark photons fired trigger
 792     if(!fIsMB)
 793       p->SetTrig(fPHOSTrigUtils->IsFiredTrigger(clu)) ;
 794
 795     if(fIsMB || !fIsMB && p->IsTrig()){
 796       FillHistogram(Form("hCluNXZM%d",mod),cellX,cellZ,1.);
 797       FillHistogram(Form("hCluEXZM%d",mod),cellX,cellZ,clu->E());
 798       if(fidArea>1){
 799         FillHistogram(Form("hCluArea2M%d",mod),cellX,cellZ,1.);
 800         if(fidArea>2){
 801           FillHistogram(Form("hCluArea3M%d",mod),cellX,cellZ,1.);
 802         }
 803       }
 804     }
 805
 806     if(fIsMC){
 807        //Look for MC particle entered PHOS
 808        FillHistogram(Form("LabelsNPrim_cent%d",fCentBin),clu->E(),float(clu->GetNLabels())) ;
 809        Int_t primLabel=clu->GetLabelAt(0) ; //FindPrimary(clu,sure) ;
 810        //Look what particle left vertex
 811        if(primLabel>-1){
 812          AliAODMCParticle * prim = (AliAODMCParticle*)fStack->At(primLabel) ;
 813          Int_t iparent=primLabel;
 814          AliAODMCParticle * parent = prim;
 815          Double_t r2=prim->Xv()*prim->Xv()+prim->Yv()*prim->Yv() ;
 816          while((r2 > rcut*rcut) && (iparent>-1)){
 817            iparent=parent->GetMother();
 818            parent=(AliAODMCParticle*)fStack->At(iparent);
 819            r2=parent->Xv()*parent->Xv()+parent->Yv()*parent->Yv() ;
 820          }
 821          p->SetPrimary(primLabel) ;
 822          p->SetPrimaryAtVertex(iparent) ;
 823          p->SetWeight(PrimaryParticleWeight(parent)) ;
 824        }
 825        else{
 826          p->SetPrimary(-1); //Primary index
 827          p->SetPrimaryAtVertex(-1) ;
 828          p->SetWeight(1.) ;
 829        }
 830     }
 831     else{
 832       p->SetPrimary(-1); //Primary index
 833       p->SetPrimaryAtVertex(-1) ;
 834       p->SetWeight(1.) ;
 835     }
 836     //PID criteria
 837 //  Cut on Core Lambdas
 838     p->SetDispBit(clu->Chi2()<2.5*2.5) ;
 839 //  Cut on FullLamdas
 840 //    p->SetDispBit(clu->GetDispersion()<2.5*2.5) ;
 841     p->SetTOFBit(TestTOF(clu->GetTOF(),clu->E())) ;
 842     p->SetCPVBit(clu->GetEmcCpvDistance()>2.5) ;
 843
 844     FillHistogram(Form("hPhotM%d_All",mod),p->Pt()) ;
 845     if(fidArea>1){
 846      FillHistogram(Form("hPhotM%dA2_All",mod),p->Pt()) ;
 847       if(fidArea>2){
 848         FillHistogram(Form("hPhotM%dA3_All",mod),p->Pt()) ;
 849       }
 850     }
 851     if(p->IsDispOK()){
 852       FillHistogram(Form("hPhotM%d_Disp",mod),p->Pt()) ;
 853       if(fidArea>1){
 854        FillHistogram(Form("hPhotM%dA2_Disp",mod),p->Pt()) ;
 855         if(fidArea>2){
 856           FillHistogram(Form("hPhotM%dA3_Disp",mod),p->Pt()) ;
 857         }
 858       }
 859       if(p->IsCPVOK()){
 860         FillHistogram(Form("hPhotM%d_Both",mod),p->Pt()) ;
 861         if(fidArea>1){
 862          FillHistogram(Form("hPhotM%dA2_Both",mod),p->Pt()) ;
 863           if(fidArea>2){
 864             FillHistogram(Form("hPhotM%dA3_Both",mod),p->Pt()) ;
 865           }
 866         }
 867       }
 868     }
 869     if(p->IsCPVOK()){
 870       FillHistogram(Form("hPhotM%d_CPV",mod),p->Pt()) ;
 871       if(fidArea>1){
 872        FillHistogram(Form("hPhotM%dA2_CPV",mod),p->Pt()) ;
 873         if(fidArea>2){
 874           FillHistogram(Form("hPhotM%dA3_CPV",mod),p->Pt()) ;
 875         }
 876       }
 877     }
 878   }
 879   FillHistogram("hPHOSCentrality",fCentrality,inList+0.5) ;
 880
 881
 882   FillMCHistos() ;
 883   FillTaggingHistos() ;
 884
 885   //Remove old events
 886   fCurrentMixedList->AddFirst(fPHOSEvent);
 887   fPHOSEvent=0x0 ;
 888   if(fCurrentMixedList->GetSize() > 10){
 889     TClonesArray *tmp = static_cast <TClonesArray*> (fCurrentMixedList->Last());
 890     fCurrentMixedList->Remove(tmp);
 891     delete tmp;
 892   }
 893
 894   PostData(1, fOutputContainer);
 895
 896 }
 897 //________________________________________________
 898 void AliAnalysisTaskTaggedPhotons::FillMCHistos(){
 899
 900   //MC info about this particle
 901   if(!fIsMC)
 902     return ;
 903   const Double_t rcut=1. ; //cut on vertex to consider particle as "primary"
 904   const Double_t phiMin=260.*TMath::Pi()/180. ;
 905   const Double_t phiMax=320.*TMath::Pi()/180. ;
 906
 907   AliVEvent* event = (AliVEvent*)InputEvent();
 908
 909   Int_t nPrim = fStack->GetEntriesFast() ;
 910   //Fill Primary particl yields
 911
 912   for(Int_t i=0;i<nPrim;i++){
 913     AliAODMCParticle * prim = (AliAODMCParticle*)fStack->At(i) ;
 914     Double_t r2=prim->Xv()*prim->Xv()+prim->Yv()*prim->Yv() ;
 915     if(r2>rcut*rcut)
 916       continue ;
 917
 918     Int_t pdg=prim->GetPdgCode() ;
 919     char partName[30] ;
 920     if(pdg == 111)
 921       snprintf(partName,30,"pi0") ;
 922     else
 923       if(pdg == 221)
 924         snprintf(partName,30,"eta") ;
 925       else
 926         if(pdg == 22)
 927            snprintf(partName,30,"gamma") ;
 928         else
 929           if(pdg == 310)
 930              snprintf(partName,30,"K0s") ;
 931           else
 932             if(abs(pdg) == 321)
 933               snprintf(partName,30,"Kpm") ;
 934             else
 935               if(abs(pdg) == 211)
 936                 snprintf(partName,30,"pipm") ;
 937               else
 938                 if(abs(pdg) == 2212)
 939                   snprintf(partName,30,"p") ;
 940                 else
 941                   if(abs(pdg) ==-2212)
 942                     snprintf(partName,30,"pbar") ;
 943                   else
 944                     if(abs(pdg) == 2112)
 945                       snprintf(partName,30,"n") ;
 946                     else
 947                       if(abs(pdg) ==-2112)
 948                         snprintf(partName,30,"nbar") ;
 949                       else
 950                         continue ;
 951
 952     //Primary particle
 953     Double_t phi=prim->Phi() ;
 954     while(phi<0.)phi+=TMath::TwoPi() ;
 955     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 956     Double_t pt=prim->Pt() ;
 957
 958     //Total number of pi0 with creation radius <1 cm
 959     Double_t w = PrimaryParticleWeight(prim) ;
 960     FillHistogram(Form("hMC_all_%s_cent%d",partName,fCentBin),pt,w) ;
 961     if(TMath::Abs(prim->Y())<0.13){
 962       FillHistogram(Form("hMC_phi_%s_cent%d",partName,fCentBin),phi,w) ;
 963       if(phi>phiMin && phi<phiMax){
 964         FillHistogram(Form("hMC_unitEta_%s_cent%d",partName,fCentBin),pt,w) ;
 965         if(prim->GetMother()>-1){
 966           AliAODMCParticle * primGM = (AliAODMCParticle*)fStack->At(prim->GetMother()) ;
 967           if(primGM->GetPdgCode()==111)
 968             FillHistogram(Form("hMC_unitEta_GammaPi0_cent%d",fCentBin),pt,w) ;
 969         }
 970       }
 971     }
 972
 973     FillHistogram(Form("hMC_rap_%s_cent%d",partName,fCentBin),prim->Y(),w) ;
 974     //Some additional QA
 975     if(pdg == 111){
 976        FillHistogram("hMCpi0_ptrap",pt,prim->Y(),w) ;
 977        FillHistogram("hMCpi0_ptphi",pt,phi,w) ;
 978     }
 979     if(pdg == 22){
 980        FillHistogram("hMCgamma_ptrap",pt,prim->Y(),w) ;
 981        FillHistogram("hMCgamma_ptphi",pt,phi,w) ;
 982     }
 983
 984   }
 985
 986
 987
 988   //Clussify reconstructed clusters
 989   //First - photons (from vertex) and contaminations
 990   //Second - photons from different sources
 991   //Third - photons from pi0s - missed for different reasons
 992
 993   const Int_t n=fPHOSEvent->GetEntriesFast() ;
 994   for(Int_t i=0;i<n;i++){
 995     AliCaloPhoton *p = static_cast<AliCaloPhoton*>(fPHOSEvent->At(i));
 996     Int_t label=p->GetPrimary() ;
 997     if(label<0){ //No label!
 998       FillHistogram("hMCRecNoLabel",p->Pt(),p->GetWeight());
 999       continue ;
1000     }
1001
1002
1003     AliAODMCParticle * prim = (AliAODMCParticle*)fStack->At(p->GetPrimary()) ;
1004     //Look what particle left virtex
1005     Int_t iparent=p->GetPrimary();
1006     AliAODMCParticle * parent = prim;
1007     while(parent->Xv()*parent->Xv()+parent->Yv()*parent->Yv() > rcut*rcut){
1008         iparent=parent->GetMother();
1009         if(iparent<0)
1010           break ;
1011         parent = (AliAODMCParticle*)fStack->At(iparent) ;
1012       }
1013       Int_t parentPDG=parent->GetPdgCode() ;
1014       Int_t iFidArea=p->GetFiducialArea();
1015       switch(parentPDG){
1016         case 22: //electron/positron conversion
1017           FillPIDHistograms("hMCRecPhoton",p);  //Reconstructed with photon from conversion primary
1018           if(iFidArea>0){
1019             FillPIDHistograms("hMCRecPhoton_Area1",p);
1020             if(iFidArea>1){
1021               FillPIDHistograms("hMCRecPhoton_Area2",p);
1022               if(iFidArea>1){
1023                 FillPIDHistograms("hMCRecPhoton_Area3",p);
1024               }
1025             }
1026           }
1027           break ;
1028         case  11:
1029         case -11: //electron/positron conversion
1030           FillPIDHistograms("hMCRecE",p);  //Reconstructed with photon from conversion primary
1031           break ;
1032         case -2212:
1033           FillPIDHistograms("hMCRecPbar",p);  //Reconstructed with photon from antibaryon annihilation
1034           break ;
1035         case -2112: //antineutron & antiproton conversion
1036           FillPIDHistograms("hMCRecNbar",p);  //Reconstructed with photon from antibaryon annihilation
1037           break ;
1038         case  211:
1039         case -211:
1040           FillPIDHistograms("hMCRecPipm",p);  //Reconstructed with photon from antibaryon annihilation
1041           break ;
1042         case 2212:
1043           FillPIDHistograms("hMCRecP",p);  //Reconstructed with photon from antibaryon annihilation
1044           break ;
1045         case  321:
1046         case -321:
1047           FillPIDHistograms("hMCRecKpm",p);  //Reconstructed with photon from conversion primary
1048           break ;
1049         case 310:
1050           FillPIDHistograms("hMCRecK0s",p);  //Reconstructed with photon from conversion primary
1051           break ;
1052         case 2112: //antineutron & antiproton conversion
1053           FillPIDHistograms("hMCRecN",p);  //Reconstructed with photon from antibaryon annihilation
1054           break ;
1055         case -1: //direct photon or no primary
1056           FillPIDHistograms("hMCRecNoPRim",p);
1057           break ;
1058         default:
1059           if(parent->Charge()!=0)
1060             FillPIDHistograms("hMCRecCharg",p);  //Reconstructed with photon from antibaryon annihilation
1061           else
1062             FillPIDHistograms("hMCRecNeutral",p);  //Reconstructed with photon from antibaryon annihilation
1063       }
1064
1065
1066       //Now classify decay photon
1067       if(parentPDG==22){
1068         Int_t iGrandParent=parent->GetMother();
1069         if(iGrandParent<0 || iGrandParent>=fStack->GetEntriesFast()){
1070           FillPIDHistograms("hMCRecPhotNoPrim",p);
1071           continue ;
1072         }
1073         AliAODMCParticle * grandParent = (AliAODMCParticle*)fStack->At(iGrandParent) ;
1074         Int_t grandParentPDG=grandParent->GetPdgCode() ;
1075         switch(grandParentPDG){
1076         case 111: //pi0
1077           FillPIDHistograms("hMCRecPhotPi0",p);
1078           break ;
1079         case 221: //eta decay
1080           FillPIDHistograms("hMCRecPhotEta",p);
1081           break ;
1082         case 223: //omega meson decay
1083           FillPIDHistograms("hMCRecPhotOmega",p);
1084           break ;
1085         default:
1086           FillPIDHistograms("hMCRecPhotOther",p);
1087         }
1088         //--------consider pi0 decays--------------------
1089         if(grandParentPDG==111){
1090           //First find which daughter is our cluster
1091           //iparent - index of curent photon
1092           Int_t ipartner=grandParent->GetDaughter(0) ;
1093           if(ipartner==iparent){//look for other
1094             if(grandParent->GetNDaughters()>1){
1095               ipartner=grandParent->GetDaughter(1);
1096             }
1097             else{
1098               ipartner=-1 ;
1099             }
1100           }
1101           //There is no partner in stack
1102           if(ipartner==-1){
1103             FillPIDHistograms("hMCDecWMisPartnStack",p) ;
1104           }
1105           else{
1106             AliAODMCParticle * partner = (AliAODMCParticle *)fStack->At(ipartner);
1107             //Check if partner is registered and made correct mass
1108             //If not - trace the reason
1109             AliCaloPhoton *pp = 0x0 ;
1110
1111             for(Int_t ii=0;ii<n;ii++){
1112               if(i==ii) continue;
1113               AliCaloPhoton * tmp = static_cast<AliCaloPhoton*>(fPHOSEvent->At(ii));
1114               Int_t ipartnPrim = tmp->GetPrimary() ;
1115               while(ipartnPrim>-1){
1116                 if(ipartnPrim==ipartner){
1117                   break ;
1118                 }
1119                 ipartnPrim = ((AliAODMCParticle *)fStack->At(ipartnPrim))->GetMother();
1120               }
1121               if(ipartnPrim==ipartner){
1122                 pp=tmp ;
1123                 break ;
1124               }
1125             }
1126             if(pp){
1127               //Partner reconstructed, but did not pass cuts
1128                 FillPIDHistograms("hMCDecWRecUniqPartn",p) ;
1129             }
1130             //Partner not found. Check if it is not dominant contributor?
1131             if(!pp){
1132               for(Int_t ii=0;(ii<n) && (!pp);ii++){
1133                 if(i==ii) continue;
1134                 AliCaloPhoton * tmp = static_cast<AliCaloPhoton*>(fPHOSEvent->At(ii));
1135                 Int_t iCaloCluster=tmp->GetBC();//index of AODCaloCluster
1136                 AliVCluster * clu = event->GetCaloCluster(iCaloCluster);
1137                 Int_t nCluPrimaries = clu->GetNLabels() ;
1138                 for(Int_t iAODLabel=0; (iAODLabel<nCluPrimaries) && (!pp); iAODLabel++){
1139                   Int_t ipartnPrim = clu->GetLabelAt(iAODLabel) ;
1140                   while(ipartnPrim>-1){
1141                     if(ipartnPrim==ipartner){
1142                       break ;
1143                     }
1144                     ipartnPrim = ((AliAODMCParticle *)fStack->At(ipartnPrim))->GetMother();
1145                   }
1146                   if(ipartnPrim==ipartner){
1147                     pp=tmp ;
1148                     break ;
1149                   }
1150                 }
1151               }
1152             }
1153             //If partner still not found, check if it is in same cluster
1154             if(!pp){
1155               Int_t iCaloCluster=p->GetBC();//index of AODCaloCluster
1156               AliVCluster * clu = event->GetCaloCluster(iCaloCluster);
1157               Int_t nCluPrimaries = clu->GetNLabels() ;
1158               for(Int_t iAODLabel=0; iAODLabel<nCluPrimaries ; iAODLabel++){
1159                 Int_t ipartnPrim = clu->GetLabelAt(iAODLabel) ;
1160                 while(ipartnPrim>-1){
1161                   if(ipartnPrim==ipartner){
1162                       break ;
1163                   }
1164                   ipartnPrim = ((AliAODMCParticle *)fStack->At(ipartnPrim))->GetMother();
1165                 }
1166                 if(ipartnPrim==ipartner){ //yes, this cluster contains both primary
1167                   if(clu->GetNExMax()<2){ //was not unfolded
1168                     FillPIDHistograms("hMCDecMerged",p) ;
1169                   }
1170                   else{
1171                     FillPIDHistograms("hMCDecUnfolded",p) ;
1172                   }
1173                 }
1174               }
1175             }
1176
1177             if(pp){
1178               //Partner reconstructed, but did not pass cuts
1179                 FillPIDHistograms("hMCDecWRecPartn",p) ;
1180                 Double_t invMass=(*p+ *pp).M() ;
1181                 FillHistogram(Form("hMCmass_cent%d",fCentBin),invMass,p->Pt(),p->GetWeight()) ;
1182                 Double_t nSigma=InPi0Band(invMass,p->Pt()) ;
1183                 // analog to Tag
1184                 for(Int_t eminType=0; eminType<3; eminType++){
1185                   if(pp->E()>0.1*(eminType+1)){
1186                     for(Int_t isigma=0; isigma<3; isigma++){
1187                       if(nSigma<1.+isigma){
1188                          Int_t iType=3*eminType+isigma ;
1189                          FillPIDHistograms(Form("hMCDecWithFoundPartnType%d",iType),p) ;
1190                       }
1191                     }
1192                   }
1193                 }
1194                 if(nSigma>3.){
1195                   FillPIDHistograms("hMCDecWithWrongMass",p) ;
1196                 }
1197             }
1198             else{//Partner not reconstructed
1199               if(partner->GetPdgCode()==22){
1200                 Bool_t isPartnerLost=kFALSE; //If partner is lost for some reason
1201
1202                 //Check if partner miss PHOS
1203                 Int_t modulenum ;
1204                 Double_t ztmp=0.,xtmp=0. ;
1205                 Double_t vtx[3]={partner->Xv(),partner->Yv(),partner->Zv()} ;
1206                 Bool_t impact=fPHOSgeom->ImpactOnEmc(vtx,partner->Theta(),partner->Phi(),modulenum,ztmp,xtmp) ;
1207
1208                 if(impact){//still check bad map
1209                   Int_t relId[4] ;
1210                   fPHOSgeom->RelPosToRelId(modulenum,xtmp,ztmp,relId) ;
1211                   if ( !IsGoodChannel(modulenum,relId[2],relId[3]) ) {
1212                      impact=kFALSE ;
1213                   }
1214                 }
1215
1216                 if(!impact){ //this photon cannot hit PHOS
1217                   FillPIDHistograms("hMCDecWMisPartnAccept",p) ;  //Spectrum of tagged with missed partner
1218                   if(iFidArea>0){
1219                     FillPIDHistograms("hMCDecWMisPartnAcceptFA1",p) ;  //Spectrum of tagged with missed partner
1220                     if(iFidArea>1){
1221                       FillPIDHistograms("hMCDecWMisPartnAcceptFA2",p) ;  //Spectrum of tagged with missed partner
1222                       if(iFidArea>2){
1223                         FillPIDHistograms("hMCDecWMisPartnAcceptFA3",p) ;  //Spectrum of tagged with missed partner
1224                       }
1225                     }
1226                   }
1227                   isPartnerLost=kTRUE;
1228                 }
1229
1230                 if(!isPartnerLost){
1231                   //this photon is converted before it is registered
1232                   if(partner->GetNDaughters()>0){
1233                     AliAODMCParticle* tmpP=(AliAODMCParticle*)fStack->At(partner->GetDaughter(0));
1234                     if(tmpP->Xv()*tmpP->Xv()+tmpP->Yv()*tmpP->Yv()<450.*450.){
1235                       FillPIDHistograms("hMCDecWMisPartnConv",p) ;  //Spectrum of tagged with missed partner
1236                       isPartnerLost=kTRUE;
1237                     }
1238                   }
1239                 }
1240                 if(!isPartnerLost &&
1241                    partner->E()<0.3){ //energy is not enough to be registered by PHOS
1242                   FillPIDHistograms("hMCDecWMisPartnEmin",p) ;  //Spectrum of tagged with missed partner
1243                   isPartnerLost=kTRUE;
1244                 }
1245                 if(!isPartnerLost){ //Reason not found!!!!!
1246                   FillPIDHistograms("hMCDecWMisPartnOther",p);
1247                   Int_t multClust = event->GetNumberOfCaloClusters();
1248                   for (Int_t iclu=0; (iclu<multClust) && (!isPartnerLost); iclu++) {
1249                     AliVCluster * clu = event->GetCaloCluster(iclu);
1250                     if(!clu->IsPHOS())
1251                       continue ;
1252                     if(clu->E()==p->E()) //same cluster as current
1253                       continue ;
1254                     Int_t nCluPrimaries = clu->GetNLabels() ;
1255                     for(Int_t iAODLabel=0; (iAODLabel<nCluPrimaries) && (!isPartnerLost); iAODLabel++){
1256                       Int_t ipartnPrim = clu->GetLabelAt(iAODLabel) ;
1257                       while(ipartnPrim>-1){
1258                         if(ipartnPrim==ipartner)
1259                           break ;
1260                         ipartnPrim = ((AliAODMCParticle *)fStack->At(ipartnPrim))->GetMother();
1261                       }
1262                       if(ipartnPrim==ipartner){
1263                         isPartnerLost=kTRUE;
1264                         break ;
1265                       }
1266                     }
1267                   }
1268                   if(isPartnerLost){//Did not pass default cuts
1269                     FillPIDHistograms("hMCDecWMisPartnDefCuts",p);
1270                   }
1271                 }
1272                 else{//Sum of all missed partners
1273                   FillPIDHistograms("hMCDecWMisPartnAll",p);
1274                 }
1275               }//Partner - photon
1276               else{//partner not photon
1277                 FillPIDHistograms("hMCDecWMisPartnNPhot",p);
1278               }
1279
1280             }//Partner not reconstructed
1281           }//Partner in stack
1282         }//photon from pi0 decay
1283       }//photon
1284     } //PHOS clusters
1285
1286 }
1287
1288 //________________________________________________
1289 void AliAnalysisTaskTaggedPhotons::FillTaggingHistos(){
1290   //Fill all necessary histograms
1291
1292
1293  //First fill invariant mass distrubutions and mark tagged photons
1294   //Invariant Mass analysis
1295   const Int_t n=fPHOSEvent->GetEntriesFast() ;
1296   for(Int_t i=0;i<n-1;i++){
1297     AliCaloPhoton *p1 = static_cast<AliCaloPhoton*>(fPHOSEvent->At(i));
1298     for(Int_t j = i+1 ; j < n ; j++) {
1299       AliCaloPhoton * p2 = static_cast<AliCaloPhoton*>(fPHOSEvent->At(j));
1300
1301       //At least one photon should be trigger in PHOS triggered events
1302       if(!fIsMB && !p1->IsTrig() &&  !p2->IsTrig() )
1303         continue ;
1304
1305       Double_t invMass = (*p1 + *p2).M();
1306
1307       Double_t nsigma1 = InPi0Band(invMass,p1->Pt()) ; //in band with n sigmas
1308       Double_t nsigma2 = InPi0Band(invMass,p2->Pt()) ; //in band with n sigmas
1309
1310       if((p1->E()>0.1) && (p2->E()>0.1)){
1311         FillPIDHistograms("hInvM_Re_Emin1",p1,p2,invMass,kTRUE) ;
1312         if((p1->E()>0.2) && (p2->E()>0.2)){
1313           FillPIDHistograms("hInvM_Re_Emin2",p1,p2,invMass,kTRUE) ;
1314           if((p1->E()>0.3) && (p2->E()>0.3)){
1315             FillPIDHistograms("hInvM_Re_Emin3",p1,p2,invMass,kTRUE) ;
1316
1317             //Fill izolated pi0s
1318             if(nsigma1<2 || nsigma2<2){ //2 sigma band
1319               TLorentzVector pi0=*p1+*p2 ;
1320               Int_t isolation=EvalIsolation(&pi0,0) ;
1321               for(Int_t kind=1; kind<33; kind*=2){
1322                  if((isolation&kind)){
1323                    FillHistogram(Form("hPi_Isolation%d_cent%d",kind,fCentBin),pi0.Pt()) ;
1324                  }
1325               }
1326             }
1327
1328           }
1329         }
1330       }
1331
1332
1333       if(p2->E()>0.1){
1334         FillPIDHistograms("hSingleInvM_Re_Emin1",p1,invMass) ;
1335         if(p2->E()>0.2){
1336           FillPIDHistograms("hSingleInvM_Re_Emin2",p1,invMass) ;
1337           if(p2->E()>0.3){
1338             FillPIDHistograms("hSingleInvM_Re_Emin3",p1,invMass) ;
1339           }
1340         }
1341       }
1342
1343       if(p1->E()>0.1){
1344         FillPIDHistograms("hSingleInvM_Re_Emin1",p2,invMass) ;
1345         if(p1->E()>0.2){
1346           FillPIDHistograms("hSingleInvM_Re_Emin2",p2,invMass) ;
1347           if(p1->E()>0.3){
1348             FillPIDHistograms("hSingleInvM_Re_Emin3",p2,invMass) ;
1349           }
1350         }
1351       }
1352       if(TestPID(3, p2)){
1353         if(p2->E()>0.1){
1354           FillPIDHistograms("hSingleInvM_Re_Emin1_Strict",p1,invMass) ;
1355           if(p2->E()>0.2){
1356             FillPIDHistograms("hSingleInvM_Re_Emin2_Strict",p1,invMass) ;
1357             if(p2->E()>0.3){
1358               FillPIDHistograms("hSingleInvM_Re_Emin3_Strict",p1,invMass) ;
1359             }
1360           }
1361         }
1362       }
1363       if(TestPID(3, p1)){
1364         if(p1->E()>0.1){
1365           FillPIDHistograms("hSingleInvM_Re_Emin1_Strict",p2,invMass) ;
1366           if(p1->E()>0.1){
1367             FillPIDHistograms("hSingleInvM_Re_Emin2_Strict",p2,invMass) ;
1368             if(p1->E()>0.3){
1369               FillPIDHistograms("hSingleInvM_Re_Emin3_Strict",p2,invMass) ;
1370             }
1371           }
1372         }
1373       }
1374       if(IsSameParent(p1,p2)==111){
1375         FillPIDHistograms("hMC_InvM_Re",p1,invMass) ;
1376         FillPIDHistograms("hMC_InvM_Re",p2,invMass) ;
1377         if(TestPID(3, p2)){
1378           FillPIDHistograms("hMC_InvM_Re_Strict",p1,invMass) ;
1379         }
1380         if(TestPID(3, p1)){
1381           FillPIDHistograms("hMC_InvM_Re_Strict",p2,invMass) ;
1382         }
1383       }
1384
1385       //Tagging: 1,2,3 sigma
1386       //Emin=100,200,300 Mev
1387       //InPi0Band(mass, Ptphoton, type Emin cut
1388       //Set Bits to 1 if tagged
1389       Int_t tag1=0 ;
1390       for(Int_t eminType=0; eminType<3; eminType++){
1391         if(p2->E()>0.1*(eminType+1)){
1392           //Set corresponding bit
1393           for(Int_t isigma=0; isigma<3; isigma++){
1394             if(nsigma1<1+isigma){
1395               tag1|= (1 << (3*eminType+isigma)) ;
1396               if(p2->IsPIDOK(3))
1397                 tag1|= (1 << (3*eminType+isigma+9)) ;
1398             }
1399           }
1400         }
1401       }
1402       Int_t oldTag1=p1->GetTagInfo() ;
1403       for(Int_t ibit=0; ibit<18; ibit++){
1404         if(((oldTag1 & (1<<ibit))!=0) && //Already tagged
1405            ((tag1 & (1<<ibit))!=0)){//Multiple tagging
1406          Int_t iFidArea = p1->GetFiducialArea();
1407          if(iFidArea>0){
1408            FillPIDHistograms(Form("hPhot_TaggedMult%d_Area1",ibit),p1) ;
1409            if(iFidArea>1){
1410              FillPIDHistograms(Form("hPhot_TaggedMult%d_Area2",ibit),p1) ;
1411              if(iFidArea>2){
1412                FillPIDHistograms(Form("hPhot_TaggedMult%d_Area3",ibit),p1) ;
1413              }
1414            }
1415          }
1416         }
1417       }
1418       tag1=tag1|oldTag1 ;
1419       p1->SetTagInfo(tag1) ;
1420       Int_t tag2=0 ;
1421       for(Int_t eminType=0; eminType<3; eminType++){
1422         if(p1->E()>0.1*(eminType+1)){
1423           //Set corresponding bit
1424           for(Int_t isigma=0; isigma<3; isigma++){
1425             if(nsigma2<1+isigma){
1426               tag2|= (1 << (3*eminType+isigma)) ;
1427               if(p1->IsPIDOK(3))
1428                 tag2|= (1 << (3*eminType+isigma+9)) ;
1429             }
1430           }
1431         }
1432       }
1433       Int_t oldTag2=p2->GetTagInfo() ;
1434       for(Int_t ibit=0; ibit<18; ibit++){
1435         if(((oldTag2 & (1<<ibit))!=0) && //Already tagged
1436            ((tag2 & (1<<ibit))!=0)){//Multiple tagging
1437          Int_t iFidArea = p2->GetFiducialArea();
1438          if(iFidArea>0){
1439            FillPIDHistograms(Form("hPhot_TaggedMult%d_Area1",ibit),p2) ;
1440            if(iFidArea>1){
1441              FillPIDHistograms(Form("hPhot_TaggedMult%d_Area2",ibit),p2) ;
1442              if(iFidArea>2){
1443                FillPIDHistograms(Form("hPhot_TaggedMult%d_Area3",ibit),p2) ;
1444              }
1445            }
1446          }
1447         }
1448       }
1449       tag2=tag2|oldTag2 ;
1450       p2->SetTagInfo(tag2) ;
1451
1452       Bool_t trueTag=(IsSameParent(p1,p2)==111); //true tag
1453       p1->SetUnfolded(p1->IsntUnfolded()|trueTag) ;
1454       p2->SetUnfolded(p2->IsntUnfolded()|trueTag) ;
1455
1456     }
1457   }
1458
1459
1460
1461   //Single particle histgams
1462   for(Int_t i=0;i<n;i++){
1463     AliCaloPhoton *p = static_cast<AliCaloPhoton*>(fPHOSEvent->At(i));
1464
1465    //photon should be trigger in PHOS triggered events
1466    if(!fIsMB && !p->IsTrig() )
1467      continue ;
1468
1469     Int_t isolation = p->GetIsolationTag();
1470
1471     //Inclusive spectra
1472     FillPIDHistograms("hPhot",p) ;
1473
1474     if(p->DistToBad()>1){
1475       FillPIDHistograms("hPhot_Dist2",p) ;
1476       if(p->DistToBad()>2){
1477         FillPIDHistograms("hPhot_Dist3",p) ;
1478       }
1479     }
1480
1481
1482     for(Int_t kind=1; kind<33; kind*=2){
1483       if((isolation&kind)){
1484         FillPIDHistograms(Form("hPhot_Isolation%d",kind),p) ;
1485       }
1486     }
1487
1488     Int_t iFidArea = p->GetFiducialArea();
1489     if(iFidArea>0){
1490       FillPIDHistograms("hPhot_Area1",p) ;
1491       for(Int_t kind=1; kind<33; kind*=2){
1492         if((isolation&kind)){
1493           FillPIDHistograms(Form("hPhot_Isolation%d_Area1",kind),p) ;
1494         }
1495       }
1496
1497       //Fill taggings with
1498       //3 Emin cuts
1499       //Default Emin, 1,2,3 sigmas
1500       //strict and loose PID cut on partner
1501       Int_t tag=p->GetTagInfo() ;
1502       for(Int_t ibit=0; ibit<18; ibit++){
1503         if((tag & (1<<ibit))!=0){
1504           FillPIDHistograms(Form("hPhot_Tagged%d_Area1",ibit),p) ;
1505           if(p->IsntUnfolded()) //true tag
1506              FillPIDHistograms(Form("hPhot_TrueTagged%d",ibit),p) ;
1507
1508         }
1509       }
1510
1511       if(iFidArea>1){
1512         FillPIDHistograms("hPhot_Area2",p) ;
1513         for(Int_t ibit=0; ibit<18; ibit++){
1514           if((tag & (1<<ibit))!=0){
1515             FillPIDHistograms(Form("hPhot_Tagged%d_Area2",ibit),p) ;
1516           }
1517         }
1518         if(iFidArea>2){
1519           FillPIDHistograms("hPhot_Area3",p) ;
1520           for(Int_t ibit=0; ibit<18; ibit++){
1521             if((tag & (1<<ibit))!=0){
1522               FillPIDHistograms(Form("hPhot_Tagged%d_Area3",ibit),p) ;
1523             }
1524           }
1525         }
1526       }
1527     }
1528   }
1529
1530    //Fill Mixed InvMass distributions:
1531   TIter nextEv(fCurrentMixedList) ;
1532   for(Int_t i=0;i<n;i++){
1533     AliCaloPhoton *p1 = static_cast<AliCaloPhoton*>(fPHOSEvent->At(i));
1534     while(TClonesArray * event2 = static_cast<TClonesArray*>(nextEv())){
1535       Int_t nPhotons2 = event2->GetEntriesFast() ;
1536       for(Int_t j=0; j < nPhotons2 ; j++){
1537         AliCaloPhoton * p2 = static_cast<AliCaloPhoton*>(event2->At(j)) ;
1538         Double_t invMass = (*p1 + *p2).M();
1539
1540         //At least one photon should be trigger in PHOS triggered events
1541         if(!fIsMB && !p1->IsTrig() &&  !p2->IsTrig() )
1542           continue ;
1543
1544         if((p1->E()>0.1) && (p2->E()>0.1)){
1545           FillPIDHistograms("hInvM_Mi_Emin1",p1,p2,invMass,kFALSE) ;
1546           if((p1->E())>0.2 && (p2->E()>0.2)){
1547             FillPIDHistograms("hInvM_Mi_Emin2",p1,p2,invMass,kFALSE) ;
1548             if((p1->E())>0.3 && (p2->E()>0.3)){
1549               FillPIDHistograms("hInvM_Mi_Emin3",p1,p2,invMass,kFALSE) ;
1550             }
1551           }
1552         }
1553
1554
1555         if(p2->E()>0.1){
1556           FillPIDHistograms("hSingleInvM_Mi_Emin1",p1,invMass) ;
1557           if(p2->E()>0.2){
1558             FillPIDHistograms("hSingleInvM_Mi_Emin2",p1,invMass) ;
1559             if(p2->E()>0.3){
1560               FillPIDHistograms("hSingleInvM_Mi_Emin3",p1,invMass) ;
1561             }
1562           }
1563         }
1564         if(TestPID(3, p2)){
1565           if(p2->E()>0.1){
1566             FillPIDHistograms("hSingleInvM_Mi_Emin1_Strict",p1,invMass) ;
1567             if(p2->E()>0.2){
1568               FillPIDHistograms("hSingleInvM_Mi_Emin2_Strict",p1,invMass) ;
1569               if(p2->E()>0.3){
1570                 FillPIDHistograms("hSingleInvM_Mi_Emin3_Strict",p1,invMass) ;
1571               }
1572             }
1573           }
1574         }
1575
1576         if(p1->E()>0.1){
1577           FillPIDHistograms("hSingleInvM_Mi_Emin1",p2,invMass) ;
1578           if(p1->E()>0.2){
1579             FillPIDHistograms("hSingleInvM_Mi_Emin2",p2,invMass) ;
1580             if(p1->E()>0.3){
1581               FillPIDHistograms("hSingleInvM_Mi_Emin3",p2,invMass) ;
1582             }
1583           }
1584         }
1585         if(TestPID(3, p1)){
1586           if(p1->E()>0.1){
1587             FillPIDHistograms("hSingleInvM_Mi_Emin1_Strict",p2,invMass) ;
1588             if(p1->E()>0.2){
1589               FillPIDHistograms("hSingleInvM_Mi_Emin2_Strict",p2,invMass) ;
1590               if(p1->E()>0.3){
1591                 FillPIDHistograms("hSingleInvM_Mi_Emin3_Strict",p2,invMass) ;
1592               }
1593             }
1594           }
1595         }
1596       }
1597     }
1598   }
1599
1600 }
1601
1602 //______________________________________________________________________________
1603 void AliAnalysisTaskTaggedPhotons::Init()
1604 {
1605   // Intialisation of parameters
1606 }
1607
1608 //______________________________________________________________________________
1609 void AliAnalysisTaskTaggedPhotons::Terminate(Option_t *)
1610 {
1611   // Processing when the event loop is ended
1612   if (fDebug > 1) Printf("Terminate()");
1613 }
1614 //______________________________________________________________________________
1615 Double_t AliAnalysisTaskTaggedPhotons::InPi0Band(Double_t m, Double_t pt)const
1616 {
1617   //Parameterization of the pi0 peak region
1618   //for LHC13bcdef
1619 //  Double_t mpi0mean =  0.13447 - 1.41259e-03 * TMath::Exp(-pt/1.30044) ;
1620   //Parameterization of data 30.08.2014
1621   Double_t mpi0mean =  0.135 ;
1622
1623   //Double_t mpi0sigma=TMath::Sqrt(5.22245e-03*5.22245e-03 +2.86851e-03*2.86851e-03/pt) + 9.09932e-05*pt ;
1624   //Parameterization of data 30.08.2014
1625   Double_t mpi0sigma=TMath::Sqrt(4.67491e-03*4.67491e-03 +3.42884e-03*3.42884e-03/pt) + 1.24324e-04*pt ;
1626
1627   return TMath::Abs(m-mpi0mean)/mpi0sigma ;
1628 }
1629 //______________________________________________________________________________
1630 Int_t AliAnalysisTaskTaggedPhotons::IsSameParent(const AliCaloPhoton *p1, const AliCaloPhoton *p2)const{
1631   //Looks through parents and finds if there was commont pi0 among ancestors
1632
1633   if(!fIsMC)
1634     return 0 ; //can not say anything
1635
1636   Int_t prim1 = p1->GetPrimary();
1637   while(prim1!=-1){
1638     Int_t prim2 = p2->GetPrimary();
1639
1640     while(prim2!=-1){
1641       if(prim1==prim2){
1642         return ((AliAODMCParticle*)fStack->At(prim1))->GetPdgCode() ;
1643       }
1644       prim2=((AliAODMCParticle*)fStack->At(prim2))->GetMother() ;
1645     }
1646     prim1=((AliAODMCParticle*)fStack->At(prim1))->GetMother() ;
1647   }
1648   return 0 ;
1649 }
1650 //______________________________________________________________________________
1651 Int_t AliAnalysisTaskTaggedPhotons::GetFiducialArea(const Float_t * position)const{
1652   //calculates in which kind of fiducial area photon hit
1653
1654   TVector3 global1(position) ;
1655   Int_t relId[4] ;
1656   fPHOSgeom->GlobalPos2RelId(global1,relId) ;
1657 //  Int_t mod  = relId[0] ;
1658   Int_t cellX = relId[2];
1659   Int_t cellZ = relId[3] ;
1660
1661   //New we are in good channel,
1662   //calculate distance to the closest group of bad channels
1663   const Int_t cut1=10;
1664   const Int_t cut2=15;
1665 //For 3-mod configuration
1666 //  if((mod==3 && cellX<=cut1) || (mod==1 && cellX>=65-cut1) || cellZ<=cut1 || cellZ>=57-cut1)
1667 //For 1+3 configuraion
1668   if( cellX<=cut1 ||  cellX>=65-cut1 || cellZ<=cut1 || cellZ>=57-cut1)
1669     return 1;
1670 //  //and from large dead area
1671 //Full configuration
1672 //    if((mod==3 && cellX<=cut2) || (mod==1 && cellX>=65-cut2) || cellZ<=cut2 || cellZ>=57-cut2)
1673 //1+3 configuration
1674   if( cellX<=cut2 || cellX>=65-cut2 || cellZ<=cut2 || cellZ>=57-cut2)
1675     return 2;
1676   //Very good channel
1677   return 3 ;
1678
1679 }
1680 //______________________________________________________________________________^M
1681 void  AliAnalysisTaskTaggedPhotons::InitGeometry(){
1682   //Rotation matrixes are set with Tender
1683
1684   if(fPHOSgeom) return ;
1685
1686
1687   fPHOSgeom = AliPHOSGeometry::GetInstance() ;
1688
1689   if(!fPHOSgeom){ //Geometry not yet constructed with Tender
1690     fPHOSgeom = AliPHOSGeometry::GetInstance("IHEP","");
1691
1692     AliOADBContainer geomContainer("phosGeo");
1693     geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
1694     TObjArray *matrixes = (TObjArray*)geomContainer.GetObject(170000,"PHOSRotationMatrixes");
1695     for(Int_t mod=0; mod<5; mod++) {
1696       if(!matrixes->At(mod)) continue;
1697       fPHOSgeom->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod) ;
1698     }
1699   }
1700
1701   //Read BadMap for MC simulations
1702   Int_t runNumber=196208 ; //LHC13bcdef
1703   AliOADBContainer badmapContainer(Form("phosBadMap"));
1704   badmapContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSBadMaps.root","phosBadMap");
1705   TObjArray *maps = (TObjArray*)badmapContainer.GetObject(runNumber,"phosBadMap");
1706   if(!maps){
1707       AliError("TaggedPhotons: Can not read Bad map\n") ;
1708   }
1709   else{
1710     AliInfo(Form("TaggedPhotons: Setting PHOS bad map with name %s \n",maps->GetName())) ;
1711     for(Int_t mod=0; mod<5;mod++){
1712       if(fPHOSBadMap[mod])
1713         delete fPHOSBadMap[mod] ;
1714       TH2I * h = (TH2I*)maps->At(mod) ;
1715       if(h)
1716         fPHOSBadMap[mod]=new TH2I(*h) ;
1717     }
1718   }
1719 }
1720 //_____________________________________________________________________________
1721 void AliAnalysisTaskTaggedPhotons::FillHistogram(const char * key,Double_t x)const{
1722   //FillHistogram
1723   TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
1724   if(tmpI){
1725     tmpI->Fill(x) ;
1726     return ;
1727   }
1728   TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
1729   if(tmpF){
1730     tmpF->Fill(x) ;
1731     return ;
1732   }
1733   TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
1734   if(tmpD){
1735     tmpD->Fill(x) ;
1736     return ;
1737   }
1738   AliInfo(Form("can not find histogram <%s> ",key)) ;
1739 }
1740 //_____________________________________________________________________________
1741 void AliAnalysisTaskTaggedPhotons::FillHistogram(const char * key,Double_t x,Double_t y)const{
1742   //FillHistogram
1743   TObject * tmp = fOutputContainer->FindObject(key) ;
1744   if(!tmp){
1745     AliInfo(Form("can not find histogram <%s> ",key)) ;
1746     return ;
1747   }
1748   if(tmp->IsA() == TClass::GetClass("TH1F")){
1749     ((TH1F*)tmp)->Fill(x,y) ;
1750     return ;
1751   }
1752   if(tmp->IsA() == TClass::GetClass("TH2F")){
1753     ((TH2F*)tmp)->Fill(x,y) ;
1754     return ;
1755   }
1756   AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
1757 }
1758
1759 //_____________________________________________________________________________
1760 void AliAnalysisTaskTaggedPhotons::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
1761   //Fills 1D histograms with key
1762   TObject * tmp = fOutputContainer->FindObject(key) ;
1763   if(!tmp){
1764     AliInfo(Form("can not find histogram <%s> ",key)) ;
1765     return ;
1766   }
1767   if(tmp->IsA() == TClass::GetClass("TH2F")){
1768     ((TH2F*)tmp)->Fill(x,y,z) ;
1769     return ;
1770   }
1771   if(tmp->IsA() == TClass::GetClass("TH3F")){
1772     ((TH3F*)tmp)->Fill(x,y,z) ;
1773     return ;
1774   }
1775 }
1776 //_____________________________________________________________________________
1777 void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliCaloPhoton * p) const{
1778
1779   FillHistogram(Form("%s_All_cent%d",name,fCentBin),p->Pt(),p->GetWeight()) ;
1780   if(p->IsDispOK())
1781     FillHistogram(Form("%s_Disp_cent%d",name,fCentBin),p->Pt(),p->GetWeight()) ;
1782   if(p->IsCPVOK())
1783     FillHistogram(Form("%s_CPV_cent%d",name,fCentBin),p->Pt(),p->GetWeight()) ;
1784   if(p->IsDispOK() && p->IsCPVOK())
1785     FillHistogram(Form("%s_Both_cent%d",name,fCentBin),p->Pt(),p->GetWeight()) ;
1786
1787 }
1788 //_____________________________________________________________________________
1789 void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliCaloPhoton * p,Double_t x) const{
1790
1791   FillHistogram(Form("%s_All_cent%d",name,fCentBin),x,p->Pt(),p->GetWeight()) ;
1792   if(p->IsDispOK())
1793     FillHistogram(Form("%s_Disp_cent%d",name,fCentBin),x,p->Pt(),p->GetWeight()) ;
1794   if(p->IsCPVOK())
1795     FillHistogram(Form("%s_CPV_cent%d",name,fCentBin),x,p->Pt(),p->GetWeight()) ;
1796   if(p->IsDispOK() && p->IsCPVOK())
1797     FillHistogram(Form("%s_Both_cent%d",name,fCentBin),x,p->Pt(),p->GetWeight()) ;
1798
1799 }
1800 //_____________________________________________________________________________
1801 void AliAnalysisTaskTaggedPhotons::FillPIDHistograms(const char * name, const AliCaloPhoton * p1,const AliCaloPhoton * p2,Double_t x, Bool_t /*isRe*/) const{
1802
1803   Double_t ptPi = (*p1 + *p2).Pt() ;
1804   Double_t w=TMath::Sqrt(p1->GetWeight()*p2->GetWeight()) ;
1805 //  if(isRe){
1806 //    Int_t pdg=IsSameParent(p1,p2) ;
1807 //    if(pdg!=0)
1808 //      w=p1->GetWeight() ;
1809 //  }
1810   FillHistogram(Form("%s_All_cent%d",name,fCentBin),x,ptPi,w) ;
1811   if(p1->IsDispOK() && p2->IsDispOK())
1812     FillHistogram(Form("%s_Disp_cent%d",name,fCentBin),x,ptPi,w) ;
1813   if(p1->IsCPVOK() && p2->IsCPVOK())
1814     FillHistogram(Form("%s_CPV_cent%d",name,fCentBin),x,ptPi,w) ;
1815   if(p1->IsDispOK() && p1->IsCPVOK() && p2->IsDispOK() && p2->IsCPVOK())
1816     FillHistogram(Form("%s_Both_cent%d",name,fCentBin),x,ptPi,w) ;
1817
1818 }
1819 //_____________________________________________________________________________
1820 Bool_t AliAnalysisTaskTaggedPhotons::TestLambda(Double_t pt,Double_t l1,Double_t l2){
1821
1822   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1823   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1824   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1825   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1826   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1827   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1828               0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1829               0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1830   return (R2<2.5*2.5) ;
1831
1832 }
1833 //_________________________________________________________________________________
1834 Int_t AliAnalysisTaskTaggedPhotons::EvalIsolation(TLorentzVector * ph, Bool_t isPhoton){
1835
1836    // Check if this particle is isolated.
1837    //We use several cone radii and epsilons.
1838    //As well we look at charged particles and EMCAL ones
1839
1840    const Double_t coneR1=0.3 ;
1841    const Double_t coneR2=0.4 ;
1842    const Double_t coneR3=0.5 ;
1843
1844    const Double_t cutEcone1=3. ;
1845    const Double_t cutEcone2=4.5 ;
1846    const Double_t cutEcone3=6.5 ;
1847
1848    if(!ph) return 0 ;
1849
1850    //Sum of energies in cones, tracks and clusters in EMCAL
1851    Double_t eCone1 = 0;
1852    Double_t eCone2 = 0;
1853    Double_t eCone3 = 0;
1854    //Cross-check, energy in cone perpr to photon
1855    Double_t eP1Cone1 = 0;
1856    Double_t eP1Cone2 = 0;
1857    Double_t eP1Cone3 = 0;
1858    Double_t eP2Cone1 = 0;
1859    Double_t eP2Cone2 = 0;
1860    Double_t eP2Cone3 = 0;
1861
1862
1863    Double_t  phiTrig = ph->Phi();
1864    Double_t  etaTrig = ph->Eta();
1865
1866    Int_t n=fTrackEvent->GetEntriesFast() ;
1867    for(Int_t itr=0; itr<n; itr++){
1868      AliCaloPhoton * track = (AliCaloPhoton*)fTrackEvent->At(itr) ;
1869
1870      Double_t deleta = etaTrig - track->Eta() ;
1871      Double_t delphi = phiTrig - track->Phi() ;
1872      while(delphi<-TMath::Pi()) delphi+=TMath::TwoPi() ;
1873      while(delphi>TMath::Pi()) delphi-=TMath::TwoPi() ;
1874
1875      //Perp cones
1876      Double_t delphiP1 = phiTrig +TMath::PiOver2()- track->Phi() ;
1877      while(delphiP1<-TMath::Pi()) delphiP1+=TMath::TwoPi() ;
1878      while(delphiP1>TMath::Pi()) delphiP1-=TMath::TwoPi() ;
1879      Double_t delphiP2 = phiTrig -TMath::PiOver2()- track->Phi() ;
1880      while(delphiP2<-TMath::Pi()) delphiP2+=TMath::TwoPi() ;
1881      while(delphiP2>TMath::Pi()) delphiP2-=TMath::TwoPi() ;
1882
1883
1884      Double_t dr    = TMath::Sqrt(deleta * deleta + delphi * delphi);
1885      Double_t drP1    = TMath::Sqrt(deleta * deleta + delphiP1 * delphiP1);
1886      Double_t drP2    = TMath::Sqrt(deleta * deleta + delphiP2 * delphiP2);
1887
1888      if(dr<coneR3){
1889        eCone3+=track->Pt() ;
1890        if(dr<coneR2){
1891          eCone2+=track->Pt() ;
1892          if(dr<coneR1){
1893            eCone1+=track->Pt() ;
1894          }
1895         }
1896       }
1897
1898     if(drP1<coneR3){
1899        eP1Cone3+=track->Pt() ;
1900        if(drP1<coneR2){
1901          eP1Cone2+=track->Pt() ;
1902          if(drP1<coneR1){
1903            eP1Cone1+=track->Pt() ;
1904          }
1905        }
1906     }
1907
1908     if(drP2<coneR3){
1909        eP2Cone3+=track->Pt() ;
1910        if(drP2<coneR2){
1911          eP2Cone2+=track->Pt() ;
1912          if(drP2<coneR1){
1913            eP2Cone1+=track->Pt() ;
1914          }
1915         }
1916       }
1917     }
1918
1919     //Fill QA histgams
1920     Double_t ptTrig=ph->Pt() ;
1921     if(isPhoton){
1922     FillHistogram(Form("QA_Cone1_Tracks_cent%d",fCentBin),ptTrig,eCone1) ;
1923     FillHistogram(Form("QA_Cone2_Tracks_cent%d",fCentBin),ptTrig,eCone2) ;
1924     FillHistogram(Form("QA_Cone3_Tracks_cent%d",fCentBin),ptTrig,eCone3) ;
1925     FillHistogram(Form("QA_PCone1_Tracks_cent%d",fCentBin),ptTrig,eP1Cone1) ;
1926     FillHistogram(Form("QA_PCone2_Tracks_cent%d",fCentBin),ptTrig,eP1Cone2) ;
1927     FillHistogram(Form("QA_PCone3_Tracks_cent%d",fCentBin),ptTrig,eP1Cone3) ;
1928     FillHistogram(Form("QA_PCone1_Tracks_cent%d",fCentBin),ptTrig,eP2Cone1) ;
1929     FillHistogram(Form("QA_PCone2_Tracks_cent%d",fCentBin),ptTrig,eP2Cone2) ;
1930     FillHistogram(Form("QA_PCone3_Tracks_cent%d",fCentBin),ptTrig,eP2Cone3) ;
1931     }
1932     else{
1933     FillHistogram(Form("QA_Pi0Cone1_Tracks_cent%d",fCentBin),ptTrig,eCone1) ;
1934     FillHistogram(Form("QA_Pi0Cone2_Tracks_cent%d",fCentBin),ptTrig,eCone2) ;
1935     FillHistogram(Form("QA_Pi0Cone3_Tracks_cent%d",fCentBin),ptTrig,eCone3) ;
1936     FillHistogram(Form("QA_Pi0PCone1_Tracks_cent%d",fCentBin),ptTrig,eP1Cone1) ;
1937     FillHistogram(Form("QA_Pi0PCone2_Tracks_cent%d",fCentBin),ptTrig,eP1Cone2) ;
1938     FillHistogram(Form("QA_Pi0PCone3_Tracks_cent%d",fCentBin),ptTrig,eP1Cone3) ;
1939     FillHistogram(Form("QA_Pi0PCone1_Tracks_cent%d",fCentBin),ptTrig,eP2Cone1) ;
1940     FillHistogram(Form("QA_Pi0PCone2_Tracks_cent%d",fCentBin),ptTrig,eP2Cone2) ;
1941     FillHistogram(Form("QA_Pi0PCone3_Tracks_cent%d",fCentBin),ptTrig,eP2Cone3) ;
1942
1943     }
1944
1945     //Fill Bits
1946     Int_t iCone1E1 = (cutEcone1 > eCone1) ;
1947     Int_t iCone2E1 = (cutEcone2 > eCone2) ;
1948     Int_t iCone3E1 = (cutEcone3 > eCone3) ;
1949
1950     Int_t iCone1E2 = (1.5*cutEcone1 > eCone1) ;
1951     Int_t iCone2E2 = (1.5*cutEcone2 > eCone2) ;
1952     Int_t iCone3E2 = (1.5*cutEcone3 > eCone3) ;
1953
1954
1955     Int_t isolation=   iCone1E1+  2*iCone2E1   +4*iCone3E1+
1956                     8*iCone1E2 +16*iCone2E2 +32*iCone3E2 ;
1957     return isolation ;
1958 }
1959 //_________________________________________________________________________________
1960 Bool_t AliAnalysisTaskTaggedPhotons::TestPID(Int_t iPID, AliCaloPhoton* part){
1961 //   //Checks PID of particle
1962
1963   if(!part) return kFALSE ;
1964   if(iPID==0) return kTRUE ;
1965   if(iPID==1) return part->IsDispOK() ;
1966   if(iPID==2) return part->IsCPVOK() ;
1967   if(iPID==3) return part->IsDispOK() && part->IsCPVOK() ;
1968   return kFALSE ;
1969
1970 }
1971 //_________________________________________________________________________________
1972 Double_t AliAnalysisTaskTaggedPhotons::PrimaryParticleWeight(AliAODMCParticle * particle){
1973   if(!fIsMC)  //This is real data
1974      return 1;
1975   //Classify parent at vertex
1976   //Introduce for eta and pi0 weights
1977
1978   Double_t r2=particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv() ;
1979   Int_t mother = particle->GetMother() ;
1980   while(mother>-1){
1981     if(r2<1.)
1982       break ;
1983     particle = (AliAODMCParticle*) fStack->At(mother);
1984     mother = particle->GetMother() ;
1985     r2=particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv() ;
1986   }
1987
1988   //Particle within 1 cm from the virtex
1989   Int_t pdg = particle->GetPdgCode() ;
1990   Double_t x = particle->Pt() ;
1991   mother = particle->GetMother() ;
1992   while(TMath::Abs(pdg)<100){//gamma, electrons, muons
1993     if(mother>-1){
1994       AliAODMCParticle * tmpP=(AliAODMCParticle*)fStack->At(mother) ;
1995       pdg=tmpP->GetPdgCode() ;
1996       x = tmpP->Pt() ;
1997       mother = tmpP->GetMother() ;
1998     }
1999     else{ //direct photon/electron....
2000       return 1.;
2001     }
2002   }
2003   if(pdg == 111){
2004   //Pi0
2005      if(x<1) return 1. ;
2006      else return fWeightParamPi0[0]*(TMath::Power(x,fWeightParamPi0[1])*
2007        (1.+fWeightParamPi0[2]*x+fWeightParamPi0[3]*x*x)/
2008        (1.+fWeightParamPi0[4]*x+fWeightParamPi0[5]*x*x) +fWeightParamPi0[6])  ;
2009   }
2010   if(pdg == 221){
2011   //Eta - same same shape, but yield 0.48 and Br(eta->2gamma)
2012      Double_t norm=0.48 * 0.3943;
2013      if(x<1) return norm ;
2014      else return norm*fWeightParamPi0[0]*(TMath::Power(x,fWeightParamPi0[1])*
2015        (1.+fWeightParamPi0[2]*x+fWeightParamPi0[3]*x*x)/
2016        (1.+fWeightParamPi0[4]*x+fWeightParamPi0[5]*x*x) +fWeightParamPi0[6]) ;
2017   }
2018   return 1. ;
2019 }
2020 //_________________________________________________________________________________
2021 void AliAnalysisTaskTaggedPhotons::SetPi0WeightParameters(TArrayD * ar){
2022
2023   for(Int_t i=0; i<7; i++){ //Array range
2024     if(ar->GetSize()>i) fWeightParamPi0[i]=ar->At(i) ;
2025     else fWeightParamPi0[i]=0.;
2026   }
2027   //normalize to obtain smooth transition at 1 GeV
2028   Double_t x=1. ;
2029   fWeightParamPi0[0]=1./(TMath::Power(x,fWeightParamPi0[1])*
2030        (1.+fWeightParamPi0[2]*x+fWeightParamPi0[3]*x*x)/
2031        (1.+fWeightParamPi0[4]*x+fWeightParamPi0[5]*x*x) +fWeightParamPi0[6]) ;
2032
2033
2034 }
2035 //___________________________________________________________________________
2036 Int_t AliAnalysisTaskTaggedPhotons::FindPrimary(AliVCluster*clu,  Bool_t&sure){
2037   //Finds primary and estimates if it unique one?
2038   //First check can it be photon/electron
2039   const Double_t emFraction=0.9; //part of energy of cluster to be assigned to EM particle
2040   Int_t n=clu->GetNLabels() ;
2041   FillHistogram(Form("LabelsNPrim_cent%d",fCentBin),clu->E(),n) ;
2042   if(n==1){
2043     sure=kTRUE;
2044     return clu->GetLabelAt(0) ;
2045   }
2046   //Number of clusters with one or more photons
2047   Bool_t hasGamma=kFALSE ;
2048   Double_t eMax=0. ;
2049   for(Int_t i=0;  i<n;  i++){
2050     AliAODMCParticle*  p=  (AliAODMCParticle*)fStack->At(clu->GetLabelAt(i)) ;
2051     Int_t pdg = p->GetPdgCode() ;
2052     if(pdg==22){
2053       hasGamma=kTRUE ;
2054       if(p->E()>eMax){
2055         eMax=p->E();
2056       }
2057     }
2058   }
2059   if(hasGamma){
2060     FillHistogram(Form("LabelsGamma_cent%d",fCentBin),clu->E()) ;
2061     FillHistogram(Form("LabelsGammaE_cent%d",fCentBin),clu->E(),eMax/clu->E()) ;
2062   }
2063
2064   for(Int_t i=0;  i<n;  i++){
2065     AliAODMCParticle*  p= (AliAODMCParticle*) fStack->At(clu->GetLabelAt(i)) ;
2066     Int_t pdg = p->GetPdgCode() ;
2067     if(pdg==22  ||  pdg==11 || pdg == -11){
2068       if(p->E()>emFraction*clu->E()){
2069         sure=kTRUE ;
2070         return clu->GetLabelAt(i);
2071       }
2072     }
2073   }
2074
2075   Double_t*  Ekin=  new  Double_t[n] ;
2076   for(Int_t i=0;  i<n;  i++){
2077     AliAODMCParticle*  p=(AliAODMCParticle*) fStack->At(clu->GetLabelAt(i)) ;
2078     Ekin[i]=p->P() ;  // estimate of kinetic energy
2079     if(p->GetPdgCode()==-2212  ||  p->GetPdgCode()==-2112){
2080       Ekin[i]+=1.8  ;  //due to annihilation
2081     }
2082   }
2083   Int_t iMax=0;
2084   Double_t eSubMax=0. ;
2085   eMax=0.;
2086   for(Int_t i=0;  i<n;  i++){
2087     if(Ekin[i]>eMax){
2088       eSubMax=eMax;
2089       eMax=Ekin[i];
2090       iMax=i;
2091     }
2092   }
2093   if(eSubMax>0.8*eMax)//not obvious primary
2094     sure=kFALSE;
2095   else
2096     sure=kTRUE;
2097   delete[]  Ekin;
2098   return  clu->GetLabelAt(iMax) ;
2099 }
2100 //___________________________________________________________________________
2101 Bool_t AliAnalysisTaskTaggedPhotons::IsGoodChannel(Int_t mod, Int_t ix, Int_t iz)
2102 {
2103   //Check if this channel belogs to the good ones
2104   //Used only for estimate of photon impact in MC simulations
2105   //For true clusters a bad map in Tender should be used
2106
2107   if(mod>4 || mod<1){
2108     return kFALSE ;
2109   }
2110   if(!fPHOSBadMap[mod]){
2111      return kTRUE ;
2112   }
2113   if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
2114     return kFALSE ;
2115   else
2116     return kTRUE ;
2117 }