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