#include <SiG4UniversalFluctuation.h>

Public Member Functions
double	SampleFluctuations (const double momentum, const double mass, double &tmax, const double length, const double meanLoss)
	SiG4UniversalFluctuation (CLHEP::HepRandomEngine &)
	~SiG4UniversalFluctuation ()
Private Attributes
double	alim
double	chargeSquare
double	e0
double	e1Fluct
double	e1LogFluct
double	e2Fluct
double	e2LogFluct
double	electronDensity
double	f1Fluct
double	f2Fluct
CLHEP::RandFlat *	flatDistribution
CLHEP::RandGaussQ *	gaussQDistribution
double	ipotFluct
double	ipotLogFluct
double	minLoss
double	minNumberInteractionsBohr
double	nmaxCont1
double	nmaxCont2
double	particleMass
CLHEP::RandPoissonQ *	poissonQDistribution
double	problim
double	rateFluct
CLHEP::HepRandomEngine &	rndEngine
double	sumalim
double	theBohrBeta2

Detailed Description

Definition at line 67 of file SiG4UniversalFluctuation.h.

Constructor & Destructor Documentation

SiG4UniversalFluctuation::SiG4UniversalFluctuation ( CLHEP::HepRandomEngine & eng )

Definition at line 76 of file SiG4UniversalFluctuation.cc.

References chargeSquare, e0, e1Fluct, e1LogFluct, e2Fluct, e2LogFluct, electronDensity, f1Fluct, f2Fluct, flatDistribution, gaussQDistribution, ipotFluct, ipotLogFluct, create_public_lumi_plots::log, poissonQDistribution, problim, rateFluct, rndEngine, and sumalim.

  :rndEngine(eng),
   gaussQDistribution(0),
   poissonQDistribution(0),
   flatDistribution(0),
   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();
  
  gaussQDistribution = new CLHEP::RandGaussQ(rndEngine);
  poissonQDistribution = new CLHEP::RandPoissonQ(rndEngine);
  flatDistribution = new CLHEP::RandFlat(rndEngine);

  //cout << " init new fluct +++++++++++++++++++++++++++++++++++++++++"<<endl;
}

SiG4UniversalFluctuation::~SiG4UniversalFluctuation ( )

Definition at line 119 of file SiG4UniversalFluctuation.cc.

References flatDistribution, gaussQDistribution, and poissonQDistribution.

{
  delete gaussQDistribution;
  delete poissonQDistribution;
  delete flatDistribution;

}

Member Function Documentation

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

Definition at line 128 of file SiG4UniversalFluctuation.cc.

References alim, funct::C, chargeSquare, e0, e1Fluct, e1LogFluct, e2Fluct, e2LogFluct, electronDensity, f1Fluct, f2Fluct, flatDistribution, gam, gaussQDistribution, i, ipotFluct, ipotLogFluct, gen::k, max(), minLoss, minNumberInteractionsBohr, nmaxCont1, nmaxCont2, p1, p2, p3, particleMass, poissonQDistribution, 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*flatDistribution->fire();
          x = (loss - meanLoss)/siga;
        } while (1.0 - 0.5*x*x < flatDistribution->fire());
      } else {
        do {
          loss = gaussQDistribution->fire(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.,gaussQDistribution->fire(a1,siga)+0.5);
      } else {
        p1 = double(poissonQDistribution->fire(a1));
      }
    
      // excitation type 2
      if (a2>alim) {
        siga=sqrt(a2) ;
        p2 = max(0.,gaussQDistribution->fire(a2,siga)+0.5);
      } else {
        p2 = double(poissonQDistribution->fire(a2));
      }
    
      loss = p1*e1Fluct+p2*e2Fluct;
 
      // smearing to avoid unphysical peaks
      if (p2 > 0.)
        loss += (1.-2.*flatDistribution->fire())*e2Fluct;
      else if (loss>0.)
        loss += (1.-2.*flatDistribution->fire())*e1Fluct;   
      if (loss < 0.) loss = 0.0;
    }

    // ionisation
    if (a3 > 0.) {
      if (a3>alim) {
        siga=sqrt(a3) ;
        p3 = max(0.,gaussQDistribution->fire(a3,siga)+0.5);
      } else {
        p3 = double(poissonQDistribution->fire(a3));
      }
      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              = gaussQDistribution->fire(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 += gaussQDistribution->fire(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*flatDistribution->fire());
        }
      }        
      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.,gaussQDistribution->fire(a3,siga)+0.5);
  } else {
    p3 = double(poissonQDistribution->fire(a3));
  }
  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*flatDistribution->fire());
    loss *= e0*corrfac;
    // smearing for losses near to e0
    if(p3 <= 2.)
      loss += e0*(1.-2.*flatDistribution->fire()) ;
   }
    
   return loss;
}