ATLAS Offline Software
Resplot.cxx
Go to the documentation of this file.
1 
20 #include <iostream>
21 #include <sstream>
22 #include <cmath>
23 
24 
25 #include "Directory.h"
26 #include "tagname.h"
27 #include "Resplot.h"
28 
29 #include "TCanvas.h"
30 #include "TStyle.h"
31 #include "TMath.h"
32 
33 #include "generate.h"
34 #include "rmsFrac.h"
35 
36 ClassImp(Resplot)
37 
38 
39 bool Resplot::s_mAddDirectoryStatus = true;
40 bool Resplot::s_interpolate_flag = true;
41 bool Resplot::s_nofit = false;
42 
44 bool Resplot::s_oldrms95 = true;
45 bool Resplot::s_scalerms95 = true;
46 
47 Resplot::ERROR Resplot::s_ErrorSet = Resplot::OK;
48 
49 
50 void binwidth(TH1D* h) {
51 
52  for (int i=1 ; i<=h->GetNbinsX() ; i++ ) {
53  double d = h->GetBinLowEdge(i+1) - h->GetBinLowEdge(i);
54 
55  h->SetBinContent(i, h->GetBinContent(i)/d );
56  h->SetBinError(i, h->GetBinError(i)/d );
57  }
58 
59 }
60 
61 
62 // set the errors on bins with zero entries to be one so they
63 // are not ignored in the fit
64 // USE WITH CARE - root screws up the mean and rms error calculation
65 // if these are used
66 
67 void ZeroErrors(TH1D* h) {
68  for (int i=1 ; i<=h->GetNbinsX() ; i++ ) if (h->GetBinContent(i)==0) h->SetBinError(i,1);
69 }
70 
71 // and put them back to zero again afterwards
72 
73 void unZeroErrors(TH1D* h) {
74  for (int i=1 ; i<=h->GetNbinsX() ; i++ ) if (h->GetBinContent(i)==0) h->SetBinError(i,0);
75 }
76 
77 
78 
79 double getWeights(TH1D* h) {
80  double sum = 0;
81  double sume2 = 0;
82  for ( int i=0 ; i<h->GetNbinsX()+1 ; i++ ) {
83  sum += h->GetBinContent(i);
84  sume2 += h->GetBinError(i)*h->GetBinError(i);
85  }
86  if ( sume2!=0 ) return sum*sum/sume2;
87  return 0;
88 }
89 
90 
91 double getWeights(TH2D* h) {
92 
93  double sum = 0;
94  double sume2 = 0;
95 
96  for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) {
97  for ( int j=1 ; j<=h->GetNbinsY() ; j++ ) {
98  sum += h->GetBinContent(i,j);
99  sume2 += h->GetBinError(i,j)*h->GetBinError(i,j);
100  }
101  }
102 
103  if ( sume2!=0 ) return sum*sum/sume2;
104 
105  return 0;
106 }
107 
108 
109 
110 
111 const std::string Resplot::s_rversion = "resplot-v29";
112 
113 void GetStats( TH1D* h, std::vector<double>& stats ) {
114 
115  stats.clear();
116  stats.resize(4);
117 
118  TAxis& fXaxis = *h->GetXaxis();
119 
120  double s = 0;
121  double sx = 0;
122  double sx2 = 0;
123  double sx4 = 0;
124 
125 
126  // if ( fXaxis.TestBit(TAxis::kAxisRange) ) {
127  // for (bin=0;bin<4;bin++) stats[bin] = 0;
128 
129  Int_t firstBinX = fXaxis.GetFirst();
130  Int_t lastBinX = fXaxis.GetLast();
131  // include underflow/overflow if TH1::StatOverflows(kTRUE) in case no range is set on the axis
132  // if (fgStatOverflows && !fXaxis.TestBit(TAxis::kAxisRange)) {
133  // if (firstBinX == 1) firstBinX = 0;
134  // if (lastBinX == fXaxis.GetNbins() ) lastBinX += 1;
135  // }
136  for (int binx = firstBinX; binx <= lastBinX; binx++) {
137  double x = fXaxis.GetBinCenter(binx);
138  double w = TMath::Abs(h->GetBinContent(binx));
139  // err = TMath::Abs(GetBinError(binx));
140  s += w;
141  // stats[1] += err*err;
142  sx += w*x;
143  sx2 += w*x*x;
144  sx4 += w*x*x*x*x;
145  }
146 
147 
148  double mu = sx/s;
149 
150  double v = sx2/s-mu*mu;
151 
152  double rms = std::sqrt(v);
153  double rmserror = std::sqrt(0.25*(sx4/v-v*s))/s;
154  // double rmserror1 = std::sqrt(0.4*(sx4/s-v*v)/s); << is this the correct one? it seems about 2* too large
155 
156  stats[0] = mu;
157  stats[1] = std::sqrt(v/s);
158 
159  stats[2] = rms;
160  stats[3] = rmserror;
161 
162  // double duff = std::sqrt(0.5*v/s);
163 
164 #if 0
165  std::cout << "GetStats() "
166  // << "\tmean " << stats[0] << " +- " << stats[1]
167  << "\trms " << stats[2] << " +- " << stats[3] << "\t(" << rmserror1 << ")" << "\t root " << duff << std::endl;
168 #endif
169 
170 }
171 
172 
173 
174 Resplot::~Resplot() {
175 
176  // std::cout << "Resplot::~Resplot() " << Name() << std::endl;
177 
178  delskip( m_Nentries );
179  delskip( m_h2d );
180  delskip( m_h1d );
181  delskip( m_mean );
182  delskip( m_sigma );
183  delskip( m_chi2 );
184 
185  // std::cout << "Resplot::~Resplot() deleteing slices" << Name() << std::endl;
186 
187  for ( unsigned i=0 ; i<m_slices.size() ; i++ ) delskip(m_slices[i]);
188 
189  // TH2D::SetDirectory(0);
190 }
191 
192 
193 
194 void Resplot::Initialise(const std::string& name,
195  int n1, double a1, double b1,
196  int n2, double a2, double b2) {
197 
198  SetPrimary(n1, a1, b1);
199  SetSecondary(n2, a2, b2);
200 
201  m_Set = true;
202  // TH2D::SetDirectory(0);
203 
204  m_name = name;
205  const std::string name_2d = "2d";
206  const std::string name_1d = "1d";
207 
208  // std::cout << "Resplot::Initialise() Name " << Name() << std::endl;
209 
210  m_h2d = new TH2D( name_2d.c_str(), name_2d.c_str(), n1, a1, b1, n2, a2, b2);
211  m_h1d = new TH1D( name_1d.c_str(), name_1d.c_str(), n2, a2, b2);
212  m_h1d->Sumw2();
213  skip(m_h2d);
214  skip(m_h1d);
215 
216  const std::string name_mean = "mean";
217  const std::string name_sigma = "sigma";
218  const std::string name_chi2 = "chi2";
219  // const std::string name_entries = "significance";
220  const std::string name_entries = "fractional uncertainty";
221 
222  m_Nentries = new TH1D( name_entries.c_str(), name_entries.c_str(), n1, a1, b1 );
223  skip(m_Nentries);
224 
225  m_mean = new TH1D( name_mean.c_str(), name_mean.c_str(), n1, a1, b1 );
226  m_sigma = new TH1D( name_sigma.c_str(), name_sigma.c_str(), n1, a1, b1 );
227  m_chi2 = new TH1D( name_chi2.c_str(), name_chi2.c_str(), n1, a1, b1 );
228  skip(m_mean);
229  skip(m_sigma);
230  skip(m_chi2);
231 
232 }
233 
234 
235 
236 void Resplot::Initialise(const std::string& name,
237  const std::vector<double>& a,
238  int n2, double a2, double b2) {
239  Initialise(name, a.size()-1, &a[0], n2, a2, b2);
240 }
241 
242 
243 void Resplot::Initialise(const std::string& name,
244  const std::vector<double>& a,
245  const std::vector<double>& b ) {
246  Initialise(name, a.size()-1, &a[0], b.size()-1, &b[0] );
247 }
248 
249 
250 
251 void Resplot::Initialise(const std::string& name,
252  int n1, const double* a1,
253  int n2, double a2, double b2) {
254  SetPrimary(n1, a1);
255  SetSecondary(n2, a2, b2);
256 
257  m_Set = true;
258 
259  m_name = name;
260  const std::string name_2d = "2d";
261  const std::string name_1d = "1d";
262 
263  // std::cout << "Resplot::Initialise() Name " << Name() << std::endl;
264 
265  // TH2D::SetDirectory(0);
266 
267  m_h2d = new TH2D( name_2d.c_str(), name_2d.c_str(), n1, a1, n2, a2, b2);
268  m_h1d = new TH1D( name_1d.c_str(), name_1d.c_str(), n2, a2, b2);
269  m_h1d->Sumw2();
270  skip(m_h2d);
271  skip(m_h1d);
272 
273  const std::string name_mean = "mean";
274  const std::string name_sigma = "sigma";
275  const std::string name_chi2 = "chi2";
276 
277  m_mean = new TH1D( name_mean.c_str(), name_mean.c_str(), n1, a1 );
278  m_sigma = new TH1D( name_sigma.c_str(), name_sigma.c_str(), n1, a1 );
279  m_chi2 = new TH1D( name_chi2.c_str(), name_chi2.c_str(), n1, a1);
280 
281  skip(m_mean);
282  skip(m_sigma);
283  skip(m_chi2);
284 
285  // const std::string name_entries = "significance";
286  const std::string name_entries = "fractional uncertainty";
287 
288  m_Nentries = new TH1D( name_entries.c_str(), name_entries.c_str(), n1, a1 );
289  skip(m_Nentries);
290 
291 }
292 
293 
294 
295 
296 void Resplot::Initialise(const std::string& name,
297  int n1, const double* a1,
298  int n2, const double* a2) {
299  SetPrimary(n1, a1);
300  SetSecondary(n2, a2[0], a2[n2+1]);
301 
302  m_Set = true;
303 
304  m_name = name;
305  const std::string name_2d = "2d";
306  const std::string name_1d = "1d";
307 
308  // std::cout << "Resplot::Initialise() Name " << Name() << std::endl;
309 
310  // TH2D::SetDirectory(0);
311 
312  m_h2d = new TH2D( name_2d.c_str(), name_2d.c_str(), n1, a1, n2, a2 );
313  m_h1d = new TH1D( name_1d.c_str(), name_1d.c_str(), n2, a2 );
314  m_h1d->Sumw2();
315  skip(m_h2d);
316  skip(m_h1d);
317 
318  const std::string name_mean = "mean";
319  const std::string name_sigma = "sigma";
320  const std::string name_chi2 = "chi2";
321 
322  m_mean = new TH1D( name_mean.c_str(), name_mean.c_str(), n1, a1 );
323  m_sigma = new TH1D( name_sigma.c_str(), name_sigma.c_str(), n1, a1 );
324  m_chi2 = new TH1D( name_chi2.c_str(), name_chi2.c_str(), n1, a1);
325 
326  skip(m_mean);
327  skip(m_sigma);
328  skip(m_chi2);
329 
330  // const std::string name_entries = "significance";
331  const std::string name_entries = "fractional uncertainty";
332 
333  m_Nentries = new TH1D( name_entries.c_str(), name_entries.c_str(), n1, a1 );
334  skip(m_Nentries);
335 
336 }
337 
338 
341 
342 Int_t Resplot::DWrite(TDirectory* g) const {
343 
344  // std::cout << "Resplot::Write() Name " << Name() << "; " << gDirectory->GetName();
345 
346  TDirectory* cwd = gDirectory;
347  if ( g ) g->cd();
348  Directory d(Name());
349  d.push();
350 
351  // std::cout << "/" << gDirectory->GetName() << std::endl;
352 
353  // TH2D::Write();
354  m_h2d->Write();
355  m_h1d->Write();
356  m_mean->Write();
357  m_sigma->Write();
358  m_chi2->Write();
359  m_Nentries->Write();
360  // std::cout << "Resplot::Write() Name() " << Name() << " " << m_slices.size() << std::endl;
361  gDirectory->mkdir("slices")->cd();
362  for ( unsigned i=0 ; i<m_slices.size() ; i++ ) m_slices[i]->Write();
363  d.pop();
364  cwd->cd();
365 
366  return 0;
367 }
368 
369 
370 
371 // Fill the helper histograms
372 
373 std::string number(const double& x) {
374  std::ostringstream s;
375  s << x;
376  return s.str();
377 }
378 
379 
380 double NCounter(TH1D* h) {
381  double N=0;
382  for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) N += h->GetBinContent(i);
383  return N;
384 }
385 
386 
387 // Get the offsets and resolutions overall and for each slice
388 int Resplot::Finalise(double a, double b, TF1* (*func)(TH1D* s, double a, double b)) {
389 
390  // gDirectory->pwd();
391 
392  s_ErrorSet = OK;
393 
394  // std::cout << "Resplot::Finalise() Name " << Name() << "\tm_n_primary " << m_n_primary << std::endl;
395 
396  gStyle->SetOptFit(1111);
397  gStyle->SetOptStat(2211);
398 
399  if ( m_finalised ) {
400  std::cerr << "Resplot::Finalise() " << m_name << " already called" << std::endl;
401  return -1;
402  }
403 
404  // std::cout << "Resplot::Finalise() " << m_name << std::endl;
405 
406  clear();
407 
408  m_mean->SetStats(false);
409  m_sigma->SetStats(false);
410 
411 #if 1
412 
413  for ( int i=1 ; i<=m_n_primary ; i++ ) {
414 
415  if ( i%1000==0 ) std::cout << "slice " << i << " of " << m_n_primary << std::endl;
416 
417  // std::string projname = tagname(m_name,i).c_str();
418  double pllim = m_mean->GetBinLowEdge(i);
419  double pulim = m_mean->GetBinLowEdge(i+1);
420 
421  // std::cout << "Finalise " << i << " : " << pllim << " - " << pulim << std::endl;
422 
423  std::string projname;
424  if ( m_xaxis=="" ) projname = "[ " + number(pllim) + " - " + number(pulim) + " ]";
425  else projname = "[ " + number(pllim) + " < " + m_xaxis + " < " + number(pulim) + " ]";
426 
427  TH1D* s;
428  if ( !m_finalised ) {
429  s = m_h2d->ProjectionY(tagname("slice",i), i, i, "e");
430 
431  // std::cout << i << " fitnull " << getWeights( s ) << std::endl;
432 
433  if ( !m_uniform ) binwidth(s);
434  s->SetTitle(projname.c_str());
435  if ( m_yaxis!="" ) s->SetXTitle(m_yaxis.c_str());
436  }
437  else s = (TH1D*)gDirectory->Get(tagname(m_name,i));
438  skip(s);
439 
440  // ZeroErrors(s);
441 
442  TF1* f1 = 0;
443 
444 
445  if ( !s_nofit ) f1 = func(s, a, b);
446 
447  // std::cout << "nofit " << s_nofit << "\tf1 " << f1 << std::endl;
448 
449 
450  // unZeroErrors(s);
451 
452 
453  if ( f1!=0 ) {
454 
455  f1->SetLineWidth(1);
456  f1->SetNpx(5000);
457  m_fitname = f1->GetName();
458 
459  // double N = NCounter(s);
460  double N = s->GetEffectiveEntries();
461 
462  if ( N!=0 ) {
463  m_Nentries->SetBinContent(i,1/sqrt(N));
464  m_Nentries->SetBinError(i,0);
465  }
466  else {
467  m_Nentries->SetBinContent(i,0);
468  m_Nentries->SetBinError(i,0);
469  }
470 
471  m_mean->SetBinContent(i, f1->GetParameter(1) );
472  m_mean->SetBinError(i, f1->GetParError(1) );
473 
474  m_sigma->SetBinContent(i, std::fabs(f1->GetParameter(2)) );
475  m_sigma->SetBinError(i, f1->GetParError(2) );
476 
477  int Ndof = f1->GetNDF();
478  if ( Ndof ) m_chi2->SetBinContent( i, std::fabs(f1->GetChisquare())/Ndof );
479  m_chi2->SetBinError(i, 0 );
480  delete f1;
481  }
482  else {
483  m_Nentries->SetBinContent(i,0);
484  m_mean->SetBinContent(i,0);
485  m_mean->SetBinError(i,0);
486  m_sigma->SetBinContent(i,0);
487  m_sigma->SetBinError(i,0);
488  m_chi2->SetBinContent(i,0);
489  m_chi2->SetBinError(i,0);
490  }
491 
492  // std::cout << "Resplot::Finalise() slice " << s->GetName() << std::endl;
493 
494  m_slices.push_back(s);
495  }
496 
497 
498 
499 #endif
500 
501 
502  m_chi2->SetMinimum(1e-5);
503 
504  // get the overall resolution and offset
505 
506  if ( !m_uniform ) binwidth(m_h1d);
507 
508  TF1* f2 = 0;
509 
510  if ( !s_nofit ) {
511  ZeroErrors(m_h1d);
512  f2 = func(m_h1d, a, b);
513  unZeroErrors(m_h1d);
514  }
515 
516  if ( f2!=0 ) {
517 
518  // std::cout << " fit = " << f2 << std::endl;
519  f2->SetLineWidth(1);
520  f2->SetNpx(5000);
521 
522  m_g_mean = StatVal( f2->GetParameter(1), f2->GetParError(1) );
523  m_g_sigma = StatVal( f2->GetParameter(2), f2->GetParError(2) );
524 
525  // std::cout << gDirectory->GetName() << "\tgmean " << m_g_mean << "\tgsigma " << g_sigma << std::endl;
526 
527  std::string mname = std::string(f2->GetParName(1));
528  std::string sname = std::string(f2->GetParName(2));
529 
530  // if ( mname!=std::string("") ) { m_mean->SetName(mname.c_str()); m_mean->SetTitle(mname.c_str()); }
531  // if ( sname!=std::string("") ) { m_sigma->SetName(sname.c_str()); m_sigma->SetTitle(sname.c_str()); }
532  if ( mname!=std::string("") ) { m_mean->SetTitle(mname.c_str()); }
533  if ( sname!=std::string("") ) { m_sigma->SetTitle(sname.c_str()); }
534 
535  delete f2;
536  }
537  else {
538  if ( !s_nofit ) std::cerr << "null overall fit Resplot:" << m_name << std::endl;
539  m_g_mean = StatVal(0,0);
540  m_g_sigma = StatVal(0,0);
541  }
542 
543 
544  if ( AddDirectoryStatus() ) SetDirectory(gDirectory);
545 
546  m_finalised = true;
547 
548  if ( s_ErrorSet!=OK ) {
549  if ( s_ErrorSet == HISTOWIDTH ) {
550  std::cerr << __FUNCTION__ << " for " << m_name
551  << " :\tdistribution wider than histogram width " << std::endl;
552  }
553  if ( s_ErrorSet == BINWIDTH ) {
554  std::cerr << __FUNCTION__ << " for " << m_name
555  << " :\tbins too wide: cannot calculate rms95 " << std::endl;
556  }
557  }
558 
559  return 0;
560 }
561 
562 
564 TH2D* Resplot::rotate(const TH2* h) {
565 
566  std::vector<double> xbins;
567  std::vector<double> ybins;
568 
569  TH1D* hy = (TH1D*)h->ProjectionX("1dx", 1, h->GetNbinsY());
570  TH1D* hx = (TH1D*)h->ProjectionY("1dy", 1, h->GetNbinsX());
571 
572  hy->SetDirectory(0);
573  hx->SetDirectory(0);
574 
575  for ( int i=0 ; i<=hx->GetNbinsX()+1 ; i++ ) xbins.push_back( hx->GetBinLowEdge(i+1) );
576  for ( int i=0 ; i<=hy->GetNbinsX()+1 ; i++ ) ybins.push_back( hy->GetBinLowEdge(i+1) );
577 
578  std::cout << "rotate:" << std::endl;
579 
580  std::cout << "x: " << xbins[0] << " - " << xbins.back() << std::endl;
581  std::cout << "y: " << ybins[0] << " - " << ybins.back() << std::endl;
582 
583  TH2D* h2 = new TH2D("duff","duff", hx->GetNbinsX(), &xbins[0], hy->GetNbinsX(), &ybins[0] );
584 
585  for ( int i=0 ; i<hx->GetNbinsX() ; i++ ) {
586  for ( int j=0 ; j<hy->GetNbinsX() ; j++ ) {
587  h2->SetBinContent(j+1,i+1, h->GetBinContent(i+1,j+1));
588  h2->SetBinError(j+1,i+1, h->GetBinError(i+1,j+1));
589  }
590  }
591 
592  h2->SetEntries( h->GetEntries() );
593 
594  delete hy;
595  delete hx;
596 
597  return h2;
598 }
599 
600 
601 
602 
606 TH2D* Resplot::combine(const TH2* h, double x, int N) {
607 
608  if ( N==0 ) return 0;
609 
610  std::vector<double> xbins;
611  std::vector<double> ybins;
612 
613  TH1D* hx = (TH1D*)h->ProjectionX("r1dx", 1, h->GetNbinsY());
614  TH1D* hy = (TH1D*)h->ProjectionY("r1dy", 1, h->GetNbinsX());
615 
616  hy->SetDirectory(0);
617  hx->SetDirectory(0);
618 
619  for ( int i=0 ; i<=hx->GetNbinsX()+1 ; i++ ) {
620  xbins.push_back( hx->GetBinLowEdge(i+1) );
621  if ( hx->GetBinLowEdge(i+1)>x ) for ( int k=1 ; k<N ; k++, i++ );
622  }
623 
624  for ( int i=0 ; i<=hy->GetNbinsX()+1 ; i++ ) ybins.push_back( hy->GetBinLowEdge(i+1) );
625 
626  if ( xbins.size()==0 || ybins.size()==0 ) {
627  std::cerr << "Resplot::combine() bad limits for histogram: N xbins: " << xbins.size() << "\tN ybins: " << ybins.size() << std::endl;
628  return 0;
629  }
630 
631  std::cout << "x: " << xbins[0] << " - " << xbins.back() << std::endl;
632  std::cout << "y: " << ybins[0] << " - " << ybins.back() << std::endl;
633 
634  TH2D* h2 = new TH2D("duff","duff", xbins.size()-1, &xbins[0], ybins.size()-1, &ybins[0] );
635  h2->SetDirectory(0);
636 
637 
638  for ( int j=0 ; j<hy->GetNbinsX() ; j++ ) {
639  int xbin = 1;
640  for ( int i=0 ; i<hx->GetNbinsX() ; i++ ) {
641  double entries = h->GetBinContent(i+1,j+1);
642  double errorsq = h->GetBinError(i+1,j+1)*h->GetBinError(i+1,j+1);
643  if ( hx->GetBinLowEdge(i+1)>x ) {
644  for ( int k=1 ; k<N ; k++, i++ ) {
645  entries += h->GetBinContent(i+2,j+1);
646  errorsq += h->GetBinError(i+2,j+1)*h->GetBinError(i+2,j+1);
647  }
648  }
649  h2->SetBinContent(xbin, j+1, entries);
650  h2->SetBinError(xbin++, j+1, std::sqrt(errorsq) );
651  }
652  }
653 
654  h2->SetEntries( h->GetEntries() );
655 
656  delete hy;
657  delete hx;
658 
659  return h2;
660 }
661 
662 
663 
667 TH1D* Resplot::combine(const TH1* h, double x, int N) {
668 
669  if ( N==0 ) return 0;
670 
671  std::vector<double> xbins;
672 
673  for ( int i=0 ; i<=h->GetNbinsX()+1 ; i++ ) {
674  xbins.push_back( h->GetBinLowEdge(i+1) );
675  if ( h->GetBinLowEdge(i+1)>x ) for ( int k=1 ; k<N ; k++, i++ );
676  }
677 
678  if ( xbins.size()==0 ) {
679  std::cerr << "Resplot::combine() bad limits for histogram: N xbins: " << xbins.size() << std::endl;
680  return 0;
681  }
682 
683  std::cout << "x: " << xbins[0] << " - " << xbins.back() << std::endl;
684 
685  TH1D* h2 = new TH1D( (std::string(h->GetName())+"-duff").c_str(),"duff", xbins.size()-1, &xbins[0] );
686  h2->SetDirectory(0);
687 
688 
689  int xbin = 1;
690  for ( int i=0 ; i<h->GetNbinsX() ; i++ ) {
691  double entries = h->GetBinContent(i+1);
692  double errorsq = h->GetBinError(i+1)*h->GetBinError(i+1);
693  if ( h->GetBinLowEdge(i+1)>x ) {
694  for ( int k=1 ; k<N ; k++, i++ ) {
695  entries += h->GetBinContent(i+2);
696  errorsq += h->GetBinError(i+2)*h->GetBinError(i+2);
697  }
698  }
699 
700  h2->SetBinContent(xbin, entries);
701  h2->SetBinError(xbin++, std::sqrt(errorsq) );
702  }
703 
704  h2->SetEntries( h->GetEntries() );
705 
706  return h2;
707 }
708 
709 
710 
714 // TH2D* Resplot::combine(const TH2* h, double epsilon) {
715 TH2D* Resplot::combine(const TH2* h, double inveps2) {
716 
717  std::cout << "combine" << std::endl;
718 
719  // if ( epsilon==0 ) return 0;
720 
721  std::vector<double> xbins;
722  std::vector<double> ybins;
723 
724  std::cout << "projection" << std::endl;
725 
726  TH1D* hx = (TH1D*)h->ProjectionX("r1dx", 1, h->GetNbinsY());
727  hx->SetDirectory(0);
728 
729 
730  // hy->SetDirectory(0);
731  // hx->SetDirectory(0);
732 
733  // double inveps2 = 1/(epsilon*epsilon);
734 
735  double N = 0;
736 
737  bool newbin = true;
738 
739  std::cout << "combining bins " << std::endl;
740 
741  for ( int i=1 ; i<=hx->GetNbinsX() ; i++ ) {
742 
743  TH1D* hy = (TH1D*)h->ProjectionY("r1dy", i, i );
744  hy->SetDirectory(0);
745  N += generate::GetEntries(hy);
746  delete hy;
747 
748  if ( xbins.size()==0 ) {
749  if ( N==0 ) continue;
750  else xbins.push_back( hx->GetBinLowEdge(i) );
751  }
752 
753  if ( newbin ) xbins.push_back( hx->GetBinLowEdge(i+1) );
754  else xbins.back() = hx->GetBinLowEdge(i+1);
755 
756  newbin = false;
757 
758  if ( xbins.size()>0 ) std::cout << i << "\tN " << N << " " << " (" << inveps2 << ")" << "\t" << xbins.back() << std::endl;
759 
760  if ( N >= inveps2 ) {
761  newbin = true;
762  N = 0;
763  }
764  }
765 
766  std::cout << "finishsed " << xbins.size() << std::endl;
767 
768  TH1D* hy = (TH1D*)h->ProjectionY("1dy", 1, h->GetNbinsX());
769  hy->SetDirectory(0);
770  for ( int i=0 ; i<=hy->GetNbinsX()+1 ; i++ ) ybins.push_back( hy->GetBinLowEdge(i+1) );
771 
772  std::cout << "combine" << std::endl;
773  std::cout << "x bins " << hx->GetNbinsX() << "\t y bins " << hy->GetNbinsX() << std::endl;
774 
775  if ( xbins.size()==0 || ybins.size()==0 ) {
776  std::cerr << "Resplot::combine() bad limits for histogram: N xbins: " << xbins.size() << "\tN ybins: " << ybins.size() << std::endl;
777  return 0;
778  }
779 
780  std::cout << "x: " << xbins.size() << " " << xbins[0] << " - " << xbins.back() << std::endl;
781  std::cout << "y: " << ybins.size() << " " << ybins[0] << " - " << ybins.back() << std::endl;
782 
783  TH2D* h2 = new TH2D( (std::string(h->GetName())+"-duff").c_str(), "duff", xbins.size()-1, &xbins[0], ybins.size()-1, &ybins[0] );
784  h2->SetDirectory(0);
785 
786 
787  for ( int j=1 ; j<=hy->GetNbinsX() ; j++ ) {
788 
789  unsigned xbin = 0;
790 
791  for ( int i=1 ; i<=hx->GetNbinsX() ; i++ ) {
792 
793  while ( hx->GetBinCenter(i)>xbins[xbin+1] && xbin<(xbins.size()-1) ) xbin++;
794 
795  if ( hx->GetBinCenter(i)>=xbins[xbin] && hx->GetBinCenter(i)<xbins[xbin+1] ) {
796 
798  double n = h->GetBinContent(i,j);
799  double nh = h2->GetBinContent(xbin+1, j);
800  h2->SetBinContent(xbin+1, j, nh+n);
801 
803  double ne = h->GetBinError(i,j);
804  double nhe = h2->GetBinError(xbin+1, j);
805  h2->SetBinError(xbin+1, j, std::sqrt( ne*ne + nhe*nhe ) );
806 
807  }
808  }
809 
810  }
811 
812  h2->SetEntries( h->GetEntries() );
813 
814  delete hy;
815  delete hx;
816 
817  return h2;
818 }
819 
820 
821 
822 
823 
824 
825 
826 
827 
828 
832 // TH2D* Resplot::combine(const TH2* h, double epsilon) {
833 TH1D* Resplot::combine(const TH1* h, double inveps2) {
834  return rebin( h, findbins(h, inveps2) );
835 }
836 
837 
838 std::vector<double> Resplot::findbins(const TH1* h, double inveps2) {
839 
840  std::cout << "combine" << std::endl;
841 
842  // if ( epsilon==0 ) return 0;
843 
844  std::vector<double> xbins;
845 
846 
847  // hy->SetDirectory(0);
848  // hx->SetDirectory(0);
849 
850  // double inveps2 = 1/(epsilon*epsilon);
851 
852  double N = 0;
853 
854  bool newbin = true;
855 
856  std::cout << "combining bins " << std::endl;
857 
858  for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) {
859 
860  N += h->GetBinContent(i);
861 
862  if ( xbins.size()==0 ) {
863  if ( N==0 ) continue;
864  else xbins.push_back( h->GetBinLowEdge(i) );
865  }
866 
867  if ( newbin ) xbins.push_back( h->GetBinLowEdge(i+1) );
868  else xbins.back() = h->GetBinLowEdge(i+1);
869 
870  newbin = false;
871 
872  if ( xbins.size()>0 ) std::cout << i << "\tN " << N << " " << " (" << inveps2 << ")" << "\t" << xbins.back() << std::endl;
873 
874  if ( N >= inveps2 ) {
875  newbin = true;
876  N = 0;
877  }
878  }
879 
880  std::cout << "finishsed " << xbins.size() << std::endl;
881 
882  std::cout << "combine" << std::endl;
883  std::cout << "x bins " << h->GetNbinsX() << std::endl;
884 
885  std::cout << "x: " << xbins.size() << " " << xbins[0] << " - " << xbins.back() << std::endl;
886 
887  return xbins;
888 }
889 
890 
891 TH1D* Resplot::rebin( const TH1* h, const std::vector<double>& xbins ) {
892 
893  if ( xbins.size()==0 ) {
894  std::cerr << "Resplot::combine() bad limits for histogram: N xbins: " << xbins.size() << std::endl;
895  return 0;
896  }
897 
898  TH1D* h2 = new TH1D("duff","duff", xbins.size()-1, &xbins[0] );
899  h2->SetDirectory(0);
900 
901  unsigned xbin = 0;
902  for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) {
903 
904  while ( h->GetBinCenter(i)>xbins[xbin+1] && xbin<(xbins.size()-1) ) xbin++;
905 
906  if ( h->GetBinCenter(i)>=xbins[xbin] && h->GetBinCenter(i)<xbins[xbin+1] ) {
907  double n = h->GetBinContent(i);
908  double nh = h2->GetBinContent(xbin+1);
909  h2->SetBinContent(xbin+1, nh+n);
910  }
911  }
912 
913  return h2;
914 }
915 
916 
917 
918 
919 
920 
921 std::vector<double> Resplot::getbins(const TH1* h) {
922  std::vector<double> bins(h->GetNbinsX()+1);
923  for ( int i=0 ; i<h->GetNbinsX()+1 ; i++ ) bins[i] = h->GetBinLowEdge(i+1);
924  return bins;
925 }
926 
927 
928 
929 
930 // Get the efficiency for a single slice of the major ordinate
931 
932 StatVal Resplot::GetEfficiency(TH1D* h, double lower, double upper) {
933  double total = 0;
934  double good = 0;
935  StatVal v;
936 
937  for (int i=0 ; i<=h->GetNbinsX()+1 ; i++ ) {
938  total += h->GetBinContent(i);
939  }
940  for ( int i=1 ; i<=h->GetNbinsX() ; i++ ) {
941  double x = h->GetBinCenter(i);
942  if ( x>lower && x<= upper ) {
943  good += h->GetBinContent(i);
944  }
945  }
946  if ( total>0 ) {
947  v.value = good/total;
948  v.error = sqrt(v.value*(1-v.value)/total);
949  }
950 
951  return v;
952 }
953 
954 
955 
956 // Get all efficiencies for all the major bins
957 
958 TH1D* Resplot::GetEfficiencies(const std::string& hname, double lower, double upper) {
959  // TDirectory* cwd = cd();
960 
961  // push();
962 
963  TH1D* e = (TH1D*) m_mean->Clone();
964  e->SetName(hname.c_str()); e->SetMarkerStyle(20);
965  if ( m_xaxis!="" ) e->SetXTitle(m_xaxis.c_str());
966 
967 
968  for ( int i=1 ; i<=m_n_primary ; i++ ) {
969  TH1D* s = m_h2d->ProjectionY(tagname(hname,i), i, i, "e");
970  StatVal v = GetEfficiency(s, lower, upper);
971  e->SetBinContent(i, v.value);
972  e->SetBinError(i, v.error);
973  delete s;
974  }
975 
976  // cwd->cd();
977 
978  // pop();
979 
980  return e;
981 }
982 
983 
984 
985 
986 
987 // Get all efficiencies for all the major bins
988 
989 TH1D* Resplot::GetEfficiencies(const std::string& hname, double Nsigma ) {
990  // TDirectory* cwd = cd();
991  // push();
992 
993  // check that the fits have been done for this
994  // plot
995  if ( !m_finalised ) Finalise();
996 
997  TH1D* e = (TH1D*) m_mean->Clone();
998  e->SetName(hname.c_str()); e->SetMarkerStyle(20);
999  if ( m_xaxis!="" ) e->SetXTitle(m_xaxis.c_str());
1000 
1001 
1002  for ( int i=1 ; i<=m_n_primary ; i++ ) {
1003 
1004  double mean = m_mean->GetBinContent(i);
1005  double sigma = m_sigma->GetBinContent(i);
1006 
1007  TH1D* s = m_h2d->ProjectionY(tagname(hname,i), i, i, "e");
1008  StatVal v = GetEfficiency(s, mean-Nsigma*sigma, mean+Nsigma*sigma);
1009  e->SetBinContent(i, v.value);
1010  e->SetBinError(i, v.error);
1011  delete s;
1012  }
1013 
1014  // cwd->cd();
1015 
1016  // pop();
1017 
1018  return e;
1019 }
1020 
1021 
1022 
1023 // fit Gaussian to histogram
1024 
1025 TF1* Resplot::FitGaussian(TH1D* s, double a, double b) {
1026  TF1* f1 = new TF1("gausr", "gaus");
1027  return FitInternal( s, a, b, f1 );
1028  // s->GetXaxis()->SetRangeUser(-1, 1);
1029 }
1030 
1031 TF1* Resplot::FitRGaussian(TH1D* s, double a, double b) {
1032  TF1* f1 = new TF1("rgausr", "[0]*x*exp((x-[1])([1]-x)/(2*[2]*[2]))");
1033  return FitInternal( s, a, b, f1 );
1034 }
1035 
1036 TF1* Resplot::FitRBreit(TH1D* s, double a, double b) {
1037  TF1* f1 = new TF1("rbreitr", "x*[0]*[2]*[2]/([2]*[2]+(x-[1])*(x-[1]))");
1038 
1039  f1->SetParName(0, "Maximum");
1040  f1->SetParName(1, "Median");
1041  f1->SetParName(2, "Gamma");
1042 
1043  f1->SetParameter(0,s->GetBinContent(s->GetMaximumBin()) );
1044  f1->SetParameter(1,s->GetMean());
1045  f1->SetParameter(2,s->GetRMS());
1046 
1047  return FitInternal( s, a, b, f1 );
1048 }
1049 
1050 
1051 TF1* Resplot::FitInternal(TH1D* s, double a, double b, TF1* f1) {
1052 
1053  if ( f1==0 ) return f1;
1054 
1055  f1->SetNpx(5000);
1056  f1->SetLineWidth(1);
1057 
1058  double mean = 0;
1059  double sigma = 0;
1060 
1061 
1062  f1->SetParameter(0,s->GetBinContent(s->GetMaximumBin()) );
1063  f1->SetParameter(1,mean=s->GetBinCenter(s->GetMaximumBin()));
1064  f1->SetParameter(2,sigma=s->GetRMS());
1065 
1066  int nbins=0;
1067  for (int j=1 ; j<=s->GetNbinsX() ; j++) if (s->GetBinContent(j)) nbins++;
1068 
1069  if (nbins>2) {
1070  if ( a==-999 || b==-999 ) s->Fit(f1,"Q");
1071  else s->Fit(f1,"Q", "", a, b);
1072  }
1073  else for ( int j=0 ; j<3 ; j++ ) f1->SetParameter(j, 0);
1074 
1075  return f1;
1076 }
1077 
1078 
1079 
1080 TF1* Resplot::FitCentralGaussian(TH1D* s, double , double ) {
1081  TF1* f1 = new TF1("gausr", "gaus");
1082 
1083  f1->SetNpx(5000);
1084  f1->SetLineWidth(1);
1085 
1086  double mean = 0;
1087  double sigma = 0;
1088 
1089 
1090  f1->SetParameter(0,s->GetBinContent(s->GetMaximumBin()) );
1091  f1->SetParameter(1,mean=s->GetBinCenter(s->GetMaximumBin()));
1092  f1->SetParameter(2,sigma=s->GetRMS());
1093 
1094  int nbins=0;
1095  for (int j=1 ; j<=s->GetNbinsX() ; j++) if (s->GetBinContent(j)) nbins++;
1096 
1097  if (nbins>2) {
1098 
1099  TH1D* s_tmp = (TH1D*)s->Clone("duff");
1100  s_tmp->SetDirectory(0);
1101  s_tmp->Fit(f1,"Q");
1102 
1103  if ( f1 ) {
1105 
1107 
1108  // double sig = std::sqrt( f1->GetParError(0)*f1->GetParError(0) + s->GetBinError(s->GetMaximumBin())*s->GetBinError(s->GetMaximumBin()) );
1109  // double diff = std::fabs( (f1->GetParameter(0) - s->GetBinContent(s->GetMaximumBin()) ) )/sig;
1110 
1111  // std::cout << "diff " << diff << "\tfit " << f1->GetParameter(0) << "\thist " << s->GetBinContent(s->GetMaximumBin()) << std::endl;
1112 
1113  // if ( diff>2 ) {
1114 
1115 
1116  // double n = s->GetBinContent(s->GetMaximumBin());
1117 
1118  double frac = 1.25;
1119 
1120  // if ( n<25 ) {
1121  // frac = 1.2*n/25.0 + 4.0*(25.0-n)/25.0;
1122  // }
1123 
1124  double llim = f1->GetParameter(1)-frac*f1->GetParameter(2);
1125  double ulim = f1->GetParameter(1)+frac*f1->GetParameter(2);
1126 
1127 
1128 
1129  int nbins=0;
1130  for (int j=1 ; j<=s->GetNbinsX() ; j++) if ( s->GetBinCenter(j)>llim && s->GetBinCenter(j)<ulim ) nbins++;
1131 
1132  if ( nbins>2 ) {
1133 
1134  // if ( frac>3 ) s->Fit(f1,"Q");
1135  // else s->Fit(f1,"Q", "", llim, ulim);
1136 
1137  if ( frac>3 ) s->Fit(f1,"Q");
1138  else s->Fit(f1,"Q", "", llim, ulim);
1139 
1140  }
1141  else for ( int j=0 ; j<3 ; j++ ) f1->SetParameter(j, 0);
1142 
1143  delete s_tmp;
1144 
1145  }
1146 
1147  // std::cout << "Resplot::Finalise() " << s->GetName() << "\trange " << lower << " - " << upper << std::endl;
1148 
1149  }
1150  else for ( int j=0 ; j<3 ; j++ ) f1->SetParameter(j, 0);
1151 
1152  return f1;
1153 }
1154 
1155 
1156 
1157 // or fit a Breit Wigner
1158 
1159 TF1* Resplot::FitBreit(TH1D* s, double a, double b) {
1160 
1161  // TF1* f1 = new TF1("breitr", "[0]*[2]*[2]/([2]*[2]+(x-[1])*(x-[1]))", -100, 100 );
1162  TF1* f1 = new TF1("breitr", "[0]*[2]*[2]/([2]*[2]+(x-[1])*(x-[1]))");
1163 
1164  f1->SetParName(0, "Maximum");
1165  f1->SetParName(1, "Median");
1166  f1->SetParName(2, "Gamma");
1167 
1168  f1->SetParameter(0,s->GetBinContent(s->GetMaximumBin()) );
1169  f1->SetParameter(1,s->GetMean());
1170  f1->SetParameter(2,s->GetRMS()*0.5);
1171 
1172  f1->SetNpx(5000);
1173  f1->SetLineWidth(1);
1174 
1175  for (int j=1 ; j<=s->GetNbinsX() ; j++) if (s->GetBinContent(j)==0) s->SetBinError(j, 1);
1176 
1177 
1178  int nbins=0;
1179  for (int j=1 ; j<=s->GetNbinsX() ; j++) if (s->GetBinContent(j)>0) nbins++;
1180 
1181  if (nbins>2) {
1182  // std::cout << "fitting " << f1 << std::endl;
1183  if ( a==-999 || b==-999 ) s->Fit(f1,"Q");
1184  else s->Fit(f1,"Q", "", a, b);
1185  }
1186  else {
1187  for ( int j=0 ; j<3 ; j++ ) f1->SetParameter(j,0);
1188  }
1189 
1190  for (int j=1 ; j<=s->GetNbinsX() ; j++) if (s->GetBinContent(j)==0) s->SetBinError(j, 0);
1191 
1192  return f1;
1193 }
1194 
1195 
1196 
1197 
1198 
1199 
1200 
1201 
1202 
1203 
1204 
1205 
1206 
1207 
1208 
1209 
1210 
1211 
1212 
1213 
1214 
1215 
1216 
1217 
1218 
1219 
1220 TF1* Resplot::FitATan(TH1D* s, double , double ) {
1221  static int counter = 0;
1222 
1223  std::cout << "Resplot::FitATan() " << s->GetName() << " " << s->GetTitle() << std::endl;
1224 
1225  char name[8];
1226  sprintf(name,"atan_%02d", counter); counter++;
1227 
1228  // Cannot use "M_PI" here as root does not replace M_PI in the string expressions for function in TF1 definition
1229  TF1* f1 = new TF1(name, "100*pow( (0.01*[3]-sqrt(0.0001*[2]*[2]))*(0.5-(1/3.14159)*atan( [1]*(x-[0]) ) )+sqrt(0.0001*[2]*[2]) , [4] ) ", -100, 100);
1230 
1231  f1->SetParName(0, "Mean");
1232  f1->SetParName(1, "Gradient");
1233  f1->SetParName(2, "Pedestal");
1234  f1->SetParName(3, "Norm");
1235  f1->SetParName(4, "Power");
1236 
1237  int maxn=1;
1238  for ( int i=1 ; i<=s->GetNbinsX() ; i++ ) if ( s->GetBinContent(maxn)<s->GetBinContent(i) ) maxn = i;
1239 
1240 
1241  // f1->SetParameter(0,s->GetBinContent(maxn) );
1242  // f1->SetParameter(1,s->GetBinCenter(s->GetNbinsX()/2) );
1243  // f1->SetParameter(2,0.25);
1244 
1245  f1->SetParameter(0, 2 );
1246  f1->SetParameter(1, 10 );
1247  f1->SetParameter(2, 0.001 );
1248  f1->SetParameter(3, 100 );
1249  f1->SetParameter(4, 1 );
1250 
1251  f1->FixParameter(2, 0);
1252  f1->FixParameter(3, 100);
1253  // f1->FixParameter(4, 1);
1254 
1255  f1->SetNpx(5000);
1256  f1->SetLineWidth(1);
1257 
1258 
1259  // int nbins=0;
1260 
1261 
1262  // ZeroErrors(s);
1263  s->SetMarkerStyle(20);
1264 
1265 
1266  // int n=0;
1267  // for ( int i=1 ; i<=s->GetNbinsX() ; i++ ) if ( s->GetBinContent() ) n++;
1268  // if ( n>2 ) s->Fit(f1,"Q");
1269 
1270  s->Fit(f1,"Q");
1271 
1272  std::cout << "par0 = " << f1->GetParameter(0) << std::endl;
1273  std::cout << "par1 = " << f1->GetParameter(1) << std::endl;
1274  std::cout << "par2 = " << f1->GetParameter(2) << std::endl;
1275  std::cout << "par3 = " << f1->GetParameter(3) << std::endl;
1276  std::cout << "par4 = " << f1->GetParameter(4) << std::endl;
1277 
1278 
1279  // unZeroErrors(s);
1280  // s->SetMaximum(f1->GetParameter(0)*1.1);
1281  // double range = s->GetBinLowEdge(s->GetNbinsX()+1) - s->GetBinLowEdge(1) ;
1282 
1283  return f1;
1284 }
1285 
1286 
1287 
1288 TF1* Resplot::FitNull(TH1D* s, double , double ) {
1289 
1290  // unZeroErrors(s);
1291 
1292  TF1* f=new TF1("null", "[0]+[1]+[2]");
1293 
1294  f->SetParName(0, "Maximum");
1295  f->SetParName(1, "Mean");
1296  f->SetParName(2, "RMS");
1297 
1298  f->FixParameter(0, s->GetMaximum()); f->SetParError(0, sqrt(s->GetMaximum()));
1299  f->FixParameter(1, s->GetMean()); f->SetParError(1, s->GetMeanError());
1300  f->FixParameter(2, s->GetRMS()); f->SetParError(2, s->GetRMSError());
1301 
1302  // std::cout << gDirectory->GetName() << " " << s->GetName() << "\tFitNull mean: " << f->GetParameter(1) << "\tsigma: " << f->GetParameter(2) << std::endl;
1303 
1304  return f;
1305 }
1306 
1307 
1308 
1309 double FindMean(TH1D* s, double frac=0.95) {
1310 
1311  // double entries = s->GetEntries() + s->GetBinContent(0) + s->GetBinContent(s->GetNbinsX()+1);
1312  double entries = generate::GetEntries(s,0,s->GetNbinsX()+1);
1313 
1314  // std::cout << s->GetName() << "\tentries: " << entries << std::endl;
1315 
1316  if ( entries==0 ) return 0;
1317 
1318  // s->GetXaxis()->SetRangeUser(-100, 100);
1319  // s->GetXaxis()->UnZoom();
1320  s->GetXaxis()->SetRange(1,s->GetNbinsX());
1321 
1322  double mean = s->GetMean();
1323  double meane = s->GetMeanError();
1324 
1325  int upperbin = 0;
1326  int lowerbin = 0;
1327 
1328 
1329  for ( int it=0 ; it<20 ; it++ ) {
1330 
1331  // std::cout << it << "\tmean " << mean << " +- " << meane << std::endl;
1332 
1333  int imax = s->GetXaxis()->FindBin(mean);
1334 
1335  double sumn = s->GetBinContent(imax);
1336 
1337  // double uppersum = 0;
1338  // double lowersum = 0;
1339 
1340  upperbin = imax;
1341  lowerbin = imax;
1342 
1343 
1344  int i=1;
1345  while ( true ) {
1346 
1347  const int upperbin_i = imax+i;
1348  const int lowerbin_i = imax-i;
1349 
1350  if ( upperbin_i>s->GetNbinsX() || lowerbin_i<1 ) break;
1351 
1352  double tsum = sumn + s->GetBinContent(upperbin_i) + s->GetBinContent(lowerbin_i);
1353 
1354  if ( tsum>entries*frac ) {
1355  // uppersum = tsum/entries;
1356  break;
1357  }
1358 
1359  // lowersum = tsum/entries;
1360 
1361  sumn = tsum;
1362 
1363  upperbin = upperbin_i;
1364  lowerbin = lowerbin_i;
1365 
1366  i++;
1367  }
1368 
1369  s->GetXaxis()->SetRange(lowerbin, upperbin);
1370 
1371  double m = s->GetMean();
1372  double me = s->GetMeanError();
1373 
1374  if ( it>0 ) {
1375  if ( mean==m ||
1376  std::fabs(mean-m)<me*1e-5 ||
1377  std::fabs(mean-m)<meane*1e-5 ) {
1378  mean = m;
1379  meane = me;
1380  break;
1381  }
1382  }
1383 
1384  mean = m;
1385  meane = me;
1386 
1387  }
1388 
1389  // std::cout << "mean " << mean << " +- " << meane << "\tentries: " << generate::GetEntries(s) << "\tlower " << lowerbin << " - " << upperbin << std::endl;
1390 
1391 
1392  return mean;
1393 }
1394 
1395 
1396 int findmax(TH1D* s) {
1397  int imax = 1;
1398  for ( int i=2 ; i<=s->GetNbinsX() ; i++ ) {
1400  if ( s->GetBinContent(i)>s->GetBinContent(imax) ) imax = i;
1401  }
1402  return imax;
1403 }
1404 
1405 
1406 
1407 // TF1* Resplot::FitNull95Obsolete(TH1D* s, double frac, bool useroot ) {
1408 TF1* Resplot::FitNull95Obsolete(TH1D* s, double frac, bool useroot ) {
1409 
1410  TF1* f=new TF1("null", "[0]+[1]+[2]");
1411 
1412  f->SetParName(0, "Maximum");
1413  f->SetParName(1, "Mean");
1414  f->SetParName(2, "RMS");
1415 
1416  f->SetParameter(0, 0);
1417  f->SetParameter(1, 0);
1418  f->SetParameter(2, 0);
1419 
1420  if ( s==0 ) { std::cerr << "Resplot::FitNull95() histogram s = " << s << std::endl; return 0; }
1421  if ( s->GetXaxis()==0) { std::cerr << "Resplot::FitNull95() histogram s->GetXaxis() not definied for histogram " << s->GetName() << std::endl; return 0; }
1422 
1423  unZeroErrors(s);
1424 
1425  s->GetXaxis()->SetRange(1,s->GetNbinsX());
1426 
1427 #if 0
1428  std::cout << "FitNull95 " << s->GetName()
1429  << " " << generate::GetEntries(s)
1430  << " " << generate::GetEntries(s, 0, s->GetNbinsX()+1)
1431  << std::endl;
1432 #endif
1433 
1434 
1435  if ( generate::GetEntries(s)<10 ) return f;
1436 
1437  // std::cout << "---------------------------------------\nFN95O frac: " << frac << " " << s << std::endl;
1438 
1439  // std::cout << __LINE__ << std::endl;
1440 
1441 
1442 
1443  // std::cout << __LINE__ << " " << s->GetXaxis() << std::endl;
1444 
1447  // s->GetXaxis()->UnZoom();
1448 
1449  // std::cout << __LINE__ << std::endl;
1450 
1451  // double m = s->GetMean();
1452  double m = FindMean(s,frac);
1453 
1454  // std::cout << __LINE__ << std::endl;
1455 
1456  int imax = s->GetXaxis()->FindBin(m);
1457 
1458  // std::cout << "mean " << m << " " << imax << " max " << s->GetBinCenter(findmax(s)) << " " << findmax(s) << std::endl;
1459 
1460  f->FixParameter(0, s->GetMaximum()); f->SetParError(0, std::sqrt(s->GetMaximum()));
1461  // f->FixParameter(1, s->GetMean()); f->SetParError(1, s->GetMeanError());
1462 
1463  // double entries = s->GetEntries() + s->GetBinContent(0) + s->GetBinContent(s->GetNbinsX()+1);
1464  double entries = generate::GetEntries(s,0,(s->GetNbinsX()+1));
1465 
1466  // std::cout << "entries: " << entries << "\t" << generate::GetEntries(s) << "\t" << generate::GetEntries(s,0,s->GetNbinsX()+1) << std::endl;
1467 
1468  // double sumx = s->GetBinContent(imax)*s->GetBinCenter(imax);
1469  // double sumx2 = s->GetBinContent(imax)*s->GetBinCenter(imax)*s->GetBinCenter(imax);
1470 
1471  double sumn = s->GetBinContent(imax);
1472 
1473  int upperbin = imax;
1474  int lowerbin = imax;
1475 
1476  double uppersum = 0;
1477  double lowersum = 0;
1478 
1479  if ( sumn>entries*frac ) {
1480  s->GetXaxis()->SetRange(1,s->GetNbinsX());
1481  s_ErrorSet = BINWIDTH;
1482  return f;
1483  }
1484 
1485  if ( entries>0 ) {
1486 
1487  int i=1;
1488  while ( true ) {
1489 
1490  const int upperbin_i = imax+i;
1491  const int lowerbin_i = imax-i;
1492 
1493  if ( upperbin_i>s->GetNbinsX() || lowerbin_i<1 ) break;
1494 
1495  double tsum = sumn + s->GetBinContent(upperbin_i) + s->GetBinContent(lowerbin_i);
1496 
1497  // std::cout << i << " frac: " << lowersum
1498  // << "\tx " << s->GetBinCenter(lowerbin)
1499  // << "\t " << s->GetBinCenter(upperbin)
1500  // << "\ttsum " << tsum
1501  // << "\tentries " << entries
1502  // << std::endl;
1503 
1504  if ( tsum>entries*frac ) { uppersum = tsum/entries; break; }
1505 
1506  lowersum = tsum/entries;
1507 
1508  sumn = tsum;
1509 
1510  upperbin = upperbin_i;
1511  lowerbin = lowerbin_i;
1512 
1513  // if ( std::fabs(lowersum-frac)<0.1 )
1514 
1515  i++;
1516  }
1517  }
1518 
1519  if ( uppersum==0 ) {
1520  s->GetXaxis()->SetRange(1,s->GetNbinsX());
1521  s_ErrorSet = HISTOWIDTH;
1522  return f;
1523  }
1524 
1525  // std::cout << "FitNull95 " << s->GetName() << "\tlower bin " << lowerbin << "\t upper bin " << upperbin << std::endl;
1526 
1527  // std::cout << "full rms " << s->GetRMS() << " +- " << s->GetRMSError() << "\t n=" << s->GetNbinsX()
1528  // << "\tfrac " << ( entries>0 ? sumn/entries : 0 ) << " - " << uppersum
1529  // << std::endl;
1530 
1531  // std::cout << "upper " << upperbin << " " << lowerbin << " " << interpolate_flag << std::endl;
1532 
1533  // if ( upperbin!=lowerbin ) {
1534  if ( true ) { // technically wrong - 95% contained in only 3 bins, so RMS of the "central" bin is not defined
1535 
1536  double llim = s->GetBinLowEdge(1);
1537  double ulim = s->GetBinLowEdge(s->GetNbinsX()+1);
1538 
1539  double rms = 0;
1540  double erms = 0;
1541 
1542  std::vector<double> stats;
1543 
1544  // double sentries = generate::GetEntries(s);
1545 
1546  // std::cout << "GetEntries " << s->GetName() << " " << sentries << std::endl;
1547 
1548  if ( s_interpolate_flag ) {
1549  s->GetXaxis()->SetRange(lowerbin, upperbin);
1550 
1551  // std::cout << "GetEntries " << s->GetName() << " " << generate::GetEntries(s,lowerbin,upperbin)/sentries << std::endl;
1552 
1553  double rlower;
1554  double erlower;
1555 
1556  double rupper;
1557  double erupper;
1558 
1559 
1561 
1562  if ( !useroot ) {
1563  GetStats( s, stats );
1564  rlower = stats[2];
1565  erlower = stats[3];
1566  }
1567  // else {
1568  rlower = s->GetRMS();
1569  erlower = s->GetRMSError();
1570  // }
1571 
1572  // std::cout << "FitNULL95() " << s->GetName() << "\tmean " << m << "\t lower rms " << rlower << " +- " << s->GetRMSError() << " (wrong)\t " << stats[2] << " +- " << stats[3] << " (correct)" << std::endl;
1573  // std::cout << s->GetName() << " " << s->GetTitle() << "\t lower 95% rms " << rlower << " +- " << erlower << "\t( " << lowerbin << " - " << upperbin << " )" << std::endl;
1574 
1575  s->GetXaxis()->SetRange(lowerbin-1, upperbin+1);
1576 
1577  // std::cout << "GetEntries " << s->GetName() << " " << generate::GetEntries(s,lowerbin-1,upperbin+1)/sentries << std::endl;
1578 
1579  // std::cout << "Range " << s->GetBinLowEdge(lowerbin) << " - " << s->GetBinLowEdge(upperbin+1) << std::endl;
1580 
1581  if ( !useroot ) {
1582  GetStats( s, stats );
1583  rupper = stats[2];
1584  erupper = stats[3];
1585  }
1586  // else {
1587  rupper = s->GetRMS();
1588  erupper = s->GetRMSError();
1589  // }
1590 
1591  // std::cout << "FitNULL95() " << s->GetName() << "\t upper rms " << rupper << " +- " << s->GetRMSError() << " (wrong)\t " << stats[2] << " +- " << stats[3] << " (correct)" << std::endl;
1592 
1593  // std::cout << s->GetName() << " " << s->GetTitle() << "\t upper 95% rms " << rupper << " +- " << erupper << "\t( " << lowerbin-1 << " - " << upperbin+1 << " )" << std::endl;
1594 
1595 
1596 
1597 
1598  rms = rlower + (frac-lowersum)*(rupper-rlower)/(uppersum-lowersum);
1599  erms = erlower + (frac-lowersum)*(erupper-erlower)/(uppersum-lowersum);
1600 
1601  }
1602  else {
1603  s->GetXaxis()->SetRange(lowerbin, upperbin);
1604  if ( !useroot ) {
1605  GetStats( s, stats );
1606  rms = stats[2];
1607  erms = stats[3];
1608  }
1609  else {
1610  rms = s->GetRMS();
1611  erms = s->GetRMSError();
1612  }
1613  }
1614 
1615 #if 0
1616 
1617  std::cout << "numbers " << s->GetName()
1618  << " " << generate::GetEntries(s,lowerbin,upperbin)
1619  << " - " << generate::GetEntries(s,lowerbin-1,upperbin+1) << "\ttotal " << generate::GetEntries(s,0,s->GetNbinsX()+1) << std::endl;
1620 
1621  std::cout << "fractions " << s->GetName()
1622  << "\t" << generate::GetEntries(s,lowerbin,upperbin)/generate::GetEntries(s,0,s->GetNbinsX()+1)
1623  << " - " << generate::GetEntries(s,lowerbin-1,upperbin+1)/generate::GetEntries(s,0,s->GetNbinsX()+1) << std::endl;
1624 
1625  std::cout << "lowersum " << s->GetName() << "\t" << lowersum << "\tuppersum " << uppersum << std::endl;
1626 
1627  std::cout << "limits " << s->GetName() << "\t" << s->GetBinLowEdge(lowerbin) << " - " << s->GetBinLowEdge(upperbin+1) << std::endl;
1628 
1629 
1630 
1631  std::cout << "FitNULL95Obsolete() " << s->GetName() << "\t rms " << rms << " +- " << erms << "\t inv " << 1/rms << std::endl;
1632 
1633  printf("rms %12.10lf inv %12.10lf\n", rms, 1/rms );
1634 
1635  printf("rms68 %12.10lf inv %12.10lf\n", rms/0.5369760683, rms*1.8622803865 );
1636 #endif
1637 
1638 
1641  double scale = 1;
1642  if ( s_scalerms95 ) scale = 1.1479538518;
1643 
1644 
1645 
1646  // f->FixParameter(2, s->GetRMS());
1647  // f->SetParError(2, s->GetRMSError());
1648  f->FixParameter(2, scale*rms);
1649  f->SetParError(2, scale*erms);
1650  f->FixParameter(1, s->GetMean());
1651  f->SetParError(1, s->GetMeanError());
1652  // std::cout << s->GetName() << " " << s->GetTitle() << "\t 95% inter rms " << rms << " +- " << erms << std::endl;
1653  s->GetXaxis()->SetRangeUser(llim, ulim);
1654 
1655  }
1656  else {
1657  f->FixParameter(1, 0);
1658  f->SetParError(1, 0);
1659  f->FixParameter(2, 0);
1660  f->SetParError(2, 0);
1661  // std::cout << " 95% rms " << 0 << " +- " << 0 << "\t( " << lowerbin << " - " << upperbin << " )" << std::endl;
1662  }
1663 
1664  s->GetXaxis()->SetRange(1,s->GetNbinsX());
1665 
1666  return f;
1667 }
1668 
1669 
1670 
1671 
1672 
1673 TF1* Resplot::FitNull95(TH1D* s, double , double ) {
1674  // if ( s_oldrms95 ) return Resplot::FitNull95Obsolete(s, 0.95, true ); /// use approximate root errors for speed
1675  if ( s_oldrms95 ) return Resplot::FitNull95Obsolete( s, 0.95 );
1676  else return Resplot::FitNull95New(s);
1677 }
1678 
1679 
1680 
1681 int nexperiments_ = 0;
1682 int nexperiments_max = 40;
1683 int nexperiments_min = 20;
1684 
1685 
1686 
1687 TF1* Resplot::FitNull95New(TH1D* s, double, bool ) { // double frac, bool useroot ) {
1688 
1689  // std::cout << "---------------------------------------\nFN95 " << std::endl;
1690 
1691  TF1* f=new TF1("null", "[0]+[1]+[2]");
1692 
1693  f->SetParName(0, "Maximum");
1694  f->SetParName(1, "Mean");
1695  f->SetParName(2, "RMS");
1696 
1697  f->SetParameter(0, 0);
1698  f->SetParameter(1, 0);
1699  f->SetParameter(2, 0);
1700 
1701  if ( s->GetEffectiveEntries()==0 ) return f;
1702 
1703 
1704  // double _entries = getWeights( s );
1705  const double entries = s->GetEffectiveEntries();
1706 
1707 
1708  // int nexperiments = 20+int(800/std::sqrt(s->GetEntries()));
1709  // int nexperiments = 20+int(800/std::sqrt(s->GetEntries()));
1710  int nexperiments = 20+int(1000/entries);
1711 
1712  if ( nexperiments_!=0 ) nexperiments = nexperiments_;
1713 
1714 
1715  if ( nexperiments>nexperiments_max ) nexperiments = nexperiments_max;
1716  if ( nexperiments<nexperiments_min ) nexperiments = nexperiments_min;
1717 
1718  // std::cout << "FitNull95 entries " << entries << "\tnexperiments " << nexperiments << std::endl;
1719 
1720  // std::cout << "experiments " << nexperiments_ << std::endl;
1721 
1722  // std::cout << "h entries " << s->GetEntries() << "\tn experiments " << nexperiments << "\tN " << s->GetEntries()*nexperiments << std::endl;
1723 
1725  generate::experiment e( s, nexperiments, entries );
1726 
1729  double scale = 1;
1730  if ( s_scalerms95 ) scale = 1.1479538518;
1731 
1732  // f->FixParameter( 1, s->GetMean() );
1733  f->FixParameter( 1, e.hmean() );
1734  f->SetParError( 1, e.mean_error() );
1735 
1736  f->FixParameter( 2, scale*e.hrms() );
1737  f->SetParError( 2, scale*e.rms_error() );
1738 
1739  // s->GetXaxis()->SetRange(1, s->GetNbinsX());
1740 
1741  // s->GetXaxis()->SetRangeUser(-1, 1);
1742 
1743  return f;
1744 }
1745 
1746 
1747 
1748 
1749 
1750 
1751 
1752 
1753 
1754 TF1* Resplot::FitNull95Central(TH1D* s) {
1755 
1756  TF1* f=new TF1("null", "[0]+[1]+[2]");
1757 
1758  f->SetParName(0, "Maximum");
1759  f->SetParName(1, "Mean");
1760  f->SetParName(2, "RMS");
1761 
1762  f->SetParameter(0, 0);
1763  f->SetParameter(1, 0);
1764  f->SetParameter(2, 0);
1765 
1766  if ( s==0 ) { std::cerr << "Resplot::FitNull95() histogram s = " << s << std::endl; return 0; }
1767  if ( s->GetXaxis()==0) { std::cerr << "Resplot::FitNull95() histogram s->GetXaxis() not definied for histogram " << s->GetName() << std::endl; return 0; }
1768 
1769  // std::cout << "FitNull " << s->GetName() << " " << generate::GetEntries(s) << std::endl;
1770 
1771  if ( generate::GetEntries(s)<20 ) return f;
1772  // if ( s->GetEffectiveEntries<20 ) return f;
1773 
1774  s->GetXaxis()->SetRange(1,s->GetNbinsX());
1775  // s->GetXaxis()->UnZoom();
1776 
1777 
1780 
1783 
1784  int imax = findmax(s);
1785 
1786  f->FixParameter(0, s->GetMaximum()); f->SetParError(0, sqrt(s->GetMaximum()));
1787  // f->FixParameter(1, s->GetMean()); f->SetParError(1, s->GetMeanError());
1788 
1789 
1790  double entries = s->GetEntries() + s->GetBinContent(0) + s->GetBinContent(s->GetNbinsX()+1);
1791 
1792  // double sumx = s->GetBinContent(imax)*s->GetBinCenter(imax);
1793  // double sumx2 = s->GetBinContent(imax)*s->GetBinCenter(imax)*s->GetBinCenter(imax);
1794  double sumn = s->GetBinContent(imax);
1795 
1796  double frac = 0.95;
1797 
1798  int upperbin = imax;
1799  int lowerbin = imax;
1800 
1801  double uppersum = 0;
1802  double lowersum = 0;
1803 
1804  if ( entries>0 ) {
1805 
1806  int i=1;
1807  while ( true ) {
1808 
1809  const int upperbin_i = imax+i;
1810  const int lowerbin_i = imax-i;
1811 
1812  if ( upperbin_i>s->GetNbinsX() || lowerbin_i<1 ) break;
1813 
1814  double tsum = sumn + s->GetBinContent(upperbin_i) + s->GetBinContent(lowerbin_i);
1815 
1816  if ( tsum>entries*frac ) { uppersum = tsum/entries; break; }
1817 
1818  lowersum = tsum/entries;
1819 
1820  sumn = tsum;
1821 
1822  upperbin = upperbin_i;
1823  lowerbin = lowerbin_i;
1824 
1825  i++;
1826  }
1827  }
1828 
1829  // std::cout << "FitNull95 " << s->GetName() << " " << upperbin << " " << lowerbin << std::endl;
1830 
1831  // std::cout << "full rms " << s->GetRMS() << " +- " << s->GetRMSError() << "\t n=" << s->GetNbinsX()
1832  // << "\tfrac " << ( entries>0 ? sumn/entries : 0 ) << " - " << uppersum
1833  // << std::endl;
1834 
1835  // std::cout << "upper " << upperbin << " " << lowerbin << " " << interpolate_flag << std::endl;
1836 
1837  if ( upperbin!=lowerbin ) {
1838 
1839  double llim = s->GetBinLowEdge(1);
1840  double ulim = s->GetBinLowEdge(s->GetNbinsX()+1);
1841 
1842  double rms = 0;
1843  double erms = 0;
1844 
1845  std::vector<double> stats;
1846 
1847  if ( s_interpolate_flag ) {
1848  s->GetXaxis()->SetRange(lowerbin, upperbin);
1849 
1850 
1851  double rlower;
1852  double erlower;
1853 
1854  double rupper;
1855  double erupper;
1856 
1857 
1859 
1860  GetStats( s, stats );
1861  rlower = stats[2];
1862  erlower = stats[3];
1863 
1864  // std::cout << "FitNULL95() " << s->GetName() << "\t lower rms " << rlower << " +- " << s->GetRMSError() << " (wrong)\t " << stats[2] << " +- " << stats[3] << " (correct)" << std::endl;
1865  // std::cout << s->GetName() << " " << s->GetTitle() << "\t lower 95% rms " << rlower << " +- " << erlower << "\t( " << lowerbin << " - " << upperbin << " )" << std::endl;
1866 
1867  s->GetXaxis()->SetRange(lowerbin-1, upperbin+1);
1868 
1869  GetStats( s, stats );
1870  rupper = stats[2];
1871  erupper = stats[3];
1872 
1873  // std::cout << "FitNULL95() " << s->GetName() << "\t upper rms " << rupper << " +- " << s->GetRMSError() << " (wrong)\t " << stats[2] << " +- " << stats[3] << " (correct)" << std::endl;
1874 
1875  // std::cout << s->GetName() << " " << s->GetTitle() << "\t upper 95% rms " << rupper << " +- " << erupper << "\t( " << lowerbin-1 << " - " << upperbin+1 << " )" << std::endl;
1876 
1877  rms = rlower + (0.95-lowersum)*(rupper-rlower)/(uppersum-lowersum);
1878  erms = erlower + (0.95-lowersum)*(erupper-erlower)/(uppersum-lowersum);
1879 
1880  }
1881  else {
1882  s->GetXaxis()->SetRange(lowerbin, upperbin);
1883  GetStats( s, stats );
1884  rms = stats[2];
1885  erms = stats[3];
1886 
1887  }
1888 
1889 
1890  // f->FixParameter(2, s->GetRMS());
1891  // f->SetParError(2, s->GetRMSError());
1892  f->FixParameter(2, rms);
1893  f->SetParError(2, erms);
1894  f->FixParameter(1, s->GetMean());
1895  f->SetParError(1, s->GetMeanError());
1896  // std::cout << s->GetName() << " " << s->GetTitle() << "\t 95% inter rms " << rms << " +- " << erms << std::endl;
1897  s->GetXaxis()->SetRangeUser(llim, ulim);
1898  }
1899  else {
1900  f->FixParameter(1, 0);
1901  f->SetParError(1, 0);
1902  f->FixParameter(2, 0);
1903  f->SetParError(2, 0);
1904  // std::cout << " 95% rms " << 0 << " +- " << 0 << "\t( " << lowerbin << " - " << upperbin << " )" << std::endl;
1905  }
1906 
1907  s->GetXaxis()->SetRange(1, s->GetNbinsX());
1908 
1909  return f;
1910 }
1911 
1912 
1913 
1914 
1915 
1916 
1917 
1918 
1919 
1920 
1921 
1922 
1923 
1924 // or fit a Poissonian
1925 
1926 TF1* Resplot::FitPoisson(TH1D* s, double a, double b) {
1927 
1928  char fstring[1024];
1929  sprintf(fstring, "sqrt([0]*[0])*pow(sqrt([1]*[1]*([2]*[2])),x*[2]-[3])*exp(-sqrt([1]*[1]*([2]*[2])))/TMath::Gamma((x*[2]-[3])+1)");
1930  TF1 *f1 = new TF1("poisson", fstring );
1931  // TF1* f1 = new TF1("poiss","[0]*TMath::Poisson(x*[2],[1])");
1932 
1933  f1->SetParName(0, "Maximum");
1934  f1->SetParName(1, "Mean");
1935  f1->SetParName(2, "Scale");
1936  f1->SetParName(3, "Offset");
1937 
1938  f1->SetNpx(5000);
1939  f1->SetLineWidth(1);
1940  f1->SetLineColor(s->GetLineColor());
1941 
1942  int nbins=0;
1943  for (int j=1 ; j<=s->GetNbinsX() ; j++) if (s->GetBinContent(j)>0) nbins++;
1944 
1945  f1->SetParameter(0,s->GetBinContent(s->GetMaximumBin()) );
1946  f1->SetParameter(1,s->GetBinCenter(s->GetMaximumBin()));
1947  f1->SetParameter(2,1);
1948  f1->FixParameter(3,0);
1949 
1950  if (nbins>2) {
1951  // std::cout << "fitting " << f1 << std::endl;
1952  // if ( a==-999 || b==-999 ) s->Fit(f1,"Q");
1953  // else s->Fit(f1,"Q", "", a, b);
1954  if ( a==-999 || b==-999 ) s->Fit(f1);
1955  else s->Fit(f1,"", "", a, b);
1956  }
1957  else {
1958  for ( int j=0 ; j<4 ; j++ ) f1->SetParameter(j,0);
1959  }
1960 
1961  return f1;
1962 }
1963 
1964 
1965 
1966 TF1* Resplot::FitXExp(TH1D* s, double a, double b) {
1967 
1968  // TF1* f1 = new TF1("breitr", "[0]*[2]*[2]/([2]*[2]+(x-[1])*(x-[1]))", -100, 100 );
1969  // TF1* f1 = new TF1("xexpr", "[0]*(pow([2]*x,[1])*exp(-[2]*x))+[3]");
1970  // TF1* f1 = new TF1("xexpr", "[0]*(pow((1.0/[2])*(x-[3]),[1])*exp(-(1.0/[2])*(x-[3])))");
1971  TF1* f1 = new TF1("xexpr", "[0]*(pow((1.0/[2])*x,[1])*exp(-(1.0/[2])*x))+[3]");
1972 
1973  f1->SetParName(0, "Normalisation");
1974  f1->SetParName(1, "Power");
1975  f1->SetParName(2, "Lifetime");
1976  f1->SetParName(3, "offset");
1977 
1978  // f1->SetParameter(0, 1e-156);
1979  f1->SetParameter(0, s->GetBinContent(s->GetMaximumBin()) );
1980  // f1->SetParameter(1, -1*(s->GetBinCenter(s->GetMaximumBin())) );
1981  f1->SetParameter(2, s->GetRMS() );
1982  // f1->SetParameter(1, 0);
1983  f1->SetParameter(1, 0);
1984  f1->SetParameter(3, 0);
1985 
1986  // f1->FixParameter(0, 1);
1987  // f1->FixParameter(1, 25);
1988  // f1->FixParameter(1, 0);
1989  // f1->FixParameter(3, 0);
1990 
1991  f1->SetNpx(5000);
1992  f1->SetLineWidth(1);
1993 
1994  int nbins=0;
1995  for (int j=1 ; j<=s->GetNbinsX() ; j++) if (s->GetBinContent(j)>0) nbins++;
1996 
1997  // std::cout << "starting exp fit p1=" << f1->GetParameter(1) << std::endl;
1998 
1999  if (nbins>2) {
2000  // std::cout << "fitting " << f1 << std::endl;
2001  if ( a==-999 || b==-999 ) s->Fit(f1, "Q");
2002  else s->Fit(f1, "Q", "", a, b);
2003  }
2004  else {
2005  for ( int j=0 ; j<3 ; j++ ) f1->SetParameter(j,0);
2006  }
2007 
2008  // std::cout << "done" << std::endl;
2009 
2010  return f1;
2011 }
2012 
2013 
2014 
2015 
2016 
2017 
2018 Double_t langaufun(Double_t *x_par, Double_t *par) {
2019 
2020  //Fit parameters:
2021  //par[0]=Total area (integral -inf to inf, normalization constant)
2022  //par[1]=Most Probable (MP, location) parameter of Landau density
2023  //par[2]=Width (scale) parameter of Landau density
2024  //par[3]=Width (sigma) of convoluted Gaussian function
2025  //
2026  //In the Landau distribution (represented by the CERNLIB approximation),
2027  //the maximum is located at x=-0.22278298 with the location parameter=0.
2028  //This shift is corrected within this function, so that the actual
2029  //maximum is identical to the MP parameter.
2030 
2031  // Numeric constants
2032  Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
2033  Double_t mpshift = -0.22278298; // Landau maximum location
2034 
2035  // Control constants
2036  Double_t np = 500; // number of convolution steps
2037  Double_t sc = 5; // convolution extends to +-sc Gaussian sigmas
2038 
2039  double& x = x_par[0];
2040 
2041  // Variables
2042  Double_t xx;
2043  Double_t mpc;
2044  Double_t fland;
2045  Double_t sum = 0.0;
2046  Double_t xlow,xupp;
2047  Double_t step;
2048  Double_t i;
2049 
2050  // if ( par[2]<0.1*par[3] ) return 0;
2051 
2052  // correct the most probable bin position as described above
2053  mpc = par[1] - mpshift*par[2];
2054 
2055  // Range of convolution integral
2056 
2057  // if ( par[3]<par[2] ) np = np*par[2]/par[3];
2058 
2059  xlow = x - sc * par[3];
2060  xupp = x + sc * par[3];
2061 
2062  step = (xupp-xlow) / np;
2063 
2064  // Convolution integral of Landau and Gaussian by sum
2065  for(i=1.0; i<=np/2; i++) {
2066  xx = xlow + (i-0.5)*step;
2067  fland = TMath::Landau(xx,mpc,par[2]) / par[2];
2068  sum += fland * TMath::Gaus(x,xx,par[3]);
2069 
2070  xx = xupp - (i-0.5)*step;
2071  fland = TMath::Landau(xx,mpc,par[2]) / par[2];
2072  sum += fland * TMath::Gaus(x,xx,par[3]);
2073  }
2074 
2075  return (par[0] * step * sum * invsq2pi / par[3]);
2076 }
2077 
2078 
2079 
2080 // fit Gaussian to histogram
2081 
2082 TF1* Resplot::FitLandauGauss(TH1D* s, double a, double b) {
2083  // TF1* f1 = new TF1("landgaus", "landgaus");
2084  // root is STUPID and you have to specify the FIT RANGE if you
2085  // want to specify the number of fit parameters!!!
2086  TF1* f1 = new TF1("landgaus", langaufun, 0, 1, 4);
2087 
2088  f1->SetNpx(5000);
2089  f1->SetLineWidth(1);
2090 
2091  // parameter starting values
2092  double start_vals[4]; // = { 0.1, 3.0, 1.0, 0.3 };
2093 
2094  start_vals[0] = s->GetEntries();
2095  start_vals[1] = s->GetBinCenter(s->GetMaximumBin());
2096  start_vals[2] = s->GetRMS();
2097  start_vals[3] = 0.5*start_vals[2];
2098 
2099  f1->SetParNames("Area", "MP", "Width", "GSigma");
2100  f1->SetParameters(start_vals);
2101 
2102  // parameter low and high values for the fit
2103  // double par_lo[4] = { 0.0, -20.0, 0.0, 0.00 };
2104  // double par_hi[4] = { 20000.0, 50000.0, 50.0, 0.01 };
2105  // for ( int i=0 ; i<4 ; i++ ) f1->SetParLimits(i, par_lo[i], par_hi[i]);
2106  f1->SetParLimits( 2, 0.1*start_vals[2], 10*start_vals[2] );
2107  f1->SetParLimits( 3, 0.0, 10*start_vals[2] );
2108 
2109  int nbins=0;
2110  for (int j=1 ; j<=s->GetNbinsX() ; j++) if (s->GetBinContent(j)) nbins++;
2111 
2112  if (nbins>2) {
2113  if ( a==-999 || b==-999 ) s->Fit(f1,"Q");
2114  else s->Fit(f1,"Q", "", a, b);
2115  // if ( a==-999 || b==-999 ) s->Fit(f1);
2116  // else s->Fit(f1,"", "", a, b);
2117  }
2118  else for ( int j=0 ; j<3 ; j++ ) f1->SetParameter(j, 0);
2119 
2120  return f1;
2121 }
2122 
2123 
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
FindMean
double FindMean(TH1D *s, double frac=0.95)
Definition: Resplot.cxx:1309
Resplot::AddDirectoryStatus
static bool AddDirectoryStatus()
Definition: Resplot.h:456
nexperiments_max
int nexperiments_max
Definition: Resplot.cxx:1682
pdg_comparison.sigma
sigma
Definition: pdg_comparison.py:324
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
Resplot::ERROR
ERROR
Definition: Resplot.h:59
mean
void mean(std::vector< double > &bins, std::vector< double > &values, const std::vector< std::string > &files, const std::string &histname, const std::string &tplotname, const std::string &label="")
Definition: dependence.cxx:254
Resplot::m_name
std::string m_name
Definition: Resplot.h:640
Resplot::m_mean
TH1D * m_mean
Definition: Resplot.h:643
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
langaufun
Double_t langaufun(Double_t *x_par, Double_t *par)
Definition: Resplot.cxx:2018
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
Resplot::DWrite
Int_t DWrite(TDirectory *g=0) const
boooo, need to use the stupid Int_t class because I foolishly gave this function the same name as TOb...
Definition: Resplot.cxx:342
Resplot::SetDirectory
void SetDirectory(TDirectory *=0)
Definition: Resplot.h:459
dqt_zlumi_pandas.hname
string hname
Definition: dqt_zlumi_pandas.py:279
Resplot::GetEfficiencies
TH1D * GetEfficiencies(const std::string &hname, double lower, double upper)
Definition: Resplot.cxx:958
hist_file_dump.d
d
Definition: hist_file_dump.py:137
Resplot::m_finalised
bool m_finalised
Definition: Resplot.h:667
ZeroErrors
void ZeroErrors(TH1D *h)
Definition: Resplot.cxx:67
number
std::string number(const double &x)
Definition: Resplot.cxx:373
Resplot::m_Set
bool m_Set
Definition: Resplot.h:638
Resplot::s_scalerms95
static bool s_scalerms95
Definition: Resplot.h:675
python.compressB64.sx
string sx
Definition: compressB64.py:96
Directory
Definition: Directory.h:29
skel.it
it
Definition: skel.GENtoEVGEN.py:396
PixelAthClusterMonAlgCfg.ybins
ybins
Definition: PixelAthClusterMonAlgCfg.py:169
Resplot.h
Resplot::FitXExp
static TF1 * FitXExp(TH1D *s, double a=-999, double b=-999)
Definition: Resplot.cxx:1966
generate::experiment
given a histogram, get a generator for the distribution in that histogram, then run some pseudo-exper...
Definition: generate.h:132
Resplot::combine
static TH2D * combine(const TH2 *h, double x, int N)
combine bins along the x axis after value x, combine each N bins into a single bin
Definition: Resplot.cxx:606
Resplot::m_h1d
TH1D * m_h1d
Definition: Resplot.h:648
Resplot::m_Nentries
TH1D * m_Nentries
Definition: Resplot.h:642
Resplot::m_sigma
TH1D * m_sigma
Definition: Resplot.h:644
upper
int upper(int c)
Definition: LArBadChannelParser.cxx:49
JetTiledMap::N
@ N
Definition: TiledEtaPhiMap.h:44
NCounter
double NCounter(TH1D *h)
Definition: Resplot.cxx:380
Resplot::FitPoisson
static TF1 * FitPoisson(TH1D *s, double a=-999, double b=-999)
Definition: Resplot.cxx:1926
PlotPulseshapeFromCool.np
np
Definition: PlotPulseshapeFromCool.py:64
yodamerge_tmp.scale
scale
Definition: yodamerge_tmp.py:138
trigbs_dumpHLTContentInBS.stats
stats
Definition: trigbs_dumpHLTContentInBS.py:91
x
#define x
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
read_hist_ntuple.f2
f2
Definition: read_hist_ntuple.py:20
generate.h
unZeroErrors
void unZeroErrors(TH1D *h)
Definition: Resplot.cxx:73
Resplot::FitCentralGaussian
static TF1 * FitCentralGaussian(TH1D *s, double a=-999, double b=-999)
Definition: Resplot.cxx:1080
ClassImp
ClassImp(xAOD::TFileChecker) namespace xAOD
Definition: TFileChecker.cxx:28
Resplot::getbins
static std::vector< double > getbins(const TH1 *h)
Definition: Resplot.cxx:921
Resplot::SetSecondary
void SetSecondary(int n, double a, double b)
Definition: Resplot.h:594
Resplot::FitInternal
static TF1 * FitInternal(TH1D *s, double a=-999, double b=-999, TF1 *f1=0)
Definition: Resplot.cxx:1051
Resplot::Initialise
void Initialise(const std::string &name, int n1, double a1, double b1, int n2, double a2, double b2)
Definition: Resplot.cxx:194
fitman.hy
def hy
Definition: fitman.py:515
GetStats
void GetStats(TH1D *h, std::vector< double > &stats)
Definition: Resplot.cxx:113
Resplot::m_uniform
bool m_uniform
Definition: Resplot.h:670
Resplot::FitNull95Obsolete
static TF1 * FitNull95Obsolete(TH1D *s, double frac=0.95, bool useroot=false)
Definition: Resplot.cxx:1408
defineDB.tagname
string tagname
Definition: JetTagCalibration/share/defineDB.py:19
convertTimingResiduals.sum
sum
Definition: convertTimingResiduals.py:55
ParseInputs.gDirectory
gDirectory
Definition: Final2012/ParseInputs.py:133
lumiFormat.i
int i
Definition: lumiFormat.py:85
beamspotman.n
n
Definition: beamspotman.py:731
python.CaloCondTools.g
g
Definition: CaloCondTools.py:15
Resplot::FitNull
static TF1 * FitNull(TH1D *s, double a=-999, double b=-999)
Definition: Resplot.cxx:1288
Resplot::FitRGaussian
static TF1 * FitRGaussian(TH1D *s, double a=-999, double b=-999)
Definition: Resplot.cxx:1031
checkxAOD.frac
frac
Definition: Tools/PyUtils/bin/checkxAOD.py:259
Resplot::FitNull95
static TF1 * FitNull95(TH1D *s, double a=0, double b=0)
Definition: Resplot.cxx:1673
Resplot::FitNull95New
static TF1 * FitNull95New(TH1D *s, double frac=0.95, bool useroot=false)
Definition: Resplot.cxx:1687
Resplot::m_slices
std::vector< TH1D * > m_slices
Definition: Resplot.h:664
Resplot::FitATan
static TF1 * FitATan(TH1D *s, double a=-999, double b=-999)
Definition: Resplot.cxx:1220
Resplot::FitNull95Central
static TF1 * FitNull95Central(TH1D *s)
Definition: Resplot.cxx:1754
Resplot::s_oldrms95
static bool s_oldrms95
temporarily allow toggeling between old and new rms95 error estimates
Definition: Resplot.h:674
Resplot::findbins
static std::vector< double > findbins(const TH1 *h, double epsilon)
Definition: Resplot.cxx:838
Resplot::Name
const std::string & Name() const
Definition: Resplot.h:340
hist_file_dump.f
f
Definition: hist_file_dump.py:135
plotting.yearwise_luminosity_vs_mu.bins
bins
Definition: yearwise_luminosity_vs_mu.py:30
Resplot
Definition: Resplot.h:50
Resplot::GetEfficiency
StatVal GetEfficiency(TH1D *h, double lower, double upper)
Definition: Resplot.cxx:932
Resplot::FitGaussian
static TF1 * FitGaussian(TH1D *s, double a=-999, double b=-999)
Definition: Resplot.cxx:1025
Resplot::rotate
static TH2D * rotate(const TH2 *h)
flip the x and y axes
Definition: Resplot.cxx:564
imax
int imax(int i, int j)
Definition: TileLaserTimingTool.cxx:33
findmax
int findmax(TH1D *s)
Definition: Resplot.cxx:1396
Resplot::m_yaxis
std::string m_yaxis
Definition: Resplot.h:662
Resplot::Write
Int_t Write(const char *=0, Int_t=0, Int_t=0) const
Hooray, this stupidity is to overwride both the const and non-const TObject Write methods Fixme: shou...
Definition: Resplot.h:452
Resplot::OK
@ OK
Definition: Resplot.h:59
Resplot::m_g_mean
StatVal m_g_mean
Definition: Resplot.h:650
tagname.h
getWeights
double getWeights(TH1D *h)
Definition: Resplot.cxx:79
Resplot::FitRBreit
static TF1 * FitRBreit(TH1D *s, double a=-999, double b=-999)
Definition: Resplot.cxx:1036
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:29
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
Resplot::skip
void skip(TH1D *t)
Definition: Resplot.h:583
Resplot::m_n_primary
int m_n_primary
Definition: Resplot.h:655
Resplot::m_h2d
TH2D * m_h2d
Definition: Resplot.h:647
Resplot::Finalise
int Finalise(double a=-999, double b=-999, TF1 *(*func)(TH1D *s, double a, double b)=Resplot::FitGaussian)
Definition: Resplot.cxx:388
SCT_CalibAlgs::nbins
@ nbins
Definition: SCT_CalibNumbers.h:10
Resplot::BINWIDTH
@ BINWIDTH
Definition: Resplot.h:59
LArCellBinning.xbins
int xbins
Definition: LArCellBinning.py:163
Resplot::s_interpolate_flag
static bool s_interpolate_flag
Definition: Resplot.h:634
Rtt_histogram.n1
n1
Definition: Rtt_histogram.py:21
nexperiments_
int nexperiments_
Definition: Resplot.cxx:1681
cwd
std::string cwd
Definition: listroot.cxx:38
python.PyAthena.v
v
Definition: PyAthena.py:154
Resplot::FitLandauGauss
static TF1 * FitLandauGauss(TH1D *s, double a=-999, double b=-999)
Definition: Resplot.cxx:2082
ReadBchFromCool.good
good
Definition: ReadBchFromCool.py:433
fitman.sx2
sx2
Definition: fitman.py:537
a
TList * a
Definition: liststreamerinfos.cxx:10
h
Resplot::clear
void clear()
Definition: Resplot.h:623
StatVal
simple struct to hold a value and it's associated uncertainty.
Definition: StatVal.h:26
Resplot::m_xaxis
std::string m_xaxis
Definition: Resplot.h:661
Resplot::m_fitname
std::string m_fitname
Definition: Resplot.h:665
Resplot::~Resplot
virtual ~Resplot()
Definition: Resplot.cxx:174
generate::GetEntries
double GetEntries(TH1D *h, int ilow, int ihi)
Definition: rmsFrac.cxx:20
Resplot::s_ErrorSet
static ERROR s_ErrorSet
Definition: Resplot.h:680
beamspotnt.rms
rms
Definition: bin/beamspotnt.py:1266
fitman.hx
def hx
Definition: fitman.py:514
entries
double entries
Definition: listroot.cxx:49
tagname
Definition: tagname.h:29
LArCellBinning.step
step
Definition: LArCellBinning.py:158
Resplot::m_chi2
TH1D * m_chi2
Definition: Resplot.h:645
Directory.h
nexperiments_min
int nexperiments_min
Definition: Resplot.cxx:1683
Resplot::s_nofit
static bool s_nofit
Definition: Resplot.h:677
Resplot::delskip
void delskip(TH1D *t)
Definition: Resplot.h:584
Resplot::s_mAddDirectoryStatus
static bool s_mAddDirectoryStatus
Definition: Resplot.h:633
Resplot::SetPrimary
void SetPrimary(int n, const double *)
Definition: Resplot.h:592
Resplot::FitBreit
static TF1 * FitBreit(TH1D *s, double a=-999, double b=-999)
Definition: Resplot.cxx:1159
Resplot::s_rversion
static const std::string s_rversion
Definition: Resplot.h:636
Resplot::rebin
static TH1D * rebin(const TH1 *h, const std::vector< double > &bins)
Definition: Resplot.cxx:891
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
Resplot::HISTOWIDTH
@ HISTOWIDTH
Definition: Resplot.h:59
binwidth
void binwidth(TH1D *h)
Definition: Resplot.cxx:50
test_pyathena.counter
counter
Definition: test_pyathena.py:15
CaloNoise_fillDB.mu
mu
Definition: CaloNoise_fillDB.py:53
fitman.k
k
Definition: fitman.py:528
Resplot::m_g_sigma
StatVal m_g_sigma
Definition: Resplot.h:651
read_hist_ntuple.f1
f1
Definition: read_hist_ntuple.py:4
rmsFrac.h
Generated on Sat Jan 25 2025 21:16:59 for ATLAS Offline Software by   doxygen 1.8.18