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

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

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

void fitSignalShapeW(int massBin,int id, int channels,int categ, int sample, 
		     /* float lumi, bool doSfLepton, */double rangeLow, double rangeHigh,
		     double bwSigma,
		     double fitValues[9], double fitErrors[9], double covQual[1]){
 // ------ root settings ---------
  gROOT->Reset();  
  gROOT->SetStyle("Plain");
  gStyle->SetPadGridX(kFALSE);
  gStyle->SetPadGridY(kFALSE);
  //gStyle->SetOptStat("kKsSiourRmMen");
  gStyle->SetOptStat("iourme");
  //gStyle->SetOptStat("rme");
  //gStyle->SetOptStat("");
  gStyle->SetOptFit(11);
  gStyle->SetPadLeftMargin(0.14);
  gStyle->SetPadRightMargin(0.06);
  // ------------------------------ 

  ROOT::Math::MinimizerOptions::SetDefaultTolerance( 1.E-7);

  stringstream FileName;
  //Insert the file here
  if(sample==1) FileName <<"root://lxcms03//data3/Higgs/150915/ZH125/ZZ4lAnalysis.root" ;
  else if(sample==2) FileName << "root://lxcms03//data3/Higgs/150915/WplusH125/ZZ4lAnalysis.root";
  else if(sample==3) FileName << "root://lxcms03//data3/Higgs/150915/WminusH125/ZZ4lAnalysis.root";
  else if(sample==4) FileName << "root://lxcms03//data3/Higgs/150915/ttH125/ZZ4lAnalysis.root";
  else {
    cout << "Wrong sample." << endl;
    return;
  }
    

  cout << "Using " << FileName.str() << endl;
  
 
  TFile* ggFile = TFile::Open(FileName.str().c_str()); 

  TTree* ggTree = (TTree*) ggFile->Get("ZZTree/candTree");

  float m4l;
  
  Short_t z1flav, z2flav; 
  float weight;

  Short_t nExtraLeptons;   
  float ZZPt;
  Short_t nJets;
  Short_t nBTaggedJets;
  std::vector<float> * jetpt = 0;
  std::vector<float> * jeteta = 0;
  std::vector<float> * jetphi = 0;
  std::vector<float> * jetmass = 0;
  float jet30pt[10];
  float jet30eta[10];
  float jet30phi[10];
  float jet30mass[10];
  float Fisher;
  
  int  nentries = ggTree->GetEntries();
 
  //--- ggTree part
  ggTree->SetBranchAddress("ZZMass",&m4l);
  ggTree->SetBranchAddress("Z1Flav",&z1flav);
  ggTree->SetBranchAddress("Z2Flav",&z2flav);
  ggTree->SetBranchAddress("genHEPMCweight",&weight);
  ggTree->SetBranchAddress("nExtraLep",&nExtraLeptons);
  ggTree->SetBranchAddress("nCleanedJets",&nJets);
  ggTree->SetBranchAddress("nCleanedJetsPt30BTagged",&nBTaggedJets);
  ggTree->SetBranchAddress("DiJetFisher",&Fisher);
  
  ggTree->SetBranchAddress("JetPt",&jetpt);
  ggTree->SetBranchAddress("JetEta",&jeteta);
  ggTree->SetBranchAddress("JetPhi",&jetphi);
  ggTree->SetBranchAddress("JetMass",&jetmass);
  ggTree->SetBranchAddress("ZZPt",&ZZPt);

  //--- rooFit part
  double xMin,xMax,xInit;
  xInit = (double) massBin;
  xMin = rangeLow;
  xMax = rangeHigh ;
  cout << "Fit range: [" << xMin << " , " << xMax << "]. Init value = " << xInit << endl;
  
  TH1F *hmass = new TH1F("hmass","hmass",200,xMin,xMax);
  //---------  
  RooRealVar x("mass","m_{4l}",xInit,xMin,xMax,"GeV");
  RooRealVar w("myW","myW",1.0,0.,1000.);
  RooArgSet ntupleVarSet(x,w);
  RooDataSet dataset("mass4l","mass4l",ntupleVarSet,WeightVar("myW"));

  for(int k=0; k<nentries; k++){
    ggTree->GetEvent(k);

    int njet30 = 0;
    for (unsigned int ijet = 0; ijet < jetpt->size(); ijet++) { 
      if ( (*jetpt)[ijet] > 30. ) {
	jet30pt[njet30] = (*jetpt)[ijet];      
	jet30eta[njet30] = (*jeteta)[ijet];
	jet30phi[njet30] = (*jetphi)[ijet];
	jet30mass[njet30] = (*jetmass)[ijet];
//.........这里部分代码省略.........

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

void crossfeeds_nondiag(TString title, 
			TString bkgfile,
			TString epsfile,
			TString txtfile,
			Double_t alpha_,
			Double_t mass_,
			Double_t n_,
			Double_t sigma_
			)
{

  RooRealVar mbc("mbc", "m_{BC}", 1.83, 1.89, "GeV");
  RooRealVar ebeam("ebeam", "Ebeam", 0., 100., "GeV");
  RooRealVar chg("chg", "Charge", -2, 2);
  RooCategory passed("passed", "Event should be used for plot");

  passed.defineType("yes", 1);
  passed.defineType("no", 0);

  RooRealVar arg_cutoff ("arg_cutoff", "Argus cutoff", 1.8865, 1.885, 1.8875,"GeV"); 
  RooRealVar arg_slope ("arg_slope", "Argus slope", -13, -100, 40);

  RooRealVar mbc_float ("mbc_float", "Floating D mass", mass_, "GeV"); 
  RooRealVar sigma ("sigma", "CB width", sigma_, "GeV"); 
  RooRealVar alpha("alpha", "CB shape cutoff", alpha_);
  RooRealVar n("n", "CB tail parameter", n_);

  RooCBShape cb_float ("cb_float", "Floating Crystal Barrel", mbc, mbc_float, sigma, alpha, n); 
  RooArgusBG argus("argus", "Argus BG", mbc, arg_cutoff, arg_slope);

  RooRealVar yld("yield", "D yield", 0, -30, 100000); 
  RooRealVar bkg("bkg", "Background", 20, 0, 40000);

  // Build pdf
  RooAddPdf sumpdf_float("sumpdf_float", "Generic D sum pdf", RooArgList(cb_float, argus),
			   RooArgList(yld, bkg));
  
  RooDataSet* dset = RooDataSet::read(bkgfile, RooArgList(mbc, ebeam, passed), "", "");

  RooPlot* xframe  = mbc.frame();

  RooDataSet* dset2 = dset->reduce("passed==1");

  dset2->plotOn(xframe);
  
  // RooFitResult* rv = sumpdf_float.fitTo(*dset2, Extended(kTRUE), Save(kTRUE),
  // 					Hesse(kTRUE), Verbose(kTRUE));
  RooFitResult* rv = sumpdf_float.fitTo(*dset2, "ermh");

  sumpdf_float.paramOn(xframe, dset2);

  if ((yld.getVal() < 0) && (-yld.getVal()/bkg.getVal() > 0.5)){
    yld.setVal(0);
    bkg.setVal(1);
  }
  
  sumpdf_float.plotOn(xframe);
  sumpdf_float.plotOn(xframe, Components(RooArgSet(argus)),
                      LineColor(kRed), LineStyle(kDashed));

  TCanvas* c1 = new TCanvas("c1","Canvas", 2);
  
  xframe->SetTitleOffset(2.2, "Y");
  xframe->SetTitleOffset(1.1, "X");
  xframe->SetTitle(title);

  c1->SetLeftMargin(0.17);
  xframe->Draw();
  
  if ( rv && rv->covQual() != 3){
    // fit has failed
    TText *txt = new TText();
    txt->SetTextSize(.08);
    txt->SetTextAlign(22);
    txt->SetTextAngle(30);
    txt->DrawTextNDC(0.5, 0.5, "FAILED");
  }
  

  c1->Update();
  c1->Print(epsfile);
  c1->Clear();

  FILE* table = fopen(txtfile.Data(), "w+");
  fprintf(table, "Name\t|| Value\t|| Error\n");
  //  fprintf(table, "yldsigma\t| %.10f\t| \n", yld.getVal()/yld.getError());
  fprintf(table, "entries\t| %.10f\t| \n", dset->numEntries());
  fprintf(table, "yld\t| %.10f\t| %.10f\n",  yld.getVal(), yld.getError());
  //  fprintf(table, "ratio\t| %.10f\t| \n",  yld.getVal()/dset->numEntries());
  //  fprintf(table, "ratioerr\t| %.10f\t| \n",  yld.getError()/dset->numEntries());
  fclose(table);

  cout << "Saved output as: " << txtfile << endl;

  rv->Delete();
}

///
/// Perform the 1d Prob scan.
/// Saves chi2 values and the prob-Scan p-values in a root tree
/// For the datasets stuff, we do not yet have a MethodDatasetsProbScan class, so we do it all in
/// MethodDatasetsProbScan
/// \param nRun Part of the root tree file name to facilitate parallel production.
///
int MethodDatasetsProbScan::scan1d(bool fast, bool reverse)
{
	if (fast) return 0; // tmp

	if ( arg->debug ) cout << "MethodDatasetsProbScan::scan1d() : starting ... " << endl;

    // Set limit to all parameters.
    this->loadParameterLimits(); /// Default is "free", if not changed by cmd-line parameter


    // Define scan parameter and scan range.
    RooRealVar *parameterToScan = w->var(scanVar1);
    float parameterToScan_min = hCL->GetXaxis()->GetXmin();
    float parameterToScan_max = hCL->GetXaxis()->GetXmax();

		// do a free fit
		RooFitResult *result = this->loadAndFit(this->pdf); // fit on data
		assert(result);
    RooSlimFitResult *slimresult = new RooSlimFitResult(result,true);
		slimresult->setConfirmed(true);
		solutions.push_back(slimresult);
		double freeDataFitValue = w->var(scanVar1)->getVal();

    // Define outputfile
    system("mkdir -p root");
    TString probResName = Form("root/scan1dDatasetsProb_" + this->pdf->getName() + "_%ip" + "_" + scanVar1 + ".root", arg->npoints1d);
    TFile* outputFile = new TFile(probResName, "RECREATE");

    // Set up toy root tree
    this->probScanTree = new ToyTree(this->pdf, arg);
    this->probScanTree->init();
    this->probScanTree->nrun = -999; //\todo: why does this branch even exist in the output tree of the prob scan?

    // Save parameter values that were active at function
    // call. We'll reset them at the end to be transparent
    // to the outside.
    RooDataSet* parsFunctionCall = new RooDataSet("parsFunctionCall", "parsFunctionCall", *w->set(pdf->getParName()));
    parsFunctionCall->add(*w->set(pdf->getParName()));

    // start scan
    cout << "MethodDatasetsProbScan::scan1d_prob() : starting ... with " << nPoints1d << " scanpoints..." << endl;
    ProgressBar progressBar(arg, nPoints1d);
    for ( int i = 0; i < nPoints1d; i++ )
    {
        progressBar.progress();
        // scanpoint is calculated using min, max, which are the hCL x-Axis limits set in this->initScan()
        // this uses the "scan" range, as expected
        // don't add half the bin size. try to solve this within plotting method

        float scanpoint = parameterToScan_min + (parameterToScan_max - parameterToScan_min) * (double)i / ((double)nPoints1d - 1);
				if (arg->debug) cout << "DEBUG in MethodDatasetsProbScan::scan1d_prob() " << scanpoint << " " << parameterToScan_min << " " << parameterToScan_max << endl;

        this->probScanTree->scanpoint = scanpoint;

        if (arg->debug) cout << "DEBUG in MethodDatasetsProbScan::scan1d_prob() - scanpoint in step " << i << " : " << scanpoint << endl;

        // don't scan in unphysical region
        // by default this means checking against "free" range
        if ( scanpoint < parameterToScan->getMin() || scanpoint > parameterToScan->getMax() + 2e-13 ) {
            cout << "it seems we are scanning in an unphysical region: " << scanpoint << " < " << parameterToScan->getMin() << " or " << scanpoint << " > " << parameterToScan->getMax() + 2e-13 << endl;
            exit(EXIT_FAILURE);
        }

        // FIT TO REAL DATA WITH FIXED HYPOTHESIS(=SCANPOINT).
        // THIS GIVES THE NUMERATOR FOR THE PROFILE LIKELIHOOD AT THE GIVEN HYPOTHESIS
        // THE RESULTING NUISANCE PARAMETERS TOGETHER WITH THE GIVEN HYPOTHESIS ARE ALSO
        // USED WHEN SIMULATING THE TOY DATA FOR THE FELDMAN-COUSINS METHOD FOR THIS HYPOTHESIS(=SCANPOINT)
        // Here the scanvar has to be fixed -> this is done once per scanpoint
        // and provides the scanner with the DeltaChi2 for the data as reference
        // additionally the nuisances are set to the resulting fit values

        parameterToScan->setVal(scanpoint);
        parameterToScan->setConstant(true);

        RooFitResult *result = this->loadAndFit(this->pdf); // fit on data
        assert(result);

        if (arg->debug) {
            cout << "DEBUG in MethodDatasetsProbScan::scan1d_prob() - minNll data scan at scan point " << scanpoint << " : " << 2 * result->minNll() << ": "<< 2 * pdf->getMinNll() << endl;
        }
        this->probScanTree->statusScanData = result->status();

        // set chi2 of fixed fit: scan fit on data
        // CAVEAT: chi2min from fitresult gives incompatible results to chi2min from pdf
        // this->probScanTree->chi2min           = 2 * result->minNll();
        this->probScanTree->chi2min           = 2 * pdf->getMinNll();
        this->probScanTree->covQualScanData   = result->covQual();
        this->probScanTree->scanbest  = freeDataFitValue;

        // After doing the fit with the parameter of interest constrained to the scanpoint,
        // we are now saving the fit values of the nuisance parameters. These values will be
        // used to generate toys according to the PLUGIN method.
        this->probScanTree->storeParsScan(); // \todo : figure out which one of these is semantically the right one
//.........这里部分代码省略.........

   void ws_constrained_profile3D( const char* wsfile = "rootfiles/ws-data-unblind.root",
                                   const char* new_poi_name = "n_M234_H4_3b",
                                   int npoiPoints = 20,
                                   double poiMinVal = 0.,
                                   double poiMaxVal = 20.,
                                   double constraintWidth = 1.5,
                                   double ymax = 10.,
                                   int verbLevel=0 ) {


     gStyle->SetOptStat(0) ;

     //--- make output directory.

     char command[10000] ;
     sprintf( command, "basename %s", wsfile ) ;
     TString wsfilenopath = gSystem->GetFromPipe( command ) ;
     wsfilenopath.ReplaceAll(".root","") ;
     char outputdirstr[1000] ;
     sprintf( outputdirstr, "outputfiles/scans-%s", wsfilenopath.Data() ) ;
     TString outputdir( outputdirstr ) ;


     printf("\n\n Creating output directory: %s\n\n", outputdir.Data() ) ;
     sprintf(command, "mkdir -p %s", outputdir.Data() ) ;
     gSystem->Exec( command ) ;


     //--- Tell RooFit to shut up about anything less important than an ERROR.
      RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR) ;



       if ( verbLevel > 0 ) { printf("\n\n Verbose level : %d\n\n", verbLevel) ; }


       TFile* wstf = new TFile( wsfile ) ;

       RooWorkspace* ws = dynamic_cast<RooWorkspace*>( wstf->Get("ws") );

       if ( verbLevel > 0 ) { ws->Print() ; }






       RooDataSet* rds = (RooDataSet*) ws->obj( "ra2b_observed_rds" ) ;

       if ( verbLevel > 0 ) {
          printf("\n\n\n  ===== RooDataSet ====================\n\n") ;
          rds->Print() ;
          rds->printMultiline(cout, 1, kTRUE, "") ;
       }





       ModelConfig* modelConfig = (ModelConfig*) ws->obj( "SbModel" ) ;
       RooAbsPdf* likelihood = modelConfig->GetPdf() ;

       RooRealVar* rrv_mu_susy_all0lep = ws->var("mu_susy_all0lep") ;
       if ( rrv_mu_susy_all0lep == 0x0 ) {
          printf("\n\n\n *** can't find mu_susy_all0lep in workspace.  Quitting.\n\n\n") ;
          return ;
       }





       //-- do BG only.
       rrv_mu_susy_all0lep->setVal(0.) ;
       rrv_mu_susy_all0lep->setConstant( kTRUE ) ;










       //-- do a prefit.

       printf("\n\n\n ====== Pre fit with unmodified nll var.\n\n") ;

       RooFitResult* dataFitResultSusyFixed = likelihood->fitTo(*rds, Save(true),Hesse(false),Minos(false),Strategy(1),PrintLevel(verbLevel));
       int dataSusyFixedFitCovQual = dataFitResultSusyFixed->covQual() ;
       if ( dataSusyFixedFitCovQual < 2 ) { printf("\n\n\n *** Failed fit!  Cov qual %d.  Quitting.\n\n", dataSusyFixedFitCovQual ) ; return ; }
       double dataFitSusyFixedNll = dataFitResultSusyFixed->minNll() ;

       if ( verbLevel > 0 ) {
          dataFitResultSusyFixed->Print("v") ;
       }

       printf("\n\n Nll value, from fit result : %.3f\n\n", dataFitSusyFixedNll ) ;

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

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

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





       printf("\n\n\n === Model values for observables\n\n") ;

       printObservables() ;



      //--- save the actual values of the observables.

       saveObservables() ;











       //--- evaluate the test stat on the data: fit with susy floating.

       rrv_mu_susy_sig->setVal( poiVal ) ;
       rrv_mu_susy_sig->setConstant( kTRUE ) ;

       printf("\n\n\n ====== Fitting the data with susy fixed.\n\n") ;

       RooFitResult* dataFitResultSusyFixed = likelihood->fitTo(*rds, Save(true));
       int dataSusyFixedFitCovQual = dataFitResultSusyFixed->covQual() ;
       if ( dataSusyFixedFitCovQual != 3 ) { printf("\n\n\n *** Failed fit!  Cov qual %d.  Quitting.\n\n", dataSusyFixedFitCovQual ) ; return ; }
       double dataFitSusyFixedNll = dataFitResultSusyFixed->minNll() ;


       rrv_mu_susy_sig->setVal( 0.0 ) ;
       rrv_mu_susy_sig->setConstant( kFALSE ) ;

       printf("\n\n\n ====== Fitting the data with susy floating.\n\n") ;

       RooFitResult* dataFitResultSusyFloat = likelihood->fitTo(*rds, Save(true));
       int dataSusyFloatFitCovQual = dataFitResultSusyFloat->covQual() ;
       if ( dataSusyFloatFitCovQual != 3 ) { printf("\n\n\n *** Failed fit!  Cov qual %d.  Quitting.\n\n", dataSusyFloatFitCovQual ) ; return ; }
       double dataFitSusyFloatNll = dataFitResultSusyFloat->minNll() ;

       double dataTestStat = 2.*( dataFitSusyFixedNll - dataFitSusyFloatNll) ;

       printf("\n\n\n Data value of test stat : %8.2f\n", dataTestStat ) ;












       printf("\n\n\n === Nuisance parameters\n\n") ;

       {