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