Better trigger histogram
[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 (Coll.)%s", 
240                                                      andSel.Data()));
241   ret->GetXaxis()->SetBinLabel(kBinA,           Form("A%s", andSel.Data()));
242   ret->GetXaxis()->SetBinLabel(kBinC,           Form("C%s", andSel.Data()));
243   ret->GetXaxis()->SetBinLabel(kBinE,           Form("E%s", andSel.Data()));
244   ret->GetXaxis()->SetBinLabel(kBinInel,        "Coll. & INEL");
245   ret->GetXaxis()->SetBinLabel(kBinInelGt0,     "Coll. & INEL>0");
246   ret->GetXaxis()->SetBinLabel(kBinNSD,         "Coll. & NSD");
247   ret->GetXaxis()->SetBinLabel(kBinV0AND,       "Coll. & V0AND");
248   ret->GetXaxis()->SetBinLabel(kBinMCNSD,       "NSD (MC truth)");
249   ret->GetXaxis()->SetBinLabel(kBinPileUp,      "w/Pileup");
250   ret->GetXaxis()->SetBinLabel(kBinOffline,     "w/Offline");
251   ret->GetXaxis()->SetBinLabel(kBinNClusterGt0, "w/N_{cluster}>1");
252   ret->GetXaxis()->SetBinLabel(kWithVertex,     "w/Vertex");
253   ret->GetXaxis()->SetBinLabel(kWithTrigger,    Form("w/Selected trigger (%s)",
254                                                      sel.Data()));
255   ret->GetXaxis()->SetBinLabel(kAccepted,       "Accepted by cut");
256   ret->GetXaxis()->SetNdivisions(kAccepted, false);
257   ret->SetStats(0);
258
259   return ret;
260 }
261 //____________________________________________________________________
262 UInt_t 
263 AliAODForwardMult::MakeTriggerMask(const char* what)
264 {
265   UShort_t    trgMask = 0;
266   TString     trgs(what);
267   trgs.ToUpper();
268   TObjArray*  parts = trgs.Tokenize("&");
269   TObjString* trg;
270   TIter       next(parts);
271   while ((trg = static_cast<TObjString*>(next()))) { 
272     TString s(trg->GetString());
273     s.Strip(TString::kBoth, ' ');
274     s.ToUpper();
275     if      (s.IsNull()) continue;
276     if      (s.CompareTo("INEL")  == 0) trgMask |= AliAODForwardMult::kInel;
277     else if (s.CompareTo("INEL>0")== 0) trgMask |= AliAODForwardMult::kInelGt0;
278     else if (s.CompareTo("NSD")   == 0) trgMask |= AliAODForwardMult::kNSD;
279     else if (s.CompareTo("V0AND") == 0) trgMask |= AliAODForwardMult::kV0AND;
280     else if (s.CompareTo("MCNSD") == 0) trgMask |= AliAODForwardMult::kMCNSD;
281     else if (s.CompareTo("B")     == 0) trgMask |= AliAODForwardMult::kB;
282     else if (s.CompareTo("A")     == 0) trgMask |= AliAODForwardMult::kA;
283     else if (s.CompareTo("C")     == 0) trgMask |= AliAODForwardMult::kC;
284     else if (s.CompareTo("E")     == 0) trgMask |= AliAODForwardMult::kE;
285     else if (s.CompareTo("NCLUSTER>0") == 0) 
286       trgMask |= AliAODForwardMult::kNClusterGt0;
287     else 
288       AliWarningGeneral("MakeTriggerMask", 
289                         Form("Unknown trigger %s", s.Data()));
290   }
291   parts->Delete();
292   return trgMask;
293 }
294
295 //____________________________________________________________________
296 Bool_t
297 AliAODForwardMult::CheckEvent(Int_t    triggerMask,
298                               Double_t vzMin, Double_t vzMax,
299                               UShort_t cMin,  UShort_t cMax, 
300                               TH1*     hist) const
301 {
302   // 
303   // Check if event meets the passses requirements.   
304   //
305   // It returns true if @e all of the following is true 
306   //
307   // - The trigger is within the bit mask passed.
308   // - The vertex is within the specified limits. 
309   // - The centrality is within the specified limits, or if lower
310   //   limit is equal to or larger than the upper limit.
311   // 
312   // If a histogram is passed in the last parameter, then that
313   // histogram is filled with the trigger bits. 
314   // 
315   // Parameters:
316   //    triggerMask  Trigger mask
317   //    vzMin        Minimum @f$ v_z@f$ (in centimeters)
318   //    vzMax        Maximum @f$ v_z@f$ (in centimeters) 
319   //    cMin         Minimum centrality (in percent)
320   //    cMax         Maximum centrality (in percent)
321   //    hist         Histogram to fill 
322   // 
323   // Return:
324   //    @c true if the event meets the requirements 
325   //
326   if (cMin < cMax && (cMin > fCentrality || cMax <= fCentrality)) return false;
327
328   if (hist) { 
329     Int_t tmp = triggerMask & ~kB;
330     hist->AddBinContent(kBinAll);
331     if (IsTriggerBits(kB|tmp))          hist->AddBinContent(kBinB);
332     if (IsTriggerBits(kA|tmp))          hist->AddBinContent(kBinA);
333     if (IsTriggerBits(kC|tmp))          hist->AddBinContent(kBinC);
334     if (IsTriggerBits(kE|tmp))          hist->AddBinContent(kBinE);
335     if (IsTriggerBits(kB|kInel))        hist->AddBinContent(kBinInel);
336     if (IsTriggerBits(kB|kInelGt0))     hist->AddBinContent(kBinInelGt0);
337     if (IsTriggerBits(kB|kNSD))         hist->AddBinContent(kBinNSD);
338     if (IsTriggerBits(kB|kV0AND))       hist->AddBinContent(kBinV0AND);
339     if (IsTriggerBits(kPileUp))         hist->AddBinContent(kBinPileUp);
340     if (IsTriggerBits(kMCNSD))          hist->AddBinContent(kBinMCNSD);
341     if (IsTriggerBits(kOffline))        hist->AddBinContent(kBinOffline);
342     if (IsTriggerBits(kNClusterGt0))    hist->AddBinContent(kBinNClusterGt0);
343     if (IsTriggerBits(triggerMask) && !IsTriggerBits(kB|tmp))
344       Warning("CheckEvent", "event: 0x%x, mask: 0x%x, tmp: 0x%x, tmp|b: 0x%x",
345              fTriggers, triggerMask, tmp, tmp|kB);
346   }
347   // Check if we have an event of interest. 
348   Int_t mask = triggerMask; //|kB
349   if (!IsTriggerBits(mask)) return false;
350   
351   // Check for pileup
352   if (IsTriggerBits(kPileUp)) return false;
353   if (hist) hist->AddBinContent(kWithTrigger);
354   
355   // Check that we have a valid vertex
356   if (vzMin < vzMax && !HasIpZ()) return false;
357   if (hist) hist->AddBinContent(kWithVertex);
358
359   // Check that vertex is within cuts 
360   if (vzMin < vzMax && !InRange(vzMin, vzMax)) return false;
361   if (hist) hist->AddBinContent(kAccepted);
362   
363   return true;
364 }
365
366 //____________________________________________________________________
367 void
368 AliAODForwardMult::Print(Option_t* option) const
369 {
370   // Print this object 
371   // 
372   // Parameters: 
373   //  option   Passed to TH1::Print 
374   fHist.Print(option);
375   UShort_t sys = GetSystem();
376   TString  str = "unknown";
377   switch (sys) { 
378   case 1:  str = "pp"; break;
379   case 2:  str = "PbPb"; break;
380   case 3:  str = "pPb" ; break;
381   }
382   std::cout << "Ipz:         " << fIpZ << "cm " << (HasIpZ() ? "" : "in") 
383             << "valid\n"
384             << "Triggers:    " << GetTriggerString(fTriggers)  << "\n"
385             << "sNN:         " << GetSNN() << "GeV\n" 
386             << "System:      " << str << "\n"
387             << "Centrality:  " << fCentrality << "%" 
388             << std::endl;
389 }
390
391 //____________________________________________________________________
392 //
393 // EOF
394 //