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