Garfield::TrackHeed Class Reference

#include <TrackHeed.hh>

Inheritance diagram for Garfield::TrackHeed:

Public Member Functions
	TrackHeed ()
	~TrackHeed ()
bool	NewTrack (const double x0, const double y0, const double z0, const double t0, const double dx0, const double dy0, const double dz0)
bool	GetCluster (double &xcls, double &ycls, double &zcls, double &tcls, int &n, double &e, double &extra)
bool	GetElectron (const int i, double &x, double &y, double &z, double &t, double &e, double &dx, double &dy, double &dz)
double	GetClusterDensity ()
double	GetStoppingPower ()
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)
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 ()
virtual	~Track ()
virtual void	SetParticle (std::string part)
void	SetEnergy (const double e)
void	SetBetaGamma (const double bg)
void	SetBeta (const double beta)
void	SetGamma (const double gamma)
void	SetMomentum (const double p)
void	SetKineticEnergy (const double ekin)
double	GetEnergy () const
double	GetBetaGamma () const
double	GetBeta () const
double	GetGamma () const
double	GetMomentum () const
double	GetKineticEnergy () const
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 ()
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	className
double	q
int	spin
double	mass
double	energy
double	beta2
bool	isElectron
std::string	particleName
Sensor *	sensor
bool	isChanged
bool	usePlotting
ViewDrift *	viewer
bool	debug
int	plotId

Detailed Description

Definition at line 30 of file TrackHeed.hh.

Constructor & Destructor Documentation

TrackHeed()

Garfield::TrackHeed::TrackHeed

(

)

Definition at line 98 of file TrackHeed.cc.

    : ready(false),
      hasActiveTrack(false),
      mediumDensity(-1.),
      mediumName(""),
      usePhotonReabsorption(true),
      usePacsOutput(false),
      useDelta(true),
      nDeltas(0),
      particle(0),
      matter(0),
      gas(0),
      material(0),
      m_atPacs(0),
      m_molPacs(0),
      emin(2.e-6),
      emax(2.e-1),
      nEnergyIntervals(200),
      energyMesh(0),
      transferCs(0),
      elScat(0),
      lowSigma(0),
      pairProd(0),
      deltaCs(0),
      chamber(0),
      lX(0.),
      lY(0.),
      lZ(0.),
      cX(0.),
      cY(0.),
      cZ(0.) {
 
  className = "TrackHeed";
 
  HeedInterface::sensor = 0;
  HeedInterface::useEfield = false;
  HeedInterface::useBfield = false;
 
  deltaElectrons.clear();
}

~TrackHeed()

Garfield::TrackHeed::~TrackHeed

(

)

Definition at line 139 of file TrackHeed.cc.

                      {
 
  if (particle != 0) delete particle;
  if (matter != 0) delete matter;
  if (gas != 0) delete gas;
  if (material != 0) delete material;
  if (m_atPacs != 0) delete m_atPacs;
  if (m_molPacs != 0) delete m_molPacs;
  if (energyMesh != 0) delete energyMesh;
  if (transferCs != 0) delete transferCs;
  if (elScat != 0) delete elScat;
  if (lowSigma != 0) delete lowSigma;
  if (pairProd != 0) delete pairProd;
  if (deltaCs != 0) delete deltaCs;
  if (chamber != 0) delete chamber;
 
  Garfield::HeedInterface::sensor = 0;
}

Member Function Documentation

DisableDeltaElectronTransport()

void Garfield::TrackHeed::DisableDeltaElectronTransport ( )

inline

Definition at line 68 of file TrackHeed.hh.

{ useDelta = false; }

DisableElectricField()

void Garfield::TrackHeed::DisableElectricField

(

)

Definition at line 944 of file TrackHeed.cc.

{ HeedInterface::useEfield = false; }

DisableMagneticField()

void Garfield::TrackHeed::DisableMagneticField

(

)

Definition at line 948 of file TrackHeed.cc.

{ HeedInterface::useBfield = false; }

DisablePhotoAbsorptionCrossSectionOutput()

void Garfield::TrackHeed::DisablePhotoAbsorptionCrossSectionOutput ( )

inline

Definition at line 74 of file TrackHeed.hh.

{ usePacsOutput = false; }

DisablePhotonReabsorption()

void Garfield::TrackHeed::DisablePhotonReabsorption ( )

inline

Definition at line 71 of file TrackHeed.hh.

{ usePhotonReabsorption = false; }

EnableDeltaElectronTransport()

void Garfield::TrackHeed::EnableDeltaElectronTransport ( )

inline

Definition at line 67 of file TrackHeed.hh.

{ useDelta = true; }

Referenced by GarfieldPhysics::InitializePhysics().

EnableElectricField()

void Garfield::TrackHeed::EnableElectricField

(

)

Definition at line 942 of file TrackHeed.cc.

{ HeedInterface::useEfield = true; }

EnableMagneticField()

void Garfield::TrackHeed::EnableMagneticField

(

)

Definition at line 946 of file TrackHeed.cc.

{ HeedInterface::useBfield = true; }

EnablePhotoAbsorptionCrossSectionOutput()

void Garfield::TrackHeed::EnablePhotoAbsorptionCrossSectionOutput ( )

inline

Definition at line 73 of file TrackHeed.hh.

{ usePacsOutput = true; }

EnablePhotonReabsorption()

void Garfield::TrackHeed::EnablePhotonReabsorption ( )

inline

Definition at line 70 of file TrackHeed.hh.

{ usePhotonReabsorption = true; }

GetCluster()

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

virtual

Implements Garfield::Track.

Definition at line 383 of file TrackHeed.cc.

                                                                           {
 
  // Initial settings.
  xcls = ycls = zcls = tcls = 0.;
  extra = 0.;
  n = 0;
  e = 0.;
 
  // Make sure NewTrack has successfully been called.
  if (!ready) {
    std::cerr << className << "::GetCluster:\n";
    std::cerr << "    Track has not been initialized.\n";
    std::cerr << "    Call NewTrack first.\n";
    return false;
  }
 
  if (!hasActiveTrack) {
    std::cerr << className << "::GetCluster:\n";
    std::cerr << "    There are no more clusters.\n";
    return false;
  }
 
  bool ok = false;
  Medium* medium = 0;
  AbsListNode<ActivePtr<gparticle> >* node = 0;
  Heed::HeedPhoton* virtualPhoton = 0;
  while (!ok) {
    // Get the first element from the particle bank.
    node = Heed::particle_bank.get_first_node();
 
    // Make sure the particle bank is not empty.
    if (node == 0) {
      hasActiveTrack = false;
      return false;
    }
 
    // Convert the particle to a (virtual) photon.
    virtualPhoton = dynamic_cast<Heed::HeedPhoton*>(node->el.get());
    if (virtualPhoton == 0) {
      std::cerr << className << "::GetCluster:\n";
      std::cerr << "    Particle is not a virtual photon.\n";
      std::cerr << "    Program bug!\n";
      // Delete the node.
      Heed::particle_bank.erase(node);
      // Try the next node.
      continue;
    }
 
    if (virtualPhoton->parent_particle_number != 0) {
      std::cerr << className << "::GetCluster:\n";
      std::cerr << "    Virtual photon has an unexpected parent.\n";
      // Delete this virtual photon.
      Heed::particle_bank.erase(node);
      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 + cX;
    ycls = virtualPhoton->currpos.pt.v.y * 0.1 + cY;
    zcls = virtualPhoton->currpos.pt.v.z * 0.1 + cZ;
    tcls = virtualPhoton->currpos.time;
    // Make sure the cluster is inside the drift area.
    if (!sensor->IsInArea(xcls, ycls, zcls)) {
      // Delete this virtual photon and proceed with the next one.
      Heed::particle_bank.erase(node);
      continue;
    }
    // Make sure the cluster is inside a medium.
    if (!sensor->GetMedium(xcls, ycls, zcls, medium)) {
      // Delete this virtual photon and proceed with the next one.
      Heed::particle_bank.erase(node);
      continue;
    }
    // Make sure the medium has not changed.
    if (medium->GetName() != mediumName ||
        fabs(medium->GetMassDensity() - mediumDensity) > 1.e-9 ||
        !medium->IsIonisable()) {
      // Delete this virtual photon and proceed with the next one.
      Heed::particle_bank.erase(node);
      continue;
    }
    // Seems to be ok.
    ok = true;
  }
 
  // Plot the cluster, if requested.
  if (usePlotting) PlotCluster(xcls, ycls, zcls);
 
  // Transport the virtual photon.
  virtualPhoton->fly();
  // Get the transferred energy (convert from MeV to eV).
  e = virtualPhoton->energy * 1.e6;
 
  // Make a list of parent particle id numbers.
  std::vector<int> ids;
  ids.clear();
  // At the beginning, there is only the virtual photon.
  ids.push_back(virtualPhoton->particle_number);
  int nIds = 1;
 
  // Look for daughter particles.
  deltaElectrons.clear();
  nDeltas = 0;
  chamber->conduction_electron_bank.allocate_block(1000);
  bool deleteNode = false;
  Heed::HeedDeltaElectron* delta = 0;
  Heed::HeedPhoton* photon = 0;
  AbsListNode<ActivePtr<gparticle> >* nextNode = node->get_next_node();
  AbsListNode<ActivePtr<gparticle> >* tempNode = 0;
  // Loop over the particle bank.
  while (nextNode != 0) {
    deleteNode = false;
    // Check if it is a delta electron.
    delta = dynamic_cast<Heed::HeedDeltaElectron*>(nextNode->el.get());
    if (delta != 0) {
      // Check if the delta electron was produced by one of the photons
      // belonging to this cluster.
      for (int i = nIds; i--;) {
        if (delta->parent_particle_number == ids[i]) {
          if (useDelta) {
            // Transport the delta electron.
            delta->fly();
          } else {
            // Add the delta electron to the list, for later use.
            deltaElectron newDeltaElectron;
            newDeltaElectron.x = delta->currpos.pt.v.x * 0.1 + cX;
            newDeltaElectron.y = delta->currpos.pt.v.y * 0.1 + cY;
            newDeltaElectron.z = delta->currpos.pt.v.z * 0.1 + 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;
            deltaElectrons.push_back(newDeltaElectron);
            ++nDeltas;
          }
          deleteNode = true;
          break;
        }
      }
    } else {
      // Check if it is a real photon.
      photon = dynamic_cast<Heed::HeedPhoton*>(nextNode->el.get());
      if (photon == 0) {
        std::cerr << className << "::GetCluster:\n";
        std::cerr << "    Particle is neither an electron nor a photon.\n";
        return false;
      }
      for (int i = nIds; i--;) {
        if (photon->parent_particle_number == ids[i]) {
          // Transport the photon and add its number to the list of ids.
          if (usePhotonReabsorption) photon->fly();
          deleteNode = true;
          ids.push_back(photon->particle_number);
          ++nIds;
          break;
        }
      }
    }
    // Proceed with the next node in the particle bank.
    if (deleteNode) {
      tempNode = nextNode->get_next_node();
      Heed::particle_bank.erase(nextNode);
      nextNode = tempNode;
    } else {
      nextNode = nextNode->get_next_node();
    }
  }
 
  // Get the total number of electrons produced in this step.
  if (useDelta) {
    n = chamber->conduction_electron_bank.get_qel();
  } else {
    n = nDeltas;
  }
 
  // Remove the virtual photon from the particle bank.
  Heed::particle_bank.erase(node);
 
  return true;
}

Referenced by GarfieldPhysics::DoIt().

GetClusterDensity()

double Garfield::TrackHeed::GetClusterDensity ( )

virtual

Reimplemented from Garfield::Track.

Definition at line 349 of file TrackHeed.cc.

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

GetElectron()

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

Definition at line 566 of file TrackHeed.cc.

                                        {
 
  // Make sure NewTrack has successfully been called.
  if (!ready) {
    std::cerr << className << "::GetElectron:\n";
    std::cerr << "    Track has not been initialized.\n";
    std::cerr << "    Call NewTrack first.\n";
    return false;
  }
 
  if (useDelta) {
    // Make sure an electron with this number exists.
    const int n = chamber->conduction_electron_bank.get_qel();
    if (i < 0 || i >= n) {
      std::cerr << className << "::GetElectron:\n";
      std::cerr << "    Electron number out of range.\n";
      return false;
    }
 
    x = chamber->conduction_electron_bank[i].ptloc.v.x * 0.1 + cX;
    y = chamber->conduction_electron_bank[i].ptloc.v.y * 0.1 + cY;
    z = chamber->conduction_electron_bank[i].ptloc.v.z * 0.1 + cZ;
    t = chamber->conduction_electron_bank[i].time;
    e = 0.;
    dx = dy = dz = 0.;
 
  } else {
    // Make sure a delta electron with this number exists.
    if (i < 0 || i >= nDeltas) {
      std::cerr << className << "::GetElectron:\n";
      std::cerr << "    Delta electron number out of range.\n";
      return false;
    }
 
    x = deltaElectrons[i].x;
    y = deltaElectrons[i].y;
    z = deltaElectrons[i].z;
    t = deltaElectrons[i].t;
    e = deltaElectrons[i].e;
    dx = deltaElectrons[i].dx;
    dy = deltaElectrons[i].dy;
    dz = deltaElectrons[i].dz;
  }
 
  return true;
}

Referenced by GarfieldPhysics::DoIt().

GetFanoFactor()

double Garfield::TrackHeed::GetFanoFactor

(

)

const

Definition at line 1371 of file TrackHeed.cc.

{ return matter->F; }

GetStoppingPower()

double Garfield::TrackHeed::GetStoppingPower ( )

virtual

Reimplemented from Garfield::Track.

Definition at line 366 of file TrackHeed.cc.

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

GetW()

double Garfield::TrackHeed::GetW

(

)

const

Definition at line 1370 of file TrackHeed.cc.

{ return 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

Implements Garfield::Track.

Definition at line 158 of file TrackHeed.cc.

                                           {
 
  hasActiveTrack = false;
  ready = false;
 
  // Make sure the sensor has been set.
  if (sensor == 0) {
    std::cerr << className << "::NewTrack:\n";
    std::cerr << "    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 (!sensor->GetArea(xmin, ymin, zmin, xmax, ymax, zmax)) {
    std::cerr << className << "::NewTrack:\n";
    std::cerr << "    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 (debug) {
    std::cout << className << "::NewTrack:\n";
    std::cout << "    Bounding box dimensions:\n";
    std::cout << "      x: " << lx << " cm\n";
    std::cout << "      y: " << ly << " cm\n";
    std::cout << "      z: " << lz << " cm\n";
  }
  if (fabs(lx - lX) > Small || fabs(ly - lY) > Small || fabs(lz - lZ) > Small) {
    lX = lx;
    lY = ly;
    lZ = lz;
    isChanged = true;
  }
  // Update the center of the bounding box.
  if (std::isinf(xmin) || std::isinf(xmax)) {
    cX = 0.;
  } else {
    cX = 0.5 * (xmin + xmax);
  }
  if (std::isinf(ymin) || std::isinf(ymax)) {
    cY = 0.;
  } else {
    cY = 0.5 * (ymin + ymax);
  }
  if (std::isinf(zmin) || std::isinf(zmax)) {
    cZ = 0.;
  } else {
    cZ = 0.5 * (zmin + zmax);
  }
  if (debug) {
    std::cout << className << "::NewTrack:\n";
    std::cout << "    Center of bounding box:\n";
    std::cout << "      x: " << cX << " cm\n";
    std::cout << "      y: " << cY << " cm\n";
    std::cout << "      z: " << cZ << " cm\n";
  }
 
  HeedInterface::sensor = sensor;
 
  // Make sure the initial position is inside an ionisable medium.
  Medium* medium;
  if (!sensor->GetMedium(x0, y0, z0, medium)) {
    std::cerr << className << "::NewTrack:\n";
    std::cerr << "    No medium at initial position.\n";
    return false;
  } else if (!medium->IsIonisable()) {
    std::cerr << "TrackHeed:NewTrack:\n";
    std::cerr << "    Medium at initial position is not ionisable.\n";
    return false;
  }
 
  // Check if the medium has changed since the last call.
  if (medium->GetName() != mediumName ||
      fabs(medium->GetMassDensity() - mediumDensity) > 1.e-9) {
    isChanged = true;
  }
 
  // If medium, particle or bounding box have changed,
  // update the cross-sections.
  if (isChanged) {
    if (!Setup(medium)) return false;
    isChanged = false;
    mediumName = medium->GetName();
    mediumDensity = medium->GetMassDensity();
  }
 
  Heed::particle_bank.clear();
  deltaElectrons.clear();
  Heed::cluster_bank.allocate_block(100);
  chamber->conduction_electron_bank.allocate_block(1000);
 
  // 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 (debug) {
      std::cout << className << "::NewTrack:\n";
      std::cout << "    Direction vector has zero norm.\n";
      std::cout << "    Initial direction is randomized.\n";
    }
    // Null vector. Sample the direction isotropically.
    const double ctheta = 1. - 2. * RndmUniform();
    const double stheta = sqrt(1. - ctheta * ctheta);
    const double phi = TwoPi * RndmUniform();
    dx = cos(phi) * stheta;
    dy = sin(phi) * stheta;
    dz = ctheta;
  } else {
    // Normalise the direction vector.
    dx /= d;
    dy /= d;
    dz /= d;
  }
  vec velocity(dx, dy, dz);
  velocity = velocity * Heed::c_light * GetBeta();
 
  if (debug) {
    std::cout << className << "::NewTrack:\n";
    std::cout << "    Track starts at (" << x0 << ", " << y0 << ", " << z0
              << ") at time " << t0 << "\n";
    std::cout << "    Initial direction: (" << dx << ", " << dy << ", " << dz
              << ")\n";
  }
 
  // Initial position (shift with respect to bounding box center and
  // convert from cm to mm).
  point p0((x0 - cX) * 10., (y0 - cY) * 10., (z0 - cZ) * 10.);
  // Setup the particle.
  Heed::last_particle_number = 0;
  if (particle != 0) {
    delete particle;
    particle = 0;
  }
 
  Heed::particle_def* particleType = &Heed::muon_minus_def;
  if (particleName == "e-") {
    particleType = &Heed::electron_def;
  } else if (particleName == "e+") {
    particleType = &Heed::positron_def;
  } else if (particleName == "mu-") {
    particleType = &Heed::muon_minus_def;
  } else if (particleName == "mu+") {
    particleType = &Heed::muon_plus_def;
  } else if (particleName == "pi-") {
    particleType = &Heed::pi_minus_meson_def;
  } else if (particleName == "pi+") {
    particleType = &Heed::pi_plus_meson_def;
  } else if (particleName == "K-") {
    particleType = &Heed::K_minus_meson_def;
  } else if (particleName == "K+") {
    particleType = &Heed::K_plus_meson_def;
  } else if (particleName == "p") {
    particleType = &Heed::proton_def;
  } else if (particleName == "pbar") {
    particleType = &Heed::anti_proton_def;
  } else if (particleName == "d") {
    particleType = &Heed::deuteron_def;
  } else if (particleName == "alpha") {
    particleType = &Heed::alpha_particle_def;
  } else if (particleName == "exotic") {
    // User defined particle
    Heed::user_particle_def.set_mass(mass * 1.e-6);
    Heed::user_particle_def.set_charge(q);
    particleType = &Heed::user_particle_def;
  } else {
    // Not a predefined particle, use muon definition.
    if (q > 0.) {
      particleType = &Heed::muon_minus_def;
    } else {
      particleType = &Heed::muon_plus_def;
    }
  }
 
  particle = new Heed::HeedParticle(chamber, p0, velocity, t0, particleType);
  // Transport the particle.
  particle->fly();
  hasActiveTrack = true;
  ready = true;
 
  // Plot the new track.
  if (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 950 of file TrackHeed.cc.

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

SetParticleUser()

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

Definition at line 973 of file TrackHeed.cc.

                                                              {
 
  if (fabs(z) < Small) {
    std::cerr << className << "::SetParticleUser:\n";
    std::cerr << "    Particle cannot have zero charge.\n";
    return;
  }
  if (m < Small) {
    std::cerr << className << "::SetParticleUser:\n";
    std::cerr << "    Particle mass must be greater than zero.\n";
  }
  q = z;
  mass = m;
  isElectron = false;
  spin = 0;
  particleName = "exotic";
}

TransportDeltaElectron()

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 615 of file TrackHeed.cc.

                                                 {
 
  nel = 0;
 
  // Check if delta electron transport was disabled.
  if (!useDelta) {
    std::cerr << className << "::TransportDeltaElectron:\n";
    std::cerr << "    Delta electron transport has been switched off.\n";
    return;
  }
 
  // Make sure the kinetic energy is positive.
  if (e0 <= 0.) {
    std::cerr << className << "::TransportDeltaElectron:\n";
    std::cerr << "    Kinetic energy must be positive.\n";
    return;
  }
 
  // Make sure the sensor has been set.
  if (sensor == 0) {
    std::cerr << className << "::TransportDeltaElectron:\n";
    std::cerr << "    Sensor is not defined.\n";
    ready = false;
    return;
  }
 
  // Get the bounding box.
  double xmin, ymin, zmin;
  double xmax, ymax, zmax;
  if (!sensor->GetArea(xmin, ymin, zmin, xmax, ymax, zmax)) {
    std::cerr << className << "::TransportDeltaElectron:\n";
    std::cerr << "    Drift area is not set.\n";
    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 - lX) > Small || fabs(ly - lY) > Small || fabs(lz - lZ) > Small) {
    lX = lx;
    lY = ly;
    lZ = lz;
    isChanged = true;
    update = true;
    hasActiveTrack = false;
  }
  // Update the center of the bounding box.
  cX = 0.5 * (xmin + xmax);
  cY = 0.5 * (ymin + ymax);
  cZ = 0.5 * (zmin + zmax);
 
  HeedInterface::sensor = sensor;
 
  // Make sure the initial position is inside an ionisable medium.
  Medium* medium;
  if (!sensor->GetMedium(x0, y0, z0, medium)) {
    std::cerr << className << "::TransportDeltaElectron:\n";
    std::cerr << "    No medium at initial position.\n";
    return;
  } else if (!medium->IsIonisable()) {
    std::cerr << "TrackHeed:TransportDeltaElectron:\n";
    std::cerr << "    Medium at initial position is not ionisable.\n";
    ready = false;
    return;
  }
 
  // Check if the medium has changed since the last call.
  if (medium->GetName() != mediumName ||
      fabs(medium->GetMassDensity() - mediumDensity) > 1.e-9) {
    isChanged = true;
    update = true;
    ready = false;
    hasActiveTrack = false;
  }
 
  // If medium or bounding box have changed, update the "chamber".
  if (update) {
    if (!Setup(medium)) return;
    ready = true;
    mediumName = medium->GetName();
    mediumDensity = medium->GetMassDensity();
  }
 
  deltaElectrons.clear();
  chamber->conduction_electron_bank.allocate_block(1000);
 
  // 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.
    const double phi = TwoPi * RndmUniform();
    const double ctheta = 1. - 2. * RndmUniform();
    const double stheta = sqrt(1. - ctheta * ctheta);
    dx = cos(phi) * stheta;
    dy = sin(phi) * stheta;
    dz = ctheta;
  } else {
    // Normalise the direction vector.
    dx /= d;
    dy /= d;
    dz /= d;
  }
  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::c_light * beta;
  velocity = velocity * speed;
 
  // Initial position (shift with respect to bounding box center and
  // convert from cm to mm).
  point p0((x0 - cX) * 10., (y0 - cY) * 10., (z0 - cZ) * 10.);
 
  // Transport the electron.
  Heed::HeedDeltaElectron delta(chamber, p0, velocity, t0, 0);
  delta.fly();
 
  nel = chamber->conduction_electron_bank.get_qel();
}

Referenced by GarfieldPhysics::DoIt().

TransportPhoton()

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 743 of file TrackHeed.cc.

                                                                              {
 
  nel = 0;
 
  // Make sure the energy is positive.
  if (e0 <= 0.) {
    std::cerr << className << "::TransportPhoton:\n";
    std::cerr << "    Photon energy must be positive.\n";
    return;
  }
 
  // Make sure the sensor has been set.
  if (sensor == 0) {
    std::cerr << className << "::TransportPhoton:\n";
    std::cerr << "    Sensor is not defined.\n";
    ready = false;
    return;
  }
 
  // Get the bounding box.
  double xmin, ymin, zmin;
  double xmax, ymax, zmax;
  if (!sensor->GetArea(xmin, ymin, zmin, xmax, ymax, zmax)) {
    std::cerr << className << "::TransportPhoton:\n";
    std::cerr << "    Drift area is not set.\n";
    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 - lX) > Small || fabs(ly - lY) > Small || fabs(lz - lZ) > Small) {
    lX = lx;
    lY = ly;
    lZ = lz;
    isChanged = true;
    update = true;
    hasActiveTrack = false;
  }
  // Update the center of the bounding box.
  cX = 0.5 * (xmin + xmax);
  cY = 0.5 * (ymin + ymax);
  cZ = 0.5 * (zmin + zmax);
 
  HeedInterface::sensor = sensor;
 
  // Make sure the initial position is inside an ionisable medium.
  Medium* medium;
  if (!sensor->GetMedium(x0, y0, z0, medium)) {
    std::cerr << className << "::TransportPhoton:\n";
    std::cerr << "    No medium at initial position.\n";
    return;
  } else if (!medium->IsIonisable()) {
    std::cerr << "TrackHeed:TransportPhoton:\n";
    std::cerr << "    Medium at initial position is not ionisable.\n";
    ready = false;
    return;
  }
 
  // Check if the medium has changed since the last call.
  if (medium->GetName() != mediumName ||
      fabs(medium->GetMassDensity() - mediumDensity) > 1.e-9) {
    isChanged = true;
    update = true;
    ready = false;
  }
 
  // If medium or bounding box have changed, update the "chamber".
  if (update) {
    if (!Setup(medium)) return;
    ready = true;
    mediumName = medium->GetName();
    mediumDensity = medium->GetMassDensity();
  }
 
  // Delete the particle bank.
  // Clusters from the current track will be lost.
  hasActiveTrack = false;
  Heed::last_particle_number = 0;
  Heed::particle_bank.clear();
  deltaElectrons.clear();
  nDeltas = 0;
  chamber->conduction_electron_bank.allocate_block(1000);
 
  // 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.
    const double phi = TwoPi * RndmUniform();
    const double ctheta = 1. - 2. * RndmUniform();
    const double stheta = sqrt(1. - ctheta * ctheta);
    dx = cos(phi) * stheta;
    dy = sin(phi) * stheta;
    dz = ctheta;
  } else {
    // Normalise the direction vector.
    dx /= d;
    dy /= d;
    dz /= d;
  }
  vec velocity(dx, dy, dz);
  velocity = velocity * Heed::c_light;
 
  // Initial position (shift with respect to bounding box center and
  // convert from cm to mm).
  point p0((x0 - cX) * 10., (y0 - cY) * 10., (z0 - cZ) * 10.);
 
  // Create and transport the photon.
  Heed::HeedPhoton photon(chamber, p0, velocity, t0, 0, e0 * 1.e-6, 0);
  photon.fly();
 
  // Make a list of parent particle id numbers.
  std::vector<int> ids;
  ids.clear();
  // At the beginning, there is only the original photon.
  ids.push_back(photon.particle_number);
  int nIds = 1;
 
  // Look for daughter particles.
  Heed::HeedDeltaElectron* delta = 0;
  Heed::HeedPhoton* fluorescencePhoton = 0;
 
  // Get the first element from the particle bank.
  AbsListNode<ActivePtr<gparticle> >* nextNode =
      Heed::particle_bank.get_first_node();
  AbsListNode<ActivePtr<gparticle> >* tempNode = 0;
  // Loop over the particle bank.
  while (nextNode != 0) {
    // Check if it is a delta electron.
    delta = dynamic_cast<Heed::HeedDeltaElectron*>(nextNode->el.get());
    if (delta != 0) {
      // Check if the delta electron was produced by one of the photons
      // belonging to this cluster.
      bool gotParent = false;
      for (int i = nIds; i--;) {
        if (delta->parent_particle_number == ids[i]) {
          gotParent = true;
          if (useDelta) {
            // Transport the delta electron.
            delta->fly();
          } else {
            // Add the delta electron to the list, for later use.
            deltaElectron newDeltaElectron;
            newDeltaElectron.x = delta->currpos.pt.v.x * 0.1 + cX;
            newDeltaElectron.y = delta->currpos.pt.v.y * 0.1 + cY;
            newDeltaElectron.z = delta->currpos.pt.v.z * 0.1 + 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;
            deltaElectrons.push_back(newDeltaElectron);
            ++nDeltas;
          }
          break;
        }
      }
      if (!gotParent) {
        std::cerr << className << "::TransportPhoton:\n";
        std::cerr << "    Delta electron with unknown parent.\n";
      }
    } else {
      // Check if it is a fluorescence photon.
      fluorescencePhoton = dynamic_cast<Heed::HeedPhoton*>(nextNode->el.get());
      if (fluorescencePhoton == 0) {
        std::cerr << className << "::TransportPhoton:\n";
        std::cerr << "    Unknown secondary particle.\n";
        return;
      }
      for (int i = nIds; i--;) {
        if (fluorescencePhoton->parent_particle_number == ids[i]) {
          // Transport the photon and add its number to the list of ids.
          fluorescencePhoton->fly();
          ids.push_back(fluorescencePhoton->particle_number);
          ++nIds;
          break;
        }
      }
    }
    // Proceed with the next node in the particle bank.
    tempNode = nextNode->get_next_node();
    Heed::particle_bank.erase(nextNode);
    nextNode = tempNode;
  }
 
  // Get the total number of electrons produced in this step.
  if (useDelta) {
    nel = chamber->conduction_electron_bank.get_qel();
  } else {
    nel = nDeltas;
  }
}

Referenced by GarfieldPhysics::DoIt().

---

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

• TrackHeed.hh
• TrackHeed.cc