//---MODELING---//

	//Polynomial background for noemu
	RooRealVar c1_noemu("c1_noemu", "c1_noemu", c1_poly_noemu, c1_poly_min_noemu, c1_poly_max_noemu);
	RooRealVar c2_noemu("c2_noemu", "c2_noemu", c2_poly_noemu, c2_poly_min_noemu,	c2_poly_max_noemu);
	RooRealVar c3_noemu("c3_noemu", "c3_noemu", c3_poly_noemu, c3_poly_min_noemu,	c3_poly_max_noemu);
	RooPolynomial bkg_noemu("bkg_noemu", "bkg_noemu", e_cm, RooArgList(c1_noemu, c2_noemu, c3_noemu));

	//Polynomial background for emu
	RooRealVar c1_emu("c1_emu", "c1_emu", c1_poly_emu, c1_poly_min_emu, c1_poly_max_emu);
	RooRealVar c2_emu("c2_emu", "c2_emu", c2_poly_emu, c2_poly_min_emu,	c2_poly_max_emu);
	RooRealVar c3_emu("c3_emu", "c3_emu", c3_poly_emu, c3_poly_min_emu,	c3_poly_max_emu);
	RooPolynomial bkg_emu("bkg_emu", "bkg_emu", e_cm, RooArgList(c1_emu, c2_emu, c3_emu));
	
	//Gaussian signal parameters
	RooRealVar m_G("mean_G", "mean_G", mean_G, mean_G_min, mean_G_max);
	RooRealVar s_G("sigma_G", "sigma_G", sigma_G, sigma_G_min, sigma_G_max);
	RooGaussian gauss("gauss", "Generic Gaussian Fit", e_cm, m_G, s_G);

	
	//Lock in constants
	m_G.setVal(mean_G);
	c1_noemu.setVal(-2.892);
	c2_noemu.setVal(-0.08);
	c3_noemu.setVal(10.0);
	
	//Create the branching ratio, common to each channel
	RooRealVar BR_A0("BR_A0","BR_A0",BR_A0_seed, BR_A0_min, BR_A0_max);			//shared parameter
	//Create components for noemu channel model
	
	RooRealVar BR_A0_noemu("BR_A0_noemu","BR_A0_noemu",seed_sig_evts_noemu, min_sig_evts_noemu, max_sig_evts_noemu);
	RooRealVar noemu_BR_coefficient("noemu_BR_coefficient","noemu_BR_coefficient",noemu_BR_coefficient_val);
	RooProduct nsig_noemu("nsig_noemu","nsig_noemu",RooArgSet(BR_A0,noemu_BR_coefficient));													//nsig_noemu common unscaled nsig
	RooRealVar nbkg_noemu("nbkg_noemu","background fraction_noemu",seed_bkg_evts_noemu,min_bkg_evts_noemu,max_bkg_evts_noemu);
	
	//Create components for emu channel model
	RooRealVar BR_A0_emu("BR_A0_emu","BR_A0_emu",seed_sig_evts_emu, min_sig_evts_emu, max_sig_evts_emu);
	RooRealVar emu_BR_coefficient("emu_BR_coefficient","emu_BR_coefficient",emu_BR_coefficient_val);
	RooProduct nsig_emu("nsig_emu","nsig_emu",RooArgSet(BR_A0,emu_BR_coefficient));															//nsig_emu common unscaled nsig
	RooRealVar nbkg_emu("nbkg_emu","background fraction_emu",seed_bkg_evts_emu,min_bkg_evts_emu,max_bkg_evts_emu);
	
	//Create models for noemu and emu channel models
	RooAddPdf model_noemu("model_noemu","model_noemu",RooArgList(gauss,bkg_noemu),RooArgList(nsig_noemu,nbkg_noemu));
	RooAddPdf model_emu("model_emu","model_emu",RooArgList(gauss,bkg_emu),RooArgList(nsig_emu,nbkg_emu));


	//---DATA---//
	//Create index category and joint data sets
	RooCategory sample("sample","sample");
	sample.defineType("noemu");
	sample.defineType("emu");

	//Construct combined dataset in (x,sample)
	RooDataSet combData("combData","combined data",e_cm,Index(sample),Import("noemu",data_noemu),Import("emu",data_emu)) ;


	//---FITTING---//
	//Construct a simultaneous pdf using category sample as index
	RooSimultaneous simPdf("simPdf","simultaneous pdf",sample) ;

	//Make a combined model for the two signal channels 
	simPdf.addPdf(model_noemu,"no_emu");
	simPdf.addPdf(model_emu,"emu");
	simPdf.fitTo(combData);
	

	
	//---PLOTTING DATA---//
	//noemu channel plots
	data_noemu.plotOn(xframe_noemu, Name("data_noemu"), MarkerColor(kBlack));	
	model_noemu.plotOn(xframe_noemu,Components(bkg_noemu),LineColor(kRed), LineStyle(kDashed)) ;
	model_noemu.plotOn(xframe_noemu,Components(gauss),LineColor(kGreen), LineStyle(kDashed)) ;
	model_noemu.plotOn(xframe_noemu, Name("model_noemu"), DrawOption("L"));
	xframe_noemu->Draw();
	
	//emu channel plots
	canvas1->cd(2);
	data_emu.plotOn(xframe_emu, Name("data_emu"), MarkerColor(kBlack));	
	model_emu.plotOn(xframe_emu,Components(bkg_emu),LineColor(kRed), LineStyle(kDashed)) ;
	model_emu.plotOn(xframe_emu,Components(gauss),LineColor(kGreen), LineStyle(kDashed)) ;
	model_emu.plotOn(xframe_emu, Name("model_emu"), DrawOption("L"));
	xframe_emu->Draw();
	
	//Put parameters on for noemu
	canvas1->cd(3);
	RooPlot* xframe3 = e_cm.frame();
	model_noemu.paramOn(xframe3);
	xframe3->SetTitle("noemu Parameters");
	xframe3->Draw();
	
	//Put parameters on for emu
	canvas1->cd(4);
	RooPlot* xframe4 = e_cm.frame();
	model_emu.paramOn(xframe4);
	xframe4->SetTitle("emu Parameters");
	xframe4->Draw();
	
	//--Fit Values--//
	cout<<endl;
	signal_events = nsig_noemu.getVal();
	cout<<"nsig noemu = "<<nsig_noemu.getVal()<<endl;
	cout<<"nsig emu = "<<nsig_emu.getVal()<<endl;
	cout<<"BR_A0 x10^5= "<<BR_A0.getVal()<<endl;
	cout<<endl;
	
	//Plot fit to data
	canvas1->SaveAs(strip+TString("ecm_")+strip2+TString("A0mass.png"));
	
	////////////////////////////
	//UPPER LIMIT CALCULATIONS//
	////////////////////////////

	//BR_A0.setRange(0, 10);		//RESETTING OF RANGE OF BR_A0! Must limit to physical values!
	//BR_A0.setRange(-10, 10);		//debugging
	//BR_A0.setRange(-100, 100);	//debugging
	//BR_A0.setRange(0, 1);			//debugging
	//BR_A0.setRange(0, 100);		//debugging
	BR_A0.setRange(0, 1000);		//debugging
	
	
	//--Bayesean Calculator--//
	RooArgSet POI_Bayes(BR_A0);												//Set BR as POI, option for other POIs
	RooUniform flat_prior("flat_prior","Flat Prior", RooArgSet(BR_A0));		//Make flat prior in BR
	BayesianCalculator bcalc(combData,simPdf,POI_Bayes, flat_prior);		//Call BayesianCalculator, create object
	bcalc.SetConfidenceLevel(.90);											//Set 90% CL
	bcalc.SetLeftSideTailFraction(0);										//Make one-sided upper limit
	SimpleInterval* interval = bcalc.GetInterval();							//Get the interval
	cout<<endl;
	std::cout << "90% CL interval: [ " << interval->LowerLimit() << " - " << interval->UpperLimit() << " ]\n"; //Read out bounds of interval
	cout<<endl;
	UL_Bayes = interval->UpperLimit();										//Read out upper limit
	RooPlot * plot = bcalc.GetPosteriorPlot();								//Make plot of likelihood and integral region.
	TCanvas* canvas4 = new TCanvas("canvas4", "Fourth canvas", 1200, 1000);
	canvas4->cd(1);
	plot->Draw();
	canvas4->SaveAs(strip+TString("ecm_")+strip2+TString("A0mass_BayesPlot.png"));
	
	//---PLC Upper Limit---//
	RooArgSet POI_PLC_UL(BR_A0);
	ProfileLikelihoodCalculator plc_UL(combData, simPdf, POI_PLC_UL);	//Call to PLC with new range on nsig
	plc_UL.SetConfidenceLevel(0.90);
	LikelihoodInterval* plInt_UL = plc_UL.GetInterval();
	//UL_PLC = plInt_UL->UpperLimit(nsig);
	UL_PLC = plInt_UL->UpperLimit(BR_A0);
	cout<<endl;
	std::cout<<"UL_PLC = "<<UL_PLC<<endl;
	cout<<endl;