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

//.........这里部分代码省略.........
int npt_a = sizeof(selptbin)/sizeof(double)-1;
for(int icent_a=0;icent_a<ncent+1;icent_a++){
for(int icent_b=0; icent_b<ncent+1; icent_b++)
if(selcentbin[icent_a] == centbin[icent_b]) break;
int xcentmin = icent_b;
for(int icent_b=0; icent_b<ncent+1; icent_b++)
if(selcentbin[icent_a+1] == centbin[icent_b]) break;
int xcentmax = icent_b;
if((xcentmin == ncent+1) || (xcentmax == ncent+1)) exit(0);
for(int ipt_a=0;ipt_a<npt_a;ipt_a++){
for(int ipt_b=0; ipt_b<npt+1; ipt_b++)
if(selptbin[ipt_a] == ptbin[ipt_b]) break;
int xptmin = ipt_b;
for(int ipt_b=0; ipt_b<npt+1; ipt_b++)
if(selptbin[ipt_a+1] == ptbin[ipt_b]) break;
int xptmax = ipt_b;
if((xptmin == npt+1) || (xptmax == npt+1)) exit(0);
cout<<xcentmin<<"\t"<<xcentmax<<"\t"<<endl;
cout<<xptmin<<"\t"<<xptmax<<"\t"<<endl;
kforebbcw_In = (TH1F*)kforebbcw[xcentmin][xptmin]->Clone();
hforebbcw_In = (TH1F*)hforebbcw[xcentmin][xptmin]->Clone();
kbackbbcw2_In = (TH1F*)kbackbbcw2[xcentmin][xptmin]->Clone();
kforebbcw_In->Reset();
hforebbcw_In->Reset();
kbackbbcw2_In->Reset();
for(int icent=xcentmin; icent<xcentmax; icent++){
for(int ipt=xptmin; ipt<xptmax; ipt++){
kforebbcw_In->Add(kforebbcw[icent][ipt]);
hforebbcw_In->Add(hforebbcw[icent][ipt]);
kbackbbcw2_In->Add(kbackbbcw2[icent][ipt]);
}
}
kforebbcw_In->Rebin(2);
hforebbcw_In->Rebin(2);
kbackbbcw2_In->Rebin(2);


TH1F* hpp;
TH1F* hbackpp;
  hpp = (TH1F*)kforebbcw_In->Clone();
  hbackpp = (TH1F*)kbackbbcw2_In->Clone();
  float nbackpp = hbackpp->Integral()/2.0/PI;
  float nforepp = hpp->Integral()/2.0/PI;
  //float ntrig0 = ptforedis_0->Integral(11,30);
   for(int i=0; i<20; i++){
     float pp_cont = 1.0*hpp->GetBinContent(i+1);
     //float pp0_err = 1.0*hpp->GetBinError(i+1);
     float weight2 = sqrt(1.0*hforebbcw_In->GetBinContent(i+1));

     float backpp_cont = 1.0*hbackpp->GetBinContent(i+1);
    
     float con = pp_cont/backpp_cont*nbackpp/nforepp;
     float err = weight2/backpp_cont*nbackpp/nforepp;

     hpp->SetBinContent(i+1, con);
     hpp->SetBinError(i+1, err);

   }
  TF1 *fun0 = new TF1("fun0","[0]*(1+2*[1]*cos(x)+2*[2]*cos(2*x)+2*[3]*cos(3*x)+2*[4]*cos(4*x))", -0.5*PI, 1.5*PI);

  fun0->SetLineColor(1);
  hpp->Fit("fun0","NORQ");

  TF1 *fun1 = new TF1("fun1","[0]*(1+2*[1]*cos(x))",   -0.5*PI, 1.5*PI);
  TF1 *fun2 = new TF1("fun2","[0]*(1+2*[1]*cos(2*x))", -0.5*PI, 1.5*PI);
  TF1 *fun3 = new TF1("fun3","[0]*(1+2*[1]*cos(3*x))", -0.5*PI, 1.5*PI);
  TF1 *fun4 = new TF1("fun4","[0]*(1+2*[1]*cos(4*x))", -0.5*PI, 1.5*PI);

const double ptmean[npt] = {0.360943, 0.691833, 1.1911, 1.69654, 2.20117, 2.70571, 3.2097, 3.71372, 4.21814, 4.72014};
  fout<<fun0->GetParameter(1)<<" "<<fun0->GetParError(1)<<" "
      <<fun0->GetParameter(2)<<" "<<fun0->GetParError(2)<<" "//<<0.05*fun0->GetParameter(2)<<" "<<endl;
      <<fun0->GetParameter(3)<<" "<<fun0->GetParError(3)<<endl;
}
}
  fout.close();

/*  
  cout<<"*************** v2 ***********"<<endl;
  float v2_0 = funvn0->GetParameter(1)/(zym_pp0 + funvn0->GetParameter(0));
  float v2_1 = funvn1->GetParameter(1)/(zym_pp1 + funvn1->GetParameter(0));
  float v2_2 = funvn2->GetParameter(1)/(zym_pp2 + funvn2->GetParameter(0));
  float v2_3 = funvn3->GetParameter(1)/(zym_pp3 + funvn3->GetParameter(0));
  

  cout<<v2_0<<" "<<v2_1<<" "<<v2_2<<" "<<v2_3<<endl;
 
  
  cout<<funvn0->GetParameter(0)*funvn0->GetParameter(1)<<" "
      <<funvn1->GetParameter(0)*funvn1->GetParameter(1)<<" "
      <<funvn2->GetParameter(0)*funvn2->GetParameter(1)<<" "
      <<funvn3->GetParameter(0)*funvn3->GetParameter(1)<<endl;
  

  cout<<"*************** v2 ***********"<<endl;
  cout<<funvn0->GetParameter(2)<<" "<<funvn1->GetParameter(2)<<" "<<funvn2->GetParameter(2)<<" "<<funvn3->GetParameter(2)<<endl;

  cout<<"*************** v3 ***********"<<endl;
  cout<<funvn0->GetParameter(3)<<" "<<funvn1->GetParameter(3)<<" "<<funvn2->GetParameter(3)<<" "<<funvn3->GetParameter(3)<<endl;
*/
}

//.........这里部分代码省略.........
    leg->AddEntry(ggH  , TString::Format("%.0f#timesH(125 GeV)#rightarrow#tau#tau", SIGNAL_SCALE) , "L" );
  }
  else{
    leg->AddEntry(ggH  , "SM H(125 GeV)#rightarrow#tau#tau" , "L" );
  }
#endif
#endif
#ifdef ASIMOV
  leg->AddEntry(data , "sum(bkg) + H(125)"              , "LP");
#else
  leg->AddEntry(data , "Observed"                       , "LP");
#endif
  leg->AddEntry(Ztt  , "Z#rightarrow#tau#tau"           , "F" );
  leg->AddEntry(EWK  , "Z#rightarrow ee"                , "F" );
  leg->AddEntry(EWK1 , "W+jets"                         , "F" );
  leg->AddEntry(ttbar, "t#bar{t}"                       , "F" );
  leg->AddEntry(Fakes, "QCD"                            , "F" );
  $ERROR_LEGEND
  leg->Draw();

  /*
    Ratio Data over MC
  */
  TCanvas *canv0 = MakeCanvas("canv0", "histograms", 600, 400);
  canv0->SetGridx();
  canv0->SetGridy();
  canv0->cd();

  TH1F* model = (TH1F*)Ztt ->Clone("model");
  TH1F* test1 = (TH1F*)data->Clone("test1"); 
  for(int ibin=0; ibin<test1->GetNbinsX(); ++ibin){
    //the small value in case of 0 entries in the model is added to prevent the chis2 test from failing
    model->SetBinContent(ibin+1, model->GetBinContent(ibin+1)>0 ? model->GetBinContent(ibin+1)*model->GetBinWidth(ibin+1) : 0.01);
    model->SetBinError  (ibin+1, CONVERVATIVE_CHI2 ? 0. : model->GetBinError  (ibin+1)*model->GetBinWidth(ibin+1));
    test1->SetBinContent(ibin+1, test1->GetBinContent(ibin+1)*test1->GetBinWidth(ibin+1));
    test1->SetBinError  (ibin+1, test1->GetBinError  (ibin+1)*test1->GetBinWidth(ibin+1));
  }
  double chi2prob = test1->Chi2Test      (model,"PUW");        std::cout << "chi2prob:" << chi2prob << std::endl;
  double chi2ndof = test1->Chi2Test      (model,"CHI2/NDFUW"); std::cout << "chi2ndf :" << chi2ndof << std::endl;
  double ksprob   = test1->KolmogorovTest(model);              std::cout << "ksprob  :" << ksprob   << std::endl;
  double ksprobpe = test1->KolmogorovTest(model,"DX");         std::cout << "ksprobpe:" << ksprobpe << std::endl;  

  std::vector<double> edges;
  TH1F* zero = (TH1F*)ref->Clone("zero"); zero->Clear();
  TH1F* rat1 = (TH1F*)data->Clone("rat1"); 
  for(int ibin=0; ibin<rat1->GetNbinsX(); ++ibin){
    rat1->SetBinContent(ibin+1, Ztt->GetBinContent(ibin+1)>0 ? data->GetBinContent(ibin+1)/Ztt->GetBinContent(ibin+1) : 0);
    rat1->SetBinError  (ibin+1, Ztt->GetBinContent(ibin+1)>0 ? data->GetBinError  (ibin+1)/Ztt->GetBinContent(ibin+1) : 0);
    zero->SetBinContent(ibin+1, 0.);
    zero->SetBinError  (ibin+1, Ztt->GetBinContent(ibin+1)>0 ? Ztt ->GetBinError  (ibin+1)/Ztt->GetBinContent(ibin+1) : 0);
  }
  for(int ibin=0; ibin<rat1->GetNbinsX(); ++ibin){
    if(rat1->GetBinContent(ibin+1)>0){
      edges.push_back(TMath::Abs(rat1->GetBinContent(ibin+1)-1.)+TMath::Abs(rat1->GetBinError(ibin+1)));
      // catch cases of 0 bins, which would lead to 0-alpha*0-1
      rat1->SetBinContent(ibin+1, rat1->GetBinContent(ibin+1)-1.);
    }
  }
  float range = 0.1;
  std::sort(edges.begin(), edges.end());
  if (edges[edges.size()-2]>0.1) { range = 0.2; }
  if (edges[edges.size()-2]>0.2) { range = 0.5; }
  if (edges[edges.size()-2]>0.5) { range = 1.0; }
  if (edges[edges.size()-2]>1.0) { range = 1.5; }
  if (edges[edges.size()-2]>1.5) { range = 2.0; }
  rat1->SetLineColor(kBlack);

//.........这里部分代码省略.........
  }

  if(ppPbPb){
    for(int i=0;i<tB->GetEntries();i++){
      tB->GetEntry(i);
      if(fabs(jtetaB)<2 && binB>=cbinlo && binB<cbinhi && abs(refparton_flavorForBB)==5)
	hB->Fill(refptB,jtptB/refptB,weightB);
    }
    for(int i=0;i<tC->GetEntries();i++){
      tC->GetEntry(i);
      if(fabs(jtetaC)<2 && binC>=cbinlo && binC<cbinhi && abs(refparton_flavorForBC)==4)
	hC->Fill(refptC,jtptC/refptC,weightC);
    }
  }

  
 
  hL->SetMinimum(0.);
  
  hL->SetLineColor(kBlue);
  hB->SetLineColor(kRed);
  hC->SetLineColor(kGreen);

  hL->SetMarkerColor(kBlue);
  hB->SetMarkerColor(kRed);
  hC->SetMarkerColor(kGreen);
  
  //hL->SetMarkerStyle(4);
  //hB->SetMarkerStyle(4);
  //hC->SetMarkerStyle(4);
  
  hL->SetXTitle("genJet p_{T} (GeV/c)");
  hL->SetYTitle("<reco p_{T} / gen p_{T} >");

  hL->GetXaxis()->SetRangeUser(50.,199.999);
  hL->GetYaxis()->SetRangeUser(0.5,1.05);
  
  TCanvas *c1=new TCanvas("c1","c1",800,600);
  c1->SetGridx(1);
  c1->SetGridy(1);

  hL->Draw("e1");
  hB->Draw("e1,same");
  hC->Draw("e1,same");

  TLegend *leg=new TLegend(0.4,0.15,0.9,0.45);
  leg->SetBorderSize(0);
  leg->SetFillStyle(0);
  if(ppPbPb&&cbinlo==0&&cbinhi==40)leg->SetHeader("Pythia+Hydjet, 0-100%");
  leg->AddEntry(hL,"Inclusive jets","pl");
  leg->AddEntry(hC,"c-jets","pl");
  leg->AddEntry(hB,"b-jets","pl");
  leg->Draw();

  TCanvas *c2=new TCanvas("c2","c2",1);
  /*
  TH1F *hL2 = (TH1F*)hL->Clone("hL2");
  TH1F *hB2 = (TH1F*)hB->Clone("hB2");
  hL2->Add(hB2,-1);
  hL2->Draw();
  */

  TH1F  *hcorr = new TH1F("hcorr","hcorr",250,50,300);
  hcorr->Sumw2();

  for(int i=0;i<hL->GetNbinsX();i++){
    cout<<" b resp "<<hB->GetBinContent(i+1)<<endl;
    cout<<" l resp "<<hL->GetBinContent(i+1)<<endl;
    cout<<" l offset "<<1.-hL->GetBinContent(i+1)<<endl;
    cout<<" corrected b resp "<<hB->GetBinContent(i+1)+1.-hL->GetBinContent(i+1)<<endl;
    float jesOffset = 1.-hL->GetBinContent(i+1);

    hcorr->SetBinContent(i+1,hB->GetBinContent(i+1)+jesOffset);

    hcorr->SetBinError(i+1,sqrt(hB->GetBinError(i+1)*hB->GetBinError(i+1)+hL->GetBinError(i+1)*hL->GetBinError(i+1)));


  }

  hcorr->SetMinimum(0.5);
  hcorr->SetMaximum(1.1);
      
  hcorr->SetLineColor(kRed);
  hcorr->SetMarkerColor(kRed);
  hcorr->SetMarkerStyle(4);
  hcorr->Draw();

  TF1 *fCorr = new TF1("fCorr","[0]+[1]*log(x)+[2]*log(x)*log(x)",50,300);
  fCorr->SetLineWidth(1);
  fCorr->SetLineColor(kBlue);
  hcorr->Fit(fCorr);

  TFile *fout;
  if(ppPbPb) fout =new TFile(Form("bJEShistos/bJetScale_PbPb_Cent_fineBin_%d_%d.root",cbinlo,cbinhi),"recreate");
  else fout =new TFile("bJEShistos/bJetScale_PP_fineBin.root","recreate");
  hcorr->Write();
  fCorr->Write();
  fout->Close();

}

void alicePlots(){
    
    TFile* alice = new TFile("~/Downloads/HEPData-ins1288320-v1-root.root");
    alice->cd("Table 16");
    TGraph* aliceData = Graph1D_y1;
    TH1F* hist = Hist1D_y1;
    TH1F* stat = Hist1D_y1_e1;
    TH1F* syst = Hist1D_y1_e2;
    TGraphAsymmErrors* graph2 = (TGraphAsymmErrors*)aliceData->Clone("graph2");

    Int_t numPts = aliceData->GetN();
    Double_t x, y;
    for(int i = 0; i<numPts; i++){
        aliceData->GetPoint(i, x, y);
        aliceData->SetPoint(i, x, (y - 0.89581));
        graph2->SetPoint(i, x, (y- 0.89581));
        hist->SetBinContent(i+1, hist->GetBinContent(i+1) - 0.89581);
        hist->SetBinError(i+1, stat->GetBinContent(i+1));
        graph2->SetPointEXhigh(i, 0.1);
        graph2->SetPointEXlow(i, 0.1);
    }

    graph2->SetLineColor(kBlue-10);
    graph2->SetLineWidth(2);
    graph2->SetMarkerColor(kBlue-10);
    graph2->SetFillColor(kBlue-10);
    hist->SetLineColor(kBlue-2);
    hist->SetLineWidth(2);

    aliceData->SetTitle("");
    aliceData->GetYaxis()->SetTitle("Mass - Vacuum Mass (GeV/c^{2})");
    aliceData->GetYaxis()->SetTitleSize(0.06);
    aliceData->GetYaxis()->SetLabelSize(0.04);
    aliceData->GetYaxis()->SetTitleOffset(1.65);
    aliceData->GetYaxis()->SetTitleFont(42);
    aliceData->GetYaxis()->SetLabelFont(42);
    aliceData->GetXaxis()->SetTitle("p_{T} (GeV/c)");
    aliceData->GetXaxis()->SetTitleSize(0.06);
    aliceData->GetXaxis()->SetLabelSize(0.05);
    aliceData->GetXaxis()->SetTitleFont(42);
    aliceData->GetXaxis()->SetLabelFont(42);
    aliceData->SetMarkerStyle(29);
    aliceData->SetMarkerSize(2.5);
    aliceData->SetMarkerColor(kBlue-2);
    aliceData->SetLineColor(kBlue-2);

    aliceData->GetYaxis()->SetRangeUser(-0.02, 0.015);
    aliceData->GetXaxis()->SetRangeUser(0, 5);

    TFile* phsd = new TFile("~/utaustin/resonancefits/finalplotting/20170721_KKbarAdded2_fixedwidth42_recon_pf100_scaled_error05.root");
    TH1D* mass = phsd->Get("kstar0mass");
    mass->SetName("mass");
    mass->SetMarkerStyle(26);
    mass->SetMarkerSize(2.5);
    mass->SetMarkerColor(2);
    mass->SetLineColor(2);

    TF1* line = new TF1("line", "[0]", 0.0, 5.0);
    line->SetParameter(0, 0.0);
    line->SetLineColor(1);
    line->SetLineStyle(7);
    line->SetLineWidth(3);

    for(int j = 0; j<mass->GetNbinsX(); j++){
        mass->SetBinContent(j+1, (mass->GetBinContent(j+1) - 0.892));
    }


    TFile* phsd2 = new TFile("~/utaustin/resonancefits/finalplotting/20170616_KKbarAdded2_fixedwidth_recon_pf100_scaled_error05.root");
    TH1D* mass2 = phsd2->Get("kstar0mass");
    mass2->SetName("mass2");
    mass2->SetMarkerStyle(22);
    mass2->SetMarkerSize(2.5);
    mass2->SetMarkerColor(2);
    mass2->SetLineColor(2);
    for(int j = 0; j<mass2->GetNbinsX(); j++){
        mass2->SetBinContent(j+1, (mass2->GetBinContent(j+1) - 0.892));
    }



    TExec *exec1 = new TExec("exec1", "gStyle->SetErrorX(0.1)");
    TExec *exec2 = new TExec("exec2", "gStyle->SetErrorX(0.5)");
    
    TCanvas *c = new TCanvas ("c", "c", 50, 50, 650, 600);
    c->cd()->SetMargin(0.1997, 0.0369, 0.1396, 0.0681);
    aliceData->Draw("APX");
    //exec1->Draw();
    graph2->Draw("SAME P2");
    //exec2->Draw();
    hist->Draw("SAME E1");
    line->Draw("SAME");
    mass2->Draw("SAME P E1");
    mass->Draw("SAME P E1");
    aliceData->Draw("SAME PX");
    

    TLegend* legend = new TLegend(0.5836, 0.1815, 0.9489, 0.3438);
    legend->SetMargin(0.2);
    legend->SetTextSizePixels(20);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            bin1 = pt_bound[j]; bin2 = pt_bound[j+1];
            if( AccJpsi == 1 && (AllCut == 1) && (vars >= bin1 && vars < bin2) && TMath::Abs(JpsiRap) <= maxRap && TMath::Abs(JpsiRap) >= minRap) {
              hRecoDiMuon[j]->Fill(JpsiMass,RecWeight);
            }
          }
          if(iSpec == 1){ 
            bin1 = rap_bound[j]; bin2 = rap_bound[j+1];
            if( AccJpsi == 1 && (AllCut == 1) && (vars >= bin1 && vars < bin2) && JpsiPt >= minPt && JpsiPt <= maxPt) {
              hRecoDiMuon[j]->Fill(JpsiMass,RecWeight);
            }
          }
        }
      }
    } //rec tree loop ends

    //======================  File loop Starts ============================

    ///////////////////////////////////////////////////////////////////
    cout<< " adding "<<endl;
    char gtmp[512], gtmp1[512];
    char rtmp[512], rtmp1[512];

    for(int i = 0; i < nBins; i++){
      if(iSpec == 1) sprintf(gtmp,"hGenDiMuon_Pt_%0.1f_%0.1f",pt_bound[i],pt_bound[i+1]);
      if(iSpec == 2) sprintf(gtmp,"hGenDiMuon_Rap_%0.1f_%0.1f",rap_bound[i],rap_bound[i+1]);
      if(iSpec == 1) sprintf(gtmp1,"plots/Gen/hGenDiMuon_Pt_%0.1f_%0.1f.png",pt_bound[i],pt_bound[i+1]);
      if(iSpec == 2) sprintf(gtmp1,"plots/Gen/hGenDiMuon_Rap_%0.1f_%0.1f.png",rap_bound[i],rap_bound[i+1]);
      if(iSpec == 1) sprintf(rtmp,"hRecDiMuon_Pt_%0.1f_%0.1f",pt_bound[i],pt_bound[i+1]);
      if(iSpec == 2) sprintf(rtmp,"hRecDiMuon_Rap_%0.1f_%0.1f",rap_bound[i],rap_bound[i+1]);
      if(iSpec == 1) sprintf(rtmp1,"plots/Rec/hRecDiMuon_Pt_%0.1f_%0.1f.png",pt_bound[i],pt_bound[i+1]);
      if(iSpec == 2) sprintf(rtmp1,"plots/Rec/hRecDiMuon_Rap_%0.1f_%0.1f.png",rap_bound[i],rap_bound[i+1]);
      hGenDiMuon[i]->SetName(gtmp);
      hRecoDiMuon[i]->SetName(rtmp);
    }

    cout<<"Starts to calculate efficiency"<<endl;
    dataFile<<""<<endl;
    //=====================Loop for eff========================================================================================//
    //define stuff here for error on weighted samples
    for(int i = 0; i < nBins; i++){
      int Gbinlow =hGenDiMuon[i]->GetXaxis()->FindBin(2.0);//2.95
      int Gbinhi=hGenDiMuon[i]->GetXaxis()->FindBin(4.0);//2.95

      int Rbinlow =hRecoDiMuon[i]->GetXaxis()->FindBin(2.95);//2.95
      int Rbinhi=hRecoDiMuon[i]->GetXaxis()->FindBin(3.25);//2.95

      genNo[i] = hGenDiMuon[i]->IntegralAndError(Gbinlow, Gbinhi, genErr[i]); 
      recoNo[i] = hRecoDiMuon[i]->IntegralAndError(Rbinlow, Rbinhi, recoErr[i]);
      //calculate Eff         
      if(genNo[i] == 0 || recoNo[i] == 0) {
        eff[i] = 0;
        effErr[i] = 0;
      }else{
        eff[i] = recoNo[i]/genNo[i]; 

        effErr[i] = recoNo[i]/genNo[i]*TMath::Sqrt(genErr[i]*genErr[i]/(genNo[i]*genNo[i]) + recoErr[i]*recoErr[i]/(recoNo[i]*recoNo[i]));

        //error without weight
        dataFile<<" Bin ["<<i<<"] - "<< " Reco Jpsi : "<< recoNo[i]  <<", Gen Jpsi : "<< genNo[i] <<endl;
        dataFile<<" Eff ["<<i<<"] - "<< eff[i] <<" Error "<< effErr[i] <<endl;
        cout<<" Bin ["<<i<<"] - "<< " Reco Jpsi : "<< recoNo[i]  <<", Gen Jpsi : "<< genNo[i] <<endl;
        cout<<" Eff ["<<i<<"] - "<< eff[i] <<" Error "<< effErr[i] <<endl;
      }
    }
    TH1F *hEff = new TH1F();
    int nEffBins;
    if(iSpec == 0) {hEff = new TH1F("hEff","hEff;p_{T} GeV/c;Efficiency",nPtBins,pt_bound); nEffBins = nPtBins;}
    if(iSpec == 1) {hEff = new TH1F("hEff","hEff;y;Efficiency",nRapBins,rap_bound); nEffBins = nRapBins;}

    dataFile<<"Efficiency"<<endl;
    for(int i = 0; i < nEffBins; i++){
      hEff->SetBinContent(i+1,eff[i]);
      hEff->SetBinError(i+1,effErr[i]);
      cout<<"Trying to measure eff vs "<<cSp[iSpec]<<endl;
      cout<<"Eff : "<<eff[i]<<", err : "<<effErr[i]<<endl;
      dataFile<<eff[i]<<endl;
    }

    dataFile<<""<<endl;
    dataFile<<"Errors"<<endl;
    for(int i = 0; i < nEffBins; i++){
      hEff->SetBinContent(i+1,eff[i]);
      hEff->SetBinError(i+1,effErr[i]);
      cout<<"Trying to measure eff vs "<<cSp[iSpec]<<endl;
      cout<<"Eff : "<<eff[i]<<", err : "<<effErr[i]<<endl;
      dataFile<<effErr[i]<<endl;
    }


    outfile->cd();
    hEff->Draw();
    char tmp_histo[512];
    sprintf(tmp_histo,"hEff_%s",cSp[iSpec]);
    hEff->SetName(tmp_histo);
    hEff->Write();
    dataFile<<""<<endl;
    dataFile.close();
  }
  outfile->Write();
}

//.........这里部分代码省略.........
    //TF1* ffff = new TF1("ffff","expo(0) + gaus(2)",0.4,1.7);

    ring_histo[i]->Sumw2();
    //ring_histo[i]->Rebin(3);

    ringfit[i] = new TF1(name,"pol1(0) + gaus(2)");
    ringfit[i]->SetParLimits(0,0,10.0*ring_histo[i]->GetBinContent(1));
    ringfit[i]->SetParLimits(1,-10000,0);
    ringfit[i]->SetParLimits(2,0,10.0*ring_histo[i]->GetMaximum());
    ringfit[i]->SetParLimits(3,0,10);
    ringfit[i]->SetParameter(0,ring_histo[i]->GetBinContent(1));
    ringfit[i]->SetParameter(1,-ring_histo[i]->GetBinContent(1)/6.0);
    ringfit[i]->SetParameter(2,ring_histo[i]->GetMaximum());
    ringfit[i]->SetParameter(3,0.95);
    ringfit[i]->SetParameter(4,0.15);
    ringfit[i]->SetParNames("constant1","Slope","constant2","Mean","Sigma");
    /*
    ringfit[i] = new TF1(name,fit_function2,0.1,2.5,8);
    ringfit[i]->SetParameter(1,1.);
    ringfit[i]->SetParameter(2,0.2);
    ringfit[i]->SetParameter(3,1.5); //relative height of peak to bg
    ringfit[i]->SetParameter(4,0.25);
    ringfit[i]->SetParameter(5,0.15);
    ringfit[i]->SetParameter(7,0.8);
    ringfit[i]->SetParNames("Constant","Mean","Sigma","Peak Ratio","Bg Mean","Bg Sigma","Bg2 constant","Bg2 decay");
    	
    */
    ringfit[i]->SetLineColor(kBlue);
    ringfit[i]->SetLineWidth(0.6);
    
    ring_histo[i]->Fit(ringfit[i],"rql","",0.2,1.7);
		
    ringprec->SetBinContent(i+1,(ringfit[i]->GetParameter(3)));
    ringprec->SetBinError(i+1,ringfit[i]->GetParameter(4));
    ringprec2->SetBinContent(i+1,(ringfit[i]->GetParameter(3)));
    ringprec2->SetBinError(i+1,ringfit[i]->GetParError(3));
    //ew[i] = 4066/(60*(fit[i]->GetParameter(2))*(fit[i]->GetParameter(2)));

    float mean = ringfit[i]->GetParameter(3);
    float merr = ringfit[i]->GetParError(3);
    cout<<"ring "<<i<<" "<<mean<<" "<<merr/mean<<" "<<ring_histo[i]->GetEntries()<<endl;
  }

  for(int i = 0; i < nrings/2; i++){
    ring2_histo[i]->Add(ring_histo[2*i]);
    ring2_histo[i]->Add(ring_histo[2*i+1]);
    sprintf(name,"ring2_fit_%i",i);
    ring2fit[i] = new TF1(name,"pol1(0) + gaus(2)",0.3,1.7);

    ring2_histo[i]->Rebin(3);

    ring2fit[i]->SetParLimits(0,0,10.0*ring2_histo[i]->GetBinContent(1));
    ring2fit[i]->SetParLimits(1,-10000,0);
    ring2fit[i]->SetParLimits(2,0,10.0*ring2_histo[i]->GetMaximum());
    ring2fit[i]->SetParLimits(3,0,10);
    ring2fit[i]->SetParLimits(4,0.17,0.175);
    ring2fit[i]->SetParameter(0,ring2_histo[i]->GetBinContent(1));
    ring2fit[i]->SetParameter(1,-ring2_histo[i]->GetBinContent(1)/6.0);
    ring2fit[i]->SetParameter(2,ring2_histo[i]->GetMaximum());
    ring2fit[i]->SetParameter(3,0.95);
    ring2fit[i]->SetParameter(4,0.11245);
    ring2fit[i]->SetParNames("constant1","Slope","constant2","Mean","Sigma");

    ring2_histo[i]->Fit(ring2fit[i],"rql","",0.3,1.7);

    ring2prec->SetBinContent(i+1,(ring2fit[i]->GetParameter(3)));

void oDependence()
{
  //set plot options
  gStyle->SetOptFit(1); 
  gStyle->SetOptStat(1);  
  TCanvas  * c1 = new TCanvas("canPRF","Pad response function",700,900);
  TPad * pad1 = new TPad("pad1THR","",0.05,0.55,0.45,0.95,21);
  pad1->Draw();
  TPad * pad2 = new TPad("pad2PRF","",0.55,0.55,0.95,0.95,21);
  pad2->Draw(); 
  TPad * pad3 = new TPad("pad3PRF","",0.55,0.05,0.95,0.45,21);
  pad3->Draw(); 
 
  pad1->cd();
  pad1->SetGridx();
  pad1->SetGridy();
  pad2->SetGridx();
  pad2->SetGridy();
  pad3->SetGridx();
  pad3->SetGridy();

  //make histogram of threshold dependence
  TH1F * hotd =new TH1F("Occupancy dependence on threshold",
			"Ocupancy at first pad row as function of threshold",
                        25,0.,25.);

  hotd->SetBinContent(5,0.625);
  hotd->SetBinError(5,0.02);
  hotd->SetBinContent(10,0.559);
  hotd->SetBinError(10,0.02); 
  hotd->SetBinContent(20,0.478);
  hotd->SetBinError(20,0.02);
  hotd->SetXTitle("Threshold   [channels]");
  hotd->SetYTitle("occupancy");
  hotd->Fit("pol1","+");  
  hotd->Draw("error");
  //make histogram of PRF  dependence
  TH1F * hoprfd =new TH1F("Occupancy dependence on PRF width",
			"Occupancy at first pad row as function of generic PRF sigma for  2.05x0.35 cm pad size ",
                        65, 0.,6.5);
  hoprfd->SetBinContent(10,0.492);
  hoprfd->SetBinError(10,0.02);

  hoprfd->SetBinContent(20,0.524);
  hoprfd->SetBinError(20,0.02); 

  hoprfd->SetBinContent(30,0.559);
  hoprfd->SetBinError(30,0.02);
  hoprfd->SetXTitle("Sigma of PRF   [mm]");
  hoprfd->SetYTitle("occupancy");
  pad2->cd();
  hoprfd->Fit("pol1","+");  
  hoprfd->Draw("error");
  pad2->Draw();
  //pad 3 histogram  
  pad3->cd();
   TH1F * hoprfd88 =new TH1F("Occupancy dependence on PRF width 08x08",
			"Occupancy at first pad row as function of generic PRF sigma for  0.8x0.8 cm pad size ",
                        65, 0.,6.5);

  hoprfd88->SetBinContent(20,0.322);
  hoprfd88->SetBinError(20,0.02);

  hoprfd88->SetBinContent(30,0.344);
  hoprfd88->SetBinError(30,0.02); 

  hoprfd88->SetBinContent(40,0.369);
  hoprfd88->SetBinError(40,0.02);
 
  hoprfd88->SetBinContent(60,0.416);
  hoprfd88->SetBinError(60,0.02);
  hoprfd88->SetXTitle("Sigma of PRF   [mm]");
  hoprfd88->SetYTitle("occupancy");
  hoprfd88->Fit("pol1","+");  
  hoprfd88->Draw("error");
  c1->cd();
  TPaveText * comment = new TPaveText(0.05,0.15,0.45,0.35,"NDC");
  comment->SetTextAlign(12);
  comment->SetFillColor(42);  
  comment->ReadFile("commentdep.txt");  
  comment->Draw();

}

//================================================
void DeltaZVsPos(const Int_t save = 0)
{
  THnSparseF *hn = (THnSparseF*)f->Get(Form("mhTrkDzDy_%s",trigName[kTrigType]));
  TList *list = new TList;

  // dz vs BL
  TH2F *hTrkDzVsBL = (TH2F*)hn->Projection(1,3);
  c = draw2D(hTrkDzVsBL,Form("%s: #Deltaz of matched track-hit pairs",trigName[kTrigType]));
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_vs_BL_%s.pdf",run_type,run_cfg_name.Data(),trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_vs_BL_%s.png",run_type,run_cfg_name.Data(),trigName[kTrigType]));
    }

  list->Clear();
  TString legName[30];
  TH1F *hTrkDzInBL[30];
  Int_t counter = 0;
  for(Int_t i=0; i<30; i++)
    {
      hTrkDzInBL[i] = (TH1F*)hTrkDzVsBL->ProjectionY(Form("hDeltaZ_BL%d",i+1),i+1,i+1);
      if(hTrkDzInBL[i]->GetEntries()>0)
	{
	  legName[counter] = Form("Module %d",i+1);
	  hTrkDzInBL[i]->SetLineColor(color[counter]);
	  list->Add(hTrkDzInBL[i]);
	  counter ++;
	}
    }
  c = drawHistos(list,"TrkDzInBL",Form("%s: #Deltaz of matched track-hit pairs in backleg;#Deltaz (cm)",trigName[kTrigType]),kTRUE,-100,100,kTRUE,0,1.2*hTrkDzInBL[1]->GetMaximum(),kFALSE,kTRUE,legName,kTRUE,"",0.15,0.25,0.2,0.88,kFALSE,0.04,0.04,kFALSE,1,kTRUE,kFALSE);
  TLine *line = GetLine(0,0,0,1.1*hTrkDzInBL[1]->GetMaximum(),1);
  line->Draw();
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_in_BL_%s.pdf",run_type,run_cfg_name.Data(),trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_in_BL_%s.png",run_type,run_cfg_name.Data(),trigName[kTrigType]));
    }

  // dz vs Mod
  TH2F *hTrkDzVsMod = (TH2F*)hn->Projection(1,4);
  c = draw2D(hTrkDzVsMod,Form("%s: #Deltaz of matched track-hit pairs",trigName[kTrigType]));
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_vs_Mod_%s.pdf",run_type,run_cfg_name.Data(),trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_vs_Mod_%s.png",run_type,run_cfg_name.Data(),trigName[kTrigType]));
    }

  TH1F *hMthMod = (TH1F*)hTrkDzVsMod->ProjectionX("hMthMod");
  hMthMod->Sumw2();
  hMthMod->Scale(1./hMthMod->Integral());
  TH2F *hMtdHitMap = (TH2F*)f->Get(Form("mhMtdHitMap_%s",trigName[kTrigType]));
  TH1F *htmp = (TH1F*)hMtdHitMap->ProjectionY("hHitMod_finebin");
  htmp->Rebin(12);
  TH1F *hMtdHitMod = new TH1F(Form("hMtdHitMod_%s",trigName[kTrigType]),"# of MTD hits per module;module",5,1,6);
  for(int i=0; i<hMtdHitMod->GetNbinsX(); i++)
    {
      hMtdHitMod->SetBinContent(i+1,htmp->GetBinContent(i+1));
      hMtdHitMod->SetBinError(i+1,htmp->GetBinError(i+1));
    }
  hMtdHitMod->Scale(1./hMtdHitMod->Integral());
  list->Clear();
  list->Add(hMthMod);
  list->Add(hMtdHitMod);
  TString legName3[2] = {"Matched good hits","All good hits"};
  c = drawHistos(list,"MtdHitMod",Form("%s: MTD hits per module;module;probability",trigName[kTrigType]),kFALSE,0,5,kTRUE,0,0.5,kFALSE,kTRUE,legName3,kTRUE,"",0.15,0.25,0.6,0.88,kTRUE);
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sCompMtdHitMod_%s.pdf",run_type,run_cfg_name.Data(),trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sCompMtdHitMod_%s.png",run_type,run_cfg_name.Data(),trigName[kTrigType]));
    }

  list->Clear();
  TString legName2[5];
  TH1F *hTrkDzInMod[5];
  for(Int_t i=0; i<5; i++)
    {
      hTrkDzInMod[i] = (TH1F*)hTrkDzVsMod->ProjectionY(Form("hDeltaZ_Mod%d",i+1),i+1,i+1);
      legName2[i] = Form("Module %d",i+1);
      list->Add(hTrkDzInMod[i]);
    }
  c = drawHistos(list,"TrkDzInMod",Form("%s: #Deltaz of matched track-hit pairs in module;#Deltaz (cm)",trigName[kTrigType]),kTRUE,-100,100,kTRUE,0,1.2*hTrkDzInMod[3]->GetMaximum(),kFALSE,kTRUE,legName2,kTRUE,"",0.15,0.25,0.6,0.88,kTRUE);
  TLine *line = GetLine(0,0,0,hTrkDzInMod[3]->GetMaximum()*1.05,1);
  line->Draw();
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_in_Mod_%s.pdf",run_type,run_cfg_name.Data(),trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_in_Mod_%s.png",run_type,run_cfg_name.Data(),trigName[kTrigType]));
    }
}

TObjArray* get_mass_v_vars(RooDataSet* ds){
  //mass variable
  RooRealVar m( "Z_M", "Z_M", 60., 120. );
  //other variables needed to make dataset
  int nvars = 4;
  string vars[] = { "mup_eta", "mum_eta", "mup_phi", "mum_phi"};
  double los[]  = { 2.0, 2.0, -TMath::Pi(), -TMath::Pi() };
  double his[]  = { 4.5, 4.5, TMath::Pi(), TMath::Pi()   };
  int binss[]   = { 10, 10, 10, 10 };

  TObjArray* hists = new TObjArray();

  //fit_mass(ds, "", "");

  for (int i = 0 ; i < nvars ; ++i){
    string var = vars[i] ;
    double lo  = los[i]  ;
    double hi  = his[i]  ;
    int bins   = binss[i];
    vector<RooDataSet*> datasets   ;
    for (int i = 0 ; i < bins ; ++i) {
      ostringstream cut;
      double binlo = lo + i*((hi - lo)/bins);
      double binhi = lo + (i+1)*((hi-lo)/bins);
      cut  <<var<<">"<<binlo<<" && "<<var<<"<"<<binhi;
      //add dataset with specific cuts
      datasets.push_back((RooDataSet*)ds->reduce(cut.str().c_str()));
      
    }
    //make mass histogram
    TH1F* hist   = new TH1F(("mass_v_"+var+"_truth").c_str()   , ("mass_v_"+var+"_truth").c_str()   , bins, lo, hi);
    //calculate latex version of variable if possible
    string latexvar;
    if (var.find("eta") != std::string::npos) latexvar = "#eta";
    else if (var.find("phi") != std::string::npos ) latexvar = "#phi";
    else latexvar = var;
    for (int i = 0 ; i < bins ; ++i) {
      ostringstream s1;
      s1<<var<<"_truth_"<<i;
      double binlo = lo + i*((hi - lo)/bins);
      double binhi = lo + (i+1)*((hi-lo)/bins);
      //cut to be presented on plot
      ostringstream label;
      label<<binlo<<" < "<<latexvar<<" < "<<binhi;
      std::pair<double, double> result = fit_mass(datasets.at(i), s1.str(), label.str());
      hist->SetBinContent(i+1, result.first);
      hist->SetBinError(i+1, result.second);
    }
    
    TCanvas* canv = new TCanvas( "canv", "canv", 800.0, 600.0 );
    hist->SetXTitle(latexvar.c_str());
    hist->SetYTitle("Fitted Z0 Mass");
    hist->Draw("e1");
    canv->SaveAs(("mass_v_"+var+"_truth.pdf").c_str());

    hists->Add(hist);
  }

  return hists;

}

//.........这里部分代码省略.........
             hDYtautau = (TH1F*)gDirectory->Get("hDYtautau");
             httbar_mc = (TH1F*)gDirectory->Get("httbar_mc");
             hWlepton = (TH1F*)gDirectory->Get("hWlepton");
             hdiboson = (TH1F*)gDirectory->Get("hdiboson");
           }
          }

          if (USE_OVERFLOWS) nmass_bins = 7;
          for(int i = 0; i < nmass_bins; i++ ) {
            int nRapbins = 24;
            if(i == 5 && !USE_OVERFLOWS) nRapbins = 12;
            if(i == 6 && USE_OVERFLOWS) nRapbins = 12;
            for( int j = 0; j < nRapbins; j++ ) {
                //take data and subtract backgrounds
                double background = 0;
                double error = 0;
                //FIXME this assumes the data-driven error is proper in the file acessed
                if (DATA_DRIVEN_BKG) { 
                   background = (hqcd->GetBinContent(i*24+j+1)+httbar->GetBinContent(i*24+j+1));
                   error = sqrt(pow(hdata->GetBinError(i*24+j+1),2)+ pow(hqcd->GetBinError(i*24+j+1),2)+pow(httbar->GetBinError(i*24+j+1),2)); 
                } else { 
                   background = (hqcd->GetBinContent(i*24+j+1)+httbar_mc->GetBinContent(i*24+j+1)+hDYtautau->GetBinContent(i*24+j+1)+hWlepton->GetBinContent(i*24+j+1)+hdiboson->GetBinContent(i*24+j+1));
                   error = sqrt(pow(hdata->GetBinError(i*24+j+1),2)+ pow(hqcd->GetBinError(i*24+j+1),2)+pow(httbar_mc->GetBinError(i*24+j+1),2)+pow(hDYtautau->GetBinError(i*24+j+1),2)+pow(hWlepton->GetBinError(i*24+j+1),2)+pow(hdiboson->GetBinError(i*24+j+1),2));
                }
                input_data.at(i*24+j) = hdata->GetBinContent(i*24+j+1) - background;
                input_err.at(i*24+j) = error;
            }
          }

  	  int numbins = _inMeasData->GetNbinsX();
	  for( int ibin = 0; ibin < numbins; ibin++ ) {
	    //cout << "input(" << ibin << ") = " << input_data.at(ibin) << " " << input_err.at(ibin) << endl;
	    _inMeasData->SetBinContent(ibin+1, input_data.at(ibin));
	    _inMeasData->SetBinError(ibin+1, input_err.at(ibin));
	  }
	  _inMeasData->SetBinContent(0, 0);
	  _inMeasData->SetBinError(0, 0);
        } ////not closure test

        //
        //SET UP CANVAS TO DRAW
        //
        can = new TCanvas("can","can",800,600);
        gStyle->SetPalette(1);
        //can->SetLogy();
        //can->SetLogx();

        //
        //STEP2: SET UP THE UNFOLDER CLASS
        //   
        _unfolder = new Unfolder("badger");
        _unfolder->setTrue(_inTrue);
        if (!CLOSURE_TEST) {
          _inMeas->Add(_inMeasData);
          _unfolder->setMeas(_inMeas);
        } else {
          _unfolder->setMeas(_inMeasMC);
        }
        _unfolder->setMigrationMatrix(_inResMat);
        _unfolder->unfold();
        _unfolder->getUnfolded()->SetLineColor(kRed);//Draw();  
        _unfolder->getTrue()->Draw("hist");
        //_unfolder->getMeas()->Draw("hist");

        //_inMeas->Print("all");        

//.........这里部分代码省略.........
       pow(histo1D_ttSig->GetBinContent(ilum+1)*SF_trigger_error, 2)+
       //*************************
       //uncertinty on lepton sel
       pow(histo1D_ttSig->GetBinContent(ilum+1)*SF_Lepton_error, 2)+
       //*************************
       //uncertinty on met sel
       pow( histo1D_ttSig->GetBinContent(ilum+1)*SF_MET_error, 2);
     
       error_all += pow(histo1D_ttSig->GetBinContent(ilum+1)*systeError_ttbar, 2);   
       if(ilum > 2) error_all += pow(histo1D_ttSig->GetBinContent(ilum+1)*0.01, 2);    
       if(ilum > 6) error_all += pow(histo1D_ttSig->GetBinContent(ilum+1)*0.02, 2);   
       error_all += pow(histo1D_VV->GetBinContent(ilum+1)*systeError_VV, 2);
       error_all += pow(histo1D_tW->GetBinContent(ilum+1)*systeError_stop, 2);
     
        
       //*************************
       //uncertinty on met sel
       if(ilum < 5) pow( histo1D_ttSig->GetBinContent(ilum+1)*0.02, 2);
       if(ilum >= 5 && ilum < 8) pow( histo1D_ttSig->GetBinContent(ilum+1)*0.03, 2);
       if(ilum >= 8)             pow( histo1D_ttSig->GetBinContent(ilum+1)*0.04, 2);
	  
	
	    
	    //frac_ee   = frac_ee/  (frac_ee+frac_mumu+frac_emu);
	    //frac_mumu = frac_mumu/(frac_ee+frac_mumu+frac_emu);
	    //frac_emu  = frac_emu/ (frac_ee+frac_mumu+frac_emu);
	    
	    //cout << "frac_ee   " << frac_ee << endl;
	    //cout << "frac_mumu " << frac_mumu << endl;
	    //cout << "frac_emu  " << frac_emu << endl;
	    
	   error_all += pow(histo1D_DY->GetBinContent(ilum+1)*((0.40/1.9281)*frac_ee + (0.37/1.82219)*frac_mumu+  (0.30/1.38872)*frac_emu), 2);
	   
     lumiband->SetBinError(ilum+1,sqrt(error_all));
     
     //modifications
     histo1D_mc->SetBinError(ilum+1,sqrt(error_all));
  }
  
  
  TGraphErrors *thegraph = new TGraphErrors(lumiband);
  thegraph->SetFillStyle(3005);
  thegraph->SetFillColor(1);
  
  histo1D_ttSig->SetFillStyle(1001);
  histo1D_ttBkg->SetFillStyle(1001);
  histo1D_DY->SetFillStyle(1001);
  histo1D_VV->SetFillStyle(1001);
  histo1D_tW->SetFillStyle(1001);
    
  histo1D_ttSig->SetFillColor(kRed+1);
  histo1D_ttBkg->SetFillColor(kRed-7);
  histo1D_DY->SetFillColor(kAzure-2);
  histo1D_VV->SetFillColor(13);
  histo1D_tW->SetFillColor(kMagenta);
  
  
  histo1D_ttSig->GetYaxis()->CenterTitle();
  histo1D_ttSig->GetYaxis()->SetTitle("");
  histo1D_ttSig->GetXaxis()->SetLabelSize(0);
  histo1D_ttSig->GetXaxis()->SetTitleSize(0);
  
  histo1D_ttBkg->GetYaxis()->CenterTitle();
  histo1D_ttBkg->GetYaxis()->SetTitle("");
  histo1D_ttBkg->GetXaxis()->SetLabelSize(0);
  histo1D_ttBkg->GetXaxis()->SetTitleSize(0);

void AnalysisBase::findCutoutRegion(TH2F* h){

  cout << "===============================================================" << endl;
  cout << "findCutoutRegion(): Looking for cutout region in D type sensor " << endl;
  cout << "===============================================================" << endl;


  int ilx = h->GetXaxis()->FindBin(xMin) + 1;
  int ihx = h->GetXaxis()->FindBin(xMax) - 1;
  int ily = 1;//h->GetYaxis()->GetBinFindBin(yMin) - 10;
  int ihy = h->GetYaxis()->FindBin(yMax) + 10;
  if(ily < 0) ily = 0;
  if(ihy > h->GetYaxis()->GetNbins()) ihy = h->GetYaxis()->GetNbins();  
  int nb = h->GetXaxis()->GetNbins(); 

  double xlow = h->GetXaxis()->GetBinLowEdge(1);
  double xhi = h->GetXaxis()->GetBinLowEdge(nb)+h->GetXaxis()->GetBinWidth(1);

  TH1F *hpr = new TH1F("hpr","Y vs X, edge",nb,xlow,xhi);

  TH1D *hpy;
  int ipeak = 0;
  for(int i=ilx; i<ihx-1; i++){
    hpy = h->ProjectionY("hpy",i,i);
    double maxcon  = 0.0;
    double maxbin  = 0.0;
    for(int j=ily;j<ihy;j++){
      maxcon = maxcon + hpy->GetBinContent(j);
      if(hpy->GetBinContent(j) > maxbin) {
        maxbin = hpy->GetBinContent(j) ;
        ipeak = j;
      }
    }
    // Find the lower "edge"
    for(int j=ipeak; j>=ily; j--){
      double r0 = (hpy->GetBinContent(j-1)+hpy->GetBinContent(j-2))/2.0;
      double r1 = hpy->GetBinContent(j);
      double r2 = (hpy->GetBinContent(j+3)+hpy->GetBinContent(j+4))/2.0;
      if(r1<=2 && r0<=1.5 && r2>5*r1 && r2>8){
        hpr->SetBinContent(i,hpy->GetBinCenter(j)+2*hpy->GetBinWidth(j));
        hpr->SetBinError(i,hpy->GetBinWidth(1)/2.0);
        //cout << "found: " << i << " " << hpy->GetEntries() << " " << hpy->GetBinCenter(j) << " " 
        //     << hpy->GetBinContent(j) << " " << r0 << " " << r1 << " " << r2 << " " << endl;
        break;
      }
    }
    
  }

  TF1* poly2 = new TF1("poly2","[0]+[1]*x+[2]*x*x",xMin,xMax);
  poly2->SetParameters(-1.5,-0.17,-0.15);
  hpr->Fit(poly2,"R0");
  hpr->SetLineColor(kRed);

  holeQuadPar[0] = poly2->GetParameter(0);
  holeQuadPar[1] = poly2->GetParameter(1);
  holeQuadPar[2] = poly2->GetParameter(2);

  //return;
  

  delete poly2;
  delete hpr;
  delete hpy;
  
  return;
}

//.........这里部分代码省略.........
	}

      }
    }



    int HV=0;

    for(int iDepth = MinDepth; iDepth <= MaxDepth; iDepth++){
      for(int iPhi = MinPhi; iPhi <= MaxPhi; iPhi++){
	for(int iEta = MinEta; iEta <= MaxEta; iEta++){

	  //cout<<iDepth<<" "<<iPhi<<" "<<iEta<<endl;

	  if(!LED[iDepth][iEta][iPhi]) continue;

	  sprintf(spehistname,"led %d %d %d",iDepth,iEta,iPhi);
	  TH1F *hspe_temp = (TH1F *)LED[iDepth][iEta][iPhi]->Clone(spehistname);
	  sprintf(spehistname,"ped %d %d %d",iDepth,iEta,iPhi);
	  TH1F *hped = (TH1F *)PED[iDepth][iEta][iPhi]->Clone(spehistname);
	  hspe->Reset();
	  sprintf (spehistname, "SumLED_Depth_%d_Eta_%d_Phi_%d",iDepth,iEta,iPhi);
	  hspe->SetTitle(spehistname);

	  //combine bins of original SPE histogram
	  for(int ib=1; ib<=hspe_temp->GetNbinsX(); ib++) {
	    double bin_center = hspe_temp->GetBinCenter(ib);
	    if(bin_center>hspe->GetXaxis()->GetXmax()) continue;
	    int newbin = hspe->FindBin(bin_center);
	    double new_content = hspe->GetBinContent(newbin) + hspe_temp->GetBinContent(ib);
	    double new_error = sqrt(pow(hspe->GetBinError(newbin),2)+pow(hspe_temp->GetBinError(ib),2));
	    hspe->SetBinContent(newbin,new_content);
	    hspe->SetBinError(newbin,new_error);
	  }
	  TH1F* hspe_unscaled = (TH1F*)hspe->Clone("hspe_unscaled");
	  //renormalize bins of new SPE histogram
	  for(int ib=1; ib<=hspe->GetNbinsX(); ib++) {
	    double new_content = hspe->GetBinContent(ib)/hspe->GetXaxis()->GetBinWidth(ib)*hspe_temp->GetXaxis()->GetBinWidth(1);
	    double new_error = hspe->GetBinError(ib)/hspe->GetXaxis()->GetBinWidth(ib)*hspe_temp->GetXaxis()->GetBinWidth(1);
	    hspe->SetBinContent(ib,new_content);
	    hspe->SetBinError(ib,new_error);
	  }
	  
	  if(hspe_temp->Integral()==0) continue;
	  else drawflag[iDepth][iPhi] = true;	  

	  Nev = hspe_temp->Integral()*hspe_temp->GetXaxis()->GetBinWidth(1); 
	  //TF1 *fped = new TF1("fped","gaus",0, 80);
	  TF1 *fped = new TF1("fped","gaus",0, 20);
	  hped->Fit(fped,"NQR");
	  double pploc = fped->GetParameter(1), ppwidth = fped->GetParameter(2);
	  cout<<"depth "<<iDepth<<" ieta "<<iEta<<" iphi "<<iPhi
	      <<" pploc "<<pploc<<" ppwidth "<<ppwidth<<endl;
	  //hspe->Fit(fped, "NQ", "", pploc - 3*ppwidth, pploc + ppwidth);  
	  hspe->Fit(fped, "NQ", "", pploc - 3*ppwidth, pploc + 3*ppwidth);  
	  
	  //estimate SPE peak location
	  int max_SPE_bin, maxbin, Nbins;
	  double max_SPE_height=0, minheight, max_SPE_location;
	  bool minflag = false;
	  maxbin=hspe->FindBin(fped->GetParameter(1)); //location of pedestal peak
	  minheight=hspe->GetBinContent(maxbin); //initialize minheight
	  Nbins = hspe->GetNbinsX();
	  for(int j=maxbin+1; j<Nbins-1; j++) { //start from pedestal peak and loop through bins
	    if(hspe->GetBinContent(j) > minheight && !minflag) minflag=true; //only look for SPE peak when minflag=true

//.........这里部分代码省略.........
         GenTemp->Draw("same");

         TPaveText *label2 = util::LabelFactory::createPaveTextWithOffset(2,1.0,0.05);
         label2->AddText("Anti-k_{T} (R=0.5) PFchs Jets");
         label2->AddText( Form("%0.1f #leq |#eta| #leq %0.1f, %3.0f #leq #bar{ p}_{T} [GeV] #leq %3.0f", eta_bins[ieta], eta_bins[ieta+1], pt_bins[ipt], pt_bins[ipt+1]) );
         label2->Draw("same");
  
         TLegend* leg2 = util::LabelFactory::createLegendWithOffset(2,0.15);
         leg2->AddEntry(extrapol_Gen,"Extrapolation (PLI)","LP");

         leg2->Draw();
         cmsPrel(-1, false , 8);

         TString name2;
         name2 = Form("ClosureTest/Extrapol_Eta%i_pt%i_gen" + suffix + ".eps", ieta, ipt);
         cb->Print(name2);

         float par_data = lin_extrapol_data->GetParameter(0);
         float par_data_err = lin_extrapol_data->GetParError(0);
         float par_mc = lin_extrapol_mc->GetParameter(0);
         float par_mc_err = lin_extrapol_mc->GetParError(0);
         float par_gen = lin_extrapol_gen->GetParameter(0);
         float par_gen_err = lin_extrapol_gen->GetParError(0);
       
         cout << "ieta : " << ieta << "  ipt : " << ipt << endl;
         cout << "Parameter data: " << par_data << endl;
         cout << "Parameter error data: " << par_data_err << endl;
         cout << "Parameter mc: " << par_mc << endl;
         cout << "Parameter error mc: " << par_mc_err << endl;
         cout << "Parameter gen: " << par_gen << endl;
         cout << "Parameter error gen: " << par_gen_err << endl;

         extrapolated_mcsmeared->SetBinContent(ipt+1, par_data);
         extrapolated_mcsmeared->SetBinError(ipt+1, par_data_err);
         extrapolated_mc->SetBinContent(ipt+1, par_mc);
         extrapolated_mc->SetBinError(ipt+1, par_mc_err);   
         extrapolated_gen->SetBinContent(ipt+1, par_gen);
         extrapolated_gen->SetBinError(ipt+1, par_gen_err);    

         float par_data_pli_corr = 0;
         float par_data_pli_corr_err = 0;
         float par_mc_pli_corr = 0;
         float par_mc_pli_corr_err = 0;

         if(par_gen > 0 && par_data > 0 && par_mc > 0 && par_data > par_gen && par_mc > par_gen) {
            par_data_pli_corr = TMath::Sqrt(pow(par_data,2) - pow(par_gen,2));
            par_data_pli_corr_err = TMath::Sqrt( pow(par_data,2)/(pow(par_data,2) - pow(par_gen,2)) * pow(par_data_err,2) +  pow(par_gen,2)/(pow(par_data,2) - pow(par_gen,2)) * pow(par_gen_err,2));
            par_mc_pli_corr = TMath::Sqrt(pow(par_mc,2) - pow(par_gen,2));
            par_mc_pli_corr_err = TMath::Sqrt( pow(par_mc,2)/(pow(par_mc,2) - pow(par_gen,2)) * pow(par_mc_err,2) +  pow(par_gen,2)/(pow(par_mc,2) - pow(par_gen,2)) * pow(par_gen_err,2));
         }
            
         extrapolated_mcsmeared_with_pli->SetBinContent(ipt+1, par_data_pli_corr);
         extrapolated_mcsmeared_with_pli->SetBinError(ipt+1, par_data_pli_corr_err);
         extrapolated_mc_with_pli->SetBinContent(ipt+1, par_mc_pli_corr);
         extrapolated_mc_with_pli->SetBinError(ipt+1, par_mc_pli_corr_err);

         cout << "Parameter data after pli: " << par_data_pli_corr << endl;
         cout << "Parameter error data after pli: " << par_data_pli_corr_err << endl;
         cout << "Parameter mc after pli: " << par_mc_pli_corr << endl;
         cout << "Parameter error mc after pli: " << par_mc_pli_corr_err << endl;
      }

      // --------------------------------------- //
      // calc data/mc ratio and fit with constant
      TH1F* ratio = new TH1F(*extrapolated_mc);
      TH1F* ratio_with_pli = new TH1F(*extrapolated_mc);

void drawDifference(TH1* iH0,TH1 *iH1,TH1* iH2, TGraphErrors* iH3, int chnl,TGraphErrors* iH4,TGraphAsymmErrors* iH5,TH1* StatErrBand,TGraphErrors* iH6){
  std::string lName = std::string(iH0->GetName());
  TH1F *lHDiff  = new TH1F((lName+"Diff").c_str(),(lName+"Diff").c_str(),nBins-1,WptLogBins);// lHDiff->Sumw2();
  TH1F *lXHDiff1 = new TH1F((lName+"XDiff1").c_str(),(lName+"XDiff1").c_str(),iH0->GetNbinsX(),iH0->GetXaxis()->GetXmin(),iH0->GetXaxis()->GetXmax());
  int i1 = 0;
  lXHDiff1->SetLineWidth(2); lXHDiff1->SetLineColor(kBlack); //lXHDiff1->SetLineStyle(2);
  
  StatErrBand->SetMarkerStyle(kFullCircle); StatErrBand->SetMarkerColor(kBlack);StatErrBand->SetMarkerSize(0.6);
  StatErrBand->SetLineWidth(2); StatErrBand->SetLineColor(kBlack);

  //lHDiff->GetYaxis()->SetRangeUser(0.2,1.8);
  lHDiff->GetYaxis()->SetRangeUser(0.6,1.5);
  if (chnl == 2)
    lHDiff->GetYaxis()->SetRangeUser(0.4,1.4);
  if (chnl == 3)
    lHDiff->GetYaxis()->SetRangeUser(0.3,1.3);
  lHDiff->GetYaxis()->SetTitleOffset(0.4);
  lHDiff->GetYaxis()->SetTitleSize(0.12);
  lHDiff->GetYaxis()->SetLabelSize(0.12);
  lHDiff->GetYaxis()->CenterTitle();
  lHDiff->GetXaxis()->SetTitleOffset(1.);
  lHDiff->GetXaxis()->SetTitleSize(0.12);
  lHDiff->GetXaxis()->SetLabelSize(0);
  if(chnl==3)
    lHDiff->GetXaxis()->SetLabelSize(0.12);
  lHDiff->GetYaxis()->SetNdivisions(405);
  if (chnl == 3)
    lHDiff->GetXaxis()->SetTitle(" W p_{T} ");
  lHDiff->GetYaxis()->SetTitle("Theory/Data");
  //lHDiff->GetYaxis()->SetTitle("Data/ResBos");
  gStyle->SetOptStat(0);
  
  for(int i0 = 0; i0 < lHDiff->GetNbinsX()+1; i0++) {
    double lXCenter = lHDiff->GetBinCenter(i0);
    double lXVal     = iH0   ->GetBinContent(i0);
    lXHDiff1->SetBinContent(i0, 1.0);
    while(iH1->GetBinCenter(i1) < lXCenter) {i1++;}
    if(iH1->GetBinContent(i0) > 0) lHDiff->SetBinContent(i0,lXVal/(iH1->GetBinContent(i0)));
    //if(iH1->GetBinContent(i0) > 0) lHDiff->SetBinError(i0,iH0->GetBinError(i0)/(iH1->GetBinContent(i0)));
    if(iH1->GetBinContent(i0) > 0) lHDiff->SetBinError(i0,0.00001);
  }
  
  TGraphErrors* ErrBand = new TGraphErrors(iH2);
  ErrBand->SetFillColor(kBlack);
  ErrBand->SetFillStyle(3354);
  ErrBand->SetLineWidth(1);
  
  if (chnl == 1)
  {
    lHDiff->SetMarkerStyle(kOpenCircle);
    lHDiff->SetMarkerColor(kBlue);
    lHDiff->SetLineColor(kBlue);
  }
  if (chnl == 2)
  {
    lHDiff->SetMarkerStyle(kOpenTriangleUp);
    lHDiff->SetMarkerColor(kRed);
    lHDiff->SetLineColor(kRed);
  }
  if (chnl == 3)
  {
    lHDiff->SetMarkerStyle(kOpenSquare);
    lHDiff->SetMarkerColor(kGreen+3);
    lHDiff->SetLineColor(kGreen+3);
  }
  
  lHDiff->SetMarkerSize(0.8);
  
  //lHDiff->SetLineColor(kBlack); lHDiff->SetMarkerColor(kBlack);
  
  lHDiff->SetTitle("");
  lHDiff->Draw("E");
  if (chnl == 1)
    iH5->Draw("2same");
  if (chnl == 2 || chnl == 3)
    iH4->Draw("2same");
  if (chnl == 2 || chnl == 3)
    iH3->Draw("2same");
  if (chnl == 3)
    iH6->Draw("2");
  lXHDiff1->Draw("histsame");
  ErrBand->Draw("2same");
  lHDiff->Draw("Esame");
  StatErrBand->Draw("E1same");
}