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 = event->GetHeader() ;
 718
 719   // Checks if we have a primary vertex
 720   // Get primary vertices form ESD
 721   const AliAODVertex *esdVertex5 = event->GetPrimaryVertex();
 722
 723  // don't rely on ESD vertex, assume (0,0,0)
 724   Double_t vtx0[3] ={0.,0.,0.};
 725
 726
 727   FillHistogram("hZvertex",esdVertex5->GetZ());
 728   if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
 729     PostData(1, fOutputContainer);
 730     return;
 731   }
 732   FillHistogram("hSelEvents",2.5) ;
 733
 734   //Vtx class z-bin
 735   //  Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
 736   //  if(zvtx<0)zvtx=0 ;
 737   //  if(zvtx>9)zvtx=9 ;
 738   Int_t zvtx=0 ;
 739
 740 //  fCentrality=header->GetCentralityP()->GetCentralityPercentile("V0M"); // returns the centrality percentile,
 741 //                                                          //a float from 0 to 100 (or to the trigger efficiency)
 742    fCentrality=header->GetZDCN2Energy() ;
 743
 744   if( fCentrality < 0. || fCentrality>80.){
 745     PostData(1, fOutputContainer);
 746     return;
 747   }
 748   FillHistogram("hSelEvents",3.5) ;
 749   Float_t bins[7]={0.,5.,10.,20.,40.,60.,80.} ;
 750   fCenBin=0 ;
 751   while(fCenBin<6 && fCentrality > bins[fCenBin+1])
 752     fCenBin++ ;
 753
 754   const Int_t nMixEvents[6]={4,4,5,10,20,20} ;
 755
 756
 757   //reaction plane
 758   fRPfull= header->GetZDCN1Energy() ;
 759   if(fRPfull==999){ //reaction plain was not defined
 760     PostData(1, fOutputContainer);
 761     return;
 762   }
 763
 764   FillHistogram("hSelEvents",4.5) ;
 765   //All event selections done
 766   FillHistogram("hCentrality",fCentrality) ;
 767   //Reaction plain is defined in the range (-pi/2;pi/2)
 768   //We have 10 bins
 769   Int_t irp=Int_t(10.*fRPfull/TMath::Pi());
 770   if(irp>9)irp=9 ;
 771
 772   if(!fPHOSEvents[zvtx][fCenBin][irp])
 773     fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
 774   TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
 775
 776   // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
 777   if(fEventCounter == 0) {
 778     for(Int_t mod=0; mod<5; mod++) {
 779       const TGeoHMatrix* m =header->GetPHOSMatrix(mod) ;
 780       fPHOSGeo->SetMisalMatrix(m,mod) ;
 781       Printf("PHOS geo matrix for module # %d is set: %p\n", mod,m);
 782     }
 783     fEventCounter++ ;
 784   }
 785
 786   TClonesArray *mcArray = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
 787
 788   ProcessMC() ;
 789
 790   if(fPHOSEvent1){
 791     fPHOSEvent1->Clear() ;
 792     fPHOSEvent2->Clear() ;
 793   }
 794   else{
 795     fPHOSEvent1 = new TClonesArray("AliCaloPhoton",200) ;
 796     fPHOSEvent2 = new TClonesArray("AliCaloPhoton",200) ;
 797   }
 798
 799   TClonesArray * clustersEmb = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
 800   AliAODCaloCells * cellsEmb = (AliAODCaloCells *)event->FindListObject("EmbeddedPHOScells") ;
 801   TClonesArray * clustersOld = event->GetCaloClusters() ;
 802   AliAODCaloCells * cellsOld = event->GetPHOSCells() ;
 803 //  TClonesArray *mcArray = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
 804
 805
 806   TVector3 vertex(vtx0);
 807   char key[55] ;
 808   //Before Embedding
 809   Int_t multClustOld = clustersOld->GetEntriesFast();
 810   Int_t multClustEmb = clustersEmb->GetEntriesFast();
 811   Int_t inPHOSold=0 ;
 812   Int_t inPHOSemb=0 ;
 813   for (Int_t i=0; i<multClustOld; i++) {
 814     AliAODCaloCluster *clu = (AliAODCaloCluster*)clustersOld->At(i);
 815     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
 816
 817     Bool_t survive=kFALSE ;
 818     for(Int_t ii=0;(ii<multClustEmb)&&(!survive);ii++){
 819        AliAODCaloCluster *clu2 = (AliAODCaloCluster*)clustersEmb->At(ii);
 820        survive=IsSameCluster(clu,clu2);
 821     }
 822
 823
 824     Float_t  position[3];
 825     clu->GetPosition(position);
 826     TVector3 global(position) ;
 827     Int_t relId[4] ;
 828     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 829     Int_t mod  = relId[0] ;
 830     Int_t cellX = relId[2];
 831     Int_t cellZ = relId[3] ;
 832     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
 833       continue ;
 834     if(clu->GetNCells()<3)
 835       continue ;
 836
 837     snprintf(key,55,"hCluM%d",mod) ;
 838     FillHistogram(key,cellX,cellZ,1.);
 839
 840     TLorentzVector pv1 ;
 841     clu->GetMomentum(pv1 ,vtx0);
 842
 843     pv1*=Scale(pv1.E()) ;
 844
 845     if(inPHOSold>=fPHOSEvent1->GetSize()){
 846       fPHOSEvent1->Expand(inPHOSold+50) ;
 847     }
 848     AliCaloPhoton * ph = new((*fPHOSEvent1)[inPHOSold]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
 849     ph->SetModule(mod) ;
 850     AliPHOSAodCluster cluPHOS1(*clu);
 851     cluPHOS1.Recalibrate(fPHOSCalibData,cellsOld); // modify the cell energies
 852     Double_t ecore=CoreEnergy(&cluPHOS1) ;
 853     ecore*=Scale(ecore) ;
 854     pv1*= ecore/pv1.E() ;
 855     ph->SetMomV2(&pv1) ;
 856     ph->SetNCells(clu->GetNCells());
 857     Double_t m02=0.,m20=0.;
 858     EvalLambdas(&cluPHOS1,0,m02, m20);
 859     ph->SetDispBit(TestLambda(clu->E(),m20,m02)) ;
 860     EvalLambdas(&cluPHOS1,1,m02, m20);
 861     ph->SetDisp2Bit(TestLambda2(clu->E(),m20,m02)) ;
 862
 863     ph->SetCPVBit(clu->GetEmcCpvDistance()>2.) ;
 864     ph->SetCPV2Bit(clu->GetEmcCpvDistance()>4.) ;
 865     if(!survive) //this cluster found in list after embedding, skipping it
 866       ph->SetTagged(1) ;
 867     ph->SetPhoton(clu->GetNExMax()<2); // Remember, if it is unfolded
 868     ph->SetWeight(1.) ; //All weights for real particles ==1.
 869
 870     if(!survive){
 871       Double_t distBC=clu->GetDistanceToBadChannel();
 872       if(distBC>2.)
 873         FillHistogram(Form("hPhotAll_DistBad2_cen%d",fCenBin),ph->Pt(),-1.) ;
 874         if(distBC>4.)
 875           FillHistogram(Form("hPhotAll_DistBad4_cen%d",fCenBin),ph->Pt(),-1.) ;
 876           if(distBC>6.)
 877             FillHistogram(Form("hPhotAll_DistBad6_cen%d",fCenBin),ph->Pt(),-1.) ;
 878     }
 879     inPHOSold++ ;
 880   }
 881
 882   for (Int_t i=0; i<multClustEmb; i++) {
 883     AliAODCaloCluster *clu = (AliAODCaloCluster*)clustersEmb->At(i);
 884     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
 885
 886     Bool_t survive=kFALSE ;
 887     for(Int_t ii=0;(ii<multClustOld)&&(!survive);ii++){
 888        AliAODCaloCluster *clu2 = (AliAODCaloCluster*)clustersOld->At(ii);
 889        survive=IsSameCluster(clu,clu2);
 890     }
 891
 892     Float_t  position[3];
 893     clu->GetPosition(position);
 894     TVector3 global(position) ;
 895     Int_t relId[4] ;
 896     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 897     Int_t mod  = relId[0] ;
 898     Int_t cellX = relId[2];
 899     Int_t cellZ = relId[3] ;
 900     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
 901       continue ;
 902     if(clu->GetNCells()<3)
 903       continue ;
 904
 905     snprintf(key,55,"hCluM%d",mod) ;
 906     FillHistogram(key,cellX,cellZ,1.);
 907
 908     TLorentzVector pv1 ;
 909     clu->GetMomentum(pv1 ,vtx0);
 910
 911     pv1*=Scale(pv1.E()) ;
 912
 913     if(inPHOSemb>=fPHOSEvent2->GetSize()){
 914       fPHOSEvent2->Expand(inPHOSemb+50) ;
 915     }
 916     AliCaloPhoton * ph = new((*fPHOSEvent2)[inPHOSemb]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
 917     ph->SetModule(mod) ;
 918     AliPHOSAodCluster cluPHOS1(*clu);
 919     cluPHOS1.Recalibrate(fPHOSCalibData,cellsEmb); // modify the cell energies
 920     Double_t ecore=CoreEnergy(&cluPHOS1) ;
 921     ecore*=Scale(ecore) ;
 922     pv1*= ecore/pv1.E() ;
 923     ph->SetMomV2(&pv1) ;
 924     ph->SetNCells(clu->GetNCells());
 925     Double_t m02=0.,m20=0.;
 926     EvalLambdas(&cluPHOS1,0,m02, m20);
 927     ph->SetDispBit(TestLambda(clu->E(),m20,m02)) ;
 928     EvalLambdas(&cluPHOS1,1,m02, m20);
 929     ph->SetDisp2Bit(TestLambda2(clu->E(),m20,m02)) ;
 930
 931     ph->SetCPVBit(clu->GetEmcCpvDistance()>2.) ;
 932     ph->SetCPV2Bit(clu->GetEmcCpvDistance()>4.) ;
 933     if(!survive) //this cluster found in list after embedding, skipping it
 934       ph->SetTagged(1) ;
 935     ph->SetPhoton(clu->GetNExMax()<2); // Remember, if it is unfolded
 936
 937     //Set weight for embedded particles
 938     Double_t w=1. ;
 939     Int_t iprim = clu->GetLabel() ;
 940     if(iprim<mcArray->GetEntriesFast() && iprim>-1){
 941       AliAODMCParticle* particle =  (AliAODMCParticle*) mcArray->At(iprim);
 942       iprim=particle->GetMother() ;
 943       while(iprim>-1){
 944         particle =  (AliAODMCParticle*) mcArray->At(iprim);
 945         iprim=particle->GetMother() ;
 946       }
 947       if(particle->GetPdgCode()==111){
 948         Double_t pt = particle->Pt() ;
 949         w=PrimaryWeight(pt) ;
 950       }
 951     }
 952     ph->SetWeight(w) ;
 953
 954
 955     if(!survive){
 956       Double_t distBC=clu->GetDistanceToBadChannel();
 957       if(distBC>2.)
 958         FillHistogram(Form("hPhotAll_DistBad2_cen%d",fCenBin),ph->Pt()) ;
 959         if(distBC>4.)
 960           FillHistogram(Form("hPhotAll_DistBad4_cen%d",fCenBin),ph->Pt()) ;
 961           if(distBC>6.)
 962             FillHistogram(Form("hPhotAll_DistBad6_cen%d",fCenBin),ph->Pt()) ;
 963     }
 964
 965     inPHOSemb++ ;
 966   }
 967
 968   //Single photon
 969   for (Int_t i1=0; i1<inPHOSold; i1++) {
 970     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent1->At(i1) ;
 971     if(!ph1->IsTagged())
 972       continue ;
 973
 974     Double_t dphi=ph1->Phi()-fRPfull ;
 975     while(dphi<0)dphi+=TMath::Pi() ;
 976     while(dphi>TMath::Pi())dphi-=TMath::Pi() ;
 977     Double_t pt = ph1->Pt() ;
 978     Double_t ptV2=ph1->GetMomV2()->Pt() ;
 979 //always 1!    Double_t w=ph1->GetWeight() ;
 980
 981     FillHistogram(Form("hPhotPhiTPCAll_cen%d",fCenBin),pt,dphi,-1.) ;
 982     FillHistogram(Form("hPhotPhiTPCAllcore_cen%d",fCenBin),ptV2,dphi,-1.) ;
 983
 984     FillHistogram(Form("hPhotAll_cen%d",fCenBin),pt,-1.) ;
 985     FillHistogram(Form("hNegPhotAll_cen%d",fCenBin),pt,-1.) ;
 986     FillHistogram(Form("hPhotAllcore_cen%d",fCenBin),ptV2,-1.) ;
 987     FillHistogram(Form("hNegPhotAllcore_cen%d",fCenBin),ptV2,-1.) ;
 988     if(ph1->IsPhoton()){
 989       FillHistogram(Form("hPhotAllwou_cen%d",fCenBin),pt,-1.) ;
 990     }
 991     if(ph1->IsCPVOK() ){
 992       FillHistogram(Form("hPhotCPV_cen%d",fCenBin),pt,-1.) ;
 993       FillHistogram(Form("hNegPhotCPV_cen%d",fCenBin),pt,-1.) ;
 994       FillHistogram(Form("hPhotCPVcore_cen%d",fCenBin),ptV2,-1.) ;
 995       FillHistogram(Form("hNegPhotCPVcore_cen%d",fCenBin),ptV2,-1.) ;
 996       FillHistogram(Form("hPhotPhiTPCCPV_cen%d",fCenBin),pt,dphi,-1.) ;
 997       FillHistogram(Form("hPhotPhiTPCCPVcore_cen%d",fCenBin),ptV2,dphi,-1.) ;
 998     }
 999     if(ph1->IsCPV2OK() ){
1000       FillHistogram(Form("hPhotCPV2_cen%d",fCenBin),pt,-1.) ;
1001       FillHistogram(Form("hPhotCPV2core_cen%d",fCenBin),ptV2,-1.) ;
1002       FillHistogram(Form("hNegPhotCPV2_cen%d",fCenBin),pt,-1.) ;
1003       FillHistogram(Form("hPhotPhiTPCCPV2_cen%d",fCenBin),pt,dphi,-1.) ;
1004       FillHistogram(Form("hPhotPhiTPCCPV2core_cen%d",fCenBin),ptV2,dphi,-1.) ;
1005     }
1006     if(ph1->IsDisp2OK()){
1007       FillHistogram(Form("hPhotDisp2_cen%d",fCenBin),pt,-1.) ;
1008       FillHistogram(Form("hPhotDisp2core_cen%d",fCenBin),ptV2,-1.) ;
1009       FillHistogram(Form("hNegPhotDisp2_cen%d",fCenBin),pt,-1.) ;
1010       FillHistogram(Form("hNegPhotDisp2core_cen%d",fCenBin),ptV2,-1.) ;
1011
1012       FillHistogram(Form("hPhotPhiTPCDisp2_cen%d",fCenBin),pt,dphi,-1.) ;
1013       FillHistogram(Form("hPhotPhiTPCDisp2core_cen%d",fCenBin),ptV2,dphi,-1.) ;
1014       if(ph1->IsCPVOK()){
1015         FillHistogram(Form("hPhotBoth2_cen%d",fCenBin),pt,-1.) ;
1016         FillHistogram(Form("hNegPhotBoth2_cen%d",fCenBin),pt,-1.) ;
1017         FillHistogram(Form("hPhotBoth2core_cen%d",fCenBin),ptV2,-1.) ;
1018         FillHistogram(Form("hNegPhotBoth2core_cen%d",fCenBin),ptV2,-1.) ;
1019         FillHistogram(Form("hPhotPhiTPCBoth2_cen%d",fCenBin),pt,dphi,-1.) ;
1020         FillHistogram(Form("hPhotPhiTPCBoth2core_cen%d",fCenBin),ptV2,dphi,-1.) ;
1021       }
1022     }
1023     if(ph1->IsDispOK()){
1024       FillHistogram(Form("hPhotDisp_cen%d",fCenBin),pt,-1.) ;
1025       FillHistogram(Form("hNegPhotDisp_cen%d",fCenBin),pt,-1.) ;
1026       FillHistogram(Form("hPhotDispcore_cen%d",fCenBin),ptV2,-1.) ;
1027       FillHistogram(Form("hNegPhotDispcore_cen%d",fCenBin),ptV2,-1.) ;
1028       if(ph1->IsPhoton()){
1029         FillHistogram(Form("hPhotDispwou_cen%d",fCenBin),pt,-1.) ;
1030       }
1031
1032       FillHistogram(Form("hPhotPhiTPCDisp_cen%d",fCenBin),pt,dphi,-1.) ;
1033       FillHistogram(Form("hPhotPhiTPCDispcore_cen%d",fCenBin),ptV2,dphi,-1.) ;
1034       if(ph1->IsCPVOK()){
1035         FillHistogram(Form("hPhotBoth_cen%d",fCenBin),pt,-1.) ;
1036         FillHistogram(Form("hNegPhotBoth_cen%d",fCenBin),pt,-1.) ;
1037         FillHistogram(Form("hPhotBothcore_cen%d",fCenBin),ptV2,-1.) ;
1038         FillHistogram(Form("hNegPhotBothcore_cen%d",fCenBin),ptV2,-1.) ;
1039
1040         FillHistogram(Form("hPhotPhiTPCBoth_cen%d",fCenBin),pt,dphi,-1.) ;
1041         FillHistogram(Form("hPhotPhiTPCBothcore_cen%d",fCenBin),ptV2,dphi,-1.) ;
1042       }
1043     }
1044   } // end of loop i1
1045
1046   for (Int_t i1=0; i1<inPHOSemb; i1++) {
1047     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
1048     if(!ph1->IsTagged())
1049       continue ;
1050
1051     Double_t dphi=ph1->Phi()-fRPfull ;
1052     while(dphi<0)dphi+=TMath::Pi() ;
1053     while(dphi>TMath::Pi())dphi-=TMath::Pi() ;
1054     Double_t pt = ph1->Pt() ;
1055     Double_t ptV2=ph1->GetMomV2()->Pt() ;
1056     Double_t w=ph1->GetWeight() ;
1057
1058     FillHistogram(Form("hPhotPhiTPCAll_cen%d",fCenBin),pt,dphi,w) ;
1059     FillHistogram(Form("hPhotPhiTPCAllcore_cen%d",fCenBin),ptV2,dphi,w) ;
1060
1061     FillHistogram(Form("hPhotAll_cen%d",fCenBin),pt,w) ;
1062     FillHistogram(Form("hPhotAllcore_cen%d",fCenBin),ptV2,w) ;
1063     if(ph1->IsPhoton()){
1064       FillHistogram(Form("hPhotAllwou_cen%d",fCenBin),pt,w) ;
1065     }
1066     if(ph1->IsCPVOK() ){
1067       FillHistogram(Form("hPhotCPV_cen%d",fCenBin),pt,w) ;
1068       FillHistogram(Form("hPhotCPVcore_cen%d",fCenBin),ptV2,w) ;
1069
1070       FillHistogram(Form("hPhotPhiTPCCPV_cen%d",fCenBin),pt,dphi,w) ;
1071       FillHistogram(Form("hPhotPhiTPCCPVcore_cen%d",fCenBin),ptV2,dphi,w) ;
1072     }
1073     if(ph1->IsCPV2OK() ){
1074       FillHistogram(Form("hPhotCPV2_cen%d",fCenBin),pt,w) ;
1075       FillHistogram(Form("hPhotCPV2core_cen%d",fCenBin),ptV2,w) ;
1076       FillHistogram(Form("hPhotPhiTPCCPV2_cen%d",fCenBin),pt,dphi,w) ;
1077       FillHistogram(Form("hPhotPhiTPCCPV2core_cen%d",fCenBin),ptV2,dphi,w) ;
1078     }
1079     if(ph1->IsDisp2OK()){
1080       FillHistogram(Form("hPhotDisp2_cen%d",fCenBin),pt,w) ;
1081       FillHistogram(Form("hPhotDisp2core_cen%d",fCenBin),ptV2,w) ;
1082
1083       FillHistogram(Form("hPhotPhiTPCDisp2_cen%d",fCenBin),pt,dphi,w) ;
1084       FillHistogram(Form("hPhotPhiTPCDisp2core_cen%d",fCenBin),ptV2,dphi,w) ;
1085       if(ph1->IsCPVOK()){
1086         FillHistogram(Form("hPhotBoth2_cen%d",fCenBin),pt,w) ;
1087         FillHistogram(Form("hPhotBoth2core_cen%d",fCenBin),ptV2,w) ;
1088
1089         FillHistogram(Form("hPhotPhiTPCBoth2_cen%d",fCenBin),pt,dphi,w) ;
1090         FillHistogram(Form("hPhotPhiTPCBoth2core_cen%d",fCenBin),ptV2,dphi,w) ;
1091       }
1092     }
1093     if(ph1->IsDispOK()){
1094       FillHistogram(Form("hPhotDisp_cen%d",fCenBin),pt,w) ;
1095       FillHistogram(Form("hPhotDispcore_cen%d",fCenBin),ptV2,w) ;
1096       if(ph1->IsPhoton()){
1097         FillHistogram(Form("hPhotDispwou_cen%d",fCenBin),pt,w) ;
1098       }
1099
1100       FillHistogram(Form("hPhotPhiTPCBoth_cen%d",fCenBin),pt,dphi,w) ;
1101       FillHistogram(Form("hPhotPhiTPCBothcore_cen%d",fCenBin),ptV2,dphi,w) ;
1102       if(ph1->IsCPVOK()){
1103         FillHistogram(Form("hPhotBoth_cen%d",fCenBin),pt,w) ;
1104         FillHistogram(Form("hPhotBothcore_cen%d",fCenBin),ptV2,w) ;
1105
1106         FillHistogram(Form("hPhotPhiTPCBoth_cen%d",fCenBin),pt,dphi,w) ;
1107         FillHistogram(Form("hPhotPhiTPCBothcore_cen%d",fCenBin),ptV2,dphi,w) ;
1108       }
1109     } // end of loop i2
1110   } // end of loop i1
1111
1112
1113   const Double_t prob[10]={0.1,0.2,0.3,1.,1.,1.,1.,1.,1.,1.} ; //Probabilities to accept Tagged+Bg pair
1114
1115
1116   // Fill Real disribution:
1117   // Disappeared clusters enter with negative contribution
1118   // In addition fill control histogam with Real before embedding
1119   for (Int_t i1=0; i1<inPHOSold-1; i1++) {
1120     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent1->At(i1) ;
1121     Double_t w1 = ph1->GetWeight() ;
1122     for (Int_t i2=i1+1; i2<inPHOSold; i2++) {
1123       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent1->At(i2) ;
1124
1125       TLorentzVector p12  = *ph1  + *ph2;
1126       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1127       Double_t w2 = ph2->GetWeight() ;
1128       Double_t w = TMath::Sqrt(w1*w2) ;
1129
1130       Double_t dphi=p12.Phi()-fRPfull ;
1131       while(dphi<0)dphi+=TMath::Pi() ;
1132       while(dphi>TMath::Pi())dphi-=TMath::Pi() ;
1133       Double_t m=p12.M() ;
1134       Double_t mV2=pv12.M() ;
1135       Double_t pt = p12.Pt() ;
1136       Double_t ptV2 = pv12.Pt() ;
1137
1138
1139       //Fill Controll histogram: Real before embedding
1140       FillHistogram(Form("hOldMassPtAll_cen%d",fCenBin),m,pt,-w) ;
1141       FillHistogram(Form("hOldMassPtAllcore_cen%d",fCenBin),mV2,ptV2,-w) ;
1142       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1143         FillHistogram(Form("hOldMassPtCPV_cen%d",fCenBin),m ,pt,-w) ;
1144         FillHistogram(Form("hOldMassPtCPV_cen%d",fCenBin),mV2, ptV2,-w) ;
1145       }
1146       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1147         FillHistogram(Form("hOldMassPtCPV2_cen%d",fCenBin),m ,pt,-w) ;
1148       }
1149       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1150         FillHistogram(Form("hOldMassPtDisp2_cen%d",fCenBin),m ,pt,-w) ;
1151         FillHistogram(Form("hOldMassPtDisp2core_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1152         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1153           FillHistogram(Form("hOldMassPtBoth2_cen%d",fCenBin),m ,pt,-w) ;
1154           FillHistogram(Form("hOldMassPtBoth2core_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1155         }
1156       }
1157       if(ph1->IsDispOK() && ph2->IsDispOK()){
1158         FillHistogram(Form("hOldMassPtDisp_cen%d",fCenBin),m ,pt,-w) ;
1159         FillHistogram(Form("hOldMassPtDispcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1160         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1161           FillHistogram(Form("hOldMassPtBoth_cen%d",fCenBin),m ,pt,-w) ;
1162           FillHistogram(Form("hOldMassPtBothcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1163         }
1164       }
1165
1166       //Now fill main histograms with negative contributions
1167       if(!(ph1->IsTagged() || ph2->IsTagged()) )
1168         continue ;
1169       if(!ph1->IsTagged() || !ph2->IsTagged()){ //Tagged + Bg combination
1170         if(gRandom->Uniform()>prob[fCenBin])
1171           continue ;
1172       }
1173       FillHistogram(Form("hMassPtAll_cen%d",fCenBin),m ,pt,-w) ;
1174       FillHistogram(Form("hMassPtAllcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1175       if(ph1->IsPhoton()&&ph2->IsPhoton()){
1176         FillHistogram(Form("hMassPtAllwou_cen%d",fCenBin),m,pt,-w) ;
1177       }
1178
1179       FillHistogram(Form("hMassPtTPCAll_cen%d",fCenBin),m,pt,dphi,-w) ;
1180       FillHistogram(Form("hMassPtTPCAllcore_cen%d",fCenBin),mV2,ptV2,dphi,-w) ;
1181
1182       FillHistogram(Form("hNegMassPtAll_cen%d",fCenBin),m ,pt,-w) ;
1183       FillHistogram(Form("hNegMassPtAllcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1184
1185       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1186         FillHistogram(Form("hMassPtCPV_cen%d",fCenBin),m ,pt,-w) ;
1187         FillHistogram(Form("hMassPtCPVcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1188
1189         FillHistogram(Form("hMassPtTPCCPV_cen%d",fCenBin),m,pt,dphi,-w) ;
1190         FillHistogram(Form("hMassPtTPCCPVcore_cen%d",fCenBin),mV2,ptV2,dphi,-w) ;
1191
1192         FillHistogram(Form("hNegMassPtCPV_cen%d",fCenBin),m ,pt,-w) ;
1193         FillHistogram(Form("hNegMassPtCPVcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1194       }
1195       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1196         FillHistogram(Form("hMassPtCPV2_cen%d",fCenBin),m ,pt,-w) ;
1197         FillHistogram(Form("hNegMassPtCPV2_cen%d",fCenBin),m ,pt,-w) ;
1198
1199         FillHistogram(Form("hMassPtTPCCPV2_cen%d",fCenBin),m,pt,dphi,-w) ;
1200         FillHistogram(Form("hMassPtTPCCPV2core_cen%d",fCenBin),mV2,ptV2,dphi,-w) ;
1201       }
1202
1203       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1204         FillHistogram(Form("hMassPtDisp2_cen%d",fCenBin),m ,pt,-w) ;
1205         FillHistogram(Form("hMassPtDisp2core_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1206
1207         FillHistogram(Form("hMassPtTPCDisp2_cen%d",fCenBin),m,pt,dphi,-w) ;
1208         FillHistogram(Form("hMassPtTPCDisp2core_cen%d",fCenBin),mV2,ptV2,dphi,-w) ;
1209
1210         FillHistogram(Form("hNegMassPtDisp2_cen%d",fCenBin),m ,pt,-w) ;
1211         FillHistogram(Form("hNegMassPtDisp2core_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1212         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1213           FillHistogram(Form("hMassPtBoth2_cen%d",fCenBin),m ,pt,-w) ;
1214           FillHistogram(Form("hMassPtBoth2core_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1215           FillHistogram(Form("hNegMassPtBoth2_cen%d",fCenBin),m ,pt,-w) ;
1216           FillHistogram(Form("hNegMassPtBoth2core_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1217
1218           FillHistogram(Form("hMassPtTPCBoth2_cen%d",fCenBin),m,pt,dphi,-w) ;
1219           FillHistogram(Form("hMassPtTPCBoth2core_cen%d",fCenBin),mV2,ptV2,dphi,-w) ;
1220         }
1221       }
1222       if(ph1->IsDispOK() && ph2->IsDispOK()){
1223         FillHistogram(Form("hMassPtDisp_cen%d",fCenBin),m ,pt,-w) ;
1224         FillHistogram(Form("hMassPtDispcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1225         if(ph1->IsPhoton()&&ph2->IsPhoton()){
1226           FillHistogram(Form("hMassPtDispwou_cen%d",fCenBin),m,pt,-w) ;
1227         }
1228         FillHistogram(Form("hNegMassPtDisp_cen%d",fCenBin),m ,pt,-w) ;
1229         FillHistogram(Form("hNegMassPtDispcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1230
1231         FillHistogram(Form("hMassPtTPCDisp_cen%d",fCenBin),m,pt,dphi,-w) ;
1232         FillHistogram(Form("hMassPtTPCDispcore_cen%d",fCenBin),mV2,ptV2,dphi,-w) ;
1233         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1234           FillHistogram(Form("hMassPtBoth_cen%d",fCenBin),m ,pt,-w) ;
1235           FillHistogram(Form("hMassPtBothcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1236           FillHistogram(Form("hNegMassPtBoth_cen%d",fCenBin),m ,pt,-w) ;
1237           FillHistogram(Form("hNegMassPtBothcore_cen%d",fCenBin),mV2 ,ptV2,-w) ;
1238
1239           FillHistogram(Form("hMassPtTPCBoth_cen%d",fCenBin),m,pt,dphi,-w) ;
1240           FillHistogram(Form("hMassPtTPCBothcore_cen%d",fCenBin),mV2,ptV2,dphi,-w) ;
1241         }
1242       }
1243     } // end of loop i2
1244   } // end of loop i1
1245
1246
1247   // Further fill Real disribution
1248   // now with positive contribution from new clusters
1249   // ass well fill controll histogram
1250   for (Int_t i1=0; i1<inPHOSemb-1; i1++) {
1251     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
1252     Double_t w1 = ph1->GetWeight() ;
1253     for (Int_t i2=i1+1; i2<inPHOSemb; i2++) {
1254       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent2->At(i2) ;
1255
1256       TLorentzVector p12  = *ph1  + *ph2;
1257       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1258       Double_t m=p12.M() ;
1259       Double_t mV2=pv12.M() ;
1260       Double_t pt = p12.Pt() ;
1261       Double_t ptV2 = pv12.Pt() ;
1262       Double_t w2 = ph2->GetWeight() ;
1263       Double_t w = TMath::Sqrt(w1*w2) ;
1264
1265       // Controll histogram: Real after embedding
1266       FillHistogram(Form("hNewMassPtAll_cen%d",fCenBin),m ,pt,w) ;
1267       FillHistogram(Form("hNewMassPtAllcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1268       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1269         FillHistogram(Form("hNewMassPtCPV_cen%d",fCenBin),m ,pt,w) ;
1270         FillHistogram(Form("hNewMassPtCPVcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1271       }
1272       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1273         FillHistogram(Form("hNewMassPtCPV2_cen%d",fCenBin),m ,pt,w) ;
1274       }
1275       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1276         FillHistogram(Form("hNewMassPtDisp2_cen%d",fCenBin),m ,pt,w) ;
1277         FillHistogram(Form("hNewMassPtDisp2core_cen%d",fCenBin),mV2 ,ptV2,w) ;
1278         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1279           FillHistogram(Form("hNewMassPtBoth2_cen%d",fCenBin),m ,pt,w) ;
1280           FillHistogram(Form("hNewMassPtBoth2core_cen%d",fCenBin),mV2,ptV2,w) ;
1281         }
1282       }
1283       if(ph1->IsDispOK() && ph2->IsDispOK()){
1284         FillHistogram(Form("hNewMassPtDisp_cen%d",fCenBin),m ,pt,w) ;
1285         FillHistogram(Form("hNewMassPtDispcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1286         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1287           FillHistogram(Form("hNewMassPtBoth_cen%d",fCenBin),m ,pt,w) ;
1288           FillHistogram(Form("hNewMassPtBothcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1289         }
1290       }
1291
1292       //Now fill main histogamm
1293       //new clusters with positive contribution
1294       if(!(ph1->IsTagged() || ph2->IsTagged()) )
1295         continue ;
1296       if(!ph1->IsTagged() || !ph2->IsTagged()){ //Tagged + Bg combination
1297         if(gRandom->Uniform()>prob[fCenBin])
1298           continue ;
1299       }
1300
1301       Double_t dphi=p12.Phi()-fRPfull ;
1302       while(dphi<0)dphi+=TMath::Pi() ;
1303       while(dphi>TMath::Pi())dphi-=TMath::Pi() ;
1304
1305       FillHistogram(Form("hMassPtAll_cen%d",fCenBin),m ,pt,w) ;
1306       FillHistogram(Form("hMassPtAllcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1307       if(ph1->IsPhoton()&&ph2->IsPhoton()){
1308         FillHistogram(Form("hMassPtAllwou_cen%d",fCenBin),m,pt,w) ;
1309       }
1310
1311       FillHistogram(Form("hMassPtTPCAll_cen%d",fCenBin),m,pt,dphi,w) ;
1312       FillHistogram(Form("hMassPtTPCAllcore_cen%d",fCenBin),mV2,ptV2,dphi,w) ;
1313
1314       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1315         FillHistogram(Form("hMassPtCPV_cen%d",fCenBin),m ,pt,w) ;
1316         FillHistogram(Form("hMassPtCPVcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1317
1318         FillHistogram(Form("hMassPtTPCCPV_cen%d",fCenBin),m,pt,dphi,w) ;
1319         FillHistogram(Form("hMassPtTPCCPVcore_cen%d",fCenBin),mV2,ptV2,dphi,w) ;
1320       }
1321       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1322         FillHistogram(Form("hMassPtCPV2_cen%d",fCenBin),m ,pt,w) ;
1323
1324         FillHistogram(Form("hMassPtTPCCPV2_cen%d",fCenBin),m,pt,dphi,w) ;
1325         FillHistogram(Form("hMassPtTPCCPV2core_cen%d",fCenBin),mV2,ptV2,dphi,w) ;
1326       }
1327       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1328         FillHistogram(Form("hMassPtDisp2_cen%d",fCenBin),m ,pt,w) ;
1329         FillHistogram(Form("hMassPtDisp2core_cen%d",fCenBin),mV2 ,ptV2,w) ;
1330
1331         FillHistogram(Form("hMassPtTPCDisp2_cen%d",fCenBin),m,pt,dphi,w) ;
1332         FillHistogram(Form("hMassPtTPCDisp2core_cen%d",fCenBin),mV2,ptV2,dphi,w) ;
1333         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1334           FillHistogram(Form("hMassPtBoth2_cen%d",fCenBin),m ,pt,w) ;
1335           FillHistogram(Form("hMassPtBoth2core_cen%d",fCenBin),mV2 ,ptV2,w) ;
1336
1337           FillHistogram(Form("hMassPtTPCBoth2_cen%d",fCenBin),m,pt,dphi,w) ;
1338           FillHistogram(Form("hMassPtTPCBoth2core_cen%d",fCenBin),mV2,ptV2,dphi,w) ;
1339         }
1340      }
1341       if(ph1->IsDispOK() && ph2->IsDispOK()){
1342         FillHistogram(Form("hMassPtDisp_cen%d",fCenBin),m ,pt,w) ;
1343         FillHistogram(Form("hMassPtDispcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1344         if(ph1->IsPhoton()&&ph2->IsPhoton()){
1345           FillHistogram(Form("hMassPtDispwou_cen%d",fCenBin),m,pt,w) ;
1346         }
1347
1348         FillHistogram(Form("hMassPtTPCDisp_cen%d",fCenBin),m,pt,dphi,w) ;
1349         FillHistogram(Form("hMassPtTPCDispcore_cen%d",fCenBin),mV2,ptV2,dphi,w) ;
1350
1351         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1352           FillHistogram(Form("hMassPtBoth_cen%d",fCenBin),m ,pt,w) ;
1353           FillHistogram(Form("hMassPtBothcore_cen%d",fCenBin),mV2 ,ptV2,w) ;
1354
1355           FillHistogram(Form("hMassPtTPCBoth_cen%d",fCenBin),m,pt,dphi,w) ;
1356           FillHistogram(Form("hMassPtTPCBothcore_cen%d",fCenBin),mV2,ptV2,dphi,w) ;
1357         }
1358       }
1359     } // end of loop i2
1360   } // end of loop i1
1361
1362
1363   //now mixed, does not really matter old or new list of clusters
1364   for (Int_t i1=0; i1<inPHOSemb; i1++) {
1365     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
1366     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
1367       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
1368       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
1369         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
1370         TLorentzVector p12  = *ph1  + *ph2;
1371         TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1372
1373         Double_t dphi=p12.Phi()-fRPfull ;
1374         while(dphi<0)dphi+=TMath::Pi() ;
1375         while(dphi>TMath::Pi())dphi-=TMath::Pi() ;
1376         Double_t m=p12.M() ;
1377         Double_t mV2=pv12.M() ;
1378         Double_t pt = p12.Pt() ;
1379         Double_t ptV2 = pv12.Pt() ;
1380
1381         FillHistogram(Form("hMiMassPtAll_cen%d",fCenBin),m ,pt,1.) ;
1382         FillHistogram(Form("hMiMassPtAllcore_cen%d",fCenBin),mV2 ,ptV2,1.) ;
1383         if(ph1->IsPhoton()&&ph2->IsPhoton()){
1384           FillHistogram(Form("hMiMassPtAllwou_cen%d",fCenBin),m,pt,1.) ;
1385         }
1386
1387         FillHistogram(Form("hMiMassPtTPCAll_cen%d",fCenBin),m,pt,dphi) ;
1388         FillHistogram(Form("hMiMassPtTPCAllcore_cen%d",fCenBin),mV2,ptV2,dphi) ;
1389
1390         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1391           FillHistogram(Form("hMiMassPtCPV_cen%d",fCenBin),m ,pt,1.) ;
1392           FillHistogram(Form("hMiMassPtCPVcore_cen%d",fCenBin),mV2 ,ptV2,1.) ;
1393
1394           FillHistogram(Form("hMiMassPtTPCCPV_cen%d",fCenBin),m,pt,dphi) ;
1395           FillHistogram(Form("hMiMassPtTPCCPVcore_cen%d",fCenBin),mV2,ptV2,dphi) ;
1396
1397         }
1398         if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1399           FillHistogram(Form("hMiMassPtCPV2_cen%d",fCenBin),m ,pt,1.) ;
1400
1401           FillHistogram(Form("hMiMassPtTPCCPV2_cen%d",fCenBin),m,pt,dphi) ;
1402           FillHistogram(Form("hMiMassPtTPCCPV2core_cen%d",fCenBin),mV2,ptV2,dphi) ;
1403         }
1404         if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1405           FillHistogram(Form("hMiMassPtDisp2_cen%d",fCenBin),m ,pt,1.) ;
1406           FillHistogram(Form("hMiMassPtDisp2core_cen%d",fCenBin),mV2 ,ptV2,1.) ;
1407
1408           FillHistogram(Form("hMiMassPtTPCDisp2_cen%d",fCenBin),m,pt,dphi) ;
1409           FillHistogram(Form("hMiMassPtTPCDisp2core_cen%d",fCenBin),mV2,ptV2,dphi) ;
1410
1411           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1412             FillHistogram(Form("hMiMassPtBoth2_cen%d",fCenBin),m ,pt,1.) ;
1413             FillHistogram(Form("hMiMassPtBoth2core_cen%d",fCenBin),mV2 ,ptV2,1.) ;
1414
1415             FillHistogram(Form("hMiMassPtTPCBoth2_cen%d",fCenBin),m,pt,dphi) ;
1416             FillHistogram(Form("hMiMassPtTPCBoth2core_cen%d",fCenBin),mV2,ptV2,dphi) ;
1417           }
1418         }
1419         if(ph1->IsDispOK() && ph2->IsDispOK()){
1420           FillHistogram(Form("hMiMassPtDisp_cen%d",fCenBin),m ,pt,1.) ;
1421           FillHistogram(Form("hMiMassPtDispcore_cen%d",fCenBin),mV2 ,ptV2,1.) ;
1422           if(ph1->IsPhoton()&&ph2->IsPhoton()){
1423             FillHistogram(Form("hMiMassPtDispwou_cen%d",fCenBin),m,pt,1.) ;
1424           }
1425
1426           FillHistogram(Form("hMiMassPtTPCDisp_cen%d",fCenBin),m,pt,dphi) ;
1427           FillHistogram(Form("hMiMassPtTPCDispcore_cen%d",fCenBin),mV2,ptV2,dphi) ;
1428
1429           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1430             FillHistogram(Form("hMiMassPtBoth_cen%d",fCenBin),m ,pt,1.) ;
1431             FillHistogram(Form("hMiMassPtBothcore_cen%d",fCenBin),mV2 ,ptV2,1.) ;
1432
1433             FillHistogram(Form("hMiMassPtTPCBoth_cen%d",fCenBin),m,pt,dphi) ;
1434             FillHistogram(Form("hMiMassPtTPCBothcore_cen%d",fCenBin),mV2,ptV2,dphi) ;
1435           }
1436         }
1437       } // end of loop i2
1438     }
1439   } // end of loop i1
1440
1441
1442
1443   //Now we either add current events to stack or remove
1444   //If no photons in current event - no need to add it to mixed
1445   if(fPHOSEvent2->GetEntriesFast()>0){
1446     prevPHOS->AddFirst(fPHOSEvent2) ;
1447     fPHOSEvent2=0;
1448     delete fPHOSEvent1;
1449     fPHOSEvent1=0;
1450     if(prevPHOS->GetSize()>nMixEvents[fCenBin]){//Remove redundant events
1451       TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
1452       prevPHOS->RemoveLast() ;
1453       delete tmp ;
1454     }
1455   }
1456   // Post output data.
1457   PostData(1, fOutputContainer);
1458   fEventCounter++;
1459 }
1460 //___________________________________________________________________________
1461 Bool_t AliAnalysisTaskPi0DiffEfficiency::IsSameCluster(AliAODCaloCluster * c1, AliAODCaloCluster * c2)const{
1462  //Compare clusters before and after embedding
1463  //clusters are the same if
1464  // - Energy changed less than 0.1%  (numerical accuracy in reconstruction)
1465  // - lists of digits are the same
1466
1467  if(c1->GetNCells() != c2->GetNCells())
1468    return kFALSE ;
1469
1470  if(TMath::Abs(c1->E()-c2->E())>0.01*c1->E())
1471    return kFALSE ;
1472
1473  UShort_t *list1 = c1->GetCellsAbsId() ;
1474  UShort_t *list2 = c2->GetCellsAbsId() ;
1475  for(Int_t i=0; i<c1->GetNCells(); i++){
1476   if(list1[i] != list2[i])
1477     return kFALSE ;
1478  }
1479  return kTRUE ;
1480
1481 }
1482 //____________________________________________________________________________
1483 void  AliAnalysisTaskPi0DiffEfficiency::EvalLambdas(AliAODCaloCluster * clu, Int_t iR,Double_t &m02, Double_t &m20){
1484   //calculate dispecrsion of the cluster in the circle with radius distanceCut around the maximum
1485
1486   Double_t rCut=0. ;
1487   if(iR==0)
1488     rCut=3.5 ;
1489   else
1490     rCut=4.5 ;
1491
1492
1493   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1494 // Calculates the center of gravity in the local PHOS-module coordinates
1495   Float_t wtot = 0;
1496   Double_t xc[100]={0} ;
1497   Double_t zc[100]={0} ;
1498   Double_t x = 0 ;
1499   Double_t z = 0 ;
1500   Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
1501   const Double_t logWeight=4.5 ;
1502   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1503     Int_t relid[4] ;
1504     Float_t xi ;
1505     Float_t zi ;
1506     fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
1507     fPHOSGeo->RelPosInModule(relid, xi, zi);
1508     xc[iDigit]=xi ;
1509     zc[iDigit]=zi ;
1510     if (clu->E()>0 && elist[iDigit]>0) {
1511       Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1512       x    += xc[iDigit] * w ;
1513       z    += zc[iDigit] * w ;
1514       wtot += w ;
1515     }
1516   }
1517   if (wtot>0) {
1518     x /= wtot ;
1519     z /= wtot ;
1520   }
1521
1522   wtot = 0. ;
1523   Double_t dxx  = 0.;
1524   Double_t dzz  = 0.;
1525   Double_t dxz  = 0.;
1526   Double_t xCut = 0. ;
1527   Double_t zCut = 0. ;
1528   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1529     if (clu->E()>0 && elist[iDigit]>0.) {
1530         Double_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1531         Double_t xi= xc[iDigit] ;
1532         Double_t zi= zc[iDigit] ;
1533         if((xi-x)*(xi-x)+(zi-z)*(zi-z) < rCut*rCut){
1534           xCut += w * xi ;
1535           zCut += w * zi ;
1536           dxx  += w * xi * xi ;
1537           dzz  += w * zi * zi ;
1538           dxz  += w * xi * zi ;
1539           wtot += w ;
1540         }
1541     }
1542
1543   }
1544   if (wtot>0) {
1545     xCut/= wtot ;
1546     zCut/= wtot ;
1547     dxx /= wtot ;
1548     dzz /= wtot ;
1549     dxz /= wtot ;
1550     dxx -= xCut * xCut ;
1551     dzz -= zCut * zCut ;
1552     dxz -= xCut * zCut ;
1553
1554     m02 =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
1555     m20 =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
1556   }
1557   else {
1558     m20=m02=0.;
1559   }
1560
1561 }
1562 //___________________________________________________________________________
1563 void AliAnalysisTaskPi0DiffEfficiency::ProcessMC(){
1564   //fill histograms for efficiensy etc. calculation
1565   const Double_t rcut = 1. ; //cut for primary particles
1566   //---------First pi0/eta-----------------------------
1567   char partName[10] ;
1568
1569   AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
1570   if(!event) return ;
1571   TClonesArray *mcArray = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
1572   for(Int_t i=0;i<mcArray->GetEntriesFast();i++){
1573      AliAODMCParticle* particle =  (AliAODMCParticle*) mcArray->At(i);
1574     if(particle->GetPdgCode() == 111)
1575       snprintf(partName,10,"pi0") ;
1576     else
1577       if(particle->GetPdgCode() == 221)
1578         snprintf(partName,10,"eta") ;
1579       else
1580         if(particle->GetPdgCode() == 22)
1581            snprintf(partName,10,"gamma") ;
1582         else
1583            continue ;
1584
1585     //Primary particle
1586     Double_t r=TMath::Sqrt(particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv());
1587     if(r >rcut)
1588       continue ;
1589
1590     Double_t pt = particle->Pt() ;
1591     Double_t w = PrimaryWeight(pt) ;
1592     //Total number of pi0 with creation radius <1 cm
1593     FillHistogram(Form("hMC_all_%s_cen%d",partName,fCenBin),pt,w) ;
1594     if(TMath::Abs(particle->Y())<0.12){
1595       FillHistogram(Form("hMC_unitEta_%s_cen%d",partName,fCenBin),pt,w) ;
1596     }
1597
1598     FillHistogram(Form("hMC_rap_%s_cen%d",partName,fCenBin),particle->Y(),w) ;
1599
1600     Double_t phi=particle->Phi() ;
1601     while(phi<0.)phi+=TMath::TwoPi() ;
1602     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
1603     FillHistogram(Form("hMC_phi_%s_cen%d",partName,fCenBin),phi,w) ;
1604
1605   }
1606
1607   //Now calculate "Real" distribution of clusters with primary
1608   TClonesArray * clusters = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
1609   AliAODCaloCells * cellsEmb = (AliAODCaloCells *)event->FindListObject("EmbeddedPHOScells") ;
1610   Int_t multClust = clusters->GetEntriesFast();
1611   TClonesArray cluPrim("AliCaloPhoton",multClust) ; //clusters with primary
1612   Int_t inPHOS=0 ;
1613   Double_t vtx0[3] = {0,0,0};
1614   for (Int_t i=0; i<multClust; i++) {
1615     AliAODCaloCluster *clu = (AliAODCaloCluster*)clusters->At(i);
1616     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
1617     if(clu->GetLabel()<0) continue ;
1618
1619     Float_t  position[3];
1620     clu->GetPosition(position);
1621     TVector3 global(position) ;
1622     Int_t relId[4] ;
1623     fPHOSGeo->GlobalPos2RelId(global,relId) ;
1624     Int_t mod  = relId[0] ;
1625     Int_t cellX = relId[2];
1626     Int_t cellZ = relId[3] ;
1627     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
1628       continue ;
1629     if(clu->GetNCells()<3)
1630       continue ;
1631
1632     TLorentzVector pv1 ;
1633     clu->GetMomentum(pv1 ,vtx0);
1634
1635     pv1*=Scale(pv1.E()) ;
1636
1637
1638
1639     if(inPHOS>=cluPrim.GetSize()){
1640       cluPrim.Expand(inPHOS+50) ;
1641     }
1642     AliCaloPhoton * ph = new(cluPrim[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
1643     //AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
1644     ph->SetModule(mod) ;
1645
1646     AliPHOSAodCluster cluPHOS1(*clu);
1647     cluPHOS1.Recalibrate(fPHOSCalibData,cellsEmb); // modify the cell energies
1648     Double_t ecore=CoreEnergy(&cluPHOS1) ;
1649     ecore*=Scale(ecore) ;
1650     pv1*= ecore/pv1.E() ;
1651     ph->SetMomV2(&pv1) ;
1652     ph->SetNCells(clu->GetNCells());
1653     Double_t m02=0.,m20=0.;
1654     EvalLambdas(&cluPHOS1,0,m02, m20);
1655     ph->SetDispBit(TestLambda(clu->E(),m20,m02)) ;
1656     EvalLambdas(&cluPHOS1,1,m02, m20);
1657     ph->SetDisp2Bit(TestLambda2(clu->E(),m20,m02)) ;
1658
1659     ph->SetCPVBit(clu->GetEmcCpvDistance()>2.) ; //radius in sigmas
1660     ph->SetCPV2Bit(clu->GetEmcCpvDistance()>4.) ;
1661     ph->SetPhoton(clu->GetNExMax()<2); // Remember, if it is unfolded
1662     Double_t w=1. ;
1663     Int_t iprim = clu->GetLabel() ;
1664     if(iprim<mcArray->GetEntriesFast() && iprim>-1){
1665       AliAODMCParticle* particle =  (AliAODMCParticle*) mcArray->At(iprim);
1666       iprim=particle->GetMother() ;
1667       while(iprim>-1){
1668         particle =  (AliAODMCParticle*) mcArray->At(iprim);
1669         iprim=particle->GetMother() ;
1670       }
1671       if(particle->GetPdgCode()==111){
1672         Double_t pt = particle->Pt() ;
1673         w=PrimaryWeight(pt) ;
1674       }
1675     }
1676     ph->SetWeight(w) ;
1677
1678     inPHOS++ ;
1679   }
1680
1681   //Single photon
1682   char key[55] ;
1683   for (Int_t i1=0; i1<inPHOS; i1++) {
1684     AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
1685     Double_t pt=ph1->Pt() ;
1686     Double_t ptV2=ph1->GetMomV2()->Pt() ;
1687     Double_t w = ph1->GetWeight() ;
1688     FillHistogram(Form("hMCPhotAll_cen%d",fCenBin),pt,w) ;
1689     FillHistogram(Form("hMCPhotAllcore_cen%d",fCenBin),ptV2,w) ;
1690     if(ph1->IsPhoton()){
1691       FillHistogram(Form("hMCPhotAllwou_cen%d",fCenBin),pt,w) ;
1692     }
1693     if(ph1->IsCPVOK() ){
1694       FillHistogram(Form("hMCPhotCPV_cen%d",fCenBin),pt,w) ;
1695       FillHistogram(Form("hMCPhotCPVcore_cen%d",fCenBin),ptV2,w) ;
1696     }
1697     if(ph1->IsCPV2OK() ){
1698       snprintf(key,55,"hMCPhotCPV2_cen%d",fCenBin) ;
1699       FillHistogram(Form("hMCPhotCPV2_cen%d",fCenBin),pt,w) ;
1700     }
1701     if(ph1->IsDisp2OK()){
1702       FillHistogram(Form("hMCPhotDisp2_cen%d",fCenBin),pt,w) ;
1703       FillHistogram(Form("hMCPhotDisp2core_cen%d",fCenBin),ptV2,w) ;
1704       if(ph1->IsCPVOK()){
1705         FillHistogram(Form("hMCPhotBoth2_cen%d",fCenBin),pt,w) ;
1706         FillHistogram(Form("hMCPhotBoth2core_cen%d",fCenBin),ptV2,w) ;
1707       }
1708     }
1709     if(ph1->IsDispOK()){
1710       FillHistogram(Form("hMCPhotDisp_cen%d",fCenBin),pt,w) ;
1711       FillHistogram(Form("hMCPhotDispcore_cen%d",fCenBin),ptV2,w) ;
1712       if(ph1->IsPhoton()){
1713         FillHistogram(Form("hMCPhotDispwou_cen%d",fCenBin),pt,w) ;
1714       }
1715       if(ph1->IsCPVOK()){
1716         FillHistogram(Form("hMCPhotBoth_cen%d",fCenBin),pt,w) ;
1717         FillHistogram(Form("hMCPhotBothcore_cen%d",fCenBin),ptV2,w) ;
1718       }
1719     } // end of loop i2
1720   } // end of loop i1
1721
1722   // Fill Real disribution
1723   for (Int_t i1=0; i1<inPHOS-1; i1++) {
1724     AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
1725     Double_t w1 = ph1->GetWeight() ;
1726     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
1727       AliCaloPhoton * ph2=(AliCaloPhoton*)cluPrim.At(i2) ;
1728       Double_t w2 = ph2->GetWeight() ;
1729       Double_t w = TMath::Sqrt(w1*w2) ;
1730       TLorentzVector p12  = *ph1  + *ph2;
1731       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1732       Double_t m=p12.M() ;
1733       Double_t pt=p12.Pt() ;
1734       Double_t mV2=pv12.M() ;
1735       Double_t ptV2=pv12.Pt() ;
1736
1737       FillHistogram(Form("hMCMassPtAll_cen%d",fCenBin),m,pt,w) ;
1738       snprintf(key,55,"hMCMassPtAllcore_cen%d",fCenBin) ;
1739       FillHistogram(Form("hMCMassPtAllcore_cen%d",fCenBin),mV2,ptV2,w) ;
1740       if(ph1->IsPhoton()&&ph2->IsPhoton() ){
1741         FillHistogram(Form("hMCMassPtAllwou_cen%d",fCenBin),m,pt,w) ;
1742       }
1743
1744       if(ph1->Module()==1 && ph2->Module()==1)
1745         FillHistogram("hMCPi0M11",m,pt,w);
1746       else if(ph1->Module()==2 && ph2->Module()==2)
1747         FillHistogram("hMCPi0M22",m,pt,w);
1748       else if(ph1->Module()==3 && ph2->Module()==3)
1749         FillHistogram("hMCPi0M33",m,pt,w);
1750       else if(ph1->Module()==1 && ph2->Module()==2)
1751         FillHistogram("hMCPi0M12",m,pt,w);
1752       else if(ph1->Module()==1 && ph2->Module()==3)
1753         FillHistogram("hMCPi0M13",m,pt,w);
1754       else if(ph1->Module()==2 && ph2->Module()==3)
1755         FillHistogram("hMCPi0M23",m,pt,w);
1756
1757       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1758         FillHistogram(Form("hMCMassPtCPV_cen%d",fCenBin),m,pt,w) ;
1759         FillHistogram(Form("hMCMassPtCPVcore_cen%d",fCenBin),mV2,ptV2,w) ;
1760       }
1761       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1762         FillHistogram(Form("hMCMassPtCPV2core_cen%d",fCenBin),mV2,ptV2,w) ;
1763       }
1764       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1765         FillHistogram(Form("hMCMassPtDisp2_cen%d",fCenBin),m,pt,w) ;
1766         FillHistogram(Form("hMCMassPtDisp2core_cen%d",fCenBin),mV2,ptV2,w) ;
1767         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1768           FillHistogram(Form("hMCMassPtBoth2_cen%d",fCenBin),m,pt,w) ;
1769           FillHistogram(Form("hMCMassPtBoth2core_cen%d",fCenBin),mV2,ptV2,w) ;
1770         }
1771       }
1772       if(ph1->IsDispOK() && ph2->IsDispOK()){
1773         FillHistogram(Form("hMCMassPtDisp_cen%d",fCenBin),m,pt,w) ;
1774         FillHistogram(Form("hMCMassPtDispcore_cen%d",fCenBin),mV2,ptV2,w) ;
1775         if(ph1->IsPhoton()&& ph2->IsPhoton()){
1776           FillHistogram(Form("hMCMassPtDispwou_cen%d",fCenBin),m,pt,w) ;
1777         }
1778         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1779           FillHistogram(Form("hMCMassPtBoth_cen%d",fCenBin),m,pt,w) ;
1780           FillHistogram(Form("hMCMassPtBothcore_cen%d",fCenBin),mV2,ptV2,w) ;
1781         }
1782       }
1783     } // end of loop i2
1784   } // end of loop i1
1785 }
1786 //____________________________________________________________________________
1787 Double_t  AliAnalysisTaskPi0DiffEfficiency::CoreEnergy(AliPHOSAodCluster * clu){
1788   //calculate energy of the cluster in the circle with radius distanceCut around the maximum
1789
1790   //Can not use already calculated coordinates?
1791   //They have incidence correction...
1792   const Double_t distanceCut =3.5 ;
1793   const Double_t logWeight=4.5 ;
1794
1795   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1796 // Calculates the center of gravity in the local PHOS-module coordinates
1797   Float_t wtot = 0;
1798   Double_t xc[100]={0} ;
1799   Double_t zc[100]={0} ;
1800   Double_t x = 0 ;
1801   Double_t z = 0 ;
1802   Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
1803   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1804     Int_t relid[4] ;
1805     Float_t xi ;
1806     Float_t zi ;
1807     fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
1808     fPHOSGeo->RelPosInModule(relid, xi, zi);
1809     xc[iDigit]=xi ;
1810     zc[iDigit]=zi ;
1811     if (clu->E()>0 && elist[iDigit]>0) {
1812       Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1813       x    += xc[iDigit] * w ;
1814       z    += zc[iDigit] * w ;
1815       wtot += w ;
1816     }
1817   }
1818   if (wtot>0) {
1819     x /= wtot ;
1820     z /= wtot ;
1821   }
1822   Double_t coreE=0. ;
1823   for(Int_t iDigit=0; iDigit < mulDigit; iDigit++) {
1824     Double_t distance = TMath::Sqrt((xc[iDigit]-x)*(xc[iDigit]-x)+(zc[iDigit]-z)*(zc[iDigit]-z)) ;
1825     if(distance < distanceCut)
1826       coreE += elist[iDigit] ;
1827   }
1828   //Apply non-linearity correction
1829   return (0.0241+1.0504*coreE+0.000249*coreE*coreE) ;
1830 }
1831 //_____________________________________________________________________________
1832 Bool_t AliAnalysisTaskPi0DiffEfficiency::TestLambda(Double_t pt,Double_t l1,Double_t l2){
1833
1834   //For R=3.5
1835    Double_t  l1Mean  = 1.170014 -0.059465/(1.+0.019343*pt+0.147544*pt*pt) ;
1836    Double_t  l2Mean = 1.626270 + 0.761554*exp(-1.213839*pt)-0.020027*pt ;
1837    Double_t  l1Sigma = 0.133409 + 0.261307*exp(-0.636874*pt)-0.002849*pt ;
1838    Double_t  l2Sigma = 0.289698 + 0.459400*exp(-1.214242*pt)-0.012578*pt ;
1839    Double_t  c=-0.124103 ;
1840 /*
1841   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1842   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1843   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1844   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1845   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1846 */
1847   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1848               0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1849               0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1850   return (R2<2.5*2.5) ;
1851
1852 }
1853 //_____________________________________________________________________________
1854 Bool_t AliAnalysisTaskPi0DiffEfficiency::TestLambda2(Double_t pt,Double_t l1,Double_t l2){
1855
1856 //For R=4.5
1857   Double_t   l1Mean  = 1.150200 + 0.097886/(1.+1.486645*pt+0.000038*pt*pt) ;
1858   Double_t   l2Mean = 1.574706 + 0.997966*exp(-0.895075*pt)-0.010666*pt ;
1859   Double_t   l1Sigma = 0.100255 + 0.337177*exp(-0.517684*pt)+0.001170*pt ;
1860   Double_t   l2Sigma = 0.232580 + 0.573401*exp(-0.735903*pt)-0.002325*pt ;
1861   Double_t   c = -0.110983 -0.017353/(1.-1.836995*pt+0.934517*pt*pt) ;
1862
1863 /*
1864   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1865   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1866   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1867   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1868   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1869 */
1870   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1871               0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1872               0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1873   return (R2<2.5*2.5) ;
1874
1875 }
1876 //____________________________________________________________________________
1877 Double_t AliAnalysisTaskPi0DiffEfficiency::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge){
1878   //Parameterization of LHC10h period
1879   //_true if neutral_
1880
1881   Double_t meanX=0;
1882   Double_t meanZ=0.;
1883   Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
1884               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);
1885   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) ;
1886   AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
1887   if(!event)AliFatal("Can not get ESD event") ;
1888   Double_t mf = event->GetMagneticField(); //Positive for ++ and negative for --
1889
1890   if(mf<0.){ //field --
1891     meanZ = -0.468318 ;
1892     if(charge>0)
1893       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)) ;
1894     else
1895       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)) ;
1896   }
1897   else{ //Field ++
1898     meanZ= -0.468318;
1899     if(charge>0)
1900       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)) ;
1901     else
1902       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)) ;
1903   }
1904
1905   Double_t rz=(dz-meanZ)/sz ;
1906   Double_t rx=(dx-meanX)/sx ;
1907   return TMath::Sqrt(rx*rx+rz*rz) ;
1908 }
1909 //____________________________________________________________________________
1910 Bool_t AliAnalysisTaskPi0DiffEfficiency::TestTOF(Double_t t, Double_t e){
1911
1912   Double_t sigma = TMath::Sqrt(2.23183e-09*2.23183e-09
1913                              +2.24611e-09*2.24611e-09/e
1914                              +5.65054e-09*5.65054e-09/e/e) ;
1915   sigma=TMath::Min(20.e-9,sigma) ; //for the soft (<400 MeV) part
1916   Double_t mean=1.1e-9 ;
1917   if(TMath::Abs(t-mean)<2.*sigma)
1918     return kTRUE ;
1919   else
1920     if(TMath::Abs(t-mean+100.e-9)<2.*sigma)
1921       return kTRUE ;
1922
1923   return kFALSE ;
1924 }
1925 //_____________________________________________________________________________
1926 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x)const{
1927   //FillHistogram
1928   TObject * tmp = fOutputContainer->FindObject(key) ;
1929   if(!tmp){
1930     AliInfo(Form("can not find histogram <%s> ",key)) ;
1931     return ;
1932   }
1933   if(tmp->IsA() == TClass::GetClass("TH1I")){
1934     ((TH1I*)tmp)->Fill(x) ;
1935     return ;
1936   }
1937   if(tmp->IsA() == TClass::GetClass("TH1F")){
1938     ((TH1F*)tmp)->Fill(x) ;
1939     return ;
1940   }
1941   if(tmp->IsA() == TClass::GetClass("TH1D")){
1942     ((TH1D*)tmp)->Fill(x) ;
1943     return ;
1944   }
1945   AliInfo(Form("can not find 1D histogram <%s> ",key)) ;
1946 }
1947 //_____________________________________________________________________________
1948 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x,Double_t y)const{
1949   //FillHistogram
1950   TObject * tmp = fOutputContainer->FindObject(key) ;
1951   if(!tmp){
1952     AliInfo(Form("can not find histogram <%s> ",key)) ;
1953     return ;
1954   }
1955   if(tmp->IsA() == TClass::GetClass("TH1F")){
1956     ((TH1F*)tmp)->Fill(x,y) ;
1957     return ;
1958   }
1959   if(tmp->IsA() == TClass::GetClass("TH2F")){
1960     ((TH2F*)tmp)->Fill(x,y) ;
1961     return ;
1962   }
1963   AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
1964 }
1965
1966 //_____________________________________________________________________________
1967 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
1968   //Fills 1D histograms with key
1969   TObject * tmp = fOutputContainer->FindObject(key) ;
1970   if(!tmp){
1971     AliInfo(Form("can not find histogram <%s> ",key)) ;
1972     return ;
1973   }
1974   if(tmp->IsA() == TClass::GetClass("TH2F")){
1975     ((TH2F*)tmp)->Fill(x,y,z) ;
1976     return ;
1977   }
1978   if(tmp->IsA() == TClass::GetClass("TH3F")){
1979     ((TH3F*)tmp)->Fill(x,y,z) ;
1980     return ;
1981   }
1982 }
1983 //_____________________________________________________________________________
1984 void AliAnalysisTaskPi0DiffEfficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z, Double_t w) const{
1985   //Fills 1D histograms with key
1986   TObject * tmp = fOutputContainer->FindObject(key) ;
1987   if(!tmp){
1988     AliInfo(Form("can not find histogram <%s> ",key)) ;
1989     return ;
1990   }
1991   if(tmp->IsA() == TClass::GetClass("TH3F")){
1992     ((TH3F*)tmp)->Fill(x,y,z,w) ;
1993     return ;
1994   }
1995 }
1996 //_____________________________________________________________________________
1997 Bool_t AliAnalysisTaskPi0DiffEfficiency::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
1998 {
1999   //Check if this channel belogs to the good ones
2000
2001   if(strcmp(det,"PHOS")==0){
2002     if(mod>5 || mod<1){
2003       AliError(Form("No bad map for PHOS module %d ",mod)) ;
2004       return kTRUE ;
2005     }
2006     if(!fPHOSBadMap[mod]){
2007       AliError(Form("No Bad map for PHOS module %d",mod)) ;
2008       return kTRUE ;
2009     }
2010     if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
2011       return kFALSE ;
2012     else
2013       return kTRUE ;
2014   }
2015   else{
2016     AliError(Form("Can not find bad channels for detector %s ",det)) ;
2017   }
2018   return kTRUE ;
2019 }
2020 //_____________________________________________________________________________
2021 Double_t AliAnalysisTaskPi0DiffEfficiency::PrimaryWeight(Double_t x){
2022
2023    Double_t w=1 ;
2024
2025
2026    //Parameterization of LHC10h data from Jan 2013 (pi0 spectrum)
2027    //Should be consistend with spectrum parameterization used in simulation
2028    if(fCenBin==0) //0-5
2029      w = (0.561741+0.332841*x-0.007082*x*x)/(1.-0.447804*x+0.157830*x*x)+0.080394*x ;
2030    if(fCenBin==1) //5-10
2031      w = (0.659096+0.101701*x+0.042395*x*x)/(1.-0.470110*x+0.154665*x*x)+0.052932*x ;
2032    if(fCenBin==2) //10-20
2033      w = (0.615575+0.005621*x+0.069263*x*x)/(1.-0.485422*x+0.160822*x*x)+0.040865*x ;
2034    if(fCenBin==3) //20-40
2035      w = (0.441240+0.158358*x+0.059458*x*x)/(1.-0.332609*x+0.147528*x*x)+0.037926*x ;
2036    if(fCenBin==4) //40-60
2037      w = (0.467895-0.001113*x+0.029610*x*x)/(1.-0.266502*x+0.065105*x*x)+0.025431*x ;
2038    if(fCenBin==5) //60-80
2039      w = (0.465204-0.139362*x+0.043500*x*x)/(1.-0.371689*x+0.067606*x*x)+0.006519*x ;
2040
2041 /*
2042   //Parameterization of photon spectrum 25.02
2043   if(fCenBin==0) //0-5
2044      w=(0.870487-0.494032*x+0.076334*x*x+0.001065*x*x*x)/(1.-0.646014*x+0.113839*x*x);
2045   if(fCenBin==1) //5-10
2046      w=(-8.310403+15.767226*x-2.167756*x*x+0.184356*x*x*x)/(1.+4.556793*x+0.980941*x*x);
2047   if(fCenBin==2) //10-20
2048      w=(-5.281594+7.477165*x-0.688609*x*x+0.097601*x*x*x)/(1.+1.102693*x+0.882454*x*x);
2049   if(fCenBin==3) //20-40
2050      w=(0.789472-4.750155*x+4.381545*x*x-0.029239*x*x*x)/(1.-3.738304*x+3.328318*x*x);
2051   if(fCenBin==4) //40-60
2052      w=(0.876792-0.491379*x+0.130872*x*x-0.001390*x*x*x)/(1.+-0.511934*x+0.112705*x*x);
2053   if(fCenBin==5) //60-80
2054      w=(0.719912-0.167292*x+0.017196*x*x+0.000861*x*x*x)/(1.-0.336146*x+0.037731*x*x);
2055 */
2056   return TMath::Max(0.,w) ;
2057 }
2058
2059
2060