CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
fitUtilities.h
Go to the documentation of this file.
1 #ifndef DQM_SiStripHistoricInfoClient_fitUtilities_H
2 #define DQM_SiStripHistoricInfoClient_fitUtilities_H
3 
4 #include <vector>
5 #include <cstring>
6 #include <iostream>
7 #include <sstream>
8 #include <string>
9 #include "TH1.h"
10 #include "TF1.h"
11 #include "TMath.h"
12 #include "DQMServices/Core/interface/MonitorElement.h"
13 #include "FWCore/MessageLogger/interface/MessageLogger.h"
14 
15 
24 double langaufun(double *x, double *par);
25 int32_t langaupro(double *params, double &maxx, double &FWHM);
26 double Gauss(double *x, double *par);
27 
28 class fitUtilities{
29 
30  public:
31 
32  fitUtilities();
33  ~fitUtilities();
34 
35  void init();
36  double doLanGaussFit(MonitorElement* ME){return doLanGaussFit(ME->getTH1F());}
37  double doLanGaussFit(TH1F*);
38 
39  double doGaussFit(MonitorElement* ME){return doGaussFit(ME->getTH1F());}
40  double doGaussFit(TH1F*);
41 
42  double getLanGaussPar(std::string s) ;
43  double getLanGaussParErr(std::string s);
44  double getLanGaussConv(std::string s) ;
45 
46  double getGaussPar(std::string s) ;
47  double getGaussParErr(std::string s);
48 
49  double getFitChi() {return chi2GausS;}
50  int getFitnDof() {return nDofGausS;}
51 
52  private:
53 
54  double pLanGausS[4], epLanGausS[4];
55  double pGausS[3], epGausS[3];
56  double pLanConv[2];
57  double chi2GausS;
58  int32_t nDofGausS;
59  TF1 *langausFit;
60  TF1 *gausFit;
61 };
62 
63 #endif // DQM_SiStripHistoricInfoClient_fitUtilities_H
64 
65 //-----------------------------------------------------------------------------------------------
66 double langaufun(double *x, double *par){
67  //Fit parameters:
68  //par[0]=Width (scale) parameter of Landau density
69  //par[1]=Most Probable (MP, location) parameter of Landau density
70  //par[2]=Total area (integral -inf to inf, normalization constant)
71  //par[3]=Width (sigma) of convoluted Gaussian function
72  //
73  //In the Landau distribution (represented by the CERNLIB approximation),
74  //the maximum is located at x=-0.22278298 with the location parameter=0.
75  //This shift is corrected within this function, so that the actual
76  //maximum is identical to the MP parameter.
77 
78  // Numeric constants
79  double invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
80  double mpshift = -0.22278298; // Landau maximum location
81 
82  // Control constants
83  double np = 100.0; // number of convolution steps
84  double sc = 5.0; // convolution extends to +-sc Gaussian sigmas
85 
86  // Variables
87  double xx;
88  double mpc;
89  double fland;
90  double sum = 0.0;
91  double xlow,xupp;
92  double step;
93  double i;
94 
95 
96  // MP shift correction
97  mpc = par[1] - mpshift * par[0];
98 
99  // Range of convolution integral
100  xlow = x[0] - sc * par[3];
101  xupp = x[0] + sc * par[3];
102 
103  step = (xupp-xlow) / np;
104 
105  // Landau Distribution Production
106  for(i=1.0; i<=np/2; i++) {
107  xx = xlow + (i-.5) * step;
108  fland = TMath::Landau(xx,mpc,par[0]) / par[0];
109  sum += fland * TMath::Gaus(x[0],xx,par[3]);
110 
111  xx = xupp - (i-.5) * step;
112  fland = TMath::Landau(xx,mpc,par[0]) / par[0];
113  sum += fland * TMath::Gaus(x[0],xx,par[3]);
114  }
115 
116  return (par[2] * step * sum * invsq2pi / par[3]);
117 }
118 
119 //-----------------------------------------------------------------------------------------------
120 int32_t langaupro(double *params, double &maxx, double &FWHM) {
121  edm::LogInfo("fitUtility") << "inside langaupro " << std::endl;
122  // Seaches for the location (x value) at the maximum of the
123  // Landau and its full width at half-maximum.
124  //
125  // The search is probably not very efficient, but it's a first try.
126 
127  double p,x,fy,fxr,fxl;
128  double step;
129  double l,lold,dl;
130  int32_t i = 0;
131  const int32_t MAXCALLS = 10000;
132  const double dlStop = 1e-3; // relative change < .001
133 
134  // Search for maximum
135  p = params[1] - 0.1 * params[0];
136  step = 0.05 * params[0];
137  lold = -2.0;
138  l = -1.0;
139 
140  dl = (l-lold)/lold; // FIXME catch divide by zero
141  while ( (TMath::Abs(dl)>dlStop ) && (i < MAXCALLS) ) {
142  i++;
143 
144  lold = l;
145  x = p + step;
146  l = langaufun(&x,params);
147  dl = (l-lold)/lold; // FIXME catch divide by zero
148 
149  if (l < lold)
150  step = -step/10;
151 
152  p += step;
153  }
154 
155  if (i == MAXCALLS)
156  return (-1);
157 
158  maxx = x;
159 
160  fy = l/2;
161 
162 
163  // Search for right x location of fy
164  p = maxx + params[0];
165  step = params[0];
166  lold = -2.0;
167  l = -1e300;
168  i = 0;
169 
170  dl = (l-lold)/lold; // FIXME catch divide by zero
171  while ( ( TMath::Abs(dl)>dlStop ) && (i < MAXCALLS) ) {
172  i++;
173 
174  lold = l;
175  x = p + step;
176  l = TMath::Abs(langaufun(&x,params) - fy);
177  dl = (l-lold)/lold; // FIXME catch divide by zero
178 
179  if (l > lold)
180  step = -step/10;
181 
182  p += step;
183  }
184 
185  if (i == MAXCALLS)
186  return (-2);
187 
188  fxr = x;
189 
190 
191  // Search for left x location of fy
192  p = maxx - 0.5 * params[0];
193  step = -params[0];
194  lold = -2.0;
195  l = -1e300;
196  i = 0;
197 
198  dl = (l-lold)/lold; // FIXME catch divide by zero
199  while ( ( TMath::Abs(dl)>dlStop ) && (i < MAXCALLS) ) {
200  i++;
201 
202  lold = l;
203  x = p + step;
204  l = TMath::Abs(langaufun(&x,params) - fy);
205  dl = (l-lold)/lold; // FIXME catch divide by zero
206 
207  if (l > lold)
208  step = -step/10;
209 
210  p += step;
211  }
212 
213  if (i == MAXCALLS)
214  return (-3);
215 
216 
217  fxl = x;
218 
219  FWHM = fxr - fxl;
220  return (0);
221 }
222 
223 
224 //-----------------------------------------------------------------------------------------------
225 double Gauss(double *x, double *par)
226 {
227  // The noise function: a gaussian
228 
229  double arg = 0;
230  if (par[2]) arg = (x[0] - par[1])/par[2];
231 
232  double noise = par[0]*TMath::Exp(-0.5*arg*arg);
233  return noise;
234 }
235 
236 
237 
238 //-----------------------------------------------------------------------------------------------
239 fitUtilities::fitUtilities():langausFit(0),gausFit(0){
240  init();
241 }
242 
243 //-----------------------------------------------------------------------------------------------
244 void fitUtilities::init(){
245  pLanGausS[0]=0 ; pLanGausS[1]=0; pLanGausS[2]=0; pLanGausS[3]=0;
246  epLanGausS[0]=0; epLanGausS[1]=0; epLanGausS[2]=0; epLanGausS[3]=0;
247  pGausS[0]=0 ; pGausS[1]=0; pGausS[2]=0;
248  epGausS[0]=0; epGausS[1]=0; epGausS[2]=0;
249  pLanConv[0]=0; pLanConv[1]=0;
250 
251  //if ( langausFit!=0 ) delete langausFit;
252  //if ( gausFit!=0 ) delete gausFit;
253 
254  chi2GausS = -9999.;
255  nDofGausS = -9999;
256 }
257 
258 //-----------------------------------------------------------------------------------------------
259 fitUtilities::~fitUtilities(){
260  if ( langausFit!=0 ) delete langausFit;
261  if ( gausFit!=0 ) delete gausFit;
262 }
263 
264 
265 
266 //-----------------------------------------------------------------------------------------------
267 double fitUtilities::doLanGaussFit(TH1F* htoFit){
268  init();
269 
270  if (htoFit->GetEntries()!=0) {
271  edm::LogInfo("fitUtility")<<"Fitting "<< htoFit->GetTitle() <<std::endl;
272  // Setting fit range and start values
273  double fr[2];
274  double sv[4], pllo[4], plhi[4];
275  fr[0]=0.5*htoFit->GetMean();
276  fr[1]=3.0*htoFit->GetMean();
277 
278  // (EM) parameters setting good for signal only
279  int32_t imax=htoFit->GetMaximumBin();
280  double xmax=htoFit->GetBinCenter(imax);
281  double ymax=htoFit->GetBinContent(imax);
282  int32_t i[2];
283  int32_t iArea[2];
284 
285  i[0]=htoFit->GetXaxis()->FindBin(fr[0]);
286  i[1]=htoFit->GetXaxis()->FindBin(fr[1]);
287 
288  iArea[0]=htoFit->GetXaxis()->FindBin(fr[0]);
289  iArea[1]=htoFit->GetXaxis()->FindBin(fr[1]);
290  double AreaFWHM=htoFit->Integral(iArea[0],iArea[1],"width");
291 
292  sv[1]=xmax;
293  sv[2]=htoFit->Integral(i[0],i[1],"width");
294  sv[3]=AreaFWHM/(4*ymax);
295  sv[0]=sv[3];
296 
297  plhi[0]=25.0; plhi[1]=200.0; plhi[2]=1000000.0; plhi[3]=50.0;
298  pllo[0]=1.5 ; pllo[1]=10.0 ; pllo[2]=1.0 ; pllo[3]= 1.0;
299 
300  Char_t FunName[100];
301  sprintf(FunName,"FitfcnLG_%s",htoFit->GetName());
302 
303  langausFit = new TF1(FunName,langaufun,fr[0],fr[1],4);
304  langausFit->SetParameters(sv);
305  langausFit->SetParNames("Width","MP","Area","GSigma");
306 
307  for (int32_t i=0; i<4; i++) {
308  langausFit->SetParLimits(i,pllo[i],plhi[i]);
309  }
310 
311  try{
312  htoFit->Fit(langausFit,"R0"); // "R" fit in a range,"0" quiet fit
313 
314  langausFit->SetRange(fr[0],fr[1]);
315  langausFit->GetParameters(pLanGausS);
316  std::memcpy((void*) epLanGausS, (void*) langausFit->GetParErrors(), 4*sizeof(double));
317 
318  chi2GausS =langausFit->GetChisquare(); // obtain chi^2
319  nDofGausS = langausFit->GetNDF(); // obtain ndf
320 
321  double sPeak, sFWHM;
322  langaupro(pLanGausS,sPeak,sFWHM);
323  pLanConv[0]=sPeak;
324  pLanConv[1]=sFWHM;
325  edm::LogInfo("fitUtility") << "langaupro: max " << sPeak << std::endl;
326  edm::LogInfo("fitUtility") << "langaupro: FWHM " << sFWHM << std::endl;
327  }
328  catch(...){
329  edm::LogError("fitUtility") << "problem in fitting " << htoFit->GetTitle() << " \n\tDefault values of the parameters will be used";
330  pLanGausS[0]=-9999; pLanGausS[1]=-9999; pLanGausS[2]=-9999; pLanGausS[3]=-9999;
331  epLanGausS[0]=-9999; epLanGausS[1]=-9999; epLanGausS[2]=-9999; epLanGausS[3]=-9999;
332  pLanConv[0]=-9999; pLanConv[1]=-9999;
333  chi2GausS=-9999; nDofGausS=-9999;
334  }
335  }
336  else {
337  pLanGausS[0]=-9999; pLanGausS[1]=-9999; pLanGausS[2]=-9999; pLanGausS[3]=-9999;
338  epLanGausS[0]=-9999; epLanGausS[1]=-9999; epLanGausS[2]=-9999; epLanGausS[3]=-9999;
339  pLanConv[0]=-9999; pLanConv[1]=-9999;
340  chi2GausS=-9999; nDofGausS=-9999;
341  }
342 
343  return htoFit->GetEntries();
344 
345 }
346 
347 
348 //-----------------------------------------------------------------------------------------------
349 double fitUtilities::doGaussFit(TH1F* htoFit){
350  init();
351  if (htoFit->GetEntries()!=0) {
352 
353  // Setting fit range and start values
354  double fr[2];
355  double sv[3], pllo[3], plhi[3];
356  fr[0]=htoFit->GetMean()-5*htoFit->GetRMS();
357  fr[1]=htoFit->GetMean()+5*htoFit->GetRMS();
358 
359  int32_t imax=htoFit->GetMaximumBin();
360  double xmax=htoFit->GetBinCenter(imax);
361  double ymax=htoFit->GetBinContent(imax);
362  int32_t i[2];
363  int32_t iArea[2];
364 
365  i[0]=htoFit->GetXaxis()->FindBin(fr[0]);
366  i[1]=htoFit->GetXaxis()->FindBin(fr[1]);
367 
368  iArea[0]=htoFit->GetXaxis()->FindBin(fr[0]);
369  iArea[1]=htoFit->GetXaxis()->FindBin(fr[1]);
370  double AreaFWHM=htoFit->Integral(iArea[0],iArea[1],"width");
371 
372  sv[2]=AreaFWHM/(4*ymax);
373  sv[1]=xmax;
374  sv[0]=htoFit->Integral(i[0],i[1],"width");
375 
376  plhi[0]=1000000.0; plhi[1]=10.0; plhi[2]=10.;
377  pllo[0]=1.5 ; pllo[1]=0.1; pllo[2]=0.3;
378  Char_t FunName[100];
379  sprintf(FunName,"FitfcnLG_%s",htoFit->GetName());
380  gausFit = new TF1(FunName,Gauss,fr[0],fr[1],3);
381  gausFit->SetParameters(sv);
382  gausFit->SetParNames("Constant","GaussPeak","Sigma");
383 
384  for (int32_t i=0; i<3; i++) {
385  gausFit->SetParLimits(i,pllo[i],plhi[i]);
386  }
387 
388  try{
389  htoFit->Fit(gausFit,"R0");
390 
391  gausFit->SetRange(fr[0],fr[1]);
392  gausFit->GetParameters(pGausS);
393  std::memcpy((void*) epGausS, (void*) gausFit->GetParErrors(), 3*sizeof(double));
394 
395  chi2GausS =langausFit->GetChisquare(); // obtain chi^2
396  nDofGausS = langausFit->GetNDF();// obtain ndf
397  }
398  catch(...){
399  edm::LogError("fitUtility") << "problem in fitting " << htoFit->GetTitle() << " \n\tDefault values of the parameters will be used";
400  pGausS[0]=-9999; pGausS[1]=-9999; pGausS[2]=-9999;
401  epGausS[0]=-9999; epGausS[1]=-9999; epGausS[2]=-9999;
402  chi2GausS=-9999; nDofGausS=-9999;
403  }
404 
405  }
406  else {
407  pGausS[0]=-9999; pGausS[1]=-9999; pGausS[2]=-9999;
408  epGausS[0]=-9999; epGausS[1]=-9999; epGausS[2]=-9999;
409  chi2GausS=-9999; nDofGausS=-9999;
410  }
411 
412  return htoFit->GetEntries();
413 }
414 
415 //-----------------------------------------------------------------------------------------------
416 double fitUtilities::getLanGaussPar(std::string s){
417  if(s=="landau_width")
418  return pLanGausS[0];
419  else if(s=="mpv")
420  return pLanGausS[1];
421  else if(s=="area")
422  return pLanGausS[2];
423  else if(s=="gauss_sigma")
424  return pLanGausS[3];
425  else
426  return -99999;
427 }
428 
429 double fitUtilities::getLanGaussParErr(std::string s){
430  if(s=="landau_width")
431  return epLanGausS[0];
432  else if(s=="mpv")
433  return epLanGausS[1];
434  else if(s=="area")
435  return epLanGausS[2];
436  else if(s=="gauss_sigma")
437  return epLanGausS[3];
438  else
439  return -99999;
440 }
441 
442 double fitUtilities::getLanGaussConv(std::string s) {
443  if(s=="mpv")
444  return pLanConv[0];
445  else if(s=="fwhm")
446  return pLanConv[1];
447  else
448  return -99999;
449 }
450 
451 double fitUtilities::getGaussPar(std::string s) {
452  if(s=="area")
453  return pGausS[0];
454  else if(s=="mean")
455  return pGausS[1];
456  else if(s=="sigma")
457  return pGausS[2];
458  else
459  return -99999;
460 }
461 
462 double fitUtilities::getGaussParErr(std::string s) {
463  if(s=="area")
464  return epGausS[0];
465  else if(s=="mean")
466  return epGausS[1];
467  else if(s=="sigma")
468  return epGausS[2];
469  else
470  return -99999;
471 }
i
int i
Definition: DBlmapReader.cc:9
fitUtilities::pGausS
double pGausS[3]
Definition: fitUtilities.h:55
fitUtilities::getFitnDof
int getFitnDof()
Definition: fitUtilities.h:50
fitUtilities::langausFit
TF1 * langausFit
Definition: fitUtilities.h:59
fitUtilities::getLanGaussParErr
double getLanGaussParErr(std::string s)
Definition: fitUtilities.h:429
fitUtilities::pLanGausS
double pLanGausS[4]
Definition: fitUtilities.h:54
fitUtilities::getFitChi
double getFitChi()
Definition: fitUtilities.h:49
fitUtilities::epLanGausS
double epLanGausS[4]
Definition: fitUtilities.h:54
SiStripMonitorClusterAlca_cfi.ymax
tuple ymax
Definition: SiStripMonitorClusterAlca_cfi.py:37
funct::arg
A arg
Definition: Factorize.h:36
fitUtilities::getGaussPar
double getGaussPar(std::string s)
Definition: fitUtilities.h:451
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
ME
Definition: ME.h:11
fitUtilities::doGaussFit
double doGaussFit(MonitorElement *ME)
Definition: fitUtilities.h:39
np
int np
Definition: AMPTWrapper.h:33
Gauss
double Gauss(double *x, double *par)
Definition: fitUtilities.h:225
fitUtilities::getLanGaussConv
double getLanGaussConv(std::string s)
Definition: fitUtilities.h:442
fitUtilities::getLanGaussPar
double getLanGaussPar(std::string s)
Definition: fitUtilities.h:416
langaupro
int32_t langaupro(double *params, double &maxx, double &FWHM)
Definition: fitUtilities.h:120
fitUtilities::getGaussParErr
double getGaussParErr(std::string s)
Definition: fitUtilities.h:462
fitUtilities::pLanConv
double pLanConv[2]
Definition: fitUtilities.h:56
MonitorElement::getTH1F
TH1F * getTH1F(void) const
Definition: MonitorElement.cc:1527
fitUtilities::doLanGaussFit
double doLanGaussFit(MonitorElement *ME)
Definition: fitUtilities.h:36
SiStripMonitorClusterAlca_cfi.xmax
tuple xmax
Definition: SiStripMonitorClusterAlca_cfi.py:35
fitUtilities::chi2GausS
double chi2GausS
Definition: fitUtilities.h:57
relval_parameters_module.step
string step
Definition: relval_parameters_module.py:36
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
fitUtilities::nDofGausS
int32_t nDofGausS
Definition: fitUtilities.h:58
x
Definition: DDAxes.h:10
langaufun
Double_t langaufun(Double_t *x, Double_t *par)
Definition: HOCalibAnalyzer.cc:116
fitUtilities::epGausS
double epGausS[3]
Definition: fitUtilities.h:55
Generated for CMSSW Reference Manual by doxygen 1.8.5