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


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.2
Committed:	Tue Apr 10 19:41:12 2012 UTC (13 years, 1 month ago) by yangyong
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V2012combv2b, V2012combv2a, V2012combv2, V2012combv1, V2011combv1, pi0etaeb_laser20111122
Changes since 1.1:	+812 -0 lines
Log Message:	LaserTag: EcalLaserAPDPNRatios_data_20111122_158851_180363 and alpha Tag ( for Endcap Only) EcalLaserAlphas_lto420-620_progr_data_20111122
#	User	Rev	Content
1	yangyong	1.2	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
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=vq(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