C++ RooDataHist类代码示例

void chi2(int xmin = 0, int xmax = 200, TString filename="../DsubMC/met2j0bIso2ewkScale_0_200.root", int nparam = 2){
RooAbsData::ErrorType errorType = RooAbsData::SumW2;

file = new TFile(filename);
RooRealVar* h = new RooRealVar("h","h",xmin,xmax); 

//Get Data
TH1D* hData = file->Get("dataih");
hData->SetName("hData");
//hData->Draw();
RooDataHist* data = new RooDataHist("data","data",*h,hData);

//Get Summed MC
TH1D* hMC = file->Get("hh");
hMC->SetName("hMC");
hMC->Draw();
RooDataHist rdhMC("MC","MC",*h,hMC);
RooHistPdf pdfMC("MCpdf","MCpdf",*h,rdhMC);

//make (plot) the curves
RooPlot* hFrame = h->frame(Name("hFrame"));
data->plotOn(hFrame,RooFit::DataError(errorType));
pdfMC.plotOn(hFrame,ProjWData(*data),Components(pdfMC),Name("h_total"));

//Determine Chi^2
double chi2fit = hFrame->chiSquare("h_total", "h_data", nparam);
cout<<"Chi 2 / dof: "<<chi2fit<<endl;

//Determine K-S
double ks = hMC->KolmogorovTest(hData);
cout<<"Kolmogorov-Smirnov: "<<ks<<endl;
}

void makeDataset( TString fname, TString tname, TString outfname ) {

  RooWorkspace *w = new RooWorkspace("w","w");

  w->factory( "Dst_M[1950.,2070.]" );
  w->factory( "D0_M[1810.,1920.]" );
  w->factory( "D0_LTIME_ps[0.00,5.]" );

  RooArgSet *observables = new RooArgSet();
  observables->add( *w->var("Dst_M") );
  observables->add( *w->var("D0_M") );
  observables->add( *w->var("D0_LTIME_ps") );
  w->defineSet("observables", *observables);

  w->var("Dst_M")->setBins(240);
  w->var("D0_M")->setBins(220);
  w->var("D0_LTIME_ps")->setBins(200);

  double Dst_M = -999.;
  double D0_M = -999.;
  double D0_LTIME_ps = -999.;

  TFile *tf = TFile::Open(fname);
  TTree *tree = (TTree*)tf->Get(tname);

  tree->SetBranchAddress( "Dst_M", &Dst_M );
  tree->SetBranchAddress( "D0_M" , &D0_M );
  tree->SetBranchAddress( "D0_LTIME_ps", &D0_LTIME_ps );

  RooDataSet *data = new RooDataSet("Data","Data",*observables);
  RooDataHist *dataH = new RooDataHist("DataHist","Data",*observables);

  for ( int ev=0; ev<tree->GetEntries(); ev++) {

    tree->GetEntry(ev);

    if ( ev%10000 == 0 ) cout << ev << " / " << tree->GetEntries() << endl;

    if ( Dst_M < w->var("Dst_M")->getMin() || Dst_M > w->var("Dst_M")->getMax() ) continue;
    if ( D0_M < w->var("D0_M")->getMin() || D0_M > w->var("D0_M")->getMax() ) continue;
    if ( D0_LTIME_ps < w->var("D0_LTIME_ps")->getMin() || D0_LTIME_ps > w->var("D0_LTIME_ps")->getMax() ) continue;

    w->var("Dst_M")->setVal(Dst_M);
    w->var("D0_M")->setVal(D0_M);
    w->var("D0_LTIME_ps")->setVal(D0_LTIME_ps);

    data->add( *observables );
    dataH->add( *observables );

  }
  tf->Close();

  w->import(*data);
  w->import(*dataH);

  w->writeToFile(outfname);

}

void rf903_numintcache(Int_t mode=0)
{
   // Mode = 0 : Run plain fit (slow)
   // Mode = 1 : Generate workspace with pre-calculated integral and store it on file (prepare for accelerated running)
   // Mode = 2 : Run fit from previously stored workspace including cached integrals (fast, requires run in mode=1 first)

   // C r e a t e ,   s a v e   o r   l o a d   w o r k s p a c e   w i t h   p . d . f .
   // -----------------------------------------------------------------------------------

   // Make/load workspace, exit here in mode 1
   RooWorkspace* w1 = getWorkspace(mode) ;
   if (mode==1) {

      // Show workspace that was created
      w1->Print() ;

      // Show plot of cached integral values
      RooDataHist* hhcache = (RooDataHist*) w1->expensiveObjectCache().getObj(1) ;
      if (hhcache) {

         new TCanvas("rf903_numintcache","rf903_numintcache",600,600) ;
         hhcache->createHistogram("a")->Draw() ;

      }
      else {
         Error("rf903_numintcache","Cached histogram is not existing in workspace");
      }
         return ;
   }

   // U s e   p . d . f .   f r o m   w o r k s p a c e   f o r   g e n e r a t i o n   a n d   f i t t i n g
   // -----------------------------------------------------------------------------------

   // This is always slow (need to find maximum function value empirically in 3D space)
   RooDataSet* d = w1->pdf("model")->generate(RooArgSet(*w1->var("x"),*w1->var("y"),*w1->var("z")),1000) ;

   // This is slow in mode 0, but fast in mode 1
   w1->pdf("model")->fitTo(*d,Verbose(kTRUE),Timer(kTRUE)) ;

   // Projection on x (always slow as 2D integral over Y,Z at fitted value of a is not cached)
   RooPlot* framex = w1->var("x")->frame(Title("Projection of 3D model on X")) ;
   d->plotOn(framex) ;
   w1->pdf("model")->plotOn(framex) ;

   // Draw x projection on canvas
   new TCanvas("rf903_numintcache","rf903_numintcache",600,600) ;
   framex->Draw() ;

   // Make workspace available on command line after macro finishes
   gDirectory->Add(w1) ;

   return ;

}

//_________________________________________________
void TestJeffreysGaussSigma(){
  // this one is VERY sensitive
  // if the Gaussian is narrow ~ range(x)/nbins(x) then the peak isn't resolved
  //   and you get really bizzare shapes
  // if the Gaussian is too wide range(x) ~ sigma then PDF gets renormalized
  //   and the PDF falls off too fast at high sigma
  RooWorkspace w("w");
  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,1,5])");
  w.factory("n[100,.1,2000]");
  w.factory("ExtendPdf::p(g,n)");
  //  w.var("sigma")->setConstant();
  w.var("mu")->setConstant();
  w.var("n")->setConstant();
  w.var("x")->setBins(301);

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;


  //  w.defineSet("poi","mu,sigma");
  //w.defineSet("poi","mu,sigma,n");
  w.defineSet("poi","sigma");
  w.defineSet("obs","x");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",3);

  const RooArgSet* temp = w.set("poi");
  pi.getParameters(*temp)->Print();
  //  return;

  //  return;
  RooGenericPdf* test = new RooGenericPdf("test","test","sqrt(2.)/sigma",*w.set("poi"));

  TCanvas* c1 = new TCanvas;
  RooPlot* plot = w.var("sigma")->frame();
  pi.plotOn(plot);
  test->plotOn(plot,LineColor(kRed),LineStyle(kDotted));
  plot->Draw();


}

void JeffreysPriorDemo(){
  RooWorkspace w("w");
  w.factory("Uniform::u(x[0,1])");
  w.factory("mu[100,1,200]");
  w.factory("ExtendPdf::p(u,mu)");

  //  w.factory("Poisson::pois(n[0,inf],mu)");

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());
  //  RooDataHist* asimov2 = w.pdf("pois")->generateBinned(*w.var("n"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;

  w.defineSet("poi","mu");
  w.defineSet("obs","x");
  //  w.defineSet("obs2","n");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  //  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",10);

  //  JeffreysPrior pi2("jeffreys2","jeffreys",*w.pdf("pois"),*w.set("poi"),*w.set("obs2"));

  //  return;
  RooGenericPdf* test = new RooGenericPdf("test","test","1./sqrt(mu)",*w.set("poi"));

  TCanvas* c1 = new TCanvas;
  RooPlot* plot = w.var("mu")->frame();
  //  pi.plotOn(plot, Normalization(1,RooAbsReal::Raw),Precision(.1));
  pi.plotOn(plot);
  //  pi2.plotOn(plot,LineColor(kGreen),LineStyle(kDotted));
  test->plotOn(plot,LineColor(kRed));
  plot->Draw();

}

void MCStudy()
{
    
    string fileName = "h4l_Template_2012_lowmass_unconstrained_new_v00.root";
    
    const char * wsName      = "combined";    //"w";
    const char * modelSBName = "ModelConfig"; //"modelConfig";
    const char * dataName    = dataname.c_str();
    
    TFile* file = TFile::Open(fileName);

    // get the workspace out of file
    RooWorkspace* w = (RooWorkspace*) file->Get(wsName);
    if (!wsName) {
        cout << "workspace not found" << endl;
        return -1.0;
    }
    
    // get the modelConfig out of the file
    ModelConfig* sbModel = (ModelConfig*) w->obj(modelSBName);
    
    // get the data out of file
    //RooAbsData* data = w->data(dataName);
    
    // get toy MC data
    TFile *fhist = new TFile("../ToyMC/toyMC_ggF125.root");
    TH1F *hsum = (TH1F*) fhist->Get("sumtoy1");
    
    RooRealVar m4l("m4l","m4l",120,130);
    RooDataHist* data = new RooDataHist("data","data",x,hsum);
    
    // get pdf
    RooAbsPdf* model = sbModel->GetPdf();
    
    RooPlot* xframe = m4l.frame();
    data->plotOn(xframe);
    model->fitTo(*data);
    model->plotOn(xframe,LineColor(kRed));
    
    TCanvas *c = new TCanvas("c","c");
    xframe->Draw();
    
}

//_________________________________________________
void TestJeffreysGaussMean(){
  RooWorkspace w("w");
  w.factory("Gaussian::g(x[0,-20,20],mu[0,-5,5],sigma[1,0,10])");
  w.factory("n[10,.1,200]");
  w.factory("ExtendPdf::p(g,n)");
  w.var("sigma")->setConstant();
  w.var("n")->setConstant();

  RooDataHist* asimov = w.pdf("p")->generateBinned(*w.var("x"),ExpectedData());

  RooFitResult* res = w.pdf("p")->fitTo(*asimov,Save(),SumW2Error(kTRUE));

  asimov->Print();
  res->Print();
  TMatrixDSym cov = res->covarianceMatrix();
  cout << "variance = " << (cov.Determinant()) << endl;
  cout << "stdev = " << sqrt(cov.Determinant()) << endl;
  cov.Invert();
  cout << "jeffreys = " << sqrt(cov.Determinant()) << endl;

  //  w.defineSet("poi","mu,sigma");
  w.defineSet("poi","mu");
  w.defineSet("obs","x");

  RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  //  pi.specialIntegratorConfig(kTRUE)->method1D().setLabel("RooAdaptiveGaussKronrodIntegrator1D")  ;
  //  pi.specialIntegratorConfig(kTRUE)->getConfigSection("RooIntegrator1D").setRealValue("maxSteps",3);

  const RooArgSet* temp = w.set("poi");
  pi.getParameters(*temp)->Print();

  //  return;
  RooGenericPdf* test = new RooGenericPdf("test","test","1",*w.set("poi"));

  TCanvas* c1 = new TCanvas;
  RooPlot* plot = w.var("mu")->frame();
  pi.plotOn(plot);
  test->plotOn(plot,LineColor(kRed),LineStyle(kDotted));
  plot->Draw();


}

void simul(){

	TFile* in = new TFile( "../bin/Kdata.root", "READ" );

	TH1D* h = in->Get("tof_8");

	RooRealVar * x = new RooRealVar( "x", "x", -10, 10 );

	RooDataHist * rdh = new RooDataHist( "data", "data", RooArgSet( *x ), h);

	RooRealVar m( "mean", "mean", 0, -1, 1 );
	RooRealVar s( "sigma", "sigma", .5, 0, 1 );
	RooGaussian g( "gauss", "gauss", )

	RooPlot * frame = x->frame();
	rdh->plotOn( frame );
	frame->Draw();
	

}

void MakeSpinPlots::DrawSpinSubTotBackground(TString mcName,bool signal){
  bool drawSM = (smName!="" && smName!=mcName);

  TCanvas cv;
  double thisN  = ws->data(mcName+"_Combined")->sumEntries();
  float norm = thisN;


  if(signal) norm = ws->var(Form("Data_%s_FULLFIT_Nsig",mcName.Data()))->getVal();
  RooPlot *frame = ws->var("cosT")->frame(0,1,10);

  RooDataSet* tmp = (RooDataSet*)ws->data(Form("Data_Combined"))->reduce("(mass>115 && mass<120) || (mass>130 && mass<135)");
  tmp->plotOn(frame,RooFit::Rescale(norm/tmp->sumEntries()));

  ws->pdf(Form("%s_FIT_cosTpdf",mcName.Data()))->plotOn(frame,RooFit::LineColor(kGreen),RooFit::Normalization(norm/tmp->sumEntries()));
  if(drawSM)  ws->pdf(Form("%s_FIT_cosTpdf",smName.Data()))->plotOn(frame,RooFit::LineColor(kRed),RooFit::Normalization(norm/tmp->sumEntries()));
  if(signal){
    RooDataHist *h = (RooDataHist*)ws->data( Form("Data_%s_Combined_bkgSub_cosT",mcName.Data()) );
    h->plotOn(frame,RooFit::MarkerStyle(4));
    std::cout << "Nsig: " << h->sumEntries() << std::endl;
  }

  
  frame->SetMaximum(frame->GetMaximum()*(signal?2.:1.2)*norm/tmp->sumEntries());
  frame->SetMinimum(-1*frame->GetMaximum());
  TLegend l(0.6,0.2,0.95,0.45);
  l.SetFillColor(0);
  l.SetBorderSize(0);
  l.SetHeader("Combined");
  l.AddEntry(frame->getObject(0),"Data m#in [115,120]#cup[130,135]","p");
  l.AddEntry(frame->getObject(1),mcName,"l");
  if(drawSM) l.AddEntry(frame->getObject(2),"SM Higgs","l");
  if(signal) l.AddEntry(frame->getObject(2+drawSM),"bkg-subtracted Data","p");
  
  frame->Draw();
  l.Draw("SAME");
  cv.SaveAs( basePath+Form("/cosThetaPlots/CosThetaDist_SimpleSub_%s%s_%s.png",outputTag.Data(),(signal ? "":"_BLIND"),mcName.Data()) );
  cv.SaveAs( basePath+Form("/cosThetaPlots/C/CosThetaDist_SimpleSub_%s%s_%s.C",outputTag.Data(),(signal ? "":"_BLIND"),mcName.Data()) );
  cv.SaveAs( basePath+Form("/cosThetaPlots/CosThetaDist_SimpleSub_%s%s_%s.pdf",outputTag.Data(),(signal ? "":"_BLIND"),mcName.Data()) );
}

RooFitResult* fitter_Zee(TH1D *hist){

  RooRealVar  Zmassvar("Zmassvar","Zmassvar", 82, 100);
  RooDataHist *datahist = new RooDataHist("data","Z Mass",Zmassvar,hist);
  RooPlot *Zmassvarframe = Zmassvar.frame(Name("Zmassvarframe"),Title(hist->GetTitle())) ;
  datahist->plotOn(Zmassvarframe);
  
  RooRealVar alpha  ("alpha"  ,        "alpha" , 0.005,0.001,0.1); 
  RooRealVar n      ("n"      ,            "n" , 1,0.001,10); 
  RooRealVar cbmean ("cbmean" ,       "cbmean" , 1, 0.8, 1.2);
  RooRealVar cbsigma("cbsigma",      "cbsigma" , 0.01, 0.001, 0.2);

  RooRealVar bwmean("bwmean","bwmean",91,85,95);
  RooRealVar bwsigma("bwsigma","bwsigma",3,2,4);
  RooRealVar expoconst("expoconst","expoconst",-0.1,-0.5,0);

  RooCBShape cball  ("cball"  , "crystal ball" , Zmassvar, cbmean, cbsigma, alpha, n);
  RooBreitWigner bw("bw","breit wigner",Zmassvar,bwmean,bwsigma);

  RooFFTConvPdf cballXbw("cballXbw","cball (X) bw",Zmassvar,bw,cball);
  RooExponential expo("expo", "exponential", Zmassvar, expoconst);

  RooRealVar frac("frac","frac",0.1,0.001,0.2);

  RooAddPdf Zshapemodel("Zshapemodel","expo + cball (X) bw",RooArgList(expo,cballXbw),frac);


  RooFitResult *fitres =Zshapemodel.fitTo(*datahist,Range(82,100),Save());   
  Zshapemodel.plotOn(Zmassvarframe,LineColor(kBlue));

  Zmassvarframe->Draw();

  fitres->Print();

  return fitres;

};

int scaleSmearFit(TString RDFile, TString MCFile, char BaseName[30])
{
  cout<<"Processing "<<BaseName<<endl;
  gStyle->SetPalette(1);
  //Data and histograms
  TFile *f_RD = new TFile(RDFile);
  TFile *f_MC = new TFile(MCFile);
  TH1D *h1_ZmassDaughEtaRD = (TH1D*)f_RD->Get("h1_ZmassDaughEtaRD")->Clone();
  TH1D *h1_ZmassDaughEtaMC[ScaleBins][ScaleBins];
  //Variables
  char histName[30];
  RooRealVar zMass("zMass","zMass",60,120);
  RooRealVar *nMC;
  RooRealVar nBRD("nBRD","nBRD",5,2000);
  RooRealVar nSRD("nSRD","nSRD",1000,0,20000);
  RooRealVar nBMC("nBMC","nBMC",5,2000);
  RooRealVar nSMC("nSMC","nSMC",1000,0,20000);
  RooPlot *zmassframe;
  RooPlot *zmassframeMC;
  CPlot *plotFit;
  CPlot *plotFitMC;
  //zMass.setBins(50);
  RooFitResult* fitRes;
  RooFitResult* fitResMC;
  TCanvas *myCan = MakeCanvas("myCan","myCan",800,600);
  RooDataHist *ZmassMC;
  //RooHistPdf *pdfMC;
  RooAddPdf *pdfMC;
  RooAddPdf* pdfRD;
  //RD fitting
  RooDataHist *ZmassRD =
    new RooDataHist("ZmassRD","ZmassRD", RooArgSet(zMass),h1_ZmassDaughEtaRD);
  //CBreitWignerConvCrystalBall ZsignalPdf("ZsigPdf",zMass);
  //pdfRD=new RooAddPdf("pdfRD","pdfRD",RooArgList(*(ZsignalPdf.model)),RooArgList(nRD));
  CVoigtian     ZsigRD("ZsigRD",zMass);
  CErfExpo ZbgRD("ZbgRD",zMass);
  //CQuadraticExp ZbgRD("ZbgRD",zMass);
  pdfRD=new RooAddPdf("pdfRD","pdfRD",RooArgList(*(ZsigRD.model),*(ZbgRD.model)),RooArgList(nSRD,nBRD));
  pdfRD=new RooAddPdf("pdfRD","pdfRD",RooArgList(*(ZsigRD.model),*(ZbgRD.model)),RooArgList(nSRD,nBRD));
  fitRes=pdfRD->fitTo(*ZmassRD,Extended(),Minos(kTRUE),Save(kTRUE));
  //fitRes=ZsignalPdf.model->fitTo(*ZmassRD,Minos(kTRUE),Save(kTRUE));
  zmassframe=zMass.frame(Bins(60));
  ZmassRD->plotOn(zmassframe,MarkerStyle(kFullCircle),MarkerSize(0.9),DrawOption("zp"));
  //ZsignalPdf.model->plotOn(zmassframe,LineColor(kBlue),DrawOption("l"));
  pdfRD->plotOn(zmassframe,LineColor(kBlack),DrawOption("l"));
  pdfRD->plotOn(zmassframe,Components(RooArgSet(*(ZsigRD.model))),LineColor(kBlue),DrawOption("l"));
  pdfRD->plotOn(zmassframe,Components(RooArgSet(*(ZbgRD.model))),LineColor(kRed),DrawOption("l"));
  sprintf(histName,"ZmassRD_%s",BaseName);
  plotFit = new CPlot(histName,zmassframe,"","","Z mass");
  plotFit->setOutDir("Plot");
  plotFit->Draw(myCan,kTRUE,"png");

  //CErfExpo *pdfZbg;
  //double nLL[41][41];

  //                                 90 0.004
  double ScaleWidth = (ScaleH-ScaleL)/(ScaleBins-1);
  double SmearWidth = (SmearH-SmearL)/(ScaleBins-1);
  sprintf(histName,"h2_NLL_%s",BaseName);
  TH2D *h2_NLL = new TH2D(histName,histName,
      ScaleBins,ScaleL-ScaleWidth/2,ScaleH+ScaleWidth/2,
      ScaleBins,SmearL-SmearWidth/2,SmearH+SmearWidth/2);
  //TH2D *h2_NLL = new TH2D("h2_NLL","NLL",41,0.97,1.05,41,0.5,1.5);
  double x,prob;
  //RooAbsReal *nll;
  double nll;
  double binContent;
  //***
  CVoigtian ZsigMC("ZsigMC",zMass);
  RooArgSet allArgset(zMass);
  RooArgSet anaVar;
  //for(int i(0);i<=0;i++)for(int j(0);j<=0;j++)
  for(int i(0);i<=ScaleBins-1;i++)for(int j(0);j<=ScaleBins-1;j++)
  {
    sprintf(histName,"h1_ZmassDaughEta_%d_%d",i,j);
    h1_ZmassDaughEtaMC[i][j] = (TH1D*)f_MC->Get(histName)->Clone(histName);
    ZmassMC =
      new RooDataHist("ZmassMC","ZmassMC",RooArgSet(zMass),h1_ZmassDaughEtaMC[i][j]);
    //                                               interpolation order
    //pdfMC =new RooHistPdf ("pdfMC", "pdfMC", zMass,*ZmassMC, 1);
    //Using fitting MC
    nSMC.setVal(285);
    ZsigMC.mass->setVal(91.37);
    ZsigMC.sigma->setVal(0.42);
    ZsigMC.width->setVal(4.4);
    pdfMC = new RooAddPdf("pdfMC","pdfMC",RooArgList(*(ZsigMC.model)),RooArgList(nSMC));
    fitResMC = pdfMC->fitTo(*ZmassMC,Extended(),Minos(kTRUE),Save(kTRUE),SumW2Error(kTRUE));//SumW2Error(kTRUE) default

    nll=0;
    //nll= h1_ZmassDaughEtaMC[i][j]->Chi2Test(h1_ZmassDaughEtaRD,"CHI2/NDF");
    // code 1000 dataHist sum
    int intCode = pdfMC->getAnalyticalIntegral(allArgset,anaVar);
    double norm = pdfMC->analyticalIntegral(intCode);
    cout<<"norm: code "<<norm<<" : "<<intCode<<"======================"<<endl;
    //double totalProb(0);
    //****
    for(int k(1);k<=60;k++)
    {
      x=h1_ZmassDaughEtaMC[i][j]->GetBinCenter(k);
    //  binContent = h1_ZmassDaughEtaRD->GetBinContent(k);
//.........这里部分代码省略.........

void Back_2dFit(void){
  TFile *ifile = TFile::Open("/home/vitaly/B0toDh0/TMVA/fil_b2dh_gen.root");
  TTree *tree = (TTree*)ifile->Get("TEvent");

  RooCategory b0f("b0f","b0f");
  b0f.defineType("comb",-1);

  RooArgSet argset;
  
  const double mbcMin = 5.20;
  const double mbcMax = 5.29;
  const double deMin = -0.15;
  const double deMax = 0.3;

  RooRealVar mbc("mbc","M_{bc}",mbcMin,mbcMax,"GeV"); argset.add(mbc);
  mbc.setRange("Signal",mbc_min,mbc_max);
  mbc.setRange("mbcSignal",mbc_min,mbc_max);
  mbc.setRange("deSignal",mbcMin,mbcMax);
  RooRealVar de("de","#DeltaE",deMin,deMax,"GeV"); argset.add(de);
  de.setRange("Signal",de_min,de_max);
  de.setRange("mbcSignal",deMin,deMax);
  de.setRange("deSignal",de_min,de_max);
  RooRealVar md("md","md",DMass-md_cut,DMass+md_cut,"GeV"); argset.add(md);
  RooRealVar mk("mk","mk",KMass-mk_cut,KMass+mk_cut,"GeV"); argset.add(mk);
  RooRealVar mpi0("mpi0","mpi0",Pi0Mass-mpi0_cut,Pi0Mass+mpi0_cut,"GeV"); argset.add(mpi0);
  RooRealVar bdtgs("bdtgs","bdtgs",bdtgs_cut,1.); argset.add(bdtgs);
  RooRealVar atckpi_max("atckpi_max","atckpi_max",0.,atckpi_cut); argset.add(atckpi_max);

  argset.add(b0f);

  RooDataSet ds("ds","ds",tree,argset);
  RooDataSet* ds0 = ds.reduce(RooArgSet(de,mbc));
  RooDataHist* dh = ds0->binnedClone();
  ds.Print();
  dh->Print();

  ////////////
  // de pdf //
  ////////////
//  RooRealVar c1("c1","c1",mc_c1,-1.,0.); if(cComb) c1.setConstant(kTRUE);
//  RooRealVar c2("c2","c2",mc_c2,0.,0.1); if(cComb) c2.setConstant(kTRUE);
  RooRealVar c3("c3","c3",mc_c1,-1.,0.); if(cComb) c3.setConstant(kTRUE);
  RooRealVar c4("c4","c4",mc_c2,0.,0.1); if(cComb) c4.setConstant(kTRUE);
  RooChebychev pdf_de("pdf_de","pdf_de",de,RooArgSet(c3,c4));

  /////////////
  // mbc pdf //
  /////////////
  RooRealVar argpar("argpar","argus shape parameter",mc_argpar,-100.,-1.); if(cComb)
  argpar.setConstant(kTRUE);
  RooRealVar argedge("argedge","argedge",mc_argedge,5.285,5.292); if(cComb) argedge.setConstant(kTRUE);
  RooArgusBG pdf_mbc("pdf_mbc","Argus PDF",mbc,argedge,argpar);

  /////////
  // pdf //
  /////////
  RooProdPdf pdf("pdf","pdf",RooArgList(pdf_de,pdf_mbc));

  pdf.fitTo(*dh,Verbose(),Timer(true));

  /////////////
  //  Plots  //
  /////////////
  // de //
  RooPlot* deFrame = de.frame();
  ds.plotOn(deFrame,DataError(RooAbsData::SumW2),MarkerSize(1),MarkerColor(kGreen));
  pdf.plotOn(deFrame,LineWidth(2),LineColor(kGreen));
  ds.plotOn(deFrame,DataError(RooAbsData::SumW2),MarkerSize(1),CutRange("mbcSignal"));
  pdf.plotOn(deFrame,LineWidth(2),ProjectionRange("mbcSignal"));

  RooHist* hdepull = deFrame->pullHist();
  RooPlot* dePull = de.frame(Title("#Delta E pull distribution"));
  dePull->addPlotable(hdepull,"P");
  dePull->GetYaxis()->SetRangeUser(-5,5);

  TCanvas* cm = new TCanvas("#Delta E, Signal","#Delta E, Signal",600,700);
  cm->cd();

  TPad *pad3 = new TPad("pad3","pad3",0.01,0.20,0.99,0.99);
  TPad *pad4 = new TPad("pad4","pad4",0.01,0.01,0.99,0.20);
  pad3->Draw();
  pad4->Draw();

  pad3->cd();
  pad3->SetLeftMargin(0.15);
  pad3->SetFillColor(0);

  deFrame->GetXaxis()->SetTitleSize(0.05);
  deFrame->GetXaxis()->SetTitleOffset(0.85);
  deFrame->GetXaxis()->SetLabelSize(0.04);
  deFrame->GetYaxis()->SetTitleOffset(1.6);
  deFrame->Draw();

  stringstream out;
  out.str("");
  out << "#chi^{2}/n.d.f = " << deFrame->chiSquare();
  TPaveText *pt = new TPaveText(0.6,0.8,0.98,0.9,"brNDC");
  pt->SetFillColor(0);
  pt->SetTextAlign(12);
  pt->AddText(out.str().c_str());
//.........这里部分代码省略.........

void makeSystPlot( TFile * f, TString oldFolder, RooWorkspace *WS,  string channel, string syst, int toMassNo, int fromMassNo, int rightBinNo, int addRightBin,  int addRightBinm1) //massNo 0-51, see xSec7TeV.h 
{


  std::cout << "oldFolder,  channel , addRightBin, addRightBinm1: " << oldFolder << " , " <<channel << " , " << addRightBin  << " , "<< addRightBinm1 << std::endl;

  RooArgList  * hobs = new RooArgList("hobs");
  //  RooRealVar BDT("BDT", "BDT", -1, 1);///OLD VARIABLE NAME HERE
  RooRealVar BDT("CMS_vhbb_BDT_Wln", "CMS_vhbb_BDT_Wln", -1, 1);///OLD VARIABLE NAME HERE
  hobs->add(*WS->var("CMS_vhbb_BDT_Wln"));  ///NEW VARIABLE NAME HERE
  RooWorkspace *tempWS =  (RooWorkspace*) f->Get(oldFolder.Data());
  TString systT(syst);
  TString chanT(channel);

  if((kount < 3) && (channel=="data_obs"))
  {
    kount++;
    std::string namen  = channel;
    
    std::cout << oldFolder.Data() << std::endl;
    std::cout << namen << std::endl;
    RooDataHist* tempRooDataHistNom = (RooDataHist*)  tempWS->data(namen.c_str());
    TH1 *tempHistNom = tempRooDataHistNom->createHistogram(namen.c_str(),BDT,RooFit::Binning(bins));

   tempHistNom->Rebin(rebin);

 if(addRightBin == 1)
    {
     float  err0 = tempHistNom->GetBinError(rightBinNo);
     float  con0 = tempHistNom->GetBinContent(rightBinNo);
 
     float  err1 = tempHistNom->GetBinError(rightBinNo-1);
     float  con1 = tempHistNom->GetBinContent(rightBinNo-1);
       
     tempHistNom->SetBinContent(rightBinNo,0);
     tempHistNom->SetBinError(rightBinNo,0);
     tempHistNom->SetBinContent(rightBinNo-1,con0+con1);
     tempHistNom->SetBinError(rightBinNo-1,sqrt(err0*err0+err1*err1));
    }


 if(addRightBinm1 == 1)
    {
     float  err0 = tempHistNom->GetBinError(rightBinNo-1);
     float  con0 = tempHistNom->GetBinContent(rightBinNo-1);
 
     float  err1 = tempHistNom->GetBinError(rightBinNo-2);
     float  con1 = tempHistNom->GetBinContent(rightBinNo-2);
       
     tempHistNom->SetBinContent(rightBinNo-1,0);
     tempHistNom->SetBinError(rightBinNo-1,0);
     tempHistNom->SetBinContent(rightBinNo-2,con0+con1);
     tempHistNom->SetBinError(rightBinNo-2,sqrt(err0*err0+err1*err1));
    }





    RooDataHist *DHnom = new RooDataHist(channel.c_str(),"",*hobs,tempHistNom);  
    WS->import(*(new RooHistPdf(channel.c_str(),"",*hobs,*DHnom)));
 
 }

 if (channel!="data_obs")
{
  std::string nameUp; 
  std::string namen; 
  std::string nameDown;


  if(syst == "stat")
  {
   if(oldFolder.Contains("Wenu"))
   { 
     nameUp  = channel + "_CMS_vhbb_stat" + channel + "_WenuUp";
     namen  = channel;
     nameDown  = channel + "_CMS_vhbb_stat" + channel + "_WenuDown";

     if(channel == "s_Top")
     {
     nameUp  = channel + "_CMS_vhbb_statsTop_WenuUp";
     namen  = channel;
     nameDown  = channel + "_CMS_vhbb_statsTop_WenuDown";
     }

   }
   else
   {
     nameUp  = channel + "_CMS_vhbb_stat" + channel + "_WmunuUp";
     namen  = channel;
     nameDown  = channel + "_CMS_vhbb_stat" + channel + "_WmunuDown";

     if(channel == "s_Top")
     {
     nameUp  = channel + "_CMS_vhbb_statsTop_WmunuUp";
     namen  = channel;
     nameDown  = channel + "_CMS_vhbb_statsTop_WmunuDown";
     }

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

void FitBias(TString CAT,TString CUT,float SIG,float BKG,int NTOYS)
{
  gROOT->ForceStyle();
  
  RooMsgService::instance().setSilentMode(kTRUE);
  RooMsgService::instance().setStreamStatus(0,kFALSE);
  RooMsgService::instance().setStreamStatus(1,kFALSE);
  
  // -----------------------------------------
  TFile *fTemplates = TFile::Open("templates_"+CUT+"_"+CAT+"_workspace.root");
  RooWorkspace *wTemplates = (RooWorkspace*)fTemplates->Get("w");
  RooRealVar *x            = (RooRealVar*)wTemplates->var("mTop");
  RooAbsPdf *pdf_signal    = (RooAbsPdf*)wTemplates->pdf("ttbar_pdf_Nominal");
  RooAbsPdf *pdf_bkg       = (RooAbsPdf*)wTemplates->pdf("qcdCor_pdf"); 
  TRandom *rnd = new TRandom();
  rnd->SetSeed(0);
  x->setBins(250);   
  RooPlot *frame;

  TFile *outf;

  if (NTOYS > 1) { 
    outf = TFile::Open("FitBiasToys_"+CUT+"_"+CAT+".root","RECREATE");
  }

  float nSigInj,nBkgInj,nSigFit,nBkgFit,eSigFit,eBkgFit,nll;

  TTree *tr = new TTree("toys","toys");
  
  tr->Branch("nSigInj",&nSigInj,"nSigInj/F");
  tr->Branch("nSigFit",&nSigFit,"nSigFit/F");
  tr->Branch("nBkgInj",&nBkgInj,"nBkgInj/F");
  tr->Branch("nBkgFit",&nBkgFit,"nBkgFit/F");
  tr->Branch("eSigFit",&eSigFit,"eSigFit/F");
  tr->Branch("eBkgFit",&eBkgFit,"eBkgFit/F");
  tr->Branch("nll"    ,&nll    ,"nll/F");

  for(int itoy=0;itoy<NTOYS;itoy++) {
    // generate pseudodataset
    nSigInj = rnd->Poisson(SIG);
    nBkgInj = rnd->Poisson(BKG);
    RooRealVar *nSig = new RooRealVar("nSig","nSig",nSigInj);
    RooRealVar *nBkg = new RooRealVar("nBkg","nBkg",nBkgInj);
    RooAddPdf *model = new RooAddPdf("model","model",RooArgList(*pdf_signal,*pdf_bkg),RooArgList(*nSig,*nBkg)); 
    RooDataSet *data = model->generate(*x,nSigInj+nBkgInj);
    
    RooDataHist *roohist = new RooDataHist("roohist","roohist",RooArgList(*x),*data);
    // build fit model
    RooRealVar *nFitSig = new RooRealVar("nFitSig","nFitSig",SIG,0,10*SIG);
    RooRealVar *nFitBkg = new RooRealVar("nFitBkg","nFitBkg",BKG,0,10*BKG);
    RooAddPdf *modelFit = new RooAddPdf("modelFit","modelFit",RooArgList(*pdf_signal,*pdf_bkg),RooArgList(*nFitSig,*nFitBkg)); 
    // fit the pseudo dataset
    RooFitResult *res = modelFit->fitTo(*roohist,RooFit::Save(),RooFit::Extended(kTRUE));
    //res->Print();
    nSigFit = nFitSig->getVal();
    nBkgFit = nFitBkg->getVal();
    eSigFit = nFitSig->getError();
    eBkgFit = nFitBkg->getError();
    nll     = res->minNll();
    tr->Fill();
    if (itoy % 100 == 0) {
      cout<<"Toy #"<<itoy<<": injected = "<<nSigInj<<", fitted = "<<nSigFit<<", error = "<<eSigFit<<endl;
    }
    if (NTOYS == 1) {
      frame = x->frame();
      roohist->plotOn(frame); 
      model->plotOn(frame);
    }
  }
  if (NTOYS == 1) {
    TCanvas *can = new TCanvas("Toy","Toy",900,600);
    frame->Draw();
  }  
  else {
    outf->cd();
    tr->Write();
    outf->Close();
    fTemplates->Close();
  }  
}

void TwoCategorySimZFitter( TH1& h_TT, TH1& h_TF, 
                            double intLumi, int removeEE )
{
  // The fit variable - lepton invariant mass
  rooMass_ = new RooRealVar("Mass","m_{ee}",60.0, 120.0, "GeV/c^{2}");
  rooMass_->setBins(20.0);
  RooRealVar Mass = *rooMass_;

  // Make the category variable that defines the two fits,
  // namely whether the probe passes or fails the eff criteria.
  RooCategory sample("sample","") ;
  sample.defineType("TT", 1) ;
  sample.defineType("TF", 2) ; 

  gROOT->cd();

  ///////// convert Histograms into RooDataHists
  RooDataHist* data_TT = new RooDataHist("data_TT","data_TT",
					  RooArgList(Mass), &h_TT);
  RooDataHist* data_TF = new RooDataHist("data_TF","data_TF",
					  RooArgList(Mass), &h_TF);

  RooDataHist* data = new RooDataHist( "fitData","fitData",
				       RooArgList(Mass),Index(sample),
				       Import("TT",*data_TT), Import("TF",*data_TF) ); 

  data->get()->Print();
  cout << "Made datahist" << endl;


  // ********** Construct signal & bkg shape PDF ********** //
  makeSignalPdf();
  cout << "Made signal pdf" << endl;
  makeBkgPdf();
  cout << "Made bkg pdf" << endl;

  // Now supply integrated luminosity in inverse picobarn
  // -->  we get this number from the CMS lumi group
  // https://twiki.cern.ch/twiki/bin/view/CMS/LumiWiki2010Data
  RooRealVar lumi("lumi","lumi", intLumi);


  // Now define Z production cross section variable (in pb) 
  RooRealVar xsec("xsec","xsec", 1300., 400.0, 2000.0);


  // Define efficiency variables  
  RooRealVar eff("eff","eff", 0.9, 0.5, 1.0);

  // Now define acceptance variables --> we get these numbers from MC   
//   RooRealVar acc_BB("acc_BB","acc_BB", 0.2253);
//   RooRealVar acc_EB("acc_EB","acc_EB", 0.1625);
//   RooRealVar acc_EE("acc_EE","acc_EE", 0.0479);

  double ACCEPTANCE=0.4357;
  if(removeEE==1) ACCEPTANCE= 0.3878;
  if(removeEE==2) ACCEPTANCE= 0.2253;
  RooRealVar acc("acc","acc", ACCEPTANCE);

  // Define background yield variables: they are not related to each other  
  RooRealVar nBkgTF("nBkgTF","nBkgTF",10.,-10.,100000.);
 
  ////////////////////////////////////////////////////////////////////////////
  ////////////////////////////////////////////////////////////////////////////
 //  Define signal yield variables.  
  // They are linked together by the total cross section:  e.g. 
  //          Nbb = sigma*L*Abb*effB

  char* formula;
  RooArgList* args;
  formula = "lumi*xsec*acc*eff*eff";
  args = new RooArgList(lumi,xsec,acc,eff);
  RooFormulaVar nSigTT("nSigTT", formula, *args);
  delete args;

  formula = "lumi*xsec*acc*eff*(1.0-eff)";
  args = new RooArgList(lumi,xsec,acc,eff);
  RooFormulaVar nSigTF("nSigTF",formula, *args);
  delete args;

  /////////////////////////////////////////////////////////////////////////
  /////////////////////////////////////////////////////////////////////////

   RooArgList componentsTF(*signalShapePdfTF_,*bkgShapePdfTF_);
   RooArgList yieldsTF(nSigTF, nBkgTF );	  

   RooExtendPdf pdfTT("pdfTT","extended sum pdf", *signalShapePdfTT_, nSigTT);
   RooAddPdf pdfTF("pdfTF","extended sum pdf",componentsTF, yieldsTF);


   // The total simultaneous fit ...
   RooSimultaneous totalPdf("totalPdf","totalPdf", sample);
   totalPdf.addPdf(pdfTT,"TT");
   totalPdf.Print();
   totalPdf.addPdf(pdfTF,"TF");
   totalPdf.Print();


  // ********* Do the Actual Fit ********** //  
   RooFitResult *fitResult = totalPdf.fitTo(*data, Save(true), 
//.........这里部分代码省略.........