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