#include <SiG4UniversalFluctuation.h>

Public Member Functions
double	SampleFluctuations (const double momentum, const double mass, double &tmax, const double length, const double meanLoss, CLHEP::HepRandomEngine *)

	SiG4UniversalFluctuation ()

	~SiG4UniversalFluctuation ()

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

double	nmaxCont1

double	nmaxCont2

double	particleMass

double	problim

double	rateFluct

double	sumalim

double	theBohrBeta2

Detailed Description

Definition at line 63 of file SiG4UniversalFluctuation.h.

Constructor & Destructor Documentation

SiG4UniversalFluctuation::SiG4UniversalFluctuation ( )

Definition at line 75 of file SiG4UniversalFluctuation.cc.

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

   :minNumberInteractionsBohr(10.0),
    theBohrBeta2(50.0*keV/proton_mass_c2),
    minLoss(10.*eV),
    problim(5.e-3),  
    alim(10.),
    nmaxCont1(4.),
    nmaxCont2(16.)
 {
   sumalim = -log(problim);
   //lastMaterial = 0;
   
   // Add these definitions d.k.
   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();
   
   //cout << " init new fluct +++++++++++++++++++++++++++++++++++++++++"<<endl;
 }

SiG4UniversalFluctuation::~SiG4UniversalFluctuation ( )

Definition at line 110 of file SiG4UniversalFluctuation.cc.

111 {

112 }

Member Function Documentation

double SiG4UniversalFluctuation::SampleFluctuations	(	const double	momentum,
		const double	mass,
		double &	tmax,
		const double	length,
		const double	meanLoss,
		CLHEP::HepRandomEngine *	engine
	)

Definition at line 115 of file SiG4UniversalFluctuation.cc.

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

 {
 // Calculate actual loss from the mean loss.
 // The model used to get the fluctuations is essentially the same
 // as in Glandz in Geant3 (Cern program library W5013, phys332).
 // L. Urban et al. NIM A362, p.416 (1995) and Geant4 Physics Reference Manual
 
   // shortcut for very very small loss (out of validity of the model)
   //
   if (meanLoss < minLoss) return meanLoss;
 
   //if(!particle) InitialiseMe(dp->GetDefinition());
   //G4double tau   = dp->GetKineticEnergy()/particleMass;
   //G4double gam   = tau + 1.0;
   //G4double gam2  = gam*gam;
   //G4double beta2 = tau*(tau + 2.0)/gam2;
 
   particleMass   = mass; // dp->GetMass();
   double gam2   = (momentum*momentum)/(particleMass*particleMass) + 1.0;
   double beta2 = 1.0 - 1.0/gam2;
   double gam = sqrt(gam2);
 
   double loss(0.), siga(0.);
   
   // Gaussian regime
   // for heavy particles only and conditions
   // for Gauusian fluct. has been changed 
   //
   if ((particleMass > electron_mass_c2) &&
       (meanLoss >= minNumberInteractionsBohr*tmax))
   {
     double massrate = electron_mass_c2/particleMass ;
     double tmaxkine = 2.*electron_mass_c2*beta2*gam2/
       (1.+massrate*(2.*gam+massrate)) ;
     if (tmaxkine <= 2.*tmax)   
     {
       //electronDensity = material->GetElectronDensity();
       siga  = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
                                 * electronDensity * chargeSquare;
       siga = sqrt(siga);
       double twomeanLoss = meanLoss + meanLoss;
       if (twomeanLoss < siga) {
         double x;
         do {
           loss = twomeanLoss*CLHEP::RandFlat::shoot(engine);
           x = (loss - meanLoss)/siga;
         } while (1.0 - 0.5*x*x < CLHEP::RandFlat::shoot(engine));
       } else {
         do {
           loss = CLHEP::RandGaussQ::shoot(engine, meanLoss, siga);
         } while (loss < 0. || loss > twomeanLoss);
       }
       return loss;
     }
   }
 
   // Glandz regime : initialisation
   //
 //   if (material != lastMaterial) {
 //     f1Fluct      = material->GetIonisation()->GetF1fluct();
 //     f2Fluct      = material->GetIonisation()->GetF2fluct();
 //     e1Fluct      = material->GetIonisation()->GetEnergy1fluct();
 //     e2Fluct      = material->GetIonisation()->GetEnergy2fluct();
 //     e1LogFluct   = material->GetIonisation()->GetLogEnergy1fluct();
 //     e2LogFluct   = material->GetIonisation()->GetLogEnergy2fluct();
 //     rateFluct    = material->GetIonisation()->GetRateionexcfluct();
 //     ipotFluct    = material->GetIonisation()->GetMeanExcitationEnergy();
 //     ipotLogFluct = material->GetIonisation()->GetLogMeanExcEnergy();
 //     lastMaterial = material;
 //   }
 
   double a1 = 0. , a2 = 0., a3 = 0. ;
   double p1,p2,p3;
   double rate = rateFluct ;
 
   double w1 = tmax/ipotFluct;
   double w2 = vdt::fast_log(2.*electron_mass_c2*beta2*gam2)-beta2;
 
   if(w2 > ipotLogFluct)
   {
     double C = meanLoss*(1.-rateFluct)/(w2-ipotLogFluct);
     a1 = C*f1Fluct*(w2-e1LogFluct)/e1Fluct;
     a2 = C*f2Fluct*(w2-e2LogFluct)/e2Fluct;
     if(a2 < 0.)
     {
       a1 = 0. ;
       a2 = 0. ;
       rate = 1. ;  
     }
   }
   else
   {
     rate = 1. ;
   }
 
   // added
   if(tmax > ipotFluct) {
     a3 = rate*meanLoss*(tmax-ipotFluct)/(ipotFluct*tmax*vdt::fast_log(w1));
   }
   double suma = a1+a2+a3;
   
   // Glandz regime
   //
   if (suma > sumalim)
   {
     p1 = 0., p2 = 0 ;
     if((a1+a2) > 0.)
     {
       // excitation type 1
       if (a1>alim) {
         siga=sqrt(a1) ;
         p1 = max(0., CLHEP::RandGaussQ::shoot(engine, a1, siga) + 0.5);
       } else {
         CLHEP::RandPoissonQ randPoissonQ(*engine, a1);
         p1 = double(randPoissonQ.fire());
       }
     
       // excitation type 2
       if (a2>alim) {
         siga=sqrt(a2) ;
         p2 = max(0., CLHEP::RandGaussQ::shoot(engine, a2, siga) + 0.5);
       } else {
         CLHEP::RandPoissonQ randPoissonQ(*engine, a2);
         p2 = double(randPoissonQ.fire());
       }
     
       loss = p1*e1Fluct+p2*e2Fluct;
  
       // smearing to avoid unphysical peaks
       if (p2 > 0.)
         loss += (1.-2.*CLHEP::RandFlat::shoot(engine))*e2Fluct;
       else if (loss>0.)
         loss += (1.-2.*CLHEP::RandFlat::shoot(engine))*e1Fluct;
       if (loss < 0.) loss = 0.0;
     }
 
     // ionisation
     if (a3 > 0.) {
       if (a3>alim) {
         siga=sqrt(a3) ;
         p3 = max(0., CLHEP::RandGaussQ::shoot(engine, a3, siga) + 0.5);
       } else {
         CLHEP::RandPoissonQ randPoissonQ(*engine, a3);
         p3 = double(randPoissonQ.fire());
       }
       double lossc = 0.;
       if (p3 > 0) {
         double na = 0.; 
         double alfa = 1.;
         if (p3 > nmaxCont2) {
           double rfac   = p3/(nmaxCont2+p3);
           double namean = p3*rfac;
           double sa     = nmaxCont1*rfac;
           na              = CLHEP::RandGaussQ::shoot(engine, namean, sa);
           if (na > 0.) {
             alfa   = w1*(nmaxCont2+p3)/(w1*nmaxCont2+p3);
             double alfa1  = alfa*vdt::fast_log(alfa)/(alfa-1.);
             double ea     = na*ipotFluct*alfa1;
             double sea    = ipotFluct*sqrt(na*(alfa-alfa1*alfa1));
             lossc += CLHEP::RandGaussQ::shoot(engine, ea, sea);
           }
         }
 
         if (p3 > na) {
           w2 = alfa*ipotFluct;
           double w  = (tmax-w2)/tmax;
           int    nb = int(p3-na);
           for (int k=0; k<nb; k++) lossc += w2/(1.-w*CLHEP::RandFlat::shoot(engine));
         }
       }        
       loss += lossc;  
     }
     return loss;
   }
   
   // suma < sumalim;  very small energy loss;  
   //
   //double e0 = material->GetIonisation()->GetEnergy0fluct();
 
   a3 = meanLoss*(tmax-e0)/(tmax*e0*vdt::fast_log(tmax/e0));
   if (a3 > alim)
   {
     siga=sqrt(a3);
     p3 = max(0., CLHEP::RandGaussQ::shoot(engine, a3, siga) + 0.5);
   } else {
     CLHEP::RandPoissonQ randPoissonQ(*engine, a3);
     p3 = double(randPoissonQ.fire());
   }
   if (p3 > 0.) {
     double w = (tmax-e0)/tmax;
     double corrfac = 1.;
     if (p3 > nmaxCont2) {
       corrfac = p3/nmaxCont2;
       p3 = nmaxCont2;
     } 
     int ip3 = (int)p3;
     for (int i=0; i<ip3; i++) loss += 1./(1.-w*CLHEP::RandFlat::shoot(engine));
     loss *= e0*corrfac;
     // smearing for losses near to e0
     if(p3 <= 2.)
       loss += e0*(1.-2.*CLHEP::RandFlat::shoot(engine)) ;
    }
     
    return loss;
 }