22 return (-
pow(residual - center, 2) / 2. / sigma / sigma) +
log(1. /
sqrt(2.*
M_PI) / sigma);
31 if (gammaoversigma == 0.)
return (-toversigma*toversigma/2.) -
log(
sqrt(2*
M_PI));
35 double integrandplus_last = 0.;
36 double integrandminus_last = 0.;
37 for (
double inc = 0.; uplus <
max; inc +=
step) {
38 double uplus_last = uplus;
39 uplus =
pow(inc, power);
41 double integrandplus =
exp(-
pow(uplus - toversigma, 2)/2.) / (uplus*uplus/gammaoversigma + gammaoversigma) /
M_PI /
sqrt(2.*
M_PI);
42 double integrandminus =
exp(-
pow(-uplus - toversigma, 2)/2.) / (uplus*uplus/gammaoversigma + gammaoversigma) /
M_PI /
sqrt(2.*
M_PI);
44 sum += integrandplus * (uplus - uplus_last);
45 sum += 0.5 * fabs(integrandplus_last - integrandplus) * (uplus - uplus_last);
47 sum += integrandminus * (uplus - uplus_last);
48 sum += 0.5 * fabs(integrandminus_last - integrandminus) * (uplus - uplus_last);
50 integrandplus_last = integrandplus;
51 integrandminus_last = integrandminus;
59 double gammaoversigma = fabs(gamma / sigma);
60 double toversigma = fabs((residual - center) / sigma);
70 if (gsisbad || tsisbad) {
71 return log(gammaoversigma/
M_PI) -
log(toversigma*toversigma + gammaoversigma*gammaoversigma) -
log(sigma);
84 return val -
log(sigma);
94 return log(TMath::Voigt(residual - center, fabs(sigma), fabs(gamma)*2.));
98 return par[0] * TMath::Voigt(xvec[0] - par[1], fabs(par[2]), fabs(par[3])*2.);
102 double x = residual-center;
103 double s = fabs(sigma);
108 if (fabs(x)<m)
return -x*x/(2*s*
s) -
log(n);
109 else return log(a) -4*
log(fabs(x)) -
log(n);
117 return (erf((high + center) /
sqrt(2.) / sigma) - erf((low + center) /
sqrt(2.) / sigma)) *
exp(0.5/sigma/sigma) / 2.;
124 std::ifstream convolution_table(
"convolution_table.txt");
125 if (convolution_table.is_open()) {
128 double tsbinsize = 0.;
129 double gsbinsize = 0.;
131 convolution_table >> numgsbins >> numtsbins >> tsbinsize >> gsbinsize;
134 throw cms::Exception(
"MuonResidualsFitter") <<
"convolution_table.txt has the wrong bin width/bin size. Throw it away and let the fitter re-create the file." << std::endl;
143 convolution_table >> read_gsbin >> read_tsbin >>
value;
144 if (read_gsbin != gsbin || read_tsbin != tsbin) {
145 throw cms::Exception(
"MuonResidualsFitter") <<
"convolution_table.txt is out of order. Throw it away and let the fitter re-create the file." << std::endl;
152 convolution_table.close();
156 std::ofstream convolution_table2(
"convolution_table.txt");
158 if (!convolution_table2.is_open()) {
159 throw cms::Exception(
"MuonResidualsFitter") <<
"Couldn't write to file convolution_table.txt" << std::endl;
164 edm::LogWarning(
"MuonResidualsFitter") <<
"Initializing convolution look-up table (takes a few minutes)..." << std::endl;
165 std::cout <<
"Initializing convolution look-up table (takes a few minutes)..." << std::endl;
170 std::cout <<
" gsbin " << gsbin <<
"/" << MuonResidualsFitter_numgsbins << std::endl;
185 convolution_table2.close();
187 edm::LogWarning(
"MuonResidualsFitter") <<
"Initialization done!" << std::endl;
188 std::cout <<
"Initialization done!" << std::endl;
192 bool MuonResidualsFitter::dofit(
void (*
fcn)(
int&,
double*,
double&,
double*,
int), std::vector<int> &parNum, std::vector<std::string> &parName, std::vector<double> &start, std::vector<double> &
step, std::vector<double> &low, std::vector<double> &high) {
201 std::vector<int>::const_iterator iNum = parNum.begin();
202 std::vector<std::string>::const_iterator iName = parName.begin();
203 std::vector<double>::const_iterator istart = start.begin();
204 std::vector<double>::const_iterator istep = step.begin();
205 std::vector<double>::const_iterator ilow = low.begin();
206 std::vector<double>::const_iterator ihigh = high.begin();
208 for (; iNum != parNum.end(); ++iNum, ++iName, ++istart, ++istep, ++ilow, ++ihigh) {
218 for (
int i = 0;
i < 10;
i++) arglist[
i] = 0.;
226 for (
int i = 0;
i < 10;
i++) arglist[
i] = 0.;
232 bool try_again =
false;
235 for (
int i = 0;
i < 10;
i++) arglist[
i] = 0.;
239 if (ierflg != 0) try_again =
true;
244 for (
int i = 0;
i < 10;
i++) arglist[
i] = 0.;
249 Double_t fmin, fedm, errdef;
250 Int_t npari, nparx, istat;
254 for (
int i = 0;
i < 10;
i++) arglist[
i] = 0.;
264 for (
int i = 0;
i <
npar();
i++) {
273 if (smierflg != 0)
return false;
280 int whichcopy = which;
282 fwrite(&rows,
sizeof(
long), 1, file);
283 fwrite(&cols,
sizeof(
int), 1, file);
284 fwrite(&whichcopy,
sizeof(
int), 1, file);
286 double *likeAChecksum =
new double[cols];
287 double *likeAChecksum2 =
new double[cols];
288 for (
int i = 0;
i < cols;
i++) {
289 likeAChecksum[
i] = 0.;
290 likeAChecksum2[
i] = 0.;
294 fwrite((*residual),
sizeof(
double), cols, file);
296 for (
int i = 0;
i < cols;
i++) {
297 if (fabs((*residual)[
i]) > likeAChecksum[i]) likeAChecksum[
i] = fabs((*residual)[i]);
298 if (fabs((*residual)[i]) < likeAChecksum2[i]) likeAChecksum2[
i] = fabs((*residual)[i]);
304 fwrite(likeAChecksum,
sizeof(
double), cols, file);
305 fwrite(likeAChecksum2,
sizeof(
double), cols, file);
307 delete [] likeAChecksum;
308 delete [] likeAChecksum2;
314 int readwhich = -100;
316 fread(&rows,
sizeof(
long), 1, file);
317 fread(&cols,
sizeof(
int), 1, file);
318 fread(&readwhich,
sizeof(
int), 1, file);
320 if (cols !=
ndata() || rows < 0 || readwhich != which)
throw cms::Exception(
"MuonResidualsFitter") <<
"temporary file is corrupted (which = " << which <<
" readwhich = " << readwhich <<
" rows = " << rows <<
" cols = " << cols <<
")" << std::endl;
322 double *likeAChecksum =
new double[cols];
323 double *likeAChecksum2 =
new double[cols];
324 for (
int i = 0;
i < cols;
i++) {
325 likeAChecksum[
i] = 0.;
326 likeAChecksum2[
i] = 0.;
329 for (
long row = 0; row <
rows; row++) {
330 double *residual =
new double[cols];
331 fread(residual,
sizeof(
double), cols, file);
334 for (
int i = 0;
i < cols;
i++) {
335 if (fabs(residual[
i]) > likeAChecksum[i]) likeAChecksum[
i] = fabs(residual[i]);
336 if (fabs(residual[i]) < likeAChecksum2[i]) likeAChecksum2[
i] = fabs(residual[i]);
340 double *readChecksum =
new double[cols];
341 double *readChecksum2 =
new double[cols];
342 fread(readChecksum,
sizeof(
double), cols, file);
343 fread(readChecksum2,
sizeof(
double), cols, file);
345 for (
int i = 0;
i < cols;
i++) {
346 if (fabs(likeAChecksum[
i] - readChecksum[
i]) > 1
e-10 || fabs(1./likeAChecksum2[i] - 1./readChecksum2[i]) > 1e10) {
347 throw cms::Exception(
"MuonResidualsFitter") <<
"temporary file is corrupted (which = " << which <<
" rows = " << rows <<
" likeAChecksum " << likeAChecksum[
i] <<
" != readChecksum " << readChecksum[
i] <<
" " <<
" likeAChecksum2 " << likeAChecksum2[
i] <<
" != readChecksum2 " << readChecksum2[
i] <<
")" << std::endl;
351 delete [] likeAChecksum;
352 delete [] likeAChecksum2;
357 if (which == 0) window = 2.*30.;
358 else if (which == 1) window = 2.*30.;
359 else if (which == 2) window = 2.*20.;
360 else if (which == 3) window = 2.*50.;
365 hist->Fill(multiplier * (*
r)[which]);
371 if (which == 0) window = 2.*30.;
372 else if (which == 1) window = 2.*30.;
373 else if (which == 2) window = 2.*20.;
374 else if (which == 3) window = 2.*50.;
379 double weight = 1./(*r)[whichredchi2];
380 if (TMath::Prob(1./weight*12, 12) < 0.99) {
381 hist->Fill(multiplier * (*
r)[which], weight);
double MuonResidualsFitter_logGaussPowerTails(double residual, double center, double sigma)
double MuonResidualsFitter_integrate_pureGaussian(double low, double high, double center, double sigma)
Double_t MuonResidualsFitter_powerLawTails_TF1(Double_t *xvec, Double_t *par)
virtual void inform(TMinuit *tMinuit)=0
const int MuonResidualsFitter_numgsbins
void write(FILE *file, int which=0)
long numResiduals() const
std::vector< double > m_error
void read(FILE *file, int which=0)
void fill(double *residual)
int residualsModel() const
Exp< T >::type exp(const T &t)
bool dofit(void(*fcn)(int &, double *, double &, double *, int), std::vector< int > &parNum, std::vector< std::string > &parName, std::vector< double > &start, std::vector< double > &step, std::vector< double > &low, std::vector< double > &high)
double MuonResidualsFitter_logROOTVoigt(double residual, double center, double sigma, double gamma)
double MuonResidualsFitter_compute_log_convolution(double toversigma, double gammaoversigma, double max=1000., double step=0.001, double power=4.)
std::vector< double > m_value
const T & max(const T &a, const T &b)
const int MuonResidualsFitter_numtsbins
double MuonResidualsFitter_logPowerLawTails(double residual, double center, double sigma, double gamma)
const double MuonResidualsFitter_gsbinsize
Double_t MuonResidualsFitter_GaussPowerTails_TF1(Double_t *xvec, Double_t *par)
const double MuonResidualsFitter_tsbinsize
static TMinuit * MuonResidualsFitter_TMinuit
Double_t MuonResidualsFitter_pureGaussian_TF1(Double_t *xvec, Double_t *par)
void plotsimple(std::string name, TFileDirectory *dir, int which, double multiplier)
Log< T >::type log(const T &t)
void fcn(int &, double *, double &, double *, int)
std::vector< double * >::const_iterator residuals_end() const
void plotweighted(std::string name, TFileDirectory *dir, int which, int whichredchi2, double multiplier)
bool MuonResidualsFitter_table_initialized
double MuonResidualsFitter_lookup_table[MuonResidualsFitter_numgsbins][MuonResidualsFitter_numtsbins]
T * make() const
make new ROOT object
Double_t MuonResidualsFitter_ROOTVoigt_TF1(Double_t *xvec, Double_t *par)
std::vector< double * >::const_iterator residuals_begin() const
double MuonResidualsFitter_logPureGaussian(double residual, double center, double sigma)
Power< A, B >::type pow(const A &a, const B &b)
const double par[8 *NPar][4]