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
#	User	Rev	Content
1	yangyong	1.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