Coverity fixes (Jens)
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_PbPbQA / AliAnalysisTaskPHOSPbPbQA.cxx
CommitLineData
6c0d310b 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#include "AliAnalysisManager.h"
18#include "AliAnalysisTaskSE.h"
19#include "AliAnalysisTaskPHOSPbPbQA.h"
20#include "AliCaloPhoton.h"
21#include "AliPHOSGeometry.h"
22#include "AliESDEvent.h"
23#include "AliESDCaloCells.h"
24#include "AliCentrality.h"
25#include "AliLog.h"
26#include "TObjArray.h"
27#include "TList.h"
28#include "TH1.h"
29#include "TH2.h"
30
31// Stripped-down version of Dmitri Peressounko' AliAnalysisTaskPi0Flow class
32// used for the fast QA of PbPb data.
33//...
34// Author: Boris Polishchuk
35// Date : 19.10.2011
36
37ClassImp(AliAnalysisTaskPHOSPbPbQA)
38
9aa24ff1 39//________________________________________________________________________
40AliAnalysisTaskPHOSPbPbQA::AliAnalysisTaskPHOSPbPbQA() : AliAnalysisTaskSE(),
41 fOutputContainer(0),fPHOSEvent(0),fCentrality(0),fCenBin(0),
42 fPHOSGeo(0),fEventCounter(0)
43{
44 //Default constructor
45
46 for(Int_t i=0;i<1;i++){
47 for(Int_t j=0;j<2;j++)
48 fPHOSEvents[i][j]=0 ;
49 }
50
17e719cf 51 // Initialize the PHOS geometry
52 fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
53
9aa24ff1 54}
55
6c0d310b 56//________________________________________________________________________
57AliAnalysisTaskPHOSPbPbQA::AliAnalysisTaskPHOSPbPbQA(const char *name)
58: AliAnalysisTaskSE(name),
59 fOutputContainer(0),
60 fPHOSEvent(0),
61 fCentrality(0),fCenBin(0),
62 fPHOSGeo(0),
63 fEventCounter(0)
64{
65 // Constructor
66 for(Int_t i=0;i<1;i++){
67 for(Int_t j=0;j<2;j++)
68 fPHOSEvents[i][j]=0 ;
69 }
70
71 // Output slots #0 write into a TH1 container
72 DefineOutput(1,TList::Class());
73
74 // Initialize the PHOS geometry
75 fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
76
77}
78
79//________________________________________________________________________
80void AliAnalysisTaskPHOSPbPbQA::UserCreateOutputObjects()
81{
82 // Create histograms
83 // Called once
84
85 // ESD histograms
86 if(fOutputContainer != NULL){
87 delete fOutputContainer;
88 }
89
90 fOutputContainer = new TList();
91 fOutputContainer->SetOwner(kTRUE);
92
93
94 fOutputContainer->Add(new TH2F("hCenPHOS","Centrality vs PHOSclusters", 100,0.,100.,200,0.,200.)) ;
95 fOutputContainer->Add(new TH2F("hCenPHOSCells","Centrality vs PHOS cells", 100,0.,100.,100,0.,1000.)) ;
96 fOutputContainer->Add(new TH2F("hCenTrack","Centrality vs tracks", 100,0.,100.,100,0.,15000.)) ;
97
98 //pi0 spectrum
99 Int_t nPtPhot = 300 ;
100 Double_t ptPhotMax = 30 ;
101 Int_t nM = 500;
102 Double_t mMin = 0.0;
103 Double_t mMax = 1.0;
104 char key[55] ;
105
106 for(Int_t cent=0; cent<2; cent++){
107
108 snprintf(key,55,"hPi0All_cen%d",cent) ;
109 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
110
111 snprintf(key,55,"hMiPi0All_cen%d",cent) ;
112 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
113 }
114
115 //per module
116 for(Int_t cent=0; cent<2; cent++){
117 for(Int_t sm=1; sm<4; sm++) {
118
119 snprintf(key,55,"hPi0AllSM%d_cen%d",sm,cent) ;
120 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
121
122 snprintf(key,55,"hMiPi0AllSM%d_cen%d",sm,cent) ;
123 fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
124 }
125 }
126
127 PostData(1, fOutputContainer);
128
129}
130
131//________________________________________________________________________
132void AliAnalysisTaskPHOSPbPbQA::UserExec(Option_t *)
133{
134 // Main loop, called for each event
135 // Analyze ESD/AOD
136
137 AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
138
139 if (!event) {
140 Printf("ERROR: Could not retrieve event");
141 PostData(1, fOutputContainer);
142 return;
143 }
144
145 // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
146 char key[55] ;
147
148 if(fEventCounter == 0) {
149 for(Int_t mod=0; mod<5; mod++) {
150 if(!event->GetPHOSMatrix(mod)) continue;
151 fPHOSGeo->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
152 }
153 fEventCounter++ ;
154 }
155
156 Int_t zvtx=0 ;
157
158 AliCentrality *centrality = event->GetCentrality();
159 fCentrality=centrality->GetCentralityPercentile("V0M");
160
4ac2548a 161 if( fCentrality < 0. ){
6c0d310b 162 PostData(1, fOutputContainer);
163 return;
164 }
165
166 if(fCentrality < 20.) fCenBin = 0;
167 else fCenBin = 1;
168
169 printf("centrality %.3f [%d]\n",fCentrality,fCenBin);
170
171 if(!fPHOSEvents[zvtx][fCenBin])
172 fPHOSEvents[zvtx][fCenBin]=new TList() ;
173
174 TList * prevPHOS = fPHOSEvents[zvtx][fCenBin] ;
175
176 if(fPHOSEvent)
177 fPHOSEvent->Clear() ;
178 else
179 fPHOSEvent = new TClonesArray("AliCaloPhoton",200) ;
180
181 Int_t multClust = event->GetNumberOfCaloClusters();
182 AliESDCaloCells * cells = event->GetPHOSCells() ;
183
184 FillHistogram("hCenPHOSCells",fCentrality,cells->GetNumberOfCells()) ;
185 FillHistogram("hCenTrack",fCentrality,event->GetNumberOfTracks()) ;
186
187 Int_t inPHOS = 0;
188 Double_t vtx0[3] = {0,0,0}; // vertex
189
190 for (Int_t i=0; i<multClust; i++) {
191
192 AliESDCaloCluster *clu = event->GetCaloCluster(i);
193
194 if ( !clu->IsPHOS() || clu->E()<0.3) continue;
195 if(clu->GetNCells()<3) continue;
196
197 Float_t position[3];
198 clu->GetPosition(position);
199 TVector3 global(position) ;
200 Int_t relId[4] ;
201 fPHOSGeo->GlobalPos2RelId(global,relId) ;
202 Int_t mod = relId[0] ;
203 Int_t cellX = relId[2];
204 Int_t cellZ = relId[3] ;
205
206 TLorentzVector pv1 ;
207 clu->GetMomentum(pv1 ,vtx0);
208
209 if(inPHOS>=fPHOSEvent->GetSize()){
210 fPHOSEvent->Expand(inPHOS+50) ;
211 }
212
213 new((*fPHOSEvent)[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
214 AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
215 ph->SetModule(mod) ;
216 ph->SetMomV2(&pv1) ;
217 ph->SetNCells(clu->GetNCells());
218
219 ph->SetEMCx(float(cellX)) ;
220 ph->SetEMCz(float(cellZ)) ;
221
222 inPHOS++ ;
223 }
224
225 FillHistogram("hCenPHOS",fCentrality,inPHOS) ;
226
227 //pi0
228 for (Int_t i1=0; i1<inPHOS-1; i1++) {
229 AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
230 Int_t sm1 = ph1->Module();
231
232 for (Int_t i2=i1+1; i2<inPHOS; i2++) {
233 AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent->At(i2) ;
234
235 Int_t sm2 = ph2->Module();
236 TLorentzVector p12 = *ph1 + *ph2;
237
238 snprintf(key,55,"hPi0All_cen%d",fCenBin) ;
239 FillHistogram(key,p12.M() ,p12.Pt()) ;
240
241 if(sm1==sm2) {
242 snprintf(key,55,"hPi0AllSM%d_cen%d",sm1,fCenBin) ;
243 FillHistogram(key,p12.M() ,p12.Pt()) ;
244 }
245
246 } // end of loop i2
247 } // end of loop i1
248
249 //now mixed
250 for (Int_t i1=0; i1<inPHOS; i1++) {
251 AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
252 Int_t sm1 = ph1->Module();
253
254 for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
255 TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
256
257 for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
258 AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
259
260 Int_t sm2 = ph2->Module();
261 TLorentzVector p12 = *ph1 + *ph2;
262
263 snprintf(key,55,"hMiPi0All_cen%d",fCenBin) ;
264 FillHistogram(key,p12.M() ,p12.Pt()) ;
265
266 if(sm1==sm2) {
1e404dfe 267 snprintf(key,55,"hMiPi0AllSM%d_cen%d",sm1,fCenBin) ;
6c0d310b 268 FillHistogram(key,p12.M() ,p12.Pt()) ;
269 }
270
271 } // end of loop i2
272 }
273 } // end of loop i1
274
275
276 //Now we either add current events to stack or remove
277 //If no photons in current event - no need to add it to mixed
278 if(fPHOSEvent->GetEntriesFast()>0){
279 prevPHOS->AddFirst(fPHOSEvent) ;
280 fPHOSEvent=0;
281 if(prevPHOS->GetSize()>100){//Remove redundant events
282 TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
283 prevPHOS->RemoveLast() ;
284 delete tmp ;
285 }
286 }
287 // Post output data.
288 PostData(1, fOutputContainer);
289 fEventCounter++;
290}
291
292//_____________________________________________________________________________
293void AliAnalysisTaskPHOSPbPbQA::FillHistogram(const char * key,Double_t x)const{
294 //FillHistogram
295 TH1I * tmpI = dynamic_cast<TH1I*>(fOutputContainer->FindObject(key)) ;
296 if(tmpI){
297 tmpI->Fill(x) ;
298 return ;
299 }
300 TH1F * tmpF = dynamic_cast<TH1F*>(fOutputContainer->FindObject(key)) ;
301 if(tmpF){
302 tmpF->Fill(x) ;
303 return ;
304 }
305 TH1D * tmpD = dynamic_cast<TH1D*>(fOutputContainer->FindObject(key)) ;
306 if(tmpD){
307 tmpD->Fill(x) ;
308 return ;
309 }
310 AliInfo(Form("can not find histogram <%s> ",key)) ;
311}
312//_____________________________________________________________________________
313void AliAnalysisTaskPHOSPbPbQA::FillHistogram(const char * key,Double_t x,Double_t y)const{
314 //FillHistogram
315 TObject * tmp = fOutputContainer->FindObject(key) ;
316 if(!tmp){
317 AliInfo(Form("can not find histogram <%s> ",key)) ;
318 return ;
319 }
320 if(tmp->IsA() == TClass::GetClass("TH1F")){
321 ((TH1F*)tmp)->Fill(x,y) ;
322 return ;
323 }
324 if(tmp->IsA() == TClass::GetClass("TH2F")){
325 ((TH2F*)tmp)->Fill(x,y) ;
326 return ;
327 }
328 AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
329}