yangyong/combineICEB/usefullcode.cc

TRandom3 *rgeng_;

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; 
}


float getGaussian(float mean, float sigma, float min,float max){
  
  if( rgeng_ == 0){
    rgeng_ = new TRandom3(0);
  }
  
  float tmp = rgeng_->Gaus(mean,sigma);
  while( tmp <= min || tmp >= max ){
    tmp = rgeng_->Gaus(mean,sigma);
  }
  return tmp; 
  
}


double interpolatedTwoBinTH1(TH1 *h1, double c){
  
  int nbin1 = h1->GetNbinsX();
  for(int b=1; b<nbin1; b++){
    double c1 = h1->GetBinCenter(b);
    double c2 = h1->GetBinCenter(b+1);
    if( c >= c1 && c <= c2){
      double v1 = h1->GetBinContent(b);
      double v2 = h1->GetBinContent(b+1);
      double k = (v2-v1)/(c2-c1);
      return v1 + k * ( c-c1);
    }
  }
  return 0; 
  
}


void NormHistogram(TH1 *h1){
  float total = h1->GetEntries(); 
  
  if( total <=0){
    cout<<" empty histogram? " << total <<endl; 
  }
  
  h1->Scale(1./total);
  
}


void emptyHistogram( TH1 *h1){
  int nbin1 = h1->GetNbinsX();
  for(int b=1; b<=nbin1; b++){
    h1->SetBinContent(b,0);
  }
  
}


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 scaleHistogramMCtoData(TH1 *hda, TH1 *hmc){
  
  double ndata = hda->Integral();
  double nmc = hmc->Integral();

  if( ndata >0 && nmc >0){
    hmc ->Scale( ndata/nmc);
  }
  

}


void scaleHistogramMCtoDatav1(TH1 *hda, TH1 *hmc, float xlow, float xhigh){
  
  double width = hda->GetBinWidth(1);
  int nbin_start =  int ( (xlow - hda->GetXaxis()->GetXmin())/ width) +1; 
  int nbin_end =  int ( (xhigh - hda->GetXaxis()->GetXmin())/ width); 
  
  double ndata = hda->Integral(nbin_start, nbin_end);
  double nmc = hmc->Integral(nbin_start,nbin_end);
  
  cout<<" nbin_start/end "<< nbin_start <<" "<< nbin_end <<" "<< ndata <<" "<<nmc <<endl; 
  
  if( ndata >0 && nmc >0){
    hmc ->Scale( ndata/nmc);
  }
  

}


void addTH2(TH2 *h1,TH2 *h2){
  
  int binx = h1->GetNbinsX();
  int biny = h1->GetNbinsY();
  
  if( binx !=  h2->GetNbinsX() || biny != h2->GetNbinsY()){
    cout<<"diff bins addTH2 ? " <<binx <<" "<< h2->GetNbinsX() <<" "<<  biny <<" "<< h2->GetNbinsY()<<endl; 
    return; 
  }
  for(int bx = 1; bx<= binx; bx++){
    for(int by = 1; by<= biny; by++){
      double tmp1 = h1->GetBinContent(bx,by);
      double tmp2 = h2->GetBinContent(bx,by);
      h2->SetBinContent(bx,by,tmp1+tmp2);
    }
  }
  
}


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;
}


TH1F *divideHistogramDataMC(TH1 *hdata, TH1 *hmc){
  
  TString hname = TString(hdata->GetName()) + TString(hmc->GetName());
  TH1F *htmp = (TH1F*)hdata->Clone(hname);
  
  for(int b=1; b<= hdata->GetNbinsX();b++){
    float r=0; 
    float rErr = 0; 
    if( hmc->GetBinContent(b) >0){
      r = hdata->GetBinContent(b) / hmc->GetBinContent(b);
      rErr = sqrt(r) / hmc->GetBinContent(b);
    }
    htmp->SetBinContent(b,r);
    htmp->SetBinError(b,rErr);
    
  }
  return htmp;
  
  
}


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);
  
}


TGraphAsymmErrors* dividingTwoHistogramNoErrBinContent(TH1 *h1, TH1 *h2){
  
  
  float x[1000];
  float exl[1000];
  float exh[1000];
    
  float y[1000];
  float yup[1000];
  float ydw[1000];
  

  if( h1->GetNbinsX()  != h2->GetNbinsX() ){
    cout<<" dividingTwoHistogramNoErrBinContent diff bins " << endl; 
    return 0 ; 
  }
  
  if( h1->GetNbinsX() > 1000){
    cout<<" dividingTwoHistogramNoErrBinContent too much bins " << endl; 
    return 0 ; 
  }

  int n =0; 
  
  for(int b=1; b<= h1->GetNbinsX(); b++){
    float b1 = h1->GetBinContent(b);
    float b2 = h2->GetBinContent(b);
    
    float r,rup,rdown;
    
    
    x[n] = h1->GetBinCenter(b);
    exl[n] = 0.5*h1->GetBinWidth(b);
    exh[n] = 0.5*h1->GetBinWidth(b);
    
    if( b1>0 && b2>0){
      if( b1<=b2){
        calcRatioBayes(b1,b2,r,rup,rdown);
      }else{
        calcRatioBayes(b2,b1,r,rup,rdown);
      }
    }else{
      r = 0; 
      rup = 0; 
      rdown = 0; 
    }
    
    
    y[n] = b1/b2; 
    yup[n] = rup;
    ydw[n] = rdown; 
        
    
    n ++; 
    
  }
  TGraphAsymmErrors *g1 = new TGraphAsymmErrors(n, x, y, exl, exh, ydw,yup);
  return g1; 
    

}


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++){
    double 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; 
  
  
}


Revision:	1.3
Committed:	Mon Aug 27 21:37:40 2012 UTC (12 years, 8 months ago) by yangyong
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V2012combv2c, HEAD
Changes since 1.2:	+0 -64 lines
Log Message:	* empty log message *
#	Content
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(tmp1tmp1 + xxxx));
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(tmp1tmp1 - xxxx));
106	}else{
107	///h1->SetBinContent(b,sqrt(tmp1tmp1 - xxxx));
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( tmp1tmp1 - tmp2tmp2);
187	diffErr = sqrt( pow( tmp1/difftmp1err,2) + pow( tmp2/difftmp2err,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.0corrxFrErryFrErr)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 += 4pow(XP[j]-X0,2)(XPErr[j]XPErr[j] + X0ErrX0Err);
488	}else if( fabs(res[1]-1)<0.01){
489	ddXP += 4pow(XM[j]-X0,2)(XMErr[j]XMErr[j] + X0ErrX0Err);
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 += 4pow(XP[j]-X0,2)(XPErr[j]XPErr[j] + X0ErrX0Err);
496	}else if( fabs(res[1]-1)<0.01){
497	ddXM += 4pow(XM[j]-X0,2)(XMErr[j]XMErr[j] + X0ErrX0Err);
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/(tmp1etmp1e) + tmp2/ (tmp2etmp2e)) / ( 1./(tmp1etmp1e) + 1./(tmp2etmp2e));
606	ave = sqrt(1/ (1./(tmp1etmp1e) + 1./(tmp2etmp2e)) );
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