hpstr
The Heavy Photon Search Toolkit for Reconstruction (hpstr) provides an interface to physics data from the HPS experiment saved in the LCIO format and converts it into an ROOT based format.
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
ZBiHistos.cxx
Go to the documentation of this file.
1#include "ZBiHistos.h"
2#include <sstream>
3#include "TMath.h"
4
7
10
11ZBiHistos::ZBiHistos(const std::string& inputName){
12 m_name = inputName;
13}
14
15void ZBiHistos::addHisto1d(std::string histoname, std::string xtitle, int nbinsX, float xmin, float xmax){
16 histos1d[m_name+"_"+histoname] = plot1D(m_name+"_"+histoname, xtitle, nbinsX, xmin, xmax);
17}
18
19void ZBiHistos::addHisto2d(std::string histoname, std::string xtitle, int nbinsX, float xmin, float xmax, std::string ytitle, int nbinsY, float ymin, float ymax){
20 histos2d[m_name+"_"+histoname] = plot2D(m_name+"_"+histoname, xtitle, nbinsX, xmin, xmax, ytitle, nbinsY, ymin, ymax);
21}
22
23void ZBiHistos::resetHistograms1d(){
24 for(it1d it=histos1d.begin(); it != histos1d.end(); it ++){
25 if(it->second != nullptr)
26 it->second->Reset();
27 }
28}
29
30void ZBiHistos::resetHistograms2d(){
31 for(it2d it=histos2d.begin(); it != histos2d.end(); it ++){
32 if(it->second != nullptr)
33 it->second->Reset();
34 }
35}
36
37double ZBiHistos::integrateHistogram1D(std::string histoname){
38 int xmax;
39 int xmin;
40 double integral;
41
42 //Check if histogram exists to integrate
43 if(!histos1d.count(histoname)){
44 std::cout << "[ZBiHistos}::ERROR::NO HISTOGRAM NAMED " << histoname <<
45 " FOUND! DEFINE HISTOGRAM IN HISTOCONFIG!" << std::endl;
46 integral = -9999.9;
47
48 }
49
50 if(histos1d[histoname] == nullptr){
51 std::cout << "ERROR: Histogram for " << histoname << " is NULLPTR" << std::endl;
52 return -9999.9;
53 }
54
55 else{
56 xmax = histos1d[histoname]->FindLastBinAbove(0.0);
57 xmin = histos1d[histoname]->FindFirstBinAbove(0.0);
58 //integral = histos1d[histoname]->Integral(0,histos1d[histoname]->GetNbinsX()+1);
59 integral = histos1d[histoname]->Integral(xmin,xmax);
60 }
61 return integral;
62}
63
64double ZBiHistos::cutFractionOfSignalVariable(std::string cutvariable, bool isCutGreaterThan, double cutFraction, double initialIntegral){
65
66 TH1F* histo = histos1d[m_name+"_"+cutvariable+"_h"];
67 int xmax = histo->FindLastBinAbove(0.0);
68 int xmin = histo->FindFirstBinAbove(0.0);
69
70 if(debug_){
71 std::cout << "Initial Integral: " << initialIntegral << std::endl;
72 std::cout << "Cut variable: " << cutvariable << std::endl;
73 std::cout << "Cut fraction: " << cutFraction << std::endl;
74 }
75
76 double cutvalue;
77 if(isCutGreaterThan){
78 cutvalue = histo->GetXaxis()->GetBinLowEdge(xmin);
79 while(histo->Integral(xmin,xmax) > initialIntegral*(1.0-cutFraction)){
80 xmin = xmin + 1;
81 cutvalue = histo->GetXaxis()->GetBinLowEdge(xmin);
82 }
83 }
84 else {
85 cutvalue = histo->GetXaxis()->GetBinUpEdge(xmax);
86 while(histo->Integral(xmin,xmax) > initialIntegral*(1.0-cutFraction)){
87 xmax = xmax - 1;
88 cutvalue = histo->GetXaxis()->GetBinUpEdge(xmax);
89 }
90 }
91
92 return cutvalue;
93}
94
95void ZBiHistos::defineZBiCutflowProcessorHistograms(){
96 //These histograms are used specifically for tracking the ZBitCutflowProcessor iterative process
97 addHisto2d("persistent_cuts_hh","pct_sig_cut", 1000,-0.5,99.5,"cut_id",50,0.5,50.5);
98 addHisto2d("test_cuts_ZBi_hh","pct_sig_cut", 1000,-0.5,99.5,"cut_id",50,0.5,50.5);
99 addHisto2d("test_cuts_values_hh","pct_sig_cut", 1000,-0.5,99.5,"cut_id",50,0.5,50.5);
100 addHisto2d("test_cuts_zcut_hh","pct_sig_cut", 1000,-0.5,99.5,"zcut",50,0.5,50.5);
101 addHisto2d("test_cuts_nbkg_hh","pct_sig_cut", 1000,-0.5,99.5,"zcut",50,0.5,50.5);
102 addHisto2d("test_cuts_nsig_hh","pct_sig_cut", 1000,-0.5,99.5,"zcut",50,0.5,50.5);
103 addHisto2d("best_test_cut_ZBi_hh","pct_sig_cut", 1000,-0.5,99.5,"ZBi",2000,0.0,20.0);
104
105 addHisto1d("best_test_cut_ZBi_h","pct_sig_cut", 1000,-0.5,99.5);
106 addHisto1d("best_test_cut_zcut_h","pct_sig_cut", 1000,-0.5,99.5);
107 addHisto1d("best_test_cut_nsig_h","pct_sig_cut", 1000,-0.5,99.5);
108 addHisto1d("best_test_cut_nbkg_h","pct_sig_cut", 1000,-0.5,99.5);
109 addHisto1d("best_test_cut_id_h","pct_sig_cut", 1000,-0.5,99.5);
110}
111
112void ZBiHistos::set2DHistoYlabel(std::string histoName, int ybin, std::string ylabel){
113 histos2d[m_name+"_"+histoName]->GetYaxis()->SetBinLabel(ybin, ylabel.c_str());
114}
115
116std::vector<double> ZBiHistos::defineImpactParameterCut(double alpha){
117 //alpha defines the % of signal we allow to be cut. Default is 15%
118 TH2F* hh = (TH2F*)histos2d[m_name+"_z0_v_recon_z_hh"];
119 addHisto1d("impact_parameter_up_h","recon_z [mm]", 450, -20.0, 70.0);
120 addHisto1d("impact_parameter_down_h","recon_z [mm]", 450, -20.0, 70.0);
121 TH1F* up_h = (TH1F*)histos1d[m_name+"_impact_parameter_up_h"];
122 TH1F* down_h = (TH1F*)histos1d[m_name+"_impact_parameter_down_h"];
123
124 for(int i=0; i < hh->GetNbinsX(); i++){
125
126 TH1F* projy = (TH1F*)hh->ProjectionY(("projy_bin_"+std::to_string(i+1)).c_str(),i+1,i+1);
127 if(projy->GetEntries() < 1) continue;
128
129 //Impact param for z0 > 0
130 int start_bin = projy->FindBin(0.0);
131 int end_bin = projy->FindLastBinAbove(0.0);
132 double refIntegral = projy->Integral(start_bin,end_bin);
133 int cutz0_bin = start_bin;
134 double testIntegral = refIntegral;
135 while(testIntegral > (1.0-alpha)*refIntegral && cutz0_bin < end_bin-1){
136 cutz0_bin = cutz0_bin + 1;
137 testIntegral = projy->Integral(cutz0_bin, end_bin);
138 }
139
140 double cutz0_up = projy->GetXaxis()->GetBinLowEdge(cutz0_bin);
141 for(int j=0; j<(int)refIntegral; j++){
142 up_h->Fill(hh->GetXaxis()->GetBinCenter(i+1),cutz0_up/refIntegral);
143 }
144
145 //impact param for z0 < 0
146 end_bin = projy->FindFirstBinAbove(0.0);
147 start_bin = projy->FindBin(0.0);
148 refIntegral = projy->Integral(end_bin, start_bin);
149 cutz0_bin = start_bin;
150 testIntegral = refIntegral;
151 while(testIntegral > (1.0-alpha)*refIntegral && cutz0_bin > end_bin+1){
152 cutz0_bin = cutz0_bin - 1;
153 testIntegral = projy->Integral(end_bin, cutz0_bin);
154 }
155
156 double cutz0_down = projy->GetXaxis()->GetBinUpEdge(cutz0_bin);
157 for(int j=0; j<(int)refIntegral; j++){
158 down_h->Fill(hh->GetXaxis()->GetBinCenter(i+1),cutz0_down/refIntegral);
159 }
160 }
161
162 TF1* fitFunc = new TF1("linear_fit","[0]*(x-[1])",5.0,70.0);
163 TFitResultPtr fitResult = (TFitResultPtr)up_h->Fit("linear_fit","QS","",5.0,70.0);
164 fitFunc->Draw();
165 double m_p = fitResult->GetParams()[0];
166 double a_p = fitResult->GetParams()[1];
167
168 fitResult = (TFitResultPtr)down_h->Fit("linear_fit","QS","",5.0,70.0);
169 fitFunc->Draw();
170 double m_d = fitResult->GetParams()[0];
171 double a_d = fitResult->GetParams()[1];
172
173 //Find the location in z0 where two lines meet
174 double diff = 999.9;
175 double x = 10.0;
176 while(std::abs((m_p*(x-a_p) - (m_d*(x-a_d)))) < diff ){
177 //while( std::abs((a_p+b_p*x) - (a_d+b_d*x)) < diff ){
178 diff = std::abs((m_p*(x-a_p)) - (m_d*(x-a_d)));
179 x = x - 0.01;
180 }
181 double beta = .5*((m_p*(x-a_p)) + (m_p*(x-a_d)));
182 std::vector<double> params {m_p,a_p,m_d,a_d, beta, x};
183
184 delete fitFunc;
185 return params;
186}
187
188/*
189TF1* ZBiHistos::fitExponentialPlusConst(std::string histogramName, double start_nevents){
190
191 //Get z vertex distribution corresponding to Test Cut
192 std::string histoname = m_name+"_"+histogramName+"_h";
193
194 //If histogram is empty, return nullptr
195 if(histos1d[histoname]->GetEntries() < 1){
196 std::cout << "[ZBiHistos]::WARNING: Background Model is NULL: " <<
197 histogramName << " distribution is empty and could not be fit!" << std::endl;
198 return nullptr;
199 }
200 //Start fit where start_nevents are left in tail
201 int lastbin = histos1d[histoname]->FindLastBinAbove(0.0);
202 int firstbin = lastbin - 1;
203 double test_integral = 0.0;
204 while(test_integral < start_nevents && histos1d[histoname]->GetBinLowEdge(firstbin) > 0.0){
205 test_integral = histos1d[histoname]->Integral(firstbin, lastbin);
206 firstbin = firstbin - 1;
207 }
208
209 double xmin = histos1d[histoname]->GetBinLowEdge(firstbin);
210 double xmax = histos1d[histoname]->GetBinLowEdge(lastbin+1);
211 histos1d[histoname]->GetXaxis()->SetRangeUser(xmin,xmax);
212 std::cout << "bkg xmin: " << xmin << std::endl;
213 std::cout << "bkg xmax: " << xmax << std::endl;
214
215 //Initially fit tail with exponential to seed params
216 double best_seed_chi2 = 9999999.9;
217 double seed_0;
218 double seed_1;
219 TF1* fitFunc_seed = new TF1("fitfunc_seed", "[0]*exp([1]*x)", xmin, xmax);
220 fitFunc_seed->SetParameter(0, (double)start_nevents);
221 fitFunc_seed->SetParameter(1, -0.5);
222 TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc_seed, "QLSIM", "", xmin,xmax);
223 seed_0 = fitResult->Parameter(0);
224 seed_1 = fitResult->Parameter(1);
225 if(fitResult->Ndf() <= 0){
226 best_seed_chi2 = fitResult->Chi2()/fitResult->Ndf();
227 }
228
229
230 std::cout << "First Seed 0: " << seed_0 << std::endl;
231 std::cout << "First Seed 1: " << seed_1 << std::endl;
232 std::cout << "First Seed Chi2: " << best_seed_chi2 << std::endl;
233
234 TRandom3* ran = new TRandom3();
235 ran->SetSeed(0);
236 TF1* fitFunc_seed = new TF1("fitfunc_seed", "[0]*exp([1]*x)", xmin, xmax);
237 //fitFunc_seed->SetParameter(0, (double)start_nevents);
238 fitFunc_seed->SetParameter(0, 10.0);
239 fitFunc_seed->SetParameter(1, -0.5);
240 TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc_seed, "QSIM", "", xmin,xmax);
241 std::cout << "Start with seed 0:10 and 1:-0.5" << std::endl;
242
243 if(fitResult->Ndf() > 0){
244 best_seed_chi2 = fitResult->Chi2()/fitResult->Ndf();
245 seed_0 = fitResult->Parameter(0);
246 seed_1 = fitResult->Parameter(1);
247 }
248 std::cout << "First Seed 0: " << seed_0 << std::endl;
249 std::cout << "First Seed 1: " << seed_1 << std::endl;
250 std::cout << "First Seed Chi2: " << best_seed_chi2 << std::endl;
251 for(int i=0; i < 50; i++){
252 //fitFunc = new TF1("fitfunc_seed", "[0]*exp([1]*x)", xmin, xmax);
253 //fitFunc_seed->SetParameters(std::abs(ran->Gaus(10.0,2.0)), -std::abs(ran->Uniform(0.1,0.8)));
254 fitFunc_seed->SetParameters(std::abs(ran->Uniform(1.0,20.0)), -std::abs(ran->Uniform(0.1,0.8)));
255 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc_seed, "QSIM", "", xmin,xmax);
256
257 if(fitResult->Ndf() <= 0)
258 continue;
259
260 std::cout << "Iter Seed 0: " << fitResult->Parameter(0) << std::endl;
261 std::cout << "Iter Seed 1: " << fitResult->Parameter(1) << std::endl;
262 std::cout << "Iter Seed Chi2: " << fitResult->Chi2()/fitResult->Ndf() << std::endl;
263
264 if(fitResult->Chi2()/fitResult->Ndf() < best_seed_chi2){
265 best_seed_chi2 = fitResult->Chi2()/fitResult->Ndf();
266 seed_0 = fitResult->Parameter(0);
267 seed_1 = fitResult->Parameter(1);
268 std::cout << "Update Iter Seed 0: " << seed_0 << std::endl;
269 std::cout << "Update Iter Seed 1: " << seed_1 << std::endl;
270 std::cout << "Update Iter Seed Chi2: " << best_seed_chi2 << std::endl;
271 }
272 }
273 std::cout << "Seed 0: " << seed_0 << std::endl;
274 std::cout << "Seed 1: " << seed_1 << std::endl;
275 std::cout << "Seed Chi2: " << best_seed_chi2 << std::endl;
276
277 double best_chi2 = best_seed_chi2;
278 double param_2;
279 TF1* fitFunc = new TF1("fitfunc", " (x<[2])*([0]*exp([1]*x)) + (x>=[2])*([0]*exp([1]*[2]))", xmin, xmax);
280 fitFunc->FixParameter(0, seed_0);
281 fitFunc->FixParameter(1, seed_1);
282 for(int i=0; i < 1000; i++){
283 fitFunc->FixParameter(2, (double)i/10.0);
284 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax);
285
286 if(fitResult->Ndf() <= 0)
287 continue;
288
289 if(fitResult->Chi2()/fitResult->Ndf() < best_chi2){
290 best_chi2 = fitResult->Chi2()/fitResult->Ndf();
291 param_2 = fitResult->Parameter(2);
292 }
293 }
294
295 std::cout << "Final Seed 0: " << seed_0 << std::endl;
296 std::cout << "Final Seed 1: " << seed_1 << std::endl;
297 std::cout << "param 2: " << param_2 << std::endl;
298 std::cout << "Best Chi2: " << best_chi2 << std::endl;
299
300 delete ran;
301
302 if(best_chi2 < best_seed_chi2){
303 fitFunc->FixParameter(0, seed_0);
304 fitFunc->FixParameter(1, seed_1);
305 fitFunc->FixParameter(2, param_2);
306 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
307 fitFunc->Draw();
308 delete fitFunc_seed;
309 return fitFunc;
310 }
311 else{
312 fitFunc->FixParameter(0, seed_0);
313 fitFunc->FixParameter(1, seed_1);
314 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc_seed, "QLSIM", "", xmin,xmax);
315 fitFunc_seed->Draw();
316 delete fitFunc;
317 return fitFunc_seed;
318 }
319 if(fitResult->Chi2()/fitResult->Ndf() < best_seed_chi2){
320 delete fitFunc_seed;
321 fitFunc->Draw();
322 return fitFunc;
323 }
324 else{
325 delete fitFunc;
326 fitFunc_seed->Draw();
327 return fitFunc_seed;
328 }
329 */
330
331 /* Fails too often. Any failure breaks the code...not very robust
332 TRandom3* ran = new TRandom3();
333 ran->SetSeed(0);
334
335 //Initially fit tail with exponential to seed params
336 double best_chi2 = 9999999.9;
337 double seed_0;
338 double seed_1;
339 TF1* fitFunc = new TF1("fitfunc", "[0]*exp([1]*x)", xmin, xmax);
340 TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
341 seed_0 = fitResult->Parameter(0);
342 seed_1 = fitResult->Parameter(1);
343
344 for(int i=0; i < 30; i++){
345 fitFunc = new TF1("fitfunc", "[0]*exp([1]*x)", xmin, xmax);
346 fitFunc->SetParameters(ran->Gaus(seed_0,1), ran->Gaus(seed_1,1));
347 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
348
349 if(fitResult->Ndf() <= 0)
350 continue;
351
352 if(fitResult->Chi2()/fitResult->Ndf() < best_chi2){
353 best_chi2 = fitResult->Chi2()/fitResult->Ndf();
354 seed_0 = fitResult->Parameter(0);
355 seed_1 = fitResult->Parameter(1);
356 }
357 }
358
359 //Exponential plus constant
360 fitFunc = new TF1("fitfunc", " (x<[2])*([0]*exp([1]*x)) + (x>=[2])*([0]*exp([1]*[2]))", xmin, xmax);
361
362 best_chi2 = 9999999.9;
363 double best_0 = seed_0;
364 double best_1 = seed_1;
365 double best_2;
366
367 int iteration = 30;
368 fitFunc->FixParameter(0, seed_0);
369 fitFunc->FixParameter(1, seed_1);
370 for(int i=0; i < iteration; i++){
371 //fitFunc->SetParameter(0, best_0);
372 //fitFunc->SetParameter(1, best_1);
373 fitFunc->SetParameter(2, ran->Uniform(50.0));
374 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
375
376 if(fitResult->Ndf() <= 0)
377 continue;
378
379 if(fitResult->Chi2()/fitResult->Ndf() < best_chi2){
380 best_chi2 = fitResult->Chi2()/fitResult->Ndf();
381 //best_0 = fitResult->Parameter(0);
382 //best_1 = fitResult->Parameter(1);
383 best_2 = fitResult->Parameter(2);
384 }
385 }
386
387 //fitFunc->SetParameters(best_0, best_1, best_2);
388 //fitFunc->SetParameters(seed_0, seed_1, best_2);
389 fitFunc->FixParameter(0, seed_0);
390 fitFunc->FixParameter(1, seed_1);
391 fitFunc->FixParameter(2, best_2);
392 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
393 fitFunc->Draw();
394
395 return fitFunc;
396}*/
397
398TF1* ZBiHistos::fitExponentialPlusExp(std::string histogramName, double start_nevents){
399
400 //Get z vertex distribution corresponding to Test Cut
401 std::string histoname = m_name+"_"+histogramName+"_h";
402
403 //If histogram is empty, return nullptr
404 if(histos1d[histoname]->GetEntries() < 1){
405 std::cout << "[ZBiHistos]::WARNING: Background Model is NULL: " <<
406 histogramName << " distribution is empty and could not be fit!" << std::endl;
407 return nullptr;
408 }
409 double originalXmin = histos1d[histoname]->GetXaxis()->GetXmin();
410 double originalXmax = histos1d[histoname]->GetXaxis()->GetXmax();
411
412 //Start fit where start_nevents are left in tail
413 int lastbin = histos1d[histoname]->FindLastBinAbove(0.0);
414 if(histos1d[histoname]->Integral() < start_nevents)
415 start_nevents = histos1d[histoname]->GetBinContent(histos1d[histoname]->GetMaximumBin());
416 int firstbin = lastbin - 1;
417 double test_integral = 0.0;
418 while(test_integral < start_nevents && histos1d[histoname]->GetBinLowEdge(firstbin) > 0.0){
419 test_integral = histos1d[histoname]->Integral(firstbin, lastbin);
420 firstbin = firstbin - 1;
421 }
422
423 double xmin = histos1d[histoname]->GetBinLowEdge(firstbin);
424 double xmax = histos1d[histoname]->GetBinLowEdge(lastbin+1);
425 histos1d[histoname]->GetXaxis()->SetRangeUser(xmin,xmax);
426
427 //Initially fit tail with exponential to seed params
428 double best_seed_chi2 = 9999999.9;
429 double seed_0;
430 double seed_1;
431
432 TF1* fitFunc_seed = new TF1("fitFunc_seed", "[0]*exp([1]*x)", xmin, xmax);
433 //TF1* fitFunc_seed = new TF1("fitFunc_seed", "expo", xmin, xmax);
434 TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc_seed, "QSIM", "", xmin,xmax);
435 seed_0 = fitResult->Parameter(0);
436 seed_1 = fitResult->Parameter(1);
437 if(fitResult->Ndf() > 0){
438 best_seed_chi2 = fitResult->Chi2()/fitResult->Ndf();
439 }
440
441 TRandom3* ran = new TRandom3();
442 ran->SetSeed(0);
443
444 best_seed_chi2 = 999999.9;
445 for(int i=0; i < 50; i++){
446 fitFunc_seed->SetParameters(std::abs(ran->Gaus((double)start_nevents,3.0)), -std::abs(ran->Uniform(0.01,1.0)));
447 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc_seed, "QSIM", "", xmin,xmax);
448
449 if(fitResult->Ndf() <= 0)
450 continue;
451
452 if(fitResult->Chi2()/fitResult->Ndf() < best_seed_chi2){
453 if(fitResult->Parameter(0) < 1000.0)
454 continue;
455 best_seed_chi2 = fitResult->Chi2()/fitResult->Ndf();
456 seed_0 = fitResult->Parameter(0);
457 seed_1 = fitResult->Parameter(1);
458 }
459 }
460
461 double best_chi2 = 99999.9;
462 double param_2;
463 double best_seed_3;
464 double best_seed_4;
465 //TF1* fitFunc = new TF1("fitFunc", " (x<[2])*([0]*exp([1]*x)) + (x>=[2])*([3]*exp([4]*x))", xmin, xmax);
466 TF1* fitFunc = new TF1("fitFunc", " (x<[2])*([0]*exp([1]*x)) + (x>=[2])*(-[3]*x)", xmin, xmax);
467 fitFunc->FixParameter(0, seed_0);
468 fitFunc->FixParameter(1, seed_1);
469 for(int i=0; i < 1000; i++){
470 fitFunc->FixParameter(2, (double)i /10.0);
471 //TFitResultPtr intermed_fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc_seed, "QS", "", (double)i /10.0 ,xmax);
472 //double seed_3 = intermed_fitResult->Parameter(0);
473 //double seed_4 = intermed_fitResult->Parameter(1);
474 //fitFunc->SetParameter(3,seed_3);
475 //fitFunc->SetParameter(4,seed_4);
476
477 fitFunc_seed->FixParameter(0, seed_0);
478 fitFunc_seed->FixParameter(1, seed_1);
479 double seed_3 = fitFunc_seed->Eval(histos1d[histoname]->GetBinCenter(histos1d[histoname]->FindLastBinAbove(0.0)));
480 //std::cout << "LOOK FUCKER: " << seed_3 << std::endl;
481 //fitFunc->SetParameter(3,seed_3);
482
483 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax);
484
485 if(fitResult->Ndf() <= 0)
486 continue;
487
488 //std::cout << "test param 2: " << (double)i/10.0 << std::endl;
489 //std::cout << "chi2: " << fitResult->Chi2()/fitResult->Ndf() << std::endl;
490 if(fitResult->Chi2()/fitResult->Ndf() < best_chi2){
491 best_chi2 = fitResult->Chi2()/fitResult->Ndf();
492 param_2 = fitResult->Parameter(2);
493 best_seed_3 = seed_3;
494 //best_seed_4 = seed_4;
495 }
496 }
497
498 //delete ran;
499
500 fitFunc->FixParameter(0, seed_0);
501 fitFunc->FixParameter(1, seed_1);
502 fitFunc->FixParameter(2, param_2);
503 fitFunc->FixParameter(3, best_seed_3);
504 //fitFunc->FixParameter(4, best_seed_4);
505 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax);
506 fitFunc->Draw();
507
508 //return histogram to original range
509 histos1d[histoname]->GetXaxis()->SetRangeUser(originalXmin,originalXmax);
510
511 delete fitFunc_seed;
512 return fitFunc;
513}
514
515TF1* ZBiHistos::fitExponentialTail(std::string histogramName, double start_nevents){
516
517 //Get z vertex distribution corresponding to Test Cut
518 std::string histoname = m_name+"_"+histogramName+"_h";
519
520 //If histogram is empty, return nullptr
521 if(histos1d[histoname]->GetEntries() < 1){
522 std::cout << "[ZBiHistos]::WARNING: Background Model is NULL: " <<
523 histogramName << " distribution is empty and could not be fit!" << std::endl;
524 return nullptr;
525 }
526 double originalXmin = histos1d[histoname]->GetXaxis()->GetXmin();
527 double originalXmax = histos1d[histoname]->GetXaxis()->GetXmax();
528
529 //Start fit where start_nevents are left in tail
530 int lastbin = histos1d[histoname]->FindLastBinAbove(0.0);
531 if(histos1d[histoname]->Integral() < start_nevents)
532 start_nevents = histos1d[histoname]->GetBinContent(histos1d[histoname]->GetMaximumBin());
533 int firstbin = lastbin - 1;
534 double test_integral = 0.0;
535 while(test_integral < start_nevents && histos1d[histoname]->GetBinLowEdge(firstbin) > 0.0){
536 test_integral = histos1d[histoname]->Integral(firstbin, lastbin);
537 firstbin = firstbin - 1;
538 }
539
540 double xmin = histos1d[histoname]->GetBinLowEdge(firstbin);
541 double xmax = histos1d[histoname]->GetBinLowEdge(lastbin+1);
542 histos1d[histoname]->GetXaxis()->SetRangeUser(xmin,xmax);
543
544 //Initially fit tail with exponential to seed params
545 double best_seed_chi2 = 9999999.9;
546 double seed_0;
547 double seed_1;
548
549 TF1* fitFunc = new TF1("fitFunc", "[0]*exp([1]*x)", xmin, xmax);
550 TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax);
551 seed_0 = fitResult->Parameter(0);
552 seed_1 = fitResult->Parameter(1);
553 if(fitResult->Ndf() > 0){
554 best_seed_chi2 = fitResult->Chi2()/fitResult->Ndf();
555 }
556
557 TRandom3* ran = new TRandom3();
558 ran->SetSeed(0);
559
560 //best_seed_chi2 = 999999.9;
561 for(int i=0; i < 30; i++){
562 fitFunc->SetParameters(std::abs(ran->Gaus((double)start_nevents,3.0)), -std::abs(ran->Uniform(0.01,1.0)));
563 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax);
564
565 if(fitResult->Ndf() <= 0)
566 continue;
567 if(fitFunc->GetProb() < 0.001)
568 continue;
569
570 if(fitResult->Chi2()/fitResult->Ndf() < best_seed_chi2){
571 if(fitResult->Parameter(0) < 1000.0)
572 continue;
573 best_seed_chi2 = fitResult->Chi2()/fitResult->Ndf();
574 seed_0 = fitResult->Parameter(0);
575 seed_1 = fitResult->Parameter(1);
576 }
577 }
578
579 double param_0 = seed_0;
580 double param_1 = seed_1;
581 for(int i=0; i < 30; i++){
582 fitFunc->SetParameters(ran->Gaus(seed_0,1.0), -std::abs(ran->Gaus(seed_1,1.0)));
583 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax);
584 if(fitResult->Ndf() <= 0)
585 continue;
586 if(fitFunc->GetProb() < 0.001)
587 continue;
588
589 if(fitResult->Chi2()/fitResult->Ndf() < best_seed_chi2){
590 if(fitResult->Parameter(0) < 1000.0)
591 continue;
592 best_seed_chi2 = fitResult->Chi2()/fitResult->Ndf();
593 param_0 = fitResult->Parameter(0);
594 param_1 = fitResult->Parameter(1);
595 }
596 }
597
598 fitFunc->FixParameter(0, param_0);
599 fitFunc->FixParameter(1, param_1);
600 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax);
601 fitFunc->Draw();
602
603 //return histogram to original range
604 histos1d[histoname]->GetXaxis()->SetRangeUser(originalXmin,originalXmax);
605
606 delete ran;
607 return fitFunc;
608}
609
610
611
612TF1* ZBiHistos::fitExponentialPlusConst(std::string histogramName, double start_nevents){
613
614 //Get z vertex distribution corresponding to Test Cut
615 std::string histoname = m_name+"_"+histogramName+"_h";
616
617 //If histogram is empty, return nullptr
618 if(histos1d[histoname]->GetEntries() < 1){
619 std::cout << "[ZBiHistos]::WARNING: Background Model is NULL: " <<
620 histogramName << " distribution is empty and could not be fit!" << std::endl;
621 return nullptr;
622 }
623 double originalXmin = histos1d[histoname]->GetXaxis()->GetXmin();
624 double originalXmax = histos1d[histoname]->GetXaxis()->GetXmax();
625
626 //Start fit where start_nevents are left in tail
627 int lastbin = histos1d[histoname]->FindLastBinAbove(0.0);
628 if(histos1d[histoname]->Integral() < start_nevents)
629 start_nevents = histos1d[histoname]->GetBinContent(histos1d[histoname]->GetMaximumBin());
630 int firstbin = lastbin - 1;
631 double test_integral = 0.0;
632 while(test_integral < start_nevents && histos1d[histoname]->GetBinLowEdge(firstbin) > 0.0){
633 test_integral = histos1d[histoname]->Integral(firstbin, lastbin);
634 firstbin = firstbin - 1;
635 }
636
637 double xmin = histos1d[histoname]->GetBinLowEdge(firstbin);
638 double xmax = histos1d[histoname]->GetBinLowEdge(lastbin+1);
639 histos1d[histoname]->GetXaxis()->SetRangeUser(xmin,xmax);
640
641 //Initially fit tail with exponential to seed params
642 double best_seed_chi2 = 9999999.9;
643 double seed_0;
644 double seed_1;
645
646 TF1* fitFunc_seed = new TF1("fitFunc_seed", "[0]*exp([1]*x)", xmin, xmax);
647 //TF1* fitFunc_seed = new TF1("fitFunc_seed", "expo", xmin, xmax);
648 TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc_seed, "QSIM", "", xmin,xmax);
649 seed_0 = fitResult->Parameter(0);
650 seed_1 = fitResult->Parameter(1);
651 if(fitResult->Ndf() > 0 && fitFunc_seed->GetProb() > 0.001){
652 best_seed_chi2 = fitResult->Chi2()/fitResult->Ndf();
653 }
654
655 TRandom3* ran = new TRandom3();
656 ran->SetSeed(0);
657
658 //best_seed_chi2 = 999999.9;
659 for(int i=0; i < 30; i++){
660 fitFunc_seed->SetParameters(std::abs(ran->Gaus((double)start_nevents,3.0)), -std::abs(ran->Uniform(0.01,1.0)));
661 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc_seed, "QSIM", "", xmin,xmax);
662
663 if(fitResult->Ndf() <= 0)
664 continue;
665 if(fitFunc_seed->GetProb() < 0.001)
666 continue;
667
668 if(fitResult->Chi2()/fitResult->Ndf() < best_seed_chi2){
669 if(fitResult->Parameter(0) < 1000.0)
670 continue;
671 best_seed_chi2 = fitResult->Chi2()/fitResult->Ndf();
672 seed_0 = fitResult->Parameter(0);
673 seed_1 = fitResult->Parameter(1);
674 }
675 }
676
677 double best_chi2 = 99999.9;
678 double param_2;
679 TF1* fitFunc = new TF1("fitFunc", " (x<[2])*([0]*exp([1]*x)) + (x>=[2])*([0]*exp([1]*[2]))", xmin, xmax);
680 fitFunc->FixParameter(0, seed_0);
681 fitFunc->FixParameter(1, seed_1);
682 for(int i=0; i < 1000; i++){
683 fitFunc->FixParameter(2, (double)i /10.0);
684 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax);
685
686 if(fitResult->Ndf() <= 0)
687 continue;
688 if(fitFunc->GetProb() < 0.001)
689 continue;
690
691 if(fitResult->Chi2()/fitResult->Ndf() < best_chi2){
692 best_chi2 = fitResult->Chi2()/fitResult->Ndf();
693 param_2 = fitResult->Parameter(2);
694 }
695 }
696
697 delete ran;
698
699 fitFunc->FixParameter(0, seed_0);
700 fitFunc->FixParameter(1, seed_1);
701 fitFunc->FixParameter(2, param_2);
702 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QSIM", "", xmin,xmax);
703 fitFunc->Draw();
704
705 //return histogram to original range
706 histos1d[histoname]->GetXaxis()->SetRangeUser(originalXmin,originalXmax);
707
708 delete fitFunc_seed;
709 return fitFunc;
710}
711
712/*
713
714TF1* ZBiHistos::fitExponentialTail(std::string histogramName, double start_nevents){
715 //Background Model Fit Function
716 //TF1* fitFunc = new TF1("fitfunc","[0]*exp([1]*x)",0.0,100.0);
717 TF1* fitFunc = new TF1("fitfunc","expo",0.0,100.0);
718
719 //Get z vertex distribution corresponding to Test Cut
720 std::string histoname = m_name+"_"+histogramName+"_h";
721
722 //If histogram is empty, return nullptr
723 if(histos1d[histoname]->GetEntries() < 1){
724 std::cout << "[ZBiHistos]::WARNING: Background Model is NULL: " <<
725 histogramName << " distribution is empty and could not be fit!" << std::endl;
726 return nullptr;
727 }
728
729 //start mod
730 double originalXmin = histos1d[histoname]->GetXaxis()->GetXmin();
731 double originalXmax = histos1d[histoname]->GetXaxis()->GetXmax();
732
733 //Start fit where start_nevents are left in tail
734 int lastbin = histos1d[histoname]->FindLastBinAbove(0.0);
735 if(histos1d[histoname]->Integral() < start_nevents)
736 start_nevents = histos1d[histoname]->GetBinContent(histos1d[histoname]->GetMaximumBin());
737 int firstbin = lastbin - 1;
738 double test_integral = 0.0;
739 int peak_bin = histos1d[histoname]->GetMaximumBin();
740 while(test_integral < start_nevents && firstbin > peak_bin){
741 test_integral = histos1d[histoname]->Integral(firstbin, lastbin);
742 firstbin = firstbin - 1;
743 }
744
745 double xmin = histos1d[histoname]->GetBinLowEdge(firstbin);
746 double xmax = histos1d[histoname]->GetBinLowEdge(lastbin+1);
747 std::cout << "fit xmin " << xmin << std::endl;
748 std::cout << "fit xmax " << xmax << std::endl;
749 histos1d[histoname]->GetXaxis()->SetRangeUser(xmin,xmax);
750
751 fitFunc->SetParameters(50*start_nevents, -0.5);
752 TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
753 double parm0 = fitResult->Parameter(0);
754 double parm1 = fitResult->Parameter(1);
755
756 TRandom3* ran = new TRandom3();
757 ran->SetSeed(0);
758 double best_chi2 = 999999999.9;
759 double best_parm0 = parm0;
760 double best_parm1 = parm1;
761
762 int iteration = 50;
763 for(int i=0; i < iteration; i++){
764 fitFunc->SetParameters(std::abs(ran->Gaus(parm0,2)),ran->Uniform(-1.0,-0.1));
765 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
766 if(fitResult->Ndf() <= 0)
767 continue;
768
769 if(fitResult->Chi2()/fitResult->Ndf() < best_chi2){
770 best_chi2 = fitResult->Chi2()/fitResult->Ndf();
771 best_parm0 = fitResult->Parameter(0);
772 best_parm1 = fitResult->Parameter(1);
773 }
774 }
775
776 std::cout << "Fit Params: [0] = " << best_parm0 <<" [1] = " << best_parm1
777 << std::endl;
778 fitFunc->SetParameter(0,best_parm0);
779 fitFunc->SetParameter(1,best_parm1);
780 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "LSIM", "", xmin, xmax);
781 fitFunc->Draw();
782 //std::cout << "Fit Params: [0] = " << best_parm0 <<" [1] = " << best_parm1
783 // << " [2] = " << best_parm2 << std::endl;
784 std::cout << "Final fit param 0: " << fitResult->Parameter(0) << std::endl;
785 std::cout << "Final fit param 1: " << fitResult->Parameter(1) << std::endl;
786 //return histogram to original range
787 fitFunc->Draw();
788 histos1d[histoname]->GetXaxis()->SetRangeUser(originalXmin,originalXmax);
789
790 return fitFunc;
791}
792*/
793/*
794TF1* ZBiHistos::fitExponentialTail(std::string histogramName, double start_nevents){
795 //Background Model Fit Function
796 TF1* fitFunc = new TF1("fitfunc","[0]*exp([1]*x)",0.0,100.0);
797
798 //Get z vertex distribution corresponding to Test Cut
799 std::string histoname = m_name+"_"+histogramName+"_h";
800
801 //If histogram is empty, return nullptr
802 if(histos1d[histoname]->GetEntries() < 1){
803 std::cout << "[ZBiHistos]::WARNING: Background Model is NULL: " <<
804 histogramName << " distribution is empty and could not be fit!" << std::endl;
805 return nullptr;
806 }
807
808 //Start fit where start_nevents are left in tail
809 int lastbin = histos1d[histoname]->FindLastBinAbove(0.0);
810 int firstbin = lastbin - 1;
811 double test_integral = 0.0;
812 while(test_integral < start_nevents && histos1d[histoname]->GetBinLowEdge(firstbin) > 0.0){
813 test_integral = histos1d[histoname]->Integral(firstbin, lastbin);
814 firstbin = firstbin - 1;
815 }
816
817 double xmin = histos1d[histoname]->GetBinLowEdge(firstbin);
818 double xmax = histos1d[histoname]->GetBinLowEdge(lastbin+1);
819
820 //start mod
821 double originalXmin = histos1d[histoname]->GetXaxis()->GetXmin();
822 double originalXmax = histos1d[histoname]->GetXaxis()->GetXmax();
823
824 //Start fit where start_nevents are left in tail
825 int lastbin = histos1d[histoname]->FindLastBinAbove(0.0);
826 if(histos1d[histoname]->Integral() < start_nevents)
827 start_nevents = histos1d[histoname]->GetBinContent(histos1d[histoname]->GetMaximumBin());
828 int firstbin = lastbin - 1;
829 double test_integral = 0.0;
830 while(test_integral < start_nevents && histos1d[histoname]->GetBinLowEdge(firstbin) > 0.0){
831 test_integral = histos1d[histoname]->Integral(firstbin, lastbin);
832 firstbin = firstbin - 1;
833 }
834
835 double xmin = histos1d[histoname]->GetBinLowEdge(firstbin);
836 double xmax = histos1d[histoname]->GetBinLowEdge(lastbin+1);
837 std::cout << "fit xmin " << xmin << std::endl;
838 std::cout << "fit xmax " << xmax << std::endl;
839 //double xmax = 100.0;
840 histos1d[histoname]->GetXaxis()->SetRangeUser(xmin,xmax);
841
842 TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
843 double parm0 = fitResult->Parameter(0);
844 double parm1 = fitResult->Parameter(1);
845
846 TRandom3* ran = new TRandom3();
847 ran->SetSeed(0);
848 double best_chi2 = 999999999.9;
849 double best_parm0 = parm0;
850 double best_parm1 = parm1;
851 double best_parm2;
852
853 //fitFunc = new TF1("fitfunc","[0]*exp([1]*x) + abs([2])",0.0,100.0);
854 fitFunc = new TF1("fitfunc","[0]*exp([1]*x)",0.0,100.0);
855 int iteration = 25;
856 for(int i=0; i < iteration; i++){
857 //fitFunc->SetParameters(ran->Gaus(parm0,2),ran->Gaus(parm1,2), ran->Uniform(10));
858 fitFunc->SetParameters(ran->Gaus(parm0,2),ran->Gaus(parm1,2));
859 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
860 if(fitResult->Ndf() <= 0)
861 continue;
862
863 if(fitResult->Chi2()/fitResult->Ndf() < best_chi2){
864 best_chi2 = fitResult->Chi2()/fitResult->Ndf();
865 best_parm0 = fitResult->Parameter(0);
866 best_parm1 = fitResult->Parameter(1);
867 //best_parm2 = fitResult->Parameter(2);
868 }
869 }
870
871 //fitFunc->SetParameters(best_parm0, best_parm1, best_parm2);
872 fitFunc->SetParameters(best_parm0, best_parm1);
873 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "LSIM", "", xmin, xmax);
874 fitFunc->Draw();
875 //std::cout << "Fit Params: [0] = " << best_parm0 <<" [1] = " << best_parm1
876 // << " [2] = " << best_parm2 << std::endl;
877 std::cout << "Fit Params: [0] = " << best_parm0 <<" [1] = " << best_parm1
878 << std::endl;
879
880 //return histogram to original range
881 histos1d[histoname]->GetXaxis()->SetRangeUser(originalXmin,originalXmax);
882
883 return fitFunc;
884}
885*/
886
887/*
888TF1* ZBiHistos::fitExponentialTail(std::string histogramName, double start_nevents){
889 //Background Model Fit Function
890 TF1* fitFunc = new TF1("fitfunc","[0]*exp([1]*x)",10.0,100.0);
891 //TF1* fitFunc = new TF1("fitfunc","[0]*exp([1]*x) + [2]*exp([3]*x)",10.0,100.0);
892 //TF1* fitFunc = new TF1("fitfunc","[0]*exp([1]*x) + [2]",10.0,100.0);
893
894 //Get z vertex distribution corresponding to Test Cut
895 std::string histoname = m_name+"_"+histogramName+"_h";
896
897 //If histogram is empty, return nullptr
898 if(histos1d[histoname]->GetEntries() < 1){
899 std::cout << "[ZBiHistos]::WARNING: Background Model is NULL: " <<
900 histogramName << " distribution is empty and could not be fit!" << std::endl;
901 return nullptr;
902 }
903 //Start fit where start_nevents are left in tail
904 int lastbin = histos1d[histoname]->FindLastBinAbove(0.0);
905 int firstbin = lastbin - 1;
906 double test_integral = 0.0;
907 while(test_integral < start_nevents && histos1d[histoname]->GetBinLowEdge(firstbin) > 0.0){
908 test_integral = histos1d[histoname]->Integral(firstbin, lastbin);
909 firstbin = firstbin - 1;
910 }
911
912 double xmin = histos1d[histoname]->GetBinLowEdge(firstbin);
913 double xmax = histos1d[histoname]->GetBinLowEdge(lastbin+1);
914
915 TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
916 double parm0 = fitResult->Parameter(0);
917 double parm1 = fitResult->Parameter(1);
918 //double parm2 = fitResult->Parameter(2);
919 //double parm3 = fitResult->Parameter(3);
920
921 double best_chi2 = 9999999.9;
922 double best_parm0, best_parm1, best_parm2, best_parm3;
923
924 int iteration = 20;
925 for(int i=0; i < iteration; i++){
926 fitFunc->SetParameters(parm0,parm1);
927 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
928 if(fitResult->Ndf() <= 0)
929 continue;
930 if(fitResult->Chi2()/fitResult->Ndf() < best_chi2){
931 best_chi2 = fitResult->Chi2()/fitResult->Ndf();
932 best_parm0 = fitResult->Parameter(0);
933 best_parm1 = fitResult->Parameter(1);
934 //best_parm2 = fitResult->Parameter(2);
935 //best_parm3 = fitResult->Parameter(3);
936 }
937 }
938
939 //fitFunc = new TF1("fitfunc","[0]*exp([1]*x) + [2]*exp([3]*x)",10.0,100.0); //TRASH
940 fitFunc = new TF1("fitfunc","[0]*exp([1]*x) + [2]",10.0,100.0);
941 double parm2 = best_parm0*0.5;
942 for(int i=0; i < iteration; i++){
943 fitFunc->SetParameters(parm0,parm1,parm2);
944 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
945 if(fitResult->Ndf() <= 0)
946 continue;
947 if(fitResult->Chi2()/fitResult->Ndf() < best_chi2){
948 best_chi2 = fitResult->Chi2()/fitResult->Ndf();
949 best_parm0 = fitResult->Parameter(0);
950 best_parm1 = fitResult->Parameter(1);
951 best_parm2 = fitResult->Parameter(2);
952 //best_parm3 = fitResult->Parameter(3);
953 }
954 }
955
956
957 //fitFunc->SetParameters(best_parm0,best_parm1, best_parm2, best_parm3);
958 fitFunc->SetParameters(best_parm0,best_parm1, best_parm2);
959 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "LSIM", "", xmin, xmax);
960 fitFunc->Draw();
961
962 return fitFunc;
963}
964*/
965
966/*
967TF1* ZBiHistos::fitExponentialTail(std::string histogramName, double start_nevents){
968 //Background Model Fit Function
969 TF1* fitFunc = new TF1("fitfunc","[0]*exp([1]*x)",10.0,100.0);
970
971 //Get z vertex distribution corresponding to Test Cut
972 std::string histoname = m_name+"_"+histogramName+"_h";
973
974 //If histogram is empty, return nullptr
975 if(histos1d[histoname]->GetEntries() < 1){
976 std::cout << "[ZBiHistos]::WARNING: Background Model is NULL: " <<
977 histogramName << " distribution is empty and could not be fit!" << std::endl;
978 return nullptr;
979 }
980 //Start fit where start_nevents are left in tail
981 int lastbin = histos1d[histoname]->FindLastBinAbove(0.0);
982 int firstbin = lastbin - 1;
983 double test_integral = 0.0;
984 while(test_integral < start_nevents && histos1d[histoname]->GetBinLowEdge(firstbin) > 0.0){
985 test_integral = histos1d[histoname]->Integral(firstbin, lastbin);
986 firstbin = firstbin - 1;
987 }
988
989 double xmin = histos1d[histoname]->GetBinLowEdge(firstbin);
990 double xmax = histos1d[histoname]->GetBinLowEdge(lastbin+1);
991
992 TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
993 double parm0 = fitResult->Parameter(0);
994 double parm1 = fitResult->Parameter(1);
995
996 double best_chi2 = fitResult->Chi2()/fitResult->Ndf();
997 double best_parm0 = parm0;
998 double best_parm1 = parm1;
999
1000 TRandom3* ran = new TRandom3();
1001 ran->SetSeed(0);
1002 int iteration = 25;
1003 for(int i=0; i < iteration; i++){
1004 fitFunc->SetParameters(ran->Gaus(best_parm0,2),ran->Gaus(best_parm1,2));
1005 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
1006 if(fitResult->Ndf() <= 0)
1007 continue;
1008
1009 if(fitResult->Chi2()/fitResult->Ndf() < best_chi2){
1010 best_chi2 = fitResult->Chi2()/fitResult->Ndf();
1011 best_parm0 = fitResult->Parameter(0);
1012 best_parm1 = fitResult->Parameter(1);
1013 }
1014 }
1015
1016 fitFunc = new TF1("fitfunc","[0]*exp([1]*x) + [2]",10.0,100.0);
1017 double best_parm2;
1018 best_chi2 = 999999999.9;
1019 for(int i=0; i < iteration; i++){
1020 fitFunc->SetParameter(2,ran->Uniform(10));
1021 fitFunc->FixParameter(0,best_parm0);
1022 fitFunc->FixParameter(1,best_parm1);
1023 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
1024 if(fitResult->Ndf() <= 0)
1025 continue;
1026
1027 if(fitResult->Chi2()/fitResult->Ndf() < best_chi2){
1028 best_chi2 = fitResult->Chi2()/fitResult->Ndf();
1029 //best_parm0 = fitResult->Parameter(0);
1030 //best_parm1 = fitResult->Parameter(1);
1031 best_parm2 = fitResult->Parameter(2);
1032 //best_parm3 = fitResult->Parameter(3);
1033 }
1034 }
1035
1036 fitFunc->FixParameter(0,best_parm0);
1037 fitFunc->FixParameter(1,best_parm1);
1038 fitFunc->FixParameter(2,best_parm2);
1039 fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "LSIM", "", xmin, xmax);
1040 fitFunc->Draw();
1041
1042 return fitFunc;
1043}*/
1044
1045TF1* ZBiHistos::getZTailFit(std::string cutname){
1046 TF1* fitFunc = new TF1("fitfunc", "[0]* exp( (x<=([1]+[3]))*(-0.5*((x-[1])^2)/([2]^2)) + (x>=[1]+[3])*(-0.5*([3]^2/[2]^2)-(x-[1]-[3])/[4]) )", -100.0, 100.0);
1047
1048 std::string histoname = m_name+"_tritrig_zVtx_"+cutname+"_h";
1049 if(histos1d[histoname]->GetEntries() < 1){
1050 std::cout << "WARNING: Background Model is NULL: " <<
1051 cutname << " tritrig zVtx distribution was empty and could not be fit!" << std::endl;
1052 return nullptr;
1053 }
1054
1055 TFitResultPtr gausResult = (TFitResultPtr)histos1d[histoname]->Fit("gaus","QS");
1056 double gaus0 = gausResult->GetParams()[0];
1057 double gaus1 = gausResult->GetParams()[1];
1058 double gaus2 = gausResult->GetParams()[2];
1059 //gausResult = histos1d[histoname]->Fit("gaus","QS","",gaus1-2.5*gaus2, gaus1+10.0*gaus2);
1060 gausResult = histos1d[histoname]->Fit("gaus","QS","",gaus1-2.5*gaus2, gaus1+10.0*gaus2);
1061 gaus0 = gausResult->GetParams()[0];
1062 gaus1 = gausResult->GetParams()[1];
1063 gaus2 = gausResult->GetParams()[2];
1064 double xmin = gaus1 - 2.5*gaus2;
1065 double xmax = histos1d[histoname]->GetBinCenter(histos1d[histoname]->FindLastBinAbove(0.0)+1);
1066 xmax = 100.0;
1067
1068 //Fit function for a few different seeds for tail start, and length, keep the best fit
1069 double tailZ = gaus1 + 3.0*gaus2;
1070 double tail_l = 50.0;
1071 TRandom3* ran = new TRandom3();
1072 ran->SetSeed(0);
1073
1074 double best_chi2 = 9999999.9;
1075 double best_tailZ;
1076 double best_tail_l;
1077 double best_gaus0, best_gaus1, best_gaus2;
1078
1079 //Having issues with this fit sometimes. It occassionally underestimates bkg model
1080 int iteration = 80;
1081 tail_l = 50.0;
1082 for(int i =10; i < iteration; i=i+2){
1083 tailZ = gaus1 +(double)(iteration/10.0)*gaus2;
1084 fitFunc->SetParameters(gaus0, gaus1, gaus2, tailZ, tail_l);
1085 TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "QLSIM", "", xmin,xmax);
1086 if(fitResult->Ndf() <= 0)
1087 continue;
1088 if(fitResult->Chi2()/fitResult->Ndf() < best_chi2){
1089 best_chi2 = fitResult->Chi2()/fitResult->Ndf();
1090 best_gaus0 = fitResult->Parameter(0);
1091 best_gaus1 = fitResult->Parameter(1);
1092 best_gaus2 = fitResult->Parameter(2);
1093 best_tailZ = fitResult->Parameter(3);
1094 best_tail_l = fitResult->Parameter(4);
1095 tail_l = fitResult->Parameter(4);
1096 }
1097 }
1098
1099 fitFunc->SetParameters(best_gaus0, best_gaus1, best_gaus2, best_tailZ, best_tail_l);
1100 TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "LSIM", "", xmin, xmax);
1101 fitFunc->Draw();
1102 std::cout << "Finished fitting ztail" << std::endl;
1103
1104 delete ran;
1105 return fitFunc;
1106}
1107
1108double ZBiHistos::shosFitZTail(std::string cutname, double max_tail_events){
1109 TF1* fitFunc = new TF1("fitfunc", "[0]* exp( (x<=([1]+[3]))*(-0.5*((x-[1])^2)/([2]^2)) + (x>=[1]+[3])*(-0.5*([3]^2/[2]^2)-(x-[1]-[3])/[4]) )", -100.0, 100.0);
1110
1111 std::string histoname = m_name+"_tritrig_zVtx_"+cutname+"_h";
1112 if(histos1d[histoname]->GetEntries() < 1)
1113 return -4.3;
1114
1115 TFitResultPtr gausResult = (TFitResultPtr)histos1d[histoname]->Fit("gaus","QS");
1116 double gaus0 = gausResult->GetParams()[0];
1117 double gaus1 = gausResult->GetParams()[1];
1118 double gaus2 = gausResult->GetParams()[2];
1119 gausResult = histos1d[histoname]->Fit("gaus","QS","",gaus1-3.0*gaus2, gaus1+10.0*gaus2);
1120 gaus0 = gausResult->GetParams()[0];
1121 gaus1 = gausResult->GetParams()[1];
1122 gaus2 = gausResult->GetParams()[2];
1123 double xmin = gaus1 - 2.5*gaus2;
1124 double xmax = histos1d[histoname]->GetBinCenter(histos1d[histoname]->FindLastBinAbove(0.0)+1);
1125 xmax = 100.0;
1126
1127 //Fit function for a few different seeds for tail start, and length, keep the best fit
1128 double tailZ = gaus1 + 3.0*gaus2;
1129 double tail_l = 50.0;
1130 TRandom3* ran = new TRandom3();
1131 ran->SetSeed(0);
1132
1133 double best_chi2 = 99999.9;
1134 double best_tailZ;
1135 double best_tail_l;
1136 double best_gaus0, best_gaus1, best_gaus2;
1137
1138 int iteration = 80;
1139 for(int i =10; i < iteration; i=i+2){
1140 tail_l = 50.0;
1141 tailZ = gaus1 +(double)(iteration/10.0)*gaus2;
1142 fitFunc->SetParameters(gaus0, gaus1, gaus2, tailZ, tail_l);
1143 TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "LSIM", "", xmin,xmax);
1144 if(fitResult->Ndf() <= 0)
1145 continue;
1146 if(fitResult->Chi2() < best_chi2){
1147 best_chi2 = fitResult->Chi2()/fitResult->Ndf();
1148 best_gaus0 = fitResult->Parameter(0);
1149 best_gaus1 = fitResult->Parameter(1);
1150 best_gaus2 = fitResult->Parameter(2);
1151 best_tailZ = fitResult->Parameter(3);
1152 best_tail_l = fitResult->Parameter(4);
1153 }
1154 }
1155
1156 fitFunc->SetParameters(best_gaus0, best_gaus1, best_gaus2, best_tailZ, best_tail_l);
1157 TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "LSIM", "", xmin, xmax);
1158
1159 double zcut = -6.0;
1160 double testIntegral = fitFunc->Integral(zcut, 100.0);
1161 while(testIntegral > max_tail_events){
1162 zcut = zcut+0.1;
1163 testIntegral = fitFunc->Integral(zcut, 100.0);
1164 }
1165 fitFunc->Draw();
1166 delete ran;
1167 return zcut;
1168}
1169
1170double ZBiHistos::fitZTail(std::string zVtxHistoname, double max_tail_events){
1171 TF1* fitFunc = new TF1("fitfunc","[0]*exp( (((x-[1])/[2])<[3])*(-0.5*(x-[1])^2/[2]^2) + (((x-[1])/[2])>=[3])*(0.5*[3]^2-[3]*(x-[1])/[2]))", -100.0, 100.0);
1172
1173 std::string histoname = m_name+"_"+zVtxHistoname;
1174 if(histos1d[histoname]->GetEntries() < 1)
1175 return -4.3;
1176 TFitResultPtr gausResult = (TFitResultPtr)histos1d[histoname]->Fit("gaus","QS");
1177 double gaus1 = gausResult->GetParams()[1];
1178 double gaus2 = gausResult->GetParams()[2];
1179 gausResult = histos1d[histoname]->Fit("gaus","QS","",gaus1-3.0*gaus2, gaus1+10.0*gaus2);
1180 gaus1 = gausResult->GetParams()[1];
1181 gaus2 = gausResult->GetParams()[2];
1182 double tailZ = gaus1 + 3.0*gaus2;
1183
1184 fitFunc->SetParameters(gausResult->GetParams()[0], gaus1, gaus2, 3.0);
1185 TFitResultPtr fitResult = (TFitResultPtr)histos1d[histoname]->Fit(fitFunc, "LSIM", "", gaus1-2.0*gaus2, gaus1+10.0*gaus2);
1186
1187 double zcut = -6.0;
1188 double testIntegral = fitFunc->Integral(zcut, 90.0);
1189 while(testIntegral > max_tail_events){
1190 zcut = zcut+0.1;
1191 testIntegral = fitFunc->Integral(zcut, 90.0);
1192 }
1193 fitFunc->Draw();
1194 return zcut;
1195}
1196
1197void ZBiHistos::writeGraphs(TFile* outF, std::string folder){
1198 if (outF) outF->cd();
1199 TDirectory* dir{nullptr};
1200 std::cout<<folder.c_str()<<std::endl;
1201 if (!folder.empty()) {
1202 dir = outF->mkdir(folder.c_str(),"",true);
1203 dir->cd();
1204 }
1205 for(std::map<std::string,TGraph*>::iterator it=graphs_.begin(); it != graphs_.end(); ++it){
1206 if (!it->second){
1207 std::cout<<it->first<<" Null ptr in saving.."<<std::endl;
1208 continue;
1209 }
1210 it->second->Write();
1211 }
1212}
1213
1214void ZBiHistos::writeHistos(TFile* outF, std::string folder) {
1215 if (outF) outF->cd();
1216 TDirectory* dir{nullptr};
1217 std::cout<<folder.c_str()<<std::endl;
1218 if (!folder.empty()) {
1219 dir = outF->mkdir(folder.c_str(),"",true);
1220 dir->cd();
1221 }
1222 for (it1d it = histos1d.begin(); it!=histos1d.end(); ++it) {
1223 if (!it->second){
1224 std::cout<<it->first<<" Null ptr in saving.."<<std::endl;
1225 continue;
1226 }
1227 it->second->Write();
1228 }
1229 for (it2d it = histos2d.begin(); it!=histos2d.end(); ++it) {
1230 std::cout << it->first << std::endl;
1231 if (!it->second){
1232 std::cout<<it->first<<" Null ptr in saving.."<<std::endl;
1233 continue;
1234 }
1235 it->second->Write();
1236 }
1237}
1238
HistoManager::histos2d
std::map< std::string, TH2F * > histos2d
description
Definition HistoManager.h:335
HistoManager::plot2D
TH2F * plot2D(std::string name, std::string xtitle, int nbinsX, float xmin, float xmax, std::string ytitle, int nbinsY, float ymin, float ymax)
description
Definition HistoManager.cxx:203
HistoManager::it1d
std::map< std::string, TH1F * >::iterator it1d
description
Definition HistoManager.h:334
HistoManager::histos1d
std::map< std::string, TH1F * > histos1d
description
Definition HistoManager.h:333
HistoManager::it2d
std::map< std::string, TH2F * >::iterator it2d
description
Definition HistoManager.h:336
HistoManager::plot1D
TH1F * plot1D(const std::string &name, const std::string &xtitle, int nbinsX, float xmin, float xmax)
description
Definition HistoManager.cxx:189
HistoManager::m_name
std::string m_name
description
Definition HistoManager.h:331
HistoManager::debug_
bool debug_
description
Definition HistoManager.h:340
ZBiHistos::shosFitZTail
double shosFitZTail(std::string cutname, double max_tail_events)
description
Definition ZBiHistos.cxx:1108
ZBiHistos::getZTailFit
TF1 * getZTailFit(std::string cutname)
description
Definition ZBiHistos.cxx:1045
ZBiHistos::fitExponentialTail
TF1 * fitExponentialTail(std::string histogramName, double start_nevents)
description
Definition ZBiHistos.cxx:515
ZBiHistos::writeGraphs
void writeGraphs(TFile *outF, std::string folder)
description
Definition ZBiHistos.cxx:1197
ZBiHistos::cutFractionOfSignalVariable
double cutFractionOfSignalVariable(std::string cutvariable, bool isCutGreaterThan, double cutFraction, double initialIntegral)
description
Definition ZBiHistos.cxx:64
ZBiHistos::fitExponentialPlusExp
TF1 * fitExponentialPlusExp(std::string histogramName, double starnt_nevents)
description
Definition ZBiHistos.cxx:398
ZBiHistos::integrateHistogram1D
double integrateHistogram1D(std::string histoname)
description
Definition ZBiHistos.cxx:37
ZBiHistos::addHisto1d
void addHisto1d(std::string histoname, std::string xtitle, int nbinsX, float xmin, float xmax)
description
Definition ZBiHistos.cxx:15
ZBiHistos::fitZTail
double fitZTail(std::string zVtxHistoname, double max_tail_events)
description
Definition ZBiHistos.cxx:1170
ZBiHistos::set2DHistoYlabel
void set2DHistoYlabel(std::string histoName, int ybin, std::string ylabel)
description
Definition ZBiHistos.cxx:112
ZBiHistos::addHisto2d
void addHisto2d(std::string histoname, std::string xtitle, int nbinsX, float xmin, float xmax, std::string ytitle, int nbinsY, float ymin, float ymax)
description
Definition ZBiHistos.cxx:19
ZBiHistos::defineImpactParameterCut
std::vector< double > defineImpactParameterCut(double alpha=0.15)
description
Definition ZBiHistos.cxx:116
ZBiHistos::resetHistograms2d
void resetHistograms2d()
description
Definition ZBiHistos.cxx:30
ZBiHistos::defineZBiCutflowProcessorHistograms
void defineZBiCutflowProcessorHistograms()
description
Definition ZBiHistos.cxx:95
ZBiHistos::graphs_
std::map< std::string, TGraph * > graphs_
hold graphs
Definition ZBiHistos.h:156
ZBiHistos::resetHistograms1d
void resetHistograms1d()
description
Definition ZBiHistos.cxx:23
ZBiHistos::fitExponentialPlusConst
TF1 * fitExponentialPlusConst(std::string histogramName, double starnt_nevents)
description
Definition ZBiHistos.cxx:612
ZBiHistos::writeHistos
void writeHistos(TFile *outF, std::string folder)
description
Definition ZBiHistos.cxx:1214
Generated by doxygen 1.9.8