C++ RooArgSet类代码示例

void AnalyzeToy::extract_signal()
{
   calculate_yield();
   MakeSpinSPlot splotter(toyData);
   splotter.addSpecies("signal",ws->pdf("model_signal_mass"),signalYield);
   splotter.addSpecies("background",ws->pdf("model_bkg_mass"),backgroundYield);
   splotter.addVariable(ws->var("mass"));
   splotter.calculate();
   RooDataSet *sweights = splotter.getSWeightDataSet();
   sweights->SetName("sweights");

   RooRealVar weight("weight","",-5.,5.);
   RooArgSet event;
   event.add(*mass);
   event.add(*cosT);
   event.add(weight);

   extractedData = new RooDataSet("extractedData","",event,WeightVar("weight"));

   Long64_t nEntries = toyData->numEntries();
   for(int i=0;i<nEntries;i++)
   {
      double weight_double=0;
      weight_double += sweights->get(i)->getRealValue("signal_sw");
//      weight_double += sweights->get(i)->getRealValue("background_sw");

      mass->setVal(toyData->get(i)->getRealValue("mass"));
      cosT->setVal(toyData->get(i)->getRealValue("cosT"));
      extractedData->add(event,weight_double);
   }
   delete toyData;
}

void reweight_eta_1d(TH1F* weight_eta,TH1F* weight_etao,TH2F*weight_etan,TH2F* weight_eta2o,TH2F* weight_etanr,TH2F*weight_eta2,RooDataSet **dset, RooDataSet *dsetdestination, int numvar){

  if (!(*dset)) return;

  TH1F *hnum = new TH1F("hnum","hnum",n_etabins_forreweighting,etabins_forreweighting);
 TH1F *hden = new TH1F("hden","hden",n_etabins_forreweighting,etabins_forreweighting);
//  TH1F *hnum = new TH1F("hnum","hnum",25,0.,2.5);
//  TH1F *hden = new TH1F("hden","hden",25,0.,2.5);
  hnum->Sumw2();
  hden->Sumw2();

  const char* etaname=Form("rooeta%d",numvar);

  for (int i=0; i<(*dset)->numEntries(); i++){
    hden->Fill(fabs((*dset)->get(i)->getRealValue(etaname)),(*dset)->store()->weight(i));
  }
  for (int i=0; i<dsetdestination->numEntries(); i++){
    hnum->Fill(fabs(dsetdestination->get(i)->getRealValue(etaname)),dsetdestination->store()->weight(i));
  }
  hnum->Scale(1.0/hnum->Integral());
  hden->Scale(1.0/hden->Integral());

  hnum->Divide(hden);
  TH1F *h = hnum;

  RooDataSet *newdset = new RooDataSet(**dset,Form("%s_etarew",(*dset)->GetName()));
  newdset->reset();
  for (int i=0; i<(*dset)->numEntries(); i++){
    RooArgSet args = *((*dset)->get(i));
    float oldw = (*dset)->store()->weight(i);
    float eta = args.getRealValue(etaname);
    float neww = oldw*h->GetBinContent(h->FindBin(fabs(eta)));
    if(debug){ 
		weight_eta->Fill(neww);	
		weight_etao->Fill(oldw);
		weight_etan->Fill(h->FindBin(fabs(eta)),neww);	
		weight_eta2o->Fill(h->FindBin(fabs(eta)),oldw);
		if(oldw!=0 && neww!=0)weight_etanr->Fill(h->FindBin(fabs(eta)),oldw/neww);
		else {weight_etanr->Fill(-10,1);}
	   // weight_pt2->Fill(pt,neww/oldw);
	   if(oldw!=0 && neww!=0)weight_eta2->Fill(fabs(eta),oldw/neww);
	   else {weight_eta2->Fill(-10,1);}
	}  
	  newdset->add(args,neww);
	  
	  }


  newdset->SetName((*dset)->GetName());
  newdset->SetTitle((*dset)->GetTitle());
  delete hnum; delete hden;

  RooDataSet *old_dset = *dset;
  *dset=newdset;
  std::cout << "Eta 1d rew: norm from " << old_dset->sumEntries() << " to " << newdset->sumEntries() << std::endl;

  delete old_dset;

};

void checkBestFitPoint(std::string workspace, std::string fitFile, bool splusb){
	
	// Open the ws file...
	TFile *fd_=0;
	TFile *fw_=0;
		
	gSystem->Load("$CMSSW_BASE/lib/$SCRAM_ARCH/libHiggsAnalysisCombinedLimit.so");
	gROOT->SetBatch(true);
	gStyle->SetOptFit(0);
	gStyle->SetOptStat(0);
	gStyle->SetPalette(1,0);

	fw_ =  TFile::Open(workspace.c_str());
	w   = (RooWorkspace*) fw_->Get("w");
  w->Print();
	RooDataSet *data = (RooDataSet*) w->data("data_obs");
  if (splusb) {
    mc_s = (RooStats::ModelConfig*)w->genobj("ModelConfig");
  } else {
    mc_s = (RooStats::ModelConfig*)w->genobj("ModelConfig_bonly");
  }
	std::cout << "make nll"<<std::endl;
	nll = mc_s->GetPdf()->createNLL(
		*data,RooFit::Constrain(*mc_s->GetNuisanceParameters())
		,RooFit::Extended(mc_s->GetPdf()->canBeExtended()));
	
	// Now get the best fit result
	fd_ =  TFile::Open(fitFile.c_str());
	RooFitResult *fit;
  if (splusb) {
    fit =(RooFitResult*)fd_->Get("fit_s");
  } else {
    fit =(RooFitResult*)fd_->Get("fit_b");
  }
	RooArgSet fitargs = fit->floatParsFinal();
	
	std::cout << "Got the best fit values" <<std::endl;		
	w->saveSnapshot("bestfitall",fitargs,true);
	
  TString filename;
  if (splusb) {
    filename = "minimum_s.pdf";
  } else {
    filename = "minimum_b.pdf";
  }
	// Now make the plots!	
	TCanvas *c = new TCanvas("c","",600,600);
	c->SaveAs((filename+"["));

	TIterator* iter(fitargs->createIterator());
        for (TObject *a = iter->Next(); a != 0; a = iter->Next()) {
                 RooRealVar *rrv = dynamic_cast<RooRealVar *>(a);      
                 std::string name = rrv->GetName();
		 TGraph *gr = graphLH(name,rrv->getError());
		 gr->Draw("ALP");
		 c->SaveAs((filename+"["));
	}
	c->SaveAs((filename+"]"));
}

///
/// Set all parameters to values found in
/// a provided fit result.
///
void PDF_Abs::loadExtParameters(RooFitResult *r)
{
	RooArgSet *tmp = new RooArgSet();
	tmp->add(r->floatParsFinal());
	tmp->add(r->constPars());
	setParameters(parameters, tmp);
	delete tmp;
}

void makeDataset( TString fname, TString tname, TString outfname ) {

  RooWorkspace *w = new RooWorkspace("w","w");

  w->factory( "Dst_M[1950.,2070.]" );
  w->factory( "D0_M[1810.,1920.]" );
  w->factory( "D0_LTIME_ps[0.00,5.]" );

  RooArgSet *observables = new RooArgSet();
  observables->add( *w->var("Dst_M") );
  observables->add( *w->var("D0_M") );
  observables->add( *w->var("D0_LTIME_ps") );
  w->defineSet("observables", *observables);

  w->var("Dst_M")->setBins(240);
  w->var("D0_M")->setBins(220);
  w->var("D0_LTIME_ps")->setBins(200);

  double Dst_M = -999.;
  double D0_M = -999.;
  double D0_LTIME_ps = -999.;

  TFile *tf = TFile::Open(fname);
  TTree *tree = (TTree*)tf->Get(tname);

  tree->SetBranchAddress( "Dst_M", &Dst_M );
  tree->SetBranchAddress( "D0_M" , &D0_M );
  tree->SetBranchAddress( "D0_LTIME_ps", &D0_LTIME_ps );

  RooDataSet *data = new RooDataSet("Data","Data",*observables);
  RooDataHist *dataH = new RooDataHist("DataHist","Data",*observables);

  for ( int ev=0; ev<tree->GetEntries(); ev++) {

    tree->GetEntry(ev);

    if ( ev%10000 == 0 ) cout << ev << " / " << tree->GetEntries() << endl;

    if ( Dst_M < w->var("Dst_M")->getMin() || Dst_M > w->var("Dst_M")->getMax() ) continue;
    if ( D0_M < w->var("D0_M")->getMin() || D0_M > w->var("D0_M")->getMax() ) continue;
    if ( D0_LTIME_ps < w->var("D0_LTIME_ps")->getMin() || D0_LTIME_ps > w->var("D0_LTIME_ps")->getMax() ) continue;

    w->var("Dst_M")->setVal(Dst_M);
    w->var("D0_M")->setVal(D0_M);
    w->var("D0_LTIME_ps")->setVal(D0_LTIME_ps);

    data->add( *observables );
    dataH->add( *observables );

  }
  tf->Close();

  w->import(*data);
  w->import(*dataH);

  w->writeToFile(outfname);

}

void PDF_GLWADS_DK_hh_Dmix::initRelations()
{
	RooArgSet *p = (RooArgSet*)parameters;
	delete theory; theory = new RooArgList("theory"); ///< the order of this list must match that of the COR matrix!
	RooRealVar& Rcab   = *((RooRealVar*)p->find("RBRdkdpi"));
	RooRealVar& rbk    = *((RooRealVar*)p->find("r_dk"));
	RooRealVar& dbk    = *((RooRealVar*)p->find("d_dk"));
	RooConstVar& kbk = RooConst(1);
	RooConstVar& rbp = RooConst(0); // assume rb(Dpi)=0 for DK only
	RooConstVar& dbp = RooConst(1); // then the phase and the
	RooConstVar& kbp = RooConst(1); // coherence factor values don't matter
	RooConstVar& kf = RooConst(1);
	RooRealVar& rf     = *((RooRealVar*)p->find("rD_kpi"));
	RooConstVar& rfGLW = RooConst(1);
	RooRealVar& df     = *((RooRealVar*)p->find("dD_kpi"));
	RooConstVar& dfGLW = RooConst(0);
	RooRealVar& g      = *((RooRealVar*)p->find("g"));
	RooRealVar& xD     = *((RooRealVar*)p->find("xD"));
	RooRealVar& yD     = *((RooRealVar*)p->find("yD"));
	RooRealVar& AcpDKK = *((RooRealVar*)p->find("AcpDKK"));
	RooRealVar& AcpDpp = *((RooRealVar*)p->find("AcpDpipi"));
	RooConstVar& AcpD = RooConst(0);
	theory->add(*(new RooGLWADSDmixRkpVar("rkp_kpi_th",      "RooGLWADSDmixRkpVar", Rcab, rbk, dbk, kbk, rbp, dbp, kbp, rf,    df,    kf, g, xD, yD, _Mxy, "fav")));
	theory->add(*(new RooGLWADSDmixRkpVar("rkp_pipi_th",     "RooGLWADSDmixRkpVar", Rcab, rbk, dbk, kbk, rbp, dbp, kbp, rfGLW, dfGLW, kf, g, xD, yD, _Mxy, "fav")));
	theory->add(*(new RooGLWADSDmixRkpVar("rkp_kk_th",       "RooGLWADSDmixRkpVar", Rcab, rbk, dbk, kbk, rbp, dbp, kbp, rfGLW, dfGLW, kf, g, xD, yD, _Mxy, "fav")));
	theory->add(*(new RooGLWADSDmixAcpVar("afav_dk_kpi_th",  "RooGLWADSDmixAcpVar", rbk,  dbk, kbk, rf,    df,    kf, g, xD, yD, _Mxy, AcpD)));
	theory->add(*(new RooGLWADSDmixAcpVar("acp_dk_kk_th",    "RooGLWADSDmixAcpVar", rbk,  dbk, kbk, rfGLW, dfGLW, kf, g, xD, yD, _Mxy, AcpDKK)));
	theory->add(*(new RooGLWADSDmixAcpVar("acp_dk_pipi_th",  "RooGLWADSDmixAcpVar", rbk,  dbk, kbk, rfGLW, dfGLW, kf, g, xD, yD, _Mxy, AcpDpp)));
	theory->add(*(new RooGLWADSDmixRpmVar("rp_dk_kpi_th",    "RooGLWADSDmixRpmVar", rbk,  dbk, kbk, rf,    df,    kf, g, xD, yD, _Mxy, "+")));
	theory->add(*(new RooGLWADSDmixRpmVar("rm_dk_kpi_th",    "RooGLWADSDmixRpmVar", rbk,  dbk, kbk, rf,    df,    kf, g, xD, yD, _Mxy, "-")));
}

void PDF_ADS_DKDpi_K3pi_Dmix_newVars::initRelations()
{
  RooArgSet *p = (RooArgSet*)parameters;
  delete theory; theory = new RooArgList("theory"); ///< the order of this list must match that of the COR matrix!
  RooRealVar& rbk = *((RooRealVar*)p->find("r_dk"));
  RooRealVar& dbk = *((RooRealVar*)p->find("d_dk"));
  RooConstVar& kbk = RooConst(1);
  RooRealVar& rbp = *((RooRealVar*)p->find("r_dpi"));
  RooRealVar& dbp = *((RooRealVar*)p->find("d_dpi"));
  RooConstVar& kbp = RooConst(1);
  RooRealVar& kf  = *((RooRealVar*)p->find("kD_k3pi"));
  RooRealVar& rf  = *((RooRealVar*)p->find("rD_k3pi"));
	RooConstVar& rfGLW  = RooConst(1);
  RooRealVar& df  = *((RooRealVar*)p->find("dD_k3pi"));
	RooConstVar& dfGLW  = RooConst(0);
  RooRealVar& g   = *((RooRealVar*)p->find("g"));
  RooRealVar& xD  = *((RooRealVar*)p->find("xD"));
  RooRealVar& yD  = *((RooRealVar*)p->find("yD"));
	RooRealVar& F4pi = *((RooRealVar*)p->find("F_pipipipi"));
  RooRealVar& Rcab = *((RooRealVar*)p->find("RBRdkdpi"));
  RooConstVar& AcpD = RooConst(0);
	RooConstVar& AcpDpppp = RooConst(0);

	theory->add(*(new RooGLWADSDmixAcpADSVar   ("aads_dk_k3pi_th",   "RooGLWADSDmixAcpADSVar", rbk, dbk, kbk, rf, df, kf, g, xD, yD, _Mxy, AcpD)));
	theory->add(*(new RooGLWADSDmixAcpADSVar   ("aads_dpi_k3pi_th",  "RooGLWADSDmixAcpADSVar", rbp, dbp, kbp, rf, df, kf, g, xD, yD, _Mxy, AcpD)));
	theory->add(*(new RooGLWADSDmixAcpVar      ("afav_dk_k3pi_th",   "RooGLWADSDmixAcpVar"   , rbk, dbk, kbk, rf, df, kf, g, xD, yD, _Mxy, AcpD)));
	theory->add(*(new RooGLWADSDmixRADSVar     ("rads_dk_k3pi_th",   "RooGLWADSDmixRADSVar"  , rbk, dbk, kbk, rf, df, kf, g, xD, yD, _Mxy, "")));
	theory->add(*(new RooGLWADSDmixRADSVar     ("rads_dpi_k3pi_th",  "RooGLWADSDmixRADSVar"  , rbp, dbp, kbp, rf, df, kf, g, xD, yD, _Mxy, "")));
	theory->add(*(new RooGLWADSDmixRkpVar      ("rkp_k3pi_th",       "RooGLWADSDmixRkpVar"   , Rcab, rbk, dbk, kbk, rbp, dbp, kbp, rf, df, kf, g, xD, yD, _Mxy, "fav")));

}

void PDF_GLWADS_DKDpi_KSKpi_DCPV_Dmix::initRelations()
{
	RooArgSet *p = (RooArgSet*)parameters;
	delete theory; theory = new RooArgList("theory"); ///< the order of this list must match that of the COR matrix!
	RooRealVar& rbk = *((RooRealVar*)p->find("r_dk"));
	RooRealVar& dbk = *((RooRealVar*)p->find("d_dk"));
	RooConstVar& kbk = RooConst(1);
	RooRealVar& rbp = *((RooRealVar*)p->find("r_dpi"));
	RooRealVar& dbp = *((RooRealVar*)p->find("d_dpi"));
	RooConstVar& kbp = RooConst(1);
	RooRealVar& kf  = *((RooRealVar*)p->find("kD_kskpi"));
	RooRealVar& rf  = *((RooRealVar*)p->find("rD_kskpi"));
	RooRealVar& df  = *((RooRealVar*)p->find("dD_kskpi"));
	RooRealVar& g   = *((RooRealVar*)p->find("g"));
	RooRealVar& xD  = *((RooRealVar*)p->find("xD"));
	RooRealVar& yD  = *((RooRealVar*)p->find("yD"));
	RooRealVar& Rcab = *((RooRealVar*)p->find("RBRdkdpi"));
	RooRealVar&  AcpDfav = *((RooRealVar*)p->find("AcpDzKstpKm"));
	RooRealVar&  AcpDsup = *((RooRealVar*)p->find("AcpDzKstmKp"));

	theory->add(*(new RooGLWADSDmixRADSVar(  "rfavsup_dpi_kskpi_th", "RooGLWADSDmixRADSVar (inverse)", rbp, dbp, kbp, rf, df, kf, g, xD, yD, _Mxy, "inverse")));
	theory->add(*(new RooGLWADSDmixRkpVar(   "rfav_dkdpi_kskpi_th",  "RooGLWADSDmixRkpVar (fav)", Rcab, rbk, dbk, kbk, rbp, dbp, kbp, rf, df, kf, g, xD, yD, _Mxy, "fav")));
	theory->add(*(new RooGLWADSDmixRkpVar(   "rsup_dkdpi_kskpi_th",  "RooGLWADSDmixRkpVar (sup)", Rcab, rbk, dbk, kbk, rbp, dbp, kbp, rf, df, kf, g, xD, yD, _Mxy, "sup")));
	theory->add(*(new RooGLWADSDmixAcpVar(   "afav_dk_kskpi_th",     "RooGLWADSDmixAcpVar",    rbk, dbk, kbk, rf, df, kf, g, xD, yD, _Mxy, AcpDfav)));
	theory->add(*(new RooGLWADSDmixAcpADSVar("asup_dk_kskpi_th",     "RooGLWADSDmixAcpADSVar", rbk, dbk, kbk, rf, df, kf, g, xD, yD, _Mxy, AcpDsup)));
	theory->add(*(new RooGLWADSDmixAcpVar(   "afav_dpi_kskpi_th",    "RooGLWADSDmixAcpVar",    rbp, dbp, kbp, rf, df, kf, g, xD, yD, _Mxy, AcpDfav)));
	theory->add(*(new RooGLWADSDmixAcpADSVar("asup_dpi_kskpi_th",    "RooGLWADSDmixAcpADSVar", rbp, dbp, kbp, rf, df, kf, g, xD, yD, _Mxy, AcpDsup)));
}

void PDF_GLW_DK_2var::initRelations()
{
  RooArgSet *p = (RooArgSet*)parameters;
  RooRealVar  &rbk = *((RooRealVar*)p->find("r_dk"));
  RooRealVar  &dbk = *((RooRealVar*)p->find("d_dk"));
  RooConstVar &kbk = RooConst(1.);
  RooRealVar  &g   = *((RooRealVar*)p->find("g"));
	theory = new RooArgList("theory"); ///< the order of this list must match that of the COR matrix!
  theory->add(*(new RooGLWAcpVar("acpp_dk_th", "Acp+ (DK)", rbk, dbk, kbk, g, "+")));
  theory->add(*(new RooGLWRcpVar("rcpp_dk_th", "Rcp+ (DK)", rbk, dbk, kbk, g, "+")));
}

void PDF_GLW_DK_pipipi0_2var::initRelations()
{
  RooArgSet *p = (RooArgSet*)parameters;
  RooRealVar    &rbk = *((RooRealVar*)p->find("r_dk"));
  RooRealVar    &dbk = *((RooRealVar*)p->find("d_dk"));
  RooFormulaVar &kbk = *(new RooFormulaVar("dilutionpp_2var","dilutionpp_2var","2*F_pipipi0 - 1", *p));
  RooRealVar    &g   = *((RooRealVar*)p->find("g"));
	theory = new RooArgList("theory"); ///< the order of this list must match that of the COR matrix!
  theory->add(*(new RooGLWAcpVar("aglw_dk_pipipi0_th", "Aglw+ (DK_pipipi0)", rbk, dbk, kbk, g, "+")));
  theory->add(*(new RooGLWRcpVar("rglw_dk_pipipi0_th", "Rglw+ (DK_pipipi0)", rbk, dbk, kbk, g, "+")));
}

void reweight_pt_1d(TH1F* weight_pt, TH1F* weight_pto,TH2F*weight_ptn,TH2F*weight_pt2o,TH2F* weight_ptnr, TH2F* weight_pt2,RooDataSet **dset, RooDataSet *dsetdestination, int numvar){

  if (!(*dset)) return;
///numerator and denominator
  TH1F *hnum = new TH1F("hnum","hnum",n_ptbins_forreweighting,ptbins_forreweighting);
  TH1F *hden = new TH1F("hden","hden",n_ptbins_forreweighting,ptbins_forreweighting);
  hnum->Sumw2();
  hden->Sumw2();

  const char* ptname=Form("roopt%d",numvar);
// RooAbsData->get*() Load a given row of data
 //getRealValue Get value of a RooAbsReal stored in set with given name. If none is found, value of defVal is returned.
  for (int i=0; i<(*dset)->numEntries(); i++){
    hden->Fill(fabs((*dset)->get(i)->getRealValue(ptname)),(*dset)->store()->weight(i));
  }
  for (int i=0; i<dsetdestination->numEntries(); i++){
    hnum->Fill(fabs(dsetdestination->get(i)->getRealValue(ptname)),dsetdestination->store()->weight(i));
  }
  
 //normalize to one 
  hnum->Scale(1.0/hnum->Integral());
  hden->Scale(1.0/hden->Integral());
  hnum->Divide(hden);
  TH1F *h = hnum;

  RooDataSet *newdset = new RooDataSet(**dset,Form("%s_ptrew",(*dset)->GetName()));
  newdset->reset();
  for (int i=0; i<(*dset)->numEntries(); i++){
    RooArgSet args = *((*dset)->get(i));
    float oldw = (*dset)->store()->weight(i);
    float pt = args.getRealValue(ptname);
    float neww = oldw*h->GetBinContent(h->FindBin(fabs(pt)));
    if(debug){
		weight_pt->Fill(neww);	
		weight_pto->Fill(oldw);
		weight_pt2o->Fill(h->FindBin(fabs(pt)),oldw);
		weight_ptn->Fill(h->FindBin(fabs(pt)),neww);	
		if(oldw!=0 && neww!=0)weight_ptnr->Fill(h->FindBin(fabs(pt)),oldw/neww);
		else {weight_ptnr->Fill(-10,1);}
		if(oldw!=0 && neww!=0)weight_pt2->Fill(pt,oldw/neww);
		else {weight_pt2->Fill(-10,1);}
     }
	 newdset->add(args,neww);
  }

  newdset->SetName((*dset)->GetName());
  newdset->SetTitle((*dset)->GetTitle());
  delete hnum; delete hden;

  RooDataSet *old_dset = *dset;
  *dset=newdset;
  std::cout << "Pt 1d rew: norm from " << old_dset->sumEntries() << " to " << newdset->sumEntries() << std::endl;
  delete old_dset;
};

//2d reweighting of rho and its sigma
void reweight_rhosigma(TH1F* weight_rho, TH1F* weight_rhoo,TH2F*weight_rhon,TH2F*weight_rho2o,TH2F* weight_rhonr, TH2F* weight_rho2,TH2F*weight_sigman,TH2F*weight_sigma2o,TH2F* weight_sigmanr, TH2F* weight_sigma2,RooDataSet **dset, RooDataSet *dsetdestination, bool deleteold){
  if (!(*dset)) return;
//  TH2F *hnum = new TH2F("hnum","hnum",n_rhobins_forreweighting,rhobins_forreweighting,n_sigmabins_forreweighting,sigmabins_forreweighting);
//  TH2F *hden = new TH2F("hden","hden",n_rhobins_forreweighting,rhobins_forreweighting,n_sigmabins_forreweighting,sigmabins_forreweighting);
  TH2F *hnum = new TH2F("hnum","hnum",100,0,100,20,0,20);
  TH2F *hden = new TH2F("hden","hden",100,0,100,20,0,20);
  hnum->Sumw2();
  hden->Sumw2();
  for (int i=0; i<(*dset)->numEntries(); i++){
    hden->Fill(fabs((*dset)->get(i)->getRealValue("roorho")),fabs((*dset)->get(i)->getRealValue("roosigma")),(*dset)->store()->weight(i));
  }
  for (int i=0; i<dsetdestination->numEntries(); i++){
    hnum->Fill(fabs(dsetdestination->get(i)->getRealValue("roorho")),fabs(dsetdestination->get(i)->getRealValue("roosigma")),dsetdestination->store()->weight(i));
  }

  hnum->Scale(1.0/hnum->Integral());
  hden->Scale(1.0/hden->Integral());
//data/MC
  hnum->Divide(hden);
  TH2F *h = hnum;

  RooDataSet *newdset = new RooDataSet(**dset,Form("%s_rhosigmarew",(*dset)->GetName()));
  newdset->reset();
  for (int i=0; i<(*dset)->numEntries(); i++){
    RooArgSet args = *((*dset)->get(i));
    float oldw = (*dset)->store()->weight(i);
    float rho = args.getRealValue("roorho");
    float sigma = args.getRealValue("roosigma");
    float neww = oldw*h->GetBinContent(h->FindBin(rho,sigma));
	if(debug){
		weight_rho->Fill(neww);
		weight_rhoo->Fill(oldw);
		weight_rho2o->Fill(h->GetXaxis()->FindBin(rho),oldw);	
		weight_rhon->Fill(h->GetXaxis()->FindBin(rho),neww);	
		if(oldw!=0)weight_rhonr->Fill(h->GetXaxis()->FindBin(rho),oldw/neww);
		else {weight_rhonr->Fill(-10,1);}//cout << "dipho weight old 0" << endl;}
		if(oldw!=0)weight_rho2->Fill(rho,oldw/neww);
		weight_sigma2o->Fill(h->GetYaxis()->FindBin(sigma),oldw);	
		weight_sigman->Fill(h->GetYaxis()->FindBin(sigma),neww);	
		if(oldw!=0)weight_sigmanr->Fill(h->GetYaxis()->FindBin(sigma),oldw/neww);
		else {weight_sigmanr->Fill(-10,1);}//cout << "dipho weight old 0" << endl;}
		if(oldw!=0)weight_sigma2->Fill(sigma,oldw/neww);
	}
		newdset->add(args,neww);
	  }
	  newdset->SetName((*dset)->GetName());
	  newdset->SetTitle((*dset)->GetTitle());
	  delete hnum; delete hden;
	  RooDataSet *old_dset = *dset;
	  *dset=newdset;
	  std::cout << "RhoSigma2D rew: norm from " << old_dset->sumEntries() << " to " << newdset->sumEntries() << std::endl;
	  if (deleteold) delete old_dset;
	};

void PDF_ADS_DK_Kpi_2var::initRelations()
{
  RooArgSet *p = (RooArgSet*)parameters;
  RooRealVar  &rbk = *((RooRealVar*)p->find("r_dk"));
  RooRealVar  &dbk = *((RooRealVar*)p->find("d_dk"));
  RooConstVar &kbk = RooConst(1.);
  RooRealVar  &rd  = *((RooRealVar*)p->find("rD_kpi"));
  RooRealVar  &dd  = *((RooRealVar*)p->find("dD_kpi"));
  RooConstVar &kd  = RooConst(1.);
  RooRealVar  &g   = *((RooRealVar*)p->find("g"));
	theory = new RooArgList("theory"); ///< the order of this list must match that of the COR matrix!
  theory->add(*(new RooADSAVar("aads_dk_th", "Aads (DK)", rbk, dbk, kbk, rd, dd, kd, g )));
  theory->add(*(new RooADSRVar("rads_dk_th", "Rads (DK)", rbk, dbk, kbk, rd, dd, kd, g )));
}

RooMultiVarGaussian *MultiVarGaussian(RooDataSet *dataset){

  RooArgSet *args = (RooArgSet *) dataset->get()->Clone();
  //   argSet.remove(*chain->GetNLLVar(), kFALSE, kTRUE);

  RooArgSet *mu = (RooArgSet *)args->Clone();
  TMatrixDSym *matrix = dataset->covarianceMatrix();
  RooMultiVarGaussian *g = new RooMultiVarGaussian("multi","",RooArgList(*args), RooArgList(*mu),*matrix);
  return g;
  
  

  return NULL;
}

//____________________________________
void DoHypothesisTest(RooWorkspace* wks){


  // Use a RooStats ProfileLikleihoodCalculator to do the hypothesis test.
  ModelConfig model;
  model.SetWorkspace(*wks);
  model.SetPdf("model");

  //plc.SetData("data");

  ProfileLikelihoodCalculator plc;
  plc.SetData( *(wks->data("data") ));

  // here we explicitly set the value of the parameters for the null.
  // We want no signal contribution, eg. mu = 0
  RooRealVar* mu = wks->var("mu");
//   RooArgSet* nullParams = new RooArgSet("nullParams");
//   nullParams->addClone(*mu);
  RooArgSet poi(*mu);
  RooArgSet * nullParams = (RooArgSet*) poi.snapshot();
  nullParams->setRealValue("mu",0);


  //plc.SetNullParameters(*nullParams);
  plc.SetModel(model);
  // NOTE: using snapshot will import nullparams
  // in the WS and merge with existing "mu"
  // model.SetSnapshot(*nullParams);

  //use instead setNuisanceParameters
  plc.SetNullParameters( *nullParams);



  // We get a HypoTestResult out of the calculator, and we can query it.
  HypoTestResult* htr = plc.GetHypoTest();
  cout << "-------------------------------------------------" << endl;
  cout << "The p-value for the null is " << htr->NullPValue() << endl;
  cout << "Corresponding to a signifcance of " << htr->Significance() << endl;
  cout << "-------------------------------------------------\n\n" << endl;


}