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

//.........这里部分代码省略.........
    if(Buffer[n]=='p') Buffer[n] = '+';
    if(Buffer[n]=='m') Buffer[n] = '-';
  }
  if(!strcmp(Iso, "32"))  sprintf(Buffer, "%s^{3/2}",  Buffer, Iso);
  if(!strcmp(Iso, "p12")) sprintf(Buffer, "%s^{p1/2}", Buffer, Iso);
  PlotsRe->GetYaxis()->SetTitle(Buffer);

  //Open text file with model values for given p pi0 multipole
  sprintf(Buffer, "model/ppi0/%s.txt", Mlp);
  InModel_0 = fopen(Buffer, "r");
  //Skip two lines with table header
  fgets(Buffer, sizeof(Buffer), InModel_0);
  fgets(Buffer, sizeof(Buffer), InModel_0);
  //Read multipole model values from file
  MoPts_0 = 0;
  while(!feof(InModel_0))
  {
    if(fscanf(InModel_0, "%lf %lf %lf", &W, &Re, &Im)==3)
    {
      MoW_0[MoPts_0]  = W;
      MoRe_0[MoPts_0] = Re;
      MoIm_0[MoPts_0] = Im;
      MoPts_0++;
    }
  }
  //Close file with model values
  fclose(InModel_0);

  //Open text file with model values for given n pi+ multipole
  sprintf(Buffer, "model/npip/%s.txt", Mlp);
  InModel_p = fopen(Buffer, "r");
  //Skip two lines with table header
  fgets(Buffer, sizeof(Buffer), InModel_p);
  fgets(Buffer, sizeof(Buffer), InModel_p);
  //Read multipole model values from file
  MoPts_p = 0;
  while(!feof(InModel_p))
  {
    if(fscanf(InModel_p, "%lf %lf %lf", &W, &Re, &Im)==3)
    {
      MoW_p[MoPts_p]  = W;
      MoRe_p[MoPts_p] = Re;
      MoIm_p[MoPts_p] = Im;
      MoPts_p++;
    }
  }
  //Close file with model values
  fclose(InModel_p);

  //Create selected isospin multipole from p pi0 and n pi+ multipoles
  for(Int_t wp=0; wp<MoPts_p; wp++)
  {
    //Find corresponding energy bin between n pi+ and p pi0 multipoles
    Int_t w0;
    for(Int_t w0=0; w0<MoPts_0; w0++)
      if(MoW_p[wp]==MoW_0[w0]) break;
    //Create isospin multipoles
    if(!strcmp(Iso, "32"))  { MoRe[wp] = A_32(MoRe_0[w0], MoRe_p[wp]); MoIm[wp] = A_32(MoIm_0[w0], MoIm_p[wp]); }
    if(!strcmp(Iso, "p12")) { MoRe[wp] = A_12(MoRe_0[w0], MoRe_p[wp]); MoIm[wp] = A_12(MoIm_0[w0], MoIm_p[wp]); }
  }

  //Create graphs for real and imaginary parts of model multipole
  ModelRe = new TGraph(MoPts_p, MoW_p, MoRe);
  ModelIm = new TGraph(MoPts_p, MoW_p, MoIm);

  //Color, line size adjustments
  ModelRe->SetLineColor(kRed);
  ModelIm->SetLineColor(kBlue);
  ModelRe->SetLineWidth(2);
  ModelIm->SetLineWidth(2);

  //Set plot titles and object names
  sprintf(Buffer, "%s_%s", Mlp, Iso);
  ModelRe->SetTitle(Buffer);
  ModelIm->SetTitle(Buffer);
  sprintf(Buffer, "Model_Re%s_%s", Mlp, Iso);
  ModelRe->SetName(Buffer);
  sprintf(Buffer, "Model_Im%s_%s", Mlp, Iso);
  ModelIm->SetName(Buffer);

  //Plot graphs
  ModelRe->Draw("L"); //Plot as line, into same frame
  ModelIm->Draw("L"); //Plot as line, into same frame

  //Adjust drawing ranges for y-axis
  if((Lo==0.0) && (Hi==0.0))
  {
    PlotsRe->SetMinimum(Min*1.05);
    PlotsRe->SetMaximum(Max*1.05);
  }
  else
  {
    PlotsRe->SetMinimum(Lo);
    PlotsRe->SetMaximum(Hi);
  }

  if(SAVE) PlotsRe->SetTitle("");
  if(SAVE) sprintf(Buffer, "isospin/%s_%s.eps", Mlp, Iso);
  if(SAVE) Canvas->SaveAs(Buffer);
}

void plotPerCategory(){
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);
  //   setTDRStyle();
   gStyle->SetPadTopMargin   (0.04);
   gStyle->SetPadBottomMargin(0.12);
   gStyle->SetPadRightMargin (0.05);
   gStyle->SetPadLeftMargin  (0.12);
   gStyle->SetTitleSize(0.04, "XYZ");
   gStyle->SetTitleXOffset(1.1);
   gStyle->SetTitleYOffset(1.45);
   gStyle->SetPalette(1);
   gStyle->SetNdivisions(505);

   TCanvas* c1;
   TH1F* framework;
   TH2F* framework2d;
   TLegend* LEG, *LEGTH;
   TGraph* Ref;

   string Directories[]={"cards_SB13TeV_cp1.00_brn0.00", "cards_SB13TeV_GGF_cp1.00_brn0.00", "cards_SB13TeV_VBF_cp1.00_brn0.00"};
   for(unsigned int D=0;D<sizeof(Directories)/sizeof(string);D++){
      string Dir = Directories[D];

     string prod = "pp";
     if(Dir.find("GGF")!=std::string::npos)prod="gg";
     if(Dir.find("VBF")!=std::string::npos)prod="qq";
     bool strengthLimit = false;
     if(prod=="pp")strengthLimit=true;

      c1 = new TCanvas("c", "c",600,600);
      c1->SetLogy(true);
      framework = new TH1F("Graph","Graph",1,190,1510);
      framework->SetStats(false);
      framework->SetTitle("");
      framework->GetXaxis()->SetTitle("Higgs boson mass [GeV]");
      if(strengthLimit){
         framework->GetYaxis()->SetTitle("#mu = #sigma_{95%} / #sigma_{th}");
         framework->GetYaxis()->SetRangeUser(1E-2,1E3);
      }else{
         framework->GetYaxis()->SetTitle((string("#sigma_{95%} (") + prod + " #rightarrow H #rightarrow ZZ) (fb)").c_str());        
         framework->GetYaxis()->SetRangeUser(1E1,1E5);
      }
      framework->GetYaxis()->SetTitleOffset(1.40);
      framework->Draw();

      LEG = new TLegend(0.70,0.70,0.95,0.94);
      LEG->SetFillStyle(0);
      LEG->SetBorderSize(0);
      LEG->SetHeader("Observed:");

      TLegend* LEGExp = NULL;
      LEGExp = new TLegend(0.45,0.70,0.70,0.94);
      LEGExp->SetFillStyle(0);
      LEGExp->SetBorderSize(0);
      LEGExp->SetHeader("Expected:");


      getGraph("=0 Jet"                       , 2, 2, 1, LEG  , NULL, 2, Dir+ "_eq0jets/Stength_LimitSummary")->Draw("C same");
      getGraph("#geq1 Jets"                   , 4, 2, 1, LEG  , NULL, 2, Dir+"_geq1jets/Stength_LimitSummary")->Draw("C same");
      getGraph("VBF"                          , 6, 2, 1, LEG  , NULL, 2, Dir+     "_vbf/Stength_LimitSummary")->Draw("C same");
      getGraph("Combined"                     , 1, 2, 1, LEG  , NULL, 2, Dir+         "/Stength_LimitSummary")->Draw("C same");

      getGraph("=0 Jet"                       , 2, 2, 2, LEGExp  , NULL, 1, Dir+ "_eq0jets/Stength_LimitSummary")->Draw("C same");
      getGraph("#geq1 Jets"                   , 4, 2, 2, LEGExp  , NULL, 1, Dir+"_geq1jets/Stength_LimitSummary")->Draw("C same");
      getGraph("VBF"                          , 6, 2, 2, LEGExp  , NULL, 1, Dir+     "_vbf/Stength_LimitSummary")->Draw("C same");
      getGraph("Combined"                     , 1, 2, 2, LEGExp  , NULL, 1, Dir+         "/Stength_LimitSummary")->Draw("C same");


   //   LEGTH->Draw("same");
      LEGExp  ->Draw("same");
      LEG  ->Draw("same");


      /*
      char LumiLabel[1024];
      sprintf(LumiLabel,"CMS preliminary,  #sqrt{s}=%.0f TeV #scale[0.5]{#int} L=%6.1ffb^{-1}",13.0,2.3);
      TPaveText *pave = new TPaveText(0.1,0.96,0.94,0.99,"NDC");
      pave->SetBorderSize(0);
      pave->SetFillStyle(0);
      pave->SetTextAlign(32);
      pave->SetTextFont(42);
      pave->AddText(LumiLabel);
      pave->Draw("same");*/

      utils::root::DrawPreliminary(2268.759, 13, gPad->GetLeftMargin(),gPad->GetBottomMargin(),gPad->GetRightMargin(),gPad->GetTopMargin()+0.025);     
      if(strengthLimit){
         TLine* SMLine = new TLine(framework->GetXaxis()->GetXmin(),1.0,framework->GetXaxis()->GetXmax(),1.0);
         SMLine->SetLineWidth(2); SMLine->SetLineStyle(1); SMLine->SetLineColor(4);      
         SMLine->Draw("same C");
      }else{
         TGraph* THXSec   = Hxswg::utils::getXSec(Dir); 
         THXSec->SetLineWidth(2); THXSec->SetLineStyle(1); THXSec->SetLineColor(4);
         scaleGraph(THXSec, 1000);  //convert cross-section to fb
         THXSec->Draw("same C");
      }

      c1->SaveAs((Dir+"/perCat_FinalPlot.png").c_str());
      c1->SaveAs((Dir+"/perCat_FinalPlot.pdf").c_str());
      c1->SaveAs((Dir+"/perCat_FinalPlot.C"  ).c_str());
//.........这里部分代码省略.........

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

//.........这里部分代码省略.........
    //Branches
    double xx = 0.0;
    double yy = 0.0;
    
    PxxTreeNu->SetBranchAddress("Xx",&xx);
    PxxTreeNu->SetBranchAddress("Ratio",&yy);
    
    Long64_t nentries = PxxTreeNu->GetEntries();
    
    TGraph * g1 = new TGraph();
    
    for (Long64_t i=0;i<nentries;i++) {
      PxxTreeNu->GetEntry(i);
      g1->SetPoint( i, xx, yy);
    }
    
    v_Graphs->Add( g1 );
      
  }
  
  TString cname = TString("Ratio") + TString("_") + TString(model) +  TString("_") + TString(config);
  
  TCanvas * c1 = new TCanvas( cname.Data(), "track/shower ratio", 206,141,722,575); 

  c1->SetBorderSize(2);
    
  TLegend * leg = new TLegend(0.18,0.64,0.44,0.87);
  
  leg->SetBorderSize(0);
  leg->SetTextFont(22);
  leg->SetTextSize(0.062);
  leg->SetLineColor(1);
  leg->SetLineStyle(1);
  leg->SetLineWidth(1);
  leg->SetFillColor(0);
  leg->SetFillStyle(0);

  int labelpos = 0;
  
  for( int k = 0; k < max; ++k ) 
  {

    TGraph * gg = (TGraph*)v_Graphs->At(k);

    gg->SetMarkerStyle(25);
    gg->SetFillColor(10);

    gg->SetLineColor(linecolor[k]);
    gg->SetLineWidth(linewidth[k]);
    gg->SetLineStyle(linestyle[k]);

    gg->SetMaximum(MAXY);
    gg->SetMinimum(1.6);

    gg->GetXaxis()->SetLimits( 0.0, 0.055 );

    gg->GetXaxis()->SetTitle("sin^{2}#theta_{13}");
    gg->GetXaxis()->CenterTitle(true);
    gg->GetXaxis()->SetLabelFont(42);
    gg->GetXaxis()->SetLabelOffset(0.006);
    gg->GetXaxis()->SetLabelSize(0.06);
    gg->GetXaxis()->SetTitleSize(0.06);
    gg->GetXaxis()->SetTickLength(0.05);
    gg->GetXaxis()->SetTitleOffset(1.07);
    gg->GetXaxis()->SetTitleFont(42);
    gg->GetXaxis()->SetNdivisions(509);

void limitSemiLep(bool only_leptonTAG=false, bool save_plots = false, float min=0.006, float max=900., TString plot_name="Limit_prova.pdf") {
  setTDRStyle();
  TCanvas* c1  = new TCanvas("c1","c1",0,0,800,600); //c1->SetLogy();
  c1->SetLogy();
  save_plots = true;
  plot_name="LimitSemiLep.pdf";

  Double_t XMass[val]    = {1000, 1500, 2000, 2500};
	
  //LIMIT - Breit-Wigner
  Double_t Limit_Obs[val]     = {0,0,0,0};
  Double_t Limit_Exp_m2s[val] = {1.8, 0.9, 0.2, 0.3};
  Double_t Limit_Exp_m1s[val] = {2.5, 1.3, 0.4, 0.4};
  Double_t Limit_Exp[val]     = {3.6, 1.9, 0.7, 0.8};
  Double_t Limit_Exp_p1s[val] = {5.3, 2.9, 1.3, 1.5};
  Double_t Limit_Exp_p2s[val] = {7.5, 4.3, 2.1, 2.4};

  Double_t XSEC[val] = {1.54, 0.22, 0.03, 0.01};
	
  //Draw limit
  TLegend *legendr = new TLegend(0.48,0.70,0.87,0.92);
  legendr->SetTextSize(0.030); 
  legendr->SetFillColor(0);
  legendr->SetBorderSize(1);
	
  //2 sigma 
  TGraph *band_2s = new TGraph(2*val+1);
  for(int i=0;i<val;i++){
    band_2s->SetPoint(i,XMass[i],Limit_Exp_p2s[i]);
    band_2s->SetPoint(i+val,XMass[val-1-i],Limit_Exp_m2s[val-1-i]);
  }
  band_2s->SetPoint(2*val,XMass[0],Limit_Exp_p2s[0]);
  band_2s->SetLineStyle(2);
  band_2s->SetFillColor(kYellow);
  band_2s->Draw("ALF2");
  band_2s->SetMinimum(min);
  band_2s->SetMaximum(max);
  band_2s->SetTitle(0);
  band_2s->GetXaxis()->SetTitle("M(Z,H) [GeV]");
  band_2s->GetYaxis()->SetTitle("#sigma(pp #rightarrow Z' #rightarrow ZH) x BR(Z #rightarrow qq) X BR(H #rightarrow #tau#tau) [fb]");
  band_2s->GetXaxis()->SetTitleOffset(1.1);
  band_2s->GetYaxis()->SetTitleOffset(1.3);
  band_2s->GetXaxis()->SetRangeUser(1000,2500);
  band_2s->GetYaxis()->SetTitleSize(0.040);
  band_2s->GetXaxis()->SetNdivisions(1005);
  band_2s->Draw("LF2");
	
  //1 sigma
  TGraph *band_1s = new TGraph(2*val+1);
  for(int i=0;i<val;i++){
    band_1s->SetPoint(i,XMass[i],Limit_Exp_p1s[i]);
    band_1s->SetPoint(i+val,XMass[val-1-i],Limit_Exp_m1s[val-1-i]);
  }
  band_1s->SetPoint(2*val,XMass[0],Limit_Exp_p1s[0]);
  band_1s->SetLineStyle(2);
  band_1s->SetFillColor(kGreen);
  band_1s->Draw("LF2");
  //Expected
  TGraph *limit_exp = new TGraph(val,XMass,Limit_Exp);
  limit_exp->SetLineWidth(2);
  limit_exp->SetLineStyle(2);
  limit_exp->SetMarkerSize(1.3);
  limit_exp->Draw("LP");
  //Observed
  TGraph *limit_obs = new TGraph(val,XMass,Limit_Obs);
  limit_obs->SetMarkerStyle(21);
  //limit_obs->Draw("LP");

  TGraph *xsec = new TGraph(val,XMass,XSEC);
  xsec->SetMarkerStyle(21);
  xsec->SetLineColor(2);
  xsec->SetMarkerColor(2);
  xsec->SetLineWidth(2);
  //xsec->Draw("LP");

  TGraph *band_1s2 = new TGraph(2*val+1);
  TGraph *band_2s2 = new TGraph(2*val+1);
  band_1s2->SetLineStyle(2);
  band_1s2->SetFillColor(kGreen);
  band_1s2->SetLineColor(kBlue);
  band_1s2->SetMarkerColor(kBlue);
  band_2s2->SetLineStyle(2);
  band_2s2->SetFillColor(kYellow);
  band_2s2->SetLineColor(kBlue);
  band_2s2->SetMarkerColor(kBlue);

  //Legend
  //legendr->AddEntry(xsec,"Theoretical cross-section","LP");
  legendr->AddEntry(band_1s, "Expected #pm 1#sigma","flP");
  legendr->AddEntry(band_2s, "Expected #pm 2#sigma","flP");
  //legendr->AddEntry(band_1s2,"Expected #pm 1#sigma - GaussianXLandau","flP");
  //legendr->AddEntry(band_2s2,"Expected #pm 2#sigma - GaussianXLandau","flP");
  legendr->Draw();
	
  //Save plot
  if(save_plots) c1->SaveAs(plot_name); 
}

void graphLE3(Long64_t entry=0,int num=2,int spot=0){
  
  gROOT->ProcessLine(".L ~/analysis/scripts/LoadStyle.C");
  LoadStyle();

  LendaEvent * event = new LendaEvent();
  
  TTree* flt =(TTree*)gDirectory->Get("flt");
  
  flt->SetBranchAddress("Event",&event);

  flt->GetEntry(entry);

  cout<<event->Traces.size()<<endl;


  int size = (int) event->Traces[spot].size();

  Double_t* x = malloc(size*sizeof(Double_t));
  Double_t* y = malloc(size*sizeof(Double_t));

  Double_t* y1= malloc(size*sizeof(Double_t));
    

  Double_t* y2=malloc(size*sizeof(Double_t));
    

    cout<<"size is "<<size<<endl;
    for (int i=0;i<size;i++){

      x[i]=i*10;
      y[i]=event->Traces[spot][i] -380;
      y1[i]=event->Filters[spot][i];
      y2[i]=event->CFDs[spot][i];
    }
    cout<<"50 "<<event->CFDs[spot][50]<<endl;
    cout<<"50 "<<event->CFDs[spot][51]<<endl;
    cout<<"50 "<<event->CFDs[spot][52]<<endl;
    cout<<"50 "<<event->CFDs[spot][53]<<endl;



    TGraph *gr = new TGraph(size,x,y);
    TGraph *gr1 = new TGraph(size,x,y1);
    TGraph *gr2 = new TGraph(size,x,y2);

    TCanvas *c = new TCanvas();

    c->cd(1);

    gr->SetFillColor(kBlack);
    gr1->SetFillColor(0);
    gr2->SetFillColor(0);

    gr->SetLineColor(kBlack);
    gr1->SetLineColor(kBlue);
    gr2->SetLineColor(kRed);
    gr->SetLineWidth(3);
    gr2->SetLineWidth(3);
    gr1->SetLineWidth(3);

    gr->SetMarkerSize(.7);
    gr1->SetMarkerSize(.7);
    gr2->SetMarkerSize(.7);


    gr->SetMarkerStyle(8);
    gr1->SetMarkerStyle(8);
    gr2->SetMarkerStyle(8);
    TMultiGraph *mg = new TMultiGraph();

    leg = new TLegend(0.7,0.7,1,1);
    leg->SetLineWidth(0);
    leg->SetHeader("");
    leg->AddEntry(gr, "Trace","l");
    leg->AddEntry(gr1,"Fast Filter","l");
    leg->AddEntry(gr2,"CFD Filter","l");

    leg->SetTextSize(.05);
    leg->SetFillColor(kWhite);
    mg->Add(gr);
    mg->Add(gr1);
    mg->Add(gr2);

    mg->SetTitle("Pixie Digital Waveform");
    mg->Draw("a L P");
    mg->GetHistogram()->GetXaxis()->SetTitle("Time [ns]");
    mg->GetHistogram()->GetYaxis()->SetTitle("ADC Channel");
    mg->GetHistogram()->GetYaxis()->SetTitleOffset(1.3);

    leg->Draw();      
    
}

void DMplot()
{

    TCanvas *canv = new TCanvas("canv", "limits canvas", 800., 680.);
	gStyle->SetCanvasDefH(600); //Height of canvas
	gStyle->SetCanvasDefW(640); //Width of canvas
	gStyle->SetCanvasDefX(0);   //POsition on screen
	gStyle->SetCanvasDefY(0);

	gStyle->SetPadLeftMargin(0.14);//0.16);
	gStyle->SetPadRightMargin(0.165);//0.02);
	gStyle->SetPadTopMargin(0.085);//0.02);
	gStyle->SetPadBottomMargin(0.12);//0.02);

	  // For g axis titles:
	gStyle->SetTitleColor(1, "XYZ");
	gStyle->SetTitleFont(42, "XYZ");
	gStyle->SetTitleSize(0.045, "Z");
	gStyle->SetTitleSize(0.055, "XY");
	gStyle->SetTitleXOffset(1.0);//0.9);
	gStyle->SetTitleYOffset(1.15); // => 1.15 if exponents

	// For g axis labels:
	gStyle->SetLabelColor(1, "XYZ");
	gStyle->SetLabelFont(42, "XYZ");
	gStyle->SetLabelOffset(0.007, "XYZ");
	gStyle->SetLabelSize(0.04, "XYZ");

	// Legends
	gStyle->SetLegendBorderSize(0);
	gStyle->SetLegendFillColor(kWhite);
	gStyle->SetLegendFont(42);
    TPad* t1d = new TPad();
    t1d = new TPad("t1d","t1d", 0.0, 0.0, 1.0, 1.0);
    t1d->Draw();
    t1d->SetTicky();
    t1d->SetTickx();
    t1d->SetRightMargin(0.03);
    t1d->cd();

    //t1d->SetGridx(1);
    //t1d->SetGridy(1);
    t1d->SetLogy();
    t1d->SetLogx();

    //double const fn = 0.629;//0.326;
    TGraph *h_scalar_min = MakeGraph(0,0.260);
    TGraph *h_scalar_lat = MakeGraph(0,0.326);
    TGraph *h_scalar_max = MakeGraph(0,0.629);
    TGraph *h_fermion_min = MakeGraph(1,0.260);
    TGraph *h_fermion_lat = MakeGraph(1,0.326);
    TGraph *h_fermion_max = MakeGraph(1,0.629);

//////////////////////////////////////////////////

    // Get LUX bound 
    //TGraph *z1 = new TGraph("LUX_90CL.dat","%lg %lg");
    //
    TFile *fi_LUX2015 = TFile::Open("LUX_latest_2015.root");
    TGraph *z12015 =(TGraph*) fi_LUX2015->Get("LUX_2015");
    //TFile *fi_LUX2016 = TFile::Open("LUX_latest_2016.root");
    //TGraph *z12016 =(TGraph*) fi_LUX2016->Get("LUX_2016");
    TFile *fi_LUX2016 = TFile::Open("LUX_2013_2014_2015_combined.root");
    TGraph *z12016 =(TGraph*) fi_LUX2016->Get("LUX_2013_2014_2015_combined");

    z12015->SetLineWidth(3);
    z12015->SetLineColor(kGreen+2);
//    z12015->SetLineStyle(2);
    z12016->SetLineWidth(3);
    //z12016->SetLineStyle(2);
    z12016->SetLineColor(kGreen+2);

//    TLegend *leg2 = new TLegend(0.65, 0.15, 0.93, 0.42);
    TLegend *leg2 = new TLegend(0.70, 0.15, 0.97, 0.42);
    leg2->SetFillColor(0);
    leg2->SetBorderSize(0);
    //leg2->AddEntry(z1,"LUX(90\%CL)","L");
//////////////////////////////////////////////////
//  CDMS/ CRESST TOO
    TFile *CDMS_2016f = TFile::Open("CDMS_2016.root");
    TGraph *SCDMS = (TGraph*)CDMS_2016f->Get("CDMS_2016");

    TFile *CRESST2_2016f = TFile::Open("CRESST_2.root");
    TGraph *CRESST2 = (TGraph*)CRESST2_2016f->Get("CRESST_2_2016");

    TFile *PANDAX_2016f = TFile::Open("PANDAX.root");
    TGraph *PANDAX = (TGraph*)PANDAX_2016f->Get("PANDAX");

    SCDMS->SetLineColor(kAzure+7); SCDMS->SetLineStyle(9); SCDMS->SetLineWidth(3);
    CRESST2->SetLineColor(kMagenta+2); CRESST2->SetLineStyle(7); CRESST2->SetLineWidth(3);
    PANDAX->SetLineColor(kMagenta+2); PANDAX->SetLineStyle(4); PANDAX->SetLineWidth(3);
	
//
//

    h_scalar_min->SetTitle("");
    h_scalar_min->SetMinimum(2.0e-47);
    h_scalar_min->SetMaximum(1.0e-39);
    //h_scalar_min->SetMinimum(0.6e-46);
    //h_scalar_min->SetMaximum(1.0e-42);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
		{
			g->GetPoint(i, pointx0, pointy0);
			gNew->SetPoint(i+2, pointx0, pointy0);
		}

		// add origin
		gNew->SetPoint(0, hCL->GetXaxis()->GetXmin(), 0);

		// add a point at first y height but at x=origin.
		g->GetPoint(0, pointx0, pointy0);
		gNew->SetPoint(1, hCL->GetXaxis()->GetXmin(), pointy0);

		// add a point at last y height but at x=xmax.
		g->GetPoint(g->GetN()-1, pointx0, pointy0);
		gNew->SetPoint(gNew->GetN()-2, hCL->GetXaxis()->GetXmax(), pointy0);

		// add a point at xmax, 0
		gNew->SetPoint(gNew->GetN()-1, hCL->GetXaxis()->GetXmax(), 0);
		g = gNew;
	}

	int color = s->getLineColor();
	if(CLsType>0 && s->getMethodName().Contains("Plugin") && !arg->plotpluginonly) {
    if (CLsType==1) color = kBlue-7;
    else if (CLsType==2) color = kBlue+2;
  }
	else if(CLsType>0) {
    if (CLsType==1) color = s->getLineColor() - 5;
    if (CLsType==2) color = s->getLineColor() - 4;
  }
	g->SetLineColor(color);

	if ( filled ){
		g->SetLineWidth(2);
    double alpha = arg->isQuickhack(12) ? 0.4 : 1.;
    if ( arg->isQuickhack(24) ) alpha = 0.;
		g->SetFillColorAlpha(color,alpha);
		g->SetLineStyle(1);
    g->SetFillStyle( s->getFillStyle() );
	}
	else{
		g->SetLineWidth(2);
		g->SetLineStyle(s->getLineStyle());
    if ( last && arg->isQuickhack(25) ) g->SetLineWidth(3);
	}

	if ( plotPoints ){
		g->SetLineWidth(1);
		g->SetMarkerColor(color);
		g->SetMarkerStyle(8);
		g->SetMarkerSize(0.6);
		if(CLsType==1) {
			g->SetMarkerStyle(33);
			g->SetMarkerSize(1);
		}
		if(CLsType==2) {
			g->SetMarkerStyle(21);
		}
	}

	// build a histogram which holds the axes
	float min = arg->scanrangeMin == arg->scanrangeMax ? hCL->GetXaxis()->GetXmin() : arg->scanrangeMin;
	float max = arg->scanrangeMin == arg->scanrangeMax ? hCL->GetXaxis()->GetXmax() : arg->scanrangeMax;
	TH1F *haxes = new TH1F("haxes"+getUniqueRootName(), "", 100, min, max);
	haxes->SetStats(0);
	haxes->GetXaxis()->SetTitle(s->getScanVar1()->GetTitle());

//.........这里部分代码省略.........
  //scale TH cross-section and limits according to scale factor 
  //this only apply to NarrowResonnance case
  for(unsigned int i=0;i<N;i++){
    if(strengthLimit){
    ObsLimit[i]  /= XSecScaleFactor;
    ExpLimitm2[i]/= XSecScaleFactor;
    ExpLimitm1[i]/= XSecScaleFactor;
    ExpLimit  [i]/= XSecScaleFactor;
    ExpLimitp1[i]/= XSecScaleFactor;
    ExpLimitp2[i]/= XSecScaleFactor;
    }
    ThXSec[i]    *= XSecScaleFactor;
  }

    
  //limits in terms of signal strength
  TCanvas* c = new TCanvas("c", "c",600,600);
  TH1F* framework = new TH1F("Graph","Graph",1,15,70);
  framework->SetStats(false);
  framework->SetTitle("");
  framework->GetXaxis()->SetTitle("A boson mass [GeV]");
  framework->GetYaxis()->SetTitleOffset(1.70);
  if(strengthLimit){
  framework->GetYaxis()->SetTitle("#mu = #sigma_{95%} / #sigma_{th}");
  framework->GetYaxis()->SetRangeUser(5E-2,5E1);
  c->SetLogy(true);
  }else{
  framework->GetYaxis()->SetTitle("#sigma_{95%} (fb)");
  framework->GetYaxis()->SetRangeUser(1E-1,1E3);
  c->SetLogy(true);
  }
  framework->Draw();

  TGraph* TGObsLimit   = new TGraph(N,MassAxis,ObsLimit);  TGObsLimit->SetLineWidth(2);
  TGraph* TGExpLimit   = new TGraph(N,MassAxis,ExpLimit);  TGExpLimit->SetLineWidth(2); TGExpLimit->SetLineStyle(2);
  TCutG* TGExpLimit1S  = GetErrorBand("1S", N, MassAxis, ExpLimitm1, ExpLimitp1);  
  TCutG* TGExpLimit2S  = GetErrorBand("2S", N, MassAxis, ExpLimitm2, ExpLimitp2);  TGExpLimit2S->SetFillColor(5);
  TGraph* THXSec        = new TGraph(N,MassAxis,ThXSec); THXSec->SetLineWidth(2); THXSec->SetLineStyle(1); THXSec->SetLineColor(4);


  TGExpLimit->SetLineColor(2);  TGExpLimit->SetLineStyle(1);
  TGObsLimit->SetLineWidth(2);  TGObsLimit->SetMarkerStyle(20);
  TGExpLimit2S->Draw("fc same");
  TGExpLimit1S->Draw("fc same");
  if(!blind) TGObsLimit->Draw("same CP");
  TGExpLimit->Draw("same C");

  if(strengthLimit){
     TLine* SMLine = new TLine(framework->GetXaxis()->GetXmin(),1.0,framework->GetXaxis()->GetXmax(),1.0);
     SMLine->SetLineWidth(2); SMLine->SetLineStyle(1); SMLine->SetLineColor(4);      
     SMLine->Draw("same C");
  }else{
     THXSec->Draw("same C");
  }


  TPaveText *pave = new TPaveText(0.1,0.96,0.99,0.99,"NDC");
  char LumiLabel[1024];
  if(energy<9){  sprintf(LumiLabel,"CMS preliminary,  #sqrt{s}=%.0f TeV, #scale[0.5]{#int} L=%6.1ffb^{-1} - %20s",energy, luminosity,legendName.Data());
  }else{         sprintf(LumiLabel,"CMS preliminary,  #sqrt{s}=%.0f TeV #scale[0.5]{#int} L=%6.1ffb^{-1}, #sqrt{s}=%.0f TeV #scale[0.5]{#int} L=%6.1ffb^{-1}",7.0,5.0,8.0,19.7);
  }
  pave->SetBorderSize(0);
  pave->SetFillStyle(0);
  pave->SetTextFont(42);
  TObjArray* tokens = (TString(LumiLabel)).Tokenize("\\\\");
  int nt = tokens->GetEntries();

void timeonaxis()
{
// This macro illustrates the use of the time mode on the axis
// with different time intervals and time formats. It's result can
// be seen begin_html <a href="gif/timeonaxis.gif">here</a> end_html
// Through all this script, the time is expressed in UTC. some
// information about this format (and others like GPS) may be found at
// begin_html <a href="http://tycho.usno.navy.mil/systime.html">http://tycho.usno.navy.mil/systime.html</a> end_html
//  or
// begin_html <a href="http://www.topology.org/sci/time.html">http://www.topology.org/sci/time.html</a> end_html
//
// The start time is : almost NOW (the time at which the script is executed)
// actualy, the nearest preceding hour beginning.
// The time is in general expressed in UTC time with the C time() function
// This will obviously most of the time not be the time displayed on your watch
// since it is universal time. See the C time functions for converting this time
// into more useful structures.
//Author: Damir Buskulic

   time_t script_time;
   script_time = time(0);
   script_time = 3600*(int)(script_time/3600);

// The time offset is the one that will be used by all graphs.
// If one changes it, it will be changed even on the graphs already defined
   gStyle->SetTimeOffset(script_time);

   TCanvas *ct = new TCanvas("ct","Time on axis",10,10,700,900);
   ct->Divide(1,3);
   ct->SetFillColor(28);

   int i;

//======= Build a signal : noisy damped sine ======
//        Time interval : 30 minutes

   gStyle->SetTitleH(0.08);
   float noise;
   TH1F *ht = new TH1F("ht","Love at first sight",3000,0.,2000.);
   for (i=1;i<3000;i++) {
      noise = gRandom->Gaus(0,120);
      if (i>700) {
         noise += 1000*sin((i-700)*6.28/30)*exp((double)(700-i)/300);
      }
      ht->SetBinContent(i,noise);
   }
   ct->cd(1);
   ct_1->SetFillColor(41);
   ct_1->SetFrameFillColor(33);
   ht->SetLineColor(2);
   ht->GetXaxis()->SetLabelSize(0.05);
   ht->Draw();
// Sets time on the X axis
// The time used is the one set as time offset added to the value
// of the axis. This is converted into day/month/year hour:min:sec and
// a reasonnable tick interval value is chosen.
   ht->GetXaxis()->SetTimeDisplay(1);

//======= Build a simple graph beginning at a different time ======
//        Time interval : 5 seconds

   float x[100], t[100];
   for (i=0;i<100;i++) {
      x[i] = sin(i*4*3.1415926/50)*exp(-(double)i/20);
      t[i] = 6000+(double)i/20;
   }
   TGraph *gt = new TGraph(100,t,x);
   gt->SetTitle("Politics");
   ct->cd(2);
   ct_2->SetFillColor(41);
   ct_2->SetFrameFillColor(33);
   gt->SetFillColor(19);
   gt->SetLineColor(5);
   gt->SetLineWidth(2);
   gt->Draw("AL");
   gt->GetXaxis()->SetLabelSize(0.05);
// Sets time on the X axis
   gt->GetXaxis()->SetTimeDisplay(1);
   gPad->Modified();

//======= Build a second simple graph for a very long time interval ======
//        Time interval : a few years

   float x2[10], t2[10];
   for (i=0;i<10;i++) {
      x2[i] = gRandom->Gaus(500,100)*i;
      t2[i] = i*365*86400;
   }
   TGraph *gt2 = new TGraph(10,t2,x2);
   gt2->SetTitle("Number of monkeys on the moon");
   ct->cd(3);
   ct_3->SetFillColor(41);
   ct_3->SetFrameFillColor(33);
   gt2->SetFillColor(19);
   gt2->SetMarkerColor(4);
   gt2->SetMarkerStyle(29);
   gt2->SetMarkerSize(1.3);
   gt2->Draw("AP");
   gt2->GetXaxis()->SetLabelSize(0.05);
// Sets time on the X axis
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  gErr2->SetName("gErr2");

  double x,y,yerrUp,yerrDn;
  double xerr = (hExp->GetBinCenter(2)-hExp->GetBinCenter(1))/2.;

  // protect against smoothing over 1
  for (int i=0; i<gExp->GetN(); i++) {
    gExp->GetPoint(i,x,y); gExp->SetPoint(i,x,TMath::Min(y,1.));
    gErr1Up->GetPoint(i,x,y); gErr1Up->SetPoint(i,x,TMath::Min(y,1.));
    gErr1Dn->GetPoint(i,x,y); gErr1Dn->SetPoint(i,x,TMath::Min(y,1.));
    gErr2Up->GetPoint(i,x,y); gErr2Up->SetPoint(i,x,TMath::Min(y,1.));
    gErr2Dn->GetPoint(i,x,y); gErr2Dn->SetPoint(i,x,TMath::Min(y,1.));
  }

  for (int i=0; i<gExp->GetN(); i++) {
    gExp->GetPoint(i,x,y);
    gErr1->SetPoint(i,x,y);
    gErr2->SetPoint(i,x,y);

    gErr1Up->GetPoint(i,x,yerrUp);
    gErr1Dn->GetPoint(i,x,yerrDn);
    gErr1->SetPointError(i, xerr, xerr, y-yerrDn, yerrUp-y );

    gErr2Up->GetPoint(i,x,yerrUp);
    gErr2Dn->GetPoint(i,x,yerrDn);
    gErr2->SetPointError(i, xerr, xerr, y-yerrDn, yerrUp-y );
  }

  gErr2->SetFillColor( TColor::GetColor("#3182bd") );
  gErr2->SetLineColor( TColor::GetColor("#3182bd") );
  gErr1->SetFillColor( TColor::GetColor("#9ecae1") );
  gErr1->SetLineColor( TColor::GetColor("#9ecae1") );
  gExp->SetLineColor(kRed);
  gExp->SetLineWidth(3);
  gObs->SetLineColor(kBlack);
  gObs->SetMarkerColor(kBlack);
  gObs->SetLineWidth(3);
  gObs->SetMarkerSize(1);
  gObs->SetMarkerStyle(20);

	float min = arg->scanrangeMin == arg->scanrangeMax ? hObs->GetXaxis()->GetXmin() : arg->scanrangeMin;
	float max = arg->scanrangeMin == arg->scanrangeMax ? hObs->GetXaxis()->GetXmax() : arg->scanrangeMax;
	TH1F *haxes = new TH1F("haxes"+getUniqueRootName(), "", 100, min, max);
	haxes->SetStats(0);
	haxes->GetXaxis()->SetTitle(s->getScanVar1()->GetTitle());
	haxes->GetYaxis()->SetTitle("CL_{S}");
	haxes->GetXaxis()->SetLabelFont(font);
	haxes->GetYaxis()->SetLabelFont(font);
	haxes->GetXaxis()->SetTitleFont(font);
	haxes->GetYaxis()->SetTitleFont(font);
	haxes->GetXaxis()->SetTitleOffset(0.9);
	haxes->GetYaxis()->SetTitleOffset(0.85);
	haxes->GetXaxis()->SetLabelSize(labelsize);
	haxes->GetYaxis()->SetLabelSize(labelsize);
	haxes->GetXaxis()->SetTitleSize(titlesize);
	haxes->GetYaxis()->SetTitleSize(titlesize);
	int xndiv = arg->ndiv==-1 ? 407 : abs(arg->ndiv);
	bool optimizeNdiv = arg->ndiv<0 ? true : false;
	haxes->GetXaxis()->SetNdivisions(xndiv, optimizeNdiv);
	haxes->GetYaxis()->SetNdivisions(407, true);
  haxes->GetYaxis()->SetRangeUser(0.,1.);

  // Legend:
	// make the legend short, the text will extend over the boundary, but the symbol will be shorter
  float legendXmin = 0.68 ;
  float legendYmin = 0.58 ;

// -------------------------
// -- Functions           --
// -------------------------
void makeDeadEcalEff() {
   
  gStyle->SetOptStat(0);  
  TCanvas *c1 = new TCanvas("c1","A Simple Graph with error bars",200,10,700,500);
  
  const Int_t nbins = 6;

  TString lifetime = "0.5";
  //   TString lifetime = "1";
   //   TString lifetime = "5";

   Float_t masses      [nbins] = {103, 164, 246, 328, 408, 488};  
   Float_t effOld0p5ns [nbins] = {0.884, 0.958, 0.889, 0.895, 0.897, 0.909}; 
   Float_t effDPG0p5ns [nbins] = {0.877, 0.958, 0.894, 0.884, 0.901, 0.888};  // Federico's DPG map 

   Float_t effOld1ns   [nbins] = {0.896, 0.93, 0.912, 0.944, 0.918, 0.914}; 
   Float_t effDPG1ns   [nbins] = {0.899, 0.92, 0.912, 0.926, 0.90, 0.901};  // Federico's DPG map 

   Float_t effOld5ns   [nbins] = {0.921, 0.938, 0.942, 0.92, 0.936, 0.942}; 
   Float_t effDPG5ns   [nbins] = {0.909, 0.913, 0.934, 0.911, 0.93, 0.936};  // Federico's DPG map 

   Float_t* effOld;
   Float_t* effDPG; 

   if (lifetime=="0p5") {
     effOld = effOld0p5ns;
     effDPG = effDPG0p5ns;
   }

   if (lifetime=="1") {
     effOld = effOld1ns;
     effDPG = effDPG1ns;
   }

   if (lifetime=="5") {
     effOld = effOld5ns;
     effDPG = effDPG5ns;
   }

   for (int i=0; i<nbins; i++) {
     effOld[i] = 1.0 - effOld[i];
     effDPG[i] = 1.0 - effDPG[i];
   }

   TGraph *grEffOld = new TGraph(nbins, masses, effOld);  
   TGraph *grEffDPG = new TGraph(nbins, masses, effDPG);  

   grEffOld->SetMarkerStyle(21); // 21:  square  
   grEffOld->SetMarkerSize(0.8);
   grEffOld->SetMarkerColor(kBlue);
   grEffOld->SetLineWidth(2);  
   grEffOld->SetMinimum(0);
   grEffOld->SetMaximum(0.25);
   //  grEffOld->SetTitle(";chargino mass;efficiency of dead ECAL veto"); 
   grEffOld->SetTitle(";chargino mass;1 - #epsilon"); 
   grEffOld->Draw("AP");
   //   grEffOld->Draw("P");

   grEffDPG->SetMarkerStyle(26); // 26:  hollow triangle 
   grEffDPG->SetMarkerSize(0.8);
   grEffDPG->SetMarkerColor(kRed);
   grEffDPG->SetLineWidth(2);  
   grEffDPG->Draw("P, same");

  TLine l;
  l.SetLineColor(kRed);
  l.DrawLine(2.5,1.0,6.5,1.0); 
  TLegend leg2(0.50,0.6,0.85,0.8);
  leg2.AddEntry(grEffOld, "old map ("      + lifetime + " ns)", "pl");  
  leg2.AddEntry(grEffDPG, "ECAL DPG map (" + lifetime + " ns)", "pl");  
  leg2.SetBorderSize(0);
  leg2.SetFillStyle(1001);
  leg2.SetFillColor(kWhite);
  leg2.SetBorderSize(0);
  leg2.SetTextFont(gStyle->GetTitleFont());
  leg2.Draw();


  TPaveText* pt = new TPaveText(0.50, 0.82, 0.90, 0.88, "NDC");
  pt->SetFillStyle(0);
  pt->SetBorderSize(0);
  pt->SetTextFont(gStyle->GetTitleFont());
  pt->AddText("CMS Preliminary, #sqrt{s} = 8 TeV");
  //  pt->Draw();


  c1->SetLogy(0);
  c1->Print("deadEcalEff" + lifetime + "ns.pdf");  
  c1->Clear();

  return;

} 

void mssm_limit_RHW_12p9_Final()
{
//=========Macro generated from canvas: c1/c1
//=========  (Wed Jul 27 09:42:19 2016) by ROOT version6.04/14
    //=========  (Wed Jul 27 09:41:49 2016) by ROOT version6.04/14
    TCanvas *c1 = new TCanvas("c1", "c1",0,45,600,600);
    gStyle->SetOptStat(0);
    c1->SetHighLightColor(2);
    c1->Range(0,0,1,1);
    c1->SetFillColor(0);
    c1->SetBorderMode(0);
    c1->SetBorderSize(2);
    c1->SetLeftMargin(0.12);
    c1->SetRightMargin(0.04);
    c1->SetTopMargin(0.06);
    c1->SetBottomMargin(0.12);
    c1->SetFrameFillStyle(0);
    c1->SetFrameLineWidth(2);
    c1->SetFrameBorderMode(0);
    
    // ------------>Primitives in pad: pad_leg
    TPad *pad_leg = new TPad("pad_leg", "pad_leg",0.4,0.65,0.9,0.9);
    pad_leg->Draw();
    pad_leg->cd();
    pad_leg->Range(0,0,1,1);
    pad_leg->SetFillColor(0);
    pad_leg->SetFillStyle(4000);
    pad_leg->SetBorderMode(0);
    pad_leg->SetBorderSize(2);
    pad_leg->SetLeftMargin(0.12);
    pad_leg->SetRightMargin(0.04);
    pad_leg->SetTopMargin(0.06);
    pad_leg->SetBottomMargin(0.12);
    pad_leg->SetFrameFillStyle(0);
    pad_leg->SetFrameLineWidth(2);
    pad_leg->SetFrameBorderMode(0);
    
    TLegend *leg = new TLegend(0.29,0.01,0.79,0.91,NULL,"NBNDC");
    leg->SetBorderSize(0);
    leg->SetTextFont(62);
    leg->SetTextSize(0.15);
    leg->SetLineColor(1);
    leg->SetLineStyle(1);
    leg->SetLineWidth(1);
    leg->SetFillColor(0);
    leg->SetFillStyle(1001);
    
    //    TLegendEntry *entry=leg->AddEntry("NULL","95% CL Excluded:","h");
    //   entry->SetLineColor(1);
    //   entry->SetLineStyle(1);
    //   entry->SetLineWidth(1);
    //   entry->SetMarkerColor(1);
    //   entry->SetMarkerStyle(21);
    //   entry->SetMarkerSize(1);
    //   entry->SetTextFont(62);
    
    TLegendEntry * entry=leg->AddEntry("Graph1","Theory (M_{Nu}=M_{W}/2)","L");
    entry->SetLineColor(4);
    entry->SetLineStyle(1);
    entry->SetLineWidth(3);
    entry->SetMarkerColor(4);
    entry->SetMarkerStyle(21);
    entry->SetMarkerSize(1);
    entry->SetTextFont(62);
    
    entry=leg->AddEntry("Graph","Observed","L");
    entry->SetLineColor(1);
    entry->SetLineStyle(1);
    entry->SetLineWidth(3);
    entry->SetMarkerColor(1);
    entry->SetMarkerStyle(21);
    entry->SetMarkerSize(1);
    entry->SetTextFont(62);
    
    entry=leg->AddEntry("Graph","Expected","L");
    
    Int_t ci;      // for color index setting
    TColor *color; // for color definition with alpha
    
    ci = TColor::GetColor("#ff0000");
    entry->SetLineColor(ci);
    entry->SetLineStyle(1);
    entry->SetLineWidth(3);
    entry->SetMarkerColor(1);
    entry->SetMarkerStyle(21);
    entry->SetMarkerSize(1);
    entry->SetTextFont(62);
    
    
    entry=leg->AddEntry("","#pm 1#sigma Expected","F");
    
    
    ci = TColor::GetColor("#00ff00");
    entry->SetFillColor(ci);
    entry->SetFillStyle(1001);
    entry->SetLineColor(1);
    entry->SetLineStyle(1);
    entry->SetLineWidth(1);
    entry->SetMarkerColor(1);
    entry->SetMarkerStyle(21);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
       double poiMinus2stdv(-1.) ;
       double poiPlus2stdv(-1.) ;
       double twoDeltalnLMin(1e9) ;

       int nscan(1000) ;
       for ( int xi=0; xi<nscan; xi++ ) {

          double x = poiVals[0] + xi*(poiVals[npoiPoints-1]-poiVals[0])/(nscan-1) ;

          double twoDeltalnL = graph -> Eval( x, 0, "S" ) ;

          if ( poiMinus1stdv < 0. && twoDeltalnL < 1.0 ) { poiMinus1stdv = x ; printf(" set m1 : %d, x=%g, 2dnll=%g\n", xi, x, twoDeltalnL) ;}
          if ( poiMinus2stdv < 0. && twoDeltalnL < 4.0 ) { poiMinus2stdv = x ; printf(" set m2 : %d, x=%g, 2dnll=%g\n", xi, x, twoDeltalnL) ;}
          if ( twoDeltalnL < twoDeltalnLMin ) { poiBest = x ; twoDeltalnLMin = twoDeltalnL ; }
          if ( twoDeltalnLMin < 0.3 && poiPlus1stdv < 0. && twoDeltalnL > 1.0 ) { poiPlus1stdv = x ; printf(" set p1 : %d, x=%g, 2dnll=%g\n", xi, x, twoDeltalnL) ;}
          if ( twoDeltalnLMin < 0.3 && poiPlus2stdv < 0. && twoDeltalnL > 4.0 ) { poiPlus2stdv = x ; printf(" set p2 : %d, x=%g, 2dnll=%g\n", xi, x, twoDeltalnL) ;}

          if ( xi%100 == 0 ) { printf( " %4d : poi=%6.2f,  2DeltalnL = %6.2f\n", xi, x, twoDeltalnL ) ; }

       }
       printf("\n\n POI estimate :  %g  +%g  -%g    [%g,%g],   two sigma errors: +%g  -%g   [%g,%g]\n\n",
               poiBest,
               (poiPlus1stdv-poiBest), (poiBest-poiMinus1stdv), poiMinus1stdv, poiPlus1stdv,
               (poiPlus2stdv-poiBest), (poiBest-poiMinus2stdv), poiMinus2stdv, poiPlus2stdv
               ) ;

       printf(" %s val,pm1sig,pm2sig: %7.2f  %7.2f  %7.2f  %7.2f  %7.2f\n",
          new_poi_name, poiBest, (poiPlus1stdv-poiBest), (poiBest-poiMinus1stdv), (poiPlus2stdv-poiBest), (poiBest-poiMinus2stdv) ) ;

       char htitle[1000] ;
       sprintf(htitle, "%s profile likelihood scan: -2ln(L/Lm)", new_poi_name ) ;
       TH1F* hscan = new TH1F("hscan", htitle, 10, poiMinVal, poiMaxVal ) ;
       hscan->SetMinimum(0.) ;
       hscan->SetMaximum(ymax) ;


       hscan->DrawCopy() ;
       graph->SetLineColor(4) ;
       graph->SetLineWidth(3) ;
       graph->Draw("CP") ;
       gPad->SetGridx(1) ;
       gPad->SetGridy(1) ;
       cscan->Update() ;

       TLine* line = new TLine() ;
       line->SetLineColor(2) ;
       line->DrawLine(poiMinVal, 1., poiPlus1stdv, 1.) ;
       line->DrawLine(poiMinus1stdv,0., poiMinus1stdv, 1.) ;
       line->DrawLine(poiPlus1stdv ,0., poiPlus1stdv , 1.) ;

       TText* text = new TText() ;
       text->SetTextSize(0.04) ;
       char tstring[1000] ;

       sprintf( tstring, "%s = %.1f +%.1f -%.1f", new_poi_name, poiBest, (poiPlus1stdv-poiBest), (poiBest-poiMinus1stdv) ) ;
       text -> DrawTextNDC( 0.15, 0.85, tstring ) ;

       sprintf( tstring, "68%% interval [%.1f,  %.1f]", poiMinus1stdv, poiPlus1stdv ) ;
       text -> DrawTextNDC( 0.15, 0.78, tstring ) ;


       char hname[1000] ;
       sprintf( hname, "hscanout_%s", new_poi_name ) ;
       TH1F* hsout = new TH1F( hname,"scan results",4,0.,4.) ;
       double obsVal(-1.) ;
       hsout->SetBinContent(1, obsVal ) ;
       hsout->SetBinContent(2, poiPlus1stdv ) ;
       hsout->SetBinContent(3, poiBest ) ;
       hsout->SetBinContent(4, poiMinus1stdv ) ;
       TAxis* xaxis = hsout->GetXaxis() ;
       xaxis->SetBinLabel(1,"Observed val.") ;
       xaxis->SetBinLabel(2,"Model+1sd") ;
       xaxis->SetBinLabel(3,"Model") ;
       xaxis->SetBinLabel(4,"Model-1sd") ;

       char outrootfile[10000] ;
       sprintf( outrootfile, "%s/scan-ff-%s.root", outputdir.Data(), new_poi_name ) ;

       char outpdffile[10000] ;
       sprintf( outpdffile, "%s/scan-ff-%s.pdf", outputdir.Data(), new_poi_name ) ;

       cscan->Update() ; cscan->Draw() ;

       printf("\n Saving %s\n", outpdffile ) ;
       cscan->SaveAs( outpdffile ) ;



     //--- save in root file

       printf("\n Saving %s\n", outrootfile ) ;
       TFile fout(outrootfile,"recreate") ;
       graph->Write() ;
       hsout->Write() ;
       fout.Close() ;

       delete ws ;
       wstf->Close() ;

   }

void
draw2DLimitContours(map<string,TList *>& m_contours,
		    const TString& par1,
		    const TString& par2,
		    const TString& plotprefix,
		    TLegend *legend,
		    float par1_bestfit,
		    float par2_bestfit)
{

  //from here we build the two-dimensional aTGC limit

  TCanvas *finalPlot = new TCanvas("final","limits",500,500);
  finalPlot->cd();

  cout << "Drawing expected 68%" << endl;

  TList *contLevel = m_contours["exp68"];
  TGraph *curv;

  std::cout << "m_contours.size() = " << m_contours.size() << std::endl;

  for (map<string,TList *>::const_iterator iter = m_contours.begin(); iter != m_contours.end(); iter++ ){
    std::cout << "iter->first = " << iter->first << std::endl;
    std::cout << "iter->second = " << iter->second << std::endl;
  }

  std::cout << "contLevel = " << contLevel << std::endl;

  assert(contLevel);

  curv = (TGraph*)(contLevel->First());

  curv->GetXaxis()->SetLimits(parmin(par1),parmax(par1));
  curv->GetYaxis()->SetRangeUser(parmin(par2),parmax(par2));

  curv->SetTitle();
  curv->GetXaxis()->SetTitle(par2latex(par1));
  curv->GetXaxis()->SetTitleFont(42);
  curv->GetYaxis()->SetTitle(par2latex(par2));
  curv->GetYaxis()->SetTitleFont(42);
  curv->GetYaxis()->SetTitleOffset(1.20);

  legend->SetNColumns(2);

  for (int i=0; i<contLevel->GetSize(); i++) {
    assert(curv);
    curv->SetLineColor(kBlue);
    curv->SetLineWidth(2);
    curv->SetLineStyle(9);
    if (!i) {
      curv->Draw("AC");
      legend->AddEntry(curv,"Expected 68% C.L.","L");
    } else 
      curv->Draw("SAME C");
    curv=(TGraph *)(contLevel->After(curv));
  }

  cout << "Drawing expected 95%" << endl;
  
  contLevel = m_contours["exp95"];

  curv = (TGraph*)(contLevel->First());

  for (int i=0; i<contLevel->GetSize(); i++) {
    curv->SetLineColor(kGreen);
    curv->SetLineWidth(2);
    curv->SetLineStyle(9);
    curv->Draw("SAME C");
    if (!i) legend->AddEntry(curv,"Expected 95% C.L.","L");
    curv=(TGraph *)(contLevel->After(curv));
  }

  cout << "Drawing expected 99%" << endl;

  contLevel = m_contours["exp99"];
  curv = (TGraph*)(contLevel->First());
  for (int i=0; i<contLevel->GetSize(); i++) {
    curv->SetLineColor(kRed);
    curv->SetLineWidth(2);
    curv->SetLineStyle(9);
    curv->Draw("SAME C");
    if (!i) legend->AddEntry(curv,"Expected 99% C.L.","L");
    curv=(TGraph *)(contLevel->After(curv));
  }

  
  contLevel = m_contours["obs95"];

  if (contLevel) {
    cout << "Drawing obs95" << endl;

    curv = (TGraph*)(contLevel->First());

    for (int i=0; i<contLevel->GetSize(); i++) {
      curv->Draw("SAME C");
      curv->SetLineWidth(3);
      if (!i) legend->AddEntry(curv,"Observed 95% C.L.","L");
      curv=(TGraph *)(contLevel->After(curv));
    }
//.........这里部分代码省略.........