C++ RooAbsReal类代码示例

   RooFitResult* GenericModel::fitTo(RooDataSet* data)
   {
      // Perform fit of the pseudo-PDF to the data
      // On multi-core machines, this automatically uses all available processor cores

      SafeDelete(fLastFit);
#ifdef WITH_MULTICORE_CPU
      fLastFit = fModelPseudoPDF->fitTo(*data, Save(), NumCPU(WITH_MULTICORE_CPU));
#else
      fLastFit = fModelPseudoPDF->fitTo(*data, Save());
#endif

      SafeDelete(fParamDataHist);
      fParamDataHist = new RooDataHist("params", "params", GetParameters());

      // store weights of component pdfs => distribution of parameters
      fWeights.removeAll();
      const RooArgList& coefs = fModelPseudoPDF->coefList();
      for (int i = 0; i < GetNumberOfDataSets(); i++) {
         RooAbsReal* coef = (RooAbsReal*)coefs.at(i);
         RooRealVar w(Form("w%d", i), Form("Fitted weight of kernel#%d", i), coef->getVal());
         if (coef->InheritsFrom(RooRealVar::Class())) {
            w.setError(((RooRealVar*)coef)->getError());
         } else {
            w.setError(coef->getPropagatedError(*fLastFit));
         }
         fWeights.addClone(w);
         fParamDataHist->set(*GetParametersForDataset(i), w.getVal(), w.getError());
      }

      SafeDelete(fParameterPDF);
      fParameterPDF = new RooHistPdf("paramPDF", "paramPDF", GetParameters(), *fParamDataHist);

      return fLastFit;
   }

void rf313_paramranges()
{

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

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

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



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

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

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

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



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

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

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



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

  return ;
}

// grab the initial normalization from a datacard converted in workspace
// with: scripts/text2workspace.py -b -o model.root datacards/hww-12.1fb.mH125.comb_0j1j2j_shape.txt
void fillInitialNorms(RooArgSet *args, std::map<std::string, std::pair<double,double> > &vals, std::string workspace){
  TFile *fw_ =  TFile::Open(workspace.c_str());
  RooWorkspace *ws = (RooWorkspace*)fw_->Get("w");
  TIterator* iter(args->createIterator());
  for (TObject *a = iter->Next(); a != 0; a = iter->Next()) {
    RooAbsReal *rar = (RooAbsReal*)ws->obj(a->GetName());
    std::string name = rar->GetName();
    std::pair<double,double> valE(rar->getVal(),0.0);
    vals.insert( std::pair<std::string,std::pair<double ,double> > (name,valE)) ;
  }
}

double NormalizedIntegral(RooAbsPdf & function, RooRealVar & integrationVar, double lowerLimit, double upperLimit){

  integrationVar.setRange("integralRange",lowerLimit,upperLimit);
  RooAbsReal* integral = function.createIntegral(integrationVar,NormSet(integrationVar),Range("integralRange"));


  double normlizedIntegralValue = integral->getVal();

  //  cout<<normlizedIntegralValue<<endl;


  return normlizedIntegralValue;


}

THSEventsPDF::THSEventsPDF(const char *name, const char *title, 
			   RooAbsReal& _x,
			   RooAbsReal& _alpha,
			   RooAbsReal& _offset,
			   RooAbsReal& _scale,
			   Int_t NAlphBins  ) :
  RooAbsPdf(name,title), 
  x("x","x",this,_x),
  offset("offset","offset",this,_offset),
  scale("scale","scale",this,_scale),
  alpha("alpha","alpha",this,_alpha),
  fx_off(0),
  falpha(0),
  fHist(0),
  fHistPdf(0),
  fRHist(0),
  fWeightHist(0),
  fTree(0)
{ 
  
  RooRealVar *rx=dynamic_cast<RooRealVar*>(&_x);
  RooRealVar *ra=dynamic_cast<RooRealVar*>(&_alpha);
  RooRealVar *rs=dynamic_cast<RooRealVar*>(&_scale);
  RooRealVar *ro=dynamic_cast<RooRealVar*>(&_offset);
  //Work out number of bins for x axis
  //Take bin width as being 1/10 of the alpha range (arbitrary!) or 100 whichever is larger
  // NbinX=40;
  Double_t rsmin=1;
  fOldScale=rs->getVal();
  if(rs->getMin()) rsmin=rs->getMin();
  else cout<<"THSEventsPDF::THSEventsPDF Warning no scale minimum set take = 1"<<endl;
  Double_t mid=(rx->getMax()+rx->getMin())/2;
  Double_t diff=(rx->getMax()-rx->getMin())/2;
  Double_t rMin=mid-diff/rsmin - ro->getMin(); //additional range for possible tranformation or scaling
  Double_t rMax=mid+diff/rsmin + ro->getMax();
  Int_t NbinX=200/rsmin;
  cout<<GetName()<<" hist ranges "<<rMin<<" to "<<rMax<<endl;
  if(NbinX<10) NbinX=10;
   fRHist=new TH2F(TString("hmc_model_")+_x.GetName()+name,TString("MC model for ")+_x.GetName(),NbinX,rMin,rMax,NAlphBins,ra->getMin(),ra->getMax());
  fRHist->Sumw2();
  //  fRHist=new TH2F(TString("hmc_model_")+_x.GetName()+name,TString("MC model for ")+_x.GetName(),NbinX,rx->getMin(),rx->getMax(),NAlphBins,ra->getMin(),ra->getMax());
  fx=new RooRealVar(_x.GetName(),"Vx",0,x.min(),x.max());
  fx_off=new RooRealVar(_x.GetName(),"Vx_off",0,rMin,rMax);
  falpha=new RooRealVar("Valpha","Valpha",0,alpha.min(),alpha.max());
  
  fNWdim=0;
} 

void fit( RooAbsReal & chi2, int numberOfBins, const char * outFileNameWithFitResult ){
  TFile * out_root_file = new TFile(outFileNameWithFitResult , "recreate");
  RooMinuit m_tot(chi2) ;
  m_tot.migrad();
  // m_tot.hesse();
  RooFitResult* r_chi2 = m_tot.save() ;
  cout << "==> Chi2 Fit results " << endl ;
  r_chi2->Print("v") ;
  //  r_chi2->floatParsFinal().Print("v") ;
  int NumberOfFreeParameters =  r_chi2->floatParsFinal().getSize() ;
  for (int i =0; i< NumberOfFreeParameters; ++i){
    r_chi2->floatParsFinal()[i].Print();
  }
  cout<<"chi2:" <<chi2.getVal() << ", numberOfBins: " << numberOfBins  << ", NumberOfFreeParameters: " << NumberOfFreeParameters << endl; 
  cout<<"Normalized Chi2   = " << chi2.getVal()/ (numberOfBins - NumberOfFreeParameters)<<endl; 
  r_chi2->Write( ) ;
  delete out_root_file;
}

///
/// Make an Asimov toy: set all observables set to truth values.
/// The Asimov point needs to be loaded in the combiner before.
/// \param c - combiner which should be set to an asimov toy
///
void GammaComboEngine::setAsimovObservables(Combiner* c)
{
    if ( !c->isCombined() ) {
        cout << "GammaComboEngine::setAsimovObservables() : ERROR : Can't set an Asimov toy before "
             "the combiner is combined. Call combine() first." << endl;
        exit(1);
    }

    // set observables to asimov values in workspace
    RooWorkspace* w = c->getWorkspace();
    TIterator* itObs = c->getObservables()->createIterator();
    while(RooRealVar* pObs = (RooRealVar*) itObs->Next()) {
        // get theory name from the observable name
        TString pThName = pObs->GetName();
        pThName.ReplaceAll("obs","th");
        // get the theory relation
        RooAbsReal* th = w->function(pThName);
        if ( th==0 ) {
            cout << "GammaComboEngine::setAsimovObservables() : ERROR : theory relation not found in workspace: " << pThName << endl;
            exit(1);
        }
        // set the observable to what the theory relation predicts
        pObs->setVal(th->getVal());
    }
    delete itObs;

    // write back the asimov values to the PDF object so that when
    // the PDF is printed, the asimov values show up
    for ( int i=0; i<c->getPdfs().size(); i++ ) {
        PDF_Abs* pdf = c->getPdfs()[i];
        pdf->setObservableSourceString("Asimov");
        TIterator* itObs = pdf->getObservables()->createIterator();
        while(RooRealVar* pObs = (RooRealVar*) itObs->Next()) {
            RooAbsReal* obs =  w->var(pObs->GetName());
            if ( obs==0 ) {
                cout << "GammaComboEngine::setAsimovObservables() : ERROR : observable not found in workspace: " << pObs->GetName() << endl;
                exit(1);
            }
            pdf->setObservable(pObs->GetName(), obs->getVal());
        }
        delete itObs;
    }
}

void Zbi_Zgamma() {

  // Make model for prototype on/off problem
  // Pois(x | s+b) * Pois(y | tau b )
  // for Z_Gamma, use uniform prior on b.
  RooWorkspace* w = new RooWorkspace("w",true);
  w->factory("Poisson::px(x[150,0,500],sum::splusb(s[0,0,100],b[100,0,300]))");
  w->factory("Poisson::py(y[100,0,500],prod::taub(tau[1.],b))");
  w->factory("Uniform::prior_b(b)");

  // construct the Bayesian-averaged model (eg. a projection pdf)
  // p'(x|s) = \int db p(x|s+b) * [ p(y|b) * prior(b) ]
  w->factory("PROJ::averagedModel(PROD::foo(px|b,py,prior_b),b)") ;

  // plot it, blue is averaged model, red is b known exactly
  RooPlot* frame = w->var("x")->frame() ;
  w->pdf("averagedModel")->plotOn(frame) ;
  w->pdf("px")->plotOn(frame,LineColor(kRed)) ;
  frame->Draw() ;

  // compare analytic calculation of Z_Bi
  // with the numerical RooFit implementation of Z_Gamma
  // for an example with x = 150, y = 100

  // numeric RooFit Z_Gamma
  w->var("y")->setVal(100);
  w->var("x")->setVal(150);
  RooAbsReal* cdf = w->pdf("averagedModel")->createCdf(*w->var("x"));
  cdf->getVal(); // get ugly print messages out of the way

  cout << "Hybrid p-value = " << cdf->getVal() << endl;
  cout << "Z_Gamma Significance  = " <<
    PValueToSignificance(1-cdf->getVal()) << endl;

  // analytic Z_Bi
  double Z_Bi = NumberCountingUtils::BinomialWithTauObsZ(150, 100, 1);
  std::cout << "Z_Bi significance estimation: " << Z_Bi << std::endl;

  // OUTPUT
  // Hybrid p-value = 0.999058
  // Z_Gamma Significance  = 3.10804
  // Z_Bi significance estimation: 3.10804
}

//#include <typeinfo.h>
void addFlatNuisances(std::string fi){
  gSystem->Load("libHiggsAnalysisCombinedLimit.so");
  TFile *fin = TFile::Open(fi.c_str());
  RooWorkspace *wspace = (RooWorkspace*)fin->Get("w_hmumu");

  wspace->Print("");

  RooStats::ModelConfig *mc = (RooStats::ModelConfig*)wspace->genobj("ModelConfig");
  RooArgSet *nuis = (RooArgSet*) mc->GetNuisanceParameters();
  std::cout << "Before...." << std::endl;
  nuis->Print();
  
  RooRealVar *mgg = (RooRealVar*)wspace->var("mmm");
  // Get all of the "flat" nuisances to be added to the nusiances:
  RooArgSet pdfs = (RooArgSet) wspace->allVars();
  RooAbsReal *pdf;
  TIterator *it_pdf = pdfs.createIterator();
  

  while ( (pdf=(RooAbsReal*)it_pdf->Next()) ){
  	  if (!(std::string(pdf->GetName()).find("zmod") != std::string::npos )) {
  	   if (!(std::string(pdf->GetName()).find("__norm") != std::string::npos )) {
	   	continue;
	   }
	  }
	  pdf->Print();
	  RooArgSet* pdfpars = (RooArgSet*)pdf->getParameters(RooArgSet(*mgg));
	  pdfpars->Print();

	  std::string newname_pdf = (std::string("unconst_")+std::string(pdf->GetName()));
	  wspace->import(*pdf,RooFit::RenameVariable(pdf->GetName(),newname_pdf.c_str()));
	  pdf->SetName(newname_pdf.c_str());
	  nuis->add(*pdf);
  }
 
  wspace->var("MH")->setVal(125.0);
  std::cout << "After..." << std::endl;
  nuis->Print();
  mc->SetNuisanceParameters(*nuis);
  //RooWorkspace *wspace_new = wspace->Clone();
  //mc->SetWorkspace(*wspace_new);
  //wspace_new->import(*mc,true);

  TFile *finew = new TFile((std::string(fin->GetName())+std::string("_unconst.root")).c_str(),"RECREATE");
  //wspace_new->SetName("w");
  finew->WriteTObject(wspace);
  finew->Close();
}

// get effective sigma from culmalative distribution function
pair<double,double> getEffSigma(RooRealVar *mass, RooAbsPdf *pdf, double wmin=110., double wmax=130., double step=0.002, double epsilon=1.e-4){

  RooAbsReal *cdf = pdf->createCdf(RooArgList(*mass));
  cout << "Computing effSigma...." << endl;
  TStopwatch sw;
  sw.Start();
  double point=wmin;
  vector<pair<double,double> > points;
  
  while (point <= wmax){
    mass->setVal(point);
    if (pdf->getVal() > epsilon){
      points.push_back(pair<double,double>(point,cdf->getVal())); 
    }
    point+=step;
  }
  double low = wmin;
  double high = wmax;
  double width = wmax-wmin;
  for (unsigned int i=0; i<points.size(); i++){
    for (unsigned int j=i; j<points.size(); j++){
      double wy = points[j].second - points[i].second;
      if (TMath::Abs(wy-0.683) < epsilon){
        double wx = points[j].first - points[i].first;
        if (wx < width){
          low = points[i].first;
          high = points[j].first;
          width=wx;
        }
      }
    }
  }
  sw.Stop();
  cout << "effSigma: [" << low << "-" << high << "] = " << width/2. << endl;
  cout << "\tTook: "; sw.Print();
  pair<double,double> result(low,high);
  return result;
}

prepDataFiles(){
//	TDirectory *theDr = (TDirectory*) myFile->Get("eleIDdir");///denom_pt/fit_eff_plots");
	//theDr->ls();
	int myIndex;	
	
	TSystemDirectory dir(thePath, thePath);
	TSystemFile *file;
	TString fname;
	TIter next(dir.GetListOfFiles());
	while ((file=(TSystemFile*)next())) {
		fname = file->GetName();
		if (fname.BeginsWith("TnP")&& fname.Contains("mc")) {
	
			ofstream myfile;

			TFile *myFile = new TFile(fname);
			TIter nextkey(myFile->GetListOfKeys());
			TKey *key;
			while (key = (TKey*)nextkey()) {
				TString theTypeClasse = key->GetClassName();
				TString theNomClasse = key->GetTitle();
				if ( theTypeClasse == "TDirectoryFile"){
					TDirectory *theDr = (TDirectory*) myFile->Get(theNomClasse);
					TIter nextkey2(theDr->GetListOfKeys());
					TKey *key2;
					while (key2 = (TKey*)nextkey2()) {
						TString theTypeClasse2 = key2->GetClassName();
						TString theNomClasse2 = key2->GetTitle();	
						myfile.open (theNomClasse2+".info");
						if ( theTypeClasse == "TDirectoryFile"){
							cout << "avant " << endl;
							TDirectory *theDr2 = (TDirectory*) myFile->Get(theNomClasse+"/"+theNomClasse2);
							cout << "apres " << endl;
							TIter nextkey3(theDr2->GetListOfKeys());
							TKey *key3;
							while (key3 = (TKey*)nextkey3()) {
								TString theTypeClasse3 = key3->GetClassName();
								TString theNomClasse3 = key3->GetTitle();	
								if ((theNomClasse3.Contains("FromMC"))) {

									TString localClasse3 = theNomClasse3;
									localClasse3.ReplaceAll("__","%");
									cout << "apres " << localClasse3 << endl;
									TObjArray* listBin = localClasse3.Tokenize('%');
									TString first = ((TObjString*)listBin->At(0))->GetString();
									TString second = ((TObjString*)listBin->At(2))->GetString();
									myfile << first;
									myfile << " " << second << " ";
									cout << "coucou la on va récupérer le rooFitResult " << endl;

									RooFitResult *theResults = (RooFitResult*) myFile->Get(theNomClasse+"/"+theNomClasse2+"/"+theNomClasse3+"/fitresults");
									theResults->Print();
									RooArgList theParam = theResults->floatParsFinal();
									int taille = theParam.getSize();
									for (int m = 0 ; m < taille ; m++){
										cout << "m=" << m << endl;
									RooAbsArg *theArg = (RooAbsArg*) theParam.at(m);
									RooAbsReal *theReal = (RooAbsReal*) theArg;
										myfile << theReal->getVal() << " " ;
									}		
															
									myfile << "\n";

								}
							}
						}
						myfile.close();

					}
			
				}
			}
			delete myFile;
		}
	
	}

}

//.........这里部分代码省略.........
    RooFIter iter = vars->fwdIterator();
    RooAbsArg* a;
    while( (a = iter.next()) ){
      if(string(a->GetName()).compare("mgg")==0) continue;
      static_cast<RooRealVar*>(a)->setConstant(kTRUE);
    }

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

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

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


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

    RooFitResult* sbres;
    RooAbsReal* nll;
    if(constrainMu) {
      fitModel.fitTo(data,RooFit::Strategy(0),RooFit::Extended(kTRUE),RooFit::ExternalConstraints(RooArgSet(HiggsMassConstraint)));
      sbres = fitModel.fitTo(data,RooFit::Strategy(2),RooFit::Save(kTRUE),RooFit::Extended(kTRUE),RooFit::ExternalConstraints(RooArgSet(HiggsMassConstraint)));
      nll = fitModel.createNLL(data,RooFit::NumCPU(4),RooFit::Extended(kTRUE),RooFit::ExternalConstraints(RooArgSet(HiggsMassConstraint)));
    } else {
      fitModel.fitTo(data,RooFit::Strategy(0),RooFit::Extended(kTRUE));
      sbres = fitModel.fitTo(data,RooFit::Strategy(2),RooFit::Save(kTRUE),RooFit::Extended(kTRUE));
      nll = fitModel.createNLL(data,RooFit::NumCPU(4),RooFit::Extended(kTRUE));
    }
    sbres->SetName(tag+"_sb_fitres");
    ws->import(*sbres);
    ws->import(fitModel);

    RooPlot *fmgg = mgg.frame();
    data.plotOn(fmgg);
    fitModel.plotOn(fmgg);
    ws->pdf("b_"+tag+"_ext")->plotOn(fmgg,RooFit::LineColor(kRed),RooFit::Range("Full"),RooFit::NormRange("Full"));
    fmgg->SetName(tag+"_frame");
    ws->import(*fmgg);
    delete fmgg;


    RooMinuit(*nll).migrad();

    RooPlot *fNs = Nsig.frame(0,25);
    fNs->SetName(tag+"_Nsig_pll");
    RooAbsReal *pll = nll->createProfile(Nsig);
    //nll->plotOn(fNs,RooFit::ShiftToZero(),RooFit::LineColor(kRed));
    pll->plotOn(fNs);
    ws->import(*fNs);

    delete fNs;

//.........这里部分代码省略.........
    error_curve[0] = new TGraph(2*nPoints);
    error_curve[1] = new TGraph(2*nPoints);
    
    error_curve[0]->SetFillStyle(1001);
    error_curve[1]->SetFillStyle(1001);
    
    error_curve[0]->SetFillColor(kGreen);
    error_curve[1]->SetFillColor(kYellow);
    
    error_curve[0]->SetLineColor(kGreen);
    error_curve[1]->SetLineColor(kYellow);
    
    if (plotBands) {
        RooDataHist pred2("pred2", "Prediction from SB", RooArgList(x), h_SR_Prediction2);

        error_curve[3]->SetFillStyle(1001);
        error_curve[4]->SetFillStyle(1001);
        
        error_curve[3]->SetFillColor(kGreen);
        error_curve[4]->SetFillColor(kYellow);
        
        error_curve[3]->SetLineColor(kGreen);
        error_curve[4]->SetLineColor(kYellow);
        
        error_curve[2]->SetLineColor(kBlue);
        error_curve[2]->SetLineWidth(3);
        
        double binSize = rebin;
        
        for (int i=0; i!=nPoints; ++i) {
            double x0,y0, x1,y1;
            error_curve[2]->GetPoint(i,x0,y0);
            
            RooAbsReal* nlim = new RooRealVar("nlim","y0",y0,-100000,100000);
            //double lowedge = x0 - (SR_hi - SR_lo)/double(2*nPoints);
            //double upedge = x0 + (SR_hi - SR_lo)/double(2*nPoints);
            
            double lowedge = x0 - binSize/2.;
            double upedge = x0 + binSize/2.;
            
            x.setRange("errRange",lowedge,upedge);
            
            RooExtendPdf* epdf = new RooExtendPdf("epdf","extpdf",bg, *nlim,"errRange");
            
            // Construct unbinned likelihood
            RooAbsReal* nll = epdf->createNLL(pred2,NumCPU(2));
            // Minimize likelihood w.r.t all parameters before making plots
            RooMinimizer* minim = new RooMinimizer(*nll);
            minim->setMinimizerType("Minuit2");
            minim->setStrategy(2);
            minim->setPrintLevel(-1);
            minim->migrad();
            
            minim->hesse();
            RooFitResult* result = minim->lastMinuitFit();
            double errm = nlim->getPropagatedError(*result);
            
            //std::cout<<x0<<" "<<lowedge<<" "<<upedge<<" "<<y0<<" "<<nlim->getVal()<<" "<<errm<<std::endl;
            
            error_curve[0]->SetPoint(i,x0,(y0-errm));
            error_curve[0]->SetPoint(2*nPoints-i-1,x0,y0+errm);
            
            error_curve[1]->SetPoint(i,x0,(y0-2*errm));
            error_curve[1]->SetPoint(2*nPoints-i-1,x0,(y0+2*errm));
            
            error_curve[3]->SetPoint(i,x0,-errm/sqrt(y0));

void rf308_normintegration2d()
{
  // S e t u p   m o d e l 
  // ---------------------

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

  // Create p.d.f. gaussx(x,-2,3), gaussy(y,2,2) 
  RooGaussian gx("gx","gx",x,RooConst(-2),RooConst(3)) ;
  RooGaussian gy("gy","gy",y,RooConst(+2),RooConst(2)) ;

  // Create gxy = gx(x)*gy(y)
  RooProdPdf gxy("gxy","gxy",RooArgSet(gx,gy)) ;



  // R e t r i e v e   r a w  &   n o r m a l i z e d   v a l u e s   o f   R o o F i t   p . d . f . s
  // --------------------------------------------------------------------------------------------------

  // Return 'raw' unnormalized value of gx
  cout << "gxy = " << gxy.getVal() << endl ;
  
  // Return value of gxy normalized over x _and_ y in range [-10,10]
  RooArgSet nset_xy(x,y) ;
  cout << "gx_Norm[x,y] = " << gxy.getVal(&nset_xy) << endl ;

  // Create object representing integral over gx
  // which is used to calculate  gx_Norm[x,y] == gx / gx_Int[x,y]
  RooAbsReal* igxy = gxy.createIntegral(RooArgSet(x,y)) ;
  cout << "gx_Int[x,y] = " << igxy->getVal() << endl ;

  // NB: it is also possible to do the following

  // Return value of gxy normalized over x in range [-10,10] (i.e. treating y as parameter)
  RooArgSet nset_x(x) ;
  cout << "gx_Norm[x] = " << gxy.getVal(&nset_x) << endl ;

  // Return value of gxy normalized over y in range [-10,10] (i.e. treating x as parameter)
  RooArgSet nset_y(y) ;
  cout << "gx_Norm[y] = " << gxy.getVal(&nset_y) << endl ;



  // I n t e g r a t e   n o r m a l i z e d   p d f   o v e r   s u b r a n g e
  // ----------------------------------------------------------------------------

  // Define a range named "signal" in x from -5,5
  x.setRange("signal",-5,5) ;
  y.setRange("signal",-3,3) ;
  
  // Create an integral of gxy_Norm[x,y] over x and y in range "signal"
  // This is the fraction of of p.d.f. gxy_Norm[x,y] which is in the
  // range named "signal"
  RooAbsReal* igxy_sig = gxy.createIntegral(RooArgSet(x,y),NormSet(RooArgSet(x,y)),Range("signal")) ;
  cout << "gx_Int[x,y|signal]_Norm[x,y] = " << igxy_sig->getVal() << endl ;




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

  // Create the cumulative distribution function of gx
  // i.e. calculate Int[-10,x] gx(x') dx'
  RooAbsReal* gxy_cdf = gxy.createCdf(RooArgSet(x,y)) ;
  
  // Plot cdf of gx versus x
  TH1* hh_cdf = gxy_cdf->createHistogram("hh_cdf",x,Binning(40),YVar(y,Binning(40))) ;
  hh_cdf->SetLineColor(kBlue) ;

  new TCanvas("rf308_normintegration2d","rf308_normintegration2d",600,600) ;
  gPad->SetLeftMargin(0.15) ; hh_cdf->GetZaxis()->SetTitleOffset(1.8) ; 
  hh_cdf->Draw("surf") ;

}

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");
//.........这里部分代码省略.........