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