38d20e6154ebc1c3ac6fafa5719d31fe71d42d45
[u/mrichter/AliRoot.git] / PWG2 / 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 pwg2_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->Projection("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 pwg2_forward 
100  * @ingroup pwg2_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   };
132   /** 
133    * Bin numbers in trigger histograms 
134    */
135   enum { 
136     kBinAll=1,
137     kBinInel, 
138     kBinInelGt0, 
139     kBinNSD, 
140     kBinA, 
141     kBinB, 
142     kBinC, 
143     kBinE,
144     kBinPileUp, 
145     kBinMCNSD,
146     kBinOffline,
147     kWithTrigger, 
148     kWithVertex, 
149     kAccepted
150   };
151   /** 
152    * Default constructor 
153    * 
154    * Used by ROOT I/O sub-system - do not use
155    */
156   AliAODForwardMult();
157   /** 
158    * Constructor 
159    * 
160    * @param isMC Whether this was from MC or not 
161    */
162   AliAODForwardMult(Bool_t isMC);
163   /** 
164    * Destructor 
165    */
166   virtual ~AliAODForwardMult() {} // Destructor 
167   /** 
168    * Initialize 
169    * 
170    * @param etaAxis  Pseudo-rapidity axis
171    */
172   void Init(const TAxis& etaAxis);
173   /** 
174    * Get the @f$ d^2N_{ch}/d\eta d\phi@f$ histogram, 
175    *
176    * @return @f$ d^2N_{ch}/d\eta d\phi@f$ histogram, 
177    */  
178   const TH2D& GetHistogram() const { return fHist; } // Get histogram 
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   TH2D& GetHistogram() { return fHist; } // Get histogram 
185   /** 
186    * Get the trigger mask 
187    * 
188    * @return Trigger mask 
189    */
190   UInt_t GetTriggerMask() const { return fTriggers; } // Get triggers
191   /** 
192    * Set the trigger mask 
193    * 
194    * @param trg Trigger mask
195    */
196   void SetTriggerMask(UInt_t trg) { fTriggers = trg; } // Set triggers 
197   /** 
198    * Set bit(s) in trigger mask 
199    * 
200    * @param bits bit(s) to set 
201    */
202   void SetTriggerBits(UInt_t bits) { fTriggers |= bits; } // Set trigger bits
203   /** 
204    * Check if bit(s) are set in the trigger mask 
205    * 
206    * @param bits Bits to test for 
207    * 
208    * @return 
209    */
210   Bool_t IsTriggerBits(UInt_t bits) const;
211   /** 
212    * Whether we have any trigger bits 
213    */
214   Bool_t HasTrigger() const { return fTriggers != 0; } // Check for triggers
215   /** 
216    * Clear all data 
217    * 
218    * @param option  Passed on to TH2::Reset verbatim
219    */
220   void Clear(Option_t* option="");
221   /** 
222    * browse this object 
223    * 
224    * @param b Browser 
225    */
226   void Browse(TBrowser* b);
227   /** 
228    * This is a folder 
229    * 
230    * @return Always true
231    */
232   Bool_t IsFolder() const { return kTRUE; } // Always true 
233   /** 
234    * Print content 
235    * 
236    * @param option Passed verbatim to TH2::Print 
237    */
238   void Print(Option_t* option="") const;
239   /** 
240    * Set the z coordinate of the interaction point
241    * 
242    * @param ipZ Interaction point z coordinate
243    */
244   void SetIpZ(Float_t ipZ) { fIpZ = ipZ; } // Set Ip's Z coordinate
245   /** 
246    * Set the center of mass energy per nucleon-pair.  This is stored 
247    * in the (0,0) of the histogram 
248    * 
249    * @param sNN Center of mass energy per nucleon pair (GeV)
250    */
251   void SetSNN(UShort_t sNN); 
252   /** 
253    * Get the collision system number
254    * - 0: Unknown 
255    * - 1: pp
256    * - 2: PbPb
257    * 
258    * @param sys Collision system number
259    */
260   void SetSystem(UShort_t sys);
261   /** 
262    * Set the event centrality 
263    * 
264    * @param c Centrality 
265    */
266   void SetCentrality(Float_t c) { fCentrality = c; }
267   /** 
268    * Set the z coordinate of the interaction point
269    * 
270    * @return Interaction point z coordinate
271    */
272   Float_t GetIpZ() const { return fIpZ; } // Get Ip's Z coordinate 
273   /** 
274    * Check if we have a valid z coordinate of the interaction point
275    *
276    * @return True if we have a valid interaction point z coordinate
277    */
278   Bool_t HasIpZ() const;
279   /** 
280    * Get the center of mass energy per nucleon pair (GeV)
281    * 
282    * @return Center of mass energy per nucleon pair (GeV)
283    */
284   UShort_t GetSNN() const;
285   /** 
286    * Get the collision system number
287    * - 0: Unknown 
288    * - 1: pp
289    * - 2: PbPb
290    * 
291    * @return Collision system number
292    */
293   UShort_t GetSystem() const;
294   /** 
295    * Check if the z coordinate of the interaction point is within the
296    * given limits.  Note that the convention used corresponds to the
297    * convention used in ROOTs TAxis.
298    * 
299    * @param low  Lower cut (inclusive)
300    * @param high Upper cut (exclusive)
301    * 
302    * @return true if @f$ low \ge ipz < high@f$ 
303    */
304   Bool_t InRange(Float_t low, Float_t high) const;
305   /** 
306    * Get the event centrality 
307    * 
308    * 
309    * @return 
310    */
311   Float_t GetCentrality() const { return fCentrality; }
312   /** 
313    * Check if we have a valid centrality 
314    * 
315    * 
316    * @return 
317    */
318   Bool_t  HasCentrality() const { return !(fCentrality  < 0); }
319   
320   /** 
321    * Get the name of the object 
322    * 
323    * @return Name of object 
324    */
325   const Char_t* GetName() const { return (fIsMC ? "ForwardMC" : "Forward"); }
326   /** 
327    * Check if event meets the passses requirements.   
328    *
329    * It returns true if @e all of the following is true 
330    *
331    * - The trigger is within the bit mask passed.
332    * - The vertex is within the specified limits. 
333    * - The centrality is within the specified limits, or if lower
334    *   limit is equal to or larger than the upper limit.
335    * 
336    * If a histogram is passed in the last parameter, then that
337    * histogram is filled with the trigger bits. 
338    * 
339    * @param triggerMask  Trigger mask
340    * @param vzMin        Minimum @f$ v_z@f$ (in centimeters)
341    * @param vzMax        Maximum @f$ v_z@f$ (in centimeters) 
342    * @param cMin         Minimum centrality (in percent)
343    * @param cMax         Maximum centrality (in percent)
344    * @param hist         Histogram to fill 
345    * 
346    * @return @c true if the event meets the requirements 
347    */
348   Bool_t CheckEvent(Int_t    triggerMask=kInel,
349                     Double_t vzMin=-10, Double_t vzMax=10,
350                     UShort_t cMin=0,    UShort_t cMax=100, 
351                     TH1*     hist=0) const;
352   /** 
353    * Get a string correspondig to the trigger mask
354    * 
355    * @param mask Trigger mask 
356    * 
357    * @return Static string (copy before use)
358    */
359   static const Char_t* GetTriggerString(UInt_t mask);
360   /** 
361    * Make a histogram to record triggers in. 
362    *
363    * The bins defined by the trigger enumeration in this class.  One
364    * can use this enumeration to retrieve the number of triggers for
365    * each class.
366    * 
367    * @param name Name of the histogram 
368    * 
369    * @return Newly allocated histogram 
370    */
371   static TH1I* MakeTriggerHistogram(const char* name="triggers");
372 protected: 
373   Bool_t  fIsMC;     // Whether this is from MC 
374   TH2D    fHist;     // Histogram of d^2N_{ch}/(deta dphi) for this event
375   UInt_t  fTriggers; // Trigger bit mask 
376   Float_t fIpZ;      // Z coordinate of the interaction point
377   Float_t fCentrality; // Event centrality 
378
379   static const Float_t fgkInvalidIpZ; // Invalid IpZ value 
380   ClassDef(AliAODForwardMult,2); // AOD forward multiplicity 
381 };
382
383 //____________________________________________________________________
384 inline Bool_t
385 AliAODForwardMult::InRange(Float_t low, Float_t high) const 
386 {
387   return HasIpZ() && fIpZ >= low && fIpZ < high;
388 }
389
390 //____________________________________________________________________
391 inline Bool_t 
392 AliAODForwardMult::IsTriggerBits(UInt_t bits) const 
393
394   return HasTrigger() && ((fTriggers & bits) == bits); 
395 }
396
397 #endif
398 // Local Variables:
399 //  mode: C++
400 // End:
401