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

static const UInt_t NFiols = 7;
Double_t VertexRes = 0.006;

// Survey, used for calibration
//Double_t targetfoils[NFiols]={0.0019+0.001565-0.2,0.0019+0.001565-0.2+0.2/3.,0.0019+0.001565-0.2+0.2/3.*2.,0.0019+0.001565,0.0019+0.001565+0.2/3.,0.0019+0.001565+0.2/3.*2.,+0.0019+0.001565+0.2};

// Fit run 3912
Double_t targetfoils[]={-1.965143e-01,    -1.296833e-01,  -6.286111e-02, 3.785412e-03,    7.027552e-02,   1.363940e-01,   2.028485e-01}; //Extracted from fitting


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

	const UInt_t VertexResN = 6;

	Double_t fArbitaryVertexShift[100]={0};

	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 +- %0.1f mm"
						,1000*(f->GetParameter(1)-targetfoils[FoilID])
						,1000*f->GetParError(1)));
		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]);

	double aver = 0;
	for(UInt_t FoilID=0; FoilID<NFiols; FoilID++)
	{
		printf("fArbitaryVertexShift[%d] = %e;\n"
				,FoilID,fArbitaryVertexShift[FoilID]);

		aver += fArbitaryVertexShift[FoilID];
	}

	printf("Double_t targetfoils[]={");
	for(UInt_t FoilID=0; FoilID<NFiols; FoilID++)
	{
		printf("%e,\t"
				,targetfoils[FoilID]+fArbitaryVertexShift[FoilID]);

		aver += fArbitaryVertexShift[FoilID];
	}
	printf("0}; //Extracted from fitting \n");

	printf("Overall Shift = %f;\n",aver/NFiols);
}


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

//	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_L.z>>CarbonFoil%d(500,-0.4,0.4)",Run),
			"abs(ExTgtCor_L.dp)<0.05&&L.tr.n==1&&abs(ExTgtCor_L.th)<.05&&abs(ExTgtCor_L.ph)<.03 && L.a1.asum_c>100 && L.prl1.asum_c>10");

	DrawVertex();
	TH1F* h= (TH1F*)gDirectory->FindObjectAny(Form("CarbonFoil%d",Run));
	assert(h);
	h -> SetTitle(Form("Vertex for multi-cabon run %d; Vertex (m)",Run));

	FitVertex_New(h);

	SaveCanvas(c1,TString("PolarityStudy_")+c1->GetName());
}



void PolarityStudy()
{
//	VertexZ_Fit(5832);
//	VertexZ_Fit(5415);
	//	VertexZ_Fit(5460);
//	VertexZ_Fit(5443);
//	VertexZ_Fit(3924);
	//	VertexZ_Fit(3916);
		//	VertexZ_Fit(4130);


	VertexZ_Fit(3924);
	VertexZ_Fit(4221);
}