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

using namespace RooFit ;
using namespace RooPlot ;

void *run()
{
  // 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", -.5, .5),			// center of distribution in x
	ymean("ymean","ymean", -.5, .5),			// center of distribution in y
	sigma("sigma","sigma", 0.0, .2),			// sigma
	nevent("nevent","nevent", Ntot*.9, Ntot*.01, Ntot*2.);	// fit to number of events

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

/*
  // Set region of map to fit (a circle)
   RooFormulaVar xlo("xlo","-sqrt(@0*@0-@1*@1)",RooArgSet(RooConst(radi), y));
   RooFormulaVar xhi("xhi","sqrt(@0*@0-@1*@1)",RooArgSet(RooConst(radi), y));   
   y.setRange("circle", RooConst(-radi), RooConst(radi));
   x.setRange("circle", xlo, xhi);
*/

  // Create test function and generate pdf and test data from it
   TString generatingFnc("exp(-(@0^2+@1^2)/(2*@2^2))");
   RooGenericPdf gauss("gauss","2D Guassian", generatingFnc, RooArgList(x, y, csigma));
   RooDataSet *testdata = gauss.generate(RooArgSet(x,y), Ntot);


  // Create PDF for fitting to distribution
   TString fittingFnc("exp(-((@0-@2)^2+(@1-@3)^2)/(2*@4^2))");   
   RooGenericPdf fitgauss("fitgauss","Fitted 2D Gaussian", fittingFnc, RooArgList(x,y,xmean,ymean,sigma));
   RooExtendPdf  extfitgauss("exfitgauss", "Extended Fitted 2D Gaussian", fitgauss, nevent);


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



//++++++++++++++++++++++++++++++++++++++++++++++++++++++++
//++  Everything beyond here is just printing & plotting
//++++++++++++++++++++++++++++++++++++++++++++++++++++++++

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


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

   testdata->plotOn(xproj);
   testdata->plotOn(yproj);
   extfitgauss.plotOn(xproj);
   extfitgauss.plotOn(yproj);
   

  // Create actual canvas to plot on
   TCanvas *c1 = new TCanvas("c1","c1",1000,500);
   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)));
   TEllipse *eli1 = new TEllipse(0.,0.,radi);		// ellipse showing cuttoff region
   eli1->SetFillStyle(0);
   c1->cd(2);		test_hist->Draw("colz");	eli1->Draw("SAME");
   return ;
}