yangyong/zSelection/usefullcode.cc


float getMaximumY(TH1 *h1, TH1 *h2){
  return h1->GetMaximum() > h2->GetMaximum() ? h1->GetMaximum() : h2->GetMaximum();
}

int getBinNumberwithLabel(TH1 *h1, string bnam){
  int nbin1 = h1->GetNbinsX();
  
   
  for(int b=1; b<=nbin1; b++){
    if( h1->GetXaxis()->GetBinLabel(b) == bnam){
      return b; 
    }
  }
  cout<<"not found " << bnam.c_str()<<endl; 
  return -1; 
}


void AddHistByQuadrature(TH1 *h1, float xx){
  
  int nbin1 = h1->GetNbinsX();

  
  for(int b=1; b<=nbin1; b++){
    
    float tmp1 =  h1->GetBinContent(b);
    h1->SetBinContent(b,sqrt(tmp1*tmp1 + xx*xx));
    
  }
    
}


void subtractHistByQuadrature(TH1 *h1, float xx){
  
  int nbin1 = h1->GetNbinsX();

  
  for(int b=1; b<=nbin1; b++){
    
    float tmp1 =  h1->GetBinContent(b);
    if(tmp1 > xx){
      h1->SetBinContent(b,sqrt(tmp1*tmp1 - xx*xx));
    }else{
      ///h1->SetBinContent(b,sqrt(tmp1*tmp1 - xx*xx));
    }
  }
    
}


void subtractTwoHistQuadrature(TH1 *h1, TH1 *h2, TH1 *h3){
  
  int nbin1 = h1->GetNbinsX();
  int nbin2 = h2->GetNbinsX();

  if( nbin1 != nbin2) {
    cout<<" subtractTwoHistQuadrature nbin " << nbin1 <<" "<<nbin2 <<endl; 
  }
  
  for(int b=1; b<=nbin1; b++){
    
    float tmp1 =  h1->GetBinContent(b);
    float tmp1err =  h1->GetBinError(b);
    float tmp2 =  h2->GetBinContent(b);
    float tmp2err =  h2->GetBinError(b);
    
    float diff = -1; 
    float diffErr = 0; 
    if( tmp1 >= tmp2){
      diff = sqrt( tmp1*tmp1 - tmp2*tmp2);
      diffErr = sqrt(  pow( tmp1/diff*tmp1err,2) +  pow( tmp2/diff*tmp2err,2) );

      //cout<<"check " << tmp1 <<" "<< tmp2 <<" "<< diff <<endl; 

    }
    h3->SetBinContent(b,diff);
    h3->SetBinError(b,diffErr);
  }
  
  
}


double one_sigma_up(std::vector<double> limits) {
 double v_one_sigma_up = 0;
 sort(limits.begin(),limits.end());
 int samplet_size = limits.size();
 int v_one_sigma_up_index = 0;
 v_one_sigma_up_index = samplet_size * 841 / 1000;
 v_one_sigma_up = limits[v_one_sigma_up_index-1];
 return v_one_sigma_up;
}

double one_sigma_down(std::vector<double> limits) {
 double v_one_sigma_down = 0;
 sort(limits.begin(),limits.end());
 int samplet_size = limits.size();
 int v_one_sigma_down_index = 0;
 v_one_sigma_down_index = samplet_size * 159 / 1000;
 v_one_sigma_down = limits[v_one_sigma_down_index];
 return v_one_sigma_down;
}

double two_sigma_up(std::vector<double> limits) {
 double v_two_sigma_up = 0;
 sort(limits.begin(),limits.end());
 int samplet_size = limits.size();
 int v_two_sigma_up_index = 0;
 v_two_sigma_up_index = samplet_size * 979 / 1000;
 v_two_sigma_up = limits[v_two_sigma_up_index-1];
 return v_two_sigma_up;
}

double two_sigma_down(std::vector<double> limits) {
 double v_two_sigma_down = 0;
 sort(limits.begin(),limits.end());
 int samplet_size = limits.size();
 int v_two_sigma_down_index = 0;
 v_two_sigma_down_index = samplet_size * 21 / 1000;
 v_two_sigma_down = limits[v_two_sigma_down_index];
 return v_two_sigma_down;
}


float getIntegralErrorHistogram(TH1F *h1){
  
  float err = 0; 
  for(int b=1; b<= h1->GetNbinsX(); b++){
    err += pow( h1->GetBinError(b),2); 
    
  }
  return sqrt(err);
}


void getIntegralAndWtErrorHistogram(TH1F *h1, int startBin, float res[]){
  
  float tot = 0; 
  float err = 0; 
  for(int b=startBin; b<= h1->GetNbinsX(); b++){
    err += pow( h1->GetBinError(b),2); 
    tot +=  h1->GetBinContent(b);
  }
  
  res[0] = tot; 
  res[1] = sqrt(err);
  
  
}


void printbinContentHistogram(TH1 *h1){
  
  for(int b=1; b<= h1->GetNbinsX(); b++){
    cout<<"bin: "<< h1->GetBinLowEdge(b)<<"to"<<  h1->GetBinLowEdge(b)+h1->GetBinWidth(b) <<" "<< h1->GetBinContent(b)<<"+/-"<< h1->GetBinError(b)<<endl; 
  }
  
}


void calcRatioBayes(float n1, float n2, float &r, float &rup, float &rlow){
  
  //compute using Bayes                                                                                                                                 
  TH1D h1("dummy1","",1,1,2);
  h1.SetBinContent(1,n1);
  
  TH1D h2("dummy2","",1,1,2);
  h2.SetBinContent(1,n2);
  
  TGraphAsymmErrors g;
  g.BayesDivide(&h1,&h2);
  r = g.GetY()[0];
  rup = g.GetErrorYhigh(0);
  rlow = g.GetErrorYlow(0);
  
}


double ErrorInProduct(double x, double errx, double y, 
                      double erry, double corr) {
   double xFrErr = errx/x;
   double yFrErr = erry/y;
   return sqrt(pow(xFrErr,2) + pow(yFrErr,2) + 2.0*corr*xFrErr*yFrErr)*x*y;
}


void  getErrorInRatio(double x, double errx, double y, 
                      double erry, double res[]) {

  if( x==0 || y==0){
    cout<<"getErrorInRatio x/y" << x <<" "<<y<<endl; 
    exit(1);
  }
  
  res[0] = x/y; 
  res[1] = x/y *sqrt( pow(errx/x,2) + pow(erry/y,2));

}


void get3Max(double x1, double x2, double x3,double res[]){
  double max = x1;
  res[1] = 0;

  if(max<x2) {
    max = x2;
    res[1] = 1;
  }
  if(max<x3){
    max = x3;
    res[1] = 2;
  }

  res[0] = max;
  //  return max;

}


float getMaxOfThree(float x1, float x2, float x3){
  
  float xmax = x1 > x2 ? x1: x2; 
  xmax = xmax > x3 ? xmax: x3; 
  return xmax; 
  
  
}


void generateToyDataFromHist(TH1* h1, TH1 *h2){
  
  int nbinx1 = h1->GetNbinsX();
  int nbinx2 = h2->GetNbinsX();
  
  if(nbinx1 != nbinx2){
    cout<<" generateToyDataFromHist " << nbinx1 <<" "<<  nbinx2 <<endl; 
    exit(1);
  }
  TRandom3 *grand = new TRandom3(0);
  for(int b=1 ; b<= nbinx1; b++){
    float tmp = h1->GetBinContent(b);
    float tmp1 = grand->Poisson(tmp);
    h2->SetBinContent(b,tmp1);
  }
  
  
}


void calcTotalErorrPDF(int npdf, double X0,double X0Err,double XP[100],double XPErr[100],double XM[100],double XMErr[100],double results[]){
  
  
  double dXP = 0; 
  double dXM = 0; 
  double simpleDX = 0; 
  
  double dsimpleDX = 0; 
  double ddXP = 0; 
  double ddXM = 0; 
  
  double res[10];
  
  if(npdf%2 !=0) {
    cout<<" calcTotalErorrPDF npdf: "<<npdf <<endl; 
    exit(1);
  }
    
  
  for( int j=1; j<= npdf/2; j++){
    
    get3Max(XP[j]-X0,XM[j]-X0,0,res);
    double dxP = res[0];
    if( fabs(res[1]-0)<0.01){
      ddXP += 4*pow(XP[j]-X0,2)*(XPErr[j]*XPErr[j] + X0Err*X0Err);
    }else if( fabs(res[1]-1)<0.01){
      ddXP += 4*pow(XM[j]-X0,2)*(XMErr[j]*XMErr[j] + X0Err*X0Err);
    }
    
    get3Max(X0-XP[j],X0-XM[j],0,res);
    double dxM = res[0];
    if( fabs(res[1]-0)<0.01){
      ddXM += 4*pow(XP[j]-X0,2)*(XPErr[j]*XPErr[j] + X0Err*X0Err);
    }else if( fabs(res[1]-1)<0.01){
      ddXM += 4*pow(XM[j]-X0,2)*(XMErr[j]*XMErr[j] + X0Err*X0Err);
    }
    
    dXP += dxP*dxP; 
    dXM += dxM*dxM; 
    
    simpleDX +=  pow(XP[j]-XM[j],2);

    dsimpleDX += pow(XP[j]-XM[j],2)*(XPErr[j]*XPErr[j]+XMErr[j]*XMErr[j]);
    
    
    ///cout<<"checkme: "<<j<<"  "<< dXP <<" "<< dXM <<endl; 
    

  }


  dXP = sqrt(dXP);
  dXM = sqrt(dXM);

  ddXP = sqrt(ddXP);

  if( dXP>0){
    ddXP = ddXP/(2*dXP);
  }
  ddXM = sqrt(ddXM);
  if( dXM >0){
    ddXM = ddXM/(2*dXM);
  }
  
  
  simpleDX = 0.5* sqrt(simpleDX);
  dsimpleDX =  sqrt(dsimpleDX);
  if(simpleDX >0){
    dsimpleDX = dsimpleDX/(2*simpleDX);
  }
  
  ////asymmetry errors
  results[0] = dXP; 
  results[1] = dXM; 
  ///master formula
  results[2] = simpleDX; 
  
  ////error of above
  results[3] = ddXP; 
  results[4] = ddXM; 
  results[5] = dsimpleDX; 
  
  
}

int mostProb(float b){
  
  int nbest = int(b+1); 
  float prob_max = 0; 
  for(int n = int(b-1) ; n<= int(b+1); n++){
    if( n<0) n =0; 


    float prob = TMath::Poisson(n,b);

    //    cout<<prob<<" "<<n<<" "<<b<<endl;
    
    if( prob > prob_max){
      prob_max = prob; 
      nbest = n; 
    }
  }
  return nbest;
  
}


float getBinContentHistogram(TH1F *hhtmp, float pt){
  
  for(int b=1; b<= hhtmp->GetNbinsX(); b++){
    if( pt >= hhtmp->GetBinLowEdge(b) && pt < hhtmp->GetBinLowEdge(b) +  hhtmp->GetBinWidth(b)){
      return  hhtmp->GetBinContent(b);
    }
  }
  return -999; 
  
  
}


void weightedAeverageTwoHistogram(TH1F *hhtmp1, TH1F *hhtmp2, TH1F *hhres){
  
  if( hhtmp1->GetNbinsX() != hhtmp2->GetNbinsX()){
    cout<<"warning  weightedAeverageTwoHistogram  "<< hhtmp1->GetNbinsX() <<" "<< hhtmp2->GetNbinsX() <<endl; 
    exit(1); 
  }

  for(int b=1; b<= hhtmp1->GetNbinsX(); b++){
    float tmp1 =  hhtmp1->GetBinContent(b);
    float tmp2 =  hhtmp2->GetBinContent(b);
    float tmp1e =  hhtmp1->GetBinError(b);
    float tmp2e =  hhtmp2->GetBinError(b);
    
    float av = 0; 
    float ave = 0; 
    if( tmp1 > 0 && tmp2 >0){
      av = (tmp1/(tmp1e*tmp1e) + tmp2/ (tmp2e*tmp2e)) / ( 1./(tmp1e*tmp1e) + 1./(tmp2e*tmp2e));
      ave =  sqrt(1/ (1./(tmp1e*tmp1e) + 1./(tmp2e*tmp2e)) ); 
    }else if (tmp1 > 0 && tmp2 ==0){
      av = tmp1; 
      ave = tmp1e; 
    }else if( tmp1 ==0 && tmp2 >0){
      av = tmp2; 
      ave = tmp2e; 
    }else{
      cout<<"no data for combine.."<< b<<endl; 
    }
    
    hhres->SetBinContent(b,av);
    hhres->SetBinError(b,ave);
    
  }

}


float getFractionOfHistogram(TH1F *hhtmp, double xmin, double xmax){
  
  double sum = hhtmp->Integral(); 
  if( sum <=0) {
    cout<<" getFractionOfHistogram zero.." << sum << " "<< hhtmp->GetEntries();
    exit(1); 
  }
  
  int nbinsx = hhtmp->GetNbinsX();
  int bin1 =1;  
  int bin2 = nbinsx; 
  
  for(int b=1; b<= nbinsx; b++){
    if( xmin >= hhtmp->GetBinLowEdge(b) && xmin -hhtmp->GetBinLowEdge(b) - hhtmp->GetBinWidth(b) < 1E-6 ){
      bin1 = b; 
    }
    if( xmax > hhtmp->GetBinLowEdge(b) && xmax - hhtmp->GetBinLowEdge(b) - hhtmp->GetBinWidth(b) < 1E-6 ){
      bin2 = b; 
    }
  }
  
  ///cout<<"checkme"<<bin1 <<" "<<bin2 <<" "<< hhtmp->GetBinLowEdge(bin1)  <<" "<< xmax - hhtmp->GetBinLowEdge(bin2) - hhtmp->GetBinWidth(bin2) <<" "<< xmax - hhtmp->GetBinLowEdge(bin2-1) - hhtmp->GetBinWidth(bin2-1)<<endl; 
  
  return hhtmp->Integral(bin1,bin2) / sum; 
  
  
}


void cloneHistogramWithWeight(TH1 *hhtmp1, TH1* hhtmp2, float wt){
  
  int nbins = hhtmp1->GetNbinsX();
  if( nbins != hhtmp2->GetNbinsX()){
    cout<<"error.cloneHistogramWithWeight  "<<nbins <<" "<< hhtmp2->GetNbinsX() << endl; 
    exit(1); 
  }
  
  for(int j=1; j<= nbins+1; j++){ ///bins+1 for overflowbins
    
    float tmp = hhtmp1->GetBinContent(j) * wt; 

    hhtmp2->SetBinContent(j,tmp);

  }
  
  
}


void cloneHistogramWithWeightv1(TH1 *hhtmp1, TH1* hhtmp2, float wt){
  
  int nbins = hhtmp1->GetNbinsX();
  if( nbins >  hhtmp2->GetNbinsX()){
    cout<<"error.cloneHistogramWithWeightv1  "<<nbins <<" "<< hhtmp2->GetNbinsX() << endl; 
    exit(1); 
  }
  
  for(int j=1; j<= nbins; j++){
    
    float tmp = hhtmp1->GetBinContent(j) * wt; 

    hhtmp2->SetBinContent(j,tmp);

  }
  
  
}


///check which bin the value in histogram
int getBinOfHistogram(TH1 *hh, float tmp){
  
  int nbins = hh->GetNbinsX();
  
  if( tmp < hh->GetXaxis()->GetXmin() || tmp >= hh->GetXaxis()->GetXmax()) return -1; 
  
  for(int b = 1; b<= nbins; b++){
    
    if( tmp >= hh->GetBinLowEdge(b) && tmp < hh->GetBinLowEdge(b) +  hh->GetBinWidth(b)) return b; 
    
  }
  
  cout<<"error.. getBinOfHistogram "<<tmp <<endl;
  
  return -1; 
  
  
}


///check which bin the value in histogram
int getBinOfHistogramv1(TH1 *hh, float tmp){
  
  int nbins = hh->GetNbinsX();
  
  if( tmp < hh->GetXaxis()->GetXmin() ) return 1; 
  if( tmp >= hh->GetXaxis()->GetXmax() ) return nbins; 
    
  for(int b = 1; b<= nbins; b++){
    
    if( tmp >= hh->GetBinLowEdge(b) && tmp < hh->GetBinLowEdge(b) +  hh->GetBinWidth(b)) return b; 
    
  }
  
  cout<<"error.. getBinOfHistogram "<<tmp <<endl;
  exit(1);
  
  return -1; 
  
  
}


double BinP(int N, double p, int x1, int x2) {
   double q=p/(1-p); 
   int k=0; 
   double v = 1; 
   double s=0; 
   double tot=0.0;
    while(k<=N) {
       tot=tot+v;
       if(k>=x1 & k<=x2) { s=s+v; }
       if(tot>1e30){s=s/1e30; tot=tot/1e30; v=v/1e30;}
       k=k+1; 
       v=v*q*(N+1-k)/k;
    }
    return s/tot;
}


void ClopperPearsonLimits(double numerator, double denominator, double &ratio,
double &lowerLimit, double &upperLimit, const double CL_low=1.0, 
const double CL_high=1.0) 
{  
  //Confidence intervals are in the units of \sigma.
    

  ratio = numerator/denominator;
   
// first get the lower limit
   if(numerator==0)   lowerLimit = 0.0; 
   else { 
      double v=ratio/2; 
      double vsL=0; 
      double vsH=ratio; 
      double p=CL_low/100;
      while((vsH-vsL)>1e-5) { 
        if(BinP(int(denominator),v,int(numerator),int(denominator))>p) 
         { vsH=v; v=(vsL+v)/2; } 
         else { vsL=v; v=(v+vsH)/2; } 
      }
      lowerLimit = v; 
   }
   
// now get the upper limit
   if(numerator==denominator) upperLimit = 1.0;
   else { 
      double v=(1+ratio)/2; 
      double vsL=ratio; 
      double vsH=1; 
      double p=CL_high/100;
      while((vsH-vsL)>1e-5) { 
        if(BinP(int(denominator),v,0,int(numerator))<p) { vsH=v; v=(vsL+v)/2; } 
         else { vsL=v; v=(v+vsH)/2; } 
      }
      upperLimit = v;
   }
}


Revision:	1.1
Committed:	Fri Sep 30 15:23:49 2011 UTC (13 years, 7 months ago) by yangyong
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V01-00-01, V01-00-00, HEAD
Log Message:	z selection code
#	Content
1
2	float getMaximumY(TH1 h1, TH1 h2){
3	return h1->GetMaximum() > h2->GetMaximum() ? h1->GetMaximum() : h2->GetMaximum();
4	}
5
6	int getBinNumberwithLabel(TH1 *h1, string bnam){
7	int nbin1 = h1->GetNbinsX();
8
9
10	for(int b=1; b<=nbin1; b++){
11	if( h1->GetXaxis()->GetBinLabel(b) == bnam){
12	return b;
13	}
14	}
15	cout<<"not found " << bnam.c_str()<<endl;
16	return -1;
17	}
18
19
20
21
22	void AddHistByQuadrature(TH1 *h1, float xx){
23
24	int nbin1 = h1->GetNbinsX();
25
26
27	for(int b=1; b<=nbin1; b++){
28
29	float tmp1 = h1->GetBinContent(b);
30	h1->SetBinContent(b,sqrt(tmp1tmp1 + xxxx));
31
32	}
33
34	}
35
36
37
38	void subtractHistByQuadrature(TH1 *h1, float xx){
39
40	int nbin1 = h1->GetNbinsX();
41
42
43	for(int b=1; b<=nbin1; b++){
44
45	float tmp1 = h1->GetBinContent(b);
46	if(tmp1 > xx){
47	h1->SetBinContent(b,sqrt(tmp1tmp1 - xxxx));
48	}else{
49	///h1->SetBinContent(b,sqrt(tmp1tmp1 - xxxx));
50	}
51	}
52
53	}
54
55
56
57
58	void subtractTwoHistQuadrature(TH1 h1, TH1 h2, TH1 *h3){
59
60	int nbin1 = h1->GetNbinsX();
61	int nbin2 = h2->GetNbinsX();
62
63	if( nbin1 != nbin2) {
64	cout<<" subtractTwoHistQuadrature nbin " << nbin1 <<" "<<nbin2 <<endl;
65	}
66
67	for(int b=1; b<=nbin1; b++){
68
69	float tmp1 = h1->GetBinContent(b);
70	float tmp1err = h1->GetBinError(b);
71	float tmp2 = h2->GetBinContent(b);
72	float tmp2err = h2->GetBinError(b);
73
74	float diff = -1;
75	float diffErr = 0;
76	if( tmp1 >= tmp2){
77	diff = sqrt( tmp1tmp1 - tmp2tmp2);
78	diffErr = sqrt( pow( tmp1/difftmp1err,2) + pow( tmp2/difftmp2err,2) );
79
80	//cout<<"check " << tmp1 <<" "<< tmp2 <<" "<< diff <<endl;
81
82	}
83	h3->SetBinContent(b,diff);
84	h3->SetBinError(b,diffErr);
85	}
86
87
88	}
89
90
91
92
93	double one_sigma_up(std::vector<double> limits) {
94	double v_one_sigma_up = 0;
95	sort(limits.begin(),limits.end());
96	int samplet_size = limits.size();
97	int v_one_sigma_up_index = 0;
98	v_one_sigma_up_index = samplet_size * 841 / 1000;
99	v_one_sigma_up = limits[v_one_sigma_up_index-1];
100	return v_one_sigma_up;
101	}
102
103	double one_sigma_down(std::vector<double> limits) {
104	double v_one_sigma_down = 0;
105	sort(limits.begin(),limits.end());
106	int samplet_size = limits.size();
107	int v_one_sigma_down_index = 0;
108	v_one_sigma_down_index = samplet_size * 159 / 1000;
109	v_one_sigma_down = limits[v_one_sigma_down_index];
110	return v_one_sigma_down;
111	}
112
113	double two_sigma_up(std::vector<double> limits) {
114	double v_two_sigma_up = 0;
115	sort(limits.begin(),limits.end());
116	int samplet_size = limits.size();
117	int v_two_sigma_up_index = 0;
118	v_two_sigma_up_index = samplet_size * 979 / 1000;
119	v_two_sigma_up = limits[v_two_sigma_up_index-1];
120	return v_two_sigma_up;
121	}
122
123	double two_sigma_down(std::vector<double> limits) {
124	double v_two_sigma_down = 0;
125	sort(limits.begin(),limits.end());
126	int samplet_size = limits.size();
127	int v_two_sigma_down_index = 0;
128	v_two_sigma_down_index = samplet_size * 21 / 1000;
129	v_two_sigma_down = limits[v_two_sigma_down_index];
130	return v_two_sigma_down;
131	}
132
133
134
135
136
137
138	float getIntegralErrorHistogram(TH1F *h1){
139
140	float err = 0;
141	for(int b=1; b<= h1->GetNbinsX(); b++){
142	err += pow( h1->GetBinError(b),2);
143
144	}
145	return sqrt(err);
146	}
147
148
149	void getIntegralAndWtErrorHistogram(TH1F *h1, int startBin, float res[]){
150
151	float tot = 0;
152	float err = 0;
153	for(int b=startBin; b<= h1->GetNbinsX(); b++){
154	err += pow( h1->GetBinError(b),2);
155	tot += h1->GetBinContent(b);
156	}
157
158	res[0] = tot;
159	res[1] = sqrt(err);
160
161
162	}
163
164
165
166	void printbinContentHistogram(TH1 *h1){
167
168	for(int b=1; b<= h1->GetNbinsX(); b++){
169	cout<<"bin: "<< h1->GetBinLowEdge(b)<<"to"<< h1->GetBinLowEdge(b)+h1->GetBinWidth(b) <<" "<< h1->GetBinContent(b)<<"+/-"<< h1->GetBinError(b)<<endl;
170	}
171
172	}
173
174
175
176
177	void calcRatioBayes(float n1, float n2, float &r, float &rup, float &rlow){
178
179	//compute using Bayes
180	TH1D h1("dummy1","",1,1,2);
181	h1.SetBinContent(1,n1);
182
183	TH1D h2("dummy2","",1,1,2);
184	h2.SetBinContent(1,n2);
185
186	TGraphAsymmErrors g;
187	g.BayesDivide(&h1,&h2);
188	r = g.GetY()[0];
189	rup = g.GetErrorYhigh(0);
190	rlow = g.GetErrorYlow(0);
191
192	}
193
194
195
196
197	double ErrorInProduct(double x, double errx, double y,
198	double erry, double corr) {
199	double xFrErr = errx/x;
200	double yFrErr = erry/y;
201	return sqrt(pow(xFrErr,2) + pow(yFrErr,2) + 2.0corrxFrErryFrErr)x*y;
202	}
203
204
205
206	void getErrorInRatio(double x, double errx, double y,
207	double erry, double res[]) {
208
209	if( x==0 \|\| y==0){
210	cout<<"getErrorInRatio x/y" << x <<" "<<y<<endl;
211	exit(1);
212	}
213
214	res[0] = x/y;
215	res[1] = x/y *sqrt( pow(errx/x,2) + pow(erry/y,2));
216
217	}
218
219
220
221	void get3Max(double x1, double x2, double x3,double res[]){
222	double max = x1;
223	res[1] = 0;
224
225	if(max<x2) {
226	max = x2;
227	res[1] = 1;
228	}
229	if(max<x3){
230	max = x3;
231	res[1] = 2;
232	}
233
234	res[0] = max;
235	// return max;
236
237	}
238
239
240
241	float getMaxOfThree(float x1, float x2, float x3){
242
243	float xmax = x1 > x2 ? x1: x2;
244	xmax = xmax > x3 ? xmax: x3;
245	return xmax;
246
247
248	}
249
250
251	void generateToyDataFromHist(TH1* h1, TH1 *h2){
252
253	int nbinx1 = h1->GetNbinsX();
254	int nbinx2 = h2->GetNbinsX();
255
256	if(nbinx1 != nbinx2){
257	cout<<" generateToyDataFromHist " << nbinx1 <<" "<< nbinx2 <<endl;
258	exit(1);
259	}
260	TRandom3 *grand = new TRandom3(0);
261	for(int b=1 ; b<= nbinx1; b++){
262	float tmp = h1->GetBinContent(b);
263	float tmp1 = grand->Poisson(tmp);
264	h2->SetBinContent(b,tmp1);
265	}
266
267
268	}
269
270
271	void calcTotalErorrPDF(int npdf, double X0,double X0Err,double XP[100],double XPErr[100],double XM[100],double XMErr[100],double results[]){
272
273
274	double dXP = 0;
275	double dXM = 0;
276	double simpleDX = 0;
277
278	double dsimpleDX = 0;
279	double ddXP = 0;
280	double ddXM = 0;
281
282	double res[10];
283
284	if(npdf%2 !=0) {
285	cout<<" calcTotalErorrPDF npdf: "<<npdf <<endl;
286	exit(1);
287	}
288
289
290	for( int j=1; j<= npdf/2; j++){
291
292	get3Max(XP[j]-X0,XM[j]-X0,0,res);
293	double dxP = res[0];
294	if( fabs(res[1]-0)<0.01){
295	ddXP += 4pow(XP[j]-X0,2)(XPErr[j]XPErr[j] + X0ErrX0Err);
296	}else if( fabs(res[1]-1)<0.01){
297	ddXP += 4pow(XM[j]-X0,2)(XMErr[j]XMErr[j] + X0ErrX0Err);
298	}
299
300	get3Max(X0-XP[j],X0-XM[j],0,res);
301	double dxM = res[0];
302	if( fabs(res[1]-0)<0.01){
303	ddXM += 4pow(XP[j]-X0,2)(XPErr[j]XPErr[j] + X0ErrX0Err);
304	}else if( fabs(res[1]-1)<0.01){
305	ddXM += 4pow(XM[j]-X0,2)(XMErr[j]XMErr[j] + X0ErrX0Err);
306	}
307
308	dXP += dxP*dxP;
309	dXM += dxM*dxM;
310
311	simpleDX += pow(XP[j]-XM[j],2);
312
313	dsimpleDX += pow(XP[j]-XM[j],2)(XPErr[j]XPErr[j]+XMErr[j]*XMErr[j]);
314
315
316	///cout<<"checkme: "<<j<<" "<< dXP <<" "<< dXM <<endl;
317
318
319	}
320
321
322	dXP = sqrt(dXP);
323	dXM = sqrt(dXM);
324
325	ddXP = sqrt(ddXP);
326
327	if( dXP>0){
328	ddXP = ddXP/(2*dXP);
329	}
330	ddXM = sqrt(ddXM);
331	if( dXM >0){
332	ddXM = ddXM/(2*dXM);
333	}
334
335
336
337
338	simpleDX = 0.5* sqrt(simpleDX);
339	dsimpleDX = sqrt(dsimpleDX);
340	if(simpleDX >0){
341	dsimpleDX = dsimpleDX/(2*simpleDX);
342	}
343
344	////asymmetry errors
345	results[0] = dXP;
346	results[1] = dXM;
347	///master formula
348	results[2] = simpleDX;
349
350	////error of above
351	results[3] = ddXP;
352	results[4] = ddXM;
353	results[5] = dsimpleDX;
354
355
356	}
357
358	int mostProb(float b){
359
360	int nbest = int(b+1);
361	float prob_max = 0;
362	for(int n = int(b-1) ; n<= int(b+1); n++){
363	if( n<0) n =0;
364
365
366	float prob = TMath::Poisson(n,b);
367
368	// cout<<prob<<" "<<n<<" "<<b<<endl;
369
370	if( prob > prob_max){
371	prob_max = prob;
372	nbest = n;
373	}
374	}
375	return nbest;
376
377	}
378
379
380	float getBinContentHistogram(TH1F *hhtmp, float pt){
381
382	for(int b=1; b<= hhtmp->GetNbinsX(); b++){
383	if( pt >= hhtmp->GetBinLowEdge(b) && pt < hhtmp->GetBinLowEdge(b) + hhtmp->GetBinWidth(b)){
384	return hhtmp->GetBinContent(b);
385	}
386	}
387	return -999;
388
389
390	}
391
392
393
394
395
396
397	void weightedAeverageTwoHistogram(TH1F hhtmp1, TH1F hhtmp2, TH1F *hhres){
398
399	if( hhtmp1->GetNbinsX() != hhtmp2->GetNbinsX()){
400	cout<<"warning weightedAeverageTwoHistogram "<< hhtmp1->GetNbinsX() <<" "<< hhtmp2->GetNbinsX() <<endl;
401	exit(1);
402	}
403
404	for(int b=1; b<= hhtmp1->GetNbinsX(); b++){
405	float tmp1 = hhtmp1->GetBinContent(b);
406	float tmp2 = hhtmp2->GetBinContent(b);
407	float tmp1e = hhtmp1->GetBinError(b);
408	float tmp2e = hhtmp2->GetBinError(b);
409
410	float av = 0;
411	float ave = 0;
412	if( tmp1 > 0 && tmp2 >0){
413	av = (tmp1/(tmp1etmp1e) + tmp2/ (tmp2etmp2e)) / ( 1./(tmp1etmp1e) + 1./(tmp2etmp2e));
414	ave = sqrt(1/ (1./(tmp1etmp1e) + 1./(tmp2etmp2e)) );
415	}else if (tmp1 > 0 && tmp2 ==0){
416	av = tmp1;
417	ave = tmp1e;
418	}else if( tmp1 ==0 && tmp2 >0){
419	av = tmp2;
420	ave = tmp2e;
421	}else{
422	cout<<"no data for combine.."<< b<<endl;
423	}
424
425	hhres->SetBinContent(b,av);
426	hhres->SetBinError(b,ave);
427
428	}
429
430	}
431
432
433
434
435	float getFractionOfHistogram(TH1F *hhtmp, double xmin, double xmax){
436
437	double sum = hhtmp->Integral();
438	if( sum <=0) {
439	cout<<" getFractionOfHistogram zero.." << sum << " "<< hhtmp->GetEntries();
440	exit(1);
441	}
442
443	int nbinsx = hhtmp->GetNbinsX();
444	int bin1 =1;
445	int bin2 = nbinsx;
446
447	for(int b=1; b<= nbinsx; b++){
448	if( xmin >= hhtmp->GetBinLowEdge(b) && xmin -hhtmp->GetBinLowEdge(b) - hhtmp->GetBinWidth(b) < 1E-6 ){
449	bin1 = b;
450	}
451	if( xmax > hhtmp->GetBinLowEdge(b) && xmax - hhtmp->GetBinLowEdge(b) - hhtmp->GetBinWidth(b) < 1E-6 ){
452	bin2 = b;
453	}
454	}
455
456	///cout<<"checkme"<<bin1 <<" "<<bin2 <<" "<< hhtmp->GetBinLowEdge(bin1) <<" "<< xmax - hhtmp->GetBinLowEdge(bin2) - hhtmp->GetBinWidth(bin2) <<" "<< xmax - hhtmp->GetBinLowEdge(bin2-1) - hhtmp->GetBinWidth(bin2-1)<<endl;
457
458	return hhtmp->Integral(bin1,bin2) / sum;
459
460
461	}
462
463
464
465
466	void cloneHistogramWithWeight(TH1 hhtmp1, TH1 hhtmp2, float wt){
467
468	int nbins = hhtmp1->GetNbinsX();
469	if( nbins != hhtmp2->GetNbinsX()){
470	cout<<"error.cloneHistogramWithWeight "<<nbins <<" "<< hhtmp2->GetNbinsX() << endl;
471	exit(1);
472	}
473
474	for(int j=1; j<= nbins+1; j++){ ///bins+1 for overflowbins
475
476	float tmp = hhtmp1->GetBinContent(j) * wt;
477
478	hhtmp2->SetBinContent(j,tmp);
479
480	}
481
482
483	}
484
485
486
487
488	void cloneHistogramWithWeightv1(TH1 hhtmp1, TH1 hhtmp2, float wt){
489
490	int nbins = hhtmp1->GetNbinsX();
491	if( nbins > hhtmp2->GetNbinsX()){
492	cout<<"error.cloneHistogramWithWeightv1 "<<nbins <<" "<< hhtmp2->GetNbinsX() << endl;
493	exit(1);
494	}
495
496	for(int j=1; j<= nbins; j++){
497
498	float tmp = hhtmp1->GetBinContent(j) * wt;
499
500	hhtmp2->SetBinContent(j,tmp);
501
502	}
503
504
505	}
506
507
508
509
510	///check which bin the value in histogram
511	int getBinOfHistogram(TH1 *hh, float tmp){
512
513	int nbins = hh->GetNbinsX();
514
515	if( tmp < hh->GetXaxis()->GetXmin() \|\| tmp >= hh->GetXaxis()->GetXmax()) return -1;
516
517	for(int b = 1; b<= nbins; b++){
518
519	if( tmp >= hh->GetBinLowEdge(b) && tmp < hh->GetBinLowEdge(b) + hh->GetBinWidth(b)) return b;
520
521	}
522
523	cout<<"error.. getBinOfHistogram "<<tmp <<endl;
524
525	return -1;
526
527
528	}
529
530
531
532
533
534	///check which bin the value in histogram
535	int getBinOfHistogramv1(TH1 *hh, float tmp){
536
537	int nbins = hh->GetNbinsX();
538
539	if( tmp < hh->GetXaxis()->GetXmin() ) return 1;
540	if( tmp >= hh->GetXaxis()->GetXmax() ) return nbins;
541
542	for(int b = 1; b<= nbins; b++){
543
544	if( tmp >= hh->GetBinLowEdge(b) && tmp < hh->GetBinLowEdge(b) + hh->GetBinWidth(b)) return b;
545
546	}
547
548	cout<<"error.. getBinOfHistogram "<<tmp <<endl;
549	exit(1);
550
551	return -1;
552
553
554	}
555
556
557
558
559
560
561	double BinP(int N, double p, int x1, int x2) {
562	double q=p/(1-p);
563	int k=0;
564	double v = 1;
565	double s=0;
566	double tot=0.0;
567	while(k<=N) {
568	tot=tot+v;
569	if(k>=x1 & k<=x2) { s=s+v; }
570	if(tot>1e30){s=s/1e30; tot=tot/1e30; v=v/1e30;}
571	k=k+1;
572	v=vq(N+1-k)/k;
573	}
574	return s/tot;
575	}
576
577
578
579
580	void ClopperPearsonLimits(double numerator, double denominator, double &ratio,
581	double &lowerLimit, double &upperLimit, const double CL_low=1.0,
582	const double CL_high=1.0)
583	{
584	//Confidence intervals are in the units of \sigma.
585
586
587	ratio = numerator/denominator;
588
589	// first get the lower limit
590	if(numerator==0) lowerLimit = 0.0;
591	else {
592	double v=ratio/2;
593	double vsL=0;
594	double vsH=ratio;
595	double p=CL_low/100;
596	while((vsH-vsL)>1e-5) {
597	if(BinP(int(denominator),v,int(numerator),int(denominator))>p)
598	{ vsH=v; v=(vsL+v)/2; }
599	else { vsL=v; v=(v+vsH)/2; }
600	}
601	lowerLimit = v;
602	}
603
604	// now get the upper limit
605	if(numerator==denominator) upperLimit = 1.0;
606	else {
607	double v=(1+ratio)/2;
608	double vsL=ratio;
609	double vsH=1;
610	double p=CL_high/100;
611	while((vsH-vsL)>1e-5) {
612	if(BinP(int(denominator),v,0,int(numerator))<p) { vsH=v; v=(vsL+v)/2; }
613	else { vsL=v; v=(v+vsH)/2; }
614	}
615	upperLimit = v;
616	}
617	}
618
619
620
621
622
623