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

int main() {

  std::vector<float> masses;
  masses.push_back( 300. );
  masses.push_back( 400. );
  masses.push_back( 500. );
  masses.push_back( 750. );
  masses.push_back( 1000. );
  masses.push_back( 1250. );
  masses.push_back( 1500. );
  masses.push_back( 2000. );
  masses.push_back( 2500. );
  masses.push_back( 3000. );
  masses.push_back( 5000. );


  for( unsigned i=0; i<masses.size(); ++i ) {

    std::cout << "-> Starting mass: " << masses[i] << std::endl;

    TFile* file  = TFile::Open( Form("toySignal_m%.0f.root", masses[i]), "recreate" );
    file->cd();

    float thisMass = masses[i];

    RooRealVar x("mass", "mass", thisMass, 0.2*thisMass, 2.*thisMass );

    // Crystal-Ball
    RooRealVar mean( "mean", "mean", thisMass );
    RooRealVar sigma( "sigma", "sigma", 0.015*thisMass );
    RooRealVar alpha( "alpha", "alpha", 1.2 );
    RooRealVar n( "n", "n", 0.81 );
    RooCBShape cb( "cb", "cb", x, mean, sigma, alpha, n );

    // Generate a toy dataset from the interpreted p.d.f
    RooDataSet* data = cb.generate(x,10000) ;
    data->SetName("data");
    data->Write();

    file->Close();
  
  }

  return 0;

}

//.........这里部分代码省略.........
  d0LifetimeErr.setBins(10);

  RooDataSet *ProtoData_bkg = (RooDataSet*)databkg->reduce(RooArgSet(d0LifetimeErr));
  TH1F *err_bkg_h = (TH1F*)d0LifetimeErr.createHistogram("err_bkg_h");
  ProtoData_bkg->fillHistogram(err_bkg_h,RooArgSet(d0LifetimeErr));
  err_bkg_h->Scale(err_sig_h->Integral()/err_bkg_h->Integral());
  RooDataHist terr_bkg("terr_bkg","terr_bkg",RooArgSet(d0LifetimeErr),err_bkg_h);
  RooHistPdf terr_bkg_pdf("terr_bkg_pdf","terr_bkg_pdf",RooArgSet(d0LifetimeErr),terr_bkg,3);

  RooDataSet *ProtoData_flav_1 = (RooDataSet*)RSdata->reduce(RooArgSet(D0flav));
  RooDataSet *ProtoData_flav = (RooDataSet*)ProtoData_flav_1->reduce(EventRange(1,(Int_t)((1508 + 791)*ratio)));
  //RooDataHist flav_sig("flav_sig","flav_sig",RooArgSet(D0flav),*ProtoData_flav);
  //RooHistPdf flav_sig_pdf("flav_sig_pdf","flav_sig_pdf",RooArgSet(D0flav),flav_sig);

  RooDataSet *ProtoData_flav_bkg_1 = (RooDataSet*)databkg->reduce(RooArgSet(D0flav));
  RooDataSet *ProtoData_flav_bkg = (RooDataSet*)ProtoData_flav_bkg_1->reduce(EventRange(1,(Int_t)((663 + 47)*ratio)));
  //RooDataHist flav_bkg("flav_bkg","flav_bkg",RooArgSet(D0flav),*ProtoData_flav_bkg);
  //RooHistPdf flav_bkg_pdf("flav_bkg_pdf","flav_bkg_pdf",RooArgSet(D0flav),flav_bkg);

  ProtoData_flav->append(*ProtoData_flav_bkg);

  RooProdPdf totsig_norm("totsig_norm","totsig_norm",RooArgSet(terr_sig_pdf),Conditional(TOTTime,observ));
  RooProdPdf totmis_norm("totmis_norm","totmis_norm",RooArgSet(terr_sig_pdf),Conditional(Tot_mis,observ));
  RooProdPdf totbkg_norm("totbkg_norm","totbkg_norm",RooArgSet(terr_bkg_pdf),Conditional(Tot_comb,observ));
  /*
  RooProdPdf totsig_norm("totsig_norm","totsig_norm",RooArgSet(terr_sig_pdf),Conditional(TOTsigD0,observ));
  RooProdPdf totmis_norm("totmis_norm","totmis_norm",RooArgSet(terr_sig_pdf),Conditional(Tot_mis,observ));
  RooProdPdf totbkg_norm("totbkg_norm","totbkg_norm",RooArgSet(terr_bkg_pdf),Conditional(Tot_comb,observ));
  */
  //Signal + background
  RooAddPdf TOTpdf("TOTpdf","TOTpdf",RooArgList(totsig_norm,totmis_norm,totbkg_norm),RooArgList(Nsig,Nmis,Ncomb));

  RooMCStudy mgr(totsig_norm,tot_var,FitOptions(Save(1),Extended(0),Minos(0)),Extended(0),ConditionalObservables(D0flav),RooFit::ProtoData(*ProtoData_flav,1));
  //RooMCStudy mgr(totsig_norm,tot_var,FitOptions(ConditionalObservables(RooArgSet(D0flav)),Save(1),Extended(0)),Extended(0),RooFit::ProtoData(*ConditDataset,1));

  TString dataname(filename);
  dataname.Append(".dat");
  TString toyname(filename);
  toyname.Append("_fixed.root");

  if(doGen && !doFit) mgr.generate(nexp,nev,kFALSE,dataname);
  if(doFit && !doGen) mgr.fit(nexp,dataname);
  if(doGen && doFit) mgr.generateAndFit(nexp,nev,kFALSE,dataname);

  if(doFit){
  TFile f(toyname,"RECREATE");
  f.cd();
  (RooFitResult*)mgr.fitResult(0)->Write();
  mgr.fitParDataSet().tree().Write();
  f.Close();
  }

  if(extract){

  RooFitResult *theRes = (RooFitResult*)mgr.fitResult(0);

    RooRealVar myc1("myc1","myc1",-10.,10.);
    RooRealVar myc2("myc2","myc2",-10.,10.);
    RooRealVar myratio("myratio","myratio",0.,0.,1.);
    RooRealVar myx("myx","myx",0.,-1.,1.);
    RooRealVar myy("myy","myy",0.,-1.,1.);

    Double_t NrsNws = 2760./1198329;
    Double_t ratioerr = sqrt(pow(96.,2.) + pow(NrsNws,2.)*pow(2760.,2.))/1198329.;
    RooDataSet *parFloat = new RooDataSet("parFloat","parFloat",RooArgList(myratio,myx,myy));

    RooTimepdf mysigD0("mysigD0","mysigD0",d0Lifetime,m2Kpi_d0mass,m2Kpi0_d0mass,gm,&dalitzSpace,tau,c1,c2,-1,dirname);
    Double_t myDenom = mysigD0.createIntegral(RooArgSet(m2Kpi_d0mass,m2Kpi0_d0mass,d0Lifetime))->getVal();

    cout << "x' " << c2.getVal()*sqrt(NrsNws/myDenom) << endl;
    cout << "y' " << c1.getVal()*sqrt(NrsNws/myDenom) << endl;
    cout << "r0^2 " << NrsNws/myDenom << endl;

    for(Int_t j=0;j<1000;j++){
      cout << "Performing step number " << j << endl;
      RooArgList floated = theRes->randomizePars();
      myc1.setVal(((RooAbsReal*)floated.find("c1"))->getVal());
      myc2.setVal(((RooAbsReal*)floated.find("c2"))->getVal());

      Double_t myDenom = mysigD0.createIntegral(RooArgSet(m2Kpi_d0mass,m2Kpi0_d0mass,d0Lifetime))->getVal();
      Double_t myNum = RooRandom::randomGenerator()->Gaus(NrsNws,ratioerr);

      myratio.setVal(myNum/myDenom);
      myx.setVal(myc2.getVal()*sqrt((myNum/myDenom)));
      myy.setVal(myc1.getVal()*sqrt((myNum/myDenom)));

      parFloat->add(RooArgSet(myratio,myx,myy));
    }

    TString floatname(filename);
    floatname.Append("_float.root");

    TFile *f1 = new TFile(floatname,"RECREATE");
    f1->cd();
    parFloat->Write();
    f1->Close();
  }

  return 0;
}

//.........这里部分代码省略.........
  RooDataSet *ProtoData_err = (RooDataSet*)RSdata->reduce(RooArgSet(d0LifetimeErr));
  TH1F *err_sig_h = (TH1F*)d0LifetimeErr.createHistogram("err_sig_h");
  ProtoData_err->fillHistogram(err_sig_h,RooArgSet(d0LifetimeErr));
  RooDataHist terr_sig("terr_sig","terr_sig",RooArgSet(d0LifetimeErr),err_sig_h);
  RooHistPdf terr_sig_pdf("terr_sig_pdf","terr_sig_pdf",RooArgSet(d0LifetimeErr),terr_sig,3);
  d0LifetimeErr.setBins(10);

  RooDataSet *ProtoData_bkg = (RooDataSet*)bkgdata->reduce(RooArgSet(d0LifetimeErr));
  TH1F *err_bkg_h = (TH1F*)d0LifetimeErr.createHistogram("err_bkg_h");
  ProtoData_bkg->fillHistogram(err_bkg_h,RooArgSet(d0LifetimeErr));
  err_bkg_h->Scale(err_sig_h->Integral()/err_bkg_h->Integral());
  RooDataHist terr_bkg("terr_bkg","terr_bkg",RooArgSet(d0LifetimeErr),err_bkg_h);
  RooHistPdf terr_bkg_pdf("terr_bkg_pdf","terr_bkg_pdf",RooArgSet(d0LifetimeErr),terr_bkg,3);

  RooProdPdf totsig_norm("totsig_norm","totsig_norm",RooArgSet(terr_sig_pdf),Conditional(TOTsigD0,observ));
  RooProdPdf totmis_norm("totmis_norm","totmis_norm",RooArgSet(terr_sig_pdf),Conditional(Tot_mis,observ));
  RooProdPdf totbkg_norm("totbkg_norm","totbkg_norm",RooArgSet(terr_bkg_pdf),Conditional(Tot_comb,observ));

  RooProdPdf totsig_norm23("totsig_norm23","totsig_norm23",RooArgSet(terr_sig_pdf),Conditional(TOTsigD023,observ23));
  RooProdPdf totmis_norm23("totmis_norm23","totmis_norm23",RooArgSet(terr_sig_pdf),Conditional(Tot_mis23,observ23));
  RooProdPdf totbkg_norm23("totbkg_norm23","totbkg_norm23",RooArgSet(terr_bkg_pdf),Conditional(Tot_comb23,observ23));

  //Signal + background
  RooAddPdf TOTpdf("TOTpdf","TOTpdf",RooArgList(totsig_norm,totmis_norm,totbkg_norm),RooArgList(Nsig,Nmis,Ncomb));
  RooAddPdf TOTpdf23("TOTpdf23","TOTpdf23",RooArgList(totsig_norm23,totmis_norm23,totbkg_norm23),RooArgList(Nsig,Nmis,Ncomb));

  if(doFit){

    RooFitResult *theRes = TOTpdf.fitTo(*dataWS,Save(1),Minos(0),Extended(0));

    string filename = "syst_root/fit_WS_" + name + ".root";
    TFile fout(filename.c_str(),"RECREATE");
    fout.cd();
    theRes->Write();
    fout.Close();
  }

  if(extract){
    TFile f("fitWS.root");
    RooFitResult* theRes = (RooFitResult*)f.Get("nll");
 
    RooRealVar myc1("myc1","myc1",-10.,10.);
    RooRealVar myc2("myc2","myc2",-10.,10.);
    RooRealVar myratio("myratio","myratio",0.,0.,1.);
    RooRealVar myx("myx","myx",0.,-1.,1.);
    RooRealVar myy("myy","myy",0.,-1.,1.);
 
    Double_t NrsNws = 2562./1132332.;
    Double_t ratioerr = sqrt(pow(89.,2.) + pow(NrsNws,2.)*pow(2862.,2.))/1132332.;
    RooDataSet *parFloat = new RooDataSet("parFloat","parFloat",RooArgList(myratio,myx,myy));
 
    for(Int_t j=0;j<400;j++){
      cout << "Performing step number " << j << endl;
      RooArgList floated = theRes->randomizePars();
      myc1.setVal(((RooAbsReal*)floated.find("c1"))->getVal());
      myc2.setVal(((RooAbsReal*)floated.find("c2"))->getVal());
 
      RooTimepdf mysigD0("mysigD0","mysigD0",d0Lifetime,m2Kpi_d0mass,m2Kpi0_d0mass,gm,&dalitzSpace,tau,myc1,myc2,-1,dirname);
      Double_t myDenom = mysigD0.createIntegral(RooArgSet(m2Kpi_d0mass,m2Kpi0_d0mass,d0Lifetime))->getVal();
      Double_t myNum = RooRandom::randomGenerator()->Gaus(NrsNws,ratioerr);
 
      myratio.setVal(myNum/myDenom);
      myx.setVal(myc2.getVal()*sqrt((myNum/myDenom)));
      myy.setVal(myc1.getVal()*sqrt((myNum/myDenom)));
 
      parFloat->add(RooArgSet(myratio,myx,myy));