4002dd6627ace717f07466b40a74cb6c064ae912
[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 
17  */
18 class AliForwardUtil : public TObject
19 {
20 public:
21   //==================================================================
22   /** 
23    * @{ 
24    * @name 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 cms  Center of mass energy per nucleon
73    * 
74    * @return String representation of the center of mass energy per nuclean
75    */
76   static const char* CenterOfMassEnergyString(UShort_t cms);
77   //__________________________________________________________________
78   /** 
79    * Parse the magnetic field (in kG) as given by a floating point number
80    * 
81    * @param field  Magnetic field in kG 
82    * 
83    * @return Short integer value of magnetic field in kG 
84    */
85   static Short_t ParseMagneticField(Float_t field);
86   /** 
87    * Get a string representation of the magnetic field
88    * 
89    * @param field Magnetic field in kG
90    * 
91    * @return String representation of the magnetic field
92    */
93   static const char* MagneticFieldString(Short_t field);
94   /* @} */
95
96   /** 
97    * @{ 
98    * @name Energy stragling functions 
99    */
100   //__________________________________________________________________
101   /**
102    * Number of steps to do in the Landau, Gaussiam convolution 
103    */
104   static Int_t fgConvolutionSteps;
105   //------------------------------------------------------------------
106   /** 
107    * How many sigma's of the Gaussian in the Landau, Gaussian
108    * convolution to integrate over
109    */
110   static Double_t fgConvolutionNSigma;
111   //------------------------------------------------------------------
112   /** 
113    * Calculate the shifted Landau
114    * @f[
115    *    f'_{L}(x;\Delta,\xi) = f_L(x;\Delta+0.22278298\xi)
116    * @f]
117    *
118    * where @f$ f_{L}@f$ is the ROOT implementation of the Landau
119    * distribution (known to have @f$ \Delta_{p}=-0.22278298@f$ for
120    * @f$\Delta=0,\xi=1@f$. 
121    *
122    * @param x      Where to evaluate @f$ f'_{L}@f$ 
123    * @param delta  Most probable value 
124    * @param xi     The 'width' of the distribution 
125    *
126    * @return @f$ f'_{L}(x;\Delta,\xi) @f$
127    */
128   static Double_t Landau(Double_t x, Double_t delta, Double_t xi);
129   
130   //------------------------------------------------------------------
131   /** 
132    * Calculate the value of a Landau convolved with a Gaussian 
133    * 
134    * @f[ 
135    * f(x;\Delta,\xi,\sigma') = \frac{1}{\sigma' \sqrt{2 \pi}}
136    *    \int_{-\infty}^{+\infty} d\Delta' f'_{L}(x;\Delta',\xi)
137    *    \exp{-\frac{(\Delta-\Delta')^2}{2\sigma'^2}}
138    * @f]
139    * 
140    * where @f$ f'_{L}@f$ is the Landau distribution, @f$ \Delta@f$ the
141    * energy loss, @f$ \xi@f$ the width of the Landau, and 
142    * @f$ \sigma'^2=\sigma^2-\sigma_n^2 @f$.  Here, @f$\sigma@f$ is the
143    * variance of the Gaussian, and @f$\sigma_n@f$ is a parameter modelling 
144    * noise in the detector.  
145    *
146    * Note that this function uses the constants fgConvolutionSteps and
147    * fgConvolutionNSigma
148    * 
149    * References: 
150    *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
151    *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
152    *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
153    * 
154    * @param x         where to evaluate @f$ f@f$
155    * @param delta     @f$ \Delta@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
156    * @param xi        @f$ \xi@f$ of @f$ f(x;\Delta,\xi,\sigma')@f$
157    * @param sigma     @f$ \sigma@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
158    * @param sigma_n   @f$ \sigma_n@f$ of @f$\sigma'^2=\sigma^2-\sigma_n^2 @f$
159    * 
160    * @return @f$ f@f$ evaluated at @f$ x@f$.  
161    */
162   static Double_t LandauGaus(Double_t x, Double_t delta, Double_t xi, 
163                              Double_t sigma, Double_t sigma_n);
164
165   //------------------------------------------------------------------
166   /** 
167    * Evaluate 
168    * @f[ 
169    *    f_i(x;\Delta,\xi,\sigma') = f(x;\Delta_i,\xi_i,\sigma_i')
170    * @f] 
171    * corresponding to @f$ i@f$ particles i.e., with the substitutions 
172    * @f{eqnarray*}{ 
173    *    \Delta    \rightarrow \Delta_i    &=& i(\Delta + \xi\log(i))\\
174    *    \xi       \rightarrow \xi_i       &=& i \xi\\
175    *    \sigma    \rightarrow \sigma_i    &=& \sqrt{i}\sigma\\
176    *    \sigma'^2 \rightarrow \sigma_i'^2 &=& \sigma_n^2 + \sigma_i^2
177    * @f} 
178    * 
179    * @param x        Where to evaluate 
180    * @param delta    @f$ \Delta@f$ 
181    * @param xi       @f$ \xi@f$ 
182    * @param sigma    @f$ \sigma@f$ 
183    * @param sigma_n  @f$ \sigma_n@f$
184    * @param i        @f$ i @f$
185    * 
186    * @return @f$ f_i @f$ evaluated
187    */  
188   static Double_t ILandauGaus(Double_t x, Double_t delta, Double_t xi, 
189                               Double_t sigma, Double_t sigma_n, Int_t i);
190
191   //------------------------------------------------------------------
192   /** 
193    * Numerically evaluate 
194    * @f[ 
195    *    \left.\frac{\partial f_i}{\partial p_i}\right|_{x}
196    * @f] 
197    * where @f$ p_i@f$ is the @f$ i^{\mbox{th}}@f$ parameter.  The mapping 
198    * of the parameters is given by 
199    *
200    * - 0: @f$\Delta@f$ 
201    * - 1: @f$\xi@f$ 
202    * - 2: @f$\sigma@f$ 
203    * - 3: @f$\sigma_n@f$ 
204    *
205    * This is the partial derivative with respect to the parameter of
206    * the response function corresponding to @f$ i@f$ particles i.e.,
207    * with the substitutions
208    * @f[ 
209    *    \Delta    \rightarrow \Delta_i    = i(\Delta + \xi\log(i))\\
210    *    \xi       \rightarrow \xi_i       = i \xi\\
211    *    \sigma    \rightarrow \sigma_i    = \sqrt{i}\sigma\\
212    *    \sigma'^2 \rightarrow \sigma_i'^2 = \sigma_n^2 + \sigma_i^2
213    * @f] 
214    * 
215    * @param x        Where to evaluate 
216    * @param ipar     Parameter number 
217    * @param dp       @f$ \epsilon\delta p_i@f$ for some value of @f$\epsilon@f$
218    * @param delta    @f$ \Delta@f$ 
219    * @param xi       @f$ \xi@f$ 
220    * @param sigma    @f$ \sigma@f$ 
221    * @param sigma_n  @f$ \sigma_n@f$
222    * @param i        @f$ i@f$
223    * 
224    * @return @f$ f_i@f$ evaluated
225    */  
226   static Double_t IdLandauGausdPar(Double_t x, UShort_t ipar, Double_t dp,
227                                    Double_t delta, Double_t xi, 
228                                    Double_t sigma, Double_t sigma_n, Int_t i);
229
230   //------------------------------------------------------------------
231   /** 
232    * Evaluate 
233    * @f[ 
234    *   f_N(x;\Delta,\xi,\sigma') = \sum_{i=1}^N a_i f_i(x;\Delta,\xi,\sigma'a)
235    * @f] 
236    * 
237    * where @f$ f(x;\Delta,\xi,\sigma')@f$ is the convolution of a
238    * Landau with a Gaussian (see LandauGaus).  Note that 
239    * @f$ a_1 = 1@f$, @f$\Delta_i = i(\Delta_1 + \xi\log(i))@f$, 
240    * @f$\xi_i=i\xi_1@f$, and @f$\sigma_i'^2 = \sigma_n^2 + i\sigma_1^2@f$. 
241    *  
242    * References: 
243    *  - <a href="http://dx.doi.org/10.1016/0168-583X(84)90472-5">Nucl.Instrum.Meth.B1:16</a>
244    *  - <a href="http://dx.doi.org/10.1103/PhysRevA.28.615">Phys.Rev.A28:615</a>
245    *  - <a href="http://root.cern.ch/root/htmldoc/tutorials/fit/langaus.C.html">ROOT implementation</a>
246    * 
247    * @param x        Where to evaluate @f$ f_N@f$
248    * @param delta    @f$ \Delta_1@f$ 
249    * @param xi       @f$ \xi_1@f$
250    * @param sigma    @f$ \sigma_1@f$ 
251    * @param sigma_n  @f$ \sigma_n@f$ 
252    * @param n        @f$ N@f$ in the sum above.
253    * @param a        Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
254    *                 @f$ i > 1@f$ 
255    * 
256    * @return @f$ f_N(x;\Delta,\xi,\sigma')@f$ 
257    */
258   static Double_t NLandauGaus(Double_t x, Double_t delta, Double_t xi, 
259                               Double_t sigma, Double_t sigma_n, Int_t n, 
260                               Double_t* a);
261   /** 
262    * Generate a TF1 object of @f$ f_I@f$ 
263    * 
264    * @param c        Constant
265    * @param delta    @f$ \Delta@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 i        @f$ i@f$ - the number of particles
270    * @param xmin     Least value of range
271    * @param xmax     Largest value of range
272    * 
273    * @return Newly allocated TF1 object
274    */
275   static TF1* MakeILandauGaus(Double_t c, 
276                               Double_t delta, Double_t xi, 
277                               Double_t sigma, Double_t sigma_n,
278                               Int_t    i, 
279                               Double_t xmin,  Double_t  xmax);
280   /** 
281    * Generate a TF1 object of @f$ f_N@f$ 
282    * 
283    * @param c         Constant                         
284    * @param delta     @f$ \Delta@f$                    
285    * @param xi        @f$ \xi_1@f$                     
286    * @param sigma     @f$ \sigma_1@f$                  
287    * @param sigma_n   @f$ \sigma_n@f$                  
288    * @param n         @f$ N@f$ - how many particles to sum to
289    * @param a         Array of size @f$ N-1@f$ of the weights @f$ a_i@f$ for 
290    *                  @f$ i > 1@f$ 
291    * @param xmin      Least value of range  
292    * @param xmax      Largest value of range
293    * 
294    * @return Newly allocated TF1 object
295    */
296   static TF1* MakeNLandauGaus(Double_t c, 
297                               Double_t delta, Double_t  xi, 
298                               Double_t sigma, Double_t  sigma_n,
299                               Int_t    n,     Double_t* a, 
300                               Double_t xmin,  Double_t  xmax);
301                                                     
302   //__________________________________________________________________
303   /** 
304    * Structure to do fits to the energy loss spectrum 
305    * 
306    * @ingroup pwg2_forward 
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 
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   /**
455    * Base class for structure holding ring specific histograms
456    * 
457    * @ingroup pwg2_forward 
458    */
459   struct RingHistos : public TObject
460   {
461     /** 
462      * Constructor
463      * 
464      */
465     RingHistos() : fDet(0), fRing('\0'), fName("") {}
466     /** 
467      * 
468      * 
469      * @param d Detector
470      * @param r Ring 
471      */
472     RingHistos(UShort_t d, Char_t r) 
473       : fDet(d), fRing(r), fName(TString::Format("FMD%d%c", d, r)) 
474     {}
475     /** 
476      * Copy constructor
477      * 
478      * @param o Object to copy from 
479      */
480     RingHistos(const RingHistos& o) 
481       : TObject(o), fDet(o.fDet), fRing(o.fRing), fName(o.fName)
482     {}
483     /** 
484      * 
485      */
486     virtual ~RingHistos() {}
487     /** 
488      * Assignement operator
489      * 
490      * @param o Object to assign from
491      * 
492      * @return Reference to this
493      */
494     RingHistos& operator=(const RingHistos& o) 
495     {
496       TObject::operator=(o);
497       fDet  = o.fDet;
498       fRing = o.fRing;
499       fName = o.fName;
500       return *this;
501     }
502     /** 
503      * 
504      * 
505      * @param d 
506      * 
507      * @return 
508      */
509     TList* DefineOutputList(TList* d) const;
510     /** 
511      * 
512      * 
513      * @param d 
514      * 
515      * @return 
516      */
517     TList* GetOutputList(TList* d) const;
518     /** 
519      * 
520      * 
521      * @param d 
522      * @param name 
523      * 
524      * @return 
525      */
526     TH1* GetOutputHist(TList* d, const char* name) const;
527     /** 
528      * 
529      * 
530      * 
531      * @return 
532      */
533     Color_t Color() const 
534     { 
535       return ((fDet == 1 ? kRed : (fDet == 2 ? kGreen : kBlue))
536               + ((fRing == 'I' || fRing == 'i') ? 2 : -2));
537     }
538     UShort_t fDet;   // Detector
539     Char_t   fRing;  // Ring
540     TString  fName;  // Name
541
542     ClassDef(RingHistos,1) 
543   };
544   /* @} */
545     
546 };
547
548 #endif
549 // Local Variables:
550 //  mode: C++
551 // End:
552