Garfield::TrackHeed Class Reference

Generate tracks using Heed++. More...

#include <TrackHeed.hh>

Inheritance diagram for Garfield::TrackHeed:

Public Member Functions
	TrackHeed ()
	Constructor.
virtual	~TrackHeed ()
	Destructor.
virtual bool	NewTrack (const double x0, const double y0, const double z0, const double t0, const double dx0, const double dy0, const double dz0)
virtual bool	GetCluster (double &xcls, double &ycls, double &zcls, double &tcls, int &n, double &e, double &extra)
bool	GetCluster (double &xcls, double &ycls, double &zcls, double &tcls, int &ne, int &ni, double &e, double &extra)
bool	GetElectron (const unsigned int i, double &x, double &y, double &z, double &t, double &e, double &dx, double &dy, double &dz)
bool	GetIon (const unsigned int i, double &x, double &y, double &z, double &t) const
virtual double	GetClusterDensity ()
virtual double	GetStoppingPower ()
	Get the stopping power (mean energy loss [eV] per cm).
double	GetW () const
double	GetFanoFactor () const
void	TransportDeltaElectron (const double x0, const double y0, const double z0, const double t0, const double e0, const double dx0, const double dy0, const double dz0, int &nel, int &ni)
void	TransportDeltaElectron (const double x0, const double y0, const double z0, const double t0, const double e0, const double dx0, const double dy0, const double dz0, int &nel)
void	TransportPhoton (const double x0, const double y0, const double z0, const double t0, const double e0, const double dx0, const double dy0, const double dz0, int &nel, int &ni)
void	TransportPhoton (const double x0, const double y0, const double z0, const double t0, const double e0, const double dx0, const double dy0, const double dz0, int &nel)
void	EnableElectricField ()
void	DisableElectricField ()
void	EnableMagneticField ()
void	DisableMagneticField ()
void	EnableDeltaElectronTransport ()
void	DisableDeltaElectronTransport ()
void	EnablePhotonReabsorption ()
void	DisablePhotonReabsorption ()
void	EnablePhotoAbsorptionCrossSectionOutput ()
void	DisablePhotoAbsorptionCrossSectionOutput ()
void	SetEnergyMesh (const double e0, const double e1, const int nsteps)
void	SetParticleUser (const double m, const double z)
Public Member Functions inherited from Garfield::Track
	Track ()
	Constructor.
virtual	~Track ()
	Destructor.
virtual void	SetParticle (const std::string &part)
	Set the type of particle.
void	SetEnergy (const double e)
	Set the particle energy.
void	SetBetaGamma (const double bg)
	Set the relative momentum of the particle.
void	SetBeta (const double beta)
	Set the speed ( ) of the particle.
void	SetGamma (const double gamma)
	Set the Lorentz factor of the particle.
void	SetMomentum (const double p)
	Set the particle momentum.
void	SetKineticEnergy (const double ekin)
	Set the kinetic energy of the particle.
double	GetEnergy () const
double	GetBetaGamma () const
double	GetBeta () const
double	GetGamma () const
double	GetMomentum () const
double	GetKineticEnergy () const
double	GetCharge () const
	Get the charge of the projectile.
double	GetMass () const
	Get the mass [eV / c2] of the projectile.
void	SetSensor (Sensor *s)
virtual bool	NewTrack (const double x0, const double y0, const double z0, const double t0, const double dx0, const double dy0, const double dz0)=0
virtual bool	GetCluster (double &xcls, double &ycls, double &zcls, double &tcls, int &n, double &e, double &extra)=0
virtual double	GetClusterDensity ()
virtual double	GetStoppingPower ()
	Get the stopping power (mean energy loss [eV] per cm).
void	EnablePlotting (ViewDrift *viewer)
void	DisablePlotting ()
void	EnableDebugging ()
void	DisableDebugging ()

Additional Inherited Members
Protected Member Functions inherited from Garfield::Track
void	PlotNewTrack (const double x0, const double y0, const double z0)
void	PlotCluster (const double x0, const double y0, const double z0)
Protected Attributes inherited from Garfield::Track
std::string	m_className
double	m_q
int	m_spin
double	m_mass
double	m_energy
double	m_beta2
bool	m_isElectron
std::string	m_particleName
Sensor *	m_sensor
bool	m_isChanged
bool	m_usePlotting
ViewDrift *	m_viewer
bool	m_debug
int	m_plotId

Detailed Description

Generate tracks using Heed++.

Definition at line 37 of file TrackHeed.hh.

Constructor & Destructor Documentation

TrackHeed()

Garfield::TrackHeed::TrackHeed

(

)

Constructor.

Definition at line 60 of file TrackHeed.cc.

    : m_ready(false),
      m_hasActiveTrack(false),
      m_mediumDensity(-1.),
      m_mediumName(""),
      m_usePhotonReabsorption(true),
      m_usePacsOutput(false),
      m_doDeltaTransport(true),
      m_matter(NULL),
      m_gas(NULL),
      m_material(NULL),
      m_atPacs(NULL),
      m_molPacs(NULL),
      m_emin(2.e-6),
      m_emax(2.e-1),
      m_nEnergyIntervals(200),
      m_energyMesh(NULL),
      m_transferCs(NULL),
      m_elScat(NULL),
      m_lowSigma(NULL),
      m_pairProd(NULL),
      m_deltaCs(NULL),
      m_chamber(NULL),
      m_lX(0.),
      m_lY(0.),
      m_lZ(0.),
      m_cX(0.),
      m_cY(0.),
      m_cZ(0.) {
 
  m_className = "TrackHeed";
  m_conductionElectrons.reserve(1000);
  m_conductionIons.reserve(1000);
}

~TrackHeed()

Garfield::TrackHeed::~TrackHeed ( )

virtual

Destructor.

Definition at line 95 of file TrackHeed.cc.

                      {
 
  if (m_matter) delete m_matter;
  if (m_gas) delete m_gas;
  if (m_material) delete m_material;
  if (m_atPacs) delete m_atPacs;
  if (m_molPacs) delete m_molPacs;
  if (m_energyMesh) delete m_energyMesh;
  if (m_transferCs) delete m_transferCs;
  if (m_elScat) delete m_elScat;
  if (m_lowSigma) delete m_lowSigma;
  if (m_pairProd) delete m_pairProd;
  if (m_deltaCs) delete m_deltaCs;
  if (m_chamber) delete m_chamber;
}

Member Function Documentation

DisableDeltaElectronTransport()

void Garfield::TrackHeed::DisableDeltaElectronTransport ( )

inline

Definition at line 100 of file TrackHeed.hh.

{ m_doDeltaTransport = false; }

DisableElectricField()

void Garfield::TrackHeed::DisableElectricField

(

)

Definition at line 848 of file TrackHeed.cc.

{ m_fieldMap.UseEfield(false); }

DisableMagneticField()

void Garfield::TrackHeed::DisableMagneticField

(

)

Definition at line 850 of file TrackHeed.cc.

{ m_fieldMap.UseBfield(false); }

DisablePhotoAbsorptionCrossSectionOutput()

void Garfield::TrackHeed::DisablePhotoAbsorptionCrossSectionOutput ( )

inline

Definition at line 106 of file TrackHeed.hh.

{ m_usePacsOutput = false; }

DisablePhotonReabsorption()

void Garfield::TrackHeed::DisablePhotonReabsorption ( )

inline

Definition at line 103 of file TrackHeed.hh.

{ m_usePhotonReabsorption = false; }

EnableDeltaElectronTransport()

void Garfield::TrackHeed::EnableDeltaElectronTransport ( )

inline

Definition at line 99 of file TrackHeed.hh.

{ m_doDeltaTransport = true; }

Referenced by GarfieldPhysics::InitializePhysics().

EnableElectricField()

void Garfield::TrackHeed::EnableElectricField

(

)

Definition at line 847 of file TrackHeed.cc.

{ m_fieldMap.UseEfield(true); }

EnableMagneticField()

void Garfield::TrackHeed::EnableMagneticField

(

)

Definition at line 849 of file TrackHeed.cc.

{ m_fieldMap.UseBfield(true); }

EnablePhotoAbsorptionCrossSectionOutput()

void Garfield::TrackHeed::EnablePhotoAbsorptionCrossSectionOutput ( )

inline

Definition at line 105 of file TrackHeed.hh.

{ m_usePacsOutput = true; }

EnablePhotonReabsorption()

void Garfield::TrackHeed::EnablePhotonReabsorption ( )

inline

Definition at line 102 of file TrackHeed.hh.

{ m_usePhotonReabsorption = true; }

GetCluster() [1/2]

bool Garfield::TrackHeed::GetCluster ( double &  xcls, double &  ycls, double &  zcls, double &  tcls, int &  n, double &  e, double &  extra  )

virtual

Get the next "cluster" (ionising collision of the charged particle).

Parameters

xcls,ycls,zcls	coordinates of the collision
tcls	time of the collision
n	number of electrons produced
e	deposited energy
extra	additional information (not always implemented)

Implements Garfield::Track.

Definition at line 327 of file TrackHeed.cc.

                                                                           {
 
  int ni = 0;
  return GetCluster(xcls, ycls, zcls, tcls, n, ni, e, extra);
}

Referenced by GarfieldPhysics::DoIt(), and GetCluster().

GetCluster() [2/2]

bool Garfield::TrackHeed::GetCluster	(	double &	xcls,
		double &	ycls,
		double &	zcls,
		double &	tcls,
		int &	ne,
		int &	ni,
		double &	e,
		double &	extra
	)

Definition at line 334 of file TrackHeed.cc.

                                          {
 
  // Initialise and reset.
  xcls = ycls = zcls = tcls = 0.;
  extra = 0.;
  ne = ni = 0;
  e = 0.;
 
  m_deltaElectrons.clear();
  m_conductionElectrons.clear();
  m_conductionIons.clear();
 
  // Make sure NewTrack has been called successfully.
  if (!m_ready) {
    std::cerr << m_className << "::GetCluster:\n"
              << "    Track has not been initialized. Call NewTrack first.\n";
    return false;
  }
 
  std::vector<Heed::gparticle*>::const_iterator end = m_particleBank.end();
  if (m_bankIterator == end) {
    std::cerr << m_className << "::GetCluster:\n"
              << "    There are no more clusters on this track.\n";
    return false;
  }
 
  // Look for the next cluster (i. e. virtual photon) in the list.
  Heed::HeedPhoton* virtualPhoton = NULL;
  for (; m_bankIterator != end; ++m_bankIterator) {
    // Convert the particle to a (virtual) photon.
    virtualPhoton = dynamic_cast<Heed::HeedPhoton*>(*m_bankIterator);
    if (!virtualPhoton) {
      std::cerr << m_className << "::GetCluster:\n"
                << "    Particle is not a virtual photon. Program bug!\n";
      // Try the next element.
      continue;
    }
 
    // Get the location of the interaction (convert from mm to cm
    // and shift with respect to bounding box center).
    xcls = virtualPhoton->currpos.pt.v.x * 0.1 + m_cX;
    ycls = virtualPhoton->currpos.pt.v.y * 0.1 + m_cY;
    zcls = virtualPhoton->currpos.pt.v.z * 0.1 + m_cZ;
    tcls = virtualPhoton->currpos.time;
    // Skip clusters outside the drift area or outside the active medium.
    if (!IsInside(xcls, ycls, zcls)) continue;
    // Add the first ion (at the position of the cluster).
    m_conductionIons.push_back(
        Heed::HeedCondElectron(virtualPhoton->currpos.pt, tcls));
    ++m_bankIterator;
    break;
  }
 
  // Stop if we did not find a virtual photon.
  if (m_bankIterator == end || !virtualPhoton) return false;
  // Plot the cluster, if requested.
  if (m_usePlotting) PlotCluster(xcls, ycls, zcls);
 
  std::vector<Heed::gparticle*> secondaries;
  // Transport the virtual photon.
  virtualPhoton->fly(secondaries);
  // Get the transferred energy (convert from MeV to eV).
  e = virtualPhoton->energy * 1.e6;
 
  while (!secondaries.empty()) {
    std::vector<Heed::gparticle*> newSecondaries;
    // Loop over the secondaries.
    std::vector<Heed::gparticle*>::iterator it;
    std::vector<Heed::gparticle*>::const_iterator send = secondaries.end();
    for (it = secondaries.begin(); it != send; ++it) {
      // Check if it is a delta electron.
      Heed::HeedDeltaElectron* delta = dynamic_cast<Heed::HeedDeltaElectron*>(*it);
      if (delta) {
        extra += delta->curr_kin_energy * 1.e6;
        const double x = delta->currpos.pt.v.x * 0.1 + m_cX;
        const double y = delta->currpos.pt.v.y * 0.1 + m_cY;
        const double z = delta->currpos.pt.v.z * 0.1 + m_cZ;
        if (!IsInside(x, y, z)) continue;
        if (m_doDeltaTransport) {
          // Transport the delta electron.
          delta->fly(newSecondaries);
          // Add the conduction electrons and ions to the list.
          m_conductionElectrons.insert(m_conductionElectrons.end(), 
                                       delta->conduction_electrons.begin(),
                                       delta->conduction_electrons.end());
          m_conductionIons.insert(m_conductionIons.end(), 
                                  delta->conduction_ions.begin(),
                                  delta->conduction_ions.end());
        } else {
          // Add the delta electron to the list, for later use.
          deltaElectron newDeltaElectron;
          newDeltaElectron.x = delta->currpos.pt.v.x * 0.1 + m_cX;
          newDeltaElectron.y = delta->currpos.pt.v.y * 0.1 + m_cY;
          newDeltaElectron.z = delta->currpos.pt.v.z * 0.1 + m_cZ;
          newDeltaElectron.t = delta->currpos.time;
          newDeltaElectron.e = delta->curr_kin_energy * 1.e6;
          newDeltaElectron.dx = delta->currpos.dir.x;
          newDeltaElectron.dy = delta->currpos.dir.y;
          newDeltaElectron.dz = delta->currpos.dir.z;
          m_deltaElectrons.push_back(newDeltaElectron);
        }
        continue;
      } 
      // Check if it is a real photon.
      Heed::HeedPhoton* photon = dynamic_cast<Heed::HeedPhoton*>(*it);
      if (!photon) {
        std::cerr << m_className << "::GetCluster:\n"
                  << "    Particle is neither an electron nor a photon.\n";
      }
      extra += photon->energy * 1.e6;
      const double x = photon->currpos.pt.v.x * 0.1 + m_cX;
      const double y = photon->currpos.pt.v.y * 0.1 + m_cY;
      const double z = photon->currpos.pt.v.z * 0.1 + m_cZ;
      if (!IsInside(x, y, z)) continue;
      // Transport the photon.
      if (m_usePhotonReabsorption) photon->fly(newSecondaries);
    }
    for (it = secondaries.begin(); it != secondaries.end(); ++it) {
      if (*it) delete *it;
    }
    secondaries.clear();
    secondaries.swap(newSecondaries);
  }
  // Get the total number of electrons produced in this step.
  ne = m_doDeltaTransport ? m_conductionElectrons.size() : m_deltaElectrons.size();
  ni = m_conductionIons.size();
  return true;
}

GetClusterDensity()

double Garfield::TrackHeed::GetClusterDensity ( )

virtual

Get the cluster density (number of ionizing collisions per cm or inverse mean free path for ionization).

Reimplemented from Garfield::Track.

Definition at line 293 of file TrackHeed.cc.

                                    {
 
  if (!m_ready) {
    std::cerr << m_className << "::GetClusterDensity:\n"
              << "    Track has not been initialized.\n";
    return 0.;
  }
 
  if (!m_transferCs) {
    std::cerr << m_className << "::GetClusterDensity:\n"
              << "    Ionisation cross-section is not available.\n";
    return 0.;
  }
 
  return m_transferCs->quanC;
}

GetElectron()

bool Garfield::TrackHeed::GetElectron	(	const unsigned int	i,
		double &	x,
		double &	y,
		double &	z,
		double &	t,
		double &	e,
		double &	dx,
		double &	dy,
		double &	dz
	)

Retrieve the properties of a conduction or delta electron in the current cluster.

Parameters

i	index of the electron
x,y,z	coordinates of the electron
t	time
e	kinetic energy (only meaningful for delta-electrons)
dx,dy,dz	direction vector (only meaningful for delta-electrons)

Definition at line 465 of file TrackHeed.cc.

                                                                           {
 
  // Make sure NewTrack has successfully been called.
  if (!m_ready) {
    std::cerr << m_className << "::GetElectron:\n"
              << "    Track has not been initialized. Call NewTrack first.\n";
    return false;
  }
 
  if (m_doDeltaTransport) {
    // Make sure an electron with this number exists.
    if (i >= m_conductionElectrons.size()) {
      std::cerr << m_className << "::GetElectron: Index out of range.\n";
      return false;
    }
 
    x = m_conductionElectrons[i].ptloc.v.x * 0.1 + m_cX;
    y = m_conductionElectrons[i].ptloc.v.y * 0.1 + m_cY;
    z = m_conductionElectrons[i].ptloc.v.z * 0.1 + m_cZ;
    t = m_conductionElectrons[i].time;
    e = 0.;
    dx = dy = dz = 0.;
 
  } else {
    // Make sure a delta electron with this number exists.
    if (i >= m_deltaElectrons.size()) {
      std::cerr << m_className << "::GetElectron:\n"
                << "    Delta electron number out of range.\n";
      return false;
    }
 
    x = m_deltaElectrons[i].x;
    y = m_deltaElectrons[i].y;
    z = m_deltaElectrons[i].z;
    t = m_deltaElectrons[i].t;
    e = m_deltaElectrons[i].e;
    dx = m_deltaElectrons[i].dx;
    dy = m_deltaElectrons[i].dy;
    dz = m_deltaElectrons[i].dz;
  }
  return true;
}

Referenced by GarfieldPhysics::DoIt().

GetFanoFactor()

double Garfield::TrackHeed::GetFanoFactor

(

)

const

Definition at line 1271 of file TrackHeed.cc.

{ return m_matter->F; }

GetIon()

bool Garfield::TrackHeed::GetIon	(	const unsigned int	i,
		double &	x,
		double &	y,
		double &	z,
		double &	t
	)		const

Retrieve the properties of an ion in the current cluster.

Parameters

i	index of the ion
x,y,z	coordinates of the ion
t	time

Definition at line 510 of file TrackHeed.cc.

                                                                         {
 
  // Make sure a "conduction" ion with this number exists.
  if (i >= m_conductionIons.size()) {
    std::cerr << m_className << "::GetIon: Index out of range.\n";
    return false;
  }
 
  x = m_conductionIons[i].ptloc.v.x * 0.1 + m_cX;
  y = m_conductionIons[i].ptloc.v.y * 0.1 + m_cY;
  z = m_conductionIons[i].ptloc.v.z * 0.1 + m_cZ;
  t = m_conductionIons[i].time;
  return true;
}

GetStoppingPower()

double Garfield::TrackHeed::GetStoppingPower ( )

virtual

Get the stopping power (mean energy loss [eV] per cm).

Reimplemented from Garfield::Track.

Definition at line 310 of file TrackHeed.cc.

                                   {
 
  if (!m_ready) {
    std::cerr << m_className << "::GetStoppingPower:\n"
              << "    Track has not been initialized.\n";
    return 0.;
  }
 
  if (!m_transferCs) {
    std::cerr << m_className << "::GetStoppingPower:\n"
              << "    Ionisation cross-section is not available.\n";
    return 0.;
  }
 
  return m_transferCs->meanC1 * 1.e6;
}

GetW()

double Garfield::TrackHeed::GetW

(

)

const

Definition at line 1270 of file TrackHeed.cc.

{ return m_matter->W * 1.e6; }

Referenced by GarfieldPhysics::DoIt().

NewTrack()

bool Garfield::TrackHeed::NewTrack ( const double  x0, const double  y0, const double  z0, const double  t0, const double  dx0, const double  dy0, const double  dz0  )

virtual

Calculate a new track starting from (x0, y0, z0) at time t0 in direction (dx0, dy0, dz0).

Implements Garfield::Track.

Definition at line 111 of file TrackHeed.cc.

                                           {
 
  m_hasActiveTrack = false;
  m_ready = false;
 
  // Make sure the sensor has been set.
  if (!m_sensor) {
    std::cerr << m_className << "::NewTrack:\n"
              << "    Sensor is not defined.\n";
    return false;
  }
 
  // Get the bounding box.
  double xmin = 0., ymin = 0., zmin = 0.;
  double xmax = 0., ymax = 0., zmax = 0.;
  if (!m_sensor->GetArea(xmin, ymin, zmin, xmax, ymax, zmax)) {
    std::cerr << m_className << "::NewTrack:\n"
              << "    Drift area is not set.\n";
    return false;
  }
  // Check if the bounding box has changed.
  const double lx = fabs(xmax - xmin);
  const double ly = fabs(ymax - ymin);
  const double lz = fabs(zmax - zmin);
  if (m_debug) {
    std::cout << m_className << "::NewTrack:\n"
              << "    Bounding box dimensions:\n"
              << "      x: " << lx << " cm\n"
              << "      y: " << ly << " cm\n"
              << "      z: " << lz << " cm\n";
  }
  if (fabs(lx - m_lX) > Small || fabs(ly - m_lY) > Small || fabs(lz - m_lZ) > Small) {
    m_lX = lx;
    m_lY = ly;
    m_lZ = lz;
    m_isChanged = true;
  }
  // Update the center of the bounding box.
  if (std::isinf(xmin) || std::isinf(xmax)) {
    m_cX = 0.;
  } else {
    m_cX = 0.5 * (xmin + xmax);
  }
  if (std::isinf(ymin) || std::isinf(ymax)) {
    m_cY = 0.;
  } else {
    m_cY = 0.5 * (ymin + ymax);
  }
  if (std::isinf(zmin) || std::isinf(zmax)) {
    m_cZ = 0.;
  } else {
    m_cZ = 0.5 * (zmin + zmax);
  }
  if (m_debug) {
    std::cout << m_className << "::NewTrack:\n"
              << "    Center of bounding box:\n"
              << "      x: " << m_cX << " cm\n"
              << "      y: " << m_cY << " cm\n"
              << "      z: " << m_cZ << " cm\n";
  }
 
  m_fieldMap.SetSensor(m_sensor);
  m_fieldMap.SetCentre(m_cX, m_cY, m_cZ);
 
  // Make sure the initial position is inside an ionisable medium.
  Medium* medium = NULL;
  if (!m_sensor->GetMedium(x0, y0, z0, medium)) {
    std::cerr << m_className << "::NewTrack:\n"
              << "    No medium at initial position.\n";
    return false;
  } else if (!medium->IsIonisable()) {
    std::cerr << m_className << "::NewTrack:\n"
              << "    Medium at initial position is not ionisable.\n";
    return false;
  }
 
  // Check if the medium has changed since the last call.
  if (medium->GetName() != m_mediumName ||
      fabs(medium->GetMassDensity() - m_mediumDensity) > 1.e-9) {
    m_isChanged = true;
  }
 
  // If medium, particle or bounding box have changed,
  // update the cross-sections.
  if (m_isChanged) {
    if (!Setup(medium)) return false;
    m_isChanged = false;
    m_mediumName = medium->GetName();
    m_mediumDensity = medium->GetMassDensity();
  }
 
  ClearParticleBank();
 
  m_deltaElectrons.clear();
  m_conductionElectrons.clear();
  m_conductionIons.clear();
 
  // Check the direction vector.
  double dx = dx0, dy = dy0, dz = dz0;
  const double d = sqrt(dx * dx + dy * dy + dz * dz);
  if (d < Small) {
    if (m_debug) {
      std::cout << m_className << "::NewTrack:\n"
                << "    Direction vector has zero norm.\n"
                << "    Initial direction is randomized.\n";
    }
    // Null vector. Sample the direction isotropically.
    RndmDirection(dx, dy, dz);
  } else {
    // Normalise the direction vector.
    dx /= d;
    dy /= d;
    dz /= d;
  }
  Heed::vec velocity(dx, dy, dz);
  velocity = velocity * Heed::CLHEP::c_light * GetBeta();
 
  if (m_debug) {
    std::cout << m_className << "::NewTrack:\n    Track starts at ("
              << x0 << ", " << y0 << ", " << z0 << ") at time " << t0 << "\n"
              << "    Direction: (" << dx << ", " << dy << ", " << dz << ")\n";
  }
 
  // Initial position (shift with respect to bounding box center and
  // convert from cm to mm).
  Heed::point p0((x0 - m_cX) * 10., (y0 - m_cY) * 10., (z0 - m_cZ) * 10.);
 
  // Setup the particle.
  Heed::particle_def* particleType = &Heed::muon_minus_def;
  if (m_particleName == "e-") {
    particleType = &Heed::electron_def;
  } else if (m_particleName == "e+") {
    particleType = &Heed::positron_def;
  } else if (m_particleName == "mu-") {
    particleType = &Heed::muon_minus_def;
  } else if (m_particleName == "mu+") {
    particleType = &Heed::muon_plus_def;
  } else if (m_particleName == "pi-") {
    particleType = &Heed::pi_minus_meson_def;
  } else if (m_particleName == "pi+") {
    particleType = &Heed::pi_plus_meson_def;
  } else if (m_particleName == "K-") {
    particleType = &Heed::K_minus_meson_def;
  } else if (m_particleName == "K+") {
    particleType = &Heed::K_plus_meson_def;
  } else if (m_particleName == "p") {
    particleType = &Heed::proton_def;
  } else if (m_particleName == "pbar") {
    particleType = &Heed::anti_proton_def;
  } else if (m_particleName == "d") {
    particleType = &Heed::deuteron_def;
  } else if (m_particleName == "alpha") {
    particleType = &Heed::alpha_particle_def;
  } else if (m_particleName == "exotic") {
    // User defined particle
    Heed::user_particle_def.set_mass(m_mass * 1.e-6);
    Heed::user_particle_def.set_charge(m_q);
    particleType = &Heed::user_particle_def;
  } else {
    // Not a predefined particle, use muon definition.
    if (m_q > 0.) {
      particleType = &Heed::muon_minus_def;
    } else {
      particleType = &Heed::muon_plus_def;
    }
  }
 
  Heed::HeedParticle particle(m_chamber, p0, velocity, t0,
                              particleType, &m_fieldMap);
  // Transport the particle.
  particle.fly(m_particleBank);
  m_bankIterator = m_particleBank.begin();
  m_hasActiveTrack = true;
  m_ready = true;
 
  // Plot the new track.
  if (m_usePlotting) PlotNewTrack(x0, y0, z0);
  return true;
}

Referenced by GarfieldPhysics::DoIt().

SetEnergyMesh()

void Garfield::TrackHeed::SetEnergyMesh	(	const double	e0,
		const double	e1,
		const int	nsteps
	)

Definition at line 852 of file TrackHeed.cc.

                                                {
 
  if (fabs(e1 - e0) < Small) {
    std::cerr << m_className << "::SetEnergyMesh:\n";
    std::cerr << "    Invalid energy range:\n";
    std::cerr << "    " << e0 << " < E [eV] < " << e1 << "\n";
    return;
  }
 
  if (nsteps <= 0) {
    std::cerr << m_className << "::SetEnergyMesh:\n";
    std::cerr << "    Number of intervals must be > 0.\n";
    return;
  }
 
  m_emin = std::min(e0, e1);
  m_emax = std::max(e0, e1);
  m_emin *= 1.e-6;
  m_emax *= 1.e-6;
  m_nEnergyIntervals = nsteps;
}

SetParticleUser()

void Garfield::TrackHeed::SetParticleUser	(	const double	m,
		const double	z
	)

Definition at line 875 of file TrackHeed.cc.

                                                              {
 
  if (fabs(z) < Small) {
    std::cerr << m_className << "::SetParticleUser:\n"
              << "    Particle cannot have zero charge.\n";
    return;
  }
  if (m < Small) {
    std::cerr << m_className << "::SetParticleUser:\n"
              << "    Particle mass must be greater than zero.\n";
  }
  m_q = z;
  m_mass = m;
  m_isElectron = false;
  m_spin = 0;
  m_particleName = "exotic";
}

TransportDeltaElectron() [1/2]

void Garfield::TrackHeed::TransportDeltaElectron	(	const double	x0,
		const double	y0,
		const double	z0,
		const double	t0,
		const double	e0,
		const double	dx0,
		const double	dy0,
		const double	dz0,
		int &	nel
	)

Definition at line 526 of file TrackHeed.cc.

                                                 {
  int ni = 0;
  return TransportDeltaElectron(x0, y0, z0, t0, e0, dx0, dy0, dz0, nel, ni);
}

TransportDeltaElectron() [2/2]

void Garfield::TrackHeed::TransportDeltaElectron	(	const double	x0,
		const double	y0,
		const double	z0,
		const double	t0,
		const double	e0,
		const double	dx0,
		const double	dy0,
		const double	dz0,
		int &	nel,
		int &	ni
	)

Definition at line 535 of file TrackHeed.cc.

                                                          {
 
  nel = 0;
  ni = 0;
 
  // Check if delta electron transport was disabled.
  if (!m_doDeltaTransport) {
    std::cerr << m_className << "::TransportDeltaElectron:\n"
              << "    Delta electron transport has been switched off.\n";
    return;
  }
 
 
  // Make sure the sensor has been set.
  if (!m_sensor) {
    std::cerr << m_className << "::TransportDeltaElectron:\n"
              << "    Sensor is not defined.\n";
    m_ready = false;
    return;
  }
 
  // Get the bounding box.
  double xmin, ymin, zmin;
  double xmax, ymax, zmax;
  if (!m_sensor->GetArea(xmin, ymin, zmin, xmax, ymax, zmax)) {
    std::cerr << m_className << "::TransportDeltaElectron:\n"
              << "    Drift area is not set.\n";
    m_ready = false;
    return;
  }
  // Check if the bounding box has changed.
  bool update = false;
  const double lx = fabs(xmax - xmin);
  const double ly = fabs(ymax - ymin);
  const double lz = fabs(zmax - zmin);
  if (fabs(lx - m_lX) > Small || fabs(ly - m_lY) > Small || fabs(lz - m_lZ) > Small) {
    m_lX = lx;
    m_lY = ly;
    m_lZ = lz;
    m_isChanged = true;
    update = true;
    m_hasActiveTrack = false;
  }
  // Update the center of the bounding box.
  m_cX = 0.5 * (xmin + xmax);
  m_cY = 0.5 * (ymin + ymax);
  m_cZ = 0.5 * (zmin + zmax);
 
  m_fieldMap.SetSensor(m_sensor);
  m_fieldMap.SetCentre(m_cX, m_cY, m_cZ);
 
  // Make sure the initial position is inside an ionisable medium.
  Medium* medium = NULL;
  if (!m_sensor->GetMedium(x0, y0, z0, medium)) {
    std::cerr << m_className << "::TransportDeltaElectron:\n"
              << "    No medium at initial position.\n";
    return;
  } else if (!medium->IsIonisable()) {
    std::cerr << "TrackHeed:TransportDeltaElectron:\n"
              << "    Medium at initial position is not ionisable.\n";
    m_ready = false;
    return;
  }
 
  // Check if the medium has changed since the last call.
  if (medium->GetName() != m_mediumName ||
      fabs(medium->GetMassDensity() - m_mediumDensity) > 1.e-9) {
    m_isChanged = true;
    update = true;
    m_ready = false;
    m_hasActiveTrack = false;
  }
 
  // If medium or bounding box have changed, update the "chamber".
  if (update) {
    if (!Setup(medium)) return;
    m_ready = true;
    m_mediumName = medium->GetName();
    m_mediumDensity = medium->GetMassDensity();
  }
 
  m_deltaElectrons.clear();
  m_conductionElectrons.clear();
  m_conductionIons.clear();
 
  // Initial position (shift with respect to bounding box center and
  // convert from cm to mm).
  Heed::point p0((x0 - m_cX) * 10., (y0 - m_cY) * 10., (z0 - m_cZ) * 10.);
 
  // Make sure the kinetic energy is positive.
  if (e0 <= 0.) {
    // Just create a conduction electron on the spot.
    m_conductionElectrons.push_back(Heed::HeedCondElectron(p0, t0));
    nel = 1;
    return;
  }
 
  // Check the direction vector.
  double dx = dx0, dy = dy0, dz = dz0;
  const double d = sqrt(dx * dx + dy * dy + dz * dz);
  if (d <= 0.) {
    // Null vector. Sample the direction isotropically.
    RndmDirection(dx, dy, dz);
  } else {
    // Normalise the direction vector.
    dx /= d;
    dy /= d;
    dz /= d;
  }
  Heed::vec velocity(dx, dy, dz);
 
  // Calculate the speed for the given kinetic energy.
  const double gamma = 1. + e0 / ElectronMass;
  const double beta = sqrt(1. - 1. / (gamma * gamma));
  double speed = Heed::CLHEP::c_light * beta;
  velocity = velocity * speed;
 
  // Transport the electron.
  std::vector<Heed::gparticle*> secondaries;
  Heed::HeedDeltaElectron delta(m_chamber, p0, velocity, t0, 0, &m_fieldMap);
  delta.fly(secondaries);
  ClearBank(secondaries);
 
  m_conductionElectrons.swap(delta.conduction_electrons);
  m_conductionIons.swap(delta.conduction_ions);
  nel = m_conductionElectrons.size();
  ni = m_conductionIons.size();
}

Referenced by GarfieldPhysics::DoIt(), and TransportDeltaElectron().

TransportPhoton() [1/2]

void Garfield::TrackHeed::TransportPhoton	(	const double	x0,
		const double	y0,
		const double	z0,
		const double	t0,
		const double	e0,
		const double	dx0,
		const double	dy0,
		const double	dz0,
		int &	nel
	)

Definition at line 668 of file TrackHeed.cc.

                                                                              {
  int ni = 0;
  TransportPhoton(x0, y0, z0, t0, e0, dx0, dy0, dz0, nel, ni);
}

TransportPhoton() [2/2]

void Garfield::TrackHeed::TransportPhoton	(	const double	x0,
		const double	y0,
		const double	z0,
		const double	t0,
		const double	e0,
		const double	dx0,
		const double	dy0,
		const double	dz0,
		int &	nel,
		int &	ni
	)

Definition at line 676 of file TrackHeed.cc.

                                                   {
 
  nel = 0;
  ni = 0;
 
  // Make sure the energy is positive.
  if (e0 <= 0.) {
    std::cerr << m_className << "::TransportPhoton:\n"
              << "    Photon energy must be positive.\n";
    return;
  }
 
  // Make sure the sensor has been set.
  if (!m_sensor) {
    std::cerr << m_className << "::TransportPhoton:\n"
              << "    Sensor is not defined.\n";
    m_ready = false;
    return;
  }
 
  // Get the bounding box.
  double xmin, ymin, zmin;
  double xmax, ymax, zmax;
  if (!m_sensor->GetArea(xmin, ymin, zmin, xmax, ymax, zmax)) {
    std::cerr << m_className << "::TransportPhoton:\n"
              << "    Drift area is not set.\n";
    m_ready = false;
    return;
  }
  // Check if the bounding box has changed.
  bool update = false;
  const double lx = fabs(xmax - xmin);
  const double ly = fabs(ymax - ymin);
  const double lz = fabs(zmax - zmin);
  if (fabs(lx - m_lX) > Small || fabs(ly - m_lY) > Small || fabs(lz - m_lZ) > Small) {
    m_lX = lx;
    m_lY = ly;
    m_lZ = lz;
    m_isChanged = true;
    update = true;
    m_hasActiveTrack = false;
  }
  // Update the center of the bounding box.
  m_cX = 0.5 * (xmin + xmax);
  m_cY = 0.5 * (ymin + ymax);
  m_cZ = 0.5 * (zmin + zmax);
 
  m_fieldMap.SetSensor(m_sensor);
  m_fieldMap.SetCentre(m_cX, m_cY, m_cZ);
 
  // Make sure the initial position is inside an ionisable medium.
  Medium* medium = NULL;
  if (!m_sensor->GetMedium(x0, y0, z0, medium)) {
    std::cerr << m_className << "::TransportPhoton:\n"
              << "    No medium at initial position.\n";
    return;
  } else if (!medium->IsIonisable()) {
    std::cerr << "TrackHeed:TransportPhoton:\n"
              << "    Medium at initial position is not ionisable.\n";
    m_ready = false;
    return;
  }
 
  // Check if the medium has changed since the last call.
  if (medium->GetName() != m_mediumName ||
      fabs(medium->GetMassDensity() - m_mediumDensity) > 1.e-9) {
    m_isChanged = true;
    update = true;
    m_ready = false;
  }
 
  // If medium or bounding box have changed, update the "chamber".
  if (update) {
    if (!Setup(medium)) return;
    m_ready = true;
    m_mediumName = medium->GetName();
    m_mediumDensity = medium->GetMassDensity();
  }
 
  // Delete the particle bank.
  // Clusters from the current track will be lost.
  m_hasActiveTrack = false;
  ClearParticleBank();
  m_deltaElectrons.clear();
  m_conductionElectrons.clear();
  m_conductionIons.clear();
 
  // Check the direction vector.
  double dx = dx0, dy = dy0, dz = dz0;
  const double d = sqrt(dx * dx + dy * dy + dz * dz);
  if (d <= 0.) {
    // Null vector. Sample the direction isotropically.
    RndmDirection(dx, dy, dz);
  } else {
    // Normalise the direction vector.
    dx /= d;
    dy /= d;
    dz /= d;
  }
  Heed::vec velocity(dx, dy, dz);
  velocity = velocity * Heed::CLHEP::c_light;
 
  // Initial position (shift with respect to bounding box center and
  // convert from cm to mm).
  Heed::point p0((x0 - m_cX) * 10., (y0 - m_cY) * 10., (z0 - m_cZ) * 10.);
 
  // Create and transport the photon.
  Heed::HeedPhoton photon(m_chamber, p0, velocity, t0, 0, e0 * 1.e-6, 
                          &m_fieldMap);
  ClearParticleBank();
  std::vector<Heed::gparticle*> secondaries;
  photon.fly(secondaries);
 
  while (!secondaries.empty()) {
    std::vector<Heed::gparticle*> newSecondaries;
    // Loop over the particle bank and look for daughter particles.
    std::vector<Heed::gparticle*>::iterator it;
    for (it = secondaries.begin(); it != secondaries.end(); ++it) {
      // Check if it is a delta electron.
      Heed::HeedDeltaElectron* delta = dynamic_cast<Heed::HeedDeltaElectron*>(*it);
      if (delta) {
        if (m_doDeltaTransport) {
          // Transport the delta electron.
          delta->fly(newSecondaries);
          // Add the conduction electrons to the list.
          m_conductionElectrons.insert(m_conductionElectrons.end(), 
                                       delta->conduction_electrons.begin(),
                                       delta->conduction_electrons.end());
          m_conductionIons.insert(m_conductionIons.end(), 
                                  delta->conduction_ions.begin(),
                                  delta->conduction_ions.end());
        } else {
          // Add the delta electron to the list, for later use.
          deltaElectron newDeltaElectron;
          newDeltaElectron.x = delta->currpos.pt.v.x * 0.1 + m_cX;
          newDeltaElectron.y = delta->currpos.pt.v.y * 0.1 + m_cY;
          newDeltaElectron.z = delta->currpos.pt.v.z * 0.1 + m_cZ;
          newDeltaElectron.t = delta->currpos.time;
          newDeltaElectron.e = delta->curr_kin_energy * 1.e6;
          newDeltaElectron.dx = delta->currpos.dir.x;
          newDeltaElectron.dy = delta->currpos.dir.y;
          newDeltaElectron.dz = delta->currpos.dir.z;
          m_deltaElectrons.push_back(newDeltaElectron);
        }
        continue;
      }
      // Check if it is a fluorescence photon.
      Heed::HeedPhoton* fluorescencePhoton = dynamic_cast<Heed::HeedPhoton*>(*it);
      if (!fluorescencePhoton) {
        std::cerr << m_className << "::TransportPhoton:\n"
                  << "    Unknown secondary particle.\n";
        ClearBank(secondaries);
        ClearBank(newSecondaries);
        return;
      }
      fluorescencePhoton->fly(newSecondaries);
    }
    secondaries.swap(newSecondaries);
    ClearBank(newSecondaries);
  }
 
  // Get the total number of electrons produced in this step.
  nel = m_doDeltaTransport ? m_conductionElectrons.size() : 
                             m_deltaElectrons.size();
  ni = m_conductionIons.size();
}

Referenced by GarfieldPhysics::DoIt(), and TransportPhoton().

---

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

• TrackHeed.hh
• TrackHeed.cc