//
// ********************************************************************
// * DISCLAIMER                                                       *
// *                                                                  *
// * The following disclaimer summarizes all the specific disclaimers *
// * of contributors to this software. The specific disclaimers,which *
// * govern, are listed with their locations in:                      *
// *   http://cern.ch/geant4/license                                  *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.                                                             *
// *                                                                  *
// * This  code  implementation is the  intellectual property  of the *
// * GEANT4 collaboration.                                            *
// * By copying,  distributing  or modifying the Program (or any work *
// * based  on  the Program)  you indicate  your  acceptance of  this *
// * statement, and all its terms.                                    *
// ********************************************************************
//
// -------------------------------------------------------------------
//
// GEANT4 Class file
//
//
// File name:     G4UniversalFluctuation
//
// Author:        Vladimir Ivanchenko 
// 
// Creation date: 03.01.2002
//
// Modifications: 
//
// 28-12-02 add method Dispersion (V.Ivanchenko)
// 07-02-03 change signature (V.Ivanchenko)
// 13-02-03 Add name (V.Ivanchenko)
// Modified for standalone use in ORCA. d.k. 6/04
//
// -------------------------------------------------------------------
 

#include "MyG4UniversalFluctuationForSi.h"
#include "CLHEP/Units/SystemOfUnits.h"
#include "CLHEP/Units/PhysicalConstants.h"
#include "CLHEP/Random/RandGaussQ.h"
#include "CLHEP/Random/RandPoisson.h"
#include "CLHEP/Random/RandFlat.h"
#include <cmath>
#include <iostream>
using std::cout;
using std::endl;
namespace CLHEP{}    // declare namespace CLHEP for backward compatibility
using namespace CLHEP ;


//#include "Randomize.hh"
//#include "G4Poisson.hh"
//#include "G4Step.hh"
//#include "G4Material.hh"
//#include "G4DynamicParticle.hh"
//#include "G4ParticleDefinition.hh"

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
// The constructor setups various constants pluc eloss parameters
// for silicon.   
MyG4UniversalFluctuationForSi::MyG4UniversalFluctuationForSi()
 :minNumberInteractionsBohr(10.0),
  theBohrBeta2(50.0*keV/proton_mass_c2),
  minLoss(0.000001*eV),
  problim(0.01),
  alim(10.),
  nmaxCont1(4),
  nmaxCont2(16)
{
  sumalim = -log(problim);

  chargeSquare   = 1.;  //Assume all particles have charge 1
  // Taken from Geant4 printout, HARDWIRED for Silicon.
  ipotFluct = 0.0001736; //material->GetIonisation()->GetMeanExcitationEnergy();
  electronDensity = 6.797E+20; // material->GetElectronDensity();
  f1Fluct      = 0.8571; // material->GetIonisation()->GetF1fluct();
  f2Fluct      = 0.1429; //material->GetIonisation()->GetF2fluct();
  e1Fluct      = 0.000116;// material->GetIonisation()->GetEnergy1fluct();
  e2Fluct      = 0.00196; //material->GetIonisation()->GetEnergy2fluct();
  e1LogFluct   = -9.063;  //material->GetIonisation()->GetLogEnergy1fluct();
  e2LogFluct   = -6.235;  //material->GetIonisation()->GetLogEnergy2fluct();
  rateFluct    = 0.4;     //material->GetIonisation()->GetRateionexcfluct();
  ipotLogFluct = -8.659;  //material->GetIonisation()->GetLogMeanExcEnergy();
  e0 = 1.E-5;             //material->GetIonisation()->GetEnergy0fluct();

}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
MyG4UniversalFluctuationForSi::~MyG4UniversalFluctuationForSi()
{}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
// The main dedx fluctuation routine.
// Arguments: momentum in MeV/c, mass in MeV, delta ray cut (tmax) in
// MeV, silicon thickness in mm, mean eloss in MeV. 
double MyG4UniversalFluctuationForSi::SampleFluctuations(const double momentum,
                                                         const double mass,
                                                         double& tmax,
                                                         const double length,
                                                         const double meanLoss)
{
//  calculate actual loss from the mean loss
//  The model used to get the fluctuation is essentially the same
// as in Glandz in Geant3.
  
  // shortcut for very very small loss 
  if(meanLoss < minLoss) return meanLoss;

  //if(dp->GetDefinition() != particle) {
  particleMass   = mass; // dp->GetMass();
  //G4double q     = dp->GetCharge(); 
  //chargeSquare   = q*q;
  //}

  //double gam   = (dp->GetKineticEnergy())/particleMass + 1.0;
  //double gam2  = gam*gam; 
  double gam2   = (momentum*momentum)/(particleMass*particleMass) + 1.0;
  double beta2 = 1.0 - 1.0/gam2;
 
  // Validity range for delta electron cross section
  double loss, siga;
  // Gaussian fluctuation 
  if(meanLoss >= minNumberInteractionsBohr*tmax || 
     tmax <= ipotFluct*minNumberInteractionsBohr)
  {
    siga  = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length 
                              * electronDensity * chargeSquare ;
    siga = sqrt(siga);
    do {
     //loss = G4RandGauss::shoot(meanLoss,siga);
     loss = RandGaussQ::shoot(meanLoss,siga);
    } while (loss < 0. || loss > 2.*meanLoss);

    return loss;
  }

  // Non Gaussian fluctuation 
  double suma,w1,w2,C,lossc,w;
  double a1,a2,a3;
  int p1,p2,p3;
  int nb;
  double corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
  double dp3;

  w1 = tmax/ipotFluct;
  w2 = log(2.*electron_mass_c2*(gam2 - 1.0));

  C = meanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);

  a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
  a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
  a3 = rateFluct*meanLoss*(tmax-ipotFluct)/(ipotFluct*tmax*log(w1));
  if(a1 < 0.) a1 = 0.;
  if(a2 < 0.) a2 = 0.;
  if(a3 < 0.) a3 = 0.;

  suma = a1+a2+a3;

  loss = 0. ;

  if(suma < sumalim)             // very small Step
    {
      //e0 = material->GetIonisation()->GetEnergy0fluct();//Hardwired in const
      if(tmax == ipotFluct)
      {
        a3 = meanLoss/e0;

        if(a3>alim)
        {
          siga=sqrt(a3) ;
          //p3 = G4std::max(0,int(G4RandGauss::shoot(a3,siga)+0.5));
          p3 = std::max(0,int(RandGaussQ::shoot(a3,siga)+0.5));
        }
        else
          p3 = RandPoisson::shoot(a3);
        //p3 = G4Poisson(a3);

        loss = p3*e0 ;

        if(p3 > 0)
          //loss += (1.-2.*G4UniformRand())*e0 ;
          loss += (1.-2.*RandFlat::shoot())*e0 ;

      }
      else
      {
        tmax = tmax-ipotFluct+e0 ;
        a3 = meanLoss*(tmax-e0)/(tmax*e0*log(tmax/e0));

        if(a3>alim)
        {
          siga=sqrt(a3) ;
          //p3 = G4std::max(0,int(G4RandGauss::shoot(a3,siga)+0.5));
          p3 = std::max(0,int(RandGaussQ::shoot(a3,siga)+0.5));
        }
        else
          p3 = RandPoisson::shoot(a3);
        //p3 = G4Poisson(a3);

        if(p3 > 0)
        {
          w = (tmax-e0)/tmax ;
          if(p3 > nmaxCont2)
          {
            dp3 = float(p3) ;
            corrfac = dp3/float(nmaxCont2) ;
            p3 = nmaxCont2 ;
          }
          else
            corrfac = 1. ;

          //for(int i=0; i<p3; i++) loss += 1./(1.-w*G4UniformRand()) ;
          for(int i=0; i<p3; i++) loss += 1./(1.-w*RandFlat::shoot()) ;
          loss *= e0*corrfac ;  
        }        
      }
    }
    
  else                              // not so small Step
    {
      // excitation type 1
      if(a1>alim)
      {
        siga=sqrt(a1) ;
        //p1 = std::max(0,int(G4RandGauss::shoot(a1,siga)+0.5));
        p1 = std::max(0,int(RandGaussQ::shoot(a1,siga)+0.5));
      }
      else
       p1 = RandPoisson::shoot(a1);
      //p1 = G4Poisson(a1);

      // excitation type 2
      if(a2>alim)
      {
        siga=sqrt(a2) ;
        //p2 = std::max(0,int(G4RandGauss::shoot(a2,siga)+0.5));
        p2 = std::max(0,int(RandGaussQ::shoot(a2,siga)+0.5));
      }
      else
        p2 = RandPoisson::shoot(a2);
      //p2 = G4Poisson(a2);

      loss = p1*e1Fluct+p2*e2Fluct;
 
      // smearing to avoid unphysical peaks
      if(p2 > 0)
        //loss += (1.-2.*G4UniformRand())*e2Fluct;   
        loss += (1.-2.*RandFlat::shoot())*e2Fluct;   
      else if (loss>0.)
        loss += (1.-2.*RandFlat::shoot())*e1Fluct;   

      // ionisation .......................................
     if(a3 > 0.)
     {
      if(a3>alim)
      {
        siga=sqrt(a3) ;
        p3 = std::max(0,int(RandGaussQ::shoot(a3,siga)+0.5));
      }
      else
        p3 = RandPoisson::shoot(a3);

      lossc = 0.;
      if(p3 > 0)
      {
        na = 0.; 
        alfa = 1.;
        if (p3 > nmaxCont2)
        {
          dp3        = float(p3);
          rfac       = dp3/(float(nmaxCont2)+dp3);
          namean     = float(p3)*rfac;
          sa         = float(nmaxCont1)*rfac;
          na         = RandGaussQ::shoot(namean,sa);
          if (na > 0.)
          {
            alfa   = w1*float(nmaxCont2+p3)/
                    (w1*float(nmaxCont2)+float(p3));
            alfa1  = alfa*log(alfa)/(alfa-1.);
            ea     = na*ipotFluct*alfa1;
            sea    = ipotFluct*sqrt(na*(alfa-alfa1*alfa1));
            lossc += RandGaussQ::shoot(ea,sea);
          }
        }

        nb = int(float(p3)-na);
        if (nb > 0)
        {
          w2 = alfa*ipotFluct;
          w  = (tmax-w2)/tmax;      
          for (int k=0; k<nb; k++) lossc += w2/(1.-w*RandFlat::shoot());
        }
      }        
      loss += lossc;  
     }
    } 

  return loss;
}