C++ TGraph::Eval方法代码示例

vector<double> getValsFromLikelihood(TGraph *graph){
 
  TGraph *grRot = new TGraph();
  TGraph *grRotUpp = new TGraph();
  TGraph *grRotLow = new TGraph();
  int pLow=0,pHigh=0;
  pair<double,double> graphMin = getGraphMin(graph);
  for (int p=0; p<graph->GetN(); p++){
    double x,y;
    graph->GetPoint(p,x,y);
    grRot->SetPoint(p,y,x);
    if (x<=graphMin.first) {
      grRotLow->SetPoint(pLow,y,x);
      pLow++;
    }
    if (x>=graphMin.first) {
      grRotUpp->SetPoint(pHigh,y,x);
      pHigh++;
    }
  }
  // return best fit and +/- 1/2 sigma errors
  vector<double> result;
  result.push_back(grRot->Eval(0.));
  result.push_back(grRotLow->Eval(1.0));
  result.push_back(grRotUpp->Eval(1.0));
  result.push_back(grRotLow->Eval(4.0));
  result.push_back(grRotUpp->Eval(4.0));
  return result; 
}

void points(TString filename) {

	
	TString cmssw;
	// 167
	
	cmssw = "$1.6.7$";
	
	TFile *f = TFile::Open(filename);

	std::vector< TString > taggers;
	taggers.push_back( "gTC2_udsg" );
	taggers.push_back( "gTC3_udsg" );
	taggers.push_back( "gTP_udsg" );
	taggers.push_back( "gJBP_udsg" );
 	taggers.push_back( "gSSV_udsg" );
 	taggers.push_back( "gCSV_udsg" );

	std::vector< TString > discriminators;
	discriminators.push_back( "discTC2_udsg" );
	discriminators.push_back( "discTC3_udsg" );
	discriminators.push_back( "discTP_udsg" );
	discriminators.push_back( "discJBP_udsg" );
 	discriminators.push_back( "discSSV_udsg" );
 	discriminators.push_back( "discCSV_udsg" );

	//TCanvas *cv_TC = new TCanvas("cv_TC","cv_TC",700,700);
	//TCanvas *cv_TP = new TCanvas("cv_TP","cv_TP",700,700);

	std::cout << "Tagger & Point & Discriminator & light mistag & b-efficiency \\\\ \\hline" << std::endl;
	for ( size_t itagger = 0; itagger < taggers.size(); ++itagger ) {

		TString tag = taggers[itagger];
		TGraphErrors *agraph = (TGraphErrors*) gDirectory->Get("Histograms/MCtruth/"+tag);

		//if (taggers == "gTC2_udsg" || taggers =
		
		TGraph *dgraph = (TGraph*) gDirectory->Get("Histograms/MCtruth/"+discriminators[itagger]);
		dgraph->Sort();

		TGraphErrors *g = new TGraphErrors(agraph->GetN(),agraph->GetY(),agraph->GetX(),agraph->GetEY(),agraph->GetEX());
		g->Sort();
		
		//cv[itagger] = new TCanvas("cv","cv",600,600);
		//g->Draw("ACP");
		TString se = " & ";
		
		
		std::cout << tag << se << "Loose" << se << std::setprecision(3) << dgraph->Eval(0.1) << se << "0.1" << se << std::setprecision(2) << g->Eval(0.1) << "\\\\" << std::endl;
		std::cout << tag << se << "Medium" << se << std::setprecision(3) << dgraph->Eval(0.01) << se << "0.01" << se << std::setprecision(2) << g->Eval(0.01) << "\\\\" << std::endl;
		std::cout << tag << se << "Tight" << se << std::setprecision(3) << dgraph->Eval(0.001) << se << "0.001" << se << std::setprecision(2) << g->Eval(0.001) << "\\\\ \\hline" << std::endl;
		
	}

	
}

  /** 
   * Create ratios to other data 
   * 
   * @param ib      Bin number  
   * @param res     Result
   * @param alice   ALICE result if any
   * @param cms     CMS result if any
   * @param all     Stack to add ratio to 
   */
  void Ratio2Stack(Int_t ib, TH1* res, TGraph* alice, TGraph* cms, THStack* all)
  {
    if (!all || !res || !(alice || cms)) return;

    Int_t        off  = 5*ib;
    TGraph*      gs[] = { (alice ? alice : cms), (alice ? cms : 0), 0 };
    TGraph**     pg   = gs;
    while (*pg) { 
      TGraph*     g = *pg;
      const char* n = (g == alice ? "ALICE" : "CMS");

      TH1*    r = static_cast<TH1*>(res->Clone(Form("ratio%s", n)));
      TString tit(r->GetTitle());
      tit.ReplaceAll("Corrected", Form("Ratio to %s", n));
      r->SetTitle(tit);
      r->SetMarkerColor(g->GetMarkerColor());
      r->SetLineColor(g->GetLineColor());

      TObject* tst = r->FindObject("legend");
      if (tst) r->GetListOfFunctions()->Remove(tst);

      for (Int_t i = 1; i <= r->GetNbinsX(); i++) {
	Double_t c = r->GetBinContent(i);
	Double_t e = r->GetBinError(i);
	Double_t o = g->Eval(r->GetBinCenter(i));
	if (o < 1e-12) { 
	  r->SetBinContent(i, 0);
	  r->SetBinError(i, 0);
	  continue;
	}
	r->SetBinContent(i, (c - o) / o + off);
	r->SetBinError(i, e / o);
      }
      all->Add(r);
      pg++;
    }
    TLegend* leg = StackLegend(all);
    if (!leg) return;
      
    TString   txt      = res->GetTitle();
    txt.ReplaceAll("Corrected P(#it{N}_{ch}) in ", "");
    if      (ib == 0) txt.Append(" "); // (#times1)");
    // else if (ib == 1) txt.Append(" (#times10)");
    else              txt.Append(Form(" (+%d)", off));

    TObject* dummy = 0;
    TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
    e->SetMarkerStyle(res->GetMarkerStyle());
    e->SetMarkerSize(res->GetMarkerSize());
    e->SetMarkerColor(kBlack);
    e->SetFillColor(0);
    e->SetFillStyle(0);
    e->SetLineColor(kBlack);
  }

void disceff(TString filename) {
	gROOT->SetStyle("Plain");
	
	TString cmssw;
	// 167
	
	cmssw = "$3.1.0_pre9$";
		
	TFile *f = TFile::Open(filename);

	std::vector< TString > taggers;
	taggers.push_back( "TC2" );
	taggers.push_back( "TC3" );
	taggers.push_back( "TP" );
	taggers.push_back( "BTP" );
	taggers.push_back( "SSV" );
	taggers.push_back( "CSV" );
	taggers.push_back( "MSV" );
	taggers.push_back( "SMT" );
	taggers.push_back( "SETbyIP3d" );
	taggers.push_back( "SETbyPt" );
	taggers.push_back( "SMTbyIP3d" );
	taggers.push_back( "SMTbyPt" );

	std::vector< TString > discriminators;
	for ( size_t itagger = 0; itagger < taggers.size(); ++itagger ) {
		discriminators.push_back( "disc"+taggers[itagger]+"_udsg" );
	}
//	discriminators.push_back( "discTC3_udsg" );
//	discriminators.push_back( "discTP_udsg" );

	const int dim=taggers.size();
        TCanvas *cv[dim];
	TMultiGraph* mg[dim];
	for ( size_t itagger = 0; itagger < taggers.size(); ++itagger ) {
		TString tag = "g"+taggers[itagger]+"_udsg";
		TString tagb = "g"+taggers[itagger]+"_b";
		TString tagc = "g"+taggers[itagger]+"_c";
		cv[itagger] = new TCanvas("cv_"+taggers[itagger],"cv_"+taggers[itagger],700,700);
		TLegend *legend0 = new TLegend(0.68,0.70,0.88,0.90);

		TGraphErrors *agraph = (TGraphErrors*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+tag);
		TGraphErrors *bgraph = (TGraphErrors*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+tagb);
		TGraphErrors *cgraph = (TGraphErrors*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+tagc);

		TGraph *dgraph = (TGraph*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+discriminators[itagger]);
		TGraph *bvsdgraph = new TGraph(dgraph->GetN(),dgraph->GetY(),bgraph->GetY());
		TGraph *cvsdgraph = new TGraph(dgraph->GetN(),dgraph->GetY(),cgraph->GetY());
		TGraph *lvsdgraph = new TGraph(dgraph->GetN(),dgraph->GetY(),agraph->GetY());
		TGraph *udsgvsdgraph = new TGraph(dgraph->GetN(),dgraph->GetY(),dgraph->GetX());
		dgraph->Sort();
//		udsgvsdgraph->Sort();
//		udsgvsdgraph->SetLineColor(1);
//              legend0 -> AddEntry(udsgvsdgraph,"control","l");
		lvsdgraph->Sort();
		lvsdgraph->SetLineColor(2);
                legend0 -> AddEntry(lvsdgraph,tag,"l");
		cvsdgraph->Sort();
		cvsdgraph->SetLineColor(3);
                legend0 -> AddEntry(cvsdgraph,tagc,"l");
		bvsdgraph->Sort();
		bvsdgraph->SetLineColor(4);
                legend0 -> AddEntry(bvsdgraph,tagb,"l");
		mg[itagger]= new TMultiGraph();
//		mg[itagger]->Add(udsgvsdgraph);
		mg[itagger]->Add(lvsdgraph);
		mg[itagger]->Add(cvsdgraph);
		mg[itagger]->Add(bvsdgraph);
//		mg[itagger]->Add(dgraph);

		cv[itagger]->cd(1);
		mg[itagger]->Draw("ALP");
		mg[itagger]->GetYaxis()->SetTitle("eff");
		mg[itagger]->GetXaxis()->SetTitle("discriminant");
		legend0 -> Draw();
		cv[itagger]->Update();
	        cv[itagger]->cd(0);
		cv[itagger]-> Print ("BTagPATeff_vs_disc"+taggers[itagger]+".eps");
		cv[itagger]-> Print ("BTagPATeff_vs_disc"+taggers[itagger]+".ps");

		TGraphErrors *g = new TGraphErrors(agraph->GetN(),agraph->GetY(),agraph->GetX(),agraph->GetEY(),agraph->GetEX());
		g->Sort();
		g->SetLineColor(itagger+1);

		std::cout << " Tagger: " << tag << std::endl;
		std::cout << " Loose(10%): " << " cut > " << std::setprecision(4) << dgraph->Eval(0.1) << " b-eff = " << g->Eval(0.1) << std::endl;
		std::cout << " Medium(1%):  " << " cut > " << std::setprecision(4) << dgraph->Eval(0.01) << " b-eff = " << g->Eval(0.01) << std::endl;
		std::cout << " Tight(0.1%) = " << " cut > " << std::setprecision(4) << dgraph->Eval(0.001) << " b-eff = " << g->Eval(0.001) << std::endl;

	}//end for
 
	
}

//.........这里部分代码省略.........
    h_scalar_min->GetXaxis()->SetTitleOffset(1.03);
    h_scalar_min->GetXaxis()->SetTitle("DM mass [GeV]");
    h_scalar_min->GetYaxis()->SetTitle("DM-nucleon cross section [cm^{2}]");
    h_scalar_lat->SetLineColor(4);
    h_scalar_lat->SetLineStyle(1);
    h_scalar_lat->SetLineWidth(3);
    h_scalar_max->SetLineColor(4);
    h_scalar_max->SetLineStyle(2);
    h_scalar_max->SetLineWidth(3);


    h_fermion_min->SetLineColor(kRed);
    h_fermion_min->SetLineStyle(2);
    h_fermion_min->SetLineWidth(3);
    h_fermion_lat->SetLineColor(kRed);
    h_fermion_lat->SetLineStyle(1);
    h_fermion_lat->SetLineWidth(3);
    h_fermion_max->SetLineColor(kRed);
    h_fermion_max->SetLineStyle(2);
    h_fermion_max->SetLineWidth(3);


    //h_scalar_lat->SetFillStyle(3005);
    //h_scalar_lat->SetLineWidth(-402);

    h_scalar_min->Draw("AL");
    h_scalar_lat->Draw("L");
    h_scalar_max->Draw("L");

    h_fermion_min->Draw("L");
    h_fermion_lat->Draw("L");
    h_fermion_max->Draw("L");
    z12016->Draw("L");
    //z12015->Draw("L");
    //CRESST2->Draw("C");
    SCDMS->Draw("C");
    PANDAX->Draw("C");

    //leg2->Draw();
    TLatex *lat = new TLatex(); lat->SetTextSize(0.025);
    lat->SetTextFont(42);
    lat->SetTextColor(kGreen+2);
    lat->SetTextAngle(15);
    //lat->DrawLatex(130,z1->Eval(130)*1.5,"LUX #it{Phys. Rev. Lett.} #bf{116} (2016)");


    //lat->DrawLatex(130,z12015->Eval(130)*1.5,"LUX (2015)");
    lat->DrawLatex(130,z12016->Eval(130)*0.5,"LUX (2013+2014-16)");
    
    lat->SetTextColor(SCDMS->GetLineColor());
    lat->SetTextAngle(344);
    //lat->SetFillColor(kWhite);
    lat->DrawLatex(5,SCDMS->Eval(5)*1.5,"CDMSlite (2015)");
    
    lat->SetTextColor(kBlue);
    lat->SetTextAngle(330);
    lat->DrawLatex(13,h_scalar_lat->Eval(13)*1.5,"Scalar DM");

    lat->SetTextColor(kRed);
    lat->SetTextAngle(4);
    lat->DrawLatex(5,h_fermion_lat->Eval(5)*1.5,"Fermion DM");


    lat->SetTextColor(kMagenta+2);
    lat->SetTextAngle(15);
    lat->DrawLatex(130,PANDAX->Eval(130)*1.25,"PandaX-II (2016)");

	TLatex * tex = new TLatex();
	tex->SetNDC();
	tex->SetTextFont(42);
	tex->SetLineWidth(2);
	tex->SetTextSize(0.04);
	tex->SetTextAlign(31);
	tex->DrawLatex(0.93,0.78,"4.9 fb^{-1} (7 TeV) + 19.7 fb^{-1} (8 TeV)");
	tex->DrawLatex(0.93,0.74,"+ 2.3 fb^{-1} (13 TeV)");
	//tex->SetTextAlign(11);
  	tex->SetTextFont(42);
	tex->SetTextSize(0.06);
        TLegend *legBOX; 
	if (isPrelim)	{
   		legBOX = new TLegend(0.16,0.83,0.38,0.89);
		legBOX->SetFillColor(kWhite);
		legBOX->SetLineWidth(0);
		legBOX->Draw();
		tex->DrawLatex(0.17, 0.84, "#bf{CMS} #it{Preliminary}");
	} else {
   		legBOX = new TLegend(0.16,0.83,0.28,0.89);
		legBOX->SetFillColor(kWhite);
		legBOX->SetLineWidth(0);
		//legBOX->Draw();
		tex->DrawLatex(0.93, 0.84, "#bf{CMS}");
	}
  tex->SetTextSize(0.045);
  tex->DrawLatex(0.93,0.68,Form("B(H #rightarrow inv.) < %.2f",BRinv));
  tex->DrawLatex(0.93,0.5,"90% CL limits");
  canv->SetTicky(1);
  canv->SetTickx(1);
  canv->SaveAs("limitsDM.pdf"); 

}

//.........这里部分代码省略.........
   graph->SetPoint(182,94.90000153,1.380649686);
   graph->SetPoint(183,94.94000244,1.434825897);
   graph->SetPoint(184,94.98000336,1.490056992);
   graph->SetPoint(185,95.01999664,1.546342969);
   graph->SetPoint(186,95.05999756,1.603683949);
   graph->SetPoint(187,95.09999847,1.66207993);
   graph->SetPoint(188,95.13999939,1.721530676);
   graph->SetPoint(189,95.18000031,1.782036185);
   graph->SetPoint(190,95.22000122,1.843595982);
   graph->SetPoint(191,95.22000122,1.843595982);
   graph->SetPoint(192,95.26000214,1.906209826);
   graph->SetPoint(193,95.30000305,1.969877124);
   graph->SetPoint(194,95.33999634,2.034597158);
   graph->SetPoint(195,95.37999725,2.100369453);
   graph->SetPoint(196,95.41999817,2.167192936);
   graph->SetPoint(197,95.45999908,2.235066652);
   graph->SetPoint(198,95.5,2.30398941);
   graph->SetPoint(199,95.54000092,2.37395978);
   graph->SetPoint(200,95.58000183,2.44497633);
   graph->SetPoint(201,95.62000275,2.517037392);
   graph->SetPoint(202,95.66000366,2.590141296);
   graph->SetPoint(203,95.69999695,2.66428566);
   graph->SetPoint(204,95.73999786,2.739468575);
   graph->SetPoint(205,95.77999878,2.815687418);
   graph->SetPoint(206,95.81999969,2.892939568);
   graph->SetPoint(207,95.86000061,2.971222401);
   graph->SetPoint(208,95.90000153,3.050532341);
   graph->SetPoint(209,95.94000244,3.130866289);
   graph->SetPoint(210,95.98000336,3.212220907);


   graph->Sort();

   cout << "4e: " << graph->Eval(91.1876) << endl;


   //2e2mu
   TGraph *graph1 = new TGraph(211);
   graph1->SetName("Graph1");
   graph1->SetTitle("Graph1");
   graph1->SetFillColor(1);
   graph1->SetLineWidth(3);
   graph1->SetLineColor(kBlue);
   graph1->SetMarkerStyle(20);
   graph1->SetPoint(0,88.01999664,8.795412064);
   graph1->SetPoint(1,88.05999756,8.61370182);
   graph1->SetPoint(2,88.09999847,8.43285656);
   graph1->SetPoint(3,88.13999939,8.252916336);
   graph1->SetPoint(4,88.18000031,8.073918343);
   graph1->SetPoint(5,88.22000122,7.895903111);
   graph1->SetPoint(6,88.26000214,7.718908787);
   graph1->SetPoint(7,88.30000305,7.542974472);
   graph1->SetPoint(8,88.33999634,7.36813879);
   graph1->SetPoint(9,88.37999725,7.194441795);
   graph1->SetPoint(10,88.41999817,7.021921158);
   graph1->SetPoint(11,88.45999908,6.850616932);
   graph1->SetPoint(12,88.5,6.680567265);
   graph1->SetPoint(13,88.54000092,6.51181221);
   graph1->SetPoint(14,88.58000183,6.344389915);
   graph1->SetPoint(15,88.62000275,6.178339481);
   graph1->SetPoint(16,88.66000366,6.013700008);
   graph1->SetPoint(17,88.69999695,5.850510597);
   graph1->SetPoint(18,88.73999786,5.688809872);
   graph1->SetPoint(19,88.77999878,5.528635979);
   graph1->SetPoint(20,88.81999969,5.370028496);
   graph1->SetPoint(21,88.81999969,5.370028496);

void plotLimit(string outputDir="./", TString inputs="", TString inputs_blinded="", TString inputXSec="", bool strengthLimit=true, bool blind=false, double energy=7, double luminosity=5.035, TString legendName="ee and #mu#mu channels")
{
   setTDRStyle();  
   gStyle->SetPadTopMargin   (0.05);
   gStyle->SetPadBottomMargin(0.12);
   gStyle->SetPadRightMargin (0.16);
   gStyle->SetPadLeftMargin  (0.14);
   gStyle->SetTitleSize(0.04, "XYZ");
   gStyle->SetTitleXOffset(1.1);
   gStyle->SetTitleYOffset(1.45);
   gStyle->SetPalette(1);
   gStyle->SetNdivisions(505);
  
  //get the limits from the tree
  TFile* file = TFile::Open(inputs);
  printf("Looping on %s\n",inputs.Data());
  if(!file) return;
  if(file->IsZombie()) return;
  TFile* file_blinded = TFile::Open(inputs_blinded);
  printf("Looping on %s\n",inputs_blinded.Data());
  if(!file_blinded) return;
  if(file_blinded->IsZombie()) return;
  TTree* tree_blinded = (TTree*)file_blinded->Get("limit");
  tree_blinded->GetBranch("mh"              )->SetAddress(&Tmh      );
  tree_blinded->GetBranch("limit"           )->SetAddress(&Tlimit   );
  tree_blinded->GetBranch("limitErr"        )->SetAddress(&TlimitErr);
  tree_blinded->GetBranch("quantileExpected")->SetAddress(&TquantExp);
  TGraph* ExpLimitm2 = getLimitGraph(tree_blinded,0.025);
  TGraph* ExpLimitm1 = getLimitGraph(tree_blinded,0.160);
  TGraph* ExpLimit   = getLimitGraph(tree_blinded,0.500);
  TGraph* ExpLimitp1 = getLimitGraph(tree_blinded,0.840);
  TGraph* ExpLimitp2 = getLimitGraph(tree_blinded,0.975);
  file_blinded->Close(); 
  TTree* tree = (TTree*)file->Get("limit");
  tree->GetBranch("mh"              )->SetAddress(&Tmh      );
  tree->GetBranch("limit"           )->SetAddress(&Tlimit   );
  tree->GetBranch("limitErr"        )->SetAddress(&TlimitErr);
  tree->GetBranch("quantileExpected")->SetAddress(&TquantExp);
  TGraph* ObsLimit   = getLimitGraph(tree,-1   ); 
  file->Close(); 

  FILE* pFileSStrenght = fopen((outputDir+"SignalStrenght").c_str(),"w");
  std::cout << "Printing Signal Strenght" << std::endl;
  for(int i=0;i<ExpLimit->GetN();i++){
     double M = ExpLimit->GetX()[i];
     std::cout << "Mass: " << M << "; ExpLimit: " << ExpLimit->Eval(M) << std::endl; 
     printf("$%8.6E$ & $%8.6E$ & $[%8.6E,%8.6E]$ & $[%8.6E,%8.6E]$ \\\\\\hline\n",M, ExpLimit->Eval(M), ExpLimitm1->Eval(M), ExpLimitp1->Eval(M), ExpLimitm2->Eval(M),  ExpLimitp2->Eval(M));
     fprintf(pFileSStrenght, "$%8.6E$ & $%8.6E$ & $[%8.6E,%8.6E]$ & $[%8.6E,%8.6E]$ & $%8.6E$ \\\\\\hline\n",M, ExpLimit->Eval(M), ExpLimitm1->Eval(M), ExpLimitp1->Eval(M), ExpLimitm2->Eval(M),  ExpLimitp2->Eval(M), ObsLimit->Eval(M));
    if(int(ExpLimit->GetX()[i])%50!=0)continue; //printf("%f ",ObsLimit->Eval(M));
  }printf("\n");
  fclose(pFileSStrenght); 
 

  //get the pValue
  inputs = inputs.ReplaceAll("/LimitTree", "/PValueTree");
  file = TFile::Open(inputs);
  
  printf("Looping on %s\n",inputs.Data());
  if(!file) return;
  if(file->IsZombie()) return;
  
  tree = (TTree*)file->Get("limit");
  
  tree->GetBranch("limit"           )->SetAddress(&Tlimit   );
  
  TGraph* pValue     = getLimitGraph(tree,-1);
  
  file->Close();

  
  //make TH Cross-sections
   string suffix = outputDir;
   TGraph* THXSec   = Hxswg::utils::getXSec(outputDir); 
   scaleGraph(THXSec, 1000);  //convert cross-section to fb
   double cprime=1.0; double  brnew=0.0;
   double XSecScaleFactor = 1.0;
   if(suffix.find("_cp")!=string::npos){
     sscanf(suffix.c_str()+suffix.find("_cp"), "_cp%lf_brn%lf", &cprime, &brnew);
     XSecScaleFactor = pow(cprime,2) * (1-brnew);
   }
  //XSecScaleFactor = 0.001; //pb to fb
  scaleGraph(THXSec, XSecScaleFactor);


  string prod = "pp_SM";
  if(outputDir.find("ggH")!=std::string::npos)prod="gg";
  if(outputDir.find("qqH")!=std::string::npos)prod="qq";
  if(outputDir.find("ppH")!=std::string::npos)prod="pp";

  
  strengthLimit = false;
  if(prod=="pp_SM")strengthLimit=true;
 
  //TGraph *XSecMELA = Hxswg::utils::getXSecMELA(cprime);

  //Hxswg::utils::multiplyGraph(   ObsLimit, XSecMELA);
  //Hxswg::utils::multiplyGraph( ExpLimitm2, XSecMELA);
  //Hxswg::utils::multiplyGraph( ExpLimitm1, XSecMELA);
  //Hxswg::utils::multiplyGraph(   ExpLimit, XSecMELA);
  //Hxswg::utils::multiplyGraph( ExpLimitp1, XSecMELA);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
             gr_true_eff->GetErrorX(i), gr_true_eff->GetErrorY(i));
      if (!save_output) continue;
      gSystem->mkdir(dir_prefix.Data());
      tout = Form("%f %f %f %f", oux, ouy, ouxe, ouye);
      if (i == 0)
        tout = Form("Printf(\"%s\"); > %s/effi_%s", tout.Data(),
                    dir_prefix.Data(), grapht.Data());
      else
        tout = Form("Printf(\"%s\"); >> %s/effi_%s", tout.Data(),
                    dir_prefix.Data(), grapht.Data());
      //      Printf(":::::: %s", tout.Data());
      gROOT->ProcessLine(tout.Data());
    }
    // ------------------
    c = new TCanvas("cfinal", "mc_effi", 1200, 450);
    c->Divide(2, 1, 0.001, 0.001); c->Modified(); c->Draw();
    c->cd(1);
    gr_true->SetMinimum(0);
    gr_true->SetTitle(Form("phi (true & raw), %s", gtitle.Data()));
    gr_true->SetMarkerColor(kGreen+1);
    gr_true->Draw("AP");
    gr->SetMarkerColor(kRed+1);
    gr->Draw("P");
    c->cd(2)->SetGrid();
    gr_true_eff->SetMinimum(0);
    gr_true_eff->SetTitle(Form("efficiency, %s", grapht.Data()));
    gr_true_eff->SetMarkerColor(kGreen+1);
    gr_true_eff->Draw("AP");
    gr_e->SetMarkerColor(kRed+1);
    gr_e->Draw("P");
    //    c->SaveAs(Form("%s_%s.eps", blabla.Data(), grapht.Data()));
    return;
  }

  //  TGraph *geff = new TGraph(Form("effi_raw_%s", grapht.Data()));
  //  TGraph *geff = new TGraph(Form("effi_true_%s", grapht.Data()));
  //  TGraph *geff = new TGraph("effi_true_Phi2010_qualityonly");
  TGraph *geff = new TGraph("effi_true_PhiNsigma_qualityonly");
  if (geff->IsZombie()) return;
  geff->SetMarkerStyle(22);
  geff->SetMarkerColor(kBlack);
  geff->GetXaxis()->SetTitle("p_{t}, GeV/c");
  geff->SetTitle(Form("efficiency, %s", grapht.Data()));
  c = new TCanvas();
  c->SetGrid();
  geff->Draw("AP");
  Double_t tpcsigma = 9999.9;
  if (ilist == 1) tpcsigma = 1.0;
  if (ilist == 2) tpcsigma = 1.5;
  if (ilist == 3) tpcsigma = 2.0;
  if (ilist == 4) tpcsigma = 2.5;
  if (ilist == 5) tpcsigma = 3.0;
  Double_t sss = TMath::Erf(tpcsigma/TMath::Sqrt(2.0));
  if (noSigma) sss = 1.0;
  Printf("sigma = %10f", sss);

  // for (Int_t i = 0; i < count; i++)
  //   geff->GetY()[i] = (sss*sss)/(geff->GetY()[i]);
  // geff->SetMaximum(1.0);
  // geff->Draw("AP");

  for (Int_t i = 0; i < count; i++) {
    Double_t deno = geff->Eval(grx[i])*sss*sss;
    if (deno < 0.00001) deno = 1;
    gry_fix[i] = gry[i]/deno;
  }
  new TCanvas;
  TGraph *gr_fix = new TGraph(count, grx, gry_fix);
  gr_fix->SetMarkerStyle(21);
  gr_fix->SetMarkerColor(hh->GetLineColor());
  gr_fix->GetXaxis()->SetTitle("p_{t}, GeV/c");
  gr_fix->SetTitle(Form("corrected phi * #sigma^{2}, %s", gtitle.Data()));
  if (noSigma)
    gr_fix->SetTitle(Form("corrected phi (no #sigma), %s", gtitle.Data()));
  gr_fix->Draw("AP");

  //---------------------
  c = new TCanvas("cfinald", "data_correct", 1200, 450);
  c->Divide(2, 1, 0.001, 0.001); c->Modified(); c->Draw();
  c->cd(1);
  gr->SetMinimum(0);
  gr->SetMarkerColor(kBlack);
  gr->Draw("AP");
  c->cd(2);
  gr_fix->SetMinimum(0);
  gr_fix->SetMarkerColor(kGreen+3);
  gr_fix->Draw("AP");
  TString bla9 = Form("qualityonly_PID2_%s", grapht.Data());
  if (noSigma) bla9 = Form("%s_noSig.C", bla9.Data());
  else         bla9 = Form("%s.C", bla9.Data());
  //  gr_fix->SaveAs(bla9.Data());
  //  TPad *cp = new TPad("cpf", "", 0.45,0.45,0.99,0.92);
  TPad *cp = new TPad("cpf", "", 0.60,0.55,0.99,0.93);
  cp->SetLogy(); cp->Draw(); cp->cd();
  TGraph *cloneg = ((TGraph *)gr_fix->Clone());
  cloneg->SetTitle(); cloneg->SetMarkerSize(0.8);
  cloneg->Draw("AP");
  //  c->SaveAs(Form("%s_%s.eps", blabla.Data(), grapht.Data()));
  f->Close();
}

//.........这里部分代码省略.........
	if (check == "Uncor10t2_nConstituent_") PerfTitle="uncorr pt>10 #tracks > 2 & nConstituents >1";

	std::vector< TString > discriminators;
        for ( size_t itagger = 0; itagger < taggers.size(); ++itagger ) {
  		discriminators.push_back( check+"disc"+taggers[itagger]+"_udsg" );
	}
//	discriminators.push_back( "discTC3_udsg" );
//	discriminators.push_back( "discTP_udsg" );

	TCanvas *cv = new TCanvas("cv","cv",900,900);
	TMultiGraph *mg =new TMultiGraph();
	TLegend *legend0 = new TLegend(0.68,0.12,0.88,0.32);
        // inverted axis
	TMultiGraph *mginv =new TMultiGraph();
	TLegend *legend1 = new TLegend(0.65,0.18,0.85,0.48);
	std::ofstream salida("BTagPATop3"+check+".txt");
	
	for ( size_t itagger = 0; itagger < taggers.size(); ++itagger ) {

		TString tag    = check+"g"+taggers[itagger]+"_udsg";
		TGraphErrors *agraph = (TGraphErrors*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+tag);

		
		TGraph *dgraph = (TGraph*) gDirectory->Get("BTagPATAnalyzer"+taggers[itagger]+"/"+taggers[itagger]+"/"+discriminators[itagger]);
//		TGraph *udsgvsdgraph = new TGraph(dgraph->GetN(),dgraph->GetY(),dgraph->GetX());
		dgraph->Sort();
//		udsgvsdgraph->Sort();

		TGraphErrors *g = new TGraphErrors(agraph->GetN(),agraph->GetY(),agraph->GetX(),agraph->GetEY(),agraph->GetEX());
		g->Sort();
		g->SetLineColor(itagger+1);
                legend0 -> AddEntry(g,taggers[itagger],"l");
		mg->Add(g);

                //inverted axis
		TGraphErrors *ginv = new TGraphErrors(agraph->GetN(),agraph->GetX(),agraph->GetY(),agraph->GetEX(),agraph->GetEY());
		ginv->Sort();
		ginv->SetLineColor(colorlines[itagger]);
		ginv->SetMarkerStyle(8);
		ginv->SetMarkerColor(colorlines[itagger]);
                legend1 -> AddEntry(ginv,labeltag[itagger],"l");
		mginv->Add(ginv);

		std::cout << " Tagger: " << tag << std::endl;
		std::cout << " Loose(10%): " << " cut > " << std::setprecision(4) << dgraph->Eval(0.1) << " b-eff = " << g->Eval(0.1) << std::endl;
		std::cout << " Medium(1%):  " << " cut > " << std::setprecision(4) << dgraph->Eval(0.01) << " b-eff = " << g->Eval(0.01) << std::endl;
		std::cout << " Tight(0.1%) = " << " cut > " << std::setprecision(4) << dgraph->Eval(0.001) << " b-eff = " << g->Eval(0.001) << std::endl;
                salida << " " << taggers[itagger] << endl;
	        salida << " #bin \t b-eff  \t err-beff \t  non-beff \t err-non-b"<< endl;
        	for ( size_t i=0;i< agraph->GetN();i++){ 
			salida	<< "  " << (i+1) << " \t "
				<< agraph->GetX()[i] <<" \t " 
				<< agraph->GetEX()[i]<<" \t "
				<< agraph->GetY()[i] <<" \t "
				<< agraph->GetEY()[i]
                       		<< " "<< endl;
		}

	}
//	cv->SetLogx();
	mg->Draw("ALP");
	mg->GetYaxis()->SetTitle("b-eff");
	mg->GetXaxis()->SetTitle("udsg-mistagging");
	legend0 -> Draw();
        cv-> SetGrid();
	cv-> Print ("BTagPATop"+check+".eps");
	cv-> Print ("BTagPATop"+check+".png");
	cv->cd(0);
	cv->SetLogx();
	mg->Draw("ALP");
	legend0 -> Draw();
	cv-> Print ("BTagPATop1"+check+".eps");
	cv-> Print ("BTagPATop1"+check+".png");
        //inverted axis
	TCanvas *cvinv = new TCanvas("cvinv","cvinv",700,700);
        cvinv->SetLogy(0);
	mginv->Draw("ALP");
	mginv->GetXaxis()->SetTitle("b-eff");
	mginv->GetYaxis()->SetTitle("udsg-mistagging");
	legend1 -> Draw();
        cvinv->SetGrid();
	cvinv-> Print ("BTagPATop2"+check+".eps");
	cvinv-> Print ("BTagPATop2"+check+".png");
	cvinv->cd(0);
//	cvinv->Update();
//	cvinv->SetLogx();
//	mginv->SetXmin(10^-4);
//      axis range using a histogram helper
    TH2D*hpx = new TH2D("hpx",PerfTitle,200,0.0,1.005,200,0.00002,1.05);
        hpx->SetStats(kFALSE);
		hpx->Draw();
	hpx->GetXaxis()->SetTitle("b-eff");
	hpx->GetYaxis()->SetTitle("udsg-mistagging");
	cvinv->SetLogy(1);
	mginv->Draw("AP");
	legend1 -> Draw();
	cvinv-> Print ("BTagPATop3"+check+".eps");
	cvinv-> Print ("BTagPATop3"+check+".png");
	//ed}
}

//.........这里部分代码省略.........
			if ( hardjets[j].getVarF("jn_jp") > maxJetBTag && fabs(jet.Eta()) < 2.5 )
				maxJetBTag = hardjets[j].getVarF("jn_jp");
			double tempDelPhiJetMet = deltaPhi(met.Phi(), jet.Phi());
			if ( tempDelPhiJetMet < minDeltaPhiJetMet )
				minDeltaPhiJetMet = tempDelPhiJetMet;
		}

		nvtx = *nvtxP;

		if (isData)
			ni = -1;
		else
			ni = *niP;

		if (isSignal) {
			const int nMC = ev.getSVV<int>("mcn");
			const int * mcID = ev.getAVA<int>("mc_id");
			int hIdx = 0;
			for (; hIdx < nMC; ++hIdx)
				if (fabs(mcID[hIdx]) == 25)
					break;
			if (hIdx == nMC)
				throw string("ERROR: Higgs not found in signal sample!");

			float Hpx = ev.getAVV<float>("mc_px", hIdx);
			float Hpy = ev.getAVV<float>("mc_py", hIdx);
			float Hpz = ev.getAVV<float>("mc_pz", hIdx);
			float Hen = ev.getAVV<float>("mc_en", hIdx);
			TLorentzVector higgs;
			higgs.SetPxPyPzE( Hpx, Hpy, Hpz, Hen );
			hmass = higgs.M();

			if (higgsW) {
				hweight = higgsW->Eval(hmass);
				if (higgsI)
					hweight *= higgsI->Eval(hmass);
			} else
				hweight = 1;
		}

		if ( opt.checkBoolOption("ADDITIONAL_LEPTON_VETO") && (type == ELE || type == MU || type == EMU) && ((nele + nmu + nsoftmu) > 2) )
			continue;
		if ( opt.checkBoolOption("ADDITIONAL_LEPTON_VETO") && (type == PHOT) && ((nele + nmu + nsoftmu) > 0) )
			continue;
		if ( opt.checkBoolOption("ZPT_CUT") && zpt < 55 )
			continue;
		// for different background estimation methods different window should be applied:
		// * sample for photons should have 76.0 < zmass < 106.0
		// * sample for non-resonant background should not have this cut applied
		if ( opt.checkBoolOption("TIGHT_ZMASS_CUT") && (type == ELE || type == MU) && (zmass < 76.0 || zmass > 106.0))
			continue;
		if ( opt.checkBoolOption("WIDE_ZMASS_CUT") && (type == ELE || type == MU) && (zmass < 76.0 || zmass > 106.0))
			continue;
		if ( opt.checkBoolOption("BTAG_CUT") && ( maxJetBTag > 0.264) )
			continue;
		if ( opt.checkBoolOption("DPHI_CUT") && ( minDeltaPhiJetMet < 0.5) )
			continue;


#ifdef PRINTEVENTS
		evPrint.setJetCollection(hardjets);
		evPrint.setMET(met);
		evPrint.setMT(mt);
		string channelType;
		if (type == ELE)
			channelType = "ee";

//.........这里部分代码省略.........
  
 std::cerr << "... I'm minimizing ... DATA analysis" << std::endl;
 std::cerr << ">>>>>>> numEvents = " << numEvents << " => " << vET_data.size() << " selected (=" << mylist->GetN() << ")" << std::endl;
 numSelectedData = vET_data.size();
 
 
 ///===== Chi2 ====
 std::cerr << " === Chi2 === " << std::endl;
 minuit->SetFunction(functorChi2);

 TGraph * grChi2 = new TGraph(iNoSteps);
 minuit->Scan(iPar_NoBG,iNoSteps,grChi2->GetX(),grChi2->GetY(),MinScan,MaxScan);

// TGraph * grChi2 = new TGraph();
// for (int iStep = 0; iStep < iNoSteps; iStep++){
//  double x = MinScan + (MaxScan - MinScan) / iNoSteps * (iStep+0.5);
//  double y = Chi2F(&x);
//  grChi2->SetPoint(iStep+1,x,y);
// }
 grChi2->Draw("AL");
 outFile->cd();
 minuit->PrintResults();
 outFile->cd();
 grChi2->SetTitle("grChi2");
 grChi2->Write();
 const double *outParametersTemp = minuit->X();
 const double *errParametersTemp = minuit->Errors();
  
 double *outParameters = new double;
 double *errParameters = new double;
 outParameters[0] = outParametersTemp[0];
 errParameters[0] = errParametersTemp[0];
   
 double minChi2 = grChi2->Eval(outParameters[0]);
 std::cerr << " numEvents = " << numEvents << " Scale = " << outParameters[0] << " +/- " << errParameters[0] << std::endl;
 ///===== end Chi2 ====
  
 ///==== likelihood ====
 std::cerr << " === LL === " << std::endl;
 minuit->SetFunction(functorLL);

 TGraph * grLL_temp = new TGraph(iNoSteps);
 minuit->Scan(iPar_NoBG,iNoSteps,grLL_temp->GetX(),grLL_temp->GetY(),MinScan,MaxScan);
 TGraph * grLL = new TGraph();
 grLL->SetName("grLL");
 int nPointLL = 0;
 for (unsigned int iStep = 0; iStep < iNoSteps; iStep++){
  double x = MinScan + (MaxScan - MinScan) / iNoSteps * (iStep+0.5);
  double y = LLFunc(&x);
//  std::cerr << " y = " << y << std::endl;
  if (y != numberDATA * numEvents) {
   std::cerr << " Ok y = " << y << std::endl;
   grLL->SetPoint(nPointLL,x,y);
   nPointLL++;
  }
 }
 
 
 std::cerr << " finito " << std::endl;
 
 grLL->Draw("AL");
 outFile->cd();
 minuit->PrintResults();
 outFile->cd();
 grLL->SetTitle("grLL");
 grLL->Write();

///**** LL ****
void doMC_LL(){
 TTree* MyTreeMC = (TTree*) fileInMC->Get(treeNameMC.c_str());

 for (nIter = 0; nIter<maxIter; nIter++){
  if (!(nIter%1)) std::cerr << ">>> nIter = " << nIter << " : " << maxIter << std::endl;
  vET_data.clear();
  outFile->cd();
  TString nameDATA = Form("hDATA_%d_%d_%.5f",Data_or_MC,nIter,ScaleTrue);
  TH1F hDATA(nameDATA,nameDATA,numBINS,minBINS,maxBINS);
  
  MyTreeMC->Draw(">> myListMC",(AdditionalCut+Form("&& (ET * (1+(%f)))>%f",ScaleTrue,minET)).Data(),"entrylist");
  TEntryList *myListMC = (TEntryList*)gDirectory->Get("myListMC");
  
  MyTreeMC->SetEntryList(0); 
  TEntryList *listMCHere = new TEntryList("listMCHere","listMCHere");    
  for (int iEvt = 0; iEvt < numSelectedData; iEvt ++){
   listMCHere->Enter(myListMC->GetEntry(gRandom->Uniform(0,myListMC->GetN())));
  }
  
  MyTreeMC->SetEntryList(listMCHere);
  MyTreeMC->Draw(Form("(1+%f) * %s >> %s",ScaleTrue,variableName.c_str(),nameDATA.Data()));
  
  ConvertStdVectDouble(vET_data,MyTreeMC->GetV1(),numSelectedData);
  
  ///==== likelihood ====
 std::cerr << " === LL === " << std::endl;
 std::cerr << " === pseudo vET_data.size() = " << vET_data.size() << std::endl;
 
 
  minuit->SetFunction(functorLL); 
  TGraph * grLL_temp = new TGraph(iNoSteps);
  minuit->Scan(iPar_NoBG,iNoSteps,grLL_temp->GetX(),grLL_temp->GetY(),MinScan,MaxScan);
  TGraph * grLL = new TGraph();
  int nPointLL = 0;
 for (unsigned int iStep = 0; iStep < iNoSteps; iStep++){
   double x = MinScan + (MaxScan - MinScan) / iNoSteps * (iStep+0.5);
   double y = LLFunc(&x);
   if (y != numberDATA * numEvents) {
    grLL->SetPoint(nPointLL,x,y);
    nPointLL++;
   }
  }
  grLL->Draw("AL");
  outFile->cd();
  minuit->PrintResults();
  const double *outParametersTemp2 = minuit->X();
  const double *errParametersTemp2 = minuit->Errors();
  
  double *outParametersLL = new double;
  double *errParametersLL = new double;
  outParametersLL[0] = outParametersTemp2[0];
  errParametersLL[0] = errParametersTemp2[0];
  
  std::cerr << " nPointLL = " << nPointLL << std::endl;
  
  double minLL = grLL->Eval(outParametersLL[0]);
  ///==== end likelihood ====
  
    ///==== Save the whole shape of LL/Chi2 ====
  for (unsigned int ii=0; ii < iNoSteps; ii++){
   double X_ii = (MaxScan - MinScan) / iNoSteps * ii + MinScan;
   
   Alpha   = X_ii;
   Chi2    = 0;
   LL      = grLL->Eval(X_ii);
   NewChi2 = 0;
   
   myTreeChi2->Fill();
  }

  ///===== Look for minima =====
  double a;
  double b;
  double c;
  
  
  double errX_low = -9999;
  double errX_up = 9999;
  int err_low = 0;
  int err_up = 0;
  for (unsigned int ii=0; ii < iNoSteps; ii++){
   double X_ii = (MaxScan - MinScan) / iNoSteps * ii + MinScan;
   double here = grLL->Eval(X_ii);
   if (err_low == 0){
    if (here < (minLL + DELTA_LL)){
     errX_low = X_ii;
     err_low = 1;
    }
   }
   else if (err_up == 0 && here > (minLL + DELTA_LL) && X_ii > outParametersLL[0]){
    errX_up = X_ii; 
    err_up = 1;
   }
  }
  
  AlphaMean = outParametersLL[0];
  AlphaMinus = errX_low;
  AlphaPlus = errX_up;
 
//.........这里部分代码省略.........

void bjetefficiency() { 	  	  	
std::vector<string> samples_;
//samples_.push_back("WPHJET.root");
//samples_.push_back("WW.root");
//samples_.push_back("WZ.root");
//samples_.push_back("ZZ.root");
samples_.push_back("SIGNALtGC.root");
samples_.push_back("SIGNALtGu.root");

//open files
std::vector<TFile*> files;
for(unsigned int idx=0; idx<samples_.size(); ++idx){
files.push_back(new TFile(samples_[idx].c_str()));}
std::vector<string> variables2_;
/*
variables2_.push_back("btageffL");
variables2_.push_back("btageffM");
variables2_.push_back("btageffT");

variables2_.push_back("ctageffL");
variables2_.push_back("ctageffM");
variables2_.push_back("ctageffT");
*/
variables2_.push_back("btageffL");
variables2_.push_back("btageffM");
variables2_.push_back("btageffT");


std::vector<double> finaleff;

Double_t bdiscriminator[35];
for (unsigned int i=0; i<15; ++i){
bdiscriminator[i] =0.244-0.1+i*0.02; }
for (unsigned int i=0; i<10; ++i){
bdiscriminator[15+i] = 0.679-0.1+i*0.02;}
for (unsigned int i=0; i<10; ++i){
bdiscriminator[25+i] = 0.898-0.1+i*0.02;}

Double_t count[3];
count[0]=15;
count[1]=10;
count[2]=10;

for(unsigned int ii=0; ii<variables2_.size(); ++ii){
std::vector<TH1F*> hists;
for(unsigned int idx=0; idx<files.size(); ++idx){hists.push_back((TH1F*)files[idx]->Get((std::string("STEP1/").append(variables2_[ii]).c_str())));}

for(unsigned int idx=1; idx<files.size(); ++idx){
hists[idx]->Add(hists[idx-1]);}

//hists[files.size()-1]->Draw();
std::vector<double> eff;
int bin; 

for ( bin=1; bin<count[ii]*2+1; ++bin){eff.push_back(hists[files.size()-1]->GetBinContent(bin));}
for (unsigned int i=0; i<count[ii]; ++i){finaleff.push_back(eff[2*i+1]/(eff[2*i]+eff[2*i+1])); 
//cout<< eff[2*i+1]<<"                "<<eff[2*i+1]/(eff[2*i]+eff[2*i+1])<<endl;
}
}

//   TGraph *gr = new TGraph(30,bdiscriminator,finaleff);

Double_t effi[35];
for (unsigned int i=0; i<35; ++i){
effi[i] =finaleff[i]; }
   TGraph *gr = new TGraph(35,effi,bdiscriminator);
   TGraph *grREV = new TGraph(35,bdiscriminator,effi);

   cout<<"the new working points for nominal SF"<<endl;
   cout<<"loose        "<<gr->Eval(1.008*grREV->Eval(0.244))<<endl;
   cout<<"medium        "<<gr->Eval(0.963*grREV->Eval(0.679))<<endl;
   cout<<"tight        "<<gr->Eval(0.947*grREV->Eval(0.898))<<endl;


   cout<<"the new working points for SF UP"<<endl;
   cout<<"loose        "<<gr->Eval((1.008+0.023)*grREV->Eval(0.244))<<endl;
   cout<<"medium        "<<gr->Eval((0.963+0.020)*grREV->Eval(0.679))<<endl;
   cout<<"tight        "<<gr->Eval((0.947+0.025)*grREV->Eval(0.898))<<endl;


   cout<<"the new working points for SF DOWN"<<endl;
   cout<<"loose        "<<gr->Eval((1.008-0.023)*grREV->Eval(0.244))<<endl;
   cout<<"medium        "<<gr->Eval((0.963-0.020)*grREV->Eval(0.679))<<endl;
   cout<<"tight        "<<gr->Eval((0.947-0.025)*grREV->Eval(0.898))<<endl;

//   cout<<gr->Eval(1.008*grREV->Eval(0.244))<<endl;
//   cout<<gr->Eval(0.966*grREV->Eval(0.679))<<endl;
//   cout<<gr->Eval(0.931*grREV->Eval(0.898))<<endl;

//   cout<<gr->Eval(0.0990)<<endl;
 //  cout<<gr->Eval(0.0142)<<endl;
//   cout<<gr->Eval(0.00155)<<endl;
//   cout<<gr->Eval(0.0142)<<endl;
 //  cout<<gr->Eval(0.0016)<<endl;
   gr->SetLineColor(4);
   gr->SetLineWidth(4);
   gr->SetMarkerColor(4);
   gr->SetMarkerStyle(21);

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

//.........这里部分代码省略.........
    if (cat==4 || cat==5 || cat==6 || cat==7 || cat==8) {
      modelBuilder.addBkgPdf("Bernstein",3,Form("pol3_cat%d",cat));
      modelBuilder.addBkgPdf("Bernstein",4,Form("pol4_cat%d",cat));
      /*
      modelBuilder.addBkgPdf("PowerLaw",1,Form("pow1_cat%d",cat));
      modelBuilder.addBkgPdf("PowerLaw",3,Form("pow3_cat%d",cat));
      modelBuilder.addBkgPdf("Exponential",1,Form("exp1_cat%d",cat));
      modelBuilder.addBkgPdf("Exponential",3,Form("exp3_cat%d",cat));
      modelBuilder.addBkgPdf("Laurent",1,Form("lau1_cat%d",cat));
      modelBuilder.addBkgPdf("Laurent",3,Form("lau3_cat%d",cat));
      */
    }
    map<string,RooAbsPdf*> bkgPdfs = modelBuilder.getBkgPdfs();
    modelBuilder.setSignalPdfFromMC(sigMC);
    modelBuilder.makeSBPdfs();
    map<string,RooAbsPdf*> sbPdfs = modelBuilder.getSBPdfs();

    modelBuilder.fitToData(dataBinned,true,true);
    modelBuilder.fitToData(dataBinned,false,true);

    modelBuilder.throwToy(Form("cat%d_toy0",cat),dataBinned->sumEntries(),true,true);

    // Profile this category using ProfileMultiplePdfs
    ProfileMultiplePdfs profiler;
    for (map<string,RooAbsPdf*>::iterator pdf=sbPdfs.begin(); pdf!=sbPdfs.end(); pdf++) {
      string bkgOnlyName = pdf->first.substr(pdf->first.find("sb_")+3,string::npos);
      if (bkgPdfs.find(bkgOnlyName)==bkgPdfs.end()){
        cerr << "ERROR -- couldn't find bkg only pdf " << bkgOnlyName << " for SB pdf " << pdf->first << endl;
        pdf->second->fitTo(*dataBinned);
        exit(1);
      }
      int nParams = bkgPdfs[bkgOnlyName]->getVariables()->getSize()-1;
      profiler.addPdf(pdf->second,2*nParams);
      //profiler.addPdf(pdf->second);
      cout << pdf->second->GetName() << " nParams=" << pdf->second->getVariables()->getSize() << " nBkgParams=" << nParams << endl;
    }
    profiler.printPdfs();
    //cout << "Continue?" << endl;
    //string bus; cin >> bus;
    profiler.plotNominalFits(dataBinned,mass,80,Form("cat%d",cat));
    pair<double,map<string,TGraph*> > minNlls = profiler.profileLikelihood(dataBinned,mass,mu,mu_low,mu_high,mu_step);
    pair<double,map<string,TGraph*> > correctedNlls = profiler.computeEnvelope(minNlls,Form("cat%d",cat),2.);
    minNlltrack.push_back(make_pair(correctedNlls.first,correctedNlls.second["envelope"]));
    //minNlls.second.insert(pair<string,TGraph*>("envelope",envelopeNll.second));
    //map<string,TGraph*> minNLLs = profiler.profileLikelihoodEnvelope(dataBinned,mu,mu_low,mu_high,mu_step);
    profiler.plot(correctedNlls.second,Form("cat%d_nlls",cat));
    //profiler.print(minNLLs,mu_low,mu_high,mu_step);
    /*
    if (minNLLs.find("envelope")==minNLLs.end()){
      cerr << "ERROR -- envelope TGraph not found in minNLLs" << endl;
      exit(1);
    }
    */
    //minNlltrack.push_back(make_pair(profiler.getGlobalMinNLL(),minNLLs["envelope"]));
  }
  //exit(1);
  TGraph *comb = new TGraph();
  for (vector<pair<double,TGraph*> >::iterator it=minNlltrack.begin(); it!=minNlltrack.end(); it++){
    if (it->second->GetN()!=minNlltrack.begin()->second->GetN()){
      cerr << "ERROR -- unequal number of points for TGraphs " << it->second->GetName() << " and " << minNlltrack.begin()->second->GetName() << endl;
      exit(1);
    }
  }
  for (int p=0; p<minNlltrack.begin()->second->GetN(); p++){
    double x,y,sumy=0;
    for (vector<pair<double,TGraph*> >::iterator it=minNlltrack.begin(); it!=minNlltrack.end(); it++){
      it->second->GetPoint(p,x,y);
      sumy += (y+it->first);
    }
    comb->SetPoint(p,x,sumy);
  }
  pair<double,double> globalMin = getGraphMin(comb);
  for (int p=0; p<comb->GetN(); p++){
    double x,y;
    comb->GetPoint(p,x,y);
    comb->SetPoint(p,x,y-globalMin.second);
  }
  vector<double> fitVal = getValsFromLikelihood(comb);

  cout << "Best fit.." << endl;
  cout << "\t mu = " << Form("%4.3f",fitVal[0]) << " +/- (1sig) = " << fitVal[2]-fitVal[0] << " / " << fitVal[0]-fitVal[1] << endl;
  cout << "\t      " << "    " << " +/- (2sig) = " << fitVal[4]-fitVal[0] << " / " << fitVal[0]-fitVal[3] << endl;

  cout << comb->Eval(fitVal[0]) << " " << comb->Eval(fitVal[1]) << " " << comb->Eval(fitVal[2]) << " " << comb->Eval(fitVal[3]) << " " << comb->Eval(fitVal[4]) << endl;

  double quadInterpVal = ProfileMultiplePdfs::quadInterpMinimum(comb);
  cout << "quadInterp: mu = " << quadInterpVal << endl;
  cout << "\t " << comb->Eval(quadInterpVal) << " " << comb->Eval(quadInterpVal-0.005) << " " << comb->Eval(quadInterpVal-0.01) << " " << comb->Eval(quadInterpVal+0.005) << " " << comb->Eval(quadInterpVal+0.01) << endl;
  
  comb->SetLineWidth(2);
  TCanvas *canv = new TCanvas();
  comb->Draw("ALP");
  canv->Print("plots/comb.pdf");
  TFile *tempOut = new TFile("tempOut.root","RECREATE");
  tempOut->cd();
  comb->SetName("comb");
  comb->Write();
  tempOut->Close();
  return 0;
}

void ExtractSignalYield(string ModelName, string decayMode, float natural_width, bool dependsOnKtilde=0)
{
  
  string dir = "/usr/users/dschaefer/root/results/";
  string outputfile = ModelName+"_"+decayMode+"_signalYield.root";
  //float natural_width =0.05;
  std::stringstream s;
      if(int(natural_width*100+0.5) == 10 or int(100*natural_width+0.5)==20 or int(natural_width*100+0.5)==30)
      {
	s << std::fixed << std::setprecision(1) << natural_width;
      }
      else
      {
	s << std::fixed << std::setprecision(2) << natural_width;
      }
  string swidth = s.str();
  float ktilde;
  string sktilde;
  if(dependsOnKtilde)
  {
   ktilde = natural_width;
   sktilde = swidth;
  }
  
  TGraph* limit = getObservedLimit(ModelName,decayMode,natural_width,dependsOnKtilde);
  TGraph* limitUP = getLimitUp(ModelName,decayMode,natural_width,1,dependsOnKtilde);
  TGraph* limitDOWN = getLimitDown(ModelName,decayMode,natural_width,1,dependsOnKtilde);
  TGraph* limitUP2 = getLimitUp(ModelName,decayMode,natural_width,2,dependsOnKtilde);
  TGraph* limitDOWN2 = getLimitDown(ModelName,decayMode,natural_width,2,dependsOnKtilde);
  
  
  int number;
  float* masses;
  if(decayMode.find("lvjj")!=string::npos)
  {
    masses = new float[38];
    number =38;
    float array[38] = {800,900,1000,1100,1200,1300,1400,1500,1600,1700,1800,1900,2000,2100,2200,2300,2400,2500,2600,2700,2800,2900,3000,3100,3200,3300,3400,3500,3600,3700,3800,3900,4000,4100,4200,4300,4400,4500};
    masses = array;
  }
  else
  {
   number = 29;
   masses = new float[29];
   float array[29] = {1200,1300,1400,1500,1600,1700,1800,1900,2000,2100,2200,2300,2400,2500,2600,2700,2800,2900,3000,3100,3200,3300,3400,3500,3600,3700,3800,3900,4000};
   masses = array;
  }
  TGraph* gryellow = new TGraph(2*number);
  TGraph* grgreen  = new TGraph(2*number);
  
  for(int m=0;m<number;m++)
  {
    gryellow->SetPoint(m,masses[m],limitUP->Eval(masses[m]));
    gryellow->SetPoint(number+m,masses[number-m-1],limitDOWN->Eval(masses[number-m-1]));
    grgreen->SetPoint(m,masses[m],limitUP2->Eval(masses[m]));
    grgreen->SetPoint(number+m,masses[number-m-1],limitDOWN2->Eval(masses[number-m-1]));
    
  }
  gryellow->SetFillColor(kYellow);
  grgreen->SetFillColor(kGreen);
  
  TCanvas* test = new TCanvas("test","test",400,400);
  test->SetLogy();
  limitUP->SetFillColor(kGreen+2);
  limitUP->SetLineColor(kGreen+2);
  limitUP->SetLineWidth(-1504);
  
  //limitUP->Draw("ACP");
  limitDOWN->SetFillColor(kBlue);
  limitDOWN->SetLineColor(kGreen+2);
  limitDOWN->SetLineWidth(-1504);
  //limitDOWN->Draw("CPsame");
  gryellow->Draw("AFsame");
  grgreen->Draw("Fsame");
  limit->SetLineColor(kBlue);
  limit->SetLineWidth(3);
  limit->Draw("Lsame");
  test->Update();
  
  float lumi = 2.10*1000;  //pb^-1
  if(decayMode.find("jjjj")!=string::npos){lumi = 2.6*1000;}//pb^-1
  
  
  
  
  //float *xs = getTheoryCrossSection(ModelName,0.1);
  
  //x values for which a theory prediction for production cross section exists
  float m_xs[10]  =  {800,900,1000,1500,1800,2000,2500,3000,3500,4500};
  float BR_G_to_WW[10];
  float BR_G_to_ZZ[10];
  float *acceptance;
  float *N_expected;
  //float xs = 0.01;
  float *xs = new float[10];
  
  if(dependsOnKtilde)
  {
    xs = getTheoryCrossSection(ModelName,ktilde);
  }
//.........这里部分代码省略.........