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

void rf602_chi2fit()
{

   // S e t u p   m o d e l
   // ---------------------

   // Declare observable x
   RooRealVar x("x","x",0,10) ;

   // Create two Gaussian PDFs g1(x,mean1,sigma) anf g2(x,mean2,sigma) and their parameters
   RooRealVar mean("mean","mean of gaussians",5) ;
   RooRealVar sigma1("sigma1","width of gaussians",0.5) ;
   RooRealVar sigma2("sigma2","width of gaussians",1) ;

   RooGaussian sig1("sig1","Signal component 1",x,mean,sigma1) ;  
   RooGaussian sig2("sig2","Signal component 2",x,mean,sigma2) ;  
   
   // Build Chebychev polynomial p.d.f.  
   RooRealVar a0("a0","a0",0.5,0.,1.) ;
   RooRealVar a1("a1","a1",0.2,0.,1.) ;
   RooChebychev bkg("bkg","Background",x,RooArgSet(a0,a1)) ;

   // Sum the signal components into a composite signal p.d.f.
   RooRealVar sig1frac("sig1frac","fraction of component 1 in signal",0.8,0.,1.) ;
   RooAddPdf sig("sig","Signal",RooArgList(sig1,sig2),sig1frac) ;

   // Sum the composite signal and background 
   RooRealVar bkgfrac("bkgfrac","fraction of background",0.5,0.,1.) ;
   RooAddPdf  model("model","g1+g2+a",RooArgList(bkg,sig),bkgfrac) ;


   // C r e a t e   b i n n e d   d a t a s e t
   // -----------------------------------------

   RooDataSet* d = model.generate(x,10000) ;
   RooDataHist* dh = d->binnedClone() ;

   // Construct a chi^2 of the data and the model.
   // When a p.d.f. is used in a chi^2 fit, the probability density scaled
   // by the number of events in the dataset to obtain the fit function
   // If model is an extended p.d.f, the expected number events is used
   // instead of the observed number of events.
   model.chi2FitTo(*dh) ;

   // NB: It is also possible to fit a RooAbsReal function to a RooDataHist
   // using chi2FitTo(). 

   // Note that entries with zero bins are _not_ allowed 
   // for a proper chi^2 calculation and will give error
   // messages  
   RooDataSet* dsmall = (RooDataSet*) d->reduce(EventRange(1,100)) ;
   RooDataHist* dhsmall = dsmall->binnedClone() ;
   RooChi2Var chi2_lowstat("chi2_lowstat","chi2",model,*dhsmall) ;
   cout << chi2_lowstat.getVal() << endl ;


}

RooDataHist * genHistFromModelPdf(const char * name, RooAbsPdf *model,  RooRealVar *var,    double ScaleLumi,  int range, int rebin, int seed ) {
  double genEvents =  model->expectedEvents(*var);
  TRandom3 *rndm = new TRandom3();
  rndm->SetSeed(seed);
  double nEvt = rndm->PoissonD( genEvents) ;
  int intEvt = ( (nEvt- (int)nEvt) >= 0.5) ? (int)nEvt +1 : int(nEvt);
  RooDataSet * data = model->generate(*var ,   intEvt   );
  cout<< " expected events for " << name << " = "<< genEvents << endl; 
  cout<< " data->numEntries() for name " << name << " == " << data->numEntries()<< endl;
  // cout<< " nEvt from PoissonD for" << name << " == " << nEvt<< endl;
  //cout<< " cast of nEvt  for" << name << " == " << intEvt<< endl; 
  RooAbsData *binned_data = data->binnedClone();
  TH1 * toy_hist = binned_data->createHistogram( name, *var, Binning(range/rebin )  );
  for(int i = 1; i <= toy_hist->GetNbinsX(); ++i) {
    toy_hist->SetBinError( i,  sqrt( toy_hist->GetBinContent(i)) );
    if(toy_hist->GetBinContent(i) == 0.00) {
      cout<< " WARNING: histo " << name << " has 0 enter in bin number " << i << endl;   
    }
    if(toy_hist->GetBinContent(i) < 0.1) {
      toy_hist->SetBinContent(i, 0.0);
      toy_hist->SetBinError(i, 0.0);
      cout<< " WARNING: setting value 0.0 to histo " << name << " for bin number " << i << endl;   
    }  
  }
  RooDataHist * toy_rooHist = new RooDataHist(name, name , RooArgList(*var), toy_hist );
  return toy_rooHist; 
}

void runUPWW() {
  
  int higgsMass=125;
  double intLumi=5.1;
  int nToys = 10;
  bool draw=true;
  
  using namespace RooFit;
  
  gROOT->ProcessLine(".L ~/tdrstyle.C");
  setTDRStyle();
  gStyle->SetPadLeftMargin(0.16);
  gROOT->ForceStyle();
  gROOT->ProcessLine(".L statsFactory.cc+");
  
  //
  // set up test kind 
  // 
  
  double sigRate;
  double bkgRate;
  
  if(higgsMass==125){
    sigRate = 7.;
    bkgRate = 66.;
  }else{
    cout << "HMMMM.... I don't know that mass point...BYE!" << endl;
    return;
  }
  
  RooRealVar* mll  = new RooRealVar("mll","dilepton mass [GeV]", 12, 80.);
  mll->setBins(17);
  
  RooArgSet* obs = new RooArgSet(*mll) ;
  
  // read signal hypothesis 1
  TChain *tsigHyp1 = new TChain("angles");
  tsigHyp1->Add(Form("datafiles/bdtpresel/%i/SMHiggsWW_%i_JHU.root",higgsMass, higgsMass));
  RooDataSet *sigHyp1Data = new RooDataSet("sigHyp1Data","sigHyp1Data",tsigHyp1,*obs);
  RooDataHist *sigHyp1Hist = sigHyp1Data->binnedClone(0);
  RooHistPdf* sigHyp1Pdf = new RooHistPdf("sigHyp1Pdf", "sigHyp1Pdf", *obs, *sigHyp1Hist);

  // read background
  TChain *bkgTree = new TChain("angles");
  bkgTree->Add(Form("datafiles/bdtpresel/%i/WW_madgraph_8TeV.root",higgsMass));
  RooDataSet *bkgData = new RooDataSet("bkgData","bkgData",bkgTree,*obs);
  RooDataHist *bkgHist = bkgData->binnedClone(0);
  RooHistPdf* bkgPdf = new RooHistPdf("bkgPdf", "bkgPdf", *obs, *bkgHist);
    
  char statResults[25];
  statsFactory *hwwuls;
  sprintf(statResults,"uls_hww125_%.0ffb.root", intLumi);
  hwwuls = new statsFactory(obs, sigHyp1Pdf, sigHyp1Pdf, statResults);
  hwwuls->runUpperLimitWithBackground(sigRate*intLumi, bkgRate*intLumi, bkgPdf, nToys);
  delete hwwuls;
  

  // draw plots 
  if(draw) {
    RooPlot* plot1 = mll->frame();
    TString plot1Name = "mll";
    TCanvas* c1 = new TCanvas("c1","c1",400,400); 
    
    bkgData->plotOn(plot1,MarkerColor(kBlack));
    bkgPdf->plotOn(plot1, LineColor(kBlack), LineStyle(kDashed));
    sigHyp1Data->plotOn(plot1,MarkerColor(kRed));
    sigHyp1Pdf->plotOn(plot1,LineColor(kRed), LineStyle(kDashed));      
    
    // draw...
    plot1->Draw();
    c1->SaveAs(Form("plots/ul/epsfiles/%s.eps", plot1Name.Data()));
    c1->SaveAs(Form("plots/ul/pngfiles/%s.png", plot1Name.Data()));
    
    delete c1;
  }
  
  
}

void rf307_fullpereventerrors()
{
  // B - p h y s i c s   p d f   w i t h   p e r - e v e n t  G a u s s i a n   r e s o l u t i o n
  // ----------------------------------------------------------------------------------------------

  // Observables
  RooRealVar dt("dt","dt",-10,10) ;
  RooRealVar dterr("dterr","per-event error on dt",0.01,10) ;

  // Build a gaussian resolution model scaled by the per-event error = gauss(dt,bias,sigma*dterr)
  RooRealVar bias("bias","bias",0,-10,10) ;
  RooRealVar sigma("sigma","per-event error scale factor",1,0.1,10) ;
  RooGaussModel gm("gm1","gauss model scaled bt per-event error",dt,bias,sigma,dterr) ;

  // Construct decay(dt) (x) gauss1(dt|dterr)
  RooRealVar tau("tau","tau",1.548) ;
  RooDecay decay_gm("decay_gm","decay",dt,tau,gm,RooDecay::DoubleSided) ;



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

  // Use landau p.d.f to get empirical distribution with long tail
  RooLandau pdfDtErr("pdfDtErr","pdfDtErr",dterr,RooConst(1),RooConst(0.25)) ;
  RooDataSet* expDataDterr = pdfDtErr.generate(dterr,10000) ;

  // Construct a histogram pdf to describe the shape of the dtErr distribution
  RooDataHist* expHistDterr = expDataDterr->binnedClone() ;
  RooHistPdf pdfErr("pdfErr","pdfErr",dterr,*expHistDterr) ;


  // C o n s t r u c t   c o n d i t i o n a l   p r o d u c t   d e c a y _ d m ( d t | d t e r r ) * p d f ( d t e r r )
  // ----------------------------------------------------------------------------------------------------------------------

  // Construct production of conditional decay_dm(dt|dterr) with empirical pdfErr(dterr)
  RooProdPdf model("model","model",pdfErr,Conditional(decay_gm,dt)) ;

  // (Alternatively you could also use the landau shape pdfDtErr)
  //RooProdPdf model("model","model",pdfDtErr,Conditional(decay_gm,dt)) ;

  

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

  // Specify external dataset with dterr values to use model_dm as conditional p.d.f.
  RooDataSet* data = model.generate(RooArgSet(dt,dterr),10000) ;

  

  // F i t   c o n d i t i o n a l   d e c a y _ d m ( d t | d t e r r )
  // ---------------------------------------------------------------------

  // Specify dterr as conditional observable
  model.fitTo(*data) ;


  
  // P l o t   c o n d i t i o n a l   d e c a y _ d m ( d t | d t e r r )
  // ---------------------------------------------------------------------


  // Make two-dimensional plot of conditional p.d.f in (dt,dterr)
  TH1* hh_model = model.createHistogram("hh_model",dt,Binning(50),YVar(dterr,Binning(50))) ;
  hh_model->SetLineColor(kBlue) ;


  // Make projection of data an dt
  RooPlot* frame = dt.frame(Title("Projection of model(dt|dterr) on dt")) ;
  data->plotOn(frame) ;
  model.plotOn(frame) ;


  // Draw all frames on canvas
  TCanvas* c = new TCanvas("rf307_fullpereventerrors","rf307_fullperventerrors",800, 400);
  c->Divide(2) ;
  c->cd(1) ; gPad->SetLeftMargin(0.20) ; hh_model->GetZaxis()->SetTitleOffset(2.5) ; hh_model->Draw("surf") ;
  c->cd(2) ; gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.6) ; frame->Draw() ;



}

void runSigSepWWSingle(int higgsMass, double intLumi, int nToys,  const TestType test, int var, int toy, bool draw, const unsigned int seed) {

    // location of data
    // for ucsd batch submission
    const char *dataLocation = "root://xrootd.unl.edu//store/user/yygao/HWWAngular/datafiles/";
    // for local tests
    // const char *dataLocation = "datafiles/";
    

    //
    // set up test kind 
    // 

    TString testName = getTestName(test);
    TString varName = getVarName(var);
    TString toyName = getToyName(toy);
    
    std::cout << "Doing " << toyName << " studies on " << testName << " separation based on " << varName << "\n";

    double lowMt(0.);
    double highMt = higgsMass;
    double sigRate;
    double bkgRate;
    
    if(higgsMass==125){
      sigRate = 13.0;
      bkgRate = 155.;
    }else{
      cout << "HMMMM.... I don't know that mass point...BYE!" << endl;
      return;
    }
    
    RooRealVar* dphill = new RooRealVar("dphill","#Delta#phi(leptons) [radian]", 0, TMath::Pi());
    dphill->setBins(20);
    RooRealVar* mt  = new RooRealVar("mt","transverse higgs mass", 60, 130);
    mt->setBins(10);
    RooRealVar* mll  = new RooRealVar("mll","dilepton mass [GeV]", 10, 90.);
    mll->setBins(10);
    
    RooArgSet* obs;

    if ( var == DPHI )
      obs = new RooArgSet(*dphill) ;
    
    if ( var == MLL ) 
      obs = new RooArgSet(*mll) ;
    
    if ( var == MLLMT ) 
      obs = new RooArgSet(*mll, *mt) ;

    if ( var == DPHIMT ) 
      obs = new RooArgSet(*dphill, *mt) ;

    //
    // read signal hypothesis 1 always SMHiggs
    // 
    TChain *tsigHyp1 = new TChain("angles");
    tsigHyp1->Add(Form("%s/%i/SMHiggsWW_%i_JHU.root", dataLocation, higgsMass, higgsMass));
    
    RooDataSet *sigHyp1Data = new RooDataSet("sigHyp1Data","sigHyp1Data",tsigHyp1,*obs);
    RooDataHist *sigHyp1Hist = sigHyp1Data->binnedClone(0);
    RooHistPdf* sigHyp1Pdf = new RooHistPdf("sigHyp1Pdf", "sigHyp1Pdf", *obs, *sigHyp1Hist);
      
    // read signal hypothesis 2
    TChain *tsigHyp2 = new TChain("angles");
    TString secondhypName = getSecondHypInputName(test, float(higgsMass));
    tsigHyp2->Add(Form("%s/%i/%s",dataLocation, higgsMass, secondhypName.Data()));
    
    std::cout << secondhypName << "\n";
    
    RooDataSet *sigHyp2Data = new RooDataSet("sigHyp2Data","sigHyp2Data",tsigHyp2,*obs);
    RooDataHist *sigHyp2Hist = sigHyp2Data->binnedClone(0);
    RooHistPdf* sigHyp2Pdf = new RooHistPdf("sigHyp2Pdf", "sigHyp2Pdf", *obs, *sigHyp2Hist);

    // read background
    TChain *bkgTree = new TChain("angles");
    bkgTree->Add(Form("%s/%i/WW_madgraph_8TeV_0j.root",dataLocation,higgsMass));
    RooDataSet *bkgData = new RooDataSet("bkgData","bkgData",bkgTree,*obs);
    RooDataHist *bkgHist = bkgData->binnedClone(0);
    RooHistPdf* bkgPdf = new RooHistPdf("bkgPdf", "bkgPdf", *obs, *bkgHist);


    char statResults[50];
    statsFactory *myHypothesisSeparation;
    sprintf(statResults,"stat_%s_%s_%s_%.0ffb_%u.root",testName.Data(), toyName.Data(), varName.Data(), intLumi, seed);
    printf(statResults);
    myHypothesisSeparation = new statsFactory(obs, sigHyp1Pdf, sigHyp2Pdf, seed, statResults);
    // running pure toys
    myHypothesisSeparation->hypothesisSeparationWithBackground(sigRate*intLumi,sigRate*intLumi,nToys,bkgPdf,bkgRate*intLumi);
    delete myHypothesisSeparation;
    std::cout << "deleted myHypothesisSeparation" << std::endl;

    
    // draw plots 
    if(draw) {
      RooPlot* plot1;
      TString plot1Name;
      TCanvas* c1 = new TCanvas("c1","c1",400,400); 
      
      if ( var == DPHIMT || var == DPHI) {
//.........这里部分代码省略.........