More code clean up.
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliFMDCorrELossFit.h
1 #ifndef ALIROOT_PWG2_FORWARD_ALIFMDCORRELOSSFIT_H
2 #define ALIROOT_PWG2_FORWARD_ALIFMDCORRELOSSFIT_H
3 #include <TObject.h>
4 #include <TAxis.h>
5 #include <TObjArray.h>
6 class TF1;
7 class TBrowser;
8
9 /** 
10  * Object holding the Energy loss fit 'correction'
11  * 
12  * These are generated from Monte-Carlo or real ESDs. 
13  *
14  * @ingroup pwg2_forward_corr
15  */
16 class AliFMDCorrELossFit : public TObject 
17 {
18 public:
19   /** 
20    * POD structure to hold data from fits 
21    * 
22    */
23   struct ELossFit : public TObject 
24   {
25     Int_t     fN;      // Number of peaks fitted
26     UShort_t  fNu;     // Number degrees of freedom
27     Double_t  fChi2;   // Chi square from fit
28     Double_t  fC;      // Scaling constant 
29     Double_t  fDelta;  // Most probable value 
30     Double_t  fXi;     // Width parameter of Landau 
31     Double_t  fSigma;  // Sigma on folded gaussian 
32     Double_t  fSigmaN; // Sigma of detector noise 
33     Double_t* fA;      // [fN] Weights 
34     Double_t  fEC;     // Error on C 
35     Double_t  fEDelta; // Error on Delta 
36     Double_t  fEXi;    // Error on Xi
37     Double_t  fESigma; // Error on sigma 
38     Double_t  fESigmaN;// Error on sigma (noise)
39     Double_t* fEA;     // [fN] Error on weights
40     Int_t     fQuality;// Assigned quality 
41     UShort_t  fDet;    // Detector 
42     Char_t    fRing;   // Ring
43     UShort_t  fBin;    // Eta bin
44
45     static Double_t fgMaxRelError;  // Global default max relative error
46     static Double_t fgLeastWeight;  // Global default least weight 
47     static Double_t fgMaxChi2nu;    // Global default maximum reduced chi^2
48     /**
49      * Default constructor 
50      * 
51      */
52     ELossFit();
53     /** 
54      * Construct from a function
55      * 
56      * @param quality Quality flag
57      * @param f       Function
58      */
59     ELossFit(Int_t quality,const TF1& f);
60     /** 
61      * Constructor with full parameter set
62      * 
63      * @param d         Detector number
64      * @param r         Ring identifier 
65      * @param quality   Quality flag
66      * @param n         @f$ N@f$ - Number of fitted peaks
67      * @param chi2      @f$ \chi^2 @f$
68      * @param nu        @f$ \nu @f$ - number degrees of freedom
69      * @param c         @f$ C2f$ - scale constant
70      * @param ec        @f$ \delta C@f$ - error on @f$ C@f$ 
71      * @param delta     @f$ \Delta2f$ - scale constant            
72      * @param edelta    @f$ \delta\Delta@f$ - error on @f$\Delta@f$ 
73      * @param xi        @f$ \xi2f$ - scale constant               
74      * @param exi       @f$ \delta\xi@f$ - error on @f$\xi@f$ 
75      * @param sigma     @f$ \sigma@f$ - scale constant             
76      * @param esigma    @f$ \delta\sigma@f$ - error on @f$\sigma@f$ 
77      * @param sigman    @f$ \sigma_n@f$ - scale constant                  
78      * @param esigman   @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$ 
79      * @param a         Array of @f$ N-1@f$ weights @f$ a_i@f$ for 
80      *                  @f$ i=2,\ldots@f$ 
81      */
82     ELossFit(Int_t     quality,UShort_t  n, 
83              Double_t  chi2,   UShort_t  nu, 
84              Double_t  c,      Double_t  ec, 
85              Double_t  delta,  Double_t  edelta, 
86              Double_t  xi,     Double_t  exi,
87              Double_t  sigma,  Double_t  esigma, 
88              Double_t  sigman, Double_t  esigman, 
89              Double_t* a,      Double_t* ea);
90     /** 
91      * Copy constructor 
92      * 
93      * @param o Object to copy from 
94      */
95     ELossFit(const ELossFit& o);
96     /** 
97      * Assignment operator 
98      * 
99      * @param o Object to assign from 
100      * 
101      * @return Reference to this object 
102      */
103     ELossFit& operator=(const ELossFit& o);
104     /** 
105      * Destructor 
106      */
107     ~ELossFit();
108     /**
109      * @{
110      * @name Access to parameters 
111      */
112     /**
113      * @return Number of peaks fitted
114      */
115     Int_t GetN() const { return fN; }
116     /**
117      * @return Number degrees of freedom
118      */
119     UShort_t GetNu() const { return fNu; }
120     /**
121      * @return Chi square from fit
122      */
123     Double_t GetChi2() const { return fChi2; }
124     /**
125      * @return Scaling constant 
126      */
127     Double_t GetC() const { return fC; }
128     /**
129      * @return Most probable value 
130      */
131     Double_t GetDelta() const { return fDelta; }
132     /**
133      * @return Width parameter of Landau 
134      */
135     Double_t GetXi() const { return fXi; }
136     /**
137      * @return Sigma on folded gaussian 
138      */
139     Double_t GetSigma() const { return fSigma; }
140     /**
141      * @return Sigma of detector noise 
142      */
143     Double_t GetSigmaN() const { return fSigmaN; }
144     /**
145      * @return Weights 
146      */
147     Double_t* GetAs() const { return fA; }
148     /**
149      * @return Weights 
150      */
151     Double_t GetA(UShort_t i) const;    
152     /**
153      * @return Error on C 
154      */
155     Double_t GetEC() const { return fEC; }
156     /**
157      * @return Error on Delta 
158      */
159     Double_t GetEDelta() const { return fEDelta; }
160     /**
161      * @return Error on Xi
162      */
163     Double_t GetEXi() const { return fEXi; }
164     /**
165      * @return Error on sigma 
166      */
167     Double_t GetESigma() const { return fESigma; }
168     /**
169      * @return Error on sigma (noise)
170      */
171     Double_t GetESigmaN() const { return fESigmaN; }
172     /**
173      * @return Error on weights
174      */
175     Double_t* GetEAs() const { return fEA; }
176     /**
177      * @return Error on weights
178      */
179     Double_t GetEA(UShort_t i) const;
180     /**
181      * @return Assigned quality 
182      */
183     Int_t GetQuality() const { return fQuality; }
184     /**
185      * @return Detector 
186      */
187     UShort_t GetDet() const { return fDet; }
188     /**
189      * @return Ring
190      */
191     Char_t GetRing() const { return fRing; }
192     /**
193      * @return Eta bin
194      */
195     UShort_t GetBin() const { return fBin; }
196     /* @} */
197
198     /** 
199      * @{ 
200      * @name Evaluation 
201      */
202     /** 
203      * Evaluate 
204      * @f[ 
205      *  f_N(x;\Delta,\xi,\sigma') = 
206      *     \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')
207      * @f] 
208      *
209      * (see AliForwardUtil::NLandauGaus) for the maximum @f$ N @f$
210      * that fulfills the requirements 
211      *
212      * @param x           Where to evaluate 
213      * @param maxN        @f$ \max{N}@f$    
214      * 
215      * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$ 
216      */
217     Double_t Evaluate(Double_t x, 
218                       UShort_t maxN=999) const;
219     /** 
220      * Evaluate 
221      * @f[ 
222      *   f_W(x;\Delta,\xi,\sigma') = 
223      *   \frac{\sum_{i=1}^{n} i a_i f_i(x;\Delta,\xi,\sigma')}{
224      *     f_N(x;\Delta,\xi,\sigma')} = 
225      *   \frac{\sum_{i=1}^{n} i a_i f(x;\Delta_i,\xi_i,\sigma_i')}{
226      *     \sum_{i=1}^{n} a_i f(x;\Delta_i,\xi_i,\sigma_i')}
227      * @f] 
228      * where @f$ n@f$ fulfills the requirements (see FindMaxWeight). 
229      *
230      * If the denominator is zero, then 1 is returned. 
231      *
232      * See also AliForwardUtil::ILandauGaus and AliForwardUtil::NLandauGaus
233      * for more information on the evaluated functions. 
234      * 
235      * @param x           Where to evaluate 
236      * @param maxN        @f$ \max{N}@f$      
237      * 
238      * @return @f$ f_W(x;\Delta,\xi,\sigma')@f$.  
239      */
240     Double_t EvaluateWeighted(Double_t x, 
241                               UShort_t maxN=9999) const;
242     /** 
243      * Find the maximum weight to use.  The maximum weight is the
244      * largest i for which 
245      * 
246      * - @f$ i \leq \max{N}@f$ 
247      * - @f$ a_i > \min{a}@f$ 
248      * - @f$ \delta a_i/a_i > \delta_{max}@f$ 
249      * 
250      * @param maxRelError @f$ \min{a}@f$ 
251      * @param leastWeight @f$ \delta_{max}@f$ 
252      * @param maxN        @f$ \max{N}@f$      
253      * 
254      * @return The largest index @f$ i@f$ for which the above
255      * conditions hold.  Will never return less than 1. 
256      */
257     Int_t FindMaxWeight(Double_t maxRelError=2*fgMaxRelError, 
258                         Double_t leastWeight=fgLeastWeight, 
259                         UShort_t maxN=999) const;
260     /* @} */
261     /** 
262      * @{
263      * @name TObject Sortable interface 
264      */
265     /** 
266      * Declare this object as sortable 
267      * 
268      * @return Always true 
269      */
270     Bool_t IsSortable() const { return kTRUE; }
271     /** 
272      * Compare to another ELossFit object. 
273      * 
274      * - +1, if this quality is better (larger) than other objects quality
275      * - -1, if this quality is worse (smaller) than other objects quality
276      * - +1, if this @f$|\chi^2/\nu-1|@f$ is smaller than the same for other
277      * - -1, if this @f$|\chi^2/\nu-1|@f$ is larger than the same for other
278      * - 0 otherwise 
279      * 
280      * @param o Other object to compare to 
281      */
282     Int_t Compare(const TObject* o) const;
283     /* @} */
284     /** 
285      * @{ 
286      * name Auxiliary member functions  
287      */
288     /** 
289      * Information to standard output 
290      * 
291      * @param option Not used 
292      */
293     void Print(Option_t* option) const; // *MENU*
294     /** 
295      * Draw this fit 
296      * 
297      * @param option Options 
298      *  - COMP  Draw components too 
299      */
300     void Draw(Option_t* option="comp"); // *MENU*
301     /** 
302      * Browse this object 
303      * 
304      * @param b Browser
305      */
306     void Browse(TBrowser* b);
307     /** 
308      * Get the name of this object 
309      * 
310      * 
311      * @return 
312      */
313     const Char_t* GetName() const;
314     /** 
315      * Calculate the quality 
316      */
317     void CalculateQuality(Double_t maxChi2nu=fgMaxChi2nu, 
318                           Double_t maxRelError=fgMaxRelError, 
319                           Double_t leastWeight=fgLeastWeight);
320     /* @} */
321     ClassDef(ELossFit,1); // Result of fit 
322   };
323
324   /** 
325    * Default constructor 
326    */
327   AliFMDCorrELossFit();
328   /** 
329    * Copy constructor 
330    * 
331    * @param o Object to copy from 
332    */
333   AliFMDCorrELossFit(const AliFMDCorrELossFit& o);
334   /** 
335    * Destructor 
336    */
337   virtual ~AliFMDCorrELossFit(); 
338   /** 
339    * Assignment operator 
340    * 
341    * @param o Object to assign from 
342    * 
343    * @return Reference to this object 
344    */
345   AliFMDCorrELossFit& operator=(const AliFMDCorrELossFit& o);
346
347   /** 
348    * @{ 
349    * @name Set fits 
350    */
351   /** 
352    * Set the fit parameters from a function 
353    * 
354    * @param d        Detector
355    * @param r        Ring 
356    * @param eta      Eta 
357    * @param quality  Quality flag
358    * @param f        Function from fit 
359    */  
360   Bool_t SetFit(UShort_t d, Char_t r, Double_t eta, Int_t quality, 
361                 const TF1& f);
362   /** 
363    * Set the fit parameters from a function 
364    * 
365    * @param d    Detector
366    * @param r    Ring 
367    * @param eta  Eta 
368    * @param f    ELoss fit result - note, the object will take ownership
369    */  
370   Bool_t SetFit(UShort_t d, Char_t r, Double_t eta, ELossFit* f);
371   /** 
372    * Set the fit parameters from a function 
373    * 
374    * @param d    Detector
375    * @param r    Ring 
376    * @param eta  Eta (bin number, 1->nBins)
377    * @param f    ELoss fit result - note, the object will take ownership
378    */  
379   Bool_t SetFit(UShort_t d, Char_t r, Int_t etaBin, ELossFit* f);
380   /** 
381    * Set the fit parameters from a function 
382    * 
383    * @param d         Detector number
384    * @param r         Ring identifier 
385    * @param quality   Quality flag
386    * @param n         @f$ N@f$ - Number of fitted peaks
387    * @param chi2      @f$ \chi^2 @f$
388    * @param nu        @f$ \nu @f$ - number degrees of freedom
389    * @param c         @f$ C2f$ - scale constant
390    * @param ec        @f$ \delta C@f$ - error on @f$ C@f$ 
391    * @param delta     @f$ \Delta2f$ - scale constant              
392    * @param edelta    @f$ \delta\Delta@f$ - error on @f$\Delta@f$ 
393    * @param xi        @f$ \xi2f$ - scale constant                 
394    * @param exi       @f$ \delta\xi@f$ - error on @f$\xi@f$ 
395    * @param sigma     @f$ \sigma@f$ - scale constant               
396    * @param esigma    @f$ \delta\sigma@f$ - error on @f$\sigma@f$ 
397    * @param sigman    @f$ \sigma_n@f$ - scale constant            
398    * @param esigman   @f$ \delta\sigma_n@f$ - error on @f$\sigma_n@f$ 
399    * @param a         Array of @f$ N-1@f$ weights @f$ a_i@f$ for 
400    *                  @f$ i=2,\ldots@f$ 
401    */
402   Bool_t SetFit(UShort_t  d,      Char_t    r, Double_t eta, 
403                 Int_t     quality,UShort_t  n, 
404                 Double_t  chi2,   UShort_t  nu, 
405                 Double_t  c,      Double_t  ec, 
406                 Double_t  delta,  Double_t  edelta, 
407                 Double_t  xi,     Double_t  exi,
408                 Double_t  sigma,  Double_t  esigma, 
409                 Double_t  sigman, Double_t  esigman, 
410                 Double_t* a,      Double_t* ea);
411   /* @} */
412   
413   /** 
414    * @{
415    * @name Set and get eta axis
416    */
417   /** 
418    * Set the eta axis to use 
419    * 
420    * @param axis Eta axis 
421    */
422   void SetEtaAxis(const TAxis& axis);
423   /** 
424    * Set the eta axis to use 
425    * 
426    * @param axis Eta axis 
427    */
428   void SetEtaAxis(Int_t nBins, Double_t min, Double_t max);
429   /** 
430    * Get the eta axis used
431    * 
432    * @return 
433    */
434   const TAxis& GetEtaAxis() const { return fEtaAxis; }
435   /** 
436    * Set the low cut used when fitting 
437    * 
438    * @param cut Cut value 
439    */
440   void SetLowCut(Double_t cut) { fLowCut = cut; }
441   /** 
442    * Get the low cut used when fitting 
443    * 
444    * @return Low cut used for fitting 
445    */
446   Double_t GetLowCut() const { return fLowCut; }
447   /** 
448    * Find the eta bin corresponding to the given eta 
449    * 
450    * @param eta  Eta value 
451    * 
452    * @return Bin (in @f$[1,N_{bins}]@f$) corresponding to the given
453    * eta, or 0 if out of range.
454    */
455   Int_t FindEtaBin(Double_t eta) const;
456   /* @} */
457
458   /**
459    * @{                                         
460    * @name Find fits 
461    */
462   /** 
463    * Find the fit corresponding to the specified parameters 
464    * 
465    * @param d   Detector 
466    * @param r   Ring 
467    * @param eta Eta value 
468    * 
469    * @return Fit parameters or null in case of problems 
470    */
471   ELossFit* FindFit(UShort_t d, Char_t r, Double_t eta) const;
472   /** 
473    * Find the fit corresponding to the specified parameters 
474    * 
475    * @param d      Detector 
476    * @param r      Ring 
477    * @param etabin Eta bin (1 based)
478    * 
479    * @return Fit parameters or null in case of problems 
480    */
481   ELossFit* FindFit(UShort_t d, Char_t r, Int_t etabin) const;
482   /* @} */
483
484   /**                                           
485    * @{ 
486    * @name Miscellaneous
487    */
488   /** 
489    * Get the ring array corresponding to the specified ring
490    * 
491    * @param d Detector 
492    * @param r Ring 
493    * 
494    * @return Pointer to ring array, or null in case of problems
495    */
496   TObjArray* GetRingArray(UShort_t d, Char_t r) const;
497   /** 
498    * Signal that this is a folder
499    * 
500    * @return Always true 
501    */
502   Bool_t IsFolder() const { return true; }
503   /** 
504    * Browse this object 
505    * 
506    * @param b 
507    */
508   void Browse(TBrowser* b);
509   /** 
510    * Draw this object 
511    *
512    * @param option Options.  Possible values are 
513    *  - err Plot error bars 
514    */
515   void Draw(Option_t* option=""); //*MENU*
516   /** 
517    * Print this object.  
518    * 
519    * @param option Options 
520    *   - R   Print recursive  
521    *
522    */
523   void Print(Option_t* option="R") const; //*MENU*
524   /* @} */
525 protected:
526   /** 
527    * Get the ring array corresponding to the specified ring
528    * 
529    * @param d Detector 
530    * @param r Ring 
531    * 
532    * @return Pointer to ring array, or newly created container 
533    */
534   TObjArray* GetOrMakeRingArray(UShort_t d, Char_t r);
535
536   TObjArray  fRings;    // Array of rings
537   TAxis      fEtaAxis;  // Eta axis used
538   Double_t   fLowCut;   // Low cut used when fitting 
539
540   ClassDef(AliFMDCorrELossFit,2); 
541 };
542
543 //____________________________________________________________________
544 inline void 
545 AliFMDCorrELossFit::SetEtaAxis(Int_t nBins, Double_t min, Double_t max)
546 {
547   fEtaAxis.Set(nBins, min, max);
548 }
549 //____________________________________________________________________
550 inline void 
551 AliFMDCorrELossFit::SetEtaAxis(const TAxis& e)
552 {
553   fEtaAxis.Set(e.GetNbins(), e.GetXmin(), e.GetXmax());
554 }
555 //____________________________________________________________________
556 inline Double_t
557 AliFMDCorrELossFit::ELossFit::GetA(UShort_t i) const
558 {
559   if (i <  1)   return 0;
560   if (i >  fN)  return 0;
561   if (i == 1)   return 1;
562   return fA[i-2];
563 }
564 //____________________________________________________________________
565 inline Double_t
566 AliFMDCorrELossFit::ELossFit::GetEA(UShort_t i) const
567 {
568   if (i <  1)   return 0;
569   if (i >  fN)  return 0;
570   if (i == 1)   return 1;
571   return fEA[i-2];
572 }
573
574
575 #endif
576 // Local Variables:
577 //   mode: C++ 
578 // End:
579