ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/yangyong/combineICEB/usefullcode.cc

Comparing UserCode/yangyong/combineICEB/usefullcode.cc (file contents):
Revision 1.1 by yangyong, Mon Dec 26 17:27:05 2011 UTC vs.
Revision 1.2 by yangyong, Tue Apr 10 19:41:12 2012 UTC

#	Line 0 | Line 1
1	+	TRandom3 *rgeng_;
2	+
3	+	float getMaximumY(TH1 *h1, TH1 *h2){
4	+	return h1->GetMaximum() > h2->GetMaximum() ? h1->GetMaximum() : h2->GetMaximum();
5	+	}
6	+
7	+	int getBinNumberwithLabel(TH1 *h1, string bnam){
8	+	int nbin1 = h1->GetNbinsX();
9	+
10	+
11	+	for(int b=1; b<=nbin1; b++){
12	+	if( h1->GetXaxis()->GetBinLabel(b) == bnam){
13	+	return b;
14	+	}
15	+	}
16	+	cout<<"not found " << bnam.c_str()<<endl;
17	+	return -1;
18	+	}
19	+
20	+
21	+	float getGaussian(float mean, float sigma, float min,float max){
22	+
23	+	if( rgeng_ == 0){
24	+	rgeng_ = new TRandom3(0);
25	+	}
26	+
27	+	float tmp = rgeng_->Gaus(mean,sigma);
28	+	while( tmp <= min || tmp >= max ){
29	+	tmp = rgeng_->Gaus(mean,sigma);
30	+	}
31	+	return tmp;
32	+
33	+	}
34	+
35	+
36	+
37	+
38	+	double interpolatedTwoBinTH1(TH1 *h1, double c){
39	+
40	+	int nbin1 = h1->GetNbinsX();
41	+	for(int b=1; b<nbin1; b++){
42	+	double c1 = h1->GetBinCenter(b);
43	+	double c2 = h1->GetBinCenter(b+1);
44	+	if( c >= c1 && c <= c2){
45	+	double v1 = h1->GetBinContent(b);
46	+	double v2 = h1->GetBinContent(b+1);
47	+	double k = (v2-v1)/(c2-c1);
48	+	return v1 + k * ( c-c1);
49	+	}
50	+	}
51	+	return 0;
52	+
53	+	}
54	+
55	+
56	+	void NormHistogram(TH1 *h1){
57	+	float total = h1->GetEntries();
58	+
59	+	if( total <=0){
60	+	cout<<" empty histogram? " << total <<endl;
61	+	}
62	+
63	+	h1->Scale(1./total);
64	+
65	+	}
66	+
67	+
68	+	void emptyHistogram( TH1 *h1){
69	+	int nbin1 = h1->GetNbinsX();
70	+	for(int b=1; b<=nbin1; b++){
71	+	h1->SetBinContent(b,0);
72	+	}
73	+
74	+	}
75	+
76	+
77	+
78	+	void AddHistByQuadrature(TH1 *h1, float xx){
79	+
80	+	int nbin1 = h1->GetNbinsX();
81	+
82	+
83	+	for(int b=1; b<=nbin1; b++){
84	+
85	+	float tmp1 =  h1->GetBinContent(b);
86	+	h1->SetBinContent(b,sqrt(tmp1*tmp1 + xx*xx));
87	+
88	+	}
89	+
90	+	}
91	+
92	+
93	+
94	+
95	+
96	+	void subtractHistByQuadrature(TH1 *h1, float xx){
97	+
98	+	int nbin1 = h1->GetNbinsX();
99	+
100	+
101	+	for(int b=1; b<=nbin1; b++){
102	+
103	+	float tmp1 =  h1->GetBinContent(b);
104	+	if(tmp1 > xx){
105	+	h1->SetBinContent(b,sqrt(tmp1*tmp1 - xx*xx));
106	+	}else{
107	+	///h1->SetBinContent(b,sqrt(tmp1*tmp1 - xx*xx));
108	+	}
109	+	}
110	+
111	+	}
112	+
113	+
114	+	void scaleHistogramMCtoData(TH1 *hda, TH1 *hmc){
115	+
116	+	double ndata = hda->Integral();
117	+	double nmc = hmc->Integral();
118	+
119	+	if( ndata >0 && nmc >0){
120	+	hmc ->Scale( ndata/nmc);
121	+	}
122	+
123	+
124	+	}
125	+
126	+
127	+	void scaleHistogramMCtoDatav1(TH1 *hda, TH1 *hmc, float xlow, float xhigh){
128	+
129	+	double width = hda->GetBinWidth(1);
130	+	int nbin_start =  int ( (xlow - hda->GetXaxis()->GetXmin())/ width) +1;
131	+	int nbin_end =  int ( (xhigh - hda->GetXaxis()->GetXmin())/ width);
132	+
133	+	double ndata = hda->Integral(nbin_start, nbin_end);
134	+	double nmc = hmc->Integral(nbin_start,nbin_end);
135	+
136	+	cout<<" nbin_start/end "<< nbin_start <<" "<< nbin_end <<" "<< ndata <<" "<<nmc <<endl;
137	+
138	+	if( ndata >0 && nmc >0){
139	+	hmc ->Scale( ndata/nmc);
140	+	}
141	+
142	+
143	+	}
144	+
145	+
146	+
147	+	void addTH2(TH2 *h1,TH2 *h2){
148	+
149	+	int binx = h1->GetNbinsX();
150	+	int biny = h1->GetNbinsY();
151	+
152	+	if( binx !=  h2->GetNbinsX() || biny != h2->GetNbinsY()){
153	+	cout<<"diff bins addTH2 ? " <<binx <<" "<< h2->GetNbinsX() <<" "<<  biny <<" "<< h2->GetNbinsY()<<endl;
154	+	return;
155	+	}
156	+	for(int bx = 1; bx<= binx; bx++){
157	+	for(int by = 1; by<= biny; by++){
158	+	double tmp1 = h1->GetBinContent(bx,by);
159	+	double tmp2 = h2->GetBinContent(bx,by);
160	+	h2->SetBinContent(bx,by,tmp1+tmp2);
161	+	}
162	+	}
163	+
164	+	}
165	+
166	+
167	+	void subtractTwoHistQuadrature(TH1 *h1, TH1 *h2, TH1 *h3){
168	+
169	+	int nbin1 = h1->GetNbinsX();
170	+	int nbin2 = h2->GetNbinsX();
171	+
172	+	if( nbin1 != nbin2) {
173	+	cout<<" subtractTwoHistQuadrature nbin " << nbin1 <<" "<<nbin2 <<endl;
174	+	}
175	+
176	+	for(int b=1; b<=nbin1; b++){
177	+
178	+	float tmp1 =  h1->GetBinContent(b);
179	+	float tmp1err =  h1->GetBinError(b);
180	+	float tmp2 =  h2->GetBinContent(b);
181	+	float tmp2err =  h2->GetBinError(b);
182	+
183	+	float diff = -1;
184	+	float diffErr = 0;
185	+	if( tmp1 >= tmp2){
186	+	diff = sqrt( tmp1*tmp1 - tmp2*tmp2);
187	+	diffErr = sqrt(  pow( tmp1/diff*tmp1err,2) +  pow( tmp2/diff*tmp2err,2) );
188	+
189	+	//cout<<"check " << tmp1 <<" "<< tmp2 <<" "<< diff <<endl;
190	+
191	+	}
192	+	h3->SetBinContent(b,diff);
193	+	h3->SetBinError(b,diffErr);
194	+	}
195	+
196	+
197	+	}
198	+
199	+
200	+
201	+
202	+	double one_sigma_up(std::vector<double> limits) {
203	+	double v_one_sigma_up = 0;
204	+	sort(limits.begin(),limits.end());
205	+	int samplet_size = limits.size();
206	+	int v_one_sigma_up_index = 0;
207	+	v_one_sigma_up_index = samplet_size * 841 / 1000;
208	+	v_one_sigma_up = limits[v_one_sigma_up_index-1];
209	+	return v_one_sigma_up;
210	+	}
211	+
212	+	double one_sigma_down(std::vector<double> limits) {
213	+	double v_one_sigma_down = 0;
214	+	sort(limits.begin(),limits.end());
215	+	int samplet_size = limits.size();
216	+	int v_one_sigma_down_index = 0;
217	+	v_one_sigma_down_index = samplet_size * 159 / 1000;
218	+	v_one_sigma_down = limits[v_one_sigma_down_index];
219	+	return v_one_sigma_down;
220	+	}
221	+
222	+	double two_sigma_up(std::vector<double> limits) {
223	+	double v_two_sigma_up = 0;
224	+	sort(limits.begin(),limits.end());
225	+	int samplet_size = limits.size();
226	+	int v_two_sigma_up_index = 0;
227	+	v_two_sigma_up_index = samplet_size * 979 / 1000;
228	+	v_two_sigma_up = limits[v_two_sigma_up_index-1];
229	+	return v_two_sigma_up;
230	+	}
231	+
232	+	double two_sigma_down(std::vector<double> limits) {
233	+	double v_two_sigma_down = 0;
234	+	sort(limits.begin(),limits.end());
235	+	int samplet_size = limits.size();
236	+	int v_two_sigma_down_index = 0;
237	+	v_two_sigma_down_index = samplet_size * 21 / 1000;
238	+	v_two_sigma_down = limits[v_two_sigma_down_index];
239	+	return v_two_sigma_down;
240	+	}
241	+
242	+
243	+
244	+	TH1F *divideHistogramDataMC(TH1 *hdata, TH1 *hmc){
245	+
246	+	TString hname = TString(hdata->GetName()) + TString(hmc->GetName());
247	+	TH1F *htmp = (TH1F*)hdata->Clone(hname);
248	+
249	+	for(int b=1; b<= hdata->GetNbinsX();b++){
250	+	float r=0;
251	+	float rErr = 0;
252	+	if( hmc->GetBinContent(b) >0){
253	+	r = hdata->GetBinContent(b) / hmc->GetBinContent(b);
254	+	rErr = sqrt(r) / hmc->GetBinContent(b);
255	+	}
256	+	htmp->SetBinContent(b,r);
257	+	htmp->SetBinError(b,rErr);
258	+
259	+	}
260	+	return htmp;
261	+
262	+
263	+	}
264	+
265	+
266	+	float getIntegralErrorHistogram(TH1F *h1){
267	+
268	+	float err = 0;
269	+	for(int b=1; b<= h1->GetNbinsX(); b++){
270	+	err += pow( h1->GetBinError(b),2);
271	+
272	+	}
273	+	return sqrt(err);
274	+	}
275	+
276	+
277	+	void getIntegralAndWtErrorHistogram(TH1F *h1, int startBin, float res[]){
278	+
279	+	float tot = 0;
280	+	float err = 0;
281	+	for(int b=startBin; b<= h1->GetNbinsX(); b++){
282	+	err += pow( h1->GetBinError(b),2);
283	+	tot +=  h1->GetBinContent(b);
284	+	}
285	+
286	+	res[0] = tot;
287	+	res[1] = sqrt(err);
288	+
289	+
290	+	}
291	+
292	+
293	+
294	+	void printbinContentHistogram(TH1 *h1){
295	+
296	+	for(int b=1; b<= h1->GetNbinsX(); b++){
297	+	cout<<"bin: "<< h1->GetBinLowEdge(b)<<"to"<<  h1->GetBinLowEdge(b)+h1->GetBinWidth(b) <<" "<< h1->GetBinContent(b)<<"+/-"<< h1->GetBinError(b)<<endl;
298	+	}
299	+
300	+	}
301	+
302	+
303	+
304	+
305	+	void calcRatioBayes(float n1, float n2, float &r, float &rup, float &rlow){
306	+
307	+	//compute using Bayes
308	+	TH1D h1("dummy1","",1,1,2);
309	+	h1.SetBinContent(1,n1);
310	+
311	+	TH1D h2("dummy2","",1,1,2);
312	+	h2.SetBinContent(1,n2);
313	+
314	+	TGraphAsymmErrors g;
315	+	g.BayesDivide(&h1,&h2);
316	+	r = g.GetY()[0];
317	+	rup = g.GetErrorYhigh(0);
318	+	rlow = g.GetErrorYlow(0);
319	+
320	+	}
321	+
322	+
323	+
324	+
325	+	TGraphAsymmErrors* dividingTwoHistogramNoErrBinContent(TH1 *h1, TH1 *h2){
326	+
327	+
328	+
329	+	float x[1000];
330	+	float exl[1000];
331	+	float exh[1000];
332	+
333	+	float y[1000];
334	+	float yup[1000];
335	+	float ydw[1000];
336	+
337	+
338	+	if( h1->GetNbinsX()  != h2->GetNbinsX() ){
339	+	cout<<" dividingTwoHistogramNoErrBinContent diff bins " << endl;
340	+	return 0 ;
341	+	}
342	+
343	+	if( h1->GetNbinsX() > 1000){
344	+	cout<<" dividingTwoHistogramNoErrBinContent too much bins " << endl;
345	+	return 0 ;
346	+	}
347	+
348	+	int n =0;
349	+
350	+	for(int b=1; b<= h1->GetNbinsX(); b++){
351	+	float b1 = h1->GetBinContent(b);
352	+	float b2 = h2->GetBinContent(b);
353	+
354	+	float r,rup,rdown;
355	+
356	+
357	+	x[n] = h1->GetBinCenter(b);
358	+	exl[n] = 0.5*h1->GetBinWidth(b);
359	+	exh[n] = 0.5*h1->GetBinWidth(b);
360	+
361	+	if( b1>0 && b2>0){
362	+	if( b1<=b2){
363	+	calcRatioBayes(b1,b2,r,rup,rdown);
364	+	}else{
365	+	calcRatioBayes(b2,b1,r,rup,rdown);
366	+	}
367	+	}else{
368	+	r = 0;
369	+	rup = 0;
370	+	rdown = 0;
371	+	}
372	+
373	+
374	+	y[n] = b1/b2;
375	+	yup[n] = rup;
376	+	ydw[n] = rdown;
377	+
378	+
379	+	n ++;
380	+
381	+	}
382	+	TGraphAsymmErrors *g1 = new TGraphAsymmErrors(n, x, y, exl, exh, ydw,yup);
383	+	return g1;
384	+
385	+
386	+	}
387	+
388	+
389	+	double ErrorInProduct(double x, double errx, double y,
390	+	double erry, double corr) {
391	+	double xFrErr = errx/x;
392	+	double yFrErr = erry/y;
393	+	return sqrt(pow(xFrErr,2) + pow(yFrErr,2) + 2.0*corr*xFrErr*yFrErr)*x*y;
394	+	}
395	+
396	+
397	+
398	+	void  getErrorInRatio(double x, double errx, double y,
399	+	double erry, double res[]) {
400	+
401	+	if( x==0 || y==0){
402	+	cout<<"getErrorInRatio x/y" << x <<" "<<y<<endl;
403	+	exit(1);
404	+	}
405	+
406	+	res[0] = x/y;
407	+	res[1] = x/y *sqrt( pow(errx/x,2) + pow(erry/y,2));
408	+
409	+	}
410	+
411	+
412	+
413	+	void get3Max(double x1, double x2, double x3,double res[]){
414	+	double max = x1;
415	+	res[1] = 0;
416	+
417	+	if(max<x2) {
418	+	max = x2;
419	+	res[1] = 1;
420	+	}
421	+	if(max<x3){
422	+	max = x3;
423	+	res[1] = 2;
424	+	}
425	+
426	+	res[0] = max;
427	+	//  return max;
428	+
429	+	}
430	+
431	+
432	+
433	+	float getMaxOfThree(float x1, float x2, float x3){
434	+
435	+	float xmax = x1 > x2 ? x1: x2;
436	+	xmax = xmax > x3 ? xmax: x3;
437	+	return xmax;
438	+
439	+
440	+	}
441	+
442	+
443	+	void generateToyDataFromHist(TH1* h1, TH1 *h2){
444	+
445	+	int nbinx1 = h1->GetNbinsX();
446	+	int nbinx2 = h2->GetNbinsX();
447	+
448	+	if(nbinx1 != nbinx2){
449	+	cout<<" generateToyDataFromHist " << nbinx1 <<" "<<  nbinx2 <<endl;
450	+	exit(1);
451	+	}
452	+	TRandom3 *grand = new TRandom3(0);
453	+	for(int b=1 ; b<= nbinx1; b++){
454	+	float tmp = h1->GetBinContent(b);
455	+	float tmp1 = grand->Poisson(tmp);
456	+	h2->SetBinContent(b,tmp1);
457	+	}
458	+
459	+
460	+	}
461	+
462	+
463	+	void calcTotalErorrPDF(int npdf, double X0,double X0Err,double XP[100],double XPErr[100],double XM[100],double XMErr[100],double results[]){
464	+
465	+
466	+	double dXP = 0;
467	+	double dXM = 0;
468	+	double simpleDX = 0;
469	+
470	+	double dsimpleDX = 0;
471	+	double ddXP = 0;
472	+	double ddXM = 0;
473	+
474	+	double res[10];
475	+
476	+	if(npdf%2 !=0) {
477	+	cout<<" calcTotalErorrPDF npdf: "<<npdf <<endl;
478	+	exit(1);
479	+	}
480	+
481	+
482	+	for( int j=1; j<= npdf/2; j++){
483	+
484	+	get3Max(XP[j]-X0,XM[j]-X0,0,res);
485	+	double dxP = res[0];
486	+	if( fabs(res[1]-0)<0.01){
487	+	ddXP += 4*pow(XP[j]-X0,2)*(XPErr[j]*XPErr[j] + X0Err*X0Err);
488	+	}else if( fabs(res[1]-1)<0.01){
489	+	ddXP += 4*pow(XM[j]-X0,2)*(XMErr[j]*XMErr[j] + X0Err*X0Err);
490	+	}
491	+
492	+	get3Max(X0-XP[j],X0-XM[j],0,res);
493	+	double dxM = res[0];
494	+	if( fabs(res[1]-0)<0.01){
495	+	ddXM += 4*pow(XP[j]-X0,2)*(XPErr[j]*XPErr[j] + X0Err*X0Err);
496	+	}else if( fabs(res[1]-1)<0.01){
497	+	ddXM += 4*pow(XM[j]-X0,2)*(XMErr[j]*XMErr[j] + X0Err*X0Err);
498	+	}
499	+
500	+	dXP += dxP*dxP;
501	+	dXM += dxM*dxM;
502	+
503	+	simpleDX +=  pow(XP[j]-XM[j],2);
504	+
505	+	dsimpleDX += pow(XP[j]-XM[j],2)*(XPErr[j]*XPErr[j]+XMErr[j]*XMErr[j]);
506	+
507	+
508	+	///cout<<"checkme: "<<j<<"  "<< dXP <<" "<< dXM <<endl;
509	+
510	+
511	+	}
512	+
513	+
514	+	dXP = sqrt(dXP);
515	+	dXM = sqrt(dXM);
516	+
517	+	ddXP = sqrt(ddXP);
518	+
519	+	if( dXP>0){
520	+	ddXP = ddXP/(2*dXP);
521	+	}
522	+	ddXM = sqrt(ddXM);
523	+	if( dXM >0){
524	+	ddXM = ddXM/(2*dXM);
525	+	}
526	+
527	+
528	+
529	+
530	+	simpleDX = 0.5* sqrt(simpleDX);
531	+	dsimpleDX =  sqrt(dsimpleDX);
532	+	if(simpleDX >0){
533	+	dsimpleDX = dsimpleDX/(2*simpleDX);
534	+	}
535	+
536	+	////asymmetry errors
537	+	results[0] = dXP;
538	+	results[1] = dXM;
539	+	///master formula
540	+	results[2] = simpleDX;
541	+
542	+	////error of above
543	+	results[3] = ddXP;
544	+	results[4] = ddXM;
545	+	results[5] = dsimpleDX;
546	+
547	+
548	+	}
549	+
550	+	int mostProb(float b){
551	+
552	+	int nbest = int(b+1);
553	+	float prob_max = 0;
554	+	for(int n = int(b-1) ; n<= int(b+1); n++){
555	+	if( n<0) n =0;
556	+
557	+
558	+	float prob = TMath::Poisson(n,b);
559	+
560	+	//    cout<<prob<<" "<<n<<" "<<b<<endl;
561	+
562	+	if( prob > prob_max){
563	+	prob_max = prob;
564	+	nbest = n;
565	+	}
566	+	}
567	+	return nbest;
568	+
569	+	}
570	+
571	+
572	+	float getBinContentHistogram(TH1F *hhtmp, float pt){
573	+
574	+	for(int b=1; b<= hhtmp->GetNbinsX(); b++){
575	+	if( pt >= hhtmp->GetBinLowEdge(b) && pt < hhtmp->GetBinLowEdge(b) +  hhtmp->GetBinWidth(b)){
576	+	return  hhtmp->GetBinContent(b);
577	+	}
578	+	}
579	+	return -999;
580	+
581	+
582	+	}
583	+
584	+
585	+
586	+
587	+
588	+
589	+	void weightedAeverageTwoHistogram(TH1F *hhtmp1, TH1F *hhtmp2, TH1F *hhres){
590	+
591	+	if( hhtmp1->GetNbinsX() != hhtmp2->GetNbinsX()){
592	+	cout<<"warning  weightedAeverageTwoHistogram  "<< hhtmp1->GetNbinsX() <<" "<< hhtmp2->GetNbinsX() <<endl;
593	+	exit(1);
594	+	}
595	+
596	+	for(int b=1; b<= hhtmp1->GetNbinsX(); b++){
597	+	float tmp1 =  hhtmp1->GetBinContent(b);
598	+	float tmp2 =  hhtmp2->GetBinContent(b);
599	+	float tmp1e =  hhtmp1->GetBinError(b);
600	+	float tmp2e =  hhtmp2->GetBinError(b);
601	+
602	+	float av = 0;
603	+	float ave = 0;
604	+	if( tmp1 > 0 && tmp2 >0){
605	+	av = (tmp1/(tmp1e*tmp1e) + tmp2/ (tmp2e*tmp2e)) / ( 1./(tmp1e*tmp1e) + 1./(tmp2e*tmp2e));
606	+	ave =  sqrt(1/ (1./(tmp1e*tmp1e) + 1./(tmp2e*tmp2e)) );
607	+	}else if (tmp1 > 0 && tmp2 ==0){
608	+	av = tmp1;
609	+	ave = tmp1e;
610	+	}else if( tmp1 ==0 && tmp2 >0){
611	+	av = tmp2;
612	+	ave = tmp2e;
613	+	}else{
614	+	cout<<"no data for combine.."<< b<<endl;
615	+	}
616	+
617	+	hhres->SetBinContent(b,av);
618	+	hhres->SetBinError(b,ave);
619	+
620	+	}
621	+
622	+	}
623	+
624	+
625	+
626	+
627	+	float getFractionOfHistogram(TH1F *hhtmp, double xmin, double xmax){
628	+
629	+	double sum = hhtmp->Integral();
630	+	if( sum <=0) {
631	+	cout<<" getFractionOfHistogram zero.." << sum << " "<< hhtmp->GetEntries();
632	+	exit(1);
633	+	}
634	+
635	+	int nbinsx = hhtmp->GetNbinsX();
636	+	int bin1 =1;
637	+	int bin2 = nbinsx;
638	+
639	+	for(int b=1; b<= nbinsx; b++){
640	+	if( xmin >= hhtmp->GetBinLowEdge(b) && xmin -hhtmp->GetBinLowEdge(b) - hhtmp->GetBinWidth(b) < 1E-6 ){
641	+	bin1 = b;
642	+	}
643	+	if( xmax > hhtmp->GetBinLowEdge(b) && xmax - hhtmp->GetBinLowEdge(b) - hhtmp->GetBinWidth(b) < 1E-6 ){
644	+	bin2 = b;
645	+	}
646	+	}
647	+
648	+	///cout<<"checkme"<<bin1 <<" "<<bin2 <<" "<< hhtmp->GetBinLowEdge(bin1)  <<" "<< xmax - hhtmp->GetBinLowEdge(bin2) - hhtmp->GetBinWidth(bin2) <<" "<< xmax - hhtmp->GetBinLowEdge(bin2-1) - hhtmp->GetBinWidth(bin2-1)<<endl;
649	+
650	+	return hhtmp->Integral(bin1,bin2) / sum;
651	+
652	+
653	+	}
654	+
655	+
656	+
657	+
658	+	void cloneHistogramWithWeight(TH1 *hhtmp1, TH1* hhtmp2, float wt){
659	+
660	+	int nbins = hhtmp1->GetNbinsX();
661	+	if( nbins != hhtmp2->GetNbinsX()){
662	+	cout<<"error.cloneHistogramWithWeight  "<<nbins <<" "<< hhtmp2->GetNbinsX() << endl;
663	+	exit(1);
664	+	}
665	+
666	+	for(int j=1; j<= nbins+1; j++){ ///bins+1 for overflowbins
667	+
668	+	float tmp = hhtmp1->GetBinContent(j) * wt;
669	+
670	+	hhtmp2->SetBinContent(j,tmp);
671	+
672	+	}
673	+
674	+
675	+	}
676	+
677	+
678	+
679	+
680	+	void cloneHistogramWithWeightv1(TH1 *hhtmp1, TH1* hhtmp2, float wt){
681	+
682	+	int nbins = hhtmp1->GetNbinsX();
683	+	if( nbins >  hhtmp2->GetNbinsX()){
684	+	cout<<"error.cloneHistogramWithWeightv1  "<<nbins <<" "<< hhtmp2->GetNbinsX() << endl;
685	+	exit(1);
686	+	}
687	+
688	+	for(int j=1; j<= nbins; j++){
689	+	double tmp = hhtmp1->GetBinContent(j) * wt;
690	+	hhtmp2->SetBinContent(j,tmp);
691	+	}
692	+
693	+
694	+	}
695	+
696	+
697	+
698	+
699	+	///check which bin the value in histogram
700	+	int getBinOfHistogram(TH1 *hh, float tmp){
701	+
702	+	int nbins = hh->GetNbinsX();
703	+
704	+	if( tmp < hh->GetXaxis()->GetXmin() || tmp >= hh->GetXaxis()->GetXmax()) return -1;
705	+
706	+	for(int b = 1; b<= nbins; b++){
707	+
708	+	if( tmp >= hh->GetBinLowEdge(b) && tmp < hh->GetBinLowEdge(b) +  hh->GetBinWidth(b)) return b;
709	+
710	+	}
711	+
712	+	cout<<"error.. getBinOfHistogram "<<tmp <<endl;
713	+
714	+	return -1;
715	+
716	+
717	+	}
718	+
719	+
720	+
721	+
722	+
723	+	///check which bin the value in histogram
724	+	int getBinOfHistogramv1(TH1 *hh, float tmp){
725	+
726	+	int nbins = hh->GetNbinsX();
727	+
728	+	if( tmp < hh->GetXaxis()->GetXmin() ) return 1;
729	+	if( tmp >= hh->GetXaxis()->GetXmax() ) return nbins;
730	+
731	+	for(int b = 1; b<= nbins; b++){
732	+
733	+	if( tmp >= hh->GetBinLowEdge(b) && tmp < hh->GetBinLowEdge(b) +  hh->GetBinWidth(b)) return b;
734	+
735	+	}
736	+
737	+	cout<<"error.. getBinOfHistogram "<<tmp <<endl;
738	+	exit(1);
739	+
740	+	return -1;
741	+
742	+
743	+	}
744	+
745	+
746	+
747	+
748	+
749	+
750	+	double BinP(int N, double p, int x1, int x2) {
751	+	double q=p/(1-p);
752	+	int k=0;
753	+	double v = 1;
754	+	double s=0;
755	+	double tot=0.0;
756	+	while(k<=N) {
757	+	tot=tot+v;
758	+	if(k>=x1 & k<=x2) { s=s+v; }
759	+	if(tot>1e30){s=s/1e30; tot=tot/1e30; v=v/1e30;}
760	+	k=k+1;
761	+	v=v*q*(N+1-k)/k;
762	+	}
763	+	return s/tot;
764	+	}
765	+
766	+
767	+
768	+
769	+	void ClopperPearsonLimits(double numerator, double denominator, double &ratio,
770	+	double &lowerLimit, double &upperLimit, const double CL_low=1.0,
771	+	const double CL_high=1.0)
772	+	{
773	+	//Confidence intervals are in the units of \sigma.
774	+
775	+
776	+	ratio = numerator/denominator;
777	+
778	+	// first get the lower limit
779	+	if(numerator==0)   lowerLimit = 0.0;
780	+	else {
781	+	double v=ratio/2;
782	+	double vsL=0;
783	+	double vsH=ratio;
784	+	double p=CL_low/100;
785	+	while((vsH-vsL)>1e-5) {
786	+	if(BinP(int(denominator),v,int(numerator),int(denominator))>p)
787	+	{ vsH=v; v=(vsL+v)/2; }
788	+	else { vsL=v; v=(v+vsH)/2; }
789	+	}
790	+	lowerLimit = v;
791	+	}
792	+
793	+	// now get the upper limit
794	+	if(numerator==denominator) upperLimit = 1.0;
795	+	else {
796	+	double v=(1+ratio)/2;
797	+	double vsL=ratio;
798	+	double vsH=1;
799	+	double p=CL_high/100;
800	+	while((vsH-vsL)>1e-5) {
801	+	if(BinP(int(denominator),v,0,int(numerator))<p) { vsH=v; v=(vsL+v)/2; }
802	+	else { vsL=v; v=(v+vsH)/2; }
803	+	}
804	+	upperLimit = v;
805	+	}
806	+	}
807	+
808	+
809	+
810	+
811	+
812	+

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines