C++ TH1F::GetBinLowEdge方法代码示例

  TH1F * getpdgs(std::string const & dr, std::string const & we) {
    
    TCanvas * canvas_temp = new TCanvas("temp");
    tree->Draw((dr+">>h").c_str(), we.c_str());
    TH1F * hist = (TH1F*)gDirectory->Get("h");
    int const min = hist->GetBinLowEdge(hist->FindFirstBinAbove(0));
    int const max = hist->GetBinLowEdge(hist->FindLastBinAbove(0)) +
      hist->GetBinWidth(hist->FindLastBinAbove(0)) + 1;
    int const binno = max - min;
    delete canvas_temp;
    delete hist;

    canvas_temp = new TCanvas("temp");
    tree->Draw(Form((dr+">>h(%d,%d,%d)").c_str(), binno, min, max), we.c_str());
    ofile->cd();
    delete canvas_temp; 
    TH1F * histb = (TH1F*)gDirectory->Get("h");

    if(histb->GetEntries() == 0) {
      delete histb;
      return nullptr;
    }    

    return histb;
    
  }

void performClosure(RooRealVar *mass, RooAbsPdf *pdf, RooDataSet *data, string closurename, double wmin=110., double wmax=130., double slow=110., double shigh=130., double step=0.002) {
  
  // plot to perform closure test
  cout << "Performing closure test..." << endl; 
  double nbins = (wmax-wmin)/step;
  TH1F *h = new TH1F("h","h",int(floor(nbins+0.5)),wmin,wmax);
  if (data){
    pdf->fillHistogram(h,RooArgList(*mass),data->sumEntries());
    h->Scale(2*h->GetNbinsX()/double(binning_));
  }
  else {
    pdf->fillHistogram(h,RooArgList(*mass));
  }
  int binLow = h->FindBin(slow);
  int binHigh = h->FindBin(shigh)-1;
  TH1F *copy = new TH1F("copy","c",binHigh-binLow,h->GetBinLowEdge(binLow),h->GetBinLowEdge(binHigh+1));
  for (int b=0; b<copy->GetNbinsX(); b++) copy->SetBinContent(b+1,h->GetBinContent(b+1+binLow));
  double areaCov = 100*h->Integral(binLow,binHigh)/h->Integral();
 
  // style
  h->SetLineColor(kBlue);
  h->SetLineWidth(3);
  h->SetLineStyle(7);
  copy->SetLineWidth(3);
  copy->SetFillColor(kGray);
  
  TCanvas *c = new TCanvas();
  if (data){
    RooPlot *plot = mass->frame(Bins(binning_),Range("higgsRange"));
    plot->addTH1(h,"hist");
    plot->addTH1(copy,"same f");
    if (data) data->plotOn(plot);
    pdf->plotOn(plot,Normalization(h->Integral(),RooAbsReal::NumEvent),NormRange("higgsRange"),Range("higgsRange"),LineWidth(1),LineColor(kRed),LineStyle(kDashed));
    plot->Draw();
    c->Print(closurename.c_str());
  }
  else {
    RooPlot *plot = mass->frame(Bins(binning_),Range("higgsRange"));
    h->Scale(plot->getFitRangeBinW()/h->GetBinWidth(1));
    copy->Scale(plot->getFitRangeBinW()/h->GetBinWidth(1));
    pdf->plotOn(plot,LineColor(kRed),LineWidth(3));
    plot->Draw();
    h->Draw("hist same");
    copy->Draw("same f");
    c->Print(closurename.c_str());
  }
  cout << "IntH: [" << h->GetBinLowEdge(binLow) << "-" << h->GetBinLowEdge(binHigh+1) << "] Area = " << areaCov << endl;
  delete c;
  delete copy;
  delete h;
}

int main(int argc, char* argv[]){

  ic::Plot::SetTdrStyle();
  TH1F data = GetFromTFile<TH1F>("datacard_m_vis_inclusive_et_2011.root", "/eleTau_inclusive", "data_obs");
  TH1F zee = GetFromTFile<TH1F>("datacard_m_vis_inclusive_et_2011.root", "/eleTau_inclusive", "ZL");

  double data_rate = data.Integral();
  double zee_init_rate = zee.Integral();

  RooRealVar x("x","mvis",data.GetBinLowEdge(1),data.GetBinLowEdge(data.GetNbinsX()+1));
  RooRealVar c1("c1","c1",-10,10);
  RooRealVar c2("c2","c2",-1,1);
  RooRealVar c3("c3","c3",-1,1);
  RooGenericPdf bkg("bkg","exp(c1*x+c2*x*x+c3*x*x*x)",RooArgSet(x,c1,c2,c3));


  RooDataHist zee_hist("zee_hist","zee_hist",RooArgSet(x),&zee);
  RooHistPdf zee_pdf("zee_pdf","zee_pdf",RooArgSet(x),zee_hist);


  ///Fit Data
  TCanvas C("canvas", "canvas", 800,800);

  RooRealVar yield("yield","yield",0.01,.9);
  RooAddPdf full_pdf("full_pdf","full_pdf",RooArgList(zee_pdf,bkg),RooArgList(yield));
  RooDataHist data_hist("data_hist","data_hist",RooArgSet(x),&data);
  RooChi2Var chi("chi","chi",full_pdf,data_hist,RooFit::DataError(RooAbsData::SumW2));
  RooMinuit minuit(chi);
  minuit.migrad();
  RooPlot* plot=x.frame();
  data_hist.plotOn(plot);
  full_pdf.plotOn(plot);
  full_pdf.plotOn(plot,RooFit::Components(zee_pdf),RooFit::LineColor(2));
  C.Clear();
  plot->SetTitle("Title");
  plot->GetYaxis()->SetTitle("");
  plot->GetXaxis()->SetTitle("m(e#tau)");
  plot->Draw();
  C.Print("plotZEE.pdf");

  double fit_frac = yield.getVal();
  double fit_frac_err = yield.getError();
  double zee_fit_rate = data_rate * fit_frac;
  double err = (data_rate * fit_frac_err) / zee_init_rate;
  std::cout << "Scale factor: " << zee_fit_rate / zee_init_rate << " +/- " << err << std::endl;

  return 0;
}

void PoissonianMagic()
{
  TFile file_PUnum("/nfs/dust/cms/user/rathjd/VBF-LS-tau/PU/DataPUFile_22Jan2013ReReco_Run2012.root", "read");
  TFile file_PUden("/nfs/dust/cms/user/rathjd/VBF-LS-tau/PU/S10MC_PUFile.root", "read");
  
  TH1F *data = (TH1F*)file_PUnum.Get("analyzeHiMassTau/NVertices_0");
  data->Scale(1/data->Integral(0,-1));
  TH1F *MC   = (TH1F*)file_PUden.Get("analyzeHiMassTau/NVertices_0");
  MC->Scale(1/MC->Integral(0,-1));
  
  //define empty histograms for the smeared versions
  TH1F *num = (TH1F*)data->Clone("ratio");
  for(int i=0; i<num->GetNbinsX(); i++) num->SetBinContent(i+1,0);
  TH1F *den = (TH1F*)data->Clone("MC");
  for(int i=0; i<den->GetNbinsX(); i++) den->SetBinContent(i+1,0);
  
  //generate the poissonian distributions
  for(unsigned int i=0; i<data->GetNbinsX(); i++){
    TF1 *Pd = new TF1("Pd", "TMath::PoissonI(x,[0])",0,100);
    Pd->SetParameter(0,data->GetBinLowEdge(i+1));
    std::cout<<data->GetBinLowEdge(i+1)<<std::endl;
    TF1 *Pm = new TF1("Pm", "TMath::PoissonI(x,[0])",0,100);
    Pm->SetParameter(0,MC->GetBinLowEdge(i+1));
    //add up the poissonians
    for(unsigned int j=0; j<data->GetNbinsX(); j++){
      num->SetBinContent(j+1,num->GetBinContent(j+1)+Pd->Eval(j)/Pd->Integral(0,100)*data->GetBinContent(i+1));
      den->SetBinContent(j+1,den->GetBinContent(j+1)+Pm->Eval(j)/Pm->Integral(0,100)*MC->GetBinContent(i+1));
    }
  }
  
  //make the ratio and save the result
  num->Divide(den);
  
  TFile *f=new TFile("PUreweightHistogram.root","RECREATE");
  num->Write();
  den->Write();
  
  f->Close();
  file_PUnum.Close();
  file_PUden.Close();
}

//.........这里部分代码省略.........
  datapredictiont->Draw();
  functions[1]->Draw("same");
  TText *title1 = write_title("Top Background Prediction (JZB<0, with osof subtr)");
  title1->Draw();
  
  c4->cd(2);
  c4->cd(2)->SetLogy(1);
  TH1F *observedt = allsamples.Draw("observedt", datajzb,  nbins,low,hi, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,  luminosity);
  observedt->Draw();
  datapredictiont->Draw("histo,same");
  functions[1]->Draw("same");
  TText *title2 = write_title("Observed and predicted background");
  title2->Draw();
  
  c4->cd(3);
  c4->cd(3)->SetLogy(1);
//  TH1F *ratio = (TH1F*)observedt->Clone();

  TH1F *analytical_background_prediction= new TH1F("analytical_background_prediction","",nbins,low,hi);
  for(int i=0;i<=nbins;i++) {
    analytical_background_prediction->SetBinContent(i+1,functions[1]->Eval(((hi-low)/((float)nbins))*(i+0.5)));
    analytical_background_prediction->SetBinError(i+1,TMath::Sqrt(functions[1]->Eval(((hi-low)/((float)nbins))*(i+0.5))));
  }
  analytical_background_prediction->GetYaxis()->SetTitle("JZB [GeV]");
  analytical_background_prediction->GetYaxis()->CenterTitle();
  TH1F *analyticaldrawonly = (TH1F*)analytical_background_prediction->Clone();
  analytical_background_prediction->SetFillColor(TColor::GetColor("#3399FF"));
  analytical_background_prediction->SetMarkerSize(0);
  analytical_background_prediction->Draw("e5");
  analyticaldrawonly->Draw("histo,same");
  functions[1]->Draw("same");
  TText *title3 = write_title("Analytical bg pred histo");
  title3->Draw();
  
  c4->cd(4);
//  c4->cd(4)->SetLogy(1);
  vector<float> ratio_binning;
  ratio_binning.push_back(0);
  ratio_binning.push_back(5);
  ratio_binning.push_back(10);
  ratio_binning.push_back(20);
  ratio_binning.push_back(50);
//  ratio_binning.push_back(60);
/*
  ratio_binning.push_back(51);
  ratio_binning.push_back(52);
  ratio_binning.push_back(53);
  ratio_binning.push_back(54);
  ratio_binning.push_back(55);
  ratio_binning.push_back(56);
  ratio_binning.push_back(57);
  ratio_binning.push_back(58);
  ratio_binning.push_back(59);
  ratio_binning.push_back(60);
//  ratio_binning.push_back(70);*/
//  ratio_binning.push_back(80);
//  ratio_binning.push_back(90);
ratio_binning.push_back(80);
//  ratio_binning.push_back(110);
  ratio_binning.push_back(500);
  
  TH1F *observedtb = allsamples.Draw("observedtb", datajzb,  ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,  luminosity);
  TH1F *datapredictiontb = allsamples.Draw("datapredictiontb",    "-"+datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets,data,  luminosity);
  TH1F *datapredictiontbo = allsamples.Draw("datapredictiontbo",    "-"+datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSOF&&cutnJets,data,  luminosity);
  datapredictiontb->Add(datapredictiontbo,-1);
  TH1F *analytical_background_predictionb = allsamples.Draw("analytical_background_predictionb",datajzb, ratio_binning, "JZB [GeV]", "events", cutmass&&cutOSSF&&cutnJets&&"mll<2",data,  luminosity);
  for(int i=0;i<=(int)ratio_binning.size();i++) {
    analytical_background_predictionb->SetBinContent(i+1,functions[1]->Eval(analytical_background_predictionb->GetBinCenter(i)));
    analytical_background_predictionb->SetBinError(i+1,TMath::Sqrt(functions[1]->Eval(analytical_background_predictionb->GetBinCenter(i))));
  }
  
  TH1F *ratio = (TH1F*) observedtb->Clone();
  ratio->Divide(datapredictiontb);
  
  for (int i=0;i<=(int)ratio_binning.size();i++) {
    dout << observedtb->GetBinLowEdge(i+1) << ";"<<observedtb->GetBinContent(i+1) << ";" << datapredictiontb->GetBinContent(i+1) << " --> " << ratio->GetBinContent(i+1) << "+/-" << ratio->GetBinError(i+1) << endl;
  }
  
//  ratio->Divide(datapredictiontb);
//  ratio->Divide(analytical_background_predictionb);
//  TGraphAsymmErrors *JZBratio= histRatio(observedtb,analytical_background_predictionb,data,ratio_binning);
//  ratio->Divide(analytical_background_prediction);
//  ratio->Divide(datapredictiont);
//  ratio->GetYaxis()->SetTitle("obs/pred");
//  JZBratio->Draw("AP");
  ratio->GetYaxis()->SetRangeUser(0,10);
  ratio->Draw();
//analytical_background_predictionb->Draw();
//  JZBratio->SetTitle("");
  TText *title4 = write_title("Ratio of observed to predicted");
  title4->Draw();
  
//  CompleteSave(c4,"test/ttbar_discovery_dataprediction___analytical_function");
  CompleteSave(c4,"test/ttbar_discovery_dataprediction__analytical__new_binning_one_huge_bin_from_80");

  
  
  
  
}

void AnalysisBase::findDeadRegions(){

  cout << "=================================================" << endl;
  cout << "findDeadRegions(): Looking for dead strip regions " << endl;
  cout << "=================================================" << endl;

  TH1F* hf = new TH1F("hf","X position of matched cluster",200,-10.0,10.0);
  TH1F* hnf = new TH1F("hnf","X position of matched cluster",200,-10.0,10.0);
  
  double nomStrip = 0, detStrip = 0;
  Long64_t nbytes = 0, nb = 0;
  Int_t nentries = fChain->GetEntriesFast();
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
      
    if(n_tp3_tracks > 1) continue;

    // Loop over TPIX tracks in event
    for(int k=0; k<n_tp3_tracks; k++){
      double x_trk = vec_trk_tx->at(k)*z_DUT+vec_trk_x->at(k);
      double y_trk = vec_trk_ty->at(k)*z_DUT+vec_trk_y->at(k);

      transformTrackToDUTFrame(k, x_trk, y_trk, nomStrip, detStrip );

      bool goodTrack = false;
      bool inFiducial = false;
      if(x_trk>xMin && x_trk<xMax && y_trk>yMin && y_trk<yMax) inFiducial = true;          
      inFiducial = inFiducial && (x_trk<xLeftHole || x_trk>xRightHole);
      double tx = 1000*vec_trk_tx->at(k);
      double ty = 1000*vec_trk_ty->at(k);
      if(tx>txMin && tx<txMax && ty>tyMin && ty<tyMax) goodTrack = true;        
      bool goodTime =  (clustersTDC >= tdcLo && clustersTDC <= tdcHi);

      bool foundHit = false;
      // Loop over clusters
      for(int j=0; j<min(clusterNumberPerEvent,10); j++){
        if(clustersPosition[j] < 0.1) continue;
        if(polarity*clustersCharge[j] < kClusterChargeMin) continue;
        //bool goodHit = (clustersPosition[j]>iLo || clustersPosition[j]<iHi);
        double x_dut = getDUTHitPosition(j);

        double dx = x_dut - x_trk;

        if(goodTrack && inFiducial && goodTime && fabs(dx)<dxWin) {
          foundHit = true;
        }   
      }
      if(inFiducial && goodTrack && goodTime) {
        hf->Fill(x_trk);    
        if(!foundHit) hnf->Fill(x_trk);          
      }
        
    }
  }

  TH1F *hineff = (TH1F*)hnf->Clone("hineff");
  hineff->Divide(hf);

  // Only for debugging
  /*
    TCanvas *cc = new TCanvas("cc","Hits",800,800);
    cc->Divide(2,2);
   
    cc->cd(1);
    hnf->Draw();
    cc->cd(2);
    hf->Draw();
    cc->cd(3);
    hineff->Draw();
  */


  for(int k = 0; k<hineff->GetNbinsX()-1; k++){
    if(nDeadRegion > 20){
      cout << "WARNING: Exceeded 20 dead regions, exiting from loop" << endl;       
      break;
    }     
    double v1 = hineff->GetBinContent(k);
    if(v1 > k_DeadStrip){
      deadRegionLo[nDeadRegion] = hf->GetBinLowEdge(k);
      for(int j = k; j<hineff->GetNbinsX(); j++){
        double v2 = hineff->GetBinContent(j);
        if(v2 < k_DeadStrip){
          deadRegionHi[nDeadRegion] = hf->GetBinLowEdge(j);
          cout << "====> Found inefficiency strip region from x :    " 
               << deadRegionLo[nDeadRegion] << " <===> " << deadRegionHi[nDeadRegion] << " mm" << endl;
          k = j + 1;
          nDeadRegion++;
          break;
        }    
      }
    }
  }
   
  delete hf;
  delete hnf;
  delete hineff;
   
//.........这里部分代码省略.........

void plotHitEnergy(string inputDir, float radius=0.4){
  setNCUStyle(true);

  TH1F* hecal;
  TH1F* hhcal;


  TString energy=gSystem->GetFromPipe(Form("file=%s; test=${file##*/}; test2=${test%%mumu*}; echo \"${test2##*tev}\"",inputDir.data()));
  cout << "energy = " << energy.Data() << endl;
  
  string inputFile = inputDir + "/radius" + Form("%0.1f",radius)+ "_rawhit.root";
  cout << "opening " << inputFile.data() << endl;


  
  TFile *f = TFile::Open(inputFile.data());
  hecal = (TH1F*)f->FindObjectAny("hecalhit_energy");
  hecal->SetLineWidth(3);
  hecal->SetFillStyle(1001);
  hecal->SetFillColor(4);
  hecal->SetLineColor(4);
  hecal->SetXTitle("EM_BARREL hit energy [GeV]");
  hecal->SetYTitle(Form("Number of hits per %d GeV",(int)hecal->GetBinWidth(1)));
  hecal->SetTitleOffset(1.2,"X");
  hecal->SetTitleOffset(1.4,"Y");
  hhcal = (TH1F*)f->FindObjectAny("hhcalhit_energy");
  hhcal->SetLineWidth(3);
  hhcal->SetFillStyle(1001);
  hhcal->SetFillColor(2);
  hhcal->SetLineColor(2);
  hhcal->SetXTitle("HAD_BARREL hit energy [GeV]");
  hhcal->SetYTitle(Form("Number of hits per %d GeV",(int)hhcal->GetBinWidth(1)));
  hhcal->SetTitleOffset(1.2,"X");
  hhcal->SetTitleOffset(1.4,"Y");


  TLegend* leg = new TLegend(0.444,0.690,0.990,0.903);
  leg->SetFillColor(0);
  leg->SetFillStyle(0);

  TCanvas* c1 = new TCanvas("c1","",500,500);
  int lastbin=hecal->FindLastBinAbove(0)+50;
  float xmax=hecal->GetBinLowEdge(lastbin);
  hecal->GetXaxis()->SetRangeUser(0,xmax);
  hecal->Draw("hist");
  c1->SetLogy(1);
  
  leg->SetHeader("rfull012");
  leg->AddEntry(hecal,Form("%s-TeV Z'#rightarrow qq",energy.Data()),"f");
  leg->Draw("same");

  c1->Print(Form("rfull012/rfull012_EM_BARREL_hit_energy_%sTeVZp.pdf",energy.Data()));
  c1->Print(Form("rfull012/rfull012_EM_BARREL_hit_energy_%sTeVZp.eps",energy.Data()));
  leg->Clear();


  lastbin=hhcal->FindLastBinAbove(0)+50;
  xmax=hhcal->GetBinLowEdge(lastbin);
  hhcal->GetXaxis()->SetRangeUser(0,xmax);
  hhcal->Draw("hist");
  c1->SetLogy(1);
  
  leg->SetHeader("rfull012");
  leg->AddEntry(hhcal,Form("%s-TeV Z'#rightarrow qq",energy.Data()),"f");
  leg->Draw("same");

  c1->Print(Form("rfull012/rfull012_HAD_BARREL_hit_energy_%sTeVZp.pdf",energy.Data()));
  c1->Print(Form("rfull012/rfull012_HAD_BARREL_hit_energy_%sTeVZp.eps",energy.Data()));


}

void Fit(TString SIGNAL,TString HISTO,double scaleSGN)
{
  gROOT->ForceStyle();
  TFile *fDat = TFile::Open("Histo_flatTree_data_tmva"+SIGNAL+".root");
  TFile *fBkg = TFile::Open("Histo_flatTree_qcd_weights_tmva"+SIGNAL+".root");
  TFile *fSgn = TFile::Open("Histo_flatTree_"+SIGNAL+"_weights_tmva"+SIGNAL+".root");
  
  TH1 *hDat = (TH1*)fDat->Get(HISTO);
  TH1 *hBkgRaw = (TH1*)fBkg->Get(HISTO);
  TH1 *hSgn = (TH1*)fSgn->Get(HISTO);
  TH1 *hDat_JESlo = (TH1*)fDat->Get(HISTO+"_JESlo");
  TH1 *hBkgRaw_JESlo = (TH1*)fBkg->Get(HISTO+"_JESlo");
  TH1 *hSgn_JESlo = (TH1*)fSgn->Get(HISTO+"_JESlo");
  TH1 *hDat_JESup = (TH1*)fDat->Get(HISTO+"_JESup");
  TH1 *hBkgRaw_JESup = (TH1*)fBkg->Get(HISTO+"_JESup");
  TH1 *hSgn_JESup = (TH1*)fSgn->Get(HISTO+"_JESup");
  
  TH1F *hBkg = (TH1F*)hBkgRaw->Clone("Bkg");
  TH1F *hBkg_JESlo = (TH1F*)hBkgRaw_JESlo->Clone("Bkg_JESlo");
  TH1F *hBkg_JESup = (TH1F*)hBkgRaw_JESup->Clone("Bkg_JESup");
  
  hBkg->Smooth(2);
  hBkg_JESlo->Smooth(2);
  hBkg_JESup->Smooth(2);
  hSgn->Smooth(2);
  hSgn_JESlo->Smooth(2);
  hSgn_JESup->Smooth(2);
  
  double lumi = 4967;
  hBkg->Scale(lumi);
  hBkg_JESlo->Scale(lumi);
  hBkg_JESup->Scale(lumi);
  double k_factor = hDat->Integral()/hBkg->Integral();
  double k_factor_JESlo = hDat->Integral()/hBkg_JESlo->Integral();
  double k_factor_JESup = hDat->Integral()/hBkg_JESup->Integral();
  hBkg->Scale(k_factor);
  cout<<"Signal entries = "<<hSgn->GetEntries()<<endl;
  hSgn->Scale(lumi/scaleSGN);
  hBkg_JESlo->Scale(k_factor_JESlo);
  hSgn_JESlo->Scale(lumi/scaleSGN);
  hBkg_JESup->Scale(k_factor_JESup);
  hSgn_JESup->Scale(lumi/scaleSGN);
  hSgn_JESlo->Scale(hSgn->Integral()/hSgn_JESlo->Integral());
  hSgn_JESup->Scale(hSgn->Integral()/hSgn_JESup->Integral());
  
  TH1 *hBkg_STATlo = (TH1*)hBkg->Clone(HISTO+"_STATlo");
  TH1 *hSgn_STATlo = (TH1*)hSgn->Clone(HISTO+"_STATlo");
  TH1 *hBkg_STATup = (TH1*)hBkg->Clone(HISTO+"_STATup");
  TH1 *hSgn_STATup = (TH1*)hSgn->Clone(HISTO+"_STATup");
  
  float y1,e1;
  for(int i=0;i<hBkg->GetNbinsX();i++) {
    y1 = hBkg->GetBinContent(i+1);
    e1 = hBkg->GetBinError(i+1);
    hBkg_STATlo->SetBinContent(i+1,y1-e1);
    hBkg_STATup->SetBinContent(i+1,y1+e1);
    y1 = hSgn->GetBinContent(i+1);
    e1 = hSgn->GetBinError(i+1);
    hSgn_STATlo->SetBinContent(i+1,y1-e1);
    hSgn_STATup->SetBinContent(i+1,y1+e1);
  }
  hBkg_STATlo->Scale(hBkg->Integral()/hBkg_STATlo->Integral());
  hBkg_STATup->Scale(hBkg->Integral()/hBkg_STATup->Integral());
  hSgn_STATlo->Scale(hSgn->Integral()/hSgn_STATlo->Integral());
  hSgn_STATup->Scale(hSgn->Integral()/hSgn_STATup->Integral());
  
  double xMIN = hBkg->GetBinLowEdge(1);
  double xMAX = hBkg->GetBinLowEdge(hBkg->GetNbinsX()+1);
  double xMIN2 = hDat->GetBinLowEdge(hDat->FindFirstBinAbove(0.5));
  double xMAX2 = hDat->GetBinLowEdge(hDat->FindLastBinAbove(0.5)+1);
  RooRealVar x("x","x",xMIN2,xMAX2);
  
  RooDataHist data("data","dataset with x",x,hDat);
  RooDataHist bkg("qcd","bkg with x",x,hBkg);
  RooDataHist sgn("signal","sgn with x",x,hSgn);
  
  RooHistPdf bkgPDF("bkgPDF","bkgPDF",x,bkg,0);
  RooHistPdf sgnPDF("sgnPDF","sgnPDF",x,sgn,0);
  
  RooRealVar f("f","f",0,0.,1.);
  
  RooAddPdf model("model","model",RooArgList(sgnPDF,bkgPDF),RooArgList(f));
  
  RooFitResult* r = model.fitTo(data,Save());
  r->Print("v");
  double N = hDat->Integral();
  double B = hBkg->Integral();
  double S = hSgn->Integral();
  double m = f.getVal();
  double e = f.getError();
  cout<<"k-factor = "<<k_factor<<endl;
  cout<<N<<" "<<B<<" "<<S<<endl;
  cout<<"Total cross section =       "<<N/lumi<<" pb"<<endl;
  cout<<"Model cross section =       "<<S/lumi<<" pb"<<endl;
  cout<<"Fitted signal strength =    "<<m*N/S<<endl;
  cout<<"Fitted signal error =       "<<e*N/S<<endl;
  double p = 0.95;
  double xup = (N/S)*(m+sqrt(2.)*e*TMath::ErfInverse((1-p)*TMath::Erf(m/e)+p));
  cout<<"Bayesian Upper limit =      "<<xup<<endl;
  RooPlot* frame1 = x.frame();  
//.........这里部分代码省略.........

void calCrystalDepo(TTree *t1041, int ThisModule = -99, bool IsBatch =0){
  if(IsBatch) gROOT->SetBatch();

  TString ThisModuleString = TString::Itoa(ThisModule,10);
  int UpOrDown = ThisModule/abs(ThisModule);
  float Scale = 1.;
  if(UpOrDown == -1) Scale = 0.5;
  cout << "Creating energy distributions of events with maximum deposition in module: " << ThisModuleString << endl;
//Define Histograms

//This histogram is a tool
  TH1F* MatrixDepo = new TH1F("MatrixDepo", "Deposition on Matrix", 68, -17, 17);

  TH1F* fracDep = new TH1F("fracDep", "Highest Amplitude Crystal/Total Deposition", 100, 0, 1.2);
  TH1F* MaxBin = new TH1F("MaxBin", "Deposition of Highest Amplitude Crystal", 100,400*Scale,7000*Scale);
  TH1F* Integral = new TH1F("Integral", "Integral", 100,4000*Scale,13000*Scale);
  TH1F* SigNoise = new TH1F("SigNoise", "Signa/Pedestal", 100, 0, 2.);

//

  gStyle->SetOptStat(0);
  

  Mapper *mapper=Mapper::Instance();

  TBEvent *event = new TBEvent();
  t1041->SetBranchAddress("tbevent", &event);
  bool haverechits=false;
  vector<TBRecHit> *rechits=0;

  if(t1041->GetListOfBranches()->FindObject("tbrechits")) {
    cout <<"found rechits"<<endl;
    t1041->SetBranchAddress("tbrechits",&rechits);
    haverechits=true;
  }

  Int_t start=0; Int_t end=t1041->GetEntries();  
  
  int nEntries=0;
//  end = 10;
  for (Int_t i=start; i<end; i++) {
    t1041->GetEntry(i);
    if (i==0) mapper->SetEpoch(event->GetTimeStamp());
    if(haverechits && rechits->size() < 1) continue;

    MatrixDepo->Reset();
    float maxDepoModuleID = -99;
    float maxDepo = 0;
    float secondMaxDepo = 0;
    float secondMaxDepoModuleID = -99;
    float totdep = 0;
    float totalnoise = 0;
    float sumrms = 0;
    for (Int_t j = 0; j < event->NPadeChan(); j++){
      if (haverechits && j>=(int)rechits->size()) break;
      double ped,sig, max, maxTime;
      int channelID;
      if (haverechits){
	TBRecHit &hit=rechits->at(j);
	ped=hit.Pedestal();
	sig=hit.NoiseRMS();
	max=hit.AMax();
	maxTime=hit.TRise();
	channelID=hit.GetChannelID();
      }

      double x,y;
      int moduleID,fiberID;
      mapper->ChannelID2ModuleFiber(channelID,moduleID,fiberID);
      mapper->ModuleXY(moduleID,x,y);

      MatrixDepo->Fill(moduleID, max);      
      if(UpOrDown > 0 && moduleID > 0)	totalnoise+=ped;
      if(UpOrDown < 0 && moduleID < 0)	totalnoise+=ped;
    }

    int minx = 0;
    int maxx = 0;
    if(UpOrDown > 0){
      minx = MatrixDepo->GetXaxis()->FindBin(0.);
      maxx = MatrixDepo->GetXaxis()->FindBin(17.);
      MatrixDepo->GetXaxis()->SetRange(minx, maxx);
    }
    if(UpOrDown < 0){
      minx = MatrixDepo->GetXaxis()->FindBin(-17.);
      maxx = MatrixDepo->GetXaxis()->FindBin(0.);	
      MatrixDepo->GetXaxis()->SetRange(minx, maxx);
    }

//    cout << "TOTAL NOISE: " << totalnoise << endl;

    float maxbin = MatrixDepo->GetBinContent(MatrixDepo->GetMaximumBin());
    float integral = MatrixDepo->Integral(minx,maxx);

    maxDepoModuleID = MatrixDepo->GetBinLowEdge(MatrixDepo->GetMaximumBin());

    if( abs(ThisModule) != 99 && maxDepoModuleID != ThisModule ) continue;
    if( abs(ThisModule) == 100 ){
	if( abs(maxDepoModuleID) != 6 &&
	    abs(maxDepoModuleID) != 7 &&
//.........这里部分代码省略.........

void plotFrac(TList* HistList, TH1F* compT, std::string name, bool fracAll){

  gROOT->SetBatch();
  system("mkdir -p plots/ada/fracs");
  system("mkdir -p plots/grad/fracs");

  std::string bdt;
  TString str = compT->GetName();
  if (str.Contains("ada")) bdt="ada";
  if (str.Contains("grad")) bdt="grad";

  int nHists=HistList->GetEntries();
  TH1F *comp = linearBin(compT);
  TH1F *uandd[nHists];
  
  gROOT->SetStyle("Plain");
  gROOT->ForceStyle();
  gStyle->SetOptStat(0);
  TCanvas *canv = new TCanvas("","",700,700);
  
  canv->Divide(1,2,0.1,0.1);
  canv->cd(1);
  gPad->SetPad(0.01,0.3,0.99,0.99);
  gPad->SetBottomMargin(0.0);
  gPad->SetLeftMargin(0.1);
  comp->SetLineColor(1);
  comp->SetFillColor(kGray);
  comp->Draw("hist");
  comp->SetTitle(("Up, down and interpolated templates for "+bdt+" "+name).c_str());
  comp->GetYaxis()->SetRangeUser(0.0,((TH1F*)HistList->At(0))->GetMaximum()+0.5);
  comp->GetYaxis()->SetTitle("Events / bin");
  comp->GetXaxis()->SetTitle("Category");
  int mass;
  if (name=="syst120") mass=120;
  if (name=="syst135") mass=135;
  TLegend *leg = new TLegend(0.6,0.5,0.88,0.88);
  leg->SetLineColor(0);
  leg->SetFillColor(0);
  if (name=="syst120" || name=="syst135"){
    leg->AddEntry(comp,Form("True %d signal binned at %d",mass,mass),"f");
  }
  
  canv->cd(2);
  gPad->SetPad(0.01,0.01,0.99,0.3);
  gPad->SetTopMargin(0.0);
  gPad->SetBottomMargin(0.2);
  gPad->SetLeftMargin(0.1);
  TF1 *line = new TF1("line","0.0",0.,comp->GetBinLowEdge(comp->GetNbinsX()+1));
  line->SetLineColor(kBlack);
  
  for (int i=0; i<nHists; i++){
    TH1F *temp = linearBin((TH1F*)HistList->At(i));
    temp->SetLineColor((i+1)*2);
    if (!fracAll && i==nHists-1) temp->SetLineColor(kBlack); 
    canv->cd(1);
    temp->DrawCopy("same e");
    uandd[i]= (TH1F*)temp->Clone();
    uandd[i]->Add(comp,-1);
    uandd[i]->Divide(temp);
    uandd[i]->Scale(-100.);
    //uandd[i]->Divide(comp);
    if (name=="syst120" || name=="syst135") {
      if (i==0) leg->AddEntry(uandd[i],Form("True %d signal binned at %d",mass+5,mass),"lep");
      if (i==1) leg->AddEntry(uandd[i],Form("True %d signal binned at %d",mass-5,mass),"lep");
      if (i==nHists-1) leg->AddEntry(uandd[i],Form("Interpolated %d signal",mass),"lep");
    }
    //uandd[i]->Scale(100.0);
    canv->cd(2);
    if (fracAll){
      if (i==0) uandd[i]->Draw("e");
      else {
        uandd[i]->Draw("same e");
        line->Draw("same");
      }
    }
    else {
      if (i==nHists-1) {
        uandd[i]->Draw("e");
        line->Draw("same");
      }
    }
    uandd[i]->SetTitle("");
    uandd[i]->GetYaxis()->SetLabelSize(0.08);
    uandd[i]->GetYaxis()->SetRangeUser(-50.,50.);
    uandd[i]->GetYaxis()->SetTitle("#frac{#Delta_{int}}{int} %");
    uandd[i]->GetYaxis()->SetTitleOffset(0.4);
    uandd[i]->GetYaxis()->SetTitleSize(0.08);
    uandd[i]->GetXaxis()->SetLabelSize(0.08);
    uandd[i]->GetXaxis()->SetTitle("BDT output");
    uandd[i]->GetXaxis()->SetTitleSize(0.08);
    canv->cd(1);
    if (name=="syst120" || name=="syst135") leg->Draw("same");
  }
  

  canv->Print(("plots/"+bdt+"/fracs/"+name+".png").c_str(),"png");
  delete canv;
  
  fracCalls++;
  
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      }
    }
    // fill data in binned TH1F
    hist= new TH1F ( variable, variable, (int)((xmax-xmin)/binwidth), xmin, xmax);
    TH1F* ori=(TH1F*)hist->Clone("original points");
    std::cout << "fine binned theory prediction has " << hist->GetNbinsX() << " bins" << std::endl;
    std::cout << "loaded values from .dat file: " << std::endl;
    // list all data values loaded and the corresponding bin
    for(int bin=1; bin<=Nbins; ++bin){
      double x=-999;
      double y=-999;
      // NB: choose if loaded data is interpreted as points with nothing 
      //     between (like kidonakis) or as integrated over the bin range (like ahrens)
      if(points){
	// check if you are still inside the array
	//std::cout << "data point " << bin-1 << "/" << sizeof(Xvalues)/sizeof(double) << std::endl;
	if(rawHist->GetPoint(bin-1, x, y)!=-1){
	  //x=Xvalues[bin-1]; // get value from data points
	  x+=0.5*binwidth;  // add half the binwidth to get the center of the bin
	}
      }
      else{
	x=hist->GetBinCenter(bin);  // get bin center
	y=rawHist->GetY()[bin-1];
      }
      if(x!=-999){
	std::cout << "data point: " << bin;
	std::cout << "(<x>=" << x << ")-> bin";
	// get bin in target (fine binned) plot
	int bin2=bin;
	if(points) bin2=hist->FindBin(x);
	//double y=Yvalues[bin2-1];
	std::cout << bin2 << " (";
	std::cout << hist->GetBinLowEdge(bin2) << ".." << hist->GetBinLowEdge(bin2+1);
	std::cout << "): " << y << std::endl;
	// fill target (fine binned) plot
	if(!points) hist->SetBinContent(bin2, y);
	// mark bins coming from the original prediction
	ori->SetBinContent(bin2, 1.);
	// -------------------------------
	// fit range without data entries
	// -------------------------------
	// NB: needed if loaded data is interpreted as points with nothing 
	//     between (like kidonakis)
	if(points){
	  // perform a linear fit wrt previous point
	  // get the two points (this bin and the previous one)
	  int binPrev= (bin==1&&variable=="topPt") ? 0 : hist->FindBin(Xvalues[bin-2]+0.5*binwidth);
	  double x2=-1;
	  double x1= 0;
	  double y2=-1;
	  double y1= 0.;
	  rawHist->GetPoint(bin-1, x2, y2);
	  x2+=0.5*binwidth;
	  if(bin==1&&variable=="topPt"){
	    y1=0;
	    x1=0;
	  }
	  else{
	    rawHist->GetPoint(bin-2, x1, y1);
	    x1+=0.5*binwidth;
	  }
	  // calculate linear funtion
	  double a=(y2-y1)/(x2-x1);
	  double b=y1-a*x1;
	  TF1* linInterpol=new TF1("linInterpol"+getTStringFromInt(bin), "[0]*x+[1]", x1, x2);

//.........这里部分代码省略.........
 //  cout<<"estimate of gain "<<gain_est<<", fit "<<fit->GetParameter(3)<<endl;
 //  cout<<"Shape Parameter "<<fit->GetParameter(4)<<endl;
 //  cout<<"SPE width "<<rms_estimate<<endl;
 //  cout<<"maxfitrange "<<maxfitrange<<endl;

            if(500<maxfitrange) maxfitrange = 500;
	    hspe->Fit(fit, "MNQL", "", minfitrange, maxfitrange);

	    //calculate NDOF of fit excluding bins with 0 entries
	    int myNDOF=-3; //three free parameters
	    for(int j=hspe->FindBin(minfitrange); j<=hspe->FindBin(maxfitrange); j++) { //loop through fitted spe bins
	      if(hspe->GetBinContent(j)) myNDOF++;
            } //loop through fitted spe bins


//	    cout<<"estimate of gain "<<gain_est<<", fit "<<fit->GetParameter(3)<<endl;
//	    cout<<"Shape Parameter "<<fit->GetParameter(4)<<endl;
//	    double SPE_rms = fit->GetParameter(3)/sqrt(fit->GetParameter(4));
//	    cout<<"SPE width "<<SPE_rms<<endl;

	    //cout<<"SPE width "<<fit->GetParameter(4)<<endl;
	    //cout<<"Fit normalization constant: estimate "<<mu<<" fit "<<fit->GetParameter(0)<<endl;
	    //cout<<spename<<endl;


	    //calculate means and integrals of the fit and data
	    double fint, fint_error, hint, favg, havg;
	    int temp_lowbin, temp_highbin;
	    temp_lowbin = hspe->FindBin(minfitrange);
	    temp_highbin = hspe->FindBin(maxfitrange);
	    hspe_unscaled->GetXaxis()->SetRangeUser(minfitrange, maxfitrange);
	    havg = hspe_unscaled->GetMean();
	    hint = hspe->Integral(temp_lowbin,temp_highbin,"width");
	    double min_frange = hspe->GetBinLowEdge(temp_lowbin);
	    favg = fit->Mean(min_frange, maxfitrange);
	    fint = fit->Integral(min_frange, maxfitrange);
	    //fint_error = fit->IntegralError(min_frange, maxfitrange);
	    
	    double PE5int = 0; //integral of events with >=5 PE
	    double PE5loc =  fped->GetParameter(1)+ 5*fit->GetParameter(3);
	    if(PE5loc>500) PE5int = 0;
	    else {
	      int PE5bin =  hspe_temp->FindBin(PE5loc);
	      temp_highbin = hspe_temp->FindBin(maxfitrange)-1;
	      PE5int =  hspe_temp->Integral(PE5bin,temp_highbin,"width");
	    }
	    int PE5flag = 0;
	    if(PE5int/hint>0.05) PE5flag = 1; //set flag if more than 5% of events in the fit correspond to >=5PE
        //=========================================    
            for(int i1=1;i1<hspe->GetNbinsX();i1++){
constants_file1<<HV<<"\t"<<iSpig<<"\t"<<i<<"\t"<<2.6*hspe->GetBinCenter(i1)<<"\t"<<hspe->GetBinContent(i1)<<"\t"<<fit->Eval(hspe->GetBinCenter(i1))<<"\n";
            }
        //=========================================    


	    //output calibrations constants
	    //constants_file<<endl<<"LED_amplitude HV Spigot Channel Ped_mean Ped_mean_err Ped_RMS  Ped_RMS_err SPEPeak_RMS SPEPeak_RMS_err Gain Gain_err Normalized_Chi2 MeanPE_fit MeanPE_fit_err MeanPE_estimate PE5flag Polya_shape Polya_shape_err Polya_mode"<<endl;
	    constants_file<<LED_amp<<" "<<HV<<" "<<iSpig<<" "<<(bb-1)*8+pmt<<" "<<scale*fped->GetParameter(1)<<" "<<scale*fped->GetParError(1)<<" "<<scale*fped->GetParameter(2)<<" "<<scale*fped->GetParError(2)<<" "<<scale*fit->GetParameter(3)/sqrt(fit->GetParameter(4))<<" "<<0<<" "<<scale*fit->GetParameter(3)*fC2electrons<<" "<<scale*fit->GetParError(3)*fC2electrons<<" "<<fit->GetChisquare()/myNDOF/*fit->GetNDF()*/<<" "<<fit->GetParameter(0)<<" "<<fit->GetParError(0)<<" "<<mu<<" "<<PE5flag<<" "<<fit->GetParameter(4)<<" "<<fit->GetParError(4)<<" "<<scale*(fit->GetParameter(4)-1.0)/fit->GetParameter(4)*fit->GetParameter(3)*fC2electrons<<endl;



//	    cout<<LED_amp<<" "<<HV<<" "<<iSpig<<" "<<QIECh[i]<<" "<<scale*fped->GetParameter(1)<<" "<<scale*fped->GetParError(1)<<" "<<scale*fped->GetParameter(2)<<" "<<scale*fped->GetParError(2)<<" "<<scale*fit->GetParameter(4)<<" "<<scale*fit->GetParError(4)<<" "<<scale*fit->GetParameter(3)*fC2electrons<<" "<<scale*fit->GetParError(3)*fC2electrons<<" "<<fit->GetChisquare()/fit->GetNDF()<<" "<<fit->GetChisquare()<<" "<<fit->GetNDF()<<" "<<myNDOF<<" "<<fit->GetParameter(0)<<" "<<fit->GetParError(0)<<" "<<mu<<endl;
	   
	    //extra_file<<endl<<"LED_amplitude HV Spigot Channel SignalAvg_inFitRange FitAvg_inFitRange SignalInt_inFitRange FitInt_inFitRange PEge5Int Gain(fC) Gain_err Normalized_Chi2 MeanPE_fit MeanPE_fit_err MeanPE_estimate PE5flag"<<endl;
	    //extra_file<<LED_amp<<" "<<HV<<" "<<iSpig<<" "<<QIECh[i]<<" "<<scale*havg<<" "<<scale*favg<<" "<<hint<<" "<<fint<<" "<<PE5int<<" "<<scale*fit->GetParameter(3)<<" "<<scale*fit->GetParError(3)<<" "<<fit->GetChisquare()/myNDOF<<" "<<fit->GetParameter(0)<<" "<<fit->GetParError(0)<<" "<<mu<<" "<<PE5flag<<endl;

//.........这里部分代码省略.........
    } else if(CYL) { // Cylindrical: The first bin with r>0 is actually the number 1 (not the 0).
      hDen1D = pData->GetH1SliceZ(pData->GetChargeFileName(1)->c_str(),"charge",1,Nbins,opth1.Data());
    } else { // 2D cartesian
      hDen1D = pData->GetH1SliceZ(pData->GetChargeFileName(1)->c_str(),"charge",-1,Nbins,opth1.Data());
    }
    hDen1D->SetName(hName); 
    
    if(opt.Contains("comov"))
      hDen1D->GetXaxis()->SetTitle("k_{p}#zeta");
    else
      hDen1D->GetXaxis()->SetTitle("k_{p}z");
  
    hDen1D->GetYaxis()->SetTitle("n_{b}/n_{0}");
  }

  // On-axis beam density vs \zeta vs time! _________________________________
  TH2F *hDen1DvsTime = NULL; 
  if(hDen1D) {
    char hName[24];
    sprintf(hName,"hDen1DvsTime");
    TH2F *hDen1DvsTimeOld = (TH2F*) ifile->Get(hName);

    Int_t nBins   = 1;
    Float_t edge0 = Time-0.5;
    Float_t edge1 = Time+0.5;
    if(hDen1DvsTimeOld!=NULL) {
      nBins = hDen1DvsTimeOld->GetNbinsX()+1;
      Float_t binwidth =  (Time - hDen1DvsTimeOld->GetXaxis()->GetBinCenter(1))/(nBins-1);
      edge0 = hDen1DvsTimeOld->GetXaxis()->GetBinCenter(1) - binwidth/2.;
      edge1 = Time + binwidth/2.;
    }
    hDen1DvsTime = new TH2F("temp","",nBins,edge0,edge1,
		       	hDen1D->GetNbinsX(),
			hDen1D->GetBinLowEdge(1),
			hDen1D->GetBinLowEdge(hDen1D->GetNbinsX()+1));
    
    for(Int_t ix=1;ix<hDen1DvsTime->GetNbinsX();ix++) {
      for(Int_t iy=1;iy<hDen1DvsTime->GetNbinsY();iy++) {
	hDen1DvsTime->SetBinContent(ix,iy,hDen1DvsTimeOld->GetBinContent(ix,iy));
      }
    }  
    delete hDen1DvsTimeOld;
  
    // Fill last bin with the newest values.
    for(Int_t iy=1;iy<=hDen1D->GetNbinsX();iy++) {
      hDen1DvsTime->SetBinContent(nBins,iy,hDen1D->GetBinContent(iy));
    }   

    hDen1DvsTime->GetZaxis()->SetTitle("n_{b}/n_{0}");
    hDen1DvsTime->GetYaxis()->SetTitle("k_{p}#zeta");
    hDen1DvsTime->GetXaxis()->SetTitle("k_{p}z");
    hDen1DvsTime->GetZaxis()->CenterTitle();
    hDen1DvsTime->GetYaxis()->CenterTitle();
    hDen1DvsTime->GetXaxis()->CenterTitle();
    hDen1DvsTime->SetName(hName);

    // Change the range of z axis 
    Float_t Denmax = hDen1DvsTime->GetMaximum();
    hDen1DvsTime->GetZaxis()->SetRangeUser(0,Denmax); 
    hDen1DvsTime->Write(hName,TObject::kOverwrite);

  }

  // RMS (vs z) of the beam's charge distribution: 
  TProfile *hDen2Dprof = NULL;
  TH1F *hRms = NULL;

void makeFitPlotSFJVF()
{
   gROOT->SetBatch();
   gROOT->SetStyle("Plain");
   
   gStyle->SetOptStat(0);

   gROOT->ProcessLine(".x def.C");
   gROOT->ProcessLine(".x common.C");

   bool doSubPlot = 0;

   std::string sel = "MV160";
   std::string var = "Sd0t1";

   const int nf = 2;
   std::string fv[nf] = {"","_JVF0"};

   double sf[100][100];
   double sferr[100][100];
   double sfX[100];
   double sfXerr[100];
   
   TCanvas *c1 = new TCanvas("c1","c1",0,0,600,500);
   c1->Draw();
   c1->cd();

   TPad *c1_1;
   if( doSubPlot )
     {	
	c1->Range(0,0,1,1);
	c1_1 = new TPad("c1_1","main",0.01,0.30,0.99,0.99);
	c1_1->Draw();
	c1_1->cd();
     }   
   
   gStyle->SetHistTopMargin(0);

   TLegend *leg = new TLegend(0.70,0.90,0.90,0.70);
   leg->SetFillColor(253);
   leg->SetBorderSize(0);

   int nbins = 0;
   int im = 0;

   TH1F *hc;

   TH1F *hsf[nf];
   
   for(int f=0;f<nf;f++)
     {
	std::string fname = "results/sfoverlay/fit_EFF_beauty_"+var+"_"+sel+fv[f]+".root";
	TFile *file = new TFile(fname.c_str());
	TH1F *h = (TH1F*)file->Get("h1c");
	
	nbins = h->GetXaxis()->GetNbins();
	for(int ib=1;ib<nbins+1;ib++)
	  {
	     double cont = h->GetBinContent(ib);
	     double err = h->GetBinError(ib);
	     sfXerr[ib-1] = h->GetBinWidth(ib)/2.0;
	     sfX[ib-1] = h->GetBinLowEdge(ib) + sfXerr[ib-1];
	     sf[ib-1][im] = cont;
	     sferr[ib-1][im] = err;
	  }	     
	
	if( f == 0 ) {h->SetMarkerStyle(20);h->SetMarkerColor(kRed);h->SetLineColor(kRed);}
	if( f == 1 ) {h->SetMarkerStyle(25);h->SetMarkerColor(kBlack);h->SetLineColor(kBlack);}
	if( f == 2 ) {h->SetMarkerStyle(22);h->SetMarkerColor(kBlue);h->SetLineColor(kBlue);}
	if( f == 3 ) {h->SetMarkerStyle(23);h->SetMarkerColor(kMagenta);h->SetLineColor(kMagenta);}
	
	std::string tit = "";
	if( f == 0 ) tit = "JVF > 0.5";
	if( f == 1 ) tit = "JVF > 0";
	if( f == 2 ) tit = "S_{d_{0}}(t_{3})";
	if( f == 3 ) tit = "S_{d_{0}}(#mu)";
	leg->AddEntry(h,tit.c_str(),"lep");

	if( f == 0 ) h->Draw("e1");
	else h->Draw("e1 same");

	ATLASLabel(0.20,0.85,"Internal",1);
	
	h->GetYaxis()->SetRangeUser(0.3,1.6);
	h->GetXaxis()->SetRangeUser(0.0,0.04);
	
	im++;
	
	std::string hname = "hsf"+std::string(Form("%d",f));
	hsf[f] = (TH1F*)h->Clone(hname.c_str());
     }   

   leg->Draw();
   
   gPad->RedrawAxis("g");
      
/*   if( doSubPlot )
     {
	std::string foutStr = "results/fit_EFF_beauty_SFOVERLAY_"+sel+".root";
	TFile *fout = new TFile(foutStr.c_str(),"RECREATE");
//.........这里部分代码省略.........

        void Draw(TCanvas *c1, const int &rebin=1, const bool &showErr=true) {

            gStyle->SetOptStat(0);
            c1->cd();
            c1->Clear();

            TPad *pad1;
            pad1 = new TPad("pad1","pad1",0,1-0.614609572,1,1);
            pad1->SetTopMargin(0.0983606557);
            pad1->SetBottomMargin(0.025);
            pad1->Draw();
            pad1->cd();

            RebinHists(rebin);
            THStack *hstack = GetStack(c1->GetLogy());
            TH1F *signal = GetSignalHist();
            TH1F *data   = GetDataHist();

            if(c1->GetLogy()) gPad->SetLogy();
            hstack->GetHistogram()->SetLabelSize(0.00,"X");
            hstack->GetHistogram()->SetLabelSize(0.06,"Y");
            hstack->GetHistogram()->SetTitleSize(0.06,"XY");
            hstack->Draw("hist");
            if(signal && ! _stackSignal) signal->Draw("hist,same");
            if(data)     data->Draw("ep,same");
            DrawLabels();
            pad1->GetFrame()->DrawClone();


            TH1F *summed = GetSummedMCHist();

            TH1F *rdat = (TH1F*)data->Clone("rdat");   
            if(gROOT->FindObject("rref")) gROOT->FindObject("rref")->Delete();
            TH1F *rref = new TH1F("rref","rref",
                summed->GetNbinsX(),
                summed->GetBinLowEdge(1),
                summed->GetBinLowEdge(summed->GetNbinsX()+1)
            );
            for (int i = 1, n = rref->GetNbinsX(); i <= n+1; ++i) {
                rref->SetBinContent(i,summed->GetBinContent(i));
                rref->SetBinError(i,summed->GetBinError(i));
            }
            rref->SetTitle("");
            rref->SetLineWidth(0);
            rref->SetFillColor(kGray+1);
            rref->SetFillStyle(1001);
            double absmax = 0;
            for (int i = 0, n = rdat->GetNbinsX(); i <= n+1; ++i) {
                double scale = rref->GetBinContent(i);
                if (scale == 0) {
                    rdat->SetBinContent(i, 0);
                    rref->SetBinContent(i, 10000);
                    rdat->SetBinError(i, 0);
                    rref->SetBinError(i, 0);
                } else {
                    rdat->SetBinContent(i, rdat->GetBinContent(i)/scale);
                    rref->SetBinContent(i, rref->GetBinContent(i)/scale);
                    rdat->SetBinError(i, rdat->GetBinError(i)/scale);
                    rref->SetBinError(i, rref->GetBinError(i)/scale);
                    double mymax = TMath::Max(1.2*fabs(rdat->GetBinContent(i)-1)+1.4*rdat->GetBinError(i), 2.0*rref->GetBinError(i));
                    absmax = TMath::Max(mymax, absmax);
                }
            }

            c1->cd();
            TPad *pad2 = new TPad("pad2","pad2",0,0,1,1-0.614609572);
            pad2->SetTopMargin(0.0261437908);
            pad2->SetBottomMargin(0.392156863);
            pad2->Draw();
            pad2->cd();

            TLine *line = new TLine(rref->GetXaxis()->GetXmin(), 1.0, rref->GetXaxis()->GetXmax(), 1.0);
            line->SetLineColor(kBlack);
            line->SetLineWidth(1);
            line->SetLineStyle(1);

            if(showErr) {
                rref->GetYaxis()->SetRangeUser(TMath::Max(-1.,1.-absmax), TMath::Min(3.,absmax+1.));
                AxisFonts(rref->GetXaxis(), "x", hstack->GetXaxis()->GetTitle());
                rref->GetYaxis()->SetTitle("data/mc");
                rref->GetYaxis()->SetLabelSize(0.10);
                rref->GetYaxis()->SetTitleSize(0.10);
                rref->GetYaxis()->SetTitleOffset(0.85);
                rref->GetXaxis()->SetLabelSize(0.10);
                rref->GetXaxis()->SetTitleSize(0.10);
                rref->GetXaxis()->SetTitleOffset(1.5);
                rref->Draw("E2"); 
                rdat->SetMarkerStyle(20);
                rdat->Draw("E SAME p");
                line->Draw("SAME"); 
                c1->Update();
                pad2->GetFrame()->DrawClone();

            } else {
                rdat->GetYaxis()->SetRangeUser(TMath::Max(-1.,1.-absmax), TMath::Min(3.,absmax+1.));
                rdat->GetYaxis()->SetTitle("data/mc");
                rdat->GetYaxis()->SetLabelSize(0.10);
                rdat->GetYaxis()->SetTitleSize(0.10);
                rdat->GetYaxis()->SetTitleOffset(0.85);
                AxisFonts(rdat->GetXaxis(), "x", hstack->GetXaxis()->GetTitle());
//.........这里部分代码省略.........