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