Various fixes to deal with centrality
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDDensityCalculator.h
1 // This class calculates the inclusive charged particle density
2 // in each for the 5 FMD rings. 
3 //
4 #ifndef ALIFMDDENSITYCALCULATOR_H
5 #define ALIFMDDENSITYCALCULATOR_H
6 #include <TNamed.h>
7 #include <TList.h>
8 #include <TArrayI.h>
9 #include "AliForwardUtil.h"
10 class AliESDFMD;
11 class TH2D;
12 class TH1D;
13 class TProfile;
14 class AliFMDCorrELossFit;
15
16 /** 
17  * This class calculates the inclusive charged particle density
18  * in each for the 5 FMD rings. 
19  *
20  * @par Input:
21  *   - AliESDFMD object possibly corrected for sharing
22  *
23  * @par Output:
24  *   - 5 RingHistos objects - each with a number of vertex dependent 
25  *     2D histograms of the inclusive charge particle density 
26  * 
27  * @par Corrections used: 
28  *   - AliFMDAnaCalibEnergyDistribution 
29  *   - AliFMDDoubleHitCorrection 
30  *   - AliFMDDeadCorrection 
31  *
32  * @ingroup pwg2_forward_algo
33  */
34 class AliFMDDensityCalculator : public TNamed
35 {
36 public:
37   /** 
38    * Constructor 
39    */
40   AliFMDDensityCalculator();
41   /** 
42    * Constructor 
43    * 
44    * @param name Name of object
45    */
46   AliFMDDensityCalculator(const char* name);
47   /** 
48    * Copy constructor 
49    * 
50    * @param o Object to copy from 
51    */
52   AliFMDDensityCalculator(const AliFMDDensityCalculator& o);
53   /** 
54    * Destructor 
55    */
56   virtual ~AliFMDDensityCalculator();
57   /** 
58    * Assignement operator
59    * 
60    * @param o Object to assign from 
61    * 
62    * @return Reference to this object
63    */
64   AliFMDDensityCalculator& operator=(const AliFMDDensityCalculator& o);
65   /** 
66    * Initialize this sub-algorithm
67    * 
68    * @param etaAxis Not used 
69    */
70   virtual void Init(const TAxis& etaAxis);
71   /** 
72    * Do the calculations 
73    * 
74    * @param fmd      AliESDFMD object (possibly) corrected for sharing
75    * @param hists    Histogram cache
76    * @param vtxBin   Vertex bin 
77    * @param lowFlux  Low flux flag. 
78    * 
79    * @return true on successs 
80    */
81   virtual Bool_t Calculate(const AliESDFMD& fmd, 
82                            AliForwardUtil::Histos& hists, 
83                            UShort_t vtxBin, Bool_t lowFlux);
84   /** 
85    * Scale the histograms to the total number of events 
86    * 
87    * @param dir     where to put the output
88    * @param nEvents Number of events 
89    */
90   virtual void ScaleHistograms(const TList* dir, Int_t nEvents);
91   /** 
92    * Output diagnostic histograms to directory 
93    * 
94    * @param dir List to write in
95    */  
96   virtual void DefineOutput(TList* dir);
97   /** 
98    * Set the debug level.  The higher the value the more output 
99    * 
100    * @param dbg Debug level 
101    */
102   void SetDebug(Int_t dbg=1) { fDebug = dbg; }
103   /** 
104    * Maximum particle weight to use 
105    * 
106    * @param m 
107    */
108   void SetMaxParticles(UShort_t m) { fMaxParticles = m; }  
109   /** 
110    * Set whether to use poisson statistics to estimate the 
111    * number of particles that has hit within a region.  If this is true, 
112    * then the average charge particle density is given by 
113    * @f[
114    *  \lamda = -\log\left(\frac{N_e}{N_t}\right)
115    * @f]
116    * where $N_e$ is the number of strips within the region that has no
117    * hits over threshold, and $N_t$ is the total number of strips in the 
118    * region/ 
119    * 
120    * @param u Whether to use poisson statistics to estimate the 
121    * number of particles that has hit within a region.
122    */
123   void SetUsePoisson(Bool_t u) { fUsePoisson = u; }
124   /** 
125    * Set the lower multiplicity cut.  This overrides the setting in
126    * the energy loss fits.
127    * 
128    * @param cut Cut to use 
129    */
130   void SetMultCut(Double_t cut) { fMultCut = cut; }
131   /** 
132    * Get the multiplicity cut.  If the user has set fMultCut (via
133    * SetMultCut) then that value is used.  If not, then the lower
134    * value of the fit range for the enery loss fits is returned.
135    * 
136    * @return Lower cut on multiplicity
137    */
138   Double_t GetMultCut() const;
139   /** 
140    * Print information 
141    * 
142    * @param option Print options 
143    *   - max  Print max weights 
144    */
145   void Print(Option_t* option="") const;
146 protected:
147   /** 
148    * Find the max weight to use for FMD<i>dr</i> in eta bin @a iEta
149    * 
150    * @param cor   Correction
151    * @param d     Detector 
152    * @param r     Ring 
153    * @param iEta  Eta bin 
154    */
155   Int_t FindMaxWeight(const AliFMDCorrELossFit* cor,
156                       UShort_t d, Char_t r, Int_t iEta) const;
157
158   /** 
159    * Find the max weights and cache them 
160    * 
161    */  
162   void CacheMaxWeights();
163   /** 
164    * Find the (cached) maximum weight for FMD<i>dr</i> in 
165    * @f$\eta@f$ bin @a iEta
166    * 
167    * @param d     Detector
168    * @param r     Ring
169    * @param iEta  Eta bin
170    * 
171    * @return max weight or <= 0 in case of problems 
172    */
173   Int_t GetMaxWeight(UShort_t d, Char_t r, Int_t iEta) const;
174   /** 
175    * Find the (cached) maximum weight for FMD<i>dr</i> iat
176    * @f$\eta@f$ 
177    * 
178    * @param d     Detector
179    * @param r     Ring
180    * @param eta   Eta bin
181    * 
182    * @return max weight or <= 0 in case of problems 
183    */
184   Int_t GetMaxWeight(UShort_t d, Char_t r, Float_t eta) const;
185
186   /** 
187    * Get the number of particles corresponding to the signal mult
188    * 
189    * @param mult     Signal
190    * @param d        Detector
191    * @param r        Ring 
192    * @param s        Sector 
193    * @param t        Strip (not used)
194    * @param v        Vertex bin 
195    * @param eta      Pseudo-rapidity 
196    * @param lowFlux  Low-flux flag 
197    * 
198    * @return The number of particles 
199    */
200   virtual Float_t NParticles(Float_t mult, 
201                              UShort_t d, Char_t r, UShort_t s, UShort_t t, 
202                              UShort_t v, Float_t eta, Bool_t lowFlux) const;
203   /** 
204    * Get the inverse correction factor.  This consist of
205    * 
206    * - acceptance correction (corners of sensors) 
207    * - double hit correction (for low-flux events) 
208    * - dead strip correction 
209    * 
210    * @param d        Detector
211    * @param r        Ring 
212    * @param s        Sector 
213    * @param t        Strip (not used)
214    * @param v        Vertex bin 
215    * @param eta      Pseudo-rapidity 
216    * @param lowFlux  Low-flux flag 
217    * 
218    * @return 
219    */
220   virtual Float_t Correction(UShort_t d, Char_t r, UShort_t s, UShort_t t, 
221                              UShort_t v, Float_t eta, Bool_t lowFlux) const;
222   /** 
223    * Get the acceptance correction for strip @a t in an ring of type @a r
224    * 
225    * @param r  Ring type ('I' or 'O')
226    * @param t  Strip number 
227    * 
228    * @return Inverse acceptance correction 
229    */
230   virtual Float_t AcceptanceCorrection(Char_t r, UShort_t t) const;
231   /** 
232    * Generate the acceptance corrections 
233    * 
234    * @param r Ring to generate for 
235    * 
236    * @return Newly allocated histogram of acceptance corrections
237    */
238   virtual TH1D*   GenerateAcceptanceCorrection(Char_t r) const;
239   /** 
240    * Internal data structure to keep track of the histograms
241    */
242   struct RingHistos : public AliForwardUtil::RingHistos
243   { 
244     /** 
245      * Default CTOR
246      */
247     RingHistos();
248     /** 
249      * Constructor
250      * 
251      * @param d detector
252      * @param r ring 
253      */
254     RingHistos(UShort_t d, Char_t r);
255     /** 
256      * Copy constructor 
257      * 
258      * @param o Object to copy from 
259      */
260     RingHistos(const RingHistos& o);
261     /** 
262      * Assignment operator 
263      * 
264      * @param o Object to assign from 
265      * 
266      * @return Reference to this 
267      */
268     RingHistos& operator=(const RingHistos& o);
269     /** 
270      * Destructor 
271      */
272     ~RingHistos();
273     /** 
274      * Initialize the object 
275      * 
276      * @param eAxis 
277      */
278     void Init(const TAxis& eAxis);
279     /** 
280      * Make output 
281      * 
282      * @param dir Where to put it 
283      */
284     void Output(TList* dir);
285     /** 
286      * Scale the histograms to the total number of events 
287      * 
288      * @param dir     Where the output is 
289      * @param nEvents Number of events 
290      */
291     void ScaleHistograms(TList* dir, Int_t nEvents);
292     /** 
293      * Create Poisson histograms 
294      */
295     void ResetPoissonHistos();
296     TH2D*     fEvsN;         // Correlation of Eloss vs uncorrected Nch
297     TH2D*     fEvsM;         // Correlation of Eloss vs corrected Nch
298     TProfile* fEtaVsN;       // Average uncorrected Nch vs eta
299     TProfile* fEtaVsM;       // Average corrected Nch vs eta
300     TProfile* fCorr;         // Average correction vs eta
301     TH2D*     fDensity;      // Distribution inclusive Nch
302     TH2D*     fELossVsPoisson; // Correlation of energy loss vs Poisson N_ch
303     TH2D*     fTotalStrips;  // Total number of strips in a region
304     TH2D*     fEmptyStrips;  // Total number of strips in a region
305     TH2D*     fBasicHits  ;  // Total number basic hits in a region
306     TH2D*     fEmptyVsTotal; // # of empty strips vs total number of strips 
307     
308     
309     ClassDef(RingHistos,3);
310   };
311   /** 
312    * Get the ring histogram container 
313    * 
314    * @param d Detector
315    * @param r Ring 
316    * 
317    * @return Ring histogram container 
318    */
319   RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
320
321   TList    fRingHistos;    //  List of histogram containers
322   Double_t fMultCut;       //  Low cut on scaled energy loss
323   TH1D*    fSumOfWeights;  //  Histogram
324   TH1D*    fWeightedSum;   //  Histogram
325   TH1D*    fCorrections;   //  Histogram
326   UShort_t fMaxParticles;  //  Maximum particle weight to use 
327   Bool_t   fUsePoisson;    //  If true, then use poisson statistics 
328   TH1D*    fAccI;          //  Acceptance correction for inner rings
329   TH1D*    fAccO;          //  Acceptance correction for outer rings
330   TArrayI  fFMD1iMax;      //  Array of max weights 
331   TArrayI  fFMD2iMax;      //  Array of max weights 
332   TArrayI  fFMD2oMax;      //  Array of max weights 
333   TArrayI  fFMD3iMax;      //  Array of max weights 
334   TArrayI  fFMD3oMax;      //  Array of max weights 
335   Int_t    fDebug;         //  Debug level 
336
337   ClassDef(AliFMDDensityCalculator,2); // Calculate Nch density 
338 };
339
340 #endif
341 // Local Variables:
342 //   mode: C++
343 // End:
344