#include <LandauFP420.h>

Public Member Functions
	LandauFP420 ()
double	SampleFluctuations (const double momentum, const double mass, double &tmax, const double length, const double meanLoss)
	~LandauFP420 ()
Private Attributes
double	alim
double	chargeSquare
double	e0
double	e1Fluct
double	e1LogFluct
double	e2Fluct
double	e2LogFluct
double	electronDensity
double	f1Fluct
double	f2Fluct
double	ipotFluct
double	ipotLogFluct
double	minLoss
double	minNumberInteractionsBohr
int	nmaxCont1
int	nmaxCont2
double	particleMass
double	problim
double	rateFluct
double	sumalim
double	theBohrBeta2

Detailed Description

Definition at line 52 of file LandauFP420.h.

Constructor & Destructor Documentation

LandauFP420::LandauFP420 ( )

Definition at line 66 of file LandauFP420.cc.

References chargeSquare, e0, e1Fluct, e1LogFluct, e2Fluct, e2LogFluct, electronDensity, f1Fluct, f2Fluct, ipotFluct, ipotLogFluct, funct::log(), problim, rateFluct, and sumalim.

 :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();

}

LandauFP420::~LandauFP420 ( )

Definition at line 93 of file LandauFP420.cc.

{}

Member Function Documentation

double LandauFP420::SampleFluctuations	(	const double	momentum,
		const double	mass,
		double &	tmax,
		const double	length,
		const double	meanLoss
	)

Definition at line 99 of file LandauFP420.cc.

References alim, funct::C, chargeSquare, e0, e1Fluct, e1LogFluct, e2Fluct, e2LogFluct, electronDensity, f1Fluct, f2Fluct, i, ipotFluct, ipotLogFluct, gen::k, funct::log(), max(), minLoss, minNumberInteractionsBohr, nmaxCont1, nmaxCont2, p1, p2, p3, particleMass, rateFluct, mathSSE::sqrt(), sumalim, and w2.

{
//  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 = CLHEP::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(CLHEP::RandGaussQ::shoot(a3,siga)+0.5));
        }
        else
          p3 = CLHEP::RandPoisson::shoot(a3);
        //p3 = G4Poisson(a3);

        loss = p3*e0 ;

        if(p3 > 0)
          //loss += (1.-2.*G4UniformRand())*e0 ;
          loss += (1.-2.*CLHEP::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(CLHEP::RandGaussQ::shoot(a3,siga)+0.5));
        }
        else
          p3 = CLHEP::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*CLHEP::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(CLHEP::RandGaussQ::shoot(a1,siga)+0.5));
      }
      else
       p1 = CLHEP::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(CLHEP::RandGaussQ::shoot(a2,siga)+0.5));
      }
      else
        p2 = CLHEP::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.*CLHEP::RandFlat::shoot())*e2Fluct;   
      else if (loss>0.)
        loss += (1.-2.*CLHEP::RandFlat::shoot())*e1Fluct;   

      // ionisation .......................................
     if(a3 > 0.)
     {
      if(a3>alim)
      {
        siga=sqrt(a3) ;
        p3 = std::max(0,int(CLHEP::RandGaussQ::shoot(a3,siga)+0.5));
      }
      else
        p3 = CLHEP::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         = CLHEP::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 += CLHEP::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*CLHEP::RandFlat::shoot());
        }
      }        
      loss += lossc;  
     }
    } 

  return loss;
}