CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MuonResidualsFitter.cc
Go to the documentation of this file.
2 #include <fstream>
3 #include "TMath.h"
4 
5 // all global variables begin with "MuonResidualsFitter_" to avoid
6 // namespace clashes (that is, they do what would ordinarily be done
7 // with a class structure, but Minuit requires them to be global)
8 
9 const double MuonResidualsFitter_gsbinsize = 0.01;
10 const double MuonResidualsFitter_tsbinsize = 0.1;
13 
16 
18 
19 // fit function
20 double MuonResidualsFitter_logPureGaussian(double residual, double center, double sigma) {
21  sigma = fabs(sigma);
22  return (-pow(residual - center, 2) / 2. / sigma / sigma) + log(1. / sqrt(2.*M_PI) / sigma);
23 }
24 
25 // TF1 interface version
26 Double_t MuonResidualsFitter_pureGaussian_TF1(Double_t *xvec, Double_t *par) {
27  return par[0] * exp(MuonResidualsFitter_logPureGaussian(xvec[0], par[1], par[2]));
28 }
29 
30 double MuonResidualsFitter_compute_log_convolution(double toversigma, double gammaoversigma, double max, double step, double power) {
31  if (gammaoversigma == 0.) return (-toversigma*toversigma/2.) - log(sqrt(2*M_PI));
32 
33  double sum = 0.;
34  double uplus = 0.;
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);
40 
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);
43 
44  sum += integrandplus * (uplus - uplus_last);
45  sum += 0.5 * fabs(integrandplus_last - integrandplus) * (uplus - uplus_last);
46 
47  sum += integrandminus * (uplus - uplus_last);
48  sum += 0.5 * fabs(integrandminus_last - integrandminus) * (uplus - uplus_last);
49 
50  integrandplus_last = integrandplus;
51  integrandminus_last = integrandminus;
52  }
53  return log(sum);
54 }
55 
56 // fit function
57 double MuonResidualsFitter_logPowerLawTails(double residual, double center, double sigma, double gamma) {
58  sigma = fabs(sigma);
59  double gammaoversigma = fabs(gamma / sigma);
60  double toversigma = fabs((residual - center) / sigma);
61 
62  int gsbin0 = int(floor(gammaoversigma / MuonResidualsFitter_gsbinsize));
63  int gsbin1 = int(ceil(gammaoversigma / MuonResidualsFitter_gsbinsize));
64  int tsbin0 = int(floor(toversigma / MuonResidualsFitter_tsbinsize));
65  int tsbin1 = int(ceil(toversigma / MuonResidualsFitter_tsbinsize));
66 
67  bool gsisbad = (gsbin0 >= MuonResidualsFitter_numgsbins || gsbin1 >= MuonResidualsFitter_numgsbins);
68  bool tsisbad = (tsbin0 >= MuonResidualsFitter_numtsbins || tsbin1 >= MuonResidualsFitter_numtsbins);
69 
70  if (gsisbad || tsisbad) {
71  return log(gammaoversigma/M_PI) - log(toversigma*toversigma + gammaoversigma*gammaoversigma) - log(sigma);
72  }
73  else {
74  double val00 = MuonResidualsFitter_lookup_table[gsbin0][tsbin0];
75  double val01 = MuonResidualsFitter_lookup_table[gsbin0][tsbin1];
76  double val10 = MuonResidualsFitter_lookup_table[gsbin1][tsbin0];
77  double val11 = MuonResidualsFitter_lookup_table[gsbin1][tsbin1];
78 
79  double val0 = val00 + ((toversigma / MuonResidualsFitter_tsbinsize) - tsbin0) * (val01 - val00);
80  double val1 = val10 + ((toversigma / MuonResidualsFitter_tsbinsize) - tsbin0) * (val11 - val10);
81 
82  double val = val0 + ((gammaoversigma / MuonResidualsFitter_gsbinsize) - gsbin0) * (val1 - val0);
83 
84  return val - log(sigma);
85  }
86 }
87 
88 // TF1 interface version
89 Double_t MuonResidualsFitter_powerLawTails_TF1(Double_t *xvec, Double_t *par) {
90  return par[0] * exp(MuonResidualsFitter_logPowerLawTails(xvec[0], par[1], par[2], par[3]));
91 }
92 
93 double MuonResidualsFitter_logROOTVoigt(double residual, double center, double sigma, double gamma) {
94  return log(TMath::Voigt(residual - center, fabs(sigma), fabs(gamma)*2.));
95 }
96 
97 Double_t MuonResidualsFitter_ROOTVoigt_TF1(Double_t *xvec, Double_t *par) {
98  return par[0] * TMath::Voigt(xvec[0] - par[1], fabs(par[2]), fabs(par[3])*2.);
99 }
100 
101 double MuonResidualsFitter_logGaussPowerTails(double residual, double center, double sigma) {
102  double x = residual-center;
103  double s = fabs(sigma);
104  double m = 2*s;
105  double a = pow(m,4)*exp(-2);
106  double n = sqrt(2*M_PI)*s*erf(sqrt(2))+2*a*pow(m,-3)/3;
107 
108  if (fabs(x)<m) return -x*x/(2*s*s) - log(n);
109  else return log(a) -4*log(fabs(x)) - log(n);
110 }
111 
112 Double_t MuonResidualsFitter_GaussPowerTails_TF1(Double_t *xvec, Double_t *par) {
113  return par[0] * exp(MuonResidualsFitter_logGaussPowerTails(xvec[0], par[1], par[2]));
114 }
115 
116 double MuonResidualsFitter_integrate_pureGaussian(double low, double high, double center, double sigma) {
117  return (erf((high + center) / sqrt(2.) / sigma) - erf((low + center) / sqrt(2.) / sigma)) * exp(0.5/sigma/sigma) / 2.;
118 }
119 
123 
124  std::ifstream convolution_table("convolution_table.txt");
125  if (convolution_table.is_open()) {
126  int numgsbins = 0;
127  int numtsbins = 0;
128  double tsbinsize = 0.;
129  double gsbinsize = 0.;
130 
131  convolution_table >> numgsbins >> numtsbins >> tsbinsize >> gsbinsize;
132  if (numgsbins != MuonResidualsFitter_numgsbins || numtsbins != MuonResidualsFitter_numtsbins ||
133  tsbinsize != MuonResidualsFitter_tsbinsize || gsbinsize != MuonResidualsFitter_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;
135  }
136 
137  for (int gsbin = 0; gsbin < MuonResidualsFitter_numgsbins; gsbin++) {
138  for (int tsbin = 0; tsbin < MuonResidualsFitter_numtsbins; tsbin++) {
139  int read_gsbin = 0;
140  int read_tsbin = 0;
141  double value = 0.;
142 
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;
146  }
147 
148  MuonResidualsFitter_lookup_table[gsbin][tsbin] = value;
149  }
150  }
151 
152  convolution_table.close();
153  }
154 
155  else {
156  std::ofstream convolution_table2("convolution_table.txt");
157 
158  if (!convolution_table2.is_open()) {
159  throw cms::Exception("MuonResidualsFitter") << "Couldn't write to file convolution_table.txt" << std::endl;
160  }
161 
163 
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;
166 
167  for (int gsbin = 0; gsbin < MuonResidualsFitter_numgsbins; gsbin++) {
168  double gammaoversigma = double(gsbin) * MuonResidualsFitter_gsbinsize;
169 
170  std::cout << " gsbin " << gsbin << "/" << MuonResidualsFitter_numgsbins << std::endl;
171 
172  for (int tsbin = 0; tsbin < MuonResidualsFitter_numtsbins; tsbin++) {
173  double toversigma = double(tsbin) * MuonResidualsFitter_tsbinsize;
174 
175  // 1e-6 errors (out of a value of ~0.01) with max=100, step=0.001, power=4 (max=1000 does a little better with the tails)
176  MuonResidualsFitter_lookup_table[gsbin][tsbin] = MuonResidualsFitter_compute_log_convolution(toversigma, gammaoversigma);
177 
178  // <10% errors with max=20, step=0.005, power=4 (faster computation for testing)
179  // MuonResidualsFitter_lookup_table[gsbin][tsbin] = MuonResidualsFitter_compute_log_convolution(toversigma, gammaoversigma, 100., 0.005, 4.);
180 
181  convolution_table2 << gsbin << " " << tsbin << " " << MuonResidualsFitter_lookup_table[gsbin][tsbin] << std::endl;
182  }
183  }
184 
185  convolution_table2.close();
186 
187  edm::LogWarning("MuonResidualsFitter") << "Initialization done!" << std::endl;
188  std::cout << "Initialization done!" << std::endl;
189  }
190 }
191 
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) {
194 
195  MuonResidualsFitter_TMinuit = new TMinuit(npar());
197  MuonResidualsFitter_TMinuit->SetObjectFit(fitinfo);
200 
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();
207 
208  for (; iNum != parNum.end(); ++iNum, ++iName, ++istart, ++istep, ++ilow, ++ihigh) {
209  MuonResidualsFitter_TMinuit->DefineParameter(*iNum, iName->c_str(), *istart, *istep, *ilow, *ihigh);
210  if (fixed(*iNum)) MuonResidualsFitter_TMinuit->FixParameter(*iNum);
211  }
212 
213  double arglist[10];
214  int ierflg;
215  int smierflg; //second MIGRAD ierflg
216 
217  // chi^2 errors should be 1.0, log-likelihood should be 0.5
218  for (int i = 0; i < 10; i++) arglist[i] = 0.;
219  arglist[0] = 0.5;
220  ierflg = 0;
221  smierflg = 0;
222  MuonResidualsFitter_TMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
223  if (ierflg != 0) { delete MuonResidualsFitter_TMinuit; delete fitinfo; return false; }
224 
225  // set strategy = 2 (more refined fits)
226  for (int i = 0; i < 10; i++) arglist[i] = 0.;
227  arglist[0] = m_strategy;
228  ierflg = 0;
229  MuonResidualsFitter_TMinuit->mnexcm("SET STR", arglist, 1, ierflg);
230  if (ierflg != 0) { delete MuonResidualsFitter_TMinuit; delete fitinfo; return false; }
231 
232  bool try_again = false;
233 
234  // minimize
235  for (int i = 0; i < 10; i++) arglist[i] = 0.;
236  arglist[0] = 50000;
237  ierflg = 0;
238  MuonResidualsFitter_TMinuit->mnexcm("MIGRAD", arglist, 1, ierflg);
239  if (ierflg != 0) try_again = true;
240 
241  // just once more, if needed (using the final Minuit parameters from the failed fit; often works)
242  if (try_again) {
243  // minimize
244  for (int i = 0; i < 10; i++) arglist[i] = 0.;
245  arglist[0] = 50000;
246  MuonResidualsFitter_TMinuit->mnexcm("MIGRAD", arglist, 1, smierflg);
247  }
248 
249  Double_t fmin, fedm, errdef;
250  Int_t npari, nparx, istat;
251  MuonResidualsFitter_TMinuit->mnstat(fmin, fedm, errdef, npari, nparx, istat);
252 
253  if (istat != 3) {
254  for (int i = 0; i < 10; i++) arglist[i] = 0.;
255  ierflg = 0;
256  MuonResidualsFitter_TMinuit->mnexcm("HESSE", arglist, 0, ierflg);
257  }
258 
259  // read-out the results
260  m_loglikelihood = -fmin;
261 
262  m_value.clear();
263  m_error.clear();
264  for (int i = 0; i < npar(); i++) {
265  double v, e;
266  MuonResidualsFitter_TMinuit->GetParameter(i, v, e);
267  m_value.push_back(v);
268  m_error.push_back(e);
269  }
270 
272  delete fitinfo;
273  if (smierflg != 0) return false;
274  return true;
275 }
276 
277 void MuonResidualsFitter::write(FILE *file, int which) {
278  long rows = numResiduals();
279  int cols = ndata();
280  int whichcopy = which;
281 
282  fwrite(&rows, sizeof(long), 1, file);
283  fwrite(&cols, sizeof(int), 1, file);
284  fwrite(&whichcopy, sizeof(int), 1, file);
285 
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.;
291  }
292 
293  for (std::vector<double*>::const_iterator residual = residuals_begin(); residual != residuals_end(); ++residual) {
294  fwrite((*residual), sizeof(double), cols, file);
295 
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]);
299  }
300  } // end loop over residuals
301 
302  // the idea is that mal-formed doubles are likely to be huge values (or tiny values)
303  // because the exponent gets screwed up; we want to check for that
304  fwrite(likeAChecksum, sizeof(double), cols, file);
305  fwrite(likeAChecksum2, sizeof(double), cols, file);
306 
307  delete [] likeAChecksum;
308  delete [] likeAChecksum2;
309 }
310 
311 void MuonResidualsFitter::read(FILE *file, int which) {
312  long rows = -100;
313  int cols = -100;
314  int readwhich = -100;
315 
316  fread(&rows, sizeof(long), 1, file);
317  fread(&cols, sizeof(int), 1, file);
318  fread(&readwhich, sizeof(int), 1, file);
319 
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;
321 
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.;
327  }
328 
329  for (long row = 0; row < rows; row++) {
330  double *residual = new double[cols];
331  fread(residual, sizeof(double), cols, file);
332  fill(residual);
333 
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]);
337  }
338  } // end loop over records in file
339 
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);
344 
345  for (int i = 0; i < cols; i++) {
346  if (fabs(likeAChecksum[i] - readChecksum[i]) > 1e-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;
348  }
349  }
350 
351  delete [] likeAChecksum;
352  delete [] likeAChecksum2;
353 }
354 
355 void MuonResidualsFitter::plotsimple(std::string name, TFileDirectory *dir, int which, double multiplier) {
356  double window = 100.;
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.;
361 
362  TH1F *hist = dir->make<TH1F>(name.c_str(), "", 200, -window, window);
363 
364  for (std::vector<double*>::const_iterator r = residuals_begin(); r != residuals_end(); ++r) {
365  hist->Fill(multiplier * (*r)[which]);
366  }
367 }
368 
369 void MuonResidualsFitter::plotweighted(std::string name, TFileDirectory *dir, int which, int whichredchi2, double multiplier) {
370  double window = 100.;
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.;
375 
376  TH1F *hist = dir->make<TH1F>(name.c_str(), "", 200, -window, window);
377 
378  for (std::vector<double*>::const_iterator r = residuals_begin(); r != residuals_end(); ++r) {
379  double weight = 1./(*r)[whichredchi2];
380  if (TMath::Prob(1./weight*12, 12) < 0.99) {
381  hist->Fill(multiplier * (*r)[which], weight);
382  }
383  }
384 }
tuple arglist
Definition: fitWZ.py:38
double MuonResidualsFitter_logGaussPowerTails(double residual, double center, double sigma)
int i
Definition: DBlmapReader.cc:9
def window
Definition: svgfig.py:642
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
virtual int npar()=0
list file
Definition: dbtoweb.py:253
std::vector< double > m_error
void read(FILE *file, int which=0)
void fill(double *residual)
double value(int parNum)
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
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)
bool fixed(int parNum)
T sqrt(T t)
Definition: SSEVec.h:28
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)
virtual int ndata()=0
#define M_PI
Definition: BFit3D.cc:3
Log< T >::type log(const T &t)
Definition: Log.h:22
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)
double a
Definition: hdecay.h:121
bool MuonResidualsFitter_table_initialized
double MuonResidualsFitter_lookup_table[MuonResidualsFitter_numgsbins][MuonResidualsFitter_numtsbins]
T * make() const
make new ROOT object
tuple cout
Definition: gather_cfg.py:41
string s
Definition: asciidump.py:422
dbl *** dir
Definition: mlp_gen.cc:35
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)
mathSSE::Vec4< T > v
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
const double par[8 *NPar][4]