reallySimplified.C

double reallySimplified(std::string modelName="signal",std::string outname="ht400Output"){
    RooArgList xlist_;
    RooArgList muA_;
    RooCategory sampleType("bin_number","Bin Number");
    RooRealVar  observation("observed","Observed Events bin",1);
    double bkg[1] = {9.0};
    double sig[1] = {3.0};
    int nbins = 1;//

    for (int b=1;b<=nbins;b++){
        sampleType.defineType(Form("%d",b-1),b-1);
        sampleType.setIndex(b-1);
    }
    RooArgSet   obsargset(observation,sampleType);
    for (int b=1;b<=nbins;b++){
        double bkgy = (double)bkg[b-1];
        RooRealVar *meanAsimov_ = new RooRealVar(Form("exp_bin_%d_In_asimov",b),Form("ASimov expected bin %d",b),bkgy);
        meanAsimov_->setConstant(true);
        RooRealVar *x_ = new RooRealVar(Form("exp_bin_%d",b),Form("bkg bin %d",b),bkgy,0.01*bkgy,bkgy*100);
        x_->setMin(0.001);
        muA_.add(*meanAsimov_);
        xlist_.add(*x_);
    }

    // Make the signal component
    RooRealVar r("r","r",0,-5,10);
    RooArgList signals_;
    for (int b=1;b<=nbins;b++) {
        RooFormulaVar *sigF = new RooFormulaVar(Form("signal_%d",b),Form("@0*%g",sig[b-1]),RooArgSet(r));
        signals_.add(*sigF);
    }

    RooArgList plist_;
    RooArgList slist_;
    sampleType.setIndex(0);
    RooSimultaneous combined_pdf("combined_pdf","combined_pdf",sampleType);
    for (int b=1;b<=nbins;b++){
        RooAddition *sum = new RooAddition(Form("splusb_bin_%d",b),Form("Signal plus background in bin %d",b),RooArgList(*((RooRealVar*)(signals_.at(b-1))),*((RooRealVar*)(xlist_.at(b-1)))));
        RooPoisson  *pois = new RooPoisson(Form("pdf_bin_%d",b),Form("Poisson in bin %d",b),observation,(*sum));
        //RooGaussian *pois = new RooGaussian(Form("pdf_bin_%d",b),Form("Poisson in bin %d",b),observation,(*sum),RooFit::RooConst(TMath::Sqrt(sum->getVal())));
        combined_pdf.addPdf(*pois,Form("%d",b-1));
        slist_.add(*sum);
        plist_.add(*pois);
    }
    RooDataSet asimovdata("AsimovData","Asimov in all Bins",obsargset);
    r.setVal(0);
    for (int b=1;b<=nbins;b++){
        sampleType.setIndex(b-1);
        double exp = (double) ((*(RooRealVar*)slist_.at(b-1)).getVal());
        ((RooRealVar*)muA_.at(b-1))->setVal(exp);
        observation.setVal(int(exp));
        // hDataAsimov->SetBinContent(b,exp);
        asimovdata.add(RooArgSet(observation,sampleType));
        std::cout << "asimov data " << observation.getVal() << " background " <<  bkg[b-1] << std::endl;
        // std::cout << " post fit Exp background At " << b << ", Simple code=" << exp << ", combine code=" << bkgcombfit->GetBinContent(b) << " Observed in the data " << data.GetBinContent(b) << " signal " << signal->GetBinContent(b) << std::endl;
        // std::cout << " Asi = "<< observation.getVal() << ", besty = " << ((RooRealVar*)muA_.at(b-1))->getVal()  << " before the fit that was " << bkg->GetBinContent(b) << std::endl;
    }
    //r.setVal(5);
    r.setConstant(false);
    RooAbsReal *nllA_ = combined_pdf.createNLL(asimovdata);//,RooFit::ExternalConstraints(RooArgList(asimov_constraint_pdf)));
    RooMinimizer mc(*nllA_);
    mc.minimize("Minuit2","minimize");
    for (int b=1;b<=nbins;b++){
        double exp = (double) ((*(RooRealVar*)slist_.at(b-1)).getVal());
        std::cout << "post fit " << exp << std::endl;
    }
    std::cout << r.getVal() << std::endl;
    return 0.;
}