C++ RooWorkspace类代码示例

int main(int argc, char* argv[]){

  string bkgFileName;
  string sigFileName;
  string sigWSName;
  string bkgWSName;
  string outFileName;
  string datFileName;
  string outDir;
  int cat;
  int ntoys;
  int jobn;
  int seed;
  float mu_low;
  float mu_high;
  float mu_step;
  float expectSignal;
  int expectSignalMass;
  bool skipPlots=false;
  int verbosity;
  bool throwHybridToys=false;
  vector<float> switchMass;
  vector<string> switchFunc;

  po::options_description desc("Allowed options");
  desc.add_options()
    ("help,h",                                                                                  "Show help")
    ("sigfilename,s", po::value<string>(&sigFileName),                                          "Signal file name")
    ("bkgfilename,b", po::value<string>(&bkgFileName),                                          "Background file name")
    ("sigwsname", po::value<string>(&sigWSName)->default_value("cms_hgg_workspace"),            "Signal workspace name")
    ("bkgwsname", po::value<string>(&bkgWSName)->default_value("cms_hgg_workspace"),            "Background workspace name")
    ("outfilename,o", po::value<string>(&outFileName)->default_value("BiasStudyOut.root"),      "Output file name")
    ("datfile,d", po::value<string>(&datFileName)->default_value("config.dat"),                 "Name of datfile containing pdf info")
    ("outDir,D", po::value<string>(&outDir)->default_value("./"),                               "Name of out directory for plots")
    ("cat,c", po::value<int>(&cat),                                                             "Category")
    ("ntoys,t", po::value<int>(&ntoys)->default_value(0),                                       "Number of toys to run")
    ("jobn,j", po::value<int>(&jobn)->default_value(0),                                         "Job number")
    ("seed,r", po::value<int>(&seed)->default_value(0),                                         "Set random seed")
    ("mulow,L", po::value<float>(&mu_low)->default_value(-3.),                                  "Value of mu to start scan")
    ("muhigh,H", po::value<float>(&mu_high)->default_value(3.),                                 "Value of mu to end scan")
    ("mustep,S", po::value<float>(&mu_step)->default_value(0.01),                               "Value of mu step size")
    ("expectSignal", po::value<float>(&expectSignal)->default_value(0.),                        "Inject signal into toy")
    ("expectSignalMass", po::value<int>(&expectSignalMass)->default_value(125),                 "Inject signal at this mass")
    ("skipPlots",                                                                               "Skip full profile and toy plots")                        
    ("verbosity,v", po::value<int>(&verbosity)->default_value(0),                               "Verbosity level")
  ;    
  
  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("skipPlots")) skipPlots=true;
  if (expectSignalMass!=110 && expectSignalMass!=115 && expectSignalMass!=120 && expectSignalMass!=125 && expectSignalMass!=130 && expectSignalMass!=135 && expectSignalMass!=140 && expectSignalMass!=145 && expectSignalMass!=150){
    cerr << "ERROR - expectSignalMass has to be integer in range (110,150,5)" << endl;
    exit(1);
  }

  vector<pair<int,pair<string,string> > > toysMap;
  vector<pair<int,pair<string,string> > > fabianMap;
  vector<pair<int,pair<string,string> > > paulMap;
  readDatFile(datFileName,cat,toysMap,fabianMap,paulMap);
  
  cout << "Toy vector.." << endl;
  printOptionsMap(toysMap);
  cout << "Fabian vector.." << endl;
  printOptionsMap(fabianMap);
  cout << "Paul vector.." << endl;
  printOptionsMap(paulMap);
  
  TStopwatch sw;
  sw.Start();
 
  if (verbosity<1) {
    RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
    RooMsgService::instance().setSilentMode(true);
  }
  
  TFile *bkgFile = TFile::Open(bkgFileName.c_str());
  TFile *sigFile = TFile::Open(sigFileName.c_str());

  //RooWorkspace *bkgWS = (RooWorkspace*)bkgFile->Get("cms_hgg_workspace");
  RooWorkspace *bkgWS = (RooWorkspace*)bkgFile->Get(bkgWSName.c_str());
  RooWorkspace *sigWS = (RooWorkspace*)sigFile->Get(sigWSName.c_str());

  if (!bkgWS || !sigWS){
    cerr << "ERROR - one of signal or background workspace is NULL" << endl;
    exit(1);
  }

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

  TFile *outFile = new TFile(outFileName.c_str(),"RECREATE");
  TTree *muTree = new TTree("muTree","muTree");
  int toyn;
  vector<string> truthModel;
  vector<double> muFab;
  vector<double> muPaul;
  vector<double> muChi2;
  vector<double> muAIC;
//.........这里部分代码省略.........

void StandardBayesianNumericalDemo(const char* infile = "",
		      const char* workspaceName = "combined",
		      const char* modelConfigName = "ModelConfig",
		      const char* dataName = "obsData"){

  /////////////////////////////////////////////////////////////
  // First part is just to access a user-defined file 
  // or create the standard example file if it doesn't exist
  ////////////////////////////////////////////////////////////
  TString filename = infile;
  if (filename.IsNull()) { 
    filename = "results/example_combined_GaussExample_model.root";
    std::cout << "Use standard file generated with HistFactory : " << filename << std::endl;
  }

  // Check if example input file exists
  TFile *file = TFile::Open(filename);

  // if input file was specified but not found, quit
  if(!file && !TString(infile).IsNull()){
     cout <<"file " << filename << " not found" << endl;
     return;
  } 

  // if default file not found, try to create it
  if(!file ){
    // Normally this would be run on the command line
    cout <<"will run standard hist2workspace example"<<endl;
    gROOT->ProcessLine(".! prepareHistFactory .");
    gROOT->ProcessLine(".! hist2workspace config/example.xml");
    cout <<"\n\n---------------------"<<endl;
    cout <<"Done creating example input"<<endl;
    cout <<"---------------------\n\n"<<endl;

    // now try to access the file again
    file = TFile::Open(filename);
  }

  if(!file){
    // if it is still not there, then we can't continue
    cout << "Not able to run hist2workspace to create example input" <<endl;
    return;
  }

  
  /////////////////////////////////////////////////////////////
  // Tutorial starts here
  ////////////////////////////////////////////////////////////

  // get the workspace out of the file
  RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);
  if(!w){
    cout <<"workspace not found" << endl;
    return;
  }

  // get the modelConfig out of the file
  ModelConfig* mc = (ModelConfig*) w->obj(modelConfigName);

  // get the modelConfig out of the file
  RooAbsData* data = w->data(dataName);

  // make sure ingredients are found
  if(!data || !mc){
    w->Print();
    cout << "data or ModelConfig was not found" <<endl;
    return;
  }

  /////////////////////////////////////////////
  // create and use the BayesianCalculator
  // to find and plot the 95% credible interval
  // on the parameter of interest as specified
  // in the model config
  
  // before we do that, we must specify our prior
  // it belongs in the model config, but it may not have
  // been specified
  RooUniform prior("prior","",*mc->GetParametersOfInterest());
  w->import(prior);
  mc->SetPriorPdf(*w->pdf("prior"));

  // do without systematics
  //mc->SetNuisanceParameters(RooArgSet() );

  
  BayesianCalculator bayesianCalc(*data,*mc);
  bayesianCalc.SetConfidenceLevel(0.95); // 95% interval

  // default of the calculator is central interval.  here use shortest , central or upper limit depending on input
  // doing a shortest interval might require a longer time since it requires a scan of the posterior function
  if (intervalType == 0)  bayesianCalc.SetShortestInterval(); // for shortest interval
  if (intervalType == 1)  bayesianCalc.SetLeftSideTailFraction(0.5); // for central interval
  if (intervalType == 2)  bayesianCalc.SetLeftSideTailFraction(0.); // for upper limit

  if (!integrationType.IsNull() ) { 
     bayesianCalc.SetIntegrationType(integrationType); // set integrationType
     bayesianCalc.SetNumIters(nToys); // set number of ietrations (i.e. number of toys for MC integrations)
  }

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

void RA2bHypoTestInvDemo(const char * fileName =0,
			 const char * wsName = "combined",
			 const char * modelSBName = "ModelConfig",
			 const char * modelBName = "",
			 const char * dataName = "obsData",                 
			 int calculatorType = 0,
			 int testStatType = 3, 
			 bool useCls = true ,  
			 int npoints = 5,   
			 double poimin = 0,  
			 double poimax = 5, 
			 int ntoys=1000,
			 int mgl = -1,
			 int mlsp = -1,
			 const char * outFileName = "test")    
{
/*

   Other Parameter to pass in tutorial
   apart from standard for filename, ws, modelconfig and data

    type = 0 Freq calculator 
    type = 1 Hybrid 

    testStatType = 0 LEP
                 = 1 Tevatron 
                 = 2 Profile Likelihood
                 = 3 Profile Likelihood one sided (i.e. = 0 if mu < mu_hat)

    useCLs          scan for CLs (otherwise for CLs+b)    

    npoints:        number of points to scan , for autoscan set npoints = -1 

    poimin,poimax:  min/max value to scan in case of fixed scans 
                    (if min >= max, try to find automatically)                           

    ntoys:         number of toys to use 

    extra options are available as global paramters of the macro. They are: 

    plotHypoTestResult   plot result of tests at each point (TS distributions) 
    useProof = true;
    writeResult = true;
    nworkers = 4;


   */

   if (fileName==0) { 
      fileName = "results/example_combined_GaussExample_model.root";
      std::cout << "Use standard file generated with HistFactory :" << fileName << std::endl;
   }
   TFile * file = new TFile(fileName); 

   RooWorkspace * w = dynamic_cast<RooWorkspace*>( file->Get(wsName) );
   HypoTestInverterResult * r = 0; 
   std::cout << w << "\t" << fileName << std::endl;
   if (w != NULL) {
      r = RunInverter(w, modelSBName, modelBName, dataName, calculatorType, testStatType, npoints, poimin, poimax,  ntoys, useCls );    
      if (!r) { 
         std::cerr << "Error running the HypoTestInverter - Exit " << std::endl;
         return;          
      }
   }
   else 
   { 
      // case workspace is not present look for the inverter result
      std::cout << "Reading an HypoTestInverterResult with name " << wsName << " from file " << fileName << std::endl;
      r = dynamic_cast<HypoTestInverterResult*>( file->Get(wsName) ); //
      if (!r) { 
         std::cerr << "File " << fileName << " does not contain a workspace or an HypoTestInverterResult - Exit " 
                   << std::endl;
         file->ls();
         return; 
      }
   }		
      		


   printf("\n\n") ;
   HypoTestResult* htr = r->GetResult(0) ;
   printf("  Data value for test stat : %7.3f\n", htr->GetTestStatisticData() ) ;
   printf("  CLsplusb : %9.4f\n", r->CLsplusb(0) ) ;
   printf("  CLb      : %9.4f\n", r->CLb(0) ) ;
   printf("  CLs      : %9.4f\n", r->CLs(0) ) ;
   printf("\n\n") ;
   cout << flush ;

   double upperLimit = r->UpperLimit();
   double ulError = r->UpperLimitEstimatedError();


   std::cout << "The computed upper limit is: " << upperLimit << " +/- " << ulError << std::endl;
 
   const int nEntries = r->ArraySize();


   const char *  typeName = (calculatorType == 0) ? "Frequentist" : "Hybrid";
   const char * resultName = (w) ? w->GetName() : r->GetName();
   TString plotTitle = TString::Format("%s CL Scan for workspace %s",typeName,resultName);
//.........这里部分代码省略.........

double StandardFrequentistDiscovery(
   const char* infile = "",
   const char* workspaceName = "channel1",
   const char* modelConfigNameSB = "ModelConfig",
   const char* dataName = "obsData",
   int toys = 1000,
   double poiValueForBackground = 0.0,
   double poiValueForSignal = 1.0
) {

   // The workspace contains the model for s+b. The b model is "autogenerated"
   // by copying s+b and setting the one parameter of interest to zero.
   // To keep the script simple, multiple parameters of interest or different
   // functional forms of the b model are not supported.

   // for now, assume there is only one parameter of interest, and these are
   // its values:

   /////////////////////////////////////////////////////////////
   // First part is just to access a user-defined file
   // or create the standard example file if it doesn't exist
   ////////////////////////////////////////////////////////////
   const char* filename = "";
   if (!strcmp(infile,"")) {
      filename = "results/example_channel1_GammaExample_model.root";
      bool fileExist = !gSystem->AccessPathName(filename); // note opposite return code
      // if file does not exists generate with histfactory
      if (!fileExist) {
#ifdef _WIN32
         cout << "HistFactory file cannot be generated on Windows - exit" << endl;
         return -1;
#endif
         // Normally this would be run on the command line
         cout <<"will run standard hist2workspace example"<<endl;
         gROOT->ProcessLine(".! prepareHistFactory .");
         gROOT->ProcessLine(".! hist2workspace config/example.xml");
         cout <<"\n\n---------------------"<<endl;
         cout <<"Done creating example input"<<endl;
         cout <<"---------------------\n\n"<<endl;
      }
      
   }
   else
      filename = infile;
   
   // Try to open the file
   TFile *file = TFile::Open(filename);
   
   // if input file was specified byt not found, quit
   if(!file ){
      cout <<"StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
      return -1;
   } 


   /////////////////////////////////////////////////////////////
   // Tutorial starts here
   ////////////////////////////////////////////////////////////

   TStopwatch *mn_t = new TStopwatch;
   mn_t->Start();

   // get the workspace out of the file
   RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);
   if (!w) {
      cout << "workspace not found" << endl;
      return -1.0;
   }

   // get the modelConfig out of the file
   ModelConfig* mc = (ModelConfig*) w->obj(modelConfigNameSB);

   // get the data out of the file
   RooAbsData* data = w->data(dataName);

   // make sure ingredients are found
   if (!data || !mc) {
      w->Print();
      cout << "data or ModelConfig was not found" << endl;
      return -1.0;
   }


   RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
   firstPOI->setVal(poiValueForSignal);
   mc->SetSnapshot(*mc->GetParametersOfInterest());
   // create null model
   ModelConfig *mcNull = mc->Clone("ModelConfigNull");
   firstPOI->setVal(poiValueForBackground);
   mcNull->SetSnapshot(*(RooArgSet*)mcNull->GetParametersOfInterest()->snapshot());



   // ----------------------------------------------------
   // Configure a ProfileLikelihoodTestStat and a SimpleLikelihoodRatioTestStat
   // to use simultaneously with ToyMCSampler
   ProfileLikelihoodTestStat* plts =  new ProfileLikelihoodTestStat(*mc->GetPdf());
   plts->SetOneSidedDiscovery(true);
   plts->SetVarName( "q_{0}/2" );
   
//.........这里部分代码省略.........

MakeAICFits::MakeAICFits() {
  //RooRandom::randomGenerator()->SetSeed(314159);
  TFile *f = TFile::Open("/mnt/hadoop/store/user/amott/Hgg2013/Hgg/workspaces/hgg_22Jan2013_R9_CIC.root");
  RooWorkspace *w = (RooWorkspace*) f->Get("cms_hgg_spin_workspace");
  RooAbsData *d = w->data("Data_Combined");
  RooCategory* cat = (RooCategory*)cms_hgg_spin_workspace->obj("evtcat");
  RooAbsData *dc = d->reduce("evtcat==evtcat::cat4");
  //RooAbsData *dc = w->data("Data_Combined")->reduce("evtcar==evtcat::cat0");
  RooRealVar *mass = w->var("mass");

  const Int_t nToys = 1;

  std::cout<<"==========  Data Set Made    ==========="<<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};
  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;
    }
    int SampleSize = dc->sumEntries();
    RooPlot* frame = mass->frame();
    TLegend *leg = new TLegend(0.55,0.55,0.9,0.9, NULL, "NDC");
    dc->plotOn(frame);
    leg->AddEntry(dc,"Hgg Background Cat 4", "lep");

    for (int type=0; type<7; type++) {
      std::cout<<type<<endl;
    }
    int display = 7;
    for (int type=0; type<display; type++) {
      RooAbsPdf* ModelShape;
      if (type<7) {
	//std::cout<<"Model Shape:    "<<type<<std::endl;
	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(0,mass);
	RooAbsPdf* Model2 = MakeAICFits::getBackgroundPdf(1,mass);
	RooAbsPdf* Model3 = MakeAICFits::getBackgroundPdf(2,mass);
	RooAbsPdf* Model4 = MakeAICFits::getBackgroundPdf(3,mass);
	RooAbsPdf* Model5 = MakeAICFits::getBackgroundPdf(4,mass);
	int k = MakeAICFits::Num_Params(3);
	k+= MakeAICFits::Num_Params(1);
	k+= MakeAICFits::Num_Params(0);
	//k+= MakeAICFits::Num_Params(3);
	//k+= MakeAICFits::Num_Params(4);
	RooRealVar* modratio1 = new RooRealVar("modrat1", "modrat1", 0.62, 0.6, 0.7);
	RooRealVar* modratio2 = new RooRealVar("modrat2", "modrat2", 0.29, 0.25, 0.35);
	RooRealVar* modratio3 = new RooRealVar("modrat3", "modrat3", 0.01);
	//RooRealVar* modratio4 = new RooRealVar("modrat4", "modrat4", 0.25);
	ModelShape = new RooAddPdf("Composite", "Background Model", RooArgList(*Model1, *Model4, *Model2), RooArgList(*modratio1, *modratio2));
      }
      RooRealVar *Nbkg = new RooRealVar("Nbkg","N Background Events", SampleSize,0,1e9);
      RooExtendPdf *BkgModel = new RooExtendPdf("BKGFIT_bkgModel", "Background Model", *ModelShape, *Nbkg);
      TH1F* Model = new TH1F("Model", "Model", 100,0,100);
 
      RooFitResult *bkg_databkg = BkgModel->fitTo(*dc, RooFit::Save(kTRUE), RooFit::Optimize(0));
      if (type == 0) {
	BkgModel->plotOn(frame, RooFit::LineColor(kBlue));
	Model->SetLineColor(kBlue);
	leg->AddEntry(Model, "Exponential Model", "l");
      }
      if (type == 4) { 
	BkgModel->plotOn(frame, RooFit::LineColor(kRed));
	Model->SetLineColor(kRed);
	leg->AddEntry(Model, "Polynomial Model", "l");
      }
      if (type == 5) { 
	BkgModel->plotOn(frame, RooFit::LineColor(kGreen));
	Model->SetLineColor(kGreen);
	leg->AddEntry(Model, "Power Model", "l");
      }
      if (type == 7) {
	BkgModel->plotOn(frame, RooFit::LineColor(kMagenta));
	Model->SetLineColor(kMagenta);
	leg->AddEntry(Model, "Composite Model", "l");
      }
      Double_t bkg_databkg_Nll = bkg_databkg->minNll();
      Int_t covariance = bkg_databkg->covQual();
      avgcov[type] += covariance;
      //assert (covariance == 3);
      // Calculate AIC for each model
      LogLikelihood[type] += -bkg_databkg_Nll;
      AICc_bkg_array[type] += 2.*(k + k*(k + 1.)/(SampleSize - (k + 1.)) + bkg_databkg_Nll);
      AIC_bkg_array[type] += 2.*(k + bkg_databkg_Nll);
      // Clean up objects
      delete bkg_databkg;
      bkg_databkg_Nll = 0.;
    }
    //delete databkg;
    TCanvas *c = new TCanvas("", "", 800, 600);
    frame->Draw();
    leg->Draw();
//.........这里部分代码省略.........

void OneSidedFrequentistUpperLimitWithBands_intermediate(const char* infile = "",
					    const char* workspaceName = "combined",
					    const char* modelConfigName = "ModelConfig",
					    const char* dataName = "obsData"){


  double confidenceLevel=0.95;
  // degrade/improve number of pseudo-experiments used to define the confidence belt.  
  // value of 1 corresponds to default number of toys in the tail, which is 50/(1-confidenceLevel)
  double additionalToysFac = 1.;  
  int nPointsToScan = 30; // number of steps in the parameter of interest 
  int nToyMC = 100; // number of toys used to define the expected limit and band

  TStopwatch t;
  t.Start();
  /////////////////////////////////////////////////////////////
  // First part is just to access a user-defined file 
  // or create the standard example file if it doesn't exist
  ////////////////////////////////////////////////////////////
  const char* filename = "";
  if (!strcmp(infile,""))
    filename = "results/example_combined_GaussExample_model.root";
  else
    filename = infile;
  // Check if example input file exists
  TFile *file = TFile::Open(filename);

  // if input file was specified byt not found, quit
  if(!file && strcmp(infile,"")){
    cout <<"file not found" << endl;
    return;
  } 

  // if default file not found, try to create it
  if(!file ){
    // Normally this would be run on the command line
    cout <<"will run standard hist2workspace example"<<endl;
    gROOT->ProcessLine(".! prepareHistFactory .");
    gROOT->ProcessLine(".! hist2workspace config/example.xml");
    cout <<"\n\n---------------------"<<endl;
    cout <<"Done creating example input"<<endl;
    cout <<"---------------------\n\n"<<endl;
  }

  // now try to access the file again
  file = TFile::Open(filename);
  if(!file){
    // if it is still not there, then we can't continue
    cout << "Not able to run hist2workspace to create example input" <<endl;
    return;
  }

  
  /////////////////////////////////////////////////////////////
  // Now get the data and workspace
  ////////////////////////////////////////////////////////////

  // get the workspace out of the file
  RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);
  if(!w){
    cout <<"workspace not found" << endl;
    return;
  }

  // get the modelConfig out of the file
  ModelConfig* mc = (ModelConfig*) w->obj(modelConfigName);

  // get the modelConfig out of the file
  RooAbsData* data = w->data(dataName);

  // make sure ingredients are found
  if(!data || !mc){
    w->Print();
    cout << "data or ModelConfig was not found" <<endl;
    return;
  }

  cout << "Found data and ModelConfig:" <<endl;
  mc->Print();

  /////////////////////////////////////////////////////////////
  // Now get the POI for convenience
  // you may want to adjust the range of your POI
  ////////////////////////////////////////////////////////////
  RooRealVar* firstPOI = (RooRealVar*) mc->GetParametersOfInterest()->first();
  //  firstPOI->setMin(0);
  //  firstPOI->setMax(10);

  /////////////////////////////////////////////
  // create and use the FeldmanCousins tool
  // to find and plot the 95% confidence interval
  // on the parameter of interest as specified
  // in the model config
  // REMEMBER, we will change the test statistic
  // so this is NOT a Feldman-Cousins interval
  FeldmanCousins fc(*data,*mc);
  fc.SetConfidenceLevel(confidenceLevel); 
  fc.AdditionalNToysFactor(additionalToysFac); // improve sampling that defines confidence belt
  //  fc.UseAdaptiveSampling(true); // speed it up a bit, but don't use for expectd limits
  fc.SetNBins(nPointsToScan); // set how many points per parameter of interest to scan
//.........这里部分代码省略.........

void new_RA4(){
  
  // let's time this challenging example
  TStopwatch t;
  t.Start();

  // set RooFit random seed for reproducible results
  RooRandom::randomGenerator()->SetSeed(4357);

  // make model
  RooWorkspace* wspace = new RooWorkspace("wspace");

  wspace->factory("Gaussian::sigCons(prime_SigEff[0,-5,5], nom_SigEff[0,-5,5], 1)");
  wspace->factory("expr::SigEff('1.0*pow(1.20,@0)',prime_SigEff)"); // // 1+-20%, 1.20=exp(20%)

  wspace->factory("Poisson::on(non[0,50], sum::splusb(prod::SigUnc(s[0,0,50],SigEff),mainb[8.8,0,50],dilep[0.9,0,20],tau[2.3,0,20],QCD[0.,0,10],MC[0.1,0,4]))");

  wspace->factory("Gaussian::mcCons(prime_rho[0,-5,5], nom_rho[0,-5,5], 1)");
  wspace->factory("expr::rho('1.0*pow(1.39,@0)',prime_rho)"); // // 1+-39%
  wspace->factory("Poisson::off(noff[0,200], prod::rhob(mainb,rho,mu_plus_e[0.74,0.01,10],1.08))");
  wspace->factory("Gaussian::mcCons2(mu_plus_enom[0.74,0.01,4], mu_plus_e, sigmatwo[.05])");

  wspace->factory("Gaussian::dilep_pred(dilep_nom[0.9,0,20], dilep, sigma3[2.2])");
  wspace->factory("Gaussian::tau_pred(tau_nom[2.3,0,20], tau, sigma4[0.5])");
  wspace->factory("Gaussian::QCD_pred(QCD_nom[0.0,0,10], QCD, sigma5[1.0])");
  wspace->factory("Gaussian::MC_pred(MC_nom[0.1,0.01,4], MC, sigma7[0.14])");

  wspace->factory("PROD::model(on,off,mcCons,mcCons2,sigCons,dilep_pred,tau_pred,QCD_pred,MC_pred)");

  RooArgSet obs(*wspace->var("non"), *wspace->var("noff"), *wspace->var("mu_plus_enom"), *wspace->var("dilep_nom"), *wspace->var("tau_nom"), "obs");
  obs.add(*wspace->var("QCD_nom"));  obs.add(*wspace->var("MC_nom"));
  RooArgSet globalObs(*wspace->var("nom_SigEff"), *wspace->var("nom_rho"), "global_obs");
  // fix global observables to their nominal values
  wspace->var("nom_SigEff")->setConstant();
  wspace->var("nom_rho")->setConstant();

  RooArgSet poi(*wspace->var("s"), "poi");
  RooArgSet nuis(*wspace->var("mainb"), *wspace->var("prime_rho"), *wspace->var("prime_SigEff"), *wspace->var("mu_plus_e"), *wspace->var("dilep"), *wspace->var("tau"), "nuis");
  nuis.add(*wspace->var("QCD"));  nuis.add(*wspace->var("MC"));


  wspace->factory("Uniform::prior_poi({s})");
  wspace->factory("Uniform::prior_nuis({mainb,mu_plus_e,dilep,tau,QCD,MC})");
  wspace->factory("PROD::prior(prior_poi,prior_nuis)");

  wspace->var("non")->setVal(8); //observed
  //wspace->var("non")->setVal(12); //expected observation
  wspace->var("noff")->setVal(7); //observed events in control region
  wspace->var("mu_plus_enom")->setVal(0.74);
  wspace->var("dilep_nom")->setVal(0.9);
  wspace->var("tau_nom")->setVal(2.3);
  wspace->var("QCD")->setVal(0.0);
  wspace->var("MC")->setVal(0.1);


  RooDataSet * data = new RooDataSet("data","",obs);
  data->add(obs);
  wspace->import(*data);


  /////////////////////////////////////////////////////
  // Now the statistical tests
  // model config
  ModelConfig* pSbModel = new ModelConfig("SbModel");
  pSbModel->SetWorkspace(*wspace);
  pSbModel->SetPdf(*wspace->pdf("model"));
  pSbModel->SetPriorPdf(*wspace->pdf("prior"));
  pSbModel->SetParametersOfInterest(poi);
  pSbModel->SetNuisanceParameters(nuis);
  pSbModel->SetObservables(obs);
  pSbModel->SetGlobalObservables(globalObs);
  wspace->import(*pSbModel);

  // set all but obs, poi and nuisance to const
  SetConstants(wspace, pSbModel);
  wspace->import(*pSbModel);


  Double_t poiValueForBModel = 0.0;
  ModelConfig* pBModel = new ModelConfig(*(RooStats::ModelConfig *)wspace->obj("SbModel"));
  pBModel->SetName("BModel");
  pBModel->SetWorkspace(*wspace);
  wspace->import(*pBModel);


  RooAbsReal * pNll = pSbModel->GetPdf()->createNLL(*data);
  RooAbsReal * pProfile = pNll->createProfile(RooArgSet());
  pProfile->getVal(); // this will do fit and set POI and nuisance parameters to fitted values
  RooArgSet * pPoiAndNuisance = new RooArgSet();
  //if(pSbModel->GetNuisanceParameters())
  //  pPoiAndNuisance->add(*pSbModel->GetNuisanceParameters());
  pPoiAndNuisance->add(*pSbModel->GetParametersOfInterest());
  cout << "\nWill save these parameter points that correspond to the fit to data" << endl;
  pPoiAndNuisance->Print("v");
  pSbModel->SetSnapshot(*pPoiAndNuisance);
  delete pProfile;
  delete pNll;
  delete pPoiAndNuisance;


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

//.........这里部分代码省略.........
  TFile *file = TFile::Open(filename);

  // if input file was specified byt not found, quit
  if(!file && strcmp(infile,"")){
    cout <<"file not found" << endl;
    return;
  } 

  // if default file not found, try to create it
  if(!file ){
    // Normally this would be run on the command line
    cout <<"will run standard hist2workspace example"<<endl;
    gROOT->ProcessLine(".! prepareHistFactory .");
    gROOT->ProcessLine(".! hist2workspace config/example.xml");
    cout <<"\n\n---------------------"<<endl;
    cout <<"Done creating example input"<<endl;
    cout <<"---------------------\n\n"<<endl;
  }

  // now try to access the file again
  file = TFile::Open(filename);
  if(!file){
    // if it is still not there, then we can't continue
    cout << "Not able to run hist2workspace to create example input" <<endl;
    return;
  }

  
  /////////////////////////////////////////////////////////////
  // Tutorial starts here
  ////////////////////////////////////////////////////////////

  // get the workspace out of the file
  RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);
  if(!w){
    cout <<"workspace not found" << endl;
    return;
  }
  w->Print();

  // get the modelConfig out of the file
  ModelConfig* sbModel = (ModelConfig*) w->obj(modelSBName);


  // get the modelConfig out of the file
  RooAbsData* data = w->data(dataName);

  // make sure ingredients are found
  if(!data || !sbModel){
    w->Print();
    cout << "data or ModelConfig was not found" <<endl;
    return;
  }
  // make b model
  ModelConfig* bModel = (ModelConfig*) w->obj(modelBName);


   // case of no systematics
   // remove nuisance parameters from model
   if (noSystematics) { 
      const RooArgSet * nuisPar = sbModel->GetNuisanceParameters();
      if (nuisPar && nuisPar->getSize() > 0) { 
         std::cout << "StandardHypoTestInvDemo" << "  -  Switch off all systematics by setting them constant to their initial values" << std::endl;
         RooStats::SetAllConstant(*nuisPar);
      }
      if (bModel) { 

void asimov() {


  std::string InputFile = "./data/example.root";

  // Create the measurement
  Measurement meas("meas", "meas");

  meas.SetOutputFilePrefix( "./results/asimov_UsingC" );
  meas.SetPOI( "SigXsecOverSM" );
  meas.AddConstantParam("alpha_syst1");
  meas.AddConstantParam("Lumi");

  meas.SetLumi( 1.0 );
  meas.SetLumiRelErr( 0.10 );
  meas.SetExportOnly( false );
  meas.SetBinHigh( 2 );

  // Create a channel

  Channel chan( "channel1" );
  chan.SetData( "data", InputFile );
  chan.SetStatErrorConfig( 0.05, "Poisson" );

  
  // Now, create some samples


  // Create the signal sample
  Sample signal( "signal", "signal", InputFile );
  signal.AddOverallSys( "syst1",  0.95, 1.05 );
  signal.AddNormFactor( "SigXsecOverSM", 1, 0, 3 );
  chan.AddSample( signal );

  // Background 1
  Sample background1( "background1", "background1", InputFile );
  background1.ActivateStatError( "background1_statUncert", InputFile );
  background1.AddOverallSys( "syst2", 0.95, 1.05  );
  chan.AddSample( background1 );


  // Background 1
  Sample background2( "background2", "background2", InputFile );
  background2.ActivateStatError();
  background2.AddOverallSys( "syst3", 0.95, 1.05  );
  chan.AddSample( background2 );
  

  // Done with this channel
  // Add it to the measurement:
  meas.AddChannel( chan );

  // Collect the histograms from their files,
  // print some output, 
  meas.CollectHistograms();
  //meas.PrintTree();

  // One can print XML code to an
  // output directory:
  // meas.PrintXML( "xmlFromCCode", meas.GetOutputFilePrefix() );

  RooWorkspace* wspace = RooStats::HistFactory::HistoToWorkspaceFactoryFast::MakeCombinedModel(meas);

  // Get the pdf, obs
  RooAbsPdf* pdf = wspace->pdf("simPdf");
  RooDataSet* obsData = (RooDataSet*) wspace->data("obsData");
  RooDataSet* asimovData = (RooDataSet*) wspace->data("asimovData");

  // fit the pdf to the asimov data
  std::cout << "Fitting to Observed Data" << std::endl;
  pdf->fitTo(*obsData);

  std::cout << "Fitting to Asimov Data" << std::endl;
  pdf->fitTo(*asimovData);

  return;

}

int main (int argc, char **argv)
{
  const char* chInFile = "ws.root";
  const char* chOutFile = "ws_data.root";
  const char* chRootFile = "BDT20.root";
  const char* chCut = "";

  char option_char;
  while ( (option_char = getopt(argc,argv, "i:o:r:c:")) != EOF )
    switch (option_char)
      {
         case 'i': chInFile = optarg; break;
         case 'o': chOutFile = optarg; break;
         case 'r': chRootFile = optarg; break;
         case 'c': chCut = optarg; break;
         case '?': fprintf (stderr,
                            "usage: %s [i<input file> o<output file>]\n", argv[0]);
      }

  cout << "In  Ws = " << chInFile << endl;
  cout << "Out Ws = " << chOutFile << endl;
  cout << "Data From = " << chRootFile << endl;
  cout << "Extra Cut = " << chCut << endl;


   TFile* in_file = new TFile(chInFile);
   RooWorkspace *rws = (RooWorkspace*) in_file->Get("rws");

   TFile* tree_file = new TFile(chRootFile);
   TTree* tree = (TTree*) tree_file->Get("tree");

   TFile* out_file = new TFile(chOutFile, "RECREATE");

   RooArgSet allVars(*rws->var("m"),*rws->var("t"),*rws->var("et"),*rws->var("cpsi"),*rws->var("ctheta"),*rws->var("phi"),*rws->var("d"),*rws->cat("dilution"));
   RooDataSet* data = new RooDataSet("data","data",allVars);
   RooDataSet* dataBkg = new RooDataSet("dataBkg","dataBkg",allVars);

   //TCut* cut = new TCut("5.17<bs_mass && bs_mass<5.57 && bs_pdl>-0.044 && bs_pdl<0.3 ");
   TCut* cut = new TCut("5.17<bs_mass && bs_mass<5.57 && bs_epdl<0.025 && bs_pdl<0.4 && bs_pdl>-0.44");
   *cut += chCut;
   tree->Draw(">>entry_list", *cut, "entrylist");
   TEntryList* event_list = (TEntryList*) out_file->Get("entry_list");

   Double_t dM, dT, dEt, dCpsi, dCtheta, dPhi, dd;
   Int_t ddDefined;
   tree->SetBranchAddress("bs_mass", &dM);
   tree->SetBranchAddress("bs_pdl", &dT);
   tree->SetBranchAddress("bs_epdl", &dEt);
   tree->SetBranchAddress("bs_angle_cpsi", &dCpsi);
   tree->SetBranchAddress("bs_angle_ctheta", &dCtheta);
   tree->SetBranchAddress("bs_angle_phi", &dPhi);
   tree->SetBranchAddress("newtag_ost", &dd);
   tree->SetBranchAddress("newtag_ost_defined", &ddDefined);


   for (Long_t i=0; i<event_list->GetN(); i++){
     tree->GetEntry(event_list->GetEntry(i));

       *rws->var("m")=dM;
       *rws->var("t")=dT/0.0299792458;
       *rws->var("et")=dEt/0.0299792458;
       *rws->var("cpsi")=dCpsi;
       *rws->var("ctheta")=dCtheta;
       *rws->var("phi")=dPhi;

       *rws->var("d")=0;
       rws->cat("dilution")->setIndex(0);
       if ( ddDefined==1 ){
    	   rws->cat("dilution")->setIndex(1);
    	   *rws->var("d")=dd;
       }

       data->add(allVars);
       if (dM<5.29 || dM>5.44)
           dataBkg->add(allVars);
   }

   rws->import(*data);
   rws->import(*dataBkg);
   rws->Write("rws");
   out_file->Close();
   in_file->Close();
   tree_file->Close();

   cout << endl << "Done." << endl;
}

void makeWorkspace2015(RooWorkspace& ws, const TString FileName, struct InputOpt opt){
  double binw=0.05;
 
  std::string finput(FileName);
  TFile *f = new TFile(finput.c_str());
  TTree* theTree       = (TTree*)gROOT->FindObject("myTree"); // OS --- all mass
 
  RooRealVar* mass       = new RooRealVar("invariantMass","#mu#mu mass", opt.dMuon.M.Min, opt.dMuon.M.Max, "GeV/c^{2}");	
  //  ws.import(*mass);
  RooRealVar* dimuPt     = new RooRealVar("dimuPt","p_{T}(#DiMuon)",0,60,"GeV/c");
  RooRealVar* dimuRapidity = new RooRealVar("dimuRapidity",  "dimuRapidity",-2.4, 2.4);
  RooRealVar* vProb      = new RooRealVar("vProb",  "vProb"  ,0.01,1.00);
  RooRealVar* QQsign     = new RooRealVar("QQsign",  "QQsign"  ,-1,5);
  RooRealVar* Centrality = new RooRealVar("Centrality","Centrality",0,200);
  RooRealVar* RunNb      = new RooRealVar("RunNb","RunNb",0,350000);
  RooRealVar* muPlusPt   = new RooRealVar("muPlusPt","muPlusPt", 0, 1000);
  RooRealVar* muMinusPt  = new RooRealVar("muMinusPt","muMinusPt", 0, 1000);
  RooRealVar* muPlusEta  = new RooRealVar("muPlusEta","muPlusEta", -2.4, 2.4);
  RooRealVar* muMinusEta = new RooRealVar("muMinusEta","muMinusEta", -2.4, 2.4);
  RooDataSet* data0, *dataOS, *dataSS;
  RooArgSet cols(*mass,*dimuPt,*dimuRapidity,*muPlusPt,*muMinusPt,*muPlusEta,*muMinusEta, *RunNb, *QQsign);
  data0 = new RooDataSet("data","data",cols); 

  // import the tree to the RooDataSet
  UInt_t          runNb;
  Int_t           centrality;
  ULong64_t       HLTriggers;
  Int_t           Reco_QQ_size;
  Int_t           Reco_QQ_sign[99];   //[Reco_QQ_size]
  TClonesArray    *Reco_QQ_4mom;
  TClonesArray    *Reco_QQ_mupl_4mom;
  TClonesArray    *Reco_QQ_mumi_4mom;
  ULong64_t       Reco_QQ_trig[99];   //[Reco_QQ_size]
  Float_t         Reco_QQ_VtxProb[99];   //[Reco_QQ_size]

  TBranch        *b_runNb;   //!
  TBranch        *b_centrality;   //!
  TBranch        *b_HLTriggers;   //!
  TBranch        *b_Reco_QQ_size;   //!
  TBranch        *b_Reco_QQ_sign;   //!
  TBranch        *b_Reco_QQ_4mom;   //!
  TBranch        *b_Reco_QQ_mupl_4mom;   //!
  TBranch        *b_Reco_QQ_mumi_4mom;   //!
  TBranch        *b_Reco_QQ_trig;   //!
  TBranch        *b_Reco_QQ_VtxProb;   //!

  Reco_QQ_4mom = 0;
  Reco_QQ_mupl_4mom = 0;
  Reco_QQ_mumi_4mom = 0;
  theTree->SetBranchAddress("runNb", &runNb, &b_runNb);
  theTree->SetBranchAddress("Centrality", &centrality, &b_centrality);
  theTree->SetBranchAddress("HLTriggers", &HLTriggers, &b_HLTriggers);
  theTree->SetBranchAddress("Reco_QQ_size", &Reco_QQ_size, &b_Reco_QQ_size);
  theTree->SetBranchAddress("Reco_QQ_sign", Reco_QQ_sign, &b_Reco_QQ_sign);
  theTree->GetBranch("Reco_QQ_4mom")->SetAutoDelete(kFALSE);
  theTree->SetBranchAddress("Reco_QQ_4mom", &Reco_QQ_4mom, &b_Reco_QQ_4mom);
  theTree->GetBranch("Reco_QQ_mupl_4mom")->SetAutoDelete(kFALSE);
  theTree->SetBranchAddress("Reco_QQ_mupl_4mom", &Reco_QQ_mupl_4mom, &b_Reco_QQ_mupl_4mom);
  theTree->GetBranch("Reco_QQ_mumi_4mom")->SetAutoDelete(kFALSE);
  theTree->SetBranchAddress("Reco_QQ_mumi_4mom", &Reco_QQ_mumi_4mom, &b_Reco_QQ_mumi_4mom);
  theTree->SetBranchAddress("Reco_QQ_trig", Reco_QQ_trig, &b_Reco_QQ_trig);
  theTree->SetBranchAddress("Reco_QQ_VtxProb", Reco_QQ_VtxProb, &b_Reco_QQ_VtxProb);

   if (theTree == 0) return;

   Long64_t nentries = theTree->GetEntriesFast();

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      // Long64_t ientry = LoadTree(jentry);
      // if (ientry < 0) break;
      nb = theTree->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;

      for (int i=0; i<Reco_QQ_size; i++) {
         TLorentzVector *qq4mom = (TLorentzVector*) Reco_QQ_4mom->At(i);
         TLorentzVector *qq4mupl = (TLorentzVector*) Reco_QQ_mupl_4mom->At(i);
         TLorentzVector *qq4mumi = (TLorentzVector*) Reco_QQ_mumi_4mom->At(i);
         mass->setVal(qq4mom->M());
         dimuPt->setVal(qq4mom->Pt());
         dimuRapidity->setVal(qq4mom->Rapidity());
         vProb->setVal(Reco_QQ_VtxProb[i]);
         QQsign->setVal(Reco_QQ_sign[i]);
         Centrality->setVal(centrality);
         muPlusPt->setVal(qq4mupl->Pt());
         muMinusPt->setVal(qq4mumi->Pt());
         muPlusEta->setVal(qq4mupl->Eta());
         muMinusEta->setVal(qq4mumi->Eta());
         RunNb->setVal(runNb);

         data0->add(cols);
      }
   }

   TString cut_ap(Form("(%.2f<invariantMass && invariantMass<%.2f) &&"
		       "(%.2f<muPlusEta && muPlusEta < %.2f) &&" 
		       "(%.2f<muMinusEta && muMinusEta < %.2f) &&" 
		       "(%.2f<dimuPt && dimuPt<%.2f) &&"
		       "(abs(dimuRapidity)>%.2f && abs(dimuRapidity)<%.2f)  &&"
		       "(muPlusPt > %.2f && muMinusPt > %.2f) &&"
//.........这里部分代码省略.........

int main(int argc, char *argv[]){
 
  OptionParser(argc,argv);

  TStopwatch sw;
  sw.Start();

  RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
  RooMsgService::instance().setSilentMode(true);

  system("mkdir -p plots/fTest");

  vector<string> procs;
  split(procs,procString_,boost::is_any_of(","));
  
  TFile *inFile = TFile::Open(filename_.c_str());
  RooWorkspace *inWS = (RooWorkspace*)inFile->Get("cms_hgg_workspace");
  
  RooRealVar *mass = (RooRealVar*)inWS->var("CMS_hgg_mass");
  //mass->setBins(320);
  //mass->setRange(mass_-10,mass_+10);
  //mass->setBins(20);
  RooRealVar *MH = new RooRealVar("MH","MH",mass_);
  MH->setVal(mass_);
  MH->setConstant(true);

  map<string,pair<int,int> > choices;
  map<string,vector<RooPlot*> > plotsRV;
  map<string,vector<RooPlot*> > plotsWV;

  for (unsigned int p=0; p<procs.size(); p++){
    vector<RooPlot*> tempRV;
    vector<RooPlot*> tempWV;
    for (int cat=0; cat<ncats_; cat++){
      RooPlot *plotRV = mass->frame(Range(mass_-10,mass_+10));
      plotRV->SetTitle(Form("%s_cat%d_RV",procs[p].c_str(),cat));
      tempRV.push_back(plotRV);
      RooPlot *plotWV = mass->frame(Range(mass_-10,mass_+10));
      plotWV->SetTitle(Form("%s_cat%d_WV",procs[p].c_str(),cat));
      tempWV.push_back(plotWV);
    }
    plotsRV.insert(pair<string,vector<RooPlot*> >(procs[p],tempRV));
    plotsWV.insert(pair<string,vector<RooPlot*> >(procs[p],tempWV));
  }

  vector<int> colors;
  colors.push_back(kBlue);
  colors.push_back(kRed);
  colors.push_back(kGreen+2);
  colors.push_back(kMagenta+1);
  
  for (int cat=0; cat<ncats_; cat++){
    for (unsigned int p=0; p<procs.size(); p++){
      string proc = procs[p];
      RooDataSet *dataRV = (RooDataSet*)inWS->data(Form("sig_%s_mass_m%d_rv_cat%d",proc.c_str(),mass_,cat));
      RooDataSet *dataWV = (RooDataSet*)inWS->data(Form("sig_%s_mass_m%d_wv_cat%d",proc.c_str(),mass_,cat));
      //mass->setBins(160);
      //RooDataHist *dataRV = dataRVtemp->binnedClone();
      //RooDataHist *dataWV = dataWVtemp->binnedClone();
      //RooDataSet *dataRVw = (RooDataSet*)dataRVtemp->reduce(Form("CMS_hgg_mass>=%3d && CMS_hgg_mass<=%3d",mass_-10,mass_+10)); 
      //RooDataSet *dataWVw = (RooDataSet*)dataWVtemp->reduce(Form("CMS_hgg_mass>=%3d && CMS_hgg_mass<=%3d",mass_-10,mass_+10));
      //RooDataHist *dataRV = new RooDataHist(Form("roohist_%s",dataRVtemp->GetName()),Form("roohist_%s",dataRVtemp->GetName()),RooArgSet(*mass),*dataRVtemp);
      //RooDataHist *dataWV = new RooDataHist(Form("roohist_%s",dataWVtemp->GetName()),Form("roohist_%s",dataWVtemp->GetName()),RooArgSet(*mass),*dataWVtemp);
      //RooDataSet *dataRV = stripWeights(dataRVweight,mass);
      //RooDataSet *dataWV = stripWeights(dataWVweight,mass);
      //RooDataSet *data = (RooDataSet*)inWS->data(Form("sig_%s_mass_m%d_cat%d",proc.c_str(),mass_,cat));
     
      int rvChoice=0;
      int wvChoice=0;
      
      // right vertex
      int order=1;
      int prev_order=0;
      int cache_order=0;
      double thisNll=0.;
      double prevNll=1.e6;
      double chi2=0.;
      double prob=0.;

      dataRV->plotOn(plotsRV[proc][cat]);
      //while (prob<0.8) {
      while (order<5) {
        RooAddPdf *pdf = buildSumOfGaussians(Form("cat%d_g%d",cat,order),mass,MH,order);
        RooFitResult *fitRes = pdf->fitTo(*dataRV,Save(true),SumW2Error(true));//,Range(mass_-10,mass_+10));
        double myNll=0.;
        thisNll = fitRes->minNll();
        //double myNll = getMyNLL(mass,pdf,dataRV);
        //thisNll = getMyNLL(mass,pdf,dataRV);
        //RooAbsReal *nll = pdf->createNLL(*dataRV);
        //RooMinuit m(*nll);
        //m.migrad();
        //thisNll = nll->getVal();
        //plot(Form("plots/fTest/%s_cat%d_g%d_rv",proc.c_str(),cat,order),mass_,mass,dataRV,pdf);
        pdf->plotOn(plotsRV[proc][cat],LineColor(colors[order-1]));
        chi2 = 2.*(prevNll-thisNll);
        //if (chi2<0. && order>1) chi2=0.;
        int diffInDof = (2*order+(order-1))-(2*prev_order+(prev_order-1));
        prob = TMath::Prob(chi2,diffInDof);
        cout << "\t RV: " << cat << " " << order << " " << diffInDof << " " << prevNll << " " << thisNll << " " << myNll << " " << chi2 << " " << prob << endl;
        prevNll=thisNll;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  sigalphavar.setConstant(false);   
  
  RooRealVar signvar("signvar","",2.);
  signvar.setConstant(false);     

  RooRealVar sigalpha2var("sigalpha2var","",1.);
  sigalpha2var.setConstant(false);   
  
  RooRealVar sign2var("sign2var","",2.);
  sign2var.setConstant(false);     
  
  
   
  RooArgList tgts;
  RooGBRFunction func("func","",condvars,4);
  RooGBRTarget sigwidtht("sigwidtht","",func,0,sigwidthtvar);
  RooGBRTarget sigmeant("sigmeant","",func,1,sigmeantvar);
  RooGBRTarget signt("signt","",func,2,signvar);
  RooGBRTarget sign2t("sign2t","",func,3,sign2var);
  
  tgts.add(sigwidtht);
  tgts.add(sigmeant);
  tgts.add(signt);
  tgts.add(sign2t);
  
  RooRealConstraint sigwidthlim("sigwidthlim","",sigwidtht,0.0002,0.5);
  RooRealConstraint sigmeanlim("sigmeanlim","",sigmeant,0.2,2.0);
  //RooRealConstraint sigmeanlim("sigmeanlim","",sigmeant,-2.0,-0.2); 
  
  RooRealConstraint signlim("signlim","",signt,1.01,110.); 

  RooRealConstraint sign2lim("sign2lim","",sign2t,1.01,110.); 
  
  RooLinearVar tgtscaled("tgtscaled","",*tgtvar,sigmeanlim,RooConst(0.));
  
  RooDoubleCBFast sigpdf("sigpdf","",tgtscaled,RooConst(1.),sigwidthlim,RooConst(2.0),signlim,RooConst(1.0),sign2lim);
  //RooDoubleCBFast sigpdf("sigpdf","",tgtscaled,RooConst(1.),sigwidthlim,RooConst(2.0),signlim,RooConst(1.0),sign2lim);
  
  //RooCBExp sigpdf("sigpdf","",tgtscaled,RooConst(-1.),sigwidthlim,sigalpha2lim,sign2lim,sigalphalim);
  
  //RooDoubleCBFast sigpdf("sigpdf","",tgtscaled,RooConst(1.),sigwidthlim,RooConst(100.),RooConst(100.),sigalpha2lim,sign2lim);
  //RooDoubleCBFast sigpdf("sigpdf","",tgtscaled,RooConst(1.),sigwidthlim,sigalphalim,signlim,RooConst(3.),sign2lim);
  //RooCBShape sigpdf("sigpdf","",tgtscaled,RooConst(1.),sigwidthlim,sigalphalim,signlim);
  
  RooConstVar etermconst("etermconst","",0.);  
  //RooFormulaVar etermconst("etermconst","","1000.*(@0-1.)*(@0-1.)",RooArgList(tgtscaled));
   
  RooRealVar r("r","",1.);
  r.setConstant();

  std::vector<RooAbsReal*> vpdf;
  vpdf.push_back(&sigpdf);  

  double minweight = 200.;

  std::vector<double> minweights;
  minweights.push_back(minweight);
  
  //ntot.setConstant();

  TFile *fres = new TFile("fres.root","RECREATE");

  if (1) {  
    std::vector<RooAbsData*> vdata;
    vdata.push_back(hdata);    
    
    RooHybridBDTAutoPdf bdtpdfdiff("bdtpdfdiff","",func,tgts,etermconst,r,vdata,vpdf);
    bdtpdfdiff.SetMinCutSignificance(5.);
    bdtpdfdiff.SetPrescaleInit(100);
   // bdtpdfdiff.SetPrescaleInit(10);
    //bdtpdfdiff.SetMaxNSpurious(300.);
    //bdtpdfdiff.SetMaxNSpurious(2400.);
    bdtpdfdiff.SetShrinkage(0.1);
    bdtpdfdiff.SetMinWeights(minweights);
    //bdtpdfdiff.SetMaxNodes(270);
    //bdtpdfdiff.SetMaxNodes(750);
    bdtpdfdiff.SetMaxNodes(500);
    //bdtpdfdiff.SetMaxDepth(8);
    bdtpdfdiff.TrainForest(1e6);  
    
  }   
     
  
  RooWorkspace *wereg = new RooWorkspace("wereg");
  wereg->import(sigpdf);
  
  if (doele && dobarrel)
    wereg->writeToFile("wereg_ele_eb.root");    
  else if (doele && !dobarrel) 
    wereg->writeToFile("wereg_ele_ee.root");    
  else if (!doele && dobarrel)
    wereg->writeToFile("wereg_ph_eb.root");    
  else if (!doele && !dobarrel)
    wereg->writeToFile("wereg_ph_ee.root");    
  
  
  return;
  
  
}

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

  ca->SetGridx();
  ca->SetGridy();
  ca->cd();
  
  gPad->SetLogx();
  gPad->SetObjectStat(1);

  frame->GetYaxis()->SetRangeUser(0,1.05);
  frame->GetXaxis()->SetRangeUser(1,100.);
  frame->GetYaxis()->SetTitle("Efficiency");
  frame->GetXaxis()->SetTitle("E_{T} [GeV]");
  frame->Draw() ;

  frame2->GetYaxis()->SetRangeUser(0,1.05);
  frame2->GetXaxis()->SetRangeUser(1,100.);
  frame2->GetYaxis()->SetTitle("Efficiency");
  frame2->GetXaxis()->SetTitle("E_{T} [GeV]");
  frame2->Draw("same") ;

  TH1F *SCeta1 = new TH1F("SCeta1","SCeta1",50,-2.5,2.5);
  TH1F *SCeta2 = new TH1F("SCeta2","SCeta2",50,-2.5,2.5);

  SCeta1->SetLineColor(color1) ;
  SCeta1->SetMarkerColor(color1);
  SCeta1->SetMarkerStyle(style1);

  SCeta2->SetLineColor(color2) ;
  SCeta2->SetMarkerColor(color2);
  SCeta2->SetMarkerStyle(style2);

  TLegend *leg = new TLegend(0.246,0.435,0.461,0.560,NULL,"brNDC"); // mid : x=353.5
  leg->SetLineColor(1);
  leg->SetTextColor(1);
  leg->SetTextFont(42);
  leg->SetTextSize(0.03);
  leg->SetShadowColor(kWhite);
  leg->SetFillColor(kWhite);
  leg->SetMargin(0.25);
  TLegendEntry *entry=leg->AddEntry("NULL","L1_SingleEG"+names[iEG],"h");
//   leg->AddEntry(SCeta1,name_leg_ecal[iECAL1]+" "+name_leg_coll[iColl1],"p");
//   leg->AddEntry(SCeta2,name_leg_ecal[iECAL2]+" "+name_leg_coll[iColl2],"p");
  leg->AddEntry(SCeta1,name_leg_ecal[iECAL1],"p");
  leg->AddEntry(SCeta2,name_leg_ecal[iECAL2],"p");
  leg->Draw();

  leg = new TLegend(0.16,0.725,0.58,0.905,NULL,"brNDC");
  leg->SetBorderSize(0);
  leg->SetTextFont(62);
  leg->SetTextSize(0.03);
  leg->SetLineColor(0);
  leg->SetLineStyle(1);
  leg->SetLineWidth(1);
  leg->SetFillColor(0);
  leg->SetFillStyle(0);
  leg->AddEntry("NULL","CMS Preliminary 2012 pp  #sqrt{s}=8 TeV","h");
  leg->AddEntry("NULL","#int L dt = "+lumi+"^{-1}","h");
  leg->AddEntry("NULL","Threshold : "+names[iEG]+" GeV","h");
  leg->Draw();

  TPaveText *pt2 = new TPaveText(0.220,0.605,0.487,0.685,"brNDC"); // mid : x=353.5                                          
  pt2->SetLineColor(1);
  pt2->SetTextColor(1);
  pt2->SetTextFont(42);
  pt2->SetTextSize(0.03);
  pt2->SetFillColor(kWhite);
  pt2->SetShadowColor(kWhite);
  pt2->AddText("L1 E/Gamma Trigger");
  pt2->AddText("Electrons from Z");
  pt2->Draw();
  
  //TString name_image="eff_EG20_2012_12fb";

  ca->Print(name_image+".cxx","cxx");
  ca->Print(name_image+".png","png");
  ca->Print(name_image+".gif","gif");
  ca->Print(name_image+".pdf","pdf");
  ca->Print(name_image+".ps","ps");
  ca->Print(name_image+".eps","eps");

  /////////////////////////////
  // SAVE THE ROO FIT RESULT //
  /////////////////////////////

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

  w->import(dataSet);
  w->import(dataSet2);
  
  w->import(*roofitres1,"roofitres1");
  w->import(*roofitres2,"roofitres2");

  cout << "CREATES WORKSPACE : " << endl;
  w->Print();
  
  w->writeToFile(name_image+"_fitres.root") ;
  //gDirectory->Add(w) ;

  //f1->Close();
}

///
/// Fit to the minimum using an improve method.
/// The idea is to kill known minima. For this, we add a multivariate
/// Gaussian to the chi2 at the position of the minimum. Ideally it should
/// be a parabola shaped function, but the Gaussian is readily available.
/// A true parabola like in the original IMPROVE algorithm doesn't work
/// because apparently it affects the other regions of
/// the chi2 too much. There are two parameters to tune: the width of the
/// Gaussian, which is controlled by the error definition of the first fit,
/// and the height of the Gaussian, which is set to 16 (=4 sigma).
/// Three fits are performed: 1) initial fit, 2) fit with improved FCN,
/// 3) fit of the non-improved FCN using the result from step 2 as the start
/// parameters.
/// So far it is only available for the Prob method, via the probimprove
/// command line flag.
///
RooFitResult* Utils::fitToMinImprove(RooWorkspace *w, TString name)
{
	TString parsName = "par_"+name;
	TString obsName  = "obs_"+name;
	TString pdfName  = "pdf_"+name;
	int printlevel = -1;
	RooMsgService::instance().setGlobalKillBelow(ERROR);

	// step 1: find a minimum to start with
	RooFitResult *r1 = 0;
	{
		RooFormulaVar ll("ll", "ll", "-2*log(@0)", RooArgSet(*w->pdf(pdfName)));
		// RooFitResult* r1 = fitToMin(&ll, printlevel);
		RooMinuit m(ll);
		m.setPrintLevel(-2);
		m.setNoWarn();
		m.setErrorLevel(4.0); ///< define 2 sigma errors. This will make the hesse PDF 2 sigma wide!
		int status = m.migrad();
		r1 = m.save();
		// if ( 102<RadToDeg(w->var("g")->getVal())&&RadToDeg(w->var("g")->getVal())<103 )
		// {
		//   cout << "step 1" << endl;
		//   r1->Print("v");
		//   gStyle->SetPalette(1);
		//   float xmin = 0.;
		//   float xmax = 3.14;
		//   float ymin = 0.;
		//   float ymax = 0.2;
		//   TH2F* histo = new TH2F("histo", "histo", 100, xmin, xmax, 100, ymin, ymax);
		//   for ( int ix=0; ix<100; ix++ )
		//   for ( int iy=0; iy<100; iy++ )
		//   {
		//     float x = xmin + (xmax-xmin)*(double)ix/(double)100;
		//     float y = ymin + (ymax-ymin)*(double)iy/(double)100;
		//     w->var("d_dk")->setVal(x);
		//     w->var("r_dk")->setVal(y);
		//     histo->SetBinContent(ix+1,iy+1,ll.getVal());
		//   }
		//   newNoWarnTCanvas("c2");
		//   histo->GetZaxis()->SetRangeUser(0,20);
		//   histo->Draw("colz");
		//   setParameters(w, parsName, r1);
		// }
	}

	// step 2: build and fit the improved fcn
	RooFitResult *r2 = 0;
	{
		// create a Hesse PDF, import both PDFs into a new workspace,
		// so that their parameters are linked
		RooAbsPdf* hessePdf = r1->createHessePdf(*w->set(parsName));
		if ( !hessePdf ) return r1;
		// RooWorkspace* wImprove = (RooWorkspace*)w->Clone();
		RooWorkspace* wImprove = new RooWorkspace();
		wImprove->import(*w->pdf(pdfName));
		wImprove->import(*hessePdf);
		hessePdf = wImprove->pdf(hessePdf->GetName());
		RooAbsPdf* fullPdf = wImprove->pdf(pdfName);

		RooFormulaVar ll("ll", "ll", "-2*log(@0) +16*@1", RooArgSet(*fullPdf, *hessePdf));
		// RooFitResult *r2 = fitToMin(&ll, printlevel);
		RooMinuit m(ll);
		m.setPrintLevel(-2);
		m.setNoWarn();
		m.setErrorLevel(1.0);
		int status = m.migrad();
		r2 = m.save();

		// if ( 102<RadToDeg(w->var("g")->getVal())&&RadToDeg(w->var("g")->getVal())<103 )
		// {
		//   cout << "step 3" << endl;
		//   r2->Print("v");
		//
		//   gStyle->SetPalette(1);
		//   float xmin = 0.;
		//   float xmax = 3.14;
		//   float ymin = 0.;
		//   float ymax = 0.2;
		//   TH2F* histo = new TH2F("histo", "histo", 100, xmin, xmax, 100, ymin, ymax);
		//   for ( int ix=0; ix<100; ix++ )
		//   for ( int iy=0; iy<100; iy++ )
		//   {
		//     float x = xmin + (xmax-xmin)*(double)ix/(double)100;
		//     float y = ymin + (ymax-ymin)*(double)iy/(double)100;
//.........这里部分代码省略.........