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

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 ComputeUpperLimit(RooAbsData *data, RooStats::ModelConfig *model, float &UpperLimit, float &signif, RooRealVar *mu, RooArgSet *nullParams,RooWorkspace *ws,REGION region,const char* tag) {

  bool StoreEverything=false; // activate if you want to store frames and all
  
  RooStats::ProfileLikelihoodCalculator *plc = new RooStats::ProfileLikelihoodCalculator(*data, *model);
  plc->SetParameters(*mu);
  plc->SetNullParameters(*nullParams);
  plc->SetTestSize(0.05);
  RooStats::LikelihoodInterval *interval = plc->GetInterval();

  bool ComputationSuccessful=false;
  UpperLimit = interval->UpperLimit(*mu,ComputationSuccessful);
  signif = 0.0; // plc->GetHypoTest()->Significance();   // deactivated significance (to make algorithm faster)

  if(!ComputationSuccessful) {
    cout << "There seems to have been a problem. Returned upper limit is " << UpperLimit << " but it will be set to -999" << endl;
    UpperLimit=-999;
    signif=-999;
  }

  if(StoreEverything) {
    // Store it all
    RooRealVar* minv = (RooRealVar*)model->GetObservables()->first();
    minv->setBins(static_cast<int>((minv->getMax()-minv->getMin())/5.));

    RooPlot* frameEE = minv->frame(RooFit::Title("ee sample"));
    frameEE->GetXaxis()->CenterTitle(1);
    frameEE->GetYaxis()->CenterTitle(1);
    
    RooPlot* frameMM = minv->frame(RooFit::Title("mm sample"));
    frameMM->GetXaxis()->CenterTitle(1);
    frameMM->GetYaxis()->CenterTitle(1);
    
    RooPlot* frameOF = minv->frame(RooFit::Title("OF sample"));
    frameOF->GetXaxis()->CenterTitle(1);
    frameOF->GetYaxis()->CenterTitle(1);
    
    data->plotOn(frameMM,RooFit::Cut("catCentral==catCentral::MMCentral"));
    model->GetPdf()->plotOn(frameMM,RooFit::Slice(*ws->cat("catCentral"), "MMCentral"),RooFit::ProjWData(*data));
    
    data->plotOn(frameEE,RooFit::Cut("catCentral==catCentral::EECentral"));
    model->GetPdf()->plotOn(frameEE,RooFit::Slice(*ws->cat("catCentral"), "EECentral"),RooFit::ProjWData(*data));
    
    data->plotOn(frameOF,RooFit::Cut("catCentral==catCentral::OFOSCentral"));
    model->GetPdf()->plotOn(frameOF,RooFit::Slice(*ws->cat("catCentral"), "OFOSCentral"),RooFit::ProjWData(*data));
    
    TFile *fout = new TFile("fout.root","UPDATE");
    frameMM->Write(Concatenate(Concatenate(data->GetName(),"_MM"),tag),TObject::kOverwrite);
    frameEE->Write(Concatenate(Concatenate(data->GetName(),"_EE"),tag),TObject::kOverwrite);
    frameOF->Write(Concatenate(Concatenate(data->GetName(),"_OF"),tag),TObject::kOverwrite);
    fout->Close();
  }

  delete plc;
  plc=0;
}

int main(int argc, char* argv[]){
 
  string fileName;
  string fileNameZee;  
  string functionName;
  string fileNameout;
  int ncats;
  int jcats;
  int bins; 
  string outfilename;
  bool is2011=false;
  bool useDoubleCB=false;  
  bool verbose=false;
  int mhLow;
  int mhHigh;



  po::options_description desc("Allowed options");
  desc.add_options()
    ("help,h",                                                                                  "Show help")
    ("infilename,i", po::value<string>(&fileName),                                              "In file name")
    ("infilenameZee,I", po::value<string>(&fileNameZee),                                              "In file name Zee")   
    ("function,f", po::value<string>(&functionName),                                              "Function to use")
    ("Outfilename,o", po::value<string>(&fileNameout),                                              "Out file name")
    ("ncats,c", po::value<int>(&ncats)->default_value(5),                                       "Number of categories")
    ("jcats,j", po::value<int>(&jcats)->default_value(0),                                       "Start number of categories")
    ("mhLow,L", po::value<int>(&mhLow)->default_value(75),                                                                                                                 "Starting point for scan") 
    ("mhHigh,H", po::value<int>(&mhHigh)->default_value(120),                                                                                                               "End point for scan") 
    ("bins,B", po::value<int>(&bins)->default_value(180),                                                                                                                 "Bins for the dataset") 
    ("is2011",                                                                                  "Run 2011 config")
    ("useDoubleCB",                                                                                  "use double crystal ball function")   
    ("verbose,v",                                                                               "Run with more output")
  ;
  po::variables_map vm;
  po::store(po::parse_command_line(argc,argv,desc),vm);
  po::notify(vm);
  if (vm.count("help")) { cout << desc << endl; exit(1); }
  if (vm.count("is2011")) is2011=true;
  if (vm.count("useDoubleCB"))  useDoubleCB=true;   
  if (vm.count("verbose")) verbose=true;

  if (!verbose) {
    RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
    RooMsgService::instance().setSilentMode(true);
  }


  TFile *outputfile;
  //RooWorkspace *outputws = new RooWorkspace("cms_hgg_workspace");  
  RooWorkspace *outputws;  
  outputfile = new TFile(fileNameout.c_str(),"RECREATE");

  
  
  TFile *inFile = TFile::Open(fileName.c_str());
  RooWorkspace *inWS = (RooWorkspace*)inFile->Get("cms_hgg_workspace");
  outputws = (RooWorkspace*)inWS->Clone("cms_hgg_workspace");   
  



  vector<string> functionClasses;
  functionClasses.push_back("Chebychev");
  functionClasses.push_back("Bernstein");
  functionClasses.push_back("Exponential");
  functionClasses.push_back("PowerLaw");
  functionClasses.push_back("Laurent");
  
  map<string,string> namingMap;
  namingMap.insert(pair<string,string>("Bernstein","pol"));  
  namingMap.insert(pair<string,string>("Exponential","exp"));
  namingMap.insert(pair<string,string>("PowerLaw","pow"));
  namingMap.insert(pair<string,string>("Laurent","lau"));
  
  vector<pair<pair<string,int> ,pair<pair<int,int>, pair<float,float> > > >  fabChoice;
  int sqrts; string ext;
  if (is2011) {
    sqrts = 7;
    ext = "7TeV";
  }
  else {
    sqrts = 8;
    ext = "8TeV";
    fabChoice.push_back(pair<pair<string,int>, pair<pair<int,int>, pair<float,float> > >(make_pair("Bernstein",-3),make_pair(make_pair(5,1), make_pair(-11.0,11.0)))); //0 
    fabChoice.push_back(pair<pair<string,int>, pair<pair<int,int>, pair<float,float> > >(make_pair("Bernstein",-3),make_pair(make_pair(5,1), make_pair(-11.0,11.0)))); //1
    fabChoice.push_back(pair<pair<string,int>, pair<pair<int,int>, pair<float,float> > >(make_pair("Chebychev",-3),make_pair(make_pair(5,1), make_pair(-11.0,11.0)))); //2
    fabChoice.push_back(pair<pair<string,int>, pair<pair<int,int>, pair<float,float> > >(make_pair("Bernstein",-3),make_pair(make_pair(5,1), make_pair(-11.0,11.0)))); //3
  }

  // store results here

  PdfModelBuilderFAN pdfsModel;
  RooRealVar *mass = (RooRealVar*)inWS->var("CMS_hgg_mass");
  mass->setRange(mhLow,mhHigh); 
  pdfsModel.setObsVar(mass);
  mass->setBins(bins); 
  ofstream outfile("Zee_Yield.log");  

cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ Initialization Done +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl; 
//.........这里部分代码省略.........

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 plotPdf_7D_XWW(double mH = 125, bool draw=true) {

    gROOT->ProcessLine(".L tdrstyle.C");
    setTDRStyle();
    TGaxis::SetMaxDigits(3);
    gROOT->ForceStyle();
    
    // Declaration of the PDFs to use
    gROOT->ProcessLine(".L  PDFs/RooSpinTwo_7D.cxx+");

    // W/Z mass and decay width constants
    double mV = 80.399;
    double gamV = 2.085;
    bool offshell = false;
    if ( mH < 2 * mV ) offshell = true;

    
    // for the pole mass and decay width of W 
    RooRealVar* mX = new RooRealVar("mX","mX", mH);
    RooRealVar* mW = new RooRealVar("mW","mW", mV);
    RooRealVar* gamW = new RooRealVar("gamW","gamW",gamV);

    //
    // Observables (7D)
    // 
    RooRealVar* wplusmass = new RooRealVar("wplusmass","m(W+)",mV,1e-09,120);
    wplusmass->setBins(50);
    RooRealVar* wminusmass = new RooRealVar("wminusmass","m(W-)",mV,1e-09,120);
    wminusmass->setBins(50);
    RooRealVar* hs = new RooRealVar("costhetastar","cos#theta*",-1,1);
    hs->setBins(20);
    RooRealVar* Phi1 = new RooRealVar("phistar1","#Phi_{1}",-TMath::Pi(),TMath::Pi());
    Phi1->setBins(20);
    RooRealVar* h1 = new RooRealVar("costheta1","cos#theta_{1}",-1,1);
    h1->setBins(20);
    RooRealVar* h2 = new RooRealVar("costheta2","cos#theta_{2}",-1,1);
    h2->setBins(20);
    RooRealVar* Phi = new RooRealVar("phi","#Phi",-TMath::Pi(),TMath::Pi());
    Phi->setBins(20);
    
    //
    // coupling constants for 2m+
    // See equation 5,6,7 in PRD 91, 075022
    //
    double s = (mH*mH-2*mV*mV)/2.;
    double c1 = 2*(1+mV*mV/s);
    c1 = c1 * 2.0; // scale up to be consistent with the generator
    // std::cout << "c1 = " << c1 << "\n"; 

    RooRealVar* c1Val = new RooRealVar("c1Val", "c1Val", c1);
    RooRealVar* c2Val = new RooRealVar("c2Val", "c2Val", -0.5);
    RooRealVar* c3Val = new RooRealVar("c3Val", "c3Val", 0.);
    RooRealVar* c4Val = new RooRealVar("c4Val", "c4Val", -1.);
    RooRealVar* c5Val = new RooRealVar("c5Val", "c5Val", 0.);
    RooRealVar* c6Val = new RooRealVar("c6Val", "c6Val", 0.);
    RooRealVar* c7Val = new RooRealVar("c7Val", "c7Val", 0.);
    
    // 
    // Alternative definition in terms of g1->g10
    // 
    RooRealVar* useGTerm = new RooRealVar("useGTerm", "useGTerm",1.); // set to 1 if using g couplings
    RooRealVar* g1Val = new RooRealVar("g1Val", "g1Val", 1);
    RooRealVar* g2Val = new RooRealVar("g2Val", "g2Val", 0.);
    RooRealVar* g3Val = new RooRealVar("g3Val", "g3Val", 0.);
    RooRealVar* g4Val = new RooRealVar("g4Val", "g4Val", 0.);
    RooRealVar* g5Val = new RooRealVar("g5Val", "g5Val", 1.);
    RooRealVar* g6Val = new RooRealVar("g6Val", "g6Val", 0.);
    RooRealVar* g7Val = new RooRealVar("g7Val", "g7Val", 0.);
    RooRealVar* g8Val = new RooRealVar("g8Val", "g8Val", 0.);
    RooRealVar* g9Val = new RooRealVar("g9Val", "g9Val", 0.);
    RooRealVar* g10Val = new RooRealVar("g10Val", "g10Val", 0.);

    // related to the gg/qq productions 
    RooRealVar* fz1Val = new RooRealVar("fz1Val", "fz1Val", 0);
    RooRealVar* fz2Val = new RooRealVar("fz2Val", "fz2Val", 1.0);

    // Even more parameters, do not have to touch, based on Z couplings
    RooRealVar* R1Val = new RooRealVar("R1Val","R1Val",1);
    RooRealVar* R2Val = new RooRealVar("R2Val","R2Val",1);
    
      
    // PDF definition SM Higgs (JP = 2+)
    RooSpinTwo_7D *myPDF;
    if ( offshell )
      myPDF = new RooSpinTwo_7D("myPDF","myPDF", *mX, *wplusmass, *wminusmass, *hs, *h1,*h2, *Phi, *Phi1, 
				  *c1Val, *c2Val, *c3Val, *c4Val, *c5Val, *c6Val, *c7Val, 
				  *useGTerm, *g1Val, *g2Val, *g3Val, *g4Val, *g5Val, *g6Val, *g7Val, *g8Val, *g9Val, *g10Val,
				  *fz1Val, *fz2Val, *R1Val, *R2Val, *mW, *gamW);
    else 
      myPDF = new RooSpinTwo_7D("myPDF","myPDF", *mX, *mW, *mW, *hs, *h1,*h2, *Phi, *Phi1, 
				  *c1Val, *c2Val, *c3Val, *c4Val, *c5Val, *c6Val, *c7Val, 
				  *useGTerm, *g1Val, *g2Val, *g3Val, *g4Val, *g5Val, *g6Val, *g7Val, *g8Val, *g9Val, *g10Val,
				  *fz1Val, *fz2Val, *R1Val, *R2Val, *mW, *gamW);
    // dataset for (JP = 2+)
    TString fileName;
    if ( useGTerm->getVal() > 0.) {
      fileName = Form("TWW_2mplus_%.0f_JHU.root", mH);
    }
    else {
      fileName = Form("TWW_%.0f_JHU_GenFromC.root", mH);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  TCut Events1_4 = "(Entry$%4==1)";
  TCut Events2_4 = "(Entry$%4==2)";
  TCut Events0_4 = "(Entry$%4==0)";

  TCut Events01_4 = "(Entry$%4<2)";
  TCut Events23_4 = "(Entry$%4>1)";

  TCut EventsTest = "(Entry$%2==1)";

  //weightvar.SetTitle(EventsTest*selcut);
  weightvar.SetTitle(selcut);
/*
  if (doele) 
    weightvar.SetTitle(prescale100alt*selcut);
  else
    weightvar.SetTitle(selcut);
  */
  //make testing dataset
  RooDataSet *hdata = RooTreeConvert::CreateDataSet("hdata",dtree,vars,weightvar);   

  if (doele) 
    weightvar.SetTitle(prescale1000alt*selcut);
  else
    weightvar.SetTitle(prescale10alt*selcut);
  //make reduced testing dataset for integration over conditional variables
  RooDataSet *hdatasmall = RooTreeConvert::CreateDataSet("hdatasmall",dtree,vars,weightvar);     
    
  //retrieve full pdf from workspace
  RooAbsPdf *sigpdf = ws->pdf("sigpdf");
  
  //input variable corresponding to sceta
  RooRealVar *scEraw = ws->var("var_0");
  scEraw->setRange(1.,2.);
  scEraw->setBins(100);
//  RooRealVar *scetavar = ws->var("var_1");
//  RooRealVar *scphivar = ws->var("var_2");
  
 
  //regressed output functions
  RooAbsReal *sigmeanlim = ws->function("sigmeanlim");
  RooAbsReal *sigwidthlim = ws->function("sigwidthlim");
  RooAbsReal *signlim = ws->function("signlim");
  RooAbsReal *sign2lim = ws->function("sign2lim");

//  RooAbsReal *sigalphalim = ws->function("sigalphalim");
  //RooAbsReal *sigalpha2lim = ws->function("sigalpha2lim");


  //formula for corrected energy/true energy ( 1.0/(etrue/eraw) * regression mean)
  RooFormulaVar ecor("ecor","","1./(@0)*@1",RooArgList(*tgtvar,*sigmeanlim));
  RooRealVar *ecorvar = (RooRealVar*)hdata->addColumn(ecor);
  ecorvar->setRange(0.,2.);
  ecorvar->setBins(800);
  
  //formula for raw energy/true energy (1.0/(etrue/eraw))
  RooFormulaVar raw("raw","","1./@0",RooArgList(*tgtvar));
  RooRealVar *rawvar = (RooRealVar*)hdata->addColumn(raw);
  rawvar->setRange(0.,2.);
  rawvar->setBins(800);

  //clone data and add regression outputs for plotting
  RooDataSet *hdataclone = new RooDataSet(*hdata,"hdataclone");
  RooRealVar *meanvar = (RooRealVar*)hdataclone->addColumn(*sigmeanlim);
  RooRealVar *widthvar = (RooRealVar*)hdataclone->addColumn(*sigwidthlim);
  RooRealVar *nvar = (RooRealVar*)hdataclone->addColumn(*signlim);
  RooRealVar *n2var = (RooRealVar*)hdataclone->addColumn(*sign2lim);

void fitMass2(int iCharge,double iEtaMin,double iEtaMax,double iPhiMin,double iPhiMax,double &iRes,double &iShift,double &iResErr,double &iShiftErr) { 
  Prep();
  RooRealVar    lXVar  ("XVar","mass(GeV/c^{2})",60,60,120); lXVar.setBins(1000);
  RooRealVar    lSPar  ("SPar","SPar", 1.,0., 2.);
  RooFormulaVar lXShift("uparshift","@0*@1",RooArgList(lXVar,lSPar));
  TFile *lMCFile = new TFile("../Efficiency/Data/ZTP8_v2.root");
  TTree *lMCTree = (TTree*) lMCFile->FindObjectAny("WNtupleIdEffNT"); 
  TH1F *lMass = new TH1F("M","M",100,60,120); 
  int   lCharge = 0;   lMCTree->SetBranchAddress("charge",&lCharge);
  float lEta    = 0;   lMCTree->SetBranchAddress("eta"   ,&lEta);
  float lPhi    = 0;   lMCTree->SetBranchAddress("phi"   ,&lPhi);
  float lMt     = 0;   lMCTree->SetBranchAddress("mt"    ,&lMt);
  float lPt     = 0;   lMCTree->SetBranchAddress("pt"    ,&lPt);
  float lOPt    = 0;   lMCTree->SetBranchAddress("jetpt" ,&lOPt);
  float lOEta   = 0;   lMCTree->SetBranchAddress("jeteta",&lOEta);
  float lOPhi   = 0;   lMCTree->SetBranchAddress("jetphi",&lOPhi);
  float lTrkIso   = 0; lMCTree->SetBranchAddress("trkiso",&lTrkIso);
  float lEcalIso  = 0; lMCTree->SetBranchAddress("ecaliso",&lEcalIso);
  float lHcalIso  = 0; lMCTree->SetBranchAddress("hcaliso",&lHcalIso);
  float lChi2     = 0; lMCTree->SetBranchAddress("chi2"   ,&lChi2);
  unsigned int   lNHit     = 0; lMCTree->SetBranchAddress("nhit"   ,&lNHit);
  Muon  lMuon;         lMCTree->SetBranchAddress("muon"   ,&lMuon);
  double lVPt   = 0; double lEt = 0; double lPx = 0; double lPy = 0;
  for(int i0 = 0; i0 < lMCTree->GetEntries(); i0++) {
    lMCTree->GetEntry(i0);
    if(lMt < 60)          continue;
    if(lChi2 > 10)        continue;
    if(lNHit < 10)        continue;
    if(lMuon.NSeg < 2)    continue;
    if(lMuon.NPixel == 0) continue;
    if(lMuon.NValid == 0) continue;
    if(lMuon.Type   != 3) continue;
    if((lTrkIso+lEcalIso+lHcalIso)/lPt > 0.15) continue;
    //if(lCharge > 0 && iCharge < 0) continue;
    //if(lCharge < 0 && iCharge > 0) continue;
    //if(lEta < iPhiMin || lEta > iPhiMax) continue;
    //if(fabs(lEta) < 1.2) continue;
    //if(lEta < iEtaMin || lEta > iEtaMax) continue;
    //if(lPhi*lOPhi > 0 || lEta*lOEta > 0) continue;
    if(lPhi < iPhiMin || lPhi > iPhiMax) continue;
    lEt = lOPt + lPt; lPx = fabs(lPt)*cos(lPhi) + fabs(lOPt)*cos(lOPhi); lPy = fabs(lPt)*sin(lPhi) + fabs(lOPt)*sin(lOPhi);
    lVPt = sqrt(lPx*lPx + lPy*lPy);
    //if(lVPt < iPhiMin || lVPt > iPhiMax) continue;
    lMass->Fill(fabs(lMt));     //lMass->Fill(lPt);
  }
  RooDataHist *lMHist = new RooDataHist("M" ,"M" ,RooArgSet(lXVar),lMass);
  RooHistPdf  *lMPdf  = new RooHistPdf ("MH","MH",lXShift,lXVar,*lMHist,5); 
  RooRealVar l1Sigma("sigma1","sigma1",0.2,0.,15.);  //l1Sigma.setConstant(kTRUE);
  RooRealVar lR0Mean("xmean","xmean",0,-10,10);    lR0Mean.setConstant(kTRUE);
  RooRealVar lExp   ("exp"   ,"exp"   ,-0.006,-15,15.); //lExp.setConstant(kTRUE);
  RooRealVar lFrac  ("frac","frac"    ,0.9,0.,1);
  RooGaussian   lGaus1("gaus1","gaus1",lXVar,lR0Mean,l1Sigma);
  RooExponential lExpF("Exp","Exp"  ,lXVar,lExp);
  RooFFTConvPdf  lConv("Conv","Conv",lXVar,*lMPdf,lGaus1); //lConv.setBufferStrategy(RooFFTConvPdf::Flat);
  RooAddPdf      lGAdd("Add","Add"  ,lConv,lExpF,lFrac);
  RooDataSet *lData = new RooDataSet("crap","crap",RooArgSet(lXVar)); 
  TFile *lFile = new TFile("../Efficiency/Data/mTPNT8_v1.root");
  TTree *lTree = (TTree*) lFile->FindObjectAny("WNtupleIdEffNT");
  lTree->SetBranchAddress("charge",&lCharge);
  lTree->SetBranchAddress("eta"   ,&lEta);
  lTree->SetBranchAddress("phi"   ,&lPhi);
  lTree->SetBranchAddress("mt"    ,&lMt);
  lTree->SetBranchAddress("pt"    ,&lPt);
  lTree->SetBranchAddress("jetpt" ,&lOPt);
  lTree->SetBranchAddress("jeteta",&lOEta);
  lTree->SetBranchAddress("jetphi",&lOPhi);
  lTree->SetBranchAddress("trkiso",&lTrkIso);
  lTree->SetBranchAddress("ecaliso",&lEcalIso);
  lTree->SetBranchAddress("hcaliso",&lHcalIso);
  lTree->SetBranchAddress("chi2"   ,&lChi2);
  lTree->SetBranchAddress("nhit"   ,&lNHit);
  lTree->SetBranchAddress("muon"   ,&lMuon);
  for(int i0 = 0; i0 < lTree->GetEntries();i0++) { 
    lTree->GetEntry(i0);
    if(lMt < 60)                   continue;
    if(lCharge > 0 && iCharge < 0) continue;
    if(lCharge < 0 && iCharge > 0) continue;
    if(lChi2 > 10) continue;
    if(lNHit < 10) continue;
    if(lMuon.NSeg < 2) continue;
    if(lMuon.NPixel == 0) continue;
    if(lMuon.NValid == 0) continue;
    if(lMuon.Type   != 3) continue;
    if((lTrkIso+lEcalIso+lHcalIso)/lPt > 0.15) continue;
    //if(fabs(lEta) < 1.2) continue;
    //if(lPhi*lOPhi > 0 || lEta*lOEta > 0) continue;
    //if(lEta < iPhiMin || lEta > iPhiMax) continue;
    if(lPhi < iPhiMin || lPhi > iPhiMax) continue;
    lEt = lOPt + lPt; lPx = fabs(lPt)*cos(lPhi) + fabs(lOPt)*cos(lOPhi); lPy = fabs(lPt)*sin(lPhi) + fabs(lOPt)*sin(lOPhi);
    lVPt = sqrt(lPx*lPx + lPy*lPy);
    //if(lVPt < iPhiMin || lVPt > iPhiMax) continue;
    lXVar.setVal(fabs(lMt));
    if(lCharge > 0) lXVar.setVal(correct(lMt,lPhi,lOPhi,lEta,lOEta));
    if(lCharge < 0) lXVar.setVal(correct(lMt,lOPhi,lPhi,lOEta,lEta));
    lData->add(RooArgSet(lXVar));
  }
  /*
  TFile *lQFile = new TFile("");
  TTree *lQTree = (TTree*) lQFile->FindObjectAny("WNtupleIdEffNT");
  lQTree->SetBranchAddress("mt"    ,&lMt);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  vector<double> muFab;
  vector<double> muPaul;
  vector<double> muChi2;
  vector<double> muAIC;
  vector<double> muFabErrLow;
  vector<double> muPaulErrLow;
  vector<double> muChi2ErrLow;
  vector<double> muAICErrLow;
  vector<double> muFabErrHigh;
  vector<double> muPaulErrHigh;
  vector<double> muChi2ErrHigh;
  vector<double> muAICErrHigh;
  muTree->Branch("jobn",&jobn);
  muTree->Branch("toyn",&toyn);
  muTree->Branch("truthModel",&truthModel);
  muTree->Branch("muFab",&muFab);
  muTree->Branch("muPaul",&muPaul);
  muTree->Branch("muChi2",&muChi2);
  muTree->Branch("muAIC",&muAIC);
  muTree->Branch("muFabErrLow",&muFabErrLow);
  muTree->Branch("muPaulErrLow",&muPaulErrLow);
  muTree->Branch("muChi2ErrLow",&muChi2ErrLow);
  muTree->Branch("muAICErrLow",&muAICErrLow);
  muTree->Branch("muFabErrHigh",&muFabErrHigh);
  muTree->Branch("muPaulErrHigh",&muPaulErrHigh);
  muTree->Branch("muChi2ErrHigh",&muChi2ErrHigh);
  muTree->Branch("muAICErrHigh",&muAICErrHigh);
  
  //TH1F *muDistFab = new TH1F("muDistFab","muDistFab",int(20*(mu_high-mu_low)),mu_low,mu_high);
  //TH1F *muDistPaul = new TH1F("muDistPaul","muDistPaul",int(20*(mu_high-mu_low)),mu_low,mu_high);
  //TH1F *muDistChi2 = new TH1F("muDistChi2","muDistChi2",int(20*(mu_high-mu_low)),mu_low,mu_high);
  //TH1F *muDistAIC = new TH1F("muDistAIC","muDistAIC",int(20*(mu_high-mu_low)),mu_low,mu_high);
  
  mass->setBins(320);
  RooDataSet *data = (RooDataSet*)bkgWS->data(Form("data_mass_cat%d",cat));
  //RooDataSet *data = (RooDataSet*)bkgWS->data(Form("data_cat%d_7TeV",cat));
  RooDataHist *dataBinned = new RooDataHist(Form("roohist_data_mass_cat%d",cat),Form("roohist_data_mass_cat%d",cat),RooArgSet(*mass),*data);
  RooDataSet *sigMC = (RooDataSet*)sigWS->data(Form("sig_ggh_mass_m%d_cat%d",expectSignalMass,cat));
  RooDataSet *sigMC_vbf = (RooDataSet*)sigWS->data(Form("sig_wzh_mass_m%d_cat%d",expectSignalMass,cat));
  RooDataSet *sigMC_wzh = (RooDataSet*)sigWS->data(Form("sig_vbf_mass_m%d_cat%d",expectSignalMass,cat));
  RooDataSet *sigMC_tth = (RooDataSet*)sigWS->data(Form("sig_tth_mass_m%d_cat%d",expectSignalMass,cat));
  sigMC->append(*sigMC_vbf);
  sigMC->append(*sigMC_wzh);
  sigMC->append(*sigMC_tth);
  //RooExtendPdf *ggh_pdf = (RooExtendPdf*)sigWS->pdf(Form("sigpdfsmrel_cat%d_7TeV_ggh",cat));
  //RooExtendPdf *vbf_pdf = (RooExtendPdf*)sigWS->pdf(Form("sigpdfsmrel_cat%d_7TeV_vbf",cat));
  //RooExtendPdf *wzh_pdf = (RooExtendPdf*)sigWS->pdf(Form("sigpdfsmrel_cat%d_7TeV_wzh",cat));
  //RooExtendPdf *tth_pdf = (RooExtendPdf*)sigWS->pdf(Form("sigpdfsmrel_cat%d_7TeV_tth",cat));
  //RooAbsPdf *sigPdf = new RooAddPdf(Form("sigpdfsmrel_cat%d_7TeV",cat),Form("sigpdfsmrel_cat%d_7TeV",cat),RooArgList(*ggh_pdf,*vbf_pdf,*wzh_pdf,*tth_pdf));
  
  if (!dataBinned || !sigMC){
    cerr << "ERROR -- one of data or signal is NULL" << endl;
    exit(1);
  }
  
  // set of truth models to throw toys from
  PdfModelBuilder toysModel;
  toysModel.setObsVar(mass);
  toysModel.setSignalModifier(mu);
  // add truth pdfs from config datfile these need to be cached
  // to throw a toy from the SB fit make sure that the cache happens at makeSBPdfs
  for (vector<pair<int,pair<string,string> > >::iterator it=toysMap.begin(); it!=toysMap.end(); it++){
    if (it->first==-1) { // this is a hyrbid toy
      throwHybridToys=true;
      vector<string> temp;
      split(temp,it->second.first,boost::is_any_of(","));

void makeDatacard(double mh, double massLow, double massHigh, double merrHigh, int ch, std::string cat, std::map<std::string, std::string> file, bool useModZ, bool doMassErr) {

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

    if (doMassErr && cat != "BB" && cat != "XX" && cat != "YY") {
        std::cout << "When using event-by-event mass errors, set cat = BB or XX or YY" << std::endl; 
        return;
    }

    /* Setting up the strings */
    std::string chstr = getChannelName(ch);

    stringstream mh_ss;
    mh_ss << mh;
    
    std::cout << "Creating datacard for " << mh_ss.str() << " GeV mass point, channel " << chstr << " and category " << cat << " ... " << std::endl;
   
    std::stringstream card_name_ss;
    card_name_ss << "card_";
    if      (!doMassErr) card_name_ss << "1D_";
    else                 card_name_ss << "2D_merr_";
    card_name_ss << "m" << mh_ss.str() << "_";
    card_name_ss << chstr << "_" << cat;
    std::string card_name = card_name_ss.str();

    std::string workspace = card_name+"_workspace.root";

    /* Higgs mass and dimuon mass variables */

    const char* massvarstr  = "CMS_hmumu_mass";
    const char* merrvarstr  = "CMS_hmumu_merr";
    const char* scalevarstr = "CMS_hmumu_scale";
    const char* resvarstr   = "CMS_hmumu_res";
    const char* mhvarstr    = "MH";

    RooRealVar rmh  ("MH"       , "MH"         , mh);
    RooRealVar m2mu (massvarstr , "Dimuon mass", mh  , massLow, massHigh, "GeV/c^{2}");
    RooRealVar merr (merrvarstr , "Mass Error" , 0.0 , 0.0    , merrHigh, ""         );
    RooRealVar scale(scalevarstr, "Scale unc. ", 0.0 , 0.0    , 1.0     , "GeV/c^{2}");
    RooRealVar res  (resvarstr  , "RFes. unc. ", 0.0 , 0.0    , 1.0);
    m2mu.setBins(200);   
    merr.setBins(200);   

    /* RooDataSet of the observed data */

    RooDataSet* data_obs = getDataSet(file["dat"].c_str(), false, ch, cat, massLow, massHigh, merrHigh, "data_obs", massvarstr, merrvarstr, false);
    
    /* Extract shape parameters */

    std::cout << "Extracting the signal fit parameters" << std::endl;
    std::map<std::string, double> sparams = doFit(file["ggH"],  true,   6, cat, 125.0, int(massHigh-massLow)*2, massLow,    140.0, false,  true, useModZ, false);
    std::cout << "Extracting the background fit parameters" << std::endl;
    std::map<std::string, double> bparams = doFit(file["dat"], false,  ch, cat, 125.0, int(massHigh-massLow) , massLow, massHigh, false, false, useModZ, true );

    std::string ecat = "XX";
    if (cat == "BB") ecat = "BB";
    if (cat == "YY") ecat = "YY";

    std::map<std::string, double> esparams;
    std::map<std::string, double> ebparams;
    if (doMassErr) {
        if (cat == "BB") {
            std::cout << "Extracting the signal mass error fit parameters" << std::endl;
            esparams = doEbEFit(file["ggH"], true, ch, ecat, 140, 0.006, 0.013, massLow, massHigh,  true,  true,  true);
            std::cout << "Extracting the background mass error fit parameters" << std::endl;
            ebparams = doEbEFit(file["DY"] , true, ch, ecat, 140, 0.006, 0.013, massLow, massHigh,  true,  true, false);
        }
        else {
            std::cout << "Extracting the signal mass error fit parameters" << std::endl;
            esparams = doEbEFit(file["ggH"], true, ch, ecat, 240, 0.005, 0.035, massLow, massHigh,  true,  true,  true);
            std::cout << "Extracting the background mass error fit parameters" << std::endl;
            ebparams = doEbEFit(file["DY"] , true, ch, ecat, 240, 0.005, 0.035, massLow, massHigh,  true,  true, false);
        }
    }

    /* Compute yields */

    double bkg_yield = computeYield(data_obs, massvarstr, false,  true);
    bkg_yield *= bparams["bkgsf"];
    std::cout << "Computing the expected background yield from the side-bands : " << bkg_yield << std::endl;

    RooDataSet* sig_gH_dset = getDataSet(file["ggH"], true, ch, cat, massLow, massHigh, merrHigh, "dset_gH", massvarstr, merrvarstr, false);
    RooDataSet* sig_qH_dset = getDataSet(file["qqH"], true, ch, cat, massLow, massHigh, merrHigh, "dset_qH", massvarstr, merrvarstr, false);
    RooDataSet* sig_PH_dset = getDataSet(file["WPH"], true, ch, cat, massLow, massHigh, merrHigh, "dset_PH", massvarstr, merrvarstr, false);
    RooDataSet* sig_MH_dset = getDataSet(file["WMH"], true, ch, cat, massLow, massHigh, merrHigh, "dset_MH", massvarstr, merrvarstr, false);
    RooDataSet* sig_ZH_dset = getDataSet(file["ZH" ], true, ch, cat, massLow, massHigh, merrHigh, "dset_ZH", massvarstr, merrvarstr, false);

    double sig_gH_yield  = computeYield(sig_gH_dset, massvarstr,  true, false);
    double sig_qH_yield  = computeYield(sig_qH_dset, massvarstr,  true, false);
    double sig_WH_yield  = computeYield(sig_PH_dset, massvarstr,  true, false);
           sig_WH_yield += computeYield(sig_MH_dset, massvarstr,  true, false);
    double sig_ZH_yield  = computeYield(sig_ZH_dset, massvarstr,  true, false);
    double sig_tH_yield  = 1e-5;

    delete sig_gH_dset;
    delete sig_qH_dset;
    delete sig_PH_dset;
    delete sig_MH_dset;
    delete sig_ZH_dset;
//.........这里部分代码省略.........

void performFit(string inputDir, string inputParameterFile, string label,
                string PassInputDataFilename, string FailInputDataFilename, 
                string PassSignalTemplateHistName, string FailSignalTemplateHistName)
{

  gBenchmark->Start("fitZCat");

  //--------------------------------------------------------------------------------------------------------------
  // Settings 
  //==============================================================================================================
  
  const Double_t mlow  = 60;
  const Double_t mhigh = 120;
  const Int_t    nbins = 24;

  TString effType = inputDir;

  // The fit variable - lepton invariant mass
  RooRealVar* rooMass_ = new RooRealVar("Mass","m_{ee}",mlow, mhigh, "GeV/c^{2}");
  RooRealVar Mass = *rooMass_;
  Mass.setBins(nbins);

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


  RooDataSet *dataPass  = RooDataSet::read((inputDir+PassInputDataFilename).c_str(),RooArgList(Mass));
  RooDataSet *dataFail  = RooDataSet::read((inputDir+FailInputDataFilename).c_str(),RooArgList(Mass));

  RooDataSet *dataCombined = new RooDataSet("dataCombined","dataCombined", RooArgList(Mass), RooFit::Index(sample), 
                                            RooFit::Import("Pass",*dataPass),
                                            RooFit::Import("Fail",*dataFail));



  //*********************************************************************************************
  //Define Free Parameters
  //*********************************************************************************************
  RooRealVar* ParNumSignal =  LoadParameters(inputParameterFile, label,"ParNumSignal");
  RooRealVar* ParNumBkgPass =  LoadParameters(inputParameterFile, label,"ParNumBkgPass");
  RooRealVar* ParNumBkgFail =  LoadParameters(inputParameterFile, label, "ParNumBkgFail");
  RooRealVar* ParEfficiency  = LoadParameters(inputParameterFile, label, "ParEfficiency");
  RooRealVar* ParPassBackgroundExpCoefficient = LoadParameters(inputParameterFile, label, "ParPassBackgroundExpCoefficient");
  RooRealVar* ParFailBackgroundExpCoefficient = LoadParameters(inputParameterFile, label, "ParFailBackgroundExpCoefficient");
  RooRealVar* ParPassSignalMassShift = LoadParameters(inputParameterFile, label, "ParPassSignalMassShift");
  RooRealVar* ParFailSignalMassShift = LoadParameters(inputParameterFile, label, "ParFailSignalMassShift");
  RooRealVar* ParPassSignalResolution = LoadParameters(inputParameterFile, label, "ParPassSignalResolution");
  RooRealVar* ParFailSignalResolution = LoadParameters(inputParameterFile, label, "ParFailSignalResolution");



// new RooRealVar  ("ParPassSignalMassShift","ParPassSignalMassShift",-2.6079e-02,-10.0, 10.0);   //ParPassSignalMassShift->setConstant(kTRUE);  
//   RooRealVar* ParFailSignalMassShift = new RooRealVar  ("ParFailSignalMassShift","ParFailSignalMassShift",7.2230e-01,-10.0, 10.0);   //ParFailSignalMassShift->setConstant(kTRUE);  
//   RooRealVar* ParPassSignalResolution = new RooRealVar ("ParPassSignalResolution","ParPassSignalResolution",6.9723e-01,0.0, 10.0);     ParPassSignalResolution->setConstant(kTRUE);  
//   RooRealVar* ParFailSignalResolution = new RooRealVar ("ParFailSignalResolution","ParFailSignalResolution",1.6412e+00,0.0, 10.0);     ParFailSignalResolution->setConstant(kTRUE);  

  //*********************************************************************************************
  //
  //Load Signal PDFs
  //
  //*********************************************************************************************

  TFile *Zeelineshape_file =  new TFile("res/photonEfffromZee.dflag1.eT1.2.gT40.mt15.root", "READ");
  TH1* histTemplatePass = (TH1D*) Zeelineshape_file->Get(PassSignalTemplateHistName.c_str());
  TH1* histTemplateFail = (TH1D*) Zeelineshape_file->Get(FailSignalTemplateHistName.c_str());

  //Introduce mass shift coordinate transformation 
//   RooFormulaVar PassShiftedMass("PassShiftedMass","@[email protected]",RooArgList(Mass,*ParPassSignalMassShift));
//   RooFormulaVar FailShiftedMass("FailShiftedMass","@[email protected]",RooArgList(Mass,*ParFailSignalMassShift));

  RooGaussian  *PassSignalResolutionFunction =  new RooGaussian("PassSignalResolutionFunction","PassSignalResolutionFunction",Mass,*ParPassSignalMassShift,*ParPassSignalResolution);
  RooGaussian  *FailSignalResolutionFunction =  new RooGaussian("FailSignalResolutionFunction","FailSignalResolutionFunction",Mass,*ParFailSignalMassShift,*ParFailSignalResolution);


  RooDataHist* dataHistPass = new RooDataHist("dataHistPass","dataHistPass", RooArgSet(Mass), histTemplatePass);
  RooDataHist* dataHistFail = new RooDataHist("dataHistFail","dataHistFail", RooArgSet(Mass), histTemplateFail);
  RooHistPdf* signalShapePassTemplatePdf = new RooHistPdf("signalShapePassTemplatePdf", "signalShapePassTemplatePdf", Mass, *dataHistPass, 1);
  RooHistPdf* signalShapeFailTemplatePdf = new RooHistPdf("signalShapeFailTemplatePdf", "signalShapeFailTemplatePdf", Mass, *dataHistFail, 1);

  RooFFTConvPdf* signalShapePassPdf = new RooFFTConvPdf("signalShapePassPdf","signalShapePassPdf"  , Mass, *signalShapePassTemplatePdf,*PassSignalResolutionFunction,2);
  RooFFTConvPdf* signalShapeFailPdf = new RooFFTConvPdf("signalShapeFailPdf","signalShapeFailPdf"  , Mass, *signalShapeFailTemplatePdf,*FailSignalResolutionFunction,2);


  // Now define some efficiency/yield variables  
  RooFormulaVar* NumSignalPass = new RooFormulaVar("NumSignalPass", "ParEfficiency*ParNumSignal", RooArgList(*ParEfficiency,*ParNumSignal));
  RooFormulaVar* NumSignalFail = new RooFormulaVar("NumSignalFail", "(1.0-ParEfficiency)*ParNumSignal", RooArgList(*ParEfficiency,*ParNumSignal));


  //*********************************************************************************************
  //
  // Create Background PDFs
  //
  //*********************************************************************************************
  RooExponential* bkgPassPdf = new RooExponential("bkgPassPdf","bkgPassPdf",Mass, *ParPassBackgroundExpCoefficient);
  RooExponential* bkgFailPdf = new RooExponential("bkgFailPdf","bkgFailPdf",Mass, *ParFailBackgroundExpCoefficient);


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

void fitPtOverMCJLST(int mass = 125, int LHCsqrts = 7, int whichtype = 1, 
		     bool correctErrors = false, /* string changeParName = "", */
		     bool showErrorPDFs = false, string systString = "Default")

// whichtype
// 0 - gg Signal
// 1 - VBF Signal
// 2 - ZZ
// 3 - ZX
// 4 - ggZZ
// 5 - WH
// 6 - ZH
// 7 - ttH

{

  string changeParName = "";
  if (systString == "Default") changeParName = "up";

  string nameSample[8] = {"gg","vbf","zz","zx","ggzz","wh","zh","tth"};
  float maxType[8] = {2.4,3.2,1.6,1.6,1.6,3.2,3.2,3.2};   
  float rebinType[8] = {1,2,1,1,4,10,10,40};
  
  for (int t = 0; t < 8; t++) {
    if (mass > 150) maxType[t] = int(117.90*maxType[t]/sqrt(mass-10.91))/10.;
  }

  char fileToOpen[200];
  sprintf(fileToOpen,"selRootFiles/PToverM_%s%d_SEL_%dTeV.root",nameSample[whichtype].c_str(),mass,LHCsqrts);
  // if (whichtype == 3) sprintf(fileToOpen,"PTOVERM_%s_SEL_allTeV.root",nameSample[whichtype].c_str());

  RooRealVar* ptoverm = new RooRealVar("ptoverm","p_{T}/M^{4l}",0.,maxType[whichtype],"GeV/c");
 
  TFile input(fileToOpen);
  // if (systString == "Mass" || systString == "Mela") {
  //  sprintf(fileToOpen,"ptovermH_%sUp",systString.c_str());
  // } else {
  sprintf(fileToOpen,"ptovermH_%s",systString.c_str());
  //}
  TH1F* ptovermH = (TH1F*)input.Get(fileToOpen);
  
  if (rebinType[whichtype] > 1) ptovermH->Rebin(rebinType[whichtype]);
  if (maxType[whichtype] < ptovermH->GetBinLowEdge(ptovermH->GetNbinsX() + 1) - ptovermH->GetBinWidth(1)) {
    int theBin = ptovermH->FindBin(maxType[whichtype]);
    ptovermH->GetXaxis()->SetRange(1,theBin-1);
  }

  gROOT->ProcessLine(".L mytdrstyle.C");
  gROOT->ProcessLine("setTDRStyle()");
  
  // cout << endl << "Signal " << endl;   
  ptoverm->setBins(ptovermH->GetNbinsX());

  RooDataHist* rdh = new RooDataHist("rdh","Some dataset",RooArgList(*ptoverm),Import(*ptovermH,kFALSE));
 
  // fit definitions
  // RooWorkspace *ws = new RooWorkspace("ws");

  RooRealVar m("m","emme", 1.,0.01, 40.);
  RooRealVar n("n","enne", 0.93, 0.05, 15.);
  RooRealVar n2("n2","enne2", 0.75, 0.5, 15.);
  RooRealVar bb("bb","bibi",0.02, 0.000005, 20.0);
  RooRealVar T("T","tti",0.2,0.00000005,1.);
  RooRealVar bb2("bb2","bibi2",0.02, 0.0005, 10.0);
  RooRealVar fexp("fexp","f_exp",0.02, 0.0, 1.0);
  if (whichtype == 0) {
    if (LHCsqrts == 8) {
      m.setVal(3.319);   // m.setConstant(kTRUE);
      n.setVal(0.7606);    if (systString != "Default" || mass != 125) n.setConstant(kTRUE); 
      n2.setVal(0.8061);   n2.setConstant(kTRUE);
      bb.setVal(3.728);   // bb.setConstant(kTRUE);
      T.setVal(0.00333);   // T.setConstant(kTRUE);
      bb2.setVal(1.7172);    // bb2.setConstant(kTRUE);
      fexp.setVal(0.002144);   if (systString != "Default" || mass != 125) fexp.setConstant(kTRUE);
    } else {
      m.setVal(0.061);   // m.setConstant(kTRUE);
      n.setVal(1.6141);   if (systString == "Resummation" || systString == "TopMass")  n.setConstant(kTRUE);
      n2.setVal(1.3294);   n2.setConstant(kTRUE);
      bb.setVal(4.2761);   // bb.setConstant(kTRUE);
      T.setVal(0.0361);   // T.setConstant(kTRUE);
      bb2.setVal(1.6643);   bb2.setConstant(kTRUE);
      fexp.setVal(0.0004);   // fexp.setConstant(kTRUE);
    }
  } else if (whichtype == 1) {
    m.setVal(1.006);   // m.setConstant(kTRUE);
    n.setVal(10.939);   n.setConstant(kTRUE);
    n2.setVal(1.1448);   n2.setConstant(kTRUE);
    bb.setVal(0.02048);   bb.setConstant(kTRUE);
    T.setVal(0.16115);   if (systString.find("Mela") != string::npos) T.setConstant(kTRUE); // T.setConstant(kTRUE);
    bb2.setVal(1.0024);   bb2.setConstant(kTRUE);
    fexp.setVal(0.005);   fexp.setConstant(kTRUE);
    if (mass > 300) {
      fexp.setVal(0.0);   fexp.setConstant(kFALSE);
    }
    if (mass > 500) {
      bb2.setVal(5.0);  //  bb2.setConstant(kFALSE);
    }
    if (mass > 500) {
      bb.setVal(15.0);  //  bb.setConstant(kFALSE);
    }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
				  sg_p8=new RooRealVar("sg_p8", "sg_p8", 0., 1.);
			  }
			  if (point[h]=="400")
			  {
				  rangeLo=300., rangeHi=600.;
				  sg_p0=new RooRealVar("sg_p0", "sg_p0", 370., 430.);
				  sg_p1=new RooRealVar("sg_p1", "sg_p1", 3., 15.);
				  sg_p2=new RooRealVar("sg_p2", "sg_p2", 370.,460.);
				  sg_p3=new RooRealVar("sg_p3", "sg_p3", 10., 100.);
				  sg_p8=new RooRealVar("sg_p8", "sg_p8", 0., 1.);
			  }
			  if (point[h]=="450")
			  {
				  rangeLo=300., rangeHi=600.;
				  sg_p0=new RooRealVar("sg_p0", "sg_p0", 420., 480.);
				  sg_p1=new RooRealVar("sg_p1", "sg_p1", 3., 15.);
				  sg_p2=new RooRealVar("sg_p2", "sg_p2", 410., 490.);
				  sg_p3=new RooRealVar("sg_p3", "sg_p3", 10., 100.);
				  sg_p8=new RooRealVar("sg_p8", "sg_p8", 0., 1.);
			  }
			  if (point[h]=="300")
			  {
				  rangeLo=SR_lo, rangeHi=550.;
				  sg_p0=new RooRealVar("sg_p0", "sg_p0", 290., 320.);
				  sg_p1=new RooRealVar("sg_p1", "sg_p1", 5., 9.);
				  sg_p2=new RooRealVar("sg_p2", "sg_p2", 250., 360.);
				  sg_p3=new RooRealVar("sg_p3", "sg_p3", 10., 130.);
				  sg_p8=new RooRealVar("sg_p8", "sg_p8", 0., 1.);
			  }


		   //	x=new RooRealVar("x", "m_{X} (GeV)", 250, 900);
		        int binning = (abs(SR_lo-SR_hi)/(rebin)); 	std::cout<<"  binning "<<binning<<std::endl;
   			x->setBins(binning);

		  RooGaussian signalCore("signalCore", "Signal Prediction", *x, *sg_p0, *sg_p1);
		  RooGaussian Cpolyn("Cpolyn", "Combinatoric", *x, *sg_p2, *sg_p3);
		  RooAddPdf signal_("signal_", "signal", RooArgList(signalCore, Cpolyn), *sg_p8);
		  RooDataHist signalHistogram("signalHistogram", "Signal Histogram", RooArgList(*x), h_mX_SR_signal);
		  signal_.fitTo(signalHistogram, RooFit::Range(SR_lo,SR_hi), RooFit::Save(),RooFit::SumW2Error(kFALSE)); 
		  RooPlot *plot=x->frame();
	          //RooRealVar *signal_p0, *signal_p1, *signal_p2, *signal_p3, *signal_p7, *signal_p8;	
		  // signal_p0 = new RooRealVar("signal_p0", "signal_p0", sg_p0->getVal());
                  //signal_p1 = new RooRealVar("signal_p1", "signal_p1", sg_p1->getVal());
                  //signal_p2 = new RooRealVar("signal_p2", "signal_p2", sg_p2->getVal());
                  //signal_p3 = new RooRealVar("signal_p3", "signal_p3", sg_p3->getVal());
                  //signal_p8 = new RooRealVar("signal_p8", "signal_p8", sg_p8->getVal());
               /*                   
                  RooGaussian g1("g1", "Signal Prediction", *x, *sg_p0, *sg_p1);
                  RooGaussian g2("g2", "Combinatoric", *x, *sg_p2,*sg_p3);
                                  
                  RooAddPdf signal_bkg("signal_bkg", "signal_bkg", RooArgList(g1,g2), *sg_p8);
		  RooRealVar * signal_bkg_norm;
		  if(point[h]=="270") signal_bkg_norm= new RooRealVar("signal_bkg_norm","signal_bkg_norm", 0., -0.27, 0.27,"");	
		  if(point[h]=="300") signal_bkg_norm= new RooRealVar("signal_bkg_norm","signal_bkg_norm", 0., -0.55, 0.55,""); 
 		  if(point[h]=="350") signal_bkg_norm= new RooRealVar("signal_bkg_norm","signal_bkg_norm", 0., -0.24, 0.24,"");
	          if(point[h]=="400") signal_bkg_norm= new RooRealVar("signal_bkg_norm","signal_bkg_norm", 0., -0.09, 0.09,"");
*/
	
		   //sg_p8.setConstant(kTRUE);
		 // sg_p2.setConstant(kTRUE);
		 // sg_p3.setConstant(kTRUE);
		 // sg_p1.setConstant(kTRUE);

         	 /* int toyy = toy;

//.........这里部分代码省略.........
  //allVals[5] = 1.;
  //allVals[6] = -0.5;
  allVals[5] = 1.;
  allVals[6] = -1;

  double allValsSorted[7];

  for(int i=0; i<7; ++i) {//make allValsSorted increase with the order
    allValsSorted[i]=allVals[i];
    int tSpot = i;
    
    for(int j = i-1; j>=0; --j) {
      if(allValsSorted[j] > allValsSorted[tSpot]) {
	// switch
	double tempVal = allValsSorted[j];	  
	allValsSorted[j] = allValsSorted[tSpot];
	allValsSorted[tSpot] = tempVal;
	tSpot--;
      } else break;
    }
  }  
  
  for(int i=6-numCat; i<=6; ++i)
    catVals[i-6+numCat] = allValsSorted[i];//catVals from -1 to 1 increase with the order
  
  //==========3.read in the input tree==============  
  TFile* modFile = new TFile(modelName.Data());
  TNtuple* tup = (TNtuple*) modFile->FindObjectAny("MITMVAtuple");

  //============fit=============
  //-----------set the independent variable----------------
  RooRealVar* mass = new RooRealVar("mass","",100.,180.);
  mass->setRange(100.,180.);
  mass->setBins(10000,"cache");//bins for numrical integration
  mass->setBins(320);//bins for fitting and plotting
  //-----------set the variable to make dataset-------------
  RooRealVar* weight = new RooRealVar("weight","",0.,100.);
  RooRealVar* proc   = new RooRealVar("proc"  ,"",0.,10.);
  RooRealVar* mva    = new RooRealVar("mva"   ,"",-10.,10.);
  RooRealVar* vbftag = new RooRealVar("vbftag"   ,"",0.,1.); 
  //----------prepare the dataset---------------------------
  TString* dataNames      = new TString[numCat+1];
  TString* sigNames       = new TString[numCat+1];
  RooDataHist** sigCat = new RooDataHist*[numCat+1];
  RooDataSet** dataCat = new RooDataSet*[numCat+1];
  RooDataHist** dataHCat = new RooDataHist*[numCat+1];
  TString baseCut = "mass > 100 && mass < 180 ";
  TString baseBG  = baseCut + TString(" && proc > 3 && proc < 11");//MC background
  TString baseSIG = baseCut + TString(" && proc < 4 ");//MC signal
  TString baseObsData = baseCut + TString(" && proc==11");//obs data for count
  int*    numObsData= new int[numCat+1];
  //----------prepare the pdf-------------------------------
  //signal model
  RooAbsPdf**   sigpdfcat_part1 = new RooAbsPdf*[numCat+1];
  RooAbsPdf**   sigpdfcat_part2 = new RooAbsPdf*[numCat+1];
  RooAbsPdf**   sigpdfcat       = new RooAbsPdf*[numCat+1];
  RooRealVar**  sigMean_part1   = new RooRealVar*[numCat+1];
  RooRealVar**  sigWidth_part1  = new RooRealVar*[numCat+1];
  RooRealVar**  sigMean_part2   = new RooRealVar*[numCat+1];
  RooRealVar**  sigWidth_part2  = new RooRealVar*[numCat+1];
  RooRealVar**  sigRatio        = new RooRealVar*[numCat+1];
  double* nominalSignal   = new double[numCat+1];
  TH1D** sigH = new TH1D*[numCat+1];//ming:used to get nominalSignal and input to RooDataHist
  TH1D** ObsDataH = new TH1D*[numCat+1];//obs data for count
  //background model
  //exp

int main(){
  
  system("mkdir -p plots");
  RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
  TFile *bkgFile = TFile::Open("comb_svn/hgg.inputbkgdata_8TeV_MVA.root");
  TFile *sigFile = TFile::Open("comb_svn/hgg.inputsig_8TeV_nosplitVH_MVA.root");
  RooWorkspace *bkgWS = (RooWorkspace*)bkgFile->Get("cms_hgg_workspace");
  RooWorkspace *sigWS = (RooWorkspace*)sigFile->Get("wsig_8TeV");

  RooRealVar *mass = (RooRealVar*)bkgWS->var("CMS_hgg_mass");
  RooRealVar *mu = new RooRealVar("mu","mu",-5.,5.);

  mass->setBins(320);
  cout << mass->getBins() << endl;
  RooDataSet *dataAll;
  
  int firstCat=1;
  int lastCat=1;
  float mu_low=-1.;
  float mu_high=3.;
  float mu_step=0.01;

  vector<pair<double,TGraph*> > minNlltrack;

  for (int cat=firstCat; cat<=lastCat; cat++){
    RooDataSet *data = (RooDataSet*)bkgWS->data(Form("data_mass_cat%d",cat));
    if (cat==firstCat) dataAll = (RooDataSet*)data->Clone("data_mass_all");
    else dataAll->append(*data);
    RooDataHist *dataBinned = new RooDataHist(Form("roohist_data_mass_cat%d",cat),Form("roohist_data_mass_cat%d",cat),RooArgSet(*mass),*data);
    RooDataSet *sigMC = (RooDataSet*)sigWS->data(Form("sig_mass_m125_cat%d",cat));

    if (!dataBinned || !sigMC){
      cerr << "ERROR -- one of data or signal is NULL" << endl;
      exit(1);
    }
    
    // Construct PDFs for this category using PdfModelBuilder
    PdfModelBuilder modelBuilder;
    modelBuilder.setObsVar(mass);
    modelBuilder.setSignalModifier(mu);
    // For Standard Analysis
    //if (cat>=0 && cat<=3) modelBuilder.addBkgPdf("Bernstein",5,Form("pol5_cat%d",cat));
    //if (cat>=4 && cat<=5) modelBuilder.addBkgPdf("Bernstein",4,Form("pol4_cat%d",cat));
    //if (cat>=6 && cat<=8) modelBuilder.addBkgPdf("Bernstein",3,Form("pol3_cat%d",cat));
    // To Profile Multiple PDFs
    if (cat==0 || cat==1 || cat==2 || cat==3){
      modelBuilder.addBkgPdf("Bernstein",4,Form("pol4_cat%d",cat));
      modelBuilder.addBkgPdf("Bernstein",5,Form("pol5_cat%d",cat));
      modelBuilder.addBkgPdf("Bernstein",6,Form("pol6_cat%d",cat));
      /*
      modelBuilder.addBkgPdf("PowerLaw",1,Form("pow1_cat%d",cat));
      modelBuilder.addBkgPdf("PowerLaw",3,Form("pow3_cat%d",cat));
      modelBuilder.addBkgPdf("PowerLaw",5,Form("pow5_cat%d",cat));
      modelBuilder.addBkgPdf("Exponential",1,Form("exp1_cat%d",cat));
      modelBuilder.addBkgPdf("Exponential",3,Form("exp3_cat%d",cat));
      modelBuilder.addBkgPdf("Exponential",5,Form("exp5_cat%d",cat));
      modelBuilder.addBkgPdf("Laurent",1,Form("lau1_cat%d",cat));
      modelBuilder.addBkgPdf("Laurent",3,Form("lau3_cat%d",cat));
      modelBuilder.addBkgPdf("Laurent",5,Form("lau5_cat%d",cat));
      */
    }
    if (cat==4 || cat==5 || cat==6 || cat==7 || cat==8) {
      modelBuilder.addBkgPdf("Bernstein",3,Form("pol3_cat%d",cat));
      modelBuilder.addBkgPdf("Bernstein",4,Form("pol4_cat%d",cat));
      /*
      modelBuilder.addBkgPdf("PowerLaw",1,Form("pow1_cat%d",cat));
      modelBuilder.addBkgPdf("PowerLaw",3,Form("pow3_cat%d",cat));
      modelBuilder.addBkgPdf("Exponential",1,Form("exp1_cat%d",cat));
      modelBuilder.addBkgPdf("Exponential",3,Form("exp3_cat%d",cat));
      modelBuilder.addBkgPdf("Laurent",1,Form("lau1_cat%d",cat));
      modelBuilder.addBkgPdf("Laurent",3,Form("lau3_cat%d",cat));
      */
    }
    map<string,RooAbsPdf*> bkgPdfs = modelBuilder.getBkgPdfs();
    modelBuilder.setSignalPdfFromMC(sigMC);
    modelBuilder.makeSBPdfs();
    map<string,RooAbsPdf*> sbPdfs = modelBuilder.getSBPdfs();

    modelBuilder.fitToData(dataBinned,true,true);
    modelBuilder.fitToData(dataBinned,false,true);

    modelBuilder.throwToy(Form("cat%d_toy0",cat),dataBinned->sumEntries(),true,true);

    // Profile this category using ProfileMultiplePdfs
    ProfileMultiplePdfs profiler;
    for (map<string,RooAbsPdf*>::iterator pdf=sbPdfs.begin(); pdf!=sbPdfs.end(); pdf++) {
      string bkgOnlyName = pdf->first.substr(pdf->first.find("sb_")+3,string::npos);
      if (bkgPdfs.find(bkgOnlyName)==bkgPdfs.end()){
        cerr << "ERROR -- couldn't find bkg only pdf " << bkgOnlyName << " for SB pdf " << pdf->first << endl;
        pdf->second->fitTo(*dataBinned);
        exit(1);
      }
      int nParams = bkgPdfs[bkgOnlyName]->getVariables()->getSize()-1;
      profiler.addPdf(pdf->second,2*nParams);
      //profiler.addPdf(pdf->second);
      cout << pdf->second->GetName() << " nParams=" << pdf->second->getVariables()->getSize() << " nBkgParams=" << nParams << endl;
    }
    profiler.printPdfs();
    //cout << "Continue?" << endl;
    //string bus; cin >> bus;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  RooAbsReal *sigmeanlim = ws->function("sigmeanlim");
  RooAbsReal *sigwidthlim = ws->function("sigwidthlim");
  RooAbsReal *signlim = ws->function("signlim");
  RooAbsReal *sign2lim = ws->function("sign2lim");
  RooAbsReal *alphalim = ws->function("sigalphalim");
  RooAbsReal *alpha2lim = ws->function("sigalpha2lim");  

  //RooFormulaVar ecor("ecor","","1./(@0*@1)",RooArgList(*tgtvar,*sigmeanlim));
  RooFormulaVar ecor("ecor","","@1/@0",RooArgList(*tgtvar,*sigmeanlim));
  //RooFormulaVar ecor("ecor","","@0/@1",RooArgList(*tgtvar,*sigmeanlim));
  //RooFormulaVar ecor("ecor","","exp(@[email protected])",RooArgList(*tgtvar,*sigmeanlim));
  
  
  RooAbsReal *condnll = sigpdf->createNLL(*hdata,ConditionalObservables(sigmeant->FuncVars()));

  double condnllval = condnll->getVal();
  
  
  //RooFormulaVar ecor("ecor","","@1/@0",RooArgList(*tgtvar,*sigmeanlim));
  //RooFormulaVar ecor("ecor","","@0/@1",RooArgList(*tgtvar,*sigmeanlim));
  //RooFormulaVar ecor("ecor","","@0",RooArgList(*tgtvar));
  //RooRealVar *ecorvar = (RooRealVar*)hdata->addColumn(ecor);
//   ecorvar->setRange(0.,2.);
//   ecorvar->setBins(800);
  
//    RooFormulaVar raw("raw","","1./@0",RooArgList(*tgtvar));
//    //RooRealVar *rawvar = (RooRealVar*)hdata->addColumn(raw);
//    rawvar->setRange(0.,2.);
//    rawvar->setBins(800);

/*  RooFormulaVar eraw("eraw","","@0",RooArgList(*tgtvar));
  RooRealVar *erawvar = (RooRealVar*)hdatasig->addColumn(eraw);
  erawvar->setRange(0.,2.);
  erawvar->setBins(400); */ 

  //RooFormulaVar ecor("ptcor","","@0/(@1)",RooArgList(*tgtvar,*sigmeanlim));
  
  RooDataSet *hdataclone = new RooDataSet(*hdata,"hdataclone");
  RooRealVar *ecorvar = (RooRealVar*)hdataclone->addColumn(ecor);
  RooRealVar *meanvar = (RooRealVar*)hdataclone->addColumn(*sigmeanlim);
  RooRealVar *widthvar = (RooRealVar*)hdataclone->addColumn(*sigwidthlim);
  RooRealVar *nvar = 0;
  if (signlim) nvar = (RooRealVar*)hdataclone->addColumn(*signlim);
  RooRealVar *n2var = 0;
  if (sign2lim) n2var = (RooRealVar*)hdataclone->addColumn(*sign2lim);
  RooRealVar *alphavar = 0;; 
  if (alphalim) alphavar = (RooRealVar*)hdataclone->addColumn(*alphalim);
  RooRealVar *alpha2var = 0;
  if (alpha2lim) alpha2var = (RooRealVar*)hdataclone->addColumn(*alpha2lim);  
  
  RooFormulaVar ecorfull("ecorfull","","@0*@1",RooArgList(*sigmeanlim,*rawevar));
  RooRealVar *ecorfullvar = (RooRealVar*)hdataclone->addColumn(ecorfull);
  
  RooFormulaVar ediff("ediff","","(@0 - @1)/@1",RooArgList(*nomevar,ecorfull));
  RooRealVar *ediffvar = (RooRealVar*)hdataclone->addColumn(ediff);  
  
  RooFormulaVar fullerr("fullerr","","@0*@1",RooArgList(*ecorvar,*sigwidthlim));
  RooRealVar *fullerrvar = (RooRealVar*)hdataclone->addColumn(fullerr);
  
  RooFormulaVar relerr("relerr","","@0/@1",RooArgList(*sigwidthlim,*sigmeanlim));
  RooRealVar *relerrvar = (RooRealVar*)hdataclone->addColumn(relerr);  
  
  ecorvar->setRange(0.,2.);
  ecorvar->setBins(800);  
  
  RooFormulaVar raw("raw","","1./@0",RooArgList(*tgtvar));