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