//*******************************************************
//* Last update : 20210803 by Masato Kimura
//*
//* Supplementary material for Phys.Rev.D100, 032002
//* Contact	: masato@kylab.sci.waseda.ac.jp (M.Kimura),
//*			: kohei.yorita@waseda.jp(K.Yorita)
//*
//* How to Run : 
//*	-1 $ root
//*	-2 > .L LeffwQy.C
//* -3 > Leff(your_input_field_in_kV/cm)
//*******************************************************

///////////Declaration
Int_t LeffwQy(Double_t field_kVcm = 0.0, TString mode = "BestFit", Double_t parkB = 0.0, Double_t parAlpha_0 = 0.0, Double_t parD_Alpha = 0.0, Double_t parGamma = 0.0, Double_t parDelta = 0.0);
Double_t funcdetailLnhrd(Double_t Er_MeV);
Double_t funcdetailBirks(Double_t Er_MeV, Double_t parkB);
Double_t funcdetailTIB(Double_t Er_MeV, Double_t parEf_kVcm, Double_t parAlpha, Double_t parvarSigma, Bool_t flagLeff);
Double_t funcLeff(Double_t *x, Double_t *par);
Double_t funcQy(Double_t *x, Double_t *par);

/////////// Main
Int_t LeffwQy(Double_t field_kVcm, TString mode, Double_t parkB, Double_t parAlpha_0, Double_t parD_Alpha, Double_t parGamma, Double_t parDelta){
	////////// Parameter seting
	if("BestFit" == mode){
		parkB		= 3.1e-4;
		parAlpha_0	= 1.0;
		parD_Alpha	= 8.9e-4;
		parGamma	= 1.2;
		parDelta	= 5.8e-1;
	} else
	if("E-calib_+03" == mode){
		parkB		= 3.7e-4;
		parAlpha_0	= 1.0;
		parD_Alpha	= 9.2e-4;
		parGamma	= 1.15;
		parDelta	= 5.75e-1;
	} else
	if("E-calib_-03" == mode){
		parkB		= 2.5e-4;
		parAlpha_0	= 1.0;
		parD_Alpha	= 8.9e-4;
		parGamma	= 1.16;
		parDelta	= 5.76e-1;
	} else
	if("TOFarm_+1" == mode){
		parkB		= 3.4e-4;
		parAlpha_0	= 1.0;
		parD_Alpha	= 8.2e-4;
		parGamma	= 1.16;
		parDelta	= 5.77e-1;
	} else
	if("TOFarm_-1" == mode){
		parkB		= 3.4e-4;
		parAlpha_0	= 1.0;
		parD_Alpha	= 8.2e-4;
		parGamma	= 1.15;
		parDelta	= 5.75e-1;
	} else
	if("TOFt0_+1" == mode){
		parkB		= 3.5e-4;
		parAlpha_0	= 1.0;
		parD_Alpha	= 8.9e-4;
		parGamma	= 1.15;
		parDelta	= 5.75e-1;
	} else
	if("TOFt0_-1" == mode){
		parkB		= 2.7e-4;
		parAlpha_0	= 1.0;
		parD_Alpha	= 8.8e-4;
		parGamma	= 1.15;
		parDelta	= 5.77e-1;
	} else
	if("g2g1_+20" == mode){
		parkB		= 3.1e-4;
		parAlpha_0	= 1.0;
		parD_Alpha	= 12.0e-4;
		parGamma	= 0.93;
		parDelta	= 5.78e-1;
	} else
	if("g2g1_-20" == mode){
		parkB		= 3.1e-4;
		parAlpha_0	= 1.0;
		parD_Alpha	= 5.4e-4;
		parGamma	= 1.62;
		parDelta	= 5.86e-1;
	} else
	if("Manual" == mode){
		std::cout << "Use your input value." << std::endl;
	} else
	return -1;

	////////// Define Function
	TF1 *fLeff = new TF1("fLeff", funcLeff, 1.0e-3, 1.0, 6);
	fLeff->SetParNames("Field_kV/cm",	// 0
					   "k_{B}",		 	// 1
					   "#alpha_{0}",	// 2
					   "D_{#alpha}",	// 3
					   "#gamma",		// 4
					   "#delta");		// 5
	fLeff->SetParameters(field_kVcm,	// 0
						 parkB,			// 1
						 parAlpha_0,	// 2
						 parD_Alpha,	// 3
						 parGamma,		// 4
						 parDelta);		// 5

	TF1 *fQy = new TF1("fQy", funcQy, 1.0e-3, 1.0, 6);
	fQy->SetParNames("Field_kV/cm",	// 0
					 "k_{B}",	 	// 1
					 "#alpha_{0}",	// 2
					 "D_{#alpha}",	// 3
					 "#gamma",		// 4
					 "#delta");		// 5
	fQy->SetParameters(field_kVcm,	// 0
					   parkB,		// 1
					   parAlpha_0,	// 2
					   parD_Alpha,	// 3
					   parGamma,	// 4
					   parDelta);	// 5

	////////// Draw
	TCanvas *c1=new TCanvas("c1","c1",900,600);
	fLeff->Draw();

	TCanvas *c2=new TCanvas("c2","c2",900,600);
	fQy->Draw();

	return 1;
}

Int_t Leff_wBand(const Double_t field_kVcm){
	////////// Define Function
	TF1 *fLeff = new TF1("fLeff", funcLeff, 1.0e-3, 1.0, 6);
	fLeff->SetLineColor(1);
	TF1 *fBand[4][2];
	TGraphAsymmErrors *gBand;
	for(Int_t ierr=0; ierr<4; ++ierr){
		for(Int_t idx=0; idx<2; ++idx){
			fBand[ierr][idx]=new TF1(Form("fBand_%d_%d",ierr,idx), funcLeff, 1.0e-3, 1.0, 6);
			fBand[ierr][idx]->SetLineColor(ierr+2);
		}
	}
	gBand=new TGraphAsymmErrors();
	gBand->SetFillStyle(1001);
	gBand->SetFillColor(kGray);
	gBand->SetLineColor(kGray);

	////////// Parameter seting
	Double_t parkB = 0.0, parAlpha_0 = 0.0, parD_Alpha = 0.0, parGamma = 0.0, parDelta = 0.0;
	// "BestFit"
	parkB		= 3.1e-4;
	parAlpha_0	= 1.0;
	parD_Alpha	= 8.9e-4;
	parGamma	= 1.2;
	parDelta	= 5.8e-1;
	fLeff->SetParameters(field_kVcm, parkB, parAlpha_0, parD_Alpha, parGamma, parDelta);

	// "E-calib_+03"
	parkB		= 3.7e-4;
	parAlpha_0	= 1.0;
	parD_Alpha	= 9.2e-4;
	parGamma	= 1.15;
	parDelta	= 5.75e-1;
	fBand[0][0]->SetParameters(field_kVcm, parkB, parAlpha_0, parD_Alpha, parGamma, parDelta);

	// "E-calib_-03"
	parkB		= 2.5e-4;
	parAlpha_0	= 1.0;
	parD_Alpha	= 8.9e-4;
	parGamma	= 1.16;
	parDelta	= 5.76e-1;
	fBand[0][1]->SetParameters(field_kVcm, parkB, parAlpha_0, parD_Alpha, parGamma, parDelta);

	// "TOFarm_+1"
	parkB		= 3.4e-4;
	parAlpha_0	= 1.0;
	parD_Alpha	= 8.2e-4;
	parGamma	= 1.16;
	parDelta	= 5.77e-1;
	fBand[1][0]->SetParameters(field_kVcm, parkB, parAlpha_0, parD_Alpha, parGamma, parDelta);

	// "TOFarm_-1"
	parkB		= 3.4e-4;
	parAlpha_0	= 1.0;
	parD_Alpha	= 8.2e-4;
	parGamma	= 1.15;
	parDelta	= 5.75e-1;
	fBand[1][1]->SetParameters(field_kVcm, parkB, parAlpha_0, parD_Alpha, parGamma, parDelta);

	// "TOFt0_+1"
	parkB		= 3.5e-4;
	parAlpha_0	= 1.0;
	parD_Alpha	= 8.9e-4;
	parGamma	= 1.15;
	parDelta	= 5.75e-1;
	fBand[2][0]->SetParameters(field_kVcm, parkB, parAlpha_0, parD_Alpha, parGamma, parDelta);

	// "TOFt0_-1"
	parkB		= 2.7e-4;
	parAlpha_0	= 1.0;
	parD_Alpha	= 8.8e-4;
	parGamma	= 1.15;
	parDelta	= 5.77e-1;
	fBand[2][1]->SetParameters(field_kVcm, parkB, parAlpha_0, parD_Alpha, parGamma, parDelta);

	// "g2g1_+20"
	parkB		= 3.1e-4;
	parAlpha_0	= 1.0;
	parD_Alpha	= 12.0e-4;
	parGamma	= 0.93;
	parDelta	= 5.78e-1;
	fBand[3][0]->SetParameters(field_kVcm, parkB, parAlpha_0, parD_Alpha, parGamma, parDelta);

	// "g2g1_-20" 
	parkB		= 3.1e-4;
	parAlpha_0	= 1.0;
	parD_Alpha	= 5.4e-4;
	parGamma	= 1.62;
	parDelta	= 5.86e-1;
	fBand[3][1]->SetParameters(field_kVcm, parkB, parAlpha_0, parD_Alpha, parGamma, parDelta);


	////////// Draw Band
	const Double_t energyRange[2] = {30.0e-3, 200.0e-3};
	Double_t currentEnergy = 0.0, currentLeff = 0.0, convolUpper = 0.0, convolLower = 0.0;
	for(Int_t ip=0; ip<1000; ++ip){
		currentEnergy = energyRange[0]*TMath::Power(10.0, ip*TMath::Log10(energyRange[1]/energyRange[0])/999.0);
		convolUpper = 0.0;
		convolLower = 0.0;
		currentLeff = fLeff->Eval(currentEnergy);
		for(Int_t ierr=0; ierr<4; ++ierr){
			if(field_kVcm==0.0 && 3==ierr) continue;
			convolUpper = TMath::Hypot(convolUpper, fBand[ierr][0]->Eval(currentEnergy)-currentLeff);
			convolLower = TMath::Hypot(convolUpper, currentLeff-fBand[ierr][1]->Eval(currentEnergy));
		}
		gBand->SetPoint(ip, currentEnergy, currentLeff);
		gBand->SetPointError(ip, 0.0, 0.0, convolLower, convolUpper);
	}

	////////// Draw Canvas
	TCanvas *c2=new TCanvas("c2","c2",900,600);
	fLeff->Draw();
	gBand->Draw("l");
	for(Int_t ierr=0; ierr<4; ++ierr){
		if(field_kVcm==0.0 && 3==ierr) continue;
		for(Int_t idx=0; idx<2; ++idx){
			fBand[ierr][idx]->Draw("same");
		}
	}
	fLeff->Draw("same"); // redraw
}

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

Double_t funcdetailLnhrd(Double_t Er_MeV){
	Double_t Er_keV	= Er_MeV*1.0e3;

	Double_t parZ	= 18.0;
	Double_t parA	= 40.0;
	Double_t parEps	= 11.5*Er_keV*TMath::Power(parZ, -7.0/3.0);
	Double_t park	= 0.133*TMath::Power(parZ, 2.0/3.0)*TMath::Power(parA, -1.0/2.0);
	Double_t parg	= 3.0*TMath::Power(parEps, 0.15) + 0.7*TMath::Power(parEps, 0.6) + parEps;
	Double_t parLndhrd	= park*parg/(1.0+park*parg);

	return parLndhrd;
}

Double_t funcdetailBirks(Double_t Er_MeV, Double_t parkB){
	Double_t pardEdx= TMath::Power(1.208 + 2.527*TMath::Power(Er_MeV, 0.1596), 1.0/0.1596); // read from Astrophys 30 12 (Mei et al) Fig.4
	Double_t parBrks= 1.0/(1.0+parkB*pardEdx);

	return parBrks;
}

Double_t funcdetailTIB(Double_t Er_MeV, Double_t parEf_kVcm, Double_t parAlpha, Double_t parvarSigma, Bool_t flagLeff){
	Double_t Er_keV		= Er_MeV*1.0e3;
	Double_t parEf_Vcm	= parEf_kVcm*1.0e3;
	Double_t wv_keV		= 19.5*1.0e-3;

	Double_t parLeff;
	if(flagLeff) parLeff = funcdetailLnhrd(Er_MeV);
	else parLeff = 1.0;

	Double_t nIon		= Er_keV/wv_keV * parLeff * (1.0/(1.0+parAlpha));
	Double_t parxi		= nIon*parvarSigma;

	return 1.0 - TMath::Log(1.0+parxi)/parxi; // return recombination probability
}

Double_t funcLeff(Double_t *x, Double_t *par){
	Double_t Er_MeV		= x[0];
	Double_t parEf_kVcm	= par[0];
	Double_t parEf_Vcm	= parEf_kVcm * 1.0e3;
	Double_t parkB		= par[1];
	Double_t parAlpha_0	= par[2];
	Double_t parD_Alpha	= par[3];
	Double_t parGamma	= par[4];
	Double_t parDelta	= par[5];

	Double_t parLndhrd	= funcdetailLnhrd(Er_MeV);
	Double_t parMei		= parLndhrd * funcdetailBirks(Er_MeV, parkB);

	if(0.0==parEf_kVcm) return parMei;

	Double_t parAlpha	= parAlpha_0 * TMath::Exp(-1.0*parEf_Vcm*parD_Alpha);
	Double_t parvarSigma	= parGamma * TMath::Power(parEf_Vcm, -1.0*parDelta);

	Double_t parTIB		= funcdetailTIB(Er_MeV, parEf_kVcm, parAlpha, parvarSigma, kTRUE);
	return parMei * (1.0-(1.0/(1.0+parAlpha)*(1.0-parTIB)));
}

Double_t funcQy(Double_t *x, Double_t *par){
	Double_t Er_MeV		= x[0];
	Double_t Er_keV		= Er_MeV*1.0e3;
	Double_t parEf_kVcm	= par[0];
	Double_t parEf_Vcm	= parEf_kVcm * 1.0e3;
	Double_t parkB		= par[1];
	Double_t parAlpha_0	= par[2];
	Double_t parD_Alpha	= par[3];
	Double_t parGamma	= par[4];
	Double_t parDelta	= par[5];

	Double_t parAlpha	= parAlpha_0 * TMath::Exp(-1.0*parEf_Vcm*parD_Alpha);
	Double_t parvarSigma	= parGamma * TMath::Power(parEf_Vcm, -1.0*parDelta);

	Double_t parLndhrd	= funcdetailLnhrd(Er_MeV);
	Double_t wv_keV		= 19.5*1.0e-3;
	Double_t parTIB		= funcdetailTIB(Er_MeV, parEf_kVcm, parAlpha, parvarSigma, true);

	Double_t parNe	= parLndhrd*(Er_keV/wv_keV)*(1.0/(1.0+parAlpha))*(1.0-parTIB);
	return parNe/Er_keV;
}