///////////////////////////////////////////////////////////////////////////
// A positive X value is to the beam left, a positive Y is too high.
// Spectrometer Slit: (DX/DY from spectrometer centerline; DZ from target)
// 		DX DY DZ
// 		-0.34 4.49 1157
Double_t SieveOffX = -(4.49)*1e-3;
Double_t SieveOffY = -(-0.34)*1e-3;
Double_t ZPos = 1157*1e-3;
Double_t SieveSpaceX = 1.025 * 25.400e-3;
Double_t SieveSpaceY = .513 * 25.400e-3;
Int_t NSieveSides = 3;

///////////////////////////////////////////////////////////////////////////

static const UInt_t NFiols = 7;
Double_t targetfoils[NFiols]={0.0019-0.2,0.0019-0.2+0.2/3.,0.0019-0.2+0.2/3.*2.,0.0019,0.0019+0.2/3.,0.0019+0.2/3.*2.,+0.0019+0.2};
Double_t VertexRes = 0.006;

///////////////////////////////////////////////////////////////////////////
#define DIPOLE_MAG2MOM(Mag) (2.702*Mag-1.6e-03*Mag*Mag*Mag)
///////////////////////////////////////////////////////////////////////////
void DrawSieve()
{
	Double_t MaxPlot = .1;
	for (Int_t x = -NSieveSides; x<=NSieveSides; x++)
	{
		TLine *l = new TLine(-MaxPlot,SieveSpaceX*x-SieveOffX,+MaxPlot,SieveSpaceX*x-SieveOffX);
		l->SetLineColor(3);
		l->Draw();
	}
	
	for (Int_t y = -NSieveSides; y<=NSieveSides; y++)
	{
		TLine *l = new TLine(SieveSpaceY*y-SieveOffY,-MaxPlot,SieveSpaceY*y-SieveOffY,MaxPlot);
		l->SetLineColor(3);
		l->Draw();
	}

	return;
}

DrawVertex()
{
    Double_t MaxPlot = 10000000000;
    for (UInt_t i; i<NFiols; i++)
    {
        TLine *l = new TLine(targetfoils[i],0,targetfoils[i],MaxPlot);
        l->SetLineColor(4);
        l->Draw();
    }
}

void FitVertex(TH1F * h)
{
	assert(h);

	const UInt_t VertexResN = 6;

	Double_t fArbitaryVertexShift[100]={0};


//	fArbitaryVertexShift[0] = 7.829146e-04;
//	fArbitaryVertexShift[1] = 7.561359e-05;
//	fArbitaryVertexShift[2] = -3.375388e-04;
//	fArbitaryVertexShift[3] = -2.566912e-04;
//	fArbitaryVertexShift[4] = -3.451242e-04;
//	fArbitaryVertexShift[5] = -7.884301e-04;
//	fArbitaryVertexShift[6] = -1.294849e-03;


	for(UInt_t FoilID=0; FoilID<NFiols; FoilID++)
// 		for(UInt_t FoilID=2; FoilID<3; FoilID++)
	{
		
		TString FitFunc = Form("Foild%d",FoilID);
		TF1 *f = new TF1(FitFunc,"gaus+[3]+[4]*x+[5]*x*x"
				,targetfoils[FoilID]-VertexRes*VertexResN,targetfoils[FoilID]+VertexRes*VertexResN);
		f->SetParameter(1,targetfoils[FoilID]);
		f->SetParameter(2,VertexRes);
		
// 		new TCanvas();
		TH1F * h1 = new TH1F(*h);
		for(UInt_t i=0; i<h1->GetNbinsX(); i++ )
		{
			if (TMath::Abs(h1->GetBinCenter(i)-targetfoils[FoilID])>=VertexRes*VertexResN)
			{
// 				cout<<"Removing Bin #"<<i<<endl;
				h1->SetBinContent(i,0);
			}
		}
// 		h1 -> Draw();
		h1 -> Fit(FitFunc,"N0");
		delete h1;
		f->Draw("SAME");
		TText *t = new TText(targetfoils[FoilID]
				,(f->GetParameter(0)
				+f->GetParameter(3)
						+f->GetParameter(4)*targetfoils[FoilID]
						+f->GetParameter(5)*targetfoils[FoilID]*targetfoils[FoilID])*1.02
				,Form("%2.1f +- %2.1f mm"
						,1000*(f->GetParameter(1)-targetfoils[FoilID])
						,1000*f->GetParameter(2)));
		t->SetTextSize(0.03);
		t->SetTextAlign(21);
		t->SetTextColor(2);
		t->Draw();

		fArbitaryVertexShift[FoilID] += f->GetParameter(1)-targetfoils[FoilID];
	}

	for(UInt_t FoilID=0; FoilID<NFiols; FoilID++)
		printf("opt->fArbitaryVertexShift[%d] = %e;\n"
				,FoilID,fArbitaryVertexShift[FoilID]);
	for(UInt_t FoilID=0; FoilID<NFiols; FoilID++)
		printf("fArbitaryVertexShift[%d] = %e;\n"
				,FoilID,fArbitaryVertexShift[FoilID]);
	
}


VertexZ4452()
{
	TCanvas * c1 = new TCanvas("VertexZ4452");
	
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_4452.root");
	T->Draw("ReactPt_L.z>>(800,-0.4,0.4)","L.tr.n==1&& abs(ReactPt_L.z)<.4 && abs(L.tr.tg_dp)<.05");
	DrawVertex();
}

VertexZ5832()
{
	TCanvas * c1 = new TCanvas("VertexZ5832");
	
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_5832.root");
	T->Draw("ReactPt_L.z>>CarbonFoil5832(500,-0.4,0.4)","abs(L.gold.dp)<0.05&&L.tr.n==1&&abs(L.tr.tg_th)<.07&&abs(L.tr.tg_ph)<.04");
	
	DrawVertex();
	FitVertex(CarbonFoil5832);

}

VertexZ2531()
{
	TCanvas * c1 = new TCanvas("VertexZ2531");
	
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_optics_L_2531.root");
	T->Draw("ReactPt_Lurb.z>>VertexZ2531(500,-0.4,0.4)","abs(L.gold.dp)<0.05&&L.tr.n==1&&abs(L.tr.tg_th)<.07&&abs(L.tr.tg_ph)<.04");
	TH1F * VertexZ2531 = gDirectory->GetObjectChecked("VertexZ2531","TH1F");
	DrawVertex();
	FitVertex(VertexZ2531);

}

VertexZ25052()
{
	TCanvas * c1 = new TCanvas("VertexZ25052");
	
// 	TFile *_file0 = TFile::Open("ROOTfiles/e06010_optics_L_2505.root");
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2505.root");
	T->Draw("ReactPt_Lurb.z>>VertexZ25052(500,-0.4,0.4)","abs(L.gold.dp)<0.05&&L.tr.n==1&&abs(L.tr.tg_th)<.07&&abs(L.tr.tg_ph)<.04");
	
	TH1F * VertexZ25052 = gDirectory->GetObjectChecked("VertexZ25052","TH1F");
	DrawVertex();
	FitVertex(VertexZ25052);

}

VertexZ(UInt_t Run)
{
	cout<<"VertexZ()"<<endl;
	TCanvas * c1 = new TCanvas(Form("VertexZ%d",Run),Form("VertexZ%d",Run));

//	TFile *_file0 = TFile::Open(Form("ROOTfiles/e06010_det_L_%d.root",Run));
	TFile *_file0 = TFile::Open(Form("ROOTfiles/e06010_L_%d.root",Run));
	T->Draw(Form("ReactPt_Lurb.z>>CarbonFoil%d(500,-0.4,0.4)",Run),"abs(L.gold.dp)<0.05&&L.tr.n==1&&abs(L.tr.tg_th)<.07&&abs(L.tr.tg_ph)<.04");
	
	DrawVertex();
	TH1F* h= (TH1F*)gDirectory->FindObjectAny(Form("CarbonFoil%d",Run));
	assert(h);
	
	FitVertex(h);

}


VertexZ2533()
{
	TCanvas * c1 = new TCanvas("VertexZ2533");
	
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2533.root");
	T->Draw("ReactPt_L.z>>CarbonFoil2533(800,-0.4,0.4)","L.tr.n==1&& abs(ReactPt_L.z)<.4 && abs(L.tr.tg_dp)<.05");
	DrawVertex();
}

SieveByFoil(TString RootFileID = "1.T400")
{
	TCanvas * c1 = new TCanvas("SieveCheck","SieveCheck",1900,1000);
	c1->Divide(3,3);
	idx=1;
	
	TFile *_file0 = TFile::Open(Form("ROOTfiles/e06010_det_L_2533.%s.root",RootFileID.Data()));

	c1->cd(idx++);
	T->Draw("-1.157*L.tr.tg_th:-L.tr.tg_y-1.157*L.tr.tg_ph>>CentralFoilMin3(100,-.1,.1,100,-.2,.2)","L.tr.n==1&& abs(ReactPt_L.z+.19)<.01 && abs(L.tr.tg_dp)<.05");
	DrawSieve();
	
	c1->cd(idx++);
	T->Draw("-1.157*L.tr.tg_th:-L.tr.tg_y-1.157*L.tr.tg_ph>>CentralFoilMin2(100,-.05,.05,100,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z+0.12)<.01 && abs(L.tr.tg_dp)<.05");
	DrawSieve();
	
	c1->cd(idx++);
	T->Draw("-1.157*L.tr.tg_th:-L.tr.tg_y-1.157*L.tr.tg_ph>>CentralFoilMin1(100,-.05,.05,100,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z+0.06)<.01 && abs(L.tr.tg_dp)<.05");
	DrawSieve();
	
	c1->cd(idx++);
	T->Draw("-1.157*L.tr.tg_th:-L.tr.tg_y-1.157*L.tr.tg_ph>>CentralFoil(100,-.05,.05,100,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z)<.015 &&abs(L.tr.tg_dp)<.05");
	DrawSieve();
	
	c1->cd(idx++);
	T->Draw("-1.157*L.tr.tg_th:-L.tr.tg_y-1.157*L.tr.tg_ph>>CentralFoilPlus1(100,-.05,.05,100,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z-0.075)<.015 && abs(L.tr.tg_dp)<.05");
	DrawSieve();
	CentralFoilPlus1->SetTitle("CentralFoilPlus1");
	
	c1->cd(idx++);
	T->Draw("-1.157*L.tr.tg_th:-L.tr.tg_y-1.157*L.tr.tg_ph>>CentralFoilPlus2(100,-.05,.05,100,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z-0.135)<.015 && abs(L.tr.tg_dp)<.05");
	DrawSieve();
	
	c1->cd(idx++);
	T->Draw("-1.157*L.tr.tg_th:-L.tr.tg_y-1.157*L.tr.tg_ph>>CentralFoilPlus3(100,-.1,.1,100,-.2,.2)","L.tr.n==1&& abs(ReactPt_L.z-0.24)<.03 && abs(L.tr.tg_dp)<.05");
	DrawSieve();

	
	CentralFoilMin3->SetTitle(RootFileID+".CentralFoilMin3");
	CentralFoilMin2->SetTitle(RootFileID+".CentralFoilMin2");
	CentralFoilMin1->SetTitle(RootFileID+".CentralFoilMin1");
	CentralFoil->SetTitle(RootFileID+".CentralFoil");
	CentralFoilPlus1->SetTitle(RootFileID+".CentralFoilPlus1");
	CentralFoilPlus2->SetTitle(RootFileID+".CentralFoilPlus2");
	CentralFoilPlus3->SetTitle(RootFileID+".CentralFoilPlus3");

	c1->Print(Form("ROOTfiles/e06010_det_L_2533.%s.root.png",RootFileID.Data()),"png");
}

const Double_t ThetaAccept = 0.062;
const Double_t PhiAccept = 0.04;

OpticsAcceptanceCorrelationCheck_MakeCut()
{
	// build R Cuts
	
	TChain * T=new TChain("T");
// 	RootFileLoader(T,"./ROOTfiles/e06010_prod4_%d.root",4100,4200);
// 	RootFileLoader(T,"./ROOTfiles/e06010_prod4_%d.root",4242,4300);
	RootFileLoader(T,"/home/jinhuang/e06010/p4_skim//e06010_phys_%d.root",4240,4240);
	
// 	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.23/3 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035/10 && L.tr.n==1 && abs(ReactPt_L.z)<.02");
// 	T->Draw(">>LooseCut","abs(ReactPt_L.z)<.3 && abs(L.tr.tg_th)<.1 && abs(L.tr.tg_ph)<.06 && abs(L.tr.tg_dp)<.075 && L.tr.n==1");
// 	T->Draw("1","abs(ReactPt_L.z)<.185 && abs(ExTgtCor_L.th)<.1 && abs(ExTgtCor_L.ph)<.06 && abs(ExTgtCor_L.dp)<.075 && L.tr.n==1");
	T->Draw(">>LooseCut","abs(ReactPt_L.z)<.185 && abs(ExTgtCor_L.th)<.1 && abs(ExTgtCor_L.ph)<.06 && abs(ExTgtCor_L.dp)<.075 && L.tr.n==1");
	
	TString CutUsed = "LooseCut";
	TString DrawOption  = "COLZ";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);
	
	TCanvas *c1 = new TCanvas("OpticsAcceptanceCorrelationCheck","OpticsAcceptanceCorrelationCheck",1900,1000);

	c1->Divide(3,2);
	UInt_t idx=1;
	
	TPad *cpad;
	TCutG* cutg;
	fstream flog("OpticsAcceptanceCorrelationCheck.C",ios_base::out|ios_base::app);
	assert(flog.is_open());

	c1->cd(idx++);
	T->Draw("ExTgtCor_L.dp:L.tr.tg_y>>dpy(400,-.1,.1,400,-.08,.08)","",DrawOption);

	cpad = gPad;
	cutg = (TCutG*) (cpad->WaitPrimitive("CUTG","CutG")); // making cut, store to CUTG
	assert(cutg);
	c1->Update();
	cutg->SetName(Form("gcut_dpy",i)); //
  // set axises' name
	cutg->SetVarX("L.tr.tg_y");
	cutg->SetVarY("ExTgtCor_L.dp");
	cutg->SavePrimitive(flog);
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");
	
	c1->cd(idx++);
	T->Draw("ExTgtCor_L.dp:ExTgtCor_L.th>>dpth(400,-.1,.1,400,-.08,.08)","",DrawOption);

	cpad = gPad;
	cutg = (TCutG*) (cpad->WaitPrimitive("CUTG","CutG")); // making cut, store to CUTG
	assert(cutg);
	c1->Update();
	cutg->SetName(Form("gcut_dpth",i)); //
  // set axises' name
	cutg->SetVarX("ExTgtCor_L.th");
	cutg->SetVarY("ExTgtCor_L.dp");
	cutg->SavePrimitive(flog);
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");
	
	c1->cd(idx++);
	T->Draw("ExTgtCor_L.dp:ExTgtCor_L.ph>>dpph(400,-.06,.06,400,-.08,.08)","",DrawOption);
	
	cpad = gPad;
	cutg = (TCutG*) (cpad->WaitPrimitive("CUTG","CutG")); // making cut, store to CUTG
	assert(cutg);
	c1->Update();
	cutg->SetName(Form("gcut_dpph",i)); //
  // set axises' name
	cutg->SetVarX("ExTgtCor_L.ph");
	cutg->SetVarY("ExTgtCor_L.dp");
	cutg->SavePrimitive(flog);
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");
	
	c1->cd(idx++);
	T->Draw("L.tr.tg_y:ExTgtCor_L.th>>yth(400,-.1,.1,400,-.1,.1)","",DrawOption);
	
	cpad = gPad;
	cutg = (TCutG*) (cpad->WaitPrimitive("CUTG","CutG")); // making cut, store to CUTG
	assert(cutg);
	c1->Update();
	cutg->SetName(Form("gcut_yth",i)); //
  // set axises' name
	cutg->SetVarX("ExTgtCor_L.th");
	cutg->SetVarY("L.tr.tg_y");
	cutg->SavePrimitive(flog);
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");
	
	c1->cd(idx++);
	T->Draw("L.tr.tg_y:ExTgtCor_L.ph>>yph(400,-.06,.06,400,-.1,.1)","",DrawOption);
	
	cpad = gPad;
	cutg = (TCutG*) (cpad->WaitPrimitive("CUTG","CutG")); // making cut, store to CUTG
	assert(cutg);
	c1->Update();
	cutg->SetName(Form("gcut_yph",i)); //
  // set axises' name
	cutg->SetVarX("ExTgtCor_L.ph");
	cutg->SetVarY("L.tr.tg_y");
	cutg->SavePrimitive(flog);
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");
	
	c1->cd(idx++);
	T->Draw("ExTgtCor_L.th:ExTgtCor_L.ph>>thph(400,-.06,.06,400,-.1,.1)","",DrawOption);
	
	cpad = gPad;
	cutg = (TCutG*) (cpad->WaitPrimitive("CUTG","CutG")); // making cut, store to CUTG
	assert(cutg);
	c1->Update();
	cutg->SetName(Form("gcut_thph",i)); //
  // set axises' name
	cutg->SetVarX("ExTgtCor_L.ph");
	cutg->SetVarY("ExTgtCor_L.th");
	cutg->SavePrimitive(flog);
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");

	flog.close();

	SaveCanvas(c1,TString(c1->GetName())+".42xx",kFALSE);
	
}

OpticsAcceptanceCorrelationCheck_DrawCut()
{
	// build R Cuts
	
	TChain * T=new TChain("T");
// 	RootFileLoader(T,"./ROOTfiles/e06010_prod4_%d.root",4100,4200);
// 	RootFileLoader(T,"./ROOTfiles/e06010_prod4_%d.root",4242,4300);
	RootFileLoader(T,"/home/jinhuang/e06010/p4_skim//e06010_phys_%d.root",4240,4299);
	
// 	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.23/3 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035/10 && L.tr.n==1 && abs(ReactPt_L.z)<.02");
// 	T->Draw(">>LooseCut","abs(ReactPt_L.z)<.3 && abs(L.tr.tg_th)<.1 && abs(L.tr.tg_ph)<.06 && abs(L.tr.tg_dp)<.075 && L.tr.n==1");
// 	T->Draw("1","abs(ReactPt_L.z)<.185 && abs(ExTgtCor_L.th)<.1 && abs(ExTgtCor_L.ph)<.06 && abs(ExTgtCor_L.dp)<.075 && L.tr.n==1");
	T->Draw(">>LooseCut","abs(ReactPt_L.z)<.185 && abs(ExTgtCor_L.th)<.1 && abs(ExTgtCor_L.ph)<.06 && abs(ExTgtCor_L.dp)<.075 && L.tr.n==1");
	
	TString CutUsed = "LooseCut";
	TString DrawOption  = "COLZ";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);
	
	TCanvas *c1 = new TCanvas("OpticsAcceptanceCorrelationCheck_DrawCut","OpticsAcceptanceCorrelationCheck_DrawCut",1900,1000);

	c1->Divide(3,2);
	UInt_t idx=1;
	
	gInterpreter->ExecuteMacro("OpticsAcceptanceCorrelationCheck.C");

	c1->cd(idx++);c1->Update();
	T->Draw("ExTgtCor_L.dp:L.tr.tg_y>>dpy(400,-.1,.1,400,-.08,.08)","",DrawOption);

	cutg = gcut_dpy; // making cut, store to CUTG
	assert(cutg);
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");
	
	c1->cd(idx++);c1->Update();
	T->Draw("ExTgtCor_L.dp:ExTgtCor_L.th>>dpth(400,-.1,.1,400,-.08,.08)","",DrawOption);

	cutg= gcut_dpth; //
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");

	c1->cd(idx++);c1->Update();
	T->Draw("ExTgtCor_L.dp:ExTgtCor_L.ph>>dpph(400,-.06,.06,400,-.08,.08)","",DrawOption);

	cutg=gcut_dpph; //
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");

	c1->cd(idx++);c1->Update();
	T->Draw("L.tr.tg_y:ExTgtCor_L.th>>yth(400,-.1,.1,400,-.1,.1)","",DrawOption);

	cutg=gcut_yth; //
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");

	c1->cd(idx++);c1->Update();
	T->Draw("L.tr.tg_y:ExTgtCor_L.ph>>yph(400,-.06,.06,400,-.1,.1)","",DrawOption);

	cutg=gcut_yph; //
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");

	c1->cd(idx++);c1->Update();
	T->Draw("ExTgtCor_L.th:ExTgtCor_L.ph>>thph(400,-.06,.06,400,-.1,.1)","",DrawOption);

	cutg=gcut_thph; //
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");


	SaveCanvas(c1,TString(c1->GetName())+".42xx",kFALSE);
	
}

void OpticsAcceptanceRCutCheck_VZCompare()
{
	//compare positive and negative he3 VZ after R-Cuts
	TFile *_fileNeg = TFile::Open("./ROOTfiles/e06010_prod4_4111.root");
	TFile *_filePos = TFile::Open("./ROOTfiles/e06010_prod4_4251.root");

	TCanvas *c1 = new TCanvas("OpticsAcceptanceRCutCheck_VZCompare","OpticsAcceptanceRCutCheck_VZCompare",1400,1100);
	gROOT->ProcessLine(".x OpticsAcceptanceCorrelationCheck.C");
	c1->Divide(2,2);
	UInt_t idx=1;

	const Double_t scaling = 1350./960.;
	
	c1->cd(idx++);
	_fileNeg->cd();
	T->Draw("ReactPt_L.z>>zNeg(100,-.25,.25)","L.tr.n==1","GOFF");
	T->Draw("ReactPt_L.z>>zNeg_cut(100,-.25,.25)",HRS_RCut&&"L.tr.n==1","GOFF");
	zNeg->Scale(scaling);
	zNeg_cut->Scale(scaling);
	zNeg->SetLineColor(1);
	zNeg_cut->SetLineColor(4);
// 	zNeg_cut->SetFillColor(4);
// 	zNeg_cut->SetFillStyle(3008);
	zNeg->SetTitle("HRS Negative Mode, Before and After R-Cuts : React Z /m");
	zNeg->Draw();
	zNeg_cut->Draw("SAME");
	
	c1->cd(idx++);
	_filePos->cd();
	T->Draw("ReactPt_L.z>>zPos(100,-.25,.25)","L.tr.n==1","GOFF");
	T->Draw("ReactPt_L.z>>zPos_cut(100,-.25,.25)",HRS_RCut&&"L.tr.n==1","GOFF");
	zPos->SetLineColor(6);
	zPos_cut->SetLineColor(2);
// 	zPos_cut->SetFillColor(2);
// 	zPos_cut->SetFillStyle(3008);
	zPos->SetTitle("HRS Positive Mode, Before and After R-Cuts : React Z /m");
	zPos->Draw();
	zPos_cut->Draw("SAME");
	
	c1->cd(idx++);
	zPos->Clone("zPos2");
	zPos2->SetTitle("Compare Before R-Cuts : React Z /m");
	zPos2->Draw();
	zNeg->Draw("SAME");
	
	c1->cd(idx++);
	zPos_cut->Clone("zPos_cut2");
	zPos_cut2->SetTitle("Compare After R-Cuts : React Z /m");
	zPos_cut2->Draw();
	zNeg_cut->Draw("SAME");
}

void OpticsAcceptanceRCutCheck()
{
	

	TChain * T=new TChain("T");
// 	RootFileLoader(T,"./ROOTfiles/e06010_prod4_%d.root",4100,4200);
// 	RootFileLoader(T,"./ROOTfiles/e06010_prod4_%d.root",4111,4111);
// 	RootFileLoader(T,"./ROOTfiles/e06010_prod4_%d.root",4251,4251);
	RootFileLoader(T,"./ROOTfiles/e06010_prod4_%d.root",4240,4240);

	T->Draw(">>LooseCut","L.tr.n==1");
	
	TString CutUsed = "LooseCut";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);
	
	TCanvas *c1 = new TCanvas("OpticsAcceptanceRCutCheck","OpticsAcceptanceRCutCheck",1900,1000);

	gROOT->ProcessLine(".x OpticsAcceptanceCorrelationCheck.C");
	
	c1->Divide(3,2);
	UInt_t idx=1;
	
	c1->cd(idx++);
	T->Draw("L.tr.tg_y>>tg_y(200,-.2,.2)","L.tr.n==1","GOFF");
	T->Draw("L.tr.tg_y>>tg_y_cut(200,-.2,.2)",HRS_RCut,"GOFF");
	tg_y_cut->SetLineColor(4);
	tg_y_cut->SetFillColor(4);
	tg_y_cut->SetFillStyle(3008);
	tg_y->SetTitle("Target Y /m");
	tg_y->Draw();
	tg_y_cut->Draw("SAME");
	
	c1->cd(idx++);
	T->Draw("ReactPt_L.z>>z(200,-.25,.25)","L.tr.n==1","GOFF");
	T->Draw("ReactPt_L.z>>z_cut(200,-.25,.25)",HRS_RCut,"GOFF");
	z_cut->SetLineColor(4);
	z_cut->SetFillColor(4);
	z_cut->SetFillStyle(3008);
	z->SetTitle("React Z /m");
	z->Draw();
	z_cut->Draw("SAME");
	
	c1->cd(idx++);
	T->Draw("ExTgtCor_L.th>>th(200,-.1,.1)","L.tr.n==1","GOFF");
	T->Draw("ExTgtCor_L.th>>th_cut(200,-.1,.1)",HRS_RCut,"GOFF");
	th_cut->SetLineColor(4);
	th_cut->SetFillColor(4);
	th_cut->SetFillStyle(3008);
	th->SetTitle("Target Theta /rad");
	th->Draw();
	th_cut->Draw("SAME");
	
	c1->cd(idx++);
	T->Draw("ExTgtCor_L.ph>>ph(200,-.07,.07)","L.tr.n==1","GOFF");
	T->Draw("ExTgtCor_L.ph>>ph_cut(200,-.07,.07)",HRS_RCut,"GOFF");
	ph_cut->SetLineColor(4);
	ph_cut->SetFillColor(4);
	ph_cut->SetFillStyle(3008);
	ph->SetTitle("Target Phi /rad");
	ph->Draw();
	ph_cut->Draw("SAME");
	
	c1->cd(idx++);
	T->Draw("ExTgtCor_L.dp>>dp(200,-.08,.08)","L.tr.n==1","GOFF");
	T->Draw("ExTgtCor_L.dp>>dp_cut(200,-.08,.08)",HRS_RCut,"GOFF");
	dp_cut->SetLineColor(4);
	dp_cut->SetFillColor(4);
	dp_cut->SetFillStyle(3008);
	dp->SetTitle("Target Delta P");
	dp->Draw();
	dp_cut->Draw("SAME");
	
	c1->cd(idx++);
	T->Draw("ExTgtCor_L.x>>x(200,-.04,.04)","L.tr.n==1","GOFF");
	T->Draw("ExTgtCor_L.x>>x_cut(200,-.04,.04)",HRS_RCut,"GOFF");
	x_cut->SetLineColor(4);
	x_cut->SetFillColor(4);
	x_cut->SetFillStyle(3008);
	x->SetTitle("Target X/m");
	x->Draw();
	x_cut->Draw("SAME");
	
	
}

OpticsAcceptanceCorrelationCheck_Old()
{
	
	TChain * T=new TChain("T");
	UInt_t startrun = 4100, endrun =4300;
	RootFileLoader(T,"./ROOTfiles/e06010_prod4_%d.root",startrun,endrun);
	
// 	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.23/3 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035/10 && L.tr.n==1 && abs(ReactPt_L.z)<.02");
// 	T->Draw(">>LooseCut","abs(ReactPt_L.z)<.3 && abs(L.tr.tg_th)<.1 && abs(L.tr.tg_ph)<.06 && abs(L.tr.tg_dp)<.075 && L.tr.n==1");
	T->Draw(">>LooseCut","abs(ReactPt_L.z)<.185 && abs(L.tr.tg_th)<.1 && abs(L.tr.tg_ph)<.06 && abs(L.tr.tg_dp)<.075 && L.tr.n==1");
	
	TString CutUsed = "LooseCut";
	TString DrawOption  = "COLZ";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);
	
	TCanvas *c1 = new TCanvas("OpticsAcceptanceCorrelationCheck","OpticsAcceptanceCorrelationCheck",1900,1000);

	c1->Divide(3,2);
	UInt_t idx=1;
	
	TPad *cpad;
	TCutG* cutg;
	fstream flog("OpticsAcceptanceCorrelationCheck.C",ios_base::out|ios_base::app);
	assert(flog.is_open());

	c1->cd(idx++);
	T->Draw("L.tr.tg_dp:L.tr.tg_y>>dpy(400,-.1,.1,400,-.08,.08)","",DrawOption);

	cpad = gPad;
	cutg = (TCutG*) (cpad->WaitPrimitive("CUTG","CutG")); // making cut, store to CUTG
	assert(cutg);
	c1->Update();
	cutg->SetName(Form("gcut_dpy",i)); //
  // set axises' name
	cutg->SetVarX("L.tr.tg_y");
	cutg->SetVarY("L.tr.tg_dp");
	cutg->SavePrimitive(flog);
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");
	
	c1->cd(idx++);
	T->Draw("L.tr.tg_dp:L.tr.tg_th>>dpth(400,-.1,.1,400,-.08,.08)","",DrawOption);

	cpad = gPad;
	cutg = (TCutG*) (cpad->WaitPrimitive("CUTG","CutG")); // making cut, store to CUTG
	assert(cutg);
	c1->Update();
	cutg->SetName(Form("gcut_dpth",i)); //
  // set axises' name
	cutg->SetVarX("L.tr.tg_th");
	cutg->SetVarY("L.tr.tg_dp");
	cutg->SavePrimitive(flog);
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");
	
	c1->cd(idx++);
	T->Draw("L.tr.tg_dp:L.tr.tg_ph>>dpph(400,-.06,.06,400,-.08,.08)","",DrawOption);
	
	cpad = gPad;
	cutg = (TCutG*) (cpad->WaitPrimitive("CUTG","CutG")); // making cut, store to CUTG
	assert(cutg);
	c1->Update();
	cutg->SetName(Form("gcut_dpph",i)); //
  // set axises' name
	cutg->SetVarX("L.tr.tg_ph");
	cutg->SetVarY("L.tr.tg_dp");
	cutg->SavePrimitive(flog);
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");
	
	c1->cd(idx++);
	T->Draw("L.tr.tg_y:L.tr.tg_th>>yth(400,-.1,.1,400,-.1,.1)","",DrawOption);
	
	cpad = gPad;
	cutg = (TCutG*) (cpad->WaitPrimitive("CUTG","CutG")); // making cut, store to CUTG
	assert(cutg);
	c1->Update();
	cutg->SetName(Form("gcut_yth",i)); //
  // set axises' name
	cutg->SetVarX("L.tr.tg_th");
	cutg->SetVarY("L.tr.tg_y");
	cutg->SavePrimitive(flog);
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");
	
	c1->cd(idx++);
	T->Draw("L.tr.tg_y:L.tr.tg_ph>>yph(400,-.06,.06,400,-.1,.1)","",DrawOption);
	
	cpad = gPad;
	cutg = (TCutG*) (cpad->WaitPrimitive("CUTG","CutG")); // making cut, store to CUTG
	assert(cutg);
	c1->Update();
	cutg->SetName(Form("gcut_yph",i)); //
  // set axises' name
	cutg->SetVarX("L.tr.tg_ph");
	cutg->SetVarY("L.tr.tg_y");
	cutg->SavePrimitive(flog);
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");
	
	c1->cd(idx++);
	T->Draw("L.tr.tg_th:L.tr.tg_ph>>thph(400,-.06,.06,400,-.1,.1)","",DrawOption);
	
	cpad = gPad;
	cutg = (TCutG*) (cpad->WaitPrimitive("CUTG","CutG")); // making cut, store to CUTG
	assert(cutg);
	c1->Update();
	cutg->SetName(Form("gcut_thph",i)); //
  // set axises' name
	cutg->SetVarX("L.tr.tg_ph");
	cutg->SetVarY("L.tr.tg_th");
	cutg->SavePrimitive(flog);
	cutg->SetLineColor(6);
	cutg->SetLineWidth(2);
	cutg->Draw("L");

	flog.close();
}

DrawAngularAcceptanceLine()
{
	const Double_t MaxPlot = 10000;
	
	TLine *l = new TLine(-MaxPlot,-ThetaAccept,+MaxPlot,-ThetaAccept);
	l->SetLineColor(3);
	l->Draw();
	TLine *l = new TLine(-MaxPlot,ThetaAccept,+MaxPlot,ThetaAccept);
	l->SetLineColor(3);
	l->Draw();
	
	TLine *l = new TLine(-PhiAccept,-MaxPlot,-PhiAccept,+MaxPlot);
	l->SetLineColor(3);
	l->Draw();
	TLine *l = new TLine(PhiAccept,-MaxPlot,PhiAccept,+MaxPlot);
	l->SetLineColor(3);
	l->Draw();
}

AngularAcceptanceCheckByFoil(TString RootFileID = "Prod2")
{
	
	gStyle->SetOptStat(11);
	
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2533.root");

// 	c1->cd(idx++);
	T->Draw("L.tr.tg_th:L.tr.tg_ph>>CentralFoilMin3(500,-.05,.05,500,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z+.2)<.01 && abs(L.tr.tg_dp)<.05","COLZ");
	
// 	c1->cd(idx++);
	T->Draw("L.tr.tg_th:L.tr.tg_ph>>CentralFoilMin2(500,-.05,.05,500,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z+0.134)<.01 && abs(L.tr.tg_dp)<.05","COLZ");
	
// 	c1->cd(idx++);
	T->Draw("L.tr.tg_th:L.tr.tg_ph>>CentralFoilMin1(500,-.05,.05,500,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z+0.067)<.01 && abs(L.tr.tg_dp)<.05","COLZ");
	
// 	c1->cd(idx++);
	T->Draw("L.tr.tg_th:L.tr.tg_ph>>CentralFoil(500,-.05,.05,500,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z)<.015 &&abs(L.tr.tg_dp)<.05","COLZ");
	
// 	c1->cd(idx++);
	T->Draw("L.tr.tg_th:L.tr.tg_ph>>CentralFoilPlus1(500,-.05,.05,500,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z-0.063)<.015 && abs(L.tr.tg_dp)<.05","COLZ");
	
// 	c1->cd(idx++);
	T->Draw("L.tr.tg_th:L.tr.tg_ph>>CentralFoilPlus2(500,-.05,.05,500,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z-0.13)<.015 && abs(L.tr.tg_dp)<.05","COLZ");
	
	TString sTitle = "Th vs Ph (Target) at ";
	CentralFoilMin3->SetTitle(sTitle+"Central Foil Upstream 3");
	CentralFoilMin2->SetTitle(sTitle+"Central Foil Upstream 2");
	CentralFoilMin1->SetTitle(sTitle+"Central Foil Upstream 1");
	CentralFoil->SetTitle(sTitle+"Central Foil");
	CentralFoilPlus1->SetTitle(sTitle+"Central Foil Downstream 1");
	CentralFoilPlus2->SetTitle(sTitle+"Central Foil Downstream 2");

	TCanvas * c1 = new TCanvas("AngularCheck","AngularCheck",1900,1000);
	c1->Divide(3,2);
	UInt_t idx=1;
	
	c1->cd(idx++);
	CentralFoilMin3->Draw("COLZ");
	DrawAngularAcceptanceLine();
	c1->cd(idx++);
	CentralFoilMin2->Draw("COLZ");
	DrawAngularAcceptanceLine();
	c1->cd(idx++);
	CentralFoilMin1->Draw("COLZ");
	DrawAngularAcceptanceLine();
	c1->cd(idx++);
	CentralFoil->Draw("COLZ");
	DrawAngularAcceptanceLine();
	c1->cd(idx++);
	CentralFoilPlus1->Draw("COLZ");
	DrawAngularAcceptanceLine();
	c1->cd(idx++);
	CentralFoilPlus2->Draw("COLZ");
	DrawAngularAcceptanceLine();
	

	c1->Print(Form("AngularAcceptanceCheck.%s.e06010_det_L_2533.root.png",RootFileID.Data()),"png");
}

AngularAcceptanceCheck(TString RootFileID = "Prod2")
{
	
	gStyle->SetOptStat(11);
	
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2533.root");

// 	c1->cd(idx++);
	T->Draw("L.tr.tg_th:L.tr.tg_ph>>CentralFoilMin3(500,-.05,.05,500,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z+.2)<.01 && abs(L.tr.tg_dp)<.05","COLZ");
	
// 	c1->cd(idx++);
	T->Draw("L.tr.tg_th:L.tr.tg_ph>>CentralFoilMin2(500,-.05,.05,500,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z+0.134)<.01 && abs(L.tr.tg_dp)<.05","COLZ");
	
// 	c1->cd(idx++);
	T->Draw("L.tr.tg_th:L.tr.tg_ph>>CentralFoilMin1(500,-.05,.05,500,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z+0.067)<.01 && abs(L.tr.tg_dp)<.05","COLZ");
	
// 	c1->cd(idx++);
	T->Draw("L.tr.tg_th:L.tr.tg_ph>>CentralFoil(500,-.05,.05,500,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z)<.015 &&abs(L.tr.tg_dp)<.05","COLZ");
	
// 	c1->cd(idx++);
	T->Draw("L.tr.tg_th:L.tr.tg_ph>>CentralFoilPlus1(500,-.05,.05,500,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z-0.063)<.015 && abs(L.tr.tg_dp)<.05","COLZ");
	
// 	c1->cd(idx++);
	T->Draw("L.tr.tg_th:L.tr.tg_ph>>CentralFoilPlus2(500,-.05,.05,500,-.1,.1)","L.tr.n==1&& abs(ReactPt_L.z-0.13)<.015 && abs(L.tr.tg_dp)<.05","COLZ");

	const TH2D & htmp = *CentralFoilMin3;
	TH2D * SieveSum = new TH2D(htmp);
	SieveSum -> Add(CentralFoilMin2);
	SieveSum -> Add(CentralFoilMin1);
	SieveSum -> Add(CentralFoil);
	SieveSum -> Add(CentralFoilPlus1);
	SieveSum -> Add(CentralFoilPlus2);
	
	SieveSum->SetTitle("Superposing all Sieve Plots");
	SieveSum->SetXTitle("Target Phi (L.tr.ph)");
	SieveSum->SetYTitle("Target Theta (L.tr.th)");

	
	TChain * t=new TChain("T");
	RootFileLoader(t,"./ROOTfiles/e06010_prod3_%d.root",5816,5817);
	RootFileLoader(t,"./ROOTfiles/e06010_prod3_%d.root",5805,5815);
	
	TCut cut = "L.tr.n==1 && BB.tr.n ==1 && abs(CT.pi.t)<.5 && L.cer.asum_c<10&&L.a1.asum_c>200&&(L.prl1.e+L.prl2.e)/2350./(1+L.tr.tg_dp)>0.04&&(L.prl1.e+L.prl2.e)/2350./(1+L.tr.tg_dp)<0.12 && !DL.edtpl&&(DL.evtypebits&1<<5)==1<<5 &&  BB.tr.chi2[BB.ts.sh_flag_tr[]]/BB.tr.ndof[BB.ts.sh_flag_tr[]]<2.4 && BB.tr.n > 0 && BB.tr.p[BB.ts.sh_flag_tr[]]>0.6 && BB.tr.p[BB.ts.sh_flag_tr[]]<2.2  && (BB.optics.charge[BB.ts.sh_flag_tr[]] ==-1 ||BB.optics.charge[BB.ts.sh_flag_tr[]] ==1) && abs(BB.tr.vz[BB.ts.sh_flag_tr[]])<0.2 && BB.ts.ps_e[]>350. && abs((BB.ts.ps_e[]/2+BB.ts.sh_e[])/BB.tr.p[BB.ts.sh_flag_tr[]]-1000)<270 && L.tr.n==1&& abs(ReactPt_L.z)<.17 && abs(ReactPt_L.z-BB.tr.vz)<.5 && abs(L.tr.tg_dp)<.05";

	const UInt_t NPassCut  = t->GetEntries(cut);
	t->Draw("L.tr.tg_th:L.tr.tg_ph>>Prod5816(500,-.05,.05,500,-.1,.1)",cut);
	Prod5816->SetTitle("Production Run with (e'pi) cuts");
	Prod5816->SetXTitle("Target Phi (L.tr.ph)");
	Prod5816->SetYTitle("Target Theta (L.tr.th)");

	////////////////////////////////////////////////////////////////////////////////
	TCanvas * c1 = new TCanvas("AngularCheck1","AngularCheck1",1200,1000);
	c1->Divide(2,1);
	UInt_t idx=1;
	
	c1->cd(idx++);
	SieveSum->Draw("COLZ");
	DrawAngularAcceptanceLine();
	c1->cd(idx++);
	Prod5816->Draw("SCAT");
	DrawAngularAcceptanceLine();
	
	c1->Print(Form("AngularAcceptanceCheck.%s.Sum.png",RootFileID.Data()),"png");

	////////////////////////////////////////////////////////////////////////////////


	Prod5816->ProjectionX("Prod5816Ph");
	Prod5816->ProjectionY("Prod5816Th");
	
	Prod5816Ph->SetTitle("Target Phi Acceptance (e'pi) Production");
	Prod5816Ph->SetXTitle("Target Phi (L.tr.ph)");
	Prod5816Th->SetTitle("Target Theta Acceptance (e'pi) Production");
	Prod5816Th->SetXTitle("Target Theta (L.tr.th)");
	
	Double_t NPassThCut  = 0;
	Double_t NPassPhCut  = 0;

	TH1D* h = 0;
	h=Prod5816Ph;
	for(UInt_t i=0;i<h->GetNbinsX(); i++)
		if (TMath::Abs(h->GetBinCenter(i))<=PhiAccept) NPassPhCut+=h->GetBinContent(i);
	h=Prod5816Th;
	for(UInt_t i=0;i<h->GetNbinsX(); i++)
		if (TMath::Abs(h->GetBinCenter(i))<=ThetaAccept) NPassThCut+=h->GetBinContent(i);

	
	Prod5816Ph->Rebin(4);
	Prod5816Th->Rebin(4);
	
	TCanvas * c1 = new TCanvas("AngularCheck2","AngularCheck2",1200,1000);
	c1->Divide(2,1);
	UInt_t idx=1;

// 	const Double_t MaxPlot=10000;
	c1->cd(idx++);
	Prod5816Ph->Draw();
	TLine *l = new TLine(-PhiAccept,-0,-PhiAccept,+Prod5816Ph->GetMaximum());
	l->SetLineColor(4);
	l->Draw();
	TLine *l = new TLine(PhiAccept,-0,PhiAccept,+Prod5816Ph->GetMaximum());
	l->SetLineColor(4);
	l->Draw();

	TText *text = new TText(0,Prod5816Ph->GetMaximum()/2,Form("Cut Acceptance ~ %2.1f%%",100*NPassPhCut/NPassCut));
	text->SetTextSize(0.05);
	text->SetTextAlign(22);
	text->SetTextColor(4);
	text->Draw();
	
	c1->cd(idx++);
	Prod5816Th->Draw();
	TLine *l = new TLine(-ThetaAccept,-0,-ThetaAccept,+Prod5816Th->GetMaximum());
	l->SetLineColor(4);
	l->Draw();
	TLine *l = new TLine(ThetaAccept,0,ThetaAccept,+Prod5816Th->GetMaximum());
	l->SetLineColor(4);
	l->Draw();

	TText *text = new TText(0,Prod5816Th->GetMaximum()/2,Form("Cut Acceptance ~ %2.1f%%",100*NPassThCut/NPassCut));
	text->SetTextSize(0.05);
	text->SetTextAlign(22);
	text->SetTextColor(4);
	text->Draw();

	
	c1->Print(Form("AngularAcceptanceCheck.%s.ThPh.png",RootFileID.Data()),"png");
}

	const Double_t VZAcceptance = 0.185;
VertexAcceptanceCheck(TString RootFileID = "Prod2")
{
	
	gStyle->SetOptStat(11);
	
	
	TChain * t=new TChain("T");
	RootFileLoader(t,"./ROOTfiles/e06010_prod3_%d.root",5816,5817);
	RootFileLoader(t,"./ROOTfiles/e06010_prod3_%d.root",5805,5815);
	
	TCut cut = "L.tr.n==1 && BB.tr.n ==1 && abs(CT.pi.t)<.5 && L.cer.asum_c<10&&L.a1.asum_c>200&&(L.prl1.e+L.prl2.e)/2350./(1+L.tr.tg_dp)>0.04&&(L.prl1.e+L.prl2.e)/2350./(1+L.tr.tg_dp)<0.12 && !DL.edtpl&&(DL.evtypebits&1<<5)==1<<5 &&  BB.tr.chi2[BB.ts.sh_flag_tr[]]/BB.tr.ndof[BB.ts.sh_flag_tr[]]<2.4 && BB.tr.n > 0 && BB.tr.p[BB.ts.sh_flag_tr[]]>0.6 && BB.tr.p[BB.ts.sh_flag_tr[]]<2.2  && (BB.optics.charge[BB.ts.sh_flag_tr[]] ==-1 ||BB.optics.charge[BB.ts.sh_flag_tr[]] ==1) && abs(BB.tr.vz[BB.ts.sh_flag_tr[]])<0.2 && BB.ts.ps_e[]>350. && abs((BB.ts.ps_e[]/2+BB.ts.sh_e[])/BB.tr.p[BB.ts.sh_flag_tr[]]-1000)<270 && L.tr.n==1 && abs(ReactPt_L.z-BB.tr.vz)<.5 && abs(L.tr.tg_dp)<.05 && abs(L.tr.tg_th)<.06 && abs(L.tr.tg_ph)<.04";

	const UInt_t NPassCut  = t->GetEntries(cut);
	t->Draw("ReactPt_L.z>>Prod5816VZ(800,-.3,.3)",cut);
	Prod5816VZ->SetTitle("LHRS React with (e'pi) Coinc Cuts");
	Prod5816VZ->SetXTitle("ReactPt_L.z (m)");

	
// 	const Double_t VZAcceptance = 0.185;
	
	Double_t NPassVzCut  = 0;
	TH1F* h = 0;
	h=Prod5816VZ;
	for(UInt_t i=0;i<h->GetNbinsX(); i++)
		if (TMath::Abs(h->GetBinCenter(i))<=VZAcceptance) NPassVzCut+=h->GetBinContent(i);
	
	////////////////////////////////////////////////////////////////////////////////
	TCanvas * c1 = new TCanvas("AngularCheck1","AngularCheck1",1900,1000);
	c1->Divide(2,1);
	UInt_t idx=1;
	
	c1->cd(idx++);
	t->Draw("ReactPt_L.z>>Prod5816VZNoCuts(800,-.3,.3)","L.tr.n==1");
	Prod5816VZNoCuts->SetTitle("Production Runs LHRS React VZ (No Cuts)");
	Prod5816VZNoCuts->SetXTitle("ReactPt_L.z (m)");
	TLine *l = new TLine(-VZAcceptance,-0,-VZAcceptance,+Prod5816VZNoCuts->GetMaximum());
	l->SetLineColor(4);
	l->Draw();
	TLine *l = new TLine(VZAcceptance,-0,VZAcceptance,+Prod5816VZNoCuts->GetMaximum());
	l->SetLineColor(4);
	l->Draw();

	c1->cd(idx++);
	Prod5816VZ->Rebin(18);
	Prod5816VZ->Draw("");
	
	TLine *l = new TLine(-VZAcceptance,-0,-VZAcceptance,+Prod5816VZ->GetMaximum());
	l->SetLineColor(4);
	l->Draw();
	TLine *l = new TLine(VZAcceptance,-0,VZAcceptance,+Prod5816VZ->GetMaximum());
	l->SetLineColor(4);
	l->Draw();

	
	TText *text = new TText(0,Prod5816VZ->GetMaximum()/2,Form("Cut Acceptance ~ %2.1f%%",100*NPassVzCut/NPassCut));
	text->SetTextSize(0.05);
	text->SetTextAlign(22);
	text->SetTextColor(4);
	text->Draw();
	
	c1->Print(Form("VertexAcceptanceCheck.%s.png",RootFileID.Data()),"png");

	////////////////////////////////////////////////////////////////////////////////

}

VertexDependancyCheck(TString RootFileID = "Prod2")
{
	
	gStyle->SetOptStat(11);
	
	
	TChain * t=new TChain("T");
	RootFileLoader(t,"./ROOTfiles/e06010_prod3_%d.root",5816,5817);
// 	RootFileLoader(t,"./ROOTfiles/e06010_prod3_%d.root",5805,5815);


	////////////////////////////////////////////////////////////////////////////////
	TCanvas * c1 = new TCanvas("AngularCheck1","AngularCheck1",1900,1000);
	c1->Divide(3,2);
	UInt_t idx=1;
	
	c1->cd(idx++);
	t->Draw("ReactPt_L.z>>Prod5816VZ(800,-.3,.3)","L.tr.n==1","");
	Prod5816VZ->SetTitle("React Z");
	Prod5816VZ->SetXTitle("ReactPt_L.z (m)");
	Prod5816VZ->SetYTitle("Delta Momentum L.tr.dp");
	
	TLine *l = new TLine(-VZAcceptance,-0,-VZAcceptance,10000);
	l->SetLineColor(6);
	l->Draw();
	TLine *l = new TLine(VZAcceptance,-0,VZAcceptance,10000);
	l->SetLineColor(6);
	l->Draw();
	
	c1->cd(idx++);
	t->Draw("L.tr.tg_dp:ReactPt_L.z>>Prod5816VZdp(800,-.3,.3,800,-.06,.06)","L.tr.n==1","COLZ");
	Prod5816VZdp->SetTitle("Target DP");
	Prod5816VZdp->SetXTitle("ReactPt_L.z (m)");
	Prod5816VZdp->SetYTitle("Delta Momentum L.tr.dp");
	
	TLine *l = new TLine(-VZAcceptance,-10000,-VZAcceptance,10000);
	l->SetLineColor(6);
	l->Draw();
	TLine *l = new TLine(VZAcceptance,-10000,VZAcceptance,10000);
	l->SetLineColor(6);
	l->Draw();

	
	c1->cd(idx++);
	t->Draw("L.tr.r_y:ReactPt_L.z>>Prod5816VZy(800,-.3,.3,800,-.1,.1)","L.tr.n==1","COLZ");
	Prod5816VZy->SetTitle("Rotated y");
	Prod5816VZy->SetXTitle("ReactPt_L.z (m)");
	Prod5816VZy->SetYTitle("Rotated y L.tr.r_y (m)");
	
	TLine *l = new TLine(-VZAcceptance,-10000,-VZAcceptance,10000);
	l->SetLineColor(6);
	l->Draw();
	TLine *l = new TLine(VZAcceptance,-10000,VZAcceptance,10000);
	l->SetLineColor(6);
	l->Draw();

	
	c1->cd(idx++);
	t->Draw("L.tr.r_ph:ReactPt_L.z>>Prod5816VZph(800,-.3,.3,800,-.06,.06)","L.tr.n==1","COLZ");
	Prod5816VZph->SetTitle("Rotated Phi");
	Prod5816VZph->SetXTitle("ReactPt_L.z (m)");
	Prod5816VZph->SetYTitle("Rotated Phi L.tr.r_ph");
	
	TLine *l = new TLine(-VZAcceptance,-10000,-VZAcceptance,10000);
	l->SetLineColor(6);
	l->Draw();
	TLine *l = new TLine(VZAcceptance,-10000,VZAcceptance,10000);
	l->SetLineColor(6);
	l->Draw();

	
	c1->cd(idx++);
	t->Draw("L.tr.tg_th:ReactPt_L.z>>Prod5816VZth(800,-.3,.3,800,-.08,.08)","L.tr.n==1","COLZ");
	Prod5816VZth->SetTitle("Target Theta");
	Prod5816VZth->SetXTitle("ReactPt_L.z (m)");
	Prod5816VZth->SetYTitle("Target Theta L.tr.tg_th");
	
	TLine *l = new TLine(-VZAcceptance,-10000,-VZAcceptance,10000);
	l->SetLineColor(6);
	l->Draw();
	TLine *l = new TLine(VZAcceptance,-10000,VZAcceptance,10000);
	l->SetLineColor(6);
	l->Draw();
	c1->cd(idx++);
	
	t->Draw("L.tr.tg_ph:ReactPt_L.z>>Prod5816VZtgph(800,-.3,.3,800,-.05,.05)","L.tr.n==1","COLZ");
	Prod5816VZtgph->SetTitle("Target Phi");
	Prod5816VZtgph->SetXTitle("ReactPt_L.z (m)");
	Prod5816VZtgph->SetYTitle("Target Theta L.tr.tg_ph");
	
	TLine *l = new TLine(-VZAcceptance,-10000,-VZAcceptance,10000);
	l->SetLineColor(6);
	l->Draw();
	TLine *l = new TLine(VZAcceptance,-10000,VZAcceptance,10000);
	l->SetLineColor(6);
	l->Draw();

	
	c1->Print(Form("VertexDependancyCheck.%s.png",RootFileID.Data()),"png");

	////////////////////////////////////////////////////////////////////////////////

}

SieveCutFile()
{
	fstream fcut("Sieve.cut",ios_base::out);
	assert(fcut.is_open());
	
	TFile * f = new TFile("./cutfiles/gcut_L_0.root","");
	assert(f);
	
	for (UInt_t FoilID=1;FoilID<=6;FoilID++)
	{
		
		TString VertexCut = Form("fcut_L_%d",FoilID);
		cout<<"Looking for "<<VertexCut<<endl;
		TCutG *cVertexCut = (TCutG*)gROOT->FindObject(VertexCut);
		assert(cVertexCut);
		for (UInt_t Col=1;Col<=7;Col++)
			for (UInt_t Row=1;Row<=7;Row++)
		{
			TString SieveCut = Form("hcut_L_%d_%d_%d",FoilID,Col,Row);
			TCutG *cSieveCut = (TCutG*)gROOT->FindObject(SieveCut);
			if (!cSieveCut)
			{
					//write an empty cut
				fcut<<"fEvtHdr.fRun==0"<<endl;
			}
			else
			{
				fcut<<VertexCut<<"&&"<<SieveCut
						<<"&&abs(L.gold.dp)<0.05&&L.tr.n==1"<<endl;
			}
		}
	}
	
	fcut.close();
}

DpCutFile()
{
	fstream fcut("Dp.cut",ios_base::out);
	assert(fcut.is_open());
	
	TFile * f = new TFile("./cutfiles/gcut_L_dpfull.root","");
	assert(f);

	const UInt_t NConfig = 5;
	TString ConfigCut[]={
		"fEvtHdr.fRun==2528&&L_dp_m4",
  "fEvtHdr.fRun==2529&&L_dp_m2",
  "(fEvtHdr.fRun==2530||fEvtHdr.fRun==2531)&&L_dp_0",
  "fEvtHdr.fRun==2532&&L_dp_p2",
  "(fEvtHdr.fRun==2533||fEvtHdr.fRun==2534)&&L_dp_p4"
	};
	TString AntiConfigCut[]={
		"fEvtHdr.fRun==2528&&!L_dp_m4",
  "fEvtHdr.fRun==2529&&!L_dp_m2",
  "(fEvtHdr.fRun==2530||fEvtHdr.fRun==2531)&&!L_dp_0",
  "fEvtHdr.fRun==2532&&!L_dp_p2",
  "(fEvtHdr.fRun==2533||fEvtHdr.fRun==2534)&&!L_dp_p4"
	};
	
	assert(gROOT->FindObject("L_dp_m4"));
	assert(gROOT->FindObject("L_dp_m2"));
	assert(gROOT->FindObject("L_dp_0"));
	assert(gROOT->FindObject("L_dp_p2"));
	assert(gROOT->FindObject("L_dp_p4"));

	// first out put major events in peaks
	for (UInt_t ConfigID = 0; ConfigID<NConfig; ConfigID++)
	{
		for (UInt_t FoilID=1;FoilID<=6;FoilID++)
		{
			TString VertexCut = Form("fcut_L_%d",FoilID);
			cout<<"Looking for "<<VertexCut<<endl;
			TCutG *cVertexCut = (TCutG*)gROOT->FindObject(VertexCut);
			assert(cVertexCut);
			for (UInt_t Col=1;Col<=7;Col++)
				for (UInt_t Row=1;Row<=7;Row++)
			{
				TString SieveCut = Form("hcut_L_%d_%d_%d",FoilID,Col,Row);
				TCutG *cSieveCut = (TCutG*)gROOT->FindObject(SieveCut);
				if (!cSieveCut)
				{
						//write an empty cut
					fcut<<"fEvtHdr.fRun==0"<<endl;
				}
				else
				{
					fcut<<VertexCut<<"&&"<<SieveCut
							<<"&&abs(L.gold.dp)<0.05&&L.tr.n==1&&"<<ConfigCut[ConfigID]<<endl;
				}
			}
		}
	}
	// then secondary events outside peaks
	for (UInt_t ConfigID = 0; ConfigID<NConfig; ConfigID++)
	{
		for (UInt_t FoilID=1;FoilID<=6;FoilID++)
		{
			TString VertexCut = Form("fcut_L_%d",FoilID);
			cout<<"Looking for "<<VertexCut<<endl;
			TCutG *cVertexCut = (TCutG*)gROOT->FindObject(VertexCut);
			assert(cVertexCut);
			for (UInt_t Col=1;Col<=7;Col++)
				for (UInt_t Row=1;Row<=7;Row++)
			{
				TString SieveCut = Form("hcut_L_%d_%d_%d",FoilID,Col,Row);
				TCutG *cSieveCut = (TCutG*)gROOT->FindObject(SieveCut);
				if (!cSieveCut)
				{
						//write an empty cut
					fcut<<"fEvtHdr.fRun==0"<<endl;
				}
				else
				{
					fcut<<VertexCut<<"&&"<<SieveCut
							<<"&&abs(L.gold.dp)<0.05&&L.tr.n==1&&"<<AntiConfigCut[ConfigID]<<endl;
				}
			}
		}
	}

	fcut.close();
}

DpCutFileNew()
{
	fstream fcut("DpNew.cut",ios_base::out);
	assert(fcut.is_open());
	
	TFile * f = new TFile("./cutfiles/gcut_L_dpfull.root","");
	assert(f);

	const UInt_t NConfig = 5;
	TString ConfigCut[]={
		"fEvtHdr.fRun==2528",
  "fEvtHdr.fRun==2529",
  "(fEvtHdr.fRun==2530||fEvtHdr.fRun==2531)",
  "fEvtHdr.fRun==2532",
  "(fEvtHdr.fRun==2533||fEvtHdr.fRun==2534)"
	};
	const TString KineName[5] = {"m4","m2","0","p2","p4"};

	
	// first out put major events in peaks
	for (UInt_t ConfigID = 0; ConfigID<NConfig; ConfigID++)
	{
		for (UInt_t FoilID=1;FoilID<=6;FoilID++)
		{
			TString VertexCut = Form("fcut_L_%d",FoilID);
// 			cout<<"Looking for "<<VertexCut<<endl;
			TCutG *cVertexCut = (TCutG*)gROOT->FindObject(VertexCut);
			assert(cVertexCut);
			
			TString DpCut = Form("L_dp_%s_%d",KineName[ConfigID].Data(),FoilID);
// 			cout<<"Looking for "<<VertexCut<<endl;
			TCutG *cVertexCut = (TCutG*)gROOT->FindObject(VertexCut);
			assert(cVertexCut);

			cout<<"Output for "<<ConfigCut[ConfigID]<<" and "<<VertexCut <<" and "<< DpCut<<endl;
			
			for (UInt_t Col=1;Col<=7;Col++)
				for (UInt_t Row=1;Row<=7;Row++)
			{
				TString SieveCut = Form("hcut_L_%d_%d_%d",FoilID,Col,Row);
				TCutG *cSieveCut = (TCutG*)gROOT->FindObject(SieveCut);
				if (!cSieveCut)
				{
						//write an empty cut
					fcut<<"fEvtHdr.fRun==0"<<endl;
				}
				else
				{
					fcut<<VertexCut<<"&&"<<SieveCut
							<<"&&abs(L.gold.dp)<0.05&&L.tr.n==1&&"<<ConfigCut[ConfigID]
							<<"&&"<<DpCut<<endl;
				}
			}
		}
	}
// 	// then secondary events outside peaks
	for (UInt_t ConfigID = 0; ConfigID<NConfig; ConfigID++)
	{
		for (UInt_t FoilID=1;FoilID<=6;FoilID++)
		{
			TString VertexCut = Form("fcut_L_%d",FoilID);
// 			cout<<"Looking for "<<VertexCut<<endl;
			TCutG *cVertexCut = (TCutG*)gROOT->FindObject(VertexCut);
			assert(cVertexCut);
			
			TString DpCut = Form("L_dp_%s_%d",KineName[ConfigID].Data(),FoilID);
// 			cout<<"Looking for "<<VertexCut<<endl;
			TCutG *cVertexCut = (TCutG*)gROOT->FindObject(VertexCut);
			assert(cVertexCut);

			cout<<"Output anti cut for "<<ConfigCut[ConfigID]<<" and "<<VertexCut <<" and "<< DpCut<<endl;
			
			for (UInt_t Col=1;Col<=7;Col++)
				for (UInt_t Row=1;Row<=7;Row++)
			{
				TString SieveCut = Form("hcut_L_%d_%d_%d",FoilID,Col,Row);
				TCutG *cSieveCut = (TCutG*)gROOT->FindObject(SieveCut);
				if (!cSieveCut)
				{
						//write an empty cut
					fcut<<"fEvtHdr.fRun==0"<<endl;
				}
				else
				{
					fcut<<VertexCut<<"&&"<<SieveCut
							<<"&&abs(L.gold.dp)<0.05&&L.tr.n==1&&"<<ConfigCut[ConfigID]
							<<"&&!("<<DpCut<<")"<<endl;
				}
			}
		}
	}

	fcut.close();
}

DrawOpticsData()
{
	TFile *fcut, *fData;
	TTree *t;
	
	fdata = new TFile("ROOTfiles/e06010_det_L_5832.root");
	t=T;
	
	fcut = new TFile("./cutfiles/gcut_L_1.root","");

	
}

DrawDeltaScan()
{
	
	TCanvas *c1 = new TCanvas("DeltaScan","DeltaScan",1900,1100);

	c1->Divide(2,5);
	UInt_t idx=1;
	TFile *_file0;
	
	_file0 = TFile::Open("ROOTfiles/e06010_det_L_2533.root");
	c1->cd(idx++);
	T->Draw("L.tr.tg_dp>>Dp_m4(500,-.05,.05)");
	c1->cd(idx++);
	T->Draw("L.tr.p>>p_m4(100,1.210,1.230)");
	
	_file0 = TFile::Open("ROOTfiles/e06010_det_L_2532.root");
	c1->cd(idx++);
	T->Draw("L.tr.tg_dp>>Dp_m2(500,-.05,.05)");
	c1->cd(idx++);
	T->Draw("L.tr.p>>p_m2(100,1.210,1.230)");
	
	_file0 = TFile::Open("ROOTfiles/e06010_det_L_2531.root");
	c1->cd(idx++);
	T->Draw("L.tr.tg_dp>>Dp_0(500,-.05,.05)");
	c1->cd(idx++);
	T->Draw("L.tr.p>>p_0(100,1.210,1.230)");
	
	_file0 = TFile::Open("ROOTfiles/e06010_det_L_2529.root");
	c1->cd(idx++);
	T->Draw("L.tr.tg_dp>>Dp_p2(500,-.05,.05)");
	c1->cd(idx++);
	T->Draw("L.tr.p>>p_p2(100,1.210,1.230)");
	
	_file0 = TFile::Open("ROOTfiles/e06010_det_L_2528.root");
	c1->cd(idx++);
	T->Draw("L.tr.tg_dp>>Dp_p4(500,-.05,.05)");
	c1->cd(idx++);
	T->Draw("L.tr.p>>p_p4(100,1.210,1.230)");

}

DrawVertexCoinc()
{

	
	TChain * T=new TChain("T");
// 	RootFileLoader(T,"./ROOTfiles/e06010_prod4_%d.root",5800,5816);
// 	RootFileLoader(T,"./ROOTfiles/e06010_prod3_%d.root",5816,5816);
	RootFileLoader(T,"/raid2/jinhuang/e06010/Production4/e06010_prod4_%d.root",3000,6000);

	Double_t VZLimit = .3;

	TCut cut = "L.tr.n==1 && BB.tr.n ==1 && abs(CT.pi.t)<.5 && L.cer.asum_c<10&&L.a1.asum_c>200&&(L.prl1.e+L.prl2.e)/2350./(1+L.tr.tg_dp)>0.04&&(L.prl1.e+L.prl2.e)/2350./(1+L.tr.tg_dp)<0.12 && !DL.edtpl&&(DL.evtypebits&1<<5)==1<<5 &&  BB.tr.chi2[BB.ts.sh_flag_tr[]]/BB.tr.ndof[BB.ts.sh_flag_tr[]]<2.4 && BB.tr.n > 0 && BB.tr.p[BB.ts.sh_flag_tr[]]>0.6 && BB.tr.p[BB.ts.sh_flag_tr[]]<2.2  && (BB.optics.charge[BB.ts.sh_flag_tr[]] ==-1 ||BB.optics.charge[BB.ts.sh_flag_tr[]] ==1) && abs(BB.tr.vz[BB.ts.sh_flag_tr[]])<0.2 && BB.ts.ps_e[]>350. && abs((BB.ts.ps_e[]/2+BB.ts.sh_e[])/BB.tr.p[BB.ts.sh_flag_tr[]]-1000)<270";


	TCanvas *c1 = new TCanvas("LOpticsBBLHRSComp","LOpticsBBLHRSComp",1000,1100);

	c1->Divide(1,2);
	UInt_t idx=1;

	c1->cd(idx++);
	T->Draw("ReactPt_L.z:BB.tr.vz-ReactPt_L.z>>VZ2D(200,-.1,.1,200,-.3,.3)"
			,cut,"SCAT");
	VZ2D->SetTitle("Compare of LHRS and Bigbite Vertex, (e'pi) Cuts");
	VZ2D->SetXTitle("Vertex Diff = BB.tr.vz-ReactPt_L.z (m)");
	VZ2D->SetYTitle("LHRS Vertex = ReactPt_L.z (m)");
	VZ2D->SetMarkerSize(2);

	
	TLine *l = new TLine(0,-VZLimit,0,VZLimit);
	l->SetLineColor(6);
	l->Draw();
	

	
	c1->cd(idx++);
	VZ2D -> ProjectionX("VZCoinc");
	VZCoinc->Rebin(2);

	VZCoinc->Draw();
	VZCoinc->SetXTitle("Vertex Diff = BB.tr.vz-ReactPt_L.z (m)");

	const Double_t fitrange = .04;
	TF1 *f = new TF1("gausl","gaus+[3]",-fitrange,fitrange);
	f->SetParameter(2,.02);

	VZCoinc->Fit("gausl","R");
	
	TText *t = new TText(f->GetParameter(2)
			,f->GetParameter(0)-f->GetParameter(3)
			,Form("BBite- LHRS = %2.1f +- %2.1f cm"
					,100*(f->GetParameter(1))
							,100*f->GetParameter(2)));
	t->SetTextSize(0.05);
	t->SetTextAlign(13);
	t->SetTextColor(4);
	t->Draw();
	
	TLine *l = new TLine(0,0,0,10000);
	l->SetLineColor(6);
	l->Draw();

	SaveCanvas(c1,c1->GetName());
	
}

DrawVertexCoincStat()
{

	
	TChain * T=new TChain("T");
	UInt_t startrun = 3000, endrun =6000;
	RootFileLoader(T,"./ROOTfiles/e06010_prod4_%d.root",startrun,endrun);
// 	RootFileLoader(T,"./ROOTfiles/e06010_prod3_%d.root",5816,5816);

	Double_t VZLimit = .3;

	TCut cut = "L.tr.n==1 && BB.tr.n ==1 && abs(CT.pi.t)<.5 && L.cer.asum_c<10&&L.a1.asum_c>200&&(L.prl1.e+L.prl2.e)/2350./(1+L.tr.tg_dp)>0.04&&(L.prl1.e+L.prl2.e)/2350./(1+L.tr.tg_dp)<0.12 && !DL.edtpl&&(DL.evtypebits&1<<5)==1<<5 &&  BB.tr.chi2[BB.ts.sh_flag_tr[]]/BB.tr.ndof[BB.ts.sh_flag_tr[]]<2.4 && BB.tr.n > 0 && BB.tr.p[BB.ts.sh_flag_tr[]]>0.6 && BB.tr.p[BB.ts.sh_flag_tr[]]<2.2  && (BB.optics.charge[BB.ts.sh_flag_tr[]] ==-1 ||BB.optics.charge[BB.ts.sh_flag_tr[]] ==1) && abs(BB.tr.vz[BB.ts.sh_flag_tr[]])<0.2 && BB.ts.ps_e[]>350. && abs((BB.ts.ps_e[]/2+BB.ts.sh_e[])/BB.tr.p[BB.ts.sh_flag_tr[]]-1000)<270";


	TCanvas *c1 = new TCanvas("LOpticsBBLHRSCompStat","LOpticsBBLHRSComp",1000,1100);

	c1->Divide(1,2);
	UInt_t idx=1;

	c1->cd(idx++);
	T->Draw(Form("BB.tr.vz-ReactPt_L.z:fEvtHdr.fRun>>VZ2D(%d,%d,%d,50,-.1,.1)"
			,endrun/10-startrun/10+1,startrun/10*10,endrun/10*10+10),cut,"COLZ");
	VZ2D->SetTitle("Compare of LHRS and Bigbite Vertex, (e'pi) Cuts");
	VZ2D->SetYTitle("Vertex Diff = BB.tr.vz-ReactPt_L.z (m)");
	VZ2D->SetXTitle("Production Run number");
	VZ2D->SetMarkerSize(2);

	
	TLine *l = new TLine(0,-VZLimit,0,VZLimit);
	l->SetLineColor(6);
	l->Draw();
	
	TLine *l = new TLine(3000,0,6000,0);
	l->SetLineColor(6);
	l->Draw();

	
	c1->cd(idx++);
	VZ2D -> ProjectionY("VZCoinc");
	VZCoinc->Rebin(2);

	VZCoinc->Draw();
	VZCoinc->SetXTitle("Vertex Diff = BB.tr.vz-ReactPt_L.z (m)");

	const Double_t fitrange = .04;
	TF1 *f = new TF1("gausl","gaus+[3]",-fitrange,fitrange);
	f->SetParameter(2,.02);

	VZCoinc->Fit("gausl","R");
	
	TText *t = new TText(f->GetParameter(2)
			,f->GetParameter(0)-f->GetParameter(3)
			,Form("BBite- LHRS = %2.1f +- %2.1f cm"
					,100*(f->GetParameter(1))
							,100*f->GetParameter(2)));
	t->SetTextSize(0.05);
	t->SetTextAlign(13);
	t->SetTextColor(4);
	t->Draw();
	
	TLine *l = new TLine(0,0,0,10000);
	l->SetLineColor(6);
	l->Draw();
	
}

DrawVertexCoincStatByBBp()
{

	
	TChain * T=new TChain("T");
	UInt_t startrun = 4000, endrun =6000;
	RootFileLoader(T,"./ROOTfiles/e06010_prod4_%d.root",startrun,endrun);
// 	RootFileLoader(T,"./ROOTfiles/e06010_prod3_%d.root",5816,5816);

	Double_t VZLimit = .3;

	TCut cut = "L.tr.n==1 && BB.tr.n ==1 && abs(CT.pi.t)<.5 && L.cer.asum_c<10&&L.a1.asum_c>200&&(L.prl1.e+L.prl2.e)/2350./(1+L.tr.tg_dp)>0.04&&(L.prl1.e+L.prl2.e)/2350./(1+L.tr.tg_dp)<0.12 && !DL.edtpl&&(DL.evtypebits&1<<5)==1<<5 &&  BB.tr.chi2[BB.ts.sh_flag_tr[]]/BB.tr.ndof[BB.ts.sh_flag_tr[]]<2.4 && BB.tr.n > 0 && BB.tr.p[BB.ts.sh_flag_tr[]]>0.6 && BB.tr.p[BB.ts.sh_flag_tr[]]<2.2  && (BB.optics.charge[BB.ts.sh_flag_tr[]] ==-1 ||BB.optics.charge[BB.ts.sh_flag_tr[]] ==1) && abs(BB.tr.vz[BB.ts.sh_flag_tr[]])<0.2 && BB.ts.ps_e[]>350. && abs((BB.ts.ps_e[]/2+BB.ts.sh_e[])/BB.tr.p[BB.ts.sh_flag_tr[]]-1000)<270";


	TCanvas *c1 = new TCanvas("LOpticsBBLHRSCompStatByDp","LOpticsBBLHRSComp",1400,1000);

	const Double_t fitrange = .04;
	TF1 *f = new TF1("gausl","gaus+[3]",-fitrange,fitrange);
	f->SetParameter(2,.02);
	
	c1->Divide(2,2);
	UInt_t idx=1;

	
	
	c1->cd(idx++);
	T->Draw("BB.tr.vz-ReactPt_L.z>>BBp600(200,-.1,.1)",cut && "abs(BB.tr.p-0.6)<0.3","");
	BBp600->SetTitle("BB 0.6 GeV/c : Coinc of LHRS and Bigbite Vertex, (e'pi) Cuts");
	BBp600->SetXTitle("Vertex Diff = BB.tr.vz-ReactPt_L.z (m)");
	BBp600->Fit("gausl","R");
	
	TText *t = new TText(f->GetParameter(2)
			,f->GetParameter(0)-f->GetParameter(3)
			,Form("BBite- LHRS = %2.1f +- %2.1f cm"
					,100*(f->GetParameter(1))
							,100*f->GetParameter(2)));
	t->SetTextSize(0.05);
	t->SetTextAlign(13);
	t->SetTextColor(4);
	t->Draw();
	
	TLine *l = new TLine(0,0,0,10000);
	l->SetLineColor(6);
	l->Draw();
	
	
	c1->cd(idx++);
	T->Draw("BB.tr.vz-ReactPt_L.z>>BBp1200(200,-.1,.1)",cut && "abs(BB.tr.p-1.2)<0.3","");
	BBp1200->SetTitle("BB 1.2 GeV/c : Coinc of LHRS and Bigbite Vertex, (e'pi) Cuts");
	BBp1200->SetXTitle("Vertex Diff = BB.tr.vz-ReactPt_L.z (m)");
	BBp1200->Fit("gausl","R");
	
	TText *t = new TText(f->GetParameter(2)
			,f->GetParameter(0)-f->GetParameter(3)
			,Form("BBite- LHRS = %2.1f +- %2.1f cm"
					,100*(f->GetParameter(1))
							,100*f->GetParameter(2)));
	t->SetTextSize(0.05);
	t->SetTextAlign(13);
	t->SetTextColor(4);
	t->Draw();
	
	TLine *l = new TLine(0,0,0,10000);
	l->SetLineColor(6);
	l->Draw();
	
	
	c1->cd(idx++);
	T->Draw("BB.tr.vz-ReactPt_L.z>>BBp1800(200,-.1,.1)",cut && "abs(BB.tr.p-1.8)<0.3","");
	BBp1800->SetTitle("BB 1.8 GeV/c : Coinc of LHRS and Bigbite Vertex, (e'pi) Cuts");
	BBp1800->SetXTitle("Vertex Diff = BB.tr.vz-ReactPt_L.z (m)");
	BBp1800->Fit("gausl","R");
	
	TText *t = new TText(f->GetParameter(2)
			,f->GetParameter(0)-f->GetParameter(3)
			,Form("BBite- LHRS = %2.1f +- %2.1f cm"
					,100*(f->GetParameter(1))
							,100*f->GetParameter(2)));
	t->SetTextSize(0.05);
	t->SetTextAlign(13);
	t->SetTextColor(4);
	t->Draw();
	
	TLine *l = new TLine(0,0,0,10000);
	l->SetLineColor(6);
	l->Draw();
	
	
	c1->cd(idx++);
	T->Draw("BB.tr.vz-ReactPt_L.z>>BBp2400(200,-.1,.1)",cut && "abs(BB.tr.p-2.4)<0.3","");
	BBp2400->SetTitle("BB 2.4 GeV/c : Coinc of LHRS and Bigbite Vertex, (e'pi) Cuts");
	BBp2400->SetXTitle("Vertex Diff = BB.tr.vz-ReactPt_L.z (m)");
	BBp2400->Fit("gausl","R");
	
	TText *t = new TText(f->GetParameter(2)
			,f->GetParameter(0)-f->GetParameter(3)
			,Form("BBite- LHRS = %2.1f +- %2.1f cm"
					,100*(f->GetParameter(1))
							,100*f->GetParameter(2)));
	t->SetTextSize(0.05);
	t->SetTextAlign(13);
	t->SetTextColor(4);
	t->Draw();
	
	TLine *l = new TLine(0,0,0,10000);
	l->SetLineColor(6);
	l->Draw();
	
}

OpticsStability()
{

	
	TChain * T=new TChain("T");
	UInt_t startrun = 3000, endrun =6000;
	RootFileLoader(T,"ProductionROOTfiles/e06010_prod4_%d.root",startrun,endrun);
	
	TCanvas *c1 = new TCanvas("OpticsStability","OpticsStability",1900,1000);

	c1->Divide(2,2);
	UInt_t idx=1;

	const Int_t sample = 1;

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_dp:fEvtHdr.fRun>>tg_dp(%d,%d,%d,100,-.1,.1)"
			,endrun/sample-startrun/sample+1,startrun/sample*sample,endrun/sample*sample+sample),"","COLZ");
	tg_dp->SetTitle("Target dp : HRS Optics Stability Check");
	tg_dp->SetYTitle("Target dp");
	tg_dp->SetXTitle("Production Run number");
	

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_y:fEvtHdr.fRun>>tg_y(%d,%d,%d,100,-.2,.2)"
			,endrun/sample-startrun/sample+1,startrun/sample*sample,endrun/sample*sample+sample),"","COLZ");
	tg_y->SetTitle("Target Y/m : HRS Optics Stability Check");
	tg_y->SetYTitle("Target Y/m");
	tg_y->SetXTitle("Production Run number");

	
	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_th:fEvtHdr.fRun>>tg_th(%d,%d,%d,100,-.12,.2)"
			,endrun/sample-startrun/sample+1,startrun/sample*sample,endrun/sample*sample+sample),"","COLZ");
	tg_th->SetTitle("Target Theta : HRS Optics Stability Check");
	tg_th->SetYTitle("Target Theta");
	tg_th->SetXTitle("Production Run number");
	

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_dp:fEvtHdr.fRun>>tg_ph(%d,%d,%d,100,-.1,.1)"
			,endrun/sample-startrun/sample+1,startrun/sample*sample,endrun/sample*sample+sample),"","COLZ");
	tg_ph->SetTitle("Target Phi : HRS Optics Stability Check");
	tg_ph->SetYTitle("Target Phi");
	tg_ph->SetXTitle("Production Run number");


	TCanvas *c1 = new TCanvas("OpticsStabilityStat","OpticsStabilityStat",1900,1000);

	c1->Divide(2,2);
	UInt_t idx=1;

	c1->cd(idx++);
	tg_dp -> ProfileX();
	tg_dp_pfx -> Draw();
	
	c1->cd(idx++);
	tg_y -> ProfileX();
	tg_y_pfx -> Draw();
	
	c1->cd(idx++);
	tg_th -> ProfileX();
	tg_th_pfx -> Draw();
	
	c1->cd(idx++);
	tg_ph -> ProfileX();
	tg_ph_pfx -> Draw();
	
	
}

OpticsCellWindowStability()
{

	
	TChain * T=new TChain("T");
	UInt_t startrun = 3000, endrun =6000;
	RootFileLoader(T,"ProductionROOTfiles/e06010_prod4_%d.root",startrun,endrun);
	
	TCanvas *c1 = new TCanvas("OpticsCellWindowStability","OpticsCellWindowStability",1000,1200);

	c1->Divide(1,4);
	UInt_t idx=1;


	c1->cd(idx++);
	T->Draw(Form("ReactPt_L.z:fEvtHdr.fRun>>upwall(%d,%d,%d,100,-.22,-.18)"
			,endrun-startrun+1,startrun,endrun+1),"","COLZ");
	upwall->SetTitle("UpStream Cell Window : HRS Optics Stability Check");
	upwall->SetYTitle("React Z/m");
	upwall->SetXTitle("Production Run number");

	gPad->SetLogz();

	c1->cd(idx++);
	T->Draw(Form("ReactPt_L.z:fEvtHdr.fRun>>downwall(%d,%d,%d,100,.18,.22)"
			,endrun-startrun+1,startrun,endrun+1),"","COLZ");
	downwall->SetTitle("DownStream Cell Window : HRS Optics Stability Check");
	downwall->SetYTitle("React Z/m");
	downwall->SetXTitle("Production Run number");

	gPad->SetLogz();
	
	c1->cd(idx++);
	T->Draw(Form("urb.x:fEvtHdr.fRun>>urbx(%d,%d,%d,100,-0.01,0.01)"
			,endrun-startrun+1,startrun,endrun+1),"","COLZ");
	downwall->SetTitle("Beam Horizontal Position");
	downwall->SetYTitle("Beam Horizontal Position/m");
	downwall->SetXTitle("Production Run number");
	c1->cd(idx++);
	T->Draw(Form("urb.y:fEvtHdr.fRun>>urby(%d,%d,%d,100,-0.01,0.01)"
			,endrun-startrun+1,startrun,endrun+1),"","COLZ");
	downwall->SetTitle("Beam Vertical Position");
	downwall->SetYTitle("Beam Vertical Position/m");
	downwall->SetXTitle("Production Run number");
}

DrawH2Test()
{
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2620.root");
// 	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2663.root");
	
// 	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1 && rb.y<0.004 && abs(rb.x+0.001)<0.00005");
	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1 && rb.y<0.004 ");
// 	T->Draw(">>TestCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1");

// 	TString CutUsed = "LooseCut";
// 	TString CutUsed = "DeepCut";
	TString CutUsed = "ProdCut";
// 	TString CutUsed = "TestCut";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);

// 	T->Draw("(( 2*0.938272*rb.e)/(2*(0.938272 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21269-1>>(100,-.5,.5)","","");
	TCanvas * c1 = new TCanvas(CutUsed+"_Canvas",CutUsed,1900,1000);
	c1->Divide(3,3);
	UInt_t idx=1;


// 	c1->cd(idx++);
	
// 	T->Draw("((( 2*0.938272*rb.e)/(2*(0.938272 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2)))))/1.21284789137319704-1>>asdf(100,-.05,.05)","","");
// 	c1->cd(idx++);
// 	T->Draw("ExTgtCor_L.dp>>asdfasdf(100,-.05,.05)","","");
// 	c1->cd(idx++);
// 	T->Draw("ExTgtCor_L.dp-((( 2*0.938272*rb.e)/(2*(0.938272 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>asdff2(100,-.01,.01)","","");
	c1->cd(idx++);
	T->Draw("rb.x:ExTgtCor_L.dp-((( 2*0.938272*rb.e)/(2*(0.938272 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>rb_x(1000,-.01,.01,1000,-.004,.004)","","COLZ");
	c1->cd(idx++);
	T->Draw("rb.y:ExTgtCor_L.dp-((( 2*0.938272*rb.e)/(2*(0.938272 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>rb_y(1000,-.01,.01,1000,-.002,.008)","","COLZ");
	c1->cd(idx++);
	T->Draw("ReactPt_L.z:ExTgtCor_L.dp-((( 2*0.938272*rb.e)/(2*(0.938272 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>ReactPt_L_z(1000,-.01,.01,1000,-.25,.25)","","COLZ");
	c1->cd(idx++);
	T->Draw("ExTgtCor_L.th:ExTgtCor_L.dp-((( 2*0.938272*rb.e)/(2*(0.938272 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>ExTgtCor_L_th(1000,-.01,.01,1000,-.07,.07)","","COLZ");
	c1->cd(idx++);
	T->Draw("L.tr.tg_ph:ExTgtCor_L.dp-((( 2*0.938272*rb.e)/(2*(0.938272 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>tg_ph(1000,-.01,.01,1000,-.04,.04)","","COLZ");
	c1->cd(idx++);
	T->Draw("L.tr.r_th:ExTgtCor_L.dp-((( 2*0.938272*rb.e)/(2*(0.938272 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>r_th(1000,-.01,.01,1000,-.04,.04)","","COLZ");
	c1->cd(idx++);
	T->Draw("L.tr.r_y:ExTgtCor_L.dp-((( 2*0.938272*rb.e)/(2*(0.938272 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>r_y(1000,-.01,.01,1000,-.05,.05)","","COLZ");
	c1->cd(idx++);
	T->Draw("L.tr.r_ph:ExTgtCor_L.dp-((( 2*0.938272*rb.e)/(2*(0.938272 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>r_ph(1000,-.01,.01,1000,-.05,.05)","","COLZ");
	c1->cd(idx++);
	T->Draw("ExTgtCor_L.x:ExTgtCor_L.dp-((( 2*0.938272*rb.e)/(2*(0.938272 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>ExTgtCor_L_x(1000,-.01,.01,1000,-.01,.01)","","COLZ");

}

DrawHe3Test()
{
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2766.root");
// 	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2663.root");
	
// 	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1 && rb.y<0.004 && abs(rb.x+0.001)<0.00005");
	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1");
// 	T->Draw(">>TestCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1");

// 	TString CutUsed = "LooseCut";
// 	TString CutUsed = "DeepCut";
	TString CutUsed = "ProdCut";
// 	TString CutUsed = "TestCut";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);

// 	T->Draw("(( 2*0.938272*rb.e)/(2*(0.938272 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21269-1>>(100,-.5,.5)","","");
	TCanvas * c1 = new TCanvas(CutUsed+"_He3Canvas",CutUsed,1900,1000);
// 	c1->Divide(3,1);
	c1->Divide(3,3);
	UInt_t idx=1;

// 	c1->cd(idx++);
// 	T->Draw("ExTgtCor_L.dp>>dp(1000,-.05,.05)","","");
// 
// 	c1->cd(idx++);
// 	T->Draw("((( 2*2.80839*rb.e)/(2*(2.80839 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2)))))/1.21284789137319704-1>>dp_exp_He3(100,-.05,.05)","","");
// 	c1->cd(idx++);
// 	T->Draw("ExTgtCor_L.dp-((( 2*2.80839*rb.e)/(2*(2.80839 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>dp_diff_He3(1000,-.05,.05)","","");
	
	c1->cd(idx++);
	T->Draw("rb.x:ExTgtCor_L.dp-((( 2*2.80839*rb.e)/(2*(2.80839 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>rb_x(1000,-.005,.005,1000,-.008,.008)","","COLZ");
	c1->cd(idx++);
	T->Draw("rb.y:ExTgtCor_L.dp-((( 2*2.80839*rb.e)/(2*(2.80839 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>rb_y(1000,-.005,.005,1000,-.008,.008)","","COLZ");
	c1->cd(idx++);
	T->Draw("ReactPt_L.z:ExTgtCor_L.dp-((( 2*2.80839*rb.e)/(2*(2.80839 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>ReactPt_L_z(1000,-.005,.005,1000,-.25,.25)","","COLZ");
	c1->cd(idx++);
	T->Draw("ExTgtCor_L.th:ExTgtCor_L.dp-((( 2*2.80839*rb.e)/(2*(2.80839 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>ExTgtCor_L_th(1000,-.005,.005,1000,-.07,.07)","","COLZ");
	c1->cd(idx++);
	T->Draw("L.tr.tg_ph:ExTgtCor_L.dp-((( 2*2.80839*rb.e)/(2*(2.80839 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>tg_ph(1000,-.005,.005,1000,-.04,.04)","","COLZ");
	c1->cd(idx++);
	T->Draw("L.tr.r_th:ExTgtCor_L.dp-((( 2*2.80839*rb.e)/(2*(2.80839 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>r_th(1000,-.005,.005,1000,-.04,.04)","","COLZ");
	c1->cd(idx++);
	T->Draw("L.tr.r_y:ExTgtCor_L.dp-((( 2*2.80839*rb.e)/(2*(2.80839 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>r_y(1000,-.005,.005,1000,-.05,.05)","","COLZ");
	c1->cd(idx++);
	T->Draw("L.tr.r_ph:ExTgtCor_L.dp-((( 2*2.80839*rb.e)/(2*(2.80839 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>r_ph(1000,-.005,.005,1000,-.05,.05)","","COLZ");
	c1->cd(idx++);
	T->Draw("ExTgtCor_L.x:ExTgtCor_L.dp-((( 2*2.80839*rb.e)/(2*(2.80839 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/1.21284789137319704-1)>>ExTgtCor_L_x(1000,-.005,.005,1000,-.01,.01)","","COLZ");

}

DrawH2()
{
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2620.root");
// 	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2663.root");
	
	T->Draw(">>LooseCut","abs(ReactPt_L.z)<.2 && abs(L.tr.tg_th)<.07 && abs(L.tr.tg_ph)<.04 && abs(L.tr.tg_dp)<.05 && L.tr.n==1");
	T->Draw(">>DeepCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.05 && abs(L.tr.tg_ph)<.03 && L.tr.n==1");
	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1 && abs(rb.y-0.0035)<0.0005 && abs(rb.x+0.001)<0.001");
// 	T->Draw(">>TestCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1");

// 	TString CutUsed = "LooseCut";
// 	TString CutUsed = "DeepCut";
	TString CutUsed = "ProdCut";
// 	TString CutUsed = "TestCut";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);

	////////////////////////////////////////////////////////////////////////////////
	
	const Double_t RealMom = DIPOLE_MAG2MOM(0.448924);
	const Double_t NMRMomCorr = RealMom-1.21272;
	cout<<"NMR Momentum Correction = "<<NMRMomCorr<<endl;

	const Double_t BeamEnergyLossOnDp
			= 0.90561*1e-3*(
							0.075//Beam window
			+0.010//He4 Before Cell
			+0.038//Cell Entrance
			+0.036//Half Cell Length
						   );
	const Double_t BeamEnergyLossOnDp_ReactZDep = 0.90561* 0.036e-3/.2;//GeV/m
	const Double_t AfterScatteringEnergyLoss
			= 1e-3*(
					3.62796/20*0.036 //in cell H2
			+0.038*10	//Exit from cell wall
			+0.0228		//He4
			+.1//Target Enclosure
			+0.109	//Air
			+.1//HRS Entrance
				   );
	TString RadCorrectionScatMom
			= Form("(%f+%f*ReactPt_L.z)"
			,BeamEnergyLossOnDp+AfterScatteringEnergyLoss
			,BeamEnergyLossOnDp_ReactZDep
				  );
	
	cout<<"Total Correction on p = "<<RadCorrectionScatMom<<endl;
	
	const Double_t RealDpKin = ( 2*0.938272013*1.2306)/(2*(0.938272013 + 1.2306 - 1.2306*TMath::Cos(16/180.*3.14159)))/RealMom-1;

	//DM = (sqrt(PriKine.W2)-0.938272013)
	const char * exp_p_kin_nocorr = "( -(((sqrt(PriKine.W2)-0.938272013)*(sqrt(PriKine.W2)-0.938272013) + 2*0.938272013*(sqrt(PriKine.W2)-0.938272013) - 2*0.938272013*1.2306)/(2*(0.938272013 + 1.2306 - 1.2306*cos(16/180.*3.14159)))) )";
	
	////////////////////////////////////////////////////////////////////////////////
	TCanvas * c1 = new TCanvas(CutUsed+"_Canvas",CutUsed,1900,1000);
	c1->Divide(3,3);
	UInt_t idx=1;


	c1->cd(idx++);
	T->Draw("L.tr.tg_dp>>dP(200,-.05,.05)","","");
	dP->SetTitle("LHRS Momentum Diff");
	dP->SetXTitle("Target dp (GeV)");

// 	c1->cd(idx++);
// 	T->Draw("L.tr.p>>P(200,1.11269,1.31269)","","");
// 	P->SetTitle("LHRS Momentum");
// 	P->SetXTitle("p (GeV)");

	c1->cd(idx++);
	T->Draw(Form("sqrt(PriKine.W2)-0.951636*(%s+%f)>>W(200,0.92,.96)",RadCorrectionScatMom.Data(),NMRMomCorr),"","");
	W->SetTitle("Inv Mass");
	W->SetXTitle("Invariant Mass (GeV)");
// 	W->SetYTitle("Target Theta L.tr.tg_th");
	TLine *l = new TLine(0.938272013,-1,0.938272013,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("%s>>PKin(200,1.16,1.18)",exp_p_kin_nocorr),"","");
	PKin->SetTitle("Simplified LHRS Momentum (Angular Independant)");
	PKin->SetXTitle("p_kin (GeV)");
// 	W->SetYTitle("Target Theta L.tr.tg_th");

	c1->cd(idx++);
	T->Draw(Form("%s/1.21269-1>>dpkin(800,-0.04,-.03)",exp_p_kin_nocorr),"","");
	dpkin->SetTitle("dp_kin (Angular Independant)");
	dpkin->SetXTitle("dp_kin");
// 	W->SetYTitle("Target Theta L.tr.tg_th");

	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();
	
	c1->cd(idx++);
	T->Draw(Form("ReactPt_L.z:(%s+%s)/1.21269-1>>dPvz(800,-0.04,-.03,800,-.3,.3)",exp_p_kin_nocorr,RadCorrectionScatMom.Data()),"","COLZ");
	dPvz->SetTitle("React Z VS dp_kin");
	dPvz->SetXTitle("dp_kin");
	dPvz->SetYTitle("React Z (m)");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_th:(%s+%s)/1.21269-1>>dPth(800,-0.04,-.03,800,-.07,.07)",exp_p_kin_nocorr,RadCorrectionScatMom.Data()),"","COLZ");
	dPth->SetTitle("Target Theta VS dp_kin");
	dPth->SetXTitle("dp_kin");
	dPth->SetYTitle("Target Theta");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_ph:(%s+%s)/1.21269-1>>dPph(800,-0.04,-.03,800,-.045,.045)",exp_p_kin_nocorr,RadCorrectionScatMom.Data()),"","COLZ");
	dPph->SetTitle("Target Phi VS dp_kin");
	dPph->SetXTitle("dp_kin");
	dPph->SetYTitle("Target Phi");

	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();
	
	c1->cd(idx++);
	T->Draw(Form("rb.x:(%s+%s)/1.21269-1>>dPx(800,-0.04,-.03,800,-.004,.004)",exp_p_kin_nocorr,RadCorrectionScatMom.Data()),"","COLZ");
	dPx->SetTitle("rb x VS dp_kin");
	dPx->SetXTitle("dp_kin");
	dPx->SetYTitle("rb x (m)");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();
	
	c1->cd(idx++);
	T->Draw(Form("rb.y:(%s+%s)/1.21269-1>>dPy(800,-0.04,-.03,800,-.0,.008)",exp_p_kin_nocorr,RadCorrectionScatMom.Data()),"","COLZ");
	dPy->SetTitle("rb y VS dp_kin");
	dPy->SetXTitle("dp_kin");
	dPy->SetYTitle("VDC y (m)");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->Print(Form("H2Check.%s.2620.png",CutUsed.Data()),"png");
	////////////////////////////////////////////////////////////////////////////////
// 	TCanvas * c1 = new TCanvas(CutUsed+"_Canvas_RadCorrection",CutUsed,1900,1000);
// 	c1->Divide(1,2);
// 	UInt_t idx=1;
	//
	//
	//
// 	c1->cd(idx++);
// 	T->Draw(Form("sqrt(PriKine.W2)+%s/0.0603824>>W(200,0.92,.96)"
// 			,RadCorrectionScatMom.Data()),"","");
// 	W->SetTitle("Inv Mass");
// 	W->SetXTitle("Invariant Mass (GeV)");
// // 	W->SetYTitle("Target Theta L.tr.tg_th");
	//
// 	c1->cd(idx++);
// 	T->Draw(Form("(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))) +%s)/1.21269-1>>dpkin(800,-0.04,-.03)",RadCorrectionScatMom.Data()),"","");
// 	dpkin->SetTitle("dp_kin (Angular Independant)");
// 	dpkin->SetXTitle("dp_kin");
// // 	W->SetYTitle("Target Theta L.tr.tg_th");
	//
// 	TLine *l = new TLine(RealDpKin,0,RealDpKin,1000000);
// 	l->SetLineColor(6);
// 	l->Draw();
	//

// 	c1->cd(idx++);
// 	T->Draw("ReactPt_L.z:(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))))/1.21269-1>>dPvz(800,-0.04,-.03,800,-.3,.3)","","COLZ");
// 	dPvz->SetTitle("React Z VS dp_kin");
// 	dPvz->SetXTitle("dp_kin");
// 	dPvz->SetYTitle("React Z (m)");
	//
// 	c1->cd(idx++);
// 	T->Draw("L.tr.tg_th:(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))))/1.21269-1>>dPth(800,-0.04,-.03,800,-.07,.07)","","COLZ");
// 	dPth->SetTitle("Target Theta VS dp_kin");
// 	dPth->SetXTitle("dp_kin");
// 	dPth->SetYTitle("Target Theta");
	//
// 	c1->cd(idx++);
// 	T->Draw("L.tr.tg_ph:(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))))/1.21269-1>>dPph(800,-0.04,-.03,800,-.045,.045)","","COLZ");
// 	dPph->SetTitle("Target Phi VS dp_kin");
// 	dPph->SetXTitle("dp_kin");
// 	dPph->SetYTitle("Target Phi");
	//
// 	c1->cd(idx++);
// 	T->Draw("L.tr.x:(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))))/1.21269-1>>dPx(800,-0.04,-.03,800,-.8,.8)","","COLZ");
// 	dPx->SetTitle("VDC x VS dp_kin");
// 	dPx->SetXTitle("dp_kin");
// 	dPx->SetYTitle("VDC x (m)");

// 	c1->Print(Form("H2Check.%s.AfterCorrection.png",CutUsed.Data()),"png");
	
}

He3_RasterCheck_WithDp()
{

	TChain * t=new TChain("T");
	RootFileLoader(t,"./ROOTfiles/e06010_det_L_%d.root",2765,2768);
	
	T->Draw(">>TestCut","abs(ReactPt_L.z)<.1 && abs(L.tr.tg_th)<.062/7 && abs(L.tr.tg_ph)<.035/10 && L.tr.n==1");
	
	TString CutUsed = "TestCut";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);
	
	TCanvas * c1 = new TCanvas("RasterCheck_He3Dp",CutUsed,1900,600);
	c1->Divide(3,1);
	UInt_t idx=1;

	c1->cd(idx++);
	T->Draw("L.tr.tg_dp>>tg_dp(100,-.01,0)");
	
	c1->cd(idx++);
	T->Draw("rb.y:L.tr.tg_dp>>tg_dpvsy(1000,-.01,0,1000,0,.008)");
	
// 	c1->cd(idx++);
// 	T->Draw("rb.y:ExTgtCor_L.dp>>ex_dpvsy(1000,-.01,0,1000,0,.008)");

	
	c1->cd(idx++);
	
	c1->cd(idx++);
	T->Draw("rb.x:L.tr.tg_dp>>tg_dpvsx(1000,-.01,0,1000,-0.004,.004)");
	
// 	c1->cd(idx++);
// 	T->Draw("rb.x:ExTgtCor_L.dp>>ex_dpvsx(1000,-.01,0,1000,-0.004,.004)");
	
	c1->Print(TString(c1->GetName())+TString(".png"),"png");
}
	
He3_RasterCheck_WithDp_2580()
{

	TChain * t=new TChain("T");
	RootFileLoader(t,"./ROOTfiles/e06010_det_L_%d.root",2580,2580);
	
	T->Draw(">>TestCut","abs(ReactPt_L.z)<.02 && abs(L.tr.tg_th)<.062/4 && abs(L.tr.tg_ph)<.035/8 && L.tr.n==1");
	
	TString CutUsed = "TestCut";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);
	
	TCanvas * c1 = new TCanvas("RasterCheck_He3Dp_2580",CutUsed,1900,600);
	c1->Divide(3,1);
	UInt_t idx=1;

	c1->cd(idx++);
	T->Draw("L.tr.tg_dp>>tg_dp(100,-.01,0.01)");
	
	c1->cd(idx++);
	T->Draw("rb.y:L.tr.tg_dp>>tg_dpvsy(1000,-.01,0.01,1000,0,.008)");
	
// 	c1->cd(idx++);
// 	T->Draw("rb.y:ExTgtCor_L.dp>>ex_dpvsy(1000,-.01,0.01,1000,0,.008)");

	
	c1->cd(idx++);
	
	c1->cd(idx++);
	T->Draw("rb.x:L.tr.tg_dp>>tg_dpvsx(1000,-.01,0.01,1000,-0.004,.004)");
	
// 	c1->cd(idx++);
// 	T->Draw("rb.x:ExTgtCor_L.dp>>ex_dpvsx(1000,-.01,0.01,1000,-0.004,.004)");
	
	c1->Print(TString(c1->GetName())+TString(".png"),"png");
}

const Double_t RefCellWall = 1.;

DrawHe3_SelfCalculation()
{
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2766.root");
// 	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2663.root");
// 	TChain * t=new TChain("T");
// 	RootFileLoader(t,"./ROOTfiles/e06010_det_L_%d.root",2765,2768);
	
	T->Draw(">>LooseCut","abs(ReactPt_L.z)<.2 && abs(L.tr.tg_th)<.07 && abs(L.tr.tg_ph)<.04 && abs(L.tr.tg_dp)<.05 && L.tr.n==1");
	T->Draw(">>DeepCut","abs(ReactPt_L.z)<.05 && abs(L.tr.tg_th)<.05 && abs(L.tr.tg_ph)<.003 && L.tr.n==1");
	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1");
	T->Draw(">>TestCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.0062 && abs(L.tr.tg_ph)<.0035 && L.tr.n==1");
	T->Draw(">>TestCut2","abs(ReactPt_L.z)<.04 && abs(L.tr.tg_th)<.0062 && abs(L.tr.tg_ph)<.0035 && L.tr.n==1");


// 	TString CutUsed = "LooseCut";
	TString CutUsed = "ProdCut";
// 	TString CutUsed = "DeepCut";
// 	TString CutUsed = "TestCut";
// 	TString CutUsed = "TestCut2";
	TString DrawOption  = "COLZ";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);

	////////////////////////////////////////////////////////////////////////////////
	
	const Double_t RealMom = DIPOLE_MAG2MOM(0.448926);
	cout<<"NMR Momentum  = "<<RealMom<<endl;
// 	const Double_t NMRMomCorr = 0;

	const Double_t BeamEnergyLossOnDp
			= 0.966718*1e-3
			*(
			  0.0895//Beam window
			+4.7435E-04*(51.6-20)//He4 Before Cell
			+0.4712*0.142//Cell Entrance
			+4.7598E-04*6.10344729*20//Half Cell Length
			 );
	const Double_t BeamEnergyLossOnDp_ReactZDep = 0.966718*1e-3*4.7598E-04*6.10344729*100 ;//GeV/m
	const Double_t AfterScatteringEnergyLoss
			= 1e-3*(
					4.7598E-04*6.10344729*1/TMath::Sin(16/180.*3.14159)	//in cell He3
			+0.4712*RefCellWall/TMath::Sin(16/180.*3.14159)	//Exit from cell wall
			+4.7435E-04*70.1548		//He4
			+.1//Target Enclosure
			+3.2433E-03*(116.078-70.1548)	//Air
			+0.0798//HRS Entrance
				   );
	TString RadCorrectionScatMom
			= Form("(%f+%f*ReactPt_L.z)"
			,BeamEnergyLossOnDp+AfterScatteringEnergyLoss
			,BeamEnergyLossOnDp_ReactZDep
				  );
	
	cout<<"Total Correction on p = "<<RadCorrectionScatMom<<endl;
	
	const Double_t RealDpKin = ( 2*(3.0160293*.931494028-.000511*2)*1.23046)/(2*((3.0160293*.931494028-.000511*2) + 1.23046 - 1.23046*TMath::Cos(16/180.*3.14159)))/RealMom-1;

	//DM = (sqrt(PriKine_He3.W2)-(3.0160293*.931494028-.000511*2))
	TString exp_dp_kin = Form("(ExTgtCor_L.dp - ( (( 2*2.80839*rb.e)/(2*(2.80839 + rb.e - (rb.e*cos(L.tr.tg_ph + 16/180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/%f-1-(%f) ) + %s)",RealMom,RealDpKin,RadCorrectionScatMom.Data());

	cout<<"exp_dp_kin = "<<exp_dp_kin<<endl;
	
	//////////////////////////////////////////////////////////////////////////////
	TCanvas * c1 = new TCanvas(CutUsed+"_DependencyCheck_WithRadCorrection_SelfCalc",CutUsed,1900,1000);
	c1->Divide(4,3);
	UInt_t idx=1;


// 	c1->cd(idx++);
// 	T->Draw("L.tr.tg_dp>>dP(200,-.05,.05)","","");
// 	dP->SetTitle("LHRS Momentum Diff");
// 	dP->SetXTitle("Target dp (GeV)");


	c1->cd(idx++);
	T->Draw(Form("%s>>dpkin1(800,-0.05,.05)",exp_dp_kin.Data()),"","");
	dpkin1->SetTitle("dp_kin (Angular Independant)");
	dpkin1->SetXTitle("dp_kin");

	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("ReactPt_L.z:%s>>dPvz(800,-0.01,.00,800,-.3,.3)",exp_dp_kin.Data()),"",DrawOption);
	dPvz->SetTitle("React Z VS dp_kin");
	dPvz->SetXTitle("dp_kin");
	dPvz->SetYTitle("React Z (m)");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("rb.x:%s>>dPx(800,-0.01,.00,800,-.004,.004)",exp_dp_kin.Data()),"",DrawOption);
	dPx->SetTitle("rb x VS dp_kin");
	dPx->SetXTitle("dp_kin");
	dPx->SetYTitle("rb x (m)");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("rb.y:%s>>dPy(800,-0.01,.00,800,-.000,.008)",exp_dp_kin.Data()),"",DrawOption);
	dPy->SetTitle("rb y VS dp_kin");
	dPy->SetXTitle("dp_kin");
	dPy->SetYTitle("rb y (m)");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_th:%s>>dPth(800,-0.01,.00,800,-.07,.07)",exp_dp_kin.Data()),"",DrawOption);
	dPth->SetTitle("Target Theta VS dp_kin");
	dPth->SetXTitle("dp_kin");
	dPth->SetYTitle("Target Theta");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_ph:%s>>dPph(800,-0.01,.00,800,-.045,.045)",exp_dp_kin.Data()),"",DrawOption);
	dPph->SetTitle("Target Phi VS dp_kin");
	dPph->SetXTitle("dp_kin");
	dPph->SetYTitle("Target Phi");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();


	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_y:%s>>dPtgy(800,-0.01,.00,800,-.06,.06)",exp_dp_kin.Data()),"",DrawOption);
	dPtgy->SetTitle("Target Y VS dp_kin");
	dPtgy->SetXTitle("dp_kin");
	dPtgy->SetYTitle("Target Y");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("ExTgtCor_L.x:%s>>dPtgx(800,-0.01,.00,800,-.01,.01)",exp_dp_kin.Data()),"",DrawOption);
	dPtgx->SetTitle("Target X VS dp_kin");
	dPtgx->SetXTitle("dp_kin");
	dPtgx->SetYTitle("Target X");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.r_x:%s>>dPr_x(800,-0.01,.00,800,-.2,.2)",exp_dp_kin.Data()),"",DrawOption);
	dPr_x->SetTitle("Rotated X VS dp_kin");
	dPr_x->SetXTitle("dp_kin");
	dPr_x->SetYTitle("Rotated X");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.r_th:%s>>dPr_th(800,-0.01,.00,800,-.03,.03)",exp_dp_kin.Data()),"",DrawOption);
	dPr_th->SetTitle("Rotated Th VS dp_kin");
	dPr_th->SetXTitle("dp_kin");
	dPr_th->SetYTitle("Rotated Th");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.r_y:%s>>dPr_y(800,-0.01,.00,800,-.06,.06)",exp_dp_kin.Data()),"",DrawOption);
	dPr_y->SetTitle("Rotated Y VS dp_kin");
	dPr_y->SetXTitle("dp_kin");
	dPr_y->SetYTitle("Rotated Y");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.r_ph:%s>>dPr_ph(800,-0.01,.00,800,-.06,.06)",exp_dp_kin.Data()),"",DrawOption);
	dPr_ph->SetTitle("Rotated Ph VS dp_kin");
	dPr_ph->SetXTitle("dp_kin");
	dPr_ph->SetYTitle("Rotated Ph");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();


	c1->Print(TString(c1->GetName())+TString(".png"),"png");
	
	
	////////////////////////////////////////////////////////////////////////////////
	
	TCanvas * c1 = new TCanvas(CutUsed+"_WithRadCorrection_SelfCalc",CutUsed,1900,1000);
	c1->Divide(2,2);
	UInt_t idx=1;



	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_dp>>dp(800,-0.05,.05)"),"","");
// 	dpkin->SetTitle("dp");
	dp->SetXTitle("dp");
	
	c1->cd(idx++);
	T->Draw(Form("ExTgtCor_L.dp>>extdp(800,-0.05,.05)",exp_dp_kin.Data()),"","");
// 	dpkin->SetTitle("dp_kin (Angular Independant)");
	extdp->SetXTitle("ExTgtCor_L.dp");

	c1->cd(idx++);
	T->Draw(Form("%s>>dpkin(800,-0.05,.05)",exp_dp_kin.Data()),"","");
	dpkin->SetTitle("dp_kin (Angular Independant)");
	dpkin->SetXTitle("dp_kin");

	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("%s>>dpkinz(800,-0.01,.01)",exp_dp_kin.Data()),"","");
	dpkinz->SetTitle("dp_kin Zoomed (Angular Independant)");
	dpkinz->SetXTitle("dp_kin");

	const Double_t VertexRes = 6e-4;
	const Double_t VertexResN = 2;
	TF1 *f = new TF1("gausl2","gaus+[3]+[4]*x"
			,RealDpKin-VertexRes*VertexResN,RealDpKin+VertexRes*VertexResN);
	f->SetParameter(1,RealDpKin);
	f->SetParameter(2,VertexRes);

	dpkinz->Fit("gausl2","R");
	
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	TLatex *t = new TLatex(f->GetParameter(1)+f->GetParameter(2)
			,f->GetParameter(0)+f->GetParameter(3)+f->GetParameter(4)*f->GetParameter(1)
			,Form("\\Delta \\pm \\sigma = %2.1f \\pm %2.1f \\times 10^{-4}"
					,10000*(f->GetParameter(1)-RealDpKin)
							,10000*f->GetParameter(2)));
	t->SetTextSize(0.05);
	t->SetTextAlign(12);
	t->SetTextColor(2);
	t->Draw();
	
	c1->Print(TString(c1->GetName())+TString(".png"),"png");
}

DrawH2_SelfCalculation()
{
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2620.root");
// 	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2663.root");
// 	TChain * t=new TChain("T");
// 	RootFileLoader(t,"./ROOTfiles/e06010_det_L_%d.root",2765,2768);
	
// 	T->Draw(">>LooseCut","abs(ReactPt_L.z)<.2 && abs(L.tr.tg_th)<.07 && abs(L.tr.tg_ph)<.04 && abs(L.tr.tg_dp)<.05 && L.tr.n==1");
// 	T->Draw(">>DeepCut","abs(ReactPt_L.z)<.05 && abs(L.tr.tg_th)<.05 && abs(L.tr.tg_ph)<.003 && L.tr.n==1");
	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1 && rb.y<0.0045");
// 	T->Draw(">>TestCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.0062 && abs(L.tr.tg_ph)<.0035 && L.tr.n==1");
// 	T->Draw(">>TestCut2","abs(ReactPt_L.z)<.04 && abs(L.tr.tg_th)<.0062 && abs(L.tr.tg_ph)<.0035 && L.tr.n==1");


// 	TString CutUsed = "LooseCut";
	TString CutUsed = "ProdCut";
// 	TString CutUsed = "DeepCut";
// 	TString CutUsed = "TestCut";
// 	TString CutUsed = "TestCut2";
	TString DrawOption  = "COLZ";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);

	////////////////////////////////////////////////////////////////////////////////
	
	const Double_t RealMom = DIPOLE_MAG2MOM(0.448924);
	cout<<"NMR Momentum  = "<<RealMom<<endl;
// 	const Double_t NMRMomCorr = 0;

	const Double_t BeamEnergyLossOnDp
			= 0.90561 *1e-3
			*(
			  0.0895//Beam window
			+4.7435E-04*(51.6-20)//He4 Before Cell
			+0.4712*0.142//Cell Entrance
			+4.7426E-04*7.80459639*20//Half Cell Length
			 );
	const Double_t BeamEnergyLossOnDp_ReactZDep = 0.90561 *1e-3*4.7426E-04*7.80459639*100 ;//GeV/m
	const Double_t AfterScatteringEnergyLoss
			= 1e-3*(
					4.7426E-04*7.80459639*1/TMath::Sin(16/180.*3.14159)	//in cell He3
			+0.4712*RefCellWall/TMath::Sin(16/180.*3.14159)	//Exit from cell wall
			+4.7435E-04*70.1548		//He4
			+.1//Target Enclosure
			+3.2433E-03*(116.078-70.1548)	//Air
			+0.0798//HRS Entrance
				   );
	TString RadCorrectionScatMom
			= Form("(%f+%f*ReactPt_L.z)"
			,BeamEnergyLossOnDp+AfterScatteringEnergyLoss
			,BeamEnergyLossOnDp_ReactZDep
				  );
	
	cout<<"Total Correction on p = "<<RadCorrectionScatMom<<endl;
	
	const Double_t RealDpKin = ( 2*(0.938272013)*1.2306)/(2*((0.938272013) + 1.2306 - 1.2306*TMath::Cos(16/180.*3.14159)))/RealMom-1;

	TString exp_dp_kin = Form("(ExTgtCor_L.dp- ( (( 2*.938272013*rb.e)/(2*(.938272013 + rb.e - (rb.e*cos(L.tr.tg_ph + 16./180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/%f-1-(%f) ) + %s)",RealMom,RealDpKin,RadCorrectionScatMom.Data());

	cout<<"exp_dp_kin = "<<exp_dp_kin<<endl;
	
	//////////////////////////////////////////////////////////////////////////////
	TCanvas * c1 = new TCanvas(TString("H2Ref.")+CutUsed+"_DependencyCheck_WithRadCorrection_SelfCalc",CutUsed,1900,1000);
	c1->Divide(4,3);
	UInt_t idx=1;


// 	c1->cd(idx++);
// 	T->Draw("L.tr.tg_dp>>dP(200,-.05,.05)","","");
// 	dP->SetTitle("LHRS Momentum Diff");
// 	dP->SetXTitle("Target dp (GeV)");


	c1->cd(idx++);
	T->Draw(Form("%s>>dpkin1(800,-0.05,.05)",exp_dp_kin.Data()),"","");
	dpkin1->SetTitle("dp_kin (Angular Independant)");
	dpkin1->SetXTitle("dp_kin");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("ReactPt_L.z:%s>>dPvz(800,-0.05,-.02,800,-.3,.3)",exp_dp_kin.Data()),"",DrawOption);
	dPvz->SetTitle("React Z VS dp_kin");
	dPvz->SetXTitle("dp_kin");
	dPvz->SetYTitle("React Z (m)");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("rb.x:%s>>dPx(800,-0.05,-.02,800,-.004,.004)",exp_dp_kin.Data()),"",DrawOption);
	dPx->SetTitle("rb x VS dp_kin");
	dPx->SetXTitle("dp_kin");
	dPx->SetYTitle("rb x (m)");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("rb.y:%s>>dPy(800,-0.05,-.02,800,-.000,.008)",exp_dp_kin.Data()),"",DrawOption);
	dPy->SetTitle("rb y VS dp_kin");
	dPy->SetXTitle("dp_kin");
	dPy->SetYTitle("rb y (m)");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_th:%s>>dPth(800,-0.05,-.02,800,-.07,.07)",exp_dp_kin.Data()),"",DrawOption);
	dPth->SetTitle("Target Theta VS dp_kin");
	dPth->SetXTitle("dp_kin");
	dPth->SetYTitle("Target Theta");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_ph:%s>>dPph(800,-0.05,-.02,800,-.045,.045)",exp_dp_kin.Data()),"",DrawOption);
	dPph->SetTitle("Target Phi VS dp_kin");
	dPph->SetXTitle("dp_kin");
	dPph->SetYTitle("Target Phi");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();


	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_y:%s>>dPtgy(800,-0.05,-.02,800,-.06,.06)",exp_dp_kin.Data()),"",DrawOption);
	dPtgy->SetTitle("Target Y VS dp_kin");
	dPtgy->SetXTitle("dp_kin");
	dPtgy->SetYTitle("Target Y");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("ExTgtCor_L.x:%s>>dPtgx(800,-0.05,-.02,800,-.01,.01)",exp_dp_kin.Data()),"",DrawOption);
	dPtgx->SetTitle("Target X VS dp_kin");
	dPtgx->SetXTitle("dp_kin");
	dPtgx->SetYTitle("Target X");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.r_x:%s>>dPr_x(800,-0.05,-.02,800,-.8,-.2)",exp_dp_kin.Data()),"",DrawOption);
	dPr_x->SetTitle("Rotated X VS dp_kin");
	dPr_x->SetXTitle("dp_kin");
	dPr_x->SetYTitle("Rotated X");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.r_th:%s>>dPr_th(800,-0.05,-.02,800,-.03,.03)",exp_dp_kin.Data()),"",DrawOption);
	dPr_th->SetTitle("Rotated Th VS dp_kin");
	dPr_th->SetXTitle("dp_kin");
	dPr_th->SetYTitle("Rotated Th");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.r_y:%s>>dPr_y(800,-0.05,-.02,800,-.06,.06)",exp_dp_kin.Data()),"",DrawOption);
	dPr_y->SetTitle("Rotated Y VS dp_kin");
	dPr_y->SetXTitle("dp_kin");
	dPr_y->SetYTitle("Rotated Y");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.r_ph:%s>>dPr_ph(800,-0.05,-.02,800,-.06,.06)",exp_dp_kin.Data()),"",DrawOption);
	dPr_ph->SetTitle("Rotated Ph VS dp_kin");
	dPr_ph->SetXTitle("dp_kin");
	dPr_ph->SetYTitle("Rotated Ph");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();


	c1->Print(TString(c1->GetName())+TString(".png"),"png");

	
	////////////////////////////////////////////////////////////////////////////////
	
	TCanvas * c1 = new TCanvas(TString("H2Ref.")+CutUsed+"_WithRadCorrection_SelfCalc",CutUsed,1900,1000);
	c1->Divide(2,2);
	UInt_t idx=1;



	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_dp>>dp(800,-0.05,.05)"),"","");
// 	dpkin->SetTitle("dp");
	dp->SetXTitle("dp");
	
	c1->cd(idx++);
	T->Draw(Form("ExTgtCor_L.dp>>extdp(800,-0.05,.05)",exp_dp_kin.Data()),"","");
// 	dpkin->SetTitle("dp_kin (Angular Independant)");
	extdp->SetXTitle("ExTgtCor_L.dp");

	c1->cd(idx++);
	T->Draw(Form("%s>>dpkin(800,-0.05,.05)",exp_dp_kin.Data()),"","");
	dpkin->SetTitle("dp_kin (Angular Independant)");
	dpkin->SetXTitle("dp_kin");

	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("%s>>dpkinz(800,-0.05,.02)",exp_dp_kin.Data()),"","");
	dpkinz->SetTitle("dp_kin Zoomed (Angular Independant)");
	dpkinz->SetXTitle("dp_kin");

	const Double_t VertexRes = 12e-4;
	const Double_t VertexResN = 2;
	TF1 *f = new TF1("gausl2","gaus+[3]+[4]*x"
			,RealDpKin-VertexRes*VertexResN,RealDpKin+VertexRes*VertexResN);
	f->SetParameter(1,RealDpKin);
	f->SetParameter(2,VertexRes);

	dpkinz->Fit("gausl2","R");

	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	TLatex *t = new TLatex(f->GetParameter(1)+f->GetParameter(2)
			,f->GetParameter(0)+f->GetParameter(3)+f->GetParameter(4)*f->GetParameter(1)
			,Form("\\Delta \\pm \\sigma = %2.1f \\pm %2.1f \\times 10^{-4}"
					,10000*(f->GetParameter(1)-RealDpKin)
							,10000*f->GetParameter(2)));
	t->SetTextSize(0.05);
	t->SetTextAlign(12);
	t->SetTextColor(2);
	t->Draw();
	
	c1->Print(TString(c1->GetName())+TString(".png"),"png");
}

DrawC_SelfCalculation()
{
// 	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2505.root");
// 	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2647.root");
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2580.root");
// 	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2705.root");
	
// 	TChain * t=new TChain("T");
// 	RootFileLoader(t,"./ROOTfiles/e06010_det_L_%d.root",2765,2768);
	
// 	T->Draw(">>LooseCut","abs(ReactPt_L.z)<.2 && abs(L.tr.tg_th)<.07 && abs(L.tr.tg_ph)<.04 && abs(L.tr.tg_dp)<.05 && L.tr.n==1");
// 	T->Draw(">>DeepCut","abs(ReactPt_L.z)<.05 && abs(L.tr.tg_th)<.05 && abs(L.tr.tg_ph)<.003 && L.tr.n==1");
// 	T->Draw(">>TestCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.0062 && abs(L.tr.tg_ph)<.0035 && L.tr.n==1");
// 	T->Draw(">>TestCut2","abs(ReactPt_L.z)<.04 && abs(L.tr.tg_th)<.0062 && abs(L.tr.tg_ph)<.0035 && L.tr.n==1");

// 	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.23 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1");
	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.23/3 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035/10 && L.tr.n==1 && abs(ReactPt_L.z)<.02");

// 	TString CutUsed = "LooseCut";
	TString CutUsed = "ProdCut";
// 	TString CutUsed = "DeepCut";
// 	TString CutUsed = "TestCut";
// 	TString CutUsed = "TestCut2";
	TString DrawOption  = "COLZ";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);

	////////////////////////////////////////////////////////////////////////////////
	
	const Double_t RealMom = DIPOLE_MAG2MOM(0.448922);
	cout<<"NMR Momentum  = "<<RealMom<<endl;
// 	const Double_t NMRMomCorr = 0;

	const Double_t BeamEnergyLossOnDp
			= 0.99152 *1e-3
			*(
			  0.0895//Beam window
			+0.4299//BeO
			+4.7435E-04*51.6//he4
			+4*0.1072//mid carbon foil
			 );
	const Double_t BeamEnergyLossOnDp_ReactZDep = 0.99152 *1e-3*0.1072/(.2/3) ;//GeV/m
	const Double_t AfterScatteringEnergyLoss
			= 1e-3*(
					+4.7435E-04*70.1548		//He4
			+.1//Target Enclosure
			+3.2433E-03*(116.078-70.1548)	//Air
			+0.0798//HRS Entrance
			-.180*1.2//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
				   );
	TString RadCorrectionScatMom
			= Form("(%f+%f*ReactPt_L.z)"
			,BeamEnergyLossOnDp+AfterScatteringEnergyLoss
			,BeamEnergyLossOnDp_ReactZDep
				  );


	const UInt_t NExcitationStates = 4; //Carbon Excitation States
	const Double_t ExcitationEnergyList[]
			={0,0.00443891,0.00765420,0.009641,0.010844,0.011160,0.011828,0.012710,0.013352,0.014083,0.015110,0.015440,0.016105,0.016570,0.017230};

	
	cout<<"Total Correction on p = "<<RadCorrectionScatMom<<endl;
	
	const Double_t RealDpKin = ( 2*(12*.931494028-.000511*6)*1.2306)/(2*((12*.931494028-.000511*6) + 1.2306 - 1.2306*TMath::Cos(16/180.*3.14159)))/RealMom-1;

	TString exp_dp_kin = Form("(ExTgtCor_L.dp- ( (( 2*(12*.931494028-.000511*6)*rb.e)/(2*((12*.931494028-.000511*6) + rb.e - (rb.e*cos(L.tr.tg_ph + 16./180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/%f-1-(%f) ) + %s)",RealMom,RealDpKin,RadCorrectionScatMom.Data());

	cout<<"exp_dp_kin = "<<exp_dp_kin<<endl;
	
	//////////////////////////////////////////////////////////////////////////////
	TCanvas * c1 = new TCanvas(TString("Carbon.")+CutUsed+"_DependencyCheck_WithRadCorrection_SelfCalc",CutUsed,1900,1000);
	c1->Divide(4,3);
	UInt_t idx=1;


// 	c1->cd(idx++);
// 	T->Draw("L.tr.tg_dp>>dP(200,-.05,.05)","","");
// 	dP->SetTitle("LHRS Momentum Diff");
// 	dP->SetXTitle("Target dp (GeV)");


	c1->cd(idx++);
	T->Draw(Form("%s>>dpkin1(800,-0.05,.05)",exp_dp_kin.Data()),"","");
	dpkin1->SetTitle("dp_kin (Angular Independant)");
	dpkin1->SetXTitle("dp_kin");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("ReactPt_L.z:%s>>dPvz(800,0,.015,800,-.3,.3)",exp_dp_kin.Data()),"",DrawOption);
	dPvz->SetTitle("React Z VS dp_kin");
	dPvz->SetXTitle("dp_kin");
	dPvz->SetYTitle("React Z (m)");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("rb.x:%s>>dPx(800,0,.015,800,-.004,.004)",exp_dp_kin.Data()),"",DrawOption);
	dPx->SetTitle("rb x VS dp_kin");
	dPx->SetXTitle("dp_kin");
	dPx->SetYTitle("rb x (m)");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("rb.y:%s>>dPy(800,0,.015,800,-.000,.008)",exp_dp_kin.Data()),"",DrawOption);
	dPy->SetTitle("rb y VS dp_kin");
	dPy->SetXTitle("dp_kin");
	dPy->SetYTitle("rb y (m)");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_th:%s>>dPth(800,0,.015,800,-.07,.07)",exp_dp_kin.Data()),"",DrawOption);
	dPth->SetTitle("Target Theta VS dp_kin");
	dPth->SetXTitle("dp_kin");
	dPth->SetYTitle("Target Theta");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_ph:%s>>dPph(800,0,.015,800,-.045,.045)",exp_dp_kin.Data()),"",DrawOption);
	dPph->SetTitle("Target Phi VS dp_kin");
	dPph->SetXTitle("dp_kin");
	dPph->SetYTitle("Target Phi");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}


	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_y:%s>>dPtgy(800,0,.015,800,-.06,.06)",exp_dp_kin.Data()),"",DrawOption);
	dPtgy->SetTitle("Target Y VS dp_kin");
	dPtgy->SetXTitle("dp_kin");
	dPtgy->SetYTitle("Target Y");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("ExTgtCor_L.x:%s>>dPtgx(800,0,.015,800,-.01,.01)",exp_dp_kin.Data()),"",DrawOption);
	dPtgx->SetTitle("Target X VS dp_kin");
	dPtgx->SetXTitle("dp_kin");
	dPtgx->SetYTitle("Target X");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("L.tr.r_x:%s>>dPr_x(800,0,.015,800,-.8,.8)",exp_dp_kin.Data()),"",DrawOption);
	dPr_x->SetTitle("Rotated X VS dp_kin");
	dPr_x->SetXTitle("dp_kin");
	dPr_x->SetYTitle("Rotated X");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}
	c1->cd(idx++);
	T->Draw(Form("L.tr.r_th:%s>>dPr_th(800,0,.015,800,-.03,.03)",exp_dp_kin.Data()),"",DrawOption);
	dPr_th->SetTitle("Rotated Th VS dp_kin");
	dPr_th->SetXTitle("dp_kin");
	dPr_th->SetYTitle("Rotated Th");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}
	c1->cd(idx++);
	T->Draw(Form("L.tr.r_y:%s>>dPr_y(800,0,.015,800,-.02,.02)",exp_dp_kin.Data()),"",DrawOption);
	dPr_y->SetTitle("Rotated Y VS dp_kin");
	dPr_y->SetXTitle("dp_kin");
	dPr_y->SetYTitle("Rotated Y");
	
	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("L.tr.r_ph:%s>>dPr_ph(800,0,.015,800,-.06,.06)",exp_dp_kin.Data()),"",DrawOption);
	dPr_ph->SetTitle("Rotated Ph VS dp_kin");
	dPr_ph->SetXTitle("dp_kin");
	dPr_ph->SetYTitle("Rotated Ph");
	
	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}


	c1->Print(TString(c1->GetName())+TString(".png"),"png");

	
	////////////////////////////////////////////////////////////////////////////////
	
	TCanvas * c1 = new TCanvas(TString("Carbon.")+CutUsed+"_WithRadCorrection_SelfCalc",CutUsed,1900,1000);
	c1->Divide(2,2);
	UInt_t idx=1;



	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_dp>>dp(800,-0.05,.05)"),"","");
// 	dpkin->SetTitle("dp");
	dp->SetXTitle("dp");
	
	c1->cd(idx++);
	T->Draw(Form("ExTgtCor_L.dp>>extdp(800,-0.05,.05)",exp_dp_kin.Data()),"","");
// 	dpkin->SetTitle("dp_kin (Angular Independant)");
	extdp->SetXTitle("ExTgtCor_L.dp");

	c1->cd(idx++);
	T->Draw(Form("%s>>dpkin(800,-0.05,.05)",exp_dp_kin.Data()),"","");
	dpkin->SetTitle("dp_kin (Angular Independant)");
	dpkin->SetXTitle("dp_kin");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}
	
	c1->cd(idx++);
	T->Draw(Form("%s>>dpkinz(800,-0.0,.015)",exp_dp_kin.Data()),"","");
	dpkinz->SetTitle("dp_kin Zoomed (Angular Independant)");
	dpkinz->SetXTitle("dp_kin");

	const Double_t VertexRes = 6e-4;
	const Double_t VertexResN = 2;
	TF1 *f = new TF1("gausl2","gaus+[3]+[4]*x"
			,RealDpKin-VertexRes*VertexResN,RealDpKin+VertexRes*VertexResN);
	f->SetParameter(1,RealDpKin);
	f->SetParameter(2,VertexRes);

	dpkinz->Fit("gausl2","R");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	TLatex *t = new TLatex(f->GetParameter(1)+f->GetParameter(2)
			,f->GetParameter(0)+f->GetParameter(3)+f->GetParameter(4)*f->GetParameter(1)
			,Form("\\Delta \\pm \\sigma = %2.1f \\pm %2.1f \\times 10^{-4}"
					,10000*(f->GetParameter(1)-RealDpKin)
							,10000*f->GetParameter(2)));
	t->SetTextSize(0.05);
	t->SetTextAlign(12);
	t->SetTextColor(2);
	t->Draw();
	
	c1->Print(TString(c1->GetName())+TString(".png"),"png");
}

DrawC_SelfCalculation2580_dp()
{
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2580.root");
// 	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2663.root");
// 	TChain * t=new TChain("T");
// 	RootFileLoader(t,"./ROOTfiles/e06010_det_L_%d.root",2765,2768);
	
// 	T->Draw(">>LooseCut","abs(ReactPt_L.z)<.2 && abs(L.tr.tg_th)<.07 && abs(L.tr.tg_ph)<.04 && abs(L.tr.tg_dp)<.05 && L.tr.n==1");
// 	T->Draw(">>DeepCut","abs(ReactPt_L.z)<.05 && abs(L.tr.tg_th)<.05 && abs(L.tr.tg_ph)<.003 && L.tr.n==1");
	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.023/3 && abs(L.tr.tg_th)<.062/4 && abs(L.tr.tg_ph)<.035/10 && L.tr.n==1");
// 	T->Draw(">>TestCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.0062 && abs(L.tr.tg_ph)<.0035 && L.tr.n==1");
// 	T->Draw(">>TestCut2","abs(ReactPt_L.z)<.04 && abs(L.tr.tg_th)<.0062 && abs(L.tr.tg_ph)<.0035 && L.tr.n==1");


// 	TString CutUsed = "LooseCut";
	TString CutUsed = "ProdCut";
// 	TString CutUsed = "DeepCut";
// 	TString CutUsed = "TestCut";
// 	TString CutUsed = "TestCut2";
	TString DrawOption  = "SCAT";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);

	////////////////////////////////////////////////////////////////////////////////
	
	const Double_t RealMom = DIPOLE_MAG2MOM(0.448922);
	cout<<"NMR Momentum  = "<<RealMom<<endl;
// 	const Double_t NMRMomCorr = 0;

	const Double_t BeamEnergyLossOnDp
			= 0.99152 *1e-3
			*(
			  0.0895//Beam window
			+0.4299//BeO
			+4.7435E-04*51.6//he4
			+4*0.1072//mid carbon foil
			 );
	const Double_t BeamEnergyLossOnDp_ReactZDep = 0.99152 *1e-3*0.1072/(.2/3) ;//GeV/m
	const Double_t AfterScatteringEnergyLoss
			= 1e-3*(
					+4.7435E-04*70.1548		//He4
			+.1//Target Enclosure
			+3.2433E-03*(116.078-70.1548)	//Air
			+0.0798//HRS Entrance
			-.180*1.2//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
				   );
	TString RadCorrectionScatMom
			= Form("(%f+%f*ReactPt_L.z)"
			,BeamEnergyLossOnDp+AfterScatteringEnergyLoss
			,BeamEnergyLossOnDp_ReactZDep
				  );


	const UInt_t NExcitationStates = 4; //Carbon Excitation States
	const Double_t ExcitationEnergyList[]
			={0,0.00443891,0.00765420,0.009641,0.010844,0.011160,0.011828,0.012710,0.013352,0.014083,0.015110,0.015440,0.016105,0.016570,0.017230};

	
	cout<<"Total Correction on p = "<<RadCorrectionScatMom<<endl;
	
	const Double_t RealDpKin = ( 2*(12*.931494028-.000511*6)*1.2306)/(2*((12*.931494028-.000511*6) + 1.2306 - 1.2306*TMath::Cos(16/180.*3.14159)))/RealMom-1;

	TString exp_dp_kin = Form("(L.tr.tg_dp-rb.y*0.2- ( (( 2*(12*.931494028-.000511*6)*rb.e)/(2*((12*.931494028-.000511*6) + rb.e - (rb.e*cos(L.tr.tg_ph + 16./180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_L.th)^2))))/%f-1-(%f) ) + %s)",RealMom,RealDpKin,RadCorrectionScatMom.Data());

	cout<<"exp_dp_kin = "<<exp_dp_kin<<endl;
	
	//////////////////////////////////////////////////////////////////////////////
	TCanvas * c1 = new TCanvas(TString("Carbon.2580.DpDependency")+CutUsed+"_DependencyCheck_WithRadCorrection_SelfCalc",CutUsed,1900,1000);
	c1->Divide(4,3);
	UInt_t idx=1;


// 	c1->cd(idx++);
// 	T->Draw("L.tr.tg_dp>>dP(200,-.05,.05)","","");
// 	dP->SetTitle("LHRS Momentum Diff");
// 	dP->SetXTitle("Target dp (GeV)");


	c1->cd(idx++);
	T->Draw(Form("%s>>dpkin1(800,-0.05,.05)",exp_dp_kin.Data()),"","");
	dpkin1->SetTitle("dp_kin (Angular Independant)");
	dpkin1->SetXTitle("dp_kin");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("ReactPt_L.z:%s>>dPvz(800,0,.015,800,-.3,.3)",exp_dp_kin.Data()),"",DrawOption);
	dPvz->SetTitle("React Z VS dp_kin");
	dPvz->SetXTitle("dp_kin");
	dPvz->SetYTitle("React Z (m)");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("rb.x:%s>>dPx(800,0,.015,800,-.004,.004)",exp_dp_kin.Data()),"",DrawOption);
	dPx->SetTitle("rb x VS dp_kin");
	dPx->SetXTitle("dp_kin");
	dPx->SetYTitle("rb x (m)");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("rb.y:%s>>dPy(800,0,.015,800,-.000,.008)",exp_dp_kin.Data()),"",DrawOption);
	dPy->SetTitle("rb y VS dp_kin");
	dPy->SetXTitle("dp_kin");
	dPy->SetYTitle("rb y (m)");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_th:%s>>dPth(800,0,.015,800,-.07,.07)",exp_dp_kin.Data()),"",DrawOption);
	dPth->SetTitle("Target Theta VS dp_kin");
	dPth->SetXTitle("dp_kin");
	dPth->SetYTitle("Target Theta");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_ph:%s>>dPph(800,0,.015,800,-.045,.045)",exp_dp_kin.Data()),"",DrawOption);
	dPph->SetTitle("Target Phi VS dp_kin");
	dPph->SetXTitle("dp_kin");
	dPph->SetYTitle("Target Phi");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}


	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_y:%s>>dPtgy(800,0,.015,800,-.06,.06)",exp_dp_kin.Data()),"",DrawOption);
	dPtgy->SetTitle("Target Y VS dp_kin");
	dPtgy->SetXTitle("dp_kin");
	dPtgy->SetYTitle("Target Y");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("ExTgtCor_L.x:%s>>dPtgx(800,0,.015,800,-.01,.01)",exp_dp_kin.Data()),"",DrawOption);
	dPtgx->SetTitle("Target X VS dp_kin");
	dPtgx->SetXTitle("dp_kin");
	dPtgx->SetYTitle("Target X");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("L.tr.r_x:%s>>dPr_x(800,0,.015,800,-.8,.8)",exp_dp_kin.Data()),"",DrawOption);
	dPr_x->SetTitle("Rotated X VS dp_kin");
	dPr_x->SetXTitle("dp_kin");
	dPr_x->SetYTitle("Rotated X");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}
	c1->cd(idx++);
	T->Draw(Form("L.tr.r_th:%s>>dPr_th(800,0,.015,800,-.03,.03)",exp_dp_kin.Data()),"",DrawOption);
	dPr_th->SetTitle("Rotated Th VS dp_kin");
	dPr_th->SetXTitle("dp_kin");
	dPr_th->SetYTitle("Rotated Th");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}
	c1->cd(idx++);
	T->Draw(Form("L.tr.r_y:%s>>dPr_y(800,0,.015,800,-.02,.02)",exp_dp_kin.Data()),"",DrawOption);
	dPr_y->SetTitle("Rotated Y VS dp_kin");
	dPr_y->SetXTitle("dp_kin");
	dPr_y->SetYTitle("Rotated Y");
	
	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("L.tr.r_ph:%s>>dPr_ph(800,0,.015,800,-.06,.06)",exp_dp_kin.Data()),"",DrawOption);
	dPr_ph->SetTitle("Rotated Ph VS dp_kin");
	dPr_ph->SetXTitle("dp_kin");
	dPr_ph->SetYTitle("Rotated Ph");
	
	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}


	c1->Print(TString(c1->GetName())+TString(".png"),"png");

	
	////////////////////////////////////////////////////////////////////////////////
	
	TCanvas * c1 = new TCanvas(TString("Carbon.2580.DpDependency.")+CutUsed+"_WithRadCorrection_SelfCalc",CutUsed,1900,1000);
	c1->Divide(2,2);
	UInt_t idx=1;



	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_dp>>dp(800,-0.05,.05)"),"","");
// 	dpkin->SetTitle("dp");
	dp->SetXTitle("dp");
	
	c1->cd(idx++);
	T->Draw(Form("ExTgtCor_L.dp>>extdp(800,-0.05,.05)",exp_dp_kin.Data()),"","");
// 	dpkin->SetTitle("dp_kin (Angular Independant)");
	extdp->SetXTitle("ExTgtCor_L.dp");

	c1->cd(idx++);
	T->Draw(Form("%s>>dpkin(800,-0.05,.05)",exp_dp_kin.Data()),"","");
	dpkin->SetTitle("dp_kin (Angular Independant)");
	dpkin->SetXTitle("dp_kin");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}
	
	c1->cd(idx++);
	T->Draw(Form("%s>>dpkinz(800,-0.0,.015)",exp_dp_kin.Data()),"","");
	dpkinz->SetTitle("dp_kin Zoomed (Angular Independant)");
	dpkinz->SetXTitle("dp_kin");

	const Double_t VertexRes = 6e-4;
	const Double_t VertexResN = 2;
	TF1 *f = new TF1("gausl2","gaus+[3]+[4]*x"
			,RealDpKin-VertexRes*VertexResN,RealDpKin+VertexRes*VertexResN);
	f->SetParameter(1,RealDpKin);
	f->SetParameter(2,VertexRes);

	dpkinz->Fit("gausl2","R");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	TLatex *t = new TLatex(f->GetParameter(1)+f->GetParameter(2)
			,f->GetParameter(0)+f->GetParameter(3)+f->GetParameter(4)*f->GetParameter(1)
			,Form("\\Delta \\pm \\sigma = %2.1f \\pm %2.1f \\times 10^{-4}"
					,10000*(f->GetParameter(1)-RealDpKin)
							,10000*f->GetParameter(2)));
	t->SetTextSize(0.05);
	t->SetTextAlign(12);
	t->SetTextColor(2);
	t->Draw();
	
	c1->Print(TString(c1->GetName())+TString(".png"),"png");
}

DrawC_SelfCalculation2531()
{
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2531.root");
// 	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2663.root");
// 	TChain * t=new TChain("T");
// 	RootFileLoader(t,"./ROOTfiles/e06010_det_L_%d.root",2765,2768);
	
// 	T->Draw(">>LooseCut","abs(ReactPt_L.z)<.2 && abs(L.tr.tg_th)<.07 && abs(L.tr.tg_ph)<.04 && abs(L.tr.tg_dp)<.05 && L.tr.n==1");
// 	T->Draw(">>DeepCut","abs(ReactPt_L.z)<.05 && abs(L.tr.tg_th)<.05 && abs(L.tr.tg_ph)<.003 && L.tr.n==1");
	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.23 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1");
// 	T->Draw(">>TestCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.0062 && abs(L.tr.tg_ph)<.0035 && L.tr.n==1");
// 	T->Draw(">>TestCut2","abs(ReactPt_L.z)<.04 && abs(L.tr.tg_th)<.0062 && abs(L.tr.tg_ph)<.0035 && L.tr.n==1");


// 	TString CutUsed = "LooseCut";
	TString CutUsed = "ProdCut";
// 	TString CutUsed = "DeepCut";
// 	TString CutUsed = "TestCut";
// 	TString CutUsed = "TestCut2";
	TString DrawOption  = "COLZ";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);

	////////////////////////////////////////////////////////////////////////////////
	
	const Double_t RealMom = DIPOLE_MAG2MOM(0.453231);
	cout<<"NMR Momentum  = "<<RealMom<<endl;
// 	const Double_t NMRMomCorr = 0;

	const Double_t BeamEnergyLossOnDp
			= 1 *1e-3
			*(
			  0.0895//Beam window
			+0.4299//BeO
			+4.7435E-04*51.6//he4
			+4*0.1072//mid carbon foil
			 );
	const Double_t BeamEnergyLossOnDp_ReactZDep = 1 *1e-3*0.1072/(.2/3) ;//GeV/m
	const Double_t AfterScatteringEnergyLoss
			= 1e-3*(
					+4.7435E-04*70.1548		//He4
			+.1//Target Enclosure
			+3.2433E-03*(116.078-70.1548)	//Air
			+0.0798//HRS Entrance
			-.180*1.2//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
				   );
	TString RadCorrectionScatMom
			= Form("(%f+%f*ReactPt_L.z)"
			,BeamEnergyLossOnDp+AfterScatteringEnergyLoss
			,BeamEnergyLossOnDp_ReactZDep
				  );


	const UInt_t NExcitationStates = 4; //Carbon Excitation States
	const Double_t ExcitationEnergyList[]
			={0,0.00443891,0.00765420,0.009641,0.010844,0.011160,0.011828,0.012710,0.013352,0.014083,0.015110,0.015440,0.016105,0.016570,0.017230};

	
	cout<<"Total Correction on p = "<<RadCorrectionScatMom<<endl;
	
	const Double_t RealDpKin = ( 2*(12*.931494028-.000511*6)*1.2306)/(2*((12*.931494028-.000511*6) + 1.2306 - 1.2306*TMath::Cos(16/180.*3.14159)))/RealMom-1;

	TString exp_dp_kin = Form("(ExTgtCor_Lurb.dp- ( (( 2*(12*.931494028-.000511*6)*rb.e)/(2*((12*.931494028-.000511*6) + rb.e - (rb.e*cos(L.tr.tg_ph + 16./180*pi)/cos(L.tr.tg_ph))/sqrt(1/cos(L.tr.tg_ph)^2 + tan(ExTgtCor_Lurb.th)^2))))/%f-1-(%f) ) + %s)",RealMom,RealDpKin,RadCorrectionScatMom.Data());

	cout<<"exp_dp_kin = "<<exp_dp_kin<<endl;
	
	//////////////////////////////////////////////////////////////////////////////
	TCanvas * c1 = new TCanvas(TString("Carbon2531.")+CutUsed+"_DependencyCheck_WithRadCorrection_SelfCalc",CutUsed,1900,1000);
	c1->Divide(4,3);
	UInt_t idx=1;


// 	c1->cd(idx++);
// 	T->Draw("L.tr.tg_dp>>dP(200,-.05,.05)","","");
// 	dP->SetTitle("LHRS Momentum Diff");
// 	dP->SetXTitle("Target dp (GeV)");


	c1->cd(idx++);
	T->Draw(Form("%s>>dpkin1(800,-0.05,.05)",exp_dp_kin.Data()),"","");
	dpkin1->SetTitle("dp_kin (Angular Independant)");
	dpkin1->SetXTitle("dp_kin");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("ReactPt_L.z:%s>>dPvz(800,-.01,.005,800,-.3,.3)",exp_dp_kin.Data()),"",DrawOption);
	dPvz->SetTitle("React Z VS dp_kin");
	dPvz->SetXTitle("dp_kin");
	dPvz->SetYTitle("React Z (m)");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("urb.x:%s>>dPx(800,-.01,.005,800,-.004,.004)",exp_dp_kin.Data()),"",DrawOption);
	dPx->SetTitle("urb x VS dp_kin");
	dPx->SetXTitle("dp_kin");
	dPx->SetYTitle("urb x (m)");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("urb.y:%s>>dPy(800,-.01,.005,800,-.000,.008)",exp_dp_kin.Data()),"",DrawOption);
	dPy->SetTitle("urb y VS dp_kin");
	dPy->SetXTitle("dp_kin");
	dPy->SetYTitle("urb y (m)");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_th:%s>>dPth(800,-.01,.005,800,-.07,.07)",exp_dp_kin.Data()),"",DrawOption);
	dPth->SetTitle("Target Theta VS dp_kin");
	dPth->SetXTitle("dp_kin");
	dPth->SetYTitle("Target Theta");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_ph:%s>>dPph(800,-.01,.005,800,-.045,.045)",exp_dp_kin.Data()),"",DrawOption);
	dPph->SetTitle("Target Phi VS dp_kin");
	dPph->SetXTitle("dp_kin");
	dPph->SetYTitle("Target Phi");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}


	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_y:%s>>dPtgy(800,-.01,.005,800,-.06,.06)",exp_dp_kin.Data()),"",DrawOption);
	dPtgy->SetTitle("Target Y VS dp_kin");
	dPtgy->SetXTitle("dp_kin");
	dPtgy->SetYTitle("Target Y");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("ExTgtCor_Lurb.x:%s>>dPtgx(800,-.01,.005,800,-.01,.01)",exp_dp_kin.Data()),"",DrawOption);
	dPtgx->SetTitle("Target X VS dp_kin");
	dPtgx->SetXTitle("dp_kin");
	dPtgx->SetYTitle("Target X");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("L.tr.r_x:%s>>dPr_x(800,-.01,.005,800,-.8,.8)",exp_dp_kin.Data()),"",DrawOption);
	dPr_x->SetTitle("Rotated X VS dp_kin");
	dPr_x->SetXTitle("dp_kin");
	dPr_x->SetYTitle("Rotated X");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}
	c1->cd(idx++);
	T->Draw(Form("L.tr.r_th:%s>>dPr_th(800,-.01,.005,800,-.03,.03)",exp_dp_kin.Data()),"",DrawOption);
	dPr_th->SetTitle("Rotated Th VS dp_kin");
	dPr_th->SetXTitle("dp_kin");
	dPr_th->SetYTitle("Rotated Th");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}
	c1->cd(idx++);
	T->Draw(Form("L.tr.r_y:%s>>dPr_y(800,-.01,.005,800,-.02,.02)",exp_dp_kin.Data()),"",DrawOption);
	dPr_y->SetTitle("Rotated Y VS dp_kin");
	dPr_y->SetXTitle("dp_kin");
	dPr_y->SetYTitle("Rotated Y");
	
	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	c1->cd(idx++);
	T->Draw(Form("L.tr.r_ph:%s>>dPr_ph(800,-.01,.005,800,-.06,.06)",exp_dp_kin.Data()),"",DrawOption);
	dPr_ph->SetTitle("Rotated Ph VS dp_kin");
	dPr_ph->SetXTitle("dp_kin");
	dPr_ph->SetYTitle("Rotated Ph");
	
	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}


	c1->Print(TString(c1->GetName())+TString(".png"),"png");

	
	////////////////////////////////////////////////////////////////////////////////
	
	TCanvas * c1 = new TCanvas(TString("Carbon2531.")+CutUsed+"_WithRadCorrection_SelfCalc",CutUsed,1900,1000);
	c1->Divide(2,2);
	UInt_t idx=1;



	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_dp>>dp(800,-0.05,.05)"),"","");
// 	dpkin->SetTitle("dp");
	dp->SetXTitle("dp");
	
	c1->cd(idx++);
	T->Draw(Form("ExTgtCor_Lurb.dp>>extdp(800,-0.05,.05)",exp_dp_kin.Data()),"","");
// 	dpkin->SetTitle("dp_kin (Angular Independant)");
	extdp->SetXTitle("ExTgtCor_Lurb.dp");

	c1->cd(idx++);
	T->Draw(Form("%s>>dpkin(800,-0.05,.05)",exp_dp_kin.Data()),"","");
	dpkin->SetTitle("dp_kin (Angular Independant)");
	dpkin->SetXTitle("dp_kin");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}
	
	c1->cd(idx++);
	T->Draw(Form("%s>>dpkinz(800,-.01,.005)",exp_dp_kin.Data()),"","");
	dpkinz->SetTitle("dp_kin Zoomed (Angular Independant)");
	dpkinz->SetXTitle("dp_kin");

	const Double_t VertexRes = 6e-4;
	const Double_t VertexResN = 2;
	TF1 *f = new TF1("gausl2","gaus+[3]+[4]*x"
			,RealDpKin-VertexRes*VertexResN,RealDpKin+VertexRes*VertexResN);
	f->SetParameter(1,RealDpKin);
	f->SetParameter(2,VertexRes);

	dpkinz->Fit("gausl2","R");

	for(UInt_t exc=0; exc<NExcitationStates;exc++)
	{
		TLine *l = new TLine(RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,-1,RealDpKin-ExcitationEnergyList[exc]*.995751/RealMom,1000000);
		l->SetLineColor(6);
		l->Draw();
	}

	TLatex *t = new TLatex(f->GetParameter(1)+f->GetParameter(2)
			,f->GetParameter(0)+f->GetParameter(3)+f->GetParameter(4)*f->GetParameter(1)
			,Form("\\Delta \\pm \\sigma = %2.1f \\pm %2.1f \\times 10^{-4}"
					,10000*(f->GetParameter(1)-RealDpKin)
							,10000*f->GetParameter(2)));
	t->SetTextSize(0.05);
	t->SetTextAlign(12);
	t->SetTextColor(2);
	t->Draw();
	
	c1->Print(TString(c1->GetName())+TString(".png"),"png");
}

DrawHe3()
{
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2766.root");
// 	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2663.root");
	
	T->Draw(">>LooseCut","abs(ReactPt_L.z)<.2 && abs(L.tr.tg_th)<.07 && abs(L.tr.tg_ph)<.04 && abs(L.tr.tg_dp)<.05 && L.tr.n==1");
	T->Draw(">>DeepCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.05 && abs(L.tr.tg_ph)<.03 && L.tr.n==1");
	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1");
// 	T->Draw(">>TestCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1");

// 	TString CutUsed = "LooseCut";
// 	TString CutUsed = "DeepCut";
	TString CutUsed = "ProdCut";
// 	TString CutUsed = "TestCut";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);

	////////////////////////////////////////////////////////////////////////////////
	
	const Double_t RealMom = DIPOLE_MAG2MOM(0.448926);
	const Double_t NMRMomCorr = RealMom-1.21271;
	cout<<"NMR Momentum Correction = "<<NMRMomCorr<<endl;

	const Double_t BeamEnergyLossOnDp
			= 0.90561*1e-3*(
							0.075//Beam window
			+0.010//He4 Before Cell
			+0.038//Cell Entrance
			+0.036//Half Cell Length
						   );
	const Double_t BeamEnergyLossOnDp_ReactZDep = 0.90561* 0.036e-3/.2;//GeV/m
	const Double_t AfterScatteringEnergyLoss
			= 1e-3*(
					3.62796/20*0.036 //in cell H2
			+0.038*10	//Exit from cell wall
			+0.0228		//He4
			+.1//Target Enclosure
			+0.109	//Air
			+.1//HRS Entrance
				   );
	TString RadCorrectionScatMom
			= Form("(%f+%f*ReactPt_L.z)"
			,BeamEnergyLossOnDp+AfterScatteringEnergyLoss
			,BeamEnergyLossOnDp_ReactZDep
				  );
	
	cout<<"Total Correction on p = "<<RadCorrectionScatMom<<endl;
	
	const Double_t RealDpKin = ( 2*(3.0160293*.931494028-.000511*2)*1.2306)/(2*((3.0160293*.931494028-.000511*2) + 1.2306 - 1.2306*TMath::Cos(16/180.*3.14159)))/RealMom-1;

	//DM = (sqrt(PriKine_He3.W2)-(3.0160293*.931494028-.000511*2))
	const char * exp_p_kin_nocorr = "( -(((sqrt(PriKine_He3.W2)-(3.0160293*.931494028-.000511*2))*(sqrt(PriKine_He3.W2)-(3.0160293*.931494028-.000511*2)) + 2*(3.0160293*.931494028-.000511*2)*(sqrt(PriKine_He3.W2)-(3.0160293*.931494028-.000511*2)) - 2*(3.0160293*.931494028-.000511*2)*1.2306)/(2*((3.0160293*.931494028-.000511*2) + 1.2306 - 1.2306*cos(16/180.*3.14159)))) )";
	
	////////////////////////////////////////////////////////////////////////////////
	TCanvas * c1 = new TCanvas(CutUsed+"_DependencyCheck_WithRadCorrection",CutUsed,1900,1000);
	c1->Divide(3,3);
	UInt_t idx=1;


	c1->cd(idx++);
	T->Draw("L.tr.tg_dp>>dP(200,-.05,.05)","","");
	dP->SetTitle("LHRS Momentum Diff");
	dP->SetXTitle("Target dp (GeV)");

// 	c1->cd(idx++);
// 	T->Draw("L.tr.p>>P(200,1.11269,1.31269)","","");
// 	P->SetTitle("LHRS Momentum");
// 	P->SetXTitle("p (GeV)");

	c1->cd(idx++);
	T->Draw(Form("sqrt(PriKine_He3.W2)-0.951636*(%s+%f)>>W(800,2.75,2.9)",RadCorrectionScatMom.Data(),NMRMomCorr),"","");
	W->SetTitle("Inv Mass");
	W->SetXTitle("Invariant Mass (GeV)");
// 	W->SetYTitle("Target Theta L.tr.tg_th");
	TLine *l = new TLine((3.0160293*.931494028-.000511*2),-1,(3.0160293*.931494028-.000511*2),1000000);
	l->SetLineColor(6);
	l->Draw();

// 	c1->cd(idx++);
// 	T->Draw(Form("%s>>PKin(200,1.16,1.18)",exp_p_kin_nocorr),"","");
// 	PKin->SetTitle("Simplified LHRS Momentum (Angular Independant)");
// 	PKin->SetXTitle("p_kin (GeV)");
// // 	W->SetYTitle("Target Theta L.tr.tg_th");

	c1->cd(idx++);
	T->Draw(Form("%s/1.21269-1>>dpkin(800,-0.05,.05)",exp_p_kin_nocorr),"","");
	dpkin->SetTitle("dp_kin (Angular Independant)");
	dpkin->SetXTitle("dp_kin");
// 	W->SetYTitle("Target Theta L.tr.tg_th");

	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();
	
	c1->cd(idx++);
	T->Draw(Form("ReactPt_L.z:(%s+%s)/1.21269-1>>dPvz(800,-0.05,.05,800,-.3,.3)",exp_p_kin_nocorr,RadCorrectionScatMom.Data()),"","COLZ");
	dPvz->SetTitle("React Z VS dp_kin");
	dPvz->SetXTitle("dp_kin");
	dPvz->SetYTitle("React Z (m)");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_th:(%s+%s)/1.21269-1>>dPth(800,-0.05,.05,800,-.07,.07)",exp_p_kin_nocorr,RadCorrectionScatMom.Data()),"","COLZ");
	dPth->SetTitle("Target Theta VS dp_kin");
	dPth->SetXTitle("dp_kin");
	dPth->SetYTitle("Target Theta");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw(Form("L.tr.tg_ph:(%s+%s)/1.21269-1>>dPph(800,-0.05,.05,800,-.045,.045)",exp_p_kin_nocorr,RadCorrectionScatMom.Data()),"","COLZ");
	dPph->SetTitle("Target Phi VS dp_kin");
	dPph->SetXTitle("dp_kin");
	dPph->SetYTitle("Target Phi");

	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();
	
	c1->cd(idx++);
	T->Draw(Form("rb.x:(%s+%s)/1.21269-1>>dPx(800,-0.05,.05,800,-.004,.004)",exp_p_kin_nocorr,RadCorrectionScatMom.Data()),"","COLZ");
	dPx->SetTitle("rb x VS dp_kin");
	dPx->SetXTitle("dp_kin");
	dPx->SetYTitle("rb x (m)");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();
	
	c1->cd(idx++);
	T->Draw(Form("rb.y:(%s+%s)/1.21269-1>>dPy(800,-0.05,.05,800,-.0,.008)",exp_p_kin_nocorr,RadCorrectionScatMom.Data()),"","COLZ");
	dPy->SetTitle("rb y VS dp_kin");
	dPy->SetXTitle("dp_kin");
	dPy->SetYTitle("rb y (m)");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->Print(TString(c1->GetName())+TString(".png"),"png");
	
}

DrawH2Wrong()
{
	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2620.root");
	
	T->Draw(">>LooseCut","abs(ReactPt_L.z)<.2 && abs(L.tr.tg_th)<.07 && abs(L.tr.tg_ph)<.04 && abs(L.tr.tg_dp)<.05 && L.tr.n==1");
	T->Draw(">>DeepCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.05 && abs(L.tr.tg_ph)<.03 && L.tr.n==1");
	T->Draw(">>ProdCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1");
	T->Draw(">>TestCut","abs(ReactPt_L.z)<.17 && abs(L.tr.tg_th)<.062 && abs(L.tr.tg_ph)<.035 && L.tr.n==1");

// 	TString CutUsed = "LooseCut";
// 	TString CutUsed = "DeepCut";
	TString CutUsed = "ProdCut";
// 	TString CutUsed = "TestCut";
	TEventList *elist = (TEventList*)gDirectory->Get(CutUsed);
	assert(elist);
	T->SetEventList(elist);

	////////////////////////////////////////////////////////////////////////////////
	
	const Double_t RealMom = DIPOLE_MAG2MOM(0.448924);
	const Double_t NMRMomCorr = RealMom-1.21272;
	cout<<"NMR Momentum Correction = "<<NMRMomCorr<<endl;

	const Double_t BeamEnergyLossOnDp
			= 0.90561*1e-3*(
							0.075//Beam window
			+0.010//He4 Before Cell
			+0.038//Cell Entrance
			+0.036//Half Cell Length
						   );
	const Double_t BeamEnergyLossOnDp_ReactZDep = 0.90561* 0.036e-3/.2;//GeV/m
	const Double_t AfterScatteringEnergyLoss
			= 1e-3*(
					3.62796/20*0.036 //in cell H2
			+0.038*10	//Exit from cell wall
			+0.0228		//He4
			+.1//Target Enclosure
			+0.109	//Air
			+.1//HRS Entrance
				   );
	TString RadCorrectionScatMom
			= Form("(%f+%f*ReactPt_L.z)"
			,BeamEnergyLossOnDp+AfterScatteringEnergyLoss
			,BeamEnergyLossOnDp_ReactZDep
				  );
	
	cout<<"Total Correction on p = "<<RadCorrectionScatMom<<endl;
	
	const Double_t RealDpKin = ( 2*0.938272013*1.2306)/(2*(0.938272013 + 1.2306 - 1.2306*TMath::Cos(16/180.*3.14159)))/RealMom-1;
	
	////////////////////////////////////////////////////////////////////////////////
	TCanvas * c1 = new TCanvas(CutUsed+"_Canvas",CutUsed,1900,1000);
	c1->Divide(3,3);
	UInt_t idx=1;


	c1->cd(idx++);
	T->Draw("L.tr.tg_dp>>dP(200,-.05,.05)","","");
	dP->SetTitle("LHRS Momentum Diff");
	dP->SetXTitle("Target dp (GeV)");

	c1->cd(idx++);
	T->Draw("L.tr.p>>P(200,1.11269,1.31269)","","");
	P->SetTitle("LHRS Momentum");
	P->SetXTitle("p (GeV)");

	c1->cd(idx++);
	T->Draw("sqrt(PriKine.W2)>>W(200,0.92,.96)","","");
	W->SetTitle("Inv Mass");
	W->SetXTitle("Invariant Mass (GeV)");
// 	W->SetYTitle("Target Theta L.tr.tg_th");
	TLine *l = new TLine(0.938272013,-1,0.938272013,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw("(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))))>>PKin(200,1.16,1.18)","","");
	PKin->SetTitle("Simplified LHRS Momentum (Angular Independant)");
	PKin->SetXTitle("p_kin (GeV)");
// 	W->SetYTitle("Target Theta L.tr.tg_th");

	c1->cd(idx++);
	T->Draw("(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))))/1.21269-1>>dpkin(800,-0.04,-.03)","","");
	dpkin->SetTitle("dp_kin (Angular Independant)");
	dpkin->SetXTitle("dp_kin");
// 	W->SetYTitle("Target Theta L.tr.tg_th");

	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();
	
	c1->cd(idx++);
	T->Draw("ReactPt_L.z:(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))))/1.21269-1>>dPvz(800,-0.04,-.03,800,-.3,.3)","","COLZ");
	dPvz->SetTitle("React Z VS dp_kin");
	dPvz->SetXTitle("dp_kin");
	dPvz->SetYTitle("React Z (m)");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw("L.tr.tg_th:(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))))/1.21269-1>>dPth(800,-0.04,-.03,800,-.07,.07)","","COLZ");
	dPth->SetTitle("Target Theta VS dp_kin");
	dPth->SetXTitle("dp_kin");
	dPth->SetYTitle("Target Theta");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

	c1->cd(idx++);
	T->Draw("L.tr.tg_ph:(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))))/1.21269-1>>dPph(800,-0.04,-.03,800,-.045,.045)","","COLZ");
	dPph->SetTitle("Target Phi VS dp_kin");
	dPph->SetXTitle("dp_kin");
	dPph->SetYTitle("Target Phi");

	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();
	
	c1->cd(idx++);
	T->Draw("L.tr.x:(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))))/1.21269-1>>dPx(800,-0.04,-.03,800,-.8,.8)","","COLZ");
	dPx->SetTitle("VDC x VS dp_kin");
	dPx->SetXTitle("dp_kin");
	dPx->SetYTitle("VDC x (m)");
	TLine *l = new TLine(RealDpKin,-1,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();

// 	c1->Print(Form("H2Check.%s.png",CutUsed.Data()),"png");
	////////////////////////////////////////////////////////////////////////////////
	TCanvas * c1 = new TCanvas(CutUsed+"_Canvas_RadCorrection",CutUsed,1900,1000);
	c1->Divide(1,2);
	UInt_t idx=1;



	c1->cd(idx++);
	T->Draw(Form("sqrt(PriKine.W2)+%s/0.0603824>>W(200,0.92,.96)"
			,RadCorrectionScatMom.Data()),"","");
	W->SetTitle("Inv Mass");
	W->SetXTitle("Invariant Mass (GeV)");
// 	W->SetYTitle("Target Theta L.tr.tg_th");

	c1->cd(idx++);
	T->Draw(Form("(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))) +%s)/1.21269-1>>dpkin(800,-0.04,-.03)",RadCorrectionScatMom.Data()),"","");
	dpkin->SetTitle("dp_kin (Angular Independant)");
	dpkin->SetXTitle("dp_kin");
// 	W->SetYTitle("Target Theta L.tr.tg_th");

	TLine *l = new TLine(RealDpKin,0,RealDpKin,1000000);
	l->SetLineColor(6);
	l->Draw();


// 	c1->cd(idx++);
// 	T->Draw("ReactPt_L.z:(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))))/1.21269-1>>dPvz(800,-0.04,-.03,800,-.3,.3)","","COLZ");
// 	dPvz->SetTitle("React Z VS dp_kin");
// 	dPvz->SetXTitle("dp_kin");
// 	dPvz->SetYTitle("React Z (m)");
// 
// 	c1->cd(idx++);
// 	T->Draw("L.tr.tg_th:(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))))/1.21269-1>>dPth(800,-0.04,-.03,800,-.07,.07)","","COLZ");
// 	dPth->SetTitle("Target Theta VS dp_kin");
// 	dPth->SetXTitle("dp_kin");
// 	dPth->SetYTitle("Target Theta");
// 
// 	c1->cd(idx++);
// 	T->Draw("L.tr.tg_ph:(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))))/1.21269-1>>dPph(800,-0.04,-.03,800,-.045,.045)","","COLZ");
// 	dPph->SetTitle("Target Phi VS dp_kin");
// 	dPph->SetXTitle("dp_kin");
// 	dPph->SetYTitle("Target Phi");
// 
// 	c1->cd(idx++);
// 	T->Draw("L.tr.x:(( 2*sqrt(PriKine.W2)*rb.e)/(2*(sqrt(PriKine.W2) + rb.e - rb.e*cos(16/180.*3.14159))))/1.21269-1>>dPx(800,-0.04,-.03,800,-.8,.8)","","COLZ");
// 	dPx->SetTitle("VDC x VS dp_kin");
// 	dPx->SetXTitle("dp_kin");
// 	dPx->SetYTitle("VDC x (m)");

// 	c1->Print(Form("H2Check.%s.AfterCorrection.png",CutUsed.Data()),"png");
	
}

RasterCheck()
{
// 	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2766.root");
	//
	TCanvas * c1 = new TCanvas("RasterCheck","Raster Checks",900,400);
// 	c1->Divide(3,1);
	c1->Divide(3,1);
	UInt_t idx=1;

	c1->cd(idx++);
	T->Draw("rb.x>>rb_x(100,-.008,.008)");
	T->Draw("urb.x>>urb_x(100,-.008,.008)","","SAME");
	rb_x->SetLineColor(2);
	urb_x->SetLineColor(4);
	c1->cd(idx++);
	T->Draw("rb.y>>rb_y(100,-.008,.008)");
	T->Draw("urb.y>>urb_y(100,-.008,.008)","","SAME");
	rb_y->SetLineColor(2);
	urb_y->SetLineColor(4);
	
	c1->cd(idx++);
	T->Draw("rb.y:rb.x>>rb_xy(400,-.008,.008,400,-.008,.008)","","COLZ");
}

RasterCorrelation()
{
// 	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2766.root");
	//
	TCanvas * c1 = new TCanvas("RasterCorrelation","Raster Correlation Checks",900,900);
// 	c1->Divide(3,1);
	c1->Divide(2,2);
	UInt_t idx=1;

	c1->cd(idx++);
	T->Draw("rb.x:rawcurx");
	c1->cd(idx++);
	T->Draw("rb.x:rawcury");
	c1->cd(idx++);
	T->Draw("rb.y:rawcurx");
	c1->cd(idx++);
	T->Draw("rb.y:rawcury");

	c1->Print(TString(c1->GetName())+".png","png");
}

RasterCorrelation2()
{
// 	TFile *_file0 = TFile::Open("ROOTfiles/e06010_det_L_2766.root");
	//
	TCanvas * c1 = new TCanvas("RasterCorrelation2.Check","Raster Correlation Checks",900,900);
// 	c1->Divide(3,1);
	c1->Divide(2,2);
	UInt_t idx=1;

	c1->cd(idx++);
	T->Draw("rb.x:urb.x");
	c1->cd(idx++);
	T->Draw("rb.x:urb.y");
	c1->cd(idx++);
	T->Draw("rb.y:urb.x");
	c1->cd(idx++);
	T->Draw("rb.y:urb.y");

	c1->Print(TString(c1->GetName())+".png","png");
}

LOpticsCalib()
{
	
	gStyle->SetOptStat(11);

// 	VertexZ5832();
	OpticsAcceptanceRCutCheck();
// 	DrawC_SelfCalculation2580_dp();
// 	AngularAcceptanceCheckByFoil();
// 	AngularAcceptanceCheck();
// 	DrawHe3();
}