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