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
  15 #include "AliAODMCParticle.h"
  16 #include "AliAnalysisManager.h"
  17 #include "AliMCEventHandler.h"
  18 #include "AliMCEvent.h"
  19 #include "AliStack.h"
  20 #include "AliAnalysisTaskSE.h"
  21 #include "AliAnalysisTaskPi0DiffEfficiency.h"
  22 #include "AliCaloPhoton.h"
  23 #include "AliPHOSGeometry.h"
  24 #include "AliPHOSAodCluster.h"
  25 #include "AliPHOSCalibData.h"
  26 #include "AliAODEvent.h"
  27 #include "AliAODCaloCluster.h"
  28 #include "AliAODVertex.h"
  29 #include "AliESDtrackCuts.h"
  30 #include "AliLog.h"
  31 #include "AliPID.h"
  32 #include "AliCDBManager.h"
  33 #include "AliCentrality.h"
  34
  35 // Analysis task to calculate pi0 registration efficiency
  36 // as a difference between spectra before and after embedding
  37 // Authors: Dmitri Peressounko
  38 // Date   : Aug.2011
  39
  40 ClassImp(AliAnalysisTaskPi0DiffEfficiency)
  41
  42 //________________________________________________________________________
  43 AliAnalysisTaskPi0DiffEfficiency::AliAnalysisTaskPi0DiffEfficiency(const char *name)
  44 : AliAnalysisTaskPi0Efficiency(name),
  45   fPHOSEvent1(0),
  46   fPHOSEvent2(0)
  47 {
  48
  49
  50
  51 }
  52
  53 //________________________________________________________________________
  54 void AliAnalysisTaskPi0DiffEfficiency::UserCreateOutputObjects()
  55 {
  56   // Create histograms
  57   // Called once
  58
  59   // ESD histograms
  60   if(fOutputContainer != NULL){
  61     delete fOutputContainer;
  62   }
  63   fOutputContainer = new TList();
  64   fOutputContainer->SetOwner(kTRUE);
  65
  66   //Event selection
  67   fOutputContainer->Add(new TH1F("hSelEvents","Event celection", 10,0.,10.)) ;
  68
  69   //vertex distribution
  70   fOutputContainer->Add(new TH1F("hZvertex","Z vertex position", 50,-25.,25.)) ;
  71
  72   //Centrality
  73   fOutputContainer->Add(new TH1F("hCentrality","Event centrality", 100,0.,100.)) ;
  74
  75   //QA histograms
  76   fOutputContainer->Add(new TH2F("hCluM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
  77   fOutputContainer->Add(new TH2F("hCluM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
  78   fOutputContainer->Add(new TH2F("hCluM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
  79
  80   Int_t nM       = 500;
  81   Double_t mMin  = 0.0;
  82   Double_t mMax  = 1.0;
  83   Int_t nPt      = 200;
  84   Double_t ptMin = 0;
  85   Double_t ptMax = 20;
  86
  87   char key[55] ;
  88   for(Int_t cent=0; cent<6; cent++){
  89     //Single photon
  90     snprintf(key,55,"hPhotAll_cen%d",cent) ;
  91     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
  92     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
  93     snprintf(key,55,"hPhotAllcore_cen%d",cent) ;
  94     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
  95     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
  96     snprintf(key,55,"hPhotAllwou_cen%d",cent) ;
  97     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
  98     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
  99     snprintf(key,55,"hPhotCPV_cen%d",cent) ;
 100     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 101     snprintf(key,55,"hPhotCPVcore_cen%d",cent) ;
 102     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 103     snprintf(key,55,"hPhotCPV2_cen%d",cent) ;
 104     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 105     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 106     snprintf(key,55,"hPhotDisp_cen%d",cent) ;
 107     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 108     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 109     snprintf(key,55,"hPhotDisp2_cen%d",cent) ;
 110     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 111     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 112     snprintf(key,55,"hPhotDispwou_cen%d",cent) ;
 113     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 114     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 115     snprintf(key,55,"hPhotBoth_cen%d",cent) ;
 116     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 117     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 118     snprintf(key,55,"hPhotBothcore_cen%d",cent) ;
 119     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 120     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 121
 122     snprintf(key,55,"hNegPhotAll_cen%d",cent) ;
 123     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 124     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 125     snprintf(key,55,"hNegPhotAllcore_cen%d",cent) ;
 126     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 127     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 128     snprintf(key,55,"hNegPhotAllwou_cen%d",cent) ;
 129     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 130     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 131     snprintf(key,55,"hNegPhotCPV_cen%d",cent) ;
 132     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 133     snprintf(key,55,"hNegPhotCPVcore_cen%d",cent) ;
 134     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 135     snprintf(key,55,"hNegPhotCPV2_cen%d",cent) ;
 136     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 137     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 138     snprintf(key,55,"hNegPhotDisp_cen%d",cent) ;
 139     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 140     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 141     snprintf(key,55,"hNegPhotDisp2_cen%d",cent) ;
 142     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 143     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 144     snprintf(key,55,"hNegPhotDispwou_cen%d",cent) ;
 145     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 146     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 147     snprintf(key,55,"hNegPhotBoth_cen%d",cent) ;
 148     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 149     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 150     snprintf(key,55,"hNegPhotBothcore_cen%d",cent) ;
 151     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 152     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 153
 154     snprintf(key,55,"hOldMassPtAll_cen%d",cent) ;
 155     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 156     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 157     snprintf(key,55,"hOldMassPtAllcore_cen%d",cent) ;
 158     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 159     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 160     snprintf(key,55,"hOldMassPtAllwou_cen%d",cent) ;
 161     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 162     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 163     snprintf(key,55,"hOldMassPtCPV_cen%d",cent) ;
 164     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 165     snprintf(key,55,"hOldMassPtCPVcore_cen%d",cent) ;
 166     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 167     snprintf(key,55,"hOldMassPtCPV2_cen%d",cent) ;
 168     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 169     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 170     snprintf(key,55,"hOldMassPtDisp_cen%d",cent) ;
 171     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 172     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 173     snprintf(key,55,"hOldMassPtDisp2_cen%d",cent) ;
 174     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 175     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 176     snprintf(key,55,"hOldMassPtDispwou_cen%d",cent) ;
 177     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 178     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 179     snprintf(key,55,"hOldMassPtBoth_cen%d",cent) ;
 180     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 181     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 182     snprintf(key,55,"hOldMassPtBothcore_cen%d",cent) ;
 183     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 184     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 185
 186     snprintf(key,55,"hNewMassPtAll_cen%d",cent) ;
 187     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 188     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 189     snprintf(key,55,"hNewMassPtAllcore_cen%d",cent) ;
 190     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 191     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 192     snprintf(key,55,"hNewMassPtAllwou_cen%d",cent) ;
 193     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 194     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 195     snprintf(key,55,"hNewMassPtCPV_cen%d",cent) ;
 196     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 197     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 198     snprintf(key,55,"hNewMassPtCPVcore_cen%d",cent) ;
 199     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 200     snprintf(key,55,"hNewMassPtCPV2_cen%d",cent) ;
 201     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 202     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 203     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 204     snprintf(key,55,"hNewMassPtDisp_cen%d",cent) ;
 205     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 206     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 207     snprintf(key,55,"hNewMassPtDisp2_cen%d",cent) ;
 208     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 209     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 210     snprintf(key,55,"hNewMassPtDispwou_cen%d",cent) ;
 211     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 212     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 213     snprintf(key,55,"hNewMassPtBoth_cen%d",cent) ;
 214     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 215     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 216     snprintf(key,55,"hNewMassPtBothcore_cen%d",cent) ;
 217     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 218     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 219
 220     snprintf(key,55,"hNegMassPtAll_cen%d",cent) ;
 221     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 222     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 223     snprintf(key,55,"hNegMassPtAllcore_cen%d",cent) ;
 224     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 225     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 226     snprintf(key,55,"hNegMassPtAllwou_cen%d",cent) ;
 227     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 228     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 229     snprintf(key,55,"hNegMassPtCPV_cen%d",cent) ;
 230     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 231     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 232     snprintf(key,55,"hNegMassPtCPVcore_cen%d",cent) ;
 233     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 234     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 235     snprintf(key,55,"hNegMassPtCPV2_cen%d",cent) ;
 236     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 237     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 238     snprintf(key,55,"hNegMassPtDisp_cen%d",cent) ;
 239     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 240     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 241     snprintf(key,55,"hNegMassPtDisp2_cen%d",cent) ;
 242     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 243     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 244     snprintf(key,55,"hNegMassPtDispwou_cen%d",cent) ;
 245     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 246     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 247     snprintf(key,55,"hNegMassPtBoth_cen%d",cent) ;
 248     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 249     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 250     snprintf(key,55,"hNegMassPtBothcore_cen%d",cent) ;
 251     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 252     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 253
 254
 255     snprintf(key,55,"hMassPtAll_cen%d",cent) ;
 256     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 257     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 258     snprintf(key,55,"hMassPtAllcore_cen%d",cent) ;
 259     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 260     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 261     snprintf(key,55,"hMassPtAllwou_cen%d",cent) ;
 262     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 263     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 264     snprintf(key,55,"hMassPtCPV_cen%d",cent) ;
 265     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 266     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 267     snprintf(key,55,"hMassPtCPVcore_cen%d",cent) ;
 268     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 269     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 270     snprintf(key,55,"hMassPtCPV2_cen%d",cent) ;
 271     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 272     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 273     snprintf(key,55,"hMassPtDisp_cen%d",cent) ;
 274     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 275     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 276     snprintf(key,55,"hMassPtDisp2_cen%d",cent) ;
 277     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 278     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 279     snprintf(key,55,"hMassPtDispwou_cen%d",cent) ;
 280     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 281     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 282     snprintf(key,55,"hMassPtBoth_cen%d",cent) ;
 283     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 284     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 285     snprintf(key,55,"hMassPtBothcore_cen%d",cent) ;
 286     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 287     ((TH2F*)fOutputContainer->Last())->Sumw2() ;
 288
 289     snprintf(key,55,"hMassPtAll_a07_cen%d",cent) ;
 290     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 291     snprintf(key,55,"hMassPtCPV_a07_cen%d",cent) ;
 292     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 293     snprintf(key,55,"hMassPtCPV2_a07_cen%d",cent) ;
 294     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 295     snprintf(key,55,"hMassPtDisp_a07_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,"hMassPtBoth_a07_cen%d",cent) ;
 298     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 299
 300     snprintf(key,55,"hMassPtAll_a08_cen%d",cent) ;
 301     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 302     snprintf(key,55,"hMassPtCPV_a08_cen%d",cent) ;
 303     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 304     snprintf(key,55,"hMassPtCPV2_a08_cen%d",cent) ;
 305     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 306     snprintf(key,55,"hMassPtDisp_a08_cen%d",cent) ;
 307     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 308     snprintf(key,55,"hMassPtBoth_a08_cen%d",cent) ;
 309     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 310
 311     snprintf(key,55,"hMassPtAll_a09_cen%d",cent) ;
 312     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 313     snprintf(key,55,"hMassPtCPV_a09_cen%d",cent) ;
 314     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 315     snprintf(key,55,"hMassPtCPV2_a09_cen%d",cent) ;
 316     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 317     snprintf(key,55,"hMassPtDisp_a09_cen%d",cent) ;
 318     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 319     snprintf(key,55,"hMassPtBoth_a09_cen%d",cent) ;
 320     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 321
 322
 323     //Mixed
 324     snprintf(key,55,"hMiMassPtAll_cen%d",cent) ;
 325     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 326     snprintf(key,55,"hMiMassPtAllcore_cen%d",cent) ;
 327     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 328     snprintf(key,55,"hMiMassPtAllwou_cen%d",cent) ;
 329     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 330     snprintf(key,55,"hMiMassPtCPV_cen%d",cent) ;
 331     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 332     snprintf(key,55,"hMiMassPtCPVcore_cen%d",cent) ;
 333     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 334     snprintf(key,55,"hMiMassPtCPV2_cen%d",cent) ;
 335     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 336     snprintf(key,55,"hMiMassPtDisp_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,"hMiMassPtDisp2_cen%d",cent) ;
 339     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 340     snprintf(key,55,"hMiMassPtDispwou_cen%d",cent) ;
 341     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 342     snprintf(key,55,"hMiMassPtBoth_cen%d",cent) ;
 343     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 344     snprintf(key,55,"hMiMassPtBothcore_cen%d",cent) ;
 345     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 346
 347     snprintf(key,55,"hMiMassPtAll_a07_cen%d",cent) ;
 348     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 349     snprintf(key,55,"hMiMassPtCPV_a07_cen%d",cent) ;
 350     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 351     snprintf(key,55,"hMiMassPtCPV2_a07_cen%d",cent) ;
 352     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 353     snprintf(key,55,"hMiMassPtDisp_a07_cen%d",cent) ;
 354     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 355     snprintf(key,55,"hMiMassPtBoth_a07_cen%d",cent) ;
 356     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 357
 358     snprintf(key,55,"hMiMassPtAll_a08_cen%d",cent) ;
 359     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 360     snprintf(key,55,"hMiMassPtCPV_a08_cen%d",cent) ;
 361     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 362     snprintf(key,55,"hMiMassPtCPV2_a08_cen%d",cent) ;
 363     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 364     snprintf(key,55,"hMiMassPtDisp_a08_cen%d",cent) ;
 365     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 366     snprintf(key,55,"hMiMassPtBoth_a08_cen%d",cent) ;
 367     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 368
 369      snprintf(key,55,"hMiMassPtAll_a09_cen%d",cent) ;
 370     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 371     snprintf(key,55,"hMiMassPtCPV_a09_cen%d",cent) ;
 372     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 373     snprintf(key,55,"hMiMassPtCPV2_a09_cen%d",cent) ;
 374     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 375     snprintf(key,55,"hMiMassPtDisp_a09_cen%d",cent) ;
 376     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 377     snprintf(key,55,"hMiMassPtBoth_a09_cen%d",cent) ;
 378     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 379
 380     snprintf(key,55,"hMCMassPtAll_cen%d",cent) ;
 381     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 382     snprintf(key,55,"hMCMassPtAllcore_cen%d",cent) ;
 383     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 384     snprintf(key,55,"hMCMassPtAllwou_cen%d",cent) ;
 385     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 386     snprintf(key,55,"hMCMassPtCPV_cen%d",cent) ;
 387     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 388     snprintf(key,55,"hMCMassPtCPVcore_cen%d",cent) ;
 389     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 390     snprintf(key,55,"hMCMassPtCPV2_cen%d",cent) ;
 391     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 392     snprintf(key,55,"hMCMassPtDisp_cen%d",cent) ;
 393     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 394     snprintf(key,55,"hMCMassPtDisp2_cen%d",cent) ;
 395     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 396     snprintf(key,55,"hMCMassPtDispwou_cen%d",cent) ;
 397     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 398     snprintf(key,55,"hMCMassPtBoth_cen%d",cent) ;
 399     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 400     snprintf(key,55,"hMCMassPtBothcore_cen%d",cent) ;
 401     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 402
 403     snprintf(key,55,"hMCMassPtAll_a07_cen%d",cent) ;
 404     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 405     snprintf(key,55,"hMCMassPtCPV_a07_cen%d",cent) ;
 406     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 407     snprintf(key,55,"hMCMassPtCPV2_a07_cen%d",cent) ;
 408     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 409     snprintf(key,55,"hMCMassPtDisp_a07_cen%d",cent) ;
 410     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 411     snprintf(key,55,"hMCMassPtBoth_a07_cen%d",cent) ;
 412     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 413
 414     snprintf(key,55,"hMCMassPtAll_a08_cen%d",cent) ;
 415     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 416     snprintf(key,55,"hMCMassPtCPV_a08_cen%d",cent) ;
 417     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 418     snprintf(key,55,"hMCMassPtCPV2_a08_cen%d",cent) ;
 419     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 420     snprintf(key,55,"hMCMassPtDisp_a08_cen%d",cent) ;
 421     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 422     snprintf(key,55,"hMCMassPtBoth_a08_cen%d",cent) ;
 423     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 424
 425     snprintf(key,55,"hMCMassPtAll_a09_cen%d",cent) ;
 426     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 427     snprintf(key,55,"hMCMassPtCPV_a09_cen%d",cent) ;
 428     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 429     snprintf(key,55,"hMCMassPtCPV2_a09_cen%d",cent) ;
 430     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 431     snprintf(key,55,"hMCMassPtDisp_a09_cen%d",cent) ;
 432     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 433     snprintf(key,55,"hMCMassPtBoth_a09_cen%d",cent) ;
 434     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 435
 436     //Single photon
 437     snprintf(key,55,"hMCPhotAll_cen%d",cent) ;
 438     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 439     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 440     snprintf(key,55,"hMCPhotAllcore_cen%d",cent) ;
 441     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 442     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 443     snprintf(key,55,"hMCPhotAllwou_cen%d",cent) ;
 444     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 445     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 446     snprintf(key,55,"hMCPhotCPV_cen%d",cent) ;
 447     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 448     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 449     snprintf(key,55,"hMCPhotCPVcore_cen%d",cent) ;
 450     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 451     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 452     snprintf(key,55,"hMCPhotCPV2_cen%d",cent) ;
 453     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 454     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 455     snprintf(key,55,"hMCPhotDisp_cen%d",cent) ;
 456     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 457     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 458     snprintf(key,55,"hMCPhotDisp2_cen%d",cent) ;
 459     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 460     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 461     snprintf(key,55,"hMCPhotDispwou_cen%d",cent) ;
 462     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 463     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 464     snprintf(key,55,"hMCPhotBoth_cen%d",cent) ;
 465     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 466     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 467     snprintf(key,55,"hMCPhotBothcore_cen%d",cent) ;
 468     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 469     ((TH1F*)fOutputContainer->Last())->Sumw2() ;
 470
 471   }
 472
 473   fOutputContainer->Add(new TH2F("hMCPi0M11","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 474   fOutputContainer->Add(new TH2F("hMCPi0M22","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 475   fOutputContainer->Add(new TH2F("hMCPi0M33","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 476   fOutputContainer->Add(new TH2F("hMCPi0M12","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 477   fOutputContainer->Add(new TH2F("hMCPi0M13","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 478   fOutputContainer->Add(new TH2F("hMCPi0M23","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 479
 480   //MC
 481   for(Int_t cent=0; cent<6; cent++){
 482     snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
 483     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
 484     snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
 485     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
 486     snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
 487     fOutputContainer->Add(new TH1F("hMC_rap_eta","Rapidity eta",200,-1.,1.)) ;
 488     snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
 489     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
 490     snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
 491     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
 492     snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
 493     fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
 494     snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
 495     fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
 496     snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
 497     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
 498     snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
 499     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 500     snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
 501     fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
 502     snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
 503     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 504     snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
 505     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 506   }
 507
 508   PostData(1, fOutputContainer);
 509
 510 }
 511
 512 //________________________________________________________________________
 513 void AliAnalysisTaskPi0DiffEfficiency::UserExec(Option_t *)
 514 {
 515   // Main loop, called for each event
 516   // Analyze ESD/AOD
 517
 518   FillHistogram("hSelEvents",0.5) ;
 519
 520   AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
 521   if (!event) {
 522     Printf("ERROR: Could not retrieve event");
 523     PostData(1, fOutputContainer);
 524     return;
 525   }
 526
 527   FillHistogram("hSelEvents",1.5) ;
 528   AliAODHeader *header = event->GetHeader() ;
 529
 530   // Checks if we have a primary vertex
 531   // Get primary vertices form ESD
 532   const AliAODVertex *esdVertex5 = event->GetPrimaryVertex();
 533
 534  // don't rely on ESD vertex, assume (0,0,0)
 535   Double_t vtx0[3] ={0.,0.,0.};
 536
 537
 538   FillHistogram("hZvertex",esdVertex5->GetZ());
 539   if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
 540     PostData(1, fOutputContainer);
 541     return;
 542   }
 543   FillHistogram("hSelEvents",2.5) ;
 544
 545   //Vtx class z-bin
 546   //  Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
 547   //  if(zvtx<0)zvtx=0 ;
 548   //  if(zvtx>9)zvtx=9 ;
 549   Int_t zvtx=0 ;
 550
 551 //  fCentrality=header->GetCentralityP()->GetCentralityPercentile("V0M"); // returns the centrality percentile,
 552 //                                                          //a float from 0 to 100 (or to the trigger efficiency)
 553    fCentrality=header->GetZDCN2Energy() ;
 554
 555   if( fCentrality < 0. || fCentrality>80.){
 556     PostData(1, fOutputContainer);
 557     return;
 558   }
 559   FillHistogram("hSelEvents",3.5) ;
 560   Float_t bins[7]={0.,5.,10.,20.,40.,60.,80.} ;
 561   fCenBin=0 ;
 562   while(fCenBin<6 && fCentrality > bins[fCenBin+1])
 563     fCenBin++ ;
 564
 565
 566   //reaction plain
 567   fRPfull= header->GetZDCN1Energy() ;
 568   if(fRPfull==999){ //reaction plain was not defined
 569     PostData(1, fOutputContainer);
 570     return;
 571   }
 572
 573   FillHistogram("hSelEvents",4.5) ;
 574   //All event selections done
 575   FillHistogram("hCentrality",fCentrality) ;
 576   //Reaction plain is defined in the range (-pi/2;pi/2)
 577   //We have 10 bins
 578   Int_t irp=Int_t(10.*fRPfull/TMath::Pi());
 579   if(irp>9)irp=9 ;
 580
 581   if(!fPHOSEvents[zvtx][fCenBin][irp])
 582     fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
 583   TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
 584
 585   // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
 586   if(fEventCounter == 0) {
 587     for(Int_t mod=0; mod<5; mod++) {
 588       const TGeoHMatrix* m =header->GetPHOSMatrix(mod) ;
 589       fPHOSGeo->SetMisalMatrix(m,mod) ;
 590       Printf("PHOS geo matrix for module # %d is set: %p\n", mod,m);
 591     }
 592     fEventCounter++ ;
 593   }
 594
 595   ProcessMC() ;
 596
 597   if(fPHOSEvent1){
 598     fPHOSEvent1->Clear() ;
 599     fPHOSEvent2->Clear() ;
 600   }
 601   else{
 602     fPHOSEvent1 = new TClonesArray("AliCaloPhoton",200) ;
 603     fPHOSEvent2 = new TClonesArray("AliCaloPhoton",200) ;
 604   }
 605
 606   TClonesArray * clustersEmb = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
 607   AliAODCaloCells * cellsEmb = (AliAODCaloCells *)event->FindListObject("EmbeddedPHOScells") ;
 608   TClonesArray * clustersOld = event->GetCaloClusters() ;
 609   AliAODCaloCells * cellsOld = event->GetPHOSCells() ;
 610   TVector3 vertex(vtx0);
 611   char key[55] ;
 612   //Before Embedding
 613   Int_t multClustOld = clustersOld->GetEntriesFast();
 614   Int_t multClustEmb = clustersEmb->GetEntriesFast();
 615   Int_t inPHOSold=0 ;
 616   Int_t inPHOSemb=0 ;
 617   for (Int_t i=0; i<multClustOld; i++) {
 618     AliAODCaloCluster *clu = (AliAODCaloCluster*)clustersOld->At(i);
 619     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
 620
 621     Bool_t survive=kFALSE ;
 622     for(Int_t ii=0;(ii<multClustEmb)&&(!survive);ii++){
 623        AliAODCaloCluster *clu2 = (AliAODCaloCluster*)clustersEmb->At(ii);
 624        survive=IsSameCluster(clu,clu2);
 625     }
 626
 627
 628     Float_t  position[3];
 629     clu->GetPosition(position);
 630     TVector3 global(position) ;
 631     Int_t relId[4] ;
 632     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 633     Int_t mod  = relId[0] ;
 634     Int_t cellX = relId[2];
 635     Int_t cellZ = relId[3] ;
 636     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
 637       continue ;
 638     if(clu->GetNCells()<3)
 639       continue ;
 640
 641     snprintf(key,55,"hCluM%d",mod) ;
 642     FillHistogram(key,cellX,cellZ,1.);
 643
 644     TLorentzVector pv1 ;
 645     clu->GetMomentum(pv1 ,vtx0);
 646
 647     if(inPHOSold>=fPHOSEvent1->GetSize()){
 648       fPHOSEvent1->Expand(inPHOSold+50) ;
 649     }
 650     AliCaloPhoton * ph = new((*fPHOSEvent1)[inPHOSold]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
 651     ph->SetModule(mod) ;
 652     AliPHOSAodCluster cluPHOS1(*clu);
 653     cluPHOS1.Recalibrate(fPHOSCalibData,cellsOld); // modify the cell energies
 654     Double_t ecore=CoreEnergy(&cluPHOS1) ;
 655     pv1*= ecore/pv1.E() ;
 656     ph->SetMomV2(&pv1) ;
 657     ph->SetNCells(clu->GetNCells());
 658     ph->SetDispBit(TestLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
 659     ph->SetDisp2Bit(TestLambda2(clu->E(),clu->GetM20(),clu->GetM02())) ;
 660     ph->SetCPVBit(clu->GetEmcCpvDistance()>2.) ;
 661     ph->SetCPV2Bit(clu->GetEmcCpvDistance()>4.) ;
 662     if(!survive) //this cluster found in list after embedding, skipping it
 663       ph->SetTagged(1) ;
 664     ph->SetPhoton(clu->GetNExMax()<2); // Remember, if it is unfolded
 665
 666     inPHOSold++ ;
 667   }
 668
 669   for (Int_t i=0; i<multClustEmb; i++) {
 670     AliAODCaloCluster *clu = (AliAODCaloCluster*)clustersEmb->At(i);
 671     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
 672
 673     Bool_t survive=kFALSE ;
 674     for(Int_t ii=0;(ii<multClustOld)&&(!survive);ii++){
 675        AliAODCaloCluster *clu2 = (AliAODCaloCluster*)clustersOld->At(ii);
 676        survive=IsSameCluster(clu,clu2);
 677     }
 678
 679     Float_t  position[3];
 680     clu->GetPosition(position);
 681     TVector3 global(position) ;
 682     Int_t relId[4] ;
 683     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 684     Int_t mod  = relId[0] ;
 685     Int_t cellX = relId[2];
 686     Int_t cellZ = relId[3] ;
 687     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
 688       continue ;
 689     if(clu->GetNCells()<3)
 690       continue ;
 691
 692     snprintf(key,55,"hCluM%d",mod) ;
 693     FillHistogram(key,cellX,cellZ,1.);
 694
 695     TLorentzVector pv1 ;
 696     clu->GetMomentum(pv1 ,vtx0);
 697
 698     if(inPHOSemb>=fPHOSEvent2->GetSize()){
 699       fPHOSEvent2->Expand(inPHOSemb+50) ;
 700     }
 701     AliCaloPhoton * ph = new((*fPHOSEvent2)[inPHOSemb]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
 702     ph->SetModule(mod) ;
 703     AliPHOSAodCluster cluPHOS1(*clu);
 704     cluPHOS1.Recalibrate(fPHOSCalibData,cellsEmb); // modify the cell energies
 705     Double_t ecore=CoreEnergy(&cluPHOS1) ;
 706     pv1*= ecore/pv1.E() ;
 707     ph->SetMomV2(&pv1) ;
 708     ph->SetNCells(clu->GetNCells());
 709     ph->SetDispBit(TestLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
 710     ph->SetDisp2Bit(TestLambda2(clu->E(),clu->GetM20(),clu->GetM02())) ;
 711     ph->SetCPVBit(clu->GetEmcCpvDistance()>2.) ;
 712     ph->SetCPV2Bit(clu->GetEmcCpvDistance()>4.) ;
 713     if(!survive) //this cluster found in list after embedding, skipping it
 714       ph->SetTagged(1) ;
 715     ph->SetPhoton(clu->GetNExMax()<2); // Remember, if it is unfolded
 716
 717     inPHOSemb++ ;
 718   }
 719
 720   //Single photon
 721   for (Int_t i1=0; i1<inPHOSold; i1++) {
 722     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent1->At(i1) ;
 723     if(!ph1->IsTagged())
 724       continue ;
 725     snprintf(key,55,"hPhotAll_cen%d",fCenBin) ;
 726     FillHistogram(key,ph1->Pt(),-1.) ;
 727     snprintf(key,55,"hNegPhotAll_cen%d",fCenBin) ;
 728     FillHistogram(key,ph1->Pt(),-1.) ;
 729     snprintf(key,55,"hPhotAllcore_cen%d",fCenBin) ;
 730     FillHistogram(key,ph1->GetMomV2()->Pt(),-1.) ;
 731     snprintf(key,55,"hNegPhotAllcore_cen%d",fCenBin) ;
 732     FillHistogram(key,ph1->GetMomV2()->Pt(),-1.) ;
 733     if(ph1->IsPhoton()){
 734       snprintf(key,55,"hPhotAllwou_cen%d",fCenBin) ;
 735       FillHistogram(key,ph1->Pt(),-1.) ;
 736       snprintf(key,55,"hNegPhotAllwou_cen%d",fCenBin) ;
 737       FillHistogram(key,ph1->Pt(),-1.) ;
 738     }
 739     if(ph1->IsCPVOK() ){
 740       snprintf(key,55,"hPhotCPV_cen%d",fCenBin) ;
 741       FillHistogram(key,ph1->Pt(),-1.) ;
 742       snprintf(key,55,"hNegPhotCPV_cen%d",fCenBin) ;
 743       FillHistogram(key,ph1->Pt(),-1.) ;
 744       snprintf(key,55,"hPhotCPVcore_cen%d",fCenBin) ;
 745       FillHistogram(key,ph1->GetMomV2()->Pt(),-1.) ;
 746       snprintf(key,55,"hNegPhotCPVcore_cen%d",fCenBin) ;
 747       FillHistogram(key,ph1->GetMomV2()->Pt(),-1.) ;
 748     }
 749     if(ph1->IsCPV2OK() ){
 750       snprintf(key,55,"hPhotCPV2_cen%d",fCenBin) ;
 751       FillHistogram(key,ph1->Pt(),-1.) ;
 752       snprintf(key,55,"hNegPhotCPV2_cen%d",fCenBin) ;
 753       FillHistogram(key,ph1->Pt(),-1.) ;
 754     }
 755     if(ph1->IsDisp2OK()){
 756       snprintf(key,55,"hPhotDisp2_cen%d",fCenBin) ;
 757       FillHistogram(key,ph1->Pt(),-1.) ;
 758       snprintf(key,55,"hNegPhotDisp2_cen%d",fCenBin) ;
 759       FillHistogram(key,ph1->Pt(),-1.) ;
 760     }
 761     if(ph1->IsDispOK()){
 762       snprintf(key,55,"hPhotDisp_cen%d",fCenBin) ;
 763       FillHistogram(key,ph1->Pt(),-1.) ;
 764       snprintf(key,55,"hNegPhotDisp_cen%d",fCenBin) ;
 765       FillHistogram(key,ph1->Pt(),-1.) ;
 766       if(ph1->IsPhoton()){
 767         snprintf(key,55,"hPhotDispwou_cen%d",fCenBin) ;
 768         FillHistogram(key,ph1->Pt(),-1.) ;
 769         snprintf(key,55,"hNegPhotDispwou_cen%d",fCenBin) ;
 770         FillHistogram(key,ph1->Pt(),-1.) ;
 771       }
 772       if(ph1->IsCPVOK()){
 773         snprintf(key,55,"hPhotBoth_cen%d",fCenBin) ;
 774         FillHistogram(key,ph1->Pt(),-1.) ;
 775         snprintf(key,55,"hNegPhotBoth_cen%d",fCenBin) ;
 776         FillHistogram(key,ph1->Pt(),-1.) ;
 777         snprintf(key,55,"hPhotBothcore_cen%d",fCenBin) ;
 778         FillHistogram(key,ph1->GetMomV2()->Pt(),-1.) ;
 779         snprintf(key,55,"hNegPhotBothcore_cen%d",fCenBin) ;
 780         FillHistogram(key,ph1->GetMomV2()->Pt(),-1.) ;
 781       }
 782     } // end of loop i2
 783   } // end of loop i1
 784
 785   for (Int_t i1=0; i1<inPHOSemb; i1++) {
 786     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
 787     if(!ph1->IsTagged())
 788       continue ;
 789     snprintf(key,55,"hPhotAll_cen%d",fCenBin) ;
 790     FillHistogram(key,ph1->Pt(),1.) ;
 791     snprintf(key,55,"hPhotAllcore_cen%d",fCenBin) ;
 792     FillHistogram(key,ph1->GetMomV2()->Pt(),1.) ;
 793     if(ph1->IsPhoton()){
 794       snprintf(key,55,"hPhotAllwou_cen%d",fCenBin) ;
 795       FillHistogram(key,ph1->Pt(),1.) ;
 796     }
 797     if(ph1->IsCPVOK() ){
 798       snprintf(key,55,"hPhotCPV_cen%d",fCenBin) ;
 799       FillHistogram(key,ph1->Pt(),1.) ;
 800       snprintf(key,55,"hPhotCPVcore_cen%d",fCenBin) ;
 801       FillHistogram(key,ph1->GetMomV2()->Pt(),1.) ;
 802     }
 803     if(ph1->IsCPV2OK() ){
 804       snprintf(key,55,"hPhotCPV2_cen%d",fCenBin) ;
 805       FillHistogram(key,ph1->Pt(),1.) ;
 806     }
 807     if(ph1->IsDisp2OK()){
 808       snprintf(key,55,"hPhotDisp2_cen%d",fCenBin) ;
 809       FillHistogram(key,ph1->Pt(),1.) ;
 810     }
 811     if(ph1->IsDispOK()){
 812       snprintf(key,55,"hPhotDisp_cen%d",fCenBin) ;
 813       FillHistogram(key,ph1->Pt(),1.) ;
 814       if(ph1->IsPhoton()){
 815         snprintf(key,55,"hPhotDispwou_cen%d",fCenBin) ;
 816         FillHistogram(key,ph1->Pt(),1.) ;
 817       }
 818       if(ph1->IsCPVOK()){
 819         snprintf(key,55,"hPhotBoth_cen%d",fCenBin) ;
 820         FillHistogram(key,ph1->Pt(),1.) ;
 821         snprintf(key,55,"hPhotBothcore_cen%d",fCenBin) ;
 822         FillHistogram(key,ph1->GetMomV2()->Pt(),1.) ;
 823       }
 824     } // end of loop i2
 825   } // end of loop i1
 826
 827
 828
 829   // Fill Real disribution:
 830   // Disappeared clusters enter with negative contribution
 831   // In addition fill control histogam with Real before embedding
 832   for (Int_t i1=0; i1<inPHOSold-1; i1++) {
 833     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent1->At(i1) ;
 834     for (Int_t i2=i1+1; i2<inPHOSold; i2++) {
 835       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent1->At(i2) ;
 836
 837       TLorentzVector p12  = *ph1  + *ph2;
 838       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 839       Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
 840
 841       //Fill Controll histogram: Real before embedding
 842       snprintf(key,55,"hOldMassPtAll_cen%d",fCenBin) ;
 843       FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 844       snprintf(key,55,"hOldMassPtAllcore_cen%d",fCenBin) ;
 845       FillHistogram(key,pv12.M() ,pv12.Pt(),-1.) ;
 846       if(ph1->IsPhoton() && ph2->IsPhoton()){
 847         snprintf(key,55,"hOldMassPtAllwou_cen%d",fCenBin) ;
 848         FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 849       }
 850       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 851         snprintf(key,55,"hOldMassPtCPV_cen%d",fCenBin) ;
 852         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 853         snprintf(key,55,"hOldMassPtCPV_cen%d",fCenBin) ;
 854         FillHistogram(key,pv12.M(), pv12.Pt(),-1) ;
 855       }
 856       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
 857         snprintf(key,55,"hOldMassPtCPV2_cen%d",fCenBin) ;
 858         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 859       }
 860       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
 861         snprintf(key,55,"hOldMassPtDisp2_cen%d",fCenBin) ;
 862         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 863       }
 864       if(ph1->IsDispOK() && ph2->IsDispOK()){
 865         snprintf(key,55,"hOldMassPtDisp_cen%d",fCenBin) ;
 866         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 867         if(ph1->IsPhoton() && ph2->IsPhoton()){
 868           snprintf(key,55,"hOldMassPtDispwou_cen%d",fCenBin) ;
 869           FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 870         }
 871         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 872           snprintf(key,55,"hOldMassPtBoth_cen%d",fCenBin) ;
 873           FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 874           snprintf(key,55,"hOldMassPtBothcore_cen%d",fCenBin) ;
 875           FillHistogram(key,pv12.M() ,pv12.Pt(),-1) ;
 876         }
 877       }
 878
 879       //Now fill main histograms with negative contributions
 880       if(!(ph1->IsTagged() || ph2->IsTagged()) )
 881         continue ;
 882       snprintf(key,55,"hMassPtAll_cen%d",fCenBin) ;
 883       FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 884       snprintf(key,55,"hMassPtAllcore_cen%d",fCenBin) ;
 885       FillHistogram(key,pv12.M() ,pv12.Pt(),-1.) ;
 886       if(ph1->IsPhoton() && ph2->IsPhoton()){
 887         snprintf(key,55,"hMassPtAllwou_cen%d",fCenBin) ;
 888         FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 889       }
 890       if(a<0.9){
 891         snprintf(key,55,"hMassPtAll_a09_cen%d",fCenBin) ;
 892         FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 893         if(a<0.8){
 894           snprintf(key,55,"hMassPtAll_a08_cen%d",fCenBin) ;
 895           FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 896           if(a<0.7){
 897             snprintf(key,55,"hMassPtAll_a07_cen%d",fCenBin) ;
 898             FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 899           }
 900         }
 901       }
 902       snprintf(key,55,"hNegMassPtAll_cen%d",fCenBin) ;
 903       FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 904       snprintf(key,55,"hNegMassPtAllcore_cen%d",fCenBin) ;
 905       FillHistogram(key,pv12.M() ,pv12.Pt(),-1.) ;
 906       if(ph1->IsPhoton() && ph2->IsPhoton()){
 907         snprintf(key,55,"hNegMassPtAllwou_cen%d",fCenBin) ;
 908         FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 909       }
 910
 911       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 912         snprintf(key,55,"hMassPtCPV_cen%d",fCenBin) ;
 913         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 914         snprintf(key,55,"hMassPtCPVcore_cen%d",fCenBin) ;
 915         FillHistogram(key,pv12.M() ,pv12.Pt(),-1) ;
 916         if(a<0.9){
 917           snprintf(key,55,"hMassPtCPV_a09_cen%d",fCenBin) ;
 918           FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 919           if(a<0.8){
 920             snprintf(key,55,"hMassPtCPV_a08_cen%d",fCenBin) ;
 921             FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 922             if(a<0.7){
 923               snprintf(key,55,"hMassPtCPV_a07_cen%d",fCenBin) ;
 924               FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 925             }
 926           }
 927         }
 928         snprintf(key,55,"hNegMassPtCPV_cen%d",fCenBin) ;
 929         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 930         snprintf(key,55,"hNegMassPtCPVcore_cen%d",fCenBin) ;
 931         FillHistogram(key,pv12.M() ,pv12.Pt(),-1) ;
 932       }
 933       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
 934         snprintf(key,55,"hMassPtCPV2_cen%d",fCenBin) ;
 935         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 936         if(a<0.9){
 937           snprintf(key,55,"hMassPtCPV2_a09_cen%d",fCenBin) ;
 938           FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 939           if(a<0.8){
 940             snprintf(key,55,"hMassPtCPV2_a08_cen%d",fCenBin) ;
 941             FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 942             if(a<0.7){
 943               snprintf(key,55,"hMassPtCPV2_a07_cen%d",fCenBin) ;
 944               FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 945             }
 946           }
 947         }
 948         snprintf(key,55,"hNegMassPtCPV2_cen%d",fCenBin) ;
 949         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 950       }
 951
 952       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
 953         snprintf(key,55,"hMassPtDisp2_cen%d",fCenBin) ;
 954         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 955       }
 956       if(ph1->IsDispOK() && ph2->IsDispOK()){
 957         snprintf(key,55,"hMassPtDisp_cen%d",fCenBin) ;
 958         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 959         if(ph1->IsPhoton() && ph2->IsPhoton()){
 960           snprintf(key,55,"hMassPtDispwou_cen%d",fCenBin) ;
 961           FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 962         }
 963         if(a<0.9){
 964           snprintf(key,55,"hMassPtDisp_a09_cen%d",fCenBin) ;
 965           FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 966           if(a<0.8){
 967             snprintf(key,55,"hMassPtDisp_a08_cen%d",fCenBin) ;
 968             FillHistogram(key,p12.M() ,p12.Pt()) ;
 969             if(a<0.7){
 970               snprintf(key,55,"hMassPtDisp_a07_cen%d",fCenBin) ;
 971               FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 972             }
 973           }
 974         }
 975         snprintf(key,55,"hNegMassPtDisp_cen%d",fCenBin) ;
 976         FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 977         if(ph1->IsPhoton() && ph2->IsPhoton()){
 978           snprintf(key,55,"hNegMassPtDispwou_cen%d",fCenBin) ;
 979           FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 980         }
 981         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 982           snprintf(key,55,"hMassPtBoth_cen%d",fCenBin) ;
 983           FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
 984           snprintf(key,55,"hMassPtBothcore_cen%d",fCenBin) ;
 985           FillHistogram(key,pv12.M() ,pv12.Pt(),-1) ;
 986           if(a<0.9){
 987             snprintf(key,55,"hMassPtBoth_a09_cen%d",fCenBin) ;
 988             FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 989             if(a<0.8){
 990               snprintf(key,55,"hMassPtBoth_a08_cen%d",fCenBin) ;
 991               FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 992               if(a<0.7){
 993                 snprintf(key,55,"hMassPtBoth_a07_cen%d",fCenBin) ;
 994                 FillHistogram(key,p12.M() ,p12.Pt(),-1.) ;
 995              }
 996             }
 997           }
 998           snprintf(key,55,"hNegMassPtBoth_cen%d",fCenBin) ;
 999           FillHistogram(key,p12.M() ,p12.Pt(),-1) ;
1000           snprintf(key,55,"hNegMassPtBothcore_cen%d",fCenBin) ;
1001           FillHistogram(key,pv12.M() ,pv12.Pt(),-1) ;
1002         }
1003       }
1004     } // end of loop i2
1005   } // end of loop i1
1006
1007
1008   // Further fill Real disribution
1009   // now with positive contribution from new clusters
1010   // ass well fill controll histogram
1011   for (Int_t i1=0; i1<inPHOSemb-1; i1++) {
1012     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
1013     for (Int_t i2=i1+1; i2<inPHOSemb; i2++) {
1014       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent2->At(i2) ;
1015
1016       TLorentzVector p12  = *ph1  + *ph2;
1017       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1018       Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
1019
1020       // Controll histogram: Real after embedding
1021       snprintf(key,55,"hNewMassPtAll_cen%d",fCenBin) ;
1022       FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
1023       snprintf(key,55,"hNewMassPtAllcore_cen%d",fCenBin) ;
1024       FillHistogram(key,pv12.M() ,pv12.Pt(),1.) ;
1025       if(ph1->IsPhoton() && ph2->IsPhoton()){
1026         snprintf(key,55,"hNewMassPtAllwou_cen%d",fCenBin) ;
1027         FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
1028       }
1029       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1030         snprintf(key,55,"hNewMassPtCPV_cen%d",fCenBin) ;
1031         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
1032         snprintf(key,55,"hNewMassPtCPVcore_cen%d",fCenBin) ;
1033         FillHistogram(key,pv12.M() ,pv12.Pt(),1) ;
1034       }
1035       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1036         snprintf(key,55,"hNewMassPtCPV2_cen%d",fCenBin) ;
1037         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
1038       }
1039       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1040         snprintf(key,55,"hNewMassPtDisp2_cen%d",fCenBin) ;
1041         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
1042       }
1043       if(ph1->IsDispOK() && ph2->IsDispOK()){
1044         snprintf(key,55,"hNewMassPtDisp_cen%d",fCenBin) ;
1045         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
1046         if(ph1->IsPhoton() && ph2->IsPhoton()){
1047           snprintf(key,55,"hNewMassPtDispwou_cen%d",fCenBin) ;
1048           FillHistogram(key,p12.M() ,p12.Pt(),1) ;
1049         }
1050         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1051           snprintf(key,55,"hNewMassPtBoth_cen%d",fCenBin) ;
1052           FillHistogram(key,p12.M() ,p12.Pt(),1) ;
1053           snprintf(key,55,"hNewMassPtBothcore_cen%d",fCenBin) ;
1054           FillHistogram(key,pv12.M() ,pv12.Pt(),1) ;
1055         }
1056       }
1057
1058       //Now fill main histogamm
1059       //new clusters with positive contribution
1060       if(!(ph1->IsTagged() || ph2->IsTagged()) )
1061         continue ;
1062       snprintf(key,55,"hMassPtAll_cen%d",fCenBin) ;
1063       FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
1064       snprintf(key,55,"hMassPtAllcore_cen%d",fCenBin) ;
1065       FillHistogram(key,pv12.M() ,pv12.Pt(),1.) ;
1066       if(ph1->IsPhoton() && ph2->IsPhoton()){
1067         snprintf(key,55,"hMassPtAllwou_cen%d",fCenBin) ;
1068         FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
1069       }
1070       if(a<0.9){
1071         snprintf(key,55,"hMassPtAll_a09_cen%d",fCenBin) ;
1072         FillHistogram(key,p12.M() ,p12.Pt()) ;
1073         if(a<0.8){
1074           snprintf(key,55,"hMassPtAll_a08_cen%d",fCenBin) ;
1075           FillHistogram(key,p12.M() ,p12.Pt()) ;
1076           if(a<0.7){
1077             snprintf(key,55,"hMassPtAll_a07_cen%d",fCenBin) ;
1078             FillHistogram(key,p12.M() ,p12.Pt()) ;
1079           }
1080         }
1081       }
1082       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1083         snprintf(key,55,"hMassPtCPV_cen%d",fCenBin) ;
1084         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
1085         snprintf(key,55,"hMassPtCPVcore_cen%d",fCenBin) ;
1086         FillHistogram(key,pv12.M() ,pv12.Pt(),1) ;
1087         if(a<0.9){
1088           snprintf(key,55,"hMassPtCPV_a09_cen%d",fCenBin) ;
1089           FillHistogram(key,p12.M() ,p12.Pt()) ;
1090           if(a<0.8){
1091             snprintf(key,55,"hMassPtCPV_a08_cen%d",fCenBin) ;
1092             FillHistogram(key,p12.M() ,p12.Pt()) ;
1093             if(a<0.7){
1094               snprintf(key,55,"hMassPtCPV_a07_cen%d",fCenBin) ;
1095               FillHistogram(key,p12.M() ,p12.Pt()) ;
1096             }
1097           }
1098         }
1099       }
1100       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1101         snprintf(key,55,"hMassPtCPV2_cen%d",fCenBin) ;
1102         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
1103          if(a<0.9){
1104           snprintf(key,55,"hMassPtCPV2_a09_cen%d",fCenBin) ;
1105           FillHistogram(key,p12.M() ,p12.Pt()) ;
1106           if(a<0.8){
1107             snprintf(key,55,"hMassPtCPV2_a08_cen%d",fCenBin) ;
1108             FillHistogram(key,p12.M() ,p12.Pt()) ;
1109             if(a<0.7){
1110               snprintf(key,55,"hMassPtCPV2_a07_cen%d",fCenBin) ;
1111               FillHistogram(key,p12.M() ,p12.Pt()) ;
1112             }
1113           }
1114         }
1115       }
1116       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1117         snprintf(key,55,"hMassPtDisp2_cen%d",fCenBin) ;
1118         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
1119       }
1120       if(ph1->IsDispOK() && ph2->IsDispOK()){
1121         snprintf(key,55,"hMassPtDisp_cen%d",fCenBin) ;
1122         FillHistogram(key,p12.M() ,p12.Pt(),1) ;
1123         if(ph1->IsPhoton() && ph2->IsPhoton()){
1124           snprintf(key,55,"hMassPtDispwou_cen%d",fCenBin) ;
1125           FillHistogram(key,p12.M() ,p12.Pt(),1) ;
1126         }
1127         if(a<0.9){
1128           snprintf(key,55,"hMassPtDisp_a09_cen%d",fCenBin) ;
1129           FillHistogram(key,p12.M() ,p12.Pt()) ;
1130           if(a<0.8){
1131             snprintf(key,55,"hMassPtDisp_a08_cen%d",fCenBin) ;
1132             FillHistogram(key,p12.M() ,p12.Pt()) ;
1133             if(a<0.7){
1134               snprintf(key,55,"hMassPtDisp_a07_cen%d",fCenBin) ;
1135               FillHistogram(key,p12.M() ,p12.Pt()) ;
1136             }
1137           }
1138         }
1139
1140         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1141           snprintf(key,55,"hMassPtBoth_cen%d",fCenBin) ;
1142           FillHistogram(key,p12.M() ,p12.Pt(),1) ;
1143           snprintf(key,55,"hMassPtBothcore_cen%d",fCenBin) ;
1144           FillHistogram(key,pv12.M() ,pv12.Pt(),1) ;
1145           if(a<0.9){
1146             snprintf(key,55,"hMassPtBoth_a09_cen%d",fCenBin) ;
1147             FillHistogram(key,p12.M() ,p12.Pt()) ;
1148             if(a<0.8){
1149               snprintf(key,55,"hMassPtBoth_a08_cen%d",fCenBin) ;
1150               FillHistogram(key,p12.M() ,p12.Pt()) ;
1151               if(a<0.7){
1152                 snprintf(key,55,"hMassPtBoth_a07_cen%d",fCenBin) ;
1153                 FillHistogram(key,p12.M() ,p12.Pt()) ;
1154               }
1155             }
1156           }
1157         }
1158       }
1159     } // end of loop i2
1160   } // end of loop i1
1161
1162
1163   //now mixed, does not really matter old or new list of clusters
1164   for (Int_t i1=0; i1<inPHOSemb; i1++) {
1165     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent2->At(i1) ;
1166     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
1167       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
1168       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
1169         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
1170         TLorentzVector p12  = *ph1  + *ph2;
1171         TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1172         Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
1173
1174         snprintf(key,55,"hMiMassPtAll_cen%d",fCenBin) ;
1175         FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
1176         snprintf(key,55,"hMiMassPtAllcore_cen%d",fCenBin) ;
1177         FillHistogram(key,pv12.M() ,pv12.Pt(),1.) ;
1178         if(ph1->IsPhoton() && ph2->IsPhoton()){
1179           snprintf(key,55,"hMiMassPtAllwou_cen%d",fCenBin) ;
1180           FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
1181         }
1182         if(a<0.9){
1183           snprintf(key,55,"hMiMassPtAll_a09_cen%d",fCenBin) ;
1184           FillHistogram(key,p12.M() ,p12.Pt()) ;
1185           if(a<0.8){
1186             snprintf(key,55,"hMiMassPtAll_a08_cen%d",fCenBin) ;
1187             FillHistogram(key,p12.M() ,p12.Pt()) ;
1188             if(a<0.7){
1189               snprintf(key,55,"hMiMassPtAll_a07_cen%d",fCenBin) ;
1190               FillHistogram(key,p12.M() ,p12.Pt()) ;
1191             }
1192           }
1193         }
1194         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1195           snprintf(key,55,"hMiMassPtCPV_cen%d",fCenBin) ;
1196           FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
1197           snprintf(key,55,"hMiMassPtCPVcore_cen%d",fCenBin) ;
1198           FillHistogram(key,pv12.M() ,pv12.Pt(),1.) ;
1199           if(a<0.9){
1200             snprintf(key,55,"hMiMassPtCPV_a09_cen%d",fCenBin) ;
1201             FillHistogram(key,p12.M() ,p12.Pt()) ;
1202             if(a<0.8){
1203               snprintf(key,55,"hMiMassPtCPV_a08_cen%d",fCenBin) ;
1204               FillHistogram(key,p12.M() ,p12.Pt()) ;
1205               if(a<0.7){
1206                 snprintf(key,55,"hMiMassPtCPV_a07_cen%d",fCenBin) ;
1207                 FillHistogram(key,p12.M() ,p12.Pt()) ;
1208               }
1209             }
1210           }
1211         }
1212         if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1213           snprintf(key,55,"hMiMassPtCPV2_cen%d",fCenBin) ;
1214           FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
1215           if(a<0.9){
1216             snprintf(key,55,"hMiMassPtCPV2_a09_cen%d",fCenBin) ;
1217             FillHistogram(key,p12.M() ,p12.Pt()) ;
1218             if(a<0.8){
1219               snprintf(key,55,"hMiMassPtCPV2_a08_cen%d",fCenBin) ;
1220               FillHistogram(key,p12.M() ,p12.Pt()) ;
1221               if(a<0.7){
1222                 snprintf(key,55,"hMiMassPtCPV2_a07_cen%d",fCenBin) ;
1223                 FillHistogram(key,p12.M() ,p12.Pt()) ;
1224               }
1225             }
1226           }
1227         }
1228         if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1229           snprintf(key,55,"hMiMassPtDisp2_cen%d",fCenBin) ;
1230           FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
1231         }
1232         if(ph1->IsDispOK() && ph2->IsDispOK()){
1233           snprintf(key,55,"hMiMassPtDisp_cen%d",fCenBin) ;
1234           FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
1235           if(ph1->IsPhoton() && ph2->IsPhoton()){
1236             snprintf(key,55,"hMiMassPtDispwou_cen%d",fCenBin) ;
1237             FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
1238           }
1239           if(a<0.9){
1240             snprintf(key,55,"hMiMassPtDisp_a09_cen%d",fCenBin) ;
1241             FillHistogram(key,p12.M() ,p12.Pt()) ;
1242             if(a<0.8){
1243               snprintf(key,55,"hMiMassPtDisp_a08_cen%d",fCenBin) ;
1244               FillHistogram(key,p12.M() ,p12.Pt()) ;
1245               if(a<0.7){
1246                 snprintf(key,55,"hMiMassPtDisp_a07_cen%d",fCenBin) ;
1247                 FillHistogram(key,p12.M() ,p12.Pt()) ;
1248               }
1249             }
1250           }
1251
1252           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1253             snprintf(key,55,"hMiMassPtBoth_cen%d",fCenBin) ;
1254             FillHistogram(key,p12.M() ,p12.Pt(),1.) ;
1255             snprintf(key,55,"hMiMassPtBothcore_cen%d",fCenBin) ;
1256             FillHistogram(key,pv12.M() ,pv12.Pt(),1.) ;
1257             if(a<0.9){
1258               snprintf(key,55,"hMiMassPtBoth_a09_cen%d",fCenBin) ;
1259               FillHistogram(key,p12.M() ,p12.Pt()) ;
1260               if(a<0.8){
1261                 snprintf(key,55,"hMiMassPtBoth_a08_cen%d",fCenBin) ;
1262                 FillHistogram(key,p12.M() ,p12.Pt()) ;
1263                 if(a<0.7){
1264                   snprintf(key,55,"hMiMassPtBoth_a07_cen%d",fCenBin) ;
1265                   FillHistogram(key,p12.M() ,p12.Pt()) ;
1266                }
1267               }
1268             }
1269           }
1270         }
1271       } // end of loop i2
1272     }
1273   } // end of loop i1
1274
1275
1276   //Now we either add current events to stack or remove
1277   //If no photons in current event - no need to add it to mixed
1278   if(fPHOSEvent2->GetEntriesFast()>0){
1279     prevPHOS->AddFirst(fPHOSEvent2) ;
1280     fPHOSEvent2=0;
1281     delete fPHOSEvent1;
1282     fPHOSEvent1=0;
1283     if(prevPHOS->GetSize()>100){//Remove redundant events
1284       TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
1285       prevPHOS->RemoveLast() ;
1286       delete tmp ;
1287     }
1288   }
1289   // Post output data.
1290   PostData(1, fOutputContainer);
1291   fEventCounter++;
1292 }
1293 //___________________________________________________________________________
1294 Bool_t AliAnalysisTaskPi0DiffEfficiency::IsSameCluster(AliAODCaloCluster * c1, AliAODCaloCluster * c2)const{
1295  //Compare clusters before and after embedding
1296  //clusters are the same if
1297  // - Energy changed less than 0.1%  (numerical accuracy in reconstruction)
1298  // - lists of digits are the same
1299
1300  if(c1->GetNCells() != c2->GetNCells())
1301    return kFALSE ;
1302
1303  if(TMath::Abs(c1->E()-c2->E())>0.01*c1->E())
1304    return kFALSE ;
1305
1306  UShort_t *list1 = c1->GetCellsAbsId() ;
1307  UShort_t *list2 = c2->GetCellsAbsId() ;
1308  for(Int_t i=0; i<c1->GetNCells(); i++){
1309   if(list1[i] != list2[i])
1310     return kFALSE ;
1311  }
1312  return kTRUE ;
1313
1314 }
1315
1316
1317