Add MCNSD flag to trigger mask
[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     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 //____________________________________________________________________
179 const Char_t*
180 AliAODForwardMult::GetTriggerString(UInt_t mask)
181 {
182   // Get a string that describes the triggers 
183   // 
184   // Parameters: 
185   //   mask  Bit pattern of triggers 
186   // Return: 
187   //   Character string representation of mask 
188   static TString trg;
189   trg = "";
190   if ((mask & kInel)        != 0x0) trg.Append("INEL ");
191   if ((mask & kInelGt0)     != 0x0) trg.Append("INEL>0 ");
192   if ((mask & kNSD)         != 0x0) trg.Append("NSD ");
193   if ((mask & kA)           != 0x0) trg.Append("A ");
194   if ((mask & kB)           != 0x0) trg.Append("B ");
195   if ((mask & kC)           != 0x0) trg.Append("C ");
196   if ((mask & kE)           != 0x0) trg.Append("E ");
197   if ((mask & kMCNSD)       != 0x0) trg.Append("MCNSD ");
198   if ((mask & kNClusterGt0) != 0x0) trg.Append("NCluster>0 ");
199   return trg.Data();
200 }
201   
202 //____________________________________________________________________
203 TH1I*
204 AliAODForwardMult::MakeTriggerHistogram(const char* name) 
205 {
206   // 
207   // Make a histogram to record triggers in. 
208   //
209   // The bins defined by the trigger enumeration in this class.  One
210   // can use this enumeration to retrieve the number of triggers for
211   // each class.
212   // 
213   // Parameters:
214   //    name Name of the histogram 
215   // 
216   // Return:
217   //    Newly allocated histogram 
218   //
219   TH1I* ret = new TH1I(name, "Triggers", kAccepted+1, -.5, kAccepted+.5);
220   ret->SetYTitle("Events");
221   ret->SetFillColor(kRed+1);
222   ret->SetFillStyle(3001);
223   ret->GetXaxis()->SetBinLabel(kBinAll,         "All events");
224   ret->GetXaxis()->SetBinLabel(kBinB,           "w/B trigger");
225   ret->GetXaxis()->SetBinLabel(kBinA,           "w/A trigger");
226   ret->GetXaxis()->SetBinLabel(kBinC,           "w/C trigger");
227   ret->GetXaxis()->SetBinLabel(kBinE,           "w/E trigger");
228   ret->GetXaxis()->SetBinLabel(kBinInel,        "Minimum Bias");
229   ret->GetXaxis()->SetBinLabel(kBinInelGt0,     "INEL>0");
230   ret->GetXaxis()->SetBinLabel(kBinNSD,         "NSD");
231   ret->GetXaxis()->SetBinLabel(kBinMCNSD,       "NSD (MC truth)");
232   ret->GetXaxis()->SetBinLabel(kBinPileUp,      "w/Pileup");
233   ret->GetXaxis()->SetBinLabel(kBinOffline,     "w/Offline");
234   ret->GetXaxis()->SetBinLabel(kBinNClusterGt0, "w/N_{cluster}>1");
235   ret->GetXaxis()->SetBinLabel(kWithVertex,     "w/Vertex");
236   ret->GetXaxis()->SetBinLabel(kWithTrigger,    "w/Selected trigger");
237   ret->GetXaxis()->SetBinLabel(kAccepted,       "Accepted by cut");
238   ret->GetXaxis()->SetNdivisions(kAccepted, false);
239   ret->SetStats(0);
240
241   return ret;
242 }
243 //____________________________________________________________________
244 UInt_t 
245 AliAODForwardMult::MakeTriggerMask(const char* what)
246 {
247   UShort_t    trgMask = 0;
248   TString     trgs(what);
249   trgs.ToUpper();
250   TObjString* trg;
251   TIter       next(trgs.Tokenize(" ,|"));
252   while ((trg = static_cast<TObjString*>(next()))) { 
253     TString s(trg->GetString());
254     if      (s.IsNull()) continue;
255     if      (s.CompareTo("INEL")  == 0) trgMask |= AliAODForwardMult::kInel;
256     else if (s.CompareTo("INEL>0")== 0) trgMask |= AliAODForwardMult::kInelGt0;
257     else if (s.CompareTo("NSD")   == 0) trgMask |= AliAODForwardMult::kNSD;
258     else if (s.CompareTo("MCNSD") == 0) trgMask |= AliAODForwardMult::kMCNSD;
259     else if (s.CompareTo("B")     == 0) trgMask |= AliAODForwardMult::kB;
260     else if (s.CompareTo("A")     == 0) trgMask |= AliAODForwardMult::kA;
261     else if (s.CompareTo("C")     == 0) trgMask |= AliAODForwardMult::kC;
262     else if (s.CompareTo("E")     == 0) trgMask |= AliAODForwardMult::kE;
263     else if (s.CompareTo("NCLUSTER>0") == 0) 
264       trgMask |= AliAODForwardMult::kNClusterGt0;
265     else 
266       AliWarningGeneral("MakeTriggerMask", 
267                         Form("Unknown trigger %s", s.Data()));
268   }
269   return trgMask;
270 }
271
272 //____________________________________________________________________
273 Bool_t
274 AliAODForwardMult::CheckEvent(Int_t    triggerMask,
275                               Double_t vzMin, Double_t vzMax,
276                               UShort_t cMin,  UShort_t cMax, 
277                               TH1*     hist) const
278 {
279   // 
280   // Check if event meets the passses requirements.   
281   //
282   // It returns true if @e all of the following is true 
283   //
284   // - The trigger is within the bit mask passed.
285   // - The vertex is within the specified limits. 
286   // - The centrality is within the specified limits, or if lower
287   //   limit is equal to or larger than the upper limit.
288   // 
289   // If a histogram is passed in the last parameter, then that
290   // histogram is filled with the trigger bits. 
291   // 
292   // Parameters:
293   //    triggerMask  Trigger mask
294   //    vzMin        Minimum @f$ v_z@f$ (in centimeters)
295   //    vzMax        Maximum @f$ v_z@f$ (in centimeters) 
296   //    cMin         Minimum centrality (in percent)
297   //    cMax         Maximum centrality (in percent)
298   //    hist         Histogram to fill 
299   // 
300   // Return:
301   //    @c true if the event meets the requirements 
302   //
303   if (cMin < cMax && (cMin > fCentrality || cMax <= fCentrality)) return false;
304
305   if (hist) { 
306     hist->AddBinContent(kBinAll);
307     if (IsTriggerBits(kB|triggerMask))  hist->AddBinContent(kBinB);
308     if (IsTriggerBits(kA|triggerMask))  hist->AddBinContent(kBinA);
309     if (IsTriggerBits(kC|triggerMask))  hist->AddBinContent(kBinC);
310     if (IsTriggerBits(kE|triggerMask))  hist->AddBinContent(kBinE);
311     if (IsTriggerBits(kB|kInel))        hist->AddBinContent(kBinInel);
312     if (IsTriggerBits(kB|kInelGt0))     hist->AddBinContent(kBinInelGt0);
313     if (IsTriggerBits(kB|kNSD))         hist->AddBinContent(kBinNSD);
314     if (IsTriggerBits(kPileUp))         hist->AddBinContent(kBinPileUp);
315     if (IsTriggerBits(kMCNSD))          hist->AddBinContent(kBinMCNSD);
316     if (IsTriggerBits(kOffline))        hist->AddBinContent(kBinOffline);
317     if (IsTriggerBits(kNClusterGt0))    hist->AddBinContent(kBinNClusterGt0);
318   }
319   // Check if we have an event of interest. 
320   Int_t mask = triggerMask; //|kB
321   if (!IsTriggerBits(mask)) return false;
322   
323   // Check for pileup
324   if (IsTriggerBits(kPileUp)) return false;
325   if (hist) hist->AddBinContent(kWithTrigger);
326   
327   // Check that we have a valid vertex
328   if (vzMin < vzMax && !HasIpZ()) return false;
329   if (hist) hist->AddBinContent(kWithVertex);
330
331   // Check that vertex is within cuts 
332   if (vzMin < vzMax && !InRange(vzMin, vzMax)) return false;
333   if (hist) hist->AddBinContent(kAccepted);
334   
335   return true;
336 }
337
338 //____________________________________________________________________
339 void
340 AliAODForwardMult::Print(Option_t* option) const
341 {
342   // Print this object 
343   // 
344   // Parameters: 
345   //  option   Passed to TH1::Print 
346   fHist.Print(option);
347   UShort_t sys = GetSystem();
348   TString  str = "unknown";
349   switch (sys) { 
350   case 1:  str = "pp"; break;
351   case 2:  str = "PbPb"; break;
352   }
353   std::cout << "Ipz:         " << fIpZ << "cm " << (HasIpZ() ? "" : "in") 
354             << "valid\n"
355             << "Triggers:    " << GetTriggerString(fTriggers) 
356             << "sNN:         " << GetSNN() << "GeV\n" 
357             << "System:      " << str 
358             << "Centrality:  " << fCentrality << "%" 
359             << std::endl;
360 }
361
362 //____________________________________________________________________
363 //
364 // EOF
365 //