MarlinTPC  1.2.0
MMPadResponseFunction.h
1 #ifndef MMPadResponseFunction_H
2 #define MMPadResponseFunction_H
3 
4 // C++
5 #include <iostream>
6 #include <string>
7 #include <vector>
8 #include <map>
9 
10 // ROOT
11 #include "Math/ParamFunctor.h"
12 #include <Minuit2/FCNBase.h>
13 #include <TMath.h>
14 
15 namespace marlintpc
16 {
17 
18  class MMPadResponseFunction : public ROOT::Minuit2::FCNBase //Make this object a subclass of FCNBase so it can be minimized
19  {
20 
21  /* Class for the Micromegas Pad Response Function. This class holds general PRF data, and has a chi2 function that can be used to
22  * fit a range of PRF functional forms to that data, using a range of normalization options. It also has functionality to evaluate a given PRF function
23  * on demand, and normalize a given set of data.
24  *
25  *
26  * Author: P. Hayman, Carleton University
27  */
28 
29  public:
30 
31  //Various fit options
32  static const int FIT_MEAN = 0; //Fit just the mean
33  static const int FIT_PARAMS = 1; //Fit just the parameters (obsolete)
34  static const int FIT_BOTH = 2; //Fit both the mean and the parameters
35  static const int SUM_FIT = 3; //For calibration: fits both mean and parameters, and uses the sum of charges as the normalization
36  static const int GAUS_FIT = 4; //For calibration: fits both mean and parameters, and uses the amplitude of a gaussian fit as the normalization
37  static const int CHI2_MIN = 5; //For calibration: fits both mean and parameters, and uses the old min-chi2 definition of the normalization
38 
39  //This is the format in which the PRF stores its data. This format allows the PRF to hold data ranging from a single hit to an entire detector.
40  typedef std::vector<std::vector<std::pair<double, float> > > PointsByRow;
41 
42  //The default settings are the classic "product" form for the PRF function, and the option to only fit the mean, as in the hit finder.
43  MMPadResponseFunction(std::string prfType = "product", int fitType = MMPadResponseFunction::FIT_MEAN)
44  {
45  _PRFType = prfType;
46  _fitType = fitType;
47  x0 = 0;
48  _restrict = false;
49 
50  //Always initialize the params array
51  int nPar = getNPar(_PRFType);
52  _params = new double[nPar];
53 
54  for(int i = 0; i < nPar; i++)
55  {
56  _params[i] = 0.;
57  }
58  }
59 
60  ~MMPadResponseFunction()
61  {
62  if(_params != 0)
63  {
64  delete _params;
65  _params = 0;
66  }
67  }
68  MMPadResponseFunction( const MMPadResponseFunction & other)
69  {
70  _theErrorDef = other._theErrorDef;
71  x0 = other.x0;
72  _fitType = other._fitType;
73  _PRFType = other._PRFType;
74  _restrict = other._restrict;
75  _dataByRow = other._dataByRow;
76 
77  int nPar = getNPar(_PRFType);
78 
79  _params = new double[nPar];
80 
81  for(int i = 0; i < nPar; i++)
82  {
83  _params[i] = other._params[i];
84  }
85  }
86 
87  MMPadResponseFunction &operator=( const MMPadResponseFunction & other )
88  {
89  _theErrorDef = other._theErrorDef;
90  x0 = other.x0;
91  _fitType = other._fitType;
92  _PRFType = other._PRFType;
93  _restrict = other._restrict;
94  _dataByRow = other._dataByRow;
95 
96  int nPar = getNPar(_PRFType);
97 
98  delete [] _params;
99  _params = new double[nPar];
100 
101  for(int i = 0; i < nPar; i++)
102  {
103  _params[i] = other._params[i];
104  }
105 
106  return *this;
107  }
108 
109  void SetData(std::vector< std::pair<double, float> >& aVec); //This will convert the given vector to a PointsByRow variable
110  void SetData(PointsByRow& aVec);
111 
112  void SetType(std::string prfType)
113  {
114  _PRFType = prfType;
115  }
116 
117  void SetFit(int fitType)
118  {
119  _fitType = fitType;
120  }
121 
122  //This option allows the user to restrict minimization of the PRF to +/- 5mm of x0. Could be useful.
123  void SetRestrictDomain(bool rest)
124  {
125  _restrict = rest;
126  }
127 
128  //These functions allow for the normalization of any given data. C.f. MMScatterProcessor
129  void normalizeData(PointsByRow& data);
130  void normalizeData(std::vector< std::pair<double, float> >& data);
131 
132  //This form is necessary for cases where the parameters are being provided by Minuit. C.f. MMRefitTool.
133  void SetParameters(const double * x, int nPar)
134  {
135  _params = new double[nPar];
136 
137  for(int i = 0; i < nPar; i++)
138  {
139  _params[i] = x[i];
140  }
141 
142  }
143 
144  // Required functions from inheritance ///////////////////////////////////////////////////
145  virtual double Up() const
146  {
147  return _theErrorDef;
148  }
149 
150  double operator()(const std::vector<double>& x) const
151  {
152  double * xx = new double[x.size()];
153  for(unsigned int i = 0; i < x.size(); i++) xx[i] = x[i];
154  return (*this)(xx);
155  };
156 
157  //The main chi2 function
158  double operator()(const double* x) const ;
159 
160  //For a 1D minimization
161  double operator()(const double xx) const
162  {
163  return (*this)(&xx);
164  };
165 
166  void SetErrorDef(double def)
167  {
168  _theErrorDef = def;
169  }
171 
172  // Important functions related to the PRF functional form ////////////////////////////////
173 
174  //This allows other classes to check the number of parameters required for a given PRF form
175  static int getNPar(const std::string & type)
176  {
177 
178  std::map<std::string, std::pair<int, int> >::const_iterator it = MMPadResponseFunction::_validPRFs.find(type);
179  if(it != _validPRFs.end())
180  {
181  return it->second.second;
182  }
183  else
184  return -1;
185 
186  }
187 
188  //This is a map of the defined PRF functional forms. The map includes an index (for faster retrieval in the prf function call)
189  //and the number of parameters required
190  static std::map<std::string, std::pair<int, int> > _validPRFs; //Variable
191  static std::map<std::string, std::pair<int, int> > generatePRFMap() //Function
192  {
193  std::map<std::string, std::pair<int, int> > retVal;
194 
195  int index = 0;
196 
197  retVal["product"] = std::pair<int, int>(index, 2); //0
198  retVal["gaus*lor"] = std::pair<int, int>(index++, 2);
199  retVal["product2"] = std::pair<int, int>(index, 3); //1
200  retVal["gaus*lor2"] = std::pair<int, int>(index++, 3);
201  retVal["gaus+lor2"] = std::pair<int, int>(index++, 2); //2
202  retVal["gaus"] = std::pair<int, int>(index++, 1); //3
203  retVal["gaus2"] = std::pair<int, int>(index++, 1); //4
204  retVal["lor"] = std::pair<int, int>(index++, 1); //5
205  retVal["sech"] = std::pair<int, int>(index++, 1); //6
206  retVal["sec"] = std::pair<int, int>(index++, 1); //7
207  retVal["sech+sec"] = std::pair<int, int>(index++, 2); //8
208  retVal["sech+sec2"] = std::pair<int, int>(index++, 2); //9
209  retVal["sec*lor"] = std::pair<int, int>(index++, 2); //10
210  retVal["sec*gaus"] = std::pair<int, int>(index++, 2); //11
211  retVal["sec+lor"] = std::pair<int, int>(index++, 2); //12
212  retVal["gaus+lor"] = std::pair<int, int>(index++, 2); //13
213  retVal["gaus+lor3"] = std::pair<int, int>(index++, 3); //14
214  retVal["gaus+gaus"] = std::pair<int, int>(index++, 2); //15
215  retVal["gaus+sech"] = std::pair<int, int>(index++, 2); //16
216  retVal["gaus+sec2"] = std::pair<int, int>(index++, 2); //17
217  retVal["gaus+sec3"] = std::pair<int, int>(index++, 3); //18
218  retVal["lor+lor"] = std::pair<int, int>(index++, 2); //19
219  retVal["sec+sec"] = std::pair<int, int>(index++, 2); //20
220 
221  return retVal;
222  }
223 
224  //The PRF function itself.
225  double prf(const double c, const double * x) const
226  {
227  //This retrieves the index of the desired PRF
228  static std::map<std::string, std::pair<int, int> >::const_iterator it = MMPadResponseFunction::_validPRFs.find(_PRFType);
229 
230  if(it != _validPRFs.end())
231  {
232  //The index is then compared against the known options to retrieve the correct function. The integer index comparison is
233  //much faster than comparing against a string
234  switch(it->second.first)
235  {
236  case 0:
237  {
238  static const double log2 = log(2.);
239  return exp(-4.*log2*(1.-x[0])*c*c/(x[1]*x[1]))/ (1. + 4. * x[0] * c*c / (x[1]*x[1]) );
240  }
241 
242  case 1:
243  {
244  static const double log2 = log(2.);
245  return exp(-4.*log2*(1.-x[0])*c*c/(x[1]*x[1]))/ (1. + 4. * x[0] * c*c / (x[2]*x[2]) );
246  }
247 
248  case 2:
249  return ( exp(-(c*c)/(x[1]*x[1]) )*(1. - x[0]) + x[0]/(1.+(c*c)/(x[1]*x[1]) ) );
250 
251  case 3:
252  return exp(-c*c/(x[0]*x[0]));
253 
254  case 4:
255  {
256  double tmp = exp(-c*c/(x[0]*x[0]));
257  return 2.*tmp/(tmp + 1.);
258  }
259 
260  case 5:
261  return 1./(1. + c*c/(x[0]*x[0]));
262 
263  case 6:
264  {
265  double tmp = exp(c*x[0]);
266  return 2.*tmp/(tmp*tmp + 1.);
267  }
268 
269  case 7:
270  {
271  double tmp = exp(c*x[0]);
272  return 4.*tmp/((tmp + 1.)*(tmp + 1.));
273  }
274 
275  case 8:
276  {
277  double tmp1 = exp(c*x[0]);
278  double tmp2 = exp(c*x[1]);
279  return (tmp1/(tmp1*tmp1 + 1.) + 2.*tmp2/((tmp2 + 1.)*(tmp2 + 1.)));
280  }
281 
282  case 9:
283  {
284  double tmp1 = exp(c/x[1]);
285  double tmp2 = exp(c*x[1]);
286  return (2.*tmp1/(tmp1*tmp1 + 1.)*(1. - x[0]) + 4.*x[0]*tmp2/((tmp2 + 1.)*(tmp2 + 1.)));
287  }
288 
289  case 10:
290  {
291  double tmp = exp(-c/x[0]);
292  return 4.*tmp/((tmp + 1.)*(tmp + 1.)*(1. + c*c/x[1]));
293  }
294 
295  case 11:
296  {
297  double tmp = exp(-c*x[0]);
298  return 4.*tmp/((tmp + 1.)*(tmp + 1.))*(exp(-c*c/(x[1]*x[1])));
299  }
300 
301  case 12:
302  {
303  double tmp = exp(-c*x[0]);
304  return 0.5*(4.*tmp/((tmp + 1.)*(tmp + 1.)) + 1./(1. + c*c/(x[1]*x[1])));
305  }
306 
307  case 13:
308  return 0.5 * ( TMath::Exp(-(c*c)/(x[0]*x[0]) ) + 1./(1.+(c*c)/(x[1]*x[1]) ) );
309 
310  case 14:
311  return ( TMath::Exp(-(c*c)/(x[1]*x[1]) )*(1. - x[0]) + x[0]/(1.+(c*c)/(x[2]*x[2]) ) );
312 
313  case 15:
314  return 0.5 * ( TMath::Exp(-(c*c)/(x[0]*x[0]) ) + TMath::Exp(-(c*c)/(x[1]*x[1]) ) );
315 
316  case 16:
317  {
318  double tmp = exp(-c*x[1]);
319  return 0.5*(exp(-(c*c)/(x[0]*x[0])) + 2.*tmp/(tmp*tmp + 1.));
320  }
321 
322  case 17:
323  {
324  double tmp = exp(-c*x[1]);
325  return 0.5*(exp(-(c*c)/(x[0]*x[0])) + 4.*tmp/((tmp + 1.)*(tmp + 1.)));
326  }
327 
328  case 18:
329  {
330  double tmp = exp(-c*x[2]);
331  return (TMath::Exp(-(c*c)/(x[1]*x[1]))*(1. - x[0]) + 4.*x[0]*tmp/((tmp + 1.)*(tmp + 1.)));
332  }
333 
334  case 19:
335  return 0.5 * ( 1./(1.+(c*c)*(x[0]*x[0])) + 1./(1.+(c*c)/(x[1]*x[1])) );
336 
337  case 20:
338  {
339  double tmp1 = exp(c*x[0]);
340  double tmp2 = exp(c*x[1]);
341  return 2.*tmp1/((tmp1 + 1.)*(tmp1 + 1.)) + 2.*tmp2/((tmp2 + 1.)*(tmp2 + 1.));
342  }
343  default:
344  return -9E9;
345  }
346  }
347  else
348  {
349  std::cout << "PRF NOT DEFINED!!!" << std::endl;
350  return -99999999.;
351  }
352 
353  }//End PRF Function
354 
355  private:
356 
357  //Needed for FCN_Base inheritance
358  double _theErrorDef;
359 
360  //Mean of the prf
361  double x0;
362 
363  //PRF Parameters
364  double* _params;
365 
366  //Fit type
367  int _fitType;
368 
369  //PRF Type
370  std::string _PRFType;
371 
372  //Restrict domain in fitting for parameters.
373  bool _restrict;
374 
375  //The data to minimize
376  PointsByRow _dataByRow;
377 
378  //Helper function to normalize the amplitudes
379  double calculateRowAmp(int row, double x0, const double *) const;
380  double calculateRowAmp(int row, double x0, const double *, const PointsByRow& data) const;
381  double calculateRowAmp(double x0, const double * x, const std::vector< std::pair<double, float> >& data) const;
382 
383 
384  };//End PRF class definition
385 
386 }//End namespace marlintpc
387 
388 #endif
marlintpc::MMPadResponseFunction
Definition: MMPadResponseFunction.h:18
Generated on Fri Dec 2 2016 12:42:17 for MarlinTPC by   doxygen 1.8.6