1 #ifndef MMPadResponseFunction_H
2 #define MMPadResponseFunction_H
11 #include "Math/ParamFunctor.h"
12 #include <Minuit2/FCNBase.h>
32 static const int FIT_MEAN = 0;
33 static const int FIT_PARAMS = 1;
34 static const int FIT_BOTH = 2;
35 static const int SUM_FIT = 3;
36 static const int GAUS_FIT = 4;
37 static const int CHI2_MIN = 5;
40 typedef std::vector<std::vector<std::pair<double, float> > > PointsByRow;
51 int nPar = getNPar(_PRFType);
52 _params =
new double[nPar];
54 for(
int i = 0; i < nPar; i++)
70 _theErrorDef = other._theErrorDef;
72 _fitType = other._fitType;
73 _PRFType = other._PRFType;
74 _restrict = other._restrict;
75 _dataByRow = other._dataByRow;
77 int nPar = getNPar(_PRFType);
79 _params =
new double[nPar];
81 for(
int i = 0; i < nPar; i++)
83 _params[i] = other._params[i];
89 _theErrorDef = other._theErrorDef;
91 _fitType = other._fitType;
92 _PRFType = other._PRFType;
93 _restrict = other._restrict;
94 _dataByRow = other._dataByRow;
96 int nPar = getNPar(_PRFType);
99 _params =
new double[nPar];
101 for(
int i = 0; i < nPar; i++)
103 _params[i] = other._params[i];
109 void SetData(std::vector< std::pair<double, float> >& aVec);
110 void SetData(PointsByRow& aVec);
112 void SetType(std::string prfType)
117 void SetFit(
int fitType)
123 void SetRestrictDomain(
bool rest)
129 void normalizeData(PointsByRow& data);
130 void normalizeData(std::vector< std::pair<double, float> >& data);
133 void SetParameters(
const double * x,
int nPar)
135 _params =
new double[nPar];
137 for(
int i = 0; i < nPar; i++)
145 virtual double Up()
const
150 double operator()(
const std::vector<double>& x)
const
152 double * xx =
new double[x.size()];
153 for(
unsigned int i = 0; i < x.size(); i++) xx[i] = x[i];
158 double operator()(
const double* x)
const ;
161 double operator()(
const double xx)
const
166 void SetErrorDef(
double def)
175 static int getNPar(
const std::string & type)
178 std::map<std::string, std::pair<int, int> >::const_iterator it = MMPadResponseFunction::_validPRFs.find(type);
179 if(it != _validPRFs.end())
181 return it->second.second;
190 static std::map<std::string, std::pair<int, int> > _validPRFs;
191 static std::map<std::string, std::pair<int, int> > generatePRFMap()
193 std::map<std::string, std::pair<int, int> > retVal;
197 retVal[
"product"] = std::pair<int, int>(index, 2);
198 retVal[
"gaus*lor"] = std::pair<int, int>(index++, 2);
199 retVal[
"product2"] = std::pair<int, int>(index, 3);
200 retVal[
"gaus*lor2"] = std::pair<int, int>(index++, 3);
201 retVal[
"gaus+lor2"] = std::pair<int, int>(index++, 2);
202 retVal[
"gaus"] = std::pair<int, int>(index++, 1);
203 retVal[
"gaus2"] = std::pair<int, int>(index++, 1);
204 retVal[
"lor"] = std::pair<int, int>(index++, 1);
205 retVal[
"sech"] = std::pair<int, int>(index++, 1);
206 retVal[
"sec"] = std::pair<int, int>(index++, 1);
207 retVal[
"sech+sec"] = std::pair<int, int>(index++, 2);
208 retVal[
"sech+sec2"] = std::pair<int, int>(index++, 2);
209 retVal[
"sec*lor"] = std::pair<int, int>(index++, 2);
210 retVal[
"sec*gaus"] = std::pair<int, int>(index++, 2);
211 retVal[
"sec+lor"] = std::pair<int, int>(index++, 2);
212 retVal[
"gaus+lor"] = std::pair<int, int>(index++, 2);
213 retVal[
"gaus+lor3"] = std::pair<int, int>(index++, 3);
214 retVal[
"gaus+gaus"] = std::pair<int, int>(index++, 2);
215 retVal[
"gaus+sech"] = std::pair<int, int>(index++, 2);
216 retVal[
"gaus+sec2"] = std::pair<int, int>(index++, 2);
217 retVal[
"gaus+sec3"] = std::pair<int, int>(index++, 3);
218 retVal[
"lor+lor"] = std::pair<int, int>(index++, 2);
219 retVal[
"sec+sec"] = std::pair<int, int>(index++, 2);
225 double prf(
const double c,
const double * x)
const
228 static std::map<std::string, std::pair<int, int> >::const_iterator it = MMPadResponseFunction::_validPRFs.find(_PRFType);
230 if(it != _validPRFs.end())
234 switch(it->second.first)
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]) );
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]) );
249 return ( exp(-(c*c)/(x[1]*x[1]) )*(1. - x[0]) + x[0]/(1.+(c*c)/(x[1]*x[1]) ) );
252 return exp(-c*c/(x[0]*x[0]));
256 double tmp = exp(-c*c/(x[0]*x[0]));
257 return 2.*tmp/(tmp + 1.);
261 return 1./(1. + c*c/(x[0]*x[0]));
265 double tmp = exp(c*x[0]);
266 return 2.*tmp/(tmp*tmp + 1.);
271 double tmp = exp(c*x[0]);
272 return 4.*tmp/((tmp + 1.)*(tmp + 1.));
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.)));
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.)));
291 double tmp = exp(-c/x[0]);
292 return 4.*tmp/((tmp + 1.)*(tmp + 1.)*(1. + c*c/x[1]));
297 double tmp = exp(-c*x[0]);
298 return 4.*tmp/((tmp + 1.)*(tmp + 1.))*(exp(-c*c/(x[1]*x[1])));
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])));
308 return 0.5 * ( TMath::Exp(-(c*c)/(x[0]*x[0]) ) + 1./(1.+(c*c)/(x[1]*x[1]) ) );
311 return ( TMath::Exp(-(c*c)/(x[1]*x[1]) )*(1. - x[0]) + x[0]/(1.+(c*c)/(x[2]*x[2]) ) );
314 return 0.5 * ( TMath::Exp(-(c*c)/(x[0]*x[0]) ) + TMath::Exp(-(c*c)/(x[1]*x[1]) ) );
318 double tmp = exp(-c*x[1]);
319 return 0.5*(exp(-(c*c)/(x[0]*x[0])) + 2.*tmp/(tmp*tmp + 1.));
324 double tmp = exp(-c*x[1]);
325 return 0.5*(exp(-(c*c)/(x[0]*x[0])) + 4.*tmp/((tmp + 1.)*(tmp + 1.)));
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.)));
335 return 0.5 * ( 1./(1.+(c*c)*(x[0]*x[0])) + 1./(1.+(c*c)/(x[1]*x[1])) );
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.));
349 std::cout <<
"PRF NOT DEFINED!!!" << std::endl;
370 std::string _PRFType;
376 PointsByRow _dataByRow;
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;
Definition: MMPadResponseFunction.h:18