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

//____________________________________
void DoSPlot(RooWorkspace* ws){
  std::cout << "Calculate sWeights" << std::endl;

  // get what we need out of the workspace to do the fit
  RooAbsPdf* model = ws->pdf("model");
  RooRealVar* zYield = ws->var("zYield");
  RooRealVar* qcdYield = ws->var("qcdYield");
  RooDataSet* data = (RooDataSet*) ws->data("data");

  // fit the model to the data.
  model->fitTo(*data, Extended() );

  // The sPlot technique requires that we fix the parameters
  // of the model that are not yields after doing the fit.
  RooRealVar* sigmaZ = ws->var("sigmaZ");
  RooRealVar* qcdMassDecayConst = ws->var("qcdMassDecayConst");  
  sigmaZ->setConstant();
  qcdMassDecayConst->setConstant();


  RooMsgService::instance().setSilentMode(true);


  // Now we use the SPlot class to add SWeights to our data set
  // based on our model and our yield variables
  RooStats::SPlot* sData = new RooStats::SPlot("sData","An SPlot",
		            *data, model, RooArgList(*zYield,*qcdYield) );


  // Check that our weights have the desired properties

  std::cout << "Check SWeights:" << std::endl;


  std::cout << std::endl <<  "Yield of Z is " 
	    << zYield->getVal() << ".  From sWeights it is "
	    << sData->GetYieldFromSWeight("zYield") << std::endl;


  std::cout << "Yield of QCD is " 
	    << qcdYield->getVal() << ".  From sWeights it is "
	    << sData->GetYieldFromSWeight("qcdYield") << std::endl
	    << std::endl;

  for(Int_t i=0; i < 10; i++)
    {
      std::cout << "z Weight   " << sData->GetSWeight(i,"zYield") 
		<< "   qcd Weight   " << sData->GetSWeight(i,"qcdYield") 
		<< "  Total Weight   " << sData->GetSumOfEventSWeight(i) 
		<< std::endl;
    }

  std::cout << std::endl;

  // import this new dataset with sWeights
 std::cout << "import new dataset with sWeights" << std::endl;
 ws->import(*data, Rename("dataWithSWeights"));


}

//____________________________________
void MakePlots(RooWorkspace* wks) {

  // Make plots of the data and the best fit model in two cases:
  // first the signal+background case
  // second the background-only case.

  // get some things out of workspace
  RooAbsPdf* model = wks->pdf("model");
  RooAbsPdf* sigModel = wks->pdf("sigModel");
  RooAbsPdf* zjjModel = wks->pdf("zjjModel");
  RooAbsPdf* qcdModel = wks->pdf("qcdModel");

  RooRealVar* mu = wks->var("mu");
  RooRealVar* invMass = wks->var("invMass");
  RooAbsData* data = wks->data("data");


  //////////////////////////////////////////////////////////
  // Make plots for the Alternate hypothesis, eg. let mu float

  mu->setConstant(kFALSE);

  model->fitTo(*data,Save(kTRUE),Minos(kFALSE), Hesse(kFALSE),PrintLevel(-1));

  //plot sig candidates, full model, and individual componenets
  new TCanvas();
  RooPlot* frame = invMass->frame() ;
  data->plotOn(frame ) ;
  model->plotOn(frame) ;
  model->plotOn(frame,Components(*sigModel),LineStyle(kDashed), LineColor(kRed)) ;
  model->plotOn(frame,Components(*zjjModel),LineStyle(kDashed),LineColor(kBlack)) ;
  model->plotOn(frame,Components(*qcdModel),LineStyle(kDashed),LineColor(kGreen)) ;

  frame->SetTitle("An example fit to the signal + background model");
  frame->Draw() ;
  //  cdata->SaveAs("alternateFit.gif");

  //////////////////////////////////////////////////////////
  // Do Fit to the Null hypothesis.  Eg. fix mu=0

  mu->setVal(0); // set signal fraction to 0
  mu->setConstant(kTRUE); // set constant

  model->fitTo(*data, Save(kTRUE), Minos(kFALSE), Hesse(kFALSE),PrintLevel(-1));

  // plot signal candidates with background model and components
  new TCanvas();
  RooPlot* xframe2 = invMass->frame() ;
  data->plotOn(xframe2, DataError(RooAbsData::SumW2)) ;
  model->plotOn(xframe2) ;
  model->plotOn(xframe2, Components(*zjjModel),LineStyle(kDashed),LineColor(kBlack)) ;
  model->plotOn(xframe2, Components(*qcdModel),LineStyle(kDashed),LineColor(kGreen)) ;

  xframe2->SetTitle("An example fit to the background-only model");
  xframe2->Draw() ;
  //  cbkgonly->SaveAs("nullFit.gif");

}

RooAbsPdf *MakeModelNoSignal(RooDataHist *data, RooRealVar *mww, char *name) {
    char *modelName = (char*)calloc(50, sizeof(char));
    strcat(modelName, name);
    strcat(modelName, "model");
    char *signalName = (char*)calloc(50, sizeof(char));
    char *backgroundName = (char*)calloc(50, sizeof(char));
    strcat(signalName, name);
    strcat(signalName, "signal");
    strcat(backgroundName, name);
    strcat(backgroundName, "background");

    char *meanName = (char*)calloc(50, sizeof(char));
    strcat(meanName, name);
    strcat(meanName, "mean");

    char *sigmaName = (char*)calloc(50, sizeof(char));
    strcat(sigmaName, name);
    strcat(sigmaName, "sigma");

    char *decayName = (char*)calloc(50, sizeof(char));
    strcat(decayName, name);
    strcat(decayName, "decay");

    char *weightName = (char*)calloc(50, sizeof(char));
    strcat(weightName, name);
    strcat(weightName, "weight");

    RooRealVar *decay = new RooRealVar(decayName, "decay", -10, 10, "GeV");
    RooExponential *background = new RooExponential(backgroundName, "background", *mww, *decay);

    background->fitTo(*data);
    decay->setConstant();
    return background;
}

void testResolution() { 
  Prep();
  RooRealVar    lXVar  ("XVar","mass(GeV/c^{2})",100,60,150); lXVar.setBins(1000);
  RooRealVar l1Sigma("sigma1","sigma1",1.6 ,0.,15.);  //l1Sigma.setConstant(kTRUE);
  RooRealVar l2Sigma("sigma2","sigma2",1.6,0.,15.);   l2Sigma.setConstant(kTRUE);
  RooRealVar l3Sigma("sigma3","sigma3",2.90,0.,35.);  l3Sigma.setConstant(kTRUE);
  RooRealVar lN     ("n"     ,"n"     ,1.00,-15,15.); lN.setConstant(kTRUE);
  RooRealVar lExp   ("exp"   ,"exp"   ,-0.003,-15,15.); //lExp.setConstant(kTRUE);
  RooRealVar lR0Mean("xmean","xmean",0,-10,10);    lR0Mean.setConstant(kTRUE);
  RooRealVar lR1Mean("mean","mean",90.8,60,150);   //lR1Mean.setConstant(kTRUE);
  RooVoigtianShape     lGAdd("Add","Add",lXVar,lR1Mean,l1Sigma,l2Sigma,lN,l3Sigma,true);
  
  RooRealVar    lSPar  ("SPar","SPar", 1.,0., 2.);
  RooFormulaVar lXShift("uparshift","@0*@1",RooArgList(lXVar,lSPar));
  TH1F *lMass = getMass(0,-5,5,-1.5,1.5);
  RooDataHist *lMHist = new RooDataHist("M" ,"M" ,RooArgSet(lXVar),lMass);
  RooHistPdf  *lMPdf  = new RooHistPdf ("MH","MH",lXShift,lXVar,*lMHist,5); 
  RooRealVar lGSigma("gsigma","gsigma",1.6 ,0.,15.);  
  RooGaussian   lGaus1("gaus1","gaus1",lXVar,lR0Mean,lGSigma);
  RooFFTConvPdf  lConv("Conv","Conv",lXVar,*lMPdf,lGaus1);

  RooDataSet *lData = new RooDataSet("crap","crap",RooArgSet(lXVar)); 
  fillData(lData,lXVar,0,-5,5,-1.5,1.5);
  lConv.fitTo(*lData,Strategy(2));
  lGSigma.setVal(lGSigma.getVal());
  lSPar.setVal(lSPar.getVal()*1.01);
  //cout << "=====> Check " << l1Sigma.getVal() << " --- " << lR1Mean.getVal() << "----" << lSPar.getVal() << endl;
  PEs(&lConv,&lGAdd,2000,500,lXVar,l1Sigma,lR1Mean,lSPar,lGSigma);
  lData = 0;
  Plot(&lConv,&lGAdd,2000,50000,lXVar,l1Sigma,lR1Mean,lSPar,lGSigma,lData);
}

RooRealVar* LoadParameters(string filename, string label, string parname )
{

  RooRealVar *newVar = 0;

  // now read in the parameters from the file.  Stored in the file as                                                                                                  
  // parameter manager   |   name   | initial val  | min  | max  | step                                                                                                
  fstream parameter_file(filename.c_str(), ios::in);
  if (! parameter_file) {
    cout << "Error:  Couldn't open parameters file " << filename << endl;
    return false;
  }
  
  string name;
  string category;
  double initial_val, min, max, step;
  
  char c;
  Bool_t foundPar = kFALSE;
  while (true) {
    
    // skip white spaces and look for a #                                                                                                                              
    while(parameter_file.get(c)) {
      if      (isspace(c) || c == '\n')
        continue;
      else if (c == '#')
        while (c != '\n') parameter_file.get(c);
      else {
        parameter_file.putback(c);
        break;
      }
    }
    if (parameter_file.fail())
      break;
    
    parameter_file >> category >> name >> initial_val >> min >> max >> step;
    if (parameter_file.fail())
      break;
    
    if (category == label && parname == name) {
      // create a new fit parameter         
      newVar = new RooRealVar(name.c_str(),name.c_str(),initial_val,min, max);
      if (step == 0) {
        newVar->setConstant(kTRUE);  
      }
      
      break;
    }
  } //end while loop

  if (!newVar) {
    cout << "Could not load parameter " << parname << " from file " << filename << " , category " << label << endl;
    assert(newVar);
  }

  return newVar;
}

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

            doofit::roofit::pdfs::DooCubicSplinePdf splinePdf("splinePdf",epdf->Var("obsEtaOS"),knots,*list,0,0.5);
            //doofit::roofit::pdfs::DooCubicSplinePdf* splinePdf = new doofit::roofit::pdfs::DooCubicSplinePdf("splinePdf", epdf->Var("obsEtaOS"), knots, *list,0,0.5);



     //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->Var("obsEtaOS"),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->Var("obsEtaOS"),epdf->Var("asym_prod"),epdf->Var("asym_det"),epdf->Var("asym_tageff"));
            DecRateCoeff *coeff_sh = new DecRateCoeff("coef_sinh","coef_sinh",DecRateCoeff::CPEven,epdf->Cat("catFinalState"),epdf->Cat("obsTag"),epdf->Var("f1_f"),epdf->Var("f1_fbar"),epdf->Var("obsEtaOS"),splinePdf,epdf->Var("tageff"),epdf->Var("obsEtaOS"),epdf->Var("asym_prod"),epdf->Var("asym_det"),epdf->Var("asym_tageff"));
            DecRateCoeff *coeff_ch = new DecRateCoeff("coef_cosh","coef_cosh",DecRateCoeff::CPEven,epdf->Cat("catFinalState"),epdf->Cat("obsTag"),epdf->Var("f0_f"),epdf->Var("f0_fbar"),epdf->Var("obsEtaOS"),splinePdf,epdf->Var("tageff"),epdf->Var("obsEtaOS"),epdf->Var("asym_prod"),epdf->Var("asym_det"),epdf->Var("asym_tageff"));

            epdf->AddRealToStore(coeff_ch);
            epdf->AddRealToStore(coeff_sh);
            epdf->AddRealToStore(coeff_c);
            epdf->AddRealToStore(coeff_s);


    ///////////////////Generiere PDF's/////////////////////
    //Zeit
    epdf->GaussModel("resTimeGauss",epdf->Var("obsTime"),epdf->Var("allTimeResMean"),epdf->Var("allTimeReso"));
    epdf->BDecay("pdfSigTime",epdf->Var("obsTime"),epdf->Var("tau"),epdf->Var("dgamma"),epdf->Real("coef_cosh"),epdf->Real("coef_sinh"),epdf->Real("coef_cos"),epdf->Real("coef_sin"),epdf->Var("deltaM"),epdf->Model("resTimeGauss"));


     //Zusammenfassen der Parameter in einem RooArgSet
     RooArgSet Observables;
     Observables.add(RooArgSet( epdf->Var("obsTime"),epdf->Cat("catFinalState"),epdf->Cat("obsTag"),epdf->Var("obsEtaOS")));


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

//.........这里部分代码省略.........
  graph_Hlt1DiMuon_Hlt2DiMuonDetached_B->SetName(tmp);
  tmp = particle;
  tmp += " Hlt1DiMuon_Hlt2DiMuonDetached UB Acceptance";
  histo_Hlt1DiMuon_Hlt2DiMuonDetached_B->SetTitle(tmp);
  graph_Hlt1DiMuon_Hlt2DiMuonDetached_B->SetTitle(tmp);
  
  tmp = prefix;
  tmp += "_";
  tmp += bins;
  tmp += "bins";
  tmp += "_Hlt1TrackAndTrackMuonExcl_Hlt2DiMuonDetached_Background";
  histo_Hlt1TrackAndTrackMuonExcl_Hlt2DiMuonDetached_B->SetName(tmp);
  tmp += "_AsymErrors";
  graph_Hlt1TrackAndTrackMuonExcl_Hlt2DiMuonDetached_B->SetName(tmp);
  tmp = particle;
  tmp += " Hlt1TrackAndTrackMuonExcl_Hlt2DiMuonDetached Acceptance";
  histo_Hlt1TrackAndTrackMuonExcl_Hlt2DiMuonDetached_B->SetTitle(tmp);
  graph_Hlt1TrackAndTrackMuonExcl_Hlt2DiMuonDetached_B->SetTitle(tmp);
  
  RooDataSet* data_sub;
  RooDataSet* data_sub2;
  RooDataSet* data_sub_A_1;
  RooDataSet* data_sub_A_2;
  RooDataSet* data_sub_B_1;
  RooDataSet* data_sub_B_2;
  
  data_sub2 = ((RooDataSet*)data->reduce("CatB>0."));
  
  n_A_1_S.setVal(TMath::Max(data_sub2->sumEntries()/2.,100.));
  n_A_1_B.setVal(TMath::Max(data_sub2->sumEntries()/2.,100.));
  
  model_A_1->fitTo(*data_sub2,RooFit::Minos(kFALSE),Optimize(kFALSE),RooFit::SumW2Error(kTRUE));
  
  rfrac.setConstant();
  frac.setConstant();
  Bmass.setConstant();
  c0.setConstant();
  Bwidth1.setConstant();
  data_sub2->Delete();
  
  for(Int_t i(0); i< bins; i++){
    tmp = "B_TAU < ";
    tmp += binning->binHigh(i);
    tmp += " && B_TAU > ";
    tmp += binning->binLow(i);
    data_sub = (RooDataSet*)data->reduce(tmp);
    data_sub2 = ((RooDataSet*)data_sub->reduce("CatB>0."));
    data_sub_A_1 = ((RooDataSet*)data_sub->reduce("CatA==CatA::CatA1"));
    data_sub_A_2 = ((RooDataSet*)data_sub->reduce("CatA==CatA::CatA2"));
    data_sub_B_1 = ((RooDataSet*)data_sub->reduce("CatB==CatB::CatB1"));
    data_sub_B_2 = ((RooDataSet*)data_sub->reduce("CatB==CatB::CatB2"));
    
    TAU.setRange("fit",binning->binLow(i),binning->binHigh(i));
    TAU.setBin(i);
    
    n_A_1_S.setVal(TMath::Max(data_sub_A_1->sumEntries()/2.,100.));
    n_A_1_B.setVal(TMath::Max(data_sub_A_1->sumEntries()/2.,100.));
    
    c0.setConstant(kFALSE);
    //Bwidth1.setConstant(kFALSE);
    
    //n_A_1_S.setMin(-1000.);
    model_A_1->fitTo(*data_sub2,RooFit::Minos(kFALSE),Optimize(kFALSE),RooFit::SumW2Error(kTRUE));
    n_A_1_S.setMin(0.0000001);
    
    //c0.setConstant();

//-------------------------------------------------------------
//Main macro for generating data and fitting
//=============================================================  
void FitMassPhotonResolutionSystematics(const string workspaceFile = "/afs/cern.ch/work/d/daan/public/releases/CMSSW_5_3_9_patch3/src/CMSAna/HHToBBGG/data/FitWorkspace_asdf.root", const string outputTree = "/afs/cern.ch/work/d/daan/public/releases/CMSSW_5_3_9_patch3/src/CMSAna/HHToBBGG/data/MassFitResults/MassFitTwoD_asdf.root", Int_t plotOption = 0, Int_t storeOption = 1, Int_t SystematicsUpOrDown = 1, Int_t NToys = 5000) {

  TRandom3 *randomnumber = new TRandom3(1200);
  TFile *wsFile = new TFile (workspaceFile.c_str(), "READ");
  RooWorkspace *ws = (RooWorkspace*)wsFile->Get("MassFitWorkspace");
  
  //Import variables from workspace
  RooAbsPdf *model2Dpdf = ws->pdf("model2Dpdf");
  RooRealVar *massBjet = ws->var("massBjet");
  RooRealVar *massPho = ws->var("massPho");
  RooRealVar *nsig = ws->var("N (Sig)"); RooRealVar constNsig(*nsig);
  RooRealVar *nres = ws->var("N (ResBkg)"); RooRealVar constNres(*nres);
  RooRealVar *nnonres = ws->var("N (NonResBkg)"); RooRealVar constNnonres(*nnonres);
  RooRealVar *expRateBjet = ws->var("expRateBjet"); RooRealVar constexpBjet(*expRateBjet);
  RooRealVar *expRatePho = ws->var("expRatePho"); RooRealVar constexpPho(*expRatePho);

  //Variables to set constant
  RooRealVar *sigMeanBjet = ws->var("sigMeanBjet"); sigMeanBjet->setConstant();
  RooRealVar *sigSigmaBjet = ws->var("sigSigmaBjet"); sigSigmaBjet->setConstant();
  RooRealVar *sigAlpha = ws->var("sigAlpha"); sigAlpha->setConstant();
  RooRealVar *sigPower = ws->var("sigPower"); sigPower->setConstant();
  RooRealVar *sigmaPho = ws->var("sigmaPho"); sigmaPho->setConstant();

  RooRealVar *resMeanBjet = ws->var("resMeanBjet"); resMeanBjet->setConstant();
  RooRealVar *resSigmaBjet = ws->var("resSigmaBjet"); resSigmaBjet->setConstant();
  RooRealVar *resAlpha = ws->var("resAlpha"); resAlpha->setConstant();
  RooRealVar *resPower = ws->var("resPower"); resPower->setConstant();
  RooRealVar *resExpo = ws->var("resExpo"); resExpo->setConstant();
  RooRealVar *nbbH = ws->var("nbbH"); nbbH->setConstant();
  RooRealVar *nOthers = ws->var("nOthers"); nOthers->setConstant();
  
  double inputSigmaPho = sigmaPho->getVal();

  //Create TTree to store the resulting yield data
  TFile *f = new TFile(outputTree.c_str(), "RECREATE");
  TTree *resultTree = new TTree("resultTree", "Parameter results from fitting");
  Float_t nsigOut, nresOut, nnonresOut;
  Float_t nsigStd, nresStd, nnonresStd;
  
  resultTree->Branch("nsigOut",&nsigOut, "nsigOut/F");
  resultTree->Branch("nresOut",&nresOut, "nresOut/F");
  resultTree->Branch("nnonresOut",&nnonresOut, "nnonresOut/F");
  resultTree->Branch("nsigStd",&nsigStd, "nsigStd/F");
  resultTree->Branch("nresStd",&nresStd, "nresStd/F");
  resultTree->Branch("nnonresStd",&nnonresStd, "nnonresStd/F");
  
  //Generate Toy MC experiment data and fits
  for(UInt_t t=0; t < NToys; ++t) {
    cout << "Toy #" << t << endl;
    nsig->setVal(constNsig.getVal()); nres->setVal(constNres.getVal()); nnonres->setVal(constNnonres.getVal());
    expRateBjet->setVal(constexpBjet.getVal()); expRatePho->setVal(constexpPho.getVal());
   
    //set jet energy resolutions to nominal
    sigmaPho->setVal(inputSigmaPho);

    cout << "Before: " << sigmaPho->getVal() << " | ";

    Float_t ran = randomnumber->Poisson(325.);
    RooDataSet *pseudoData2D = model2Dpdf->generate(RooArgList(*massBjet,*massPho), ran);

    //move jet energy resolution up/down
    if (SystematicsUpOrDown == 1) {
      sigmaPho->setVal(inputSigmaPho*1.15);
    } else if (SystematicsUpOrDown == -1) {
      sigmaPho->setVal(inputSigmaPho/1.15);
    }

    cout << "After: " << sigmaPho->getVal() << " \n";

    RooFitResult *fitResult2D = model2Dpdf->fitTo(*pseudoData2D, RooFit::Save(), RooFit::Extended(kTRUE), RooFit::Strategy(2));
//     if (t == 1763) {
//     	ws->import(*pseudoData2D, kTRUE);
//     	ws->import(*pseudoData2D, Rename("pseudoData2D"));
//     }
//     if (plotOption == 1) MakePlots(ws, fitResult2D);
    

    cout << "DEBUG: " << constexpBjet.getVal() << " , " << constexpPho.getVal() << " | " << expRateBjet->getVal() << " " << expRatePho->getVal() << "\n";


    //Store fit parameters into ROOT file
    if (storeOption == 1) {
      //Save variables into separate branches of root tree
  		nsigOut = nsig->getVal();
  		nresOut = nres->getVal();
  		nnonresOut = nnonres->getVal();
  		nsigStd = nsig->getPropagatedError(*fitResult2D);
  		nresStd = nres->getPropagatedError(*fitResult2D);
  		nnonresStd = nnonres->getPropagatedError(*fitResult2D);
  		//cout << "\n\n\n\n\n\n" << nsigOut << " | " << nresOut << " | " << nnonresOut << " | " << ran  << "\n\n\n\n\n" << endl;
  		resultTree->Fill();
    }
  }
  
  //Write to the TTree and close it
  resultTree->Write();
  f->Close();
//.........这里部分代码省略.........

MakeAICFits::MakeAICFits() {
  //RooRandom::randomGenerator()->SetSeed(314159);
  // Create Toy Data Set
  // Single Exponential
  RooRealVar *mass   = new RooRealVar("mass","mass", 50., 35., 65.); 
  //RooRealVar *alpha1 = new RooRealVar("Exp_alpha","Exp_alpha",-0.1,-1.,0.);                     
  //alpha1->setVal(-0.1);
  RooRealVar *alpha1 = new RooRealVar("Exp_alpha","Exp_alpha", -0.1);
  RooExponential expon("exp","Background Exponential",*mass,*alpha1);             
  // KPower 
  //RooRealVar *alpha2 = new RooRealVar("Power_alpha","Power_alpha",-3.,-10.,0.);                             
  //alpha2->setVal(-3.);
  RooRealVar *alpha2 = new RooRealVar("Power_alpha", "Power_alpha", -3);
  RooGenericPdf bkg("pow","Background Power","@0^@1",RooArgList(*mass,*alpha2));           
  //RooRealVar ratio("ratio","Background Ratio", 0.5, 0., 1.);
  //ratio.setVal(0.5);
  RooRealVar ratio("ratio", "Background Ratio", 0.5);


  //Create Background pdf
  //  RooAddPdf bkg("bkg","bkg",RooArgList(pow,expon),ratio);

  std::cout<<"==========  Data Model Made  ==========="<<std::endl;

  const Int_t nToys = 1;
  //ratio.setVal(0);
  RooDataSet* data = bkg.generate(RooArgSet(*mass), 1000000);
  //ratio.setVal(0.20);


  std::cout<<"==========  Data Set Made    ==========="<<std::endl;
  
  // Make plain projection of data and pdf on mass
  //bkg.fitTo(*data);
  
  RooFitResult* bkg_data = bkg.fitTo(*data, RooFit::Save(kTRUE), RooFit::Optimize(0));
  Double_t bkg_data_Nll = bkg_data->minNll();
  std::cout<<" ======== Data fitted ==========="<<std::endl;
  RooArgSet* bkgParams = bkg.getParameters(*data);
  const RooArgList& fitbkgParams = bkg_data->floatParsFinal();
  std::cout<< "=======================  parameters done  ========================"<<std::endl;   

  Double_t LogLikelihood[8] = {0,0,0,0,0,0,0,0};
  Double_t AIC_bkg_array[8] = {0,0,0,0,0,0,0,0};
  Double_t AICc_bkg_array[8] = {0,0,0,0,0,0,0,0};
  //Double_t BIC_bkg_array[7] = {0,0,0,0,0,0,0};
  
  Double_t LogLikelihood_data = bkg_data_Nll;
  Int_t avgcov[8] = {0,0,0,0,0,0,0};
  std::cout<<"======================   Starting Toys  ==============="<<std::endl;
  for (Int_t i=0; i<nToys; i++) {
    if (i%10==0) {
      std::cout<<">> Processing Toy Trial " << i << "/" << nToys << std::endl;
    }

    RooPlot* frame = mass->frame();
    leg = new TLegend(0.55,0.55,0.9,0.9);
    mass->setConstant(kFALSE);
    alpha1->setConstant(kFALSE);
    alpha2->setConstant(kFALSE);
    ratio.setConstant(kFALSE);
    const RooArgList ranbkgParams = bkg_data->randomizePars();
    *bkgParams = ranbkgParams;
    
    Int_t SampleSize = 100000;
    RooDataSet* toybkg = bkg.generate(RooArgSet(*mass),SampleSize);
    toybkg->plotOn(frame);
    leg->AddEntry(toybkg,"Toy Background", "lep");
    *bkgParams = fitbkgParams;
    mass->setConstant(kTRUE);
    alpha1->setConstant(kTRUE);
    alpha2->setConstant(kTRUE);
    ratio.setConstant(kTRUE);
    for (int type=0; type<7; type++) {
      std::cout<<type<<endl;
    }
    int display = 8;
    for (int type=0; type<display; type++) {
      if (type<7) {
	//std::cout<<"Model Shape:    "<<type<<std::endl;
	RooAbsPdf* ModelShape = MakeAICFits::getBackgroundPdf(type,mass);
	//std::cout<<"Model Shape made"<<std::endl;
	int k = MakeAICFits::Num_Params(type);
	//std::cout<<"Params counted"<<std::endl;
      }
      if (type==7) {
	RooAbsPdf* Model1 = MakeAICFits::getBackgroundPdf(5,mass);
	RooAbsPdf* Model2 = MakeAICFits::getBackgroundPdf(6,mass);
	RooAbsPdf* Model3 = MakeAICFits::getBackgroundPdf(2,mass);
	RooAbsPdf* Model4 = MakeAICFits::getBackgroundPdf(3,mass);
	RooAbsPdf* Model5 = MakeAICFits::getBackgroundPdf(4,mass);
	int k = MakeAICFits::Num_Params(5);
	k+= MakeAICFits::Num_Params(6);
	//k+= MakeAICFits::Num_Params(0);
	//k+= MakeAICFits::Num_Params(3);
	//k+= MakeAICFits::Num_Params(4);
	RooRealVar* modratio1 = new RooRealVar("modrat1", "modrat1", 0.82, 0.81, 0.83);
	RooRealVar* modratio2 = new RooRealVar("modrat2", "modrat2", 0.17, 0.25, 0.35);
	RooRealVar* modratio3 = new RooRealVar("modrat3", "modrat3", 0.01);
	//RooRealVar* modratio4 = new RooRealVar("modrat4", "modrat4", 0.25);
//.........这里部分代码省略.........

void CreateDataTemplates(double dX,int BRN_ORDER)
{
  gROOT->ForceStyle();
  RooMsgService::instance().setSilentMode(kTRUE);
  for(int i=0;i<2;i++) {
    RooMsgService::instance().setStreamStatus(i,kFALSE);
  }
  double XMIN = 80;
  double XMAX = 200;
  const int NSEL(2);
  const int NCAT[NSEL] = {4,3};
  const double MVA_BND[NSEL][NCAT[0]+1] = {{-0.6,0.0,0.7,0.84,1},{-0.1,0.4,0.8,1}};
  char name[1000];
  TString SELECTION[NSEL] = {"NOM","VBF"};
  TString SELTAG[NSEL]    = {"NOM","PRK"};
  TString MASS_VAR[NSEL]  = {"mbbReg[1]","mbbReg[2]"};

  TFile *fBKG  = TFile::Open("limit_BRN5+4_dX0p1_80-200_CAT0-6/output/bkg_shapes_workspace.root");
  RooWorkspace *wBkg = (RooWorkspace*)fBKG->Get("w");
  RooWorkspace *w = new RooWorkspace("w","workspace");
  //RooRealVar x(*(RooRealVar*)wBkg->var("mbbReg"));
  TTree *tr;
  TH1F *h,*hBlind;
  TCanvas *canFit[5]; 
  RooDataHist *roohist[5],*roohist_blind[5];

  TFile *fTransfer = TFile::Open("limit_BRN5+4_dX0p1_80-200_CAT0-6/output/transferFunctions.root");
  TF1 *transFunc;

  int counter(0);
  int NPAR = BRN_ORDER;

  for(int isel=0;isel<NSEL;isel++) {
    TFile *fDATA = TFile::Open("flat/Fit_data_sel"+SELECTION[isel]+".root");
    RooRealVar *brn[8];
    RooArgSet brn_params;
    if (isel == 1) {
      NPAR = 4;
    } 
    for(int ib=0;ib<=NPAR;ib++) {
      brn[ib] = new RooRealVar("b"+TString::Format("%d",ib)+"_sel"+SELECTION[isel],"b"+TString::Format("%d",ib)+"_sel"+SELECTION[isel],0.5,0,10.);
      brn_params.add(*brn[ib]);
    }
    for(int icat=0;icat<NCAT[isel];icat++) {
      RooRealVar x("mbbReg_"+TString::Format("CAT%d",counter),"mbbReg_"+TString::Format("CAT%d",counter),XMIN,XMAX);

      sprintf(name,"fitRatio_sel%s_CAT%d",SELTAG[isel].Data(),counter);
      transFunc = (TF1*)fTransfer->Get(name);
		transFunc->Print();
      // --- The error on the tranfer function parameters is shrinked because the correlations are ingored. 
      // --- Must be consistent with TransferFunctions.C
      float p0 = transFunc->GetParameter(0);
      float e0 = transFunc->GetParError(0);
      float p1 = transFunc->GetParameter(1);
      float e1 = transFunc->GetParError(1);
      float p2 = transFunc->GetParameter(2);
      float e2 = transFunc->GetParError(2);
      RooRealVar trans_p2(TString::Format("trans_p2_CAT%d",counter),TString::Format("trans_p2_CAT%d",counter),p2);
      RooRealVar trans_p1(TString::Format("trans_p1_CAT%d",counter),TString::Format("trans_p1_CAT%d",counter),p1);
      RooRealVar trans_p0(TString::Format("trans_p0_CAT%d",counter),TString::Format("trans_p0_CAT%d",counter),p0);
		printf("%.2f %.2f %.2f\n",p0,p1,p2);
      RooGenericPdf *transfer;
      if (isel == 0) {
        trans_p2.setError(0.5*e2);
        trans_p1.setError(0.5*e1);
        trans_p0.setError(0.5*e0);
        transfer = new RooGenericPdf(TString::Format("transfer_CAT%d",counter),"@2*@[email protected]",RooArgList(x,trans_p0,trans_p1)); 
      }
      else {
        trans_p2.setError(0.05*e2);
        trans_p1.setError(0.05*e1);
        trans_p0.setError(0.05*e0);
        transfer = new RooGenericPdf(TString::Format("transfer_CAT%d",counter),"@3*@0*@[email protected]*@[email protected]",RooArgList(x,trans_p0,trans_p1,trans_p2));
      }
      trans_p2.setConstant(kTRUE);
      trans_p1.setConstant(kTRUE);
      trans_p0.setConstant(kTRUE);

		transfer->Print();
            
      sprintf(name,"FitData_sel%s_CAT%d",SELECTION[isel].Data(),icat);
      canFit[icat] = new TCanvas(name,name,900,600);
      canFit[icat]->cd(1)->SetBottomMargin(0.4);
      sprintf(name,"Hbb/events");
      tr = (TTree*)fDATA->Get(name); 
      sprintf(name,"hMbb_%s_CAT%d",SELECTION[isel].Data(),icat);
      int NBINS = (XMAX[isel][icat]-XMIN[isel][icat])/dX;
      
      h = new TH1F(name,name,NBINS,XMIN[isel][icat],XMAX[isel][icat]);

      sprintf(name,"hMbb_blind_%s_CAT%d",SELECTION[isel].Data(),icat);
      hBlind = new TH1F(name,name,NBINS,XMIN[isel][icat],XMAX[isel][icat]);

      sprintf(name,"mva%s>%1.2f && mva%s<=%1.2f",SELECTION[isel].Data(),MVA_BND[isel][icat],SELECTION[isel].Data(),MVA_BND[isel][icat+1]);
      TCut cut(name);
      sprintf(name,"mva%s>%1.2f && mva%s<=%1.2f && %s>100 && %s<150",SELECTION[isel].Data(),MVA_BND[isel][icat],SELECTION[isel].Data(),MVA_BND[isel][icat+1],MASS_VAR[isel].Data(),MASS_VAR[isel].Data());
      TCut cutBlind(name);
      tr->Draw(MASS_VAR[isel]+">>"+h->GetName(),cut); 
      tr->Draw(MASS_VAR[isel]+">>"+hBlind->GetName(),cutBlind); 
      sprintf(name,"yield_data_CAT%d",counter);
//.........这里部分代码省略.........

void forData(string channel, string catcut, bool removeMinor=true){

  // Suppress all the INFO message

  RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);

  // Input files and sum all backgrounds

  TChain* treeData  = new TChain("tree");
  TChain* treeZjets = new TChain("tree");

  if( channel == "ele" ){

    treeData->Add(Form("%s/data/SingleElectron-Run2015D-05Oct2015-v1_toyMCnew.root",  channel.data()));
    treeData->Add(Form("%s/data/SingleElectron-Run2015D-PromptReco-V4_toyMCnew.root", channel.data()));

  }

  else if( channel == "mu" ){

    treeData->Add(Form("%s/data/SingleMuon-Run2015D-05Oct2015-v1_toyMCnew.root",  channel.data()));
    treeData->Add(Form("%s/data/SingleMuon-Run2015D-PromptReco-V4_toyMCnew.root", channel.data()));

  }

  else return;

  treeZjets->Add(Form("%s/Zjets/DYJetsToLL_M-50_HT-100to200_13TeV_toyMCnew.root", channel.data()));
  treeZjets->Add(Form("%s/Zjets/DYJetsToLL_M-50_HT-200to400_13TeV_toyMCnew.root", channel.data()));
  treeZjets->Add(Form("%s/Zjets/DYJetsToLL_M-50_HT-400to600_13TeV_toyMCnew.root", channel.data()));
  treeZjets->Add(Form("%s/Zjets/DYJetsToLL_M-50_HT-600toInf_13TeV_toyMCnew.root", channel.data()));

  // To remove minor background contribution in data set (weight is -1)

  if( removeMinor ){

    treeData->Add(Form("%s/VV/WW_TuneCUETP8M1_13TeV_toyMCnew.root", channel.data()));
    treeData->Add(Form("%s/VV/WZ_TuneCUETP8M1_13TeV_toyMCnew.root", channel.data()));
    treeData->Add(Form("%s/VV/ZZ_TuneCUETP8M1_13TeV_toyMCnew.root", channel.data()));
    treeData->Add(Form("%s/TT/TT_TuneCUETP8M1_13TeV_toyMCnew.root", channel.data()));

  }

  // Define all the variables from the trees

  RooRealVar cat ("cat", "", 0, 2);
  RooRealVar mJet("prmass", "M_{jet}",  30.,  300., "GeV");
  RooRealVar mZH ("mllbb",   "M_{ZH}", 900., 3000., "GeV");
  RooRealVar evWeight("evweight", "", -1.e3, 1.e3);

  // Set the range in jet mass

  mJet.setRange("allRange", 30., 300.);
  mJet.setRange("lowSB",    30.,  65.);
  mJet.setRange("highSB",  135., 300.);
  mJet.setRange("signal",  105., 135.);

  RooBinning binsmJet(54, 30, 300);

  RooArgSet variables(cat, mJet, mZH, evWeight);

  TCut catCut = Form("cat==%s", catcut.c_str());
  TCut sbCut  = "prmass>30 && !(prmass>65 && prmass<135) && prmass<300";
  TCut sigCut = "prmass>105 && prmass<135";

  // Create a dataset from a tree -> to process an unbinned likelihood fitting

  RooDataSet dataSetData   ("dataSetData",    "dataSetData",    variables, Cut(catCut),           WeightVar(evWeight), Import(*treeData));
  RooDataSet dataSetDataSB ("dataSetDataSB",  "dataSetDataSB",  variables, Cut(catCut && sbCut),  WeightVar(evWeight), Import(*treeData));
  RooDataSet dataSetZjets  ("dataSetZjets",   "dataSetZjets",   variables, Cut(catCut),           WeightVar(evWeight), Import(*treeZjets));
  RooDataSet dataSetZjetsSB("dataSetZjetsSB", "dataSetZjetsSB", variables, Cut(catCut && sbCut),  WeightVar(evWeight), Import(*treeZjets));  
  RooDataSet dataSetZjetsSG("dataSetZjetsSG", "dataSetZjetsSG", variables, Cut(catCut && sigCut), WeightVar(evWeight), Import(*treeZjets));
  
  // Total events number

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

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

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

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

  // Signal region jet mass

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

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

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

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

void compute_p0(const char* inFileName,
	    const char* wsName = "combined",
	    const char* modelConfigName = "ModelConfig",
	    const char* dataName = "obsData",
	    const char* asimov1DataName = "asimovData_1",
	    const char* conditional1Snapshot = "conditionalGlobs_1",
	    const char* nominalSnapshot = "nominalGlobs",
	    string smass = "130",
	    string folder = "test")
{
  double mass;
  stringstream massStr;
  massStr << smass;
  massStr >> mass;

  double mu_profile_value = 1; // mu value to profile the obs data at wbefore generating the expected
  bool doConditional      = 1; // do conditional expected data
  bool remakeData         = 0; // handle unphysical pdf cases in H->ZZ->4l
  bool doUncap            = 1; // uncap p0
  bool doInj              = 0; // setup the poi for injection study (zero is faster if you're not)
  bool doObs              = 1; // compute median significance
  bool doMedian           = 1; // compute observed significance

  TStopwatch timer;
  timer.Start();

  TFile f(inFileName);
  RooWorkspace* ws = (RooWorkspace*)f.Get(wsName);
  if (!ws)
  {
    cout << "ERROR::Workspace: " << wsName << " doesn't exist!" << endl;
    return;
  }
  ModelConfig* mc = (ModelConfig*)ws->obj(modelConfigName);
  if (!mc)
  {
    cout << "ERROR::ModelConfig: " << modelConfigName << " doesn't exist!" << endl;
    return;
  }
  RooDataSet* data = (RooDataSet*)ws->data(dataName);
  if (!data)
  {
    cout << "ERROR::Dataset: " << dataName << " doesn't exist!" << endl;
    return;
  }

  mc->GetNuisanceParameters()->Print("v");

  ROOT::Math::MinimizerOptions::SetDefaultMinimizer("Minuit2");
  ROOT::Math::MinimizerOptions::SetDefaultStrategy(0);
  ROOT::Math::MinimizerOptions::SetDefaultPrintLevel(1);
  cout << "Setting max function calls" << endl;

  ws->loadSnapshot("conditionalNuis_0");
  RooArgSet nuis(*mc->GetNuisanceParameters());

  RooRealVar* mu = (RooRealVar*)mc->GetParametersOfInterest()->first();

  RooAbsPdf* pdf = mc->GetPdf();

  string condSnapshot(conditional1Snapshot);
  RooArgSet nuis_tmp2 = *mc->GetNuisanceParameters();
  RooNLLVar* obs_nll = doObs ? (RooNLLVar*)pdf->createNLL(*data, Constrain(nuis_tmp2)) : NULL;

  RooDataSet* asimovData1 = (RooDataSet*)ws->data(asimov1DataName);
  RooRealVar* emb = (RooRealVar*)mc->GetNuisanceParameters()->find("ATLAS_EMB");
  if (!asimovData1 || (string(inFileName).find("ic10") != string::npos && emb))
  {
    if (emb) emb->setVal(0.7);
    cout << "Asimov data doesn't exist! Please, allow me to build one for you..." << endl;
    string mu_str, mu_prof_str;
    asimovData1 = makeAsimovData(mc, doConditional, ws, obs_nll, 1, &mu_str, &mu_prof_str, mu_profile_value, true);
    condSnapshot="conditionalGlobs"+mu_prof_str;
  }
  
  if (!doUncap) mu->setRange(0, 40);
  else mu->setRange(-40, 40);

  RooAbsPdf* pdf = mc->GetPdf();

  RooArgSet nuis_tmp1 = *mc->GetNuisanceParameters();
  RooNLLVar* asimov_nll = (RooNLLVar*)pdf->createNLL(*asimovData1, Constrain(nuis_tmp1));

  //do asimov
  mu->setVal(1);
  mu->setConstant(0);
  if (!doInj) mu->setConstant(1);

  int status,sign;
  double med_sig=0,obs_sig=0,asimov_q0=0,obs_q0=0;

  if (doMedian)
  {
    ws->loadSnapshot(condSnapshot.c_str());
    if (doInj) ws->loadSnapshot("conditionalNuis_inj");
    else ws->loadSnapshot("conditionalNuis_1");
    mc->GetGlobalObservables()->Print("v");
    mu->setVal(0);
    mu->setConstant(1);
    status = minimize(asimov_nll, ws);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    RooNLLVar::SetIgnoreZeroEntries(0);


    //gRandom->SetSeed(1);
    //RooRandom::randomGenerator()->SetSeed(1);

    //RooMinimizerFcn::SetOverrideEverything(false);

    cout << "med test stat: " << med_testStat_val << endl;
    ws->loadSnapshot("nominalGlobs");

    ws->loadSnapshot("conditionalNuis_0");
    mu->setVal(0);


    testStat.SetWorkspace(ws);
    testStat.SetLoadUncondSnapshot("ucmles");
    double obsTestStat_val = doObs ? 2*testStat.Evaluate(*data, poi) : 0;
    cout << "obs test stat: " << obsTestStat_val << endl;
//   obsTestStat_val = 2*testStat.Evaluate(*data, poi);
//   cout << "obs test stat: " << obsTestStat_val << endl;
//   obsTestStat_val = 2*testStat.Evaluate(*data, poi);
//   cout << "obs test stat: " << obsTestStat_val << endl;

    double obs_sig;
    int sign = obsTestStat_val == 0 ? 0 : obsTestStat_val / fabs(obsTestStat_val);
    if (!doRightSided && (obsTestStat_val < 0 && obsTestStat_val > -0.1 || mu->getVal() < 0.001)) obs_sig = 0; 
    else obs_sig = sign*sqrt(fabs(obsTestStat_val));
    if (obs_sig != obs_sig) //nan, do by hand
    {
      cout << "Obs test stat gave nan: try by hand" << endl;

      mu->setVal(0);
      mu->setConstant(1);
      mc->GetPdf()->fitTo(*data, Hesse(0), Minos(0), PrintLevel(-1), Constrain(*mc->GetNuisanceParameters()));
      mu->setConstant(0);

      double L_0 = mc->GetPdf()->getVal();

      //mu->setVal(0);
      //mu->setConstant(1);
      mc->GetPdf()->fitTo(*data, Hesse(0), Minos(0), PrintLevel(-1), Constrain(*mc->GetNuisanceParameters()));
      //mu->setConstant(0);
      double L_muhat = mc->GetPdf()->getVal();

      cout << "L_0: " << L_0 << ", L_muhat: " << L_muhat << endl;
      obs_sig = sqrt(-2*TMath::Log(L_0/L_muhat));

//still nan
      if (obs_sig != obs_sig && fabs(L_0 - L_muhat) < 0.000001) obs_sig = 0;
    }
    cout << "obs: " << obs_sig << endl;

    cout << "Observed significance: " << obs_sig << endl;
    if (med_sig)
    {
      cout << "Median test stat val: " << med_testStat_val << endl;
      cout << "Median significance:   " << med_sig << endl;
    }


    f.Close();

    stringstream fileName;
    fileName << "root_files/" << folder << "/" << mass << ".root";
    system(("mkdir -vp root_files/" + folder).c_str());

void RooToyMCFit(){
    stringstream out;
   out.str("");
   out << "toyMC_" << _sigma_over_tau << "_" << _purity << "_" << mistag_rate << ".root";
   TFile* file = TFile::Open(out.str().c_str());
   TTree* tree = (TTree*)file->Get("ToyTree");
   tree->Print();

//    TFile* file = TFile::Open("ToyData/toyMC_453_0.912888_0.3_0.8_0.root");
//    TFile* file = TFile::Open("ToyData/toyMC_10000_0.9_0.3_0.8_0.root");
//    TTree* tree = (TTree*)file->Get("ToyTree1");

    RooRealVar tau("tau","tau",_tau,"ps"); tau.setConstant(kTRUE);
    RooRealVar dm("dm","dm",_dm,"ps^{-1}"); dm.setConstant(kTRUE);

    RooAbsReal* sin2beta;
    RooAbsReal* cos2beta;
    if(softlimit){
      RooRealVar x("x","x",_cos2beta,-3.,30000.); if(constBeta) x.setConstant(kTRUE);
      RooRealVar y("y","y",_sin2beta,-3.,30000.); if(constBeta) y.setConstant(kTRUE);
      cos2beta = new RooFormulaVar("cos2beta","cos2beta","@0/sqrt(@0*@[email protected]*@1)",RooArgSet(x,y));
      sin2beta = new RooFormulaVar("sin2beta","sin2beta","@0/sqrt(@0*@[email protected]*@1)",RooArgSet(y,x));
//      sin2beta = new RooFormulaVar("sin2beta","sin2beta","2/TMath::Pi()*TMath::ATan(@0)",RooArgSet(y));
//      cos2beta = new RooFormulaVar("cos2beta","cos2beta","2/TMath::Pi()*TMath::ATan(@0)",RooArgSet(x));
    }
    else{ 
      sin2beta = new RooRealVar("sin2beta","sin2beta",_sin2beta,-3.,3.); if(constBeta) ((RooRealVar*)sin2beta)->setConstant(kTRUE);
      cos2beta = new RooRealVar("cos2beta","cos2beta",_cos2beta,-3.,3.); if(constBeta) ((RooRealVar*)cos2beta)->setConstant(kTRUE);
    }
    RooRealVar dt("dt","#Deltat",-5.,5.,"ps");
    RooRealVar avgMisgat("avgMisgat","avgMisgat",mistag_rate,0.0,0.5); if(constMistag) avgMisgat.setConstant(kTRUE);
    RooRealVar delMisgat("delMisgat","delMisgat",0); delMisgat.setConstant(kTRUE);
    RooRealVar mu("mu","mu",0); mu.setConstant(kTRUE);

    tau.Print();
    dm.Print();
    sin2beta->Print();
    cos2beta->Print();
    avgMisgat.Print();

    RooRealVar moment("moment","moment",0.); moment.setConstant(kTRUE);
    RooRealVar parity("parity","parity",-1.); parity.setConstant(kTRUE);

    RooRealVar* K[8];
    RooAbsReal* Kb[8];
    RooArgList Kset;
    RooRealVar* C[8];
    RooRealVar* S[8];
    RooFormulaVar* a1[2][8];
    RooFormulaVar* b1[2][8];
    RooFormulaVar* a2[2][8];
    RooFormulaVar* b2[2][8];

    for(int i=0; i<8; i++){
        out.str("");
        out << "K" << i+1;
        K[i] = new RooRealVar(out.str().c_str(),out.str().c_str(),_K[i],0.,1.); Kset.add(*K[i]); if(constK) K[i]->setConstant(kTRUE);
        K[i]->Print();
        if(i!=7){
          out.str("");
          out << "Kb" << i+1;
          Kb[i] = new RooRealVar(out.str().c_str(),out.str().c_str(),_Kb[i],0.,1.); Kset.add(*Kb[i]); if(constK) ((RooRealVar*)Kb[i])->setConstant(kTRUE);
          Kb[i]->Print();
        }
        out.str("");
        out << "C" << i+1;
        C[i] = new RooRealVar(out.str().c_str(),out.str().c_str(),_C[i]); C[i]->setConstant(kTRUE);
        C[i]->Print();
        out.str("");
        out << "S" << i+1;
        S[i] = new RooRealVar(out.str().c_str(),out.str().c_str(),_S[i]); S[i]->setConstant(kTRUE);
        S[i]->Print();
    }
    Kb[7] = new RooFormulaVar("K8b","K8b","[email protected]@[email protected]@[email protected]@[email protected]@[email protected]@[email protected]@[email protected]@[email protected]",Kset);

    for(int i=0; i<8; i++){
        out.str("");
        out << "a10_" << i+1;
        a1[0][i] = new RooFormulaVar(out.str().c_str(),out.str().c_str(),"-(@[email protected])/(@[email protected])",RooArgList(*K[i],*Kb[i]));
        a1[0][i]->Print();
        out.str("");
        out << "a11_" << i+1;
        a1[1][i] = new RooFormulaVar(out.str().c_str(),out.str().c_str(),"(@[email protected])/(@[email protected])",RooArgList(*K[i],*Kb[i]));
        a1[1][i]->Print();
        out.str("");
        out << "a20_" << i+1;
        a2[0][i] = new RooFormulaVar(out.str().c_str(),out.str().c_str(),"(@[email protected])/(@[email protected])",RooArgList(*K[i],*Kb[i]));
        a2[0][i]->Print();
        out.str("");
        out << "a21_" << i+1;
        a2[1][i] = new RooFormulaVar(out.str().c_str(),out.str().c_str(),"-(@[email protected])/(@[email protected])",RooArgList(*K[i],*Kb[i]));
        a2[1][i]->Print();

        out.str("");
        out << "b10_" << i+1;
        b1[0][i] = new RooFormulaVar(out.str().c_str(),out.str().c_str(),"2.*(@2*@[email protected]*@5)*TMath::Sqrt(@0*@1)/(@[email protected])",RooArgList(*K[i],*Kb[i],*C[i],*S[i],*sin2beta,*cos2beta));
        b1[0][i]->Print();
        out.str("");
        out << "b11_" << i+1;
        b1[1][i] = new RooFormulaVar(out.str().c_str(),out.str().c_str(),"2.*(@2*@[email protected]*@5)*TMath::Sqrt(@0*@1)/(@[email protected])",RooArgList(*K[i],*Kb[i],*C[i],*S[i],*sin2beta,*cos2beta));
//.........这里部分代码省略.........

double* fitHist(TCanvas *iC, bool iDoMu,int iPlot, std::string iName,TH1D *iData,TH1D *iW,TH1D *iEWK,TH1D *iAntiData,TH1D *iAntiW,TH1D *iAntiEWK,const Double_t METMAX,const Int_t NBINS,const Double_t lumi,const Int_t Ecm,int iAltQCD) {
    //
    // Declare fit parameters for signal and background yields
    // Note: W signal and EWK+top PDFs are constrained to the ratio described in MC
    //
    RooRealVar nSig(("nSig"+iName).c_str(),("nSig"+iName).c_str(),0.7*(iData->Integral()),0,iData->Integral());
    RooRealVar nQCD(("nQCD"+iName).c_str(),("nQCD"+iName).c_str(),0.3*(iData->Integral()),0,iData->Integral());
    RooRealVar cewk(("cewk"+iName).c_str(),("cewk"+iName).c_str(),0.1,0,5) ;
    cewk.setVal(iEWK->Integral()/iW->Integral());
    cewk.setConstant(kTRUE);
    RooFormulaVar nEWK(("nEWK"+iName).c_str(),("nEWK"+iName).c_str(),("cewk"+iName+"*nSig"+iName).c_str(),RooArgList(nSig,cewk));

    RooRealVar nAntiSig(("nAntiSig"+iName).c_str(),("nAntiSig"+iName).c_str(),0.05*(iAntiData->Integral()),0,iAntiData->Integral());
    RooRealVar nAntiQCD(("nAntiQCD"+iName).c_str(),("nAntiQCD"+iName).c_str(),0.9 *(iAntiData->Integral()),0,iAntiData->Integral());
    RooRealVar dewk    (("dewk"    +iName).c_str(),("dewk"    +iName).c_str(),0.1,0,5) ;
    dewk.setVal(iAntiEWK->Integral()/iAntiW->Integral());
    dewk.setConstant(kTRUE);
    RooFormulaVar nAntiEWK(("nAntiEWK"+iName).c_str(),("nAntiEWK"+iName).c_str(),("dewk"+iName+"*nAntiSig"+iName).c_str(),RooArgList(nAntiSig,dewk));

    //
    // Construct PDFs for fitting
    //
    RooRealVar pfmet(("pfmet"+iName).c_str(),("pfmet"+iName).c_str(),0,METMAX);
    pfmet.setBins(NBINS);

    // Signal PDFs
    RooDataHist wMet  (("wMET" +iName).c_str(), ("wMET" +iName).c_str(),     RooArgSet(pfmet),iW);
    RooHistPdf  pdfW(( "w"+iName).c_str(), ( "w"+iName).c_str(), pfmet, wMet, 1);
    RooDataHist awMet (("awMET"+iName).c_str(), ("awMET"+iName).c_str(),     RooArgSet(pfmet),iAntiW);
    RooHistPdf apdfW(("aw"+iName).c_str(), ("aw"+iName).c_str(), pfmet,awMet, 1);

    // EWK+top PDFs
    RooDataHist ewkMet (("ewkMET" +iName).c_str(),( "ewkMET"+iName).c_str(), RooArgSet(pfmet),iEWK);
    RooHistPdf  pdfEWK (( "ewk"+iName).c_str(),( "ewk"+iName).c_str(), pfmet,ewkMet, 1);
    RooDataHist aewkMet(("aewkMET"+iName).c_str(),("aewkMET"+iName).c_str(), RooArgSet(pfmet),iAntiEWK);
    RooHistPdf apdfEWK (("aewk"+iName).c_str(),("aewk"+iName).c_str(), pfmet,aewkMet, 1);

    // QCD Pdfs
    CPepeModel0 qcd0 (("qcd0" +iName).c_str(),pfmet);
    CPepeModel1 qcd1 (("qcd1" +iName).c_str(),pfmet);
    CPepeModel2 qcd2 (("qcd2" +iName).c_str(),pfmet);
    CPepeModel0 aqcd0(("aqcd0"+iName).c_str(),pfmet);
    CPepeModel1 aqcd1(("aqcd1"+iName).c_str(),pfmet,qcd1.sigma);
    CPepeModel2 aqcd2(("aqcd2"+iName).c_str(),pfmet);
    RooGenericPdf *lQCD  =  qcd0.model;
    RooGenericPdf *lAQCD = aqcd0.model;
    if(iDoMu) lQCD  = qcd1.model;
    if(iDoMu) lAQCD = aqcd1.model;
    if(iAltQCD == 0) lQCD = qcd0.model;
    if(iAltQCD == 1) lQCD = qcd1.model;
    if(iAltQCD == 2) lQCD = qcd2.model;
    if(iAltQCD == 0) lAQCD = aqcd0.model;
    if(iAltQCD == 1) lAQCD = aqcd1.model;
    if(iAltQCD == 2) lAQCD = aqcd2.model;

    // Signal + Background PDFs
    RooAddPdf pdfMet (("pdfMet"+iName).c_str(), ("pdfMet" +iName).c_str(), RooArgList(pdfW ,pdfEWK ,*lQCD), RooArgList(nSig,nEWK,nQCD));
    RooAddPdf apdfMet(("apdfMet"+iName).c_str(),("apdfMet"+iName).c_str(), RooArgList(apdfW,apdfEWK,*lAQCD), RooArgList(nAntiSig,nAntiEWK,nAntiQCD));

    // PDF for simultaneous fit
    RooCategory rooCat("rooCat","rooCat");
    rooCat.defineType("Select");
    rooCat.defineType("Anti");

    RooSimultaneous pdfTotal("pdfTotal","pdfTotal",rooCat);
    pdfTotal.addPdf(pdfMet, "Select");
    if(iDoMu) pdfTotal.addPdf(apdfMet,"Anti");

    // Perform fits
    RooDataHist dataMet  (("dataMet"+iName).c_str(),("dataMet"+iName).c_str(), RooArgSet(pfmet),iData);
    RooDataHist antiMet  (("antiMet"+iName).c_str(),("antiMet"+iName).c_str(), RooArgSet(pfmet),iAntiData);
    RooDataHist dataTotal(("data"   +iName).c_str(),("data"   +iName).c_str(), RooArgList(pfmet), Index(rooCat),
                          Import("Select", dataMet),
                          Import("Anti",   antiMet));

    RooFitResult *fitRes = 0;
    bool runMinos = kTRUE;
    if(iPlot == 0 || iPlot == 3) runMinos = kFALSE; //Remove Minos when running toys (too damn slow)
    if(!iDoMu) fitRes = pdfMet  .fitTo(dataMet  ,Extended(),Minos(runMinos),Save(kTRUE));
    if( iDoMu) fitRes = pdfTotal.fitTo(dataTotal,Extended(),Minos(runMinos),Save(kTRUE));

    double *lResults = new double[16];
    lResults[0]  = nSig.getVal();
    lResults[1]  = nEWK.getVal();
    lResults[2]  = nQCD.getVal();
    lResults[3]  = nAntiSig.getVal();
    lResults[4]  = nAntiEWK.getVal();
    lResults[5]  = nAntiQCD.getVal();
    if(!iDoMu) lResults[6]  = double(qcd0.sigma->getVal());
    if( iDoMu) lResults[6]  = double(qcd1.sigma->getVal());

    lResults[7]  = 0.;
    if(!iDoMu) lResults[7]  = qcd1.a1->getVal();
    lResults[8]   = nSig    .getPropagatedError(*fitRes);
    lResults[9]   = nEWK    .getPropagatedError(*fitRes);
    lResults[10]  = nQCD    .getPropagatedError(*fitRes);
    lResults[11]  = nAntiSig.getPropagatedError(*fitRes);
    lResults[12]  = nAntiEWK.getPropagatedError(*fitRes);
    lResults[13]  = nAntiQCD.getPropagatedError(*fitRes);
    if( iDoMu) lResults[14]  = qcd0.sigma->getError();
//.........这里部分代码省略.........