Fixes
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliAODForwardMult.cxx
1 //
2 // Class that contains the forward multiplicity data per event 
3 //
4 // This class contains a histogram of 
5 // @f[
6 //   \frac{d^2N_{ch}}{d\eta d\phi}\quad,
7 // @f]
8 // as well as a trigger mask for each analysed event.  
9 // 
10 // The eta acceptance of the event is stored in the underflow bins of
11 // the histogram.  So to build the final histogram, one needs to
12 // correct for this acceptance (properly weighted by the events), and
13 // the vertex efficiency.  This simply boils down to defining a 2D
14 // histogram and summing the event histograms in that histogram.  One
15 // should of course also do proper book-keeping of the accepted event.
16 //
17 #include "AliAODForwardMult.h"
18 #include <TBrowser.h>
19 #include <iostream>
20 #include <TMath.h>
21 #include <TObjString.h>
22 #include <TObjArray.h>
23 #include "AliLog.h"
24 ClassImp(AliAODForwardMult)
25 #ifdef DOXY_INPUT
26 ; // For Emacs 
27 #endif
28
29 //____________________________________________________________________
30 const Float_t AliAODForwardMult::fgkInvalidIpZ = 1e6;
31
32 //____________________________________________________________________
33 AliAODForwardMult::AliAODForwardMult()
34   : fIsMC(false),
35     fHist(),
36     fTriggers(0),
37     fIpZ(fgkInvalidIpZ), 
38     fCentrality(-1)
39 {
40   // 
41   // Constructor 
42   // 
43 }
44
45 //____________________________________________________________________
46 AliAODForwardMult::AliAODForwardMult(Bool_t isMC) 
47   : fIsMC(isMC),
48     fHist("forwardMult", "d^{2}N_{ch}/d#etad#varphi in the forward regions", 
49           200, -4, 6, 20, 0, 2*TMath::Pi()),
50     fTriggers(0),
51     fIpZ(fgkInvalidIpZ), 
52     fCentrality(-1)
53 {
54   // 
55   // Constructor 
56   // 
57   // Parameters: 
58   //  isMC   If set to true this is for MC data (effects branch name)
59   // 
60   fHist.SetXTitle("#eta");
61   fHist.SetYTitle("#varphi [radians]");
62   fHist.SetZTitle("#frac{d^{2}N_{ch}}{d#etad#varphi}");
63   fHist.SetDirectory(0);
64   fHist.Sumw2();
65 }
66
67 //____________________________________________________________________
68 void
69 AliAODForwardMult::Init(const TAxis& etaAxis)
70 {
71   // Initialize the histogram with an eta axis 
72   // 
73   // Parameters: 
74   //   etaAxis       Eta axis to use 
75   // 
76   fHist.SetBins(etaAxis.GetNbins(), etaAxis.GetXmin(), etaAxis.GetXmax(), 
77                 20, 0, 2*TMath::Pi());
78 }
79
80 //____________________________________________________________________
81 void
82 AliAODForwardMult::Clear(Option_t* option)
83 {
84   // Clear (or reset) internal values 
85   // 
86   // Parameters: 
87   //  option   Passed to TH1::Reset 
88   // 
89   fHist.Reset(option);
90   fTriggers = 0;
91   fIpZ      = fgkInvalidIpZ;
92 }
93 //____________________________________________________________________
94 void
95 AliAODForwardMult::SetSNN(UShort_t snn)
96 {
97   // set the center of mass energy per nucleon pair (GeV). 
98   // This is stored in bin (0,0) of the histogram 
99   // 
100   // Parameters: 
101   //   sNN   Center of mass energy per nuclean 
102   fHist.SetBinContent(0,0,snn);
103 }
104 //____________________________________________________________________
105 void
106 AliAODForwardMult::SetSystem(UShort_t sys)
107 {
108   // set the center of mass energy per nucleon pair (GeV). 
109   // This is stored in bin (N+1,0) of the histogram 
110   // 
111   // Parameters: 
112   //   sys   Collision system number 
113   fHist.SetBinContent(fHist.GetNbinsX()+1,0,sys);
114 }
115
116 //____________________________________________________________________
117 Bool_t
118 AliAODForwardMult::HasIpZ() const
119 {
120   // Check if we have valid z coordinate of the interaction point 
121   // 
122   // Return:
123   //   true if the z coordinate of the interaction point is valid 
124   // 
125   return TMath::Abs(fIpZ - fgkInvalidIpZ) > 1;
126 }
127
128 //____________________________________________________________________
129 UShort_t
130 AliAODForwardMult::GetSNN() const
131 {
132   // set the center of mass energy per nucleon pair (GeV). 
133   // This is stored in bin (0,0) of the histogram 
134   // 
135   // Parameters: 
136   //   sNN   Center of mass energy per nuclean 
137   return UShort_t(fHist.GetBinContent(0,0));
138 }
139
140 //____________________________________________________________________
141 UShort_t
142 AliAODForwardMult::GetSystem() const
143 {
144   // set the center of mass energy per nucleon pair (GeV). 
145   // This is stored in bin (N+1,0) of the histogram 
146   // 
147   // Parameters: 
148   //   sNN   Center of mass energy per nuclean 
149   return UShort_t(fHist.GetBinContent(fHist.GetNbinsX()+1,0));
150 }
151
152 //____________________________________________________________________
153 void
154 AliAODForwardMult::Browse(TBrowser* b)
155 {
156   // Browse this object 
157   // 
158   // Parameters: 
159   //   b   Browser to use 
160   static TObjString ipz;
161   static TObjString trg;
162   static TObjString cnt;
163   ipz = Form("ip_z=%fcm", fIpZ);
164   trg = GetTriggerString(fTriggers);
165   cnt = Form("%+6.1f%%", fCentrality);
166   b->Add(&fHist);
167   b->Add(&ipz);
168   b->Add(&trg);
169   b->Add(&cnt);
170 }
171
172 //____________________________________________________________________
173 const Char_t*
174 AliAODForwardMult::GetTriggerString(UInt_t mask)
175 {
176   // Get a string that describes the triggers 
177   // 
178   // Parameters: 
179   //   mask  Bit pattern of triggers 
180   // Return: 
181   //   Character string representation of mask 
182   static TString trg;
183   trg = "";
184   if ((mask & kInel)        != 0x0) trg.Append("INEL ");
185   if ((mask & kInelGt0)     != 0x0) trg.Append("INEL>0 ");
186   if ((mask & kNSD)         != 0x0) trg.Append("NSD ");
187   if ((mask & kA)           != 0x0) trg.Append("A ");
188   if ((mask & kB)           != 0x0) trg.Append("B ");
189   if ((mask & kC)           != 0x0) trg.Append("C ");
190   if ((mask & kE)           != 0x0) trg.Append("E ");
191   if ((mask & kMCNSD)       != 0x0) trg.Append("MCNSD ");
192   return trg.Data();
193 }
194   
195 //____________________________________________________________________
196 TH1I*
197 AliAODForwardMult::MakeTriggerHistogram(const char* name) 
198 {
199   // 
200   // Make a histogram to record triggers in. 
201   //
202   // The bins defined by the trigger enumeration in this class.  One
203   // can use this enumeration to retrieve the number of triggers for
204   // each class.
205   // 
206   // Parameters:
207   //    name Name of the histogram 
208   // 
209   // Return:
210   //    Newly allocated histogram 
211   //
212   TH1I* ret = new TH1I(name, "Triggers", kAccepted+1, -.5, kAccepted+.5);
213   ret->SetYTitle("Events");
214   ret->SetFillColor(kRed+1);
215   ret->SetFillStyle(3001);
216   ret->GetXaxis()->SetBinLabel(kBinAll,         "All events");
217   ret->GetXaxis()->SetBinLabel(kBinB,           "w/B trigger");
218   ret->GetXaxis()->SetBinLabel(kBinA,           "w/A trigger");
219   ret->GetXaxis()->SetBinLabel(kBinC,           "w/C trigger");
220   ret->GetXaxis()->SetBinLabel(kBinE,           "w/E trigger");
221   ret->GetXaxis()->SetBinLabel(kBinInel,        "Minimum Bias");
222   ret->GetXaxis()->SetBinLabel(kBinInelGt0,     "INEL>0");
223   ret->GetXaxis()->SetBinLabel(kBinNSD,         "NSD");
224   ret->GetXaxis()->SetBinLabel(kBinMCNSD,       "NSD (MC truth)");
225   ret->GetXaxis()->SetBinLabel(kBinPileUp,      "w/Pileup");
226   ret->GetXaxis()->SetBinLabel(kBinOffline,     "w/Offline");
227   ret->GetXaxis()->SetBinLabel(kWithVertex,     "w/Vertex");
228   ret->GetXaxis()->SetBinLabel(kWithTrigger,    "w/Selected trigger");
229   ret->GetXaxis()->SetBinLabel(kAccepted,       "Accepted by cut");
230   ret->GetXaxis()->SetNdivisions(kAccepted, false);
231   ret->SetStats(0);
232
233   return ret;
234 }
235 //____________________________________________________________________
236 UInt_t 
237 AliAODForwardMult::MakeTriggerMask(const char* what)
238 {
239   UShort_t    trgMask = 0;
240   TString     trgs(what);
241   trgs.ToUpper();
242   TObjString* trg;
243   TIter       next(trgs.Tokenize(" ,|"));
244   while ((trg = static_cast<TObjString*>(next()))) { 
245     TString s(trg->GetString());
246     if      (s.IsNull()) continue;
247     if      (s.CompareTo("INEL")  == 0) trgMask |= AliAODForwardMult::kInel;
248     else if (s.CompareTo("INEL>0")== 0) trgMask |= AliAODForwardMult::kInelGt0;
249     else if (s.CompareTo("NSD")   == 0) trgMask |= AliAODForwardMult::kNSD;
250     else if (s.CompareTo("B")     == 0) trgMask |= AliAODForwardMult::kB;
251     else if (s.CompareTo("A")     == 0) trgMask |= AliAODForwardMult::kA;
252     else if (s.CompareTo("C")     == 0) trgMask |= AliAODForwardMult::kC;
253     else if (s.CompareTo("E")     == 0) trgMask |= AliAODForwardMult::kE;
254     else 
255       AliWarningGeneral("MakeTriggerMask", 
256                         Form("Unknown trigger %s", s.Data()));
257   }
258   return trgMask;
259 }
260
261 //____________________________________________________________________
262 Bool_t
263 AliAODForwardMult::CheckEvent(Int_t    triggerMask,
264                               Double_t vzMin, Double_t vzMax,
265                               UShort_t cMin,  UShort_t cMax, 
266                               TH1*     hist) const
267 {
268   // 
269   // Check if event meets the passses requirements.   
270   //
271   // It returns true if @e all of the following is true 
272   //
273   // - The trigger is within the bit mask passed.
274   // - The vertex is within the specified limits. 
275   // - The centrality is within the specified limits, or if lower
276   //   limit is equal to or larger than the upper limit.
277   // 
278   // If a histogram is passed in the last parameter, then that
279   // histogram is filled with the trigger bits. 
280   // 
281   // Parameters:
282   //    triggerMask  Trigger mask
283   //    vzMin        Minimum @f$ v_z@f$ (in centimeters)
284   //    vzMax        Maximum @f$ v_z@f$ (in centimeters) 
285   //    cMin         Minimum centrality (in percent)
286   //    cMax         Maximum centrality (in percent)
287   //    hist         Histogram to fill 
288   // 
289   // Return:
290   //    @c true if the event meets the requirements 
291   //
292   if (cMin < cMax && (cMin > fCentrality || cMax <= fCentrality)) return false;
293
294   if (hist) { 
295     hist->AddBinContent(kBinAll);
296     if (IsTriggerBits(kB|triggerMask))  hist->AddBinContent(kBinB);
297     if (IsTriggerBits(kA|triggerMask))  hist->AddBinContent(kBinA);
298     if (IsTriggerBits(kC|triggerMask))  hist->AddBinContent(kBinC);
299     if (IsTriggerBits(kE|triggerMask))  hist->AddBinContent(kBinE);
300     if (IsTriggerBits(kB|kInel))        hist->AddBinContent(kBinInel);
301     if (IsTriggerBits(kB|kInelGt0))     hist->AddBinContent(kBinInelGt0);
302     if (IsTriggerBits(kB|kNSD))         hist->AddBinContent(kBinNSD);
303     if (IsTriggerBits(kPileUp))         hist->AddBinContent(kBinPileUp);
304     if (IsTriggerBits(kMCNSD))          hist->AddBinContent(kBinMCNSD);
305     if (IsTriggerBits(kOffline))        hist->AddBinContent(kBinOffline);
306   }
307   // Check if we have an event of interest. 
308   if (!IsTriggerBits(triggerMask|kB)) return false;
309   
310   // Check for pileup
311   if (IsTriggerBits(kPileUp)) return false;
312   if (hist) hist->AddBinContent(kWithTrigger);
313   
314   // Check that we have a valid vertex
315   if (vzMin < vzMax && !HasIpZ()) return false;
316   if (hist) hist->AddBinContent(kWithVertex);
317
318   // Check that vertex is within cuts 
319   if (vzMin < vzMax && !InRange(vzMin, vzMax)) return false;
320   if (hist) hist->AddBinContent(kAccepted);
321   
322   return true;
323 }
324
325 //____________________________________________________________________
326 void
327 AliAODForwardMult::Print(Option_t* option) const
328 {
329   // Print this object 
330   // 
331   // Parameters: 
332   //  option   Passed to TH1::Print 
333   fHist.Print(option);
334   UShort_t sys = GetSystem();
335   TString  str = "unknown";
336   switch (sys) { 
337   case 1:  str = "pp"; break;
338   case 2:  str = "PbPb"; break;
339   }
340   std::cout << "Ipz:         " << fIpZ << "cm " << (HasIpZ() ? "" : "in") 
341             << "valid\n"
342             << "Triggers:    " << GetTriggerString(fTriggers) 
343             << "sNN:         " << GetSNN() << "GeV\n" 
344             << "System:      " << str 
345             << "Centrality:  " << fCentrality << "%" 
346             << std::endl;
347 }
348
349 //____________________________________________________________________
350 //
351 // EOF
352 //