C++ RooDataSet::get方法代码示例

Roo2DKeysPdf *SmoothKeys(TH2F *h)
{
  RooDataSet *dataset = th22dataset(h);
  RooRealVar *v1 = (RooRealVar *) dataset->get()->find("constTerm");
  RooRealVar *v2 =  (RooRealVar *) dataset->get()->find("alpha");


  Roo2DKeysPdf *myPdf = new Roo2DKeysPdf("mypdf","", *v1, *v2, *dataset);
  return myPdf;

}

void AnalyzeToy::extract_signal()
{
   calculate_yield();
   MakeSpinSPlot splotter(toyData);
   splotter.addSpecies("signal",ws->pdf("model_signal_mass"),signalYield);
   splotter.addSpecies("background",ws->pdf("model_bkg_mass"),backgroundYield);
   splotter.addVariable(ws->var("mass"));
   splotter.calculate();
   RooDataSet *sweights = splotter.getSWeightDataSet();
   sweights->SetName("sweights");

   RooRealVar weight("weight","",-5.,5.);
   RooArgSet event;
   event.add(*mass);
   event.add(*cosT);
   event.add(weight);

   extractedData = new RooDataSet("extractedData","",event,WeightVar("weight"));

   Long64_t nEntries = toyData->numEntries();
   for(int i=0;i<nEntries;i++)
   {
      double weight_double=0;
      weight_double += sweights->get(i)->getRealValue("signal_sw");
//      weight_double += sweights->get(i)->getRealValue("background_sw");

      mass->setVal(toyData->get(i)->getRealValue("mass"));
      cosT->setVal(toyData->get(i)->getRealValue("cosT"));
      extractedData->add(event,weight_double);
   }
   delete toyData;
}

void glbToId_eta()
{
  TCanvas *myCan=new TCanvas("myCan","myCan");
  myCan->SetGrid();
  TFile *f_MC= new TFile("TnP_GlbToID_MCetaplus_WptTight2012_eta.root","read");
  RooDataSet *datasetMC = (RooDataSet*)f_MC->Get("tpTree/WptTight2012_eta/fit_eff");
  //RooDataSet *datasetMC = (RooDataSet*)f_MC->Get("tpTree/WptTight2012_eta/cnt_eff");
  cout<<"ntry: "<<datasetMC->numEntries()<<endl;
  double XMC[Nbin],XMCerrL[Nbin],XMCerrH[Nbin],YMC[Nbin],YMCerrLo[Nbin],YMCerrHi[Nbin];
  for(int i(0); i<datasetMC->numEntries();i++)
  {
    const RooArgSet &pointMC=*datasetMC->get(i);
    RooRealVar &etaMC=pointMC["eta"],&effMC = pointMC["efficiency"];
    XMC[i]=etaMC.getVal();
    XMCerrL[i]=-etaMC.getAsymErrorLo();
    XMCerrH[i]=etaMC.getAsymErrorHi();
    YMC[i]=effMC.getVal();
    YMCerrLo[i]=-effMC.getAsymErrorLo();
    YMCerrHi[i]=effMC.getAsymErrorHi();
  }
  grMC=new TGraphAsymmErrors(Nbin,XMC,YMC,XMCerrL,XMCerrH,YMCerrLo,YMCerrHi);
  grMC->SetLineColor(kRed);
  grMC->SetMarkerColor(kRed);
  grMC->Draw("AP");
  //grMC->Draw("psame");
  myCan->SaveAs("glbToId_MCplus_eta.png");
  myCan->SaveAs("glbToId_MCplus_eta.eps");

}

void trig_pt()
{
  TCanvas *myCan=new TCanvas("myCan","myCan");
  myCan->SetGrid();
  TFile *f_MC= new TFile("TnP_WptCutToTrig_MCptminus.root","read");
  RooDataSet *datasetMC = (RooDataSet*)f_MC->Get("tpTree/Tnp_WptCut_to_Mu15_eta2p1_pt/cnt_eff");
  cout<<"ntry: "<<datasetMC->numEntries()<<endl;
  double XMC[Nbin],XMCerrL[Nbin],XMCerrH[Nbin],YMC[Nbin],YMCerrLo[Nbin],YMCerrHi[Nbin];
  for(int i(0); i<datasetMC->numEntries();i++)
  {
    const RooArgSet &pointMC=*datasetMC->get(i);
    RooRealVar &ptMC=pointMC["pt"],&effMC = pointMC["efficiency"];
    XMC[i]=ptMC.getVal();
    XMCerrL[i]=-ptMC.getAsymErrorLo();
    XMCerrH[i]=ptMC.getAsymErrorHi();
    YMC[i]=effMC.getVal();
    YMCerrLo[i]=-effMC.getAsymErrorLo();
    YMCerrHi[i]=effMC.getAsymErrorHi();
  }
  grMC=new TGraphAsymmErrors(11,XMC,YMC,XMCerrL,XMCerrH,YMCerrLo,YMCerrHi);
  grMC->SetLineColor(kRed);
  grMC->SetMarkerColor(kRed);
  grMC->Draw("AP");
  //grMC->Draw("psame");
  myCan->SaveAs("trig_McMinus_pt.png");
  myCan->SaveAs("trig_McMinus_pt.eps");

}

void sc_wptCut_P_et()
{
  TCanvas *myCan=new TCanvas("myCan","myCan");
  myCan->SetGrid();
  /************************
  TFile *f_RD= new TFile("TnP_Z_Trigger_RDpt.root","read");
  RooDataSet *dataset = (RooDataSet*)f_RD->Get("tpTree/Track_To_TightCombRelIso_Mu15_eta2p1_pt/fit_eff");
  cout<<"ntry: "<<dataset->numEntries()<<endl;
  double X[11],XerrL[11],XerrH[11],Y[11],YerrLo[11],YerrHi[11];
  for(int i(0); i<dataset->numEntries();i++)
  {
    const RooArgSet &point=*dataset->get(i);
    RooRealVar &pt=point["pt"],&eff = point["efficiency"];
    X[i]=pt.getVal();
    XerrL[i]=-pt.getAsymErrorLo();
    XerrH[i]=pt.getAsymErrorHi();
    Y[i]=eff.getVal();
    YerrLo[i]=-eff.getAsymErrorLo();
    YerrHi[i]=eff.getAsymErrorHi();
  }
  gr=new TGraphAsymmErrors(11,X,Y,XerrL,XerrH,YerrLo,YerrHi);
  gr->Draw("AP");
***************************/
  TFile *f_MC= new TFile("efficiency-mc-SCToPfElectron_et_P.root","read");
  RooDataSet *datasetMC = (RooDataSet*)f_MC->Get("SuperClusterToPFElectron/SCtoWptCut_efficiency/cnt_eff");
  //RooDataSet *datasetMC = (RooDataSet*)f_MC->Get("tpTree/Track_with_TightCombRelIso_to_Mu15_eta2p1_pt/fit_eff");
  cout<<"ntry: "<<datasetMC->numEntries()<<endl;

  double XMC[binSize],XMCerrL[binSize],XMCerrH[binSize],YMC[binSize],YMCerrLo[binSize],YMCerrHi[binSize];
  for(int i(0); i<datasetMC->numEntries();i++)
  {
    const RooArgSet &pointMC=*datasetMC->get(i);
    RooRealVar &ptMC=pointMC["probe_sc_et"],&effMC = pointMC["efficiency"];
    XMC[i]=ptMC.getVal();
    XMCerrL[i]=-ptMC.getAsymErrorLo();
    XMCerrH[i]=ptMC.getAsymErrorHi();
    YMC[i]=effMC.getVal();
    YMCerrLo[i]=-effMC.getAsymErrorLo();
    YMCerrHi[i]=effMC.getAsymErrorHi();
  }
  grMC=new TGraphAsymmErrors(binSize,XMC,YMC,XMCerrL,XMCerrH,YMCerrLo,YMCerrHi);
  grMC->SetLineColor(kRed);
  grMC->SetMarkerColor(kRed);
  //myCan->SetLogx();
  grMC->Draw("AP");
  //grMC->Draw("psame");
  TLine *myLine=new TLine(25,0,25,1);
  myLine->Draw("same");
  myCan->SaveAs("sc_wptCut_P_et.png");
  myCan->SaveAs("sc_wptCut_P_et.eps");

}

void wspaceread_backgrounds(int channel = 1)
{

	gSystem->AddIncludePath("-I$ROOFITSYS/include");
	gROOT->ProcessLine(".L ~/tdrstyle.C");
	
	string schannel;
	if (channel == 1) schannel = "4mu";
	if (channel == 2) schannel = "4e";
	if (channel == 3) schannel = "2mu2e";
	std::cout << "schannel = " << schannel << std::endl;

	// R e a d   w o r k s p a c e   f r o m   f i l e
	// -----------------------------------------------
	
	// Open input file with workspace (generated by rf14_wspacewrite)
	char infile[192];
	sprintf(infile,"/scratch/hep/ntran/dataFiles/HZZ4L/datasets/datasets/%s/ZZAnalysisTree_ZZTo4L_lowmass.root",schannel.c_str());
	TFile *f = new TFile(infile) ;
	char outfile[192];
	sprintf( outfile, "figs/pdf_%s_bkg_highmass.eps", schannel.c_str() );
	//f->ls();
	
	
	RooDataSet* set = (RooDataSet*) f->Get("data");
	RooArgSet* obs = set->get() ;
	obs->Print();
	RooRealVar* CMS_zz4l_mass = (RooRealVar*) obs->find("CMS_zz4l_mass") ;
	
	for (int i=0 ; i<set->numEntries() ; i++) { 
		set->get(i) ; 
		//cout << CMS_zz4l_mass->getVal() << " = " << set->weight() << endl ; 
	} 
	
	gSystem->Load("PDFs/RooqqZZPdf_cxx.so");
	//gSystem->Load("PDFs/RooggZZPdf_cxx.so");

	// LO contribution

	//RooRealVar m4l("m4l","m4l",100.,1000.);
	RooRealVar a1("a1","a1",224.,100.,1000.);
	RooRealVar a2("a2","a2",-209.,-1000.,1000.);
	RooRealVar a3("a3","a3",121.,20.,1000.);
	RooRealVar a4("a4","a4",-0.022,-10.,10.);
	RooRealVar b1("b1","b1",181.,100.,1000.);
	RooRealVar b2("b2","b2",707.,0.,1000.);
	RooRealVar b3("b3","b3",60.,20.,1000.);
	RooRealVar b4("b4","b4",0.04,-10.,10.);
	RooRealVar b5("b5","b5",5.,0.,1000.);
	RooRealVar b6("b6","b6",0.,-10.,10.);
	RooRealVar frac_bkg("frac_bkg","frac_bkg",0.5,0.,1.);
	
	//a1.setConstant(kTRUE);
	//a2.setConstant(kTRUE);
	//a3.setConstant(kTRUE);
	//a4.setConstant(kTRUE);
	//b1.setConstant(kTRUE);
	//b2.setConstant(kTRUE);
	//b3.setConstant(kTRUE);
	//b4.setConstant(kTRUE);
	//b5.setConstant(kTRUE);
	//b6.setConstant(kTRUE);

	RooqqZZPdf bkg_qqzz("bkg_qqzz","bkg_qqzz",*CMS_zz4l_mass,a1,a2,a3,a4,b1,b2,b3,b4,b5,b6,frac_bkg);

	RooFitResult *r = bkg_qqzz.fitTo( *set, SumW2Error(kTRUE) );//, Save(kTRUE), SumW2Error(kTRUE)) ;

	// Plot Y
	RooPlot* frameM4l = CMS_zz4l_mass->frame(Title("M4L"),Bins(100)) ;
	set->plotOn(frameM4l) ;
	bkg_qqzz.plotOn(frameM4l) ;
	
	TCanvas *c = new TCanvas("c","c",800,600);
	c->cd();
	frameM4l->Draw();
	
	
	
	/*
	// Retrieve workspace from file
	RooWorkspace* w = (RooWorkspace*) f->Get("workspace") ;
	
	w->Print();
	
	///*
	RooRealVar* CMS_zz4l_mass = w->var("CMS_zz4l_mass") ;
	RooAbsPdf* background_nonorm = w->pdf("background_nonorm") ;
	//RooAbsData* backgroundData = w->data("backgroundData") ;
	RooAbsData* data_bkg_red = w->data("data_bkg_red") ;
	
	RooArgSet* obs = data_bkg_red->get() ; 
	RooRealVar* xdata = obs->find(CMS_zz4l_mass.GetName()) ; 
	for (int i=0 ; i<data_bkg_red->numEntries() ; i++) { 
		data_bkg_red->get(i) ; 
		cout << xdata->getVal() << " = " << data_bkg_red->weight() << endl ; 
	} 
	std::cout << "nEntries = " << data_bkg_red->numEntries() << std::endl;
	obs->Print();

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

void wspaceread_signals2e2mu(int channel = 3)
{

	gSystem->AddIncludePath("-I$ROOFITSYS/include");
	gROOT->ProcessLine(".L ~/tdrstyle.C");
	setTDRStyle();
	//gSystem->Load("PDFs/RooRelBW1_cxx.so");
	//gSystem->Load("PDFs/RooRelBW2_cxx.so");
	gSystem->Load("../PDFs/HZZ4LRooPdfs_cc.so");

	string schannel;
	if (channel == 1) schannel = "4mu";
	if (channel == 2) schannel = "4e";
	if (channel == 3) schannel = "2mu2e";
	std::cout << "schannel = " << schannel << std::endl;

	const int nPoints = 17.;
	int masses[nPoints] = {120,130,140,150,160,170,180,190,200,250,300,350,400,450,500,550,600};
	double mHVal[nPoints] = {120,130,140,150,160,170,180,190,200,250,300,350,400,450,500,550,600};
	double widths[nPoints] = {3.48e-03,4.88e-03,8.14e-03,1.73e-02,8.30e-02,3.80e-01,6.31e-01,1.04e+00,1.43e+00,4.04e+00,8.43e+00,1.52e+01,2.92e+01,46.95,6.80e+01,93.15,1.23e+02};
	// R e a d   w o r k s p a c e   f r o m   f i l e
	// -----------------------------------------------
	
	double a_meanBW[nPoints];
	double a_gammaBW[nPoints];
	double a_meanCB[nPoints];
	double a_sigmaCB[nPoints];
	double a_alphaCB[nPoints];
	double a_nCB[nPoints];
		
	for (int i = 0; i < nPoints; i++){
	//for (int i = 0; i < 1; i++){
		
		// Open input file with workspace (generated by rf14_wspacewrite)
		char infile[192];
		sprintf(infile,"/scratch/hep/ntran/dataFiles/HZZ4L/datasets/datasets_baseline/%s/ZZAnalysisTree_H%i%s.root",schannel.c_str(),masses[i],schannel.c_str());
		TFile *f = new TFile(infile) ;
		char outfile[192];
		sprintf( outfile, "figs/pdf_%s_bkg_highmass.eps", schannel.c_str() );
		//f->ls();
		
		double windowVal = max( widths[i], 1. );
		if (mHVal[i] >= 275){ lowside = 180.; }
		else { lowside = 100.; }
		double low_M = max( (mHVal[i] - 20.*windowVal), lowside) ;
		double high_M = min( (mHVal[i] + 15.*windowVal), 900.) ;

		//double windowVal = max( widths[i], 1.);
		//double windowVal = max ( widths[i], 1. );
		//low_M = max( (mHVal[i] - 25.*windowVal), 100.) ;
		//high_M = min( (mHVal[i] + 20.*windowVal), 1000.) ;
		//low_M = max( (mHVal[i] - 15.*windowVal), 100.) ;
		//high_M = min( (mHVal[i] + 10.*windowVal), 1000.) ;
		std::cout << "lowM = " << low_M << ", highM = " << high_M << std::endl;
			
		RooDataSet* set = (RooDataSet*) f->Get("data");
		RooArgSet* obs = set->get() ;
		obs->Print();
		RooRealVar* CMS_zz4l_mass = (RooRealVar*) obs->find("CMS_zz4l_mass") ;
		CMS_zz4l_mass->setRange(low_M,high_M);
		for (int a=0 ; a<set->numEntries() ; a++) { 
			set->get(a) ; 
			//cout << CMS_zz4l_mass->getVal() << " = " << set->weight() << endl ; 
		} 
		
		// constraining parameters...
		double l_sigmaCB = 0., s_sigmaCB = 3.;
		if (mHVal[i] >= 500.){ l_sigmaCB = 10.; s_sigmaCB = 12.; }
		
		double s_n_CB = 2.6+(-1.1/290.)*(mHVal[i]-110.);
		if (mHVal[i] >= 400){ s_n_CB = 1.5; }

		
		RooRealVar mean_CB("mean_CB","mean_CB",0.,-25.,25);
		RooRealVar sigma_CB("sigma_CB","sigma_CB",s_sigmaCB,l_sigmaCB,30.);
		RooRealVar alpha_CB("alpha_CB","alpha_CB",0.95,0.8,1.2);
		RooRealVar n_CB("n_CB","n_CB",s_n_CB,1.5,2.8);
		RooCBShape signalCB("signalCB","signalCB",*CMS_zz4l_mass,mean_CB,sigma_CB,alpha_CB,n_CB);
		
		RooRealVar mean_BW("mean_BW","mean_BW", mHVal[i] ,100.,1000.);
		RooRealVar gamma_BW("gamma_BW","gamma_BW",widths[i],0.,200.);
		//RooBreitWigner signalBW("signalBW", "signalBW",*CMS_zz4l_mass,mean_BW,gamma_BW);
		//RooRelBW1 signalBW("signalBW", "signalBW",*CMS_zz4l_mass,mean_BW,gamma_BW);
		
		RooRelBWUF signalBW("signalBW", "signalBW",*CMS_zz4l_mass,mean_BW);
		//RooRelBW1 signalBW("signalBW", "signalBW",*CMS_zz4l_mass,mean_BW,gamma_BW);
		RooBreitWigner signalBW1("signalBW1", "signalBW1",*CMS_zz4l_mass,mean_BW,gamma_BW);
		RooRelBW1 signalBW2("signalBW2", "signalBW2",*CMS_zz4l_mass,mean_BW,gamma_BW);
		
		//Set #bins to be used for FFT sampling to 10000
		CMS_zz4l_mass->setBins(100000,"fft") ;

		//Construct BW (x) CB
		RooFFTConvPdf* sig_ggH = new RooFFTConvPdf("sig_ggH","BW (X) CB",*CMS_zz4l_mass,signalBW,signalCB, 2);
		// Buffer fraction for cyclical behavior
		sig_ggH->setBufferFraction(0.2);
		
		mean_BW.setConstant(kTRUE);
		gamma_BW.setConstant(kTRUE);
		n_CB.setConstant(kTRUE);
//.........这里部分代码省略.........

void LeptonPreselectionCMG( PreselType type, RooWorkspace * w ) {
	const Options & opt = Options::getInstance(); 
	if (type == ELE)
		cout << "Running Electron Preselection :" << endl;
	else if (type == MU)
		cout << "Running Muon Preselection :" << endl;
	else if (type == EMU)
		cout << "Running Electron-Muon Preselection() ..." << endl;
	else if (type == PHOT)
		cout << "Running Photon Preselection :" << endl;

	string systVar;
	try {
		systVar = opt.checkStringOption("SYSTEMATIC_VAR");
	} catch (const std::string & exc) {
		cout << exc << endl;
	}
	if (systVar == "NONE")
		systVar.clear();

#ifdef CMSSWENV
	JetCorrectionUncertainty jecUnc("Summer13_V4_MC_Uncertainty_AK5PFchs.txt");
#endif

	string inputDir = opt.checkStringOption("INPUT_DIR");
	string outputDir = opt.checkStringOption("OUTPUT_DIR");
	string sampleName = opt.checkStringOption("SAMPLE_NAME");
	string inputFile = inputDir + '/' + sampleName + ".root";
	cout << "\tInput file: " << inputFile << endl;

	bool isSignal = opt.checkBoolOption("SIGNAL");
	TGraph * higgsW = 0;
	TGraph * higgsI = 0;
	if (isSignal) {
		double higgsM = opt.checkDoubleOption("HIGGS_MASS");
		if (higgsM >= 400) {
			string dirName = "H" + double2string(higgsM);
			bool isVBF = opt.checkBoolOption("VBF");
			string lshapeHistName = "cps";
			string intHistName = "nominal";
			
			if (systVar == "LSHAPE_UP") {
				intHistName = "up";
			} else if (systVar == "LSHAPE_DOWN") {
				intHistName = "down";
			}

			if (isVBF) {
				TFile weightFile("VBF_LineShapes.root");
				higgsW = (TGraph *) ( (TDirectory *) weightFile.Get(dirName.c_str()))->Get( lshapeHistName.c_str() )->Clone();

			} else {
				TFile weightFile("GG_LineShapes.root");
				higgsW = (TGraph *) ( (TDirectory *) weightFile.Get(dirName.c_str()))->Get( lshapeHistName.c_str() )->Clone();
				TFile interfFile("newwgts_interf.root");
				higgsI = (TGraph *) ( (TDirectory *) interfFile.Get(dirName.c_str()))->Get( intHistName.c_str() )->Clone();
			}
		}
	}

	TFile * file = new TFile( inputFile.c_str() );
	if (!file->IsOpen())
		throw string("ERROR: Can't open the file: " + inputFile + "!");
	TDirectory * dir = (TDirectory *) file->Get("dataAnalyzer");
	TH1D * nEvHisto = (TH1D *) dir->Get("cutflow");
	TH1D * puHisto = (TH1D *) dir->Get("pileup");
	TTree * tree = ( TTree * ) dir->Get( "data" );
	Event ev( tree );
	const int * runP = ev.getSVA<int>("run"); 
	const int * lumiP = ev.getSVA<int>("lumi"); 
	const int * eventP = ev.getSVA<int>("event"); 
	const bool * trigBits = ev.getAVA<bool>("t_bits");
	const int * trigPres = ev.getAVA<int>("t_prescale");
	const float * metPtA = ev.getAVA<float>("met_pt");
	const float * metPhiA = ev.getAVA<float>("met_phi");
	const float * rhoP = ev.getSVA<float>("rho");
	const float * rho25P = ev.getSVA<float>("rho25");
	const int * nvtxP = ev.getSVA<int>("nvtx"); 
	const int * niP = ev.getSVA<int>("ngenITpu"); 
	
#ifdef PRINTEVENTS
	string eventFileName;
	if (type == ELE)
		eventFileName = "events_ele.txt";
	else if (type == MU)
		eventFileName = "events_mu.txt";
	else if (type == EMU)
		eventFileName = "events_emu.txt";

	EventPrinter evPrint(ev, type, eventFileName);
	evPrint.readInEvents("diff.txt");
	evPrint.printElectrons();
	evPrint.printMuons();
	evPrint.printZboson();
	evPrint.printJets();
	evPrint.printHeader();
#endif

	string outputFile = outputDir + '/' + sampleName;

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

//.........这里部分代码省略.........
  // We can use PROOF to speed things along in parallel
  ProofConfig pc(*w, 4, "",false); 
  if(mc->GetGlobalObservables()){
    cout << "will use global observables for unconditional ensemble"<<endl;
    mc->GetGlobalObservables()->Print();
    toymcsampler->SetGlobalObservables(*mc->GetGlobalObservables());
  }
  toymcsampler->SetProofConfig(&pc);	// enable proof


  // Now get the interval
  PointSetInterval* interval = fc.GetInterval();
  ConfidenceBelt* belt = fc.GetConfidenceBelt();
 
  // print out the iterval on the first Parameter of Interest
  cout << "\n95% interval on " <<firstPOI->GetName()<<" is : ["<<
    interval->LowerLimit(*firstPOI) << ", "<<
    interval->UpperLimit(*firstPOI) <<"] "<<endl;

  // get observed UL and value of test statistic evaluated there
  RooArgSet tmpPOI(*firstPOI);
  double observedUL = interval->UpperLimit(*firstPOI);
  firstPOI->setVal(observedUL);
  double obsTSatObsUL = fc.GetTestStatSampler()->EvaluateTestStatistic(*data,tmpPOI);


  // Ask the calculator which points were scanned
  RooDataSet* parameterScan = (RooDataSet*) fc.GetPointsToScan();
  RooArgSet* tmpPoint;

  // make a histogram of parameter vs. threshold
  TH1F* histOfThresholds = new TH1F("histOfThresholds","",
				    parameterScan->numEntries(),
				    firstPOI->getMin(),
				    firstPOI->getMax());
  histOfThresholds->GetXaxis()->SetTitle(firstPOI->GetName());
  histOfThresholds->GetYaxis()->SetTitle("Threshold");

  // loop through the points that were tested and ask confidence belt
  // what the upper/lower thresholds were.
  // For FeldmanCousins, the lower cut off is always 0
  for(Int_t i=0; i<parameterScan->numEntries(); ++i){
    tmpPoint = (RooArgSet*) parameterScan->get(i)->clone("temp");
    double arMax = belt->GetAcceptanceRegionMax(*tmpPoint);
    double poiVal = tmpPoint->getRealValue(firstPOI->GetName()) ;
    histOfThresholds->Fill(poiVal,arMax);
  }
  TCanvas* c1 = new TCanvas();
  c1->Divide(2);
  c1->cd(1);
  histOfThresholds->SetMinimum(0);
  histOfThresholds->Draw();
  c1->cd(2);

  /////////////////////////////////////////////////////////////
  // Now we generate the expected bands and power-constriant
  ////////////////////////////////////////////////////////////

  // First: find parameter point for mu=0, with conditional MLEs for nuisance parameters
  RooAbsReal* nll = mc->GetPdf()->createNLL(*data);
  RooAbsReal* profile = nll->createProfile(*mc->GetParametersOfInterest());
  firstPOI->setVal(0.);
  profile->getVal(); // this will do fit and set nuisance parameters to profiled values
  RooArgSet* poiAndNuisance = new RooArgSet();
  if(mc->GetNuisanceParameters())
    poiAndNuisance->add(*mc->GetNuisanceParameters());

// The actual job
void backgroundFits_qqzz_1Dw(int channel, int sqrts, int VBFtag)
{
  if(sqrts==7)return;
  TString schannel;
  if      (channel == 1) schannel = "4mu";
  else if (channel == 2) schannel = "4e";
  else if (channel == 3) schannel = "2e2mu";
  else cout << "Not a valid channel: " << schannel << endl;

  TString ssqrts = (long) sqrts + TString("TeV");

  cout << "schannel = " << schannel << "  sqrts = " << sqrts << " VBFtag = " << VBFtag << endl;

  TString outfile;
  if(VBFtag<2) outfile = "CardFragments/qqzzBackgroundFit_" + ssqrts + "_" + schannel + "_" + Form("%d",int(VBFtag)) + ".txt";
  if(VBFtag==2) outfile = "CardFragments/qqzzBackgroundFit_" + ssqrts + "_" + schannel + ".txt";
  ofstream of(outfile,ios_base::out);
  of << "### background functions ###" << endl;


  gSystem->AddIncludePath("-I$ROOFITSYS/include");
  gROOT->ProcessLine(".L ../CreateDatacards/include/tdrstyle.cc");
  setTDRStyle(false);
  gStyle->SetPadLeftMargin(0.16);

  TString filepath;
  if (sqrts==7) {
    filepath = filePath7TeV;
  } else if (sqrts==8) {
    filepath = filePath8TeV;
  }

  TChain* tree = new TChain("SelectedTree");
  tree->Add( filepath+ "/" + (schannel=="2e2mu"?"2mu2e":schannel) + "/HZZ4lTree_ZZTo*.root");


  RooRealVar* MC_weight = new RooRealVar("MC_weight","MC_weight",0.,2.) ; 
  RooRealVar* ZZMass = new RooRealVar("ZZMass","ZZMass",100.,1000.);
  RooRealVar* NJets30 = new RooRealVar("NJets30","NJets30",0.,100.);
  RooArgSet ntupleVarSet(*ZZMass,*NJets30,*MC_weight);
  RooDataSet *set = new RooDataSet("set","set",ntupleVarSet,WeightVar("MC_weight"));

  Float_t myMC,myMass;
  Short_t myNJets;
  int nentries = tree->GetEntries();

  tree->SetBranchAddress("ZZMass",&myMass);
  tree->SetBranchAddress("MC_weight",&myMC);
  tree->SetBranchAddress("NJets30",&myNJets);

  for(int i =0;i<nentries;i++) {
    tree->GetEntry(i);
    if(VBFtag==1 && myNJets<2)continue;
    if(VBFtag==0 && myNJets>1)continue;

    ntupleVarSet.setRealValue("ZZMass",myMass);
    ntupleVarSet.setRealValue("MC_weight",myMC);
    ntupleVarSet.setRealValue("NJets30",(double)myNJets);

    set->add(ntupleVarSet, myMC);
  }

  double totalweight = 0.;
  double totalweight_z = 0.;
  for (int i=0 ; i<set->numEntries() ; i++) { 
    //set->get(i) ; 
    RooArgSet* row = set->get(i) ;
    //row->Print("v");
    totalweight += set->weight();
    if (row->getRealValue("ZZMass") < 200) totalweight_z += set->weight();
  } 
  cout << "nEntries: " << set->numEntries() << ", totalweight: " << totalweight << ", totalweight_z: " << totalweight_z << endl;

  gSystem->Load("libHiggsAnalysisCombinedLimit.so");
	
  //// ---------------------------------------
  //Background
  RooRealVar CMS_qqzzbkg_a0("CMS_qqzzbkg_a0","CMS_qqzzbkg_a0",115.3,0.,200.);
  RooRealVar CMS_qqzzbkg_a1("CMS_qqzzbkg_a1","CMS_qqzzbkg_a1",21.96,0.,200.);
  RooRealVar CMS_qqzzbkg_a2("CMS_qqzzbkg_a2","CMS_qqzzbkg_a2",122.8,0.,200.);
  RooRealVar CMS_qqzzbkg_a3("CMS_qqzzbkg_a3","CMS_qqzzbkg_a3",0.03479,0.,1.);
  RooRealVar CMS_qqzzbkg_a4("CMS_qqzzbkg_a4","CMS_qqzzbkg_a4",185.5,0.,200.);
  RooRealVar CMS_qqzzbkg_a5("CMS_qqzzbkg_a5","CMS_qqzzbkg_a5",12.67,0.,200.);
  RooRealVar CMS_qqzzbkg_a6("CMS_qqzzbkg_a6","CMS_qqzzbkg_a6",34.81,0.,100.);
  RooRealVar CMS_qqzzbkg_a7("CMS_qqzzbkg_a7","CMS_qqzzbkg_a7",0.1393,0.,1.);
  RooRealVar CMS_qqzzbkg_a8("CMS_qqzzbkg_a8","CMS_qqzzbkg_a8",66.,0.,200.);
  RooRealVar CMS_qqzzbkg_a9("CMS_qqzzbkg_a9","CMS_qqzzbkg_a9",0.07191,0.,1.);
  RooRealVar CMS_qqzzbkg_a10("CMS_qqzzbkg_a10","CMS_qqzzbkg_a10",94.11,0.,200.);
  RooRealVar CMS_qqzzbkg_a11("CMS_qqzzbkg_a11","CMS_qqzzbkg_a11",-5.111,-100.,100.);
  RooRealVar CMS_qqzzbkg_a12("CMS_qqzzbkg_a12","CMS_qqzzbkg_a12",4834,0.,10000.);
  RooRealVar CMS_qqzzbkg_a13("CMS_qqzzbkg_a13","CMS_qqzzbkg_a13",0.2543,0.,1.);
	
  if (channel == 1){
    ///* 4mu
    CMS_qqzzbkg_a0.setVal(103.854);
    CMS_qqzzbkg_a1.setVal(10.0718);
    CMS_qqzzbkg_a2.setVal(117.551);
    CMS_qqzzbkg_a3.setVal(0.0450287);
    CMS_qqzzbkg_a4.setVal(185.262);
//.........这里部分代码省略.........

void MakePlots(RooWorkspace* ws){

  // Here we make plots of the discriminating variable (invMass) after the fit
  // and of the control variable (isolation) after unfolding with sPlot.
  std::cout << "make plots" << std::endl;

  // make our canvas
  TCanvas* cdata = new TCanvas("sPlot","sPlot demo", 400, 600);
  cdata->Divide(1,3);

  // get what we need out of the workspace
  RooAbsPdf* model = ws->pdf("model");
  RooAbsPdf* zModel = ws->pdf("zModel");
  RooAbsPdf* qcdModel = ws->pdf("qcdModel");

  RooRealVar* isolation = ws->var("isolation");
  RooRealVar* invMass = ws->var("invMass");

  // note, we get the dataset with sWeights
  RooDataSet* data = (RooDataSet*) ws->data("dataWithSWeights");

  // this shouldn't be necessary, need to fix something with workspace
  // do this to set parameters back to their fitted values.
  model->fitTo(*data, Extended() );

  //plot invMass for data with full model and individual componenets overlayed
  //  TCanvas* cdata = new TCanvas();
  cdata->cd(1);
  RooPlot* frame = invMass->frame() ; 
  data->plotOn(frame ) ; 
  model->plotOn(frame) ;   
  model->plotOn(frame,Components(*zModel),LineStyle(kDashed), LineColor(kRed)) ;   
  model->plotOn(frame,Components(*qcdModel),LineStyle(kDashed),LineColor(kGreen)) ;   
    
  frame->SetTitle("Fit of model to discriminating variable");
  frame->Draw() ;
 

  // Now use the sWeights to show isolation distribution for Z and QCD.  
  // The SPlot class can make this easier, but here we demonstrait in more
  // detail how the sWeights are used.  The SPlot class should make this 
  // very easy and needs some more development.

  // Plot isolation for Z component.  
  // Do this by plotting all events weighted by the sWeight for the Z component.
  // The SPlot class adds a new variable that has the name of the corresponding
  // yield + "_sw".
  cdata->cd(2);

  // create weightfed data set 
  RooDataSet * dataw_z = new RooDataSet(data->GetName(),data->GetTitle(),data,*data->get(),0,"zYield_sw") ;

  RooPlot* frame2 = isolation->frame() ; 
  dataw_z->plotOn(frame2, DataError(RooAbsData::SumW2) ) ; 
    
  frame2->SetTitle("isolation distribution for Z");
  frame2->Draw() ;

  // Plot isolation for QCD component.  
  // Eg. plot all events weighted by the sWeight for the QCD component.
  // The SPlot class adds a new variable that has the name of the corresponding
  // yield + "_sw".
  cdata->cd(3);
  RooDataSet * dataw_qcd = new RooDataSet(data->GetName(),data->GetTitle(),data,*data->get(),0,"qcdYield_sw") ;
  RooPlot* frame3 = isolation->frame() ; 
  dataw_qcd->plotOn(frame3,DataError(RooAbsData::SumW2) ) ; 
    
  frame3->SetTitle("isolation distribution for QCD");
  frame3->Draw() ;

  //  cdata->SaveAs("SPlot.gif");

}

//.........这里部分代码省略.........
  //RooRealVar *m1 = new RooRealVar("m1","m1",5320,5380);
  //RooRealVar *s1 = new RooRealVar("s1","s1",1,20);
  //RooGaussian *sig = new RooGaussian("sig","sig",*mass,*m1,*s1);
  RooRealVar *m2 = new RooRealVar("m2","m2",5320,5380);
  RooRealVar *s2 = new RooRealVar("s2","s2",1,20);
  RooGaussian *sig_bd = new RooGaussian("sig_bd","sig_bd",*mass,*m2,*s2);

  //
  RooRealVar *bs2kstkst_l       = new RooRealVar( "bs2kstkst_l"    ,"", -5, -20, -1.);
  RooConstVar *bs2kstkst_zeta   = new RooConstVar( "bs2kstkst_zeta","",0.              );
  RooConstVar *bs2kstkst_fb     = new RooConstVar( "bs2kstkst_fb"  ,"",0.              );
  RooRealVar *bs2kstkst_sigma   = new RooRealVar( "bs2kstkst_sigma","",15    ,10   ,20 );
  RooRealVar *bs2kstkst_mu      = new RooRealVar( "bs2kstkst_mu"   ,"",5350  ,5380     );
  RooRealVar *bs2kstkst_a       = new RooRealVar( "bs2kstkst_a"    ,"",2.5  , 0    ,10 );
  RooRealVar *bs2kstkst_n       = new RooRealVar( "bs2kstkst_n"    ,"",2.5  , 0    ,10 );
  RooRealVar *bs2kstkst_a2      = new RooRealVar( "bs2kstkst_a2"   ,"",2.5  , 0    ,10 );
  RooRealVar *bs2kstkst_n2      = new RooRealVar( "bs2kstkst_n2"   ,"",2.5  , 0    ,10 );

  RooIpatia2 *sig = new RooIpatia2("sig","sig",*mass,*bs2kstkst_l,*bs2kstkst_zeta,*bs2kstkst_fb,*bs2kstkst_sigma,*bs2kstkst_mu,*bs2kstkst_a,*bs2kstkst_n,*bs2kstkst_a2,*bs2kstkst_n2);

  RooRealVar *bkg_y = new RooRealVar("bkg_y","bkg_y",10e3,10e5);
  RooRealVar *sig_y = new RooRealVar("sig_y","sig_y",0,20e3);
  RooRealVar *sig_bd_y = new RooRealVar("sig_bd_y","sig_bd_y",0,3000);

  RooArgList *pdfs = new RooArgList();
  RooArgList *yields = new RooArgList();

  pdfs->add( *exp );
  pdfs->add( *sig );
  pdfs->add( *sig_bd );

  yields->add( *bkg_y );
  yields->add( *sig_y );
  yields->add( *sig_bd_y );

  RooAddPdf *pdf = new RooAddPdf("pdf","pdf",*pdfs,*yields);

  pdf->fitTo(*Data, Extended() );

  RooPlot *plot = mass->frame();
  Data->plotOn(plot);
    // set fit params constant;
  pdf->plotOn(plot);

  TCanvas *c = new TCanvas();
  plot->Draw();
  c->Print("tmp/mass.pdf");

  // Plots Kst Ms with no sweights
  TCanvas *c1 = new TCanvas("c1","c1",800,1200);
  c1->Divide(1,2);
  c1->cd(1);
  RooPlot *c1p1 = w->var("B_s0_DTF_KST1_M")->frame();
  Data->plotOn(c1p1);
  c1p1->Draw();
  c1->cd(2);
  RooPlot *c1p2 = w->var("B_s0_DTF_KST2_M")->frame();
  Data->plotOn(c1p2);
  c1p2->Draw();
  c1->Print("tmp/nosw.pdf");

  // set fit params constant
  p1->setConstant(true);
  //m1->setConstant(true);
  //s1->setConstant(true);
  bs2kstkst_l->setConstant(true);
  //bs2kstkst_zeta->setConstant(true);
  //bs2kstkst_fb->setConstant(true);
  bs2kstkst_sigma->setConstant(true);
  bs2kstkst_mu->setConstant(true);
  bs2kstkst_a->setConstant(true);
  bs2kstkst_n->setConstant(true);
  bs2kstkst_a2->setConstant(true);
  bs2kstkst_n2->setConstant(true);
  m2->setConstant(true);
  s2->setConstant(true);

  RooStats::SPlot *sData = new RooStats::SPlot("sData","sData", *Data, pdf, *yields);

  w->import(*sData);
  w->import(*Data,Rename("Data_wsweights"));

  RooDataSet *swdata = new RooDataSet("Data_wsweights", "Data", Data, *Data->get(), 0 , "sig_y_sw");
  // Plots Kst Ms with no sweights
  TCanvas *c2 = new TCanvas("c2","c2",800,1200);
  c2->Divide(1,2);
  c2->cd(1);
  RooPlot *c2p1 = w->var("B_s0_DTF_KST1_M")->frame();
  swdata->plotOn(c2p1);
  c2p1->Draw();
  c2->cd(2);
  RooPlot *c2p2 = w->var("B_s0_DTF_KST2_M")->frame();
  swdata->plotOn(c2p2);
  c2p2->Draw();
  c2->Print("tmp/withsw.pdf");


  tf->Close();
  return 0;
}

// The actual job
void backgroundFits_ggzz_1Dw(int channel, int sqrts, int VBFtag)
{
  TString schannel;
  if      (channel == 1) schannel = "4mu";
  else if (channel == 2) schannel = "4e";
  else if (channel == 3) schannel = "2e2mu";
  else cout << "Not a valid channel: " << schannel << endl;

  TString ssqrts = (long) sqrts + TString("TeV");

  cout << "schannel = " << schannel << "  sqrts = " << sqrts << " VBFtag = "<< VBFtag << endl;

  TString outfile;
  outfile = "CardFragments/ggzzBackgroundFit_" + ssqrts + "_" + schannel + "_" + Form("%d",int(VBFtag)) + ".txt";
  ofstream of(outfile,ios_base::out);

  gSystem->AddIncludePath("-I$ROOFITSYS/include");
  gROOT->ProcessLine(".L ../CreateDatacards/include/tdrstyle.cc");
  setTDRStyle(false);
  gStyle->SetPadLeftMargin(0.16);
	
  TString filepath;filepath.Form("AAAOK/ZZ%s/ZZ4lAnalysis.root",schannel.Data());
  TFile *f = TFile::Open(filepath);
  TTree *tree = f->Get("ZZTree/candTree");

  RooRealVar* MC_weight = new RooRealVar("MC_weight","MC_weight",0.,2.) ; 
  RooRealVar* ZZMass = new RooRealVar("ZZMass","ZZMass",100,100.,1000.);
  RooRealVar* NJets30 = new RooRealVar("NJets30","NJets30",0.,5.);
  RooArgSet ntupleVarSet(*ZZMass,*NJets30,*MC_weight);
  RooDataSet *set = new RooDataSet("set","set",ntupleVarSet,WeightVar("MC_weight"));
  //RooArgSet ntupleVarSet(*ZZMass,*NJets30);  
  //RooDataSet *set = new RooDataSet("set","set",ntupleVarSet);

  Float_t myMC,myMass;
  Int_t myNJets;
  int nentries = tree->GetEntries();

  Float_t myPt,myJetPt,myJetEta,myJetPhi,myJetMass,myFisher;
  Int_t myExtralep,myBJets;
  tree->SetBranchAddress("ZZMass",&myMass);
  tree->SetBranchAddress("genHEPMCweight",&myMC);
  tree->SetBranchAddress("nCleanedJetsPt30",&myNJets);
  tree->SetBranchAddress("ZZPt",&myPt);
  tree->SetBranchAddress("nExtraLep",&myExtralep);
  tree->SetBranchAddress("nCleanedJetsPt30BTagged",&myBJets);
  tree->SetBranchAddress("DiJetDEta",&myFisher);

  for(int i =0;i<nentries;i++) {
    tree->GetEntry(i);
    if(myMass<100.)continue;
    int cat = category(myExtralep,myPt, myMass,myNJets, myBJets,/* jetpt, jeteta, jetphi, jetmass,*/myFisher);
    if(VBFtag != cat )continue;

    ntupleVarSet.setRealValue("ZZMass",myMass);
    ntupleVarSet.setRealValue("MC_weight",myMC);
    ntupleVarSet.setRealValue("NJets30",(double)cat);

    set->add(ntupleVarSet, myMC);
  }

  //RooRealVar* ZZLD = new RooRealVar("ZZLD","ZZLD",0.,1.);
  //char cut[10];
  //sprintf(cut,"ZZLD>0.5");
  //RooDataSet* set = new RooDataSet("set","set",tree,RooArgSet(*ZZMass,*MC_weight,*ZZLD),cut,"MC_weight");

  double totalweight = 0.;
  for (int i=0 ; i<set->numEntries() ; i++) { 
    set->get(i) ; 
    totalweight += set->weight();
    //cout << CMS_zz4l_mass->getVal() << " = " << set->weight() << endl ; 
  } 
  cout << "nEntries: " << set->numEntries() << ", totalweight: " << totalweight << endl;
		
  gSystem->Load("libHiggsAnalysisCombinedLimit.so");
	
  //// ---------------------------------------
  //Background
  RooRealVar CMS_qqzzbkg_a0("CMS_qqzzbkg_a0","CMS_qqzzbkg_a0",115.3,0.,200.);
  RooRealVar CMS_qqzzbkg_a1("CMS_qqzzbkg_a1","CMS_qqzzbkg_a1",21.96,0.,200.);
  RooRealVar CMS_qqzzbkg_a2("CMS_qqzzbkg_a2","CMS_qqzzbkg_a2",122.8,0.,200.);
  RooRealVar CMS_qqzzbkg_a3("CMS_qqzzbkg_a3","CMS_qqzzbkg_a3",0.03479,0.,1.);
  RooRealVar CMS_qqzzbkg_a4("CMS_qqzzbkg_a4","CMS_qqzzbkg_a4",185.5,0.,200.);
  RooRealVar CMS_qqzzbkg_a5("CMS_qqzzbkg_a5","CMS_qqzzbkg_a5",12.67,0.,200.);
  RooRealVar CMS_qqzzbkg_a6("CMS_qqzzbkg_a6","CMS_qqzzbkg_a6",34.81,0.,100.);
  RooRealVar CMS_qqzzbkg_a7("CMS_qqzzbkg_a7","CMS_qqzzbkg_a7",0.1393,0.,1.);
  RooRealVar CMS_qqzzbkg_a8("CMS_qqzzbkg_a8","CMS_qqzzbkg_a8",66.,0.,200.);
  RooRealVar CMS_qqzzbkg_a9("CMS_qqzzbkg_a9","CMS_qqzzbkg_a9",0.07191,0.,1.);
	
  RooggZZPdf_v2* bkg_ggzz = new RooggZZPdf_v2("bkg_ggzz","bkg_ggzz",*ZZMass,
					      CMS_qqzzbkg_a0,CMS_qqzzbkg_a1,CMS_qqzzbkg_a2,CMS_qqzzbkg_a3,CMS_qqzzbkg_a4,
					      CMS_qqzzbkg_a5,CMS_qqzzbkg_a6,CMS_qqzzbkg_a7,CMS_qqzzbkg_a8,CMS_qqzzbkg_a9);
	
  //// ---------------------------------------
	
  RooFitResult *r1 = bkg_ggzz->fitTo( *set, Save(kTRUE), SumW2Error(kTRUE) );//, Save(kTRUE), SumW2Error(kTRUE)) ;

  cout << endl;
  cout << "------- Parameters for " << schannel << " sqrts=" << sqrts << endl;
  cout << "  a0_bkgd = " << CMS_qqzzbkg_a0.getVal() << endl;
//.........这里部分代码省略.........

void makeEffToys(Int_t seed, TString veto="D") {
//	r.SetSeed(seed);
	r = RooRandom::randomGenerator();
	r->SetSeed(seed);

	TFile* fin = TFile::Open(veto+"veto_200.root");

	TString fName("toys/"); fName+=seed; fName+="/"+veto+"veto_200.root";

	TFile* fout = new TFile(fName,"RECREATE");

	for(Int_t j=0; j<21; ++j) {
		TString hName  = "efficiency_"; hName+=j;
		TString hName2 = "efficiencyHist_"; hName2+=j;
		TEfficiency* hin = dynamic_cast<TEfficiency*>(fin->Get(hName));

		TH1* hout = dynamic_cast<TH1*>(hin->GetTotalHistogram()->Clone(hName2));

		std::vector<Int_t> corrBins;

		Int_t n = hout->GetNbinsX();
		for(Int_t i=0; i<n; ++i) {
			Double_t eff = hin->GetEfficiency(i+1);
			Double_t erm = hin->GetEfficiencyErrorLow(i+1);
			Double_t erp = hin->GetEfficiencyErrorUp(i+1);

			Bool_t fluctuate = kTRUE;

			// don't fluctuate if the veto hasn't affected this bin
			// also ignore the odd missing entry - not sure what causes these but they don't seem reasonable
			if((eff > 0.99 && eff + erp > 0.999) //efficiency close to 1 and not significantly different
			|| ((i<1 || hin->GetEfficiency(i) == 1) && (i>n-1 || hin->GetEfficiency(i+2) == 1))) { //single bin dip (careful with this one)
				eff = 1;
				fluctuate = kFALSE;
			}
			if(eff < 0.01 && eff - erm < 0.001) {//efficiency close to 0 and not significantly different
				eff = 0;
				fluctuate = kFALSE;
			}

			//otherwise we're fluctuating the bin
			if(fluctuate) {
				//if the errors are roughly symmetric then we can symmetrise them and introduce some correlation between neighbouring bins
				//this is difficult to do with asymmetric errors so if asymmetry > 10% lets just ignore correlations
				//note that a large asymmetry in neighbouring bins will also lead to same-sign fluctuations anyway
				if((erm - erp) / (erm + erp) < 0.1) {
					//correlation is more important than asymmetry

					//add bin to the list to be fluctuated later
					corrBins.push_back(i+1);

				} else {
					//asymmetry is more important than correlation

					//first catch any cases on a limit (the previous checks for eff > 0.99 and eff < 0.01 should catch these but play it safe)
					if(erm <= 0) {
						//vary with a half Gaussian
						eff += TMath::Abs(r->Gaus(0.,erp));
					} else if(erp <= 0) {
						//vary with a half Gaussian
						eff -= TMath::Abs(r->Gaus(0.,erm));
					} else {
						//vary with a bifurcated Gaussian
						RooRealVar effVar( "effVar", "",-1.,2.);
						RooRealVar muVar(  "muVar",  "",eff);
						RooRealVar sigmVar("sigmVar","",erm);
						RooRealVar sigpVar("sigpVar","",erp);

						RooBifurGauss pdf("pdf","",effVar,muVar,sigmVar,sigpVar);
						RooDataSet* ds = pdf.generate(RooArgSet(effVar),1);
						eff = ds->get(0)->getRealValue("effVar");
						delete ds;
					}
				}
			}
			if(eff > 1.0) eff = 1.0; //std::cout << i << "\t" << eff << "\t" << erp << "\t" << erm << std::endl;
//			std::cout << hin->GetEfficiency(i+1) << "\t" << eff << std::endl;

			hout->SetBinContent(i+1, eff);
		}

		//now deal with the correlated efficiencies
		Double_t corrFactor(0.01);

		while(!doCorrelatedBinFluctuation(hin,hout,corrBins,corrFactor)) {
			corrFactor /= 2.;
		}

//		std::cout << std::endl;

		TCanvas c;
		hin->Draw();
		hout->SetMarkerColor(kRed);
		hout->SetMarkerStyle(4);
		hout->Draw("Psame");
		TString pName = "plots/toys/"; pName+=seed; pName+="/"+veto+"veto_Q"; pName+=j; pName+=".pdf";
		c.SaveAs(pName);
	}

	hout->Write();
//.........这里部分代码省略.........

///
/// Perform the 1d Prob scan.
/// Saves chi2 values and the prob-Scan p-values in a root tree
/// For the datasets stuff, we do not yet have a MethodDatasetsProbScan class, so we do it all in
/// MethodDatasetsProbScan
/// \param nRun Part of the root tree file name to facilitate parallel production.
///
int MethodDatasetsProbScan::scan1d(bool fast, bool reverse)
{
	if (fast) return 0; // tmp

	if ( arg->debug ) cout << "MethodDatasetsProbScan::scan1d() : starting ... " << endl;

    // Set limit to all parameters.
    this->loadParameterLimits(); /// Default is "free", if not changed by cmd-line parameter


    // Define scan parameter and scan range.
    RooRealVar *parameterToScan = w->var(scanVar1);
    float parameterToScan_min = hCL->GetXaxis()->GetXmin();
    float parameterToScan_max = hCL->GetXaxis()->GetXmax();

		// do a free fit
		RooFitResult *result = this->loadAndFit(this->pdf); // fit on data
		assert(result);
    RooSlimFitResult *slimresult = new RooSlimFitResult(result,true);
		slimresult->setConfirmed(true);
		solutions.push_back(slimresult);
		double freeDataFitValue = w->var(scanVar1)->getVal();

    // Define outputfile
    system("mkdir -p root");
    TString probResName = Form("root/scan1dDatasetsProb_" + this->pdf->getName() + "_%ip" + "_" + scanVar1 + ".root", arg->npoints1d);
    TFile* outputFile = new TFile(probResName, "RECREATE");

    // Set up toy root tree
    this->probScanTree = new ToyTree(this->pdf, arg);
    this->probScanTree->init();
    this->probScanTree->nrun = -999; //\todo: why does this branch even exist in the output tree of the prob scan?

    // Save parameter values that were active at function
    // call. We'll reset them at the end to be transparent
    // to the outside.
    RooDataSet* parsFunctionCall = new RooDataSet("parsFunctionCall", "parsFunctionCall", *w->set(pdf->getParName()));
    parsFunctionCall->add(*w->set(pdf->getParName()));

    // start scan
    cout << "MethodDatasetsProbScan::scan1d_prob() : starting ... with " << nPoints1d << " scanpoints..." << endl;
    ProgressBar progressBar(arg, nPoints1d);
    for ( int i = 0; i < nPoints1d; i++ )
    {
        progressBar.progress();
        // scanpoint is calculated using min, max, which are the hCL x-Axis limits set in this->initScan()
        // this uses the "scan" range, as expected
        // don't add half the bin size. try to solve this within plotting method

        float scanpoint = parameterToScan_min + (parameterToScan_max - parameterToScan_min) * (double)i / ((double)nPoints1d - 1);
				if (arg->debug) cout << "DEBUG in MethodDatasetsProbScan::scan1d_prob() " << scanpoint << " " << parameterToScan_min << " " << parameterToScan_max << endl;

        this->probScanTree->scanpoint = scanpoint;

        if (arg->debug) cout << "DEBUG in MethodDatasetsProbScan::scan1d_prob() - scanpoint in step " << i << " : " << scanpoint << endl;

        // don't scan in unphysical region
        // by default this means checking against "free" range
        if ( scanpoint < parameterToScan->getMin() || scanpoint > parameterToScan->getMax() + 2e-13 ) {
            cout << "it seems we are scanning in an unphysical region: " << scanpoint << " < " << parameterToScan->getMin() << " or " << scanpoint << " > " << parameterToScan->getMax() + 2e-13 << endl;
            exit(EXIT_FAILURE);
        }

        // FIT TO REAL DATA WITH FIXED HYPOTHESIS(=SCANPOINT).
        // THIS GIVES THE NUMERATOR FOR THE PROFILE LIKELIHOOD AT THE GIVEN HYPOTHESIS
        // THE RESULTING NUISANCE PARAMETERS TOGETHER WITH THE GIVEN HYPOTHESIS ARE ALSO
        // USED WHEN SIMULATING THE TOY DATA FOR THE FELDMAN-COUSINS METHOD FOR THIS HYPOTHESIS(=SCANPOINT)
        // Here the scanvar has to be fixed -> this is done once per scanpoint
        // and provides the scanner with the DeltaChi2 for the data as reference
        // additionally the nuisances are set to the resulting fit values

        parameterToScan->setVal(scanpoint);
        parameterToScan->setConstant(true);

        RooFitResult *result = this->loadAndFit(this->pdf); // fit on data
        assert(result);

        if (arg->debug) {
            cout << "DEBUG in MethodDatasetsProbScan::scan1d_prob() - minNll data scan at scan point " << scanpoint << " : " << 2 * result->minNll() << ": "<< 2 * pdf->getMinNll() << endl;
        }
        this->probScanTree->statusScanData = result->status();

        // set chi2 of fixed fit: scan fit on data
        // CAVEAT: chi2min from fitresult gives incompatible results to chi2min from pdf
        // this->probScanTree->chi2min           = 2 * result->minNll();
        this->probScanTree->chi2min           = 2 * pdf->getMinNll();
        this->probScanTree->covQualScanData   = result->covQual();
        this->probScanTree->scanbest  = freeDataFitValue;

        // After doing the fit with the parameter of interest constrained to the scanpoint,
        // we are now saving the fit values of the nuisance parameters. These values will be
        // used to generate toys according to the PLUGIN method.
        this->probScanTree->storeParsScan(); // \todo : figure out which one of these is semantically the right one
//.........这里部分代码省略.........