C++ RooRealVar::setVal方法代码示例

void PEs(RooAbsPdf *iGen,RooAbsPdf *iFit,int iN,int iNEvents,RooRealVar &iVar,RooRealVar &iSig,RooRealVar &iMean,
	 RooRealVar &iScale,RooRealVar &iRes) { 
  double iM0 = iMean.getVal(); double iS0 = iSig.getVal();
  //iScale.setVal(iMeanScale); iRes.setVal(iSigScale);
  TRandom1 *lRand = new TRandom1(0xDEADBEEF);
  TNtuple * lDN= new TNtuple( "xxx","xxx","ntot:m_r:m:merr:sig_r:sig:sigerr");
  for(int i0=0;i0<iN;i0++){
    if(i0 % 10 == 0) cout << "+++++++++++++++++++++++++++ running ======> " << i0 << endl;
    int lN    = lRand->Poisson(iNEvents);
    RooDataSet * lSignal  = iGen->generate(iVar,lN);
    iMean.setVal(iM0); iSig.setVal(iS0);
    iFit->fitTo(*lSignal,Strategy(1));//,Save(kTRUE),PrintLevel(1));
    if(iMean.getError() < 0.05) iFit->fitTo(*lSignal,Strategy(2));
    Float_t values[]={
      (Float_t) lN,
      (Float_t) 90.78/iScale.getVal(),
      (Float_t) iMean.getVal(),
      (Float_t) iMean.getError(),
      (Float_t) iSig.getVal(),
      (Float_t) iSig.getVal(),
      (Float_t) iSig.getError()
    };
    lDN->Fill(values);
  }
  TFile *lF = new TFile("XXX.root","RECREATE");
  lDN->Write();
  lF->Close();
}

TGraph *graphLH(std::string nuisname, double err ){

	w->loadSnapshot("bestfitall"); // SetTo BestFit values as start

	// Get The parameter we want 
	RooRealVar *nuis =(RooRealVar*) w->var(nuisname.c_str());
	double bf = nuis->getVal();
	double nll_0=nll->getVal();


	TGraph *gr = new TGraph(2*npoints+1);
	for (int i=-1*npoints;i<=npoints;i++){
		nuis->setVal(bf+err*( ((float)i)*nsigma/npoints));
		double nll_v = nll->getVal();
		gr->SetPoint(i+npoints,nuis->getVal(),nll_v-nll_0);
	}

	gr->SetTitle("");
	gr->GetYaxis()->SetTitle("NLL - obs data");
	gr->GetYaxis()->SetTitleOffset(1.1);
	gr->GetXaxis()->SetTitleSize(0.05);
	gr->GetYaxis()->SetTitleSize(0.05);
	gr->GetXaxis()->SetTitle(nuisname.c_str());
	gr->SetLineColor(4);
	gr->SetLineWidth(2);
	gr->SetMarkerStyle(21);
	gr->SetMarkerSize(0.6);
	
	return gr;
	
}

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); 
}

RooStats::ModelConfig * Tprime::SetBModel( void ) {
    //
    // Define model config and parameter snapshot to describe the b model.
    // Import to workspace.
    //

    std::string legend = "[Tprime::SetBModel]: ";

    // full signal+background model
    //RooStats::ModelConfig * pSBModel = (RooStats::ModelConfig *)pWs->genobj("ModelConfig");

    // let's make the b model (bg-only) from the alt model (s+b) with xsec=0
    //RooStats::ModelConfig * pBModel =
    //  new RooStats::ModelConfig(*(RooStats::ModelConfig *)pWs->genobj("ModelConfig"));
    RooStats::ModelConfig * _sbModel = (RooStats::ModelConfig *)pWs->genobj("ModelConfig");
    RooStats::ModelConfig * pBModel = _sbModel->Clone("BModel");

    //pBModel->SetName("BModel");

    pBModel->SetWorkspace(*pWs);

    //pBModel->SetParametersOfInterest(RooArgSet());
    pWs->import(*pBModel);

    // set POI to the b model value and take snapshot
    RooRealVar * pPoi = (RooRealVar *)pBModel->GetParametersOfInterest()->first();
    pPoi->setVal(0.0);
    pBModel->SetSnapshot(*pPoi);

    pBModel->Print();

    return pBModel;
}

void testResolution() { 
  Prep();
  RooRealVar    lXVar  ("XVar","mass(GeV/c^{2})",100,60,150); lXVar.setBins(1000);
  RooRealVar l1Sigma("sigma1","sigma1",1.6 ,0.,15.);  //l1Sigma.setConstant(kTRUE);
  RooRealVar l2Sigma("sigma2","sigma2",1.6,0.,15.);   l2Sigma.setConstant(kTRUE);
  RooRealVar l3Sigma("sigma3","sigma3",2.90,0.,35.);  l3Sigma.setConstant(kTRUE);
  RooRealVar lN     ("n"     ,"n"     ,1.00,-15,15.); lN.setConstant(kTRUE);
  RooRealVar lExp   ("exp"   ,"exp"   ,-0.003,-15,15.); //lExp.setConstant(kTRUE);
  RooRealVar lR0Mean("xmean","xmean",0,-10,10);    lR0Mean.setConstant(kTRUE);
  RooRealVar lR1Mean("mean","mean",90.8,60,150);   //lR1Mean.setConstant(kTRUE);
  RooVoigtianShape     lGAdd("Add","Add",lXVar,lR1Mean,l1Sigma,l2Sigma,lN,l3Sigma,true);
  
  RooRealVar    lSPar  ("SPar","SPar", 1.,0., 2.);
  RooFormulaVar lXShift("uparshift","@0*@1",RooArgList(lXVar,lSPar));
  TH1F *lMass = getMass(0,-5,5,-1.5,1.5);
  RooDataHist *lMHist = new RooDataHist("M" ,"M" ,RooArgSet(lXVar),lMass);
  RooHistPdf  *lMPdf  = new RooHistPdf ("MH","MH",lXShift,lXVar,*lMHist,5); 
  RooRealVar lGSigma("gsigma","gsigma",1.6 ,0.,15.);  
  RooGaussian   lGaus1("gaus1","gaus1",lXVar,lR0Mean,lGSigma);
  RooFFTConvPdf  lConv("Conv","Conv",lXVar,*lMPdf,lGaus1);

  RooDataSet *lData = new RooDataSet("crap","crap",RooArgSet(lXVar)); 
  fillData(lData,lXVar,0,-5,5,-1.5,1.5);
  lConv.fitTo(*lData,Strategy(2));
  lGSigma.setVal(lGSigma.getVal());
  lSPar.setVal(lSPar.getVal()*1.01);
  //cout << "=====> Check " << l1Sigma.getVal() << " --- " << lR1Mean.getVal() << "----" << lSPar.getVal() << endl;
  PEs(&lConv,&lGAdd,2000,500,lXVar,l1Sigma,lR1Mean,lSPar,lGSigma);
  lData = 0;
  Plot(&lConv,&lGAdd,2000,50000,lXVar,l1Sigma,lR1Mean,lSPar,lGSigma,lData);
}

//
// set value and range for a variable in the workspace
//
void setValRange (RooWorkspace* workspace, const char* name, double val, double vmin, double vmax)
{
  RooRealVar* var = workspace->var(name);
  if ( var ) {
    if ( vmax>vmin )  var->setRange(vmin,vmax);
    var->setVal(val);
  }
}

void UL_significance_Hybrid(Model* model,int n_toys,int random_seed=0,double sig=1){
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  cout<<"Calculating significance with the Hybrid method"<<endl;
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  
  //set the random seed
  RooRandom::randomGenerator()->SetSeed(random_seed);
    
   //get the calculator
  HybridCalculatorOriginal myhc(*model->get_data(),*model->get_sb_likelihood(),*model->get_b_likelihood());
  
  //for numbercounting experiments
  myhc.PatchSetExtended(false);
  
  //set likelihood ratio as the test statistics
  myhc.SetTestStatistic(1);
  
  //define the systematics to be used
  if (model->get_nuisance_set()) {
    myhc.UseNuisance(true);
    myhc.SetNuisancePdf(*model->get_nuisance_prior_pdf());
    myhc.SetNuisanceParameters(*model->get_nuisance_set());
  } else {
    myhc.UseNuisance(false);                            
  }
  
  //define the number of toys to be done
  myhc.SetNumberOfToys(n_toys);

  //RooArgSet* poi= model->get_POI_set();
  //  poi->first()->Set(sig);
  RooRealVar* firstPOI = (RooRealVar*) model->get_POI_set()->first();
  firstPOI->setVal(sig);

  //get the Hypotestresult
  HybridResult* hcResult = myhc.GetHypoTest();
  double significance = hcResult->Significance();
  double CLS= hcResult->CLs();
  double CLB= hcResult->CLb();
  double CLsplusb= hcResult->CLsplusb();
  double CLSerror= hcResult->CLsError(); 
  



  cout<<"POI: "<<sig<<endl;
  cout <<"significance:" << significance<<endl;
  cout<<"CLs: "<<CLS<<endl;
  cout<<"CLb: "<<CLB<<endl;
  cout<<"CLsplusb: "<<CLsplusb<<endl;
  cout<<"CLserror: "<<CLSerror<<endl;

  HybridPlot* plot=hcResult->GetPlot("hcPlot","p Values Plot",100);
  TCanvas *c1=new TCanvas;
  plot->Draw();
  c1->SaveAs("hybrid_REsult");
  
}

//____________________________________
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");

}

///
/// Set all observables to 'truth' values computed from the
/// current parameters.
///
void PDF_Abs::setObservablesTruth()
{
	obsValSource = "truth";
	for ( int i=0; i<nObs; i++ )
	{
		RooRealVar* pObs = (RooRealVar*)((RooArgList*)observables)->at(i);
		pObs->setVal(((RooRealVar*)((RooArgList*)theory)->at(i))->getVal());
	}
}

double effRmsFromSigmaCB(double meanCB, double meanCB_err, double sigmaCB, double sigmaCB_err) {

  RooRealVar mass("zmass","m(e^{+}e^{-})",75,105,"GeV/c^{2}");

  float zwidth=2.4952;
  RooRealVar bwMean("m_{Z}","BW Mean", 91.1876, "GeV/c^{2}");
  RooRealVar bwWidth("#Gamma_{Z}", "BW Width", zwidth, "GeV/c^{2}");

  RooRealVar cbBias ("#Deltam_{CB}", "CB Bias", 0.00, -10, 10, "GeV/c^{2}"); // typical value
  RooRealVar cbSigma("#sigma_{CB}", "CB Width", 1.5, 0.8, 5.0, "GeV/c^{2}");
  RooRealVar cbCut  ("a_{CB}","CB Cut", 1.0, 1.0, 3.0); // typical value
  RooRealVar cbPower("n_{CB}","CB Order", 2.5, 0.1, 20.0); // typical value

  float mean_lo = meanCB - meanCB_err;
  float mean_hi = meanCB + meanCB_err;
  float sigma_lo = sigmaCB - sigmaCB_err;
  float sigma_hi = sigmaCB + sigmaCB_err;

  int nstep_mean = 10;
  int nstep_sigma = 10;
  
  float step_mean = 2*meanCB_err/float(nstep_mean);
  float step_sigma = 2*sigmaCB_err/float(nstep_sigma);

  double max_rms = 0.0;
  double min_rms = 999.;
  
  RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);

  int nsteps=0;
  for(float mean=mean_lo; mean<mean_hi; mean+=step_mean) {
    for(float sigma=sigma_lo; sigma<sigma_hi; sigma+=step_sigma) {

      cbBias.setVal(meanCB);
      cbSigma.setVal(sigma);

      RooBreitWigner bw("bw", "bw", mass, bwMean, bwWidth);
      RooCBShape  cball("cball", "Crystal Ball", mass, cbBias, cbSigma, cbCut, cbPower);
      RooFFTConvPdf BWxCB("BWxCB", "bw X crystal ball", mass, bw, cball);
      
      RooDataSet *gendata = BWxCB.generate(mass,1000);
      TH1D *genh = (TH1D*) gendata->createHistogram("genh",mass,RooFit::Binning(100,75,105));
      
      math::HistogramTools ht(genh);
      double rms = ht.effectiveSigma();
      rms = sqrt(pow(rms,2)-pow(zwidth,2));
      if(rms > max_rms) max_rms=rms;
      if(rms < min_rms) min_rms=rms;
      delete genh;
      nsteps++;
    }
  }
  return (max_rms-min_rms)/sqrt(12)/sqrt(nsteps);
}

void FitterUtils::initiateParams(int nGenSignalZeroGamma, int nGenSignalOneGamma, int nGenSignalTwoGamma, RooRealVar const& expoConstGen, RooRealVar& nSignal, RooRealVar& nPartReco, 
      RooRealVar& nComb, RooRealVar& fracZero, RooRealVar& fracOne, RooRealVar& expoConst, RooRealVar&  nJpsiLeak, bool constPartReco, RooRealVar const& fracPartRecoSigma)
{
   TRandom rand;
   rand.SetSeed();

   int nGenSignal = nGenSignalZeroGamma + nGenSignalOneGamma + nGenSignalTwoGamma;

   double nGenSignal2;
   double nGenPartReco2;
   if(!constPartReco)
   {
      nGenSignal2 = rand.Uniform(nGenSignal-5*sqrt(nGenSignal), nGenSignal+5*sqrt(nGenSignal));
      nGenPartReco2 = rand.Uniform(nGenPartReco-5*sqrt(nGenPartReco), nGenPartReco+5*sqrt(nGenPartReco));
   }
   if(constPartReco)
   { 
      double nGenSigPartReco( nGenSignal+nGenPartReco );
      double nGenSigPartReco2( rand.Uniform( nGenSigPartReco-5*sqrt(nGenSigPartReco), nGenSigPartReco+5*sqrt(nGenSigPartReco) ) );
      double fracPartReco1( nGenPartReco/(1.*nGenSignal));
      double fracPartReco2( rand.Uniform(fracPartReco1-5*fracPartRecoSigma.getVal(), fracPartReco1+5*fracPartRecoSigma.getVal()) ); 

      nGenPartReco2 = fracPartReco2*nGenSigPartReco2 / (1+fracPartReco2); 
      nGenSignal2 = nGenSigPartReco2 / (1+fracPartReco2); 
   }
   double nGenComb2 = rand.Uniform(nGenComb-5*sqrt(nGenComb), nGenComb+5*sqrt(nGenComb));
   double nGenJpsiLeak2 = rand.Uniform(nGenJpsiLeak-5*sqrt(nGenJpsiLeak), nGenJpsiLeak+5*sqrt(nGenJpsiLeak));


   nSignal.setVal(nGenSignal2);
   nSignal.setRange(TMath::Max(0.,nGenSignal2-10.*sqrt(nGenSignal)) , nGenSignal2+10*sqrt(nGenSignal));

   nPartReco.setVal(nGenPartReco2);
   nPartReco.setRange(TMath::Max(0.,nGenPartReco2-10.*sqrt(nGenPartReco)), nGenPartReco2+10*sqrt(nGenPartReco));


   nComb.setVal(nGenComb2);
   nComb.setRange(TMath::Max(0.,nGenComb2-10.*sqrt(nGenComb)), nGenComb2+10*sqrt(nGenComb));

   nJpsiLeak.setVal(nGenJpsiLeak2);
   nJpsiLeak.setRange(TMath::Max(0., nGenJpsiLeak2-10*sqrt(nGenJpsiLeak)), nGenJpsiLeak2+10*sqrt(nGenJpsiLeak));

   double fracGenZero(nGenSignalZeroGamma/(1.*nGenSignal));
   double fracGenOne(nGenSignalOneGamma/(1.*nGenSignal));

   fracZero.setVal(rand.Gaus(fracGenZero, sqrt(nGenSignalZeroGamma)/(1.*nGenSignal))) ;
   fracZero.setRange(0., 1.);
   fracOne.setVal(rand.Gaus(fracGenOne, sqrt(nGenSignalOneGamma)/(1.*nGenSignal))) ;
   fracOne.setRange(0., 1.);

   expoConst.setVal(rand.Uniform( expoConstGen.getVal() - 5*expoConstGen.getError(), expoConstGen.getVal() + 5*expoConstGen.getError() ) );
   expoConst.setRange( expoConstGen.getVal() - 10*expoConstGen.getError(), expoConstGen.getVal() + 10*expoConstGen.getError() );
}

void fillData(RooDataSet *iData,RooRealVar &lXVar,int iCharge,float iPhiMin,float iPhiMax,float iEtaMin,float iEtaMax) { 
  TFile *lFile = new TFile("ZTP.root");
  TTree *lTree = (TTree*) lFile->FindObjectAny("WNtupleIdEffNT");
  int   lCharge = 0; lTree->SetBranchAddress("charge",&lCharge);
  float lEta    = 0; lTree->SetBranchAddress("eta"   ,&lEta);
  float lPhi    = 0; lTree->SetBranchAddress("phi"   ,&lPhi);
  float lMt     = 0; lTree->SetBranchAddress("mt"    ,&lMt);
  float lPt     = 0; lTree->SetBranchAddress("pt"    ,&lPt);
  float lOPt    = 0; lTree->SetBranchAddress("jetpt" ,&lOPt);
  float lOEta   = 0; lTree->SetBranchAddress("jeteta",&lOEta);
  float lOPhi   = 0; lTree->SetBranchAddress("jetphi",&lOPhi);
  for(int i0 = 0; i0 < lTree->GetEntries();i0++) { 
    lTree->GetEntry(i0);
    if(lMt < 60)                   continue;
    if(lCharge > 0 && iCharge < 0) continue;
    if(lCharge < 0 && iCharge > 0) continue;
    if(lPhi < iPhiMin || lPhi > iPhiMax) continue;
    if(lEta < iEtaMin || lEta > iEtaMax) continue;
    lXVar.setVal(fabs(lMt));
    if(lCharge > 0) lXVar.setVal(correct(lMt,lPhi,lOPhi,lEta,lOEta));
    if(lCharge < 0) lXVar.setVal(correct(lMt,lOPhi,lPhi,lOEta,lEta));
    iData->add(RooArgSet(lXVar));
  }
}

///
/// Set each parameter in workspace to the values found
/// in the fit result
///
void Utils::setParameters(RooWorkspace* w, RooFitResult* values){
	RooArgList list = values->floatParsFinal();
	list.add(values->constPars());
	TIterator* it = list.createIterator();
	while(RooRealVar* p = (RooRealVar*)it->Next()){
		RooRealVar* var = dynamic_cast<RooRealVar*>(w->allVars().find(p->GetName()));
		if(!(var)){
			std::cout <<  "WARNING in Utils::setParameters(RooWorkspace,RooFitResult) -- no Var found with name "
				<< p->GetName() << " in Workspace!" << endl;
		}
		else{
			var->setVal(p->getVal());
		}
	}
	return;
};

///
/// Set all observables to 'toy' values drawn from the
/// PDF using the current parameter values. A certain number
/// of toys is pregenerated to speed up when doing mulitple toy fits.
///
void PDF_Abs::setObservablesToy()
{
	obsValSource = "toy";
	if( !pdf ){ cout<< "PDF_Abs::setObservables(): ERROR: pdf not initialized."<<endl; exit(1); }
	if ( toyObservables==0 || iToyObs==nToyObs )
	{
		RooRandom::randomGenerator()->SetSeed(0);
		if ( iToyObs==nToyObs ) delete toyObservables;
		toyObservables = pdf->generate(*(RooArgSet*)observables, nToyObs);
		iToyObs=0;
	}
	for ( int i=0; i<nObs; i++ )
	{
		RooRealVar* pObs = (RooRealVar*)((RooArgList*)observables)->at(i);
		pObs->setVal(((RooRealVar*)toyObservables->get(iToyObs)->find(pObs->GetName()))->getVal());
	}
	iToyObs+=1;
}

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;
}