Fixes for coverity checks.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardQATask.cxx
1 // 
2 // Calculate the multiplicity in the forward regions event-by-event 
3 // 
4 // Inputs: 
5 //   - AliESDEvent 
6 //
7 // Outputs: 
8 //   - AliAODForwardMult 
9 // 
10 // Histograms 
11 //   
12 // Corrections used 
13 //
14 #include "AliForwardQATask.h"
15 #include "AliForwardUtil.h"
16 #include "AliTriggerAnalysis.h"
17 #include "AliPhysicsSelection.h"
18 #include "AliLog.h"
19 #include "AliESDEvent.h"
20 #include "AliAODHandler.h"
21 #include "AliMultiplicity.h"
22 #include "AliInputEventHandler.h"
23 #include "AliForwardCorrectionManager.h"
24 #include "AliAnalysisManager.h"
25 #include "AliAODForwardMult.h"
26 #include <TH1.h>
27 #include <TDirectory.h>
28 #include <TTree.h>
29 #include <TROOT.h>
30 #include <TStopwatch.h>
31
32 //====================================================================
33 AliForwardQATask::AliForwardQATask()
34   : AliAnalysisTaskSE(),
35     fEnableLowFlux(false), 
36     fFirstEvent(true),
37     fCorrManager(0),
38     fESDFMD(),
39     fHistos(),
40     fEventInspector(),
41     fEnergyFitter(),
42     fSharingFilter(),
43     fDensityCalculator(),
44     fList(0),
45     fDebug(0)
46 {
47   // 
48   // Constructor
49   //
50 }
51
52 //____________________________________________________________________
53 AliForwardQATask::AliForwardQATask(const char* name)
54   : AliAnalysisTaskSE(name),
55     fEnableLowFlux(false), 
56     fFirstEvent(true),
57     fCorrManager(0),
58     fESDFMD(),
59     fHistos(),
60     fEventInspector("event"),
61     fEnergyFitter("energy"),
62     fSharingFilter("sharing"), 
63     fDensityCalculator("density"),
64     fList(0),
65     fDebug(0)
66 {
67   // 
68   // Constructor 
69   // 
70   // Parameters:
71   //    name Name of task 
72   //
73   DefineOutput(1, TList::Class());
74   DefineOutput(2, TList::Class());
75   fCorrManager = &AliForwardCorrectionManager::Instance(); 
76   fEnergyFitter.SetNParticles(1); // Just find the 1st peak 
77   fEnergyFitter.SetDoMakeObject(false); 
78   fEnergyFitter.SetUseIncreasingBins(true);
79   fEnergyFitter.SetDoFits(kTRUE);
80   fEnergyFitter.SetLowCut(0.4);
81   fEnergyFitter.SetFitRangeBinWidth(4);
82   fEnergyFitter.SetMinEntries(1000);
83 }
84
85 //____________________________________________________________________
86 AliForwardQATask::AliForwardQATask(const AliForwardQATask& o)
87   : AliAnalysisTaskSE(o),
88     fEnableLowFlux(o.fEnableLowFlux), 
89     fFirstEvent(o.fFirstEvent),
90     fCorrManager(o.fCorrManager),
91     fESDFMD(o.fESDFMD),
92     fHistos(o.fHistos),
93     fEventInspector(o.fEventInspector),
94     fEnergyFitter(o.fEnergyFitter),
95     fSharingFilter(o.fSharingFilter),
96     fDensityCalculator(o.fDensityCalculator),
97     fList(o.fList),
98     fDebug(o.fDebug) 
99 {
100   // 
101   // Copy constructor 
102   // 
103   // Parameters:
104   //    o Object to copy from 
105   //
106   DefineOutput(1, TList::Class());
107   DefineOutput(2, TList::Class());
108 }
109
110 //____________________________________________________________________
111 AliForwardQATask&
112 AliForwardQATask::operator=(const AliForwardQATask& o)
113 {
114   // 
115   // Assignment operator 
116   // 
117   // Parameters:
118   //    o Object to assign from 
119   // 
120   // Return:
121   //    Reference to this object 
122   //
123   if (&o == this) return *this;
124   AliAnalysisTaskSE::operator=(o);
125
126   fEnableLowFlux     = o.fEnableLowFlux;
127   fFirstEvent        = o.fFirstEvent;
128   fCorrManager       = o.fCorrManager;
129   fEventInspector    = o.fEventInspector;
130   fEnergyFitter      = o.fEnergyFitter;
131   fSharingFilter     = o.fSharingFilter;
132   fDensityCalculator = o.fDensityCalculator;
133   fHistos            = o.fHistos;
134   fList              = o.fList;
135   fDebug             = o.fDebug;
136
137   return *this;
138 }
139
140 //____________________________________________________________________
141 void
142 AliForwardQATask::SetDebug(Int_t dbg)
143 {
144   // 
145   // Set debug level 
146   // 
147   // Parameters:
148   //    dbg Debug level
149   //
150   fDebug = dbg;
151   fEventInspector.SetDebug(dbg);
152   fEnergyFitter.SetDebug(dbg);
153   fSharingFilter.SetDebug(dbg);
154   fDensityCalculator.SetDebug(dbg);
155 }
156
157 //____________________________________________________________________
158 Bool_t 
159 AliForwardQATask::CheckCorrections(UInt_t what) const
160 {
161   // 
162   // Check if all needed corrections are there and accounted for.  If not,
163   // do a Fatal exit 
164   // 
165   // Parameters:
166   //    what Which corrections is needed
167   // 
168   // Return:
169   //    true if all present, false otherwise
170   //  
171
172   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
173   // Check that we have the energy loss fits, needed by 
174   //   AliFMDSharingFilter 
175   //   AliFMDDensityCalculator 
176   if (what & AliForwardCorrectionManager::kELossFits && !fcm.GetELossFit()) { 
177     AliWarning("No energy loss fits");
178     return false;
179   }
180   return true;
181 }
182
183 //____________________________________________________________________
184 Bool_t
185 AliForwardQATask::ReadCorrections(const TAxis*& pe, 
186                                   const TAxis*& pv, 
187                                   Bool_t        mc)
188 {
189   //
190   // Read corrections
191   //
192   //
193   UInt_t what = AliForwardCorrectionManager::kAll;
194   what ^= AliForwardCorrectionManager::kDoubleHit;
195   what ^= AliForwardCorrectionManager::kVertexBias;
196   what ^= AliForwardCorrectionManager::kAcceptance;
197   what ^= AliForwardCorrectionManager::kMergingEfficiency;
198
199   AliForwardCorrectionManager& fcm = AliForwardCorrectionManager::Instance();
200   if (!fcm.Init(GetEventInspector().GetCollisionSystem(),
201                 GetEventInspector().GetEnergy(),
202                 GetEventInspector().GetField(),
203                 mc,
204                 what)) return false;
205   if (!CheckCorrections(what)) {
206     return false;
207   }
208
209   // Sett our persistency pointer 
210   // fCorrManager = &fcm;
211
212   // Get the eta axis from the secondary maps - if read in
213   if (!pe) {
214     pe = fcm.GetEtaAxis();
215     if (!pe) AliFatal("No eta axis defined");
216   }
217   // Get the vertex axis from the secondary maps - if read in
218   if (!pv) {
219     pv = fcm.GetVertexAxis();
220     if (!pv) AliFatal("No vertex axis defined");
221   }
222
223   return true;
224 }
225
226 //____________________________________________________________________
227 AliESDEvent*
228 AliForwardQATask::GetESDEvent()
229 {
230   //
231   // Get the ESD event. IF this is the first event, initialise
232   //
233   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
234   if (!esd) {
235     AliWarning("No ESD event found for input event");
236     return 0;
237   }
238
239   // On the first event, initialize the parameters
240   if (fFirstEvent && esd->GetESDRun()) {
241     GetEventInspector().ReadRunDetails(esd);
242
243     AliInfoF("Initializing with parameters from the ESD:\n"
244              "         AliESDEvent::GetBeamEnergy()   ->%f\n"
245              "         AliESDEvent::GetBeamType()     ->%s\n"
246              "         AliESDEvent::GetCurrentL3()    ->%f\n"
247              "         AliESDEvent::GetMagneticField()->%f\n"
248              "         AliESDEvent::GetRunNumber()    ->%d\n",
249              esd->GetBeamEnergy(),
250              esd->GetBeamType(),
251              esd->GetCurrentL3(),
252              esd->GetMagneticField(),
253              esd->GetRunNumber());
254
255
256     if (!InitializeSubs()) {
257       AliWarning("Initialisation of sub algorithms failed!");
258       return 0;
259     }
260     AliInfoF("Clearing first event flag from %s to false", 
261              fFirstEvent ? "true" : "false");
262     fFirstEvent = false;
263   }
264   return esd;
265 }
266 //____________________________________________________________________
267 Bool_t
268 AliForwardQATask::InitializeSubs()
269 {
270   // 
271   // Initialise the sub objects and stuff.  Called on first event 
272   // 
273   //
274   const TAxis* pe = 0;
275   const TAxis* pv = 0;
276
277   if (!ReadCorrections(pe,pv)) { 
278     AliWarning("Failed to read corrections");
279     return false;
280   }
281
282   fHistos.Init(*pe);
283
284   fEventInspector.Init(*pv);
285   fEnergyFitter.Init(*pe);
286   fSharingFilter.Init();
287   fDensityCalculator.Init(*pe);
288
289   this->Print();
290
291   return true;
292 }
293
294 //____________________________________________________________________
295 void
296 AliForwardQATask::UserCreateOutputObjects()
297 {
298   // 
299   // Create output objects 
300   // 
301   //
302   fList = new TList;
303   fList->SetOwner();
304   
305   fEventInspector.DefineOutput(fList);
306   fEnergyFitter.DefineOutput(fList);
307   fSharingFilter.DefineOutput(fList);
308   fDensityCalculator.DefineOutput(fList);
309
310   PostData(1, fList);
311 }
312 //____________________________________________________________________
313 void
314 AliForwardQATask::UserExec(Option_t*)
315 {
316   // 
317   // Process each event 
318   // 
319   // Parameters:
320   //    option Not used
321   //  
322
323   // static Int_t cnt = 0;
324   // cnt++;
325   // Get the input data 
326   AliESDEvent* esd = GetESDEvent();
327   if (!esd) { 
328     AliWarning("Got no ESD event");
329     return;
330   }
331   if (fFirstEvent) { 
332     // If the first event flag wasn't cleared in the above call to
333     // GetESDEvent, we should not do anything, since nothing has been
334     // initialised yet, so we opt out here (with a warning) 
335     AliWarning("Nothing has been initialized yet, opt'ing out");
336     return;
337   }
338
339   // Clear stuff 
340   fHistos.Clear();
341   fESDFMD.Clear();
342   
343   Bool_t   lowFlux   = kFALSE;
344   UInt_t   triggers  = 0;
345   UShort_t ivz       = 0;
346   Double_t vz        = 0;
347   Double_t cent      = -1;
348   UShort_t nClusters = 0;
349   UInt_t   found     = fEventInspector.Process(esd, triggers, lowFlux, 
350                                                ivz, vz, cent, nClusters);
351   
352   if (found & AliFMDEventInspector::kNoEvent)    return;
353   if (found & AliFMDEventInspector::kNoTriggers) return;
354   if (found & AliFMDEventInspector::kNoSPD)      return;
355   if (found & AliFMDEventInspector::kNoFMD)      return;
356   if (found & AliFMDEventInspector::kNoVertex)   return;
357   if (triggers & AliAODForwardMult::kPileUp)     return;
358   if (found & AliFMDEventInspector::kBadVertex)  return;
359
360   // We we do not want to use low flux specific code, we disable it here. 
361   if (!fEnableLowFlux) lowFlux = false;
362
363   // Get FMD data 
364   AliESDFMD* esdFMD = esd->GetFMDData();
365   
366   // Run the energy loss fitter 
367   if (!fEnergyFitter.Accumulate(*esdFMD, cent, 
368                                 triggers & AliAODForwardMult::kEmpty)) {
369     AliWarning("Energy fitter failed");
370     return;
371   }
372     
373   //  // Apply the sharing filter (or hit merging or clustering if you like)
374   if (!fSharingFilter.Filter(*esdFMD, lowFlux, fESDFMD)) { 
375     AliWarning("Sharing filter failed!");
376     return;
377   }
378
379   // Calculate the inclusive charged particle density 
380   if (!fDensityCalculator.Calculate(fESDFMD, fHistos, ivz, lowFlux)) { 
381     // if (!fDensityCalculator.Calculate(*esdFMD, fHistos, ivz, lowFlux)) { 
382     AliWarning("Density calculator failed!");
383     return;
384   }
385   PostData(1, fList);
386 }
387
388 //____________________________________________________________________
389 void
390 AliForwardQATask::Terminate(Option_t*)
391 {
392   // 
393   // End of job
394   // 
395   // Parameters:
396   //    option Not used 
397   //
398   if (fDebug) AliInfo("In Forwards terminate");
399   TStopwatch swt;
400   swt.Start();
401
402   TList* list = dynamic_cast<TList*>(GetOutputData(1));
403   if (!list) {
404     AliError(Form("No output list defined (%p)", GetOutputData(1)));
405     if (GetOutputData(1)) GetOutputData(1)->Print();
406     return;
407   }
408   
409   // Get our histograms from the container 
410   TH1I* hEventsTr    = 0;//static_cast<TH1I*>(list->FindObject("nEventsTr"));
411   TH1I* hEventsTrVtx = 0;//static_cast<TH1I*>(list->FindObject("nEventsTrVtx"));
412   TH1I* hTriggers    = 0;
413   if (!fEventInspector.FetchHistograms(list, hEventsTr, 
414                                        hEventsTrVtx, hTriggers)) { 
415     AliError(Form("Didn't get histograms from event selector "
416                   "(hEventsTr=%p,hEventsTrVtx=%p)", 
417                   hEventsTr, hEventsTrVtx));
418     return;
419   }
420
421   TStopwatch swf;
422   swf.Start();
423   fEnergyFitter.Fit(list);
424   swf.Stop();
425   AliInfoF("Fitting took %d real-time seconds, and %f CPU seconds", 
426            Int_t(swf.RealTime()), swf.CpuTime());
427
428   fSharingFilter.ScaleHistograms(list,Int_t(hEventsTr->Integral()));
429   fDensityCalculator.ScaleHistograms(list,Int_t(hEventsTrVtx->Integral()));
430
431   // Make a deep copy and post that as output 2 
432   TList* list2 = static_cast<TList*>(list->Clone(Form("%sResults", 
433                                                       list->GetName())));
434   if (fDebug) AliInfoF("Posting post processing results to %s", 
435                        list2->GetName());
436   list2->SetOwner();
437   PostData(2, list2);
438
439   swt.Stop();
440   AliInfoF("Terminate took %d real-time seconds, and %f CPU seconds", 
441            Int_t(swt.RealTime()), swt.CpuTime());
442
443 }
444
445 //____________________________________________________________________
446 void
447 AliForwardQATask::Print(Option_t* option) const
448 {
449   // 
450   // Print information 
451   // 
452   // Parameters:
453   //    option Not used
454   //
455   
456   std::cout << ClassName() << ": " << GetName() << "\n" 
457             << "  Enable low flux code:   " << (fEnableLowFlux ? "yes" : "no") 
458             << "\n"
459             << "  Off-line trigger mask:  0x" 
460             << std::hex     << std::setfill('0') 
461             << std::setw (8) << fOfflineTriggerMask 
462             << std::dec     << std::setfill (' ') << std::endl;
463   gROOT->IncreaseDirLevel();
464   if (fCorrManager) fCorrManager->Print();
465   else  
466     std::cout << "  Correction manager not set yet" << std::endl;
467   GetEventInspector()   .Print(option);
468   GetEnergyFitter()     .Print(option);
469   GetSharingFilter()    .Print(option);
470   gROOT->DecreaseDirLevel();
471 }
472
473 //
474 // EOF
475 //