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

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


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

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)~xs_m, kFALSE eq (15)~xs_p
// 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);

  // 1~xs_m, 2~xs_p
  // 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 Test_gA_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)

  cout << " Computing case y = <-4.0,-2.5> " << endl;
  gA_XS(0.894622,0.0809586,0.0243349,0.625904,0.279219,0.0940511,
	163.916,13.8525,2.32794+5*(2.32794/100),0.368769-5*(0.368769/10),0.118036-0*(0.118036/10));
 

}