Coverity fixes (Jens)
[u/mrichter/AliRoot.git] / PWGGA / PHOSTasks / PHOS_PbPbQA / AliAnalysisTaskPHOSPbPbQA.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 #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
37 ClassImp(AliAnalysisTaskPHOSPbPbQA)
38
39 //________________________________________________________________________
40 AliAnalysisTaskPHOSPbPbQA::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   
51   // Initialize the PHOS geometry 
52   fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
53   
54 }
55
56 //________________________________________________________________________
57 AliAnalysisTaskPHOSPbPbQA::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 //________________________________________________________________________
80 void 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 //________________________________________________________________________
132 void 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
161   if( fCentrality < 0. ){
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) {
267           snprintf(key,55,"hMiPi0AllSM%d_cen%d",sm1,fCenBin) ;
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 //_____________________________________________________________________________
293 void 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 //_____________________________________________________________________________
313 void 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 }