TH1F* addhists(TH1* h1, TH1*h2, double a, double b)
{
 TH1F *h3 = new TH1F("h3", "h3", h1->GetNbinsX(), h1->GetXaxis()->GetXmin(), h1->GetXaxis()->GetXmax());
 h3->SetDirectory(0);

 for (int i=0;i<h1->GetNbinsX(); i++) 
 {
  double value = a*h1->GetBinContent(i)+b*h2->GetBinContent(i);
  double error = sqrt((a*a*h1->GetBinError(i)*h1->GetBinError(i))+(b*b*h2->GetBinError(i)*h2->GetBinError(i)));
  h3->SetBinContent(i,value);
  h3->SetBinError(i,error);
 }

 return h3;
}


// this program is to calculate and return the error of asymmetry histogram.
// Here, hAsym = (h1 - h2)/(h1 + h2).
 TH1F* GetAsymHist(TH1F* h1, TH1F *h2)
 {
  double error;
  TH1F *hasym = new TH1F("hasym", "hasym", h1->GetNbinsX(), h1->GetXaxis()->GetXmin(), h1->GetXaxis()->GetXmax());
  hasym->SetDirectory(0);
  for (int i=1;i<h1->GetNbinsX()+1; i++) 
  {
   double value1 = h1->GetBinContent(i);
   double value2 = h2->GetBinContent(i);
   double error1 = h1->GetBinError(i);
   double error2 = h2->GetBinError(i);
   if((value1+value2) == 0)
   {
    hasym->SetBinContent(i, 0);
    hasym->SetBinError(i, 0);
   }
   else
   {
    error  = 2*sqrt( pow(value2, 2)*pow(error1, 2) + pow(value1, 2)*pow(error2, 2) )/pow(value1+value2, 2);

    hasym->SetBinContent(i, (value1 - value2)/(value1 + value2));
    hasym->SetBinError(i, error);
   }
  }

  return hasym;
 }

// this program is to calculate and return the error of asymmetry histogram.
// Here, hAsym = (h1 - h2)/(h3 + h4).
 TH1F* GetAsymHist(TH1F* h1, TH1F *h2, TH1F *h3, TH1F *h4)
 {
  double error;
  TH1F *hasym = new TH1F("hasym", "hasym", h1->GetNbinsX(), h1->GetXaxis()->GetXmin(), h1->GetXaxis()->GetXmax()
);
  hasym->SetDirectory(0);
  for (int i=1;i<h1->GetNbinsX()+1; i++)
  {
   double value1 = h1->GetBinContent(i);
   double value2 = h2->GetBinContent(i);
   double value3 = h3->GetBinContent(i);
   double value4 = h4->GetBinContent(i);

   double error1 = h1->GetBinError(i);
   double error2 = h2->GetBinError(i);
   double error3 = h3->GetBinError(i);
   double error4 = h4->GetBinError(i);

   if((value3+value4) == 0)
   {
    hasym->SetBinContent(i, 0);
    hasym->SetBinError(i, 0);
   }
   else
   {
    error = sqrt( (value1-value2)*(value1-value2)*(error3*error3+error4*error4) + (value3+value4)*(value3+value4)*(error1*error1+error2*error2) )/( (value3+value4)*(value3+value4) );
//    error  = 2*sqrt( pow(value2, 2)*pow(error1, 2) + pow(value1, 2)*pow(error2, 2) )/pow(value1+value2, 2);

    hasym->SetBinContent(i, (value1 - value2)/(value3 + value4));
    hasym->SetBinError(i, error);
   }
  }

  return hasym;
 }