C++ RooFitResult::edm方法代码示例

///
/// Test PDF implementation.
/// Performs a fit to the minimum.
///
bool PDF_Abs::test()
{
	bool quiet = false;
	if(quiet) RooMsgService::instance().setGlobalKillBelow(ERROR);
	fixParameters(observables);
	floatParameters(parameters);
	setLimit(parameters, "free");
	RooFormulaVar ll("ll", "ll", "-2*log(@0)", RooArgSet(*pdf));
	RooMinuit m(ll);
	if(quiet) m.setPrintLevel(-2);
	m.setNoWarn();
	m.setLogFile("/dev/zero");
	m.setErrorLevel(1.0);
	m.setStrategy(2);
	// m.setProfile(1);
	m.migrad();
	RooFitResult *f = m.save();
	bool status = !(f->edm()<1 && f->status()==0);
	if(!quiet) f->Print("v");
	delete f;
	if(quiet) RooMsgService::instance().setGlobalKillBelow(INFO);
	if(!quiet) cout << "pdf->getVal() = " << pdf->getVal() << endl;
	return status;
}

//.........这里部分代码省略.........

		RooAbsPdf * modelLL = new RooAddPdf("modelLL","modelLL",RooArgSet(*corrPdfLL,*bkgLL),*fracLL);
		RooAbsPdf * modelDD = new RooAddPdf("modelDD","modelDD",RooArgSet(*corrPdfDD,*bkgDD),*fracDD);
		RooAbsPdf * modelLLB = new RooAddPdf("modelLLB","modelLLB",RooArgSet(*corrPdfLLB,*bkgLLB),*fracLL);
		RooAbsPdf * modelDDB = new RooAddPdf("modelDDB","modelDDB",RooArgSet(*corrPdfDDB,*bkgDDB),*fracDD);

		// CREATE COMBINED DATASET
		RooDataSet * combData = new RooDataSet(Form("combData_%i",i),"combined data",RooArgSet(*MM,*cosThetaL,*cosThetaB),Index(*samples),Import("DD",*dataDD),Import("LL",*dataLL));

		Str2VarMap params;
		params["fL"] = fL;
		params["afb"] = afb;	
		Str2VarMap paramsB;
		paramsB["afbB"] = afbB;

		// FIT COS LEPTON
		RooSimultaneous * combModel = new RooSimultaneous(Form("combModel_%i",i),"",*samples);
		combModel->addPdf(*modelLL,"LL");
		combModel->addPdf(*modelDD,"DD");

		RooFitResult * res = safeFit(combModel,combData,params,&isInAllowedArea);	
	
		// FIT COS HADRON
		RooSimultaneous * combModelB = new RooSimultaneous(Form("combModelB_%i",i),"",*samples);
		combModelB->addPdf(*modelLLB,"LL");
		combModelB->addPdf(*modelDDB,"DD");

		RooFitResult * resB = safeFit(combModelB,combData,paramsB,&isInAllowedAreaB);

		cout << endl << fixed << setprecision(6) << "AfbB = " << afbB->getVal() << " +/- " << afbB->getError() << endl;
		cout << "Afb = " << afb->getVal() << " +/- " << afb->getError() << endl;
		cout << "fL = " << fL->getVal() << " +/- " << fL->getError() << endl;
		cout << endl;
		cout << "lepton:  " << res->edm() << "   "  << res->covQual() << endl;
		cout << "baryon:  " << resB->edm() << "   "  << resB->covQual() << endl;
		cout << endl;

		TH1F * fLsys = new TH1F(Form("fLsys_%i",i),"fLsys",40,-1,1);
		TH1F * afbsys = new TH1F(Form("afbsys_%i",i),"afbsys",40,-1,1);
		TH1F * afbBsys = new TH1F(Form("afbBsys_%i",i),"afbBsys",40,-1,1);
		TH1F * fLsys_frac = new TH1F(Form("fLsys_frac%i",i),"fLsys",40,-1,1);
		TH1F * afbsys_frac = new TH1F(Form("afbsys_frac%i",i),"afbsys",40,-1,1);
		TH1F * afbBsys_frac = new TH1F(Form("afbBsys_frac%i",i),"afbBsys",40,-1,1);


		RooAbsPdf * mybkgDD_2 = NULL, * mybkgDDB_2 = NULL;
		buildBkgPdfs(q2min[i],q2max[i],"DD",CutsDef::DDcut,&mybkgDD_2,&mybkgDDB_2,"RooKeyPdf");

		//cout << nevts << endl;
		//TRandom3 r(0);

		for(int e = 0; e < nexp; e++)
		{
			histFile->cd();
			RooAbsPdf * toypdf = (RooAbsPdf *)modelDD->Clone();
			Analysis * toy = new Analysis("toy",cosThetaL,modelDD,nevts);
			RooAbsPdf * toypdfB = (RooAbsPdf *)modelDDB->Clone();
			Analysis * toyB = new Analysis("toyB",cosThetaB,modelDDB,nevts);
			
			afb->setVal(0);
			afbB->setVal(-0.37);
			fL->setVal(0.6);

			safeFit(toypdf,toy->GetDataSet("-recalc"),params,&isInAllowedArea);
			safeFit(toypdfB,toyB->GetDataSet("-recalc"),paramsB,&isInAllowedAreaB);
			double def_afb = afb->getVal();

//.........这里部分代码省略.........
	}
	RooDataSet *startPars = new RooDataSet("startParsForce", "startParsForce", *w->set(parsName));
	startPars->add(*w->set(parsName));

	// set up parameters and ranges
	RooArgList *varyPars = new RooArgList();
	TIterator* it = w->set(parsName)->createIterator();
	while ( RooRealVar* p = (RooRealVar*)it->Next() )
	{
		if ( p->isConstant() ) continue;
		if ( forceVariables=="" && ( false
					|| TString(p->GetName()).BeginsWith("d") ///< use these variables
					// || TString(p->GetName()).BeginsWith("r")
					|| TString(p->GetName()).BeginsWith("k")
					|| TString(p->GetName()) == "g"
					) && ! (
						TString(p->GetName()) == "rD_k3pi"  ///< don't use these
						|| TString(p->GetName()) == "rD_kpi"
						// || TString(p->GetName()) == "dD_kpi"
						|| TString(p->GetName()) == "d_dk"
						|| TString(p->GetName()) == "d_dsk"
						))
		{
			varyPars->add(*p);
		}
		else if ( forceVariables.Contains(TString(p->GetName())+",") )
		{
			varyPars->add(*p);
		}
	}
	delete it;
	int nPars = varyPars->getSize();
	if ( debug ) cout << "Utils::fitToMinForce() : nPars = " << nPars << " => " << pow(2.,nPars) << " fits" << endl;
	if ( debug ) cout << "Utils::fitToMinForce() : varying ";
	if ( debug ) varyPars->Print();

	//////////

	r = fitToMinBringBackAngles(w->pdf(pdfName), false, printlevel);

	//////////

	int nErrors = 0;

	// We define a binary mask where each bit corresponds
	// to parameter at max or at min.
	for ( int i=0; i<pow(2.,nPars); i++ )
	{
		if ( debug ) cout << "Utils::fitToMinForce() : fit " << i << "        \r" << flush;
		setParameters(w, parsName, startPars->get(0));

		for ( int ip=0; ip<nPars; ip++ )
		{
			RooRealVar *p = (RooRealVar*)varyPars->at(ip);
			float oldMin = p->getMin();
			float oldMax = p->getMax();
			setLimit(w, p->GetName(), "force");
			if ( i/(int)pow(2.,ip) % 2==0 ) { p->setVal(p->getMin()); }
			if ( i/(int)pow(2.,ip) % 2==1 ) { p->setVal(p->getMax()); }
			p->setRange(oldMin, oldMax);
		}

		// check if start parameters are sensible, skip if they're not
		double startParChi2 = getChi2(w->pdf(pdfName));
		if ( startParChi2>2000 ){
			nErrors += 1;
			continue;
		}

		// refit
		RooFitResult *r2 = fitToMinBringBackAngles(w->pdf(pdfName), false, printlevel);

		// In case the initial fit failed, accept the second one.
		// If both failed, still select the second one and hope the
		// next fit succeeds.
		if ( !(r->edm()<1 && r->covQual()==3) ){
			delete r;
			r = r2;
		}
		else if ( r2->edm()<1 && r2->covQual()==3 && r2->minNll()<r->minNll() ){
			// better minimum found!
			delete r;
			r = r2;
		}
		else{
			delete r2;
		}
	}

	if ( debug ) cout << endl;
	if ( debug ) cout << "Utils::fitToMinForce() : nErrors = " << nErrors << endl;

	RooMsgService::instance().setGlobalKillBelow(INFO);

	// (re)set to best parameters
	setParameters(w, parsName, r);

	delete startPars;
	return r;
}

//.........这里部分代码省略.........
    else fin = TFile::Open(fname.c_str());

    //-- normalisation in the MC --//
    float mc_norm_1ltop = 0;
    float mc_norm_tt2l = 0;
    float mc_norm_Wjets = 0;
    //float mc_norm_rare = 0;

    // C r e a t e   m o d e l   f o r  CR1_peak_lowM3b
    // -------------------------------------------------------------
    // Construct pdfs for 1ltop, tt2l, Wjets and rare
    TH1F* histo_1ltop = 0;
    TH1F* histo_tt2l = 0;
    TH1F* histo_Wjets = 0;
    
    RooHistPdf *pdf_1ltop  = GetRooHistPdf(fin,region,PROCESS_NAME_TT_1L,varname,&var,mc_norm_1ltop, histo_1ltop, setup.do_mcstat, DO_NORM, setup.Ndata*setup.rel_norm_1ltop);
    RooHistPdf *pdf_tt2l   = GetRooHistPdf(fin,region,PROCESS_NAME_TT_2L,varname,&var,mc_norm_tt2l, histo_tt2l, setup.do_mcstat, DO_NORM, setup.Ndata*setup.rel_norm_tt2l);
    RooHistPdf *pdf_Wjets  = GetRooHistPdf(fin,region,PROCESS_NAME_WJETS,varname,&var,mc_norm_Wjets, histo_Wjets, setup.do_mcstat, DO_NORM, setup.Ndata*setup.rel_norm_Wjets);
    //RooHistPdf *pdf_rare   = GetRooHistPdf(fin,region,PROCESS_NAME_RARE,varname,&var,mc_norm_rare, setup.do_mcstat);

    //cerr<<"TT_1L: "<<mc_norm_1ltop<<endl;
    //cerr<<"TT_2l: "<<mc_norm_tt2l<<endl;
    //cerr<<"WJets: "<<mc_norm_Wjets<<endl;
    //cerr<<"OTHER: "<<mc_norm_rare<<endl;

    // normalization factors (RooRealVar)
    float val_1ltop = mc_norm_1ltop;
    float val_Wjets = mc_norm_Wjets;
    if(setup.do_init_uncert)
    {
        val_1ltop = setup.init_1ltop*mc_norm_1ltop;
        val_Wjets = setup.init_Wjets*mc_norm_Wjets;
    }
    RooRealVar norm_1ltop("norm_1ltop","norm_1ltop",val_1ltop,0.25*mc_norm_1ltop,10.*mc_norm_1ltop);
    RooRealVar norm_Wjets("norm_Wjets","norm_Wjets",val_Wjets,0.25*mc_norm_Wjets,10.*mc_norm_Wjets);
    RooRealVar norm_tt2l("norm_tt2l","norm_tt2l",mc_norm_tt2l,0.25*mc_norm_tt2l,2*mc_norm_tt2l);
    //RooRealVar norm_rare("norm_rare","norm_rare",mc_norm_rare,0.25*mc_norm_rare,2*mc_norm_rare);
    // possibility to study a systematic on it
    if(setup.do_xs_tt2l_sys) mc_norm_tt2l*=setup.xs_sysfactor;
    //if(setup.do_xs_rare_sys) mc_norm_rare*=setup.xs_sysfactor;
    //RooConstVar norm_rare("norm_rare","norm_rare",mc_norm_rare);

    /*
    RooAddPdf model("model","model",
            RooArgList(*pdf_1ltop,*pdf_tt2l,*pdf_Wjets,*pdf_rare),
            RooArgList(norm_1ltop,norm_tt2l,norm_Wjets,norm_rare)) ;
    */
    RooAddPdf model("model","model",
            RooArgList(*pdf_1ltop,*pdf_tt2l,*pdf_Wjets),
            RooArgList(norm_1ltop,norm_tt2l,norm_Wjets)) ;


    //RooDataHist *data_CR1_peak_lowM3b = GetRooData(fin,region,varname,&var);
    RooDataHist *data_CR1_peak_lowM3b = GetRooData(histo_1ltop,histo_Wjets, histo_tt2l,&var);

    fin->Close();


    //--  Constraints on single top and rare --//
    float RelUncert = 0.2;
    // Construct another Gaussian constraint p.d.f on "rare" bkg
    //RooGaussian constr_rare("constr_rare","constr_rare",norm_rare,RooConst(mc_norm_rare),RooConst(RelUncert*mc_norm_rare)) ;
    // Construct another Gaussian constraint p.d.f on "tt2l" bkg
    RooGaussian constr_tt2l("constr_tt2l","constr_tt2l",norm_tt2l,RooConst(mc_norm_tt2l),RooConst(RelUncert*mc_norm_tt2l)) ;

    // P e r f o r m   t em p l a t e   f i t
    // ---------------------------------------------------

    //Minimizer(type,algo) -- Choose minimization package and algorithm to use. Default is MINUIT/MIGRAD through the RooMinimizer
    //                       interface, but rare can be specified (through RooMinimizer interface). Select OldMinuit to use
    //                       MINUIT through the old RooMinuit interface
    //
    //     Type         Algorithm
    //     ------       ---------
    //     OldMinuit    migrad, simplex, minimize (=migrad+simplex), migradimproved (=migrad+improve)
    //     Minuit       migrad, simplex, minimize (=migrad+simplex), migradimproved (=migrad+improve)
    //     Minuit2      migrad, simplex, minimize, scan
    //     GSLMultiMin  conjugatefr, conjugatepr, bfgs, bfgs2, steepestdescent
    //     GSLSimAn     -


    // ---  Perform simultaneous fit of model to data and model_ctl to data_ctl --//
    //RooFitResult* res = model.fitTo(*data_CR1_peak_lowM3b,Save());
    //RooFitResult* res = model.fitTo(*data_CR1_peak_lowM3b,ExternalConstraints(constr_rare),ExternalConstraints(constr_tt2l),PrintLevel(-1),Save(),
    RooFitResult* res = model.fitTo(*data_CR1_peak_lowM3b,ExternalConstraints(constr_tt2l),PrintLevel(-1),Save(),
            Minimizer(setup.type.c_str(),setup.algo.c_str()),Verbose(0));

    //--- Writing the results ---///
    FitResult fitRes;
    fitRes.Reset();
    fitRes.norm_1ltop = mc_norm_1ltop;
    fitRes.SF_1ltop = GetSF(res,"norm_1ltop");
    fitRes.SF_Wjets = GetSF(res,"norm_Wjets");
    fitRes.edm = res->edm();
    fitRes.correlation = res->correlationMatrix()[0][1];
    fitRes.conditions = conditions;

    return fitRes;

}

double final4_D0::doFit(bool usePixel, bool isMC)
{
  gROOT->SetBatch(kTRUE);
  gROOT->SetStyle("Plain");

  //setTDRStyle();

  RooRealVar x("","",1.7,2.05);
  x.SetTitle("M(K#pi) [GeV/c^{ 2}]");

  RooRealVar mean("mean", "mean", 1.86484,1.5, 2.2);
  //RooRealVar mean("mean", "mean", 1.865116);
  RooRealVar sigma("sigma", "sigma", 0.017, 0.0002, 0.02);
  //RooRealVar sigma("sigma", "sigma", 0.015332);
  RooGaussian gauss("gauss","gaussian PDF", x, mean, sigma);

  RooRealVar alpha("alpha", "alpha", -1.0, -10.0, 10.0);
  RooRealVar power("power", "power", 3.0, 0.0, 50.0);
  RooCBShape cball("cball", "crystal ball PDF", x, mean, sigma, alpha, power);

  RooRealVar dm0("dm0", "dm0", 0.13957);
  dm0.setConstant(kTRUE);  
  RooRealVar shape("shape","shape",0.,-100.,100.);
  RooRealVar dstp1("p1","p1",0.,-500.,500.);
  RooRealVar dstp2("p2","p2",0.,-500.,500.);

  shape.setRange(0.000001,10.0);//was 0.02
  shape.setVal(0.0017);
  dstp1.setVal(0.45);
  dstp2.setVal(13.0);

  RooDstD0BG bkg("bkg","bkg",x,dm0,shape,dstp1,dstp2);

  RooRealVar c0("c0","c0",10.0,-10.0,11.0);
  RooRealVar c1("c1","c1",10.0,-10.0,11.0);
  RooRealVar c2("c2","c2",10.0,-10.0,11.0);
  RooRealVar c3("c3","c3",10.0,-10.0,11.0);
  RooRealVar c4("c4","c4",10.0,-10.0,11.0);
  RooRealVar c5("c5","c5",10.0,-10.0,11.0);
  RooRealVar c6("c6","c6",10.0,-10.0,11.0);
  RooRealVar c7("c7","c7",10.0,-10.0,11.0);
  RooRealVar c8("c8","c8",10.0,-10.0,11.0);
  RooGenericPdf cutoff("cutoff","cutoff","(@0 > @1)*(@2*abs(@[email protected]) + @3*pow(abs(@[email protected]),2) + @4*pow(abs(@[email protected]),3) + @5*pow(abs(@[email protected]),4) + @6*pow(abs(@[email protected]),5) + @7*pow(abs(@[email protected]),6) + @8*pow(abs(@[email protected]),7))",RooArgSet(x,dm0,c0,c1,c2,c3,c4,c5,c6));

  RooRealVar poly1("poly1","poly1",0.,-5000.0,5000.0);
  RooRealVar poly2("poly2","poly2",1.0,-5000.0,5000.0);
  RooRealVar poly3("poly3","poly3",1.0,-5000.0,5000.0);
  RooRealVar poly4("poly4","poly4",1.0,-5000.0,5000.0);

  RooPolynomial polybkg("polybkg","polybkg",x,RooArgSet(poly1));


  RooRealVar cheby0("cheby0","cheby0",1.0,-500.0,500.0);
  RooRealVar cheby1("cheby1","cheby1",1.0,-500.0,500.0);
  RooRealVar cheby2("cheby2","cheby2",1.0,-500.0,500.0);
  RooRealVar cheby3("cheby3","cheby3",1.0,-500.0,500.0);

  RooChebychev chebybkg("chebybkg","chebybkg",x,RooArgSet(cheby0,cheby1,cheby2,cheby3));

  RooRealVar mean2("mean2", "mean2", 0.14548,0.144, 0.147);
  RooRealVar sigma2("sigma2", "sigma2", 0.00065, 0.0002, 0.005);
  RooGaussian gauss2("gauss2","gaussian PDF 2", x, mean2, sigma2);

  RooRealVar S("S", "Signal Yield", 1100, 0, 300000);
  //RooRealVar S("S", "Signal Yield", 0, 0, 300000);
  RooRealVar SS("SS", "Signal Yield #2", 100, 0, 100000);
  RooRealVar S2("S2", "Signal2 Yield (MC only)", 0, 0, 200);
  RooRealVar B("B", "Background Yield", 4000, 0, 30000000);
  //RooRealVar B("B", "Background Yield", 0, 0, 30000000);

  //RooAddPdf sum("sum", "gaussian plus threshold PDF",RooArgList(gauss, bkg), RooArgList(S, B));
  RooAddPdf sum("sum", "gaussian plus linear PDF", RooArgList(gauss, polybkg), RooArgList(S,B));
  //RooAddPdf sum("sum", "background PDF",RooArgList(polybkg), RooArgList(B));
  //RooAddPdf sum("sum", "background PDF",RooArgList(chebybkg), RooArgList(B));
  //RooAddPdf sum("sum", "background PDF",RooArgList(cutoff), RooArgList(B));
  // RooAddPdf sum("sum", "gaussians plus threshold PDF",RooArgList(gauss, gauss2, bkg), RooArgList(S, SS, B));
  //RooAddPdf sum("sum", "crystal ball plus threshold PDF",RooArgList(cball, bkg), RooArgList(S, B));
  RooAddPdf sumMC("sumMC","double gaussian",RooArgList(gauss, gauss2), RooArgList(S, S2));

  fstream file;

  char filename[50];
  double cut=5.5;
  sprintf(filename,"D0Mass.dat");
  
  RooDataSet* data = RooDataSet::read(filename,RooArgList(x));
  RooFitResult* fit = 0;
  if (isMC == 0)
  {
    fit = sum.fitTo(*data,RooFit::Extended(),PrintLevel(1),Save(true),RooFit::NumCPU(8),RooFit::Strategy(2));
    file << "cut: " << cut << "GeV" << endl;
    file << "status: " << fit->status() << endl;
    file << "covQual: " << fit->covQual() << endl;
    file << "edm: " << fit->edm() <<  endl;
    file << "Yield: " <<  S.getVal() << " " << S.getError() << endl;
    file << "Bkg: " << B.getVal() << " " << B.getError() << endl; 
    file << "sigma: " << sigma.getVal() <<  " " << sigma.getError() << endl;
    file << "mean: " << mean.getVal() << " " << mean.getError() << endl;
    file << "shape: " << shape.getVal() << " " << shape.getError() << endl;
    file << "dstp1: " << dstp1.getVal() << " " << dstp1.getError() << endl;
//.........这里部分代码省略.........

RooFitResult * safeFit(RooAbsPdf * pdf, RooDataSet * data, Str2VarMap p, ISVALIDF_PTR isValid, string opt = "", int nfree = -1, RooArgSet * cons = NULL, RooAbsReal * nll = NULL)
{
	RooFitResult * res = NULL;

	RooRealVar cosThetaL("cosThetaL","cosThetaL",0.,-1.,1.);
	RooRealVar cosThetaB("cosThetaB","cosThetaB",0.,-1.,1.);
	
	RooArgSet obs(cosThetaL,cosThetaB);
	
	//if(opt.find("-scan")==string::npos) res = pdf->fitTo(*data,PrintLevel(-1),Save(),Extended(true)); 
	if(p.size()==1 && p.find("afb") != p.end())     p["fL"]  = GetParam(pdf,"fL");
	else if(p.size()==1 && p.find("fL") != p.end()) p["afb"] = GetParam(pdf,"afb");
	RooArgSet * nuisances = NULL;
	/*
	bool afb_iscost = false, fL_iscost = false, afbB_iscost = false;
	if (p.find("afb") != p.end())  { afb_iscost  = ((RooRealVar*)p["afb"])->getAttribute("Constant");  ((RooRealVar*)p["afb"])->setConstant();  }
	if (p.find("fL") != p.end())   { fL_iscost   = ((RooRealVar*)p["fL"])->getAttribute("Constant");   ((RooRealVar*)p["fL"])->setConstant();   }
	if (p.find("afbB") != p.end()) { afbB_iscost = ((RooRealVar*)p["afbB"])->getAttribute("Constant"); ((RooRealVar*)p["afbB"])->setConstant(); }
	RooArgSet * nuisances = copyFreePars(pdf,obs);
	if (p.find("afb") != p.end())  ((RooRealVar*)p["afb"])->setConstant(afb_iscost);
	if (p.find("afbB") != p.end()) ((RooRealVar*)p["afbB"])->setConstant(afbB_iscost);
	if (p.find("fL") != p.end())   ((RooRealVar*)p["fL"])->setConstant(fL_iscost);
	*/
	int np = 20;
	if((!res || res->covQual()!=3 || res->edm() > 0.1) && opt.find("-noscan")==string::npos)
	{
		if(!nll) nll = pdf->createNLL(*data);

		vector < double > mins, maxs, r;

		Str2VarMap::iterator iter; int pp = 0;
		for (iter = p.begin(); iter != p.end(); iter++) 
		{
			RooRealVar * curp = (RooRealVar *)iter->second;
			maxs.push_back(curp->getMax());
			mins.push_back(curp->getMin());
			r.push_back((maxs.back() - mins.back())/(double)np);
			pp++;
		}
		
		findMin(pdf,data,nll,p,mins,maxs,np,isValid,nfree,opt+"-nofit",cons,nuisances);
		
		double prec = 1e6;
		while (prec > 0.001)
		{
			double maxr = 0;
			maxs.clear(); mins.clear(); pp=0;
			for (iter = p.begin(); iter != p.end(); iter++) 
			{
				RooRealVar * curp = (RooRealVar *)iter->second;
				if((curp->getVal() + r[pp]) < curp->getMax()) maxs.push_back(curp->getVal() + r[pp]);
				else maxs.push_back(curp->getMax());
				if((curp->getVal() - r[pp]) > curp->getMin()) mins.push_back(curp->getVal() - r[pp]);
				else mins.push_back(curp->getMin());
				r[pp] = (maxs.back() - mins.back())/(double)np;
				if(r[pp] > maxr) maxr = r[pp];
				pp++;
			}
		
			prec = maxr;
			res = findMin(pdf,data,nll,p,mins,maxs,np,isValid,nfree,opt,cons,nuisances);
		}
		
		//if(!mynll) delete nll;
	}

	return res;
}