Some fixes
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliCentralMultiplicityTask.cxx
1 //====================================================================
2 // 
3 // Base class for classes that calculate the multiplicity in the
4 // central region event-by-event
5 // 
6 // Inputs: 
7 //   - AliESDEvent 
8 //
9 // Outputs: 
10 //   - AliAODCentralMult 
11 // 
12 // Histograms 
13 //   
14 // Corrections used 
15 #include "AliCentralMultiplicityTask.h"
16 #include "AliForwardCorrectionManager.h"
17 #include "AliForwardUtil.h"
18 #include "AliLog.h"
19 #include "AliAODHandler.h"
20 #include "AliInputEventHandler.h"
21 #include "AliESDInputHandler.h"
22 #include "AliAnalysisManager.h"
23 #include "AliESDEvent.h"
24 #include "AliMultiplicity.h"
25 #include "AliFMDEventInspector.h"
26 #include <TROOT.h>
27 #include <TFile.h>
28 #include <TError.h>
29 #include <iostream>
30 #include <iomanip>
31
32 //====================================================================
33 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const char* name) 
34   : AliAnalysisTaskSE(name),
35     fData(0),
36     fList(0),
37     fAODCentral(kFALSE),
38     fManager(),
39     fUseSecondary(true),
40     fUseAcceptance(true),
41     fFirstEventSeen(false)
42 {
43   // 
44   // Constructor 
45   //   
46   DefineOutput(1, TList::Class());
47   fBranchNames = 
48     "ESD:AliESDRun.,AliESDHeader.,AliMultiplicity.,"
49     "SPDVertex.,PrimaryVertex.";
50 }
51 //____________________________________________________________________
52 AliCentralMultiplicityTask::AliCentralMultiplicityTask() 
53   : AliAnalysisTaskSE(),
54     fData(0),
55     fList(0),
56     fAODCentral(),
57     fManager(),
58     fUseSecondary(true),
59     fUseAcceptance(true),
60     fFirstEventSeen(false)
61 {
62   // 
63   // Constructor 
64   // 
65 }
66 //____________________________________________________________________
67 AliCentralMultiplicityTask::AliCentralMultiplicityTask(const AliCentralMultiplicityTask& o)
68   : AliAnalysisTaskSE(o),
69     fData(o.fData),
70     fList(o.fList),
71     fAODCentral(o.fAODCentral),
72     fManager(o.fManager),
73     fUseSecondary(o.fUseSecondary),
74     fUseAcceptance(o.fUseAcceptance),
75     fFirstEventSeen(o.fFirstEventSeen)
76 {
77   //
78   // Copy constructor 
79   // 
80 }
81 //____________________________________________________________________
82 AliCentralMultiplicityTask&
83 AliCentralMultiplicityTask::operator=(const AliCentralMultiplicityTask& o)
84 {
85   // 
86   // Assignment operator 
87   //
88   fData           = o.fData;
89   fList           = o.fList;
90   fAODCentral     = o.fAODCentral;
91   fManager        = o.fManager;
92   fUseSecondary   = o.fUseSecondary;
93   fUseAcceptance  = o.fUseAcceptance;
94   fFirstEventSeen = o.fFirstEventSeen;
95   return *this;
96 }
97 //____________________________________________________________________
98 void AliCentralMultiplicityTask::UserCreateOutputObjects() 
99 {
100   // 
101   // Create output objects 
102   // 
103   //
104
105   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
106   AliAODHandler*      ah = 
107     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
108   if (!ah) AliFatal("No AOD output handler set in analysis manager");
109   
110   
111   TObject* obj = &fAODCentral;
112   ah->AddBranch("AliAODCentralMult", &obj);
113   
114   fList = new TList();
115   fList->SetOwner();
116   PostData(1,fList);
117   
118 }
119 //____________________________________________________________________
120 void AliCentralMultiplicityTask::UserExec(Option_t* /*option*/) 
121 {
122   // 
123   // Process each event 
124   // 
125   // Parameters:
126   //    option Not used
127   //  
128   
129   AliESDInputHandler* eventHandler = 
130     dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()
131                                        ->GetInputEventHandler());
132   if (!eventHandler) {
133     AliWarning("No inputhandler found for this event!");
134     return;}
135   
136   AliESDEvent* esd = eventHandler->GetEvent();
137   
138   if(!GetManager().IsInit() && !fFirstEventSeen) {
139     AliFMDEventInspector inspector;
140     inspector.ReadRunDetails(esd);
141     GetManager().Init(inspector.GetCollisionSystem(),
142                       inspector.GetEnergy(),
143                       inspector.GetField());
144     
145     AliInfo("Manager of corrections in AliCentralMultiplicityTask init");
146     fFirstEventSeen = kTRUE;
147   }
148   
149   //Selecting only events with |valid vertex| < 10 cm
150   const AliESDVertex* vertex = esd->GetPrimaryVertexSPD();
151   if (!vertex) return;
152   if(!(vertex->GetStatus())) return;
153   if(vertex->GetNContributors() <= 0) return ;
154   if(vertex->GetZRes() > 0.1 ) return;
155   Double_t vertexXYZ[3]={0,0,0};
156   vertex->GetXYZ(vertexXYZ);
157   if(TMath::Abs(vertexXYZ[2]) > 10) return;
158   
159   Double_t delta           = 2 ;
160   Double_t vertexBinDouble = (vertexXYZ[2] + 10) / delta;
161   //HHD: The vtxbins are 1-10, not 0-9
162   Int_t    vtxbin          = Int_t(vertexBinDouble + 1) ; 
163   
164   // Make sure AOD is filled
165   AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
166   AliAODHandler*      ah = 
167     dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
168   if (!ah)  
169     AliFatal("No AOD output handler set in analysis manager");
170   
171   ah->SetFillAOD(kTRUE);
172   
173   //Doing analysis
174   fAODCentral.Clear("");
175   TH2D *aodHist = &(fAODCentral.GetHistogram());
176   
177   const AliMultiplicity* spdmult = esd->GetMultiplicity();
178   //Filling clusters in layer 1 used for tracklets...
179   for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
180     aodHist->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
181
182   //...and then the unused ones in layer 1 
183   for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) 
184     aodHist->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),
185                   spdmult->GetPhiSingle(j));
186   
187   
188   // Corrections
189   TH1D* hAcceptance = fManager.GetAcceptanceCorrection(vtxbin);
190   TH2D* hSecMap     = fManager.GetSecMapCorrection(vtxbin);
191   
192   if (!hSecMap)     AliFatal("No secondary map!");
193   if (!hAcceptance) AliFatal("No acceptance!");
194     
195   if (fUseSecondary && hSecMap) aodHist->Divide(hSecMap);
196   
197   for(Int_t nx = 1; nx <= aodHist->GetNbinsX(); nx++) {
198     Float_t accCor = hAcceptance->GetBinContent(nx);
199     Float_t accErr = hAcceptance->GetBinError(nx);
200     
201     Bool_t etabinSeen = kFALSE;  
202     for(Int_t ny = 1; ny <= aodHist->GetNbinsY(); ny++) {
203       // Get currrent value 
204       Float_t aodValue = aodHist->GetBinContent(nx,ny);
205       Float_t aodErr   = aodHist->GetBinError(nx,ny);
206
207       // Set underflow bin
208       Float_t secCor   = 0;
209       if(hSecMap) secCor   = hSecMap->GetBinContent(nx,ny);
210       if (secCor > 0.5) etabinSeen = kTRUE;
211       if (aodValue < 0.000001) { aodHist->SetBinContent(nx,ny, 0); continue; }
212
213       if (!fUseAcceptance) continue; 
214
215       // Acceptance correction 
216       if (accCor   < 0.000001) accCor = 1;
217       Float_t aodNew   = aodValue / accCor ;
218       Float_t error    = aodNew*TMath::Sqrt(TMath::Power(aodErr/aodValue,2) +
219                                             TMath::Power(accErr/accCor,2) );
220       aodHist->SetBinContent(nx,ny, aodNew);
221       //test
222       aodHist->SetBinError(nx,ny,error);
223       aodHist->SetBinError(nx,ny,aodErr);
224       
225     }
226     //Filling underflow bin if we eta bin is in range
227     if(etabinSeen) aodHist->SetBinContent(nx,0, 1.);
228   }  
229
230   PostData(1,fList);
231 }
232 //____________________________________________________________________
233 void AliCentralMultiplicityTask::Terminate(Option_t* /*option*/) 
234 {
235   // 
236   // End of job
237   // 
238   // Parameters:
239   //    option Not used 
240   //
241 }
242 //____________________________________________________________________
243 void
244 AliCentralMultiplicityTask::Print(Option_t* /*option*/) const
245 {
246   // 
247   // Print information 
248   // 
249   // Parameters:
250   //    option Not used
251   //
252 }
253 //====================================================================
254 AliCentralMultiplicityTask::Manager::Manager() :
255   fAcceptancePath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralAcceptance"),
256   fSecMapPath("$ALICE_ROOT/PWG2/FORWARD/corrections/CentralSecMap"),
257   fAcceptance(),
258   fSecmap(),
259   fAcceptanceName("centralacceptance"),
260   fSecMapName("centralsecmap"),
261   fIsInit(kFALSE)
262 {
263   //
264   // Constructor 
265   // 
266 }
267 //____________________________________________________________________
268 AliCentralMultiplicityTask::Manager::Manager(const Manager& o) 
269   :fAcceptancePath(o.fAcceptancePath),
270    fSecMapPath(o.fSecMapPath),
271    fAcceptance(o.fAcceptance),
272    fSecmap(o.fSecmap),
273    fAcceptanceName(o.fAcceptanceName),
274    fSecMapName(o.fSecMapName),
275    fIsInit(o.fIsInit)
276 {
277   //
278   // Copy Constructor 
279   // 
280 }
281 //____________________________________________________________________
282 AliCentralMultiplicityTask::Manager&
283 AliCentralMultiplicityTask::Manager::operator=(const Manager& o)
284 {
285   //
286   // Assignment operator  
287   // 
288   fAcceptancePath = o.fAcceptancePath;
289   fSecMapPath     = o.fSecMapPath;
290   fAcceptance     = o.fAcceptance;
291   fSecmap         = o.fSecmap;
292   fAcceptanceName = o.fAcceptanceName;
293   fSecMapName     = o.fSecMapName;
294   fIsInit         = o.fIsInit;
295   return *this;
296 }
297
298 //____________________________________________________________________
299 const char* 
300 AliCentralMultiplicityTask::Manager::GetFullFileName(UShort_t what, 
301                                                      UShort_t sys, 
302                                                      UShort_t sNN,  
303                                                      Short_t  field) const
304 {
305   // 
306   // Get full path name to object file 
307   // 
308   // Parameters:
309   //    what   What to get 
310   //    sys    Collision system
311   //    sNN    Center of mass energy 
312   //    field  Magnetic field 
313   // 
314   // Return:
315   //    
316   //
317   return Form("%s/%s",
318               what == 0 ? GetSecMapPath() : GetAcceptancePath(), 
319               GetFileName(what, sys, sNN, field));
320 }
321
322 //____________________________________________________________________
323 const char* 
324 AliCentralMultiplicityTask::Manager::GetFileName(UShort_t  what ,
325                                                  UShort_t  sys, 
326                                                  UShort_t  sNN,
327                                                  Short_t   field) const
328 {
329   // 
330   // Get the full path name 
331   // 
332   // Parameters:
333   //    what   What to get
334   //    sys    Collision system
335   //    sNN    Center of mass energy 
336   //    field  Magnetic field 
337   // 
338   // Return:
339   //    
340   //
341   // Must be static - otherwise the data may disappear on return from
342   // this member function
343   static TString fname = "";
344   
345   fname = "";
346   switch(what) {
347   case 0:  fname.Append(fSecMapName.Data());     break;
348   case 1:  fname.Append(fAcceptanceName.Data()); break;
349   default:
350     ::Error("GetFileName", 
351             "Invalid indentifier %d for central object, must be 0 or 1!", what);
352     break;
353   }
354   fname.Append(Form("_%s_%04dGeV_%c%1dkG.root", 
355                     AliForwardUtil::CollisionSystemString(sys), 
356                     sNN, (field < 0 ? 'm' : 'p'), TMath::Abs(field)));
357   
358   return fname.Data();
359 }
360
361 //____________________________________________________________________
362 TH2D* 
363 AliCentralMultiplicityTask::Manager::GetSecMapCorrection(UShort_t vtxbin) const
364 {
365   // 
366   // Get the secondary map
367   // 
368   // Parameters:
369   //    vtxbin 
370   // 
371   // Return:
372   //    
373   //
374   if (!fSecmap) { 
375     ::Warning("GetSecMapCorrection","No secondary map defined");
376     return 0;
377   }
378   return fSecmap->GetCorrection(vtxbin);
379 }
380 //____________________________________________________________________
381 TH1D* 
382 AliCentralMultiplicityTask::Manager::GetAcceptanceCorrection(UShort_t vtxbin) 
383   const 
384 {
385   // 
386   // Get the acceptance correction 
387   // 
388   // Parameters:
389   //    vtxbin 
390   // 
391   // Return:
392   //    
393   //
394   if (!fAcceptance) { 
395     ::Warning("GetAcceptanceCorrection","No acceptance map defined");
396     return 0;
397   }
398   return fAcceptance->GetCorrection(vtxbin);
399 }
400
401 //____________________________________________________________________
402 void 
403 AliCentralMultiplicityTask::Manager::Init(UShort_t  sys, 
404                                           UShort_t  sNN,
405                                           Short_t   field) 
406 {
407   // 
408   // Initialize 
409   // 
410   // Parameters:
411   //    sys    Collision system (1: pp, 2: PbPb)
412   //    sNN    Center of mass energy per nucleon pair [GeV]
413   //    field  Magnetic field [kG]
414   //
415   if(fIsInit) ::Warning("Init","Already initialised - overriding...");
416   
417   TFile fsec(GetFullFileName(0,sys,sNN,field));
418   fSecmap = 
419     dynamic_cast<AliCentralCorrSecondaryMap*>(fsec.Get(fSecMapName.Data()));  
420   if(!fSecmap) {
421     ::Error("Init", "no central Secondary Map found!") ;
422     return;
423   }
424   TFile facc(GetFullFileName(1,sys,sNN,field));
425   fAcceptance = 
426     dynamic_cast<AliCentralCorrAcceptance*>(facc.Get(fAcceptanceName.Data()));
427   if(!fAcceptance) {
428     ::Error("Init", "no central Acceptance found!") ;
429     return;
430   }
431   
432   if(fSecmap && fAcceptance) {
433     fIsInit = kTRUE;
434     ::Info("Init", 
435            "Central Manager initialised for sys %d, energy %d, field %d",sys,sNN,field);
436
437   }
438   
439 }
440 //
441 // EOF
442 //