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

RooFitResult *breakDownFit(RooSimultaneous *m, RooAbsData *d, RooRealVar *mass, bool precondition = false){
	if(precondition){
		 const char *catsName = m->indexCat().GetName();
		 TIterator *it = m->indexCat().typeIterator();
		 while(RooCatType* ci = dynamic_cast<RooCatType*>(it->Next())) {
			 const Text_t *catLabel = ci->GetName();
			 RooAbsPdf *pdf = m->getPdf(Form("%s",catLabel));
			 RooAbsData *reduced = d->reduce(SelectVars(*mass),Cut(Form("%s==%s::%s",catsName, catsName, catLabel)));
			 RooFitResult *r = pdf->fitTo(*reduced,PrintLevel(-1),Save(),
					 Minimizer("Minuit2","migrad"),Strategy(0),Hesse(false),Minos(false),Optimize(false)
					 );
			 cout << catsName << " " << catLabel << " M2migrad0 " << r->status() << endl;
			 if(r->status()!=0){
				 RooFitResult *r = pdf->fitTo(*reduced, PrintLevel(-1), Save());
				 cout << catsName << " " << catLabel << " Mmigrad1 " << r->status() << endl;
			 }
		 }
	}

	RooFitResult *r = m->fitTo(*d, Save(), PrintLevel(-1),
			Strategy(0));
	cout << "Global fit Mmigrad0 " << r->status() << endl;
	if(r->status()!=0){
	 RooFitResult *r = m->fitTo(*d, PrintLevel(-1), Save(),
			 Minimizer("Minuit","minimize"),Strategy(2));
	 cout << "Global fit Mminimize2 " << r->status() << endl;
	 return r;
	}

	return r;
}

void runFit(RooAbsPdf *pdf, RooDataSet *data, double *NLL, int *stat_t, int MaxTries, int mhLow, int mhHigh){

	int ntries=0;
	int stat=1;
	double minnll=10e8;
	while (stat!=0){
	  if (ntries>=MaxTries) break;
	  RooFitResult *fitTest = pdf->fitTo(*data,RooFit::Save(1),Range(mhLow,mhHigh));
	  //RooFitResult *fitTest = pdf->fitTo(*data,RooFit::Save(1),Range(85,110));
	  //RooFitResult *fitTest = pdf->fitTo(*data,RooFit::Save(1),SumW2Error(kTRUE)
          stat = fitTest->status();
	  minnll = fitTest->minNll();
	  ntries++; 
	}
	*stat_t = stat;
	*NLL = minnll;
}

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

//.........这里部分代码省略.........
   const RooArgSet * poiSet = sbModel->GetParametersOfInterest();
   RooRealVar *poi = (RooRealVar*)poiSet->first();
  
   std::cout << "StandardHypoTestInvDemo : POI initial value:   " << poi->GetName() << " = " << poi->getVal()   << std::endl;  
  
   // fit the data first (need to use constraint )
   TStopwatch tw; 

   bool doFit = initialFit;
   if (testStatType == 0 && initialFit == -1) doFit = false;  // case of LEP test statistic
   if (type == 3  && initialFit == -1) doFit = false;         // case of Asymptoticcalculator with nominal Asimov
   double poihat = 0;

   if (minimizerType.size()==0) minimizerType = ROOT::Math::MinimizerOptions::DefaultMinimizerType();
   else 
      ROOT::Math::MinimizerOptions::SetDefaultMinimizer(minimizerType.c_str());
    
   Info("StandardHypoTestInvDemo","Using %s as minimizer for computing the test statistic",
        ROOT::Math::MinimizerOptions::DefaultMinimizerType().c_str() );
   
   if (doFit)  { 

      // do the fit : By doing a fit the POI snapshot (for S+B)  is set to the fit value
      // and the nuisance parameters nominal values will be set to the fit value. 
      // This is relevant when using LEP test statistics

      Info( "StandardHypoTestInvDemo"," Doing a first fit to the observed data ");
      RooArgSet constrainParams;
      if (sbModel->GetNuisanceParameters() ) constrainParams.add(*sbModel->GetNuisanceParameters());
      RooStats::RemoveConstantParameters(&constrainParams);
      tw.Start(); 
      RooFitResult * fitres = sbModel->GetPdf()->fitTo(*data,InitialHesse(false), Hesse(false),
                                                       Minimizer(minimizerType.c_str(),"Migrad"), Strategy(0), PrintLevel(mPrintLevel), Constrain(constrainParams), Save(true) );
      if (fitres->status() != 0) { 
         Warning("StandardHypoTestInvDemo","Fit to the model failed - try with strategy 1 and perform first an Hesse computation");
         fitres = sbModel->GetPdf()->fitTo(*data,InitialHesse(true), Hesse(false),Minimizer(minimizerType.c_str(),"Migrad"), Strategy(1), PrintLevel(mPrintLevel+1), Constrain(constrainParams), Save(true) );
      }
      if (fitres->status() != 0) 
         Warning("StandardHypoTestInvDemo"," Fit still failed - continue anyway.....");
  
  
      poihat  = poi->getVal();
      std::cout << "StandardHypoTestInvDemo - Best Fit value : " << poi->GetName() << " = "  
                << poihat << " +/- " << poi->getError() << std::endl;
      std::cout << "Time for fitting : "; tw.Print(); 
  
      //save best fit value in the poi snapshot 
      sbModel->SetSnapshot(*sbModel->GetParametersOfInterest());
      std::cout << "StandardHypoTestInvo: snapshot of S+B Model " << sbModel->GetName() 
                << " is set to the best fit value" << std::endl;
  
   }

   // print a message in case of LEP test statistics because it affects result by doing or not doing a fit 
   if (testStatType == 0) {
      if (!doFit) 
         Info("StandardHypoTestInvDemo","Using LEP test statistic - an initial fit is not done and the TS will use the nuisances at the model value");
      else 
         Info("StandardHypoTestInvDemo","Using LEP test statistic - an initial fit has been done and the TS will use the nuisances at the best fit value");
   }


   // build test statistics and hypotest calculators for running the inverter 
  
   SimpleLikelihoodRatioTestStat slrts(*sbModel->GetPdf(),*bModel->GetPdf());

void LL(){

  //y0 = 0.000135096401209 sigma_y0 = 0.000103896581837 x0 = 0.000446013873443 sigma_x0 =1.81384394011e-06
  //0.014108652249 0.0168368471049 0.0219755396247 0.000120423865262 1.5575931164 1.55759310722 3.41637854038
  //0.072569437325 0.084063541977 0.0376693978906 0.000284216132439 0.51908074913 0.519080758095 1.12037749267
 // double d = 0.014108652249;
 //  double sd = 0.0168368471049;
 //  double mc = 0.0219755396247;
 //  double smc = 0.000120423865262;
 //  double r0 = d/mc;

  double d = 0.072569437325;
  double sd =  0.084063541977;
  double mc =  0.0376693978906;
  double smc =  0.00028421613243;
  double r0 = d/mc;

  RooRealVar x("x","x",mc*0.9,mc*1.1);
  RooRealVar x0("x0","x0",mc);
  RooRealVar sx("sx","sx",smc);

  RooRealVar r("r","r",r0,0.,5.);
  RooRealVar y0("y0","y0",d); 
  RooRealVar sy("sy","sy",sd); 
  
  RooProduct rx("rx","rx",RooArgList(r,x));

  RooGaussian g1("g1","g1",x,x0,sx);
  RooGaussian g2("g2","g2",rx,y0,sy);

  RooProdPdf LL("LL","LL",g1,g2);

  RooArgSet obs(x0,y0); //observables
  RooArgSet poi(r); //parameters of interest
  RooDataSet data("data", "data", obs);
  data.add(obs); //actually add the data


  RooFitResult* res = LL.fitTo(data,RooFit::Minos(poi),RooFit::Save(),RooFit::Hesse(false));
  if(res->status()==0) {
    r.Print();
    x.Print();
    cout << r.getErrorLo() << " " << r.getErrorHi() << endl;
  } else {
    cout << "Likelihood maximization failed" << endl;
  }
  
  RooAbsReal* nll = LL.createNLL(data); 
  RooPlot* frame = r.frame();
  RooAbsReal* pll = nll->createProfile(poi);
  pll->plotOn(frame);//,RooFit::LineColor(ROOT::kRed));
  frame->Draw();

  r.setVal(0.);
  cout << pll->getVal() << endl; 

  return;
    
    


}

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

//.........这里部分代码省略.........
	h_ks_c.Fill(ks_c);

	delete h1_c;
      }


      //while ( chi2_l>0.5 || std::isnan(chi2_l) ){
      while ( ks_l<0.1 ){

	MultiGaus(parMeans_l, covMatrix_l, genPars_l);
	//genPars_l.Print();

	l1.setVal(genPars_l[0]);
	l2.setVal(genPars_l[1]);
	l3.setVal(genPars_l[2]);
	l4.setVal(genPars_l[3]);
	l5.setVal(genPars_l[4]);
	l6.setVal(genPars_l[5]);

	TH1* h1_l = Pdf_l->createHistogram("h1_l",x); 
	h1_l->SetEntries((int) f_l_0*generated_events);
	h1_l->Scale(f_l_0*generated_events/h1_l->Integral());
	//h1_l->Draw("SAME");
	chi2_l = h1_l->Chi2Test(h0_l,"WW CHI2/NDF");
	h_chi2_l.Fill(chi2_l);
	ks_l = h1_l->KolmogorovTest(h0_l);
	h_ks_l.Fill(ks_l);

	delete h1_l;
      }

      RooFitResult * fitRes = model.fitTo(*dataHist,Save(),SumW2Error(kFALSE),Minimizer("Minuit2"),PrintLevel(1),Verbose(0));

      if (fitRes->status() != 0 ) continue;
      
      frame = x.frame();
      dataHist->plotOn(frame);
      model.plotOn(frame);
       
      //chi2 = frame->chiSquare("model_Norm[x]","h_dataHist",2);
      chi2 = frame->chiSquare(2);
      prob = TMath::Prob(chi2,2);
      
      h_chi2.Fill(chi2);
      h_prob.Fill(prob);
      

      //if  ( prob < 0.6 ) continue;
      if ( chi2 > 5. ) continue;
      if ( f_c.getVal()<0.01 ||
      	   f_l.getVal()<0.01 ||
	   (1.-f_c.getVal()-f_l.getVal())<0.01 ) continue;
      
      if ( f_c.getVal()>0.99 ||
      	   f_l.getVal()>0.99 ||
      	   (1.-f_c.getVal()-f_l.getVal())>0.99 ) continue;

      
      Pdf_b->plotOn(b_frame);
      Pdf_c->plotOn(c_frame);
      Pdf_l->plotOn(l_frame);
      
      
      h_b1.Fill(genPars_b[0]);
      h_b2.Fill(genPars_b[1]);
      h_b3.Fill(genPars_b[2]);

void fit(Int_t i, Double_t va=-0.43, Double_t vb=0.){
    sprintf(namestr,"%d_%.2f_%.2f",i,va,vb);

    Double_t g1, g2;
    g1 = 0.1*r.Rndm() - 0.05;
    g2 =     r.Rndm() - 0.5;

    //a->setVal(va);
    //b->setVal(vb);

    //G001->setVal(g1);
    //G002->setVal(g2);

    //RooRealVar * G000 = new RooRealVar("G000", "G000",  0.5);
    //RooRealVar * G001 = new RooRealVar("G001", "G001",  g1, -5., 5.);
    //RooRealVar * G002 = new RooRealVar("G002", "G002",  g2, -5., 5.);
    //RooRealVar * G003 = new RooRealVar("G003", "G003", 0.0);//, -1., 1.);
    //RooRealVar * G004 = new RooRealVar("G004", "G004", 0.0);//, -1., 1.);
    //RooRealVar * a    = new RooRealVar("a",    "a",    va);
    //RooRealVar * b    = new RooRealVar("b",    "b",    vb);
    //RooRealVar * n    = new RooRealVar("n",    "n",    1000., -100., 5000.);
    RooRealVar G000("G000", "G000",  0.5);
    RooRealVar G001("G001", "G001",  g1, -5., 5.);
    RooRealVar G002("G002", "G002",  g2, -5., 5.);
    RooRealVar G003("G003", "G003", 0.0);//, -1., 1.);
    RooRealVar G004("G004", "G004", 0.0);//, -1., 1.);
    RooRealVar    a("a",    "a",    va);
    RooRealVar    b("b",    "b",    vb);
    RooRealVar    n("n",    "n",    1000., -100., 5000.);

    RooB2Kll pdf("pdf", "pdf", *cosTheta, G000, G001, G002, G003, G004, a, b, n);

    RooFitResult * fitresult = pdf.fitTo(*data,Save(kTRUE), Minos(kFALSE), NumCPU(4), SumW2Error(kTRUE));
    
    RooPlot * frame = cosTheta->frame();
    data->plotOn(frame);
    pdf.plotOn(frame);

    if(fitresult->minNll() == fitresult->minNll() && fitresult->minNll()>0 ) {
    fout << i << "\t" << a.getVal() << "\t" << b.getVal() << "\t" << fitresult->minNll() << "\t" << fitresult->status() << "\t" 
         << G000.getVal() << "\t" << G001.getVal() << "\t" << G002.getVal() << "\t" << G003.getVal() << "\t" << G004.getVal() << endl;
    }

    gROOT->ProcessLine(".x ~/lhcb/lhcbStyle.C");
    TCanvas c("fit","fit", 800, 800);
    frame->Draw();
    c.SaveAs("fits/fit"+TString(namestr)+".png");
    c.SaveAs("fits/fit"+TString(namestr)+".pdf");
}

//.........这里部分代码省略.........
//    fitResult->Print("v");




//   // ********* Fix Mass Shift and Fit For Resolution ********** //  
//    ParPassSignalMassShift->setConstant(kTRUE); 
//    ParFailSignalMassShift->setConstant(kTRUE); 
//    ParPassSignalResolution->setConstant(kFALSE); 
//    ParFailSignalResolution->setConstant(kFALSE); 
//    fitResult = totalPdf.fitTo(*dataCombined, RooFit::Save(true), 
//    RooFit::Extended(true), RooFit::PrintLevel(-1));
//    fitResult->Print("v");


//   // ********* Do Final Fit ********** //  
//    ParPassSignalMassShift->setConstant(kFALSE); 
//    ParFailSignalMassShift->setConstant(kFALSE); 
//    ParPassSignalResolution->setConstant(kTRUE); 
//    ParFailSignalResolution->setConstant(kTRUE); 
//    fitResult = totalPdf.fitTo(*dataCombined, RooFit::Save(true), 
//                                             RooFit::Extended(true), RooFit::PrintLevel(-1));
//    fitResult->Print("v");





  double nSignalPass = NumSignalPass->getVal();
  double nSignalFail    = NumSignalFail->getVal();
  double denominator = nSignalPass + nSignalFail;

  printf("\nFit results:\n");
  if( fitResult->status() != 0 ){
    std::cout<<"ERROR: BAD FIT STATUS"<<std::endl;
  }

  printf("    Efficiency = %.4f +- %.4f\n", 
	 ParEfficiency->getVal(), ParEfficiency->getPropagatedError(*fitResult));  
  cout << "Signal Pass: " << nSignalPass << endl;
  cout << "Signal Fail: " << nSignalFail << endl;

  cout << "*********************************************************************\n";
  cout << "Final Parameters\n";
  cout << "*********************************************************************\n";
  PrintParameter(ParNumSignal, label,"ParNumSignal");
  PrintParameter(ParNumBkgPass, label,"ParNumBkgPass");
  PrintParameter(ParNumBkgFail, label, "ParNumBkgFail");
  PrintParameter(ParEfficiency  , label, "ParEfficiency");
  PrintParameter(ParPassBackgroundExpCoefficient , label, "ParPassBackgroundExpCoefficient");
  PrintParameter(ParFailBackgroundExpCoefficient , label, "ParFailBackgroundExpCoefficient");
  PrintParameter(ParPassSignalMassShift , label, "ParPassSignalMassShift");
  PrintParameter(ParFailSignalMassShift , label, "ParFailSignalMassShift");
  PrintParameter(ParPassSignalResolution , label, "ParPassSignalResolution");
  PrintParameter(ParFailSignalResolution , label, "ParFailSignalResolution");
  cout << endl << endl;


  //--------------------------------------------------------------------------------------------------------------
  // Make plots 
  //==============================================================================================================  
  TFile *canvasFile = new TFile("Efficiency_FitResults.root", "UPDATE");


  RooAbsData::ErrorType errorType = RooAbsData::Poisson;

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