PWGGA/PHOSTasks/PHOS_embedding/AliAnalysisTaskPi0DiffEfficiency.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TObjArray.h"
   4 #include "TF1.h"
   5 #include "TFile.h"
   6 #include "TH1F.h"
   7 #include "TH2F.h"
   8 #include "TH2I.h"
   9 #include "TH3F.h"
  10 #include "TParticle.h"
  11 #include "TCanvas.h"
  12 #include "TStyle.h"
  13 #include "TRandom.h"
  14 #include "THashList.h"
  15
  16 #include "AliAODMCParticle.h"
  17 #include "AliAnalysisManager.h"
  18 #include "AliMCEventHandler.h"
  19 #include "AliMCEvent.h"
  20 #include "AliStack.h"
  21 #include "AliAnalysisTaskSE.h"
  22 #include "AliAnalysisTaskPi0DiffEfficiency.h"
  23 #include "AliCaloPhoton.h"
  24 #include "AliPHOSGeometry.h"
  25 #include "AliPHOSAodCluster.h"
  26 #include "AliPHOSCalibData.h"
  27 #include "AliAODEvent.h"
  28 #include "AliAODCaloCluster.h"
  29 #include "AliAODVertex.h"
  30 #include "AliESDtrackCuts.h"
  31 #include "AliLog.h"
  32 #include "AliPID.h"
  33 #include "AliCDBManager.h"
  34 #include "AliCentrality.h"
  35
  36 // Analysis task to calculate pi0 registration efficiency
  37 // as a difference between spectra before and after embedding
  38 // Authors: Dmitri Peressounko
  39 // Date   : Aug.2011
  40
  41 Double_t Scale(Double_t x){
  42
  43 //  return 1./1.008/1.015*(1.+0.017/(1.+x*x/2./2.)+0.03/(1.+x*x/0.6/0.6)) ;
  44   return 1./1.015/1.015*(1.+0.017/(1.+x*x/2./2.)+0.03/(1.+x*x/0.6/0.6)) ;
  45
  46
  47 }
  48
  49 ClassImp(AliAnalysisTaskPi0DiffEfficiency)
  50
  51 //________________________________________________________________________
  52 AliAnalysisTaskPi0DiffEfficiency::AliAnalysisTaskPi0DiffEfficiency(const char *name)
  53 : AliAnalysisTaskSE(name),
  54   fStack(0),
  55   fOutputContainer(0),
  56   fPHOSEvent(0),
  57   fPHOSEvent1(0),
  58   fPHOSEvent2(0),
  59   fPHOSCalibData(0),
  60   fNonLinCorr(0),
  61   fRPfull(0),
  62   fRPA(0),
  63   fRPC(0),
  64   fRPFar(0),
  65   fRPAFar(0),
  66   fRPCFar(0),
  67   fCentrality(0),
  68   fCenBin(0),
  69   fPHOSGeo(0),
  70   fEventCounter(0)
  71 {
  72   // Constructor
  73   for(Int_t i=0;i<1;i++){
  74     for(Int_t j=0;j<10;j++)
  75       for(Int_t k=0;k<11;k++)
  76         fPHOSEvents[i][j][k]=0 ;
  77   }
  78
  79   // Output slots #0 write into a TH1 container
  80   DefineOutput(1,TList::Class());
  81
  82   // Set bad channel map
  83   char key[55] ;
  84   for(Int_t i=0; i<6; i++){
  85     snprintf(key,55,"PHOS_BadMap_mod%d",i) ;
  86     fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
  87   }
  88   // Initialize the PHOS geometry
  89   fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
  90
  91   fPHOSCalibData = new AliPHOSCalibData();
  92   for(Int_t module=1; module<=5; module++) {
  93     for(Int_t column=1; column<=56; column++) {
  94       for(Int_t row=1; row<=64; row++) {
  95         fPHOSCalibData->SetADCchannelEmc(module,column,row,1.);
  96       }
  97     }
  98   }
  99
 100
 101
 102 }
 103
 104 //________________________________________________________________________
 105 void AliAnalysisTaskPi0DiffEfficiency::UserCreateOutputObjects()
 106 {
 107   // Create histograms
 108   // Called once
 109
 110   // ESD histograms
 111   if(fOutputContainer != NULL){
 112     delete fOutputContainer;
 113   }
 114   fOutputContainer = new THashList();
 115   fOutputContainer->SetOwner(kTRUE);
 116
 117   //Event selection
 118   fOutputContainer->Add(new TH1F("hSelEvents","Event celection", 10,0.,10.)) ;
 119
 120   //vertex distribution
 121   fOutputContainer->Add(new TH1F("hZvertex","Z vertex position", 50,-25.,25.)) ;
 122
 123   //Centrality
 124   fOutputContainer->Add(new TH1F("hCentrality","Event centrality", 100,0.,100.)) ;
 125
 126   //QA histograms
 127   fOutputContainer->Add(new TH2F("hCluM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 128   fOutputContainer->Add(new TH2F("hCluM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 129   fOutputContainer->Add(new TH2F("hCluM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 130
 131   Int_t nM       = 500;
 132   Double_t mMin  = 0.0;
 133   Double_t mMax  = 1.0;
 134   Int_t nPt      = 250;
 135   Double_t ptMin = 0;
 136   Double_t ptMax = 25;
 137
 138   Int_t nPtf      = 23;
 139   Double_t xPt[24]={0.6,1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.,5.5,6.,7.,8.,9.,10.,12.,14.,16.,18.,20.,22.,24.,25.} ;
 140   Int_t nPhi=10 ;
 141   Double_t xPhi[11] ;
 142   for(Int_t i=0;i<=10;i++)xPhi[i]=i*0.1*TMath::Pi() ;
 143   Int_t nMm=150 ;
 144   Double_t xM[201] ;
 145   for(Int_t i=0;i<=200;i++){xM[i]=0.0025*i;}
 146
 147
 148   char key[55] ;
 149   for(Int_t cent=0; cent<6; cent++){
 150     //Single photon
 151
 152     fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad2_cen%d",cent),"All clusters",nPt,ptMin,ptMax));
 153     fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad4_cen%d",cent),"All clusters",nPt,ptMin,ptMax));
 154     fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad6_cen%d",cent),"All clusters",nPt,ptMin,ptMax));
 155     snprintf(key,55,"hPhotAll_cen%d",cent) ;
 156     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 157     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 158     snprintf(key,55,"hPhotAllcore_cen%d",cent) ;
 159     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 160     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 161     snprintf(key,55,"hPhotAllwou_cen%d",cent) ;
 162     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 163     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 164     snprintf(key,55,"hPhotCPV_cen%d",cent) ;
 165     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 166     snprintf(key,55,"hPhotCPVcore_cen%d",cent) ;
 167     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 168     snprintf(key,55,"hPhotCPV2_cen%d",cent) ;
 169     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 170     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 171     snprintf(key,55,"hPhotCPV2core_cen%d",cent) ;
 172     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 173     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 174     snprintf(key,55,"hPhotDisp_cen%d",cent) ;
 175     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 176     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 177     snprintf(key,55,"hPhotDispwou_cen%d",cent) ;
 178     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 179     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 180     snprintf(key,55,"hPhotDisp2_cen%d",cent) ;
 181     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 182     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 183     snprintf(key,55,"hPhotDispcore_cen%d",cent) ;
 184     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 185     snprintf(key,55,"hPhotDisp2core_cen%d",cent) ;
 186     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 187     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 188     snprintf(key,55,"hPhotBoth_cen%d",cent) ;
 189     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 190     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 191     snprintf(key,55,"hPhotBothcore_cen%d",cent) ;
 192     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 193     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 194     snprintf(key,55,"hPhotBoth2_cen%d",cent) ;
 195     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 196     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 197     snprintf(key,55,"hPhotBoth2core_cen%d",cent) ;
 198     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 199     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 200
 201     snprintf(key,55,"hNegPhotAll_cen%d",cent) ;
 202     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 203     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 204     snprintf(key,55,"hNegPhotAllcore_cen%d",cent) ;
 205     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 206     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 207     snprintf(key,55,"hNegPhotCPV_cen%d",cent) ;
 208     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 209     snprintf(key,55,"hNegPhotCPVcore_cen%d",cent) ;
 210     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 211     snprintf(key,55,"hNegPhotCPV2_cen%d",cent) ;
 212     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 213     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 214     snprintf(key,55,"hNegPhotDisp_cen%d",cent) ;
 215     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 216     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 217     snprintf(key,55,"hNegPhotDisp2_cen%d",cent) ;
 218     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 219     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 220     snprintf(key,55,"hNegPhotDispcore_cen%d",cent) ;
 221     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 222     snprintf(key,55,"hNegPhotDisp2core_cen%d",cent) ;
 223     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 224     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 225     snprintf(key,55,"hNegPhotBoth_cen%d",cent) ;
 226     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 227     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 228     snprintf(key,55,"hNegPhotBothcore_cen%d",cent) ;
 229     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 230     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 231     snprintf(key,55,"hNegPhotBoth2_cen%d",cent) ;
 232     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 233     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 234     snprintf(key,55,"hNegPhotBoth2core_cen%d",cent) ;
 235     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 236     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 237
 238     snprintf(key,55,"hOldMassPtAll_cen%d",cent) ;
 239     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 240     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 241     snprintf(key,55,"hOldMassPtAllcore_cen%d",cent) ;
 242     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 243     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 244     snprintf(key,55,"hOldMassPtCPV_cen%d",cent) ;
 245     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 246     snprintf(key,55,"hOldMassPtCPVcore_cen%d",cent) ;
 247     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 248     snprintf(key,55,"hOldMassPtCPV2_cen%d",cent) ;
 249     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 250     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 251     snprintf(key,55,"hOldMassPtDisp_cen%d",cent) ;
 252     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 253     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 254     snprintf(key,55,"hOldMassPtDisp2_cen%d",cent) ;
 255     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 256     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 257     snprintf(key,55,"hOldMassPtDispcore_cen%d",cent) ;
 258     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 259     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 260     snprintf(key,55,"hOldMassPtDisp2core_cen%d",cent) ;
 261     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 262     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 263     snprintf(key,55,"hOldMassPtBoth_cen%d",cent) ;
 264     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 265     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 266     snprintf(key,55,"hOldMassPtBothcore_cen%d",cent) ;
 267     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 268     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 269     snprintf(key,55,"hOldMassPtBoth2_cen%d",cent) ;
 270     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 271     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 272     snprintf(key,55,"hOldMassPtBoth2core_cen%d",cent) ;
 273     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 274     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 275
 276     snprintf(key,55,"hNewMassPtAll_cen%d",cent) ;
 277     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 278     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 279     snprintf(key,55,"hNewMassPtAllcore_cen%d",cent) ;
 280     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 281     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 282     snprintf(key,55,"hNewMassPtCPV_cen%d",cent) ;
 283     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 284     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 285     snprintf(key,55,"hNewMassPtCPVcore_cen%d",cent) ;
 286     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 287     snprintf(key,55,"hNewMassPtCPV2_cen%d",cent) ;
 288     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 289     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 290     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 291     snprintf(key,55,"hNewMassPtDisp_cen%d",cent) ;
 292     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 293     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 294     snprintf(key,55,"hNewMassPtDisp2_cen%d",cent) ;
 295     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 296     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 297     snprintf(key,55,"hNewMassPtDispcore_cen%d",cent) ;
 298     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 299     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 300     snprintf(key,55,"hNewMassPtDisp2core_cen%d",cent) ;
 301     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 302     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 303     snprintf(key,55,"hNewMassPtBoth_cen%d",cent) ;
 304     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 305     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 306     snprintf(key,55,"hNewMassPtBothcore_cen%d",cent) ;
 307     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 308     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 309     snprintf(key,55,"hNewMassPtBoth2_cen%d",cent) ;
 310     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 311     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 312     snprintf(key,55,"hNewMassPtBoth2core_cen%d",cent) ;
 313     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 314     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 315
 316     snprintf(key,55,"hNegMassPtAll_cen%d",cent) ;
 317     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 318     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 319     snprintf(key,55,"hNegMassPtAllcore_cen%d",cent) ;
 320     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 321     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 322     snprintf(key,55,"hNegMassPtCPV_cen%d",cent) ;
 323     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 324     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 325     snprintf(key,55,"hNegMassPtCPVcore_cen%d",cent) ;
 326     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 327     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 328     snprintf(key,55,"hNegMassPtCPV2_cen%d",cent) ;
 329     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 330     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 331     snprintf(key,55,"hNegMassPtDisp_cen%d",cent) ;
 332     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 333     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 334     snprintf(key,55,"hNegMassPtDisp2_cen%d",cent) ;
 335     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 336     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 337     snprintf(key,55,"hNegMassPtDispcore_cen%d",cent) ;
 338     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 339     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 340     snprintf(key,55,"hNegMassPtDisp2core_cen%d",cent) ;
 341     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 342     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 343     snprintf(key,55,"hNegMassPtBoth_cen%d",cent) ;
 344     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 345     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 346     snprintf(key,55,"hNegMassPtBothcore_cen%d",cent) ;
 347     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 348     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 349     snprintf(key,55,"hNegMassPtBoth2_cen%d",cent) ;
 350     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 351     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 352     snprintf(key,55,"hNegMassPtBoth2core_cen%d",cent) ;
 353     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 354     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 355
 356
 357     snprintf(key,55,"hMassPtAll_cen%d",cent) ;
 358     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 359     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 360     snprintf(key,55,"hMassPtAllwou_cen%d",cent) ;
 361     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 362     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 363     snprintf(key,55,"hMassPtAllcore_cen%d",cent) ;
 364     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 365     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 366     snprintf(key,55,"hMassPtCPV_cen%d",cent) ;
 367     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 368     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 369     snprintf(key,55,"hMassPtCPVcore_cen%d",cent) ;
 370     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 371     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 372     snprintf(key,55,"hMassPtCPV2_cen%d",cent) ;
 373     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 374     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 375     snprintf(key,55,"hMassPtDisp_cen%d",cent) ;
 376     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 377     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 378     snprintf(key,55,"hMassPtDispwou_cen%d",cent) ;
 379     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 380     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 381     snprintf(key,55,"hMassPtDisp2_cen%d",cent) ;
 382     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 383     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 384     snprintf(key,55,"hMassPtDispcore_cen%d",cent) ;
 385     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 386     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 387     snprintf(key,55,"hMassPtDisp2core_cen%d",cent) ;
 388     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 389     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 390     snprintf(key,55,"hMassPtBoth_cen%d",cent) ;
 391     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 392     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 393     snprintf(key,55,"hMassPtBothcore_cen%d",cent) ;
 394     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 395     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 396     snprintf(key,55,"hMassPtBoth2_cen%d",cent) ;
 397     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 398     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 399     snprintf(key,55,"hMassPtBoth2core_cen%d",cent) ;
 400     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 401     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 402
 403
 404     //Mixed
 405     snprintf(key,55,"hMiMassPtAll_cen%d",cent) ;
 406     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 407     snprintf(key,55,"hMiMassPtAllwou_cen%d",cent) ;
 408     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 409     snprintf(key,55,"hMiMassPtAllcore_cen%d",cent) ;
 410     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 411     snprintf(key,55,"hMiMassPtCPV_cen%d",cent) ;
 412     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 413     snprintf(key,55,"hMiMassPtCPVcore_cen%d",cent) ;
 414     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 415     snprintf(key,55,"hMiMassPtCPV2_cen%d",cent) ;
 416     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 417     snprintf(key,55,"hMiMassPtDisp_cen%d",cent) ;
 418     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 419     snprintf(key,55,"hMiMassPtDispwou_cen%d",cent) ;
 420     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 421     snprintf(key,55,"hMiMassPtDisp2_cen%d",cent) ;
 422     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 423     snprintf(key,55,"hMiMassPtDispcore_cen%d",cent) ;
 424     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 425     snprintf(key,55,"hMiMassPtDisp2core_cen%d",cent) ;
 426     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 427     snprintf(key,55,"hMiMassPtBoth_cen%d",cent) ;
 428     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 429     snprintf(key,55,"hMiMassPtBothcore_cen%d",cent) ;
 430     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 431     snprintf(key,55,"hMiMassPtBoth2_cen%d",cent) ;
 432     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 433     snprintf(key,55,"hMiMassPtBoth2core_cen%d",cent) ;
 434     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 435
 436     snprintf(key,55,"hMCMassPtAll_cen%d",cent) ;
 437     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 438     snprintf(key,55,"hMCMassPtAllwou_cen%d",cent) ;
 439     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 440     snprintf(key,55,"hMCMassPtAllcore_cen%d",cent) ;
 441     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 442     snprintf(key,55,"hMCMassPtCPV_cen%d",cent) ;
 443     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 444     snprintf(key,55,"hMCMassPtCPVcore_cen%d",cent) ;
 445     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 446     snprintf(key,55,"hMCMassPtCPV2_cen%d",cent) ;
 447     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 448     snprintf(key,55,"hMCMassPtCPV2core_cen%d",cent) ;
 449     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 450     snprintf(key,55,"hMCMassPtDisp_cen%d",cent) ;
 451     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 452     snprintf(key,55,"hMCMassPtDispwou_cen%d",cent) ;
 453     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 454     snprintf(key,55,"hMCMassPtDispcore_cen%d",cent) ;
 455     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 456     snprintf(key,55,"hMCMassPtDisp2_cen%d",cent) ;
 457     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 458     snprintf(key,55,"hMCMassPtDisp2core_cen%d",cent) ;
 459     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 460     snprintf(key,55,"hMCMassPtBoth_cen%d",cent) ;
 461     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 462     snprintf(key,55,"hMCMassPtBoth2_cen%d",cent) ;
 463     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 464     snprintf(key,55,"hMCMassPtBothcore_cen%d",cent) ;
 465     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 466     snprintf(key,55,"hMCMassPtBoth2core_cen%d",cent) ;
 467     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 468
 469     //Single photon
 470     snprintf(key,55,"hMCPhotAll_cen%d",cent) ;
 471     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 472     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 473     snprintf(key,55,"hMCPhotAllwou_cen%d",cent) ;
 474     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 475     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 476     snprintf(key,55,"hMCPhotAllcore_cen%d",cent) ;
 477     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 478     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 479     snprintf(key,55,"hMCPhotCPV_cen%d",cent) ;
 480     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 481     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 482     snprintf(key,55,"hMCPhotCPVcore_cen%d",cent) ;
 483     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 484     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 485     snprintf(key,55,"hMCPhotCPV2_cen%d",cent) ;
 486     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 487     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 488     snprintf(key,55,"hMCPhotDisp_cen%d",cent) ;
 489     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 490     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 491     snprintf(key,55,"hMCPhotDispwou_cen%d",cent) ;
 492     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 493     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 494     snprintf(key,55,"hMCPhotDisp2_cen%d",cent) ;
 495     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 496     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 497     snprintf(key,55,"hMCPhotDispcore_cen%d",cent) ;
 498     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 499     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 500     snprintf(key,55,"hMCPhotDisp2core_cen%d",cent) ;
 501     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 502     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 503     snprintf(key,55,"hMCPhotBoth_cen%d",cent) ;
 504     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 505     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 506     snprintf(key,55,"hMCPhotBothcore_cen%d",cent) ;
 507     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 508     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 509     snprintf(key,55,"hMCPhotBoth2_cen%d",cent) ;
 510     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 511     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 512     snprintf(key,55,"hMCPhotBoth2core_cen%d",cent) ;
 513     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 514     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 515
 516     char phiTitle[15]={"TPC"};
 517     TH2F * tmp = 0;
 518     tmp=new TH2F(Form("hPhotPhi%sAll_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPtf,xPt,nPhi,xPhi) ;
 519     tmp->Sumw2() ; fOutputContainer->Add(tmp) ;
 520     tmp=new TH2F(Form("hPhotPhi%sAllcore_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPtf,xPt,nPhi,xPhi) ;
 521     tmp->Sumw2() ; fOutputContainer->Add(tmp) ;
 522     tmp=new TH2F(Form("hPhotPhi%sDisp_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPtf,xPt,nPhi,xPhi) ;
 523     tmp->Sumw2() ; fOutputContainer->Add(tmp) ;
 524     tmp=new TH2F(Form("hPhotPhi%sDispcore_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPtf,xPt,nPhi,xPhi) ;
 525     tmp->Sumw2() ; fOutputContainer->Add(tmp) ;
 526     tmp=new TH2F(Form("hPhotPhi%sDisp2_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPtf,xPt,nPhi,xPhi) ;
 527     tmp->Sumw2() ; fOutputContainer->Add(tmp) ;
 528     tmp=new TH2F(Form("hPhotPhi%sDisp2core_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPtf,xPt,nPhi,xPhi) ;
 529     tmp->Sumw2() ; fOutputContainer->Add(tmp) ;
 530     tmp=new TH2F(Form("hPhotPhi%sCPV_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPtf,xPt,nPhi,xPhi) ;
 531     tmp->Sumw2() ; fOutputContainer->Add(tmp) ;
 532     tmp=new TH2F(Form("hPhotPhi%sCPVcore_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPtf,xPt,nPhi,xPhi) ;
 533     tmp->Sumw2() ; fOutputContainer->Add(tmp) ;
 534     tmp=new TH2F(Form("hPhotPhi%sCPV2_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPtf,xPt,nPhi,xPhi) ;
 535     tmp->Sumw2() ; fOutputContainer->Add(tmp) ;
 536     tmp=new TH2F(Form("hPhotPhi%sCPV2core_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPtf,xPt,nPhi,xPhi) ;
 537     tmp->Sumw2() ; fOutputContainer->Add(tmp) ;
 538     tmp=new TH2F(Form("hPhotPhi%sBoth_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPtf,xPt,nPhi,xPhi) ;
 539     tmp->Sumw2() ; fOutputContainer->Add(tmp) ;
 540     tmp=new TH2F(Form("hPhotPhi%sBothcore_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPtf,xPt,nPhi,xPhi) ;
 541     tmp->Sumw2() ; fOutputContainer->Add(tmp) ;
 542     tmp=new TH2F(Form("hPhotPhi%sBoth2_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPtf,xPt,nPhi,xPhi) ;
 543     tmp->Sumw2() ; fOutputContainer->Add(tmp) ;
 544     tmp=new TH2F(Form("hPhotPhi%sBoth2core_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}" ,nPtf,xPt,nPhi,xPhi) ;
 545     tmp->Sumw2() ; fOutputContainer->Add(tmp) ;
 546
 547
 548     //Pions for flow - with weight 1/Nclu
 549     TH3F * tmp3 = 0;
 550     tmp3 = new TH3F(Form("hMassPt%sAll_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 551     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 552     tmp3 = new TH3F(Form("hMassPt%sAllcore_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 553     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 554     tmp3 = new TH3F(Form("hMassPt%sCPV_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 555     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 556     tmp3 = new TH3F(Form("hMassPt%sCPVcore_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 557     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 558     tmp3 = new TH3F(Form("hMassPt%sCPV2_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 559     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 560     tmp3 = new TH3F(Form("hMassPt%sCPV2core_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 561     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 562     tmp3 = new TH3F(Form("hMassPt%sDisp_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 563     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 564     tmp3 = new TH3F(Form("hMassPt%sDispcore_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 565     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 566     tmp3 = new TH3F(Form("hMassPt%sDisp2_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 567     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 568     tmp3 = new TH3F(Form("hMassPt%sDisp2core_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 569     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 570     tmp3 = new TH3F(Form("hMassPt%sBoth_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 571     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 572     tmp3 = new TH3F(Form("hMassPt%sBothcore_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 573     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 574     tmp3 = new TH3F(Form("hMassPt%sBoth2_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 575     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 576     tmp3 = new TH3F(Form("hMassPt%sBoth2core_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 577     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 578
 579 //mixed
 580     tmp3 = new TH3F(Form("hMiMassPt%sAll_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 581     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 582     tmp3 = new TH3F(Form("hMiMassPt%sAllcore_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 583     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 584     tmp3 = new TH3F(Form("hMiMassPt%sCPV_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 585     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 586     tmp3 = new TH3F(Form("hMiMassPt%sCPVcore_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 587     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 588     tmp3 = new TH3F(Form("hMiMassPt%sCPV2_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 589     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 590     tmp3 = new TH3F(Form("hMiMassPt%sCPV2core_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 591     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 592     tmp3 = new TH3F(Form("hMiMassPt%sDisp_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 593     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 594     tmp3 = new TH3F(Form("hMiMassPt%sDispcore_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 595     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 596     tmp3 = new TH3F(Form("hMiMassPt%sDisp2_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 597     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 598     tmp3 = new TH3F(Form("hMiMassPt%sDisp2core_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 599     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 600     tmp3 = new TH3F(Form("hMiMassPt%sBoth_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 601     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 602     tmp3 = new TH3F(Form("hMiMassPt%sBothcore_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 603     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 604     tmp3 = new TH3F(Form("hMiMassPt%sBoth2_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 605     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 606     tmp3 = new TH3F(Form("hMiMassPt%sBoth2core_cen%d",phiTitle,cent),"(M,p_{T},d#phi)_{#gamma#gamma}",nMm,xM,nPtf,xPt,nPhi,xPhi);
 607     tmp3->Sumw2() ; fOutputContainer->Add(tmp3) ;
 608
 609   }
 610
 611   fOutputContainer->Add(new TH2F("hMCPi0M11","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 612   fOutputContainer->Add(new TH2F("hMCPi0M22","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 613   fOutputContainer->Add(new TH2F("hMCPi0M33","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 614   fOutputContainer->Add(new TH2F("hMCPi0M12","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 615   fOutputContainer->Add(new TH2F("hMCPi0M13","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 616   fOutputContainer->Add(new TH2F("hMCPi0M23","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 617
 618   //MC
 619   for(Int_t cent=0; cent<6; cent++){
 620     snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
 621     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
 622     snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
 623     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
 624     snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
 625     fOutputContainer->Add(new TH1F("hMC_rap_eta","Rapidity eta",200,-1.,1.)) ;
 626     snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
 627     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
 628     snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
 629     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
 630     snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
 631     fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
 632     snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
 633     fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
 634     snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
 635     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
 636     snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
 637     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 638     snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
 639     fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
 640     snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
 641     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 642     snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
 643     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 644   }
 645
 646   PostData(1, fOutputContainer);
 647
 648 }
 649
 650 //________________________________________________________________________
 651 void AliAnalysisTaskPi0DiffEfficiency::UserExec(Option_t *)
 652 {
 653   // Main loop, called for each event
 654   // Analyze ESD/AOD
 655
 656   FillHistogram("hSelEvents",0.5) ;
 657
 658   AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
 659   if (!event) {
 660     Printf("ERROR: Could not retrieve event");
 661     PostData(1, fOutputContainer);
 662     return;
 663   }
 664
 665   FillHistogram("hSelEvents",1.5) ;
 666   AliAODHeader *header = event->GetHeader() ;
 667
 668   // Checks if we have a primary vertex
 669   // Get primary vertices form ESD
 670   const AliAODVertex *esdVertex5 = event->GetPrimaryVertex();
 671
 672  // don't rely on ESD vertex, assume (0,0,0)
 673   Double_t vtx0[3] ={0.,0.,0.};
 674
 675
 676   FillHistogram("hZvertex",esdVertex5->GetZ());
 677   if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
 678     PostData(1, fOutputContainer);
 679     return;
 680   }
 681   FillHistogram("hSelEvents",2.5) ;
 682
 683   //Vtx class z-bin
 684   //  Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
 685   //  if(zvtx<0)zvtx=0 ;
 686   //  if(zvtx>9)zvtx=9 ;
 687   Int_t zvtx=0 ;
 688
 689 //  fCentrality=header->GetCentralityP()->GetCentralityPercentile("V0M"); // returns the centrality percentile,
 690 //                                                          //a float from 0 to 100 (or to the trigger efficiency)
 691    fCentrality=header->GetZDCN2Energy() ;
 692
 693   if( fCentrality < 0. || fCentrality>80.){
 694     PostData(1, fOutputContainer);
 695     return;
 696   }
 697   FillHistogram("hSelEvents",3.5) ;
 698   Float_t bins[7]={0.,5.,10.,20.,40.,60.,80.} ;
 699   fCenBin=0 ;
 700   while(fCenBin<6 && fCentrality > bins[fCenBin+1])
 701     fCenBin++ ;
 702
 703   const Int_t nMixEvents[6]={4,4,5,10,20,20} ;
 704
 705
 706   //reaction plane
 707   fRPfull= header->GetZDCN1Energy() ;
 708   if(fRPfull==999){ //reaction plain was not defined
 709     PostData(1, fOutputContainer);
 710     return;
 711   }
 712
 713   FillHistogram("hSelEvents",4.5) ;
 714   //All event selections done
 715   FillHistogram("hCentrality",fCentrality) ;
 716   //Reaction plain is defined in the range (-pi/2;pi/2)
 717   //We have 10 bins
 718   Int_t irp=Int_t(10.*fRPfull/TMath::Pi());
 719   if(irp>9)irp=9 ;
 720
 721   if(!fPHOSEvents[zvtx][fCenBin][irp])
 722     fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
 723   TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
 724
 725   // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
 726   if(fEventCounter == 0) {
 727     for(Int_t mod=0; mod<5; mod++) {
 728       const TGeoHMatrix* m =header->GetPHOSMatrix(mod) ;
 729       fPHOSGeo->SetMisalMatrix(m,mod) ;
 730       Printf("PHOS geo matrix for module # %d is set: %p\n", mod,m);
 731     }
 732     fEventCounter++ ;
 733   }
 734
 735   TClonesArray *mcArray = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
 736
 737   ProcessMC() ;
 738
 739   if(fPHOSEvent1){
 740     fPHOSEvent1->Clear() ;
 741     fPHOSEvent2->Clear() ;
 742   }
 743   else{
 744     fPHOSEvent1 = new TClonesArray("AliCaloPhoton",200) ;
 745     fPHOSEvent2 = new TClonesArray("AliCaloPhoton",200) ;
 746   }
 747
 748   TClonesArray * clustersEmb = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
 749   AliAODCaloCells * cellsEmb = (AliAODCaloCells *)event->FindListObject("EmbeddedPHOScells") ;
 750   TClonesArray * clustersOld = event->GetCaloClusters() ;
 751   AliAODCaloCells * cellsOld = event->GetPHOSCells() ;
 752 //  TClonesArray *mcArray = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
 753
 754
 755   TVector3 vertex(vtx0);
 756   char key[55] ;
 757   //Before Embedding
 758   Int_t multClustOld = clustersOld->GetEntriesFast();
 759   Int_t multClustEmb = clustersEmb->GetEntriesFast();
 760   Int_t inPHOSold=0 ;
 761   Int_t inPHOSemb=0 ;
 762   for (Int_t i=0; i<multClustOld; i++) {
 763     AliAODCaloCluster *clu = (AliAODCaloCluster*)clustersOld->At(i);
 764     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
 765
 766     Bool_t survive=kFALSE ;
 767     for(Int_t ii=0;(ii<multClustEmb)&&(!survive);ii++){
 768        AliAODCaloCluster *clu2 = (AliAODCaloCluster*)clustersEmb->At(ii);
 769        survive=IsSameCluster(clu,clu2);
 770     }
 771
 772
 773     Float_t  position[3];
 774     clu->GetPosition(position);
 775     TVector3 global(position) ;
 776     Int_t relId[4] ;
 777     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 778     Int_t mod  = relId[0] ;
 779     Int_t cellX = relId[2];
 780     Int_t cellZ = relId[3] ;
 781     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
 782       continue ;
 783     if(clu->GetNCells()<3)
 784       continue ;
 785
 786     snprintf(key,55,"hCluM%d",mod) ;
 787     FillHistogram(key,cellX,cellZ,1.);
 788
 789     TLorentzVector pv1 ;
 790     clu->GetMomentum(pv1 ,vtx0);
 791
 792     pv1*=Scale(pv1.E()) ;
 793
 794     if(inPHOSold>=fPHOSEvent1->GetSize()){
 795       fPHOSEvent1->Expand(inPHOSold+50) ;
 796     }
 797     AliCaloPhoton * ph = new((*fPHOSEvent1)[inPHOSold]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
 798     ph->SetModule(mod) ;
 799     AliPHOSAodCluster cluPHOS1(*clu);
 800     cluPHOS1.Recalibrate(fPHOSCalibData,cellsOld); // modify the cell energies
 801     Double_t ecore=CoreEnergy(&cluPHOS1) ;
 802     ecore*=Scale(ecore) ;
 803     pv1*= ecore/pv1.E() ;
 804     ph->SetMomV2(&pv1) ;
 805     ph->SetNCells(clu->GetNCells());
 806     Double_t m02=0.,m20=0.;
 807     EvalLambdas(&cluPHOS1,0,m02, m20);
 808     ph->SetDispBit(TestLambda(clu->E(),m20,m02)) ;
 809     EvalLambdas(&cluPHOS1,1,m02, m20);
 810     ph->SetDisp2Bit(TestLambda2(clu->E(),m20,m02)) ;
 811
 812     ph->SetCPVBit(clu->GetEmcCpvDistance()>2.) ;
 813     ph->SetCPV2Bit(clu->GetEmcCpvDistance()>4.) ;
 814     if(!survive) //this cluster found in list after embedding, skipping it
 815       ph->SetTagged(1) ;
 816     ph->SetPhoton(clu->GetNExMax()<2); // Remember, if it is unfolded
 817     ph->SetWeight(1.) ; //All weights for real particles ==1.
 818
 819     if(!survive){
 820       Double_t distBC=clu->GetDistanceToBadChannel();
 821       if(distBC>2.)
 822         FillHistogram(Form("hPhotAll_DistBad2_cen%d",fCenBin),ph->Pt(),-1.) ;
 823         if(distBC>4.)
 824           FillHistogram(Form("hPhotAll_DistBad4_cen%d",fCenBin),ph->Pt(),-1.) ;
 825           if(distBC>6.)
 826             FillHistogram(Form("hPhotAll_DistBad6_cen%d",fCenBin),ph->Pt(),-1.) ;
 827     }
 828     inPHOSold++ ;
 829   }
 830
 831   for (Int_t i=0; i<multClustEmb; i++) {
 832     AliAODCaloCluster *clu = (AliAODCaloCluster*)clustersEmb->At(i);
 833     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
 834
 835     Bool_t survive=kFALSE ;
 836     for(Int_t ii=0;(ii<multClustOld)&&(!survive);ii++){
 837        AliAODCaloCluster *clu2 = (AliAODCaloCluster*)clustersOld->At(ii);
 838        survive=IsSameCluster(clu,clu2);
 839     }
 840
 841     Float_t  position[3];
 842     clu->GetPosition(position);
 843     TVector3 global(position) ;
 844     Int_t relId[4] ;
 845     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 846     Int_t mod  = relId[0] ;
 847     Int_t cellX = relId[2];
 848     Int_t cellZ = relId[3] ;
 849     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
 850       continue ;
 851     if(clu->GetNCells()<3)
 852       continue ;
 853
 854     snprintf(key,55,"hCluM%d",mod) ;
 855     FillHistogram(key,cellX,cellZ,1.);
 856
 857     TLorentzVector pv1 ;
 858     clu->GetMomentum(pv1 ,vtx0);
 859
 860     pv1*=Scale(pv1.E()) ;
 861
 862     if(inPHOSemb>=fPHOSEvent2->GetSize()){
 863       fPHOSEvent2->Expand(inPHOSemb+50) ;
 864     }
 865     AliCaloPhoton * ph = new((*fPHOSEvent2)[inPHOSemb]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
 866     ph->SetModule(mod) ;
 867     AliPHOSAodCluster cluPHOS1(*clu);
 868     cluPHOS1.Recalibrate(fPHOSCalibData,cellsEmb); // modify the cell energies
 869     Double_t ecore=CoreEnergy(&cluPHOS1) ;
 870     ecore*=Scale(ecore) ;
 871     pv1*= ecore/pv1.E() ;
 872     ph->SetMomV2(&pv1) ;
 873     ph->SetNCells(clu->GetNCells());
 874     Double_t m02=0.,m20=0.;
 875     EvalLambdas(&cluPHOS1,0,m02, m20);
 876     ph->SetDispBit(TestLambda(clu->E(),m20,m02)) ;
 877     EvalLambdas(&cluPHOS1,1,m02, m20);
 878     ph->SetDisp2Bit(TestLambda2(clu->E(),m20,m02)) ;
 879
 880     ph->SetCPVBit(clu->GetEmcCpvDistance()>2.) ;
 881     ph->SetCPV2Bit(clu->GetEmcCpvDistance()>4.) ;
 882     if(!survive) //this cluster found in list after embedding, skipping it
 883       ph->SetTagged(1) ;
 884     ph->SetPhoton(clu->GetNExMax()<2); // Remember, if it is unfolded
 885
 886     //Set weight for embedded particles
 887     Double_t w=1. ;
 888     Int_t iprim = clu->GetLabel() ;
 889     if(iprim<mcArray->GetEntriesFast() && iprim>-1){
 890       AliAODMCParticle* particle =  (AliAODMCParticle*) mcArray->At(iprim);
 891       iprim=particle->GetMother() ;
 892       while(iprim>-1){
 893         particle =  (AliAODMCParticle*) mcArray->At(iprim);
 894         iprim=particle->GetMother() ;
 895       }
 896       if(particle->GetPdgCode()==111){
 897         Double_t pt = particle->Pt() ;
 898         w=PrimaryWeight(pt) ;
 899       }
 900     }
 901     ph->SetWeight(w) ;
 902
 903
 904     if(!survive){
 905       Double_t distBC=clu->GetDistanceToBadChannel();
 906       if(distBC>2.)
 907         FillHistogram(Form("hPhotAll_DistBad2_cen%d",fCenBin),ph->Pt()) ;
 908         if(distBC>4.)
 909           FillHistogram(Form("hPhotAll_DistBad4_cen%d",fCenBin),ph->Pt()) ;
 910           if(distBC>6.)
 911             FillHistogram(Form("hPhotAll_DistBad6_cen%d",fCenBin),ph->Pt()) ;
 912     }
 913
 914     inPHOSemb++ ;
 915   }
 916
 917   //Single photon
 918   for (Int_t i1=0; i1<inPHOSold; i1++) {
 919     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent1->At(i1) ;
 920     if(!ph1->IsTagged())
 921       continue ;
 922
 923     Double_t dphi=ph1->Phi()-fRPfull ;
 924     while(dphi<0)dphi+=TMath::Pi() ;
 925     while(dphi>TMath::Pi())dphi-=TMath::Pi() ;
 926     Double_t pt = ph1->Pt() ;
 927     Double_t ptV2=ph1->GetMomV2()->Pt() ;
 928 //always 1!    Double_t w=ph1->GetWeight() ;
 929
 930     FillHistogram(Form("hPhotPhiTPCAll_cen%d",fCenBin),pt,dphi,-1.) ;
 931     FillHistogram(Form("hPhotPhiTPCAllcore_cen%d",fCenBin),ptV2,dphi,-1.) ;
 932
 933     FillHistogram(Form("hPhotAll_cen%d",fCenBin),pt,-1.) ;
 934     FillHistogram(Form("hNegPhotAll_cen%d",fCenBin),pt,-1.) ;
 935     FillHistogram(Form("hPhotAllcore_cen%d",fCenBin),ptV2,-1.) ;
 936     FillHistogram(Form("hNegPhotAllcore_cen%d",fCenBin),ptV2,-1.) ;
 937     if(ph1->IsPhoton()){
 938       FillHistogram(Form("hPhotAllwou_cen%d",fCenBin),pt,-1.) ;
 939     }
 940     if(ph1->IsCPVOK() ){
 941       FillHistogram(Form("hPhotCPV_cen%d",fCenBin),pt,-1.) ;
 942       FillHistogram(Form("hNegPhotCPV_cen%d",fCenBin),pt,-1.) ;
 943       FillHistogram(Form("hPhotCPVcore_cen%d",fCenBin),ptV2,-1.) ;
 944       FillHistogram(Form("hNegPhotCPVcore_cen%d",fCenBin),ptV2,-1.) ;
 945       FillHistogram(Form("hPhotPhiTPCCPV_cen%d",fCenBin),pt,dphi,-1.) ;
 946       FillHistogram(Form("hPhotPhiTPCCPVcore_cen%d",fCenBin),ptV2,dphi,-1.) ;
 947     }
 948     if(ph1->IsCPV2OK() ){
 949       FillHistogram(Form("hPhotCPV2_cen%d",fCenBin),pt,-1.) ;
 950       FillHistogram(Form("hPhotCPV2core_cen%d",fCenBin),ptV2,-1.) ;
 951       FillHistogram(Form("hNegPhotCPV2_cen%d",fCenBin),pt,-1.) ;
 952       FillHistogram(Form("hPhotPhiTPCCPV2_cen%d",fCenBin),pt,dphi,-1.) ;
 953       FillHistogram(Form("hPhotPhiTPCCPV2core_cen%d",fCenBin),ptV2,dphi,-1.) ;
 954     }
 955     if(ph1->IsDisp2OK()){
 956       FillHistogram(Form("hPhotDisp2_cen%d",fCenBin),pt,-1.) ;
 957       FillHistogram(Form("hPhotDisp2core_cen%d",fCenBin),ptV2,-1.) ;
 958       FillHistogram(Form("hNegPhotDisp2_cen%d",fCenBin),pt,-1.) ;
 959       FillHistogram(Form("hNegPhotDisp2core_cen%d",fCenBin),ptV2,-1.) ;
 960
 961       FillHistogram(Form("hPhotPhiTPCDisp2_cen%d",fCenBin),pt,dphi,-1.) ;
 962       FillHistogram(Form("hPhotPhiTPCDisp2core_cen%d",fCenBin),ptV2,dphi,-1.) ;
 963       if(ph1->IsCPVOK()){
 964         FillHistogram(Form("hPhotBoth2_cen%d",fCenBin),pt,-1.) ;
 965         FillHistogram(Form("hNegPhotBoth2_cen%d",fCenBin),pt,-1.) ;
 966         FillHistogram(Form("hPhotBoth2core_cen%d",fCenBin),ptV2,-1.) ;
 967         FillHistogram(Form("hNegPhotBoth2core_cen%d",fCenBin),ptV2,-1.) ;
 968         FillHistogram(Form("hPhotPhiTPCBoth2_cen%d",fCenBin),pt,dphi,-1.) ;
 969         FillHistogram(Form("hPhotPhiTPCBoth2core_cen%d",fCenBin),ptV2,dphi,-1.) ;
 970       }
 971     }
 972     if(ph1->IsDispOK()){
 973       FillHistogram(Form("hPhotDisp_cen%d",fCenBin),pt,-1.) ;
 974       FillHistogram(Form("hNegPhotDisp_cen%d",fCenBin),pt,-1.) ;
 975       FillHistogram(Form("hPhotDispcore_cen%d",fCenBin),ptV2,-1.) ;
 976       FillHistogram(Form("hNegPhotDispcore_cen%d",fCenBin),ptV2,-1.) ;
 977       if(ph1->IsPhoton()){
 978         FillHistogram(Form("hPhotDispwou_cen%d",fCenBin),pt,-1.) ;
 979       }
 980
 981       FillHistogram(Form("hPhotPhiTPCDisp_cen%d",fCenBin),pt,dphi,-1.) ;
 982       FillHistogram(Form("hPhotPhiTPCDispcore_cen%d",fCenBin),ptV2,dphi,-1.) ;
 983       if(ph1->IsCPVOK()){
 984         FillHistogram(Form("hPhotBoth_cen%d",fCenBin),pt,-1.) ;
 985         FillHistogram(Form("hNegPhotBoth_cen%d",fCenBin),pt,-1.) ;
 986         FillHistogram(Form("hPhotBothcore_cen%d",fCenBin),ptV2,-1.) ;
 987         FillHistogram(Form("hNegPhotBothcore_cen%d",fCenBin),ptV2,-1.) ;
 988
 989         FillHistogram(Form("hPhotPhiTPCBoth_cen%d",fCenBin),pt,dphi,-1.) ;
 990         FillHistogram(Form("hPhotPhiTPCBothcore_cen%d",fCenBin),ptV2,dphi,-1.) ;
 991       }
 992     }
 993   } // end of loop i1
 994
 995   for (Int_t i1=0; i1<inPHOSemb; i1++) {
 996     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
 997     if(!ph1->IsTagged())
 998       continue ;
 999
1000     Double_t dphi=ph1->Phi()-fRPfull ;
1001     while(dphi<0)dphi+=TMath::Pi() ;
1002     while(dphi>TMath::Pi())dphi-=TMath::Pi() ;
1003     Double_t pt = ph1->Pt() ;
1004     Double_t ptV2=ph1->GetMomV2()->Pt() ;
1005     Double_t w=ph1->GetWeight() ;
1006
1007     FillHistogram(Form("hPhotPhiTPCAll_cen%d",fCenBin),pt,dphi,w) ;
1008     FillHistogram(Form("hPhotPhiTPCAllcore_cen%d",fCenBin),ptV2,dphi,w) ;
1009
1010     FillHistogram(Form("hPhotAll_cen%d",fCenBin),pt,w) ;
1011     FillHistogram(Form("hPhotAllcore_cen%d",fCenBin),ptV2,w) ;
1012     if(ph1->IsPhoton()){
1013       FillHistogram(Form("hPhotAllwou_cen%d",fCenBin),pt,w) ;
1014     }
1015     if(ph1->IsCPVOK() ){
1016       FillHistogram(Form("hPhotCPV_cen%d",fCenBin),pt,w) ;
1017       FillHistogram(Form("hPhotCPVcore_cen%d",fCenBin),ptV2,w) ;
1018
1019       FillHistogram(Form("hPhotPhiTPCCPV_cen%d",fCenBin),pt,dphi,w) ;
1020       FillHistogram(Form("hPhotPhiTPCCPVcore_cen%d",fCenBin),ptV2,dphi,w) ;
1021     }
1022     if(ph1->IsCPV2OK() ){
1023       FillHistogram(Form("hPhotCPV2_cen%d",fCenBin),pt,w) ;
1024       FillHistogram(Form("hPhotCPV2core_cen%d",fCenBin),ptV2,w) ;
1025       FillHistogram(Form("hPhotPhiTPCCPV2_cen%d",fCenBin),pt,dphi,w) ;
1026       FillHistogram(Form("hPhotPhiTPCCPV2core_cen%d",fCenBin),ptV2,dphi,w) ;
1027     }
1028     if(ph1->IsDisp2OK()){
1029       FillHistogram(Form("hPhotDisp2_cen%d",fCenBin),pt,w) ;
1030       FillHistogram(Form("hPhotDisp2core_cen%d",fCenBin),ptV2,w) ;
1031
1032       FillHistogram(Form("hPhotPhiTPCDisp2_cen%d",fCenBin),pt,dphi,w) ;
1033       FillHistogram(Form("hPhotPhiTPCDisp2core_cen%d",fCenBin),ptV2,dphi,w) ;
1034       if(ph1->IsCPVOK()){
1035         FillHistogram(Form("hPhotBoth2_cen%d",fCenBin),pt,w) ;
1036         FillHistogram(Form("hPhotBoth2core_cen%d",fCenBin),ptV2,w) ;
1037
1038         FillHistogram(Form("hPhotPhiTPCBoth2_cen%d",fCenBin),pt,dphi,w) ;
1039         FillHistogram(Form("hPhotPhiTPCBoth2core_cen%d",fCenBin),ptV2,dphi,w) ;
1040       }
1041     }
1042     if(ph1->IsDispOK()){
1043       FillHistogram(Form("hPhotDisp_cen%d",fCenBin),pt,w) ;
1044       FillHistogram(Form("hPhotDispcore_cen%d",fCenBin),ptV2,w) ;
1045       if(ph1->IsPhoton()){
1046         FillHistogram(Form("hPhotDispwou_cen%d",fCenBin),pt,w) ;
1047       }
1048
1049       FillHistogram(Form("hPhotPhiTPCBoth_cen%d",fCenBin),pt,dphi,w) ;
1050       FillHistogram(Form("hPhotPhiTPCBothcore_cen%d",fCenBin),ptV2,dphi,w) ;
1051       if(ph1->IsCPVOK()){
1052         FillHistogram(Form("hPhotBoth_cen%d",fCenBin),pt,w) ;
1053         FillHistogram(Form("hPhotBothcore_cen%d",fCenBin),ptV2,w) ;
1054
1055         FillHistogram(Form("hPhotPhiTPCBoth_cen%d",fCenBin),pt,dphi,w) ;
1056         FillHistogram(Form("hPhotPhiTPCBothcore_cen%d",fCenBin),ptV2,dphi,w) ;
1057       }
1058     } // end of loop i2
1059   } // end of loop i1
1060
1061
1062   const Double_t prob[10]={0.1,0.2,0.3,1.,1.,1.,1.,1.,1.,1.} ; //Probabilities to accept Tagged+Bg pair
1063
1064
1065   // Fill Real disribution:
1066   // Disappeared clusters enter with negative contribution
1067   // In addition fill control histogam with Real before embedding
1068   for (Int_t i1=0; i1<inPHOSold-1; i1++) {
1069     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent1->At(i1) ;
1070     Double_t w1 = ph1->GetWeight() ;
1071     for (Int_t i2=i1+1; i2<inPHOSold; i2++) {
1072       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent1->At(i2) ;
1073
1074       TLorentzVector p12  = *ph1  + *ph2;
1075       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1076       Double_t w2 = ph2->GetWeight() ;
1077       Double_t w = TMath::Sqrt(w1*w2) ;
1078
1079       Double_t dphi=p12.Phi()-fRPfull ;
1080       while(dphi<0)dphi+=TMath::Pi() ;
1081       while(dphi>TMath::Pi())dphi-=TMath::Pi() ;
1082       Double_t m=p12.M() ;
1083       Double_t mV2=pv12.M() ;
1084       Double_t pt = p12.Pt() ;
1085       Double_t ptV2 = pv12.Pt() ;
1086
1087
1088       //Fill Controll histogram: Real before embedding
1089       FillHistogram(Form("hOldMassPtAll_cen%d",fCenBin),m,pt,-w) ;
1090       FillHistogram(Form("hOldMassPtAllcore_cen%d",fCenBin),mV2,ptV2,-w) ;
1091       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1092         FillHistogram(Form("hOldMassPtCPV_cen%d",fCenBin),m ,pt,-w) ;
1093         FillHistogram(Form("hOldMassPtCPV_cen%d",fCenBin),mV2, ptV2,-w) ;
1094       }
1095       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1096         FillHistogram(Form("hOldMassPtCPV2_cen%d",fCenBin),m ,pt,-w) ;
1097       }
1098       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1099         FillHistogram(Form("hOldMassPtDisp2_cen%d",fCenBin),m ,pt,-w) ;
1100         FillHistogram(Form("hOldMassPtDisp2core_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1101         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1102           FillHistogram(Form("hOldMassPtBoth2_cen%d",fCenBin),m ,pt,-w) ;
1103           FillHistogram(Form("hOldMassPtBoth2core_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1104         }
1105       }
1106       if(ph1->IsDispOK() && ph2->IsDispOK()){
1107         FillHistogram(Form("hOldMassPtDisp_cen%d",fCenBin),m ,pt,-w) ;
1108         FillHistogram(Form("hOldMassPtDispcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1109         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1110           FillHistogram(Form("hOldMassPtBoth_cen%d",fCenBin),m ,pt,-w) ;
1111           FillHistogram(Form("hOldMassPtBothcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1112         }
1113       }
1114
1115       //Now fill main histograms with negative contributions
1116       if(!(ph1->IsTagged() || ph2->IsTagged()) )
1117         continue ;
1118       if(!ph1->IsTagged() || !ph2->IsTagged()){ //Tagged + Bg combination
1119         if(gRandom->Uniform()>prob[fCenBin])
1120           continue ;
1121       }
1122       FillHistogram(Form("hMassPtAll_cen%d",fCenBin),m ,pt,-w) ;
1123       FillHistogram(Form("hMassPtAllcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1124       if(ph1->IsPhoton()&&ph2->IsPhoton()){
1125         FillHistogram(Form("hMassPtAllwou_cen%d",fCenBin),m,pt,-w) ;
1126       }
1127
1128       FillHistogram(Form("hMassPtTPCAll_cen%d",fCenBin),m,pt,dphi,-w) ;
1129       FillHistogram(Form("hMassPtTPCAllcore_cen%d",fCenBin),mV2,ptV2,dphi,-w) ;
1130
1131       FillHistogram(Form("hNegMassPtAll_cen%d",fCenBin),m ,pt,-w) ;
1132       FillHistogram(Form("hNegMassPtAllcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1133
1134       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1135         FillHistogram(Form("hMassPtCPV_cen%d",fCenBin),m ,pt,-w) ;
1136         FillHistogram(Form("hMassPtCPVcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1137
1138         FillHistogram(Form("hMassPtTPCCPV_cen%d",fCenBin),m,pt,dphi,-w) ;
1139         FillHistogram(Form("hMassPtTPCCPVcore_cen%d",fCenBin),mV2,ptV2,dphi,-w) ;
1140
1141         FillHistogram(Form("hNegMassPtCPV_cen%d",fCenBin),m ,pt,-w) ;
1142         FillHistogram(Form("hNegMassPtCPVcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1143       }
1144       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1145         FillHistogram(Form("hMassPtCPV2_cen%d",fCenBin),m ,pt,-w) ;
1146         FillHistogram(Form("hNegMassPtCPV2_cen%d",fCenBin),m ,pt,-w) ;
1147
1148         FillHistogram(Form("hMassPtTPCCPV2_cen%d",fCenBin),m,pt,dphi,-w) ;
1149         FillHistogram(Form("hMassPtTPCCPV2core_cen%d",fCenBin),mV2,ptV2,dphi,-w) ;
1150       }
1151
1152       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1153         FillHistogram(Form("hMassPtDisp2_cen%d",fCenBin),m ,pt,-w) ;
1154         FillHistogram(Form("hMassPtDisp2core_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1155
1156         FillHistogram(Form("hMassPtTPCDisp2_cen%d",fCenBin),m,pt,dphi,-w) ;
1157         FillHistogram(Form("hMassPtTPCDisp2core_cen%d",fCenBin),mV2,ptV2,dphi,-w) ;
1158
1159         FillHistogram(Form("hNegMassPtDisp2_cen%d",fCenBin),m ,pt,-w) ;
1160         FillHistogram(Form("hNegMassPtDisp2core_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1161         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1162           FillHistogram(Form("hMassPtBoth2_cen%d",fCenBin),m ,pt,-w) ;
1163           FillHistogram(Form("hMassPtBoth2core_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1164           FillHistogram(Form("hNegMassPtBoth2_cen%d",fCenBin),m ,pt,-w) ;
1165           FillHistogram(Form("hNegMassPtBoth2core_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1166
1167           FillHistogram(Form("hMassPtTPCBoth2_cen%d",fCenBin),m,pt,dphi,-w) ;
1168           FillHistogram(Form("hMassPtTPCBoth2core_cen%d",fCenBin),mV2,ptV2,dphi,-w) ;
1169         }
1170       }
1171       if(ph1->IsDispOK() && ph2->IsDispOK()){
1172         FillHistogram(Form("hMassPtDisp_cen%d",fCenBin),m ,pt,-w) ;
1173         FillHistogram(Form("hMassPtDispcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1174         if(ph1->IsPhoton()&&ph2->IsPhoton()){
1175           FillHistogram(Form("hMassPtDispwou_cen%d",fCenBin),m,pt,-w) ;
1176         }
1177         FillHistogram(Form("hNegMassPtDisp_cen%d",fCenBin),m ,pt,-w) ;
1178         FillHistogram(Form("hNegMassPtDispcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1179
1180         FillHistogram(Form("hMassPtTPCDisp_cen%d",fCenBin),m,pt,dphi,-w) ;
1181         FillHistogram(Form("hMassPtTPCDispcore_cen%d",fCenBin),mV2,ptV2,dphi,-w) ;
1182         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1183           FillHistogram(Form("hMassPtBoth_cen%d",fCenBin),m ,pt,-w) ;
1184           FillHistogram(Form("hMassPtBothcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1185           FillHistogram(Form("hNegMassPtBoth_cen%d",fCenBin),m ,pt,-w) ;
1186           FillHistogram(Form("hNegMassPtBothcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1187
1188           FillHistogram(Form("hMassPtTPCBoth_cen%d",fCenBin),m,pt,dphi,-w) ;
1189           FillHistogram(Form("hMassPtTPCBothcore_cen%d",fCenBin),mV2,ptV2,dphi,-w) ;
1190         }
1191       }
1192     } // end of loop i2
1193   } // end of loop i1
1194
1195
1196   // Further fill Real disribution
1197   // now with positive contribution from new clusters
1198   // ass well fill controll histogram
1199   for (Int_t i1=0; i1<inPHOSemb-1; i1++) {
1200     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
1201     Double_t w1 = ph1->GetWeight() ;
1202     for (Int_t i2=i1+1; i2<inPHOSemb; i2++) {
1203       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent2->At(i2) ;
1204
1205       TLorentzVector p12  = *ph1  + *ph2;
1206       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1207       Double_t m=p12.M() ;
1208       Double_t mV2=pv12.M() ;
1209       Double_t pt = p12.Pt() ;
1210       Double_t ptV2 = pv12.Pt() ;
1211       Double_t w2 = ph2->GetWeight() ;
1212       Double_t w = TMath::Sqrt(w1*w2) ;
1213
1214       // Controll histogram: Real after embedding
1215       FillHistogram(Form("hNewMassPtAll_cen%d",fCenBin),m ,pt,w) ;
1216       FillHistogram(Form("hNewMassPtAllcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1217       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1218         FillHistogram(Form("hNewMassPtCPV_cen%d",fCenBin),m ,pt,w) ;
1219         FillHistogram(Form("hNewMassPtCPVcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1220       }
1221       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1222         FillHistogram(Form("hNewMassPtCPV2_cen%d",fCenBin),m ,pt,w) ;
1223       }
1224       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1225         FillHistogram(Form("hNewMassPtDisp2_cen%d",fCenBin),m ,pt,w) ;
1226         FillHistogram(Form("hNewMassPtDisp2core_cen%d",fCenBin),mV2 ,ptV2,w) ;
1227         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1228           FillHistogram(Form("hNewMassPtBoth2_cen%d",fCenBin),m ,pt,w) ;
1229           FillHistogram(Form("hNewMassPtBoth2core_cen%d",fCenBin),mV2,ptV2,w) ;
1230         }
1231       }
1232       if(ph1->IsDispOK() && ph2->IsDispOK()){
1233         FillHistogram(Form("hNewMassPtDisp_cen%d",fCenBin),m ,pt,w) ;
1234         FillHistogram(Form("hNewMassPtDispcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1235         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1236           FillHistogram(Form("hNewMassPtBoth_cen%d",fCenBin),m ,pt,w) ;
1237           FillHistogram(Form("hNewMassPtBothcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1238         }
1239       }
1240
1241       //Now fill main histogamm
1242       //new clusters with positive contribution
1243       if(!(ph1->IsTagged() || ph2->IsTagged()) )
1244         continue ;
1245       if(!ph1->IsTagged() || !ph2->IsTagged()){ //Tagged + Bg combination
1246         if(gRandom->Uniform()>prob[fCenBin])
1247           continue ;
1248       }
1249
1250       Double_t dphi=p12.Phi()-fRPfull ;
1251       while(dphi<0)dphi+=TMath::Pi() ;
1252       while(dphi>TMath::Pi())dphi-=TMath::Pi() ;
1253
1254       FillHistogram(Form("hMassPtAll_cen%d",fCenBin),m ,pt,w) ;
1255       FillHistogram(Form("hMassPtAllcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1256       if(ph1->IsPhoton()&&ph2->IsPhoton()){
1257         FillHistogram(Form("hMassPtAllwou_cen%d",fCenBin),m,pt,w) ;
1258       }
1259
1260       FillHistogram(Form("hMassPtTPCAll_cen%d",fCenBin),m,pt,dphi,w) ;
1261       FillHistogram(Form("hMassPtTPCAllcore_cen%d",fCenBin),mV2,ptV2,dphi,w) ;
1262
1263       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1264         FillHistogram(Form("hMassPtCPV_cen%d",fCenBin),m ,pt,w) ;
1265         FillHistogram(Form("hMassPtCPVcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1266
1267         FillHistogram(Form("hMassPtTPCCPV_cen%d",fCenBin),m,pt,dphi,w) ;
1268         FillHistogram(Form("hMassPtTPCCPVcore_cen%d",fCenBin),mV2,ptV2,dphi,w) ;
1269       }
1270       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1271         FillHistogram(Form("hMassPtCPV2_cen%d",fCenBin),m ,pt,w) ;
1272
1273         FillHistogram(Form("hMassPtTPCCPV2_cen%d",fCenBin),m,pt,dphi,w) ;
1274         FillHistogram(Form("hMassPtTPCCPV2core_cen%d",fCenBin),mV2,ptV2,dphi,w) ;
1275       }
1276       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1277         FillHistogram(Form("hMassPtDisp2_cen%d",fCenBin),m ,pt,w) ;
1278         FillHistogram(Form("hMassPtDisp2core_cen%d",fCenBin),mV2 ,ptV2,w) ;
1279
1280         FillHistogram(Form("hMassPtTPCDisp2_cen%d",fCenBin),m,pt,dphi,w) ;
1281         FillHistogram(Form("hMassPtTPCDisp2core_cen%d",fCenBin),mV2,ptV2,dphi,w) ;
1282         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1283           FillHistogram(Form("hMassPtBoth2_cen%d",fCenBin),m ,pt,w) ;
1284           FillHistogram(Form("hMassPtBoth2core_cen%d",fCenBin),mV2 ,ptV2,w) ;
1285
1286           FillHistogram(Form("hMassPtTPCBoth2_cen%d",fCenBin),m,pt,dphi,w) ;
1287           FillHistogram(Form("hMassPtTPCBoth2core_cen%d",fCenBin),mV2,ptV2,dphi,w) ;
1288         }
1289      }
1290       if(ph1->IsDispOK() && ph2->IsDispOK()){
1291         FillHistogram(Form("hMassPtDisp_cen%d",fCenBin),m ,pt,w) ;
1292         FillHistogram(Form("hMassPtDispcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1293         if(ph1->IsPhoton()&&ph2->IsPhoton()){
1294           FillHistogram(Form("hMassPtDispwou_cen%d",fCenBin),m,pt,w) ;
1295         }
1296
1297         FillHistogram(Form("hMassPtTPCDisp_cen%d",fCenBin),m,pt,dphi,w) ;
1298         FillHistogram(Form("hMassPtTPCDispcore_cen%d",fCenBin),mV2,ptV2,dphi,w) ;
1299
1300         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1301           FillHistogram(Form("hMassPtBoth_cen%d",fCenBin),m ,pt,w) ;
1302           FillHistogram(Form("hMassPtBothcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1303
1304           FillHistogram(Form("hMassPtTPCBoth_cen%d",fCenBin),m,pt,dphi,w) ;
1305           FillHistogram(Form("hMassPtTPCBothcore_cen%d",fCenBin),mV2,ptV2,dphi,w) ;
1306         }
1307       }
1308     } // end of loop i2
1309   } // end of loop i1
1310
1311
1312   //now mixed, does not really matter old or new list of clusters
1313   for (Int_t i1=0; i1<inPHOSemb; i1++) {
1314     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
1315     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
1316       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
1317       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
1318         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
1319         TLorentzVector p12  = *ph1  + *ph2;
1320         TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1321
1322         Double_t dphi=p12.Phi()-fRPfull ;
1323         while(dphi<0)dphi+=TMath::Pi() ;
1324         while(dphi>TMath::Pi())dphi-=TMath::Pi() ;
1325         Double_t m=p12.M() ;
1326         Double_t mV2=pv12.M() ;
1327         Double_t pt = p12.Pt() ;
1328         Double_t ptV2 = pv12.Pt() ;
1329
1330         FillHistogram(Form("hMiMassPtAll_cen%d",fCenBin),m ,pt,1.) ;
1331         FillHistogram(Form("hMiMassPtAllcore_cen%d",fCenBin),mV2 ,ptV2,1.) ;
1332         if(ph1->IsPhoton()&&ph2->IsPhoton()){
1333           FillHistogram(Form("hMiMassPtAllwou_cen%d",fCenBin),m,pt,1.) ;
1334         }
1335
1336         FillHistogram(Form("hMiMassPtTPCAll_cen%d",fCenBin),m,pt,dphi) ;
1337         FillHistogram(Form("hMiMassPtTPCAllcore_cen%d",fCenBin),mV2,ptV2,dphi) ;
1338
1339         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1340           FillHistogram(Form("hMiMassPtCPV_cen%d",fCenBin),m ,pt,1.) ;
1341           FillHistogram(Form("hMiMassPtCPVcore_cen%d",fCenBin),mV2 ,ptV2,1.) ;
1342
1343           FillHistogram(Form("hMiMassPtTPCCPV_cen%d",fCenBin),m,pt,dphi) ;
1344           FillHistogram(Form("hMiMassPtTPCCPVcore_cen%d",fCenBin),mV2,ptV2,dphi) ;
1345
1346         }
1347         if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1348           FillHistogram(Form("hMiMassPtCPV2_cen%d",fCenBin),m ,pt,1.) ;
1349
1350           FillHistogram(Form("hMiMassPtTPCCPV2_cen%d",fCenBin),m,pt,dphi) ;
1351           FillHistogram(Form("hMiMassPtTPCCPV2core_cen%d",fCenBin),mV2,ptV2,dphi) ;
1352         }
1353         if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1354           FillHistogram(Form("hMiMassPtDisp2_cen%d",fCenBin),m ,pt,1.) ;
1355           FillHistogram(Form("hMiMassPtDisp2core_cen%d",fCenBin),mV2 ,ptV2,1.) ;
1356
1357           FillHistogram(Form("hMiMassPtTPCDisp2_cen%d",fCenBin),m,pt,dphi) ;
1358           FillHistogram(Form("hMiMassPtTPCDisp2core_cen%d",fCenBin),mV2,ptV2,dphi) ;
1359
1360           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1361             FillHistogram(Form("hMiMassPtBoth2_cen%d",fCenBin),m ,pt,1.) ;
1362             FillHistogram(Form("hMiMassPtBoth2core_cen%d",fCenBin),mV2 ,ptV2,1.) ;
1363
1364             FillHistogram(Form("hMiMassPtTPCBoth2_cen%d",fCenBin),m,pt,dphi) ;
1365             FillHistogram(Form("hMiMassPtTPCBoth2core_cen%d",fCenBin),mV2,ptV2,dphi) ;
1366           }
1367         }
1368         if(ph1->IsDispOK() && ph2->IsDispOK()){
1369           FillHistogram(Form("hMiMassPtDisp_cen%d",fCenBin),m ,pt,1.) ;
1370           FillHistogram(Form("hMiMassPtDispcore_cen%d",fCenBin),mV2 ,ptV2,1.) ;
1371           if(ph1->IsPhoton()&&ph2->IsPhoton()){
1372             FillHistogram(Form("hMiMassPtDispwou_cen%d",fCenBin),m,pt,1.) ;
1373           }
1374
1375           FillHistogram(Form("hMiMassPtTPCDisp_cen%d",fCenBin),m,pt,dphi) ;
1376           FillHistogram(Form("hMiMassPtTPCDispcore_cen%d",fCenBin),mV2,ptV2,dphi) ;
1377
1378           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1379             FillHistogram(Form("hMiMassPtBoth_cen%d",fCenBin),m ,pt,1.) ;
1380             FillHistogram(Form("hMiMassPtBothcore_cen%d",fCenBin),mV2 ,ptV2,1.) ;
1381
1382             FillHistogram(Form("hMiMassPtTPCBoth_cen%d",fCenBin),m,pt,dphi) ;
1383             FillHistogram(Form("hMiMassPtTPCBothcore_cen%d",fCenBin),mV2,ptV2,dphi) ;
1384           }
1385         }
1386       } // end of loop i2
1387     }
1388   } // end of loop i1
1389
1390
1391
1392   //Now we either add current events to stack or remove
1393   //If no photons in current event - no need to add it to mixed
1394   if(fPHOSEvent2->GetEntriesFast()>0){
1395     prevPHOS->AddFirst(fPHOSEvent2) ;
1396     fPHOSEvent2=0;
1397     delete fPHOSEvent1;
1398     fPHOSEvent1=0;
1399     if(prevPHOS->GetSize()>nMixEvents[fCenBin]){//Remove redundant events
1400       TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
1401       prevPHOS->RemoveLast() ;
1402       delete tmp ;
1403     }
1404   }
1405   // Post output data.
1406   PostData(1, fOutputContainer);
1407   fEventCounter++;
1408 }
1409 //___________________________________________________________________________
1410 Bool_t AliAnalysisTaskPi0DiffEfficiency::IsSameCluster(AliAODCaloCluster * c1, AliAODCaloCluster * c2)const{
1411  //Compare clusters before and after embedding
1412  //clusters are the same if
1413  // - Energy changed less than 0.1%  (numerical accuracy in reconstruction)
1414  // - lists of digits are the same
1415
1416  if(c1->GetNCells() != c2->GetNCells())
1417    return kFALSE ;
1418
1419  if(TMath::Abs(c1->E()-c2->E())>0.01*c1->E())
1420    return kFALSE ;
1421
1422  UShort_t *list1 = c1->GetCellsAbsId() ;
1423  UShort_t *list2 = c2->GetCellsAbsId() ;
1424  for(Int_t i=0; i<c1->GetNCells(); i++){
1425   if(list1[i] != list2[i])
1426     return kFALSE ;
1427  }
1428  return kTRUE ;
1429
1430 }
1431 //____________________________________________________________________________
1432 void  AliAnalysisTaskPi0DiffEfficiency::EvalLambdas(AliAODCaloCluster * clu, Int_t iR,Double_t &m02, Double_t &m20){
1433   //calculate dispecrsion of the cluster in the circle with radius distanceCut around the maximum
1434
1435   Double_t rCut=0. ;
1436   if(iR==0)
1437     rCut=3.5 ;
1438   else
1439     rCut=4.5 ;
1440
1441
1442   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1443 // Calculates the center of gravity in the local PHOS-module coordinates
1444   Float_t wtot = 0;
1445   Double_t xc[100]={0} ;
1446   Double_t zc[100]={0} ;
1447   Double_t x = 0 ;
1448   Double_t z = 0 ;
1449   Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
1450   const Double_t logWeight=4.5 ;
1451   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1452     Int_t relid[4] ;
1453     Float_t xi ;
1454     Float_t zi ;
1455     fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
1456     fPHOSGeo->RelPosInModule(relid, xi, zi);
1457     xc[iDigit]=xi ;
1458     zc[iDigit]=zi ;
1459     if (clu->E()>0 && elist[iDigit]>0) {
1460       Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1461       x    += xc[iDigit] * w ;
1462       z    += zc[iDigit] * w ;
1463       wtot += w ;
1464     }
1465   }
1466   if (wtot>0) {
1467     x /= wtot ;
1468     z /= wtot ;
1469   }
1470
1471   wtot = 0. ;
1472   Double_t dxx  = 0.;
1473   Double_t dzz  = 0.;
1474   Double_t dxz  = 0.;
1475   Double_t xCut = 0. ;
1476   Double_t zCut = 0. ;
1477   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1478     if (clu->E()>0 && elist[iDigit]>0.) {
1479         Double_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1480         Double_t xi= xc[iDigit] ;
1481         Double_t zi= zc[iDigit] ;
1482         if((xi-x)*(xi-x)+(zi-z)*(zi-z) < rCut*rCut){
1483           xCut += w * xi ;
1484           zCut += w * zi ;
1485           dxx  += w * xi * xi ;
1486           dzz  += w * zi * zi ;
1487           dxz  += w * xi * zi ;
1488           wtot += w ;
1489         }
1490     }
1491
1492   }
1493   if (wtot>0) {
1494     xCut/= wtot ;
1495     zCut/= wtot ;
1496     dxx /= wtot ;
1497     dzz /= wtot ;
1498     dxz /= wtot ;
1499     dxx -= xCut * xCut ;
1500     dzz -= zCut * zCut ;
1501     dxz -= xCut * zCut ;
1502
1503     m02 =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
1504     m20 =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
1505   }
1506   else {
1507     m20=m02=0.;
1508   }
1509
1510 }
1511 //___________________________________________________________________________
1512 void AliAnalysisTaskPi0DiffEfficiency::ProcessMC(){
1513   //fill histograms for efficiensy etc. calculation
1514   const Double_t rcut = 1. ; //cut for primary particles
1515   //---------First pi0/eta-----------------------------
1516   char partName[10] ;
1517
1518   AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
1519   if(!event) return ;
1520   TClonesArray *mcArray = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
1521   for(Int_t i=0;i<mcArray->GetEntriesFast();i++){
1522      AliAODMCParticle* particle =  (AliAODMCParticle*) mcArray->At(i);
1523     if(particle->GetPdgCode() == 111)
1524       snprintf(partName,10,"pi0") ;
1525     else
1526       if(particle->GetPdgCode() == 221)
1527         snprintf(partName,10,"eta") ;
1528       else
1529         if(particle->GetPdgCode() == 22)
1530            snprintf(partName,10,"gamma") ;
1531         else
1532            continue ;
1533
1534     //Primary particle
1535     Double_t r=TMath::Sqrt(particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv());
1536     if(r >rcut)
1537       continue ;
1538
1539     Double_t pt = particle->Pt() ;
1540     Double_t w = PrimaryWeight(pt) ;
1541     //Total number of pi0 with creation radius <1 cm
1542     FillHistogram(Form("hMC_all_%s_cen%d",partName,fCenBin),pt,w) ;
1543     if(TMath::Abs(particle->Y())<0.12){
1544       FillHistogram(Form("hMC_unitEta_%s_cen%d",partName,fCenBin),pt,w) ;
1545     }
1546
1547     FillHistogram(Form("hMC_rap_%s_cen%d",partName,fCenBin),particle->Y(),w) ;
1548
1549     Double_t phi=particle->Phi() ;
1550     while(phi<0.)phi+=TMath::TwoPi() ;
1551     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
1552     FillHistogram(Form("hMC_phi_%s_cen%d",partName,fCenBin),phi,w) ;
1553
1554   }
1555
1556   //Now calculate "Real" distribution of clusters with primary
1557   TClonesArray * clusters = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
1558   AliAODCaloCells * cellsEmb = (AliAODCaloCells *)event->FindListObject("EmbeddedPHOScells") ;
1559   Int_t multClust = clusters->GetEntriesFast();
1560   TClonesArray cluPrim("AliCaloPhoton",multClust) ; //clusters with primary
1561   Int_t inPHOS=0 ;
1562   Double_t vtx0[3] = {0,0,0};
1563   for (Int_t i=0; i<multClust; i++) {
1564     AliAODCaloCluster *clu = (AliAODCaloCluster*)clusters->At(i);
1565     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
1566     if(clu->GetLabel()<0) continue ;
1567
1568     Float_t  position[3];
1569     clu->GetPosition(position);
1570     TVector3 global(position) ;
1571     Int_t relId[4] ;
1572     fPHOSGeo->GlobalPos2RelId(global,relId) ;
1573     Int_t mod  = relId[0] ;
1574     Int_t cellX = relId[2];
1575     Int_t cellZ = relId[3] ;
1576     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
1577       continue ;
1578     if(clu->GetNCells()<3)
1579       continue ;
1580
1581     TLorentzVector pv1 ;
1582     clu->GetMomentum(pv1 ,vtx0);
1583
1584     pv1*=Scale(pv1.E()) ;
1585
1586
1587
1588     if(inPHOS>=cluPrim.GetSize()){
1589       cluPrim.Expand(inPHOS+50) ;
1590     }
1591     AliCaloPhoton * ph = new(cluPrim[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
1592     //AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
1593     ph->SetModule(mod) ;
1594
1595     AliPHOSAodCluster cluPHOS1(*clu);
1596     cluPHOS1.Recalibrate(fPHOSCalibData,cellsEmb); // modify the cell energies
1597     Double_t ecore=CoreEnergy(&cluPHOS1) ;
1598     ecore*=Scale(ecore) ;
1599     pv1*= ecore/pv1.E() ;
1600     ph->SetMomV2(&pv1) ;
1601     ph->SetNCells(clu->GetNCells());
1602     Double_t m02=0.,m20=0.;
1603     EvalLambdas(&cluPHOS1,0,m02, m20);
1604     ph->SetDispBit(TestLambda(clu->E(),m20,m02)) ;
1605     EvalLambdas(&cluPHOS1,1,m02, m20);
1606     ph->SetDisp2Bit(TestLambda2(clu->E(),m20,m02)) ;
1607
1608     ph->SetCPVBit(clu->GetEmcCpvDistance()>2.) ; //radius in sigmas
1609     ph->SetCPV2Bit(clu->GetEmcCpvDistance()>4.) ;
1610     ph->SetPhoton(clu->GetNExMax()<2); // Remember, if it is unfolded
1611     Double_t w=1. ;
1612     Int_t iprim = clu->GetLabel() ;
1613     if(iprim<mcArray->GetEntriesFast() && iprim>-1){
1614       AliAODMCParticle* particle =  (AliAODMCParticle*) mcArray->At(iprim);
1615       iprim=particle->GetMother() ;
1616       while(iprim>-1){
1617         particle =  (AliAODMCParticle*) mcArray->At(iprim);
1618         iprim=particle->GetMother() ;
1619       }
1620       if(particle->GetPdgCode()==111){
1621         Double_t pt = particle->Pt() ;
1622         w=PrimaryWeight(pt) ;
1623       }
1624     }
1625     ph->SetWeight(w) ;
1626
1627     inPHOS++ ;
1628   }
1629
1630   //Single photon
1631   char key[55] ;
1632   for (Int_t i1=0; i1<inPHOS; i1++) {
1633     AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
1634     Double_t pt=ph1->Pt() ;
1635     Double_t ptV2=ph1->GetMomV2()->Pt() ;
1636     Double_t w = ph1->GetWeight() ;
1637     FillHistogram(Form("hMCPhotAll_cen%d",fCenBin),pt,w) ;
1638     FillHistogram(Form("hMCPhotAllcore_cen%d",fCenBin),ptV2,w) ;
1639     if(ph1->IsPhoton()){
1640       FillHistogram(Form("hMCPhotAllwou_cen%d",fCenBin),pt,w) ;
1641     }
1642     if(ph1->IsCPVOK() ){
1643       FillHistogram(Form("hMCPhotCPV_cen%d",fCenBin),pt,w) ;
1644       FillHistogram(Form("hMCPhotCPVcore_cen%d",fCenBin),ptV2,w) ;
1645     }
1646     if(ph1->IsCPV2OK() ){
1647       snprintf(key,55,"hMCPhotCPV2_cen%d",fCenBin) ;
1648       FillHistogram(Form("hMCPhotCPV2_cen%d",fCenBin),pt,w) ;
1649     }
1650     if(ph1->IsDisp2OK()){
1651       FillHistogram(Form("hMCPhotDisp2_cen%d",fCenBin),pt,w) ;
1652       FillHistogram(Form("hMCPhotDisp2core_cen%d",fCenBin),ptV2,w) ;
1653       if(ph1->IsCPVOK()){
1654         FillHistogram(Form("hMCPhotBoth2_cen%d",fCenBin),pt,w) ;
1655         FillHistogram(Form("hMCPhotBoth2core_cen%d",fCenBin),ptV2,w) ;
1656       }
1657     }
1658     if(ph1->IsDispOK()){
1659       FillHistogram(Form("hMCPhotDisp_cen%d",fCenBin),pt,w) ;
1660       FillHistogram(Form("hMCPhotDispcore_cen%d",fCenBin),ptV2,w) ;
1661       if(ph1->IsPhoton()){
1662         FillHistogram(Form("hMCPhotDispwou_cen%d",fCenBin),pt,w) ;
1663       }
1664       if(ph1->IsCPVOK()){
1665         FillHistogram(Form("hMCPhotBoth_cen%d",fCenBin),pt,w) ;
1666         FillHistogram(Form("hMCPhotBothcore_cen%d",fCenBin),ptV2,w) ;
1667       }
1668     } // end of loop i2
1669   } // end of loop i1
1670
1671   // Fill Real disribution
1672   for (Int_t i1=0; i1<inPHOS-1; i1++) {
1673     AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
1674     Double_t w1 = ph1->GetWeight() ;
1675     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
1676       AliCaloPhoton * ph2=(AliCaloPhoton*)cluPrim.At(i2) ;
1677       Double_t w2 = ph2->GetWeight() ;
1678       Double_t w = TMath::Sqrt(w1*w2) ;
1679       TLorentzVector p12  = *ph1  + *ph2;
1680       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1681       Double_t m=p12.M() ;
1682       Double_t pt=p12.Pt() ;
1683       Double_t mV2=pv12.M() ;
1684       Double_t ptV2=pv12.Pt() ;
1685
1686       FillHistogram(Form("hMCMassPtAll_cen%d",fCenBin),m,pt,w) ;
1687       snprintf(key,55,"hMCMassPtAllcore_cen%d",fCenBin) ;
1688       FillHistogram(Form("hMCMassPtAllcore_cen%d",fCenBin),mV2,ptV2,w) ;
1689       if(ph1->IsPhoton()&&ph2->IsPhoton() ){
1690         FillHistogram(Form("hMCMassPtAllwou_cen%d",fCenBin),m,pt,w) ;
1691       }
1692
1693       if(ph1->Module()==1 && ph2->Module()==1)
1694         FillHistogram("hMCPi0M11",m,pt,w);
1695       else if(ph1->Module()==2 && ph2->Module()==2)
1696         FillHistogram("hMCPi0M22",m,pt,w);
1697       else if(ph1->Module()==3 && ph2->Module()==3)
1698         FillHistogram("hMCPi0M33",m,pt,w);
1699       else if(ph1->Module()==1 && ph2->Module()==2)
1700         FillHistogram("hMCPi0M12",m,pt,w);
1701       else if(ph1->Module()==1 && ph2->Module()==3)
1702         FillHistogram("hMCPi0M13",m,pt,w);
1703       else if(ph1->Module()==2 && ph2->Module()==3)
1704         FillHistogram("hMCPi0M23",m,pt,w);
1705
1706       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1707         FillHistogram(Form("hMCMassPtCPV_cen%d",fCenBin),m,pt,w) ;
1708         FillHistogram(Form("hMCMassPtCPVcore_cen%d",fCenBin),mV2,ptV2,w) ;
1709       }
1710       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1711         FillHistogram(Form("hMCMassPtCPV2core_cen%d",fCenBin),mV2,ptV2,w) ;
1712       }
1713       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1714         FillHistogram(Form("hMCMassPtDisp2_cen%d",fCenBin),m,pt,w) ;
1715         FillHistogram(Form("hMCMassPtDisp2core_cen%d",fCenBin),mV2,ptV2,w) ;
1716         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1717           FillHistogram(Form("hMCMassPtBoth2_cen%d",fCenBin),m,pt,w) ;
1718           FillHistogram(Form("hMCMassPtBoth2core_cen%d",fCenBin),mV2,ptV2,w) ;
1719         }
1720       }
1721       if(ph1->IsDispOK() && ph2->IsDispOK()){
1722         FillHistogram(Form("hMCMassPtDisp_cen%d",fCenBin),m,pt,w) ;
1723         FillHistogram(Form("hMCMassPtDispcore_cen%d",fCenBin),mV2,ptV2,w) ;
1724         if(ph1->IsPhoton()&& ph2->IsPhoton()){
1725           FillHistogram(Form("hMCMassPtDispwou_cen%d",fCenBin),m,pt,w) ;
1726         }
1727         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1728           FillHistogram(Form("hMCMassPtBoth_cen%d",fCenBin),m,pt,w) ;
1729           FillHistogram(Form("hMCMassPtBothcore_cen%d",fCenBin),mV2,ptV2,w) ;
1730         }
1731       }
1732     } // end of loop i2
1733   } // end of loop i1
1734 }
1735 //____________________________________________________________________________
1736 Double_t  AliAnalysisTaskPi0DiffEfficiency::CoreEnergy(AliPHOSAodCluster * clu){
1737   //calculate energy of the cluster in the circle with radius distanceCut around the maximum
1738
1739   //Can not use already calculated coordinates?
1740   //They have incidence correction...
1741   const Double_t distanceCut =3.5 ;
1742   const Double_t logWeight=4.5 ;
1743
1744   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1745 // Calculates the center of gravity in the local PHOS-module coordinates
1746   Float_t wtot = 0;
1747   Double_t xc[100]={0} ;
1748   Double_t zc[100]={0} ;
1749   Double_t x = 0 ;
1750   Double_t z = 0 ;
1751   Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
1752   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1753     Int_t relid[4] ;
1754     Float_t xi ;
1755     Float_t zi ;
1756     fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
1757     fPHOSGeo->RelPosInModule(relid, xi, zi);
1758     xc[iDigit]=xi ;
1759     zc[iDigit]=zi ;
1760     if (clu->E()>0 && elist[iDigit]>0) {
1761       Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1762       x    += xc[iDigit] * w ;
1763       z    += zc[iDigit] * w ;
1764       wtot += w ;
1765     }
1766   }
1767   if (wtot>0) {
1768     x /= wtot ;
1769     z /= wtot ;
1770   }
1771   Double_t coreE=0. ;
1772   for(Int_t iDigit=0; iDigit < mulDigit; iDigit++) {
1773     Double_t distance = TMath::Sqrt((xc[iDigit]-x)*(xc[iDigit]-x)+(zc[iDigit]-z)*(zc[iDigit]-z)) ;
1774     if(distance < distanceCut)
1775       coreE += elist[iDigit] ;
1776   }
1777   //Apply non-linearity correction
1778   return (0.0241+1.0504*coreE+0.000249*coreE*coreE) ;
1779 }
1780 //_____________________________________________________________________________
1781 Bool_t AliAnalysisTaskPi0DiffEfficiency::TestLambda(Double_t pt,Double_t l1,Double_t l2){
1782
1783   //For R=3.5
1784    Double_t  l1Mean  = 1.170014 -0.059465/(1.+0.019343*pt+0.147544*pt*pt) ;
1785    Double_t  l2Mean = 1.626270 + 0.761554*exp(-1.213839*pt)-0.020027*pt ;
1786    Double_t  l1Sigma = 0.133409 + 0.261307*exp(-0.636874*pt)-0.002849*pt ;
1787    Double_t  l2Sigma = 0.289698 + 0.459400*exp(-1.214242*pt)-0.012578*pt ;
1788    Double_t  c=-0.124103 ;
1789 /*
1790   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1791   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1792   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1793   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1794   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1795 */
1796   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1797               0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1798               0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1799   return (R2<2.5*2.5) ;
1800
1801 }
1802 //_____________________________________________________________________________
1803 Bool_t AliAnalysisTaskPi0DiffEfficiency::TestLambda2(Double_t pt,Double_t l1,Double_t l2){
1804
1805 //For R=4.5
1806   Double_t   l1Mean  = 1.150200 + 0.097886/(1.+1.486645*pt+0.000038*pt*pt) ;
1807   Double_t   l2Mean = 1.574706 + 0.997966*exp(-0.895075*pt)-0.010666*pt ;
1808   Double_t   l1Sigma = 0.100255 + 0.337177*exp(-0.517684*pt)+0.001170*pt ;
1809   Double_t   l2Sigma = 0.232580 + 0.573401*exp(-0.735903*pt)-0.002325*pt ;
1810   Double_t   c = -0.110983 -0.017353/(1.-1.836995*pt+0.934517*pt*pt) ;
1811
1812 /*
1813   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1814   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1815   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1816   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1817   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1818 */
1819   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1820               0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1821               0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1822   return (R2<2.5*2.5) ;
1823
1824 }
1825 //____________________________________________________________________________
1826 Double_t AliAnalysisTaskPi0DiffEfficiency::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge){
1827   //Parameterization of LHC10h period
1828   //_true if neutral_
1829
1830   Double_t meanX=0;
1831   Double_t meanZ=0.;
1832   Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
1833               6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219);
1834   Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ;
1835   AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
1836   Double_t mf = event->GetMagneticField(); //Positive for ++ and negative for --
1837
1838   if(mf<0.){ //field --
1839     meanZ = -0.468318 ;
1840     if(charge>0)
1841       meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
1842     else
1843       meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
1844   }
1845   else{ //Field ++
1846     meanZ= -0.468318;
1847     if(charge>0)
1848       meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
1849     else
1850       meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
1851   }
1852
1853   Double_t rz=(dz-meanZ)/sz ;
1854   Double_t rx=(dx-meanX)/sx ;
1855   return TMath::Sqrt(rx*rx+rz*rz) ;
1856 }
1857 //____________________________________________________________________________
1858 Bool_t AliAnalysisTaskPi0DiffEfficiency::TestTOF(Double_t t, Double_t e){
1859
1860   Double_t sigma = TMath::Sqrt(2.23183e-09*2.23183e-09
1861                              +2.24611e-09*2.24611e-09/e
1862                              +5.65054e-09*5.65054e-09/e/e) ;
1863   sigma=TMath::Min(20.e-9,sigma) ; //for the soft (<400 MeV) part
1864   Double_t mean=1.1e-9 ;
1865   if(TMath::Abs(t-mean)<2.*sigma)
1866     return kTRUE ;
1867   else
1868     if(TMath::Abs(t-mean+100.e-9)<2.*sigma)
1869       return kTRUE ;
1870
1871   return kFALSE ;
1872 }
1873 //_____________________________________________________________________________
1874 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x)const{
1875   //FillHistogram
1876   TObject * tmp = fOutputContainer->FindObject(key) ;
1877   if(!tmp){
1878     AliInfo(Form("can not find histogram <%s> ",key)) ;
1879     return ;
1880   }
1881   if(tmp->IsA() == TClass::GetClass("TH1I")){
1882     ((TH1I*)tmp)->Fill(x) ;
1883     return ;
1884   }
1885   if(tmp->IsA() == TClass::GetClass("TH1F")){
1886     ((TH1F*)tmp)->Fill(x) ;
1887     return ;
1888   }
1889   if(tmp->IsA() == TClass::GetClass("TH1D")){
1890     ((TH1D*)tmp)->Fill(x) ;
1891     return ;
1892   }
1893   AliInfo(Form("can not find 1D histogram <%s> ",key)) ;
1894 }
1895 //_____________________________________________________________________________
1896 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x,Double_t y)const{
1897   //FillHistogram
1898   TObject * tmp = fOutputContainer->FindObject(key) ;
1899   if(!tmp){
1900     AliInfo(Form("can not find histogram <%s> ",key)) ;
1901     return ;
1902   }
1903   if(tmp->IsA() == TClass::GetClass("TH1F")){
1904     ((TH1F*)tmp)->Fill(x,y) ;
1905     return ;
1906   }
1907   if(tmp->IsA() == TClass::GetClass("TH2F")){
1908     ((TH2F*)tmp)->Fill(x,y) ;
1909     return ;
1910   }
1911   AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
1912 }
1913
1914 //_____________________________________________________________________________
1915 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
1916   //Fills 1D histograms with key
1917   TObject * tmp = fOutputContainer->FindObject(key) ;
1918   if(!tmp){
1919     AliInfo(Form("can not find histogram <%s> ",key)) ;
1920     return ;
1921   }
1922   if(tmp->IsA() == TClass::GetClass("TH2F")){
1923     ((TH2F*)tmp)->Fill(x,y,z) ;
1924     return ;
1925   }
1926   if(tmp->IsA() == TClass::GetClass("TH3F")){
1927     ((TH3F*)tmp)->Fill(x,y,z) ;
1928     return ;
1929   }
1930 }
1931 //_____________________________________________________________________________
1932 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z, Double_t w) const{
1933   //Fills 1D histograms with key
1934   TObject * tmp = fOutputContainer->FindObject(key) ;
1935   if(!tmp){
1936     AliInfo(Form("can not find histogram <%s> ",key)) ;
1937     return ;
1938   }
1939   if(tmp->IsA() == TClass::GetClass("TH3F")){
1940     ((TH3F*)tmp)->Fill(x,y,z,w) ;
1941     return ;
1942   }
1943 }
1944 //_____________________________________________________________________________
1945 Bool_t AliAnalysisTaskPi0DiffEfficiency::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
1946 {
1947   //Check if this channel belogs to the good ones
1948
1949   if(strcmp(det,"PHOS")==0){
1950     if(mod>5 || mod<1){
1951       AliError(Form("No bad map for PHOS module %d ",mod)) ;
1952       return kTRUE ;
1953     }
1954     if(!fPHOSBadMap[mod]){
1955       AliError(Form("No Bad map for PHOS module %d",mod)) ;
1956       return kTRUE ;
1957     }
1958     if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
1959       return kFALSE ;
1960     else
1961       return kTRUE ;
1962   }
1963   else{
1964     AliError(Form("Can not find bad channels for detector %s ",det)) ;
1965   }
1966   return kTRUE ;
1967 }
1968 //_____________________________________________________________________________
1969 Double_t AliAnalysisTaskPi0DiffEfficiency::PrimaryWeight(Double_t x){
1970
1971    Double_t w=1 ;
1972
1973
1974    //Parameterization of LHC10h data from Jan 2013 (pi0 spectrum)
1975    //Should be consistend with spectrum parameterization used in simulation
1976    if(fCenBin==0) //0-5
1977      w = (0.561741+0.332841*x-0.007082*x*x)/(1.-0.447804*x+0.157830*x*x)+0.080394*x ;
1978    if(fCenBin==1) //5-10
1979      w = (0.659096+0.101701*x+0.042395*x*x)/(1.-0.470110*x+0.154665*x*x)+0.052932*x ;
1980    if(fCenBin==2) //10-20
1981      w = (0.615575+0.005621*x+0.069263*x*x)/(1.-0.485422*x+0.160822*x*x)+0.040865*x ;
1982    if(fCenBin==3) //20-40
1983      w = (0.441240+0.158358*x+0.059458*x*x)/(1.-0.332609*x+0.147528*x*x)+0.037926*x ;
1984    if(fCenBin==4) //40-60
1985      w = (0.467895-0.001113*x+0.029610*x*x)/(1.-0.266502*x+0.065105*x*x)+0.025431*x ;
1986    if(fCenBin==5) //60-80
1987      w = (0.465204-0.139362*x+0.043500*x*x)/(1.-0.371689*x+0.067606*x*x)+0.006519*x ;
1988
1989 /*
1990   //Parameterization of photon spectrum 25.02
1991   if(fCenBin==0) //0-5
1992      w=(0.870487-0.494032*x+0.076334*x*x+0.001065*x*x*x)/(1.-0.646014*x+0.113839*x*x);
1993   if(fCenBin==1) //5-10
1994      w=(-8.310403+15.767226*x-2.167756*x*x+0.184356*x*x*x)/(1.+4.556793*x+0.980941*x*x);
1995   if(fCenBin==2) //10-20
1996      w=(-5.281594+7.477165*x-0.688609*x*x+0.097601*x*x*x)/(1.+1.102693*x+0.882454*x*x);
1997   if(fCenBin==3) //20-40
1998      w=(0.789472-4.750155*x+4.381545*x*x-0.029239*x*x*x)/(1.-3.738304*x+3.328318*x*x);
1999   if(fCenBin==4) //40-60
2000      w=(0.876792-0.491379*x+0.130872*x*x-0.001390*x*x*x)/(1.+-0.511934*x+0.112705*x*x);
2001   if(fCenBin==5) //60-80
2002      w=(0.719912-0.167292*x+0.017196*x*x+0.000861*x*x*x)/(1.-0.336146*x+0.037731*x*x);
2003 */
2004   return TMath::Max(0.,w) ;
2005 }
2006
2007
2008