#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, create_public_lumi_plots::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.

94 {}

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, create_public_lumi_plots::log, max(), minLoss, minNumberInteractionsBohr, nmaxCont1, nmaxCont2, p1, p2, p3, particleMass, rateFluct, mathSSE::sqrt(), sumalim, w(), 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;
 }