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