{
  //rebin to rb
  int rb=10;
  //scale signal to sc
  double sc=1.;
  TFile f150("ZHH-BP1-85k_1.5TeV.root"); //signal
  TFile f200("ZHH-BP2-85k_1.5TeV.root");
  TFile f250("ZHH-BP3-85k_1.5TeV.root");
  TFile f300("ZHH-BP4-85k_1.5TeV.root");
  TFile fbtt("tT-ZHH-85k_1.5TeV.root"); //backgrounds
  //TFile fbww("WW-ZHH-85k_1.5TeV.root");
  TFile fbzz("ZZ-ZHH-85k_1.5TeV.root");
  TFile fbzbb("ZbB-ZHH-85k_1.5TeV.root");
  TFile fbzzz("ZZZ-ZHH-85k_1.5TeV.root");


  TH1F *h150=(TH1F*)f150.Get("Higgsmass1");
  TH1F *h200=(TH1F*)f200.Get("Higgsmass1");
  TH1F *h250=(TH1F*)f250.Get("Higgsmass1");
  TH1F *h300=(TH1F*)f300.Get("Higgsmass1");
  TH1F *hbtt=(TH1F*)fbtt.Get("Higgsmass1");
 // TH1F *hbww=(TH1F*)fbww.Get("Higgsmass1");
  TH1F *hbzz=(TH1F*)fbzz.Get("Higgsmass1");
  TH1F *hbzbb=(TH1F*)fbzbb.Get("Higgsmass1");
  TH1F *hbzzz=(TH1F*)fbzzz.Get("Higgsmass1");
 

   // double sigmaxBR[9]={0.062,0.050,0.040,0.031,80.440,1511.9,84.287,32.818,0.739};
	   
   /*  for(int i=0;i<9;i++) {
     sigmaxBR[i]*=680.;// * 1000 fb * 0.49 (Z->jj)
       cout<<"sigma BR ====  "<<sigmaxBR[i]<<endl;
      }      
    sigmaxBR[6]*=0.68;//ttbar has two tops so BR(W->jj)^2 should be multiplied
       cout<<"sigma BR ====  "<<sigmaxBR[6]<<endl;
    */

   double sigmaxBR[8]={0.} ;
   double sigma[8]={0.062,0.050,0.040,0.031,80.440,84.287,32.818,0.739}; //just values of crossection first for signal and then for background
   double Bratio[8]={0.0007,0.0002,0.0001,0.00003,0.444,0.294,0.294,1.626};

  for(int i=0;i<8;i++) {
     sigmaxBR[i] =  sigma[i]*Bratio[i] ;          
       cout<<"sigma BR ====  "<<sigmaxBR[i]<<endl;
    }
     
  //double eff[10]={36905./1000000,22872./1000000,14866./1000000,10095./1000000,14206./1000000,5154/1000000,4491./1000000,375/2000000,4595/3000000,11800/3000000};
  double eff[8]={74.3,81.0,84.9,87.3,4.919,0.005,0.433,0.167};
  double lumi=1000.;

  h150->Rebin(rb);
  h200->Rebin(rb);
  h250->Rebin(rb);
  h300->Rebin(rb);
  hbtt->Rebin(rb);
 // hbww->Rebin(rb);
  hbzz->Rebin(rb);
  hbzbb->Rebin(rb);
  hbzzz->Rebin(rb);

  double nevents[8]={0.};
  for(int i=0;i<8;i++){
    nevents[i]=sigmaxBR[i]*eff[i]*lumi;
    if(i<4)nevents[i]*=sc;
  }


  h150->Scale(nevents[0]/h150->Integral());
  h200->Scale(nevents[1]/h200->Integral());
  h250->Scale(nevents[2]/h250->Integral());
  h300->Scale(nevents[3]/h300->Integral());
  hbtt->Scale(nevents[4]/hbtt->Integral());
  //hbww->Scale(nevents[5]/hbww->Integral());
  hbzz->Scale(nevents[5]/hbzz->Integral());
  hbzbb->Scale(nevents[6]/hbzbb->Integral());
  hbzzz->Scale(nevents[7]/hbzzz->Integral());
  


  TH1F *hb=(TH1F*)hbtt->Clone();
  hb->Add(hbtt);
 // hb->Add(hbww);
  hb->Add(hbzz);
  hb->Add(hbzbb);
  hb->Add(hbzzz);  

  TH1F *hs150=(TH1F*)hbtt->Clone();
  hs150->Add(hbtt);
  //hs150->Add(hbww);
  hs150->Add(hbzz);
  hs150->Add(hbzbb);
  hs150->Add(hbzzz);
  hs150->Add(h150);

  TH1F *hs200=(TH1F*)hbtt->Clone();
  hs200->Add(hbtt);
  //h200->Add(hbww);
  hs200->Add(hbzz);
  hs200->Add(hbzbb);
  hs200->Add(hbzzz); 
  hs200->Add(h200); 

  TH1F *hs250=(TH1F*)hbtt->Clone();
  hs250->Add(hbtt);
  //hs250->Add(hbww);
  hs250->Add(hbzz);
  hs250->Add(hbzbb);
  hs250->Add(hbzzz); 
  hs250->Add(h250);

  TH1F *hs300=(TH1F*)hbtt->Clone();
  hs250->Add(hbtt);
  //hs250->Add(hbww);
  hs250->Add(hbzz);
  hs250->Add(hbzbb);
  hs250->Add(hbzzz); 
  hs250->Add(h300);

 

  hb->SetFillColor(2);
  hs150->SetFillColor(3);
  hs200->SetFillColor(4);
  hs250->SetFillColor(5);
  hs300->SetFillColor(7);

  hs150->Draw("hist same");
  hs150->GetXaxis()->SetTitle("m_{H} [GeV]");
  hs150->GetXaxis()->CenterTitle(true);
  hs150->GetYaxis()->SetTitle("Number of Events");
  hs150->GetYaxis()->CenterTitle(true);
  hs200->Draw("hist same");
  hs250->Draw("hist same");
  hs300->Draw("hist same");
  hb->Draw("hist same");
  hbtt->Draw("hist same");
 // hbww->Draw("hist same");
  hbzz->Draw("hist same");
  hbzbb->Draw("hist same");
  hbzzz->Draw("hist same");

  TLegend *lg=new TLegend(0.3,0.3,0.8,0.8);
  lg->AddEntry(h150,"m_{H}=150 GeV","F");
  lg->AddEntry(h200,"m_{H}=200 GeV","F");
  lg->AddEntry(h250,"m_{H}=250 GeV","F");
  lg->AddEntry(h300,"m_{H}=300 GeV","F");
  lg->AddEntry(hb,"Total background","F");
  lg->Draw();

int m=0;
  double TotalEffS[5]={0.};
  int nbins=h250->GetNbinsX();
  maxsig=0;
  int bin1=0;
  int bin2=0;
  double S=0.;
  double B=0.;
  for(int i=0;i<nbins;i++){
    for(int j=i+1;j<nbins;j++){
      double s2b=(hb->Integral(i,j)==0 ? 0 : double(h150->Integral(i,j))/hb->Integral(i,j));//S/B
      double significance=(hb->Integral(i,j)==0 ? 0 : double(h150->Integral(i,j))/sqrt(hb->Integral(i,j)));//S/sqrt(B)
      if(significance>maxsig){maxsig=significance;bin1=(i*rb);bin2=(j*rb);S=h150->Integral(i,j);B=hb->Integral(i,j);TotalEffS[m]=S*eff[m]/h150->Integral();}
    }
  }
  cout<<"============================"<<endl;
  cout<<"Max significance with bins "<<bin1<<"-"<<bin2<<" is: "<<maxsig<<endl;
  cout<<"S and B : "<<S<<", "<<B<<endl;
  cout<<"S/B : "<<S/B<<endl;
  cout<<"S/sqrt(B) : "<<S/sqrt(B)<<endl;
  cout<<"Total signal eff.: "<<TotalEffS[0]<<endl;
  

  m++;
  maxsig=bin1=bin2=S=B=0;
  for(int i=0;i<nbins;i++){
    for(int j=i+1;j<nbins;j++){
      double s2b=(hb->Integral(i,j)==0 ? 0 : double(h200->Integral(i,j))/hb->Integral(i,j));
      double significance=(hb->Integral(i,j)==0 ? 0 : double(h200->Integral(i,j))/sqrt(hb->Integral(i,j)));
      if(significance>maxsig){maxsig=significance;bin1=(i*rb);bin2=(j*rb);S=h200->Integral(i,j);B=hb->Integral(i,j);TotalEffS[m]=S*eff[m]/h200->Integral();}
      //      cout<<"Mass Window, S, B, S/B, S/sqrt(b) : "<<i*rb<<"-"<<j*rb<<", "<<h500->Integral(i,j)<<", "<<hb->Integral(i,j)<<", "
      //	  <<s2b<<", "<<significance<<endl; 
    }
  }
  cout<<"============================"<<endl;
  cout<<"Max significance with bins "<<bin1<<"-"<<bin2<<" is: "<<maxsig<<endl;
  cout<<"S and B are: "<<S<<", "<<B<<endl;
  cout<<"S/B : "<<S/B<<endl;
  cout<<"S/sqrt(B) : "<<S/sqrt(B)<<endl;
  cout<<"Total signal eff.: "<<TotalEffS[1]<<endl;

  m++;
  maxsig=bin1=bin2=S=B=0;
  for(int i=0;i<nbins;i++){
    for(int j=i+1;j<nbins;j++){
      double s2b=(hb->Integral(i,j)==0 ? 0 : double(h250->Integral(i,j))/hb->Integral(i,j));
      double significance=(hb->Integral(i,j)==0 ? 0 : double(h250->Integral(i,j))/sqrt(hb->Integral(i,j)));
      if(significance>maxsig){maxsig=significance;bin1=(i*rb);bin2=(j*rb);S=h250->Integral(i,j);B=hb->Integral(i,j);TotalEffS[m]=S*eff[m]/h250->Integral();}
      //      cout<<"Mass Window, S, B, S/B, S/sqrt(b) : "<<i*rb<<"-"<<j*rb<<", "<<h600->Integral(i,j)<<", "<<hb->Integral(i,j)<<", "
      //	  <<s2b<<", "<<significance<<endl; 
    }
  }
  cout<<"============================"<<endl;
  cout<<"Max significance with bins "<<bin1<<"-"<<bin2<<" is: "<<maxsig<<endl;
  cout<<"S and B are: "<<S<<", "<<B<<endl;
  cout<<"S/B : "<<S/B<<endl;
  cout<<"S/sqrt(B) : "<<S/sqrt(B)<<endl;
  cout<<"Total signal eff.: "<<TotalEffS[2]<<endl;

  m++;
  maxsig=bin1=bin2=S=B=0;
  for(int i=0;i<nbins;i++){
    for(int j=i+1;j<nbins;j++){
      double s2b=(hb->Integral(i,j)==0 ? 0 : double(h300->Integral(i,j))/hb->Integral(i,j));
      double significance=(hb->Integral(i,j)==0 ? 0 : double(h300->Integral(i,j))/sqrt(hb->Integral(i,j)));
      if(significance>maxsig){maxsig=significance;bin1=(i*rb);bin2=(j*rb);S=h300->Integral(i,j);B=hb->Integral(i,j);TotalEffS[m]=S*eff[m]/h300->Integral();}
      //      cout<<"Mass Window, S, B, S/B, S/sqrt(b) : "<<i*rb<<"-"<<j*rb<<", "<<h700->Integral(i,j)<<", "<<hb->Integral(i,j)<<", "
      //	  <<s2b<<", "<<significance<<endl; 
    }
  }
  cout<<"============================"<<endl;
  cout<<"Max significance with bins "<<bin1<<"-"<<bin2<<" is: "<<maxsig<<endl;
  cout<<"S and B are: "<<S<<", "<<B<<endl;
  cout<<"S/B : "<<S/B<<endl;
  cout<<"S/sqrt(B) : "<<S/sqrt(B)<<endl;
  cout<<"Total signal eff.: "<<TotalEffS[3]<<endl;

}