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