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 = event->GetHeader() ;
 422
 423   // Checks if we have a primary vertex
 424   // Get primary vertices form ESD
 425   const AliAODVertex *esdVertex5 = event->GetPrimaryVertex();
 426
 427  // don't rely on ESD vertex, assume (0,0,0)
 428   Double_t vtx0[3] ={0.,0.,0.};
 429
 430
 431   FillHistogram("hZvertex",esdVertex5->GetZ());
 432   if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
 433     PostData(1, fOutputContainer);
 434     return;
 435   }
 436   FillHistogram("hSelEvents",2.5) ;
 437
 438   //Vtx class z-bin
 439   //  Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
 440   //  if(zvtx<0)zvtx=0 ;
 441   //  if(zvtx>9)zvtx=9 ;
 442   Int_t zvtx=0 ;
 443
 444 //  fCentrality=header->GetCentralityP()->GetCentralityPercentile("V0M"); // returns the centrality percentile,
 445 //                                                          //a float from 0 to 100 (or to the trigger efficiency)
 446    fCentrality=header->GetZDCN2Energy() ;
 447
 448   if( fCentrality < 0. || fCentrality>80.){
 449     PostData(1, fOutputContainer);
 450     return;
 451   }
 452   FillHistogram("hSelEvents",3.5) ;
 453   Float_t bins[7]={0.,5.,10.,20.,40.,60.,80.} ;
 454   fCenBin=0 ;
 455   while(fCenBin<6 && fCentrality > bins[fCenBin+1])
 456     fCenBin++ ;
 457
 458
 459   //reaction plain
 460   fRPfull= header->GetZDCN1Energy() ;
 461   if(fRPfull==999){ //reaction plain was not defined
 462     PostData(1, fOutputContainer);
 463     return;
 464   }
 465
 466   FillHistogram("hSelEvents",4.5) ;
 467   //All event selections done
 468   FillHistogram("hCentrality",fCentrality) ;
 469   //Reaction plain is defined in the range (-pi/2;pi/2)
 470   //We have 10 bins
 471   Int_t irp=Int_t(10.*fRPfull/TMath::Pi());
 472   if(irp>9)irp=9 ;
 473
 474   if(!fPHOSEvents[zvtx][fCenBin][irp])
 475     fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
 476   TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
 477
 478   // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
 479   if(fEventCounter == 0) {
 480     for(Int_t mod=0; mod<5; mod++) {
 481       const TGeoHMatrix* m =header->GetPHOSMatrix(mod) ;
 482       fPHOSGeo->SetMisalMatrix(m,mod) ;
 483       Printf("PHOS geo matrix for module # %d is set: %p\n", mod,m);
 484     }
 485     fEventCounter++ ;
 486   }
 487
 488   ProcessMC() ;
 489
 490   if(fPHOSEvent)
 491     fPHOSEvent->Clear() ;
 492   else
 493     fPHOSEvent = new TClonesArray("AliCaloPhoton",200) ;
 494
 495
 496   char key[55] ;
 497   Int_t inPHOS=0 ;
 498   TVector3 vertex(vtx0);
 499   TClonesArray * clusters = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
 500   AliAODCaloCells * cells = (AliAODCaloCells*)event->FindListObject("EmbeddedPHOScells") ;
 501   Int_t multClust = clusters->GetEntriesFast();
 502   for (Int_t i=0; i<multClust; i++) {
 503     AliAODCaloCluster *clu = (AliAODCaloCluster*) clusters->At(i);
 504     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
 505
 506     Float_t  position[3];
 507     clu->GetPosition(position);
 508     TVector3 global(position) ;
 509     Int_t relId[4] ;
 510     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 511     Int_t mod  = relId[0] ;
 512     Int_t cellX = relId[2];
 513     Int_t cellZ = relId[3] ;
 514     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
 515       continue ;
 516     if(clu->GetNCells()<3)
 517       continue ;
 518     if(clu->GetM02()<0.2)
 519       continue ;
 520
 521     snprintf(key,55,"hCluM%d",mod) ;
 522     FillHistogram(key,cellX,cellZ,1.);
 523
 524     TLorentzVector pv1 ;
 525     clu->GetMomentum(pv1 ,vtx0);
 526
 527     if(inPHOS>=fPHOSEvent->GetSize()){
 528       fPHOSEvent->Expand(inPHOS+50) ;
 529     }
 530     new((*fPHOSEvent)[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
 531     AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
 532     ph->SetModule(mod) ;
 533     AliPHOSAodCluster cluPHOS1(*clu);
 534     cluPHOS1.Recalibrate(fPHOSCalibData,cells); // modify the cell energies
 535     Double_t ecore=CoreEnergy(&cluPHOS1) ;
 536     pv1*= ecore/pv1.E() ;
 537     ph->SetMomV2(&pv1) ;
 538     ph->SetNCells(clu->GetNCells());
 539     ph->SetDispBit(TestLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
 540     ph->SetDisp2Bit(TestLambda2(clu->E(),clu->GetM20(),clu->GetM02())) ;
 541     ph->SetCPVBit(clu->GetEmcCpvDistance()>2.) ;
 542     ph->SetCPV2Bit(clu->GetEmcCpvDistance()>4.) ;
 543     ph->SetPhoton(clu->GetNExMax()<2); // Remember, if it is unfolded
 544
 545     inPHOS++ ;
 546   }
 547   //Single photon
 548   for (Int_t i1=0; i1<inPHOS; i1++) {
 549     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 550     snprintf(key,55,"hPhotAll_cen%d",fCenBin) ;
 551     FillHistogram(key,ph1->Pt()) ;
 552     snprintf(key,55,"hPhotAllcore_cen%d",fCenBin) ;
 553     FillHistogram(key,ph1->GetMomV2()->Pt()) ;
 554     if(ph1->IsPhoton()){
 555       snprintf(key,55,"hPhotAllwou_cen%d",fCenBin) ;
 556       FillHistogram(key,ph1->Pt()) ;
 557     }
 558     if(ph1->IsCPVOK() ){
 559       snprintf(key,55,"hPhotCPV_cen%d",fCenBin) ;
 560       FillHistogram(key,ph1->Pt()) ;
 561       snprintf(key,55,"hPhotCPVcore_cen%d",fCenBin) ;
 562       FillHistogram(key,ph1->GetMomV2()->Pt()) ;
 563     }
 564     if(ph1->IsCPV2OK() ){
 565       snprintf(key,55,"hPhotCPV2_cen%d",fCenBin) ;
 566       FillHistogram(key,ph1->Pt()) ;
 567     }
 568     if(ph1->IsDisp2OK()){
 569       snprintf(key,55,"hPhotDisp2_cen%d",fCenBin) ;
 570       FillHistogram(key,ph1->Pt()) ;
 571     }
 572     if(ph1->IsDispOK()){
 573       snprintf(key,55,"hPhotDisp_cen%d",fCenBin) ;
 574       FillHistogram(key,ph1->Pt()) ;
 575       if(ph1->IsPhoton()){
 576         snprintf(key,55,"hPhotDispwou_cen%d",fCenBin) ;
 577         FillHistogram(key,ph1->Pt()) ;
 578       }
 579       if(ph1->IsCPVOK()){
 580         snprintf(key,55,"hPhotBoth_cen%d",fCenBin) ;
 581         FillHistogram(key,ph1->Pt()) ;
 582         snprintf(key,55,"hPhotBothcore_cen%d",fCenBin) ;
 583         FillHistogram(key,ph1->GetMomV2()->Pt()) ;
 584       }
 585     } // end of loop i2
 586   } // end of loop i1
 587
 588   // Fill Real disribution
 589   for (Int_t i1=0; i1<inPHOS-1; i1++) {
 590     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 591     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
 592       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent->At(i2) ;
 593       TLorentzVector p12  = *ph1  + *ph2;
 594       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 595       Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
 596
 597       snprintf(key,55,"hMassPtAll_cen%d",fCenBin) ;
 598       FillHistogram(key,p12.M() ,p12.Pt()) ;
 599       snprintf(key,55,"hMassPtAllcore_cen%d",fCenBin) ;
 600       FillHistogram(key,pv12.M(), pv12.Pt()) ;
 601       if(ph1->IsPhoton() && ph2->IsPhoton()){
 602         snprintf(key,55,"hMassPtAllwou_cen%d",fCenBin) ;
 603         FillHistogram(key,p12.M() ,p12.Pt()) ;
 604       }
 605       if(a<0.9){
 606         snprintf(key,55,"hMassPtAll_a09_cen%d",fCenBin) ;
 607         FillHistogram(key,p12.M() ,p12.Pt()) ;
 608         if(a<0.8){
 609           snprintf(key,55,"hMassPtAll_a08_cen%d",fCenBin) ;
 610           FillHistogram(key,p12.M() ,p12.Pt()) ;
 611           if(a<0.7){
 612             snprintf(key,55,"hMassPtAll_a07_cen%d",fCenBin) ;
 613             FillHistogram(key,p12.M() ,p12.Pt()) ;
 614           }
 615         }
 616       }
 617       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 618         snprintf(key,55,"hMassPtCPV_cen%d",fCenBin) ;
 619         FillHistogram(key,p12.M() ,p12.Pt()) ;
 620         snprintf(key,55,"hMassPtCPVcore_cen%d",fCenBin) ;
 621         FillHistogram(key,pv12.M(), pv12.Pt()) ;
 622         if(a<0.9){
 623           snprintf(key,55,"hMassPtCPV_a09_cen%d",fCenBin) ;
 624           FillHistogram(key,p12.M() ,p12.Pt()) ;
 625           if(a<0.8){
 626             snprintf(key,55,"hMassPtCPV_a08_cen%d",fCenBin) ;
 627             FillHistogram(key,p12.M() ,p12.Pt()) ;
 628             if(a<0.7){
 629               snprintf(key,55,"hMassPtCPV_a07_cen%d",fCenBin) ;
 630               FillHistogram(key,p12.M() ,p12.Pt()) ;
 631             }
 632           }
 633         }
 634       }
 635       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
 636         snprintf(key,55,"hMassPtCPV2_cen%d",fCenBin) ;
 637         FillHistogram(key,p12.M() ,p12.Pt()) ;
 638         if(a<0.9){
 639           snprintf(key,55,"hMassPtCPV2_a09_cen%d",fCenBin) ;
 640           FillHistogram(key,p12.M() ,p12.Pt()) ;
 641           if(a<0.8){
 642             snprintf(key,55,"hMassPtCPV2_a08_cen%d",fCenBin) ;
 643             FillHistogram(key,p12.M() ,p12.Pt()) ;
 644             if(a<0.7){
 645               snprintf(key,55,"hMassPtCPV2_a07_cen%d",fCenBin) ;
 646               FillHistogram(key,p12.M() ,p12.Pt()) ;
 647             }
 648           }
 649         }
 650       }
 651       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
 652         snprintf(key,55,"hMassPtDisp2_cen%d",fCenBin) ;
 653         FillHistogram(key,p12.M() ,p12.Pt()) ;
 654       }
 655       if(ph1->IsDispOK() && ph2->IsDispOK()){
 656         snprintf(key,55,"hMassPtDisp_cen%d",fCenBin) ;
 657         FillHistogram(key,p12.M() ,p12.Pt()) ;
 658         if(ph1->IsPhoton() && ph2->IsPhoton()){
 659           snprintf(key,55,"hMassPtDispwou_cen%d",fCenBin) ;
 660           FillHistogram(key,p12.M() ,p12.Pt()) ;
 661         }
 662         if(a<0.9){
 663           snprintf(key,55,"hMassPtDisp_a09_cen%d",fCenBin) ;
 664           FillHistogram(key,p12.M() ,p12.Pt()) ;
 665           if(a<0.8){
 666             snprintf(key,55,"hMassPtDisp_a08_cen%d",fCenBin) ;
 667             FillHistogram(key,p12.M() ,p12.Pt()) ;
 668             if(a<0.7){
 669               snprintf(key,55,"hMassPtDisp_a07_cen%d",fCenBin) ;
 670               FillHistogram(key,p12.M() ,p12.Pt()) ;
 671             }
 672           }
 673         }
 674
 675         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 676           snprintf(key,55,"hMassPtBoth_cen%d",fCenBin) ;
 677           FillHistogram(key,p12.M() ,p12.Pt()) ;
 678           snprintf(key,55,"hMassPtBothcore_cen%d",fCenBin) ;
 679           FillHistogram(key,pv12.M(), pv12.Pt()) ;
 680           if(a<0.9){
 681             snprintf(key,55,"hMassPtBoth_a09_cen%d",fCenBin) ;
 682             FillHistogram(key,p12.M() ,p12.Pt()) ;
 683             if(a<0.8){
 684               snprintf(key,55,"hMassPtBoth_a08_cen%d",fCenBin) ;
 685               FillHistogram(key,p12.M() ,p12.Pt()) ;
 686               if(a<0.7){
 687                 snprintf(key,55,"hMassPtBoth_a07_cen%d",fCenBin) ;
 688                 FillHistogram(key,p12.M() ,p12.Pt()) ;
 689              }
 690             }
 691           }
 692         }
 693       }
 694     } // end of loop i2
 695   } // end of loop i1
 696
 697   //now mixed
 698   for (Int_t i1=0; i1<inPHOS; i1++) {
 699     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 700     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
 701       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
 702       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
 703         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
 704         TLorentzVector p12  = *ph1  + *ph2;
 705         TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 706         Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
 707
 708         snprintf(key,55,"hMiMassPtAll_cen%d",fCenBin) ;
 709         FillHistogram(key,p12.M() ,p12.Pt()) ;
 710         snprintf(key,55,"hMiMassPtAllcore_cen%d",fCenBin) ;
 711         FillHistogram(key,pv12.M(), pv12.Pt()) ;
 712         if(ph1->IsPhoton() && ph2->IsPhoton()){
 713           snprintf(key,55,"hMiMassPtAllwou_cen%d",fCenBin) ;
 714           FillHistogram(key,p12.M() ,p12.Pt()) ;
 715         }
 716         if(a<0.9){
 717           snprintf(key,55,"hMiMassPtAll_a09_cen%d",fCenBin) ;
 718           FillHistogram(key,p12.M() ,p12.Pt()) ;
 719           if(a<0.8){
 720             snprintf(key,55,"hMiMassPtAll_a08_cen%d",fCenBin) ;
 721             FillHistogram(key,p12.M() ,p12.Pt()) ;
 722             if(a<0.7){
 723               snprintf(key,55,"hMiMassPtAll_a07_cen%d",fCenBin) ;
 724               FillHistogram(key,p12.M() ,p12.Pt()) ;
 725             }
 726           }
 727         }
 728         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 729           snprintf(key,55,"hMiMassPtCPV_cen%d",fCenBin) ;
 730           FillHistogram(key,p12.M() ,p12.Pt()) ;
 731           snprintf(key,55,"hMiMassPtCPVcore_cen%d",fCenBin) ;
 732           FillHistogram(key,pv12.M(), pv12.Pt()) ;
 733           if(a<0.9){
 734             snprintf(key,55,"hMiMassPtCPV_a09_cen%d",fCenBin) ;
 735             FillHistogram(key,p12.M() ,p12.Pt()) ;
 736             if(a<0.8){
 737               snprintf(key,55,"hMiMassPtCPV_a08_cen%d",fCenBin) ;
 738               FillHistogram(key,p12.M() ,p12.Pt()) ;
 739               if(a<0.7){
 740                 snprintf(key,55,"hMiMassPtCPV_a07_cen%d",fCenBin) ;
 741                 FillHistogram(key,p12.M() ,p12.Pt()) ;
 742               }
 743             }
 744           }
 745         }
 746         if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
 747           snprintf(key,55,"hMiMassPtCPV2_cen%d",fCenBin) ;
 748           FillHistogram(key,p12.M() ,p12.Pt()) ;
 749           if(a<0.9){
 750             snprintf(key,55,"hMiMassPtCPV2_a09_cen%d",fCenBin) ;
 751             FillHistogram(key,p12.M() ,p12.Pt()) ;
 752             if(a<0.8){
 753               snprintf(key,55,"hMiMassPtCPV2_a08_cen%d",fCenBin) ;
 754               FillHistogram(key,p12.M() ,p12.Pt()) ;
 755               if(a<0.7){
 756                 snprintf(key,55,"hMiMassPtCPV2_a07_cen%d",fCenBin) ;
 757                 FillHistogram(key,p12.M() ,p12.Pt()) ;
 758               }
 759             }
 760           }
 761         }
 762         if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
 763           snprintf(key,55,"hMiMassPtDisp2_cen%d",fCenBin) ;
 764           FillHistogram(key,p12.M() ,p12.Pt()) ;
 765         }
 766         if(ph1->IsDispOK() && ph2->IsDispOK()){
 767           snprintf(key,55,"hMiMassPtDisp_cen%d",fCenBin) ;
 768           FillHistogram(key,p12.M() ,p12.Pt()) ;
 769           if(ph1->IsPhoton() && ph2->IsPhoton()){
 770             snprintf(key,55,"hMiMassPtDispwou_cen%d",fCenBin) ;
 771             FillHistogram(key,p12.M() ,p12.Pt()) ;
 772           }
 773           if(a<0.9){
 774             snprintf(key,55,"hMiMassPtDisp_a09_cen%d",fCenBin) ;
 775             FillHistogram(key,p12.M() ,p12.Pt()) ;
 776             if(a<0.8){
 777               snprintf(key,55,"hMiMassPtDisp_a08_cen%d",fCenBin) ;
 778               FillHistogram(key,p12.M() ,p12.Pt()) ;
 779               if(a<0.7){
 780                 snprintf(key,55,"hMiMassPtDisp_a07_cen%d",fCenBin) ;
 781                 FillHistogram(key,p12.M() ,p12.Pt()) ;
 782               }
 783             }
 784           }
 785           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 786             snprintf(key,55,"hMiMassPtBoth_cen%d",fCenBin) ;
 787             FillHistogram(key,p12.M() ,p12.Pt()) ;
 788             snprintf(key,55,"hMiMassPtBothcore_cen%d",fCenBin) ;
 789             FillHistogram(key,pv12.M(), pv12.Pt()) ;
 790             if(a<0.9){
 791               snprintf(key,55,"hMiMassPtBoth_a09_cen%d",fCenBin) ;
 792               FillHistogram(key,p12.M() ,p12.Pt()) ;
 793               if(a<0.8){
 794                 snprintf(key,55,"hMiMassPtBoth_a08_cen%d",fCenBin) ;
 795                 FillHistogram(key,p12.M() ,p12.Pt()) ;
 796                 if(a<0.7){
 797                   snprintf(key,55,"hMiMassPtBoth_a07_cen%d",fCenBin) ;
 798                   FillHistogram(key,p12.M() ,p12.Pt()) ;
 799                }
 800               }
 801             }
 802           }
 803         }
 804       } // end of loop i2
 805     }
 806   } // end of loop i1
 807
 808
 809   //Now we either add current events to stack or remove
 810   //If no photons in current event - no need to add it to mixed
 811   if(fPHOSEvent->GetEntriesFast()>0){
 812     prevPHOS->AddFirst(fPHOSEvent) ;
 813     fPHOSEvent=0;
 814     if(prevPHOS->GetSize()>100){//Remove redundant events
 815       TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
 816       prevPHOS->RemoveLast() ;
 817       delete tmp ;
 818     }
 819   }
 820   // Post output data.
 821   PostData(1, fOutputContainer);
 822   fEventCounter++;
 823 }
 824
 825 //________________________________________________________________________
 826 void AliAnalysisTaskPi0Efficiency::Terminate(Option_t *)
 827 {
 828   // Draw result to the screen
 829   // Called once at the end of the query
 830
 831 }
 832
 833 //________________________________________________________________________
 834 Bool_t AliAnalysisTaskPi0Efficiency::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
 835 {
 836   //Check if this channel belogs to the good ones
 837
 838   if(strcmp(det,"PHOS")==0){
 839     if(mod>5 || mod<1){
 840       AliError(Form("No bad map for PHOS module %d ",mod)) ;
 841       return kTRUE ;
 842     }
 843     if(!fPHOSBadMap[mod]){
 844       AliError(Form("No Bad map for PHOS module %d",mod)) ;
 845       return kTRUE ;
 846     }
 847     if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
 848       return kFALSE ;
 849     else
 850       return kTRUE ;
 851   }
 852   else{
 853     AliError(Form("Can not find bad channels for detector %s ",det)) ;
 854   }
 855   return kTRUE ;
 856 }
 857 //_____________________________________________________________________________
 858 void AliAnalysisTaskPi0Efficiency::FillHistogram(const char * key,Double_t x)const{
 859   //FillHistogram
 860   TH1 * hist = dynamic_cast<TH1*>(fOutputContainer->FindObject(key)) ;
 861   if(hist)
 862     hist->Fill(x) ;
 863   else
 864     AliError(Form("can not find histogram (of instance TH1) <%s> ",key)) ;
 865 }
 866 //_____________________________________________________________________________
 867 void AliAnalysisTaskPi0Efficiency::FillHistogram(const char * key,Double_t x,Double_t y)const{
 868   //FillHistogram
 869   TH1 * th1 = dynamic_cast<TH1*> (fOutputContainer->FindObject(key));
 870   if(th1)
 871     th1->Fill(x, y) ;
 872   else
 873     AliError(Form("can not find histogram (of instance TH1) <%s> ",key)) ;
 874 }
 875
 876 //_____________________________________________________________________________
 877 void AliAnalysisTaskPi0Efficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
 878   //Fills 1D histograms with key
 879   TObject * obj = fOutputContainer->FindObject(key);
 880
 881   TH2 * th2 = dynamic_cast<TH2*> (obj);
 882   if(th2) {
 883     th2->Fill(x, y, z) ;
 884     return;
 885   }
 886
 887   TH3 * th3 = dynamic_cast<TH3*> (obj);
 888   if(th3) {
 889     th3->Fill(x, y, z) ;
 890     return;
 891   }
 892
 893   AliError(Form("can not find histogram (of instance TH2) <%s> ",key)) ;
 894 }
 895 //_____________________________________________________________________________
 896 Bool_t AliAnalysisTaskPi0Efficiency::TestLambda(Double_t pt,Double_t l1,Double_t l2){
 897
 898   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
 899   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
 900   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
 901   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
 902   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
 903   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
 904               0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
 905               0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
 906   return (R2<2.5*2.5) ;
 907
 908 }
 909 //_____________________________________________________________________________
 910 Bool_t AliAnalysisTaskPi0Efficiency::TestLambda2(Double_t pt,Double_t l1,Double_t l2){
 911
 912   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
 913   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
 914   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
 915   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
 916   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
 917   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
 918               0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
 919               0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
 920   return (R2<1.5*1.5) ;
 921
 922 }
 923 //___________________________________________________________________________
 924 void AliAnalysisTaskPi0Efficiency::ProcessMC(){
 925   //fill histograms for efficiensy etc. calculation
 926   const Double_t rcut = 1. ; //cut for primary particles
 927   //---------First pi0/eta-----------------------------
 928   char partName[10] ;
 929   char hkey[55] ;
 930
 931   AliAODEvent *event = dynamic_cast<AliAODEvent*>(InputEvent());
 932   if(!event) return ;
 933   TClonesArray *mcArray = (TClonesArray*)event->FindListObject(AliAODMCParticle::StdBranchName());
 934   for(Int_t i=0;i<mcArray->GetEntriesFast();i++){
 935      AliAODMCParticle* particle =  (AliAODMCParticle*) mcArray->At(i);
 936     if(particle->GetPdgCode() == 111)
 937       snprintf(partName,10,"pi0") ;
 938     else
 939       if(particle->GetPdgCode() == 221)
 940         snprintf(partName,10,"eta") ;
 941       else
 942         if(particle->GetPdgCode() == 22)
 943            snprintf(partName,10,"gamma") ;
 944         else
 945            continue ;
 946
 947     //Primary particle
 948     Double_t r=TMath::Sqrt(particle->Xv()*particle->Xv()+particle->Yv()*particle->Yv());
 949     if(r >rcut)
 950       continue ;
 951
 952     Double_t pt = particle->Pt() ;
 953     //Total number of pi0 with creation radius <1 cm
 954     snprintf(hkey,55,"hMC_all_%s_cen%d",partName,fCenBin) ;
 955     FillHistogram(hkey,pt) ;
 956     if(TMath::Abs(particle->Y())<0.12){
 957       snprintf(hkey,55,"hMC_unitEta_%s_cen%d",partName,fCenBin) ;
 958       FillHistogram(hkey,pt) ;
 959     }
 960
 961     snprintf(hkey,55,"hMC_rap_%s_cen%d",partName,fCenBin) ;
 962     FillHistogram(hkey,particle->Y()) ;
 963
 964     Double_t phi=particle->Phi() ;
 965     while(phi<0.)phi+=TMath::TwoPi() ;
 966     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 967     snprintf(hkey,55,"hMC_phi_%s_cen%d",partName,fCenBin) ;
 968     FillHistogram(hkey,phi) ;
 969
 970
 971     //Check if one of photons converted
 972     if(particle->GetNDaughters()!=2)
 973      continue ; //Do not account Dalitz decays
 974
 975 /*
 976     TParticle * gamma1 = fStack->Particle(particle->GetFirstDaughter());
 977     TParticle * gamma2 = fStack->Particle(particle->GetLastDaughter());
 978     //Number of pi0s decayed into acceptance
 979     Int_t mod1,mod2 ;
 980     Double_t x=0.,z=0. ;
 981     Bool_t hitPHOS1 = fPHOSGeo->ImpactOnEmc(gamma1, mod1, z,x) ;
 982     Bool_t hitPHOS2 = fPHOSGeo->ImpactOnEmc(gamma2, mod2, z,x) ;
 983
 984     Bool_t goodPair=kFALSE ;
 985     if( hitPHOS1 && hitPHOS2){
 986       sprintf(hkey,"hMC_PHOSacc_%s",partName) ;
 987       FillHistogram(hkey,pt) ;
 988       goodPair=kTRUE ;
 989     }
 990
 991 */
 992   }
 993
 994   //Now calculate "Real" distribution of clusters with primary
 995   TClonesArray cluPrim("AliCaloPhoton",200) ; //clusters with primary
 996   TClonesArray * clusters = (TClonesArray*)event->FindListObject("EmbeddedCaloClusters") ;
 997   AliAODCaloCells * cells = (AliAODCaloCells *)event->FindListObject("EmbeddedPHOScells") ;
 998   Int_t multClust = clusters->GetEntriesFast();
 999   Int_t inPHOS=0 ;
1000   Double_t vtx0[3] = {0,0,0};
1001   for (Int_t i=0; i<multClust; i++) {
1002     AliAODCaloCluster *clu = (AliAODCaloCluster*)clusters->At(i);
1003     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
1004     if(clu->GetLabel()<0) continue ;
1005
1006     Float_t  position[3];
1007     clu->GetPosition(position);
1008     TVector3 global(position) ;
1009     Int_t relId[4] ;
1010     fPHOSGeo->GlobalPos2RelId(global,relId) ;
1011     Int_t mod  = relId[0] ;
1012     Int_t cellX = relId[2];
1013     Int_t cellZ = relId[3] ;
1014     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
1015       continue ;
1016     if(clu->GetNCells()<3)
1017       continue ;
1018
1019     TLorentzVector pv1 ;
1020     clu->GetMomentum(pv1 ,vtx0);
1021
1022     if(inPHOS>=cluPrim.GetSize()){
1023       cluPrim.Expand(inPHOS+50) ;
1024     }
1025     AliCaloPhoton * ph = new(cluPrim[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
1026     //AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
1027     ph->SetModule(mod) ;
1028     AliPHOSAodCluster cluPHOS1(*clu);
1029     cluPHOS1.Recalibrate(fPHOSCalibData,cells); // modify the cell energies
1030     Double_t ecore=CoreEnergy(&cluPHOS1) ;
1031     pv1*= ecore/pv1.E() ;
1032     ph->SetMomV2(&pv1) ;
1033     ph->SetNCells(clu->GetNCells());
1034     ph->SetDispBit(TestLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
1035     ph->SetDisp2Bit(TestLambda2(clu->E(),clu->GetM20(),clu->GetM02())) ;
1036     ph->SetCPVBit(clu->GetEmcCpvDistance()>2.) ; //radius in sigmas
1037     ph->SetCPV2Bit(clu->GetEmcCpvDistance()>4.) ;
1038     ph->SetPhoton(clu->GetNExMax()<2); // Remember, if it is unfolded
1039
1040
1041     inPHOS++ ;
1042
1043   }
1044
1045   //Single photon
1046   char key[55] ;
1047   for (Int_t i1=0; i1<inPHOS; i1++) {
1048     AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
1049     snprintf(key,55,"hMCPhotAll_cen%d",fCenBin) ;
1050     FillHistogram(key,ph1->Pt()) ;
1051     snprintf(key,55,"hMCPhotAllcore_cen%d",fCenBin) ;
1052     FillHistogram(key,ph1->GetMomV2()->Pt()) ;
1053     if(ph1->IsPhoton()){
1054       snprintf(key,55,"hMCPhotAllwou_cen%d",fCenBin) ;
1055       FillHistogram(key,ph1->Pt()) ;
1056     }
1057     if(ph1->IsCPVOK() ){
1058       snprintf(key,55,"hMCPhotCPV_cen%d",fCenBin) ;
1059       FillHistogram(key,ph1->Pt()) ;
1060       snprintf(key,55,"hMCPhotCPVcore_cen%d",fCenBin) ;
1061       FillHistogram(key,ph1->GetMomV2()->Pt()) ;
1062
1063     }
1064     if(ph1->IsCPV2OK() ){
1065       snprintf(key,55,"hMCPhotCPV2_cen%d",fCenBin) ;
1066       FillHistogram(key,ph1->Pt()) ;
1067
1068     }
1069     if(ph1->IsDisp2OK()){
1070       snprintf(key,55,"hMCPhotDisp2_cen%d",fCenBin) ;
1071       FillHistogram(key,ph1->Pt()) ;
1072     }
1073     if(ph1->IsDispOK()){
1074       snprintf(key,55,"hMCPhotDisp_cen%d",fCenBin) ;
1075       FillHistogram(key,ph1->Pt()) ;
1076       if(ph1->IsPhoton()){
1077         snprintf(key,55,"hMCPhotDispwou_cen%d",fCenBin) ;
1078         FillHistogram(key,ph1->Pt()) ;
1079       }
1080       if(ph1->IsCPVOK()){
1081         snprintf(key,55,"hMCPhotBoth_cen%d",fCenBin) ;
1082         FillHistogram(key,ph1->Pt()) ;
1083         snprintf(key,55,"hMCPhotBothcore_cen%d",fCenBin) ;
1084         FillHistogram(key,ph1->GetMomV2()->Pt()) ;
1085       }
1086     } // end of loop i2
1087   } // end of loop i1
1088
1089   // Fill Real disribution
1090   for (Int_t i1=0; i1<inPHOS-1; i1++) {
1091     AliCaloPhoton * ph1=(AliCaloPhoton*)cluPrim.At(i1) ;
1092     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
1093       AliCaloPhoton * ph2=(AliCaloPhoton*)cluPrim.At(i2) ;
1094       TLorentzVector p12  = *ph1  + *ph2;
1095       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1096       Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
1097
1098       snprintf(key,55,"hMCMassPtAll_cen%d",fCenBin) ;
1099       FillHistogram(key,p12.M() ,p12.Pt()) ;
1100       snprintf(key,55,"hMCMassPtAllcore_cen%d",fCenBin) ;
1101       FillHistogram(key,pv12.M(), pv12.Pt()) ;
1102       if(ph1->IsPhoton()&& ph2->IsPhoton()){
1103         snprintf(key,55,"hMCMassPtAllwou_cen%d",fCenBin) ;
1104         FillHistogram(key,p12.M() ,p12.Pt()) ;
1105       }
1106       if(a<0.9){
1107         snprintf(key,55,"hMCMassPtAll_a09_cen%d",fCenBin) ;
1108         FillHistogram(key,p12.M() ,p12.Pt()) ;
1109         if(a<0.8){
1110           snprintf(key,55,"hMCMassPtAll_a08_cen%d",fCenBin) ;
1111           FillHistogram(key,p12.M() ,p12.Pt()) ;
1112           if(a<0.7){
1113             snprintf(key,55,"hMCMassPtAll_a07_cen%d",fCenBin) ;
1114             FillHistogram(key,p12.M() ,p12.Pt()) ;
1115           }
1116         }
1117       }
1118
1119
1120            if(ph1->Module()==1 && ph2->Module()==1)
1121             FillHistogram("hMCPi0M11",p12.M(),p12.Pt() );
1122           else if(ph1->Module()==2 && ph2->Module()==2)
1123             FillHistogram("hMCPi0M22",p12.M(),p12.Pt() );
1124           else if(ph1->Module()==3 && ph2->Module()==3)
1125             FillHistogram("hMCPi0M33",p12.M(),p12.Pt() );
1126           else if(ph1->Module()==1 && ph2->Module()==2)
1127             FillHistogram("hMCPi0M12",p12.M(),p12.Pt() );
1128           else if(ph1->Module()==1 && ph2->Module()==3)
1129             FillHistogram("hMCPi0M13",p12.M(),p12.Pt() );
1130           else if(ph1->Module()==2 && ph2->Module()==3)
1131             FillHistogram("hMCPi0M23",p12.M(),p12.Pt() );
1132
1133
1134
1135
1136
1137       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1138         snprintf(key,55,"hMCMassPtCPV_cen%d",fCenBin) ;
1139         FillHistogram(key,p12.M() ,p12.Pt()) ;
1140         snprintf(key,55,"hMCMassPtCPVcore_cen%d",fCenBin) ;
1141         FillHistogram(key,pv12.M(), pv12.Pt()) ;
1142         if(a<0.9){
1143           snprintf(key,55,"hMCMassPtCPV_a09_cen%d",fCenBin) ;
1144           FillHistogram(key,p12.M() ,p12.Pt()) ;
1145           if(a<0.8){
1146             snprintf(key,55,"hMCMassPtCPV_a08_cen%d",fCenBin) ;
1147             FillHistogram(key,p12.M() ,p12.Pt()) ;
1148             if(a<0.7){
1149               snprintf(key,55,"hMCMassPtCPV_a07_cen%d",fCenBin) ;
1150               FillHistogram(key,p12.M() ,p12.Pt()) ;
1151             }
1152           }
1153         }
1154       }
1155       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1156         snprintf(key,55,"hMCMassPtCPV2_cen%d",fCenBin) ;
1157         FillHistogram(key,p12.M() ,p12.Pt()) ;
1158         if(a<0.9){
1159           snprintf(key,55,"hMCMassPtCPV2_a09_cen%d",fCenBin) ;
1160           FillHistogram(key,p12.M() ,p12.Pt()) ;
1161           if(a<0.8){
1162             snprintf(key,55,"hMCMassPtCPV2_a08_cen%d",fCenBin) ;
1163             FillHistogram(key,p12.M() ,p12.Pt()) ;
1164             if(a<0.7){
1165               snprintf(key,55,"hMCMassPtCPV2_a07_cen%d",fCenBin) ;
1166               FillHistogram(key,p12.M() ,p12.Pt()) ;
1167             }
1168           }
1169         }
1170       }
1171       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1172         snprintf(key,55,"hMCMassPtDisp2_cen%d",fCenBin) ;
1173         FillHistogram(key,p12.M() ,p12.Pt()) ;
1174       }
1175       if(ph1->IsDispOK() && ph2->IsDispOK()){
1176         snprintf(key,55,"hMCMassPtDisp_cen%d",fCenBin) ;
1177         FillHistogram(key,p12.M() ,p12.Pt()) ;
1178         if(ph1->IsPhoton()&& ph2->IsPhoton()){
1179           snprintf(key,55,"hMCMassPtDispwou_cen%d",fCenBin) ;
1180           FillHistogram(key,p12.M() ,p12.Pt()) ;
1181         }
1182         if(a<0.9){
1183           snprintf(key,55,"hMCMassPtDisp_a09_cen%d",fCenBin) ;
1184           FillHistogram(key,p12.M() ,p12.Pt()) ;
1185           if(a<0.8){
1186             snprintf(key,55,"hMCMassPtDisp_a08_cen%d",fCenBin) ;
1187             FillHistogram(key,p12.M() ,p12.Pt()) ;
1188             if(a<0.7){
1189               snprintf(key,55,"hMCMassPtDisp_a07_cen%d",fCenBin) ;
1190               FillHistogram(key,p12.M() ,p12.Pt()) ;
1191             }
1192           }
1193         }
1194
1195         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1196           snprintf(key,55,"hMCMassPtBoth_cen%d",fCenBin) ;
1197           FillHistogram(key,p12.M() ,p12.Pt()) ;
1198           snprintf(key,55,"hMCMassPtBothcore_cen%d",fCenBin) ;
1199           FillHistogram(key,pv12.M(), pv12.Pt()) ;
1200           if(a<0.9){
1201             snprintf(key,55,"hMCMassPtBoth_a09_cen%d",fCenBin) ;
1202             FillHistogram(key,p12.M() ,p12.Pt()) ;
1203             if(a<0.8){
1204               snprintf(key,55,"hMCMassPtBoth_a08_cen%d",fCenBin) ;
1205               FillHistogram(key,p12.M() ,p12.Pt()) ;
1206               if(a<0.7){
1207                 snprintf(key,55,"hMCMassPtBoth_a07_cen%d",fCenBin) ;
1208                 FillHistogram(key,p12.M() ,p12.Pt()) ;
1209               }
1210             }
1211           }
1212         }
1213       }
1214     } // end of loop i2
1215   } // end of loop i1
1216 }
1217 //____________________________________________________________________________
1218 Double_t  AliAnalysisTaskPi0Efficiency::CoreEnergy(AliPHOSAodCluster * clu){
1219   //calculate energy of the cluster in the circle with radius distanceCut around the maximum
1220
1221   //Can not use already calculated coordinates?
1222   //They have incidence correction...
1223   const Double_t distanceCut =3.5 ;
1224   const Double_t logWeight=4.5 ;
1225
1226   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1227 // Calculates the center of gravity in the local PHOS-module coordinates
1228   Float_t wtot = 0;
1229   Double_t xc[100]={0} ;
1230   Double_t zc[100]={0} ;
1231   Double_t x = 0 ;
1232   Double_t z = 0 ;
1233   Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
1234   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1235     Int_t relid[4] ;
1236     Float_t xi ;
1237     Float_t zi ;
1238     fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
1239     fPHOSGeo->RelPosInModule(relid, xi, zi);
1240     xc[iDigit]=xi ;
1241     zc[iDigit]=zi ;
1242     if (clu->E()>0 && elist[iDigit]>0) {
1243       Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1244       x    += xc[iDigit] * w ;
1245       z    += zc[iDigit] * w ;
1246       wtot += w ;
1247     }
1248   }
1249   if (wtot>0) {
1250     x /= wtot ;
1251     z /= wtot ;
1252   }
1253   Double_t coreE=0. ;
1254   for(Int_t iDigit=0; iDigit < mulDigit; iDigit++) {
1255     Double_t distance = TMath::Sqrt((xc[iDigit]-x)*(xc[iDigit]-x)+(zc[iDigit]-z)*(zc[iDigit]-z)) ;
1256     if(distance < distanceCut)
1257       coreE += elist[iDigit] ;
1258   }
1259   //Apply non-linearity correction
1260   return (0.0241+1.0504*coreE+0.000249*coreE*coreE) ;
1261 }
1262
1263