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