C++ RooRealVar::getError方法代码示例

void PEs(RooAbsPdf *iGen,RooAbsPdf *iFit,int iN,int iNEvents,RooRealVar &iVar,RooRealVar &iSig,RooRealVar &iMean,
	 RooRealVar &iScale,RooRealVar &iRes) { 
  double iM0 = iMean.getVal(); double iS0 = iSig.getVal();
  //iScale.setVal(iMeanScale); iRes.setVal(iSigScale);
  TRandom1 *lRand = new TRandom1(0xDEADBEEF);
  TNtuple * lDN= new TNtuple( "xxx","xxx","ntot:m_r:m:merr:sig_r:sig:sigerr");
  for(int i0=0;i0<iN;i0++){
    if(i0 % 10 == 0) cout << "+++++++++++++++++++++++++++ running ======> " << i0 << endl;
    int lN    = lRand->Poisson(iNEvents);
    RooDataSet * lSignal  = iGen->generate(iVar,lN);
    iMean.setVal(iM0); iSig.setVal(iS0);
    iFit->fitTo(*lSignal,Strategy(1));//,Save(kTRUE),PrintLevel(1));
    if(iMean.getError() < 0.05) iFit->fitTo(*lSignal,Strategy(2));
    Float_t values[]={
      (Float_t) lN,
      (Float_t) 90.78/iScale.getVal(),
      (Float_t) iMean.getVal(),
      (Float_t) iMean.getError(),
      (Float_t) iSig.getVal(),
      (Float_t) iSig.getVal(),
      (Float_t) iSig.getError()
    };
    lDN->Fill(values);
  }
  TFile *lF = new TFile("XXX.root","RECREATE");
  lDN->Write();
  lF->Close();
}

void Plot(RooAbsPdf *iGen,RooAbsPdf *iFit,int iN,int iNEvents,RooRealVar &iVar,RooRealVar &iSig,RooRealVar &iMean,
	  RooRealVar &iScale,RooRealVar &iRes,RooDataSet *iData=0) { 
  TRandom1 *lRand = new TRandom1(0xDEADBEEF);
  RooDataSet * lSignal  = iGen->generate(iVar,iNEvents);
  if(iData == 0) {
    iFit->fitTo(*lSignal,Strategy(1));
    if(iMean.getError() < 0.05) iFit->fitTo(*lSignal,Strategy(2));
  } else {
    iFit->fitTo(*iData,Strategy(1));
    if(iMean.getError() < 0.05) iFit->fitTo(*iData,Strategy(2));
  }    
  iVar.setBins(30);
  RooPlot *lFrame1 = iVar.frame(RooFit::Title("XXX")) ;
  if(iData == 0) lSignal->plotOn(lFrame1);
  if(iData != 0) iData->plotOn(lFrame1);
  iFit->plotOn(lFrame1);
  TCanvas *iC =new TCanvas("A","A",800,600);
  iC->cd(); lFrame1->Draw();
  iC->SaveAs("Crap.png");
  if(iData != 0) { 
    RooPlot *lFrame2 = iVar.frame(RooFit::Title("XXX")) ;
    iData->plotOn(lFrame2);
    iGen->plotOn(lFrame2);
    TCanvas *iC1 =new TCanvas("B","B",800,600);
    iC1->cd(); lFrame2->Draw();
    iC1->SaveAs("Crap.png");
  }
}

void FitterUtils::initiateParams(int nGenSignalZeroGamma, int nGenSignalOneGamma, int nGenSignalTwoGamma, RooRealVar const& expoConstGen, RooRealVar& nSignal, RooRealVar& nPartReco, 
      RooRealVar& nComb, RooRealVar& fracZero, RooRealVar& fracOne, RooRealVar& expoConst, RooRealVar&  nJpsiLeak, bool constPartReco, RooRealVar const& fracPartRecoSigma)
{
   TRandom rand;
   rand.SetSeed();

   int nGenSignal = nGenSignalZeroGamma + nGenSignalOneGamma + nGenSignalTwoGamma;

   double nGenSignal2;
   double nGenPartReco2;
   if(!constPartReco)
   {
      nGenSignal2 = rand.Uniform(nGenSignal-5*sqrt(nGenSignal), nGenSignal+5*sqrt(nGenSignal));
      nGenPartReco2 = rand.Uniform(nGenPartReco-5*sqrt(nGenPartReco), nGenPartReco+5*sqrt(nGenPartReco));
   }
   if(constPartReco)
   { 
      double nGenSigPartReco( nGenSignal+nGenPartReco );
      double nGenSigPartReco2( rand.Uniform( nGenSigPartReco-5*sqrt(nGenSigPartReco), nGenSigPartReco+5*sqrt(nGenSigPartReco) ) );
      double fracPartReco1( nGenPartReco/(1.*nGenSignal));
      double fracPartReco2( rand.Uniform(fracPartReco1-5*fracPartRecoSigma.getVal(), fracPartReco1+5*fracPartRecoSigma.getVal()) ); 

      nGenPartReco2 = fracPartReco2*nGenSigPartReco2 / (1+fracPartReco2); 
      nGenSignal2 = nGenSigPartReco2 / (1+fracPartReco2); 
   }
   double nGenComb2 = rand.Uniform(nGenComb-5*sqrt(nGenComb), nGenComb+5*sqrt(nGenComb));
   double nGenJpsiLeak2 = rand.Uniform(nGenJpsiLeak-5*sqrt(nGenJpsiLeak), nGenJpsiLeak+5*sqrt(nGenJpsiLeak));


   nSignal.setVal(nGenSignal2);
   nSignal.setRange(TMath::Max(0.,nGenSignal2-10.*sqrt(nGenSignal)) , nGenSignal2+10*sqrt(nGenSignal));

   nPartReco.setVal(nGenPartReco2);
   nPartReco.setRange(TMath::Max(0.,nGenPartReco2-10.*sqrt(nGenPartReco)), nGenPartReco2+10*sqrt(nGenPartReco));


   nComb.setVal(nGenComb2);
   nComb.setRange(TMath::Max(0.,nGenComb2-10.*sqrt(nGenComb)), nGenComb2+10*sqrt(nGenComb));

   nJpsiLeak.setVal(nGenJpsiLeak2);
   nJpsiLeak.setRange(TMath::Max(0., nGenJpsiLeak2-10*sqrt(nGenJpsiLeak)), nGenJpsiLeak2+10*sqrt(nGenJpsiLeak));

   double fracGenZero(nGenSignalZeroGamma/(1.*nGenSignal));
   double fracGenOne(nGenSignalOneGamma/(1.*nGenSignal));

   fracZero.setVal(rand.Gaus(fracGenZero, sqrt(nGenSignalZeroGamma)/(1.*nGenSignal))) ;
   fracZero.setRange(0., 1.);
   fracOne.setVal(rand.Gaus(fracGenOne, sqrt(nGenSignalOneGamma)/(1.*nGenSignal))) ;
   fracOne.setRange(0., 1.);

   expoConst.setVal(rand.Uniform( expoConstGen.getVal() - 5*expoConstGen.getError(), expoConstGen.getVal() + 5*expoConstGen.getError() ) );
   expoConst.setRange( expoConstGen.getVal() - 10*expoConstGen.getError(), expoConstGen.getVal() + 10*expoConstGen.getError() );
}

MakeBiasStudy::MakeBiasStudy() {
  int Nmodels = 8;
  //RooRealVar* mass = ws->var("mass");
  RooRealVar *mass = new RooRealVar("mass","mass", 100,180);
  RooRealVar *nBkgTruth = new RooRealVar("TruthNBkg","", 0,1e9);
  //  RooAbsData* realData = ws->data("Data_Combined")->reduce( Form("evtcat==evtcat::%s",cat.Data()) );

  double Bias[Nmodels][Nmodels];
  double BiasE[Nmodels][Nmodels];
  MakeAICFits MakeAIC_Fits;
  for(int truthType = 0; truthType < Nmodels; truthType++){

    RooAbsPdf *truthPdf = MakeAIC_Fits.getBackgroundPdf(truthType,mass);
    RooExtendPdf *truthExtendedPdf = new RooExtendPdf("truthExtendedPdf","",*truthPdf,*nBkgTruth);
    //truthExtendedPdf.fitTo(*realData,RooFit::Strategy(0),RooFit::NumCPU(NUM_CPU),RooFit::Minos(kFALSE),RooFit::Extended(kTRUE));
    //truthExtendedPdf.fitTo(*realData,RooFit::Strategy(2),RooFit::NumCPU(NUM_CPU),RooFit::Minos(kFALSE),RooFit::Extended(kTRUE));
    
    double BiasWindow = 2.00;
    mass->setRange("biasRegion", mh-BiasWindow, mh+BiasWindow);
    double TruthFrac = truthExtendedPdf->createIntegral(mass,RooFit::Range("biasRegion"),RooFit::NormSet(*mass))->getVal();
    double NTruth = TruthFrac * nBkgTruth->getVal();
    double NTruthE = TruthFrac * nBkgTruth->getError();
    
    RooDataSet* truthbkg = truthPdf->generate(RooArgSet(*mass),nBkgTruth);

    for(int modelType = 0; modelType < Nmodels; modelType++){
      RooAbsPdf* ModelShape = MakeAIC_Fits.getBackgroundPdf(modelType,mass);
      RooRealVar *nBkgFit = new RooRealVar("FitNBkg", "", 0, 1e9);
      RooExtendPdf ModelExtendedPdf = new RooExtendPdf("ModelExtendedPdf", "",*ModelShape, *nBkgFit);
      ModelExtendedPdf.fitTo(truthbkg, RooFit::Strategy(0),RooFit::NumCPU(NUM_CPU),RooFit::Minos(kFALSE),RooFit::Extended(kTRUE));
      ModelExtendedPdf.fitTo(truthbkg, RooFit::Strategy(2),RooFit::NumCPU(NUM_CPU),RooFit::Minos(kFALSE),RooFit::Extended(kTRUE));

      double FitFrac = ModelExtendedPdf.createIntegral(mass,RooFit::Range("biasRegion"),RooFit::NormSet(mass))->getVal();
      double NFit = FitFrac * nBkgFit->getVal();
      double NFitE = FitFrac * nBkgFit->getError();

      Bias[truthType][modelType] = fabs(NFit - NTruth);
      BiasE[truthType][modelType] = fabs(NFitE - NTruthE);
    }

    
  }
  
  for(int i = 0; i < Nmodels; i++) { 
    std::cout <<  "===== Truth Model : " << MakeBiasStudy::Category(i) << " ===== " << std::endl; 
    for (int j = 0; j < Nmodels; j++) { 
      std::cout << "Fit Model: " << MakeBiasStudy::Category(j) << "  , Bias = " << Bias[i][j] << " +/- " << BiasE[i][j] << std::endl; 
    }   
  }
  

}

double getSigmaError(RooWorkspace *w) 
{
  using namespace RooFit;
  // Access saved Sigma CB values and error from WorkSpace
  RooRealVar *mySavedSigma = w->var("sigma_{CB}");
  RooRealVar *mySavedDeltaM = w->var("#Deltam_{CB}");
  Double_t sigma = mySavedSigma->getVal();
  Double_t sigmaerror = mySavedSigma->getError();
  Double_t DeltaM = mySavedDeltaM->getVal();
  Double_t DeltaMerror = mySavedDeltaM->getError();
  Double_t Sigpercent = 100 * sigma / (91.187 + DeltaM);
  Double_t Sigpercenterror = 100 * pow(pow(sigmaerror / (91.198 + DeltaM), 2) + pow((sigma * DeltaMerror) / pow(DeltaM + 91.198, 2), 2), .5);
  return Sigpercenterror;
} // plotFitOn(..)

void simplePrintResults() {
    vector<string> filenames;
    filenames.push_back("Output/Test2/result/FIT_DATA_Psi2SJpsi_PPPrompt_Bkg_SecondOrderChebychev_pt65300_rap016_cent0200_262620_263757.root");

    const char* parname = "N_Jpsi_PP";

    vector<string>::iterator it = filenames.begin();
    for (it; it<filenames.end(); it++) {
        TFile *f = new TFile(it->c_str());
        if (!f) {
            cout << "Error, " << *it << " not found" << endl;
            continue;
        }
        RooWorkspace *ws = (RooWorkspace*) f->Get("workspace");
        if (!ws) {
            cout << "Error, workspace not found in " << *it << endl;
            continue;
        }

        RooRealVar *var = ws->var(parname);
        if (!ws) {
            cout << "Error, variable " << parname << " not found in " << *it << endl;
            continue;
        }
        cout << *it << " " << var->getVal() << " +- " << var->getError() << endl;
    }
}

void PDF_GLWADS_Dpi_K3pi::setUncertainties(config c)
{
    switch(c)
    {
    case lumi1fb:
    {
        obsErrSource = "1fb-1, ExpNll/sept2012K3PIResult.root";
        TString File = this->dir+"/ExpNll/sept2012K3PIResult.root";
        TFile *fr = TFile::Open(File);
        RooFitResult *r = (RooFitResult*)fr->Get("fitresult_model_reducedData_binned");
        assert(r);
        for ( int i=0; i<nObs; i++ )
        {
            RooRealVar* pObs = (RooRealVar*)((RooArgList*)observables)->at(i);
            RooRealVar* pRes = (RooRealVar*)r->floatParsFinal().find(obsTmkToMalcolm(pObs->GetName()));
            assert(pRes);
            StatErr[i] = pRes->getError();
        }
        SystErr[0] = 0.010;   // afav_dpi_obs
        SystErr[1] = 0.00011; // rp_dpi_obs
        SystErr[2] = 0.00011; // rm_dpi_obs
        fr->Close();
        delete r;
        delete fr;
        break;
    }
    default:
        cout << "PDF_GLWADS_Dpi_K3pi::setUncertainties() : ERROR : config "+ConfigToTString(c)+" not found." << endl;
        exit(1);
    }
}

void checkBestFitPoint(std::string workspace, std::string fitFile, bool splusb){
	
	// Open the ws file...
	TFile *fd_=0;
	TFile *fw_=0;
		
	gSystem->Load("$CMSSW_BASE/lib/$SCRAM_ARCH/libHiggsAnalysisCombinedLimit.so");
	gROOT->SetBatch(true);
	gStyle->SetOptFit(0);
	gStyle->SetOptStat(0);
	gStyle->SetPalette(1,0);

	fw_ =  TFile::Open(workspace.c_str());
	w   = (RooWorkspace*) fw_->Get("w");
  w->Print();
	RooDataSet *data = (RooDataSet*) w->data("data_obs");
  if (splusb) {
    mc_s = (RooStats::ModelConfig*)w->genobj("ModelConfig");
  } else {
    mc_s = (RooStats::ModelConfig*)w->genobj("ModelConfig_bonly");
  }
	std::cout << "make nll"<<std::endl;
	nll = mc_s->GetPdf()->createNLL(
		*data,RooFit::Constrain(*mc_s->GetNuisanceParameters())
		,RooFit::Extended(mc_s->GetPdf()->canBeExtended()));
	
	// Now get the best fit result
	fd_ =  TFile::Open(fitFile.c_str());
	RooFitResult *fit;
  if (splusb) {
    fit =(RooFitResult*)fd_->Get("fit_s");
  } else {
    fit =(RooFitResult*)fd_->Get("fit_b");
  }
	RooArgSet fitargs = fit->floatParsFinal();
	
	std::cout << "Got the best fit values" <<std::endl;		
	w->saveSnapshot("bestfitall",fitargs,true);
	
  TString filename;
  if (splusb) {
    filename = "minimum_s.pdf";
  } else {
    filename = "minimum_b.pdf";
  }
	// Now make the plots!	
	TCanvas *c = new TCanvas("c","",600,600);
	c->SaveAs((filename+"["));

	TIterator* iter(fitargs->createIterator());
        for (TObject *a = iter->Next(); a != 0; a = iter->Next()) {
                 RooRealVar *rrv = dynamic_cast<RooRealVar *>(a);      
                 std::string name = rrv->GetName();
		 TGraph *gr = graphLH(name,rrv->getError());
		 gr->Draw("ALP");
		 c->SaveAs((filename+"["));
	}
	c->SaveAs((filename+"]"));
}

// grab the initial parameters and errors for making pull distributions:
// Take these from a fit file to the data themselves 
void fillInitialParams(RooArgSet *args, std::map<std::string, std::pair<double,double> > &vals){
	
	 TIterator* iter(args->createIterator());
         for (TObject *a = iter->Next(); a != 0; a = iter->Next()) {
                 RooRealVar *rrv = dynamic_cast<RooRealVar *>(a);      
                 std::string name = rrv->GetName();
		 std::pair<double,double> valE(rrv->getVal(),rrv->getError());
		 vals.insert( std::pair<std::string,std::pair<double ,double> > (name,valE)) ;
	 }
	
}

void FitterUtilsSimultaneousExpOfPolyTimesX::initiateParams(int nGenSignalZeroGamma, int nGenSignalOneGamma, int nGenSignalTwoGamma,
      RooRealVar& nKemu, RooRealVar& nSignal, RooRealVar& nPartReco,
      RooRealVar& nComb, RooRealVar& fracZero, RooRealVar& fracOne,
      RooRealVar&  nJpsiLeak, bool constPartReco, RooRealVar const& fracPartRecoSigma,
      RooRealVar& l1Kee, RooRealVar& l2Kee, RooRealVar& l3Kee, RooRealVar& l4Kee, RooRealVar& l5Kee,
      RooRealVar& l1Kemu, RooRealVar& l2Kemu, RooRealVar& l3Kemu, RooRealVar& l4Kemu, RooRealVar& l5Kemu,
      RooRealVar const& l1KeeGen, RooRealVar const& l2KeeGen, RooRealVar const& l3KeeGen, RooRealVar const& l4KeeGen, RooRealVar const& l5KeeGen 

      )
{
  FitterUtilsExpOfPolyTimesX::initiateParams(nGenSignalZeroGamma, nGenSignalOneGamma, nGenSignalTwoGamma,
           nSignal, nPartReco, nComb, fracZero, fracOne, nJpsiLeak, constPartReco, fracPartRecoSigma,
           l1Kee, l2Kee, l3Kee, l4Kee, l5Kee,
           l1KeeGen, l2KeeGen, l3KeeGen, l4KeeGen, l5KeeGen ); 



  TRandom rand;
  rand.SetSeed();

  nKemu.setVal(rand.Uniform(nGenKemu-5*sqrt(nGenKemu), nGenKemu+5*sqrt(nGenKemu)));
  nKemu.setRange(nGenKemu-10*sqrt(nGenKemu), nGenKemu+10*sqrt(nGenKemu));

  l1Kemu.setVal(rand.Uniform( l1KeeGen.getVal() - 5*l1KeeGen.getError(), l1KeeGen.getVal() + 5*l1KeeGen.getError() ) );
  l1Kemu.setRange( l1KeeGen.getVal() - 10*l1KeeGen.getError(), l1KeeGen.getVal() + 10*l1KeeGen.getError() );

  l2Kemu.setVal(rand.Uniform( l2KeeGen.getVal() - 5*l2KeeGen.getError(), l2KeeGen.getVal() + 5*l2KeeGen.getError() ) );
  l2Kemu.setRange( l2KeeGen.getVal() - 10*l2KeeGen.getError(), l2KeeGen.getVal() + 10*l2KeeGen.getError() );

  l3Kemu.setVal(rand.Uniform( l3KeeGen.getVal() - 5*l3KeeGen.getError(), l3KeeGen.getVal() + 5*l3KeeGen.getError() ) );
  l3Kemu.setRange( l3KeeGen.getVal() - 10*l3KeeGen.getError(), l3KeeGen.getVal() + 10*l3KeeGen.getError() );

  l4Kemu.setVal(rand.Uniform( l4KeeGen.getVal() - 5*l4KeeGen.getError(), l4KeeGen.getVal() + 5*l4KeeGen.getError() ) );
  l4Kemu.setRange( l4KeeGen.getVal() - 10*l4KeeGen.getError(), l4KeeGen.getVal() + 10*l4KeeGen.getError() );

  l5Kemu.setVal(rand.Uniform( l5KeeGen.getVal() - 5*l5KeeGen.getError(), l5KeeGen.getVal() + 5*l5KeeGen.getError() ) );
  l5Kemu.setRange( l5KeeGen.getVal() - 10*l5KeeGen.getError(), l5KeeGen.getVal() + 10*l5KeeGen.getError() );

}

void MakeSpinPlots::MakeChannelComp(const char* mcType){
  TGraphErrors graph(catNames.size());

  RooRealVar mu("mu","",-50,50);

  float exp = ws->data(Form("%s_Combined",mcType))->sumEntries();///total * 607*lumi/12.;

  TH1F frame("frame","",catNames.size(),0,catNames.size());

  //graph.GetXaxis()->SetNdivisions(catNames.size());

  float min=99999,max=-99999;
  for(int i=0;i<catNames.size();i++){
    RooRealVar *ind = ws->var( Form("Data_%s_INDFIT_%s_Nsig",mcType, catNames.at(i).Data()) );
    RooRealVar *f = ws->var( Form("Data_%s_FULLFIT_%s_fsig",mcType, catNames.at(i).Data()) );
    float mu = ind->getVal()/exp/f->getVal();
    float muE = ind->getError()/exp/f->getVal();
    graph.SetPoint(i,i+0.5,mu);
    graph.SetPointError(i,0,muE);

    if(mu-muE < min) min = mu-muE;
    if(mu+muE > max) max = mu+muE;
    
    //graph.GetXaxis()->SetBinLabel(i+1,catNames.at(i));
  }

  TF1 fit("fit","[0]",0+0.5,catNames.size()+0.5);
  graph.Fit(&fit,"MNE");
  float mean  = fit.GetParameter(0);
  float meanE = fit.GetParError(0);
  
  frame.SetAxisRange(min-1.5,max+0.5,"Y");
  frame.SetYTitle("Fitted #sigma/#sigma_{SM}");
  frame.SetXTitle("Category");
  TCanvas cv;
  frame.Draw();

  TBox err(0,mean-meanE,catNames.size(),mean+meanE);
  err.SetFillColor(kGreen);

  frame.Draw();
  err.Draw("SAME");
  TLine mLine(0,mean,catNames.size(),mean);
  mLine.Draw("SAME");
  graph.Draw("PSAME");
  

  cv.SaveAs(basePath+Form("/ChannelComp_%s_%s.png",mcType,outputTag.Data()));
  cv.SaveAs(basePath+Form("/ChannelComp_%s_%s.pdf",mcType,outputTag.Data()));
}

void MakeSpinPlots::printYields(const char* mcType){
  RooRealVar * tot = ws->var(Form("Data_%s_FULLFIT_Nsig",mcType));
  if(tot==0) return;

  cout << "Total Yield:  " << tot->getVal() << "  +-  " << tot->getError() <<endl;
  cout << "Category Yields: CONSTRAINED FIT " << endl;
  for(int i=0;i<catNames.size();i++){
    RooRealVar *f = ws->var( Form("Data_%s_FULLFIT_%s_fsig",mcType, catNames.at(i).Data()) );
    cout << "\t" << catNames.at(i) <<":   " << tot->getVal()*f->getVal() << "  +-  " << tot->getError()*f->getVal() <<endl;
  }
  cout << "\nCategory Yields: INDEPENDENT FIT " << endl;
  for(int i=0;i<catNames.size();i++){
    RooRealVar *ind = ws->var( Form("Data_%s_INDFIT_%s_Nsig",mcType, catNames.at(i).Data()) );
    if(ind==0) continue;
    cout << "\t" << catNames.at(i) <<":   " << ind->getVal() << "  +-  " << ind->getError() <<endl;
  }

  float exp = ws->data(Form("%s_Combined",mcType))->sumEntries();///total * 607*lumi/12.;
  cout << endl << "Expected Events:  "  << exp << endl;
  cout << "Expected Yields Per Category: " <<endl;
  for(int i=0;i<catNames.size();i++){ 
    RooRealVar *f = ws->var( Form("Data_%s_FULLFIT_%s_fsig",mcType, catNames.at(i).Data()) );
    cout << "\t" << catNames.at(i) <<":   " << exp*f->getVal() <<endl;
  }

  cout << "mu:  " << tot->getVal()/exp << "  +-  "
       << tot->getError()/exp <<endl;
  
  for(int i=0;i<catNames.size();i++){
    RooRealVar *ind = ws->var( Form("Data_%s_INDFIT_%s_Nsig",mcType, catNames.at(i).Data()) );
    if(ind==0) continue;
    RooRealVar *f = ws->var( Form("Data_%s_FULLFIT_%s_fsig",mcType, catNames.at(i).Data()) );
    cout << "\t" << catNames.at(i) <<":   " << ind->getVal()/(exp*f->getVal()) << "  +-  " << ind->getError()/(exp*f->getVal()) <<endl;
  }
  MakeChannelComp(mcType);
}

void logStatisticsPar(std::ostream& out, RooDataSet *dataSet, RooRealVar *realVar, int nBins, double chi2, const RooArgList &variables)
{
    TH1 *histogram = dataSet->createHistogram(Form("h%s", dataSet->GetName()), *realVar, RooFit::Binning(nBins));
    
    // Create the TeX file
    out << "\\documentclass[10pt]{article}" << std::endl;
    out << "\\usepackage[usenames]{color} %used for font color" << std::endl;
    out << "\\usepackage{fontspec}" << std::endl;
    out << "\\usepackage{xunicode}" << std::endl;
    out << "\\usepackage{xltxtra}" << std::endl;
    out << "\\defaultfontfeatures{Scale=MatchLowercase}" << std::endl;
    out << "\\setromanfont[Mapping=tex-text]{Myriad Pro}" << std::endl;
    out << "\\setsansfont[Mapping=tex-text]{Myriad Pro}" << std::endl;
    out << "\\setmonofont{Monaco}" << std::endl;
    out << "\\begin{document}" << std::endl;
    out << "\\thispagestyle{empty}" << std::endl;
    out << "\\setlength{\\tabcolsep}{1ex}" << std::endl;
    out << "\\setlength{\\fboxsep}{0ex}" << std::endl;
    out << "{\\fontsize{7pt}{0.9em}\\selectfont" << std::endl;
    out << "\\framebox{\\begin{tabular*}{60pt}{[email protected]{\\extracolsep{\\fill}}r}" << std::endl;
    
    // This is the particular info for the histogram
    out << "Entries & " ;
    formatNumber(histogram->GetEntries(), out) << " \\\\" << std::endl;
    out << "Mean & " ;
    formatNumber(histogram->GetMean(), out) << " \\\\" << std::endl;
    out << "RMS & " ;
    formatNumber(histogram->GetRMS(), out) << " \\\\" << std::endl;
    if (chi2 > 0.0) {
        out << "Fit $\\chi^{2}$ & " ;
        formatNumber(chi2, out) << " \\\\" << std::endl;
    }
    RooRealVar *theVariable;
    for (int index = 0; index < variables.getSize(); index++) {
        theVariable = dynamic_cast<RooRealVar*>(variables.find(variables[index].GetName()));
        out << theVariable->GetTitle() << " & $\\textrm{" ;
        formatNumber(theVariable->getValV(), out) << "} \\pm \\textrm{" ;
        formatNumber(theVariable->getError(), out) << "}$ \\\\" << std::endl;
    }
    out << "\\end{tabular*}}}" << std::endl;
    out << "\\end{document}" << std::endl;
    histogram->Delete();
}

bool extractParameter(string fileName, const char* parName, pair<double,double>& value) 
{
  TFile *f = new TFile( fileName.c_str() );
  if (!f) {
    cout << "[Error] " << fileName << " not found" << endl; return false;
  }
  RooWorkspace *ws = (RooWorkspace*) f->Get("workspace");
  if (!ws) {
    cout << "[ERROR] Workspace not found in " << fileName << endl; return false;
  }    
  RooRealVar *var = ws->var(parName);
  if (!var) {
    value = make_pair( -999.9 , -999.9 ); return true;
  }
  
  value = make_pair( var->getValV() , var->getError() );
        
  return true;
}

RooDataSet * getDataAndFrac(TString name, TString q2name, TreeReader * mydata, TCut cut, RooRealVar * MM, double * frac, Str2VarMap jpsiPars, double *outnsig)
{
	RooRealVar * cosThetaL = new RooRealVar("cosThetaL","cosThetaL",0.,-1.,1.);
	RooRealVar * cosThetaB = new RooRealVar("cosThetaB","cosThetaB",0.,-1.,1.);
	TCut massCut = "Lb_MassConsLambda > 5590 && Lb_MassConsLambda < 5650";

	Analysis * ana = new Analysis(name+"_mass"+q2name,"Lb",mydata,&cut,MM);
	ana->AddVariable("J_psi_1S_MM");
	ana->AddVariable(cosThetaL);
	ana->AddVariable(cosThetaB);
	RooAbsPdf * mysig = stringToPdf("Gauss","sig",MM,jpsiPars);
	RooAbsPdf * mybkg = stringToPdf("Exp","bkgM",MM);
	RooRealVar * mynsig = new RooRealVar("mynsig","mynsig",50,0,100000);
	RooRealVar * mynbkg = new RooRealVar("mynbkg","mynbkg",10,0,100000);
	RooAbsPdf * Mmodel = new RooAddPdf("MassModel","MassModel",RooArgSet(*mysig,*mybkg),RooArgSet(*mynsig,*mynbkg));
	ana->applyCuts(&cut);
	RooDataSet * data = ana->GetDataSet("-recalc");
	Mmodel->fitTo(*data,Extended(kTRUE));
	
	double sigBkg = mybkg->createIntegral(*MM,NormSet(*MM),Range("Signal"))->getVal();
	double sig = mysig->createIntegral(*MM,NormSet(*MM),Range("Signal"))->getVal();
	double nsig = mynsig->getVal();
	double nbkg = mynbkg->getVal();
	if(frac)
	{
		frac[0] = nsig*sig/(nsig*sig+nbkg*sigBkg);
		frac[1] = frac[0]*TMath::Sqrt( TMath::Power(mynsig->getError()/nsig,2) + TMath::Power(mynbkg->getError()/nbkg,2) );
	}
	TCut mycut = cut + massCut;
	ana->applyCuts(&mycut);

	TCanvas * cc = new TCanvas();
	GetFrame(MM,Mmodel,data,"-nochi2-plotAllComp",30,NULL,0,"cos#theta_{#Lambda}")->Draw();
	cc->Print("M_"+name+"_"+q2name+".pdf");
	if(*outnsig) *outnsig = nsig;
	return ana->GetDataSet("-recalc");
}