#include "H2FS_test.h"

using namespace std;

TKratFSCalculator::TKratFSCalculator()
{
  //pa = {
  //      {1.02947, 119.19287, 0.217256, 49.043086, 1.23644, 1.210464, 1.29972, 3.310218},
  //      {0.965069, 10.031723, 0.355287, 3.988292, 0.717569, 0.672956, 0.862026, 2.24877}
  //     };
  pa[0][0] = 1.02947;
  pa[0][1] = 119.19287;
  pa[0][2] = 0.217256;
  pa[0][3] = 49.043086;
  pa[0][4] = 1.23644;
  pa[0][5] = 1.210464;
  pa[0][6] = 1.29972;
  pa[0][7] = 3.310218;

  pa[1][0] = 0.965069;
  pa[1][1] = 10.031723;
  pa[1][2] = 0.355287;
  pa[1][3] = 3.988292;
  pa[1][4] = 0.717569;
  pa[1][5] = 0.672956;
  pa[1][6] = 0.862026;
  pa[1][7] = 2.24877;
  
  
}

TKratFSCalculator::~TKratFSCalculator()
{
}
  
  
//-----------------------------------------------------------------------------
double TKratFSCalculator::E2H(double e_in,double e_out, int iz, int ia, int icsi){

  // konwersja energii (MeV) na swiatlo (CHANNELS)
  // zawiara parametry kalibracyjne dla modulu 7
  // e_in  -> energia czastki wchodzacej do krysztalu
  // e_out -> energia pozostala po przebiciu krysztalu
  // e_out = 0 dla czastek stoTMath::Poweranych w krysztale

  // icsi = 0  -> CSI1 ( 25 mm)
  // icsi = 1  -> CSI2 (125 mm)

  if(icsi<0 || icsi>1) return 0.f;
  if(iz<=0  ||  ia<=0) return 0.f;

  double a1, a2, fct_si;

  if(icsi == 0){
    a1     = 6.285211;
    a2     = 0.316901;
    fct_si = 96.197182;
    if(iz==1){
      if     (ia==1)      a1 *= 1.077039;
      else if(ia==2)      a1 *= 1.048811;
      else if(ia==3)      a1 *= 1.033668;
    }
    else if(iz==2&&ia==3) a1 *= 1.014352;
    else			a1 *= 1.009463;
  }
  else if(icsi == 1){
    a1     = 6.285211;
    a2     = 0.316901;
    fct_si = 44.021126;
    if(iz==1){
      if     (ia==1)      a1 *= 1.079550;
      else if(ia==2)      a1 *= 1.053214;
      else if(ia==3)      a1 *= 1.038230;
    }
    else if(iz==2&&ia==3) a1 *= 0.995979;
    else			a1 *= 0.991244;
  }

  double   H;
  double   a, z;
  double   az2, fbe;

  H = 0.;

  z     = double(iz);
  a     = double(ia);

  // if(iz==1 && ia==-1) a=139.56995/931.4943; //Pion case
  // if(iz==1 && ia==-2) a=105.65841/931.4943; //Muon case
  // if(iz==1 && ia==-3) a=493.67700/931.4943; //Kaon case
        
        
  az2 = a*z*z;
  fbe = 1.;
  if(ia == 8 && iz == 4){           // Be8 = 2* He4
    az2 = az2/8.;
    fbe = 2.;
  }


  double a2az2 = a2*az2;

  e_in /= fbe;
  e_out /= fbe;

  H = (a2 != 0.) ? 
      e_in - e_out - a2az2*TMath::Log((e_in + a2az2)/(e_out + a2az2))
      :
      e_in - e_out;
      
      if(H<0) H=0;


  H *= fbe*fct_si/a1;

  return H;

}

//-----------------------------------------------------------------------------
double TKratFSCalculator::slow(Double_t fast, Int_t iz, Int_t ia){

  // parameters for CSI2 of MODULE 07
  // parametrization a'la Parlog

  if(iz<=0 || ia<=0) return 0.;
  double a0,a1,a2,a3,a4,a5,a6,a7;
  double  fbe=1.;
  if(ia == 8 && iz == 4){           // Be8 = 2* He4
    ia = 4;
    iz = 2;
    fbe = 2.;
    fast /=2.;
  }
  if(iz<=2){
    a0  = pa[0][0];
    a1  = pa[0][1];
    a2  = pa[0][2];
    a3  = pa[0][3];
    a4  = pa[0][4];
    a5  = pa[0][5];
    a6  = pa[0][6];
    a7  = pa[0][7];
  }
  else{
    a0  = pa[1][0];
    a1  = pa[1][1];
    a2  = pa[1][2];
    a3  = pa[1][3];
    a4  = pa[1][4];
    a5  = pa[1][5];
    a6  = pa[1][6];
    a7  = pa[1][7];
  }

  double a   = TMath::Power(double(ia),a6);
  double F   = TMath::Power(fast,a4);
  double z2  = TMath::Power(double(iz),a7);
  double a1az2  = a1*a*z2;

  double arg = F-a1az2*TMath::Log(1.+F/a1az2)+a2*a1az2*TMath::Log((F+a1az2)/(a3*a+a1az2));

  double S=0;
  if(arg>0) S= a0*TMath::Power(arg,1./a5);
  else      S=-a0*TMath::Power(-arg,1./a5);
  if(S<0) S=0;

  return fbe*S;
}

//-----------------------------------------------------------------------------
void TKratFSCalculator::H2FS(Double_t H, Int_t iz, Int_t ia, Double_t &Fast, Double_t &Slow){
  // parameters for CSI2 of MODULE 07
  // parametrization a'la Parlog
  // returns Fast and Slow components in channels for a given total light
  int    niter=0;
  int    iblag=0;
  bool   cond;
  double eps=1.e-6;

  double fF,fpF,dF,F,F4,F5;

  double  fbe=1.;
  if(ia == 8 && iz == 4){           // Be8 = 2* He4
    ia = 4;
    iz = 2;
    fbe = 2.;
  }
  double a0,a1,a2,a3,a4,a5,a6,a7;
  if(iz<=2){
    a0  = pa[0][0];
    a1  = pa[0][1];
    a2  = pa[0][2];
    a3  = pa[0][3];
    a4  = pa[0][4];
    a5  = pa[0][5];
    a6  = pa[0][6];
    a7  = pa[0][7];
  }
  else{
    a0  = pa[1][0];
    a1  = pa[1][1];
    a2  = pa[1][2];
    a3  = pa[1][3];
    a4  = pa[1][4];
    a5  = pa[1][5];
    a6  = pa[1][6];
    a7  = pa[1][7];
  }


  double a   = TMath::Power(double(ia),a6);
  double z2  = TMath::Power(double(iz),a7);
  double a1az2  = a1*a*z2;
  H /= fbe;
  F = 0.5*H;
  cond=1;
  niter=0;
  while(cond){
    niter++;
    if(F<=0 || H-F<=0) break;
    F4  = TMath::Power(F,a4);
    F5  = TMath::Power((H-F)/a0,a5);
    fF  = F5-F4+a1az2*TMath::Log(1.+F4/a1az2)-a2*a1az2*TMath::Log((F4+a1az2)/(a3*a+a1az2));
    fpF = -((a4*F4/F*(F4+a1az2*a2))/(F4+a1az2))-(a5*F5)/(H-F);

    dF   = -fF/fpF;
    F    = F+dF;
    cond = (TMath::Abs(dF/F) >= eps && F > 0. && niter <= 1000);
  }
  if(F<=0.) iblag=13;
  else if(niter > 1000) iblag=14;
  //if(iblag>0)printf("iblag %d\n",iblag);
  Fast = F;
  Slow = H-F;
  Fast *= fbe;
  Slow *= fbe;
  return;
}

//-----------------------------------------------------------------------------

void TKratFSCalculator::E2FS(Double_t e_in, Double_t e_out, Int_t iz, Int_t ia, Double_t &Fast, Double_t &Slow){
  // parameters for CSI2 of MODULE 07
  // parametrization a'la ParTMath::Log
  // returns Fast and Slow components in channels for a given e_in e_out in MeV

  double Fout, Sout, Hout;
  Fout=Sout=Hout=0.;
  if(e_out>0){
    Hout = E2H(e_out, 0., iz, ia, 1);
    H2FS(Hout, iz, ia, Fout, Sout);
  }
  double Ftot, Stot, Htot;
  Ftot=Stot=Htot=0.;
  Htot=E2H(e_in, 0., iz, ia, 1);
  H2FS(Htot, iz, ia, Ftot, Stot);
  Fast = Ftot-Fout;
  Slow = Stot-Sout;
  return;
}



//______________________________________________________________________
void TKratFSCalculator::GetAZformPdgCode( Int_t A_pdg_code, Int_t &Ref_z, Int_t &Ref_a, Int_t debug )
{
  /// Some information: http://www3.nd.edu/~avillano/geant4/geant4_pid.html
  
  Ref_a = 0;
  Ref_z = 0;
  
  ///PDG Code ion code: 1000491150
  
  Int_t M = 1000000000;
  
  if (A_pdg_code == 2212 ){ ///proton
    Ref_a = 1;
    Ref_z = 1;
    
  }
  if (A_pdg_code == 2112 ){
    Ref_a = 1;
    Ref_z = 0;
    
  }
  if( A_pdg_code > M) {
    Ref_z = ((Int_t) (A_pdg_code / 10000)) % 1000;
    Ref_a = ((Int_t) (A_pdg_code / 10)) % 1000;
    
      
  }
  if( debug > 0 ){
    cout << "[GetAZformPdgCode:] PDG code=" << A_pdg_code;
    cout << "\ta=" << Ref_a;
    cout << "\tz=" << Ref_z;
    cout << endl;
  }
  return;
}


void TKratFSCalculator::pdg_test( void )
{
  
  int a = 0;
  int z = 0;
  int pdg_code = 100010010; /// 1000491150
  
  GetAZformPdgCode (pdg_code, z, a);
  cout << "PDG code=" << pdg_code;
  cout << "\ta=" << a;
  cout << "\tz=" << z;
  cout << endl;
  
  pdg_code = 2112;
  
  GetAZformPdgCode (pdg_code, z, a);
    cout << "PDG code=" << pdg_code;
  cout << "\ta=" << a;
  cout << "\tz=" << z;
  cout << endl;
  
  pdg_code = 2212;
  GetAZformPdgCode (pdg_code, z, a);
    cout << "PDG code=" << pdg_code;
  cout << "\ta=" << a;
  cout << "\tz=" << z;
  cout << endl;
  
  pdg_code = 21120;
  GetAZformPdgCode (pdg_code, z, a);
    cout << "PDG code=" << pdg_code;
  cout << "\ta=" << a;
  cout << "\tz=" << z;
  cout << endl;
  
    pdg_code = 1000491150;
  GetAZformPdgCode (pdg_code, z, a);
    cout << "PDG code=" << pdg_code;
  cout << "\ta=" << a;
  cout << "\tz=" << z;
  cout << endl;
  
   pdg_code = 1009986321;
  GetAZformPdgCode (pdg_code, z, a );
    cout << "PDG code=" << pdg_code;
  cout << "\ta=" << a;
  cout << "\tz=" << z;
  cout << endl;
  
  return;
  
}



ClassImp( TKratFSCalculator )


/// End of file