C++ RooAbsReal::getVal方法代码示例

   RooFitResult* GenericModel::fitTo(RooDataSet* data)
   {
      // Perform fit of the pseudo-PDF to the data
      // On multi-core machines, this automatically uses all available processor cores

      SafeDelete(fLastFit);
#ifdef WITH_MULTICORE_CPU
      fLastFit = fModelPseudoPDF->fitTo(*data, Save(), NumCPU(WITH_MULTICORE_CPU));
#else
      fLastFit = fModelPseudoPDF->fitTo(*data, Save());
#endif

      SafeDelete(fParamDataHist);
      fParamDataHist = new RooDataHist("params", "params", GetParameters());

      // store weights of component pdfs => distribution of parameters
      fWeights.removeAll();
      const RooArgList& coefs = fModelPseudoPDF->coefList();
      for (int i = 0; i < GetNumberOfDataSets(); i++) {
         RooAbsReal* coef = (RooAbsReal*)coefs.at(i);
         RooRealVar w(Form("w%d", i), Form("Fitted weight of kernel#%d", i), coef->getVal());
         if (coef->InheritsFrom(RooRealVar::Class())) {
            w.setError(((RooRealVar*)coef)->getError());
         } else {
            w.setError(coef->getPropagatedError(*fLastFit));
         }
         fWeights.addClone(w);
         fParamDataHist->set(*GetParametersForDataset(i), w.getVal(), w.getError());
      }

      SafeDelete(fParameterPDF);
      fParameterPDF = new RooHistPdf("paramPDF", "paramPDF", GetParameters(), *fParamDataHist);

      return fLastFit;
   }

void fit( RooAbsReal & chi2, int numberOfBins, const char * outFileNameWithFitResult ){
  TFile * out_root_file = new TFile(outFileNameWithFitResult , "recreate");
  RooMinuit m_tot(chi2) ;
  m_tot.migrad();
  // m_tot.hesse();
  RooFitResult* r_chi2 = m_tot.save() ;
  cout << "==> Chi2 Fit results " << endl ;
  r_chi2->Print("v") ;
  //  r_chi2->floatParsFinal().Print("v") ;
  int NumberOfFreeParameters =  r_chi2->floatParsFinal().getSize() ;
  for (int i =0; i< NumberOfFreeParameters; ++i){
    r_chi2->floatParsFinal()[i].Print();
  }
  cout<<"chi2:" <<chi2.getVal() << ", numberOfBins: " << numberOfBins  << ", NumberOfFreeParameters: " << NumberOfFreeParameters << endl; 
  cout<<"Normalized Chi2   = " << chi2.getVal()/ (numberOfBins - NumberOfFreeParameters)<<endl; 
  r_chi2->Write( ) ;
  delete out_root_file;
}

void Zbi_Zgamma() {

  // Make model for prototype on/off problem
  // Pois(x | s+b) * Pois(y | tau b )
  // for Z_Gamma, use uniform prior on b.
  RooWorkspace* w = new RooWorkspace("w",true);
  w->factory("Poisson::px(x[150,0,500],sum::splusb(s[0,0,100],b[100,0,300]))");
  w->factory("Poisson::py(y[100,0,500],prod::taub(tau[1.],b))");
  w->factory("Uniform::prior_b(b)");

  // construct the Bayesian-averaged model (eg. a projection pdf)
  // p'(x|s) = \int db p(x|s+b) * [ p(y|b) * prior(b) ]
  w->factory("PROJ::averagedModel(PROD::foo(px|b,py,prior_b),b)") ;

  // plot it, blue is averaged model, red is b known exactly
  RooPlot* frame = w->var("x")->frame() ;
  w->pdf("averagedModel")->plotOn(frame) ;
  w->pdf("px")->plotOn(frame,LineColor(kRed)) ;
  frame->Draw() ;

  // compare analytic calculation of Z_Bi
  // with the numerical RooFit implementation of Z_Gamma
  // for an example with x = 150, y = 100

  // numeric RooFit Z_Gamma
  w->var("y")->setVal(100);
  w->var("x")->setVal(150);
  RooAbsReal* cdf = w->pdf("averagedModel")->createCdf(*w->var("x"));
  cdf->getVal(); // get ugly print messages out of the way

  cout << "Hybrid p-value = " << cdf->getVal() << endl;
  cout << "Z_Gamma Significance  = " <<
    PValueToSignificance(1-cdf->getVal()) << endl;

  // analytic Z_Bi
  double Z_Bi = NumberCountingUtils::BinomialWithTauObsZ(150, 100, 1);
  std::cout << "Z_Bi significance estimation: " << Z_Bi << std::endl;

  // OUTPUT
  // Hybrid p-value = 0.999058
  // Z_Gamma Significance  = 3.10804
  // Z_Bi significance estimation: 3.10804
}

///
/// Make an Asimov toy: set all observables set to truth values.
/// The Asimov point needs to be loaded in the combiner before.
/// \param c - combiner which should be set to an asimov toy
///
void GammaComboEngine::setAsimovObservables(Combiner* c)
{
    if ( !c->isCombined() ) {
        cout << "GammaComboEngine::setAsimovObservables() : ERROR : Can't set an Asimov toy before "
             "the combiner is combined. Call combine() first." << endl;
        exit(1);
    }

    // set observables to asimov values in workspace
    RooWorkspace* w = c->getWorkspace();
    TIterator* itObs = c->getObservables()->createIterator();
    while(RooRealVar* pObs = (RooRealVar*) itObs->Next()) {
        // get theory name from the observable name
        TString pThName = pObs->GetName();
        pThName.ReplaceAll("obs","th");
        // get the theory relation
        RooAbsReal* th = w->function(pThName);
        if ( th==0 ) {
            cout << "GammaComboEngine::setAsimovObservables() : ERROR : theory relation not found in workspace: " << pThName << endl;
            exit(1);
        }
        // set the observable to what the theory relation predicts
        pObs->setVal(th->getVal());
    }
    delete itObs;

    // write back the asimov values to the PDF object so that when
    // the PDF is printed, the asimov values show up
    for ( int i=0; i<c->getPdfs().size(); i++ ) {
        PDF_Abs* pdf = c->getPdfs()[i];
        pdf->setObservableSourceString("Asimov");
        TIterator* itObs = pdf->getObservables()->createIterator();
        while(RooRealVar* pObs = (RooRealVar*) itObs->Next()) {
            RooAbsReal* obs =  w->var(pObs->GetName());
            if ( obs==0 ) {
                cout << "GammaComboEngine::setAsimovObservables() : ERROR : observable not found in workspace: " << pObs->GetName() << endl;
                exit(1);
            }
            pdf->setObservable(pObs->GetName(), obs->getVal());
        }
        delete itObs;
    }
}

// grab the initial normalization from a datacard converted in workspace
// with: scripts/text2workspace.py -b -o model.root datacards/hww-12.1fb.mH125.comb_0j1j2j_shape.txt
void fillInitialNorms(RooArgSet *args, std::map<std::string, std::pair<double,double> > &vals, std::string workspace){
  TFile *fw_ =  TFile::Open(workspace.c_str());
  RooWorkspace *ws = (RooWorkspace*)fw_->Get("w");
  TIterator* iter(args->createIterator());
  for (TObject *a = iter->Next(); a != 0; a = iter->Next()) {
    RooAbsReal *rar = (RooAbsReal*)ws->obj(a->GetName());
    std::string name = rar->GetName();
    std::pair<double,double> valE(rar->getVal(),0.0);
    vals.insert( std::pair<std::string,std::pair<double ,double> > (name,valE)) ;
  }
}

double NormalizedIntegral(RooAbsPdf & function, RooRealVar & integrationVar, double lowerLimit, double upperLimit){

  integrationVar.setRange("integralRange",lowerLimit,upperLimit);
  RooAbsReal* integral = function.createIntegral(integrationVar,NormSet(integrationVar),Range("integralRange"));


  double normlizedIntegralValue = integral->getVal();

  //  cout<<normlizedIntegralValue<<endl;


  return normlizedIntegralValue;


}

// get effective sigma from culmalative distribution function
pair<double,double> getEffSigma(RooRealVar *mass, RooAbsPdf *pdf, double wmin=110., double wmax=130., double step=0.002, double epsilon=1.e-4){

  RooAbsReal *cdf = pdf->createCdf(RooArgList(*mass));
  cout << "Computing effSigma...." << endl;
  TStopwatch sw;
  sw.Start();
  double point=wmin;
  vector<pair<double,double> > points;
  
  while (point <= wmax){
    mass->setVal(point);
    if (pdf->getVal() > epsilon){
      points.push_back(pair<double,double>(point,cdf->getVal())); 
    }
    point+=step;
  }
  double low = wmin;
  double high = wmax;
  double width = wmax-wmin;
  for (unsigned int i=0; i<points.size(); i++){
    for (unsigned int j=i; j<points.size(); j++){
      double wy = points[j].second - points[i].second;
      if (TMath::Abs(wy-0.683) < epsilon){
        double wx = points[j].first - points[i].first;
        if (wx < width){
          low = points[i].first;
          high = points[j].first;
          width=wx;
        }
      }
    }
  }
  sw.Stop();
  cout << "effSigma: [" << low << "-" << high << "] = " << width/2. << endl;
  cout << "\tTook: "; sw.Print();
  pair<double,double> result(low,high);
  return result;
}

//.........这里部分代码省略.........
  RooDataSet dataSetZjets  ("dataSetZjets",   "dataSetZjets",   variables, Cut(catCut),           WeightVar(evWeight), Import(*treeZjets));
  RooDataSet dataSetZjetsSB("dataSetZjetsSB", "dataSetZjetsSB", variables, Cut(catCut && sbCut),  WeightVar(evWeight), Import(*treeZjets));  
  RooDataSet dataSetZjetsSG("dataSetZjetsSG", "dataSetZjetsSG", variables, Cut(catCut && sigCut), WeightVar(evWeight), Import(*treeZjets));
  
  // Total events number

  float totalMcEv   = dataSetZjetsSB.sumEntries() + dataSetZjetsSG.sumEntries();
  float totalDataEv = dataSetData.sumEntries();

  RooRealVar nMcEvents("nMcEvents", "nMcEvents", 0., 99999.);
  RooRealVar nDataEvents("nDataEvents", "nDataEvents", 0., 99999.);

  nMcEvents.setVal(totalMcEv);
  nMcEvents.setConstant(true);

  nDataEvents.setVal(totalDataEv);
  nDataEvents.setConstant(true);

  // Signal region jet mass

  RooRealVar constant("constant", "constant", -0.02,  -1.,   0.);
  RooRealVar offset  ("offset",   "offset",     30., -50., 200.);
  RooRealVar width   ("width",    "width",     100.,   0., 200.);

  if( catcut == "1" ) offset.setConstant(true);
  
  RooErfExpPdf model_mJet("model_mJet", "model_mJet", mJet, constant, offset, width);
  RooExtendPdf ext_model_mJet("ext_model_mJet", "ext_model_mJet", model_mJet, nMcEvents);

  RooFitResult* mJet_result = ext_model_mJet.fitTo(dataSetZjets, SumW2Error(true), Extended(true), Range("allRange"), Strategy(2), Minimizer("Minuit2"), Save(1));

  // Side band jet mass

  RooRealVar constantSB("constantSB", "constantSB", constant.getVal(),  -1.,   0.);
  RooRealVar offsetSB  ("offsetSB",   "offsetSB",   offset.getVal(),   -50., 200.);
  RooRealVar widthSB   ("widthSB",    "widthSB",    width.getVal(),      0., 200.);

  offsetSB.setConstant(true);

  RooErfExpPdf model_mJetSB("model_mJetSB", "model_mJetSB", mJet, constantSB, offsetSB, widthSB);
  RooExtendPdf ext_model_mJetSB("ext_model_mJetSB", "ext_model_mJetSB", model_mJetSB, nMcEvents);

  RooFitResult* mJetSB_result = ext_model_mJetSB.fitTo(dataSetZjetsSB, SumW2Error(true), Extended(true), Range("lowSB,highSB"), Strategy(2), Minimizer("Minuit2"), Save(1));

  RooAbsReal* nSIGFit = ext_model_mJetSB.createIntegral(RooArgSet(mJet), NormSet(mJet), Range("signal"));

  float normFactor = nSIGFit->getVal() * totalMcEv;
  
  // Plot the results on a frame

  RooPlot* mJetFrame = mJet.frame();

  dataSetZjetsSB.  plotOn(mJetFrame, Binning(binsmJet));  
  ext_model_mJetSB.plotOn(mJetFrame, Range("allRange"), VisualizeError(*mJetSB_result), FillColor(kYellow));
  dataSetZjetsSB.  plotOn(mJetFrame, Binning(binsmJet));  
  ext_model_mJetSB.plotOn(mJetFrame, Range("allRange"));
  mJetFrame->SetTitle("M_{jet} distribution in Z+jets MC");

  // Alpha ratio part

  mZH.setRange("fullRange", 900., 3000.);

  RooBinning binsmZH(21, 900, 3000);

  RooRealVar a("a", "a",  0., -1.,    1.);
  RooRealVar b("b", "b", 1000,  0., 4000.);

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

   // make a histogram of parameter vs. threshold
   TH1F* histOfThresholds = new TH1F("histOfThresholds","",
                                       parameterScan->numEntries(),
                                       firstPOI->getMin(),
                                       firstPOI->getMax());
   histOfThresholds->GetXaxis()->SetTitle(firstPOI->GetName());
   histOfThresholds->GetYaxis()->SetTitle("Threshold");

   // loop through the points that were tested and ask confidence belt
   // what the upper/lower thresholds were.
   // For FeldmanCousins, the lower cut off is always 0
   for(Int_t i=0; i<parameterScan->numEntries(); ++i){
      tmpPoint = (RooArgSet*) parameterScan->get(i)->clone("temp");
      //cout <<"get threshold"<<endl;
      double arMax = belt->GetAcceptanceRegionMax(*tmpPoint);
      double poiVal = tmpPoint->getRealValue(firstPOI->GetName()) ;
      histOfThresholds->Fill(poiVal,arMax);
   }
   TCanvas* c1 = new TCanvas();
   c1->Divide(2);
   c1->cd(1);
   histOfThresholds->SetMinimum(0);
   histOfThresholds->Draw();
   c1->cd(2);

   // -------------------------------------------------------
   // Now we generate the expected bands and power-constraint

   // First: find parameter point for mu=0, with conditional MLEs for nuisance parameters
   RooAbsReal* nll = mc->GetPdf()->createNLL(*data);
   RooAbsReal* profile = nll->createProfile(*mc->GetParametersOfInterest());
   firstPOI->setVal(0.);
   profile->getVal(); // this will do fit and set nuisance parameters to profiled values
   RooArgSet* poiAndNuisance = new RooArgSet();
   if(mc->GetNuisanceParameters())
      poiAndNuisance->add(*mc->GetNuisanceParameters());
   poiAndNuisance->add(*mc->GetParametersOfInterest());
   w->saveSnapshot("paramsToGenerateData",*poiAndNuisance);
   RooArgSet* paramsToGenerateData = (RooArgSet*) poiAndNuisance->snapshot();
   cout << "\nWill use these parameter points to generate pseudo data for bkg only" << endl;
   paramsToGenerateData->Print("v");


   RooArgSet unconditionalObs;
   unconditionalObs.add(*mc->GetObservables());
   unconditionalObs.add(*mc->GetGlobalObservables()); // comment this out for the original conditional ensemble

   double CLb=0;
   double CLbinclusive=0;

   // Now we generate background only and find distribution of upper limits
   TH1F* histOfUL = new TH1F("histOfUL","",100,0,firstPOI->getMax());
   histOfUL->GetXaxis()->SetTitle("Upper Limit (background only)");
   histOfUL->GetYaxis()->SetTitle("Entries");
   for(int imc=0; imc<nToyMC; ++imc){

      // set parameters back to values for generating pseudo data
      //    cout << "\n get current nuis, set vals, print again" << endl;
      w->loadSnapshot("paramsToGenerateData");
      //    poiAndNuisance->Print("v");

      RooDataSet* toyData = 0;
      // now generate a toy dataset
      if(!mc->GetPdf()->canBeExtended()){
         if(data->numEntries()==1)

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

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

       delete dataFitResultSusyFixed ;






       //-- Construct the new POI parameter.
       RooAbsReal* new_poi_rar(0x0) ;

       new_poi_rar = ws->var( new_poi_name ) ;
       if ( new_poi_rar == 0x0 ) {
          printf("\n\n New POI %s is not a variable.  Trying function.\n\n", new_poi_name ) ;
          new_poi_rar = ws->function( new_poi_name ) ;
          if ( new_poi_rar == 0x0 ) {
             printf("\n\n New POI %s is not a function.  I quit.\n\n", new_poi_name ) ;
             return ;
          }
       } else {
          printf("\n\n     New POI %s is a variable with current value %.1f.\n\n", new_poi_name, new_poi_rar->getVal() ) ;
       }








       if ( npoiPoints <=0 ) {
          printf("\n\n Quitting now.\n\n" ) ;
          return ;
       }


       double startPoiVal = new_poi_rar->getVal() ;



      //--- The RooNLLVar is NOT equivalent to what minuit uses.
  //   RooNLLVar* nll = new RooNLLVar("nll","nll", *likelihood, *rds ) ;
  //   printf("\n\n Nll value, from construction : %.3f\n\n", nll->getVal() ) ;

      //--- output of createNLL IS what minuit uses, so use that.
       RooAbsReal* nll = likelihood -> createNLL( *rds, Verbose(true) ) ;

       RooRealVar* rrv_poiValue = new RooRealVar( "poiValue", "poiValue", 0., -10000., 10000. ) ;
   /// rrv_poiValue->setVal( poiMinVal ) ;
   /// rrv_poiValue->setConstant(kTRUE) ;

       RooRealVar* rrv_constraintWidth = new RooRealVar("constraintWidth","constraintWidth", 0.1, 0.1, 1000. ) ;
       rrv_constraintWidth -> setVal( constraintWidth ) ;

vector<double> FitInvMass(TH1D* histo){

	vector<double> vec;

	gROOT->ProcessLine(".x ~/rootlogon.C");

	int n = histo->GetEntries();
	double w = histo->GetXaxis()->GetBinWidth(1);
	int ndf;

	RooPlot* frame;

	double hmin0 = histo->GetXaxis()->GetXmin();
	double hmax0 = histo->GetXaxis()->GetXmax();

	histo->GetXaxis()->SetRangeUser(hmin0,hmax0);

	// Declare observable x
	RooRealVar x("x","x",hmin0,hmax0) ;
	RooDataHist dh("dh","dh",x,Import(*histo)) ;

	frame = x.frame(Title(histo->GetName())) ;
	dh.plotOn(frame,DataError(RooAbsData::SumW2), MarkerColor(1),MarkerSize(0.9),MarkerStyle(7));  //this will show histogram data points on canvas
	dh.statOn(frame);  //this will display hist stat on canvas


	x.setRange("R0",90.5,91) ;
	x.setRange("R1",70,110) ;
	x.setRange("R2",60,120) ;
	x.setRange("R3",50,130) ;

	RooRealVar mean("mean","mean",91.186/*histo->GetMean()*/, 70.0, 120.0);
	RooRealVar width("width","width",7.5, 0, 30.0);
	RooRealVar sigma("sigma","sigma",0, 0.0, 120.0);


	mean.setRange(88,94);
	width.setRange(0,20);
	sigma.setRange(0,10);

	//Choose the fitting here
	//RooGaussian gauss("gauss","gauss",x,mean,sigma);ndf = 2;
	RooBreitWigner gauss("gauss","gauss",x,mean,width);ndf = 2;
	//RooVoigtian gauss("gauss","gauss",x,mean,width,sigma); ndf = 3;
	
	RooFitResult* filters = gauss.fitTo(dh,Range("R1"),"qr");
	gauss.plotOn(frame,LineColor(4));//this will show fit overlay on canvas
	gauss.paramOn(frame); //this will display the fit parameters on canvas

	//TCanvas* b1 = new TCanvas("b1","b1",1200,800);

	//gPad->SetLeftMargin(0.15);

	//frame->GetXaxis()->SetTitle("Z mass (in GeV/c^{2})");  
	//frame->GetXaxis()->SetTitleOffset(1.2);
	//float binsize = histo->GetBinWidth(1); 
	//frame->Draw() ;
	cout<<"The chi2 is:"<<endl;
	cout<<frame->chiSquare(ndf)<<endl; 
	cout<<" "<<endl;

	//Do the integral

	//Store result in .root file
	frame->Write(histo->GetTitle());

	RooAbsReal* integral = gauss.createIntegral(x, NormSet(x), Range("R1")) ;

	vec.push_back(n*integral->getVal());
	//vec.push_back((double)n);
	vec.push_back((double)frame->chiSquare(ndf));
	
	return vec;
}

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;
    
    


}

//.........这里部分代码省略.........
   // fixed_again.add( *ws1->var("sigma1") );
   //  fixed_again.add( *ws1->var("nbkg_pp") );
   // fixed_again.add( *ws1->var("npow") );
   // fixed_again.add( *ws1->var("muPlusPt") );
   // fixed_again.add( *ws1->var("muMinusPt") );
   // fixed_again.add( *ws1->var("mscale_hi") );
   // fixed_again.add( *ws1->var("mscale_pp") );
   //  
   // ws1->Print();
   cout << "99999" << endl;

   // create signal+background Model Config
   RooStats::ModelConfig sbHypo("SbHypo");
   sbHypo.SetWorkspace( *ws1 );
   sbHypo.SetPdf( *ws1->pdf("joint") );
   sbHypo.SetObservables( obs );
   sbHypo.SetGlobalObservables( globalObs );
   sbHypo.SetParametersOfInterest( poi );
   sbHypo.SetNuisanceParameters( nuis );
   sbHypo.SetPriorPdf( *ws1->pdf("step") ); // this is optional

   // ws1->Print();
   /////////////////////////////////////////////////////////////////////
   RooAbsReal * pNll = sbHypo.GetPdf()->createNLL( *data,NumCPU(10) );
   cout << "111111" << endl;
   RooMinuit(*pNll).migrad(); // minimize likelihood wrt all parameters before making plots
   cout << "444444" << endl;
   RooPlot *framepoi = ((RooRealVar *)poi.first())->frame(Bins(10),Range(0.,0.2),Title("LL and profileLL in raa3"));
   cout << "222222" << endl;
   pNll->plotOn(framepoi,ShiftToZero());
   cout << "333333" << endl;
   
   RooAbsReal * pProfile = pNll->createProfile( globalObs ); // do not profile global observables
   pProfile->getVal(); // this will do fit and set POI and nuisance parameters to fitted values
   pProfile->plotOn(framepoi,LineColor(kRed));
   framepoi->SetMinimum(0);
   framepoi->SetMaximum(3);
   TCanvas *cpoi = new TCanvas();
   cpoi->cd(); framepoi->Draw();
   cpoi->SaveAs("cpoi.pdf");

   ((RooRealVar *)poi.first())->setMin(0.);
   RooArgSet * pPoiAndNuisance = new RooArgSet("poiAndNuisance");
   // pPoiAndNuisance->add(*sbHypo.GetNuisanceParameters());
   // pPoiAndNuisance->add(*sbHypo.GetParametersOfInterest());
   pPoiAndNuisance->add( nuis );
   pPoiAndNuisance->add( poi );
   sbHypo.SetSnapshot(*pPoiAndNuisance);

   RooPlot* xframeSB = pObs->frame(Title("SBhypo"));
   data->plotOn(xframeSB,Cut("dataCat==dataCat::hi"));
   RooAbsPdf *pdfSB = sbHypo.GetPdf();
   RooCategory *dataCat = ws1->cat("dataCat");
   pdfSB->plotOn(xframeSB,Slice(*dataCat,"hi"),ProjWData(*dataCat,*data));
   TCanvas *c1 = new TCanvas();
   c1->cd(); xframeSB->Draw();
   c1->SaveAs("c1.pdf");

   delete pProfile;
   delete pNll;
   delete pPoiAndNuisance;
   ws1->import( sbHypo );
   /////////////////////////////////////////////////////////////////////
   RooStats::ModelConfig bHypo = sbHypo;
   bHypo.SetName("BHypo");
   bHypo.SetWorkspace(*ws1);

//.........这里部分代码省略.........
         printf( "  %s\n", pname ) ;
         rv_R_msigmsb[mbi] = new RooRealVar( pname, pname, R_msigmsb_initialval, 0., 3. ) ;
         rv_R_msigmsb[mbi] -> setConstant( kFALSE ) ;
         rv_R_msigmsb[mbi] -> Print() ;

      } // mbi.

      printf("\n") ;

      sprintf( pname, "sig_strength" ) ;
      RooRealVar* rv_sig_strength = new RooRealVar( pname, pname, 1.0, 0., 10. ) ;
      rv_sig_strength -> setConstant(kFALSE) ;
      rv_sig_strength -> Print() ;
      printf("  %s\n\n", pname ) ;

    //-------------------------------------------------------------------------

     //-- Define all mu parameters.

      printf("\n\n Defining mu parameters.\n\n") ;

      RooAbsReal* rv_mu_bg_msig[bins_of_nb][max_bins_of_met] ;  // first index is number of btags, second is met bin.
      RooAbsReal* rv_mu_bg_msb[bins_of_nb][max_bins_of_met]  ;  // first index is number of btags, second is met bin.

      RooAbsReal* rv_mu_sig_msig[bins_of_nb][max_bins_of_met] ; // first index is number of btags, second is met bin.
      RooAbsReal* rv_mu_sig_msb[bins_of_nb][max_bins_of_met]  ; // first index is number of btags, second is met bin.

      for ( int nbi=0; nbi<bins_of_nb; nbi++ ) {

         for ( int mbi=0; mbi<bins_of_met; mbi++ ) {

            sprintf( pname, "mu_bg_%db_msb_met%d", nbi+2, mbi+1 ) ;
            printf( "  %s\n", pname ) ;
            rv_mu_bg_msb[nbi][mbi] = new RooRealVar( pname, pname, rv_N_msb[nbi][mbi] -> getVal(), 0., 1.e6 ) ;
            rv_mu_bg_msb[nbi][mbi] -> Print() ;



            sprintf( formula, "@0 * @1 * @2" ) ;
            sprintf( pname, "mu_bg_%db_msig_met%d", nbi+2, mbi+1 ) ;
            printf( "  %s\n", pname ) ;
            rv_mu_bg_msig[nbi][mbi] = new RooFormulaVar( pname, formula, RooArgSet( *rv_Rsigsb_corr[nbi][mbi], *rv_R_msigmsb[mbi], *rv_mu_bg_msb[nbi][mbi] ) ) ;
            rv_mu_bg_msig[nbi][mbi] -> Print() ;

            sprintf( formula, "@0 * @1" ) ;
            sprintf( pname, "mu_sig_%db_msig_met%d", nbi+2, mbi+1 ) ;
            printf( "  %s\n", pname ) ;
            rv_mu_sig_msig[nbi][mbi] = new RooFormulaVar( pname, formula, RooArgSet( *rv_sig_strength, *rv_smc_msig[nbi][mbi] ) ) ;
            rv_mu_sig_msig[nbi][mbi] -> Print() ;

            sprintf( formula, "@0 * @1" ) ;
            sprintf( pname, "mu_sig_%db_msb_met%d", nbi+2, mbi+1 ) ;
            printf( "  %s\n", pname ) ;
            rv_mu_sig_msb[nbi][mbi] = new RooFormulaVar( pname, formula, RooArgSet( *rv_sig_strength, *rv_smc_msb[nbi][mbi] ) ) ;
            rv_mu_sig_msb[nbi][mbi] -> Print() ;


         } // mbi.

      } // nbi.

     //-- Finished defining mu parameters.

    //-------------------------------------------------------------------------

     //-- Defining small n's

prepDataFiles(){
//	TDirectory *theDr = (TDirectory*) myFile->Get("eleIDdir");///denom_pt/fit_eff_plots");
	//theDr->ls();
	int myIndex;	
	
	TSystemDirectory dir(thePath, thePath);
	TSystemFile *file;
	TString fname;
	TIter next(dir.GetListOfFiles());
	while ((file=(TSystemFile*)next())) {
		fname = file->GetName();
		if (fname.BeginsWith("TnP")&& fname.Contains("mc")) {
	
			ofstream myfile;

			TFile *myFile = new TFile(fname);
			TIter nextkey(myFile->GetListOfKeys());
			TKey *key;
			while (key = (TKey*)nextkey()) {
				TString theTypeClasse = key->GetClassName();
				TString theNomClasse = key->GetTitle();
				if ( theTypeClasse == "TDirectoryFile"){
					TDirectory *theDr = (TDirectory*) myFile->Get(theNomClasse);
					TIter nextkey2(theDr->GetListOfKeys());
					TKey *key2;
					while (key2 = (TKey*)nextkey2()) {
						TString theTypeClasse2 = key2->GetClassName();
						TString theNomClasse2 = key2->GetTitle();	
						myfile.open (theNomClasse2+".info");
						if ( theTypeClasse == "TDirectoryFile"){
							cout << "avant " << endl;
							TDirectory *theDr2 = (TDirectory*) myFile->Get(theNomClasse+"/"+theNomClasse2);
							cout << "apres " << endl;
							TIter nextkey3(theDr2->GetListOfKeys());
							TKey *key3;
							while (key3 = (TKey*)nextkey3()) {
								TString theTypeClasse3 = key3->GetClassName();
								TString theNomClasse3 = key3->GetTitle();	
								if ((theNomClasse3.Contains("FromMC"))) {

									TString localClasse3 = theNomClasse3;
									localClasse3.ReplaceAll("__","%");
									cout << "apres " << localClasse3 << endl;
									TObjArray* listBin = localClasse3.Tokenize('%');
									TString first = ((TObjString*)listBin->At(0))->GetString();
									TString second = ((TObjString*)listBin->At(2))->GetString();
									myfile << first;
									myfile << " " << second << " ";
									cout << "coucou la on va récupérer le rooFitResult " << endl;

									RooFitResult *theResults = (RooFitResult*) myFile->Get(theNomClasse+"/"+theNomClasse2+"/"+theNomClasse3+"/fitresults");
									theResults->Print();
									RooArgList theParam = theResults->floatParsFinal();
									int taille = theParam.getSize();
									for (int m = 0 ; m < taille ; m++){
										cout << "m=" << m << endl;
									RooAbsArg *theArg = (RooAbsArg*) theParam.at(m);
									RooAbsReal *theReal = (RooAbsReal*) theArg;
										myfile << theReal->getVal() << " " ;
									}		
															
									myfile << "\n";

								}
							}
						}
						myfile.close();

					}
			
				}
			}
			delete myFile;
		}
	
	}

}