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