PWGGA/PHOSTasks/PHOS_embedding/AliAnalysisTaskPi0Efficiency.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 #include "TChain.h"
  17 #include "TTree.h"
  18 #include "TObjArray.h"
  19 #include "THashList.h"
  20 #include "TF1.h"
  21 #include "TFile.h"
  22 #include "TH1F.h"
  23 #include "TH2F.h"
  24 #include "TH2I.h"
  25 #include "TH3F.h"
  26 #include "TParticle.h"
  27 #include "TCanvas.h"
  28 #include "TStyle.h"
  29 #include "TRandom.h"
  30
  31 #include "AliAODMCParticle.h"
  32 #include "AliAnalysisManager.h"
  33 #include "AliMCEventHandler.h"
  34 #include "AliMCEvent.h"
  35 #include "AliStack.h"
  36 #include "AliAnalysisTaskSE.h"
  37 #include "AliAnalysisTaskPi0Efficiency.h"
  38 #include "AliCaloPhoton.h"
  39 #include "AliPHOSGeometry.h"
  40 #include "AliPHOSAodCluster.h"
  41 #include "AliPHOSCalibData.h"
  42 #include "AliAODEvent.h"
  43 #include "AliAODCaloCluster.h"
  44 #include "AliAODVertex.h"
  45 #include "AliESDtrackCuts.h"
  46 #include "AliLog.h"
  47 #include "AliPID.h"
  48 #include "AliCDBManager.h"
  49 #include "AliCentrality.h"
  50
  51 // Analysis task to fill histograms with PHOS ESD clusters and cells
  52 // Authors: Yuri Kharlov
  53 // Date   : 28.05.2009
  54
  55 ClassImp(AliAnalysisTaskPi0Efficiency)
  56
  57 //________________________________________________________________________
  58 AliAnalysisTaskPi0Efficiency::AliAnalysisTaskPi0Efficiency(const char *name)
  59 : AliAnalysisTaskSE(name),
  60   fStack(0),
  61   fOutputContainer(0),
  62   fPHOSEvent(0),
  63   fPHOSCalibData(0),
  64   fNonLinCorr(0),
  65   fRPfull(0),
  66   fRPA(0),
  67   fRPC(0),
  68   fRPFar(0),
  69   fRPAFar(0),
  70   fRPCFar(0),
  71   fCentrality(0),
  72   fCenBin(0),
  73   fPHOSGeo(0),
  74   fEventCounter(0)
  75 {
  76   // Constructor
  77   for(Int_t i=0;i<1;i++){
  78     for(Int_t j=0;j<10;j++)
  79       for(Int_t k=0;k<11;k++)
  80         fPHOSEvents[i][j][k]=0 ;
  81   }
  82
  83   // Output slots #0 write into a TH1 container
  84   DefineOutput(1,THashList::Class());
  85
  86   // Set bad channel map
  87   char key[55] ;
  88   for(Int_t i=0; i<6; i++){
  89     snprintf(key,55,"PHOS_BadMap_mod%d",i) ;
  90     fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
  91   }
  92   // Initialize the PHOS geometry
  93   fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
  94
  95   fPHOSCalibData = new AliPHOSCalibData();
  96   for(Int_t module=1; module<=5; module++) {
  97     for(Int_t column=1; column<=56; column++) {
  98       for(Int_t row=1; row<=64; row++) {
  99         fPHOSCalibData->SetADCchannelEmc(module,column,row,1.);
 100       }
 101     }
 102   }
 103
 104
 105 }
 106
 107 //________________________________________________________________________
 108 void AliAnalysisTaskPi0Efficiency::UserCreateOutputObjects()
 109 {
 110   // Create histograms
 111   // Called once
 112
 113   // ESD histograms
 114   if(fOutputContainer != NULL){
 115     delete fOutputContainer;
 116   }
 117   fOutputContainer = new THashList();
 118   fOutputContainer->SetOwner(kTRUE);
 119
 120   //Event selection
 121   fOutputContainer->Add(new TH1F("hSelEvents","Event celection", 10,0.,10.)) ;
 122
 123   //vertex distribution
 124   fOutputContainer->Add(new TH1F("hZvertex","Z vertex position", 50,-25.,25.)) ;
 125
 126   //Centrality
 127   fOutputContainer->Add(new TH1F("hCentrality","Event centrality", 100,0.,100.)) ;
 128
 129   //QA histograms
 130   fOutputContainer->Add(new TH2F("hCluM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 131   fOutputContainer->Add(new TH2F("hCluM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 132   fOutputContainer->Add(new TH2F("hCluM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 133
 134   Int_t nM       = 500;
 135   Double_t mMin  = 0.0;
 136   Double_t mMax  = 1.0;
 137   Int_t nPt      = 200;
 138   Double_t ptMin = 0;
 139   Double_t ptMax = 20;
 140
 141   char key[55] ;
 142   for(Int_t cent=0; cent<6; cent++){
 143     //Single photon
 144     snprintf(key,55,"hPhotAll_cen%d",cent) ;
 145     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 146     snprintf(key,55,"hPhotAllwou_cen%d",cent) ;
 147     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 148     snprintf(key,55,"hPhotAllcore_cen%d",cent) ;
 149     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 150     snprintf(key,55,"hPhotCPV_cen%d",cent) ;
 151     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 152     snprintf(key,55,"hPhotCPVcore_cen%d",cent) ;
 153     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 154     snprintf(key,55,"hPhotCPV2_cen%d",cent) ;
 155     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 156     snprintf(key,55,"hPhotDisp_cen%d",cent) ;
 157     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 158     snprintf(key,55,"hPhotDispwou_cen%d",cent) ;
 159     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 160     snprintf(key,55,"hPhotDisp2_cen%d",cent) ;
 161     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 162     snprintf(key,55,"hPhotBoth_cen%d",cent) ;
 163     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 164     snprintf(key,55,"hPhotBothcore_cen%d",cent) ;
 165     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 166
 167     snprintf(key,55,"hMassPtAll_cen%d",cent) ;
 168     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 169     snprintf(key,55,"hMassPtAllwou_cen%d",cent) ;
 170     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 171     snprintf(key,55,"hMassPtAllcore_cen%d",cent) ;
 172     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 173     snprintf(key,55,"hMassPtCPV_cen%d",cent) ;
 174     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 175     snprintf(key,55,"hMassPtCPVcore_cen%d",cent) ;
 176     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 177     snprintf(key,55,"hMassPtCPV2_cen%d",cent) ;
 178     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 179     snprintf(key,55,"hMassPtDisp_cen%d",cent) ;
 180     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 181     snprintf(key,55,"hMassPtDispwou_cen%d",cent) ;
 182     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 183     snprintf(key,55,"hMassPtDisp2_cen%d",cent) ;
 184     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 185     snprintf(key,55,"hMassPtBoth_cen%d",cent) ;
 186     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 187     snprintf(key,55,"hMassPtBothcore_cen%d",cent) ;
 188     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 189
 190     snprintf(key,55,"hMassPtAll_a07_cen%d",cent) ;
 191     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 192     snprintf(key,55,"hMassPtCPV_a07_cen%d",cent) ;
 193     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 194     snprintf(key,55,"hMassPtCPV2_a07_cen%d",cent) ;
 195     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 196     snprintf(key,55,"hMassPtDisp_a07_cen%d",cent) ;
 197     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 198     snprintf(key,55,"hMassPtBoth_a07_cen%d",cent) ;
 199     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 200
 201     snprintf(key,55,"hMassPtAll_a08_cen%d",cent) ;
 202     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 203     snprintf(key,55,"hMassPtCPV_a08_cen%d",cent) ;
 204     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 205     snprintf(key,55,"hMassPtCPV2_a08_cen%d",cent) ;
 206     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 207     snprintf(key,55,"hMassPtDisp_a08_cen%d",cent) ;
 208     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 209     snprintf(key,55,"hMassPtBoth_a08_cen%d",cent) ;
 210     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 211
 212     snprintf(key,55,"hMassPtAll_a09_cen%d",cent) ;
 213     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 214     snprintf(key,55,"hMassPtCPV_a09_cen%d",cent) ;
 215     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 216     snprintf(key,55,"hMassPtCPV2_a09_cen%d",cent) ;
 217     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 218     snprintf(key,55,"hMassPtDisp_a09_cen%d",cent) ;
 219     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 220     snprintf(key,55,"hMassPtBoth_a09_cen%d",cent) ;
 221     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 222
 223
 224     //Mixed
 225     snprintf(key,55,"hMiMassPtAll_cen%d",cent) ;
 226     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 227     snprintf(key,55,"hMiMassPtAllwou_cen%d",cent) ;
 228     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 229     snprintf(key,55,"hMiMassPtAllcore_cen%d",cent) ;
 230     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 231     snprintf(key,55,"hMiMassPtCPV_cen%d",cent) ;
 232     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 233     snprintf(key,55,"hMiMassPtCPVcore_cen%d",cent) ;
 234     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 235     snprintf(key,55,"hMiMassPtCPV2_cen%d",cent) ;
 236     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 237     snprintf(key,55,"hMiMassPtDisp_cen%d",cent) ;
 238     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 239     snprintf(key,55,"hMiMassPtDispwou_cen%d",cent) ;
 240     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 241     snprintf(key,55,"hMiMassPtDisp2_cen%d",cent) ;
 242     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 243     snprintf(key,55,"hMiMassPtBoth_cen%d",cent) ;
 244     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 245     snprintf(key,55,"hMiMassPtBothcore_cen%d",cent) ;
 246     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 247
 248     snprintf(key,55,"hMiMassPtAll_a07_cen%d",cent) ;
 249     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 250     snprintf(key,55,"hMiMassPtCPV_a07_cen%d",cent) ;
 251     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 252     snprintf(key,55,"hMiMassPtCPV2_a07_cen%d",cent) ;
 253     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 254     snprintf(key,55,"hMiMassPtDisp_a07_cen%d",cent) ;
 255     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 256     snprintf(key,55,"hMiMassPtBoth_a07_cen%d",cent) ;
 257     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 258
 259     snprintf(key,55,"hMiMassPtAll_a08_cen%d",cent) ;
 260     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 261     snprintf(key,55,"hMiMassPtCPV_a08_cen%d",cent) ;
 262     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 263     snprintf(key,55,"hMiMassPtCPV2_a08_cen%d",cent) ;
 264     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 265     snprintf(key,55,"hMiMassPtDisp_a08_cen%d",cent) ;
 266     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 267     snprintf(key,55,"hMiMassPtBoth_a08_cen%d",cent) ;
 268     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 269
 270      snprintf(key,55,"hMiMassPtAll_a09_cen%d",cent) ;
 271     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 272     snprintf(key,55,"hMiMassPtCPV_a09_cen%d",cent) ;
 273     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 274     snprintf(key,55,"hMiMassPtCPV2_a09_cen%d",cent) ;
 275     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 276     snprintf(key,55,"hMiMassPtDisp_a09_cen%d",cent) ;
 277     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 278     snprintf(key,55,"hMiMassPtBoth_a09_cen%d",cent) ;
 279     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 280
 281     //MC
 282     snprintf(key,55,"hMCMassPtAll_cen%d",cent) ;
 283     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 284     snprintf(key,55,"hMCMassPtAllwou_cen%d",cent) ;
 285     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 286     snprintf(key,55,"hMCMassPtAllcore_cen%d",cent) ;
 287     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 288     snprintf(key,55,"hMCMassPtCPV_cen%d",cent) ;
 289     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 290     snprintf(key,55,"hMCMassPtCPVcore_cen%d",cent) ;
 291     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 292     snprintf(key,55,"hMCMassPtCPV2_cen%d",cent) ;
 293     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 294     snprintf(key,55,"hMCMassPtDisp_cen%d",cent) ;
 295     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 296     snprintf(key,55,"hMCMassPtDispwou_cen%d",cent) ;
 297     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 298     snprintf(key,55,"hMCMassPtDisp2_cen%d",cent) ;
 299     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 300     snprintf(key,55,"hMCMassPtBoth_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,"hMCMassPtBothcore_cen%d",cent) ;
 303     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 304
 305
 306     snprintf(key,55,"hMCMassPtAll_a07_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,"hMCMassPtCPV_a07_cen%d",cent) ;
 309     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 310     snprintf(key,55,"hMCMassPtCPV2_a07_cen%d",cent) ;
 311     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 312     snprintf(key,55,"hMCMassPtDisp_a07_cen%d",cent) ;
 313     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 314     snprintf(key,55,"hMCMassPtBoth_a07_cen%d",cent) ;
 315     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 316
 317     snprintf(key,55,"hMCMassPtAll_a08_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,"hMCMassPtCPV_a08_cen%d",cent) ;
 320     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 321     snprintf(key,55,"hMCMassPtCPV2_a08_cen%d",cent) ;
 322     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 323     snprintf(key,55,"hMCMassPtDisp_a08_cen%d",cent) ;
 324     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 325     snprintf(key,55,"hMCMassPtBoth_a08_cen%d",cent) ;
 326     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 327
 328     snprintf(key,55,"hMCMassPtAll_a09_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,"hMCMassPtCPV_a09_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,"hMCMassPtCPV2_a09_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,"hMCMassPtDisp_a09_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,"hMCMassPtBoth_a09_cen%d",cent) ;
 337     fOutputContainer->Add(new TH2F(key,"(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 338
 339     //Single photon
 340     snprintf(key,55,"hMCPhotAll_cen%d",cent) ;
 341     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 342     snprintf(key,55,"hMCPhotAllwou_cen%d",cent) ;
 343     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 344     snprintf(key,55,"hMCPhotAllcore_cen%d",cent) ;
 345     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 346     snprintf(key,55,"hMCPhotCPV_cen%d",cent) ;
 347     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 348     snprintf(key,55,"hMCPhotCPVcore_cen%d",cent) ;
 349     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 350     snprintf(key,55,"hMCPhotCPV2_cen%d",cent) ;
 351     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 352     snprintf(key,55,"hMCPhotDisp_cen%d",cent) ;
 353     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 354     snprintf(key,55,"hMCPhotDispwou_cen%d",cent) ;
 355     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 356     snprintf(key,55,"hMCPhotDisp2_cen%d",cent) ;
 357     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 358     snprintf(key,55,"hMCPhotBoth_cen%d",cent) ;
 359     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 360     snprintf(key,55,"hMCPhotBothcore_cen%d",cent) ;
 361     fOutputContainer->Add(new TH1F(key,"dN/dpt" ,nPt,ptMin,ptMax));
 362
 363   }
 364
 365   fOutputContainer->Add(new TH2F("hMCPi0M11","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 366   fOutputContainer->Add(new TH2F("hMCPi0M22","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 367   fOutputContainer->Add(new TH2F("hMCPi0M33","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 368   fOutputContainer->Add(new TH2F("hMCPi0M12","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 369   fOutputContainer->Add(new TH2F("hMCPi0M13","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 370   fOutputContainer->Add(new TH2F("hMCPi0M23","(M,p_{T},d#phi)_{#gamma#gamma}" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 371
 372
 373   //MC
 374   for(Int_t cent=0; cent<6; cent++){
 375     snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
 376     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
 377     snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
 378     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
 379     snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
 380     fOutputContainer->Add(new TH1F("hMC_rap_eta","Rapidity eta",200,-1.,1.)) ;
 381     snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
 382     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
 383     snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
 384     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
 385     snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
 386     fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
 387     snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
 388     fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
 389     snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
 390     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
 391     snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
 392     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 393     snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
 394     fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
 395     snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
 396     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 397     snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
 398     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 399   }
 400
 401   PostData(1, fOutputContainer);
 402
 403 }
 404
 405 //________________________________________________________________________
 406 void AliAnalysisTaskPi0Efficiency::UserExec(Option_t *)
 407 {
 408   // Main loop, called for each event
 409   // Analyze ESD/AOD
 410
 411   FillHistogram("hSelEvents",0.5) ;
 412
 413   AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
 414   if (!event) {
 415     Printf("ERROR: Could not retrieve event");
 416     PostData(1, fOutputContainer);
 417     return;
 418   }
 419
 420   FillHistogram("hSelEvents",1.5) ;
 421   AliAODHeader *header = dynamic_cast<AliAODHeader*>(event->GetHeader()) ;
 422   if(!header) AliFatal("Not a standard AOD");
 423
 424   // Checks if we have a primary vertex
 425   // Get primary vertices form ESD
 426   const AliAODVertex *esdVertex5 = event->GetPrimaryVertex();
 427
 428  // don't rely on ESD vertex, assume (0,0,0)
 429   Double_t vtx0[3] ={0.,0.,0.};
 430
 431
 432   FillHistogram("hZvertex",esdVertex5->GetZ());
 433   if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
 434     PostData(1, fOutputContainer);
 435     return;
 436   }
 437   FillHistogram("hSelEvents",2.5) ;
 438
 439   //Vtx class z-bin
 440   //  Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
 441   //  if(zvtx<0)zvtx=0 ;
 442   //  if(zvtx>9)zvtx=9 ;
 443   Int_t zvtx=0 ;
 444
 445 //  fCentrality=header->GetCentralityP()->GetCentralityPercentile("V0M"); // returns the centrality percentile,
 446 //                                                          //a float from 0 to 100 (or to the trigger efficiency)
 447    fCentrality=header->GetZDCN2Energy() ;
 448
 449   if( fCentrality < 0. || fCentrality>80.){
 450     PostData(1, fOutputContainer);
 451     return;
 452   }
 453   FillHistogram("hSelEvents",3.5) ;
 454   Float_t bins[7]={0.,5.,10.,20.,40.,60.,80.} ;
 455   fCenBin=0 ;
 456   while(fCenBin<6 && fCentrality > bins[fCenBin+1])
 457     fCenBin++ ;
 458
 459
 460   //reaction plain
 461   fRPfull= header->GetZDCN1Energy() ;
 462   if(fRPfull==999){ //reaction plain was not defined
 463     PostData(1, fOutputContainer);
 464     return;
 465   }
 466
 467   FillHistogram("hSelEvents",4.5) ;
 468   //All event selections done
 469   FillHistogram("hCentrality",fCentrality) ;
 470   //Reaction plain is defined in the range (-pi/2;pi/2)
 471   //We have 10 bins
 472   Int_t irp=Int_t(10.*fRPfull/TMath::Pi());
 473   if(irp>9)irp=9 ;
 474
 475   if(!fPHOSEvents[zvtx][fCenBin][irp])
 476     fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
 477   TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
 478
 479   // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
 480   if(fEventCounter == 0) {
 481     for(Int_t mod=0; mod<5; mod++) {
 482       const TGeoHMatrix* m =header->GetPHOSMatrix(mod) ;
 483       fPHOSGeo->SetMisalMatrix(m,mod) ;
 484       Printf("PHOS geo matrix for module # %d is set: %p\n", mod,m);
 485     }
 486     fEventCounter++ ;
 487   }
 488
 489   ProcessMC() ;
 490
 491   if(fPHOSEvent)
 492     fPHOSEvent->Clear() ;
 493   else
 494     fPHOSEvent = new TClonesArray("AliCaloPhoton",200) ;
 495
 496
 497   char key[55] ;
 498   Int_t inPHOS=0 ;
 499   TVector3 vertex(vtx0);
 500   TClonesArray * clusters = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
 501   AliAODCaloCells * cells = (AliAODCaloCells*)event->FindListObject("EmbeddedPHOScells") ;
 502   Int_t multClust = clusters->GetEntriesFast();
 503   for (Int_t i=0; i<multClust; i++) {
 504     AliAODCaloCluster *clu = (AliAODCaloCluster*) clusters->At(i);
 505     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
 506
 507     Float_t  position[3];
 508     clu->GetPosition(position);
 509     TVector3 global(position) ;
 510     Int_t relId[4] ;
 511     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 512     Int_t mod  = relId[0] ;
 513     Int_t cellX = relId[2];
 514     Int_t cellZ = relId[3] ;
 515     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
 516       continue ;
 517     if(clu->GetNCells()<3)
 518       continue ;
 519     if(clu->GetM02()<0.2)
 520       continue ;
 521
 522     snprintf(key,55,"hCluM%d",mod) ;
 523     FillHistogram(key,cellX,cellZ,1.);
 524
 525     TLorentzVector pv1 ;
 526     clu->GetMomentum(pv1 ,vtx0);
 527
 528     if(inPHOS>=fPHOSEvent->GetSize()){
 529       fPHOSEvent->Expand(inPHOS+50) ;
 530     }
 531     new((*fPHOSEvent)[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
 532     AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
 533     ph->SetModule(mod) ;
 534     AliPHOSAodCluster cluPHOS1(*clu);
 535     cluPHOS1.Recalibrate(fPHOSCalibData,cells); // modify the cell energies
 536     Double_t ecore=CoreEnergy(&cluPHOS1) ;
 537     pv1*= ecore/pv1.E() ;
 538     ph->SetMomV2(&pv1) ;
 539     ph->SetNCells(clu->GetNCells());
 540     ph->SetDispBit(TestLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
 541     ph->SetDisp2Bit(TestLambda2(clu->E(),clu->GetM20(),clu->GetM02())) ;
 542     ph->SetCPVBit(clu->GetEmcCpvDistance()>2.) ;
 543     ph->SetCPV2Bit(clu->GetEmcCpvDistance()>4.) ;
 544     ph->SetPhoton(clu->GetNExMax()<2); // Remember, if it is unfolded
 545
 546     inPHOS++ ;
 547   }
 548   //Single photon
 549   for (Int_t i1=0; i1<inPHOS; i1++) {
 550     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 551     snprintf(key,55,"hPhotAll_cen%d",fCenBin) ;
 552     FillHistogram(key,ph1->Pt()) ;
 553     snprintf(key,55,"hPhotAllcore_cen%d",fCenBin) ;
 554     FillHistogram(key,ph1->GetMomV2()->Pt()) ;
 555     if(ph1->IsPhoton()){
 556       snprintf(key,55,"hPhotAllwou_cen%d",fCenBin) ;
 557       FillHistogram(key,ph1->Pt()) ;
 558     }
 559     if(ph1->IsCPVOK() ){
 560       snprintf(key,55,"hPhotCPV_cen%d",fCenBin) ;
 561       FillHistogram(key,ph1->Pt()) ;
 562       snprintf(key,55,"hPhotCPVcore_cen%d",fCenBin) ;
 563       FillHistogram(key,ph1->GetMomV2()->Pt()) ;
 564     }
 565     if(ph1->IsCPV2OK() ){
 566       snprintf(key,55,"hPhotCPV2_cen%d",fCenBin) ;
 567       FillHistogram(key,ph1->Pt()) ;
 568     }
 569     if(ph1->IsDisp2OK()){
 570       snprintf(key,55,"hPhotDisp2_cen%d",fCenBin) ;
 571       FillHistogram(key,ph1->Pt()) ;
 572     }
 573     if(ph1->IsDispOK()){
 574       snprintf(key,55,"hPhotDisp_cen%d",fCenBin) ;
 575       FillHistogram(key,ph1->Pt()) ;
 576       if(ph1->IsPhoton()){
 577         snprintf(key,55,"hPhotDispwou_cen%d",fCenBin) ;
 578         FillHistogram(key,ph1->Pt()) ;
 579       }
 580       if(ph1->IsCPVOK()){
 581         snprintf(key,55,"hPhotBoth_cen%d",fCenBin) ;
 582         FillHistogram(key,ph1->Pt()) ;
 583         snprintf(key,55,"hPhotBothcore_cen%d",fCenBin) ;
 584         FillHistogram(key,ph1->GetMomV2()->Pt()) ;
 585       }
 586     } // end of loop i2
 587   } // end of loop i1
 588
 589   // Fill Real disribution
 590   for (Int_t i1=0; i1<inPHOS-1; i1++) {
 591     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 592     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
 593       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent->At(i2) ;
 594       TLorentzVector p12  = *ph1  + *ph2;
 595       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 596       Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
 597
 598       snprintf(key,55,"hMassPtAll_cen%d",fCenBin) ;
 599       FillHistogram(key,p12.M() ,p12.Pt()) ;
 600       snprintf(key,55,"hMassPtAllcore_cen%d",fCenBin) ;
 601       FillHistogram(key,pv12.M(), pv12.Pt()) ;
 602       if(ph1->IsPhoton() && ph2->IsPhoton()){
 603         snprintf(key,55,"hMassPtAllwou_cen%d",fCenBin) ;
 604         FillHistogram(key,p12.M() ,p12.Pt()) ;
 605       }
 606       if(a<0.9){
 607         snprintf(key,55,"hMassPtAll_a09_cen%d",fCenBin) ;
 608         FillHistogram(key,p12.M() ,p12.Pt()) ;
 609         if(a<0.8){
 610           snprintf(key,55,"hMassPtAll_a08_cen%d",fCenBin) ;
 611           FillHistogram(key,p12.M() ,p12.Pt()) ;
 612           if(a<0.7){
 613             snprintf(key,55,"hMassPtAll_a07_cen%d",fCenBin) ;
 614             FillHistogram(key,p12.M() ,p12.Pt()) ;
 615           }
 616         }
 617       }
 618       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 619         snprintf(key,55,"hMassPtCPV_cen%d",fCenBin) ;
 620         FillHistogram(key,p12.M() ,p12.Pt()) ;
 621         snprintf(key,55,"hMassPtCPVcore_cen%d",fCenBin) ;
 622         FillHistogram(key,pv12.M(), pv12.Pt()) ;
 623         if(a<0.9){
 624           snprintf(key,55,"hMassPtCPV_a09_cen%d",fCenBin) ;
 625           FillHistogram(key,p12.M() ,p12.Pt()) ;
 626           if(a<0.8){
 627             snprintf(key,55,"hMassPtCPV_a08_cen%d",fCenBin) ;
 628             FillHistogram(key,p12.M() ,p12.Pt()) ;
 629             if(a<0.7){
 630               snprintf(key,55,"hMassPtCPV_a07_cen%d",fCenBin) ;
 631               FillHistogram(key,p12.M() ,p12.Pt()) ;
 632             }
 633           }
 634         }
 635       }
 636       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
 637         snprintf(key,55,"hMassPtCPV2_cen%d",fCenBin) ;
 638         FillHistogram(key,p12.M() ,p12.Pt()) ;
 639         if(a<0.9){
 640           snprintf(key,55,"hMassPtCPV2_a09_cen%d",fCenBin) ;
 641           FillHistogram(key,p12.M() ,p12.Pt()) ;
 642           if(a<0.8){
 643             snprintf(key,55,"hMassPtCPV2_a08_cen%d",fCenBin) ;
 644             FillHistogram(key,p12.M() ,p12.Pt()) ;
 645             if(a<0.7){
 646               snprintf(key,55,"hMassPtCPV2_a07_cen%d",fCenBin) ;
 647               FillHistogram(key,p12.M() ,p12.Pt()) ;
 648             }
 649           }
 650         }
 651       }
 652       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
 653         snprintf(key,55,"hMassPtDisp2_cen%d",fCenBin) ;
 654         FillHistogram(key,p12.M() ,p12.Pt()) ;
 655       }
 656       if(ph1->IsDispOK() && ph2->IsDispOK()){
 657         snprintf(key,55,"hMassPtDisp_cen%d",fCenBin) ;
 658         FillHistogram(key,p12.M() ,p12.Pt()) ;
 659         if(ph1->IsPhoton() && ph2->IsPhoton()){
 660           snprintf(key,55,"hMassPtDispwou_cen%d",fCenBin) ;
 661           FillHistogram(key,p12.M() ,p12.Pt()) ;
 662         }
 663         if(a<0.9){
 664           snprintf(key,55,"hMassPtDisp_a09_cen%d",fCenBin) ;
 665           FillHistogram(key,p12.M() ,p12.Pt()) ;
 666           if(a<0.8){
 667             snprintf(key,55,"hMassPtDisp_a08_cen%d",fCenBin) ;
 668             FillHistogram(key,p12.M() ,p12.Pt()) ;
 669             if(a<0.7){
 670               snprintf(key,55,"hMassPtDisp_a07_cen%d",fCenBin) ;
 671               FillHistogram(key,p12.M() ,p12.Pt()) ;
 672             }
 673           }
 674         }
 675
 676         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 677           snprintf(key,55,"hMassPtBoth_cen%d",fCenBin) ;
 678           FillHistogram(key,p12.M() ,p12.Pt()) ;
 679           snprintf(key,55,"hMassPtBothcore_cen%d",fCenBin) ;
 680           FillHistogram(key,pv12.M(), pv12.Pt()) ;
 681           if(a<0.9){
 682             snprintf(key,55,"hMassPtBoth_a09_cen%d",fCenBin) ;
 683             FillHistogram(key,p12.M() ,p12.Pt()) ;
 684             if(a<0.8){
 685               snprintf(key,55,"hMassPtBoth_a08_cen%d",fCenBin) ;
 686               FillHistogram(key,p12.M() ,p12.Pt()) ;
 687               if(a<0.7){
 688                 snprintf(key,55,"hMassPtBoth_a07_cen%d",fCenBin) ;
 689                 FillHistogram(key,p12.M() ,p12.Pt()) ;
 690              }
 691             }
 692           }
 693         }
 694       }
 695     } // end of loop i2
 696   } // end of loop i1
 697
 698   //now mixed
 699   for (Int_t i1=0; i1<inPHOS; i1++) {
 700     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 701     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
 702       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
 703       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
 704         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
 705         TLorentzVector p12  = *ph1  + *ph2;
 706         TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 707         Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
 708
 709         snprintf(key,55,"hMiMassPtAll_cen%d",fCenBin) ;
 710         FillHistogram(key,p12.M() ,p12.Pt()) ;
 711         snprintf(key,55,"hMiMassPtAllcore_cen%d",fCenBin) ;
 712         FillHistogram(key,pv12.M(), pv12.Pt()) ;
 713         if(ph1->IsPhoton() && ph2->IsPhoton()){
 714           snprintf(key,55,"hMiMassPtAllwou_cen%d",fCenBin) ;
 715           FillHistogram(key,p12.M() ,p12.Pt()) ;
 716         }
 717         if(a<0.9){
 718           snprintf(key,55,"hMiMassPtAll_a09_cen%d",fCenBin) ;
 719           FillHistogram(key,p12.M() ,p12.Pt()) ;
 720           if(a<0.8){
 721             snprintf(key,55,"hMiMassPtAll_a08_cen%d",fCenBin) ;
 722             FillHistogram(key,p12.M() ,p12.Pt()) ;
 723             if(a<0.7){
 724               snprintf(key,55,"hMiMassPtAll_a07_cen%d",fCenBin) ;
 725               FillHistogram(key,p12.M() ,p12.Pt()) ;
 726             }
 727           }
 728         }
 729         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 730           snprintf(key,55,"hMiMassPtCPV_cen%d",fCenBin) ;
 731           FillHistogram(key,p12.M() ,p12.Pt()) ;
 732           snprintf(key,55,"hMiMassPtCPVcore_cen%d",fCenBin) ;
 733           FillHistogram(key,pv12.M(), pv12.Pt()) ;
 734           if(a<0.9){
 735             snprintf(key,55,"hMiMassPtCPV_a09_cen%d",fCenBin) ;
 736             FillHistogram(key,p12.M() ,p12.Pt()) ;
 737             if(a<0.8){
 738               snprintf(key,55,"hMiMassPtCPV_a08_cen%d",fCenBin) ;
 739               FillHistogram(key,p12.M() ,p12.Pt()) ;
 740               if(a<0.7){
 741                 snprintf(key,55,"hMiMassPtCPV_a07_cen%d",fCenBin) ;
 742                 FillHistogram(key,p12.M() ,p12.Pt()) ;
 743               }
 744             }
 745           }
 746         }
 747         if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
 748           snprintf(key,55,"hMiMassPtCPV2_cen%d",fCenBin) ;
 749           FillHistogram(key,p12.M() ,p12.Pt()) ;
 750           if(a<0.9){
 751             snprintf(key,55,"hMiMassPtCPV2_a09_cen%d",fCenBin) ;
 752             FillHistogram(key,p12.M() ,p12.Pt()) ;
 753             if(a<0.8){
 754               snprintf(key,55,"hMiMassPtCPV2_a08_cen%d",fCenBin) ;
 755               FillHistogram(key,p12.M() ,p12.Pt()) ;
 756               if(a<0.7){
 757                 snprintf(key,55,"hMiMassPtCPV2_a07_cen%d",fCenBin) ;
 758                 FillHistogram(key,p12.M() ,p12.Pt()) ;
 759               }
 760             }
 761           }
 762         }
 763         if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
 764           snprintf(key,55,"hMiMassPtDisp2_cen%d",fCenBin) ;
 765           FillHistogram(key,p12.M() ,p12.Pt()) ;
 766         }
 767         if(ph1->IsDispOK() && ph2->IsDispOK()){
 768           snprintf(key,55,"hMiMassPtDisp_cen%d",fCenBin) ;
 769           FillHistogram(key,p12.M() ,p12.Pt()) ;
 770           if(ph1->IsPhoton() && ph2->IsPhoton()){
 771             snprintf(key,55,"hMiMassPtDispwou_cen%d",fCenBin) ;
 772             FillHistogram(key,p12.M() ,p12.Pt()) ;
 773           }
 774           if(a<0.9){
 775             snprintf(key,55,"hMiMassPtDisp_a09_cen%d",fCenBin) ;
 776             FillHistogram(key,p12.M() ,p12.Pt()) ;
 777             if(a<0.8){
 778               snprintf(key,55,"hMiMassPtDisp_a08_cen%d",fCenBin) ;
 779               FillHistogram(key,p12.M() ,p12.Pt()) ;
 780               if(a<0.7){
 781                 snprintf(key,55,"hMiMassPtDisp_a07_cen%d",fCenBin) ;
 782                 FillHistogram(key,p12.M() ,p12.Pt()) ;
 783               }
 784             }
 785           }
 786           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 787             snprintf(key,55,"hMiMassPtBoth_cen%d",fCenBin) ;
 788             FillHistogram(key,p12.M() ,p12.Pt()) ;
 789             snprintf(key,55,"hMiMassPtBothcore_cen%d",fCenBin) ;
 790             FillHistogram(key,pv12.M(), pv12.Pt()) ;
 791             if(a<0.9){
 792               snprintf(key,55,"hMiMassPtBoth_a09_cen%d",fCenBin) ;
 793               FillHistogram(key,p12.M() ,p12.Pt()) ;
 794               if(a<0.8){
 795                 snprintf(key,55,"hMiMassPtBoth_a08_cen%d",fCenBin) ;
 796                 FillHistogram(key,p12.M() ,p12.Pt()) ;
 797                 if(a<0.7){
 798                   snprintf(key,55,"hMiMassPtBoth_a07_cen%d",fCenBin) ;
 799                   FillHistogram(key,p12.M() ,p12.Pt()) ;
 800                }
 801               }
 802             }
 803           }
 804         }
 805       } // end of loop i2
 806     }
 807   } // end of loop i1
 808
 809
 810   //Now we either add current events to stack or remove
 811   //If no photons in current event - no need to add it to mixed
 812   if(fPHOSEvent->GetEntriesFast()>0){
 813     prevPHOS->AddFirst(fPHOSEvent) ;
 814     fPHOSEvent=0;
 815     if(prevPHOS->GetSize()>100){//Remove redundant events
 816       TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
 817       prevPHOS->RemoveLast() ;
 818       delete tmp ;
 819     }
 820   }
 821   // Post output data.
 822   PostData(1, fOutputContainer);
 823   fEventCounter++;
 824 }
 825
 826 //________________________________________________________________________
 827 void AliAnalysisTaskPi0Efficiency::Terminate(Option_t *)
 828 {
 829   // Draw result to the screen
 830   // Called once at the end of the query
 831
 832 }
 833
 834 //________________________________________________________________________
 835 Bool_t AliAnalysisTaskPi0Efficiency::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
 836 {
 837   //Check if this channel belogs to the good ones
 838
 839   if(strcmp(det,"PHOS")==0){
 840     if(mod>5 || mod<1){
 841       AliError(Form("No bad map for PHOS module %d ",mod)) ;
 842       return kTRUE ;
 843     }
 844     if(!fPHOSBadMap[mod]){
 845       AliError(Form("No Bad map for PHOS module %d",mod)) ;
 846       return kTRUE ;
 847     }
 848     if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
 849       return kFALSE ;
 850     else
 851       return kTRUE ;
 852   }
 853   else{
 854     AliError(Form("Can not find bad channels for detector %s ",det)) ;
 855   }
 856   return kTRUE ;
 857 }
 858 //_____________________________________________________________________________
 859 void AliAnalysisTaskPi0Efficiency::FillHistogram(const char * key,Double_t x)const{
 860   //FillHistogram
 861   TH1 * hist = dynamic_cast<TH1*>(fOutputContainer->FindObject(key)) ;
 862   if(hist)
 863     hist->Fill(x) ;
 864   else
 865     AliError(Form("can not find histogram (of instance TH1) <%s> ",key)) ;
 866 }
 867 //_____________________________________________________________________________
 868 void AliAnalysisTaskPi0Efficiency::FillHistogram(const char * key,Double_t x,Double_t y)const{
 869   //FillHistogram
 870   TH1 * th1 = dynamic_cast<TH1*> (fOutputContainer->FindObject(key));
 871   if(th1)
 872     th1->Fill(x, y) ;
 873   else
 874     AliError(Form("can not find histogram (of instance TH1) <%s> ",key)) ;
 875 }
 876
 877 //_____________________________________________________________________________
 878 void AliAnalysisTaskPi0Efficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
 879   //Fills 1D histograms with key
 880   TObject * obj = fOutputContainer->FindObject(key);
 881
 882   TH2 * th2 = dynamic_cast<TH2*> (obj);
 883   if(th2) {
 884     th2->Fill(x, y, z) ;
 885     return;
 886   }
 887
 888   TH3 * th3 = dynamic_cast<TH3*> (obj);
 889   if(th3) {
 890     th3->Fill(x, y, z) ;
 891     return;
 892   }
 893
 894   AliError(Form("can not find histogram (of instance TH2) <%s> ",key)) ;
 895 }
 896 //_____________________________________________________________________________
 897 Bool_t AliAnalysisTaskPi0Efficiency::TestLambda(Double_t pt,Double_t l1,Double_t l2){
 898
 899   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
 900   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
 901   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
 902   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
 903   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
 904   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
 905               0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
 906               0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
 907   return (R2<2.5*2.5) ;
 908
 909 }
 910 //_____________________________________________________________________________
 911 Bool_t AliAnalysisTaskPi0Efficiency::TestLambda2(Double_t pt,Double_t l1,Double_t l2){
 912
 913   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
 914   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
 915   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
 916   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
 917   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
 918   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
 919               0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
 920               0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
 921   return (R2<1.5*1.5) ;
 922
 923 }
 924 //___________________________________________________________________________
 925 void AliAnalysisTaskPi0Efficiency::ProcessMC(){
 926   //fill histograms for efficiensy etc. calculation
 927   const Double_t rcut = 1. ; //cut for primary particles
 928   //---------First pi0/eta-----------------------------
 929   char partName[10] ;
 930   char hkey[55] ;
 931
 932   AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
 933   if(!event) return ;
 934   TClonesArray *mcArray = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
 935   for(Int_t i=0;i<mcArray->GetEntriesFast();i++){
 936      AliAODMCParticle* particle =  (AliAODMCParticle*) mcArray->At(i);
 937     if(particle->GetPdgCode() == 111)
 938       snprintf(partName,10,"pi0") ;
 939     else
 940       if(particle->GetPdgCode() == 221)
 941         snprintf(partName,10,"eta") ;
 942       else
 943         if(particle->GetPdgCode() == 22)
 944            snprintf(partName,10,"gamma") ;
 945         else
 946            continue ;
 947
 948     //Primary particle
 949     Double_t r=TMath::Sqrt(particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv());
 950     if(r >rcut)
 951       continue ;
 952
 953     Double_t pt = particle->Pt() ;
 954     //Total number of pi0 with creation radius <1 cm
 955     snprintf(hkey,55,"hMC_all_%s_cen%d",partName,fCenBin) ;
 956     FillHistogram(hkey,pt) ;
 957     if(TMath::Abs(particle->Y())<0.12){
 958       snprintf(hkey,55,"hMC_unitEta_%s_cen%d",partName,fCenBin) ;
 959       FillHistogram(hkey,pt) ;
 960     }
 961
 962     snprintf(hkey,55,"hMC_rap_%s_cen%d",partName,fCenBin) ;
 963     FillHistogram(hkey,particle->Y()) ;
 964
 965     Double_t phi=particle->Phi() ;
 966     while(phi<0.)phi+=TMath::TwoPi() ;
 967     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 968     snprintf(hkey,55,"hMC_phi_%s_cen%d",partName,fCenBin) ;
 969     FillHistogram(hkey,phi) ;
 970
 971
 972     //Check if one of photons converted
 973     if(particle->GetNDaughters()!=2)
 974      continue ; //Do not account Dalitz decays
 975
 976 /*
 977     TParticle * gamma1 = fStack->Particle(particle->GetFirstDaughter());
 978     TParticle * gamma2 = fStack->Particle(particle->GetLastDaughter());
 979     //Number of pi0s decayed into acceptance
 980     Int_t mod1,mod2 ;
 981     Double_t x=0.,z=0. ;
 982     Bool_t hitPHOS1 = fPHOSGeo->ImpactOnEmc(gamma1, mod1, z,x) ;
 983     Bool_t hitPHOS2 = fPHOSGeo->ImpactOnEmc(gamma2, mod2, z,x) ;
 984
 985     Bool_t goodPair=kFALSE ;
 986     if( hitPHOS1 && hitPHOS2){
 987       sprintf(hkey,"hMC_PHOSacc_%s",partName) ;
 988       FillHistogram(hkey,pt) ;
 989       goodPair=kTRUE ;
 990     }
 991
 992 */
 993   }
 994
 995   //Now calculate "Real" distribution of clusters with primary
 996   TClonesArray cluPrim("AliCaloPhoton",200) ; //clusters with primary
 997   TClonesArray * clusters = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
 998   AliAODCaloCells * cells = (AliAODCaloCells *)event->FindListObject("EmbeddedPHOScells") ;
 999   Int_t multClust = clusters->GetEntriesFast();
1000   Int_t inPHOS=0 ;
1001   Double_t vtx0[3] = {0,0,0};
1002   for (Int_t i=0; i<multClust; i++) {
1003     AliAODCaloCluster *clu = (AliAODCaloCluster*)clusters->At(i);
1004     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
1005     if(clu->GetLabel()<0) continue ;
1006
1007     Float_t  position[3];
1008     clu->GetPosition(position);
1009     TVector3 global(position) ;
1010     Int_t relId[4] ;
1011     fPHOSGeo->GlobalPos2RelId(global,relId) ;
1012     Int_t mod  = relId[0] ;
1013     Int_t cellX = relId[2];
1014     Int_t cellZ = relId[3] ;
1015     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
1016       continue ;
1017     if(clu->GetNCells()<3)
1018       continue ;
1019
1020     TLorentzVector pv1 ;
1021     clu->GetMomentum(pv1 ,vtx0);
1022
1023     if(inPHOS>=cluPrim.GetSize()){
1024       cluPrim.Expand(inPHOS+50) ;
1025     }
1026     AliCaloPhoton * ph = new(cluPrim[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
1027     //AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
1028     ph->SetModule(mod) ;
1029     AliPHOSAodCluster cluPHOS1(*clu);
1030     cluPHOS1.Recalibrate(fPHOSCalibData,cells); // modify the cell energies
1031     Double_t ecore=CoreEnergy(&cluPHOS1) ;
1032     pv1*= ecore/pv1.E() ;
1033     ph->SetMomV2(&pv1) ;
1034     ph->SetNCells(clu->GetNCells());
1035     ph->SetDispBit(TestLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
1036     ph->SetDisp2Bit(TestLambda2(clu->E(),clu->GetM20(),clu->GetM02())) ;
1037     ph->SetCPVBit(clu->GetEmcCpvDistance()>2.) ; //radius in sigmas
1038     ph->SetCPV2Bit(clu->GetEmcCpvDistance()>4.) ;
1039     ph->SetPhoton(clu->GetNExMax()<2); // Remember, if it is unfolded
1040
1041
1042     inPHOS++ ;
1043
1044   }
1045
1046   //Single photon
1047   char key[55] ;
1048   for (Int_t i1=0; i1<inPHOS; i1++) {
1049     AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
1050     snprintf(key,55,"hMCPhotAll_cen%d",fCenBin) ;
1051     FillHistogram(key,ph1->Pt()) ;
1052     snprintf(key,55,"hMCPhotAllcore_cen%d",fCenBin) ;
1053     FillHistogram(key,ph1->GetMomV2()->Pt()) ;
1054     if(ph1->IsPhoton()){
1055       snprintf(key,55,"hMCPhotAllwou_cen%d",fCenBin) ;
1056       FillHistogram(key,ph1->Pt()) ;
1057     }
1058     if(ph1->IsCPVOK() ){
1059       snprintf(key,55,"hMCPhotCPV_cen%d",fCenBin) ;
1060       FillHistogram(key,ph1->Pt()) ;
1061       snprintf(key,55,"hMCPhotCPVcore_cen%d",fCenBin) ;
1062       FillHistogram(key,ph1->GetMomV2()->Pt()) ;
1063
1064     }
1065     if(ph1->IsCPV2OK() ){
1066       snprintf(key,55,"hMCPhotCPV2_cen%d",fCenBin) ;
1067       FillHistogram(key,ph1->Pt()) ;
1068
1069     }
1070     if(ph1->IsDisp2OK()){
1071       snprintf(key,55,"hMCPhotDisp2_cen%d",fCenBin) ;
1072       FillHistogram(key,ph1->Pt()) ;
1073     }
1074     if(ph1->IsDispOK()){
1075       snprintf(key,55,"hMCPhotDisp_cen%d",fCenBin) ;
1076       FillHistogram(key,ph1->Pt()) ;
1077       if(ph1->IsPhoton()){
1078         snprintf(key,55,"hMCPhotDispwou_cen%d",fCenBin) ;
1079         FillHistogram(key,ph1->Pt()) ;
1080       }
1081       if(ph1->IsCPVOK()){
1082         snprintf(key,55,"hMCPhotBoth_cen%d",fCenBin) ;
1083         FillHistogram(key,ph1->Pt()) ;
1084         snprintf(key,55,"hMCPhotBothcore_cen%d",fCenBin) ;
1085         FillHistogram(key,ph1->GetMomV2()->Pt()) ;
1086       }
1087     } // end of loop i2
1088   } // end of loop i1
1089
1090   // Fill Real disribution
1091   for (Int_t i1=0; i1<inPHOS-1; i1++) {
1092     AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
1093     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
1094       AliCaloPhoton * ph2=(AliCaloPhoton*)cluPrim.At(i2) ;
1095       TLorentzVector p12  = *ph1  + *ph2;
1096       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1097       Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
1098
1099       snprintf(key,55,"hMCMassPtAll_cen%d",fCenBin) ;
1100       FillHistogram(key,p12.M() ,p12.Pt()) ;
1101       snprintf(key,55,"hMCMassPtAllcore_cen%d",fCenBin) ;
1102       FillHistogram(key,pv12.M(), pv12.Pt()) ;
1103       if(ph1->IsPhoton()&& ph2->IsPhoton()){
1104         snprintf(key,55,"hMCMassPtAllwou_cen%d",fCenBin) ;
1105         FillHistogram(key,p12.M() ,p12.Pt()) ;
1106       }
1107       if(a<0.9){
1108         snprintf(key,55,"hMCMassPtAll_a09_cen%d",fCenBin) ;
1109         FillHistogram(key,p12.M() ,p12.Pt()) ;
1110         if(a<0.8){
1111           snprintf(key,55,"hMCMassPtAll_a08_cen%d",fCenBin) ;
1112           FillHistogram(key,p12.M() ,p12.Pt()) ;
1113           if(a<0.7){
1114             snprintf(key,55,"hMCMassPtAll_a07_cen%d",fCenBin) ;
1115             FillHistogram(key,p12.M() ,p12.Pt()) ;
1116           }
1117         }
1118       }
1119
1120
1121            if(ph1->Module()==1 && ph2->Module()==1)
1122             FillHistogram("hMCPi0M11",p12.M(),p12.Pt() );
1123           else if(ph1->Module()==2 && ph2->Module()==2)
1124             FillHistogram("hMCPi0M22",p12.M(),p12.Pt() );
1125           else if(ph1->Module()==3 && ph2->Module()==3)
1126             FillHistogram("hMCPi0M33",p12.M(),p12.Pt() );
1127           else if(ph1->Module()==1 && ph2->Module()==2)
1128             FillHistogram("hMCPi0M12",p12.M(),p12.Pt() );
1129           else if(ph1->Module()==1 && ph2->Module()==3)
1130             FillHistogram("hMCPi0M13",p12.M(),p12.Pt() );
1131           else if(ph1->Module()==2 && ph2->Module()==3)
1132             FillHistogram("hMCPi0M23",p12.M(),p12.Pt() );
1133
1134
1135
1136
1137
1138       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1139         snprintf(key,55,"hMCMassPtCPV_cen%d",fCenBin) ;
1140         FillHistogram(key,p12.M() ,p12.Pt()) ;
1141         snprintf(key,55,"hMCMassPtCPVcore_cen%d",fCenBin) ;
1142         FillHistogram(key,pv12.M(), pv12.Pt()) ;
1143         if(a<0.9){
1144           snprintf(key,55,"hMCMassPtCPV_a09_cen%d",fCenBin) ;
1145           FillHistogram(key,p12.M() ,p12.Pt()) ;
1146           if(a<0.8){
1147             snprintf(key,55,"hMCMassPtCPV_a08_cen%d",fCenBin) ;
1148             FillHistogram(key,p12.M() ,p12.Pt()) ;
1149             if(a<0.7){
1150               snprintf(key,55,"hMCMassPtCPV_a07_cen%d",fCenBin) ;
1151               FillHistogram(key,p12.M() ,p12.Pt()) ;
1152             }
1153           }
1154         }
1155       }
1156       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1157         snprintf(key,55,"hMCMassPtCPV2_cen%d",fCenBin) ;
1158         FillHistogram(key,p12.M() ,p12.Pt()) ;
1159         if(a<0.9){
1160           snprintf(key,55,"hMCMassPtCPV2_a09_cen%d",fCenBin) ;
1161           FillHistogram(key,p12.M() ,p12.Pt()) ;
1162           if(a<0.8){
1163             snprintf(key,55,"hMCMassPtCPV2_a08_cen%d",fCenBin) ;
1164             FillHistogram(key,p12.M() ,p12.Pt()) ;
1165             if(a<0.7){
1166               snprintf(key,55,"hMCMassPtCPV2_a07_cen%d",fCenBin) ;
1167               FillHistogram(key,p12.M() ,p12.Pt()) ;
1168             }
1169           }
1170         }
1171       }
1172       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1173         snprintf(key,55,"hMCMassPtDisp2_cen%d",fCenBin) ;
1174         FillHistogram(key,p12.M() ,p12.Pt()) ;
1175       }
1176       if(ph1->IsDispOK() && ph2->IsDispOK()){
1177         snprintf(key,55,"hMCMassPtDisp_cen%d",fCenBin) ;
1178         FillHistogram(key,p12.M() ,p12.Pt()) ;
1179         if(ph1->IsPhoton()&& ph2->IsPhoton()){
1180           snprintf(key,55,"hMCMassPtDispwou_cen%d",fCenBin) ;
1181           FillHistogram(key,p12.M() ,p12.Pt()) ;
1182         }
1183         if(a<0.9){
1184           snprintf(key,55,"hMCMassPtDisp_a09_cen%d",fCenBin) ;
1185           FillHistogram(key,p12.M() ,p12.Pt()) ;
1186           if(a<0.8){
1187             snprintf(key,55,"hMCMassPtDisp_a08_cen%d",fCenBin) ;
1188             FillHistogram(key,p12.M() ,p12.Pt()) ;
1189             if(a<0.7){
1190               snprintf(key,55,"hMCMassPtDisp_a07_cen%d",fCenBin) ;
1191               FillHistogram(key,p12.M() ,p12.Pt()) ;
1192             }
1193           }
1194         }
1195
1196         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1197           snprintf(key,55,"hMCMassPtBoth_cen%d",fCenBin) ;
1198           FillHistogram(key,p12.M() ,p12.Pt()) ;
1199           snprintf(key,55,"hMCMassPtBothcore_cen%d",fCenBin) ;
1200           FillHistogram(key,pv12.M(), pv12.Pt()) ;
1201           if(a<0.9){
1202             snprintf(key,55,"hMCMassPtBoth_a09_cen%d",fCenBin) ;
1203             FillHistogram(key,p12.M() ,p12.Pt()) ;
1204             if(a<0.8){
1205               snprintf(key,55,"hMCMassPtBoth_a08_cen%d",fCenBin) ;
1206               FillHistogram(key,p12.M() ,p12.Pt()) ;
1207               if(a<0.7){
1208                 snprintf(key,55,"hMCMassPtBoth_a07_cen%d",fCenBin) ;
1209                 FillHistogram(key,p12.M() ,p12.Pt()) ;
1210               }
1211             }
1212           }
1213         }
1214       }
1215     } // end of loop i2
1216   } // end of loop i1
1217 }
1218 //____________________________________________________________________________
1219 Double_t  AliAnalysisTaskPi0Efficiency::CoreEnergy(AliPHOSAodCluster * clu){
1220   //calculate energy of the cluster in the circle with radius distanceCut around the maximum
1221
1222   //Can not use already calculated coordinates?
1223   //They have incidence correction...
1224   const Double_t distanceCut =3.5 ;
1225   const Double_t logWeight=4.5 ;
1226
1227   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1228 // Calculates the center of gravity in the local PHOS-module coordinates
1229   Float_t wtot = 0;
1230   Double_t xc[100]={0} ;
1231   Double_t zc[100]={0} ;
1232   Double_t x = 0 ;
1233   Double_t z = 0 ;
1234   Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
1235   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1236     Int_t relid[4] ;
1237     Float_t xi ;
1238     Float_t zi ;
1239     fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
1240     fPHOSGeo->RelPosInModule(relid, xi, zi);
1241     xc[iDigit]=xi ;
1242     zc[iDigit]=zi ;
1243     if (clu->E()>0 && elist[iDigit]>0) {
1244       Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1245       x    += xc[iDigit] * w ;
1246       z    += zc[iDigit] * w ;
1247       wtot += w ;
1248     }
1249   }
1250   if (wtot>0) {
1251     x /= wtot ;
1252     z /= wtot ;
1253   }
1254   Double_t coreE=0. ;
1255   for(Int_t iDigit=0; iDigit < mulDigit; iDigit++) {
1256     Double_t distance = TMath::Sqrt((xc[iDigit]-x)*(xc[iDigit]-x)+(zc[iDigit]-z)*(zc[iDigit]-z)) ;
1257     if(distance < distanceCut)
1258       coreE += elist[iDigit] ;
1259   }
1260   //Apply non-linearity correction
1261   return (0.0241+1.0504*coreE+0.000249*coreE*coreE) ;
1262 }
1263
1264