Better trigger histogram
[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliAODForwardMult.h
1 //
2 // See implementation or Doxygen comments for more information
3 //
4 #ifndef ALIAODFORWARDMULT_H
5 #define ALIAODFORWARDMULT_H
6 /**
7  * @file   AliAODForwardMult.h
8  * @author Christian Holm Christensen <cholm@dalsgaard.hehi.nbi.dk>
9  * @date   Wed Mar 23 13:58:00 2011
10  * 
11  * @brief  
12  * 
13  * @ingroup pwglf_forward_aod
14  * 
15  */
16 #include <TObject.h>
17 #include <TH2D.h>
18 class TBrowser;
19 class TH1I;
20 /**
21  * Class that contains the forward multiplicity data per event 
22  *
23  * This class contains a histogram of 
24  * @f[
25  *   \frac{d^2N_{ch}}{d\eta d\phi}\quad,
26  * @f]
27  * as well as a trigger mask for each analysed event.  
28  * 
29  * The eta acceptance of the event is stored in the underflow bins of
30  * the histogram.  So to build the final histogram, one needs to
31  * correct for this acceptance (properly weighted by the events), and
32  * the vertex efficiency.  This simply boils down to defining a 2D
33  * histogram and summing the event histograms in that histogram.  One
34  * should of course also do proper book-keeping of the accepted event.
35  *
36  * @code 
37  * TTree* GetAODTree()
38  * { 
39  *    TFile* file = TFile::Open("AliAODs.root","READ");
40  *    TTree* tree = static_cast<TTree*>(file->Get("aodTree"));
41  *    return tree;
42  * }
43  * 
44  * void Analyse()
45  * { 
46  *   TH2D*              sum        = 0;                  // Summed hist
47  *   TTree*             tree       = GetAODTree();       // AOD tree
48  *   AliAODForwardMult* mult       = 0;                  // AOD object
49  *   Int_t              nTriggered = 0;                  // # of triggered ev.
50  *   Int_t              nWithVertex= 0;                  // # of ev. w/vertex
51  *   Int_t              nAccepted  = 0;                  // # of ev. used
52  *   Int_t              nAvailable = tree->GetEntries(); // How many entries
53  *   Float_t            vzLow      = -10;                // Lower ip cut
54  *   Float_t            vzHigh     =  10;                // Upper ip cut
55  *   Int_t              mask       = AliAODForwardMult::kInel;// Trigger mask
56  *   tree->SetBranchAddress("forward", &forward);        // Set the address
57  * 
58  *   for (int i = 0; i < nAvailable; i++) { 
59  *     // Create sum histogram on first event - to match binning to input
60  *     if (!sum) sum = static_cast<TH2D*>(mult->Clone("d2ndetadphi"));
61  * 
62  *     tree->GetEntry(i);
63  * 
64  *     // Other trigger/event requirements could be defined 
65  *     if (!mult->IsTriggerBits(mask)) continue; 
66  *     nTriggered++;
67  *
68  *     // Check if we have vertex 
69  *     if (!mult->HasIpZ()) continue;
70  *     nWithVertex++;
71  * 
72  *     // Select vertex range (in centimeters) 
73  *     if (!mult->InRange(vzLow, vzHigh) continue; 
74  *     nAccepted++;
75  * 
76  *     // Add contribution from this event
77  *     sum->Add(&(mult->GetHistogram()));
78  *   }
79  * 
80  *   // Get acceptance normalisation from underflow bins 
81  *   TH1D* norm   = sum->ProjectionX("norm", 0, 1, "");
82  *   // Project onto eta axis - _ignoring_underflow_bins_!
83  *   TH1D* dndeta = sum->Projection("dndeta", 1, -1, "e");
84  *   // Normalize to the acceptance 
85  *   dndeta->Divide(norm);
86  *   // Scale by the vertex efficiency 
87  *   dndeta->Scale(Double_t(nWithVertex)/nTriggered, "width");
88  *   // And draw the result
89  *   dndeta->Draw();
90  * }
91  * @endcode   
92  *     
93  * The above code will draw the final @f$ dN_{ch}/d\eta@f$ for the
94  * selected event class and vertex range
95  *
96  * The histogram can be used as input for other kinds of analysis too, 
97  * like flow, event-plane, centrality, and so on. 
98  *
99  * @ingroup pwglf_forward 
100  * @ingroup pwglf_forward_aod
101  */
102 class AliAODForwardMult : public TObject
103 {
104 public:
105   /** 
106    * Bits of the trigger pattern
107    */
108   enum { 
109     /** In-elastic collision */
110     kInel     = 0x001, 
111     /** In-elastic collision with at least one SPD tracklet */
112     kInelGt0  = 0x002, 
113     /** Non-single diffractive collision */
114     kNSD      = 0x004, 
115     /** Empty bunch crossing */
116     kEmpty    = 0x008, 
117     /** A-side trigger */
118     kA        = 0x010, 
119     /** B(arrel) trigger */
120     kB        = 0x020, 
121     /** C-side trigger */
122     kC        = 0x080,  
123     /** Empty trigger */
124     kE        = 0x100,
125     /** pileup from SPD */
126     kPileUp   = 0x200,    
127     /** true NSD from MC */
128     kMCNSD    = 0x400,    
129     /** Offline MB triggered */
130     kOffline  = 0x800,
131     /** At least one SPD cluster */ 
132     kNClusterGt0 = 0x1000,
133     /** V0-AND trigger */
134     kV0AND       = 0x2000
135   };
136   /** 
137    * Bin numbers in trigger histograms 
138    */
139   enum { 
140     kBinAll=1,
141     kBinInel, 
142     kBinInelGt0, 
143     kBinNSD, 
144     kBinV0AND,
145     kBinA, 
146     kBinB, 
147     kBinC, 
148     kBinE,
149     kBinPileUp, 
150     kBinMCNSD,
151     kBinOffline,
152     kBinNClusterGt0,
153     kWithTrigger, 
154     kWithVertex, 
155     kAccepted
156   };
157   /** 
158    * Default constructor 
159    * 
160    * Used by ROOT I/O sub-system - do not use
161    */
162   AliAODForwardMult();
163   /** 
164    * Constructor 
165    * 
166    * @param isMC Whether this was from MC or not 
167    */
168   AliAODForwardMult(Bool_t isMC);
169   /** 
170    * Destructor 
171    */
172   virtual ~AliAODForwardMult() {} // Destructor 
173   /** 
174    * Initialize 
175    * 
176    * @param etaAxis  Pseudo-rapidity axis
177    */
178   void Init(const TAxis& etaAxis);
179   /** 
180    * Get the @f$ d^2N_{ch}/d\eta d\phi@f$ histogram, 
181    *
182    * @return @f$ d^2N_{ch}/d\eta d\phi@f$ histogram, 
183    */  
184   const TH2D& GetHistogram() const { return fHist; } // Get histogram 
185   /** 
186    * Get the @f$ d^2N_{ch}/d\eta d\phi@f$ histogram, 
187    *
188    * @return @f$ d^2N_{ch}/d\eta d\phi@f$ histogram, 
189    */  
190   TH2D& GetHistogram() { return fHist; } // Get histogram 
191   /** 
192    * Get the trigger bits 
193    * 
194    * @return Trigger bits 
195    */
196   UInt_t GetTriggerBits() const { return fTriggers; } // Get triggers
197   /** 
198    * Set the trigger mask 
199    * 
200    * @param trg Trigger mask
201    */
202   void SetTriggerMask(UInt_t trg) { fTriggers = trg; } // Set triggers 
203   /** 
204    * Set bit(s) in trigger mask 
205    * 
206    * @param bits bit(s) to set 
207    */
208   void SetTriggerBits(UInt_t bits) { fTriggers |= bits; } // Set trigger bits
209   /** 
210    * Check if all bit(s) are set in the trigger mask.  Note, this is
211    * an @e and between the bits.  If you need an @e or you should use
212    * the member function IsTriggerOrBits
213    * 
214    * @param bits Bits to test for 
215    * 
216    * @return true if all enabled bits in the argument is also set in
217    * the trigger word
218    */
219   Bool_t IsTriggerBits(UInt_t bits) const;
220   /** 
221    * Check if any of bit(s) are enabled in the trigger word.  This is
222    * an @e or between the selected bits.  If you need and @a and you
223    * should use the member function IsTriggerBits;
224    * 
225    * @param bits Bits to check for 
226    * 
227    * @return true if any of the enabled bits in the arguments are also
228    * enabled in the trigger mask
229    */
230   Bool_t IsTriggerOrBits(UInt_t bits) const;
231   /** 
232    * Whether we have any trigger bits 
233    *
234    * @return true if we have some trigger 
235    */
236   Bool_t HasTrigger() const { return fTriggers != 0; } // Check for triggers
237   /** 
238    * Clear all data 
239    * 
240    * @param option  Passed on to TH2::Reset verbatim
241    */
242   void Clear(Option_t* option="");
243   /** 
244    * browse this object 
245    * 
246    * @param b Browser 
247    */
248   void Browse(TBrowser* b);
249   /** 
250    * This is a folder 
251    * 
252    * @return Always true
253    */
254   Bool_t IsFolder() const { return kTRUE; } // Always true 
255   /** 
256    * Print content 
257    * 
258    * @param option Passed verbatim to TH2::Print 
259    */
260   void Print(Option_t* option="") const;
261   /** 
262    * Set the z coordinate of the interaction point
263    * 
264    * @param ipZ Interaction point z coordinate
265    */
266   void SetIpZ(Float_t ipZ) { fIpZ = ipZ; } // Set Ip's Z coordinate
267   /** 
268    * Set the center of mass energy per nucleon-pair.  This is stored 
269    * in the (0,0) of the histogram 
270    * 
271    * @param sNN Center of mass energy per nucleon pair (GeV)
272    */
273   void SetSNN(UShort_t sNN); 
274   /** 
275    * Get the collision system number
276    * - 0: Unknown 
277    * - 1: pp
278    * - 2: PbPb
279    * 
280    * @param sys Collision system number
281    */
282   void SetSystem(UShort_t sys);
283   /** 
284    * Set the event centrality 
285    * 
286    * @param c Centrality 
287    */
288   void SetCentrality(Float_t c) { fCentrality = c; }
289   /** 
290    * Set the z coordinate of the interaction point
291    * 
292    * @return Interaction point z coordinate
293    */
294   Float_t GetIpZ() const { return fIpZ; } // Get Ip's Z coordinate 
295   /** 
296    * Check if we have a valid z coordinate of the interaction point
297    *
298    * @return True if we have a valid interaction point z coordinate
299    */
300   Bool_t HasIpZ() const;
301   /** 
302    * Get the center of mass energy per nucleon pair (GeV)
303    * 
304    * @return Center of mass energy per nucleon pair (GeV)
305    */
306   UShort_t GetSNN() const;
307   /** 
308    * Get the collision system number
309    * - 0: Unknown 
310    * - 1: pp
311    * - 2: PbPb
312    * 
313    * @return Collision system number
314    */
315   UShort_t GetSystem() const;
316   /** 
317    * Check if the z coordinate of the interaction point is within the
318    * given limits.  Note that the convention used corresponds to the
319    * convention used in ROOTs TAxis.
320    * 
321    * @param low  Lower cut (inclusive)
322    * @param high Upper cut (exclusive)
323    * 
324    * @return true if @f$ low \ge ipz < high@f$ 
325    */
326   Bool_t InRange(Float_t low, Float_t high) const;
327   /** 
328    * Get the event centrality 
329    * 
330    * 
331    * @return 
332    */
333   Float_t GetCentrality() const { return fCentrality; }
334   /** 
335    * Check if we have a valid centrality 
336    * 
337    * 
338    * @return 
339    */
340   Bool_t  HasCentrality() const { return !(fCentrality  < 0); }
341   /** 
342    * Get the number of SPD clusters seen in @f$ |\eta|<1@f$ 
343    * 
344    * @return Number of SPD clusters seen
345    */
346   UShort_t GetNClusters() const { return fNClusters; }
347   /** 
348    * Set the number of SPD clusters seen in @f$ |\eta|<1@f$ 
349    * 
350    * @param n Number of SPD clusters 
351    */
352   void SetNClusters(UShort_t n) { fNClusters = n; }
353   /** 
354    * Get the name of the object 
355    * 
356    * @return Name of object 
357    */
358   const Char_t* GetName() const { return (fIsMC ? "ForwardMC" : "Forward"); }
359   /** 
360    * Check if event meets the passses requirements.   
361    *
362    * It returns true if @e all of the following is true 
363    *
364    * - The trigger is within the bit mask passed.
365    * - The vertex is within the specified limits. 
366    * - The centrality is within the specified limits, or if lower
367    *   limit is equal to or larger than the upper limit.
368    * 
369    * Note, for data with out a centrality estimate (e.g., pp), one
370    * must pass equal centrality cuts, or no data will be accepted.  In
371    * other words, for pp data, always pass cMin=0, cMax=0
372    *
373    * If a histogram is passed in the last parameter, then that
374    * histogram is filled with the trigger bits. 
375    * 
376    * @param triggerMask  Trigger mask
377    * @param vzMin        Minimum @f$ v_z@f$ (in centimeters)
378    * @param vzMax        Maximum @f$ v_z@f$ (in centimeters) 
379    * @param cMin         Minimum centrality (in percent)
380    * @param cMax         Maximum centrality (in percent)
381    * @param hist         Histogram to fill 
382    * 
383    * @return @c true if the event meets the requirements 
384    */
385   Bool_t CheckEvent(Int_t    triggerMask=kInel,
386                     Double_t vzMin=-10, Double_t vzMax=10,
387                     UShort_t cMin=0,    UShort_t cMax=100, 
388                     TH1*     hist=0) const;
389   /** 
390    * Get a string correspondig to the trigger mask
391    * 
392    * @param mask Trigger mask 
393    * 
394    * @return Static string (copy before use)
395    */
396   static const Char_t* GetTriggerString(UInt_t mask);
397   /** 
398    * Make a histogram to record triggers in. 
399    *
400    * The bins defined by the trigger enumeration in this class.  One
401    * can use this enumeration to retrieve the number of triggers for
402    * each class.
403    * 
404    * @param name Name of the histogram 
405    * 
406    * @return Newly allocated histogram 
407    */
408   static TH1I* MakeTriggerHistogram(const char* name="triggers",
409                                     Int_t mask=0);
410   /** 
411    * Utility function to make a trigger mask from the passed string. 
412    * 
413    * The string is a comma or space seperated list of case-insensitive
414    * strings
415    * 
416    * - INEL 
417    * - INEL>0
418    * - NSD 
419    * 
420    * @param what Which triggers to put in the mask. 
421    * 
422    * @return The generated trigger mask. 
423    */
424   static UInt_t MakeTriggerMask(const char* what);
425 protected: 
426   /** From MC or not */
427   Bool_t   fIsMC;       // Whether this is from MC 
428   /** Histogram of @f$d^2N_{ch}/(d\eta d\phi)@f$ for this event */
429   TH2D     fHist;       // Histogram of d^2N_{ch}/(deta dphi) for this event
430   /** Trigger bits */
431   UInt_t   fTriggers;   // Trigger bit mask 
432   /** Interaction point @f$z@f$ coordinate */
433   Float_t  fIpZ;        // Z coordinate of the interaction point
434   /** Centrality */
435   Float_t  fCentrality; // Event centrality 
436   /** Number of clusters in @f$|\eta|<1@f$ */
437   UShort_t fNClusters;  // Number of SPD clusters in |eta|<1
438   /** Invalid value for interaction point @f$z@f$ coordiante */
439   static const Float_t fgkInvalidIpZ; // Invalid IpZ value 
440   ClassDef(AliAODForwardMult,3); // AOD forward multiplicity 
441 };
442
443 //____________________________________________________________________
444 inline Bool_t
445 AliAODForwardMult::InRange(Float_t low, Float_t high) const 
446 {
447   return HasIpZ() && fIpZ >= low && fIpZ < high;
448 }
449
450 //____________________________________________________________________
451 inline Bool_t 
452 AliAODForwardMult::IsTriggerBits(UInt_t bits) const 
453
454   return HasTrigger() && ((fTriggers & bits) == bits); 
455 }
456 //____________________________________________________________________
457 inline Bool_t 
458 AliAODForwardMult::IsTriggerOrBits(UInt_t bits) const 
459
460   return HasTrigger() && ((fTriggers & bits) != 0);
461 }
462
463 #endif
464 // Local Variables:
465 //  mode: C++
466 // End:
467