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

TGraph *MarcosExclusionLine(TH2F *exclusionshape, int scantype) {
//  write_warning(__FUNCTION__,"USING MARIAS ALGORITHM...");
//  return get_exclusion_line(exclusionshape);
  
  TH2F *fakehisto = flipth2f(exclusionshape);
  TGraph *fakegraph = get_mSUGRA_exclusion_line(fakehisto, scantype);
  TGraph *realgraph = new TGraph(fakegraph->GetN());
  double x,y;
  float last_x=0;
  float last_y=0;
  int counter=0;
  for(int i=0;i<fakegraph->GetN();i++) {
	fakegraph->GetPoint(i,x,y);
	if(scantype==PlottingSetup::SMS) {
	  if(y-x<75) {
	    realgraph->SetPoint(counter,last_x,last_y);
	    counter++;
	    continue;
	  }
	}
	realgraph->SetPoint(counter,y,x);
	last_x=y;
	last_y=x;
	counter++;
  }
  realgraph->SetLineColor(TColor::GetColor("#151515")); //nice black
  realgraph->SetLineWidth(2);
  //realgraph->SetLineWidth(4);//paper style
  
  delete fakegraph;

  return realgraph;
}

void ButkevichDrawer::adjust_graph(TGraph& graph) {
  for (int iPoint=0;iPoint<graph.GetN();iPoint++) {
    double x,y;
    graph.GetPoint(iPoint,x,y);
    graph.SetPoint(iPoint,x,y*1.0e-38);
  }
}

std::vector<double> GetCrossings(TGraph const& g, double cross) {
  std::vector<double> result;
  unsigned n = g.GetN();
  double x1 = 0;
  double y1 = 0;
  double x2 = 0;
  double y2 = 0;
  for (unsigned i = 0; i < (n-1); ++i) {
    g.GetPoint(i, x1, y1);
    g.GetPoint(i+1, x2, y2);
    if ( (y2-cross)*(y1-cross) < 0.0 ) {
      double m = (y2-y1)/(x2-x1);
      double c = (y1 - m * x1);
      double xc = (cross - c) / m;
      result.push_back(xc);
      std::cout << "Crossing at " << xc << std::endl;
    }
  }
  return result;
}

void KVIDentifier::CopyGraph(const TGraph& graph)
{
   // Copy coordinates of points from the TGraph

   Double_t x, y;
   //change number of points
   Set(graph.GetN());
   for (int i = 0; i < GetN(); i++) {
      graph.GetPoint(i, x, y);
      SetPoint(i, x, y);
   }
}

TGraph* ContourGraph( TH2F* hist,double xmin=16, double xmax=90) {

    //temporary canvas
    TCanvas* MOO = new TCanvas( TString::Format("dummy_canvas_%s", hist->GetName()), "A scan of m_{0} versus m_{12}", 0, 0, 650,640);
    MOO->cd();

    TGraph* gr0 = new TGraph();
    TH2F* h = (TH2F*)hist->Clone();
    TGraph* gr = (TGraph*)gr0->Clone(TString::Format("gr_%s", h->GetName()));

    cout << "==> Will dumb histogram: " << h->GetName() << " into a graph" <<endl;

    h->SetContour( 1 );
    //h->GetXaxis()->SetRangeUser(250,1200);
    h->GetXaxis()->SetRangeUser(xmin, xmax);
    //h->GetYaxis()->SetRangeUser(2,50);

    double pval = CombinationGlob::cl_percent[1];
    std::cout << pval << std::endl; 
    double signif = TMath::NormQuantile(1-pval);
    h->SetContourLevel( 0, signif );
    h->Draw("CONT LIST");
    h->SetDirectory(0);
    gPad->Update();

    TObjArray *contours = (TObjArray*) gROOT->GetListOfSpecials()->FindObject("contours");
    Int_t ncontours     = contours->GetSize();
    cout << "Found " << ncontours << " contours " << endl;

    TList *list = (TList*)contours->At(0);
    contours->Print("v");
    if(!list) return NULL;

    gr = (TGraph*)list->First();
    if(!gr) return NULL;

    gr->SetName(TString::Format("gr_%s", hist->GetName()));
    //gr->SetName(hist->GetName());
    int N = gr->GetN();
    double x0, y0;
    for(int j=0; j<N; j++) {
        gr->GetPoint(j,x0,y0);
        cout << j << " : " << x0 << " : "<<y0 << endl;
    }
    //  //  gr->SetMarkerSize(2.0);    
    //gr->Draw("ALP");

    delete MOO;

    cout << "Generated graph " << gr << " with name " << gr->GetName() << endl;
    return gr;
}

void drawtext()
{
   Int_t i,n;
   Double_t x,y;
   TLatex *l;

   TGraph *g = (TGraph*)gPad->GetListOfPrimitives()->FindObject("Graph");
   n = g->GetN();
   for (i=1; i<n; i++) {
      g->GetPoint(i,x,y);
      l = new TLatex(x,y+0.2,Form("%4.2f",y));
      l->SetTextSize(0.025);
      l->SetTextFont(42);
      l->SetTextAlign(21);
      l->Paint();
   }
}

vector< pair <double,double> > getMinMaxOfGraphs (vector<GraphInfo> inputGraphs) {

  cout << "-----------------Inside minMax -------------------" << endl;
  pair<double,double> yValMinMax;
  pair<double,double> xValMinMax;
  
  yValMinMax.first = 1000.0; //min
  yValMinMax.second = -1000.0; //max

  xValMinMax.first = 1000.0; //min
  xValMinMax.second = -1000.0; //max
  
  
  for (unsigned iGraph =0; iGraph < inputGraphs.size(); iGraph++){
    TGraph * thisGraph = inputGraphs[iGraph].theGraph;
    
    double iMinY = 1000.0;
    double iMaxY = -1000.0;
    double iMinX = 1000.0;
    double iMaxX = -1000.0;

    
    for (int iPoint =0; iPoint < thisGraph->GetN(); iPoint++){
      double xVal, yVal;
      thisGraph->GetPoint(iPoint, xVal, yVal);
      cout <<"n = " << iPoint << "x= " << xVal << "y = " << yVal << endl;
      if (yVal < iMinY )
        iMinY = yVal;
      if (yVal > iMaxY)
        iMaxY = yVal;

      if (xVal < iMinX )
        iMinX = xVal;
      if (xVal > iMaxX)
        iMaxX = xVal;
      
    }
    
    cout << "New plot minY = " << iMinY << ", maxY = " << iMaxY << endl
         << "         minX = " << iMinX << ", maxX = " << iMaxX << endl ;

    if (iMinY < yValMinMax.first)
      yValMinMax.first = iMinY;
    if (iMaxY > yValMinMax.second)
      yValMinMax.second = iMaxY;

    if (iMinX < xValMinMax.first)
      xValMinMax.first = iMinX;
    if (iMaxX > xValMinMax.second)
      xValMinMax.second = iMaxX;
    
  }

  //------------------------

  cout << "found min = " << yValMinMax.first
       << "   found max = " << yValMinMax.second
       << endl;

  if (yValMinMax.first > 0 )
    yValMinMax.first = 0;

  vector<pair<double,double> >  returnValMinMax;
  returnValMinMax.push_back(xValMinMax);
  returnValMinMax.push_back(yValMinMax);
  
  return returnValMinMax;
                                        
}

void draw_mSUGRA_exclusion(TH2F *ocrosssection, TH2F *oFilterEfficiency, TH2F *oabsXS, TH2F *limitmap, TH2F *expmap, TH2F *expplusmap, TH2F *expminusmap, TH2F *exp2plusmap, TH2F *exp2minusmap, bool isobserved) {
  TH2F *crosssection = (TH2F*)ocrosssection->Clone("crosssection");
//  TH2F *limitmap = (TH2F*)olimitmap->Clone(((string)olimitmap->GetName()+"clone").c_str());
  TH2F *cleanhisto = (TH2F*)limitmap->Clone("clean");
  for(int ix=1;ix<=cleanhisto->GetNbinsX();ix++) {
    for(int iy=1;iy<=cleanhisto->GetNbinsY();iy++) {
      cleanhisto->SetBinContent(ix,iy,0);
    }
  }
  
  
  TH2F *FilterEfficiency;
  TH2F *absXS;

  
  
  
  write_warning(__FUNCTION__,"You'll want to switch off 'wrongwaytodothis')");
  
  if(wrongwaytodothis) {
    //this part is the one you want to remove.
    TFile *Efficiencies = new TFile("FilterEfficiencyv3.root");
    FilterEfficiency = cast_into_shape((TH2F*) Efficiencies->Get("FilterEfficiency"),limitmap);
    assert(FilterEfficiency);
    assert(crosssection);
    absXS=(TH2F*)crosssection->Clone("absXS");
    crosssection->Multiply(FilterEfficiency);
  } else {
    //this part is the one you want to keep!
    FilterEfficiency=(TH2F*)oFilterEfficiency->Clone("FilterEfficiency");
    absXS=(TH2F*)oabsXS->Clone("absXS");
  }
    
  

  TH2F *limits = (TH2F*)limitmap->Clone("limits");
  set_range(limits,true,false);
  limitmap->Divide(crosssection);
  expminusmap->Divide(crosssection);
  expplusmap->Divide(crosssection);
  exp2minusmap->Divide(crosssection);
  exp2plusmap->Divide(crosssection);
  expmap->Divide(crosssection);
  TGraph *observed = get_mSUGRA_exclusion_line(limitmap, PlottingSetup::mSUGRA);
  observed->SetLineColor(kRed);
  TGraph *expminus = get_mSUGRA_exclusion_line(expminusmap, PlottingSetup::mSUGRA);
  TGraph *expplus  = get_mSUGRA_exclusion_line(expplusmap, PlottingSetup::mSUGRA);
  TGraph *exp2minus;
  if(draw2sigma) exp2minus = get_mSUGRA_exclusion_line(exp2minusmap, PlottingSetup::mSUGRA);
  TGraph *exp2plus;
  if(draw2sigma) exp2plus = get_mSUGRA_exclusion_line(exp2plusmap, PlottingSetup::mSUGRA);

  TGraph *expected = new TGraph(expminus->GetN()+expplus->GetN());
  TGraph *expected2; 
  if(draw2sigma) expected2 = new TGraph(exp2minus->GetN()+exp2plus->GetN());
  for(int i=0;i<=expminus->GetN();i++) {
    Double_t x,y;
    expminus->GetPoint(i,x,y);
    expected->SetPoint(i,x,y);
  }
    
  for(int i=0;i<=exp2minus->GetN();i++) {
    Double_t x,y;
    exp2minus->GetPoint(i,x,y);
    expected2->SetPoint(i,x,y);
  }
  for(int i=exp2plus->GetN()-1;i>=0;i--) {
    Double_t x,y;
    exp2plus->GetPoint(i,x,y);
    expected2->SetPoint(exp2minus->GetN()+(exp2plus->GetN()-i),x,y);
  }
  for(int i=expplus->GetN()-1;i>=0;i--) {
    Double_t x,y;
    expplus->GetPoint(i,x,y);
    expected->SetPoint(expminus->GetN()+(expplus->GetN()-i),x,y);
  }
  expected->SetFillColor(TColor::GetColor("#9FF781"));
  if(draw2sigma) expected2->SetFillColor(TColor::GetColor("#F3F781"));
  
  smooth_line(observed);
  smooth_line(expected);
  if(draw2sigma) smooth_line(expected2);
  
  TCanvas *te = new TCanvas("te","te");
  te->SetRightMargin(standardmargin);
//  decorate_mSUGRA(cleanhisto,te,expected,expected2,observed);
  TH2F *noh = new TH2F("noh","noh",1,1,2,1,1,2);
  SugarCoatThis(te,10,noh,observed);
//  expected->Draw("c");
//  observed->Draw("c");
  stringstream saveas;
  if((int)((string)limitmap->GetName()).find("limitmap")>0) {
    saveas << "Limits/";
    if(!isobserved) saveas << "expected/expected_";
    saveas << "final_exclusion_for_JZB_geq_" << ((string)limitmap->GetName()).substr(((string)limitmap->GetName()).find("limitmap")+8,10);
  } else {
    saveas << "Limits/";
    if(!isobserved) saveas << "expected/expected";
    saveas << "final_exclusion_for_bestlimits";
  }
//.........这里部分代码省略.........

void slopetest(){

   string dist = "maos210";
   double nevts = 2*50000.0;
   string dim = "mt";

   double p1=0, p2=0, p3=0;
   double dx = 0;
   double eps = 0.2;
   if( dim == "mt" ){
      p1 = 172.0; p2 = 172.5, p3 = 173.0;
      //p1 = 172.5; p2 = 173, p3 = 173.5;
      //p1 = 166.0; p2 = 166.5, p3 = 167.0;
      //p1 = 178.0; p2 = 178.5, p3 = 179.0;
   }
   if( dim == "jes" ){
      p1 = 0.999; p2 = 1.000, p3 = 1.001;
   }
   if( dist == "mbl" ){
      dx = 2.8;
   }
   if( dist == "mt2_221" ){
      dx = 1.9;
   }

   TFile *f = new TFile("../results/plotsTemplates.root");
   TDirectory *d = (TDirectory*)f->Get((dim+"shape_"+dist+"_gp").c_str());

   /*
   TCanvas *ctemp = (TCanvas*)f->Get(("c_"+dist+"_gp_signal_JSF1000").c_str());
   TH1 *h = (TH1*)(ctemp->GetListOfPrimitives()->At(1));
   dx = h->GetBinWidth(1);

   delete ctemp;
   */
   dx = 10;

   TGraph *gslope = new TGraph();
   TGraph *gint = new TGraph();

   double sig = 0;
   double d2sig = 0;

   TIter nextkey(d->GetListOfKeys());
   TKey *key;
   while( (key = (TKey*)nextkey()) ){

      string name = key->GetName();
      string classname = key->GetClassName();

      TCanvas *c = (TCanvas*)d->Get(name.c_str());
      TGraph *g = (TGraph*)(c->GetListOfPrimitives()->At(1));

      double y1=0, y2=0;
      double ycent = 0;
      for(int i=0; i < g->GetN(); i++){
         double x=0, y=0;
         g->GetPoint(i, x, y);
         if( x == p1 ) y1 = y;
         if( x == p2 ) ycent = y;
         if( x == p3 ) y2 = y;
      }
      double slope = y2-y1;

      string sbin = name;
      string ntemp = "c"+dist+"_gp";
      sbin.erase(sbin.begin(), sbin.begin()+ntemp.length());
      if( dist == "mt2_221" ) sbin.erase(sbin.end(), sbin.end()+1);
      double dbin = atof(sbin.c_str());

      double integrand = slope*slope/ycent;

      gint->SetPoint(gint->GetN(), dbin, integrand);
      gslope->SetPoint(gslope->GetN(), dbin, slope);

      sig += slope*slope*eps*eps*-0.5*dx/ycent;
      d2sig += integrand*dx;

   }

   cout << "Projected S(" << eps << ") = " << sig-2.0*nevts << endl;
   cout << "Projected sigma = " << sqrt(1.0/(d2sig*nevts)) << endl;

   TCanvas *c1 = new TCanvas();
   gslope->SetLineWidth(2);
   gslope->GetXaxis()->SetTitle("(GeV)");
   gslope->GetXaxis()->SetTitleSize(0.05);
   gslope->GetXaxis()->SetTitleOffset(0.8);
   gslope->Draw("AC");

   TCanvas *c2 = new TCanvas();
   gint->SetLineWidth(2);
   gint->GetXaxis()->SetTitle("(GeV)");
   gint->GetXaxis()->SetTitleSize(0.05);
   gint->GetXaxis()->SetTitleOffset(0.8);
   gint->Draw("AC");

   return;
}

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);
//.........这里部分代码省略.........

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

  //..................................................................................................

  idx = 1;

  double xp1 = 0.0;
  double xp2 = 0.0;
    
  for( int k = 3; k < nGraphs; ++k) 
  {
    
    if ( idx >= 4 ) idx = 1;
    
    TGraph * g1 = (TGraph*)PhiGraphs->At(k);
    
    int maxpts = g1->GetN();
    
    xx = 0.0;
    double y1    = 0.0;
    double y2    = 0.0;
    double ratio = 1.0;
    
    TGraph * r1 = new TGraph();
    int np = 0;
    
    if( idx == 1 ) 
    {
      
      TGraph * denom = (TGraph*)VacuumGraphs->At(0);
      np = 0;
      for( int j = 0; j < maxpts; ++j) 
      {
        
        g1->GetPoint(j, xp1, y1);
        denom->GetPoint(j, xp2, y2);
        if ( y2 <= 1.0e-3 ) continue;
        ratio = y1/y2;

        if (ratio > 5 ) std::cout << " nu_e> x1= " << xp1 << " x2= " << xp2 << " num= " 
                                  << y1 << " den= " << y2 << " num/den= " << ratio << std::endl;

        r1->SetPoint( np, xp1, ratio );
        ++np;
      }
      
      RatioGraphs->Add(r1);
      std::cout << " ratio added " << std::endl;
      
    } else if ( idx == 2 ) 
    {

      TGraph * denom = (TGraph*)VacuumGraphs->At(1);
      np = 0;
      for( int j = 0; j < maxpts; ++j) 
      {

        g1->GetPoint(j, xp1, y1);
        denom->GetPoint(j, xp2, y2);
        if ( y2 <= 1.0e-3 ) continue;
        ratio = y1/y2;
        if (ratio > 5 ) std::cout << " nu_e> x1= " << xp1 << " x2= " << xp2 << " num= " 
                                  << y1 << " den= " << y2 << " num/den= " << ratio << std::endl;
        r1->SetPoint( np, xp1, ratio );
        ++np;
      }
      

TGraph* getContour(TH2D* h, TString name) {

  TGraph* graph = new TGraph(100);
  graph->SetName(name);
  int ip = 0;
  int nx = h->GetXaxis()->GetNbins();
  int ny = h->GetYaxis()->GetNbins();

  // for y>x
  int ix = -1;
  for(int j=ny;true; j--) {
    int k = -1;
    for(int i=2; i<nx-1; i++) {
      if(h->GetBinContent(i,j) < 0) {
	std::cout << "i,j,z : " << i << ", " << j << ", " << h->GetBinContent(i,j) << std::endl;
        k = i;
        break;
      }
    }// for i
    if(k<0) continue;
    double y = h->GetYaxis()->GetBinCenter(j);
    double x1 = h->GetXaxis()->GetBinCenter(k-1);
    double x2 = h->GetXaxis()->GetBinCenter(k);
    double z1 = h->GetBinContent(k-1,j);
    double z2 = h->GetBinContent(k,j);
    double x = x1 + (x2-x1)*fabs(z1)/fabs(z2-z1);
    std::cout << "y, x1, x2, z1, z2, x : " << y << ", " << x1 << ", " << x2 << ", " << x << ", " << z1 << ", " << z2 << std::endl;
    graph->SetPoint(ip++,x,y);

    if(h->GetYaxis()->GetBinCenter(j) < h->GetXaxis()->GetBinCenter(k)) {
      ix = k;
      break;
    }
  }// for j

  if(ix < 0) std::cout << "Something wrong...." << std::endl;

  // for y<x
  for(int i=ix; i<=nx; i++) {
    int k = -1;
    for(int j=2; j<ny-1; j++) {
      if(h->GetBinContent(i,j) < 0) {
        k = j;
        break;
      }
    }// for j
    if(k<0) continue;
    double x = h->GetXaxis()->GetBinCenter(i);
    double y1 = h->GetYaxis()->GetBinCenter(k-1);
    double y2 = h->GetYaxis()->GetBinCenter(k);
    double z1 = h->GetBinContent(i,k-1);
    double z2 = h->GetBinContent(i,k);
    double y = y1 + (y2-y1)*fabs(z1)/fabs(z2-z1);
    std::cout << "x, y1, y2, z1, z2, y : " << x << ", " << y1 << ", " << y2 << ", " << y << ", " << z1 << ", " << z2 << std::endl;
    graph->SetPoint(ip++,x,y);
  }// for i

  ip = graph->GetN()-1;
  while(1) {
    double x, y;
    graph->GetPoint(ip,x,y);
    if(x>1) break;
    else graph->RemovePoint(ip);
    ip--;
  }

  return graph;
}

//.........这里部分代码省略.........
    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;
}

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

  //return;
  // Get Contours
   TObjArray *conts = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
   TList* contLevel = NULL;
   TGraph* curv     = NULL;
   TGraph* gc       = NULL;

   Int_t nGraphs    = 0;
   Int_t TotalConts = 0;

   if (conts == NULL){
      printf("*** No Contours Were Extracted!\n");
      TotalConts = 0;
      return;
   } else {
      TotalConts = conts->GetSize();
   }

   printf("TotalConts = %d\n", TotalConts);

   for(i = 0; i < TotalConts; i++){
      contLevel = (TList*)conts->At(i);
      printf("Contour %d has %d Graphs\n", i, contLevel->GetSize());
      nGraphs += contLevel->GetSize();
   }

   nGraphs = 0;

   TH2F *hr = new TH2F("hr", ";m_{A};tan#beta", 2, 225, 600, 2, 0.1, 100);

   hr->GetXaxis()->SetTitleOffset(1.1);
   hr->GetXaxis()->SetRangeUser(200,650);
   hr->GetYaxis()->SetTitleOffset(1.2);
   hr->GetYaxis()->SetNdivisions(110,kFALSE);
   hr->GetXaxis()->SetNdivisions(20205,kFALSE);
   hr->Draw();
   Double_t x0, y0, z0;
   TLatex l;
   l.SetTextSize(0.03);
   l.SetTextAlign(32);
   char val[20];

   for(i = 0; i < TotalConts; i++){
      contLevel = (TList*)conts->At(i);
      z0 = level[i];
      printf("Z-Level Passed in as:  Z = %f\n", z0);

      // Get first graph from list on curves on this level
      curv = (TGraph*)contLevel->First();
      for(j = 0; j < contLevel->GetSize(); j++){
        // last point
         //curv->GetPoint(curv->GetN()-1, x0, y0);
        // first point
         curv->GetPoint(2, x0, y0);

         // if (z0<0) curv->SetLineColor(kRed);
         // if (z0>0) curv->SetLineColor(kBlue);
         nGraphs ++;
         printf("\tGraph: %d  -- %d Elements\n", nGraphs,curv->GetN());

	 // Draw clones of the graphs to avoid deletions in case the 1st
	 // pad is redrawn.
         gc = (TGraph*)curv->Clone();
         gc->Draw("C");

         if (z0>=.01) sprintf(val,"%0.2f",z0);
         if (z0>=.1) sprintf(val,"%0.2f",z0);
         if (z0>=1) sprintf(val,"%0.0f",z0);
         l.DrawLatex(x0*0.99,y0,val);
         curv = (TGraph*)contLevel->After(curv); // Get Next graph
      }
   }
   gPad->SetLogy();
   gPad->SetGridx();
   gPad->SetGridy();
   gPad->SetRightMargin(0.05);
   gPad->SetTopMargin(0.10);
   c2->Update();
   printf("\n\n\tExtracted %d Contours and %d Graphs \n", TotalConts, nGraphs );



  tl.SetTextAlign(31);
  tl.DrawLatex(0.8,0.85,goodName);

  tl.SetTextAlign(31);
  tl.DrawLatex(0.8,0.77,goodType);


  pCan(c2,cname+"_BW");

   c1->cd(1);
   gPad->SetLogy();
   gPad->SetLogz();
   c1->cd(2);
   gPad->SetLogy();
   c1->Update();
  pCan(c1,cname);
}

//.........这里部分代码省略.........
   TH2D *HistStreamFn = new TH2D("HstreamFn",
   "#splitline{Histogram with negative and positive contents. Six contours are defined.}{It is plotted with options CONT LIST to retrieve the contours points in TGraphs}",
   nZsamples, z[0], z[nZsamples-1], nPhiSamples, phi[0], phi[nPhiSamples-1]);

   // Load Histogram Data
   for (Int_t i = 0; i < nZsamples; i++) {
      for(Int_t j = 0; j < nPhiSamples; j++){
         HistStreamFn->SetBinContent(i,j, HofZ[i]*FofPhi[j]);
      }
   }

   gStyle->SetPalette(1);
   gStyle->SetOptStat(0);
   gStyle->SetTitleW(0.99);
   gStyle->SetTitleH(0.08);

   Double_t contours[6];
   contours[0] = -0.7;
   contours[1] = -0.5;
   contours[2] = -0.1;
   contours[3] =  0.1;
   contours[4] =  0.4;
   contours[5] =  0.8;

   HistStreamFn->SetContour(6, contours);

   // Draw contours as filled regions, and Save points
   HistStreamFn->Draw("CONT Z LIST");
   c->Update(); // Needed to force the plotting and retrieve the contours in TGraphs

   // Get Contours
   TObjArray *conts = (TObjArray*)gROOT->GetListOfSpecials()->FindObject("contours");
   TList* contLevel = NULL;
   TGraph* curv     = NULL;
   TGraph* gc       = NULL;

   Int_t nGraphs    = 0;
   Int_t TotalConts = 0;

   if (conts == NULL){
      printf("*** No Contours Were Extracted!\n");
      TotalConts = 0;
      return;
   } else {
      TotalConts = conts->GetSize();
   }

   printf("TotalConts = %d\n", TotalConts);

   for(i = 0; i < TotalConts; i++){
      contLevel = (TList*)conts->At(i);
      printf("Contour %d has %d Graphs\n", i, contLevel->GetSize());
      nGraphs += contLevel->GetSize();
   }

   nGraphs = 0;

   TCanvas* c1 = new TCanvas("c1","Contour List",610,0,600,600);
   c1->SetTopMargin(0.15);
   TH2F *hr = new TH2F("hr",
   "#splitline{Negative contours are returned first (highest to lowest). Positive contours are returned from}{lowest to highest. On this plot Negative contours are drawn in red and positive contours in blue.}",
   2, -2, 2, 2, 0, 6.5);

   hr->Draw();
   Double_t x0, y0, z0;
   TLatex l;
   l.SetTextSize(0.03);
   char val[20];

   for(i = 0; i < TotalConts; i++){
      contLevel = (TList*)conts->At(i);
      if (i<3) z0 = contours[2-i];
      else     z0 = contours[i];
      printf("Z-Level Passed in as:  Z = %f\n", z0);

      // Get first graph from list on curves on this level
      curv = (TGraph*)contLevel->First();
      for(j = 0; j < contLevel->GetSize(); j++){
         curv->GetPoint(0, x0, y0);
         if (z0<0) curv->SetLineColor(kRed);
         if (z0>0) curv->SetLineColor(kBlue);
         nGraphs ++;
         printf("\tGraph: %d  -- %d Elements\n", nGraphs,curv->GetN());

	 // Draw clones of the graphs to avoid deletions in case the 1st
	 // pad is redrawn.
         gc = (TGraph*)curv->Clone();
         gc->Draw("C");

         sprintf(val,"%g",z0);
         l.DrawLatex(x0,y0,val);
         curv = (TGraph*)contLevel->After(curv); // Get Next graph
      }
   }
   c1->Update();
   printf("\n\n\tExtracted %d Contours and %d Graphs \n", TotalConts, nGraphs );
   gStyle->SetTitleW(0.);
   gStyle->SetTitleH(0.);
   return c1;
}