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