//
// this program tests the computation of gA cross sections from
// UPC cross sections for different classes of fwd neutrons
//

// c++ headers
#include <iostream>
#include <fstream>

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

void UPC_XS(Double_t prob_1_0n0n, Double_t prob_1_Xn0n, Double_t prob_1_XnXn,
	    Double_t prob_2_0n0n, Double_t prob_2_Xn0n, Double_t prob_2_XnXn,
	    Double_t flux_1, Double_t flux_2, Double_t xs_1, Double_t xs_2)
// in this case, one gives the fluxes, probabilities of neutrons classes
// and gA cross sections at both rapidities (denoted by 1 and 2 here) to
// obtain the UPC cross sections
// NOTE: flux_2 corresponds to the higher WgA
// NOTE: Xn0n stands for (Xn0n+0nXn)
{
  // check probs
  Double_t prob_1 = prob_1_0n0n+prob_1_Xn0n+prob_1_XnXn;
  Double_t prob_2 = prob_2_0n0n+prob_2_Xn0n+prob_2_XnXn;

  cout << endl << "   tot prob_1 = " << prob_1 << " ... tot prob_2 = " << prob_2 << endl;

  // compute UPC cross sections
  Double_t UPC_XS_0n0n = flux_1*prob_1_0n0n*xs_1+flux_2*prob_2_0n0n*xs_2;
  Double_t UPC_XS_Xn0n = flux_1*prob_1_Xn0n*xs_1+flux_2*prob_2_Xn0n*xs_2;  
  Double_t UPC_XS_XnXn = flux_1*prob_1_XnXn*xs_1+flux_2*prob_2_XnXn*xs_2;

  // print out
  cout << "   UPC XSs: " << endl;
  cout << "     total = " << (UPC_XS_0n0n+UPC_XS_Xn0n+UPC_XS_XnXn) << endl;
  cout << "     0n0n = " << UPC_XS_0n0n  << endl;
  cout << "     Xn0n = " << UPC_XS_Xn0n  << endl;
  cout << "     XnXn = " << UPC_XS_XnXn  << endl << endl;;  
}


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

void Test_UPC_XS()
{
  // data for xs and flux from Michal Krelina (hs-model for rho)
  // data for probabilities from Michal Broz

  // y=0
  cout << " Computing case y = 0 " << endl;
  UPC_XS(0.82942,0.130901,0.0394531,0.82942,0.130901,0.0394531,
	 128.076,128.076,2.26123652,2.26123652);

  // |y| = 0.8
  cout << " Computing case |y| = 0.8 " << endl;
  UPC_XS(0.856358,0.110316,0.0331726,0.791281,0.159858,0.0485461,
	 153.063,103.118,2.171556597,2.37201377);
 

}

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

Double_t get_gA_XS(Double_t flux_1m, Double_t flux_2m,
		   Double_t flux_1p, Double_t flux_2p,
		   Double_t xs_1, Double_t xs_2, Bool_t opt)
// opt = kTRUE, eq (14), kFALSE eq (15)
// Phys.Rev. C96 (2017) no.1, 015203
// https://inspirehep.net/record/1491190?ln=en
{
  Double_t D = flux_2p*flux_1m - flux_1p*flux_2m;
  Double_t N  = (opt ? (flux_2p*xs_1 - flux_1p*xs_2) : (flux_1m*xs_2 - flux_2m*xs_1));
  return N/D;
}

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

void gA_XS(Double_t prob_1_0n0n, Double_t prob_1_Xn0n, Double_t prob_1_XnXn,
	   Double_t prob_2_0n0n, Double_t prob_2_Xn0n, Double_t prob_2_XnXn,
	   Double_t flux_1, Double_t flux_2,
	   Double_t xs_0n0n, Double_t xs_Xn0n, Double_t xs_XnXn)
// in this case, one gives the fluxes, probabilities of neutrons classes
// and UPC cross sections for neutron classes to
// obtain the gA cross sections from three pairs of data
// NOTE: flux_2 corresponds to the higher WgA
// NOTE: Xn0n stands for (Xn0n+0nXn)
{
  // check probs
  Double_t prob_1 = prob_1_0n0n+prob_1_Xn0n+prob_1_XnXn;
  Double_t prob_2 = prob_2_0n0n+prob_2_Xn0n+prob_2_XnXn;

  cout << endl << "   tot prob_1 = " << prob_1 << " ... tot prob_2 = " << prob_2 << endl;

  // get cross sections for flux1 != flux2
  Double_t gA_0n0n_Xn0n_1 = get_gA_XS(prob_1_0n0n*flux_1,prob_1_Xn0n*flux_1,
				      prob_2_0n0n*flux_2,prob_2_Xn0n*flux_2,
				      xs_0n0n,xs_Xn0n, kTRUE);
  Double_t gA_0n0n_Xn0n_2 = get_gA_XS(prob_1_0n0n*flux_1,prob_1_Xn0n*flux_1,
				      prob_2_0n0n*flux_2,prob_2_Xn0n*flux_2,
				      xs_0n0n,xs_Xn0n, kFALSE);
  
  Double_t gA_0n0n_XnXn_1 = get_gA_XS(prob_1_0n0n*flux_1,prob_1_XnXn*flux_1,
				      prob_2_0n0n*flux_2,prob_2_XnXn*flux_2,
				      xs_0n0n,xs_XnXn, kTRUE);
  Double_t gA_0n0n_XnXn_2 = get_gA_XS(prob_1_0n0n*flux_1,prob_1_XnXn*flux_1,
				      prob_2_0n0n*flux_2,prob_2_XnXn*flux_2,
				      xs_0n0n,xs_XnXn, kFALSE);

  Double_t gA_Xn0n_XnXn_1 = get_gA_XS(prob_1_Xn0n*flux_1,prob_1_XnXn*flux_1,
				      prob_2_Xn0n*flux_2,prob_2_XnXn*flux_2,
				      xs_Xn0n,xs_XnXn, kTRUE);
  Double_t gA_Xn0n_XnXn_2 = get_gA_XS(prob_1_Xn0n*flux_1,prob_1_XnXn*flux_1,
				      prob_2_Xn0n*flux_2,prob_2_XnXn*flux_2,
				      xs_Xn0n,xs_XnXn, kFALSE);

  // print out
  cout << "   gA XSs: " << endl;
  cout << "     from (0n0n,Xn0n), = (" << gA_0n0n_Xn0n_1 <<","<< gA_0n0n_Xn0n_2 <<")"<< endl;
  cout << "     from (0n0n,XnXn), = (" << gA_0n0n_XnXn_1 <<","<< gA_0n0n_XnXn_2 <<")"<< endl;
  cout << "     from (Xn0n,XnXn), = (" << gA_Xn0n_XnXn_1 <<","<< gA_Xn0n_XnXn_2 <<")"<< endl << endl;
 
}

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

void gA_XS_zero(Double_t prob_1_0n0n, Double_t prob_1_Xn0n, Double_t prob_1_XnXn,
		Double_t flux_1, Double_t xs_0n0n, Double_t xs_Xn0n, Double_t xs_XnXn)
// assuming both fluxes are equal (ie y=0)
// in this case, one gives the fluxes, probabilities of neutrons classes
// and UPC cross sections for neutron classes to
// obtain the gA cross sections from three pairs of data
// NOTE: Xn0n stands for (Xn0n+0nXn)
{
  // check probs
  Double_t prob_1 = prob_1_0n0n+prob_1_Xn0n+prob_1_XnXn;

  cout << endl << "   tot prob_1 = " << prob_1  << endl;

  Double_t gA_0n0n = xs_0n0n/(2.0*flux_1*prob_1_0n0n);
  Double_t gA_Xn0n = xs_Xn0n/(2.0*flux_1*prob_1_Xn0n);
  Double_t gA_XnXn = xs_XnXn/(2.0*flux_1*prob_1_XnXn);

  // print out
  cout << "   gA XSs: " << endl;
  cout << "     from 0n0n = " << gA_0n0n << endl;
  cout << "     from Xn0n = " << gA_Xn0n << endl;
  cout << "     from XnXn = " << gA_XnXn << endl << endl;
 
}

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

void oldTest_XS()
{
  // data for flux from Michal Krelina
  // data for probabilities from Michal Broz
  // data for UPC from test above (that is hot-spot model for rho)

  // y=0
  // cout << " Computing case y = 0 " << endl;
  // gA_XS_zero(0.82942,0.130901,0.0394531,128.076,480.417,75.8205,22.852);
  
  // |y| = 0.8
  cout << " Computing case |y| = 0.8 " << endl;
  gA_XS(0.856358,0.110316,0.0331726,0.791281,0.159858,0.0485461,
	153.063,103.118,478.186,75.7682,22.9003);
 

}

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

void David_gA_XS()
{
  // data for flux from Michal Krelina
  // data for probabilities from Michal Broz
  // data for UPC from David

  // y=0
  cout << " Computing case y = 0 " << endl;
  gA_XS_zero(0.82942,0.130901,0.0394531,128.076,440.8,54.5,13.0);
  gA_XS_zero(0.82942,0.130901,0.0394531,128.076,440.8*0.993,54.5*1.054,13.0*1.014);  
  //  gA_XS_zero(0.82942,0.130901,0.0394531,128.076,436.5,58.6,13.3);  
  /*
  gA_XS_zero(0.82942,0.130901,0.0394531,128.076,440.4,54.8,12.27);  
  gA_XS_zero(0.822,0.134,0.044,128.076,444.4,53.8,10.27);
  gA_XS_zero(0.807,0.141,0.052,128.076,444.4,53.8,10.27);  
  // Evg. 82.2, 13.4, 4.4 %, read off by eye from David's plots
  // starlight: 440 80.7, 14.1, 5.2%

  // |y| = 0.325
  cout << " Computing case |y| = 0.625 " << endl;
  gA_XS(0.841451,0.121721,0.0366379,0.815578,0.141448,0.0427203,
	138,118,447.2,51.73,10.47);
   
  // |y| = 0.625
  cout << " Computing case |y| = 0.625 " << endl;
  gA_XS(0.851192,0.114272,0.0343711,0.800864,0.152616,0.0462316,
	148,108,450.1,53.39,11.06);
  */

}