// Script is used by calling PSFstudy(<see parameters below>)

using namespace RooFit ;
using namespace RooPlot ;

void GaussFit2D_CircularRegion_v2()
{
  // Some default ploting parameters
   gStyle->SetOptStat(0);			// Gets rid of stats box in all plots	
   //gStyle->SetFuncWidth(1);

  // Define some variables to be used throughout the script
   float radi    = 1.7;				// cuttoff distance from camera center
   float range   = 4.;				// distance along one side of plotted map
   float binsize = 0.025;			// bin size
   int   Nbins   = range/binsize;		// number of bins along a side
   int   Ntot = 10000;				// number of events to create


  // Create variables for fitting function (these will be determined in the actual fit)
  // These numbers represent values obtained from an actual fit
   RooRealVar
	x("x","x", -range/2., range/2.),			// x coordinate
	y("y","y", -range/2., range/2.),			// y coordinate
	xmean("xmean","xmean", -1.6, -range/2, range/2),	// center of distribution in x
	ymean("ymean","ymean", 0.0, -range/2, range/2),		// center of distribution in y
	sigma("sigma","sigma", 0.1, 0.0, .15),			// sigma
	nevent("nevent","nevent", Ntot*.9, Ntot*.01, Ntot*2.);	// fit to number of events

   RooConstVar
	csigma("sigma","sigma", .09);				// Sigma for generating pdf


  // Create test function (2D symmetric Gaussian) and generate pdf and test data from it
   TString gaussFnc("exp(-((@0-@2)**2+(@1-@3)**2)/(2*@4**2))");
   RooGenericPdf *gauss = new RooGenericPdf("gauss","2D Guassian", gaussFnc, RooArgList(x,y,xmean,ymean,csigma));
   RooDataSet *testdata = gauss->generate(RooArgSet(x,y), Ntot);


  // Create PDF for fitting to distribution
   RooGenericPdf fitgauss("fitgauss","Fit to 2D Gaussian", gaussFnc, RooArgList(x,y,xmean,ymean,sigma));
   RooExtendPdf extfitgauss("extfitgauss", "Extended Fit to 2D Gaussian", fitgauss, nevent);


  // Set region of map to fit (a circle)
   ostringstream loformula, hiformula;
   loformula << "-sqrt(" << radi*radi << "-@0*@0)";
   hiformula << "sqrt(" << radi*radi << "-@0*@0)";
   RooFormulaVar xlo("xlo",loformula.str().c_str(), y);
   RooFormulaVar xhi("xhi",hiformula.str().c_str(), y);
   x.setRange("circle", xlo, xhi);


  // fit to the distribution
//   RooFitResult *fitresult = extfitgauss.fitTo(*testdata, PrintLevel(1));
   RooFitResult *fitresult = extfitgauss.fitTo(*testdata, PrintLevel(1), Range("circle"));


//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//++  THIS SECTION PRODUCES THE DATA AND PDF PROJECTIONS
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  // Plot projections of generated data and fitted pdf
   RooPlot *xproj = new RooPlot(x, xmean.getVal()-1.0, xmean.getVal()+1.0, 50);
   RooPlot *yproj = new RooPlot(y, ymean.getVal()-1.0, ymean.getVal()+1.0, 50);
   xproj->SetTitle("Projection in X;X;integrated counts");
   yproj->SetTitle("Projection in Y;Y;integrated counts");

   ostringstream datacut;					// define cut range explicitly
   datacut << "(x*x+y*y)<(" << radi*radi << ")";

   testdata->plotOn(xproj, Cut(datacut.str().c_str()));		// works
   testdata->plotOn(yproj, Cut(datacut.str().c_str()));		// works
//  testdata->plotOn(xproj, CutRange("circle"));		// doesnt seem to work
//  testdata->plotOn(yproj, CutRange("circle"));		// doesnt seem to work

   RooFormulaVar ylo("ylo",loformula.str().c_str(), x);		// definition of second circular
   RooFormulaVar yhi("yhi",hiformula.str().c_str(), x);		//    range which is the same as 
   y.setRange("circle2", ylo, yhi);				//    "circle"

//   extfitgauss.plotOn(xproj, ProjectionRange("circle"), Normalization(1.0,RooAbsReal::RelativeExpected));
   extfitgauss.plotOn(xproj, ProjectionRange("circle2"), Normalization(1.0,RooAbsReal::RelativeExpected));
   extfitgauss.plotOn(yproj, ProjectionRange("circle"), Normalization(1.0,RooAbsReal::RelativeExpected));


//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//++ IGNORE EVERYTHING BEYOND THIS POINT. IT IS JUST CODE TO PRODUCE THE PLOTS AND PRINT OUT
//++ THE FITTED PARAMETER VALUES.
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

   extfitgauss.paramOn(xproj, Layout(.7,.95,.95));
   extfitgauss.paramOn(yproj, Layout(.7,.95,.95));
   

  // Create actual canvas to plot on
   TCanvas *c1 = new TCanvas("c1","c1",1500,750);
   TPad *pad1 = new TPad("pad1","pad1", .05, .55, .95, .95);
   TPad *pad2 = new TPad("pad2","pad2", .05, .05, .95, .45);
   c1->Divide(2,1);

   pad1->cd();		gPad->SetLogy();	xproj->Draw();
   pad2->cd();		gPad->SetLogy();	yproj->Draw();
   c1->cd(1);		pad1->Draw();		pad2->Draw();


  // Plot Distribution
   TH2F *test_hist = (TH2F*)testdata->createHistogram("data hist", x,Binning(Nbins), YVar(y,Binning(Nbins)));
   test_hist->SetMinimum(1);
   TEllipse *eli1 = new TEllipse(0.,0.,radi);		// ellipse showing cuttoff region
   eli1->SetFillStyle(0);
   c1->cd(2);		test_hist->Draw("COLZ");	eli1->Draw("SAME");	gPad->SetLogz();


  // Print values so we know how good the fit is
   cout << endl;
   cout << "\txmean : " << xmean.getVal() << " +/- " << xmean.getError() << endl;
   cout << "\tymean : " << ymean.getVal() << " +/- " << ymean.getError() << endl;
   cout << "\tsigma : " << sigma.getVal() << " +/- " << sigma.getError() << endl;
   cout << "\tnevent: " << nevent.getVal() << " +/- " << nevent.getError() << endl;
   cout << endl;

}