PWG4/UserTasks/PHOS_pp_pi0/AliAnalysisTaskPi0.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TObjArray.h"
   4 #include "TF1.h"
   5 #include "TH1F.h"
   6 #include "TH2F.h"
   7 #include "TH2I.h"
   8 #include "TH3F.h"
   9 #include "TParticle.h"
  10 #include "TCanvas.h"
  11 #include "TStyle.h"
  12
  13 #include "AliAnalysisTaskSE.h"
  14 #include "AliAnalysisTaskPi0.h"
  15 #include "AliCaloPhoton.h"
  16 #include "AliPHOSGeometry.h"
  17 #include "AliESDEvent.h"
  18 #include "AliESDCaloCells.h"
  19 #include "AliESDCaloCluster.h"
  20 #include "AliESDVertex.h"
  21 #include "AliESDtrackCuts.h"
  22 #include "AliLog.h"
  23 #include "AliPID.h"
  24 #include "AliTriggerAnalysis.h"
  25
  26 // Analysis task to fill histograms with PHOS ESD clusters and cells
  27 // Authors: Yuri Kharlov
  28 // Date   : 28.05.2009
  29
  30 ClassImp(AliAnalysisTaskPi0)
  31
  32 //________________________________________________________________________
  33 AliAnalysisTaskPi0::AliAnalysisTaskPi0(const char *name)
  34 : AliAnalysisTaskSE(name),
  35   fESDtrackCuts(0),
  36   fOutputContainer(0),
  37   fPHOSEvent(0),
  38   fnCINT1B(0),
  39   fnCINT1A(0),
  40   fnCINT1C(0),
  41   fnCINT1E(0),
  42   fPHOSGeo(0),
  43   fEventCounter(0),
  44   fTriggerAnalysis(new AliTriggerAnalysis)
  45 {
  46   // Constructor
  47   Int_t nBin=10 ;
  48   for(Int_t i=0;i<nBin;i++){
  49     for(Int_t j=0;j<2;j++)
  50       fPHOSEvents[i][j]=0 ;
  51   }
  52
  53   // Output slots #0 write into a TH1 container
  54   DefineOutput(1,TList::Class());
  55
  56   // Set bad channel map
  57   char key[55] ;
  58   for(Int_t i=0; i<6; i++){
  59     snprintf(key,55,"PHOS_BadMap_mod%d",i) ;
  60     fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
  61   }
  62   // Initialize the PHOS geometry
  63   fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
  64
  65 }
  66
  67 //________________________________________________________________________
  68 void AliAnalysisTaskPi0::UserCreateOutputObjects()
  69 {
  70   // Create histograms
  71   // Called once
  72
  73   // ESD histograms
  74   if(fOutputContainer != NULL){
  75     delete fOutputContainer;
  76   }
  77   fOutputContainer = new TList();
  78   fOutputContainer->SetOwner(kTRUE);
  79
  80   fOutputContainer->Add(new TH1I("hCellMultEvent"  ,"PHOS cell multiplicity per event"    ,2000,0,2000));
  81   fOutputContainer->Add(new TH1I("hCellMultEventM1","PHOS cell multiplicity per event, M1",2000,0,2000));
  82   fOutputContainer->Add(new TH1I("hCellMultEventM2","PHOS cell multiplicity per event, M2",2000,0,2000));
  83   fOutputContainer->Add(new TH1I("hCellMultEventM3","PHOS cell multiplicity per event, M3",2000,0,2000));
  84   fOutputContainer->Add(new TH1I("hClusterMult"      ,"CaloCluster multiplicity"     ,100,0,100));
  85   fOutputContainer->Add(new TH1I("hPHOSClusterMult"  ,"PHOS cluster multiplicity"    ,100,0,100));
  86   fOutputContainer->Add(new TH1I("hPHOSClusterMultM1","PHOS cluster multiplicity, M1",100,0,100));
  87   fOutputContainer->Add(new TH1I("hPHOSClusterMultM2","PHOS cluster multiplicity, M2",100,0,100));
  88   fOutputContainer->Add(new TH1I("hPHOSClusterMultM3","PHOS cluster multiplicity, M3",100,0,100));
  89   fOutputContainer->Add(new TH1F("hCellEnergy"  ,"Cell energy"            ,3000,0.,30.));
  90   fOutputContainer->Add(new TH1F("hCellEnergyM1","Cell energy in module 1",3000,0.,30.));
  91   fOutputContainer->Add(new TH1F("hCellEnergyM2","Cell energy in module 2",3000,0.,30.));
  92   fOutputContainer->Add(new TH1F("hCellEnergyM3","Cell energy in module 3",3000,0.,30.));
  93   fOutputContainer->Add(new TH1F("hClusterEnergy"  ,"Cluster energy"      ,3000,0.,30.));
  94   fOutputContainer->Add(new TH1F("hClusterEnergyM1","Cluster energy, M1"  ,3000,0.,30.));
  95   fOutputContainer->Add(new TH1F("hClusterEnergyM2","Cluster energy, M2"  ,3000,0.,30.));
  96   fOutputContainer->Add(new TH1F("hClusterEnergyM3","Cluster energy, M3"  ,3000,0.,30.));
  97   fOutputContainer->Add(new TH2F("hClusterEvsN"  ,"Cluster energy vs digit multiplicity"    ,3000,0.,30.,40,0.,40.));
  98   fOutputContainer->Add(new TH2F("hClusterEvsNM1","Cluster energy vs digit multiplicity, M1",3000,0.,30.,40,0.,40.));
  99   fOutputContainer->Add(new TH2F("hClusterEvsNM2","Cluster energy vs digit multiplicity, M2",3000,0.,30.,40,0.,40.));
 100   fOutputContainer->Add(new TH2F("hClusterEvsNM3","Cluster energy vs digit multiplicity, M3",3000,0.,30.,40,0.,40.));
 101   fOutputContainer->Add(new TH1I("hCellMultClu"  ,"Cell multiplicity per cluster"    ,200,0,200));
 102   fOutputContainer->Add(new TH1I("hCellMultCluM1","Cell multiplicity per cluster, M1",200,0,200));
 103   fOutputContainer->Add(new TH1I("hCellMultCluM2","Cell multiplicity per cluster, M3",200,0,200));
 104   fOutputContainer->Add(new TH1I("hCellMultCluM3","Cell multiplicity per cluster, M3",200,0,200));
 105   fOutputContainer->Add(new TH1I("hModule","Module events",5,0.,5.));
 106   fOutputContainer->Add(new TH1F("hSelEvents","Selected events",7,0.5,7.5));
 107
 108   fOutputContainer->Add(new TH2F("hCellNXZM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 109   fOutputContainer->Add(new TH2F("hCellNXZM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 110   fOutputContainer->Add(new TH2F("hCellNXZM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 111   fOutputContainer->Add(new TH2F("hCellEXZM1","Cell E(X,Z), M1",64,0.5,64.5, 56,0.5,56.5));
 112   fOutputContainer->Add(new TH2F("hCellEXZM2","Cell E(X,Z), M2",64,0.5,64.5, 56,0.5,56.5));
 113   fOutputContainer->Add(new TH2F("hCellEXZM3","Cell E(X,Z), M3",64,0.5,64.5, 56,0.5,56.5));
 114   fOutputContainer->Add(new TH2F("hCluNXZM1","Clu (X,Z), M1"   ,64,0.5,64.5, 56,0.5,56.5));
 115   fOutputContainer->Add(new TH2F("hCluNXZM2","Clu (X,Z), M2"   ,64,0.5,64.5, 56,0.5,56.5));
 116   fOutputContainer->Add(new TH2F("hCluNXZM3","Clu (X,Z), M3"   ,64,0.5,64.5, 56,0.5,56.5));
 117   fOutputContainer->Add(new TH2F("hCluEXZM1","Clu E(X,Z), M1"  ,64,0.5,64.5, 56,0.5,56.5));
 118   fOutputContainer->Add(new TH2F("hCluEXZM2","Clu E(X,Z), M2"  ,64,0.5,64.5, 56,0.5,56.5));
 119   fOutputContainer->Add(new TH2F("hCluEXZM3","Clu E(X,Z), M3"  ,64,0.5,64.5, 56,0.5,56.5));
 120
 121   Int_t nM       = 750;
 122   Double_t mMin  = 0.0;
 123   Double_t mMax  = 1.5;
 124   Int_t nPt      = 400;
 125   Double_t ptMin = 0;
 126   Double_t ptMax = 40;
 127   fOutputContainer->Add(new TH2F("hAsymPtPi0","(A,p_{T})_{#gamma#gamma} #pi^{0}"     ,50,0.,1.,    nPt,ptMin,ptMax));
 128   fOutputContainer->Add(new TH2F("hAsymPtEta","(A,p_{T})_{#gamma#gamma} #eta"        ,50,0.,1.,    nPt,ptMin,ptMax));
 129
 130   fOutputContainer->Add(new TH2F("hMassPtA10" ,"(M,p_{T})_{#gamma#gamma}, 0<A<1.0"   ,nM,mMin,mMax,nPt,ptMin,ptMax));
 131   fOutputContainer->Add(new TH2F("hMassPtA08" ,"(M,p_{T})_{#gamma#gamma}, 0<A<0.8"   ,nM,mMin,mMax,nPt,ptMin,ptMax));
 132   fOutputContainer->Add(new TH2F("hMassPtA07" ,"(M,p_{T})_{#gamma#gamma}, 0<A<0.7"   ,nM,mMin,mMax,nPt,ptMin,ptMax));
 133   fOutputContainer->Add(new TH2F("hMassPtA01" ,"(M,p_{T})_{#gamma#gamma}, 0<A<0.1"   ,nM,mMin,mMax,nPt,ptMin,ptMax));
 134
 135   fOutputContainer->Add(new TH2F("hMassPtA10nvtx" ,"(M,p_{T})_{#gamma#gamma}, 0<A<1.0, no vtx" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 136   fOutputContainer->Add(new TH2F("hMassPtA07nvtx" ,"(M,p_{T})_{#gamma#gamma}, 0<A<0.7, no vtx" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 137
 138   fOutputContainer->Add(new TH2F("hMassPtA10vtx" ,"(M,p_{T})_{#gamma#gamma}, 0<A<1.0, vtx"     ,nM,mMin,mMax,nPt,ptMin,ptMax));
 139   fOutputContainer->Add(new TH2F("hMassPtA07vtx" ,"(M,p_{T})_{#gamma#gamma}, 0<A<0.7, vtx"     ,nM,mMin,mMax,nPt,ptMin,ptMax));
 140
 141   fOutputContainer->Add(new TH2F("hMassPtA10V0AND" ,"(M,p_{T})_{#gamma#gamma}, 0<A<1.0, V0AND",nM,mMin,mMax,nPt,ptMin,ptMax));
 142   fOutputContainer->Add(new TH2F("hMassPtA07V0AND" ,"(M,p_{T})_{#gamma#gamma}, 0<A<0.7, V0AND",nM,mMin,mMax,nPt,ptMin,ptMax));
 143
 144   fOutputContainer->Add(new TH2F("hMassPtA10V0PU" ,"(M,p_{T})_{#gamma#gamma}, 0<A<1.0, pileup",nM,mMin,mMax,nPt,ptMin,ptMax));
 145   fOutputContainer->Add(new TH2F("hMassPtA07V0PU" ,"(M,p_{T})_{#gamma#gamma}, 0<A<0.7, pileup",nM,mMin,mMax,nPt,ptMin,ptMax));
 146
 147   fOutputContainer->Add(new TH2F("hMassPtvA10","(M,p_{T})_{#gamma#gamma}, 0<A<1.0, ESD vertex",nM,mMin,mMax,nPt,ptMin,ptMax));
 148   fOutputContainer->Add(new TH2F("hMassPtvA07","(M,p_{T})_{#gamma#gamma}, 0<A<0.7, ESD vertex",nM,mMin,mMax,nPt,ptMin,ptMax));
 149
 150   fOutputContainer->Add(new TH3F("hMassPtCA10","(M,p_{T},C)_{#gamma#gamma}, 0<A<1.0" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 151   fOutputContainer->Add(new TH3F("hMassPtCA10_cpv","(M,p_{T},C)_{#gamma#gamma}, 0<A<1.0" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 152   fOutputContainer->Add(new TH3F("hMassPtCA10_disp","(M,p_{T},C)_{#gamma#gamma}, 0<A<1.0" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 153   fOutputContainer->Add(new TH3F("hMassPtCA10_both","(M,p_{T},C)_{#gamma#gamma}, 0<A<1.0" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 154   fOutputContainer->Add(new TH3F("hMassPtCA07","(M,p_{T},C)_{#gamma#gamma}, 0<A<0.7" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 155   fOutputContainer->Add(new TH3F("hMassPtCA07_cpv","(M,p_{T},C)_{#gamma#gamma}, 0<A<1.0" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 156   fOutputContainer->Add(new TH3F("hMassPtCA07_disp","(M,p_{T},C)_{#gamma#gamma}, 0<A<1.0" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 157   fOutputContainer->Add(new TH3F("hMassPtCA07_both","(M,p_{T},C)_{#gamma#gamma}, 0<A<1.0" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 158
 159   fOutputContainer->Add(new TH2F("hMassSingle_all","(M,p_{T})_{#gamma#gamma}, no PID" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 160   fOutputContainer->Add(new TH2F("hMassSingle_cpv","(M,p_{T})_{#gamma#gamma}, no PID" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 161   fOutputContainer->Add(new TH2F("hMassSingle_disp","(M,p_{T})_{#gamma#gamma}, no PID" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 162   fOutputContainer->Add(new TH2F("hMassSingle_both","(M,p_{T})_{#gamma#gamma}, no PID" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 163
 164   fOutputContainer->Add(new TH2F("hMassPtM1","(M,p_{T})_{#gamma#gamma}, module 1"    ,nM,mMin,mMax,nPt,ptMin,ptMax));
 165   fOutputContainer->Add(new TH2F("hMassPtM2","(M,p_{T})_{#gamma#gamma}, module 2"    ,nM,mMin,mMax,nPt,ptMin,ptMax));
 166   fOutputContainer->Add(new TH2F("hMassPtM3","(M,p_{T})_{#gamma#gamma}, module 3"    ,nM,mMin,mMax,nPt,ptMin,ptMax));
 167   fOutputContainer->Add(new TH2F("hMassPtM12","(M,p_{T})_{#gamma#gamma}, modules 1,2",nM,mMin,mMax,nPt,ptMin,ptMax));
 168   fOutputContainer->Add(new TH2F("hMassPtM13","(M,p_{T})_{#gamma#gamma}, modules 1,3",nM,mMin,mMax,nPt,ptMin,ptMax));
 169   fOutputContainer->Add(new TH2F("hMassPtM23","(M,p_{T})_{#gamma#gamma}, modules 2,3",nM,mMin,mMax,nPt,ptMin,ptMax));
 170
 171   fOutputContainer->Add(new TH2F("hMassPtN3","(M,p_{T})_{#gamma#gamma}, N_{cell}>2"  ,nM,mMin,mMax,nPt,ptMin,ptMax));
 172   fOutputContainer->Add(new TH2F("hMassPtN4","(M,p_{T})_{#gamma#gamma}, N_{cell}>3"  ,nM,mMin,mMax,nPt,ptMin,ptMax));
 173   fOutputContainer->Add(new TH2F("hMassPtN5","(M,p_{T})_{#gamma#gamma}, N_{cell}>4"  ,nM,mMin,mMax,nPt,ptMin,ptMax));
 174   fOutputContainer->Add(new TH2F("hMassPtN6","(M,p_{T})_{#gamma#gamma}, N_{cell}>5"  ,nM,mMin,mMax,nPt,ptMin,ptMax));
 175
 176   fOutputContainer->Add(new TH2F("hMiAsymPt","(A,p_{T})_{#gamma#gamma}"                ,50,0.,1.,    nPt,ptMin,ptMax));
 177
 178   fOutputContainer->Add(new TH2F("hMiMassPtA10" ,"(M,p_{T})_{#gamma#gamma}, 0<A<1.0"   ,nM,mMin,mMax,nPt,ptMin,ptMax));
 179   fOutputContainer->Add(new TH2F("hMiMassPtA08" ,"(M,p_{T})_{#gamma#gamma}, 0<A<0.8"   ,nM,mMin,mMax,nPt,ptMin,ptMax));
 180   fOutputContainer->Add(new TH2F("hMiMassPtA07" ,"(M,p_{T})_{#gamma#gamma}, 0<A<0.7"   ,nM,mMin,mMax,nPt,ptMin,ptMax));
 181   fOutputContainer->Add(new TH2F("hMiMassPtA01" ,"(M,p_{T})_{#gamma#gamma}, 0<A<0.1"   ,nM,mMin,mMax,nPt,ptMin,ptMax));
 182
 183   fOutputContainer->Add(new TH2F("hMiMassPtA10nvtx" ,"(M,p_{T})_{#gamma#gamma}, 0<A<1.0, no vtx" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 184   fOutputContainer->Add(new TH2F("hMiMassPtA07nvtx" ,"(M,p_{T})_{#gamma#gamma}, 0<A<0.7, no vtx" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 185
 186   fOutputContainer->Add(new TH2F("hMiMassPtA10vtx" ,"(M,p_{T})_{#gamma#gamma}, 0<A<1.0, vtx"     ,nM,mMin,mMax,nPt,ptMin,ptMax));
 187   fOutputContainer->Add(new TH2F("hMiMassPtA07vtx" ,"(M,p_{T})_{#gamma#gamma}, 0<A<0.7, vtx"     ,nM,mMin,mMax,nPt,ptMin,ptMax));
 188
 189   fOutputContainer->Add(new TH2F("hMiMassPtA10V0AND" ,"(M,p_{T})_{#gamma#gamma}, 0<A<1.0, V0AND",nM,mMin,mMax,nPt,ptMin,ptMax));
 190   fOutputContainer->Add(new TH2F("hMiMassPtA07V0AND" ,"(M,p_{T})_{#gamma#gamma}, 0<A<0.7, V0AND",nM,mMin,mMax,nPt,ptMin,ptMax));
 191
 192   fOutputContainer->Add(new TH2F("hMiMassPtA10V0PU" ,"(M,p_{T})_{#gamma#gamma}, 0<A<1.0, pileup",nM,mMin,mMax,nPt,ptMin,ptMax));
 193   fOutputContainer->Add(new TH2F("hMiMassPtA07V0PU" ,"(M,p_{T})_{#gamma#gamma}, 0<A<0.7, pileup",nM,mMin,mMax,nPt,ptMin,ptMax));
 194
 195   fOutputContainer->Add(new TH2F("hMiMassPtvA10","(M,p_{T})_{#gamma#gamma}, 0<A<1.0, ESD vertex",nM,mMin,mMax,nPt,ptMin,ptMax));
 196   fOutputContainer->Add(new TH2F("hMiMassPtvA07","(M,p_{T})_{#gamma#gamma}, 0<A<0.7, ESD vertex",nM,mMin,mMax,nPt,ptMin,ptMax));
 197
 198   fOutputContainer->Add(new TH3F("hMiMassPtCA10","(M,p_{T},C)_{#gamma#gamma}, 0<A<1.0" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 199   fOutputContainer->Add(new TH3F("hMiMassPtCA10_cpv","(M,p_{T},C)_{#gamma#gamma}, 0<A<1.0" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 200   fOutputContainer->Add(new TH3F("hMiMassPtCA10_disp","(M,p_{T},C)_{#gamma#gamma}, 0<A<1.0" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 201   fOutputContainer->Add(new TH3F("hMiMassPtCA10_both","(M,p_{T},C)_{#gamma#gamma}, 0<A<1.0" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 202   fOutputContainer->Add(new TH3F("hMiMassPtCA07","(M,p_{T},C)_{#gamma#gamma}, 0<A<0.7" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 203   fOutputContainer->Add(new TH3F("hMiMassPtCA07_cpv","(M,p_{T},C)_{#gamma#gamma}, 0<A<1.0" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 204   fOutputContainer->Add(new TH3F("hMiMassPtCA07_disp","(M,p_{T},C)_{#gamma#gamma}, 0<A<1.0" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 205   fOutputContainer->Add(new TH3F("hMiMassPtCA07_both","(M,p_{T},C)_{#gamma#gamma}, 0<A<1.0" ,nM,mMin,mMax,nPt,ptMin,ptMax,8,0.,8.));
 206
 207   fOutputContainer->Add(new TH2F("hMiMassSingle_all","(M,p_{T})_{#gamma#gamma}, no PID" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 208   fOutputContainer->Add(new TH2F("hMiMassSingle_cpv","(M,p_{T})_{#gamma#gamma}, no PID" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 209   fOutputContainer->Add(new TH2F("hMiMassSingle_disp","(M,p_{T})_{#gamma#gamma}, no PID" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 210   fOutputContainer->Add(new TH2F("hMiMassSingle_both","(M,p_{T})_{#gamma#gamma}, no PID" ,nM,mMin,mMax,nPt,ptMin,ptMax));
 211
 212   fOutputContainer->Add(new TH2F("hMiMassPtM1","(M,p_{T})_{#gamma#gamma}, module 1"    ,nM,mMin,mMax,nPt,ptMin,ptMax));
 213   fOutputContainer->Add(new TH2F("hMiMassPtM2","(M,p_{T})_{#gamma#gamma}, module 2"    ,nM,mMin,mMax,nPt,ptMin,ptMax));
 214   fOutputContainer->Add(new TH2F("hMiMassPtM3","(M,p_{T})_{#gamma#gamma}, module 3"    ,nM,mMin,mMax,nPt,ptMin,ptMax));
 215   fOutputContainer->Add(new TH2F("hMiMassPtM12","(M,p_{T})_{#gamma#gamma}, modules 1,2",nM,mMin,mMax,nPt,ptMin,ptMax));
 216   fOutputContainer->Add(new TH2F("hMiMassPtM13","(M,p_{T})_{#gamma#gamma}, modules 1,3",nM,mMin,mMax,nPt,ptMin,ptMax));
 217   fOutputContainer->Add(new TH2F("hMiMassPtM23","(M,p_{T})_{#gamma#gamma}, modules 2,3",nM,mMin,mMax,nPt,ptMin,ptMax));
 218
 219   fOutputContainer->Add(new TH2F("hMiMassPtN3","(M,p_{T})_{#gamma#gamma}, N_{cell}>2"  ,nM,mMin,mMax,nPt,ptMin,ptMax));
 220   fOutputContainer->Add(new TH2F("hMiMassPtN4","(M,p_{T})_{#gamma#gamma}, N_{cell}>3"  ,nM,mMin,mMax,nPt,ptMin,ptMax));
 221   fOutputContainer->Add(new TH2F("hMiMassPtN5","(M,p_{T})_{#gamma#gamma}, N_{cell}>4"  ,nM,mMin,mMax,nPt,ptMin,ptMax));
 222   fOutputContainer->Add(new TH2F("hMiMassPtN6","(M,p_{T})_{#gamma#gamma}, N_{cell}>5"  ,nM,mMin,mMax,nPt,ptMin,ptMax));
 223
 224   fOutputContainer->Add(new TH1F("hPhotonKappa","#kappa(#gamma)",200,0.,20.));
 225   fOutputContainer->Add(new TH1F("hPhotonPt","p_{T}(#gamma)",200,0.,20.));
 226   fOutputContainer->Add(new TH1F("hPhotonPx","p_{x}(#gamma)",200,0.,20.));
 227   fOutputContainer->Add(new TH1F("hPhotonPy","p_{y}(#gamma)",200,0.,20.));
 228
 229   fOutputContainer->Add(new TH1F("hTrigClass","Trigger class",5,0.5,5.5));
 230
 231   fOutputContainer->Add(new TH1F("hZvertex","Z vertex",200,-50.,+50.));
 232   fOutputContainer->Add(new TH1F("hNvertexTracks","N of primary tracks from the primary vertex",150,0.,150.));
 233   fOutputContainer->Add(new TH1F("hTrackMult","Charged track multiplicity",150,0.,150.));
 234
 235   // Create ESD track cut
 236
 237   fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
 238   fESDtrackCuts ->SetMaxDCAToVertexZ(2);
 239   fESDtrackCuts ->SetEtaRange(-0.8, 0.8);
 240   fESDtrackCuts ->SetPtRange(0.15);
 241
 242   PostData(1, fOutputContainer);
 243
 244 }
 245
 246 //________________________________________________________________________
 247 void AliAnalysisTaskPi0::UserExec(Option_t *)
 248 {
 249   // Main loop, called for each event
 250   // Analyze ESD
 251
 252   // Event selection flags
 253
 254   Bool_t eventVtxExist    = kFALSE;
 255   Bool_t eventVtxZ10cm    = kFALSE;
 256   Bool_t eventPileup      = kFALSE;
 257   Bool_t eventV0AND       = kFALSE;
 258
 259   AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
 260   if (!event) {
 261      Printf("ERROR: Could not retrieve event");
 262      return;
 263   }
 264
 265   Int_t eventNumberInFile = event->GetEventNumberInFile();
 266   if(fPHOSEvent)
 267     fPHOSEvent->Clear() ;
 268   else
 269     fPHOSEvent = new TClonesArray("AliCaloPhoton",50) ;
 270
 271   // Checks if we have a primary vertex
 272   // Get primary vertices form ESD
 273
 274   if      (event->GetPrimaryVertexTracks()->GetNContributors()>0)
 275     eventVtxExist    = kTRUE;
 276   else if (event->GetPrimaryVertexSPD()   ->GetNContributors()>0)
 277     eventVtxExist    = kTRUE;
 278
 279   const AliESDVertex *esdVertex5 = event->GetPrimaryVertex();
 280
 281   Double_t vtx0[3] = {0,0,0}; // don't rely on ESD vertex, assume (0,0,0)
 282   Double_t vtx5[3];
 283   vtx5[0] = esdVertex5->GetX();
 284   vtx5[1] = esdVertex5->GetY();
 285   vtx5[2] = esdVertex5->GetZ();
 286
 287   FillHistogram("hNvertexTracks",esdVertex5->GetNContributors());
 288   FillHistogram("hZvertex"      ,esdVertex5->GetZ());
 289   if (TMath::Abs(esdVertex5->GetZ()) < 10. )
 290     eventVtxZ10cm = kTRUE;
 291
 292   if (event->IsPileupFromSPD())
 293     eventPileup = kTRUE;
 294
 295   eventV0AND = fTriggerAnalysis->IsOfflineTriggerFired(event, AliTriggerAnalysis::kV0AND);
 296
 297   // Fill event statistics for different selection criteria
 298
 299   FillHistogram("hSelEvents",1) ;
 300   if (eventVtxExist)
 301     FillHistogram("hSelEvents",2) ;
 302   if (eventVtxExist && eventVtxZ10cm)
 303     FillHistogram("hSelEvents",3) ;
 304   if (eventVtxExist && eventVtxZ10cm && eventV0AND)
 305     FillHistogram("hSelEvents",4) ;
 306   if (eventVtxExist && eventVtxZ10cm && eventV0AND && eventPileup)
 307     FillHistogram("hSelEvents",5) ;
 308   if (eventPileup)
 309     FillHistogram("hSelEvents",6) ;
 310   if(eventV0AND){
 311     FillHistogram("hSelEvents",7) ;
 312   }
 313
 314   //Vtx class z-bin
 315   Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
 316   if(zvtx<0)zvtx=0 ;
 317   if(zvtx>9)zvtx=9 ;
 318
 319   TString trigClasses = event->GetFiredTriggerClasses();
 320
 321   if (trigClasses.Contains("CINT1B")) fnCINT1B++;
 322   if (trigClasses.Contains("CINT1A")) fnCINT1A++;
 323   if (trigClasses.Contains("CINT1C")) fnCINT1C++;
 324   if (trigClasses.Contains("CINT1-E")) fnCINT1E++;
 325
 326   //Calculate charged multiplicity
 327   Int_t trackMult = 0;
 328   for (Int_t i=0;i<event->GetNumberOfTracks();++i) {
 329     AliESDtrack *track = new AliESDtrack(*event->GetTrack(i)) ;
 330     if(fESDtrackCuts->AcceptTrack(track) &&  TMath::Abs(track->Eta())< 0.8)
 331       trackMult++;
 332     delete track;
 333   }
 334   FillHistogram("hTrackMult",trackMult+0.5) ;
 335
 336   Int_t centr=0 ;
 337   //always zero centrality
 338   if(!fPHOSEvents[zvtx][centr]) fPHOSEvents[zvtx][centr]=new TList() ;
 339   TList * prevPHOS = fPHOSEvents[zvtx][centr] ;
 340
 341   if(trackMult<=2)
 342     centr=0 ;
 343   else
 344     if(trackMult<=5)
 345       centr=1 ;
 346     else
 347       if(trackMult<=9)
 348         centr=2 ;
 349       else
 350         if(trackMult<=14)
 351           centr=3 ;
 352         else
 353           if(trackMult<=22)
 354             centr=4 ;
 355           else
 356             if(trackMult<=35)
 357               centr=5 ;
 358             else
 359               if(trackMult<=50)
 360                 centr=6 ;
 361               else
 362                 centr=7 ;
 363
 364
 365   AliESDCaloCluster *clu1;
 366   TLorentzVector p1,p2,p12, pv1,pv2,pv12;
 367   AliESDCaloCells *cells      = event->GetPHOSCells();
 368
 369   Int_t multClust = event->GetNumberOfCaloClusters();
 370   Int_t multCells = cells->GetNumberOfCells();
 371   FillHistogram("hClusterMult",multClust);
 372   FillHistogram("hCellMultEvent",multCells);
 373
 374   Printf("Event %d, trig.class %s, period %d, bc %d, orbit %d",
 375      eventNumberInFile,trigClasses.Data(),event->GetPeriodNumber(),
 376      event->GetBunchCrossNumber(),event->GetOrbitNumber());
 377   Printf("\tthere are %d caloclusters and %d calocells",
 378      multClust,multCells);
 379
 380   // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
 381   if(fEventCounter == 0) {
 382     for(Int_t mod=0; mod<5; mod++) {
 383       if(!event->GetPHOSMatrix(mod)) continue;
 384       fPHOSGeo->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
 385       Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
 386     }
 387   }
 388
 389   Float_t  energy;
 390   Int_t    mod1, relId[4], cellAbsId, cellX, cellZ;
 391
 392   // Single loop over cells
 393
 394   Int_t nCellModule[3] = {0,0,0};
 395   for (Int_t iCell=0; iCell<multCells; iCell++) {
 396     cellAbsId = cells->GetCellNumber(iCell);
 397     fPHOSGeo->AbsToRelNumbering(cellAbsId,relId);
 398     mod1  = relId[0];
 399     cellX = relId[2];
 400     cellZ = relId[3] ;
 401     energy = cells->GetAmplitude(iCell);
 402     FillHistogram("hCellEnergy",energy);
 403     if      (mod1==1) {
 404       nCellModule[0]++;
 405       FillHistogram("hCellEnergyM1",cells->GetAmplitude(iCell));
 406       FillHistogram("hCellNXZM1",cellX,cellZ,1.);
 407       FillHistogram("hCellEXZM1",cellX,cellZ,energy);
 408     }
 409     else if (mod1==2) {
 410       nCellModule[1]++;
 411       FillHistogram("hCellEnergyM2",cells->GetAmplitude(iCell));
 412       FillHistogram("hCellNXZM2",cellX,cellZ,1.);
 413       FillHistogram("hCellEXZM2",cellX,cellZ,energy);
 414     }
 415     else if (mod1==3) {
 416       nCellModule[2]++;
 417       FillHistogram("hCellEnergyM3",cells->GetAmplitude(iCell));
 418       FillHistogram("hCellNXZM3",cellX,cellZ,1.);
 419       FillHistogram("hCellEXZM3",cellX,cellZ,energy);
 420     }
 421   }
 422   FillHistogram("hCellMultEventM1",nCellModule[0]);
 423   FillHistogram("hCellMultEventM2",nCellModule[1]);
 424   FillHistogram("hCellMultEventM3",nCellModule[2]);
 425
 426   // Single loop over clusters fills cluster histograms
 427
 428   Int_t    digMult;
 429   Int_t    multPHOSClust[4]  = {0,0,0,0};
 430   Float_t  position[3];
 431
 432   for (Int_t i1=0; i1<multClust; i1++) {
 433     clu1 = event->GetCaloCluster(i1);
 434     if ( !clu1->IsPHOS() ) continue;
 435
 436     digMult = clu1->GetNCells();
 437     clu1->GetPosition(position);
 438     TVector3 global1(position) ;
 439     fPHOSGeo->GlobalPos2RelId(global1,relId) ;
 440     mod1  = relId[0] ;
 441     cellX = relId[2];
 442     cellZ = relId[3] ;
 443     if ( !IsGoodChannel("PHOS",mod1,cellX,cellZ) ) continue ;
 444
 445     cellAbsId = clu1->GetCellAbsId(0);
 446     fPHOSGeo->AbsToRelNumbering(cellAbsId,relId);
 447     mod1   = relId[0];
 448     energy = clu1->E();
 449
 450     multPHOSClust[0]++;
 451     FillHistogram("hClusterEnergy",energy);
 452     FillHistogram("hClusterEvsN",energy,digMult);
 453     FillHistogram("hCellMultClu",digMult);
 454     if      (mod1==1) {
 455       multPHOSClust[1]++;
 456       FillHistogram("hClusterEvsNM1",energy,digMult);
 457       FillHistogram("hCellMultCluM1",digMult);
 458       FillHistogram("hClusterEnergyM1",energy);
 459       if (digMult == 1 && energy > 3.) {
 460     FillHistogram("hCluNXZM1",cellX,cellZ,1.);
 461     FillHistogram("hCluEXZM1",cellX,cellZ,energy);
 462       }
 463     }
 464     else if (mod1==2) {
 465       multPHOSClust[2]++;
 466       FillHistogram("hClusterEvsNM2",energy,digMult);
 467       FillHistogram("hCellMultCluM2",digMult);
 468       FillHistogram("hClusterEnergyM2",energy);
 469       if (digMult == 1 && energy > 3.) {
 470     FillHistogram("hCluNXZM2",cellX,cellZ,1.);
 471     FillHistogram("hCluEXZM2",cellX,cellZ,energy);
 472       }
 473     }
 474     else if (mod1==3) {
 475       multPHOSClust[3]++;
 476       FillHistogram("hClusterEvsNM3",energy,digMult);
 477       FillHistogram("hCellMultCluM3",digMult);
 478       FillHistogram("hClusterEnergyM3",energy);
 479       if (digMult == 1 && energy > 3.) {
 480     FillHistogram("hCluNXZM3",cellX,cellZ,1.);
 481     FillHistogram("hCluEXZM3",cellX,cellZ,energy);
 482       }
 483     }
 484
 485     if (digMult > 2) {
 486       clu1 ->GetMomentum(p1 ,vtx0);
 487       Double_t pAbs = p1.P();
 488       Double_t pT   = p1.Pt();
 489       Double_t pX   = p1.Px();
 490       Double_t pY   = p1.Py();
 491       Double_t kappa = pAbs - TMath::Power(0.135,2)/4./pAbs;
 492
 493       FillHistogram("hPhotonKappa",kappa);
 494       FillHistogram("hPhotonPt",pT);
 495       FillHistogram("hPhotonPx",pX);
 496       FillHistogram("hPhotonPy",pY);
 497     }
 498   }
 499   FillHistogram("hPHOSClusterMult",multPHOSClust[0]);
 500   FillHistogram("hPHOSClusterMultM1",multPHOSClust[1]);
 501   FillHistogram("hPHOSClusterMultM2",multPHOSClust[2]);
 502   FillHistogram("hPHOSClusterMultM3",multPHOSClust[3]);
 503
 504   //Select photons for inv mass calculation
 505   Int_t inPHOS=0 ;
 506   for (Int_t i1=0; i1<multClust; i1++) {
 507     clu1 = event->GetCaloCluster(i1);
 508     if ( !clu1->IsPHOS() || clu1->E()<0.3) continue;
 509
 510     clu1->GetPosition(position);
 511     TVector3 global1(position) ;
 512     fPHOSGeo->GlobalPos2RelId(global1,relId) ;
 513     mod1  = relId[0] ;
 514     cellX = relId[2];
 515     cellZ = relId[3] ;
 516     if ( !IsGoodChannel("PHOS",mod1,cellX,cellZ) ) continue ;
 517
 518     clu1 ->GetMomentum(p1 ,vtx0);
 519     clu1 ->GetMomentum(pv1,vtx5);
 520     digMult   = clu1->GetNCells();
 521     new((*fPHOSEvent)[inPHOS]) AliCaloPhoton(p1.X(),p1.Py(),p1.Z(),p1.E()) ;
 522     AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
 523     ph->SetModule(mod1) ;
 524     ph->SetMomV2(&pv1) ;
 525     ph->SetNCells(clu1->GetNCells());
 526     ph->SetDispBit(TestLambda(clu1->GetM20(),clu1->GetM02())) ;
 527     ph->SetCPVBit(clu1->GetEmcCpvDistance()>10.) ;
 528
 529     inPHOS++ ;
 530   }
 531
 532   // Fill Real disribution
 533
 534   for (Int_t i1=0; i1<inPHOS-1; i1++) {
 535     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 536     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
 537       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent->At(i2) ;
 538       p12  = *ph1  + *ph2;
 539       pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 540       Double_t asym  = TMath::Abs((ph1->Energy()-ph2->Energy())/(ph1->Energy()+ph2->Energy()));
 541
 542       if (ph1->GetNCells()>2 && ph2->GetNCells()>2) {
 543         FillHistogram("hMassPtA10",p12.M() ,p12.Pt() );
 544         FillHistogram("hMassPtvA10",pv12.M(),pv12.Pt());
 545         FillHistogram("hMassPtCA10",p12.M() ,p12.Pt(), centr+0.5);
 546         FillHistogram("hMassSingle_all",p12.M(),ph1->Pt()) ;
 547         FillHistogram("hMassSingle_all",p12.M(),ph2->Pt()) ;
 548
 549         if(!eventVtxExist)
 550           FillHistogram("hMassPtA10nvtx",p12.M() ,p12.Pt() );
 551         if(eventVtxExist)
 552           FillHistogram("hMassPtA10vtx"  ,p12.M() ,p12.Pt() );
 553         if(eventVtxExist && eventV0AND)
 554           FillHistogram("hMassPtA10V0AND",p12.M() ,p12.Pt() );
 555         if(eventPileup)
 556           FillHistogram("hMassPtA10PU"   ,p12.M() ,p12.Pt() );
 557
 558         if(ph1->IsCPVOK())
 559           FillHistogram("hMassSingle_cpv",p12.M(),ph1->Pt()) ;
 560         if(ph2->IsCPVOK())
 561           FillHistogram("hMassSingle_cpv",p12.M(),ph2->Pt()) ;
 562         if(ph1->IsDispOK())
 563           FillHistogram("hMassSingle_disp",p12.M(),ph1->Pt()) ;
 564         if(ph2->IsDispOK())
 565           FillHistogram("hMassSingle_disp",p12.M(),ph2->Pt()) ;
 566         if(ph1->IsCPVOK() && ph1->IsDispOK())
 567           FillHistogram("hMassSingle_both",p12.M(),ph1->Pt()) ;
 568         if(ph2->IsCPVOK() && ph2->IsDispOK())
 569           FillHistogram("hMassSingle_both",p12.M(),ph2->Pt()) ;
 570
 571
 572         if(ph1->IsCPVOK() && ph2->IsCPVOK())
 573           FillHistogram("hMassPtCA10_cpv",p12.M() ,p12.Pt(), centr+0.5);
 574         if(ph1->IsDispOK() && ph2->IsDispOK()){
 575           FillHistogram("hMassPtCA10_disp",p12.M() ,p12.Pt(), centr+0.5);
 576           if(ph1->IsCPVOK() && ph2->IsCPVOK())
 577             FillHistogram("hMassPtCA10_both",p12.M() ,p12.Pt(), centr+0.5);
 578         }
 579         if (asym<0.8) {
 580           FillHistogram("hMassPtA08",p12.M(),p12.Pt());
 581         }
 582         if (asym<0.7) {
 583           FillHistogram("hMassPtA07",p12.M(),p12.Pt());
 584           FillHistogram("hMassPtvA07",pv12.M(),pv12.Pt());
 585           FillHistogram("hMassPtCA07",p12.M() ,p12.Pt(), centr+0.5);
 586           if(!eventVtxExist)
 587             FillHistogram("hMassPtA07nvtx",p12.M() ,p12.Pt() );
 588           if(eventVtxExist)
 589             FillHistogram("hMassPtA07vtx"  ,p12.M() ,p12.Pt() );
 590           if(eventVtxExist && eventV0AND)
 591             FillHistogram("hMassPtA07V0AND",p12.M() ,p12.Pt() );
 592           if(eventPileup)
 593             FillHistogram("hMassPtA07PU"   ,p12.M() ,p12.Pt() );
 594           if(ph1->IsCPVOK() && ph2->IsCPVOK())
 595             FillHistogram("hMassPtCA07_cpv",p12.M() ,p12.Pt(), centr+0.5);
 596           if(ph1->IsDispOK() && ph2->IsDispOK()){
 597             FillHistogram("hMassPtCA07_disp",p12.M() ,p12.Pt(), centr+0.5);
 598             if(ph1->IsCPVOK() && ph2->IsCPVOK())
 599               FillHistogram("hMassPtCA07_both",p12.M() ,p12.Pt(), centr+0.5);
 600         }
 601
 602         }
 603         if (asym<0.1) {
 604           FillHistogram("hMassPtA01",p12.M(),p12.Pt());
 605         }
 606         if (TMath::Abs(p12.M()-0.135)<0.03)
 607           FillHistogram("hAsymPtPi0",asym   ,p12.Pt());
 608         if (TMath::Abs(p12.M()-0.547)<0.09)
 609           FillHistogram("hAsymPtEta",asym   ,p12.Pt());
 610
 611         if (ph1->Module()==1 && ph2->Module()==1) FillHistogram("hMassPtM1",p12.M() ,p12.Pt() );
 612         if (ph1->Module()==2 && ph2->Module()==2) FillHistogram("hMassPtM2",p12.M() ,p12.Pt() );
 613         if (ph1->Module()==3 && ph2->Module()==3) FillHistogram("hMassPtM3",p12.M() ,p12.Pt() );
 614         if (ph1->Module()!=3 && ph2->Module()!=3) FillHistogram("hMassPtM12",p12.M() ,p12.Pt() );
 615         if (ph1->Module()!=2 && ph2->Module()!=2) FillHistogram("hMassPtM13",p12.M() ,p12.Pt() );
 616         if (ph1->Module()!=1 && ph2->Module()!=1) FillHistogram("hMassPtM23",p12.M() ,p12.Pt() );
 617
 618       }
 619
 620       if (ph1->GetNCells()>3 && ph2->GetNCells()>3) {
 621         FillHistogram("hMassPtN3",p12.M() ,p12.Pt() );
 622       }
 623       if (ph1->GetNCells()>4 && ph2->GetNCells()>4) {
 624         FillHistogram("hMassPtN4",p12.M() ,p12.Pt() );
 625       }
 626       if (ph1->GetNCells()>5 && ph2->GetNCells()>5) {
 627         FillHistogram("hMassPtN5",p12.M() ,p12.Pt() );
 628       }
 629       if (ph1->GetNCells()>6 && ph2->GetNCells()>6) {
 630         FillHistogram("hMassPtN6",p12.M() ,p12.Pt() );
 631       }
 632
 633     } // end of loop i2
 634   } // end of loop i1
 635
 636   //now mixed
 637   for (Int_t i1=0; i1<inPHOS; i1++) {
 638     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 639     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
 640       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
 641       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
 642       AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
 643       p12  = *ph1  + *ph2;
 644       pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 645       Double_t asym  = TMath::Abs((ph1->Energy()-ph2->Energy())/(ph1->Energy()+ph2->Energy()));
 646
 647       if (ph1->GetNCells()>2 && ph2->GetNCells()>2) {
 648         FillHistogram("hMiMassPtA10",p12.M() ,p12.Pt() );
 649         FillHistogram("hMiMassPtvA10",pv12.M(),pv12.Pt());
 650         FillHistogram("hMiMassPtCA10",p12.M() ,p12.Pt(), centr+0.5);
 651         FillHistogram("hMiMassSingle_all",p12.M(),ph1->Pt()) ;
 652         FillHistogram("hMiMassSingle_all",p12.M(),ph2->Pt()) ;
 653
 654         if(!eventVtxExist)
 655           FillHistogram("hMiMassPtA10nvtx",p12.M() ,p12.Pt() );
 656         if(eventVtxExist)
 657           FillHistogram("hMiMassPtA10vtx"  ,p12.M() ,p12.Pt() );
 658         if(eventVtxExist && eventV0AND)
 659           FillHistogram("hMiMassPtA10V0AND",p12.M() ,p12.Pt() );
 660         if(eventPileup)
 661           FillHistogram("hMiMassPtA10PU"   ,p12.M() ,p12.Pt() );
 662
 663         if(ph1->IsCPVOK())
 664           FillHistogram("hMiMassSingle_cpv",p12.M(),ph1->Pt()) ;
 665         if(ph2->IsCPVOK())
 666           FillHistogram("hMiMassSingle_cpv",p12.M(),ph2->Pt()) ;
 667         if(ph1->IsDispOK())
 668           FillHistogram("hMiMassSingle_disp",p12.M(),ph1->Pt()) ;
 669         if(ph2->IsDispOK())
 670           FillHistogram("hMiMassSingle_disp",p12.M(),ph2->Pt()) ;
 671         if(ph1->IsCPVOK() && ph1->IsDispOK())
 672           FillHistogram("hMiMassSingle_both",p12.M(),ph1->Pt()) ;
 673         if(ph2->IsCPVOK() && ph2->IsDispOK())
 674           FillHistogram("hMiMassSingle_both",p12.M(),ph2->Pt()) ;
 675
 676
 677         if(ph1->IsCPVOK() && ph2->IsCPVOK())
 678           FillHistogram("hMiMassPtCA10_cpv",p12.M() ,p12.Pt(), centr+0.5);
 679         if(ph1->IsDispOK() && ph2->IsDispOK()){
 680           FillHistogram("hMiMassPtCA10_disp",p12.M() ,p12.Pt(), centr+0.5);
 681           if(ph1->IsCPVOK() && ph2->IsCPVOK())
 682             FillHistogram("hMiMassPtCA10_both",p12.M() ,p12.Pt(), centr+0.5);
 683         }
 684         if (asym<0.8) {
 685           FillHistogram("hMiMassPtA08",p12.M(),p12.Pt());
 686         }
 687         if (asym<0.7) {
 688           FillHistogram("hMiMassPtA07",p12.M(),p12.Pt());
 689           FillHistogram("hMiMassPtvA07",pv12.M(),pv12.Pt());
 690           FillHistogram("hMiMassPtCA07",p12.M() ,p12.Pt(), centr+0.5);
 691           if(!eventVtxExist)
 692             FillHistogram("hMiMassPtA07nvtx",p12.M() ,p12.Pt() );
 693           if(eventVtxExist)
 694             FillHistogram("hMiMassPtA07vtx"  ,p12.M() ,p12.Pt() );
 695           if(eventVtxExist && eventV0AND)
 696             FillHistogram("hMiMassPtA07V0AND",p12.M() ,p12.Pt() );
 697           if(eventPileup)
 698             FillHistogram("hMiMassPtA07PU"   ,p12.M() ,p12.Pt() );
 699           if(ph1->IsCPVOK() && ph2->IsCPVOK())
 700             FillHistogram("hMiMassPtCA07_cpv",p12.M() ,p12.Pt(), centr+0.5);
 701           if(ph1->IsDispOK() && ph2->IsDispOK()){
 702             FillHistogram("hMiMassPtCA07_disp",p12.M() ,p12.Pt(), centr+0.5);
 703             if(ph1->IsCPVOK() && ph2->IsCPVOK())
 704               FillHistogram("hMiMassPtCA07_both",p12.M() ,p12.Pt(), centr+0.5);
 705           }
 706         }
 707         if (asym<0.1) {
 708           FillHistogram("hMiMassPtA01",p12.M(),p12.Pt());
 709         }
 710         FillHistogram("hMiAsymPt",asym   ,p12.Pt());
 711
 712         if (ph1->Module()==1 && ph2->Module()==1) FillHistogram("hMiMassPtM1",p12.M() ,p12.Pt() );
 713         if (ph1->Module()==2 && ph2->Module()==2) FillHistogram("hMiMassPtM2",p12.M() ,p12.Pt() );
 714         if (ph1->Module()==3 && ph2->Module()==3) FillHistogram("hMiMassPtM3",p12.M() ,p12.Pt() );
 715         if (ph1->Module()!=3 && ph2->Module()!=3) FillHistogram("hMiMassPtM12",p12.M() ,p12.Pt() );
 716         if (ph1->Module()!=2 && ph2->Module()!=2) FillHistogram("hMiMassPtM13",p12.M() ,p12.Pt() );
 717         if (ph1->Module()!=1 && ph2->Module()!=1) FillHistogram("hMiMassPtM23",p12.M() ,p12.Pt() );
 718
 719       }
 720
 721       if (ph1->GetNCells()>3 && ph2->GetNCells()>3) {
 722         FillHistogram("hMiMassPtN3",p12.M() ,p12.Pt() );
 723       }
 724       if (ph1->GetNCells()>4 && ph2->GetNCells()>4) {
 725         FillHistogram("hMiMassPtN4",p12.M() ,p12.Pt() );
 726       }
 727       if (ph1->GetNCells()>5 && ph2->GetNCells()>5) {
 728         FillHistogram("hMiMassPtN5",p12.M() ,p12.Pt() );
 729       }
 730       if (ph1->GetNCells()>6 && ph2->GetNCells()>6) {
 731         FillHistogram("fMihMassPtN6",p12.M() ,p12.Pt() );
 732       }
 733
 734       } // end of loop i2
 735     }
 736   } // end of loop i1
 737
 738
 739   //Now we either add current events to stack or remove
 740   //If no photons in current event - no need to add it to mixed
 741   if(fPHOSEvent->GetEntriesFast()>0){
 742     prevPHOS->AddFirst(fPHOSEvent) ;
 743     fPHOSEvent=0;
 744     if(prevPHOS->GetSize()>100){//Remove redundant events
 745       TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
 746       prevPHOS->RemoveLast() ;
 747       delete tmp ;
 748     }
 749   }
 750   // Post output data.
 751   PostData(1, fOutputContainer);
 752   fEventCounter++;
 753 }
 754
 755 //________________________________________________________________________
 756 void AliAnalysisTaskPi0::Terminate(Option_t *)
 757 {
 758   // Draw result to the screen
 759   // Called once at the end of the query
 760
 761   Int_t nPP = fnCINT1B - fnCINT1A - fnCINT1C + 2*fnCINT1E;
 762   FillHistogram("hTrigClass",1,fnCINT1B);
 763   FillHistogram("hTrigClass",2,fnCINT1A);
 764   FillHistogram("hTrigClass",3,fnCINT1C);
 765   FillHistogram("hTrigClass",4,fnCINT1E);
 766   FillHistogram("hTrigClass",5,nPP);
 767   Printf("fnCINT1B=%d, fnCINT1A=%d ,fnCINT1C=%d ,fnCINT1E=%d, nPP=%d",
 768      fnCINT1B,fnCINT1A,fnCINT1C,fnCINT1E,nPP);
 769
 770 }
 771
 772 //________________________________________________________________________
 773 Bool_t AliAnalysisTaskPi0::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
 774 {
 775   //Check if this channel belogs to the good ones
 776
 777   if(strcmp(det,"PHOS")==0){
 778     if(mod>5 || mod<1){
 779       AliError(Form("No bad map for PHOS module %d ",mod)) ;
 780       return kTRUE ;
 781     }
 782     if(!fPHOSBadMap[mod]){
 783       AliError(Form("No Bad map for PHOS module %d",mod)) ;
 784       return kTRUE ;
 785     }
 786     if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
 787       return kFALSE ;
 788     else
 789       return kTRUE ;
 790   }
 791   else{
 792     AliError(Form("Can not find bad channels for detector %s ",det)) ;
 793   }
 794   return kTRUE ;
 795 }
 796 //_____________________________________________________________________________
 797 void AliAnalysisTaskPi0::FillHistogram(const char * key,Double_t x)const{
 798   //FillHistogram
 799   TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
 800   if(tmpI){
 801     tmpI->Fill(x) ;
 802     return ;
 803   }
 804   TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
 805   if(tmpF){
 806     tmpF->Fill(x) ;
 807     return ;
 808   }
 809   TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
 810   if(tmpD){
 811     tmpD->Fill(x) ;
 812     return ;
 813   }
 814   AliInfo(Form("can not find histogram <%s> ",key)) ;
 815 }
 816 //_____________________________________________________________________________
 817 void AliAnalysisTaskPi0::FillHistogram(const char * key,Double_t x,Double_t y)const{
 818   //FillHistogram
 819   TObject * tmp = fOutputContainer->FindObject(key) ;
 820   if(!tmp){
 821     AliInfo(Form("can not find histogram <%s> ",key)) ;
 822     return ;
 823   }
 824   if(tmp->IsA() == TClass::GetClass("TH1F")){
 825     ((TH1F*)tmp)->Fill(x,y) ;
 826     return ;
 827   }
 828   if(tmp->IsA() == TClass::GetClass("TH2F")){
 829     ((TH2F*)tmp)->Fill(x,y) ;
 830     return ;
 831   }
 832   AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
 833 }
 834
 835 //_____________________________________________________________________________
 836 void AliAnalysisTaskPi0::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
 837   //Fills 1D histograms with key
 838   TObject * tmp = fOutputContainer->FindObject(key) ;
 839   if(!tmp){
 840     AliInfo(Form("can not find histogram <%s> ",key)) ;
 841     return ;
 842   }
 843   if(tmp->IsA() == TClass::GetClass("TH2F")){
 844     ((TH2F*)tmp)->Fill(x,y,z) ;
 845     return ;
 846   }
 847   if(tmp->IsA() == TClass::GetClass("TH3F")){
 848     ((TH3F*)tmp)->Fill(x,y,z) ;
 849     return ;
 850   }
 851 }
 852 //_____________________________________________________________________________
 853 Bool_t AliAnalysisTaskPi0::TestLambda(Double_t l1,Double_t l2){
 854   Double_t l1Mean=1.22 ;
 855   Double_t l2Mean=2.0 ;
 856   Double_t l1Sigma=0.42 ;
 857   Double_t l2Sigma=0.71 ;
 858   Double_t c=-0.59 ;
 859   Double_t R2=(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma+(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma-c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
 860   return (R2<9.) ;
 861
 862 }