Various improvements
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliForwardUtil.h
1 // 
2 // Utilities used in the forward multiplcity analysis 
3 // 
4 //
5 #ifndef ALIFORWARDUTIL_H
6 #define ALIFORWARDUTIL_H
7 #include <TObject.h>
8 #include <TString.h>
9 #include <TObjArray.h>
10 class TH2D;
11 class TH1I;
12 class TH1;
13 class TF1;
14 class TAxis;
15 class AliESDEvent;
16
17 /** 
18  * Utilities used in the forward multiplcity analysis 
19  * 
20  * @ingroup pwg2_forward 
21  */
22 class AliForwardUtil : public TObject
23 {
24 public:
25     /** 
26      * Get the standard color for a ring  
27      *
28      * @param d Detector
29      * @param r Ring 
30      * 
31      * @return 
32      */
33   static Color_t RingColor(UShort_t d, Char_t r)
34   { 
35     return ((d == 1 ? kRed : (d == 2 ? kGreen : kBlue))
36             + ((r == 'I' || r == 'i') ? 2 : -2));
37   }
38   //==================================================================
39   /** 
40    * @{ 
41    * @name Collision/run parameters 
42    */
43   /**                                           
44    * Defined collision types 
45    */
46   enum ECollisionSystem {
47     kUnknown, 
48     kPP, 
49     kPbPb
50   };
51   //__________________________________________________________________
52   /** 
53    * Parse a collision system spec given in a string.   Known values are 
54    * 
55    *  - "pp", "p-p" which returns kPP 
56    *  - "PbPb", "Pb-Pb", "A-A", which returns kPbPb 
57    *  - Everything else gives kUnknown 
58    * 
59    * @param sys Collision system spec 
60    * 
61    * @return Collision system id 
62    */
63   static UShort_t ParseCollisionSystem(const char* sys);
64   /** 
65    * Get a string representation of the collision system 
66    * 
67    * @param sys  Collision system 
68    * - kPP -> "pp"
69    * - kPbPb -> "PbPb" 
70    * - anything else gives "unknown"
71    * 
72    * @return String representation of the collision system 
73    */
74   static const char* CollisionSystemString(UShort_t sys);
75   //__________________________________________________________________
76   /** 
77    * Parse the center of mass energy given as a float and return known 
78    * values as a unsigned integer
79    * 
80    * @param sys  Collision system (needed for AA)
81    * @param cms  Center of mass energy * total charge 
82    * 
83    * @return Center of mass energy per nucleon
84    */
85   static UShort_t ParseCenterOfMassEnergy(UShort_t sys, Float_t cms);
86   /** 
87    * Get a string representation of the center of mass energy per nuclean
88    * 
89    * @param cms  Center of mass energy per nucleon
90    * 
91    * @return String representation of the center of mass energy per nuclean
92    */
93   static const char* CenterOfMassEnergyString(UShort_t cms);
94   //__________________________________________________________________
95   /** 
96    * Parse the magnetic field (in kG) as given by a floating point number
97    * 
98    * @param field  Magnetic field in kG 
99    * 
100    * @return Short integer value of magnetic field in kG 
101    */
102   static Short_t ParseMagneticField(Float_t field);
103   /** 
104    * Get a string representation of the magnetic field
105    * 
106    * @param field Magnetic field in kG
107    * 
108    * @return String representation of the magnetic field
109    */
110   static const char* MagneticFieldString(Short_t field);
111   /* @} */
112
113   /** 
114    * @{ 
115    * @name Energy stragling functions 
116    */
117   //__________________________________________________________________
118   /**
119    * Number of steps to do in the Landau, Gaussiam convolution 
120    */
121   static Int_t fgConvolutionSteps;
122   //------------------------------------------------------------------
123   /** 
124    * How many sigma's of the Gaussian in the Landau, Gaussian
125    * convolution to integrate over
126    */
127   static Double_t fgConvolutionNSigma;
128   //------------------------------------------------------------------
129   /** 
130    * Calculate the shifted Landau
131    * @f[
132    *    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
133    * @f]
134    *
135    * where @f$ f_{L}@f$ is the ROOT implementation of the Landau
136    * distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
137    * @f$\Delta=0,\xi=1@f$. 
138    *
139    * @param x      Where to evaluate @f$ f'_{L}@f$ 
140    * @param delta  Most probable value 
141    * @param xi     The 'width' of the distribution 
142    *
143    * @return @f$ f'_{L}(x;\Delta,\xi) @f$
144    */
145   static Double_t Landau(Double_t x, Double_t delta, Double_t xi);
146   
147   //------------------------------------------------------------------
148   /** 
149    * Calculate the value of a Landau convolved with a Gaussian 
150    * 
151    * @f[ 
152    * f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
153    *    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
154    *    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
155    * @f]
156    * 
157    * where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
158    * energy loss, @f$ \xi@f$ the width of the Landau, and 
159    * @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
160    * variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling 
161    * noise in the detector.  
162    *
163    * Note that this function uses the constants fgConvolutionSteps and
164    * fgConvolutionNSigma
165    * 
166    * References: 
167    *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
168    *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
169    *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
170    * 
171    * @param x         where to evaluate @f$ f@f$
172    * @param delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
173    * @param xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
174    * @param sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
175    * @param sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
176    * 
177    * @return @f$ f@f$ evaluated at @f$ x@f$.  
178    */
179   static Double_t LandauGaus(Double_t x, Double_t delta, Double_t xi, 
180                              Double_t sigma, Double_t sigma_n);
181
182   //------------------------------------------------------------------
183   /** 
184    * Evaluate 
185    * @f[ 
186    *    f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
187    * @f] 
188    * corresponding to @f$ i@f$ particles i.e., with the substitutions 
189    * @f{eqnarray*}{ 
190    *    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))\\
191    *    \xi       \rightarrow \xi_i       &=& i \xi\\
192    *    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma\\
193    *    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
194    * @f} 
195    * 
196    * @param x        Where to evaluate 
197    * @param delta    @f$ \Delta@f$ 
198    * @param xi       @f$ \xi@f$ 
199    * @param sigma    @f$ \sigma@f$ 
200    * @param sigma_n  @f$ \sigma_n@f$
201    * @param i        @f$ i @f$
202    * 
203    * @return @f$ f_i @f$ evaluated
204    */  
205   static Double_t ILandauGaus(Double_t x, Double_t delta, Double_t xi, 
206                               Double_t sigma, Double_t sigma_n, Int_t i);
207
208   //------------------------------------------------------------------
209   /** 
210    * Numerically evaluate 
211    * @f[ 
212    *    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
213    * @f] 
214    * where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping 
215    * of the parameters is given by 
216    *
217    * - 0: @f$\Delta@f$ 
218    * - 1: @f$\xi@f$ 
219    * - 2: @f$\sigma@f$ 
220    * - 3: @f$\sigma_n@f$ 
221    *
222    * This is the partial derivative with respect to the parameter of
223    * the response function corresponding to @f$ i@f$ particles i.e.,
224    * with the substitutions
225    * @f[ 
226    *    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))\\
227    *    \xi       \rightarrow \xi_i       = i \xi\\
228    *    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma\\
229    *    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
230    * @f] 
231    * 
232    * @param x        Where to evaluate 
233    * @param ipar     Parameter number 
234    * @param dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
235    * @param delta    @f$ \Delta@f$ 
236    * @param xi       @f$ \xi@f$ 
237    * @param sigma    @f$ \sigma@f$ 
238    * @param sigma_n  @f$ \sigma_n@f$
239    * @param i        @f$ i@f$
240    * 
241    * @return @f$ f_i@f$ evaluated
242    */  
243   static Double_t IdLandauGausdPar(Double_t x, UShort_t ipar, Double_t dp,
244                                    Double_t delta, Double_t xi, 
245                                    Double_t sigma, Double_t sigma_n, Int_t i);
246
247   //------------------------------------------------------------------
248   /** 
249    * Evaluate 
250    * @f[ 
251    *   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
252    * @f] 
253    * 
254    * where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
255    * Landau with a Gaussian (see LandauGaus).  Note that 
256    * @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$, 
257    * @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$. 
258    *  
259    * References: 
260    *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
261    *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
262    *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
263    * 
264    * @param x        Where to evaluate @f$ f_N@f$
265    * @param delta    @f$ \Delta_1@f$ 
266    * @param xi       @f$ \xi_1@f$
267    * @param sigma    @f$ \sigma_1@f$ 
268    * @param sigma_n  @f$ \sigma_n@f$ 
269    * @param n        @f$ N@f$ in the sum above.
270    * @param a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
271    *                 @f$ i > 1@f$ 
272    * 
273    * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$ 
274    */
275   static Double_t NLandauGaus(Double_t x, Double_t delta, Double_t xi, 
276                               Double_t sigma, Double_t sigma_n, Int_t n, 
277                               Double_t* a);
278   /** 
279    * Generate a TF1 object of @f$ f_I@f$ 
280    * 
281    * @param c        Constant
282    * @param delta    @f$ \Delta@f$ 
283    * @param xi       @f$ \xi_1@f$              
284    * @param sigma    @f$ \sigma_1@f$           
285    * @param sigma_n  @f$ \sigma_n@f$           
286    * @param i        @f$ i@f$ - the number of particles
287    * @param xmin     Least value of range
288    * @param xmax     Largest value of range
289    * 
290    * @return Newly allocated TF1 object
291    */
292   static TF1* MakeILandauGaus(Double_t c, 
293                               Double_t delta, Double_t xi, 
294                               Double_t sigma, Double_t sigma_n,
295                               Int_t    i, 
296                               Double_t xmin,  Double_t  xmax);
297   /** 
298    * Generate a TF1 object of @f$ f_N@f$ 
299    * 
300    * @param c         Constant                         
301    * @param delta     @f$ \Delta@f$                    
302    * @param xi        @f$ \xi_1@f$                     
303    * @param sigma     @f$ \sigma_1@f$                  
304    * @param sigma_n   @f$ \sigma_n@f$                  
305    * @param n         @f$ N@f$ - how many particles to sum to
306    * @param a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
307    *                  @f$ i > 1@f$ 
308    * @param xmin      Least value of range  
309    * @param xmax      Largest value of range
310    * 
311    * @return Newly allocated TF1 object
312    */
313   static TF1* MakeNLandauGaus(Double_t c, 
314                               Double_t delta, Double_t  xi, 
315                               Double_t sigma, Double_t  sigma_n,
316                               Int_t    n,     Double_t* a, 
317                               Double_t xmin,  Double_t  xmax);
318                                                     
319   //__________________________________________________________________
320   /** 
321    * Structure to do fits to the energy loss spectrum 
322    * 
323    * @ingroup pwg2_forward 
324    */
325   struct ELossFitter 
326   {
327     enum { 
328       kC     = 0,
329       kDelta, 
330       kXi, 
331       kSigma, 
332       kSigmaN, 
333       kN, 
334       kA
335     };
336     /** 
337      * Constructor 
338      * 
339      * @param lowCut     Lower cut of spectrum - data below this cuts is ignored
340      * @param maxRange   Maximum range to fit to 
341      * @param minusBins  The number of bins below maximum to use 
342      */
343     ELossFitter(Double_t lowCut, Double_t maxRange, UShort_t minusBins); 
344     /** 
345      * Destructor
346      * 
347      */
348     virtual ~ELossFitter();
349     /** 
350      * Clear internal arrays 
351      * 
352      */
353     void Clear();
354     /** 
355      * Fit a 1-particle signal to the passed energy loss distribution 
356      * 
357      * Note that this function clears the internal arrays first 
358      *
359      * @param dist    Data to fit the function to 
360      * @param sigman If larger than zero, the initial guess of the
361      *               detector induced noise. If zero or less, then this 
362      *               parameter is ignored in the fit (fixed at 0)
363      * 
364      * @return The function fitted to the data 
365      */
366     TF1* Fit1Particle(TH1* dist, Double_t sigman=-1);
367     /** 
368      * Fit a N-particle signal to the passed energy loss distribution 
369      *
370      * If there's no 1-particle fit present, it does that first 
371      *
372      * @param dist   Data to fit the function to 
373      * @param n      Number of particle signals to fit 
374      * @param sigman If larger than zero, the initial guess of the
375      *               detector induced noise. If zero or less, then this 
376      *               parameter is ignored in the fit (fixed at 0)
377      * 
378      * @return The function fitted to the data 
379      */
380     TF1* FitNParticle(TH1* dist, UShort_t n, Double_t sigman=-1);
381      
382
383     const Double_t fLowCut;     // Lower cut on data 
384     const Double_t fMaxRange;   // Maximum range to fit 
385     const UShort_t fMinusBins;  // Number of bins from maximum to fit 1st peak
386     TObjArray fFitResults;      // Array of fit results 
387     TObjArray fFunctions;       // Array of functions 
388   };
389   /* @} */
390       
391
392   //==================================================================
393   /** 
394    * @{
395    * @name Convenience containers 
396    */
397   /** 
398    * Structure to hold histograms 
399    *
400    * @ingroup pwg2_forward 
401    */
402   struct Histos : public TObject
403   {     
404     /** 
405      * Constructor 
406      * 
407      * 
408      */
409     Histos() : fFMD1i(0), fFMD2i(0), fFMD2o(0), fFMD3i(0), fFMD3o(0) {}
410     /** 
411      * Copy constructor 
412      * 
413      * @param o Object to copy from 
414      */
415     Histos(const Histos& o) 
416       : TObject(o), 
417         fFMD1i(o.fFMD1i), 
418         fFMD2i(o.fFMD2i), 
419         fFMD2o(o.fFMD2o), 
420         fFMD3i(o.fFMD3i), 
421         fFMD3o(o.fFMD3o)
422     {}
423     /** 
424      * Assignement operator 
425      * 
426      * @return Reference to this 
427      */
428     Histos& operator=(const Histos&) { return *this;}
429     /** 
430      * Destructor
431      */
432     ~Histos();
433     /** 
434      * Initialize the object 
435      * 
436      * @param etaAxis Eta axis to use 
437      */
438     void Init(const TAxis& etaAxis);
439     /** 
440      * Make a histogram 
441      * 
442      * @param d        Detector
443      * @param r        Ring 
444      * @param etaAxis  Eta axis to use
445      * 
446      * @return Newly allocated histogram 
447      */
448     TH2D* Make(UShort_t d, Char_t r, const TAxis& etaAxis) const;
449     /** 
450      * Clear data 
451      * 
452      * @param option Not used 
453      */
454     void  Clear(Option_t* option="");
455     // const TH2D* Get(UShort_t d, Char_t r) const;
456     /** 
457      * Get the histogram for a particular detector,ring
458      * 
459      * @param d Detector 
460      * @param r Ring 
461      * 
462      * @return Histogram for detector,ring or nul 
463      */
464     TH2D* Get(UShort_t d, Char_t r) const;
465     TH2D* fFMD1i; // Histogram for FMD1i
466     TH2D* fFMD2i; // Histogram for FMD2i
467     TH2D* fFMD2o; // Histogram for FMD2o
468     TH2D* fFMD3i; // Histogram for FMD3i
469     TH2D* fFMD3o; // Histogram for FMD3o
470
471     ClassDef(Histos,1) 
472   };
473
474   //__________________________________________________________________
475   /**
476    * Base class for structure holding ring specific histograms
477    * 
478    * @ingroup pwg2_forward 
479    */
480   struct RingHistos : public TObject
481   {
482     /** 
483      * Constructor
484      * 
485      */
486     RingHistos() : fDet(0), fRing('\0'), fName("") {}
487     /** 
488      * 
489      * 
490      * @param d Detector
491      * @param r Ring 
492      */
493     RingHistos(UShort_t d, Char_t r) 
494       : fDet(d), fRing(r), fName(TString::Format("FMD%d%c", d, r)) 
495     {}
496     /** 
497      * Copy constructor
498      * 
499      * @param o Object to copy from 
500      */
501     RingHistos(const RingHistos& o) 
502       : TObject(o), fDet(o.fDet), fRing(o.fRing), fName(o.fName)
503     {}
504     /** 
505      * 
506      */
507     virtual ~RingHistos() {}
508     /** 
509      * Assignement operator
510      * 
511      * @param o Object to assign from
512      * 
513      * @return Reference to this
514      */
515     RingHistos& operator=(const RingHistos& o) 
516     {
517       TObject::operator=(o);
518       fDet  = o.fDet;
519       fRing = o.fRing;
520       fName = o.fName;
521       return *this;
522     }
523     /** 
524      * Define the outout list in @a d
525      * 
526      * @param d Where to put the output list
527      * 
528      * @return Newly allocated TList object or null
529      */
530     TList* DefineOutputList(TList* d) const;
531     /** 
532      * Get our output list from the container @a d
533      * 
534      * @param d where to get the output list from 
535      * 
536      * @return The found TList or null
537      */
538     TList* GetOutputList(TList* d) const;
539     /** 
540      * Find a specific histogram in the source list @a d
541      * 
542      * @param d     (top)-container 
543      * @param name  Name of histogram
544      * 
545      * @return Found histogram or null
546      */
547     TH1* GetOutputHist(TList* d, const char* name) const;
548     /** 
549      * 
550      * 
551      * 
552      * @return 
553      */
554     Color_t Color() const 
555     { 
556       return AliForwardUtil::RingColor(fDet, fRing);
557     }
558     UShort_t fDet;   // Detector
559     Char_t   fRing;  // Ring
560     TString  fName;  // Name
561
562     ClassDef(RingHistos,1) 
563   };
564   /* @} */
565     
566 };
567
568 #endif
569 // Local Variables:
570 //  mode: C++
571 // End:
572