C++ RooDataSet::Print方法代码示例

void test01()
{
  // S e t u p   m o d e l 
  // ---------------------

	TFile *openFile = new TFile("/tmp/kyolee/dimuonTree_upsiMiniTree_pA5tev_14nb_Run210498-210909_trigBit1_allTriggers0_pt4.root");
	TTree* tree = (TTree*)openFile->Get("UpsilonTree");	

  // Declare variables with associated name, title, initial value and allowed range
  RooRealVar mass("invariantMass","M_{#mu#mu} [GeV/c]",9.4,7,14) ;
  RooRealVar muPlusPt("muPlusPt","muPlusPt",0,40) ;
  RooRealVar muMinusPt("muMinusPt","muMinusPt",0,40) ;
  RooRealVar mean("mean","mean of gaussian",9.4,8,11) ;
  RooRealVar sigma("sigma","width of gaussian",1,0.1,10) ;

	RooDataSet* data = new RooDataSet("data", "data", RooArgSet(mass,muPlusPt,muMinusPt), Import(*tree), Cut("muPlusPt>4 && muMinusPt>4"));

 	data->Print();

  // Build gaussian p.d.f in terms of x,mean and sigma
  RooGaussian gauss("gauss","gaussian PDF",mass,mean,sigma) ;  

  // Construct plot frame in 'x'
  RooPlot* xframe = mass.frame(Title("Gaussian p.d.f.")) ;

	data->plotOn(xframe);

	/*
  // P l o t   m o d e l   a n d   c h a n g e   p a r a m e t e r   v a l u e s
  // ---------------------------------------------------------------------------

  // Plot gauss in frame (i.e. in x) 
  gauss.plotOn(xframe) ;

  // Change the value of sigma to 3
  sigma.setVal(3) ;

  // Plot gauss in frame (i.e. in x) and draw frame on canvas
  gauss.plotOn(xframe,LineColor(kRed)) ;
*/  

  // F i t   m o d e l   t o   d a t a
  // -----------------------------

  // Fit pdf to data
  gauss.fitTo(*data,Range(9,10)) ;
  gauss.plotOn(xframe) ;

  // Print values of mean and sigma (that now reflect fitted values and errors)
  mean.Print() ;
  sigma.Print() ;


  // Draw all frames on a canvas
  TCanvas* c = new TCanvas("rf101_basics","rf101_basics",800,400) ;
  c->cd(1) ; gPad->SetLeftMargin(0.15) ; xframe->GetYaxis()->SetTitleOffset(1.6) ; xframe->Draw() ;
 
}

void read() {

  // Open File
  
  //TFile *datafile = TFile::Open("~/Documents/uni/LHCb_CharmSummerProj/Gedcode/baryon-lifetimes-2015/data/run-II-data/datafileLambda_TAUmin200fs_max2200fs_Mmin2216_max2356.root");
 
  //TFile *datafile = TFile::Open("~/Documents/uni/LHCb_CharmSummerProj/Gedcode/baryon-lifetimes-2015/data/run-II-data/datafileLambda_TAUmin200fs_max2200fs_Mmin2216_max2356_CutIPCHI2lt3.root"); 

  TFile *datafile = TFile::Open("~/Documents/uni/LHCb_CharmSummerProj/learning_root/DataSetLambda_TAUmin0002_max0022.root"); 

  //TFile *datafile = TFile::Open("~/Documents/uni/LHCb_CharmSummerProj/learning_root/DataSetfromGed_fit_data_wSWeight.root"); 

  // Define dataset
  RooDataSet* ds = (RooDataSet*)datafile->Get("ds") ;

  // Define TAU variable, get limits.
  RooRealVar Lambda_cplus_TAU("Lambda_cplus_TAU","Lambda_cplus_TAU",0.0002 ,0.0022 ,"ns") ;  //real range of interest is [0.00025, 0.002], this is defined later.
  double highestTAU;
  double lowestTAU;
  ds->getRange(Lambda_cplus_TAU, lowestTAU, highestTAU);

  // Define Mass variable, get limits.
  RooRealVar Lambda_cplus_M("Lambda_cplus_M","Lambda_cplus_M",2216 ,2356, "GeV") ; 
  double highestM;
  double lowestM;
  ds->RooAbsData::getRange(Lambda_cplus_M, lowestM, highestM) ;
  
  // Define IPCHI2 variable
  RooRealVar Lambda_cplus_IPCHI2_OWNPV("Lambda_cplus_IPCHI2_OWNPV","Lambda_cplus_IPCHI2_OWNPV",-100 ,100) ; 
  double highestIPCHI2_OWNPV;
  double lowestIPCHI2_OWNPV;
  ds->RooAbsData::getRange(Lambda_cplus_IPCHI2_OWNPV, lowestIPCHI2_OWNPV, highestIPCHI2_OWNPV) ;

  //Lambda_cplus_TAU.setRange("R1",0.00018, 0.0012);




  // Print to Screen
  cout<<endl<<endl<<"************info************"<<endl;

  cout<< "Lowest lifetime value in dataset = "<<lowestTAU<<endl;
  cout<< "Highest lifetime value in dataset = "<<highestTAU<<endl;  
  cout<< "Lowest mass value in dataset (should be: 2216)= "<<lowestM<<endl;
  cout<< "Highest mass value in dataset (should be: 2356)= "<<highestM<<endl; 
  cout<< "Lowest IPCHI2_OWNPV value in dataset = "<<lowestIPCHI2_OWNPV<<endl;
  cout<< "Highest IPCHI2_OWNPV value in dataset = "<<highestIPCHI2_OWNPV<<endl;
 
  cout<<"number of events: "<<endl ; 
  ds->Print();  //number of events in dataset    
  //rmodel->Print(); //results  
  //cout<<"Chi squared ="<<chi2<<endl ;
  
}

///////////////////////////////////////////////////////////////////////////////
//weightDS 
///////////////////////////////////////////////////////////////////////////////
RooDataSet* weightDS(RooDataSet* set, double wgt){

  RooRealVar wt("wt","wt",wgt);
  set->addColumn(wt);

  RooDataSet* wSet = new RooDataSet(set->GetName(),set->GetTitle(),set,*set->get(),0,wt.GetName()) ;
  std::cout << "Returning weighted d.s." << std::endl;
  std::cout << std::endl;
  wSet->Print();
  return wSet;
    
}

//
// single measurement (LM0 or LM1)
//
void RA4Single (StatMethod method, double* sig, double* bkg) {

  // RooWorkspace* wspace = createWorkspace();
  RA4WorkSpace ra4WSpace("wspace",true,true,true);
  ra4WSpace.addChannel(RA4WorkSpace::MuChannel);
  ra4WSpace.finalize();

  double lm0_mc[4] = { 1.09, 7.68, 3.78, 21.13 };
  double lm1_mc[4] = { 0.05 , 0.34 , 1.12 , 3.43 };
  double* lm_mc = sig ? sig : lm0_mc;

  double bkg_mc[4] = {  14.73, 18.20, 8.48, 10.98 };

  ra4WSpace.setBackground(RA4WorkSpace::MuChannel,bkg_mc[0],bkg_mc[1],bkg_mc[2],bkg_mc[3]);
  ra4WSpace.setSignal(RA4WorkSpace::MuChannel,lm_mc[0],lm_mc[1],lm_mc[2],lm_mc[3],1.,1.,1.,1.);
  
  // setBackgrounds(wspace,bkg);
  // setSignal(wspace,lm_mc);

  RooWorkspace* wspace = ra4WSpace.workspace();
  // wspace->Print("v");
  // RooArgSet allVars = wspace->allVars();
  // // allVars.printLatex(std::cout,1);
  // TIterator* it = allVars.createIterator();
  // RooRealVar* var;
  // while ( var=(RooRealVar*)it->Next() ) {
  //   var->Print("v");
  //   var->printValue(std::cout);
  // }

  ////////////////////////////////////////////////////////////
  // Generate toy data
  // generate toy data assuming current value of the parameters
  // import into workspace. 
  // add Verbose() to see how it's being generated
  // wspace->var("s")->setVal(0.);
  // RooDataSet* data =   wspace->pdf("model")->generate(*wspace->set("obs"),1);
  // data->Print("v");
  // // wspace->import(*data);
  // wspace->var("s")->setVal(lm_mc[3]);
  RooDataSet* data = new RooDataSet("data","data",*wspace->set("obs"));
  data->add(*wspace->set("obs"));
  data->Print("v");

  MyLimit limit = computeLimit(wspace,data,method,true);
  std::cout << "Limit [ " << limit.lowerLimit << " , "
	    << limit.upperLimit << " ] ; isIn = " << limit.isInInterval << std::endl;
}

exampleScript()
{
  gSystem->CompileMacro("betaHelperFunctions.h"      ,"kO") ;
  gSystem->CompileMacro("RooNormalFromFlatPdf.cxx"      ,"kO") ;
  gSystem->CompileMacro("RooBetaInverseCDF.cxx"      ,"kO") ;
  gSystem->CompileMacro("RooBetaPrimeInverseCDF.cxx" ,"kO") ;
  gSystem->CompileMacro("RooCorrelatedBetaGeneratorHelper.cxx"  ,"kO") ;
  gSystem->CompileMacro("RooCorrelatedBetaPrimeGeneratorHelper.cxx"  ,"kO") ;
  gSystem->CompileMacro("rooFitBetaHelperFunctions.h","kO") ;

  TFile betaTest("betaTest.root","RECREATE");
  betaTest.cd();
  
  RooWorkspace workspace("workspace");
  TString correlatedName("testVariable");
  TString observables("observables");
  TString nuisances("nuisances");

  RooAbsArg* betaOne = getCorrelatedBetaConstraint(workspace,"betaOne","",
						   0.5 , 0.1 ,
						   observables, nuisances,
						   correlatedName );

  printf("\n\n *** constraint name is %s from betaOne and %s\n\n", betaOne->GetName(), correlatedName.Data() ) ;

  RooAbsArg* betaTwo = getCorrelatedBetaConstraint(workspace,"betaTwo","",
						   0 , 0 ,
						   observables, nuisances,
						   correlatedName );

  RooAbsArg* betaThree = getCorrelatedBetaConstraint(workspace,"betaThree","",
						     0.2 , 0.01 ,
						     observables, nuisances,
						     correlatedName );

  RooAbsArg* betaFour = getCorrelatedBetaConstraint(workspace,"betaFour","",
						    0.7 , 0.1 ,
						    observables, nuisances,
						    correlatedName );

  RooAbsArg* betaFourC = getCorrelatedBetaConstraint(workspace,"betaFourC","",
						    0.7 , 0.1 ,
						    observables, nuisances,
						    correlatedName, kTRUE );

  RooAbsArg* betaPrimeOne = getCorrelatedBetaPrimeConstraint(workspace,"betaPrimeOne","",
							     1.0 , 0.5 ,
							     observables, nuisances,
							     correlatedName );

  RooAbsArg* betaPrimeOneC = getCorrelatedBetaPrimeConstraint(workspace,"betaPrimeOneC","",
							     1.0 , 0.5 ,
							     observables, nuisances,
							     correlatedName, kTRUE );

  RooAbsArg* betaPrimeTwo = getCorrelatedBetaPrimeConstraint(workspace,"betaPrimeTwo","",
							     0.7 , 0.5 ,
							     observables, nuisances,
							     correlatedName );

  RooAbsArg* betaPrimeThree = getCorrelatedBetaPrimeConstraint(workspace,"betaPrimeThree","",
							       0.1 , 0.05 ,
							       observables, nuisances,
							       correlatedName );

  RooAbsArg* betaPrimeFour = getCorrelatedBetaPrimeConstraint(workspace,"betaPrimeFour","",
							      7 , 1 ,
							      observables, nuisances,
							      correlatedName );

  RooRealVar* correlatedParameter = workspace.var(correlatedName);

  RooAbsPdf* normalFromFlat = workspace.pdf(correlatedName+"_Constraint");

  RooDataSet* data = normalFromFlat->generate(RooArgSet(*correlatedParameter),1e5);

  data->addColumn(*normalFromFlat);

  data->addColumn(*betaOne);
  data->addColumn(*betaTwo);
  data->addColumn(*betaThree);
  data->addColumn(*betaFour);
  data->addColumn(*betaFourC);
  
  data->addColumn(*betaPrimeOne);
  data->addColumn(*betaPrimeTwo);
  data->addColumn(*betaPrimeThree);
  data->addColumn(*betaPrimeFour);
  data->addColumn(*betaPrimeOneC);

  data->Print("v");

  workspace.Print() ;

  //Setup Plotting Kluges:

  RooRealVar normalPlotter  (correlatedName+"_Constraint" , correlatedName+"_Constraint"  ,0,1);
  RooPlot* normalPlot = normalPlotter.frame();
  data->plotOn(normalPlot);

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

//.........这里部分代码省略.........
        etaBins.push_back(1.05);
        etaBins.push_back(1.37);
        etaBins.push_back(1.52);
        etaBins.push_back(1.74);
        etaBins.push_back(2.1);
 
    }
    etaBins.push_back(2.4);
    const int nEtaBins = etaBins.size()-1;

    // --- declare cut variables --- //
    RooRealVar muonPt("muonPt","p_{T}",0.0,350.0,"GeV");
    RooRealVar  missPt("missPt","p_{T}^{miss}",0.0,350.0,"GeV");
    RooRealVar  muonMt("muonMt","m_{T}",0.0,350.0,"GeV");
    RooRealVar  muonCharge("muonCharge","charge",-2.0,+2.0);
    RooRealVar  centrality("centrality","centrality",0.,1.0);
    RooRealVar  muonEta("muonEta","muonEta",-3.0,+3.0);
    RooRealVar  w("w","w",0.0,10.0);

    TString sCutsSig = "";
    sCutsSig = "abs(muonEta)>0.1&&abs(muonEta)<2.4&&muonPt>25.0&&missPt>25.0&&missPt<9000.0&&muonMt>40.0&&centrality>0.&&centrality<0.8";
    //cross-check w/o mpt cut
    //sCutsSig = "abs(muonEta)>0.1&&abs(muonEta)<2.4&&muonPt>25.0&&missPt>0.0&&missPt<9000.0&&muonMt>40.0&&centrality>0.&&centrality<0.8";
    if(doSystematic) sCutsSig = "muonPt>25.0&&abs(muonEta)>0.1&&abs(muonEta)<2.4&&centrality>0.&&centrality<0.8";
    std::cout << " Signals cuts: "<< sCutsSig << std::endl;

    RooArgList muonTauArgList(muonEta,muonPt,missPt,muonMt,centrality);

	RooFormulaVar cutsSig("cutsSig", "cutsSig", sCutsSig, muonTauArgList);

    RooCategory chargeCategory("chargeCategory","chargeCategory") ;
    chargeCategory.defineType("muMinus",-1) ;
    chargeCategory.defineType("muPlus",1) ;
    chargeCategory.Print();

    RooArgSet muonArgSet(muonPt,missPt,muonMt,muonEta,centrality,chargeCategory,w);

    ///Fill dataset
    RooDataSet* mcTauSet = fillHIWTauDataSet(baseString,fileNameIn+".root",muonArgSet); mcTauSet->Print(); 
    ///apply selection cuts
	RooDataSet* mcTauCutSet = (RooDataSet*)mcTauSet->reduce(Cut(cutsSig)); mcTauCutSet->Print();
    ///Create subsets
    RooDataSet* mcTauSubSet[nPtBins][nEtaBins][nCentralityBins];
    TH1F* hMcTauSubSet[nPtBins][nEtaBins][nCentralityBins];
    RooDataSet* mcTauCutSubSet[nPtBins][nEtaBins][nCentralityBins];
    TH1F* hMcTauCutSubSet[nPtBins][nEtaBins][nCentralityBins];

  	RooBinning b = RooBinning(170,0.0,510.0); // 3 GeV per bin
    std::cout << "Creating subsets..." << std::endl;
    char cTauGen[50],cTauCut[50];
    float cutEff;
	for ( int i = 0; i < nPtBins; i++ ) {
	  for ( int j = 0; j < nEtaBins; j++ ) {
	    for ( int k = 0; k < nCentralityBins; k++ ){

            sprintf(cTauGen,"hTauGen_Eta%i_Cent_%i",j,k);
            sprintf(cTauCut,"hTauCut_Eta%i_Cent_%i",j,k);

            TString sCutEff = "effForTauStudy_cent"; sCutEff+= k;
            TGraphAsymmErrors* grCutEff = (TGraphAsymmErrors*)_fCutEff->Get(sCutEff);
            ///Differential in centrality and eta
            if(doCentrality&&doEta) cutEff = grCutEff->GetY()[j];
            ///Averaged over all centrality and eta
            else cutEff =0.738509 ;
            ///bin in pt,eta,and centrality
            mcTauSubSet[i][j][k] = selectPtEtaCentrality(mcTauSet,ptBins[i], ptBins[i+1], etaBins[j], etaBins[j+1],centralityBins[k], centralityBins[k+1], true);

void rf303_conditional()
{
   // S e t u p   c o m p o s e d   m o d e l   g a u s s ( x , m ( y ) , s )
   // -----------------------------------------------------------------------

   // Create observables
   RooRealVar x("x","x",-10,10) ;
   RooRealVar y("y","y",-10,10) ;

   // Create function f(y) = a0 + a1*y
   RooRealVar a0("a0","a0",-0.5,-5,5) ;
   RooRealVar a1("a1","a1",-0.5,-1,1) ;
   RooPolyVar fy("fy","fy",y,RooArgSet(a0,a1)) ;

   // Create gauss(x,f(y),s)
   RooRealVar sigma("sigma","width of gaussian",0.5,0.1,2.0) ;
   RooGaussian model("model","Gaussian with shifting mean",x,fy,sigma) ;


   // Obtain fake external experimental dataset with values for x and y
   RooDataSet* expDataXY = makeFakeDataXY() ;



   // G e n e r a t e   d a t a   f r o m   c o n d i t i o n a l   p . d . f   m o d e l ( x | y )
   // ---------------------------------------------------------------------------------------------

   // Make subset of experimental data with only y values
   RooDataSet* expDataY= (RooDataSet*) expDataXY->reduce(y) ;

   // Generate 10000 events in x obtained from _conditional_ model(x|y) with y values taken from experimental data
   RooDataSet *data = model.generate(x,ProtoData(*expDataY)) ;
   data->Print() ;



   // F i t   c o n d i t i o n a l   p . d . f   m o d e l ( x | y )   t o   d a t a
   // ---------------------------------------------------------------------------------------------

   model.fitTo(*expDataXY,ConditionalObservables(y)) ;



   // P r o j e c t   c o n d i t i o n a l   p . d . f   o n   x   a n d   y   d i m e n s i o n s
   // ---------------------------------------------------------------------------------------------

   // Plot x distribution of data and projection of model on x = 1/Ndata sum(data(y_i)) model(x;y_i)
   RooPlot* xframe = x.frame() ;
   expDataXY->plotOn(xframe) ;
   model.plotOn(xframe,ProjWData(*expDataY)) ;


   // Speed up (and approximate) projection by using binned clone of data for projection
   RooAbsData* binnedDataY = expDataY->binnedClone() ;
   model.plotOn(xframe,ProjWData(*binnedDataY),LineColor(kCyan),LineStyle(kDotted)) ;


   // Show effect of projection with too coarse binning
   ((RooRealVar*)expDataY->get()->find("y"))->setBins(5) ;
   RooAbsData* binnedDataY2 = expDataY->binnedClone() ;
   model.plotOn(xframe,ProjWData(*binnedDataY2),LineColor(kRed)) ;


   // Make canvas and draw RooPlots
   new TCanvas("rf303_conditional","rf303_conditional",600, 460);
   gPad->SetLeftMargin(0.15) ; xframe->GetYaxis()->SetTitleOffset(1.2) ; xframe->Draw() ;

}

void rf403_weightedevts()
{
  // C r e a t e   o b s e r v a b l e   a n d   u n w e i g h t e d   d a t a s e t 
  // -------------------------------------------------------------------------------

  // Declare observable
  RooRealVar x("x","x",-10,10) ;
  x.setBins(40) ;

  // Construction a uniform pdf
  RooPolynomial p0("px","px",x) ;

  // Sample 1000 events from pdf
  RooDataSet* data = p0.generate(x,1000) ;

 

  // C a l c u l a t e   w e i g h t   a n d   m a k e   d a t a s e t   w e i g h t e d 
  // -----------------------------------------------------------------------------------

  // Construct formula to calculate (fake) weight for events
  RooFormulaVar wFunc("w","event weight","(x*x+10)",x) ;

  // Add column with variable w to previously generated dataset
  RooRealVar* w = (RooRealVar*) data->addColumn(wFunc) ;

  // Dataset d is now a dataset with two observable (x,w) with 1000 entries
  data->Print() ;

  // Instruct dataset wdata in interpret w as event weight rather than as observable
  RooDataSet wdata(data->GetName(),data->GetTitle(),data,*data->get(),0,w->GetName()) ;

  // Dataset d is now a dataset with one observable (x) with 1000 entries and a sum of weights of ~430K
  wdata.Print() ;



  // U n b i n n e d   M L   f i t   t o   w e i g h t e d   d a t a 
  // ---------------------------------------------------------------

  // Construction quadratic polynomial pdf for fitting
  RooRealVar a0("a0","a0",1) ;
  RooRealVar a1("a1","a1",0,-1,1) ;
  RooRealVar a2("a2","a2",1,0,10) ;
  RooPolynomial p2("p2","p2",x,RooArgList(a0,a1,a2),0) ;

  // Fit quadratic polynomial to weighted data

  // NOTE: A plain Maximum likelihood fit to weighted data does in general 
  //       NOT result in correct error estimates, unless individual
  //       event weights represent Poisson statistics themselves.
  //       
  // Fit with 'wrong' errors
  RooFitResult* r_ml_wgt = p2.fitTo(wdata,Save()) ;
  
  // A first order correction to estimated parameter errors in an 
  // (unbinned) ML fit can be obtained by calculating the
  // covariance matrix as
  //
  //    V' = V C-1 V
  //
  // where V is the covariance matrix calculated from a fit
  // to -logL = - sum [ w_i log f(x_i) ] and C is the covariance
  // matrix calculated from -logL' = -sum [ w_i^2 log f(x_i) ] 
  // (i.e. the weights are applied squared)
  //
  // A fit in this mode can be performed as follows:

  RooFitResult* r_ml_wgt_corr = p2.fitTo(wdata,Save(),SumW2Error(kTRUE)) ;



  // P l o t   w e i g h e d   d a t a   a n d   f i t   r e s u l t 
  // ---------------------------------------------------------------

  // Construct plot frame
  RooPlot* frame = x.frame(Title("Unbinned ML fit, binned chi^2 fit to weighted data")) ;

  // Plot data using sum-of-weights-squared error rather than Poisson errors
  wdata.plotOn(frame,DataError(RooAbsData::SumW2)) ;

  // Overlay result of 2nd order polynomial fit to weighted data
  p2.plotOn(frame) ;



  // M L  F i t   o f   p d f   t o   e q u i v a l e n t  u n w e i g h t e d   d a t a s e t
  // -----------------------------------------------------------------------------------------
  
  // Construct a pdf with the same shape as p0 after weighting
  RooGenericPdf genPdf("genPdf","x*x+10",x) ;

  // Sample a dataset with the same number of events as data
  RooDataSet* data2 = genPdf.generate(x,1000) ;

  // Sample a dataset with the same number of weights as data
  RooDataSet* data3 = genPdf.generate(x,43000) ;

  // Fit the 2nd order polynomial to both unweighted datasets and save the results for comparison
  RooFitResult* r_ml_unw10 = p2.fitTo(*data2,Save()) ;
//.........这里部分代码省略.........

//put very small data entries in a binned dataset to avoid unphysical pdfs, specifically for H->ZZ->4l
RooDataSet* makeData(RooDataSet* orig, RooSimultaneous* simPdf, const RooArgSet* observables, RooRealVar* firstPOI, double mass, double& mu_min)
{

  double max_soverb = 0;

  mu_min = -10e9;

  map<string, RooDataSet*> data_map;
  firstPOI->setVal(0);
  RooCategory* cat = (RooCategory*)&simPdf->indexCat();
  TList* datalist = orig->split(*(RooAbsCategory*)cat, true);
  TIterator* dataItr = datalist->MakeIterator();
  RooAbsData* ds;
  RooRealVar* weightVar = new RooRealVar("weightVar","weightVar",1);
  while ((ds = (RooAbsData*)dataItr->Next()))
  {
    string typeName(ds->GetName());
    cat->setLabel(typeName.c_str());
    RooAbsPdf* pdf = simPdf->getPdf(typeName.c_str());
    cout << "pdf: " << pdf << endl;
    RooArgSet* obs = pdf->getObservables(observables);
    cout << "obs: " << obs << endl;

    RooArgSet obsAndWeight(*obs, *weightVar);
    obsAndWeight.add(*cat);
    stringstream datasetName;
    datasetName << "newData_" << typeName;
    RooDataSet* thisData = new RooDataSet(datasetName.str().c_str(),datasetName.str().c_str(), obsAndWeight, WeightVar(*weightVar));

    RooRealVar* firstObs = (RooRealVar*)obs->first();
    //int ibin = 0;
    int nrEntries = ds->numEntries();
    for (int ib=0;ib<nrEntries;ib++)
    {
      const RooArgSet* event = ds->get(ib);
      const RooRealVar* thisObs = (RooRealVar*)event->find(firstObs->GetName());
      firstObs->setVal(thisObs->getVal());

      firstPOI->setVal(0);
      double b = pdf->expectedEvents(*firstObs)*pdf->getVal(obs);
      firstPOI->setVal(1);
      double s = pdf->expectedEvents(*firstObs)*pdf->getVal(obs) - b;

      if (s > 0)
      {
	mu_min = max(mu_min, -b/s);
	double soverb = s/b;
	if (soverb > max_soverb)
	{
	  max_soverb = soverb;
	  cout << "Found new max s/b: " << soverb << " in pdf " << pdf->GetName() << " at m = " << thisObs->getVal() << endl;
	}
      }

      if (b == 0 && s != 0)
      {
	cout << "Expecting non-zero signal and zero bg at m=" << firstObs->getVal() << " in pdf " << pdf->GetName() << endl;
      }
      if (s+b <= 0) 
      {
	cout << "expecting zero" << endl;
	continue;
      }


      double weight = ds->weight();
      if ((typeName.find("ATLAS_H_4mu") != string::npos || 
	   typeName.find("ATLAS_H_4e") != string::npos ||
	   typeName.find("ATLAS_H_2mu2e") != string::npos ||
	   typeName.find("ATLAS_H_2e2mu") != string::npos) && fabs(firstObs->getVal() - mass) < 10 && weight == 0)
      {
	cout << "adding event: " << firstObs->getVal() << endl;
	thisData->add(*event, pow(10., -9.));
      }
      else
      {
	//weight = max(pow(10.0, -9), weight);
	thisData->add(*event, weight);
      }
    }



    data_map[string(ds->GetName())] = (RooDataSet*)thisData;
  }

  
  RooDataSet* newData = new RooDataSet("newData","newData",RooArgSet(*observables, *weightVar), 
				       Index(*cat), Import(data_map), WeightVar(*weightVar));

  orig->Print();
  newData->Print();
  //newData->tree()->Scan("*");
  return newData;

}

//.........这里部分代码省略.........
  RooRealVar BDT_response("BDT_response","BDT_response",-1.0,1.0);
  RooRealVar gamma_CL("gamma_CL","gamma_CL",0.1,1.0);
  RooArgSet Args(MCBMass,MCEtaMass,BDT_response,gamma_CL);

  RooDataSet* MCData12 = new RooDataSet("MCData12","MCData12",Args,Import(*MCFlatInputTree12));
  
  std::cout <<" Data File 12 Loaded"<<std::endl;
  
  RooDataSet* MCData11 = new RooDataSet("MCData11","MCData11",Args,Import(*MCFlatInputTree11));

  std::cout<<" Data File 11 loaded"<<std::endl;

  RooDataSet* MCDataAll= new RooDataSet("MCDataAll","MCDataAll",Args);

  MCDataAll->append(*MCData12);
  MCDataAll->append(*MCData11);
  
  RooPlot* massFrame = MCBMass.frame(Title("Data Import Check"),Bins(50));
  MCDataAll->plotOn(massFrame);
  
  RooPlot *BDTFrame = BDT_response.frame(Title("BDT Cut Check"),Bins(50));
  MCDataAll->plotOn(BDTFrame);
  TCanvas C;
  C.Divide(2,1);
  C.cd(1);
  massFrame->Draw();
  C.cd(2);
  BDTFrame->Draw();
  C.SaveAs("ImportChecks.eps");

  //________________________________MAKE MCROODATACATEGORIES__________________________________

  RooDataSet* MCBData=(RooDataSet*)MCDataAll->reduce(RooArgSet(MCBMass));
  MCBData->Print("v");
  
  RooDataSet* MCEtaData=(RooDataSet*)MCDataAll->reduce(RooArgSet(MCEtaMass));
  MCEtaData->Print("v");

  RooCategory MCMassType("MCMassType","MCMassType") ;
  MCMassType.defineType("B") ;
  MCMassType.defineType("Eta") ;
  
  // Construct combined dataset in (x,sample)
  RooDataSet MCcombData("MCcombData","MC combined data",Args,Index(MCMassType),Import("B",*MCBData),Import("Eta",*MCEtaData));

  
  //=============================================== MC FIT MODEL===================================
  
  RooRealVar Mean("Mean","Mean",5279.29,5276.0,5284.00);
  RooRealVar Sigma("Sigma","Sigma",20.54,17.0,24.8);
  RooRealVar LAlpha("LAlpha","LAlpha",-1.064,-2.5,0.0);
  RooRealVar RAlpha("RAlpha","RAlpha",1.88,0.0,5.0);
  RooRealVar LN("LN","LN",13.0,0.0,40.0);
  RooRealVar RN("RN","RN",2.56,0.0,6.0);

  RooCBShape CBLeft("CBLeft","CBLeft",MCBMass,Mean,Sigma,LAlpha,LN);
  
  RooCBShape CBRight("CBRight","CBRight",MCBMass,Mean,Sigma,RAlpha,RN);

  RooRealVar FitFraction("FitFraction","FitFraction",0.5,0.0,1.0);
  RooAddPdf DCB("DCB","DCB",RooArgList(CBRight,CBLeft),FitFraction);

  RooRealVar SignalYield("SignalYield","SignalYield",4338.0,500.0,10000.0);
  //  RooExtendPdf ExtDCB("ExtDCB","ExtDCB",DCB,SignalYield);
  
  //==============================ETA DCB ++++++++++++++++++++++++++++++

int main(int argc, char* argv[]) {

 doofit::builder::EasyPdf *epdf = new doofit::builder::EasyPdf();

    

 epdf->Var("sig_yield");
 epdf->Var("sig_yield").setVal(153000);
 epdf->Var("sig_yield").setConstant(false);
 //decay time
 epdf->Var("obsTime");
 epdf->Var("obsTime").SetTitle("t_{#kern[-0.2]{B}_{#kern[-0.1]{ d}}^{#kern[-0.1]{ 0}}}");
 epdf->Var("obsTime").setUnit("ps");
 epdf->Var("obsTime").setRange(0.,16.);

 // tag, respectively the initial state of the produced B meson
 epdf->Cat("obsTag");
 epdf->Cat("obsTag").defineType("B_S",1);
 epdf->Cat("obsTag").defineType("Bbar_S",-1);

  //finalstate
  epdf->Cat("catFinalState");
  epdf->Cat("catFinalState").defineType("f",1);
  epdf->Cat("catFinalState").defineType("fbar",-1);

  epdf->Var("obsEtaOS");
  epdf->Var("obsEtaOS").setRange(0.0,0.5);


  std::vector<double> knots;
    knots.push_back(0.07);
    knots.push_back(0.10);
    knots.push_back(0.138);
    knots.push_back(0.16);
    knots.push_back(0.23);
    knots.push_back(0.28);
    knots.push_back(0.35);
    knots.push_back(0.42);
    knots.push_back(0.44);
    knots.push_back(0.48);
    knots.push_back(0.5);

  // empty arg list for coefficients
  RooArgList* list = new RooArgList();

  // create first coefficient
  RooRealVar* coeff_first = &(epdf->Var("parCSpline1"));
  coeff_first->setRange(0,10000);
  coeff_first->setVal(1);
  coeff_first->setConstant(false);
  list->add( *coeff_first );

  for (unsigned int i=1; i <= knots.size(); ++i){
    std::string number = boost::lexical_cast<std::string>(i);
    RooRealVar* coeff = &(epdf->Var("parCSpline"+number));
    coeff->setRange(0,10000);
    coeff->setVal(1);
    coeff->setConstant(false);
    list->add( *coeff );
    }
  
  // create last coefficient
  RooRealVar* coeff_last = &(epdf->Var("parCSpline"+boost::lexical_cast<std::string>(knots.size())));
  coeff_last->setRange(0,10000);
  coeff_last->setVal(1);
  coeff_last->setConstant(false);
  list->add( *coeff_last );
  list->Print();
  // define Eta PDF
  doofit::roofit::pdfs::DooCubicSplinePdf splinePdf("splinePdf",epdf->Var("obsEtaOS"),knots,*list,0,0.5);
  
  //Berechne die Tagging Assymetrie
  epdf->Var("p0");
  epdf->Var("p0").setVal(0.369);
  epdf->Var("p0").setConstant(true);

  epdf->Var("p1");
  epdf->Var("p1").setVal(0.952);
  epdf->Var("p1").setConstant(true);

  epdf->Var("delta_p0");
  epdf->Var("delta_p0").setVal(0.019);
  epdf->Var("delta_p0").setConstant(true);

  epdf->Var("delta_p1");
  epdf->Var("delta_p1").setVal(-0.012);
  epdf->Var("delta_p1").setConstant(true);

  epdf->Var("etamean");
  epdf->Var("etamean").setVal(0.365);
  epdf->Var("etamean").setConstant(true);

  epdf->Formula("omega","@0 [email protected]/2 +(@[email protected]/2)*(@[email protected])", RooArgList(epdf->Var("p0"),epdf->Var("delta_p0"),epdf->Var("p1"),epdf->Var("delta_p1"),epdf->Var("obsEtaOS"),epdf->Var("etamean")));
  epdf->Formula("omegabar","@0 [email protected]/2 +(@[email protected]/2)*(@[email protected])", RooArgList(epdf->Var("p0"),epdf->Var("delta_p0"),epdf->Var("p1"),epdf->Var("delta_p1"),epdf->Var("obsEtaOS"),epdf->Var("etamean")));
      


  //Koeffizienten
  DecRateCoeff *coeff_c = new DecRateCoeff("coef_cos","coef_cos",DecRateCoeff::CPOdd,epdf->Cat("catFinalState"),epdf->Cat("obsTag"),epdf->Var("C_f"),epdf->Var("C_fbar"),epdf->Var("obsEtaOS"),splinePdf,epdf->Var("tageff"),epdf->Real("omega"),epdf->Real("omegabar"),epdf->Var("asym_prod"),epdf->Var("asym_det"),epdf->Var("asym_tageff"));
  DecRateCoeff *coeff_s = new DecRateCoeff("coef_sin","coef_sin",DecRateCoeff::CPOdd,epdf->Cat("catFinalState"),epdf->Cat("obsTag"),epdf->Var("S_f"),epdf->Var("S_fbar"),epdf->Var("obsEtaOS"),splinePdf,epdf->Var("tageff"),epdf->Real("omega"),epdf->Real("omegabar"),epdf->Var("asym_prod"),epdf->Var("asym_det"),epdf->Var("asym_tageff"));
//.........这里部分代码省略.........

   void constrained_scan( const char* wsfile = "outputfiles/ws.root",
                          const char* new_poi_name="mu_bg_4b_msig_met1",
                          double constraintWidth=1.5,
                          int npoiPoints = 20,
                          double poiMinVal = 0.,
                          double poiMaxVal = 10.0,
                          double ymax = 9.,
                          int verbLevel=1  ) {

      TString outputdir("outputfiles") ;

      gStyle->SetOptStat(0) ;

      TFile* wstf = new TFile( wsfile ) ;
      RooWorkspace* ws = dynamic_cast<RooWorkspace*>( wstf->Get("ws") );
      ws->Print() ;

      RooDataSet* rds = (RooDataSet*) ws->obj( "hbb_observed_rds" ) ;
      cout << "\n\n\n  ===== RooDataSet ====================\n\n" << endl ;
      rds->Print() ;
      rds->printMultiline(cout, 1, kTRUE, "") ;

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

      RooAbsPdf* likelihood = ws->pdf("likelihood") ;
      if ( likelihood == 0x0 ) { printf("\n\n *** can't find likelihood in workspace.\n\n" ) ; return ; }
      printf("\n\n Likelihood:\n") ;
      likelihood -> Print() ;



      /////rv_sig_strength -> setConstant( kFALSE ) ;
      rv_sig_strength -> setVal(0.) ;
      rv_sig_strength -> setConstant( kTRUE ) ;

      likelihood->fitTo( *rds, Save(false), PrintLevel(0), Hesse(true), Strategy(1) ) ;
      //RooFitResult* fitResult = likelihood->fitTo( *rds, Save(true), PrintLevel(0), Hesse(true), Strategy(1) ) ;
      //double minNllSusyFloat = fitResult->minNll() ;
      //double susy_ss_atMinNll = rv_sig_strength -> getVal() ;

      RooMsgService::instance().getStream(1).removeTopic(Minimization) ;
      RooMsgService::instance().getStream(1).removeTopic(Fitting) ;



     //-- 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 Found it.\n\n") ;
         }
      } else {
         printf("\n\n     New POI %s is a variable with current value %.1f.\n\n", new_poi_name, new_poi_rar->getVal() ) ;
      }

       double startPoiVal = new_poi_rar->getVal() ;


       RooAbsReal* nll = likelihood -> createNLL( *rds, Verbose(true) ) ;

       RooRealVar* rrv_poiValue = new RooRealVar( "poiValue", "poiValue", 0., -10000., 10000. ) ;

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


       RooMinuit* rminuit( 0x0 ) ;


       RooMinuit* rminuit_uc = new RooMinuit( *nll  ) ;

       rminuit_uc->setPrintLevel(verbLevel-1) ;
       rminuit_uc->setNoWarn() ;

       rminuit_uc->migrad() ;
       rminuit_uc->hesse() ;

       RooFitResult* rfr_uc = rminuit_uc->fit("mr") ;

       double floatParInitVal[10000] ;
       char   floatParName[10000][100] ;
       int nFloatParInitVal(0) ;
       RooArgList ral_floats = rfr_uc->floatParsFinal() ;
       TIterator* floatParIter = ral_floats.createIterator() ;
       {
          RooRealVar* par ;
          while ( (par = (RooRealVar*) floatParIter->Next()) ) {
             sprintf( floatParName[nFloatParInitVal], "%s", par->GetName() ) ;
             floatParInitVal[nFloatParInitVal] = par->getVal() ;
             nFloatParInitVal++ ;
          }
//.........这里部分代码省略.........

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

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

   void fitqual_plots( const char* wsfile = "outputfiles/ws.root", const char* plottitle="" ) {

      TText* tt_title = new TText() ;
      tt_title -> SetTextAlign(33) ;

      gStyle -> SetOptStat(0) ;
      gStyle -> SetLabelSize( 0.06, "y" ) ;
      gStyle -> SetLabelSize( 0.08, "x" ) ;
      gStyle -> SetLabelOffset( 0.010, "y" ) ;
      gStyle -> SetLabelOffset( 0.010, "x" ) ;
      gStyle -> SetTitleSize( 0.07, "y" ) ;
      gStyle -> SetTitleSize( 0.05, "x" ) ;
      gStyle -> SetTitleOffset( 1.50, "x" ) ;
      gStyle -> SetTitleH( 0.07 ) ;
      gStyle -> SetPadLeftMargin( 0.15 ) ;
      gStyle -> SetPadBottomMargin( 0.15 ) ;
      gStyle -> SetTitleX( 0.10 ) ;

      gDirectory->Delete("h*") ;

      TFile* wstf = new TFile( wsfile ) ;

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

      int bins_of_met = TMath::Nint( ws->var("bins_of_met")->getVal()  ) ;
      printf("\n\n Bins of MET : %d\n\n", bins_of_met ) ;

      int bins_of_nb = TMath::Nint( ws->var("bins_of_nb")->getVal()  ) ;
      printf("\n\n Bins of nb : %d\n\n", bins_of_nb ) ;

      int nb_lookup[10] ;
      if ( bins_of_nb == 2 ) {
         nb_lookup[0] = 2 ;
         nb_lookup[1] = 4 ;
      } else if ( bins_of_nb == 3 ) {
         nb_lookup[0] = 2 ;
         nb_lookup[1] = 3 ;
         nb_lookup[2] = 4 ;
      }

      TCanvas* cfq1 = (TCanvas*) gDirectory->FindObject("cfq1") ;
      if ( cfq1 == 0x0 ) {
         if ( bins_of_nb == 3 ) {
            cfq1 = new TCanvas("cfq1","hbb fit", 700, 1000 ) ;
         } else if ( bins_of_nb == 2 ) {
            cfq1 = new TCanvas("cfq1","hbb fit", 700, 750 ) ;
         } else {
            return ;
         }
      }

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

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

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

      rds->Print() ;
      rds->printMultiline(cout, 1, kTRUE, "") ;

      RooAbsPdf* likelihood = modelConfig->GetPdf() ;

      ///RooFitResult* fitResult = likelihood->fitTo( *rds, Save(true), PrintLevel(0) ) ;
      RooFitResult* fitResult = likelihood->fitTo( *rds, Save(true), PrintLevel(3) ) ;
      fitResult->Print() ;


      char hname[1000] ;
      char htitle[1000] ;
      char pname[1000] ;




     //-- unpack observables.

      int obs_N_msig[10][50] ; // first index is n btags, second is met bin.
      int obs_N_msb[10][50]  ; // first index is n btags, second is met bin.

      const RooArgSet* dsras = rds->get() ;
      TIterator* obsIter = dsras->createIterator() ;
      while ( RooRealVar* obs = (RooRealVar*) obsIter->Next() ) {
         for ( int nbi=0; nbi<bins_of_nb; nbi++ ) {
            for ( int mbi=0; mbi<bins_of_met; mbi++ ) {
               sprintf( pname, "N_%db_msig_met%d", nb_lookup[nbi], mbi+1 ) ;
               if ( strcmp( obs->GetName(), pname ) == 0 ) { obs_N_msig[nbi][mbi] = TMath::Nint( obs -> getVal() ) ; }
               sprintf( pname, "N_%db_msb_met%d", nb_lookup[nbi], mbi+1 ) ;
               if ( strcmp( obs->GetName(), pname ) == 0 ) { obs_N_msb[nbi][mbi] = TMath::Nint( obs -> getVal() ) ; }
            } // mbi.
         } // nbi.
      } // obs iterator.


      printf("\n\n") ;
      for ( int nbi=0; nbi<bins_of_nb; nbi++ ) {
         printf(" nb=%d :  ", nb_lookup[nbi] ) ;
         for ( int mbi=0; mbi<bins_of_met; mbi++ ) {
            printf("  sig=%3d, sb=%3d  |", obs_N_msig[nbi][mbi], obs_N_msb[nbi][mbi] ) ;
//.........这里部分代码省略.........

float ComputeTestStat(TString wsfile, double mu_susy_sig_val) {

  gROOT->Reset();

  TFile* wstf = new TFile( wsfile ) ;

  RooWorkspace* ws = dynamic_cast<RooWorkspace*>( wstf->Get("ws") );
  ws->Print() ;
  
  ModelConfig* modelConfig = (ModelConfig*) ws->obj( "SbModel" ) ;
  
  modelConfig->Print() ;

  RooDataSet* rds = (RooDataSet*) ws->obj( "ra2b_observed_rds" ) ;
  
  rds->Print() ;
  rds->printMultiline(cout, 1, kTRUE, "") ;
  
  RooAbsPdf* likelihood = modelConfig->GetPdf() ;
  
  RooRealVar* rrv_mu_susy_sig = ws->var("mu_susy_all0lep") ;
  if ( rrv_mu_susy_sig == 0x0 ) {
    printf("\n\n\n *** can't find mu_susy_all0lep in workspace.  Quitting.\n\n\n") ;
    return ;
  } else {
    printf(" current value is : %8.3f\n", rrv_mu_susy_sig->getVal() ) ; cout << flush ;
    rrv_mu_susy_sig->setConstant(kFALSE) ;
  }

  /*
  // check the impact of varying the qcd normalization:

  RooRealVar *rrv_qcd_0lepLDP_ratioH1 = ws->var("qcd_0lepLDP_ratio_H1");
  RooRealVar *rrv_qcd_0lepLDP_ratioH2 = ws->var("qcd_0lepLDP_ratio_H2");
  RooRealVar *rrv_qcd_0lepLDP_ratioH3 = ws->var("qcd_0lepLDP_ratio_H3");
  
  rrv_qcd_0lepLDP_ratioH1->setVal(0.3);
  rrv_qcd_0lepLDP_ratioH2->setVal(0.3);
  rrv_qcd_0lepLDP_ratioH3->setVal(0.3);
  
  rrv_qcd_0lepLDP_ratioH1->setConstant(kTRUE);
  rrv_qcd_0lepLDP_ratioH2->setConstant(kTRUE);
  rrv_qcd_0lepLDP_ratioH3->setConstant(kTRUE);
  */
  
  printf("\n\n\n  ===== Doing a fit with SUSY component floating ====================\n\n") ;

  RooFitResult* fitResult = likelihood->fitTo( *rds, Save(true), PrintLevel(0) ) ;
  double logLikelihoodSusyFloat = fitResult->minNll() ;
  
  double logLikelihoodSusyFixed(0.) ;
  double testStatVal(-1.) ;
  if ( mu_susy_sig_val >= 0. ) {
    printf("\n\n\n  ===== Doing a fit with SUSY fixed ====================\n\n") ;
    printf(" fixing mu_susy_sig to %8.2f.\n", mu_susy_sig_val ) ;
    rrv_mu_susy_sig->setVal( mu_susy_sig_val ) ;
    rrv_mu_susy_sig->setConstant(kTRUE) ;
    
    fitResult = likelihood->fitTo( *rds, Save(true), PrintLevel(0) ) ;
    logLikelihoodSusyFixed = fitResult->minNll() ;
    testStatVal = 2.*(logLikelihoodSusyFixed - logLikelihoodSusyFloat) ;
    printf("\n\n\n ======= test statistic : -2 * ln (L_fixed / ln L_max) = %8.3f\n\n\n", testStatVal ) ;
  }


  return testStatVal ;

}