C++ RooAbsData类代码示例

RooAbsData * TwoBody::SetObservableRange( double peak, Double_t window_width, UInt_t minEvents ){
  //
  // Reduce the observable range so ~400 events are used
  // for the full combined dataset
  //

  std::string legend = "[TwoBody::SetObservableRange]: ";

  std::map<std::string, double> _range;

  _range = GetDataRange( data, peak, minEvents, window_width );
  char buf[256];

  // Should not change the range!
  //ws->var("mass")->setRange(_range["low"], _range["high"]);
  //ws->var("mass")->Print();

  // replace data
  //int iTotal = (int)data->sumEntries();
  sprintf(buf, "mass>%f && mass<%f", _range["low"], _range["high"]);
  RooAbsData * _data = data->reduce( RooFit::Cut(buf) );
  // correct the nbkg constraint accordingly
  //double _nbkg = ws->var("nbkg_est_dimuon")->getVal();
  //ws->var("nbkg_est_dimuon")->setVal(_nbkg*_data->sumEntries()/(double)(iTotal));
  //ws->var("nbkg_est_dimuon")->Print();

  data->Print();
  _data->Print();
  delete data;
  return _data;
}

pair<double,double> bkgEvPerGeV(RooWorkspace *work, int m_hyp, int cat, int spin=false){
  
  RooRealVar *mass = (RooRealVar*)work->var("CMS_hgg_mass");
  if (spin) mass = (RooRealVar*)work->var("mass");
  mass->setRange(100,180);
  RooAbsPdf *pdf = (RooAbsPdf*)work->pdf(Form("pdf_data_pol_model_8TeV_cat%d",cat));
  RooAbsData *data = (RooDataSet*)work->data(Form("data_mass_cat%d",cat));
  RooPlot *tempFrame = mass->frame();
  data->plotOn(tempFrame,Binning(80));
  pdf->plotOn(tempFrame);
  RooCurve *curve = (RooCurve*)tempFrame->getObject(tempFrame->numItems()-1);
  double nombkg = curve->Eval(double(m_hyp));
 
  RooRealVar *nlim = new RooRealVar(Form("nlim%d",cat),"",0.,0.,1.e5);
  //double lowedge = tempFrame->GetXaxis()->GetBinLowEdge(FindBin(double(m_hyp)));
  //double upedge  = tempFrame->GetXaxis()->GetBinUpEdge(FindBin(double(m_hyp)));
  //double center  = tempFrame->GetXaxis()->GetBinUpCenter(FindBin(double(m_hyp)));

  nlim->setVal(nombkg);
  mass->setRange("errRange",m_hyp-0.5,m_hyp+0.5);
  RooAbsPdf *epdf = 0;
  epdf = new RooExtendPdf("epdf","",*pdf,*nlim,"errRange");
		
  RooAbsReal *nll = epdf->createNLL(*data,Extended(),NumCPU(4));
  RooMinimizer minim(*nll);
  minim.setStrategy(0);
  minim.setPrintLevel(-1);
  minim.migrad();
  minim.minos(*nlim);
  
  double error = (nlim->getErrorLo(),nlim->getErrorHi())/2.;
  data->Print(); 
  return pair<double,double>(nombkg,error); 
}

 virtual Double_t Evaluate(RooAbsData& data, RooArgSet& /*nullPOI*/) {
   // This is the main method in the interface
   Double_t value = 0.0;
   for(int i=0; i < data.numEntries(); i++) {
     value += data.get(i)->getRealValue(fColumnName.c_str());
   }
   return value;
     }

//____________________________________
void MakePlots(RooWorkspace* wks) {

  // Make plots of the data and the best fit model in two cases:
  // first the signal+background case
  // second the background-only case.

  // get some things out of workspace
  RooAbsPdf* model = wks->pdf("model");
  RooAbsPdf* sigModel = wks->pdf("sigModel");
  RooAbsPdf* zjjModel = wks->pdf("zjjModel");
  RooAbsPdf* qcdModel = wks->pdf("qcdModel");

  RooRealVar* mu = wks->var("mu");
  RooRealVar* invMass = wks->var("invMass");
  RooAbsData* data = wks->data("data");


  //////////////////////////////////////////////////////////
  // Make plots for the Alternate hypothesis, eg. let mu float

  mu->setConstant(kFALSE);

  model->fitTo(*data,Save(kTRUE),Minos(kFALSE), Hesse(kFALSE),PrintLevel(-1));

  //plot sig candidates, full model, and individual componenets
  new TCanvas();
  RooPlot* frame = invMass->frame() ;
  data->plotOn(frame ) ;
  model->plotOn(frame) ;
  model->plotOn(frame,Components(*sigModel),LineStyle(kDashed), LineColor(kRed)) ;
  model->plotOn(frame,Components(*zjjModel),LineStyle(kDashed),LineColor(kBlack)) ;
  model->plotOn(frame,Components(*qcdModel),LineStyle(kDashed),LineColor(kGreen)) ;

  frame->SetTitle("An example fit to the signal + background model");
  frame->Draw() ;
  //  cdata->SaveAs("alternateFit.gif");

  //////////////////////////////////////////////////////////
  // Do Fit to the Null hypothesis.  Eg. fix mu=0

  mu->setVal(0); // set signal fraction to 0
  mu->setConstant(kTRUE); // set constant

  model->fitTo(*data, Save(kTRUE), Minos(kFALSE), Hesse(kFALSE),PrintLevel(-1));

  // plot signal candidates with background model and components
  new TCanvas();
  RooPlot* xframe2 = invMass->frame() ;
  data->plotOn(xframe2, DataError(RooAbsData::SumW2)) ;
  model->plotOn(xframe2) ;
  model->plotOn(xframe2, Components(*zjjModel),LineStyle(kDashed),LineColor(kBlack)) ;
  model->plotOn(xframe2, Components(*qcdModel),LineStyle(kDashed),LineColor(kGreen)) ;

  xframe2->SetTitle("An example fit to the background-only model");
  xframe2->Draw() ;
  //  cbkgonly->SaveAs("nullFit.gif");

}

//#include "/uscms_data/d3/cvernier/DiH_13TeV/CMSSW_7_1_5/src/HiggsAnalysis/CombinedLimit/interface/RooMultiPdf.h"
//#include "HiggsAnalysis/CombinedLimit/interface/RooMultiPdf.h"
void makeRooMultiPdfWorkspace(){

   // Load the combine Library 
   gSystem->Load("libHiggsAnalysisCombinedLimit.so");

   // Open the dummy H->gg workspace 
   TFile *f_hgg = TFile::Open("w_background_Bern.root");
   RooWorkspace *w_hgg = (RooWorkspace*)f_hgg->Get("HbbHbb");
   // The observable (CMS_hgg_mass in the workspace)
   RooRealVar *mass =  w_hgg->var("x");

   // Get three of the functions inside, exponential, linear polynomial, power law
   RooAbsPdf *pdf_exp = w_hgg->pdf("bg_exp");
   RooAbsPdf *pdf_pol = w_hgg->pdf("bg");


   // Fit the functions to the data to set the "prefit" state (note this can and should be redone with combine when doing 
   // bias studies as one typically throws toys from the "best-fit"
   RooAbsData *data = w_hgg->data("data_obs");
   pdf_exp->fitTo(*data);  // index 0
   pdf_pol->fitTo(*data);   // index 2

   // Make a plot (data is a toy dataset)
   RooPlot *plot = mass->frame();   data->plotOn(plot);
   pdf_exp->plotOn(plot,RooFit::LineColor(kBlue));
   pdf_pol->plotOn(plot,RooFit::LineColor(kRed));
   plot->SetTitle("PDF fits to toy data");
   plot->Draw();

   // Make a RooCategory object. This will control which of the pdfs is "active"
   RooCategory cat("pdf_index","Index of Pdf which is active");

   // Make a RooMultiPdf object. The order of the pdfs will be the order of their index, ie for below 
   // 0 == exponential
   // 1 == linear function
   // 2 == powerlaw
   RooArgList mypdfs;
   mypdfs.add(*pdf_exp);
   mypdfs.add(*pdf_pol);
   
   RooMultiPdf multipdf("roomultipdf","All Pdfs",cat,mypdfs);
   
   // As usual make an extended term for the background with _norm for freely floating yield
   RooRealVar norm("roomultipdf_norm","Number of background events",0,10000);
   
   // Save to a new workspace
   TFile *fout = new TFile("background_pdfs.root","RECREATE");
   RooWorkspace wout("backgrounds","backgrounds");
   wout.import(cat);
   wout.import(norm);
   wout.import(multipdf);
   wout.Print();
   wout.Write();

}

void rs500b_PrepareWorkspace_Poisson_withSystematics( TString fileName = "WS_Poisson_withSystematics.root", int type = 1 )
{

  // use a RooWorkspace to store the pdf models, prior informations, list of parameters,...
  RooWorkspace myWS("myWS");

  // Observable
  myWS.factory("x[0,0,1]") ;

  // Pdf in observable,
  myWS.factory("Uniform::sigPdf(x)") ;
  myWS.factory("Uniform::bkgPdf(x)") ;
  myWS.factory("SUM::model(S[100,0,1500]*sigPdf,B[1000,0,3000]*bkgPdf)") ;

  // Background only pdf
  myWS.factory("ExtendPdf::modelBkg(bkgPdf,B)") ;

  // Priors
  myWS.factory("Gaussian::priorNuisance(B,1000,200)") ;
  myWS.factory("Uniform::priorPOI(S)") ;

  // Definition of observables and parameters of interest
  myWS.defineSet("observables","x");
  myWS.defineSet("POI","S");
  myWS.defineSet("parameters","B") ;

  // Generate data
  RooAbsData* data = 0;
  // binned data with fixed number of events
  if (type ==0) data = myWS.pdf("model")->generateBinned(*myWS.set("observables"),myWS.var("S")->getVal(),Name("data"));
  // binned data with Poisson fluctuations
  if (type ==1) data = myWS.pdf("model")->generateBinned(*myWS.set("observables"),Extended(),Name("data"));
  // Asimov data: binned data without any fluctuations (average case)
  if (type == 2)  data = myWS.pdf("model")->generateBinned(*myWS.set("observables"),Name("data"),ExpectedData());
  myWS.import(*data) ;

  myWS.writeToFile(fileName);
  std::cout << "\nRooFit model initialized and stored in " << fileName << std::endl;

  // control plot of the generated data
  RooPlot* plot = myWS.var("x")->frame();
  data->plotOn(plot);
  plot->DrawClone();

}

void MakeSpinPlots::PlotSignalFits(TString tag, TString mcName,TString cosThetaBin){
  TCanvas cv;
  TString cat=tag;
  if(cosThetaBin!="") tag = tag+"_"+cosThetaBin;

  float mean = ws->var(Form("%s_FIT_%s_mean",mcName.Data(),tag.Data()))->getVal();
  RooPlot *frame = ws->var("mass")->frame(105,140,70);//mean-10,mean+10,40);
  RooAbsData *d = ws->data(mcName+"_Combined")->reduce(TString("evtcat==evtcat::")+cat);
  if(cosThetaBin!=""){
    TObjArray *arr = cosThetaBin.Tokenize("_");
    float low  = atof(arr->At(1)->GetName());
    float high = atof(arr->At(2)->GetName());
    d = d->reduce( Form("cosT < %0.2f && cosT >= %0.2f",high,low) );
    delete arr;
  }

  d->plotOn(frame);
  RooFitResult *res = (RooFitResult*)ws->obj(Form("%s_FIT_%s_fitResult",mcName.Data(),tag.Data()));
  RooAbsPdf * pdf = ws->pdf(Form("%s_FIT_%s",mcName.Data(),tag.Data())); //signal model
  std::cout << pdf << "\t" << res << std::endl;
  pdf->plotOn(frame,RooFit::FillColor(kGreen),RooFit::VisualizeError(*res,2.0));
  pdf->plotOn(frame,RooFit::FillColor(kYellow),RooFit::VisualizeError(*res,1.0));
  pdf->plotOn(frame,RooFit::LineColor(kRed));
  d->plotOn(frame); //data
  
  tPair lbl(mcName,tag);

  TLatex *prelim = new TLatex(0.18,0.9,"CMS Preliminary Simulation");
  TLatex *sigL  = new TLatex(0.18,0.6,Form("#sigma_{eff} = %0.2f GeV",fitSigEff[lbl].first,fitSigEff[lbl].second));
  prelim->SetNDC();
  sigL->SetNDC();
  prelim->SetTextSize(0.05);
  sigL->SetTextSize(0.05);
  
  frame->addObject(prelim);
  frame->addObject(sigL);
  frame->Draw();
  cv.SaveAs(basePath+Form("/signalModels/sig_%s_%s_%s.png",mcName.Data(),outputTag.Data(),tag.Data()));
  cv.SaveAs(basePath+Form("/signalModels/C/sig_%s_%s_%s.C",mcName.Data(),outputTag.Data(),tag.Data()));
  cv.SaveAs(basePath+Form("/signalModels/sig_%s_%s_%s.pdf",mcName.Data(),outputTag.Data(),tag.Data()));

}

RooAbsData * Tprime::GetPseudoData( void ) {
    //
    // Generate pseudo data, return a pointer.
    // Class member pointer data set to point to the dataset.
    // Caller does not take ownership.
    //

    static int n_toys = 0;

    // legend for printouts
    std::string legend = "[Tprime::GetPseudoData()]: ";

    delete data;

    // We will use ToyMCSampler to generate pseudo-data (and test statistic, eventually)
    // We are responsible for randomizing nuisances and global observables,
    // ToyMCSampler only generates observables (as of ROOT 5.30.00-rc1 and before)

    // MC sampler and test statistic
    if(n_toys == 0) { // on first entry
        // get B model config from workspace
        RooStats::ModelConfig * pBModel = (RooStats::ModelConfig *)pWs->obj("BModel");
        pBModel->SetWorkspace(*pWs);

        // get parameter of interest set
        //const RooArgSet * poi = pSbModel->GetParametersOfInterest();

        //RooStats::TestStatistic * pTestStatistic = new RooStats::ProfileLikelihoodTestStat(*pBModel->GetPdf());
        //RooStats::ToyMCSampler toymcs(*pTestStatistic, 1);
        pTestStatistic = new RooStats::ProfileLikelihoodTestStat(*pBModel->GetPdf());
        pToyMcSampler = new RooStats::ToyMCSampler(*pTestStatistic, 1);
        pToyMcSampler->SetPdf(*pBModel->GetPdf());
        pToyMcSampler->SetObservables(*pBModel->GetObservables());
        pToyMcSampler->SetParametersForTestStat(*pBModel->GetParametersOfInterest()); // just POI
        pToyMcSampler->SetGlobalObservables(*pBModel->GetGlobalObservables());
    }

    // load parameter point
    pWs->loadSnapshot("parametersToGenerateData");

    RooArgSet dummySet;
    data = pToyMcSampler->GenerateToyData(dummySet);
    std::cout << legend << "generated the following background-only pseudo-data:" << std::endl;
    data->Print();

    // count number of generated toys
    ++n_toys;

    return data;
}

double getChisq(RooAbsData &dat, RooAbsPdf &pdf, RooRealVar &var, bool prt=false) {

    // Find total number of events
    double nEvt;
    double nTot=0.0;

    for(int j=0; j<dat.numEntries(); j++) {
        dat.get(j);
        nEvt=dat.weight();
        nTot+=nEvt;
    }

    // Find chi-squared equivalent 2NLL
    //RooRealVar *var=(RooRealVar*)(pdf.getParameters(*dat)->find("CMS_hgg_mass"));
    double totNLL=0.0;
    double prbSum=0.0;

    for(int j=0; j<dat.numEntries(); j++) {
        double m=dat.get(j)->getRealValue(var.GetName());
        if ( m < var.getMin() || m > var.getMax())  continue;
        // Find probability density and hence probability
        var.setVal(m);
        double prb = var.getBinWidth(0)*pdf.getVal(var);
        prbSum+=prb;

        dat.get(j);
        nEvt=dat.weight();

        double mubin=nTot*prb;
        double contrib(0.);
        if (nEvt < 1) contrib = mubin;
        else contrib=mubin-nEvt+nEvt*log(nEvt/mubin);
        totNLL+=contrib;

        if(prt) cout << "Bin " << j << " prob = " << prb << " nEvt = " << nEvt << ", mu = " << mubin << " contribution " << contrib << endl;
    }

    totNLL*=2.0;
    if(prt) cout << pdf.GetName() << " nTot = " << nTot << " 2NLL constant = " << totNLL << endl;

    return totNLL;
}

void draw_data_mgg(TString folderName,bool blind=true,float min=103,float max=160)
{
  TFile inputFile(folderName+"/data.root");
  
  const int nCat = 5;
  TString cats[5] = {"HighPt","Hbb","Zbb","HighRes","LowRes"};

  TCanvas cv;

  for(int iCat=0; iCat < nCat; iCat++) {

    RooWorkspace *ws  = (RooWorkspace*)inputFile.Get(cats[iCat]+"_mgg_workspace");
    RooFitResult* res = (RooFitResult*)ws->obj("fitresult_pdf_data");

    RooRealVar * mass = ws->var("mgg");
    mass->setRange("all",min,max);
    mass->setRange("blind",121,130);
    mass->setRange("low",106,121);
    mass->setRange("high",130,160);

    mass->setUnit("GeV");
    mass->SetTitle("m_{#gamma#gamma}");
    
    RooAbsPdf * pdf = ws->pdf("pdf");
    RooPlot *plot = mass->frame(min,max,max-min);
    plot->SetTitle("");
    
    RooAbsData* data = ws->data("data")->reduce(Form("mgg > %f && mgg < %f",min,max));
    double nTot = data->sumEntries();
    if(blind) data = data->reduce("mgg < 121 || mgg>130");
    double nBlind = data->sumEntries();
    double norm = nTot/nBlind; //normalization for the plot
    
    data->plotOn(plot);
    pdf->plotOn(plot,RooFit::NormRange( "low,high" ),RooFit::Range("Full"),RooFit::LineWidth(0.1) );
    plot->Print();

    //add the fix error band
    RooCurve* c = plot->getCurve("pdf_Norm[mgg]_Range[Full]_NormRange[Full]");
    const int Nc = c->GetN();
    //TGraphErrors errfix(Nc);
    //TGraphErrors errfix2(Nc);
    TGraphAsymmErrors errfix(Nc);
    TGraphAsymmErrors errfix2(Nc);
    Double_t *x = c->GetX();
    Double_t *y = c->GetY();
    double NtotalFit = ws->var("Nbkg1")->getVal()*ws->var("Nbkg1")->getVal() + ws->var("Nbkg2")->getVal()*ws->var("Nbkg2")->getVal();
    for( int i = 0; i < Nc; i++ )
      {
	errfix.SetPoint(i,x[i],y[i]);
	errfix2.SetPoint(i,x[i],y[i]);
	mass->setVal(x[i]);      
	double shapeErr = pdf->getPropagatedError(*res)*NtotalFit;
	//double totalErr = TMath::Sqrt( shapeErr*shapeErr + y[i] );
	//total normalization error
	double totalErr = TMath::Sqrt( shapeErr*shapeErr + y[i]*y[i]/NtotalFit ); 
	if ( y[i] - totalErr > .0 )
	  {
	    errfix.SetPointError(i, 0, 0, totalErr, totalErr );
	  }
	else
	  {
	    errfix.SetPointError(i, 0, 0, y[i] - 0.01, totalErr );
	  }
	//2sigma
	if ( y[i] -  2.*totalErr > .0 )
	  {
	    errfix2.SetPointError(i, 0, 0, 2.*totalErr,  2.*totalErr );
	  }
	else
	  {
	    errfix2.SetPointError(i, 0, 0, y[i] - 0.01,  2.*totalErr );
	  }
	/*
	std::cout << x[i] << " " << y[i] << " "
		  << " ,pdf get Val: " << pdf->getVal()
		  << " ,pdf get Prop Err: " << pdf->getPropagatedError(*res)*NtotalFit
		  << " stat uncertainty: " << TMath::Sqrt(y[i]) << " Ntot: " << NtotalFit <<  std::endl;
	*/
      }
    errfix.SetFillColor(kYellow);
    errfix2.SetFillColor(kGreen);


    //pdf->plotOn(plot,RooFit::NormRange( "low,high" ),RooFit::FillColor(kGreen),RooFit::Range("Full"), RooFit::VisualizeError(*res,2.0,kFALSE));
    //pdf->plotOn(plot,RooFit::NormRange( "low,high" ),RooFit::FillColor(kYellow),RooFit::Range("Full"), RooFit::VisualizeError(*res,1.0,kFALSE));
    //pdf->plotOn(plot,RooFit::NormRange( "low,high" ),RooFit::FillColor(kGreen),RooFit::Range("Full"), RooFit::VisualizeError(*res,2.0,kTRUE));
    //pdf->plotOn(plot,RooFit::NormRange( "low,high" ),RooFit::FillColor(kYellow),RooFit::Range("Full"), RooFit::VisualizeError(*res,1.0,kTRUE));
    plot->addObject(&errfix,"4");
    plot->addObject(&errfix2,"4");
    plot->addObject(&errfix,"4");
    data->plotOn(plot);
    TBox blindBox(121,plot->GetMinimum()-(plot->GetMaximum()-plot->GetMinimum())*0.015,130,plot->GetMaximum());
    blindBox.SetFillColor(kGray);
    if(blind) {
      plot->addObject(&blindBox);
      pdf->plotOn(plot,RooFit::NormRange( "low,high" ),RooFit::FillColor(kGreen),RooFit::Range("Full"), RooFit::VisualizeError(*res,2.0,kTRUE));
      pdf->plotOn(plot,RooFit::NormRange( "low,high" ),RooFit::FillColor(kYellow),RooFit::Range("Full"), RooFit::VisualizeError(*res,1.0,kTRUE));
    }
    //plot->addObject(&errfix,"4");
//.........这里部分代码省略.........

//put very small data entries in a binned dataset to avoid unphysical pdfs, specifically for H->ZZ->4l
RooDataSet* makeData(RooDataSet* orig, RooSimultaneous* simPdf, const RooArgSet* observables, RooRealVar* firstPOI, double mass, double& mu_min)
{

  double max_soverb = 0;

  mu_min = -10e9;

  map<string, RooDataSet*> data_map;
  firstPOI->setVal(0);
  RooCategory* cat = (RooCategory*)&simPdf->indexCat();
  TList* datalist = orig->split(*(RooAbsCategory*)cat, true);
  TIterator* dataItr = datalist->MakeIterator();
  RooAbsData* ds;
  RooRealVar* weightVar = new RooRealVar("weightVar","weightVar",1);
  while ((ds = (RooAbsData*)dataItr->Next()))
  {
    string typeName(ds->GetName());
    cat->setLabel(typeName.c_str());
    RooAbsPdf* pdf = simPdf->getPdf(typeName.c_str());
    cout << "pdf: " << pdf << endl;
    RooArgSet* obs = pdf->getObservables(observables);
    cout << "obs: " << obs << endl;

    RooArgSet obsAndWeight(*obs, *weightVar);
    obsAndWeight.add(*cat);
    stringstream datasetName;
    datasetName << "newData_" << typeName;
    RooDataSet* thisData = new RooDataSet(datasetName.str().c_str(),datasetName.str().c_str(), obsAndWeight, WeightVar(*weightVar));

    RooRealVar* firstObs = (RooRealVar*)obs->first();
    //int ibin = 0;
    int nrEntries = ds->numEntries();
    for (int ib=0;ib<nrEntries;ib++)
    {
      const RooArgSet* event = ds->get(ib);
      const RooRealVar* thisObs = (RooRealVar*)event->find(firstObs->GetName());
      firstObs->setVal(thisObs->getVal());

      firstPOI->setVal(0);
      double b = pdf->expectedEvents(*firstObs)*pdf->getVal(obs);
      firstPOI->setVal(1);
      double s = pdf->expectedEvents(*firstObs)*pdf->getVal(obs) - b;

      if (s > 0)
      {
	mu_min = max(mu_min, -b/s);
	double soverb = s/b;
	if (soverb > max_soverb)
	{
	  max_soverb = soverb;
	  cout << "Found new max s/b: " << soverb << " in pdf " << pdf->GetName() << " at m = " << thisObs->getVal() << endl;
	}
      }

      if (b == 0 && s != 0)
      {
	cout << "Expecting non-zero signal and zero bg at m=" << firstObs->getVal() << " in pdf " << pdf->GetName() << endl;
      }
      if (s+b <= 0) 
      {
	cout << "expecting zero" << endl;
	continue;
      }


      double weight = ds->weight();
      if ((typeName.find("ATLAS_H_4mu") != string::npos || 
	   typeName.find("ATLAS_H_4e") != string::npos ||
	   typeName.find("ATLAS_H_2mu2e") != string::npos ||
	   typeName.find("ATLAS_H_2e2mu") != string::npos) && fabs(firstObs->getVal() - mass) < 10 && weight == 0)
      {
	cout << "adding event: " << firstObs->getVal() << endl;
	thisData->add(*event, pow(10., -9.));
      }
      else
      {
	//weight = max(pow(10.0, -9), weight);
	thisData->add(*event, weight);
      }
    }



    data_map[string(ds->GetName())] = (RooDataSet*)thisData;
  }

  
  RooDataSet* newData = new RooDataSet("newData","newData",RooArgSet(*observables, *weightVar), 
				       Index(*cat), Import(data_map), WeightVar(*weightVar));

  orig->Print();
  newData->Print();
  //newData->tree()->Scan("*");
  return newData;

}

void OneSidedFrequentistUpperLimitWithBands_intermediate(const char* infile = "",
					    const char* workspaceName = "combined",
					    const char* modelConfigName = "ModelConfig",
					    const char* dataName = "obsData"){


  double confidenceLevel=0.95;
  // degrade/improve number of pseudo-experiments used to define the confidence belt.  
  // value of 1 corresponds to default number of toys in the tail, which is 50/(1-confidenceLevel)
  double additionalToysFac = 1.;  
  int nPointsToScan = 30; // number of steps in the parameter of interest 
  int nToyMC = 100; // number of toys used to define the expected limit and band

  TStopwatch t;
  t.Start();
  /////////////////////////////////////////////////////////////
  // First part is just to access a user-defined file 
  // or create the standard example file if it doesn't exist
  ////////////////////////////////////////////////////////////
  const char* filename = "";
  if (!strcmp(infile,""))
    filename = "results/example_combined_GaussExample_model.root";
  else
    filename = infile;
  // Check if example input file exists
  TFile *file = TFile::Open(filename);

  // if input file was specified byt not found, quit
  if(!file && strcmp(infile,"")){
    cout <<"file not found" << endl;
    return;
  } 

  // if default file not found, try to create it
  if(!file ){
    // Normally this would be run on the command line
    cout <<"will run standard hist2workspace example"<<endl;
    gROOT->ProcessLine(".! prepareHistFactory .");
    gROOT->ProcessLine(".! hist2workspace config/example.xml");
    cout <<"\n\n---------------------"<<endl;
    cout <<"Done creating example input"<<endl;
    cout <<"---------------------\n\n"<<endl;
  }

  // now try to access the file again
  file = TFile::Open(filename);
  if(!file){
    // if it is still not there, then we can't continue
    cout << "Not able to run hist2workspace to create example input" <<endl;
    return;
  }

  
  /////////////////////////////////////////////////////////////
  // Now get the data and workspace
  ////////////////////////////////////////////////////////////

  // get the workspace out of the file
  RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);
  if(!w){
    cout <<"workspace not found" << endl;
    return;
  }

  // get the modelConfig out of the file
  ModelConfig* mc = (ModelConfig*) w->obj(modelConfigName);

  // get the modelConfig out of the file
  RooAbsData* data = w->data(dataName);

  // make sure ingredients are found
  if(!data || !mc){
    w->Print();
    cout << "data or ModelConfig was not found" <<endl;
    return;
  }

  cout << "Found data and ModelConfig:" <<endl;
  mc->Print();

  /////////////////////////////////////////////////////////////
  // Now get the POI for convenience
  // you may want to adjust the range of your POI
  ////////////////////////////////////////////////////////////
  RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
  //  firstPOI->setMin(0);
  //  firstPOI->setMax(10);

  /////////////////////////////////////////////
  // create and use the FeldmanCousins tool
  // to find and plot the 95% confidence interval
  // on the parameter of interest as specified
  // in the model config
  // REMEMBER, we will change the test statistic
  // so this is NOT a Feldman-Cousins interval
  FeldmanCousins fc(*data,*mc);
  fc.SetConfidenceLevel(confidenceLevel); 
  fc.AdditionalNToysFactor(additionalToysFac); // improve sampling that defines confidence belt
  //  fc.UseAdaptiveSampling(true); // speed it up a bit, but don't use for expectd limits
  fc.SetNBins(nPointsToScan); // set how many points per parameter of interest to scan
//.........这里部分代码省略.........

void StandardFeldmanCousinsDemo(const char* infile = "",
                                const char* workspaceName = "combined",
                                const char* modelConfigName = "ModelConfig",
                                const char* dataName = "obsData"){

   // -------------------------------------------------------
   // First part is just to access a user-defined file
   // or create the standard example file if it doesn't exist
   const char* filename = "";
   if (!strcmp(infile,"")) {
      filename = "results/example_combined_GaussExample_model.root";
      bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
      // if file does not exists generate with histfactory
      if (!fileExist) {
#ifdef _WIN32
         cout << "HistFactory file cannot be generated on Windows - exit" << endl;
         return;
#endif
         // Normally this would be run on the command line
         cout <<"will run standard hist2workspace example"<<endl;
         gROOT->ProcessLine(".! prepareHistFactory .");
         gROOT->ProcessLine(".! hist2workspace config/example.xml");
         cout <<"\n\n---------------------"<<endl;
         cout <<"Done creating example input"<<endl;
         cout <<"---------------------\n\n"<<endl;
      }

   }
   else
      filename = infile;

   // Try to open the file
   TFile *file = TFile::Open(filename);

   // if input file was specified byt not found, quit
   if(!file ){
      cout <<"StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
      return;
   }


   // -------------------------------------------------------
   // Tutorial starts here
   // -------------------------------------------------------

   // get the workspace out of the file
   RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);
   if(!w){
      cout <<"workspace not found" << endl;
      return;
   }

   // get the modelConfig out of the file
   ModelConfig* mc = (ModelConfig*) w->obj(modelConfigName);

   // get the modelConfig out of the file
   RooAbsData* data = w->data(dataName);

   // make sure ingredients are found
   if(!data || !mc){
      w->Print();
      cout << "data or ModelConfig was not found" <<endl;
      return;
   }

   // -------------------------------------------------------
   // create and use the FeldmanCousins tool
   // to find and plot the 95% confidence interval
   // on the parameter of interest as specified
   // in the model config
   FeldmanCousins fc(*data,*mc);
   fc.SetConfidenceLevel(0.95); // 95% interval
   //fc.AdditionalNToysFactor(0.1); // to speed up the result
   fc.UseAdaptiveSampling(true); // speed it up a bit
   fc.SetNBins(10); // set how many points per parameter of interest to scan
   fc.CreateConfBelt(true); // save the information in the belt for plotting

   // Since this tool needs to throw toy MC the PDF needs to be
   // extended or the tool needs to know how many entries in a dataset
   // per pseudo experiment.
   // In the 'number counting form' where the entries in the dataset
   // are counts, and not values of discriminating variables, the
   // datasets typically only have one entry and the PDF is not
   // extended.
   if(!mc->GetPdf()->canBeExtended()){
      if(data->numEntries()==1)
         fc.FluctuateNumDataEntries(false);
      else
         cout <<"Not sure what to do about this model" <<endl;
   }

   // We can use PROOF to speed things along in parallel
   //  ProofConfig pc(*w, 1, "workers=4", kFALSE);
   //  ToyMCSampler*  toymcsampler = (ToyMCSampler*) fc.GetTestStatSampler();
   //  toymcsampler->SetProofConfig(&pc); // enable proof


   // Now get the interval
   PointSetInterval* interval = fc.GetInterval();
   ConfidenceBelt* belt = fc.GetConfidenceBelt();
//.........这里部分代码省略.........

void Raa3S_Workspace(const char* name_pbpb="chad_ws_fits/centFits/ws_PbPbData_262548_263757_0cent10_0.0pt50.0_0.0y2.4.root", const char* name_pp="chad_ws_fits/centFits/ws_PPData_262157_262328_-1cent1_0.0pt50.0_0.0y2.4.root", const char* name_out="fitresult_combo.root"){

   //TFile File(filename);

   //RooWorkspace * ws = test_combine(name_pbpb, name_pp);

   TFile *f = new TFile("fitresult_combo_333.root") ;
   RooWorkspace * ws1 = (RooWorkspace*) f->Get("wcombo");

   //File.GetObject("wcombo", ws);
   ws1->Print();
   RooAbsData * data = ws1->data("data"); //dataOS, dataSS

   // RooDataSet * US_data = (RooDataSet*) data->reduce( "QQsign == QQsign::PlusMinus");
   // US_data->SetName("US_data");
   // ws->import(* US_data);
   // RooDataSet * hi_data = (RooDataSet*) US_data->reduce("dataCat == dataCat::hi");
   // hi_data->SetName("hi_data");
   // ws->import(* hi_data);
   // hi_data->Print();

   RooRealVar* raa3 = new RooRealVar("raa3","R_{AA}(#Upsilon (3S))",0.5,-1,1);
   RooRealVar* leftEdge = new RooRealVar("leftEdge","leftEdge",0);
   RooRealVar* rightEdge = new RooRealVar("rightEdge","rightEdge",1);
   RooGenericPdf step("step", "step", "(@0 >= @1) && (@0 < @2)", RooArgList(*raa3, *leftEdge, *rightEdge));
   ws1->import(step);
   ws1->factory( "Uniform::flat(raa3)" );

   //pp Luminosities, Taa and efficiency ratios Systematics

   ws1->factory( "Taa_hi[5.662e-9]" );
   ws1->factory( "Taa_kappa[1.062]" ); // was 1.057
   ws1->factory( "expr::alpha_Taa('pow(Taa_kappa,beta_Taa)',Taa_kappa,beta_Taa[0,-5,5])" );
   ws1->factory( "prod::Taa_nom(Taa_hi,alpha_Taa)" );
   ws1->factory( "Gaussian::constr_Taa(beta_Taa,glob_Taa[0,-5,5],1)" );

   ws1->factory( "lumipp_hi[5.4]" );
   ws1->factory( "lumipp_kappa[1.037]" ); // was 1.06
   ws1->factory( "expr::alpha_lumipp('pow(lumipp_kappa,beta_lumipp)',lumipp_kappa,beta_lumipp[0,-5,5])" );
   ws1->factory( "prod::lumipp_nom(lumipp_hi,alpha_lumipp)" );
   ws1->factory( "Gaussian::constr_lumipp(beta_lumipp,glob_lumipp[0,-5,5],1)" );

   // ws->factory( "effRat1[1]" );
   // ws->factory( "effRat2[1]" );
   ws1->factory( "effRat3_hi[0.95]" );
   ws1->factory( "effRat_kappa[1.054]" );
   ws1->factory( "expr::alpha_effRat('pow(effRat_kappa,beta_effRat)',effRat_kappa,beta_effRat[0,-5,5])" );
   // ws->factory( "prod::effRat1_nom(effRat1_hi,alpha_effRat)" );
   ws1->factory( "Gaussian::constr_effRat(beta_effRat,glob_effRat[0,-5,5],1)" );
   // ws->factory( "prod::effRat2_nom(effRat2_hi,alpha_effRat)" );
   ws1->factory( "prod::effRat3_nom(effRat3_hi,alpha_effRat)" );
   //  
   ws1->factory("Nmb_hi[1.161e9]");
   ws1->factory("prod::denominator(Taa_nom,Nmb_hi)");
   ws1->factory( "expr::lumiOverTaaNmbmodified('lumipp_nom/denominator',lumipp_nom,denominator)");
   RooAbsReal *lumiOverTaaNmbmodified = ws1->function("lumiOverTaaNmbmodified"); //RooFormulaVar *lumiOverTaaNmbmodified = ws->function("lumiOverTaaNmbmodified");
   //  
   //  RooRealVar *raa1 = ws->var("raa1");
   //  RooRealVar* nsig1_pp = ws->var("nsig1_pp");
   //  RooRealVar* effRat1 = ws->function("effRat1_nom");
   //  RooRealVar *raa2 = ws->var("raa2");
   //  RooRealVar* nsig2_pp = ws->var("nsig2_pp");
   //  RooRealVar* effRat2 = ws->function("effRat2_nom");
   RooRealVar* nsig3_pp = ws1->var("R_{#frac{3S}{1S}}_pp"); //RooRealVar* nsig3_pp = ws->var("N_{#Upsilon(3S)}_pp");
   cout << nsig3_pp << endl;
   RooAbsReal* effRat3 = ws1->function("effRat3_nom"); //RooRealVar* effRat3 = ws->function("effRat3_nom");
   //  
   //  RooFormulaVar nsig1_hi_modified("nsig1_hi_modified", "@0*@1*@3/@2", RooArgList(*raa1, *nsig1_pp, *lumiOverTaaNmbmodified, *effRat1));
   //  ws->import(nsig1_hi_modified);
   //  RooFormulaVar nsig2_hi_modified("nsig2_hi_modified", "@0*@1*@3/@2", RooArgList(*raa2, *nsig2_pp, *lumiOverTaaNmbmodified, *effRat2));
   //  ws->import(nsig2_hi_modified);
   RooFormulaVar nsig3_hi_modified("nsig3_hi_modified", "@0*@1*@3/@2", RooArgList(*raa3, *nsig3_pp, *lumiOverTaaNmbmodified, *effRat3));
   ws1->import(nsig3_hi_modified);

   //  // background yield with systematics
   ws1->factory( "nbkg_hi_kappa[1.10]" );
   ws1->factory( "expr::alpha_nbkg_hi('pow(nbkg_hi_kappa,beta_nbkg_hi)',nbkg_hi_kappa,beta_nbkg_hi[0,-5,5])" );
   ws1->factory( "SUM::nbkg_hi_nom(alpha_nbkg_hi*bkgPdf_hi)" );
   ws1->factory( "Gaussian::constr_nbkg_hi(beta_nbkg_hi,glob_nbkg_hi[0,-5,5],1)" );
   RooAbsPdf* sig1S_hi = ws1->pdf("sig1S_hi"); //RooAbsPdf* sig1S_hi = ws->pdf("cbcb_hi");
   RooAbsPdf* sig2S_hi = ws1->pdf("sig2S_hi");
   RooAbsPdf* sig3S_hi = ws1->pdf("sig3S_hi");
   RooAbsPdf* LSBackground_hi = ws1->pdf("nbkg_hi_nom");
   RooRealVar* nsig1_hi = ws1->var("N_{#Upsilon(1S)}_hi");
   RooRealVar* nsig2_hi = ws1->var("R_{#frac{2S}{1S}}_hi");
   RooAbsReal* nsig3_hi = ws1->function("nsig3_hi_modified"); //RooFormulaVar* nsig3_hi = ws->function("nsig3_hi_modified");
   cout << nsig1_hi << " " << nsig2_hi << " " << nsig3_pp << endl;
   RooRealVar* norm_nbkg_hi = ws1->var("n_{Bkgd}_hi");

   RooArgList pdfs_hi( *sig1S_hi,*sig2S_hi,*sig3S_hi, *LSBackground_hi);
   RooArgList norms_hi(*nsig1_hi,*nsig2_hi,*nsig3_hi, *norm_nbkg_hi);

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

   ws1->factory( "nbkg_pp_kappa[1.03]" );
   ws1->factory( "expr::alpha_nbkg_pp('pow(nbkg_pp_kappa,beta_nbkg_pp)',nbkg_pp_kappa,beta_nbkg_pp[0,-5,5])" );
   ws1->factory( "SUM::nbkg_pp_nom(alpha_nbkg_pp*bkgPdf_pp)" );
   ws1->factory( "Gaussian::constr_nbkg_pp(beta_nbkg_pp,glob_nbkg_pp[0,-5,5],1)" );
   RooAbsPdf* sig1S_pp = ws1->pdf("sig1S_pp"); //RooAbsPdf* sig1S_pp = ws1->pdf("cbcb_pp");
   RooAbsPdf* sig2S_pp = ws1->pdf("sig2S_pp");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    gROOT->ProcessLine(".! prepareHistFactory .");
    gROOT->ProcessLine(".! hist2workspace config/example.xml");
    cout <<"\n\n---------------------"<<endl;
    cout <<"Done creating example input"<<endl;
    cout <<"---------------------\n\n"<<endl;
  }

  // now try to access the file again
  file = TFile::Open(filename);
  if(!file){
    // if it is still not there, then we can't continue
    cout << "Not able to run hist2workspace to create example input" <<endl;
    return;
  }

  
  /////////////////////////////////////////////////////////////
  // Tutorial starts here
  ////////////////////////////////////////////////////////////

  // get the workspace out of the file
  RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);
  if(!w){
    cout <<"workspace not found" << endl;
    return;
  }
  w->Print();

  // get the modelConfig out of the file
  ModelConfig* sbModel = (ModelConfig*) w->obj(modelSBName);


  // get the modelConfig out of the file
  RooAbsData* data = w->data(dataName);

  // make sure ingredients are found
  if(!data || !sbModel){
    w->Print();
    cout << "data or ModelConfig was not found" <<endl;
    return;
  }
  // make b model
  ModelConfig* bModel = (ModelConfig*) w->obj(modelBName);


   // case of no systematics
   // remove nuisance parameters from model
   if (noSystematics) { 
      const RooArgSet * nuisPar = sbModel->GetNuisanceParameters();
      if (nuisPar && nuisPar->getSize() > 0) { 
         std::cout << "StandardHypoTestInvDemo" << "  -  Switch off all systematics by setting them constant to their initial values" << std::endl;
         RooStats::SetAllConstant(*nuisPar);
      }
      if (bModel) { 
         const RooArgSet * bnuisPar = bModel->GetNuisanceParameters();
         if (bnuisPar) 
            RooStats::SetAllConstant(*bnuisPar);
      }
   }


  if (!bModel ) {
      Info("StandardHypoTestInvDemo","The background model %s does not exist",modelBName);
      Info("StandardHypoTestInvDemo","Copy it from ModelConfig %s and set POI to zero",modelSBName);
      bModel = (ModelConfig*) sbModel->Clone();
      bModel->SetName(TString(modelSBName)+TString("B_only"));