C++ RooAbsReal::plotOn方法代码示例

void rf313_paramranges()
{

  // C r e a t e   3 D   p d f 
  // -------------------------

  // Define observable (x,y,z)
  RooRealVar x("x","x",0,10) ;
  RooRealVar y("y","y",0,10) ;
  RooRealVar z("z","z",0,10) ;

  // Define 3 dimensional pdf
  RooRealVar z0("z0","z0",-0.1,1) ;
  RooPolynomial px("px","px",x,RooConst(0)) ;
  RooPolynomial py("py","py",y,RooConst(0)) ;
  RooPolynomial pz("pz","pz",z,z0) ;
  RooProdPdf pxyz("pxyz","pxyz",RooArgSet(px,py,pz)) ;



  // D e f i n e d   n o n - r e c t a n g u l a r   r e g i o n   R   i n   ( x , y , z ) 
  // -------------------------------------------------------------------------------------

  //
  // R = Z[0 - 0.1*Y^2] * Y[0.1*X - 0.9*X] * X[0 - 10]
  //

  // Construct range parameterized in "R" in y [ 0.1*x, 0.9*x ]
  RooFormulaVar ylo("ylo","0.1*x",x) ;
  RooFormulaVar yhi("yhi","0.9*x",x) ;
  y.setRange("R",ylo,yhi) ;

  // Construct parameterized ranged "R" in z [ 0, 0.1*y^2 ]
  RooFormulaVar zlo("zlo","0.0*y",y) ;
  RooFormulaVar zhi("zhi","0.1*y*y",y) ;
  z.setRange("R",zlo,zhi) ;



  // C a l c u l a t e   i n t e g r a l   o f   n o r m a l i z e d   p d f   i n   R 
  // ----------------------------------------------------------------------------------

  // Create integral over normalized pdf model over x,y,z in "R" region
  RooAbsReal* intPdf = pxyz.createIntegral(RooArgSet(x,y,z),RooArgSet(x,y,z),"R") ;

  // Plot value of integral as function of pdf parameter z0
  RooPlot* frame = z0.frame(Title("Integral of pxyz over x,y,z in region R")) ;
  intPdf->plotOn(frame) ;



  new TCanvas("rf313_paramranges","rf313_paramranges",600,600) ;
  gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.6) ; frame->Draw() ;

  return ;
}

//.........这里部分代码省略.........
   // fixed_again.add( *ws1->var("decay_hi") );
   // fixed_again.add( *ws1->var("raa1") );
   // fixed_again.add( *ws1->var("raa2") );
   //  fixed_again.add( *ws1->var("nsig2_pp") );
   // fixed_again.add( *ws1->var("sigma1") );
   //  fixed_again.add( *ws1->var("nbkg_pp") );
   // fixed_again.add( *ws1->var("npow") );
   // fixed_again.add( *ws1->var("muPlusPt") );
   // fixed_again.add( *ws1->var("muMinusPt") );
   // fixed_again.add( *ws1->var("mscale_hi") );
   // fixed_again.add( *ws1->var("mscale_pp") );
   //  
   // ws1->Print();
   cout << "99999" << endl;

   // create signal+background Model Config
   RooStats::ModelConfig sbHypo("SbHypo");
   sbHypo.SetWorkspace( *ws1 );
   sbHypo.SetPdf( *ws1->pdf("joint") );
   sbHypo.SetObservables( obs );
   sbHypo.SetGlobalObservables( globalObs );
   sbHypo.SetParametersOfInterest( poi );
   sbHypo.SetNuisanceParameters( nuis );
   sbHypo.SetPriorPdf( *ws1->pdf("step") ); // this is optional

   // ws1->Print();
   /////////////////////////////////////////////////////////////////////
   RooAbsReal * pNll = sbHypo.GetPdf()->createNLL( *data,NumCPU(10) );
   cout << "111111" << endl;
   RooMinuit(*pNll).migrad(); // minimize likelihood wrt all parameters before making plots
   cout << "444444" << endl;
   RooPlot *framepoi = ((RooRealVar *)poi.first())->frame(Bins(10),Range(0.,0.2),Title("LL and profileLL in raa3"));
   cout << "222222" << endl;
   pNll->plotOn(framepoi,ShiftToZero());
   cout << "333333" << endl;
   
   RooAbsReal * pProfile = pNll->createProfile( globalObs ); // do not profile global observables
   pProfile->getVal(); // this will do fit and set POI and nuisance parameters to fitted values
   pProfile->plotOn(framepoi,LineColor(kRed));
   framepoi->SetMinimum(0);
   framepoi->SetMaximum(3);
   TCanvas *cpoi = new TCanvas();
   cpoi->cd(); framepoi->Draw();
   cpoi->SaveAs("cpoi.pdf");

   ((RooRealVar *)poi.first())->setMin(0.);
   RooArgSet * pPoiAndNuisance = new RooArgSet("poiAndNuisance");
   // pPoiAndNuisance->add(*sbHypo.GetNuisanceParameters());
   // pPoiAndNuisance->add(*sbHypo.GetParametersOfInterest());
   pPoiAndNuisance->add( nuis );
   pPoiAndNuisance->add( poi );
   sbHypo.SetSnapshot(*pPoiAndNuisance);

   RooPlot* xframeSB = pObs->frame(Title("SBhypo"));
   data->plotOn(xframeSB,Cut("dataCat==dataCat::hi"));
   RooAbsPdf *pdfSB = sbHypo.GetPdf();
   RooCategory *dataCat = ws1->cat("dataCat");
   pdfSB->plotOn(xframeSB,Slice(*dataCat,"hi"),ProjWData(*dataCat,*data));
   TCanvas *c1 = new TCanvas();
   c1->cd(); xframeSB->Draw();
   c1->SaveAs("c1.pdf");

   delete pProfile;
   delete pNll;
   delete pPoiAndNuisance;
   ws1->import( sbHypo );

void rf105_funcbinding()
{

   // B i n d   T M a t h : : E r f   C   f u n c t i o n
   // ---------------------------------------------------

   // Bind one-dimensional TMath::Erf function as RooAbsReal function
   RooRealVar x("x","x",-3,3) ;
   RooAbsReal* errorFunc = bindFunction("erf",TMath::Erf,x) ;

   // Print erf definition
   errorFunc->Print() ;

   // Plot erf on frame 
   RooPlot* frame1 = x.frame(Title("TMath::Erf bound as RooFit function")) ;
   errorFunc->plotOn(frame1) ;


   // B i n d   R O O T : : M a t h : : b e t a _ p d f   C   f u n c t i o n
   // -----------------------------------------------------------------------

   // Bind pdf ROOT::Math::Beta with three variables as RooAbsPdf function
   RooRealVar x2("x2","x2",0,0.999) ;
   RooRealVar a("a","a",5,0,10) ;
   RooRealVar b("b","b",2,0,10) ;
   RooAbsPdf* beta = bindPdf("beta",ROOT::Math::beta_pdf,x2,a,b) ;

   // Perf beta definition
   beta->Print() ;

   // Generate some events and fit
   RooDataSet* data = beta->generate(x2,10000) ;
   beta->fitTo(*data) ;

   // Plot data and pdf on frame
   RooPlot* frame2 = x2.frame(Title("ROOT::Math::Beta bound as RooFit pdf")) ;
   data->plotOn(frame2) ;
   beta->plotOn(frame2) ;



   // B i n d   R O O T   T F 1   a s   R o o F i t   f u n c t i o n
   // ---------------------------------------------------------------

   // Create a ROOT TF1 function
   TF1 *fa1 = new TF1("fa1","sin(x)/x",0,10);

   // Create an observable 
   RooRealVar x3("x3","x3",0.01,20) ;

   // Create binding of TF1 object to above observable
   RooAbsReal* rfa1 = bindFunction(fa1,x3) ;

   // Print rfa1 definition
   rfa1->Print() ;

   // Make plot frame in observable, plot TF1 binding function
   RooPlot* frame3 = x3.frame(Title("TF1 bound as RooFit function")) ;
   rfa1->plotOn(frame3) ;



   TCanvas* c = new TCanvas("rf105_funcbinding","rf105_funcbinding",1200,400) ;
   c->Divide(3) ;
   c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame1->GetYaxis()->SetTitleOffset(1.6) ; frame1->Draw() ;
   c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.6) ; frame2->Draw() ;
   c->cd(3) ; gPad->SetLeftMargin(0.15) ; frame3->GetYaxis()->SetTitleOffset(1.6) ; frame3->Draw() ;

}

void LL(){

  //y0 = 0.000135096401209 sigma_y0 = 0.000103896581837 x0 = 0.000446013873443 sigma_x0 =1.81384394011e-06
  //0.014108652249 0.0168368471049 0.0219755396247 0.000120423865262 1.5575931164 1.55759310722 3.41637854038
  //0.072569437325 0.084063541977 0.0376693978906 0.000284216132439 0.51908074913 0.519080758095 1.12037749267
 // double d = 0.014108652249;
 //  double sd = 0.0168368471049;
 //  double mc = 0.0219755396247;
 //  double smc = 0.000120423865262;
 //  double r0 = d/mc;

  double d = 0.072569437325;
  double sd =  0.084063541977;
  double mc =  0.0376693978906;
  double smc =  0.00028421613243;
  double r0 = d/mc;

  RooRealVar x("x","x",mc*0.9,mc*1.1);
  RooRealVar x0("x0","x0",mc);
  RooRealVar sx("sx","sx",smc);

  RooRealVar r("r","r",r0,0.,5.);
  RooRealVar y0("y0","y0",d); 
  RooRealVar sy("sy","sy",sd); 
  
  RooProduct rx("rx","rx",RooArgList(r,x));

  RooGaussian g1("g1","g1",x,x0,sx);
  RooGaussian g2("g2","g2",rx,y0,sy);

  RooProdPdf LL("LL","LL",g1,g2);

  RooArgSet obs(x0,y0); //observables
  RooArgSet poi(r); //parameters of interest
  RooDataSet data("data", "data", obs);
  data.add(obs); //actually add the data


  RooFitResult* res = LL.fitTo(data,RooFit::Minos(poi),RooFit::Save(),RooFit::Hesse(false));
  if(res->status()==0) {
    r.Print();
    x.Print();
    cout << r.getErrorLo() << " " << r.getErrorHi() << endl;
  } else {
    cout << "Likelihood maximization failed" << endl;
  }
  
  RooAbsReal* nll = LL.createNLL(data); 
  RooPlot* frame = r.frame();
  RooAbsReal* pll = nll->createProfile(poi);
  pll->plotOn(frame);//,RooFit::LineColor(ROOT::kRed));
  frame->Draw();

  r.setVal(0.);
  cout << pll->getVal() << endl; 

  return;
    
    


}

int DiagnosisMacro(int Nbins = 10, int Nsigma = 10, int CPUused = 1, TString Filename = "FIT_DATA_Psi2SJpsi_PPPrompt_Bkg_SecondOrderChebychev_pt65300_rap016_cent0200_262620_263757.root", TString Outputdir = "./")
//Nbins: Number of points for which to calculate profile likelihood. Time required is about (1/CPU) minutes per point per parameter. 0 means do plain likelihood only
//Nsigma: The range in which the scan is performed (value-Nsigma*sigma, value+Nsigma*sigma)
//CPUused: anything larger than 1 causes weird fit results on my laptop, runs fine on lxplus with more (16)

{
    // R e a d   w o r k s p a c e   f r o m   f i l e
    // -----------------------------------------------
    // Open input file with workspace

    //Filename = "FIT_DATA_Psi2SJpsi_PP_Jpsi_DoubleCrystalBall_Psi2S_DoubleCrystalBall_Bkg_Chebychev2_pt6590_rap016_cent0200.root";
    //Filename = "FIT_DATA_Psi2SJpsi_PbPb_Jpsi_DoubleCrystalBall_Psi2S_DoubleCrystalBall_Bkg_Chebychev1_pt6590_rap016_cent0200.root";

    TFile *f = new TFile(Filename);
    // Retrieve workspace from file
    RooWorkspace* w = (RooWorkspace*)f->Get("workspace");

    // Retrieve x,model and data from workspace
    RooRealVar* x = w->var("invMass");
    RooAbsPdf* model = w->pdf("simPdf_syst");
    if (model == 0) {
        model = w->pdf("simPdf");
    }
    if (model == 0) {
        model = w->pdf("pdfMASS_Tot_PP");
    }
    if (model == 0) {
        model = w->pdf("pdfMASS_Tot_PbPb");
    }
    if (model == 0) {
        cout << "[ERROR] pdf failed to load from the workspace" << endl;
        return false;
    }

    RooAbsData* data = w->data("dOS_DATA");
    if (data == 0) {
        data = w->data("dOS_DATA_PP");
    }
    if (data == 0) {
        data = w->data("dOS_DATA_PbPb");
    }
    if (data == 0) {
        cout << "[ERROR] data failed to load from the workspace" << endl;
        return false;
    }

    // Print structure of composite p.d.f.
    model->Print("t");

    /*
    // P l o t   m o d e l
    // ---------------------------------------------------------
    // Plot data and PDF overlaid
    RooPlot* xframe = x->frame(Title("J/psi Model and Data"));
    data->plotOn(xframe);
    model->plotOn(xframe);

    // Draw the frame on the canvas
    TCanvas* c2 = new TCanvas("PlotModel", "PlotModel", 1000, 1000);
    gPad->SetLeftMargin(0.15);
    xframe->GetYaxis()->SetTitleOffset(2.0);
    xframe->Draw();//*/

    ///// Check parameters

    RooArgSet* paramSet1 = model->getDependents(data);
    paramSet1->Print("v");  // Just check
    RooArgSet* paramSet2 = model->getParameters(data);
    paramSet2->Print("v");
    int Nparams = paramSet2->getSize();
    cout << "Number of parameters: " << Nparams<<endl<<endl;


    // C o n s t r u c t   p l a i n   l i k e l i h o o d
    // ---------------------------------------------------
    // Construct unbinned likelihood
    RooAbsReal* nll = model->createNLL(*data, NumCPU(CPUused));
    // Minimize likelihood w.r.t all parameters before making plots
    RooMinuit(*nll).migrad();


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

    ///////////////////   L O O P   O V E R   P A R A M E T E R S

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

    /// Set up loop over parameters
    TString ParamName;
    double ParamValue;
    double ParamError;
    double ParamLimitLow;
    double ParamLimitHigh;
    double FitRangeLow;
    double FitRangeHigh;
    RooRealVar* vParam;
    int counter = 0;

    // Loop start
    TIterator* iter = paramSet2->createIterator();
//.........这里部分代码省略.........

RooWorkspace* makeInvertedANFit(TTree* tree, float forceSigma=-1, bool constrainMu=false, float forceMu=-1) {
  RooWorkspace *ws = new RooWorkspace("ws","");

  std::vector< TString (*)(TString, RooRealVar&, RooWorkspace&) > bkgPdfList;
  bkgPdfList.push_back(makeSingleExp);
  bkgPdfList.push_back(makeDoubleExp);
#if DEBUG==0
  //bkgPdfList.push_back(makeTripleExp);
  bkgPdfList.push_back(makeModExp);
  bkgPdfList.push_back(makeSinglePow);
  bkgPdfList.push_back(makeDoublePow);
  bkgPdfList.push_back(makePoly2);
  bkgPdfList.push_back(makePoly3);
#endif



  RooRealVar mgg("mgg","m_{#gamma#gamma}",103,160,"GeV");
  mgg.setBins(38);

  mgg.setRange("sideband_low", 103,120);
  mgg.setRange("sideband_high",131,160);
  mgg.setRange("signal",120,131);

  RooRealVar MR("MR","",0,3000,"GeV");
  MR.setBins(60);
  
  RooRealVar Rsq("t1Rsq","",0,1,"GeV");
  Rsq.setBins(20);

  RooRealVar hem1_M("hem1_M","",-1,2000,"GeV");
  hem1_M.setBins(40);

  RooRealVar hem2_M("hem2_M","",-1,2000,"GeV");
  hem2_M.setBins(40);

  RooRealVar ptgg("ptgg","p_{T}^{#gamma#gamma}",0,500,"GeV");
  ptgg.setBins(50);

  RooDataSet data("data","",tree,RooArgSet(mgg,MR,Rsq,hem1_M,hem2_M,ptgg));

  RooDataSet* blind_data = (RooDataSet*)data.reduce("mgg<121 || mgg>130");

  std::vector<TString> tags;
  //fit many different background models
  for(auto func = bkgPdfList.begin(); func != bkgPdfList.end(); func++) {
    TString tag = (*func)("bonly",mgg,*ws);
    tags.push_back(tag);
    ws->pdf("bonly_"+tag+"_ext")->fitTo(data,RooFit::Strategy(0),RooFit::Extended(kTRUE),RooFit::Range("sideband_low,sideband_high"));
    RooFitResult* bres = ws->pdf("bonly_"+tag+"_ext")->fitTo(data,RooFit::Strategy(2),RooFit::Save(kTRUE),RooFit::Extended(kTRUE),RooFit::Range("sideband_low,sideband_high"));
    bres->SetName(tag+"_bonly_fitres");
    ws->import(*bres);

    //make blinded fit
    RooPlot *fmgg_b = mgg.frame();
    blind_data->plotOn(fmgg_b,RooFit::Range("sideband_low,sideband_high"));
    TBox blindBox(121,fmgg_b->GetMinimum()-(fmgg_b->GetMaximum()-fmgg_b->GetMinimum())*0.015,130,fmgg_b->GetMaximum());
    blindBox.SetFillColor(kGray);
    fmgg_b->addObject(&blindBox);
    ws->pdf("bonly_"+tag+"_ext")->plotOn(fmgg_b,RooFit::LineColor(kRed),RooFit::Range("Full"),RooFit::NormRange("sideband_low,sideband_high"));
    fmgg_b->SetName(tag+"_blinded_frame");
    ws->import(*fmgg_b);
    delete fmgg_b;
    

    //set all the parameters constant
    RooArgSet* vars = ws->pdf("bonly_"+tag)->getVariables();
    RooFIter iter = vars->fwdIterator();
    RooAbsArg* a;
    while( (a = iter.next()) ){
      if(string(a->GetName()).compare("mgg")==0) continue;
      static_cast<RooRealVar*>(a)->setConstant(kTRUE);
    }

    //make the background portion of the s+b fit
    (*func)("b",mgg,*ws);

    RooRealVar sigma(tag+"_s_sigma","",5,0,100);
    if(forceSigma!=-1) {
      sigma.setVal(forceSigma);
      sigma.setConstant(true);
    }
    RooRealVar mu(tag+"_s_mu","",126,120,132);
    if(forceMu!=-1) {
      mu.setVal(forceMu);
      mu.setConstant(true);
    }
    RooGaussian sig(tag+"_sig_model","",mgg,mu,sigma);
    RooRealVar Nsig(tag+"_sb_Ns","",5,0,100);
    RooRealVar Nbkg(tag+"_sb_Nb","",100,0,100000);
    

    RooRealVar HiggsMass("HiggsMass","",125.1);
    RooRealVar HiggsMassError("HiggsMassError","",0.24);
    RooGaussian HiggsMassConstraint("HiggsMassConstraint","",mu,HiggsMass,HiggsMassError);


    RooAddPdf fitModel(tag+"_sb_model","",RooArgList( *ws->pdf("b_"+tag), sig ),RooArgList(Nbkg,Nsig));

    RooFitResult* sbres;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    }
    PKPhiMCFitResult->Print("v");
    RooRealVar* PkPhiMean=NULL;
    RooRealVar* PkPhiSigma=NULL;
    RooRealVar* PkPhiLAlpha=NULL;
    RooRealVar* PkPhiRAlpha=NULL;
    RooRealVar* PkPhiLN=NULL;
    RooRealVar* PkPhiRN=NULL;
    try {
        PkPhiMean=SafeGetVar(PKPhiMCFitResult,"PkPhiMean");
        PkPhiSigma=SafeGetVar(PKPhiMCFitResult,"PkPhiSigma");
        PkPhiLAlpha=SafeGetVar(PKPhiMCFitResult,"PkPhiLAlpha");
        PkPhiLN=SafeGetVar(PKPhiMCFitResult,"PkPhiLN");
    } catch(std::exception &e) {
        std::cout<<e.what()<<std::endl;
        return 1;
    }

    RooRealVar RarePkPhiMean("RarePkPhiMean","RarePkPhiMean",PkPhiMean->getVal());
    RooRealVar RarePkPhiSigma("RarePkPhiSigma","RarePkPhiSigma",PkPhiSigma->getVal());
    RooRealVar RarePkPhiLN("RarePkPhiLN","RarePkPhiLN",PkPhiLN->getVal());
    RooRealVar RarePkPhiLAlpha("RarePkPhiLAlpha","RarePkPhiLAlpha",PkPhiLAlpha->getVal());

    RooCBShape RarePkPhiModel("RarePkPhiModel","RarePkPhiModel",LbMass,RarePkPhiMean,RarePkPhiSigma,RarePkPhiLAlpha,RarePkPhiLN);

    RooRealVar RarePkPhiYield("RarePkPhiYield","RarePkPhiYield",50.0,1.0,150.0);
    RooExtendPdf RarePkPhiPdf("RarePkPhiPdf","RarePkPhiPdf",RarePkPhiModel,RarePkPhiYield);

    RooAddPdf RarePdf("RarePdf","RarePdf",RooArgList(RarePkPhiPdf,RareBkgPdf,RareSignalPdf));

    /*RarePdf.fitTo(*RareData,Extended(kTRUE));
    RooPlot* RareFrame=LbMass.frame(Bins(35),Range(5200.0,6100.0));
    TCanvas RareCanvas;
    RareData->plotOn(RareFrame);
    RarePdf.plotOn(RareFrame);
    RareFrame->Draw();
    RareCanvas.SaveAs("RareCanvas.pdf");*/

    //________________________________ Fit Rare 2_______________________________
    TFile* Rare2InputFile = new TFile("Rare2SingleFile.root");
    TTree* Rare2InputTree=NULL;
    try {
        Rare2InputTree=SafeGetTree(Rare2InputFile,"DecayTree");
    } catch(std::exception &e) {
        std::cout<<e.what()<<std::endl;
        return 1;
    }
    RooRealVar Rare2Lambda_b0_PE("Lambda_b0_PE","Lambda_b0_PE",-RooNumber::infinity(),RooNumber::infinity());
    RooRealVar Rare2Lambda_b0_PX("Lambda_b0_PX","Lambda_b0_PX",-RooNumber::infinity(),RooNumber::infinity());
    RooRealVar Rare2Lambda_b0_PY("Lambda_b0_PY","Lambda_b0_PY",-RooNumber::infinity(),RooNumber::infinity());
    RooRealVar Rare2Lambda_b0_PZ("Lambda_b0_PZ","Lambda_b0_PZ",-RooNumber::infinity(),RooNumber::infinity());
    RooRealVar Rare2Proton_PE("Proton_PE","Proton_PE",-RooNumber::infinity(),RooNumber::infinity());
    RooRealVar Rare2Proton_PX("Proton_PX","Proton_PX",-RooNumber::infinity(),RooNumber::infinity());
    RooRealVar Rare2Proton_PY("Proton_PY","Proton_PY",-RooNumber::infinity(),RooNumber::infinity());
    RooRealVar Rare2Proton_PZ("Proton_PZ","Proton_PZ",-RooNumber::infinity(),RooNumber::infinity());
    RooRealVar Rare2Kaon_PE("Kaon_PE","Kaon_PE",-RooNumber::infinity(),RooNumber::infinity());
    RooRealVar Rare2Kaon_PX("Kaon_PX","Kaon_PX",-RooNumber::infinity(),RooNumber::infinity());
    RooRealVar Rare2Kaon_PY("Kaon_PY","Kaon_PY",-RooNumber::infinity(),RooNumber::infinity());
    RooRealVar Rare2Kaon_PZ("Kaon_PZ","Kaon_PZ",-RooNumber::infinity(),RooNumber::infinity());
    RooRealVar Rare2eta_prime_PE("eta_prime_PE","eta_prime_PE",-RooNumber::infinity(),RooNumber::infinity());
    RooRealVar Rare2eta_prime_PX("eta_prime_PX","eta_prime_PX",-RooNumber::infinity(),RooNumber::infinity());
    RooRealVar Rare2eta_prime_PY("eta_prime_PY","eta_prime_PY",-RooNumber::infinity(),RooNumber::infinity());
    RooRealVar Rare2eta_prime_PZ("eta_prime_PZ","eta_prime_PZ",-RooNumber::infinity(),RooNumber::infinity());
    RooArgSet Rare2Args(LbMass,Rare2Lambda_b0_PE,Rare2Lambda_b0_PX,Rare2Lambda_b0_PY,Rare2Lambda_b0_PZ,Rare2Proton_PE,Rare2Proton_PX,Rare2Proton_PY,Rare2Proton_PZ);
    Rare2Args.add(Rare2Kaon_PE);
    Rare2Args.add(Rare2Kaon_PX);

void plot_pll(TString fname="monoh_withsm_SRCR_bg11.7_bgslop-0.0_nsig0.0.root")
{
  SetAtlasStyle();



  TFile* file =  TFile::Open(fname);
  RooWorkspace* wspace = (RooWorkspace*) file->Get("wspace");

  cout << "\n\ncheck that eff and reco terms included in BSM component to make fiducial cross-section" <<endl;
  wspace->function("nsig")->Print();
  RooRealVar* reco = wspace->var("reco");
  if(  wspace->function("nsig")->dependsOn(*reco) ) {
    cout << "all good." <<endl;
  } else {
    cout << "need to rerun fit_withsm using DO_FIDUCIAL_LIMIT true" <<endl;
    return;
  }

  /*
  // DANGER
  // TEST WITH EXAGGERATED UNCERTAINTY
  wspace->var("unc_theory")->setMax(1);
  wspace->var("unc_theory")->setVal(1);
  wspace->var("unc_theory")->Print();
  */

  // this was for making plot about decoupling/recoupling approach
  TCanvas* tc = new TCanvas("tc","",400,400);
  RooPlot *frame = wspace->var("xsec_bsm")->frame();
  RooAbsPdf* pdfc = wspace->pdf("jointModeld");
  RooAbsData* data = wspace->data("data");
  RooAbsReal *nllJoint = pdfc->createNLL(*data, RooFit::Constrained()); // slice with fixed xsec_bsm
  RooAbsReal *profileJoint = nllJoint->createProfile(*wspace->var("xsec_bsm"));

  wspace->allVars().Print("v");
  pdfc->fitTo(*data);
  wspace->allVars().Print("v");
  wspace->var("xsec_bsm")->Print();
  double nllmin = 2*nllJoint->getVal();
  wspace->var("xsec_bsm")->setVal(0);
  double nll0 = 2*nllJoint->getVal();
  cout << Form("nllmin = %f, nll0 = %f, Z=%f", nllmin, nll0, sqrt(nll0-nllmin)) << endl;
  nllJoint->plotOn(frame, RooFit::LineColor(kGreen), RooFit::LineStyle(kDotted), RooFit::ShiftToZero(), RooFit::Name("nll_statonly")); // no error
  profileJoint->plotOn(frame,RooFit::Name("pll") );
  wspace->var("xsec_sm")->Print();
  wspace->var("theory")->Print();
  wspace->var("theory")->setConstant();
  profileJoint->plotOn(frame, RooFit::LineColor(kRed), RooFit::LineStyle(kDashed), RooFit::Name("pll_smfixed") );

  frame->GetXaxis()->SetTitle("#sigma_{BSM, fid} [fb]");
  frame->GetYaxis()->SetTitle("-log #lambda  ( #sigma_{BSM, fid} )");
  double temp = frame->GetYaxis()->GetTitleOffset();
  frame->GetYaxis()->SetTitleOffset( 1.1* temp );

  frame->SetMinimum(1e-7);
  frame->SetMaximum(4);


  // Legend
  double x1,y1,x2,y2;
  GetX1Y1X2Y2(tc,x1,y1,x2,y2);
  TLegend *legend_sr=FastLegend(x2-0.75,y2-0.3,x2-0.25,y2-0.5,0.045);
  legend_sr->AddEntry(frame->findObject("pll"),"with #sigma_{SM} uncertainty","L");
  legend_sr->AddEntry(frame->findObject("pll_smfixed"),"with #sigma_{SM} constant","L");
  legend_sr->AddEntry(frame->findObject("nll_statonly"),"no systematics","L");
  frame->Draw();
  legend_sr->Draw("SAME");



  // descriptive text
  vector<TString> pavetext11;
  pavetext11.push_back("#bf{#it{ATLAS Internal}}");
  pavetext11.push_back("#sqrt{#it{s}} = 8 TeV #scale[0.6]{#int}Ldt = 20.3 fb^{-1}");
  pavetext11.push_back("#it{H}+#it{E}_{T}^{miss} , #it{H #rightarrow #gamma#gamma}, #it{m}_{#it{H}} = 125.4 GeV");

  TPaveText* text11=CreatePaveText(x2-0.75,y2-0.25,x2-0.25,y2-0.05,pavetext11,0.045);
  text11->Draw();

  tc->SaveAs("pll.pdf");



  /*
  wspace->var("xsec_bsm")->setConstant(true);
  wspace->var("eff"     )->setConstant(true);
  wspace->var("mh"      )->setConstant(true);
  wspace->var("sigma_h" )->setConstant(true);
  wspace->var("lumi"    )->setConstant(true);
  wspace->var("xsec_sm" )->setVal(v_xsec_sm);
  wspace->var("eff"     )->setVal(1.0);
  wspace->var("lumi"    )->setVal(v_lumi);
  TH1* nllHist = profileJoint->createHistogram("xsec_bsm",100);
  TFile* out = new TFile("nllHist.root","REPLACE");
  nllHist->Write()
  out->Write();
  out->Close();
  */

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

//.........这里部分代码省略.........
  
  
  //RooFormulaVar RareYield2("RareYield2","RareYield2","@0*@1",RooArgSet(ControlYield,YieldRatio2));
  //RooFormulaVar RareYield2("RareYield2","RareYield2","(ControlYield*(1/Mode2EfficiencyRatio))*((ObservableBFRatio/(SubBRRatio*fdFl))-(YieldRatio*ModeEfficiencyRatio))",RooArgSet(ControlYield,Mode2EfficiencyRatio,ObservableBFRatio,SubBRRatio,fdFl,YieldRatio,ModeEfficiencyRatio));
  RooRealVar RareYield2("RareYield2","RareYield2",1.0,0.0,100.0);
  RooExtendPdf ExtRare2("ExtRare2","ExtRare2",RareMode2,RareYield2);
  
  RooRealVar BMass("BMass","BMass",5000.0,5500.0);					
  
  //  RooFormulaVar ControlMean("ControlMean","ControlMean","@0-339.72",RooArgSet(RareMean));
  RooRealVar MCControlSigma("MCControlSigma","MCControlSigma",17.0);
  RooFormulaVar ControlSigma("ControlSigma","ControlSigma","@0*@1",RooArgSet(MCControlSigma,SigmaCorrection));
  //RooRealVar ControlSigma("ControlSigma","ControlSigma",20.0,10.0,40.0);
  RooGaussian ControlMode("ControlMode","ControlMode",BMass,ControlMean,ControlSigma);
  RooFormulaVar ControlYield("ControlYield","ControlYield","(1/ObservableBFRatio)*((ModeEfficiencyRatio*RareYield)+(Mode2EfficiencyRatio*RareYield2))*SubBRRatio*fdFl",RooArgSet(ObservableBFRatio,ModeEfficiencyRatio,RareYield,Mode2EfficiencyRatio,RareYield2,SubBRRatio,fdFl));
  RooExtendPdf ExtControl("ExtControl","ExtControl",ControlMode,ControlYield);

    
  RooDataSet* ControlData=ControlMode.generate(RooArgSet(BMass),GenControl);
  
  RooCategory Mode("Mode","Mode");
  Mode.defineType("Rare");
  Mode.defineType("Rare2");
  Mode.defineType("Control");

  RooDataSet CombData("CombData","CombData",RooArgSet(BMass,LambdaMass),Index(Mode),Import("Rare2",*RareData2),Import("Rare",*RareData),Import("Control",*ControlData));

  RooSimultaneous SimPdf("SimPdf","SimPdf",Mode);
  SimPdf.addPdf(ExtRare,"Rare");
  SimPdf.addPdf(ExtRare2,"Rare2");
  SimPdf.addPdf(ExtControl,"Control");

  RooFitResult* SimResult=SimPdf.fitTo(CombData,Save(kTRUE),Minos(kTRUE));
  /*  double FreeYield=-1*SimResult->minNll();
  std::cout<<"With free yield = "<<SimResult->minNll()<<std::endl;
  RareYield.setVal(0);
  RareYield.setConstant(kTRUE);
  RooFitResult* Rare1Fixed=SimPdf.fitTo(CombData,Save(kTRUE),Minos(kTRUE));
  double NullYield=-1*Rare1Fixed->minNll();
  std::cout<<"With not free yield = "<<Rare1Fixed->minNll()<<std::endl;

  double DeltaLogLikelihood=NullYield-FreeYield;
  std::cout<<"DeltaNll= "<<DeltaLogLikelihood<<std::endl;
  double Significance=TMath::Sqrt(-2*DeltaLogLikelihood); 
  std::cout<<"Significance= "<<Significance<<std::endl;
  SimPdf.fitTo(CombData,Save(kTRUE),Minos(kTRUE));*/

  
  RooPlot* BFrame= BMass.frame(Bins(50),Title("B Mass"),Range(5200.0,5400.0));
  ControlData->plotOn(BFrame);
  ControlMode.plotOn(BFrame);

  RooPlot* LambdaFrame = LambdaMass.frame(Bins(50),Title("Lambda mass"),Range(5200.0,6100.0));
  RareData->plotOn(LambdaFrame);
  ExtRare.plotOn(LambdaFrame);

  RooPlot* LambdaFrame2 = LambdaMass.frame(Bins(50),Title("Lambda mass"),Range(5200.0,6100.0));
  RareData2->plotOn(LambdaFrame2);
  RareMode2.plotOn(LambdaFrame2);
  
  TCanvas BCanvas;
  BFrame->Draw();
  BCanvas.SaveAs("BCanvas.pdf");

  TCanvas LambdaCanvas;
  LambdaFrame->Draw();
  LambdaCanvas.SaveAs("LambdaCanvas.pdf");

  TCanvas LambdaCanvas2;
  LambdaFrame2->Draw();
  LambdaCanvas2.SaveAs("LambdaCanvas2.pdf");
  

  SimResult->Print("v");

  // S imPdf.graphVizTree("model.dot");
  std::cout<<"Real Ratio = "<<GenRare/(double)GenControl<<std::endl;
  ObservableBFRatio.Print("v");
  

  std::cout<<"Lb BF = "<<ObservableBFRatio.getVal()*7.06E-5<<" + "<<ObservableBFRatio.getErrorHi()*7.06E-5<<" - "<<ObservableBFRatio.getErrorLo()*7.06E-5<<std::endl;
  

  //________________________________________________ATTEMPT TO SWEIGHT____________________________________________
  
  RooStats::SPlot* sDataMass = new RooStats::SPlot("sData","An SPlot",*RareData,&ExtRare,RooArgList(RareYield,RareBkgYield));
  std::cout << std::endl <<  "Yield of signal is " << RareYield.getVal() << ".  From sWeights it is " << sDataMass->GetYieldFromSWeight("RareYield") << std::endl;
  std::cout << "Yield of background is " << RareBkgYield.getVal() << ".  From sWeights it is " << sDataMass->GetYieldFromSWeight("RareBkgYield") << std::endl << std::endl;

  RooAbsReal* nll = SimPdf.createNLL(CombData);
  RooMinuit(*nll).migrad();
  
  RooPlot* LLFrame=ObservableBFRatio.frame(Title("Some Title"),Range(0.005,0.03));
  nll->plotOn(LLFrame,ShiftToZero());
  LLFrame->GetYaxis()->SetRangeUser(0.0,1000.0);

  TCanvas LLCanvas;
  LLFrame->Draw();
  LLCanvas.SaveAs("LLCanvas.pdf");
}

void rf605_profilell()
{
  // C r e a t e   m o d e l   a n d   d a t a s e t 
  // -----------------------------------------------

  // Observable
  RooRealVar x("x","x",-20,20) ;

  // Model (intentional strong correlations)
  RooRealVar mean("mean","mean of g1 and g2",0,-10,10) ;
  RooRealVar sigma_g1("sigma_g1","width of g1",3) ;
  RooGaussian g1("g1","g1",x,mean,sigma_g1) ;

  RooRealVar sigma_g2("sigma_g2","width of g2",4,3.0,6.0) ;
  RooGaussian g2("g2","g2",x,mean,sigma_g2) ;

  RooRealVar frac("frac","frac",0.5,0.0,1.0) ;
  RooAddPdf model("model","model",RooArgList(g1,g2),frac) ;

  // Generate 1000 events
  RooDataSet* data = model.generate(x,1000) ;
  


  // C o n s t r u c t   p l a i n   l i k e l i h o o d
  // ---------------------------------------------------

  // Construct unbinned likelihood
  RooAbsReal* nll = model.createNLL(*data,NumCPU(2)) ;

  // Minimize likelihood w.r.t all parameters before making plots
  RooMinuit(*nll).migrad() ;

  // Plot likelihood scan frac 
  RooPlot* frame1 = frac.frame(Bins(10),Range(0.01,0.95),Title("LL and profileLL in frac")) ;
  nll->plotOn(frame1,ShiftToZero()) ;

  // Plot likelihood scan in sigma_g2
  RooPlot* frame2 = sigma_g2.frame(Bins(10),Range(3.3,5.0),Title("LL and profileLL in sigma_g2")) ;
  nll->plotOn(frame2,ShiftToZero()) ;



  // C o n s t r u c t   p r o f i l e   l i k e l i h o o d   i n   f r a c
  // -----------------------------------------------------------------------

  // The profile likelihood estimator on nll for frac will minimize nll w.r.t
  // all floating parameters except frac for each evaluation

  RooAbsReal* pll_frac = nll->createProfile(frac) ;

  // Plot the profile likelihood in frac
  pll_frac->plotOn(frame1,LineColor(kRed)) ;

  // Adjust frame maximum for visual clarity
  frame1->SetMinimum(0) ;
  frame1->SetMaximum(3) ;



  // C o n s t r u c t   p r o f i l e   l i k e l i h o o d   i n   s i g m a _ g 2 
  // -------------------------------------------------------------------------------

  // The profile likelihood estimator on nll for sigma_g2 will minimize nll
  // w.r.t all floating parameters except sigma_g2 for each evaluation
  RooAbsReal* pll_sigmag2 = nll->createProfile(sigma_g2) ;

  // Plot the profile likelihood in sigma_g2
  pll_sigmag2->plotOn(frame2,LineColor(kRed)) ;

  // Adjust frame maximum for visual clarity
  frame2->SetMinimum(0) ;
  frame2->SetMaximum(3) ;



  // Make canvas and draw RooPlots
  TCanvas *c = new TCanvas("rf605_profilell","rf605_profilell",800, 400);
  c->Divide(2);
  c->cd(1) ; gPad->SetLeftMargin(0.15) ; frame1->GetYaxis()->SetTitleOffset(1.4) ; frame1->Draw() ;
  c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame2->GetYaxis()->SetTitleOffset(1.4) ; frame2->Draw() ;

  delete pll_frac ;
  delete pll_sigmag2 ;
  delete nll ;
}