void ProjectSieveHoleH(void) {

	TFile *_filetmp = TFile::Open(
			"ArchivedRootFiles/SieveCentralLHRS/SieveData.f51_1170.Sieve.root");
	assert(_filetmp);

	TCanvas * result = _filetmp->GetObjectChecked("SieveCheck", "TCanvas");
	assert(result);

	result->Draw();

	TH2D
			* Sieve_Foil0 =
					new TH2D(
							*(TH2D *) ((TPad *) (result->GetListOfPrimitives()->Last())->FindObject(
									"Sieve_Foil0")));
	assert(Sieve_Foil0);
	//	Sieve_Foil0 -> Draw("COLZ");


	TH2D * hHori = new TH2D(*Sieve_Foil0);
	//	hHori -> Draw("COLZ");

	const Double_t vlim1 = -0.004, vlim2 = 0.0065;// center hole
	//	const Double_t vlim1 = -0.029, vlim2 = -0.023;// 2nd big hole
	//	const Double_t vlim1 = 0.0137, vlim2 = 0.0171;// row 5


	for (int i = 1; i <= hHori->GetNbinsY(); i++) {
		for (int j = 1; j <= hHori->GetNbinsX(); j++) {
			const double binc = hHori->GetYaxis()->GetBinCenter(i);

			if (binc < vlim1 || binc > vlim2) {
				//				cout << i << binc << endl;
				hHori -> SetBinContent(j, i, 0);
			}
		}
	}

	hHori -> Draw("COLZ");
	TH1D* hHoriProj = hHori -> ProjectionX("hHoriProj");

	c1 = new TCanvas("ProjectSieveHole", "ProjectSieveHole", 800, 500);
	//	c1 -> Divid(2,2);
	//	int idx=1;

	//	c1->cd(idx++);
	hHoriProj -> Draw();
}

void ProjectSieveHoleV(void) {

	TFile
			*_filetmp =
					TFile::Open(
							"ArchivedRootFiles/SieveCentralLHRS/db_L.vdc.dat.21.Ph.PlusJulyData.replay/SieveData.f51_2Pass.Sieve.root");
	assert(_filetmp);

	TCanvas * result = _filetmp->GetObjectChecked("SieveCheck", "TCanvas");
	assert(result);

	result->Draw();

	TH2D
			* Sieve_Foil0 =
					new TH2D(
							*(TH2D *) ((TPad *) (result->GetListOfPrimitives()->Last())->FindObject(
									"Sieve_Foil0")));
	assert(Sieve_Foil0);
	//	Sieve_Foil0 -> Draw("COLZ");


	TH2D * hHori = new TH2D(*Sieve_Foil0);
	//	hHori -> Draw("COLZ");

	//	const Double_t vlim1 = -0.0001, vlim2 = 0.0026;// center hole
	const Double_t vlim1 = 0.006, vlim2 = 0.0085;// col 4


	for (int i = 1; i <= hHori->GetNbinsX(); i++) {
		for (int j = 1; j <= hHori->GetNbinsY(); j++) {
			const double binc = hHori->GetXaxis()->GetBinCenter(i);

			if (binc < vlim1 || binc > vlim2) {
				//				cout << i << binc << endl;
				hHori -> SetBinContent(i, j, 0);
			}
		}
	}

	//	hHori -> Draw("COLZ");
	TH1D* hHoriProj = hHori -> ProjectionY("hHoriProj");

	c1 = new TCanvas("ProjectSieveHoleV", "ProjectSieveHoleV", 800, 500);
	//	c1 -> Divid(2,2);
	//	int idx=1;

	//	c1->cd(idx++);
	//	hHori -> Draw("COLZ");
	hHoriProj -> Draw();
}

void ProjectSieveHoleRHRS2Pass(void) {

	// --------------------------------------------------------------------
	TFile
			*_filetmp =
					TFile::Open(
							"ArchivedRootFiles/SieveCentralRHRS/db_R.vdc.dat.8.Ph.SortTerms.replay//SieveData.f51_2Pass_sieveout.Sieve.root");
	assert(_filetmp);

	TCanvas * result = _filetmp->GetObjectChecked("SieveCheck", "TCanvas");
	assert(result);

	result->Draw();

	TH2D
			* Sieve_Foil0_sieveout =
					new TH2D(
							*(TH2D *) ((TPad *) (result->GetListOfPrimitives()->Last())->FindObject(
									"Sieve_Foil0")));
	assert(Sieve_Foil0_sieveout);
	//	Sieve_Foil0 -> Draw("COLZ");


	TH1D* hHoriProj_sieveout = Sieve_Foil0_sieveout -> ProjectionX("hHoriProj");
	// --------------------------------------------------------------------
	TFile
			*_filetmp =
					TFile::Open(
							"ArchivedRootFiles/SieveCentralRHRS/db_R.vdc.dat.8.Ph.SortTerms.replay//SieveData.f51_2Pass_sievehalfin.Sieve.root");
	assert(_filetmp);

	TCanvas * result = _filetmp->GetObjectChecked("SieveCheck", "TCanvas");
	assert(result);

	result->Draw();

	TH2D
			* Sieve_Foil0_sievehalfin =
					new TH2D(
							*(TH2D *) ((TPad *) (result->GetListOfPrimitives()->Last())->FindObject(
									"Sieve_Foil0")));
	assert(Sieve_Foil0_sievehalfin);
	//	Sieve_Foil0 -> Draw("COLZ");


	TH1D* hHoriProj_sievehalfin = Sieve_Foil0_sievehalfin -> ProjectionX(
			"hHoriProj");
	// --------------------------------------------------------------------
	TFile
			*_filetmp =
					TFile::Open(
							"ArchivedRootFiles/SieveCentralRHRS/db_R.vdc.dat.8.Ph.SortTerms.replay//SieveData.f51_2Pass_sievein.Sieve.root");
	assert(_filetmp);

	TCanvas * result = _filetmp->GetObjectChecked("SieveCheck", "TCanvas");
	assert(result);

	result->Draw();

	TH2D
			* Sieve_Foil0_sievein =
					new TH2D(
							*(TH2D *) ((TPad *) (result->GetListOfPrimitives()->Last())->FindObject(
									"Sieve_Foil0")));
	assert(Sieve_Foil0_sievein);
	//	Sieve_Foil0 -> Draw("COLZ");


	TH1D* hHoriProj_sievein = Sieve_Foil0_sievein -> ProjectionX("hHoriProj");

	hHoriProj_sieveout -> SetLineColor(kRed);
	hHoriProj_sievehalfin -> SetLineColor(kGreen);
	hHoriProj_sievein -> SetLineColor(kBlue);

	hHoriProj_sievein -> SetYTitle("Arbitrary Unit");
	hHoriProj_sieveout -> Scale(1. / hHoriProj_sieveout->GetMaximum());
	hHoriProj_sievehalfin -> Scale(1. / hHoriProj_sievehalfin->GetMaximum());
	hHoriProj_sievein -> Scale(1. / hHoriProj_sievein->GetMaximum());

	//	hHoriProj_sievein -> SetYTitle("Yield / Coulomb of beam");
	//	hHoriProj_sieveout -> Scale (1./0.02818);
	//	hHoriProj_sievehalfin -> Scale (1./0.02214);
	//	hHoriProj_sievein -> Scale (1./(0.01998+0.02));

	hHoriProj_sievein -> SetTitle(
			"Horizontal Sieve Plane Projection: Sieve in, half-in, out");
	hHoriProj_sievein -> SetXTitle("X [m], +X = hori towards beamline");

	c1 = new TCanvas("ProjectSieveHoleRHRS2Pass", "ProjectSieveHoleRHRS2Pass",
			800, 800);
	//	c1 -> Divid(2,2);
	//	int idx=1;

	//	c1->cd(idx++);
	//	hHori -> Draw("COLZ");
	hHoriProj_sievein -> Draw();
	hHoriProj_sievehalfin -> Draw("same");
	hHoriProj_sieveout -> Draw("same");

	TLegend *l1 = new TLegend(0.1, 0.7, 0.8, 0.85);
	l1->AddEntry(hHoriProj_sieveout, "Sieve Out", "l");
	l1->AddEntry(hHoriProj_sievehalfin, "Sieve Half-In", "l");
	l1->AddEntry(hHoriProj_sievein, "Sieve In", "l");

	l1 -> Draw();
}

void ProjectSieveHoleH_LHRSCenterAngle(void) {

	TFile
			*_filetmp =
					TFile::Open(
							"Production5/db_L.vdc.dat.Prod5.1.2GeV.sieve.root");
	assert(_filetmp);

	TCanvas * result = _filetmp->GetObjectChecked("SieveCheck",	"TCanvas");
	assert(result);

	result->Draw();

	TPad * pad = (TPad *) (result->GetListOfPrimitives()->FindObject("SieveCheck_4"));
	assert(pad);

	TH2D * Sieve_Foil0 = (TH2D *) (pad->GetListOfPrimitives()->FindObject(
			"Sieve_Foil3"));

	assert(Sieve_Foil0);
	Sieve_Foil0 -> Draw("COLZ");

	TH2D * hHori = new TH2D(*Sieve_Foil3);
	//	hHori -> Draw("COLZ");

	const Double_t vlim1 = -0.015, vlim2 = 0.006;// center hole
	//	const Double_t vlim1 = -0.029, vlim2 = -0.023;// 2nd big hole
	//	const Double_t vlim1 = 0.0137, vlim2 = 0.0171;// row 5


	for (int i = 1; i <= hHori->GetNbinsY(); i++) {
		for (int j = 1; j <= hHori->GetNbinsX(); j++) {
			const double binc = hHori->GetYaxis()->GetBinCenter(i);

			if (binc < vlim1 || binc > vlim2) {
				//				cout << i << binc << endl;
				hHori -> SetBinContent(j, i, 0);
			}
		}
	}
	//
	hHori -> Draw("COLZ");
	TH1D* hHoriProj = hHori -> ProjectionX("hHoriProj");

	c1 = new TCanvas("ProjectSieveHole", "ProjectSieveHole", 800, 500);
	//	c1 -> Divid(2,2);
	//	int idx=1;

	//	c1->cd(idx++);
	hHoriProj -> Draw();

	FitHole(hHoriProj);
}

void ProjectSieveHoleV_LHRSCenterAngle(void) {

	TFile
			*_filetmp =
					TFile::Open(
							"Production5/db_L.vdc.dat.Prod5.1.2GeV.sieve.root");
	assert(_filetmp);

	TCanvas * result = _filetmp->GetObjectChecked("SieveCheck",	"TCanvas");
	assert(result);

	result->Draw();

	TPad * pad = (TPad *) (result->GetListOfPrimitives()->FindObject("SieveCheck_4"));
	assert(pad);

	TH2D * Sieve_Foil0 = (TH2D *) (pad->GetListOfPrimitives()->FindObject(
			"Sieve_Foil3"));

	assert(Sieve_Foil0);
	Sieve_Foil0 -> Draw("COLZ");

	TH2D * hHori = new TH2D(*Sieve_Foil0);
	//	hHori -> Draw("COLZ");

	const Double_t vlim1 = -0.006, vlim2 = 0.006;// center hole
//		const Double_t vlim1 = -0.017, vlim2 = -0.01;// center hole to left


	for (int i = 1; i <= hHori->GetNbinsX(); i++) {
		for (int j = 1; j <= hHori->GetNbinsY(); j++) {
			const double binc = hHori->GetXaxis()->GetBinCenter(i);

			if (binc < vlim1 || binc > vlim2) {
				//				cout << i << binc << endl;
				hHori -> SetBinContent(i, j, 0);
			}
		}
	}

	//	hHori -> Draw("COLZ");
	TH1D* hHoriProj = hHori -> ProjectionY("hHoriProj");
	//

	c1 = new TCanvas("ProjectSieveHole", "ProjectSieveHole", 800, 500);
	//	c1 -> Divid(2,2);
	//	int idx=1;

	//	c1->cd(idx++);
	hHoriProj -> Draw();

	FitHole(hHoriProj);
}

void FitHole(TH1D * hHoriProj) {
	TVirtualFitter::SetDefaultFitter("Minuit2");


	RooRealVar xCoord("xCoord", "Phi", 16, 30);
	xCoord.setUnit("mrad");

	hHoriProj->GetXaxis()->Set(
			hHoriProj->GetXaxis()->GetNbins(),
			hHoriProj->GetXaxis()->GetXmin()*1000,
			hHoriProj->GetXaxis()->GetXmax()*1000
	);
	RooDataHist data("data","dataset with x",xCoord,hHoriProj) ;

	RooPlot* xframe = xCoord.frame() ;
	data.plotOn(xframe) ;

//	RooRealVar R("R", "R", 3.9878/2); // Big Hole
//	RooRealVar R("R", "R", 2/2); // Small Hole
	RooRealVar meang("meang","Center of Sieve Hole",22.1062) ;

	RooGenericPdf aPdf("aPdf", "aPdf",
			"(abs(xCoord- meang)<R)*(sqrt(R**2-(xCoord- meang)**2))", RooArgList(xCoord, R, meang));
//	aPdf.plotOn(xframe) ;

	RooRealVar meang0("meang0","mean of Gaussian",0) ;
	RooRealVar sigmag("sigmag","sigma of Gaussian",2.18292 ) ;
	RooGaussian gauss("gauss","gauss",xCoord,meang0,sigmag) ;
//	gauss.plotOn(xframe) ;

//	aPdf.setConvolutionWindow(meang,sigmag,5);

	RooNumConvPdf model("model","model",xCoord,gauss,aPdf) ;



	R.setConstant(kTRUE);
//	R.setConstant(kFALSE);
	meang.setConstant(kFALSE);
	meang0.setConstant(kTRUE);
	sigmag.setConstant(kFALSE);

	model.fitTo(data,RooFit::PrintLevel(3));
	R.Print();
	meang.Print();
	sigmag.Print();


	model.plotOn(xframe) ;
//
	aPdf.plotOn(xframe) ;

	new TCanvas("HoleFit.1PassCenterHole","");
	xframe->Draw();


}

void ProjectSieveHole(void) {
//	ProjectSieveHoleH_LHRSCenterAngle();
	ProjectSieveHoleV_LHRSCenterAngle();
	//	ProjectSieveHoleRHRS2Pass();
}