405057eb63362771875fd9092764bd6a8ce4420b
[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    */
307   struct ELossFitter 
308   {
309     enum { 
310       kC     = 0,
311       kDelta, 
312       kXi, 
313       kSigma, 
314       kSigmaN, 
315       kN, 
316       kA
317     };
318     /** 
319      * Constructor 
320      * 
321      * @param lowCut     Lower cut of spectrum - data below this cuts is ignored
322      * @param maxRange   Maximum range to fit to 
323      * @param minusBins  The number of bins below maximum to use 
324      */
325     ELossFitter(Double_t lowCut, Double_t maxRange, UShort_t minusBins); 
326     virtual ~ELossFitter();
327     /** 
328      * Clear internal arrays 
329      * 
330      */
331     void Clear();
332     /** 
333      * Fit a 1-particle signal to the passed energy loss distribution 
334      * 
335      * Note that this function clears the internal arrays first 
336      *
337      * @param dist    Data to fit the function to 
338      * @param sigman If larger than zero, the initial guess of the
339      *               detector induced noise. If zero or less, then this 
340      *               parameter is ignored in the fit (fixed at 0)
341      * 
342      * @return The function fitted to the data 
343      */
344     TF1* Fit1Particle(TH1* dist, Double_t sigman=-1);
345     /** 
346      * Fit a N-particle signal to the passed energy loss distribution 
347      *
348      * If there's no 1-particle fit present, it does that first 
349      *
350      * @param dist   Data to fit the function to 
351      * @param n      Number of particle signals to fit 
352      * @param sigman If larger than zero, the initial guess of the
353      *               detector induced noise. If zero or less, then this 
354      *               parameter is ignored in the fit (fixed at 0)
355      * 
356      * @return The function fitted to the data 
357      */
358     TF1* FitNParticle(TH1* dist, UShort_t n, Double_t sigman=-1);
359      
360
361     const Double_t fLowCut;     // Lower cut on data 
362     const Double_t fMaxRange;   // Maximum range to fit 
363     const UShort_t fMinusBins;  // Number of bins from maximum to fit 1st peak
364     TObjArray fFitResults;      // Array of fit results 
365     TObjArray fFunctions;       // Array of functions 
366   };
367   /* @} */
368       
369
370   //==================================================================
371   /** 
372    * @{
373    * @name Convenience containers 
374    */
375   /** 
376    * Structure to hold histograms 
377    *
378    * @ingroup pwg2_forward 
379    */
380   struct Histos : public TObject
381   {     
382     /** 
383      * Constructor 
384      * 
385      * 
386      */
387     Histos() : fFMD1i(0), fFMD2i(0), fFMD2o(0), fFMD3i(0), fFMD3o(0) {}
388     /** 
389      * Copy constructor 
390      * 
391      * @param o Object to copy from 
392      */
393     Histos(const Histos& o) 
394       : TObject(o), 
395         fFMD1i(o.fFMD1i), 
396         fFMD2i(o.fFMD2i), 
397         fFMD2o(o.fFMD2o), 
398         fFMD3i(o.fFMD3i), 
399         fFMD3o(o.fFMD3o)
400     {}
401     /** 
402      * Assignement operator 
403      * 
404      * @return Reference to this 
405      */
406     Histos& operator=(const Histos&) { return *this;}
407     /** 
408      * Destructor
409      */
410     ~Histos();
411     /** 
412      * Initialize the object 
413      * 
414      * @param etaAxis Eta axis to use 
415      */
416     void Init(const TAxis& etaAxis);
417     /** 
418      * Make a histogram 
419      * 
420      * @param d        Detector
421      * @param r        Ring 
422      * @param etaAxis  Eta axis to use
423      * 
424      * @return Newly allocated histogram 
425      */
426     TH2D* Make(UShort_t d, Char_t r, const TAxis& etaAxis) const;
427     /** 
428      * Clear data 
429      * 
430      * @param option Not used 
431      */
432     void  Clear(Option_t* option="");
433     // const TH2D* Get(UShort_t d, Char_t r) const;
434     /** 
435      * Get the histogram for a particular detector,ring
436      * 
437      * @param d Detector 
438      * @param r Ring 
439      * 
440      * @return Histogram for detector,ring or nul 
441      */
442     TH2D* Get(UShort_t d, Char_t r) const;
443     TH2D* fFMD1i; // Histogram for FMD1i
444     TH2D* fFMD2i; // Histogram for FMD2i
445     TH2D* fFMD2o; // Histogram for FMD2o
446     TH2D* fFMD3i; // Histogram for FMD3i
447     TH2D* fFMD3o; // Histogram for FMD3o
448
449     ClassDef(Histos,1) 
450   };
451
452   //__________________________________________________________________
453   struct RingHistos : public TObject
454   {
455     RingHistos() : fDet(0), fRing('\0'), fName("") {}
456     RingHistos(UShort_t d, Char_t r) 
457       : fDet(d), fRing(r), fName(TString::Format("FMD%d%c", d, r)) 
458     {}
459     RingHistos(const RingHistos& o) 
460       : TObject(o), fDet(o.fDet), fRing(o.fRing), fName(o.fName)
461     {}
462     virtual ~RingHistos() {}
463     RingHistos& operator=(const RingHistos& o) 
464     {
465       TObject::operator=(o);
466       fDet  = o.fDet;
467       fRing = o.fRing;
468       fName = o.fName;
469       return *this;
470     }
471     TList* DefineOutputList(TList* d) const;
472     TList* GetOutputList(TList* d) const;
473     TH1* GetOutputHist(TList* d, const char* name) const;
474     Color_t Color() const 
475     { 
476       return ((fDet == 1 ? kRed : (fDet == 2 ? kGreen : kBlue))
477               + ((fRing == 'I' || fRing == 'i') ? 2 : -2));
478     }
479
480     UShort_t fDet; 
481     Char_t   fRing;
482     TString  fName;
483
484     ClassDef(RingHistos,1) 
485   };
486   /* @} */
487     
488 };
489
490 #endif
491 // Local Variables:
492 //  mode: C++
493 // End:
494