G4ChannelingFastSimCrystalData Class Reference

#include <G4ChannelingFastSimCrystalData.hh>

Inheritance diagram for G4ChannelingFastSimCrystalData:

Public Member Functions
	G4ChannelingFastSimCrystalData ()=default
	~G4ChannelingFastSimCrystalData ()=default
void	SetMaterialProperties (const G4Material *crystal, const G4String &lattice, const G4String &filePath)
G4ThreeVector	CoordinatesFromBoxToLattice (const G4ThreeVector &pos0)
G4ThreeVector	CoordinatesFromLatticeToBox (const G4ThreeVector &pos)
G4ThreeVector	ChannelChange (G4double &x, G4double &y, G4double &z)
	change the channel if necessary, recalculate x o y
G4double	AngleXFromBoxToLattice (G4double tx, G4double z)
G4double	AngleXFromLatticeToBox (G4double tx, G4double z)
G4double	AngleXShift (G4double z)
	auxialiary function to transform the horizontal angle
G4double	GetChannelWidthX ()
	get channel width in x and y
G4double	GetChannelWidthY ()
Public Member Functions inherited from G4VChannelingFastSimCrystalData
	G4VChannelingFastSimCrystalData ()
virtual	~G4VChannelingFastSimCrystalData ()
G4double	Ex (G4double x, G4double y)
	electric fields produced by crystal lattice
G4double	Ey (G4double x, G4double y)
G4double	ElectronDensity (G4double x, G4double y)
	electron density function
G4double	MinIonizationEnergy (G4double x, G4double y)
	minimum energy of ionization function
G4double	NuclearDensity (G4double x, G4double y, G4int ielement)
	nuclear density function (normalized to average nuclear density)
G4double	GetLindhardAngle (G4double etotal, G4double mass, G4double charge)
	Calculate the value of the Lindhard angle (!!! the value for a straight crystal).
G4double	GetLindhardAngle ()
	Calculate the value of the Lindhard angle (!!! the value for a straight crystal).
G4double	GetSimulationStep (G4double tx, G4double ty)
G4double	GetMaxSimulationStep (G4double etotal, G4double mass, G4double charge)
	Calculate maximal simulation step (standard value for channeling particles).
G4double	GetBeta ()
	get particle velocity/c
G4int	GetNelements ()
G4int	GetModel ()
G4double	GetBendingAngle ()
G4double	GetMiscutAngle ()
G4double	GetCurv (G4double z)
G4double	GetCUx (G4double z)
	get crystalline undulator wave function
G4double	GetCUtetax (G4double z)
	get crystalline undulator wave 1st derivative function
void	SetGeometryParameters (const G4LogicalVolume *crystallogic)
	set geometry parameters from current logical volume
void	SetBendingAngle (G4double tetab, const G4LogicalVolume *crystallogic)
void	SetMiscutAngle (G4double tetam, const G4LogicalVolume *crystallogic)
void	SetCrystallineUndulatorParameters (G4double amplitude, G4double period, G4double phase, const G4LogicalVolume *crystallogic)
void	SetCrystallineUndulatorParameters (const G4LogicalVolume *crystallogic, const G4String &filename="CUgeometry.dat")
	set importing geometry of a crystalline undulator from file for specific volume
void	SetCUParameters (const G4ThreeVector &amplitudePeriodPhase, const G4LogicalVolume *crystallogic)
void	SetParticleProperties (G4double etotal, G4double mp, G4double charge, const G4String &particleName)
G4double	GetCorrectionZ ()
G4ThreeVector	CoulombAtomicScattering (G4double effectiveStep, G4double step, G4int ielement)
	multiple and single scattering on screened potential
G4ThreeVector	CoulombElectronScattering (G4double eMinIonization, G4double electronDensity, G4double step)
	multiple and single scattering on electrons
G4double	IonizationLosses (G4double dz, G4int ielement)
	ionization losses
void	SetVerbosity (G4int ver)

Additional Inherited Members
Protected Attributes inherited from G4VChannelingFastSimCrystalData
G4ChannelingFastSimInterpolation *	fElectricFieldX {nullptr}
	classes containing interpolation coefficients
G4ChannelingFastSimInterpolation *	fElectricFieldY {nullptr}
G4ChannelingFastSimInterpolation *	fElectronDensity {nullptr}
G4ChannelingFastSimInterpolation *	fMinIonizationEnergy {nullptr}
std::vector< G4ChannelingFastSimInterpolation * >	fNucleiDensity
G4ThreeVector	fHalfDimBoundingBox
	values related to the crystal geometry
G4int	fBent =0
G4double	fBendingAngle =0.
G4double	fBendingR = 0.
G4double	fBending2R =0.
G4double	fBendingRsquare =0.
G4double	fCurv =0.
G4double	fMiscutAngle = 0.
G4double	fCosMiscutAngle =1.
G4double	fSinMiscutAngle =0.
G4double	fCorrectionZ = 1.
G4bool	fCU = false
G4double	fCUAmplitude =0.
G4double	fCUK =0.
G4double	fCUPhase =0.
G4double	fCUAmplitudeK =0.
G4double	fCUK2 =0.
G4int	fNelements =1
	values related to the crystal lattice
G4int	iModel =1
G4double	fVmax =0
G4double	fVmax2 =0
G4double	fVMinCrystal =0
G4double	fChangeStep =0
std::vector< G4double >	fI0
std::vector< G4double >	fRF
std::vector< G4double >	fTeta10
	angles necessary for multiple and single coulomb scattering
std::vector< G4double >	fTetamax0
std::vector< G4double >	fTetamax2
std::vector< G4double >	fTetamax12
std::vector< G4double >	fTeta12
std::vector< G4double >	fK20
	coefficients necessary for multiple and single coulomb scattering
std::vector< G4double >	fK2
std::vector< G4double >	fK40
G4double	fK30 =0
G4double	fK3 =0
std::vector< G4double >	fKD
std::vector< G4double >	fLogPlasmaEdI0
std::vector< G4double >	fPu11
	coefficients for multiple scattering suppression
std::vector< G4double >	fPzu11
std::vector< G4double >	fBB
std::vector< G4double >	fE1XBbb
std::vector< G4double >	fBBDEXP
G4int	fVerbosity = 1

Detailed Description

Definition at line 50 of file G4ChannelingFastSimCrystalData.hh.

Constructor & Destructor Documentation

G4ChannelingFastSimCrystalData()

G4ChannelingFastSimCrystalData::G4ChannelingFastSimCrystalData ( )

default

~G4ChannelingFastSimCrystalData()

G4ChannelingFastSimCrystalData::~G4ChannelingFastSimCrystalData ( )

default

Member Function Documentation

AngleXFromBoxToLattice()

G4double G4ChannelingFastSimCrystalData::AngleXFromBoxToLattice ( G4double tx, G4double z )

inlinevirtual

calculate the horizontal angle in the co-rotating reference system within a channel (periodic cell) (connected with crystal planes/axes either bent or straight)

Implements G4VChannelingFastSimCrystalData.

Definition at line 78 of file G4ChannelingFastSimCrystalData.hh.

    {return tx-AngleXShift(z)-GetCUtetax(z);}

AngleXFromLatticeToBox()

G4double G4ChannelingFastSimCrystalData::AngleXFromLatticeToBox ( G4double tx, G4double z )

inlinevirtual

calculate the horizontal angle in the Box reference system (connected with the bounding box of the volume)

Implements G4VChannelingFastSimCrystalData.

Definition at line 83 of file G4ChannelingFastSimCrystalData.hh.

    {return tx+AngleXShift(z)+GetCUtetax(z);}

AngleXShift()

G4double G4ChannelingFastSimCrystalData::AngleXShift ( G4double z )

inlinevirtual

auxialiary function to transform the horizontal angle

Implements G4VChannelingFastSimCrystalData.

Definition at line 87 of file G4ChannelingFastSimCrystalData.hh.

{return fMiscutAngle + z*fCurv;}

Referenced by AngleXFromBoxToLattice(), and AngleXFromLatticeToBox().

ChannelChange()

G4ThreeVector G4ChannelingFastSimCrystalData::ChannelChange ( G4double & x, G4double & y, G4double & z )

virtual

change the channel if necessary, recalculate x o y

Implements G4VChannelingFastSimCrystalData.

Definition at line 399 of file G4ChannelingFastSimCrystalData.cc.

{
 
    //test of enter in other channel
    if (x<0)
    {
        fNChannelx-=1;
        x+=fDx; //enter in other channel
        //correction of the longitudinal coordinate
        if (fBent) {fCorrectionZ = fBendingR/(fBendingR-fNChannelx*fDx);}
    }
    else if (x>=fDx)
    {
        fNChannelx+=1;
        x-=fDx; //enter in other channel
        //correction of the longitudinal coordinate
        if (fBent) {fCorrectionZ = fBendingR/(fBendingR-fNChannelx*fDx);}
    }
 
    //test of enter in other channel
    if (y<0)
    {
        fNChannely-=1;
        y+=fDy; //enter in other channel
    }
    else if (y>=fDy)
    {
        fNChannely+=1;
        y-=fDy; //enter in other channel
    }
 
    return G4ThreeVector(x,y,z);
}

CoordinatesFromBoxToLattice()

G4ThreeVector G4ChannelingFastSimCrystalData::CoordinatesFromBoxToLattice ( const G4ThreeVector & pos0 )

virtual

calculate the coordinates in the co-rotating reference system within a channel (periodic cell) (connected with crystal planes/axes either bent or straight)

Implements G4VChannelingFastSimCrystalData.

Definition at line 316 of file G4ChannelingFastSimCrystalData.cc.

{
   G4double x0=pos0.x(),y0=pos0.y(),z0=pos0.z();
   z0+=fHalfDimBoundingBox.z();
   G4double x,y,z;
 
   if (fBent)
   {
       // for bent crystal
       G4double rsqrt = std::sqrt(fBendingRsquare -
                                  fBending2R*(x0*fCosMiscutAngle - z0*fSinMiscutAngle) +
                                  x0*x0 + z0*z0);
       //transform to co-rotating reference system connected with "central plane/axis"
       x = fBendingR - rsqrt;
       y = y0;
       z = fBendingR*std::asin((z0*fCosMiscutAngle + x0*fSinMiscutAngle)/rsqrt);
   }
   else
   {
       //for straight crystal
       x = x0*fCosMiscutAngle - z0*fSinMiscutAngle;
       y = y0;
       z = x0*fSinMiscutAngle + z0*fCosMiscutAngle;
 
       //for crystalline undulator
       if(fCU){x-=GetCUx(z);}
   }
 
   //calculation of coordinates within a channel (periodic cell)
   fNChannelx=std::floor(x/fDx); //remember the horizontal channel number
                                 //to track the particle
   x-=fNChannelx*fDx;
   fNChannely=std::floor(y/fDy);//remember the vertical channel number
                                //to track the particle (=0 for planar case)
   y-=fNChannely*fDy;
   //correction of the longitudinal coordinate
   if (fBent) {fCorrectionZ = fBendingR/(fBendingR-fNChannelx*fDx);}
 
   return G4ThreeVector(x,y,z);
}

CoordinatesFromLatticeToBox()

G4ThreeVector G4ChannelingFastSimCrystalData::CoordinatesFromLatticeToBox ( const G4ThreeVector & pos )

virtual

calculate the coordinates in the Box reference system (connected with the bounding box of the volume)

Implements G4VChannelingFastSimCrystalData.

Definition at line 360 of file G4ChannelingFastSimCrystalData.cc.

{
   G4double x=pos.x(),y=pos.y(),z=pos.z();
 
   //transform to co-rotating reference system connected with "central plane/axis"
   x+=fNChannelx*fDx;
   y+=fNChannely*fDy;
 
   G4double x0,y0,z0;
 
   if (fBent)
   {
       // for bent crystal
       G4double rcos = (fBendingR - x)*(1. - std::cos(z/fBendingR));
       G4double a = x + rcos;
       G4double b = std::sqrt(x*x + fBending2R*rcos - a*a);
 
       //transform to Box coordinates
       x0 = a*fCosMiscutAngle + b*fSinMiscutAngle;
       y0 = y;
       z0 = b*fCosMiscutAngle - a*fSinMiscutAngle;
   }
   else
   {
       //for crystalline undulator
       if(fCU){x+=GetCUx(z);}
 
       //for straight crystal
       x0 = x*fCosMiscutAngle + z*fSinMiscutAngle;
       y0 = y;
       z0 =-x*fSinMiscutAngle + z*fCosMiscutAngle;
   }
 
   return G4ThreeVector(x0,y0,z0-fHalfDimBoundingBox.z());
}

GetChannelWidthX()

G4double G4ChannelingFastSimCrystalData::GetChannelWidthX ( )

inline

get channel width in x and y

Definition at line 90 of file G4ChannelingFastSimCrystalData.hh.

{return fDx;}

GetChannelWidthY()

G4double G4ChannelingFastSimCrystalData::GetChannelWidthY ( )

inline

Definition at line 91 of file G4ChannelingFastSimCrystalData.hh.

{return fDy;}

SetMaterialProperties()

void G4ChannelingFastSimCrystalData::SetMaterialProperties ( const G4Material * crystal, const G4String & lattice, const G4String & filePath )

virtual

find and upload crystal lattice input files, calculate all the basic values (to do only once)

Implements G4VChannelingFastSimCrystalData.

Definition at line 31 of file G4ChannelingFastSimCrystalData.cc.

{
    G4String filename=crystal->GetName(); //input file
    filename.erase(0,3);
 
    if (fVerbosity)
      {
    G4cout << 
      "======================================================================="  
           << G4endl;
    G4cout <<
      "======                 Crystal lattice data                    ========"
           << G4endl;
    G4cout << 
      "======================================================================="  
           << G4endl;    
    G4cout << "Crystal material: " << filename << G4endl;
      }
 
    //choice between planes (1D model) and axes (2D model)
    if (lattice.compare(0,1,"(")==0)
    {
      iModel=1; //planes
      filename = filename + "_planes_"; //temporary name
      if (fVerbosity)
    G4cout << "Crystal planes: " << lattice << G4endl;
    }
    else if (lattice.compare(0,1,"<")==0)
    {
      iModel=2; //axes
      filename = filename + "_axes_"; //temporary name
      if (fVerbosity)
    G4cout << "Crystal axes: " << lattice << G4endl;
    }
 
    //input file:
    filename = filename + lattice.substr(1,(lattice.length())-2) + ".dat";
 
    if(filePath=="")
    {
      //standard file path if another one is not set
      filename = "/" + filename;  
      filename = G4FindDataDir("G4CHANNELINGDATA") + filename;
    }
    else
    {
      //custom file path
      filename = filePath + filename;
    }
 
    fNelements=(G4int)crystal->GetNumberOfElements();
    for(G4int i=0; i<fNelements; i++)
    {
      fZ1.push_back(crystal->GetElement(i)->GetZ());
      fAN.push_back(crystal->GetElement(i)->GetAtomicMassAmu());
      fI0.push_back(crystal->GetElement(i)->GetIonisation()->GetMeanExcitationEnergy());
    }
 
    G4double var;//just variable
    G4double unitIF;//unit of interpolation function
 
    std::ifstream vfilein;
    vfilein.open(filename);
    //check if the input file was found, otherwise return an exception
    if(!vfilein.is_open())
    {
        G4String outputMessage="Input file " +
                               filename +
                               " is not found!";
        G4Exception("SetMaterialProperties",
                    "001",
                     FatalException,
                    outputMessage);
    }
 
    //read nuclear concentration
    for(G4int i=0; i<fNelements; i++)
    {
      vfilein >> var;
      fN0.push_back(var/CLHEP::cm3);
    }
 
    //read amplitude of thermal oscillations
    for(G4int i=0; i<fNelements; i++)
    {
      vfilein >> var;
      fU1.push_back(var*CLHEP::cm);
    }
 
        if (iModel==1)
    {
      //  read channel dimensions
      vfilein >> fDx;
      fDx*=CLHEP::cm;
      //  read interpolation step size
      vfilein >> fNpointsx;
 
      fDy = fDx;
      fNpointsy = 0;
    }
    else if (iModel==2)
    {
      //  read channel dimensions
      vfilein >> fDx >> fDy;
      fDx*=CLHEP::cm;
      fDy*=CLHEP::cm;
      //  read the number of nodes of interpolation
      vfilein >> fNpointsx >> fNpointsy;
    }
 
    //read the height of the potential well, necessary only for step length calculation
    vfilein >> fVmax;
    fVmax*=CLHEP::eV;
    fVmax2=2.*fVmax;
 
    //read the on-zero minimal potential inside the crystal,
    //necessary for angle recalculation for entrance/exit through
    //the crystal lateral surface
    vfilein >> fVMinCrystal;
    fVMinCrystal*=CLHEP::eV;
 
    // to create the class of interpolation for any function
    fElectricFieldX =
            new G4ChannelingFastSimInterpolation(fDx,fDy,fNpointsx,fNpointsy,iModel);
    if(iModel==2) {fElectricFieldY =
            new G4ChannelingFastSimInterpolation(fDx,fDy,fNpointsx,fNpointsy,iModel);}
    fElectronDensity =
            new G4ChannelingFastSimInterpolation(fDx,fDy,fNpointsx,fNpointsy,iModel);
    fMinIonizationEnergy =
            new G4ChannelingFastSimInterpolation(fDx,fDy,fNpointsx,fNpointsy,iModel);
 
    // do it for any element of crystal material
    for(G4int i=0; i<fNelements; i++)
    {
        fNucleiDensity.push_back(
            new G4ChannelingFastSimInterpolation(fDx,fDy,fNpointsx,fNpointsy,iModel));
    }
 
    if (iModel==1)
    {
      G4double ai, bi, ci, di;
      for(G4int i=0; i<fNpointsx; i++)
      {
        //reading the coefficients of cubic spline
        vfilein >> ai >> bi >> ci >> di;
        //setting spline coefficients for electric field
        unitIF=CLHEP::eV/CLHEP::cm;
        fElectricFieldX->SetCoefficients1D(ai*unitIF, bi*unitIF,
                                           ci*unitIF, di*unitIF, i);
 
        //reading the coefficients of cubic spline
        vfilein >> ai >> bi >> ci >> di;
        //setting spline coefficients for nuclear density (first element)
        unitIF=1.;
        fNucleiDensity[0]->SetCoefficients1D(ai*unitIF, bi*unitIF,
                                             ci*unitIF, di*unitIF, i);
 
        //reading the coefficients of cubic spline
        vfilein >> ai >> bi >> ci >> di;
        //setting spline coefficients for electron density
        unitIF=1./CLHEP::cm3;
        fElectronDensity->SetCoefficients1D(ai*unitIF, bi*unitIF,
                                            ci*unitIF, di*unitIF, i);
 
        //reading the coefficients of cubic spline
        vfilein >> ai >> bi >> ci >> di;
        //setting spline coefficients for minimal ionization energy
        unitIF=CLHEP::eV;
        fMinIonizationEnergy->SetCoefficients1D(ai*unitIF, bi*unitIF,
                                                ci*unitIF, di*unitIF, i);
 
        for(G4int ii=1; ii<fNelements; ii++)
        {
            //reading the coefficients of cubic spline
            vfilein >> ai >> bi >> ci >> di;
            //setting spline coefficients for nuclear density (other elements if any)
            unitIF=1.;
            fNucleiDensity[ii]->SetCoefficients1D(ai*unitIF, bi*unitIF,
                                                  ci*unitIF, di*unitIF, i);
        }
      }
    }
    else if (iModel==2)
    {
      G4double ai3D, bi3D, ci3D;
      for(G4int j=0; j<fNpointsy; j++)
      {
        for(G4int i=0; i<fNpointsx+1; i++)
        {
          for(G4int k=0; k<2; k++)
          {
            //reading the coefficients of cubic spline
            vfilein >> ai3D >> bi3D >> ci3D;
            unitIF=CLHEP::eV;
            //setting spline coefficients for minimal ionization energy
            fMinIonizationEnergy->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF, ci3D*unitIF,
                                                    i, j, k);
 
            //reading the coefficients of cubic spline
            vfilein >> ai3D >> bi3D >> ci3D;
            //setting spline coefficients for horizontal electric field
            unitIF=CLHEP::eV/CLHEP::cm;
            fElectricFieldX->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF, ci3D*unitIF,
                                               i, j, k);
 
            //reading the coefficients of cubic spline
            vfilein >> ai3D >> bi3D >> ci3D;
            //setting spline coefficients for vertical electric field
            unitIF=CLHEP::eV/CLHEP::cm;
            fElectricFieldY->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF, ci3D*unitIF,
                                               i, j, k);
 
            //reading the coefficients of cubic spline
            vfilein >> ai3D >> bi3D >> ci3D;
            //setting spline coefficients for nuclear density (first element)
            unitIF=1.;
            fNucleiDensity[0]->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF, ci3D*unitIF,
                                                 i, j, k);
 
            //reading the coefficients of cubic spline
            vfilein >> ai3D >> bi3D >> ci3D;
            //setting spline coefficients for electron density
            unitIF=1./CLHEP::cm3;
            fElectronDensity->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF, ci3D*unitIF,
                                                i, j, k);
 
            for(G4int ii=1; ii<fNelements; ii++)
            {
                //reading the coefficients of cubic spline
                vfilein >> ai3D >> bi3D >> ci3D;
                //setting spline coefficients for nuclear density (other elements if any)
                unitIF=1.;
                fNucleiDensity[ii]->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF,
                                                      ci3D*unitIF,
                                                      i, j, k);
            }
 
          }
        }
      }
    }
 
    vfilein.close();
 
    //set special values and coefficients
    G4double alphahbarc2=std::pow(CLHEP::fine_structure_const*CLHEP::hbarc ,2.);
    fK30=2.*CLHEP::pi*alphahbarc2/CLHEP::electron_mass_c2;
 
    for(G4int i=0; i<fNelements; i++)
    {
      fRF.push_back((std::pow(9*CLHEP::pi*CLHEP::pi/128/fZ1[i],1/3.))
                    *0.5291772109217*CLHEP::angstrom);//Thomas-Fermi screening radius
 
      fTetamax0.push_back(CLHEP::hbarc/(fR0*std::pow(fAN[i],1./3.)));
      fTeta10.push_back(CLHEP::hbarc/fRF[i]);
      fPu11.push_back(std::pow(fU1[i]/CLHEP::hbarc,2.));
 
      fK20.push_back(alphahbarc2*4*CLHEP::pi*fN0[i]*fZ1[i]*fZ1[i]);
 
      fK40.push_back(3.76*std::pow(CLHEP::fine_structure_const*fZ1[i],2.));
 
      fKD.push_back(fK30*fZ1[i]*fN0[i]);
 
      fLogPlasmaEdI0.push_back(G4Log((crystal->GetIonisation()->GetPlasmaEnergy())/fI0[i]));
    }
 
    fBB.resize(fNelements);
    fE1XBbb.resize(fNelements);
    fBBDEXP.resize(fNelements);
    fPzu11.resize(fNelements);
    fTeta12.resize(fNelements);
    fTetamax2.resize(fNelements);
    fTetamax12.resize(fNelements);
    fK2.resize(fNelements);
 
    fChangeStep = CLHEP::pi*std::min(fDx,fDy)/fNsteps;//necessary to define simulation
                                         //step = fChannelingStep =
                                         // = fChannelingStep0*sqrt(pv)
}

---

The documentation for this class was generated from the following files:

• G4ChannelingFastSimCrystalData.hh
• G4ChannelingFastSimCrystalData.cc