Garfield++ 4.0
A toolkit for the detailed simulation of particle detectors based on ionisation measurement in gases and semiconductors
Loading...
Searching...
No Matches
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.
 
bool NewTrack (const double x0, const double y0, const double z0, const double t0, const double dx0, const double dy0, const double dz0) override
 
bool GetCluster (double &xcls, double &ycls, double &zcls, double &tcls, int &n, double &e, double &extra) override
 
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
 
double GetClusterDensity () override
 
double GetStoppingPower () override
 Get the stopping power (mean energy loss [eV] per cm).
 
double GetW () const
 Return the W value of the medium (of the last simulated track).
 
double GetFanoFactor () const
 Return the Fano factor of the medium (of the last simulated track).
 
double GetPhotoAbsorptionCrossSection (const double e) const
 Return the photoabsorption cross-section at a given energy.
 
bool Initialise (Medium *medium)
 Compute the differential cross-section for a given medium.
 
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 &ne, 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 &ne)
 
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 &ne, 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 &ne)
 
void EnableElectricField ()
 Take the electric field into account in the stepping algorithm.
 
void DisableElectricField ()
 Do not take the electric field into account in the stepping algorithm.
 
void EnableMagneticField ()
 Take the magnetic field into account in the stepping algorithm.
 
void DisableMagneticField ()
 Do not take the magnetic field into account in the stepping algorithm.
 
void SetSteppingLimits (const double maxStep, const double radStraight, const double stepAngleStraight, const double stepAngleCurved)
 
void GetSteppingLimits (double &maxStep, double &radStraight, double &stepAngleStraight, double &stepAngleCurved)
 
void EnableDeltaElectronTransport ()
 Switch simulation of delta electrons on.
 
void DisableDeltaElectronTransport ()
 Switch simulation of delta electrons off.
 
void EnablePhotonReabsorption (const bool on=true)
 Simulate (or not) the photons produced in the atomic relaxation cascade.
 
void EnablePhotoAbsorptionCrossSectionOutput (const bool on)
 Write the photoabsorption cross-sections used to a text file.
 
void SetEnergyMesh (const double e0, const double e1, const int nsteps)
 
void SetParticleUser (const double m, const double z)
 
void EnableOneStepFly (const bool on)
 
- Public Member Functions inherited from Garfield::Track
 Track ()
 Constructor.
 
virtual ~Track ()
 Destructor.
 
virtual void SetParticle (const std::string &part)
 
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 ( $\beta = v/c$) 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
 Return the particle energy.
 
double GetBetaGamma () const
 Return the $\beta\gamma$ of the projectile.
 
double GetBeta () const
 Return the speed ( $\beta = v/c$) of the projectile.
 
double GetGamma () const
 Return the Lorentz factor of the projectile.
 
double GetMomentum () const
 Return the particle momentum.
 
double GetKineticEnergy () const
 Return the kinetic energy of the projectile.
 
double GetCharge () const
 Get the charge of the projectile.
 
double GetMass () const
 Get the mass [eV / c2] of the projectile.
 
void SetSensor (Sensor *s)
 Set the sensor through which to transport the particle.
 
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)
 Switch on plotting.
 
void DisablePlotting ()
 Switch off plotting.
 
void EnableDebugging ()
 Switch on debugging messages.
 
void DisableDebugging ()
 Switch off debugging messages.
 

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 = "Track"
 
double m_q = -1.
 
int m_spin = 1
 
double m_mass
 
double m_energy = 0.
 
double m_beta2
 
bool m_isElectron = false
 
std::string m_particleName = "mu-"
 
Sensorm_sensor = nullptr
 
bool m_isChanged = true
 
ViewDriftm_viewer = nullptr
 
bool m_debug = false
 
size_t m_plotId = 0
 

Detailed Description

Generate tracks using Heed++.

Definition at line 35 of file TrackHeed.hh.

Constructor & Destructor Documentation

◆ TrackHeed()

Garfield::TrackHeed::TrackHeed ( )

Constructor.

Definition at line 49 of file TrackHeed.cc.

49 : Track() {
50 m_className = "TrackHeed";
51 m_conductionElectrons.reserve(1000);
52 m_conductionIons.reserve(1000);
53
54 m_fieldMap.reset(new Heed::HeedFieldMap());
55}
std::string m_className
Definition: Track.hh:102
Track()
Constructor.
Definition: Track.cc:13
Retrieve electric and magnetic field from Sensor.
Definition: HeedFieldMap.h:15

◆ ~TrackHeed()

Garfield::TrackHeed::~TrackHeed ( )
virtual

Destructor.

Definition at line 57 of file TrackHeed.cc.

57{}

Member Function Documentation

◆ DisableDeltaElectronTransport()

void Garfield::TrackHeed::DisableDeltaElectronTransport ( )
inline

Switch simulation of delta electrons off.

Definition at line 161 of file TrackHeed.hh.

161{ m_doDeltaTransport = false; }

◆ DisableElectricField()

void Garfield::TrackHeed::DisableElectricField ( )

Do not take the electric field into account in the stepping algorithm.

Definition at line 677 of file TrackHeed.cc.

677{ m_fieldMap->UseEfield(false); }

◆ DisableMagneticField()

void Garfield::TrackHeed::DisableMagneticField ( )

Do not take the magnetic field into account in the stepping algorithm.

Definition at line 679 of file TrackHeed.cc.

679{ m_fieldMap->UseBfield(false); }

◆ EnableDeltaElectronTransport()

void Garfield::TrackHeed::EnableDeltaElectronTransport ( )
inline

Switch simulation of delta electrons on.

Definition at line 159 of file TrackHeed.hh.

159{ m_doDeltaTransport = true; }

Referenced by GarfieldPhysics::InitializePhysics().

◆ EnableElectricField()

void Garfield::TrackHeed::EnableElectricField ( )

Take the electric field into account in the stepping algorithm.

Definition at line 676 of file TrackHeed.cc.

676{ m_fieldMap->UseEfield(true); }

◆ EnableMagneticField()

void Garfield::TrackHeed::EnableMagneticField ( )

Take the magnetic field into account in the stepping algorithm.

Definition at line 678 of file TrackHeed.cc.

678{ m_fieldMap->UseBfield(true); }

Referenced by main().

◆ EnableOneStepFly()

void Garfield::TrackHeed::EnableOneStepFly ( const bool  on)
inline

Definition at line 182 of file TrackHeed.hh.

182{ m_useOneStepFly=on; }

◆ EnablePhotoAbsorptionCrossSectionOutput()

void Garfield::TrackHeed::EnablePhotoAbsorptionCrossSectionOutput ( const bool  on)
inline

Write the photoabsorption cross-sections used to a text file.

Definition at line 169 of file TrackHeed.hh.

169 {
170 m_usePacsOutput = on;
171 }

◆ EnablePhotonReabsorption()

void Garfield::TrackHeed::EnablePhotonReabsorption ( const bool  on = true)
inline

Simulate (or not) the photons produced in the atomic relaxation cascade.

Definition at line 164 of file TrackHeed.hh.

164 {
165 m_usePhotonReabsorption = on;
166 }

◆ GetCluster() [1/2]

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

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

Parameters
xcls,ycls,zclscoordinates of the collision
tclstime of the collision
nnumber of electrons produced
edeposited energy
extraadditional information (not always implemented)

Implements Garfield::Track.

Definition at line 222 of file TrackHeed.cc.

223 {
224 int ni = 0;
225 return GetCluster(xcls, ycls, zcls, tcls, n, ni, e, extra);
226}
bool GetCluster(double &xcls, double &ycls, double &zcls, double &tcls, int &n, double &e, double &extra) override
Definition: TrackHeed.cc:222

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

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

230 {
231 // Initialise and reset.
232 xcls = ycls = zcls = tcls = 0.;
233 extra = 0.;
234 ne = ni = 0;
235 e = 0.;
236
237 m_deltaElectrons.clear();
238 m_conductionElectrons.clear();
239 m_conductionIons.clear();
240
241 // Make sure NewTrack has been called successfully.
242 if (!m_ready) {
243 std::cerr << m_className << "::GetCluster:\n"
244 << " Track has not been initialized. Call NewTrack first.\n";
245 return false;
246 }
247
248 if (m_particleBank.empty()) return false;
249 std::vector<Heed::gparticle*>::const_iterator end = m_particleBank.end();
250 if (m_bankIterator == end) return false;
251
252 // Look for the next cluster (i. e. virtual photon) in the list.
253 Heed::HeedPhoton* virtualPhoton = nullptr;
254 for (; m_bankIterator != end; ++m_bankIterator) {
255 // Convert the particle to a (virtual) photon.
256 virtualPhoton = dynamic_cast<Heed::HeedPhoton*>(*m_bankIterator);
257 if (!virtualPhoton) {
258 std::cerr << m_className << "::GetCluster:\n"
259 << " Particle is not a virtual photon. Program bug!\n";
260 // Try the next element.
261 continue;
262 }
263
264 // Get the location of the interaction (convert from mm to cm
265 // and shift with respect to bounding box center).
266 xcls = virtualPhoton->position().x * 0.1 + m_cX;
267 ycls = virtualPhoton->position().y * 0.1 + m_cY;
268 zcls = virtualPhoton->position().z * 0.1 + m_cZ;
269 tcls = virtualPhoton->time();
270 // Skip clusters outside the drift area or outside the active medium.
271 if (!IsInside(xcls, ycls, zcls)) continue;
272 // Add the first ion (at the position of the cluster).
273 m_conductionIons.emplace_back(
274 Heed::HeedCondElectron(Heed::point(virtualPhoton->position()), tcls));
275 ++m_bankIterator;
276 break;
277 }
278
279 // Stop if we did not find a virtual photon.
280 if (!virtualPhoton) return false;
281 // Plot the cluster, if requested.
282 if (m_viewer) PlotCluster(xcls, ycls, zcls);
283
284 std::vector<Heed::gparticle*> secondaries;
285 // Transport the virtual photon.
286 virtualPhoton->fly(secondaries);
287 // Get the transferred energy (convert from MeV to eV).
288 e = virtualPhoton->m_energy * 1.e6;
289
290 while (!secondaries.empty()) {
291 std::vector<Heed::gparticle*> newSecondaries;
292 // Loop over the secondaries.
293 for (auto secondary : secondaries) {
294 // Check if it is a delta electron.
295 auto delta = dynamic_cast<Heed::HeedDeltaElectron*>(secondary);
296 if (delta) {
297 extra += delta->kinetic_energy() * 1.e6;
298 const double x = delta->position().x * 0.1 + m_cX;
299 const double y = delta->position().y * 0.1 + m_cY;
300 const double z = delta->position().z * 0.1 + m_cZ;
301 if (!IsInside(x, y, z)) continue;
302 if (m_doDeltaTransport) {
303 // Transport the delta electron.
304 delta->fly(newSecondaries);
305 // Add the conduction electrons and ions to the list.
306 m_conductionElectrons.insert(m_conductionElectrons.end(),
307 delta->conduction_electrons.begin(),
308 delta->conduction_electrons.end());
309 m_conductionIons.insert(m_conductionIons.end(),
310 delta->conduction_ions.begin(),
311 delta->conduction_ions.end());
312 } else {
313 // Add the delta electron to the list, for later use.
314 deltaElectron newDeltaElectron;
315 newDeltaElectron.x = delta->position().x * 0.1 + m_cX;
316 newDeltaElectron.y = delta->position().y * 0.1 + m_cY;
317 newDeltaElectron.z = delta->position().z * 0.1 + m_cZ;
318 newDeltaElectron.t = delta->time();
319 newDeltaElectron.e = delta->kinetic_energy() * 1.e6;
320 newDeltaElectron.dx = delta->direction().x;
321 newDeltaElectron.dy = delta->direction().y;
322 newDeltaElectron.dz = delta->direction().z;
323 m_deltaElectrons.push_back(std::move(newDeltaElectron));
324 }
325 continue;
326 }
327 // Check if it is a real photon.
328 auto photon = dynamic_cast<Heed::HeedPhoton*>(secondary);
329 if (!photon) {
330 std::cerr << m_className << "::GetCluster:\n"
331 << " Particle is neither an electron nor a photon.\n";
332 continue;
333 }
334 extra += photon->m_energy * 1.e6;
335 const double x = photon->position().x * 0.1 + m_cX;
336 const double y = photon->position().y * 0.1 + m_cY;
337 const double z = photon->position().z * 0.1 + m_cZ;
338 if (!IsInside(x, y, z)) continue;
339 // Transport the photon.
340 if (m_usePhotonReabsorption) photon->fly(newSecondaries);
341 }
342 for (auto secondary : secondaries)
343 if (secondary) delete secondary;
344 secondaries.clear();
345 secondaries.swap(newSecondaries);
346 }
347 // Get the total number of electrons produced in this step.
348 ne = m_doDeltaTransport ? m_conductionElectrons.size()
349 : m_deltaElectrons.size();
350 ni = m_conductionIons.size();
351 return true;
352}
ViewDrift * m_viewer
Definition: Track.hh:116
void PlotCluster(const double x0, const double y0, const double z0)
Definition: Track.cc:192
double m_energy
Photon energy [MeV].
Definition: HeedPhoton.h:36
vfloat time() const
Get the current time of the particle.
Definition: gparticle.h:188
const vec & position() const
Get the current position of the particle.
Definition: gparticle.h:186
virtual void fly(std::vector< gparticle * > &secondaries)
Transport the particle.
Definition: gparticle.h:155
double kinetic_energy() const
Get the current kinetic energy.
Definition: mparticle.h:34
Point.
Definition: vec.h:368
vfloat x
Definition: vec.h:192
vfloat z
Definition: vec.h:194
vfloat y
Definition: vec.h:193

◆ GetClusterDensity()

double Garfield::TrackHeed::GetClusterDensity ( )
overridevirtual

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

Reimplemented from Garfield::Track.

Definition at line 200 of file TrackHeed.cc.

200 {
201
202 if (!m_transferCs) {
203 std::cerr << m_className << "::GetClusterDensity:\n"
204 << " Ionisation cross-section is not available.\n";
205 return 0.;
206 }
207
208 return m_transferCs->quanC;
209}

◆ 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
iindex of the electron
x,y,zcoordinates of the electron
ttime
ekinetic energy (only meaningful for delta-electrons)
dx,dy,dzdirection vector (only meaningful for delta-electrons)

Definition at line 354 of file TrackHeed.cc.

356 {
357 // Make sure NewTrack has successfully been called.
358 if (!m_ready) {
359 std::cerr << m_className << "::GetElectron:\n"
360 << " Track has not been initialized. Call NewTrack first.\n";
361 return false;
362 }
363
364 if (m_doDeltaTransport) {
365 // Make sure an electron with this number exists.
366 if (i >= m_conductionElectrons.size()) {
367 std::cerr << m_className << "::GetElectron: Index out of range.\n";
368 return false;
369 }
370
371 x = m_conductionElectrons[i].x * 0.1 + m_cX;
372 y = m_conductionElectrons[i].y * 0.1 + m_cY;
373 z = m_conductionElectrons[i].z * 0.1 + m_cZ;
374 t = m_conductionElectrons[i].time;
375 e = 0.;
376 dx = dy = dz = 0.;
377
378 } else {
379 // Make sure a delta electron with this number exists.
380 if (i >= m_deltaElectrons.size()) {
381 std::cerr << m_className << "::GetElectron:\n"
382 << " Delta electron number out of range.\n";
383 return false;
384 }
385
386 x = m_deltaElectrons[i].x;
387 y = m_deltaElectrons[i].y;
388 z = m_deltaElectrons[i].z;
389 t = m_deltaElectrons[i].t;
390 e = m_deltaElectrons[i].e;
391 dx = m_deltaElectrons[i].dx;
392 dy = m_deltaElectrons[i].dy;
393 dz = m_deltaElectrons[i].dz;
394 }
395 return true;
396}

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

◆ GetFanoFactor()

double Garfield::TrackHeed::GetFanoFactor ( ) const

Return the Fano factor of the medium (of the last simulated track).

Definition at line 1048 of file TrackHeed.cc.

1048{ 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
iindex of the ion
x,y,zcoordinates of the ion
ttime

Definition at line 398 of file TrackHeed.cc.

399 {
400 // Make sure a "conduction" ion with this number exists.
401 if (i >= m_conductionIons.size()) {
402 std::cerr << m_className << "::GetIon: Index out of range.\n";
403 return false;
404 }
405
406 x = m_conductionIons[i].x * 0.1 + m_cX;
407 y = m_conductionIons[i].y * 0.1 + m_cY;
408 z = m_conductionIons[i].z * 0.1 + m_cZ;
409 t = m_conductionIons[i].time;
410 return true;
411}

◆ GetPhotoAbsorptionCrossSection()

double Garfield::TrackHeed::GetPhotoAbsorptionCrossSection ( const double  e) const

Return the photoabsorption cross-section at a given energy.

Definition at line 1050 of file TrackHeed.cc.

1050 {
1051
1052 if (!m_matter) return 0.;
1053 // Convert eV to MeV.
1054 const double e = 1.e-6 * en;
1055 double cs = 0.;
1056 const auto n = m_matter->apacs.size();
1057 for (size_t i = 0; i < n; ++i) {
1058 const double w = m_matter->matter->weight_quan(i);
1059 cs += m_matter->apacs[i]->get_ACS(e) * w;
1060 }
1061 // Convert Mbarn to cm-2.
1062 return cs * 1.e-18;
1063}

◆ GetSteppingLimits()

void Garfield::TrackHeed::GetSteppingLimits ( double &  maxStep,
double &  radStraight,
double &  stepAngleStraight,
double &  stepAngleCurved 
)
inline

Definition at line 150 of file TrackHeed.hh.

151 {
152 maxStep = m_maxStep;
153 radStraight = m_radStraight;
154 stepAngleStraight = m_stepAngleStraight;
155 stepAngleCurved = m_stepAngleCurved;
156 }

◆ GetStoppingPower()

double Garfield::TrackHeed::GetStoppingPower ( )
overridevirtual

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

Reimplemented from Garfield::Track.

Definition at line 211 of file TrackHeed.cc.

211 {
212
213 if (!m_transferCs) {
214 std::cerr << m_className << "::GetStoppingPower:\n"
215 << " Ionisation cross-section is not available.\n";
216 return 0.;
217 }
218
219 return m_transferCs->meanC1 * 1.e6;
220}

◆ GetW()

double Garfield::TrackHeed::GetW ( ) const

Return the W value of the medium (of the last simulated track).

Definition at line 1047 of file TrackHeed.cc.

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

Referenced by GarfieldPhysics::DoIt().

◆ Initialise()

bool Garfield::TrackHeed::Initialise ( Medium medium)

Compute the differential cross-section for a given medium.

Definition at line 720 of file TrackHeed.cc.

720 {
721 // Make sure the path to the Heed database is known.
722 std::string databasePath;
723 char* dbPath = std::getenv("HEED_DATABASE");
724 if (dbPath == NULL) {
725 // Try GARFIELD_HOME.
726 dbPath = std::getenv("GARFIELD_HOME");
727 if (dbPath == NULL) {
728 std::cerr << m_className << "::Initialise:\n"
729 << " Cannot retrieve database path "
730 << "(environment variables HEED_DATABASE and GARFIELD_HOME "
731 << "are not defined).\n Cannot proceed.\n";
732 return false;
733 }
734 databasePath = std::string(dbPath) + "/Heed/heed++/database";
735 } else {
736 databasePath = dbPath;
737 }
738 if (databasePath[databasePath.size() - 1] != '/') {
739 databasePath.append("/");
740 }
741
742 // Check once more that the medium exists.
743 if (!medium) {
744 std::cerr << m_className << "::Initialise: Null pointer.\n";
745 return false;
746 }
747
748 // Setup the energy mesh.
749 m_energyMesh.reset(new Heed::EnergyMesh(m_emin, m_emax, m_nEnergyIntervals));
750
751 if (medium->IsGas()) {
752 if (!SetupGas(medium)) return false;
753 } else {
754 if (!SetupMaterial(medium)) return false;
755 }
756
757 // Energy transfer cross-section
758 // Set a flag indicating whether the primary particle is an electron.
759 m_transferCs.reset(new Heed::EnTransfCS(1.e-6 * m_mass, GetGamma() - 1.,
760 m_isElectron, m_matter.get(),
761 long(m_q)));
762
763 if (!SetupDelta(databasePath)) return false;
764
765 if (m_debug) {
766 const double nc = m_transferCs->quanC;
767 const double dedx = m_transferCs->meanC * 1.e3;
768 const double dedx1 = m_transferCs->meanC1 * 1.e3;
769 const double w = m_matter->W * 1.e6;
770 const double f = m_matter->F;
771 const double minI = m_matter->min_ioniz_pot * 1.e6;
772 std::cout << m_className << "::Initialise:\n";
773 std::cout << " Cluster density: " << nc << " cm-1\n";
774 std::cout << " Stopping power (restricted): " << dedx << " keV/cm\n";
775 std::cout << " Stopping power (incl. tail): " << dedx1 << " keV/cm\n";
776 std::cout << " W value: " << w << " eV\n";
777 std::cout << " Fano factor: " << f << "\n";
778 std::cout << " Min. ionization potential: " << minI << " eV\n";
779 }
780
781 Heed::fixsyscoor primSys(Heed::point(0., 0., 0.), Heed::basis("primary"),
782 "primary");
783 m_chamber.reset(new HeedChamber(primSys, m_lX, m_lY, m_lZ,
784 *m_transferCs.get(), *m_deltaCs.get()));
785 m_fieldMap->SetSensor(m_sensor);
786 return true;
787}
Sensor * m_sensor
Definition: Track.hh:112
bool m_debug
Definition: Track.hh:118
double GetGamma() const
Return the Lorentz factor of the projectile.
Definition: Track.hh:56
bool m_isElectron
Definition: Track.hh:109
double m_q
Definition: Track.hh:104
double m_mass
Definition: Track.hh:106
Basis.
Definition: vec.h:313

Referenced by NewTrack(), TransportDeltaElectron(), and TransportPhoton().

◆ 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 
)
overridevirtual

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

Implements Garfield::Track.

Definition at line 59 of file TrackHeed.cc.

61 {
62 m_hasActiveTrack = false;
63 m_ready = false;
64
65 // Make sure the sensor has been set.
66 if (!m_sensor) {
67 std::cerr << m_className << "::NewTrack: Sensor is not defined.\n";
68 return false;
69 }
70
71 bool update = false;
72 if (!UpdateBoundingBox(update)) return false;
73
74 // Make sure the initial position is inside an ionisable medium.
75 Medium* medium = nullptr;
76 if (!m_sensor->GetMedium(x0, y0, z0, medium)) {
77 std::cerr << m_className << "::NewTrack:\n"
78 << " No medium at initial position.\n";
79 return false;
80 } else if (!medium->IsIonisable()) {
81 std::cerr << m_className << "::NewTrack:\n"
82 << " Medium at initial position is not ionisable.\n";
83 return false;
84 }
85
86 // Check if the medium has changed since the last call.
87 if (medium->GetName() != m_mediumName ||
88 fabs(medium->GetMassDensity() - m_mediumDensity) > 1.e-9) {
89 m_isChanged = true;
90 }
91
92 // If medium, particle or bounding box have changed,
93 // update the cross-sections.
94 if (m_isChanged) {
95 if (!Initialise(medium)) return false;
96 m_isChanged = false;
97 m_mediumName = medium->GetName();
98 m_mediumDensity = medium->GetMassDensity();
99 }
100
101 ClearParticleBank();
102
103 m_deltaElectrons.clear();
104 m_conductionElectrons.clear();
105 m_conductionIons.clear();
106
107 // Check the direction vector.
108 double dx = dx0, dy = dy0, dz = dz0;
109 const double d = sqrt(dx * dx + dy * dy + dz * dz);
110 if (d < Small) {
111 if (m_debug) {
112 std::cout << m_className << "::NewTrack:\n"
113 << " Direction vector has zero norm.\n"
114 << " Initial direction is randomized.\n";
115 }
116 // Null vector. Sample the direction isotropically.
117 RndmDirection(dx, dy, dz);
118 } else {
119 // Normalise the direction vector.
120 dx /= d;
121 dy /= d;
122 dz /= d;
123 }
124 Heed::vec velocity(dx, dy, dz);
125 velocity = velocity * Heed::CLHEP::c_light * GetBeta();
126
127 if (m_debug) {
128 std::cout << m_className << "::NewTrack:\n Track starts at (" << x0
129 << ", " << y0 << ", " << z0 << ") at time " << t0 << "\n"
130 << " Direction: (" << dx << ", " << dy << ", " << dz << ")\n";
131 }
132
133 // Initial position (shift with respect to bounding box center and
134 // convert from cm to mm).
135 Heed::point p0((x0 - m_cX) * 10., (y0 - m_cY) * 10., (z0 - m_cZ) * 10.);
136
137 // Setup the particle.
139 if (m_particleName == "e-") {
140 particleType = &Heed::electron_def;
141 } else if (m_particleName == "e+") {
142 particleType = &Heed::positron_def;
143 } else if (m_particleName == "mu-") {
144 particleType = &Heed::muon_minus_def;
145 } else if (m_particleName == "mu+") {
146 particleType = &Heed::muon_plus_def;
147 } else if (m_particleName == "pi-") {
148 particleType = &Heed::pi_minus_meson_def;
149 } else if (m_particleName == "pi+") {
150 particleType = &Heed::pi_plus_meson_def;
151 } else if (m_particleName == "K-") {
152 particleType = &Heed::K_minus_meson_def;
153 } else if (m_particleName == "K+") {
154 particleType = &Heed::K_plus_meson_def;
155 } else if (m_particleName == "p") {
156 particleType = &Heed::proton_def;
157 } else if (m_particleName == "pbar") {
158 particleType = &Heed::anti_proton_def;
159 } else if (m_particleName == "d") {
160 particleType = &Heed::deuteron_def;
161 } else if (m_particleName == "alpha") {
162 particleType = &Heed::alpha_particle_def;
163 } else if (m_particleName == "exotic") {
164 // User defined particle
167 particleType = &Heed::user_particle_def;
168 } else {
169 // Not a predefined particle, use muon definition.
170 if (m_q > 0.) {
171 particleType = &Heed::muon_minus_def;
172 } else {
173 particleType = &Heed::muon_plus_def;
174 }
175 }
176
177 Heed::HeedParticle particle(m_chamber.get(), p0, velocity, t0, particleType,
178 m_fieldMap.get());
179 // Set the step limits.
180 particle.set_step_limits(m_maxStep * Heed::CLHEP::cm,
181 m_radStraight * Heed::CLHEP::cm,
182 m_stepAngleStraight * Heed::CLHEP::rad,
183 m_stepAngleCurved * Heed::CLHEP::rad);
184 // Transport the particle.
185 if (m_useOneStepFly) {
186 particle.fly(m_particleBank, true);
187 } else {
188 particle.fly(m_particleBank);
189 }
190
191 m_bankIterator = m_particleBank.begin();
192 m_hasActiveTrack = true;
193 m_ready = true;
194
195 // Plot the new track.
196 if (m_viewer) PlotNewTrack(x0, y0, z0);
197 return true;
198}
bool GetMedium(const double x, const double y, const double z, Medium *&medium)
Get the medium at (x, y, z).
Definition: Sensor.cc:159
bool Initialise(Medium *medium)
Compute the differential cross-section for a given medium.
Definition: TrackHeed.cc:720
double GetBeta() const
Return the speed ( ) of the projectile.
Definition: Track.hh:54
bool m_isChanged
Definition: Track.hh:114
std::string m_particleName
Definition: Track.hh:110
void PlotNewTrack(const double x0, const double y0, const double z0)
Definition: Track.cc:187
void set_mass(const double m)
void set_charge(const double z)
Definition: vec.h:179
void RndmDirection(double &dx, double &dy, double &dz, const double length=1.)
Draw a random (isotropic) direction vector.
Definition: Random.hh:107
particle_def pi_minus_meson_def("pi_minus_meson", "pi-", 139.56755 *MeV/c_squared, -eplus, 0, 0, 0.0, spin_def(1.0, -1.0))
Definition: particle_def.h:121
particle_def pi_plus_meson_def("pi_plus_meson", "pi+", 139.56755 *MeV/c_squared, eplus, 0, 0, 0.0, spin_def(1.0, 1.0))
Definition: particle_def.h:119
particle_def muon_minus_def("muon_minus", "mu-", 105.658367 *MeV/c_squared, electron_charge, 1, 0, 0.5, spin_def(0.0, 0.0))
Definition: particle_def.h:108
particle_def alpha_particle_def("alpha_particle", "alpha", 3727.417 *MeV/c_squared, 2 *eplus, 0, 4, 0.0, spin_def(0.0, 0.0))
Definition: particle_def.h:127
particle_def anti_proton_def("", "p-", proton_def)
Definition: particle_def.h:111
particle_def deuteron_def("deuteron", "dtr", 1875.613 *MeV/c_squared, eplus, 0, 2, 0.0, spin_def(0.0, 0.0))
Definition: particle_def.h:126
particle_def proton_def("proton", "p+", proton_mass_c2/c_squared, eplus, 0, 1, 0.5, spin_def(0.5, 0.5))
Definition: particle_def.h:110
particle_def K_minus_meson_def("K_minus_meson_def", "K-", K_plus_meson_def)
Definition: particle_def.h:124
particle_def K_plus_meson_def("K_plus_meson_def", "K+", 493.677 *MeV/c_squared, 1, 0, 0, 0.0, spin_def(0.5, -0.5))
Definition: particle_def.h:123
particle_def user_particle_def("user_particle", "X", 139.56755 *MeV/c_squared, eplus, 0, 0, 0.0, spin_def(0.0, 0.0))
Definition: particle_def.h:130
DoubleAc fabs(const DoubleAc &f)
Definition: DoubleAc.h:615
particle_def muon_plus_def("muon_plus", "mu+", muon_minus_def)
Definition: particle_def.h:109
particle_def positron_def("positron", "e+", electron_def)
Definition: particle_def.h:107
particle_def electron_def("electron", "e-", electron_mass_c2/c_squared, electron_charge, 1, 0, 0.5, spin_def(0.0, 0.0))
Definition: particle_def.h:106
DoubleAc sqrt(const DoubleAc &f)
Definition: DoubleAc.cpp:314

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

◆ SetEnergyMesh()

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

Specify the energy mesh to be used.

Parameters
e0,e1lower/higher limit of the energy range [eV]
nstepsnumber of intervals

Definition at line 681 of file TrackHeed.cc.

682 {
683 if (fabs(e1 - e0) < Small) {
684 std::cerr << m_className << "::SetEnergyMesh:\n"
685 << " Invalid energy range:\n"
686 << " " << e0 << " < E [eV] < " << e1 << "\n";
687 return;
688 }
689
690 if (nsteps <= 0) {
691 std::cerr << m_className << "::SetEnergyMesh:\n"
692 << " Number of intervals must be > 0.\n";
693 return;
694 }
695
696 m_emin = std::min(e0, e1);
697 m_emax = std::max(e0, e1);
698 m_emin *= 1.e-6;
699 m_emax *= 1.e-6;
700 m_nEnergyIntervals = nsteps;
701}

◆ SetParticleUser()

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

Define particle mass and charge (for exotic particles). For standard particles Track::SetParticle should be used.

Definition at line 703 of file TrackHeed.cc.

703 {
704 if (fabs(z) < Small) {
705 std::cerr << m_className << "::SetParticleUser:\n"
706 << " Particle cannot have zero charge.\n";
707 return;
708 }
709 if (m < Small) {
710 std::cerr << m_className << "::SetParticleUser:\n"
711 << " Particle mass must be greater than zero.\n";
712 }
713 m_q = z;
714 m_mass = m;
715 m_isElectron = false;
716 m_spin = 0;
717 m_particleName = "exotic";
718}

◆ SetSteppingLimits()

void Garfield::TrackHeed::SetSteppingLimits ( const double  maxStep,
const double  radStraight,
const double  stepAngleStraight,
const double  stepAngleCurved 
)
inline

Set parameters for calculating the particle trajectory.

Parameters
maxStepmaximum step length
radStraightradius beyond which to approximate circles by polylines.
stepAngleStraightmax. angular step (in radian) when using polyline steps.
stepAngleCurvedmax. angular step (in radian) when using circular steps.

Definition at line 142 of file TrackHeed.hh.

144 {
145 m_maxStep = maxStep;
146 m_radStraight = radStraight;
147 m_stepAngleStraight = stepAngleStraight;
148 m_stepAngleCurved = stepAngleCurved;
149 }

◆ 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 &  ne 
)

Simulate a delta electron.

Parameters
x0,y0,z0initial position of the delta electron
t0initial time
e0initial kinetic energy of the delta electron
dx0,dy0,dz0initial direction of the delta electron
nenumber of electrons produced by the delta electron

Definition at line 413 of file TrackHeed.cc.

417 {
418 int ni = 0;
419 return TransportDeltaElectron(x0, y0, z0, t0, e0, dx0, dy0, dz0, nel, ni);
420}
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 &ne, int &ni)
Definition: TrackHeed.cc:422

◆ 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 &  ne,
int &  ni 
)

Simulate a delta electron.

Parameters
x0,y0,z0initial position of the delta electron
t0initial time
e0initial kinetic energy of the delta electron
dx0,dy0,dz0initial direction of the delta electron
ne,ninumber of electrons/ions produced by the delta electron

Definition at line 422 of file TrackHeed.cc.

426 {
427 nel = 0;
428 ni = 0;
429
430 // Check if delta electron transport was disabled.
431 if (!m_doDeltaTransport) {
432 std::cerr << m_className << "::TransportDeltaElectron:\n"
433 << " Delta electron transport has been switched off.\n";
434 return;
435 }
436
437 // Make sure the sensor has been set.
438 if (!m_sensor) {
439 std::cerr << m_className << "::TransportDeltaElectron:\n"
440 << " Sensor is not defined.\n";
441 m_ready = false;
442 return;
443 }
444
445 bool update = false;
446 if (!UpdateBoundingBox(update)) return;
447
448 // Make sure the initial position is inside an ionisable medium.
449 Medium* medium = nullptr;
450 if (!m_sensor->GetMedium(x0, y0, z0, medium)) {
451 std::cerr << m_className << "::TransportDeltaElectron:\n"
452 << " No medium at initial position.\n";
453 return;
454 } else if (!medium->IsIonisable()) {
455 std::cerr << "TrackHeed:TransportDeltaElectron:\n"
456 << " Medium at initial position is not ionisable.\n";
457 m_ready = false;
458 return;
459 }
460
461 // Check if the medium has changed since the last call.
462 if (medium->GetName() != m_mediumName ||
463 fabs(medium->GetMassDensity() - m_mediumDensity) > 1.e-9) {
464 m_isChanged = true;
465 update = true;
466 m_ready = false;
467 m_hasActiveTrack = false;
468 }
469
470 // If medium or bounding box have changed, update the "chamber".
471 if (update) {
472 if (!Initialise(medium)) return;
473 m_ready = true;
474 m_mediumName = medium->GetName();
475 m_mediumDensity = medium->GetMassDensity();
476 }
477
478 m_deltaElectrons.clear();
479 m_conductionElectrons.clear();
480 m_conductionIons.clear();
481
482 // Initial position (shift with respect to bounding box center and
483 // convert from cm to mm).
484 Heed::point p0((x0 - m_cX) * 10., (y0 - m_cY) * 10., (z0 - m_cZ) * 10.);
485
486 // Make sure the kinetic energy is positive.
487 if (e0 <= 0.) {
488 // Just create a conduction electron on the spot.
489 m_conductionElectrons.emplace_back(Heed::HeedCondElectron(p0, t0));
490 nel = 1;
491 return;
492 }
493
494 // Check the direction vector.
495 double dx = dx0, dy = dy0, dz = dz0;
496 const double d = sqrt(dx * dx + dy * dy + dz * dz);
497 if (d <= 0.) {
498 // Null vector. Sample the direction isotropically.
499 RndmDirection(dx, dy, dz);
500 } else {
501 // Normalise the direction vector.
502 dx /= d;
503 dy /= d;
504 dz /= d;
505 }
506 Heed::vec velocity(dx, dy, dz);
507
508 // Calculate the speed for the given kinetic energy.
509 const double gamma = 1. + e0 / ElectronMass;
510 const double beta = sqrt(1. - 1. / (gamma * gamma));
511 double speed = Heed::CLHEP::c_light * beta;
512 velocity = velocity * speed;
513
514 // Transport the electron.
515 std::vector<Heed::gparticle*> secondaries;
516 Heed::HeedDeltaElectron delta(m_chamber.get(), p0, velocity, t0, 0,
517 m_fieldMap.get());
518 delta.fly(secondaries);
519 ClearBank(secondaries);
520
521 m_conductionElectrons.swap(delta.conduction_electrons);
522 m_conductionIons.swap(delta.conduction_ions);
523 nel = m_conductionElectrons.size();
524 ni = m_conductionIons.size();
525}

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 &  ne 
)

Simulate a photon.

Parameters
x0,y0,z0initial position of the photon
t0initial time
e0initial energy of the photon
dx0,dy0,dz0initial direction of the photon
nenumber of electrons produced by the photon

Definition at line 527 of file TrackHeed.cc.

530 {
531 int ni = 0;
532 TransportPhoton(x0, y0, z0, t0, e0, dx0, dy0, dz0, nel, ni);
533}
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 &ne, int &ni)
Definition: TrackHeed.cc:535

◆ 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 &  ne,
int &  ni 
)

Simulate a photon.

Parameters
x0,y0,z0initial position of the photon
t0initial time
e0initial energy of the photon
dx0,dy0,dz0initial direction of the photon
ne,ninumber of electrons/ions produced by the photon

Definition at line 535 of file TrackHeed.cc.

539 {
540 nel = 0;
541 ni = 0;
542
543 // Make sure the energy is positive.
544 if (e0 <= 0.) {
545 std::cerr << m_className << "::TransportPhoton:\n"
546 << " Photon energy must be positive.\n";
547 return;
548 }
549
550 // Make sure the sensor has been set.
551 if (!m_sensor) {
552 std::cerr << m_className << "::TransportPhoton: Sensor is not defined.\n";
553 m_ready = false;
554 return;
555 }
556
557 bool update = false;
558 if (!UpdateBoundingBox(update)) return;
559
560 // Make sure the initial position is inside an ionisable medium.
561 Medium* medium = nullptr;
562 if (!m_sensor->GetMedium(x0, y0, z0, medium)) {
563 std::cerr << m_className << "::TransportPhoton:\n"
564 << " No medium at initial position.\n";
565 return;
566 } else if (!medium->IsIonisable()) {
567 std::cerr << "TrackHeed:TransportPhoton:\n"
568 << " Medium at initial position is not ionisable.\n";
569 m_ready = false;
570 return;
571 }
572
573 // Check if the medium has changed since the last call.
574 if (medium->GetName() != m_mediumName ||
575 fabs(medium->GetMassDensity() - m_mediumDensity) > 1.e-9) {
576 m_isChanged = true;
577 update = true;
578 m_ready = false;
579 }
580
581 // If medium or bounding box have changed, update the "chamber".
582 if (update) {
583 if (!Initialise(medium)) return;
584 m_ready = true;
585 m_mediumName = medium->GetName();
586 m_mediumDensity = medium->GetMassDensity();
587 }
588
589 // Delete the particle bank.
590 // Clusters from the current track will be lost.
591 m_hasActiveTrack = false;
592 ClearParticleBank();
593 m_deltaElectrons.clear();
594 m_conductionElectrons.clear();
595 m_conductionIons.clear();
596
597 // Check the direction vector.
598 double dx = dx0, dy = dy0, dz = dz0;
599 const double d = sqrt(dx * dx + dy * dy + dz * dz);
600 if (d <= 0.) {
601 // Null vector. Sample the direction isotropically.
602 RndmDirection(dx, dy, dz);
603 } else {
604 // Normalise the direction vector.
605 dx /= d;
606 dy /= d;
607 dz /= d;
608 }
609 Heed::vec velocity(dx, dy, dz);
610 velocity = velocity * Heed::CLHEP::c_light;
611
612 // Initial position (shift with respect to bounding box center and
613 // convert from cm to mm).
614 Heed::point p0((x0 - m_cX) * 10., (y0 - m_cY) * 10., (z0 - m_cZ) * 10.);
615
616 // Create and transport the photon.
617 Heed::HeedPhoton photon(m_chamber.get(), p0, velocity, t0, 0, e0 * 1.e-6,
618 m_fieldMap.get());
619 std::vector<Heed::gparticle*> secondaries;
620 photon.fly(secondaries);
621
622 while (!secondaries.empty()) {
623 std::vector<Heed::gparticle*> newSecondaries;
624 // Loop over the particle bank and look for daughter particles.
625 std::vector<Heed::gparticle*>::iterator it;
626 for (it = secondaries.begin(); it != secondaries.end(); ++it) {
627 // Check if it is a delta electron.
628 auto delta = dynamic_cast<Heed::HeedDeltaElectron*>(*it);
629 if (delta) {
630 if (m_doDeltaTransport) {
631 // Transport the delta electron.
632 delta->fly(newSecondaries);
633 // Add the conduction electrons to the list.
634 m_conductionElectrons.insert(m_conductionElectrons.end(),
635 delta->conduction_electrons.begin(),
636 delta->conduction_electrons.end());
637 m_conductionIons.insert(m_conductionIons.end(),
638 delta->conduction_ions.begin(),
639 delta->conduction_ions.end());
640 } else {
641 // Add the delta electron to the list, for later use.
642 deltaElectron newDeltaElectron;
643 newDeltaElectron.x = delta->position().x * 0.1 + m_cX;
644 newDeltaElectron.y = delta->position().y * 0.1 + m_cY;
645 newDeltaElectron.z = delta->position().z * 0.1 + m_cZ;
646 newDeltaElectron.t = delta->time();
647 newDeltaElectron.e = delta->kinetic_energy() * 1.e6;
648 newDeltaElectron.dx = delta->direction().x;
649 newDeltaElectron.dy = delta->direction().y;
650 newDeltaElectron.dz = delta->direction().z;
651 m_deltaElectrons.push_back(std::move(newDeltaElectron));
652 }
653 continue;
654 }
655 // Check if it is a fluorescence photon.
656 auto fluorescencePhoton = dynamic_cast<Heed::HeedPhoton*>(*it);
657 if (!fluorescencePhoton) {
658 std::cerr << m_className << "::TransportPhoton:\n"
659 << " Unknown secondary particle.\n";
660 ClearBank(secondaries);
661 ClearBank(newSecondaries);
662 return;
663 }
664 fluorescencePhoton->fly(newSecondaries);
665 }
666 secondaries.swap(newSecondaries);
667 ClearBank(newSecondaries);
668 }
669 ClearBank(secondaries);
670 // Get the total number of electrons produced in this step.
671 nel = m_doDeltaTransport ? m_conductionElectrons.size()
672 : m_deltaElectrons.size();
673 ni = m_conductionIons.size();
674}

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


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