Doc fixes
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDEnergyFitter.h
1 #ifndef ALIFMDENERGYFITTER_H
2 #define ALIFMDENERGYFITTER_H
3 #include <TNamed.h>
4 #include <TH1D.h>
5 #include <TAxis.h>
6 #include <TList.h>
7 #include <TObjArray.h>
8 #include <TClonesArray.h>
9 #include "AliFMDCorrELossFit.h"
10 #include "AliForwardUtil.h"
11 class AliESDFMD;
12 class TFitResult;
13 class TF1;
14 class TArrayD;
15
16 /**
17  * Class to fit the energy distribution.  
18  *
19  * @par Input: 
20  *    - AliESDFMD object  - from reconstruction
21  *
22  * @par Output: 
23  *    - Lists of histogram - one per ring.  Each list has a number of 
24  *      histograms corresponding to the number of eta bins defined.  
25  *
26  * @par Corrections used: 
27  *    - None
28  *
29  *
30  * @ingroup pwg2_forward_algo
31  */
32 class AliFMDEnergyFitter : public TNamed
33 {
34 public: 
35     enum { 
36       kC        = AliForwardUtil::ELossFitter::kC,
37       kDelta    = AliForwardUtil::ELossFitter::kDelta, 
38       kXi       = AliForwardUtil::ELossFitter::kXi, 
39       kSigma    = AliForwardUtil::ELossFitter::kSigma, 
40       kSigmaN   = AliForwardUtil::ELossFitter::kSigmaN,
41       kN        = AliForwardUtil::ELossFitter::kN, 
42       kA        = AliForwardUtil::ELossFitter::kA
43     };
44
45   /** 
46    * Destructor
47    */
48   virtual ~AliFMDEnergyFitter();
49   /** 
50    * Default Constructor - do not use 
51    */
52   AliFMDEnergyFitter();
53   /** 
54    * Constructor 
55    * 
56    * @param title Title of object  - not significant 
57    */
58   AliFMDEnergyFitter(const char* title);
59   /** 
60    * Copy constructor 
61    * 
62    * @param o Object to copy from 
63    */
64   AliFMDEnergyFitter(const AliFMDEnergyFitter& o);
65   /** 
66    * Assignment operator 
67    * 
68    * @param o Object to assign from 
69    * 
70    * @return Reference to this 
71    */
72   AliFMDEnergyFitter& operator=(const AliFMDEnergyFitter& o);
73
74   /** 
75    * Initialise the task
76    * 
77    * @param etaAxis The eta axis to use.  Note, that if the eta axis
78    * has already been set (using SetEtaAxis), then this parameter will be 
79    * ignored
80    */
81   void Init(const TAxis& etaAxis);
82   /** 
83    * Set the eta axis to use.  This will force the code to use this
84    * eta axis definition - irrespective of whatever axis is passed to
85    * the Init member function.  Therefore, this member function can be
86    * used to force another eta axis than one found in the correction
87    * objects. 
88    * 
89    * @param nBins  Number of bins 
90    * @param etaMin Minimum of the eta axis 
91    * @param etaMax Maximum of the eta axis 
92    */
93   void SetEtaAxis(Int_t nBins, Double_t etaMin, Double_t etaMax);
94   /** 
95    * Set the eta axis to use.  This will force the code to use this
96    * eta axis definition - irrespective of whatever axis is passed to
97    * the Init member function.  Therefore, this member function can be
98    * used to force another eta axis than one found in the correction
99    * objects. 
100    * 
101    * @param etaAxis Eta axis to use 
102    */
103   void SetEtaAxis(const TAxis& etaAxis);
104   /** 
105    * Set the low cut used for energy 
106    * 
107    * @param lowCut Low cut
108    */
109   void SetLowCut(Double_t lowCut=0.3) { fLowCut = lowCut; }
110   /** 
111    * Set the number of bins to subtract 
112    * 
113    * @param n 
114    */
115   void SetFitRangeBinWidth(UShort_t n=4) { fFitRangeBinWidth = n; }
116   /** 
117    * Whether or not to enable fitting of the final merged result.  
118    * Note, fitting takes quite a while and one should be careful not to do 
119    * this needlessly 
120    * 
121    * @param doFit Whether to do the fits or not 
122    */
123   void SetDoFits(Bool_t doFit=kTRUE) { fDoFits = doFit; }
124   /** 
125    * Set whether to make the corrections object on the output.  Note,
126    * fits should be enable for this to have any effect.
127    * 
128    * @param doMake If true (false is default), do make the corrections object. 
129    */
130   void SetDoMakeObject(Bool_t doMake=kTRUE) { fDoMakeObject = doMake; }
131   /** 
132    * Set how many particles we will try to fit at most to the data
133    * 
134    * @param n Max number of particle to try to fit 
135    */
136   void SetNParticles(UShort_t n) { fNParticles = (n < 1 ? 1 : (n > 5 ? 5 : n)); }
137   /** 
138    * Set the minimum number of entries each histogram must have 
139    * before we try to fit our response function to it
140    * 
141    * @param n Minimum number of entries
142    */
143   void SetMinEntries(UShort_t n) { fMinEntries = (n < 1 ? 1 : n); }
144   /**
145    * Set maximum energy loss to consider 
146    *
147    * @param x Maximum energy loss to consider 
148    */
149   void SetMaxE(Double_t x) { fMaxE = x; }
150   /**
151    * Set number of energy loss bins 
152    *
153    * @param x Number of energy loss bins 
154    */
155   void SetNEbins(Int_t x) { fNEbins = x; }
156   /** 
157    * Set the maximum relative error 
158    * 
159    * @param e Maximum relative error 
160    */
161   void SetMaxRelativeParameterError(Double_t e=0.2) { fMaxRelParError = e; }
162   /** 
163    * Set the maximum @f$ \chi^2/\nu@f$ 
164    * 
165    * @param c Maximum @f$ \chi^2/\nu@f$ 
166    */
167   void SetMaxChi2PerNDF(Double_t c=10) { fMaxChi2PerNDF = c; }
168   /** 
169    * Set the least weight
170    * 
171    * @param c Least weight
172    */
173   void SetMinWeight(Double_t c=1e-7) { fMinWeight = c; }
174   /**
175    * Set wheter to use increasing bin sizes 
176    *
177    * @param x Wheter to use increasing bin sizes 
178    */
179   void SetUseIncreasingBins(Bool_t x) { fUseIncreasingBins = x; }
180   /** 
181    * Fitter the input AliESDFMD object
182    * 
183    * @param input     Input 
184    * @param empty     Whether the event is 'empty'
185    * 
186    * @return True on success, false otherwise 
187    */
188   Bool_t Accumulate(const AliESDFMD& input, 
189                     Bool_t           empty);
190   /** 
191    * Scale the histograms to the total number of events 
192    * 
193    * @param dir Where the histograms are  
194    */
195   void Fit(const TList* dir);
196   /** 
197    * Generate the corrections object 
198    * 
199    * @param dir List to analyse 
200    */
201   void MakeCorrectionsObject(TList* dir);
202   
203   /** 
204    * Define the output histograms.  These are put in a sub list of the
205    * passed list.   The histograms are merged before the parent task calls 
206    * AliAnalysisTaskSE::Terminate 
207    * 
208    * @param dir Directory to add to 
209    */
210   void DefineOutput(TList* dir);
211   /** 
212    * Set the debug level.  The higher the value the more output 
213    * 
214    * @param dbg Debug level 
215    */
216   void SetDebug(Int_t dbg=1);
217   /** 
218    * Print information
219    * 
220    * @param option Not used 
221    */
222   void Print(Option_t* option="") const;
223 protected:
224   /** 
225    * Internal data structure to keep track of the histograms
226    */
227   struct RingHistos : public AliForwardUtil::RingHistos
228   { 
229     /** 
230      * Default CTOR
231      */
232     RingHistos();
233     /** 
234      * Constructor
235      * 
236      * @param d detector
237      * @param r ring 
238      */
239     RingHistos(UShort_t d, Char_t r);
240     /** 
241      * Copy constructor 
242      * 
243      * @param o Object to copy from 
244      */
245     RingHistos(const RingHistos& o);
246     /** 
247      * Assignment operator 
248      * 
249      * @param o Object to assign from 
250      * 
251      * @return Reference to this 
252      */
253     RingHistos& operator=(const RingHistos& o);
254     /** 
255      * Destructor 
256      */
257     ~RingHistos();
258     /** 
259      * Define outputs
260      * 
261      * @param dir 
262      */
263     void Output(TList* dir);
264     /** 
265      * Initialise object 
266      * 
267      * @param eAxis      Eta axis
268      * @param maxDE      Max energy loss to consider 
269      * @param nDEbins    Number of bins 
270      * @param useIncrBin Whether to use an increasing bin size 
271      */
272     void Init(const TAxis& eAxis, 
273               Double_t     maxDE=10, 
274               Int_t        nDEbins=300, 
275               Bool_t       useIncrBin=true);
276     /** 
277      * Fill histogram 
278      * 
279      * @param empty  True if event is empty
280      * @param ieta   Eta bin
281      * @param mult   Signal 
282      */
283     void Fill(Bool_t empty, Int_t ieta, Double_t mult);
284     /** 
285      * Fit each histogram to up to @a nParticles particle responses.
286      * 
287      * @param dir         Output list 
288      * @param eta         Eta axis 
289      * @param lowCut      Lower cut 
290      * @param nParticles  Max number of convolved landaus to fit
291      * @param minEntries  Minimum number of entries 
292      * @param minusBins   Number of bins from peak to subtract to 
293      *                    get the fit range 
294      * @param relErrorCut Cut applied to relative error of parameter. 
295      *                    Note, for multi-particle weights, the cut 
296      *                    is loosend by a factor of 2 
297      * @param chi2nuCut   Cut on @f$ \chi^2/\nu@f$ - 
298      *                    the reduced @f$\chi^2@f$ 
299      */
300     TObjArray* Fit(TList* dir, 
301                    const TAxis& eta,
302                    Double_t     lowCut, 
303                    UShort_t     nParticles,
304                    UShort_t     minEntries,
305                    UShort_t     minusBins,
306                    Double_t     relErrorCut, 
307                    Double_t     chi2nuCut) const;
308     /** 
309      * Fit a signal histogram.  First, the bin @f$ b_{min}@f$ with
310      * maximum bin content in the range @f$ [E_{min},\infty]@f$ is
311      * found.  Then the fit range is set to the bin range 
312      * @f$ [b_{min}-\Delta b,b_{min}+2\Delta b]@f$, and a 1 
313      * particle signal is fitted to that.  The parameters of that fit 
314      * is then used as seeds for a fit of the @f$ N@f$ particle response 
315      * to the data in the range 
316      * @f$ [b_{min}-\Delta b,N(\Delta_1+\xi_1\log(N))+2N\xi@f$
317      * 
318      * @param dist        Histogram to fit 
319      * @param lowCut      Lower cut @f$ E_{min}@f$ on signal 
320      * @param nParticles  Max number @f$ N@f$ of convolved landaus to fit
321      * @param minusBins   Number of bins @f$ \Delta b@f$ from peak to 
322      *                    subtract to get the fit range 
323      * @param relErrorCut Cut applied to relative error of parameter. 
324      *                    Note, for multi-particle weights, the cut 
325      *                    is loosend by a factor of 2 
326      * @param chi2nuCut   Cut on @f$ \chi^2/\nu@f$ - 
327      *                    the reduced @f$\chi^2@f$ 
328      * 
329      * @return The best fit function 
330      */
331     TF1* FitHist(TH1*     dist,
332                  Double_t lowCut, 
333                  UShort_t nParticles,
334                  UShort_t minusBins,
335                  Double_t relErrorCut, 
336                  Double_t chi2nuCut) const;
337     /** 
338      * Find the best fits 
339      * 
340      * @param d              Parent list
341      * @param obj            Object to add fits to
342      * @param eta            Eta axis 
343      * @param relErrorCut    Cut applied to relative error of parameter. 
344      *                       Note, for multi-particle weights, the cut 
345      *                       is loosend by a factor of 2 
346      * @param chi2nuCut      Cut on @f$ \chi^2/\nu@f$ - 
347      *                       the reduced @f$\chi^2@f$ 
348      * @param minWeightCut   Least valid @f$ a_i@f$ 
349      */
350     void FindBestFits(TList*              d, 
351                       AliFMDCorrELossFit& obj,
352                       const TAxis&        eta,     
353                       Double_t            relErrorCut, 
354                       Double_t            chi2nuCut,
355                       Double_t            minWeightCut);
356     /** 
357      * Find the best fit 
358      * 
359      * @param dist           Histogram 
360      * @param relErrorCut    Cut applied to relative error of parameter. 
361      *                       Note, for multi-particle weights, the cut 
362      *                       is loosend by a factor of 2 
363      * @param chi2nuCut      Cut on @f$ \chi^2/\nu@f$ - 
364      *                       the reduced @f$\chi^2@f$ 
365      * @param minWeightCut   Least valid @f$ a_i@f$ 
366      * 
367      * @return Best fit 
368      */
369     AliFMDCorrELossFit::ELossFit* FindBestFit(TH1* dist,
370                                               Double_t relErrorCut, 
371                                               Double_t chi2nuCut,
372                                               Double_t minWeightCut);
373     /** 
374      * Check the result of the fit. Returns true if 
375      * - @f$ \chi^2/\nu < \max{\chi^2/\nu}@f$
376      * - @f$ \Delta p_i/p_i < \delta_e@f$ for all parameters.  Note, 
377      *   for multi-particle fits, this requirement is relaxed by a 
378      *   factor of 2
379      * - @f$ a_{n} > 10^{-7}@f$ when fitting to an @f$ n@f$ 
380      *   particle response 
381      * 
382      * @param r           Result to check
383      * @param relErrorCut Cut @f$ \delta_e@f$ applied to relative error 
384      *                    of parameter.  
385      * @param chi2nuCut   Cut @f$ \max{\chi^2/\nu}@f$ 
386      * 
387      * @return true if fit is good. 
388      */
389     Bool_t CheckResult(TFitResult* r,
390                        Double_t    relErrorCut, 
391                        Double_t    chi2nuCut) const;
392     /** 
393      * Make an axis with increasing bins 
394      * 
395      * @param n    Number of bins 
396      * @param min  Minimum 
397      * @param max  Maximum
398      * 
399      * @return An axis with quadratically increasing bin size 
400      */
401     TArrayD MakeIncreasingAxis(Int_t n, Double_t min, Double_t max) const;
402     /** 
403      * Make E/E_mip histogram 
404      * 
405      * @param ieta    Eta bin
406      * @param eMin    Least signal
407      * @param eMax    Largest signal 
408      * @param deMax   Maximum energy loss 
409      * @param nDeBins Number energy loss bins 
410      * @param incr    Whether to make bins of increasing size
411      */
412     void Make(Int_t ieta, Double_t eMin, Double_t eMax, 
413               Double_t deMax=12, Int_t nDeBins=300, Bool_t incr=true);
414     /** 
415      * Make a parameter histogram
416      * 
417      * @param name   Name of histogram.
418      * @param title  Title of histogram. 
419      * @param eta    Eta axis 
420      * 
421      * @return 
422      */
423     TH1D* MakePar(const char* name, const char* title, const TAxis& eta) const;
424     /** 
425      * Make a histogram that contains the results of the fit over the full ring 
426      * 
427      * @param name  Name 
428      * @param title Title
429      * @param eta   Eta axis 
430      * @param low   Least bin
431      * @param high  Largest bin
432      * @param val   Value of parameter 
433      * @param err   Error on parameter 
434      * 
435      * @return The newly allocated histogram 
436      */
437     TH1D* MakeTotal(const char* name, 
438                     const char* title, 
439                     const TAxis& eta, 
440                     Int_t low, 
441                     Int_t high, 
442                     Double_t val, 
443                     Double_t err) const;
444     TH1D*        fEDist;        // Ring energy distribution 
445     TH1D*        fEmpty;        // Ring energy distribution for empty events
446     TList        fEtaEDists;    // Energy distributions per eta bin. 
447     TList*       fList;
448     TClonesArray fFits;
449     Int_t        fDebug;
450     ClassDef(RingHistos,1);
451   };
452   /** 
453    * Get the ring histogram container 
454    * 
455    * @param d Detector
456    * @param r Ring 
457    * 
458    * @return Ring histogram container 
459    */
460   RingHistos* GetRingHistos(UShort_t d, Char_t r) const;
461
462   TList    fRingHistos;    // List of histogram containers
463   Double_t fLowCut;        // Low cut on energy
464   UShort_t fNParticles;    // Number of landaus to try to fit 
465   UShort_t fMinEntries;    // Minimum number of entries
466   UShort_t fFitRangeBinWidth;// Number of bins to subtract from found max
467   Bool_t   fDoFits;        // Whether to actually do the fits 
468   Bool_t   fDoMakeObject;  // Whether to make corrections object
469   TAxis    fEtaAxis;       // Eta axis 
470   Double_t fMaxE;          // Maximum energy loss to consider 
471   Int_t    fNEbins;        // Number of energy loss bins 
472   Bool_t   fUseIncreasingBins; // Wheter to use increasing bin sizes 
473   Double_t fMaxRelParError;// Relative error cut
474   Double_t fMaxChi2PerNDF; // chi^2/nu cit
475   Double_t fMinWeight;     // Minimum weight value 
476   Int_t    fDebug;         // Debug level 
477   
478
479   ClassDef(AliFMDEnergyFitter,1); //
480 };
481
482 #endif
483 // Local Variables:
484 //  mode: C++ 
485 // End: