C++ RooDataSet::addColumn方法代码示例

void rf403_weightedevts()
{
  // C r e a t e   o b s e r v a b l e   a n d   u n w e i g h t e d   d a t a s e t 
  // -------------------------------------------------------------------------------

  // Declare observable
  RooRealVar x("x","x",-10,10) ;
  x.setBins(40) ;

  // Construction a uniform pdf
  RooPolynomial p0("px","px",x) ;

  // Sample 1000 events from pdf
  RooDataSet* data = p0.generate(x,1000) ;

 

  // C a l c u l a t e   w e i g h t   a n d   m a k e   d a t a s e t   w e i g h t e d 
  // -----------------------------------------------------------------------------------

  // Construct formula to calculate (fake) weight for events
  RooFormulaVar wFunc("w","event weight","(x*x+10)",x) ;

  // Add column with variable w to previously generated dataset
  RooRealVar* w = (RooRealVar*) data->addColumn(wFunc) ;

  // Dataset d is now a dataset with two observable (x,w) with 1000 entries
  data->Print() ;

  // Instruct dataset wdata in interpret w as event weight rather than as observable
  RooDataSet wdata(data->GetName(),data->GetTitle(),data,*data->get(),0,w->GetName()) ;

  // Dataset d is now a dataset with one observable (x) with 1000 entries and a sum of weights of ~430K
  wdata.Print() ;



  // U n b i n n e d   M L   f i t   t o   w e i g h t e d   d a t a 
  // ---------------------------------------------------------------

  // Construction quadratic polynomial pdf for fitting
  RooRealVar a0("a0","a0",1) ;
  RooRealVar a1("a1","a1",0,-1,1) ;
  RooRealVar a2("a2","a2",1,0,10) ;
  RooPolynomial p2("p2","p2",x,RooArgList(a0,a1,a2),0) ;

  // Fit quadratic polynomial to weighted data

  // NOTE: A plain Maximum likelihood fit to weighted data does in general 
  //       NOT result in correct error estimates, unless individual
  //       event weights represent Poisson statistics themselves.
  //       
  // Fit with 'wrong' errors
  RooFitResult* r_ml_wgt = p2.fitTo(wdata,Save()) ;
  
  // A first order correction to estimated parameter errors in an 
  // (unbinned) ML fit can be obtained by calculating the
  // covariance matrix as
  //
  //    V' = V C-1 V
  //
  // where V is the covariance matrix calculated from a fit
  // to -logL = - sum [ w_i log f(x_i) ] and C is the covariance
  // matrix calculated from -logL' = -sum [ w_i^2 log f(x_i) ] 
  // (i.e. the weights are applied squared)
  //
  // A fit in this mode can be performed as follows:

  RooFitResult* r_ml_wgt_corr = p2.fitTo(wdata,Save(),SumW2Error(kTRUE)) ;



  // P l o t   w e i g h e d   d a t a   a n d   f i t   r e s u l t 
  // ---------------------------------------------------------------

  // Construct plot frame
  RooPlot* frame = x.frame(Title("Unbinned ML fit, binned chi^2 fit to weighted data")) ;

  // Plot data using sum-of-weights-squared error rather than Poisson errors
  wdata.plotOn(frame,DataError(RooAbsData::SumW2)) ;

  // Overlay result of 2nd order polynomial fit to weighted data
  p2.plotOn(frame) ;



  // M L  F i t   o f   p d f   t o   e q u i v a l e n t  u n w e i g h t e d   d a t a s e t
  // -----------------------------------------------------------------------------------------
  
  // Construct a pdf with the same shape as p0 after weighting
  RooGenericPdf genPdf("genPdf","x*x+10",x) ;

  // Sample a dataset with the same number of events as data
  RooDataSet* data2 = genPdf.generate(x,1000) ;

  // Sample a dataset with the same number of weights as data
  RooDataSet* data3 = genPdf.generate(x,43000) ;

  // Fit the 2nd order polynomial to both unweighted datasets and save the results for comparison
  RooFitResult* r_ml_unw10 = p2.fitTo(*data2,Save()) ;
//.........这里部分代码省略.........

exampleScript()
{
  gSystem->CompileMacro("betaHelperFunctions.h"      ,"kO") ;
  gSystem->CompileMacro("RooNormalFromFlatPdf.cxx"      ,"kO") ;
  gSystem->CompileMacro("RooBetaInverseCDF.cxx"      ,"kO") ;
  gSystem->CompileMacro("RooBetaPrimeInverseCDF.cxx" ,"kO") ;
  gSystem->CompileMacro("RooCorrelatedBetaGeneratorHelper.cxx"  ,"kO") ;
  gSystem->CompileMacro("RooCorrelatedBetaPrimeGeneratorHelper.cxx"  ,"kO") ;
  gSystem->CompileMacro("rooFitBetaHelperFunctions.h","kO") ;

  TFile betaTest("betaTest.root","RECREATE");
  betaTest.cd();
  
  RooWorkspace workspace("workspace");
  TString correlatedName("testVariable");
  TString observables("observables");
  TString nuisances("nuisances");

  RooAbsArg* betaOne = getCorrelatedBetaConstraint(workspace,"betaOne","",
						   0.5 , 0.1 ,
						   observables, nuisances,
						   correlatedName );

  printf("\n\n *** constraint name is %s from betaOne and %s\n\n", betaOne->GetName(), correlatedName.Data() ) ;

  RooAbsArg* betaTwo = getCorrelatedBetaConstraint(workspace,"betaTwo","",
						   0 , 0 ,
						   observables, nuisances,
						   correlatedName );

  RooAbsArg* betaThree = getCorrelatedBetaConstraint(workspace,"betaThree","",
						     0.2 , 0.01 ,
						     observables, nuisances,
						     correlatedName );

  RooAbsArg* betaFour = getCorrelatedBetaConstraint(workspace,"betaFour","",
						    0.7 , 0.1 ,
						    observables, nuisances,
						    correlatedName );

  RooAbsArg* betaFourC = getCorrelatedBetaConstraint(workspace,"betaFourC","",
						    0.7 , 0.1 ,
						    observables, nuisances,
						    correlatedName, kTRUE );

  RooAbsArg* betaPrimeOne = getCorrelatedBetaPrimeConstraint(workspace,"betaPrimeOne","",
							     1.0 , 0.5 ,
							     observables, nuisances,
							     correlatedName );

  RooAbsArg* betaPrimeOneC = getCorrelatedBetaPrimeConstraint(workspace,"betaPrimeOneC","",
							     1.0 , 0.5 ,
							     observables, nuisances,
							     correlatedName, kTRUE );

  RooAbsArg* betaPrimeTwo = getCorrelatedBetaPrimeConstraint(workspace,"betaPrimeTwo","",
							     0.7 , 0.5 ,
							     observables, nuisances,
							     correlatedName );

  RooAbsArg* betaPrimeThree = getCorrelatedBetaPrimeConstraint(workspace,"betaPrimeThree","",
							       0.1 , 0.05 ,
							       observables, nuisances,
							       correlatedName );

  RooAbsArg* betaPrimeFour = getCorrelatedBetaPrimeConstraint(workspace,"betaPrimeFour","",
							      7 , 1 ,
							      observables, nuisances,
							      correlatedName );

  RooRealVar* correlatedParameter = workspace.var(correlatedName);

  RooAbsPdf* normalFromFlat = workspace.pdf(correlatedName+"_Constraint");

  RooDataSet* data = normalFromFlat->generate(RooArgSet(*correlatedParameter),1e5);

  data->addColumn(*normalFromFlat);

  data->addColumn(*betaOne);
  data->addColumn(*betaTwo);
  data->addColumn(*betaThree);
  data->addColumn(*betaFour);
  data->addColumn(*betaFourC);
  
  data->addColumn(*betaPrimeOne);
  data->addColumn(*betaPrimeTwo);
  data->addColumn(*betaPrimeThree);
  data->addColumn(*betaPrimeFour);
  data->addColumn(*betaPrimeOneC);

  data->Print("v");

  workspace.Print() ;

  //Setup Plotting Kluges:

  RooRealVar normalPlotter  (correlatedName+"_Constraint" , correlatedName+"_Constraint"  ,0,1);
  RooPlot* normalPlot = normalPlotter.frame();
  data->plotOn(normalPlot);

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

void rf316_llratioplot()
{

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

    // Create observables
    RooRealVar x("x","x",-5,5) ;
    RooRealVar y("y","y",-5,5) ;
    RooRealVar z("z","z",-5,5) ;

    // Create signal pdf gauss(x)*gauss(y)*gauss(z)
    RooGaussian gx("gx","gx",x,RooConst(0),RooConst(1)) ;
    RooGaussian gy("gy","gy",y,RooConst(0),RooConst(1)) ;
    RooGaussian gz("gz","gz",z,RooConst(0),RooConst(1)) ;
    RooProdPdf sig("sig","sig",RooArgSet(gx,gy,gz)) ;

    // Create background pdf poly(x)*poly(y)*poly(z)
    RooPolynomial px("px","px",x,RooArgSet(RooConst(-0.1),RooConst(0.004))) ;
    RooPolynomial py("py","py",y,RooArgSet(RooConst(0.1),RooConst(-0.004))) ;
    RooPolynomial pz("pz","pz",z) ;
    RooProdPdf bkg("bkg","bkg",RooArgSet(px,py,pz)) ;

    // Create composite pdf sig+bkg
    RooRealVar fsig("fsig","signal fraction",0.1,0.,1.) ;
    RooAddPdf model("model","model",RooArgList(sig,bkg),fsig) ;

    RooDataSet* data = model.generate(RooArgSet(x,y,z),20000) ;



    // P r o j e c t   p d f   a n d   d a t a   o n   x
    // -------------------------------------------------

    // Make plain projection of data and pdf on x observable
    RooPlot* frame = x.frame(Title("Projection of 3D data and pdf on X"),Bins(40)) ;
    data->plotOn(frame) ;
    model.plotOn(frame) ;



    // D e f i n e   p r o j e c t e d   s i g n a l   l i k e l i h o o d   r a t i o
    // ----------------------------------------------------------------------------------

    // Calculate projection of signal and total likelihood on (y,z) observables
    // i.e. integrate signal and composite model over x
    RooAbsPdf* sigyz = sig.createProjection(x) ;
    RooAbsPdf* totyz = model.createProjection(x) ;

    // Construct the log of the signal / signal+background probability
    RooFormulaVar llratio_func("llratio","log10(@0)-log10(@1)",RooArgList(*sigyz,*totyz)) ;



    // P l o t   d a t a   w i t h   a   L L r a t i o   c u t
    // -------------------------------------------------------

    // Calculate the llratio value for each event in the dataset
    data->addColumn(llratio_func) ;

    // Extract the subset of data with large signal likelihood
    RooDataSet* dataSel = (RooDataSet*) data->reduce(Cut("llratio>0.7")) ;

    // Make plot frame
    RooPlot* frame2 = x.frame(Title("Same projection on X with LLratio(y,z)>0.7"),Bins(40)) ;

    // Plot select data on frame
    dataSel->plotOn(frame2) ;



    // M a k e   M C   p r o j e c t i o n   o f   p d f   w i t h   s a m e   L L r a t i o   c u t
    // ---------------------------------------------------------------------------------------------

    // Generate large number of events for MC integration of pdf projection
    RooDataSet* mcprojData = model.generate(RooArgSet(x,y,z),10000) ;

    // Calculate LL ratio for each generated event and select MC events with llratio)0.7
    mcprojData->addColumn(llratio_func) ;
    RooDataSet* mcprojDataSel = (RooDataSet*) mcprojData->reduce(Cut("llratio>0.7")) ;

    // Project model on x, integrating projected observables (y,z) with Monte Carlo technique
    // on set of events with the same llratio cut as was applied to data
    model.plotOn(frame2,ProjWData(*mcprojDataSel)) ;



    TCanvas* c = new TCanvas("rf316_llratioplot","rf316_llratioplot",800,400) ;
    c->Divide(2) ;
    c->cd(1) ;
    gPad->SetLeftMargin(0.15) ;
    frame->GetYaxis()->SetTitleOffset(1.4) ;
    frame->Draw() ;
    c->cd(2) ;
    gPad->SetLeftMargin(0.15) ;
    frame2->GetYaxis()->SetTitleOffset(1.4) ;
    frame2->Draw() ;



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

void rf208_convolution()
{
  // S e t u p   c o m p o n e n t   p d f s 
  // ---------------------------------------

  // Construct observable
  RooRealVar t("t","t",-10,30) ;

  // Construct landau(t,ml,sl) ;
  RooRealVar ml("ml","mean bw",5.,-20,20) ;
  RooRealVar sl("sl","sigma bw",1,0.1,10) ;
  RooBreitWigner bw("bw","bw",t,ml,sl) ;
  
  // Construct gauss(t,mg,sg)
  RooRealVar mg("mg","mg",0) ;
  RooRealVar sg("sg","sg",2,0.1,10) ;
  RooGaussian gauss("gauss","gauss",t,mg,sg) ;


  // C o n s t r u c t   c o n v o l u t i o n   p d f 
  // ---------------------------------------

  // Set #bins to be used for FFT sampling to 10000
  t.setBins(10000,"cache") ; 

  // Construct landau (x) gauss
  RooFFTConvPdf lxg("lxg","bw (X) gauss",t,bw,gauss) ;



  // S a m p l e ,   f i t   a n d   p l o t   c o n v o l u t e d   p d f 
  // ----------------------------------------------------------------------

  // Sample 1000 events in x from gxlx
  RooDataSet* data = lxg.generate(t,10000) ;

  // Fit gxlx to data
  lxg.fitTo(*data) ;

  // Plot data, landau pdf, landau (X) gauss pdf
  RooPlot* frame = t.frame(Title("landau (x) gauss convolution")) ;
  data->plotOn(frame) ;
  lxg.plotOn(frame) ;
  bw.plotOn(frame,LineStyle(kDashed)) ;


  // Draw frame on canvas
  new TCanvas("rf208_convolution","rf208_convolution",600,600) ;
  gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.4) ; frame->Draw() ;

  //add a variable to the dataset
  RooFormulaVar *r_formula     = new RooFormulaVar("r_formula","","@0",t);
  RooRealVar* r = (RooRealVar*) data->addColumn(*r_formula);
  r->SetName("r");
  r->SetTitle("r");

  RooDataSet* data_r =(RooDataSet*) data->reduce(*r, "");
  r->setRange("sigrange",-10.,30.);
  RooPlot* r_frame = r->frame(Range("sigRange"),Title(" r (x) gauss convolution")) ;
  data_r->plotOn(r_frame, MarkerColor(kRed));
  r_frame->GetXaxis()->SetRangeUser(-10., 30.);
  r_frame->Draw() ;
}

//.........这里部分代码省略.........
    RooWorkspace* wsp = (RooWorkspace*) Fworkspace->Get("wsp");
    wsp->Print();

    RooRealVar* B_postcalib_M = wsp -> var("B_postcalib_M");
    RooRealVar* nsig_sw = wsp-> var("nsig_sw");
    RooRealVar* nbkg_sw = wsp-> var("nbkg_sw");
    RooRealVar* B_M13_Subst3_gamma2pi0 = wsp-> var("B_M13_Subst3_gamma2pi0");
    RooRealVar* B_M023 = wsp -> var("B_M023");
    RooRealVar* K_1_1270_plus_M = wsp -> var("K_1_1270_plus_M");
    RooRealVar* K_1_1270_plus_SMALLESTDELTACHI2 = wsp -> var("K_1_1270_plus_SMALLESTDELTACHI2");
    RooRealVar* gamma_CL = wsp -> var("gamma_CL");
    RooRealVar* piminus_PIDK = wsp -> var("piminus_PIDK");
    RooRealVar* piplus_PIDK = wsp -> var("piplus_PIDK");
    RooRealVar* Kplus_PIDp = wsp -> var("Kplus_PIDp");
    RooRealVar* Kplus_PIDK = wsp -> var("Kplus_PIDK");
    RooRealVar* B_M02 = wsp -> var("B_M02");
    RooRealVar* L_nsig = wsp -> var("L_nsig");
    RooRealVar* L_nbkg = wsp -> var("L_nbkg");

    RooArgSet arg(*B_postcalib_M,*gamma_CL,*B_M13_Subst3_gamma2pi0,*B_M023,*piminus_PIDK,*piplus_PIDK,*Kplus_PIDK,*Kplus_PIDp);
    arg.add(*K_1_1270_plus_M);
    arg.add(*K_1_1270_plus_SMALLESTDELTACHI2);
    arg.add(*B_M02);
    arg.add(*nsig_sw);
    arg.add(*L_nsig);
    arg.add(*nbkg_sw);
    arg.add(*L_nbkg);

    arg.add(*m_Kpipi);


    RooDataSet* DataSWeights = (RooDataSet*) wsp -> data("DataSWeights");
    RooFormulaVar newMass("m_Kpipi", "m_Kpipi", "K_1_1270_plus_M", RooArgList(*(wsp->var("K_1_1270_plus_M"))));
    DataSWeights->addColumn(newMass);
    RooDataSet* splot = new RooDataSet(DataSWeights->GetName(),DataSWeights->GetTitle(),DataSWeights,RooArgSet(arg),"","nsig_sw");


    cout << "\n\n >>>> Defining components and fitting  \n\n" << endl;

    // Defining here pdfs for other resonances to be fitted

    // K1(1270)

    Double_t R = 0.0015; //  was 3.1 GeV-1
    RooRealVar mean_K1_1270("mean_K1_1270","",1272.,1262.,1282.);
    RooRealVar width_K1_1270("width_K1_1270","",90.,70.,110.);

    // K1(1270) -> K rho
    /*RooBreitWigner totalPdf_K1_1270("totalPdf_K1_1270","",*m_Kpipi,mean_K1_1270,width_K1_1270);*/
    /*RooRelBreitWigner totalPdf_K1_1270("totalPdf_K1_1270","totalPdf_K1_1270_rho",*m_Kpipi,mean_K1_1270,width_K1_1270,RooConst(0),RooConst(R),RooConst(770.),RooConst(493.7));*/

    // K1(1270) -> K*0(1430) pi
    /*RooBreitWigner totalPdf_K1_1270_Kst0_1430("totalPdf_K1_1270_Kst0_1430","totalPdf_K1_1270_Kst0_1430",*m_Kpipi,mean_K1_1270,width_K1_1270);*/
    /*RooRelBreitWigner totalPdf_K1_1270_Kst0_1430("totalPdf_K1_1270_Kst0_1430","totalPdf_K1_1270_Kst0_1430",*m_Kpipi,mean_K1_1270,width_K1_1270,RooConst(0),RooConst(R),RooConst(1425.),RooConst(139.6));*/

    // K1(1270) -> K*0(892) pi
    /*RooBreitWigner totalPdf_K1_1270_Kst0_892("totalPdf_K1_1270_Kst0_892","totalPdf_K1_1270_Kst0_892",*m_Kpipi,mean_K1_1270,width_K1_1270);*/
    /*RooRelBreitWigner totalPdf_K1_1270_Kst0_892("totalPdf_K1_1270_Kst0_892","totalPdf_K1_1270_Kst0_892",*m_Kpipi,mean_K1_1270,width_K1_1270,RooConst(0),RooConst(R),RooConst(895.5),RooConst(139.6));*/

    // K1(1400)

    RooRealVar mean_K1_1400("mean_K1_1400","",1403./*,1396.,1410.*/);
    RooRealVar width_K1_1400("width_K1_1400","",174./*,151.,197.*/);
    RooBreitWigner totalPdf_K1_1400("totalPdf_K1_1400","",*m_Kpipi,mean_K1_1400,width_K1_1400);
    /*[>RooRelBreitWigner totalPdf_K1_1400("totalPdf_K1_1400","totalPdf_K1_1400",*m_Kpipi,mean_K1_1400,width_K1_1400,RooConst(0),RooConst(R),RooConst(895.),RooConst(139.6)); //K*(892) pi is 93%<]*/

//.........这里部分代码省略.........
  //make testing dataset
  RooDataSet *hdata = RooTreeConvert::CreateDataSet("hdata",dtree,vars,weightvar);   

  if (doele) 
    weightvar.SetTitle(prescale1000alt*selcut);
  else
    weightvar.SetTitle(prescale10alt*selcut);
  //make reduced testing dataset for integration over conditional variables
  RooDataSet *hdatasmall = RooTreeConvert::CreateDataSet("hdatasmall",dtree,vars,weightvar);     
    
  //retrieve full pdf from workspace
  RooAbsPdf *sigpdf = ws->pdf("sigpdf");
  
  //input variable corresponding to sceta
  RooRealVar *scEraw = ws->var("var_0");
  scEraw->setRange(1.,2.);
  scEraw->setBins(100);
//  RooRealVar *scetavar = ws->var("var_1");
//  RooRealVar *scphivar = ws->var("var_2");
  
 
  //regressed output functions
  RooAbsReal *sigmeanlim = ws->function("sigmeanlim");
  RooAbsReal *sigwidthlim = ws->function("sigwidthlim");
  RooAbsReal *signlim = ws->function("signlim");
  RooAbsReal *sign2lim = ws->function("sign2lim");

//  RooAbsReal *sigalphalim = ws->function("sigalphalim");
  //RooAbsReal *sigalpha2lim = ws->function("sigalpha2lim");


  //formula for corrected energy/true energy ( 1.0/(etrue/eraw) * regression mean)
  RooFormulaVar ecor("ecor","","1./(@0)*@1",RooArgList(*tgtvar,*sigmeanlim));
  RooRealVar *ecorvar = (RooRealVar*)hdata->addColumn(ecor);
  ecorvar->setRange(0.,2.);
  ecorvar->setBins(800);
  
  //formula for raw energy/true energy (1.0/(etrue/eraw))
  RooFormulaVar raw("raw","","1./@0",RooArgList(*tgtvar));
  RooRealVar *rawvar = (RooRealVar*)hdata->addColumn(raw);
  rawvar->setRange(0.,2.);
  rawvar->setBins(800);

  //clone data and add regression outputs for plotting
  RooDataSet *hdataclone = new RooDataSet(*hdata,"hdataclone");
  RooRealVar *meanvar = (RooRealVar*)hdataclone->addColumn(*sigmeanlim);
  RooRealVar *widthvar = (RooRealVar*)hdataclone->addColumn(*sigwidthlim);
  RooRealVar *nvar = (RooRealVar*)hdataclone->addColumn(*signlim);
  RooRealVar *n2var = (RooRealVar*)hdataclone->addColumn(*sign2lim);
 
//  RooRealVar *alphavar = (RooRealVar*)hdataclone->addColumn(*sigalphalim);
//  RooRealVar *alpha2var = (RooRealVar*)hdataclone->addColumn(*sigalpha2lim);
  
  
  //plot target variable and weighted regression prediction (using numerical integration over reduced testing dataset)
  TCanvas *craw = new TCanvas;
  //RooPlot *plot = tgtvar->frame(0.6,1.2,100);
  RooPlot *plot = tgtvar->frame(0.6,2.0,100);
  hdata->plotOn(plot);
  sigpdf->plotOn(plot,ProjWData(*hdatasmall));
  plot->Draw();
  craw->SaveAs("RawE.pdf");
  craw->SaveAs("RawE.png");
  craw->SetLogy();
  plot->SetMinimum(0.1);
  craw->SaveAs("RawElog.pdf");

//.........这里部分代码省略.........
  bdtpdfid.SetMinCutSignificance(5.0);
  bdtpdfid.SetShrinkage(0.1);
  bdtpdfid.SetMinWeightTotal(200.);
  bdtpdfid.SetMaxNodes(200);
  bdtpdfid.TrainForest(1e6);    
  
  
  RooAbsReal *finalcdferr = eerrpdf.createCDF(eerr);
  
  RooFormulaVar transerr("transerr","","sqrt(2.)*TMath::ErfInverse(2.*@0-1.)",*finalcdferr);
  


  RooAbsReal *finalcdfid = idpdf.createCDF(idmva);
  
  RooFormulaVar transid("transid","","sqrt(2.)*TMath::ErfInverse(2.*@0-1.)",*finalcdfid);
  
  
  RooWorkspace *wsout = new RooWorkspace("wsfiteerr");
  wsout->import(*hdataSingle);
  
  wsout->import(eerrpdf,RecycleConflictNodes());
  wsout->import(idpdf,RecycleConflictNodes());
//   wsout->import(transerr,RecycleConflictNodes());
//   wsout->import(transid,RecycleConflictNodes());
  
  wsout->defineSet("datavars",vars,true);
    
  wsout->writeToFile("hggfiteerr.root");  
  
  
    
  
  RooRealVar *cdfidvar = (RooRealVar*)hdataSingle->addColumn(*finalcdfid);    
  RooRealVar *transidvar = (RooRealVar*)hdataSingle->addColumn(transid);
    
  RooGaussianFast unormpdfid("unormpdfid","",*transidvar,RooConst(0.),RooConst(1.));    

  RooRealVar *cdferrvar = (RooRealVar*)hdataSingle->addColumn(*finalcdferr);    
  RooRealVar *transerrvar = (RooRealVar*)hdataSingle->addColumn(transerr);
    
  RooGaussianFast unormpdferr("unormpdferr","",*transerrvar,RooConst(0.),RooConst(1.));    
  
  
  //RooDataSet *testdata = (RooDataSet*)hdataSingle->reduce("abs(sceta)>1.3 && abs(sceta)<1.4");
  RooDataSet *testdata = hdataSingle;
    
  
  new TCanvas;
  RooPlot *eerrplot = eerr.frame(0.,0.1,200);
  testdata->plotOn(eerrplot);
  eerrpdf.plotOn(eerrplot,ProjWData(*testdata));
  eerrplot->Draw();    
  

  
  new TCanvas;
  RooPlot *transplot = transerrvar->frame(-5.,5.,100);
  hdataSingle->plotOn(transplot);
  unormpdferr.plotOn(transplot);
  transplot->Draw();
  //return;
  
  new TCanvas;
  RooPlot *cdfploterr = cdferrvar->frame(0.,1.,100);
  hdataSingle->plotOn(cdfploterr);

void ws_v03()
{
    gROOT->ProcessLine(".x ./mystyle.C");

    TFile *f1 = new TFile("K1_1270/ws_K1_1270.root");
    TFile *f2 = new TFile("K1_1400/ws_K1_1400.root");
    TFile *f3 = new TFile("K2_1430/ws_K2_1430.root");

    RooWorkspace* ws_K1_1270 = (RooWorkspace*) f1->Get("ws_K1_1270");
    RooWorkspace* ws_K1_1400 = (RooWorkspace*) f2->Get("ws_K1_1400");
    RooWorkspace* ws_K2_1430 = (RooWorkspace*) f3->Get("ws_K2_1430");
    ws_K1_1270->Print();
    ws_K1_1400->Print();
    ws_K2_1430->Print();

    // Importing variables from workspaces

    RooRealVar* m_Kpipi = ws_K1_1270 -> var("m_Kpipi");
    RooAbsPdf* totalPdf_K1_1270 = ws_K1_1270 -> pdf("totalPdf_K1_1270");
    RooAbsData* data_K1_1270 = ws_K1_1270 -> data("totalPdf_K1_1270Data");
    RooHistPdf* pdf_K1_1270_to_Krho = ws_K1_1270 -> pdf("histPdf_K1toKrho");

    /*RooRealVar* m_Kpipi = ws_K1_1400 -> var("m_Kpipi");*/
    /*RooAbsPdf* totalPdf_K1_1400 = ws_K1_1400 -> pdf("totalPdf_K1_1400");*/
    /*RooAbsData* data_K1_1400 = ws_K1_1400 -> data("totalPdf_K1_1400Data");*/
    /*RooHistPdf* pdf_K1_1400_to_Krho = ws_K1_1400 -> pdf("histPdf_K1_1400toKrho");*/

    /*RooRealVar* m_Kpipi = ws_K2_1430 -> var("m_Kpipi");*/
    /*RooAbsPdf* totalPdf_K2_1430 = ws_K2_1430 -> pdf("totalPdf_K2_1430");*/
    /*RooAbsData* data_K2_1430 = ws_K2_1430 -> data("totalPdf_K2_1430Data");*/
    /*RooHistPdf* pdf_K2_1430_to_Krho = ws_K2_1430 -> pdf("histPdf_K2_1430toKrho");*/

    // Plotting pdf

    /*TCanvas* c1 = new TCanvas("c1","canvas",20,20,1200,600);*/
    /*c1->Divide(3,1);*/
    /*c1->cd(1);*/
    /*RooPlot* frame_K1_1270 = m_Kpipi -> frame(Bins(200),Title("K1(1270) -> K#pi#pi"));*/
    /*data_K1_1270->plotOn(frame_K1_1270,DrawOption("C"));*/
    /*frame_K1_1270->Draw();*/
    /*c1->cd(2);*/
    /*RooPlot* frame_K1_1400 = m_Kpipi -> frame(Bins(200),Title("K1(1400) -> K#pi#pi"));*/
    /*data_K1_1400->plotOn(frame_K1_1400,DrawOption("C"));*/
    /*frame_K1_1400->Draw();*/
    /*c1->cd(3);*/
    /*RooPlot* frame_K2_1430= m_Kpipi -> frame(Bins(200),Title("K2*(1430) -> K#pi#pi"));*/
    /*data_K2_1430->plotOn(frame_K2_1430,DrawOption("C"));*/
    /*frame_K2_1430->Draw();*/

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

    TFile *Fworkspace = new TFile("workspace.root");

    RooWorkspace* wsp = (RooWorkspace*) Fworkspace->Get("wsp");
    wsp->Print();

    RooRealVar* B_postcalib_M = wsp -> var("B_postcalib_M");
    RooRealVar* nsig_sw = wsp-> var("nsig_sw");
    RooRealVar* nbkg_sw = wsp-> var("nbkg_sw");
    RooRealVar* B_M13_Subst3_gamma2pi0 = wsp-> var("B_M13_Subst3_gamma2pi0");
    RooRealVar* B_M023 = wsp -> var("B_M023");
    RooRealVar* K_1_1270_plus_M = wsp -> var("K_1_1270_plus_M");
    RooRealVar* K_1_1270_plus_SMALLESTDELTACHI2 = wsp -> var("K_1_1270_plus_SMALLESTDELTACHI2");
    RooRealVar* gamma_CL = wsp -> var("gamma_CL");
    RooRealVar* piminus_PIDK = wsp -> var("piminus_PIDK");
    RooRealVar* piplus_PIDK = wsp -> var("piplus_PIDK");
    RooRealVar* Kplus_PIDp = wsp -> var("Kplus_PIDp");
    RooRealVar* Kplus_PIDK = wsp -> var("Kplus_PIDK");
    RooRealVar* B_M02 = wsp -> var("B_M02");
    RooRealVar* L_nsig = wsp -> var("L_nsig");
    RooRealVar* L_nbkg = wsp -> var("L_nbkg");

    RooArgSet arg(*B_postcalib_M,*gamma_CL,*B_M13_Subst3_gamma2pi0,*B_M023,*piminus_PIDK,*piplus_PIDK,*Kplus_PIDK,*Kplus_PIDp);
    arg.add(*K_1_1270_plus_M);
    arg.add(*K_1_1270_plus_SMALLESTDELTACHI2);
    arg.add(*B_M02);
    arg.add(*nsig_sw);
    arg.add(*L_nsig);
    arg.add(*nbkg_sw);
    arg.add(*L_nbkg);

    arg.add(*m_Kpipi);


    RooDataSet* DataSWeights = (RooDataSet*) wsp -> data("DataSWeights");
    RooFormulaVar newMass("m_Kpipi", "m_Kpipi", "K_1_1270_plus_M", RooArgList(*(wsp->var("K_1_1270_plus_M"))));
    DataSWeights->addColumn(newMass);
    RooDataSet* splot = new RooDataSet(DataSWeights->GetName(),DataSWeights->GetTitle(),DataSWeights,RooArgSet(arg),"","nsig_sw");


    // Defining here pdfs for other resonances to be fitted

// TRY TO ADD 2 SIMPLE BWs FOR K1 1270 AND K1 1400 AS WELL
// also, binned clone and try to fit a toy dataset

        // K1(1270)

    /*Double_t R = 0.0015; //  was 3.1 GeV-1*/
    /*RooRealVar mean_K1_1270("mean_K1_1270","",1272.,1262.,1282.);*/
    /*RooRealVar width_K1_1270("width_K1_1270","",90.,70.,110.);*/
//.........这里部分代码省略.........

void plot_CL_chi2_roofit(char * filename, double min, double max, double initial, double ndof_min, double ndof_max, char * plot, char * var = "chi2")
{

	//gStyle->SetOptStat(0);
	//gStyle->SetOptFit(1);
	//gStyle->SetStatFontSize(0.02);
	TFile * _file0 = TFile::Open(filename);
	TTree * t = (TTree*)_file0->Get("tuple");

	RooRealVar * chi2 = new RooRealVar(var, "#chi^{2}", min, max);
	RooRealVar * ndof = new RooRealVar("ndof", "ndof", initial, ndof_min, ndof_max);
	RooChiSquarePdf * pdf = new RooChiSquarePdf("pdf", "pdf", *chi2, *ndof);

	RooDataSet * data = new RooDataSet("data", "data", RooArgSet(*chi2), RooFit::Import(*t));

	pdf->fitTo(*data);

	char formula[30];
	sprintf(formula, "TMath::Prob(%s,ndof)", var);

	RooFormulaVar * CL_ndof_eff_formula = new RooFormulaVar("CL","CL(#chi^{2})",formula, RooArgList(*chi2, *ndof));
	RooRealVar * CL_ndof_eff = (RooRealVar*) data->addColumn(*CL_ndof_eff_formula);
	CL_ndof_eff->setRange(0, 1);
	RooUniform * uniform = new RooUniform("uniform", "uniform", *CL_ndof_eff);
	uniform->fitTo(*data);

	//RooFormulaVar * CL_ndof_min_formula = new RooFormulaVar("CL","CL(#chi^{2})","TMath::Prob(chi2,39)", RooArgList(*chi2));
	//RooRealVar * CL_ndof_min = (RooRealVar*) data->addColumn(*CL_ndof_min_formula);
	//CL_ndof_min->setRange(0, 1);

	RooPlot * frame0 = chi2->frame(RooFit::Bins(25));
	data->plotOn(frame0);
	pdf->plotOn(frame0);
	pdf->paramOn(frame0, RooFit::Format("NELU", RooFit::AutoPrecision(2)), RooFit::Layout(0.6,0.95,0.75)); 
	data->statOn(frame0, RooFit::Format("NELU", RooFit::AutoPrecision(2)), RooFit::Layout(0.6,0.95,0.95));

	RooPlot * frame1 = CL_ndof_eff->frame(RooFit::Bins(10));
	data->plotOn(frame1);
	uniform->plotOn(frame1);	

        TCanvas * c = new TCanvas("c","c",1200, 600);
        c->Divide(2,1);
        c->cd(1);
	frame0->Draw();
        c->cd(2);
	frame1->Draw();

/*
        char buf[30];
        sprintf(buf, "TMath::Prob(chi2,%f)>>h1", f1->GetParameter(0));
        cout << buf << endl;
        c->Modified();
        c->Update();

        c->cd(2);
        t->Draw("TMath::Prob(chi2,ndof-8)>>h0");
        t->Draw(buf);
        h1->Draw();
        h1->Fit("pol0");
        h0->Draw("same");
        h1->GetXaxis()->SetTitle("CL(#chi^{2})");
        h1->GetYaxis()->SetTitle("Number of toys / 0.1");
        h1->SetMinimum(0);
        h1->SetMaximum(2*t->GetEntries()/nbins);
*/
        c->SaveAs(plot);



}

//.........这里部分代码省略.........
    /*RooArgSet param(*m_Kpipi);*/
    /*RooDataSet MCdataset("MCdataset","MCdataset",param);*/

    /*Double_t K_1_1270_plus_TRUEM = 0.;*/

    /*t_tree->SetBranchAddress("K_1_1270_plus_TRUEM",&K_1_1270_plus_TRUEM);*/

    /*Int_t n = 0;*/
    /*for (int i=0;i<t_tree->GetEntries();i++)*/
    /*{*/
        /*t_tree->GetEntry(i);*/
        /**m_Kpipi_K1_1270 = K_1_1270_plus_TRUEM;*/

        /*MCdataset.add(param);*/
        /*n++;*/
    /*}*/

    /*pdf_K1_1270_to_Krho.fitTo(MCdataset,Extended(true));*/

    /*TCanvas* canvas = new TCanvas("canvas","MC data for K1(1270) -> #rho K",20,20,800,600);*/
    /*RooPlot* frame_MCdataset = m_Kpipi.frame(Bins(200));*/
    /*pdf_K1_1270_to_Krho.paramOn(frame_MCdataset);*/
    /*MCdataset.plotOn(frame_MCdataset,DrawOption("C"));*/
    /*pdf_K1_1270_to_Krho->plotOn(frame_MCdataset);*/
    /*frame_MCdataset->Draw();*/

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

    TFile *Fworkspace = new TFile("workspace.root");

    RooWorkspace* wsp = (RooWorkspace*) Fworkspace->Get("wsp");
    wsp->Print();

    RooRealVar* B_postcalib_M = wsp -> var("B_postcalib_M");
    RooRealVar* nsig_sw = wsp-> var("nsig_sw");
    RooRealVar* nbkg_sw = wsp-> var("nbkg_sw");
    RooRealVar* B_M13_Subst3_gamma2pi0 = wsp-> var("B_M13_Subst3_gamma2pi0");
    RooRealVar* B_M023 = wsp -> var("B_M023");
    RooRealVar* K_1_1270_plus_M = wsp -> var("K_1_1270_plus_M");
    RooRealVar* K_1_1270_plus_SMALLESTDELTACHI2 = wsp -> var("K_1_1270_plus_SMALLESTDELTACHI2");
    RooRealVar* gamma_CL = wsp -> var("gamma_CL");
    RooRealVar* piminus_PIDK = wsp -> var("piminus_PIDK");
    RooRealVar* piplus_PIDK = wsp -> var("piplus_PIDK");
    RooRealVar* Kplus_PIDp = wsp -> var("Kplus_PIDp");
    RooRealVar* Kplus_PIDK = wsp -> var("Kplus_PIDK");
    RooRealVar* B_M02 = wsp -> var("B_M02");
    RooRealVar* L_nsig = wsp -> var("L_nsig");
    RooRealVar* L_nbkg = wsp -> var("L_nbkg");

    RooArgSet arg(*B_postcalib_M,*gamma_CL,*B_M13_Subst3_gamma2pi0,*B_M023,*piminus_PIDK,*piplus_PIDK,*Kplus_PIDK,*Kplus_PIDp);
    arg.add(*K_1_1270_plus_M);
    arg.add(*K_1_1270_plus_SMALLESTDELTACHI2);
    arg.add(*B_M02);
    arg.add(*nsig_sw);
    arg.add(*L_nsig);
    arg.add(*nbkg_sw);
    arg.add(*L_nbkg);

    arg.add(*m_Kpipi);


    RooDataSet* DataSWeights = (RooDataSet*) wsp -> data("DataSWeights");
    RooFormulaVar newMass("m_Kpipi", "m_Kpipi", "K_1_1270_plus_M", RooArgList(*(wsp->var("K_1_1270_plus_M"))));
    DataSWeights->addColumn(newMass);
    RooDataSet* splot = new RooDataSet(DataSWeights->GetName(),DataSWeights->GetTitle(),DataSWeights,RooArgSet(arg),"","nsig_sw");

    // here add pdfs and FIT sum of pdf to splot

    RooRealVar K1_1270_y("K1_1270_y","K1_1270_y",1000.,0.,10000.);
    RooRealVar K1_1400_y("K1_1400_y","K1_1400_y",100.,0.,10000.);
    /*RooRealVar K2_1430_y("K2_1430_y","K2_1430_y",300.,0.,10000.);*/
    RooFormulaVar K2_1430_y("K2_1430_y","K2_1430_y","K1_1270_y/3.",K1_1270_y); // K2_1430 yield is fixed to 1/3 of the K1_1270 yield (as found from Belle ... TO CHECK)

    RooArgList shapes;
    RooArgList yields;
    shapes.add(*totalPdf_K1_1270);
    shapes.add(*totalPdf_K1_1400);
    shapes.add(*totalPdf_K2_1430);
    yields.add(K1_1270_y);
    yields.add(K1_1400_y);
    yields.add(K2_1430_y);

    RooAddPdf PDF("PDF","total Pdf for the resonances considered", shapes,yields);

    PDF->fitTo(*splot,Extended(),SumW2Error(kTRUE),Range(1000,2000));

    // Plotting

    TCanvas* canvas_sPlot = new TCanvas("canvas_sPlot","sPlot with weights",40,20,800,600);
    RooPlot* frame_splot = m_Kpipi->frame(1000,2000,80);
    splot->plotOn(frame_splot);

    PDF->paramOn(frame_splot);
    PDF->plotOn(frame_splot,Components(*totalPdf_K1_1270),LineColor(kRed),LineStyle(kDashed));
    PDF->plotOn(frame_splot,Components(*totalPdf_K1_1400),LineColor(1),LineStyle(kDashed));
    PDF->plotOn(frame_splot,Components(*totalPdf_K2_1430),LineColor(51),LineStyle(kDashed));
    PDF->plotOn(frame_splot);
    frame_splot->Draw();

}

//.........这里部分代码省略.........
  else
    weightvar.SetTitle(selcut);
  
  //make testing dataset
  RooDataSet *hdata = RooTreeConvert::CreateDataSet("hdata",dtree,vars,weightvar);   

  if (doele) 
    weightvar.SetTitle(prescale1000alt*selcut);
  else
    weightvar.SetTitle(prescale10alt*selcut);
  //make reduced testing dataset for integration over conditional variables
  RooDataSet *hdatasmall = RooTreeConvert::CreateDataSet("hdatasmall",dtree,vars,weightvar);     
    
  //retrieve full pdf from workspace
  RooAbsPdf *sigpdf = ws->pdf("sigpdf");
  
  //input variable corresponding to sceta
  RooRealVar *scetavar = ws->var("var_1");
  RooRealVar *scphivar = ws->var("var_2");
  
 
  //regressed output functions
  RooAbsReal *sigmeanlim = ws->function("sigmeanlim");
  RooAbsReal *sigwidthlim = ws->function("sigwidthlim");
  RooAbsReal *signlim = ws->function("signlim");
  RooAbsReal *sign2lim = ws->function("sign2lim");

  RooAbsReal *sigalphalim = ws->function("sigalphalim");
  RooAbsReal *sigalpha2lim = ws->function("sigalpha2lim");


  //formula for corrected energy/true energy ( 1.0/(etrue/eraw) * regression mean)
  RooFormulaVar ecor("ecor","","1./(@0)*@1",RooArgList(*tgtvar,*sigmeanlim));
  RooRealVar *ecorvar = (RooRealVar*)hdata->addColumn(ecor);
  ecorvar->setRange(0.,2.);
  ecorvar->setBins(800);
  
  //formula for raw energy/true energy (1.0/(etrue/eraw))
  RooFormulaVar raw("raw","","1./@0",RooArgList(*tgtvar));
  RooRealVar *rawvar = (RooRealVar*)hdata->addColumn(raw);
  rawvar->setRange(0.,2.);
  rawvar->setBins(800);

  //clone data and add regression outputs for plotting
  RooDataSet *hdataclone = new RooDataSet(*hdata,"hdataclone");
  RooRealVar *meanvar = (RooRealVar*)hdataclone->addColumn(*sigmeanlim);
  RooRealVar *widthvar = (RooRealVar*)hdataclone->addColumn(*sigwidthlim);
  RooRealVar *nvar = (RooRealVar*)hdataclone->addColumn(*signlim);
  RooRealVar *n2var = (RooRealVar*)hdataclone->addColumn(*sign2lim);
 
  RooRealVar *alphavar = (RooRealVar*)hdataclone->addColumn(*sigalphalim);
  RooRealVar *alpha2var = (RooRealVar*)hdataclone->addColumn(*sigalpha2lim);
  
  
  //plot target variable and weighted regression prediction (using numerical integration over reduced testing dataset)
  TCanvas *craw = new TCanvas;
  //RooPlot *plot = tgtvar->frame(0.6,1.2,100);
  RooPlot *plot = tgtvar->frame(0.6,2.0,100);
  hdata->plotOn(plot);
  sigpdf->plotOn(plot,ProjWData(*hdatasmall));
  plot->Draw();
  craw->SaveAs("RawE.pdf");
  craw->SaveAs("RawE.png");
  craw->SetLogy();
  plot->SetMinimum(0.1);
  craw->SaveAs("RawElog.pdf");

void rf603_multicpu()
{

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

    // Create observables
    RooRealVar x("x","x",-5,5) ;
    RooRealVar y("y","y",-5,5) ;
    RooRealVar z("z","z",-5,5) ;

    // Create signal pdf gauss(x)*gauss(y)*gauss(z)
    RooGaussian gx("gx","gx",x,RooConst(0),RooConst(1)) ;
    RooGaussian gy("gy","gy",y,RooConst(0),RooConst(1)) ;
    RooGaussian gz("gz","gz",z,RooConst(0),RooConst(1)) ;
    RooProdPdf sig("sig","sig",RooArgSet(gx,gy,gz)) ;

    // Create background pdf poly(x)*poly(y)*poly(z)
    RooPolynomial px("px","px",x,RooArgSet(RooConst(-0.1),RooConst(0.004))) ;
    RooPolynomial py("py","py",y,RooArgSet(RooConst(0.1),RooConst(-0.004))) ;
    RooPolynomial pz("pz","pz",z) ;
    RooProdPdf bkg("bkg","bkg",RooArgSet(px,py,pz)) ;

    // Create composite pdf sig+bkg
    RooRealVar fsig("fsig","signal fraction",0.1,0.,1.) ;
    RooAddPdf model("model","model",RooArgList(sig,bkg),fsig) ;

    // Generate large dataset
    RooDataSet* data = model.generate(RooArgSet(x,y,z),200000) ;



    // P a r a l l e l   f i t t i n g
    // -------------------------------

    // In parallel mode the likelihood calculation is split in N pieces,
    // that are calculated in parallel and added a posteriori before passing
    // it back to MINUIT.

    // Use four processes and time results both in wall time and CPU time
    model.fitTo(*data,NumCPU(4),Timer(kTRUE)) ;



    // P a r a l l e l   M C   p r o j e c t i o n s
    // ----------------------------------------------

    // Construct signal, total likelihood projection on (y,z) observables and likelihood ratio
    RooAbsPdf* sigyz = sig.createProjection(x) ;
    RooAbsPdf* totyz = model.createProjection(x) ;
    RooFormulaVar llratio_func("llratio","log10(@0)-log10(@1)",RooArgList(*sigyz,*totyz)) ;

    // Calculate likelihood ratio for each event, define subset of events with high signal likelihood
    data->addColumn(llratio_func) ;
    RooDataSet* dataSel = (RooDataSet*) data->reduce(Cut("llratio>0.7")) ;

    // Make plot frame and plot data
    RooPlot* frame = x.frame(Title("Projection on X with LLratio(y,z)>0.7"),Bins(40)) ;
    dataSel->plotOn(frame) ;

    // Perform parallel projection using MC integration of pdf using given input dataSet.
    // In this mode the data-weighted average of the pdf is calculated by splitting the
    // input dataset in N equal pieces and calculating in parallel the weighted average
    // one each subset. The N results of those calculations are then weighted into the
    // final result

    // Use four processes
    model.plotOn(frame,ProjWData(*dataSel),NumCPU(4)) ;


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

}

void rf405_realtocatfuncs()
{

  // D e f i n e   p d f   i n   x ,   s a m p l e   d a t a s e t   i n   x 
  // ------------------------------------------------------------------------


  // Define a dummy PDF in x 
  RooRealVar x("x","x",0,10) ;
  RooArgusBG a("a","argus(x)",x,RooConst(10),RooConst(-1)) ;

  // Generate a dummy dataset 
  RooDataSet *data = a.generate(x,10000) ;



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

  // A RooThresholdCategory is a category function that maps regions in a real-valued 
  // input observable observables to state names. At construction time a 'default'
  // state name must be specified to which all values of x are mapped that are not
  // otherwise assigned
  RooThresholdCategory xRegion("xRegion","region of x",x,"Background") ;

  // Specify thresholds and state assignments one-by-one. 
  // Each statement specifies that all values _below_ the given value 
  // (and above any lower specified threshold) are mapped to the
  // category state with the given name
  //
  // Background | SideBand | Signal | SideBand | Background
  //           4.23       5.23     8.23       9.23 
  xRegion.addThreshold(4.23,"Background") ;
  xRegion.addThreshold(5.23,"SideBand") ;
  xRegion.addThreshold(8.23,"Signal") ;
  xRegion.addThreshold(9.23,"SideBand") ; 



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

  // Add values of threshold function to dataset so that it can be used as observable
  data->addColumn(xRegion) ;

  // Make plot of data in x
  RooPlot* xframe = x.frame(Title("Demo of threshold and binning mapping functions")) ;
  data->plotOn(xframe) ;

  // Use calculated category to select sideband data
  data->plotOn(xframe,Cut("xRegion==xRegion::SideBand"),MarkerColor(kRed),LineColor(kRed)) ;



  // C r e a t e   a   b i n n i n g    r e a l - > c a t   f u n c t i o n
  // ----------------------------------------------------------------------

  // A RooBinningCategory is a category function that maps bins of a (named) binning definition 
  // in a real-valued input observable observables to state names. The state names are automatically
  // constructed from the variable name, the binning name and the bin number. If no binning name
  // is specified the default binning is mapped

  x.setBins(10,"coarse") ;
  RooBinningCategory xBins("xBins","coarse bins in x",x,"coarse") ;



  // U s e   b i n n i n g   f u n c t i o n   f o r   t a b u l a t i o n   a n d   p l o t t i n g
  // -----------------------------------------------------------------------------------------------

  // Print table of xBins state multiplicity. Note that xBins does not need to be an observable in data
  // it can be a function of observables in data as well
  Roo1DTable* xbtable = data->table(xBins) ;
  xbtable->Print("v") ;

  // Add values of xBins function to dataset so that it can be used as observable
  RooCategory* xb = (RooCategory*) data->addColumn(xBins) ;

  // Define range "alt" as including bins 1,3,5,7,9 
  xb->setRange("alt","x_coarse_bin1,x_coarse_bin3,x_coarse_bin5,x_coarse_bin7,x_coarse_bin9") ;
  
  // Construct subset of data matching range "alt" but only for the first 5000 events and plot it on the frame
  RooDataSet* dataSel = (RooDataSet*) data->reduce(CutRange("alt"),EventRange(0,5000)) ;
  dataSel->plotOn(xframe,MarkerColor(kGreen),LineColor(kGreen)) ;



  new TCanvas("rf405_realtocatfuncs","rf405_realtocatfuncs",600,600) ;
  xframe->SetMinimum(0.01) ;
  gPad->SetLeftMargin(0.15) ; xframe->GetYaxis()->SetTitleOffset(1.4) ; xframe->Draw() ;


}

void eregtestingExample(bool dobarrel=true, bool doele=true) {
  
  //output dir
  TString dirname = "/data/bendavid/eregexampletest/eregexampletest_test/"; 
  gSystem->mkdir(dirname,true);
  gSystem->cd(dirname);    
  
  //read workspace from training
  TString fname;
  if (doele && dobarrel) 
    fname = "wereg_ele_eb.root";
  else if (doele && !dobarrel) 
    fname = "wereg_ele_ee.root";
  else if (!doele && dobarrel) 
    fname = "wereg_ph_eb.root";
  else if (!doele && !dobarrel) 
    fname = "wereg_ph_ee.root";
  
  TString infile = TString::Format("/data/bendavid/eregexampletest/%s",fname.Data());
  
  TFile *fws = TFile::Open(infile); 
  RooWorkspace *ws = (RooWorkspace*)fws->Get("wereg");
  
  //read variables from workspace
  RooGBRTargetFlex *meantgt = static_cast<RooGBRTargetFlex*>(ws->arg("sigmeant"));  
  RooRealVar *tgtvar = ws->var("tgtvar");
  
  
  RooArgList vars;
  vars.add(meantgt->FuncVars());
  vars.add(*tgtvar);
   
  //read testing dataset from TTree
  RooRealVar weightvar("weightvar","",1.);

  TTree *dtree;
  
  if (doele) {
    //TFile *fdin = TFile::Open("root://eoscms.cern.ch//eos/cms/store/cmst3/user/bendavid/regTreesAug1/hgg-2013Final8TeV_reg_s12-zllm50-v7n_noskim.root");
    TFile *fdin = TFile::Open("/data/bendavid/regTreesAug1/hgg-2013Final8TeV_reg_s12-zllm50-v7n_noskim.root");

    TDirectory *ddir = (TDirectory*)fdin->FindObjectAny("PhotonTreeWriterSingleInvert");
    dtree = (TTree*)ddir->Get("hPhotonTreeSingle");       
  }
  else {
    TFile *fdin = TFile::Open("root://eoscms.cern.ch///eos/cms/store/cmst3/user/bendavid/idTreesAug1/hgg-2013Final8TeV_ID_s12-h124gg-gf-v7n_noskim.root");
    TDirectory *ddir = (TDirectory*)fdin->FindObjectAny("PhotonTreeWriterPreselNoSmear");
    dtree = (TTree*)ddir->Get("hPhotonTreeSingle");       
  }
  
  //selection cuts for testing
  TCut selcut;
  if (dobarrel) 
    selcut = "ph.genpt>25. && ph.isbarrel && ph.ispromptgen"; 
  else
    selcut = "ph.genpt>25. && !ph.isbarrel && ph.ispromptgen"; 
  
  TCut selweight = "xsecweight(procidx)*puweight(numPU,procidx)";
  TCut prescale10 = "(evt%10==0)";
  TCut prescale10alt = "(evt%10==1)";
  TCut prescale25 = "(evt%25==0)";
  TCut prescale100 = "(evt%100==0)";  
  TCut prescale1000 = "(evt%1000==0)";  
  TCut evenevents = "(evt%2==0)";
  TCut oddevents = "(evt%2==1)";
  TCut prescale100alt = "(evt%100==1)";
  TCut prescale1000alt = "(evt%1000==1)";
  TCut prescale50alt = "(evt%50==1)";
  
  if (doele) 
    weightvar.SetTitle(prescale100alt*selcut);
  else
    weightvar.SetTitle(selcut);
  
  //make testing dataset
  RooDataSet *hdata = RooTreeConvert::CreateDataSet("hdata",dtree,vars,weightvar);   

  if (doele) 
    weightvar.SetTitle(prescale1000alt*selcut);
  else
    weightvar.SetTitle(prescale10alt*selcut);
  //make reduced testing dataset for integration over conditional variables
  RooDataSet *hdatasmall = RooTreeConvert::CreateDataSet("hdatasmall",dtree,vars,weightvar);     
    
  //retrieve full pdf from workspace
  RooAbsPdf *sigpdf = ws->pdf("sigpdf");
  
  //input variable corresponding to sceta
  RooRealVar *scetavar = ws->var("var_1");
  
  //regressed output functions
  RooAbsReal *sigmeanlim = ws->function("sigmeanlim");
  RooAbsReal *sigwidthlim = ws->function("sigwidthlim");
  RooAbsReal *signlim = ws->function("signlim");
  RooAbsReal *sign2lim = ws->function("sign2lim");

  //formula for corrected energy/true energy ( 1.0/(etrue/eraw) * regression mean)
  RooFormulaVar ecor("ecor","","1./(@0)*@1",RooArgList(*tgtvar,*sigmeanlim));
  RooRealVar *ecorvar = (RooRealVar*)hdata->addColumn(ecor);
  ecorvar->setRange(0.,2.);
//.........这里部分代码省略.........

void mvaPUPPETEvaluation() {
  
  //output dir
  TString dirname = "."; 
  gSystem->mkdir(dirname,true);
  gSystem->cd(dirname);    
  
  //read workspace from training
  TString fname;
  fname = "mvaPUPPET.root";
  
  TString infile = fname.Data();
  
  TFile *fws = TFile::Open(infile); 
  RooWorkspace *ws = (RooWorkspace*)fws->Get("wereg");
  
  //read variables from workspace
  RooGBRTargetFlex *perpwidth = static_cast<RooGBRTargetFlex*>(ws->arg("perpwidth"));  
  RooGBRTargetFlex *perpmean  = static_cast<RooGBRTargetFlex*>(ws->arg("perpmean"));  
  RooGBRTargetFlex *parpwidth = static_cast<RooGBRTargetFlex*>(ws->arg("parwidth"));  
  RooGBRTargetFlex *parmean   = static_cast<RooGBRTargetFlex*>(ws->arg("parmean"));  
  RooRealVar *tgtvarX = ws->var("tgtX");
  RooRealVar *tgtvarY = ws->var("tgtY");
  
  
  RooArgList vars;
  vars.add(perpwidth->FuncVars());
  vars.add(perpmean->FuncVars());
  vars.add(perpwidth->FuncVars());
  vars.add(parmean->FuncVars());
  vars.add(*tgtvarX);
  vars.add(*tgtvarY);
   
  //read testing dataset from TTree
  RooRealVar weightvar("weightvar","",1.);

  TChain *dtree = new TChain("tree");
  
  dtree->Add("root://eoscms.cern.ch//store//group/dpg_ecal/alca_ecalcalib/ecalMIBI/rgerosa/PUPPETAnalysis/DYJetsToLL_M-50_TuneCUETP8M1_13TeV-amcatnloFXFX-pythia8_Asympt50ns_MCRUN2_74_V9A_forMVATraining/DYJetsToLL_M-50_TuneCUETP8M1_13TeV-amcatnloFXFX-pythia8/crab_20150724_111858/150724_091912/0000/output_mc_1.root/PUPPET/t");
  //TFile *fdin = TFile::Open("/data/bendavid/regTreesAug1/hgg-2013Final8TeV_reg_s12-zllm50-v7n_noskim.root");

//  TDirectory *ddir = (TDirectory*)fdin->FindObjectAny("");
//  dtree = (TTree*)ddir->Get("t");       
  
  //selection cuts for testing
  TCut selcut;
  selcut = "Boson_daughter==13"; 
  
  TCut selweight = "1";
  TCut prescale100alt = "(evt%100==1)";
  
  weightvar.SetTitle(selcut);
  
  //make testing dataset
  RooDataSet *hdata = RooTreeConvert::CreateDataSet("hdata",dtree,vars,weightvar);   

  weightvar.SetTitle(selcut);
  //make reduced testing dataset for integration over conditional variables
  RooDataSet *hdatasmall = RooTreeConvert::CreateDataSet("hdatasmall",dtree,vars,weightvar);     
    
  //retrieve full pdf from workspace
  RooAbsPdf *sigpdfX = ws->pdf("sigpdfX");
  RooAbsPdf *sigpdfY = ws->pdf("sigpdfY");
  
  //input variable corresponding to sceta
  RooRealVar *scetavar = ws->var("var_1");
  
  //regressed output functions
  RooAbsReal *perpwidthlim = ws->function("perpwidthlim");
  RooAbsReal *perpmeanlim  = ws->function("perpmeanlim");
  RooAbsReal *parwidthlim  = ws->function("parwidthlim");
  RooAbsReal *parmeanlim   = ws->function("parmeanlim");

//////////////////////////////////////////////////////////////////////////////////////
//
// formula for corrected recoil?
//
/*
  //formula for corrected energy/true energy ( 1.0/(etrue/eraw) * regression mean)
  RooFormulaVar ecor("ecor","","1./(@0)*@1",RooArgList(*tgtvar,*sigmeanlim));
  RooRealVar *ecorvar = (RooRealVar*)hdata->addColumn(ecor);
  ecorvar->setRange(0.,2.);
  ecorvar->setBins(800);
  
  //formula for raw energy/true energy (1.0/(etrue/eraw))
  RooFormulaVar raw("raw","","1./@0",RooArgList(*tgtvar));
  RooRealVar *rawvar = (RooRealVar*)hdata->addColumn(raw);
  rawvar->setRange(0.,2.);
  rawvar->setBins(800);
*/
//////////////////////////////////////////////////////////////////////////////////////
  //clone data and add regression outputs for plotting
  RooDataSet *hdataclone = new RooDataSet(*hdata,"hdataclone");

  RooRealVar *perpwidthvar = (RooRealVar*)hdataclone->addColumn(*perpwidthlim);
  RooRealVar *perpmeanvar  = (RooRealVar*)hdataclone->addColumn(*perpmeanlim);
  RooRealVar *parwidthvar  = (RooRealVar*)hdataclone->addColumn(*parwidthlim);
  RooRealVar *parmeanvar   = (RooRealVar*)hdataclone->addColumn(*parmeanlim);
  
  
//.........这里部分代码省略.........