3838
3939#include < fairlogger/Logger.h>
4040
41- #include " TG4RunConfiguration.h"
4241#include " TG4ComposedPhysicsList.h"
4342#include " SimSetup/G4RunConfiguration.h"
4443
4847#include < G4PhysicsListHelper.hh>
4948#include < G4ProcessManager.hh>
5049#include < G4ProcessVector.hh>
51- #include < G4VEnergyLossProcess.hh>
5250#include < G4mplIonisation.hh>
5351
5452#include < G4ChargeState.hh>
5553#include < G4ChordFinder.hh>
5654#include < G4ClassicalRK4.hh>
5755#include < G4EquationOfMotion.hh>
58- #include < G4EventManager.hh>
5956#include < G4FieldManager.hh>
60- #include < G4MagIntegratorStepper.hh>
6157#include < G4MagneticField.hh>
58+ #include < G4PropagatorInField.hh>
6259#include < G4Track.hh>
63- #include < G4TrackingManager.hh>
6460#include < G4CoupledTransportation.hh>
6561#include < G4Transportation.hh>
6662#include < G4TransportationManager.hh>
7066
7167#include < array>
7268#include < cmath>
73- #include < vector>
7469
7570namespace o2
7671{
@@ -86,8 +81,6 @@ namespace
8681constexpr std::array<int , 4 > gMonopolePDGs = {o2::sim::MonopolePdgSymm, -o2::sim::MonopolePdgSymm,
8782 o2::sim::MonopolePdgAsymm, -o2::sim::MonopolePdgAsymm};
8883
89- inline bool isMonopolePDG (int pdg) { return o2::sim::isMonopole (pdg); }
90-
9184// / GEANT4 only queries the magnetic field when a particle has non-zero electric
9285// / charge or non-zero magnetic moment (μ) and the monopole has neither. So this is a workaround
9386// / so that G4Transportation queries the field for it
@@ -157,27 +150,16 @@ class O2MonopoleEquation : public G4EquationOfMotion
157150 {
158151 }
159152
153+ // / Sign of the magnetic charge of the monopole about to be transported (+1 or -1),
154+ // / set by O2MonopoleTransportation before each step
155+ void setMagneticChargeSign (int sign) { mMagneticChargeSign = sign; }
156+
160157 void SetChargeMomentumMass (G4ChargeState particleChargeState,
161158 G4double /* momentum*/ , G4double particleMass) override
162159 {
163160 mElCharge = CLHEP ::eplus * particleChargeState.GetCharge () * CLHEP ::c_light;
164161 mMassCof = particleMass * particleMass;
165-
166- // Non-zero only while an O2 monopole is being transported. The sign follows
167- // the PDG sign, so monopole and anti-monopole are pushed in opposite
168- // directions along B.
169- // TODO: both species are electrically neutral today, so their magnetic sign is
170- // the same function of the PDG sign; once they carry electric charge the
171- // asymmetric species (4120000) needs the opposite relative sign.
172- double signedMagneticCharge = 0 .;
173- int pdg = 0 ;
174- if (const G4Track* track = currentTrack ()) {
175- pdg = track->GetDefinition ()->GetPDGEncoding ();
176- if (isMonopolePDG (pdg)) {
177- signedMagneticCharge = (pdg > 0 ) ? mMagneticChargeEplus : -mMagneticChargeEplus ;
178- }
179- }
180- mMagCharge = CLHEP ::eplus * signedMagneticCharge * CLHEP ::c_light;
162+ mMagCharge = CLHEP ::eplus * mMagneticChargeSign * mMagneticChargeEplus * CLHEP ::c_light;
181163 }
182164
183165 void EvaluateRhsGivenB (const G4double y[], const G4double B[3 ], G4double dydx[]) const override
@@ -205,8 +187,7 @@ class O2MonopoleEquation : public G4EquationOfMotion
205187
206188 // Coupling of the magnetic charge to an electric field: the full dual
207189 // Lorentz force is F = g*(B - v x E/c^2).
208- // mMagCharge is zero for every non-monopole, so nothing else is affected.
209- if (mMagCharge != 0 . && mTPCDriftField > 0 .) {
190+ if (mTPCDriftField > 0 .) {
210191 G4double E[3 ] = {0 ., 0 ., 0 .};
211192 if (tpcDriftField (y, mTPCDriftField , E)) {
212193 // Relative to cofMag this carries 1/(c*energy), so the term is of order
@@ -227,38 +208,32 @@ class O2MonopoleEquation : public G4EquationOfMotion
227208 }
228209
229210 private:
230- static const G4Track* currentTrack ()
231- {
232- auto * eventManager = G4EventManager::GetEventManager ();
233- if (eventManager == nullptr ) {
234- return nullptr ;
235- }
236- auto * trackingManager = eventManager->GetTrackingManager ();
237- return trackingManager != nullptr ? trackingManager->GetTrack () : nullptr ;
238- }
239-
240211 double mMagneticChargeEplus ; // /< |g| in eplus units (1 g_D = 1/(2*alpha) ~ 68.5)
212+ int mMagneticChargeSign = 1 ; // /< sign of the magnetic charge of the current monopole
241213 G4double mTPCDriftField = 0 .; // /< TPC drift field, Geant4 units; 0 disables the coupling
242214 G4double mElCharge = 0 .;
243215 G4double mMagCharge = 0 .;
244216 G4double mMassCof = 0 .;
245217};
246218
247219// ____________________________________________________________________________
248- // / The field manager of the run, together with a second chord finder that
249- // / integrates O2MonopoleEquation
220+ // / The chord finder integrating O2MonopoleEquation, and that equation
250221struct MonopoleFieldSetup {
251- G4FieldManager* fieldManager = nullptr ;
252222 G4ChordFinder* chordFinder = nullptr ;
223+ O2MonopoleEquation* equation = nullptr ; // /< to pass the sign of the magnetic charge at each step
253224};
254225
255226// ____________________________________________________________________________
256227// / Build a chord finder that integrates O2MonopoleEquation instead of the
257- // / default electric-charge-only equation.
228+ // / default electric-charge-only equation, in the field of the run .
258229// /
259- // / It is deliberately NOT installed on the field manager here. The default
260- // / chord finder is what Geant4-VMC configured for this run (NystromRK4 unless
261- // / the macro says otherwise) and every ordinary particle keeps it
230+ // / It is deliberately NOT installed on a field manager here. The default
231+ // / chord finder is what Geant4-VMC configured for the field (NystromRK4 unless
232+ // / the macro says otherwise) and every ordinary particle keeps it.
233+ // / O2MonopoleTransportation installs it only for the duration of a monopole step, on the
234+ // / field manager of the current volume: the global one, or one of the local ones created for
235+ // / the volumes with their own field parameters (G4LocalFieldConstruction). The local field
236+ // / managers carry the same field as the global one, so one chord finder serves them all.
262237// /
263238// / SetUserEquationOfMotion() in Geant4-VMC is not usable here: it
264239// / registers the object with TG4GeometryManager, and the field integrator is
@@ -308,58 +283,90 @@ inline MonopoleFieldSetup buildMonopoleFieldSetup(double magneticChargeEplusUnit
308283 LOG (info) << " O2MonopolePhysics: monopole equation of motion built (F = g*B), magnetic charge = "
309284 << magneticChargeEplusUnits << " eplus (TPC drift field coupling disabled)" ;
310285 }
311- return {fieldManager, chordFinder };
286+ return {chordFinder, equation };
312287}
313288
314289// ____________________________________________________________________________
315290// / Transportation for the monopole species only.
316291// /
317292// / These must be true for the monopoles and false for every other particle:
318293// /
319- // / - the field manager has to use the monopole chord finder, otherwise O2MonopoleEquation
320- // / might not be evaluated at all (such as it happens with the default NystromRK4)
294+ // / - the field manager of the current volume (the global one, or a local one) has to use the
295+ // / monopole chord finder, otherwise O2MonopoleEquation might not be evaluated at all
296+ // / (such as it happens with the default NystromRK4)
321297// / - G4Transportation has to consider the magnetic moment
298+ // / - the field manager has to declare that the field changes the energy: the g*B force does
299+ // / work on the monopole. Otherwise G4Transportation would reset the kinetic energy at the
300+ // / end of every step to the one at its start, keeping only the change of direction
322301// /
323302// / Ordinary tracks keep the stepper the run was configured with, and neutral particles that
324303// / happen to carry a magnetic moment (neutrons above all) keep their
325304// / straight-line transport.
326305// /
327- // / This mirrors G4MonopoleTransportation from the Geant4 monopole example, but
328- // / by derives directly from G4Transportation, so that it follows the Geant4 versions
306+ // / This mirrors G4MonopoleTransportation from the Geant4 monopole example, which takes the
307+ // / kinetic energy and the time of flight from the integrated track, but derives directly
308+ // / from G4Transportation, so that it follows the Geant4 versions
329309class O2MonopoleTransportation : public G4Transportation
330310{
331311 public:
332- O2MonopoleTransportation (G4FieldManager* fieldManager, G4ChordFinder* monopoleChordFinder )
333- : G4Transportation(0 ), mFieldManager (fieldManager ), mMonopoleChordFinder (monopoleChordFinder )
312+ O2MonopoleTransportation (const MonopoleFieldSetup& fieldSetup )
313+ : G4Transportation(0 ), mChordFinder (fieldSetup.chordFinder ), mEquation (fieldSetup.equation )
334314 {
335315 }
336316
317+ void StartTracking (G4Track* track) override
318+ {
319+ G4Transportation::StartTracking (track);
320+ // G4Transportation::StartTracking() clears the integration state only of the chord finders
321+ // installed on the field managers. Without this the step estimate of the previous monopole
322+ // would carry over, and the result would depend on which events the worker simulated before
323+ mChordFinder ->ResetStepEstimate ();
324+ }
325+
337326 G4double AlongStepGetPhysicalInteractionLength (const G4Track& track, G4double previousStepSize,
338327 G4double currentMinimumStep, G4double& currentSafety,
339328 G4GPILSelection* selection) override
340329 {
330+ // the field manager the propagator uses in this volume, as G4Transportation finds it
331+ auto * fieldManager = GetPropagatorInField ()->FindAndSetFieldManager (track.GetVolume ());
332+ if (fieldManager->GetDetectorField () == nullptr ) {
333+ // volume without field: straight-line transport
334+ return G4Transportation::AlongStepGetPhysicalInteractionLength (track, previousStepSize, currentMinimumStep,
335+ currentSafety, selection);
336+ }
337+ // The sign of the magnetic charge follows the PDG sign, so monopole and anti-monopole
338+ // are pushed in opposite directions along B.
339+ // TODO: both species are electrically neutral currently, so their magnetic sign is
340+ // the same function of the PDG sign; once they carry electric charge the
341+ // asymmetric species (4120000) needs the opposite relative sign.
342+ mEquation ->setMagneticChargeSign (track.GetDefinition ()->GetPDGEncoding () > 0 ? 1 : -1 );
343+
341344 const G4bool previousMoment = G4Transportation::EnableMagneticMoment (true );
342- auto * previousChordFinder = mFieldManager ->GetChordFinder ();
343- mFieldManager ->SetChordFinder (mMonopoleChordFinder );
345+ auto * previousChordFinder = fieldManager->GetChordFinder ();
346+ const G4bool previousChangesEnergy = fieldManager->DoesFieldChangeEnergy ();
347+ fieldManager->SetChordFinder (mChordFinder );
348+ // the end kinetic energy and time of flight are then taken from the integrated track
349+ fieldManager->SetFieldChangesEnergy (true );
344350
345351 const G4double length = G4Transportation::AlongStepGetPhysicalInteractionLength (
346352 track, previousStepSize, currentMinimumStep, currentSafety, selection);
347353
348- mFieldManager ->SetChordFinder (previousChordFinder);
354+ fieldManager->SetFieldChangesEnergy (previousChangesEnergy);
355+ fieldManager->SetChordFinder (previousChordFinder);
349356 G4Transportation::EnableMagneticMoment (previousMoment);
350357 return length;
351358 }
352359
353360 private:
354- G4FieldManager* mFieldManager ; // /< field manager of the run , not owned
355- G4ChordFinder* mMonopoleChordFinder ; // /< installed only for the duration of a monopole step , not owned
361+ G4ChordFinder* mChordFinder ; // /< installed only for the duration of a monopole step , not owned
362+ O2MonopoleEquation* mEquation ; // /< the equation integrated by mChordFinder , not owned
356363};
357364
358365// ____________________________________________________________________________
359366// / Replace the transport of one monopole species with
360367// / O2MonopoleTransportation, keeping it first in the DoIt vectors exactly as
361368// / G4VUserPhysicsList::AddTransportation() left it.
362- inline bool installMonopoleTransport (G4ProcessManager* pmanager, const G4ParticleDefinition* particle,
369+ inline void installMonopoleTransport (G4ProcessManager* pmanager, const G4ParticleDefinition* particle,
363370 const MonopoleFieldSetup& fieldSetup)
364371{
365372 G4VProcess* existing = nullptr ;
@@ -374,7 +381,7 @@ inline bool installMonopoleTransport(G4ProcessManager* pmanager, const G4Particl
374381 LOG (fatal) << " O2MonopolePhysics: " << particle->GetParticleName ()
375382 << " has no transportation process to replace; the monopole could not be "
376383 " coupled to the field" ;
377- return false ;
384+ return ;
378385 }
379386 if (dynamic_cast <G4CoupledTransportation*>(existing) != nullptr ) {
380387 // Parallel worlds are in use. O2MonopoleTransportation derives from plain
@@ -383,15 +390,15 @@ inline bool installMonopoleTransport(G4ProcessManager* pmanager, const G4Particl
383390 LOG (fatal) << " O2MonopolePhysics: " << particle->GetParticleName ()
384391 << " uses G4CoupledTransportation (parallel worlds); the monopole transportation "
385392 " does not support that. Rerun with G4.monopole=0 or without parallel worlds" ;
386- return false ;
393+ return ;
387394 }
388395
389396 // One G4Transportation instance is shared by every particle
390397 // (G4VUserPhysicsList::AddTransportation creates a single one), so it is
391398 // detached from this particle only and must not be deleted.
392399 pmanager->RemoveProcess (existing);
393400
394- auto * transportation = new O2MonopoleTransportation (fieldSetup. fieldManager , fieldSetup. chordFinder );
401+ auto * transportation = new O2MonopoleTransportation (fieldSetup);
395402 pmanager->AddProcess (transportation);
396403 // Transportation has to be first in the DoIt vectors
397404 //
@@ -402,17 +409,17 @@ inline bool installMonopoleTransport(G4ProcessManager* pmanager, const G4Particl
402409 // the warning in the stdout is expected and harmless.
403410 pmanager->SetProcessOrderingToFirst (transportation, idxAlongStep);
404411 pmanager->SetProcessOrderingToFirst (transportation, idxPostStep);
405- return true ;
406412}
407413
408414// ____________________________________________________________________________
409- // / Minimal physics list whose only job is to attach the magnetic-monopole
410- // / ionisation process to the already-defined O2 monopole particles.
415+ // / Minimal physics list attaching the magnetic-monopole ionisation and transport
416+ // / to the already-defined O2 monopole particles.
411417// /
412418// / It is implemented as a G4VUserPhysicsList so that it can be handed to VMC's
413419// / TG4ComposedPhysicsList::AddPhysicsList(). It intentionally does NOT create
414- // / any particle and does NOT add transportation (both are handled by the
415- // / reference physics list), it only adds one extra process.
420+ // / any particle (they are defined by O2MCApplication::AddParticles()). For each
421+ // / monopole it adds mplIoni and replaces the transportation of the reference
422+ // / physics list with O2MonopoleTransportation.
416423class O2MonopolePhysics : public G4VUserPhysicsList
417424{
418425 public:
@@ -431,7 +438,7 @@ class O2MonopolePhysics : public G4VUserPhysicsList
431438 auto * table = G4ParticleTable::GetParticleTable ();
432439
433440 // Built on the first monopole species actually found, then shared by the rest:
434- // O2MonopoleEquation reads the sign of the magnetic charge off the track , so
441+ // O2MonopoleTransportation sets the sign of the magnetic charge at each step , so
435442 // one chord finder serves monopoles and anti-monopoles
436443 // No monopoles in the particles table == no monopole ionisation attached
437444 MonopoleFieldSetup fieldSetup;
@@ -451,28 +458,9 @@ class O2MonopolePhysics : public G4VUserPhysicsList
451458 LOG (warning) << " O2MonopolePhysics: no process manager for PDG " << pdg << " , skipping" ;
452459 continue ;
453460 }
454- // The reference EM physics list attaches the standard hadron ionisation
455- // (hIoni) to the monopole, which is wrong and makes G4LossTableManager crash
456- // while merging the two dE/dx tables.
457- // Any pre-existing energy-loss process is removed so that mplIoni is the
458- // monopole single, correct ionisation process. This is taken from Geant4
459- // exoticphysics/monopole example
460- {
461- std::vector<G4VProcess*> toRemove;
462- G4ProcessVector* plist = pmanager->GetProcessList ();
463- for (G4int ip = 0 ; ip < static_cast <G4int>(plist->size ()); ++ip) {
464- G4VProcess* proc = (*plist)[ip];
465- if (dynamic_cast <G4VEnergyLossProcess*>(proc) != nullptr ) {
466- toRemove.push_back (proc);
467- }
468- }
469- for (auto * proc : toRemove) {
470- LOG (info) << " O2MonopolePhysics: removing pre-existing energy-loss process "
471- << proc->GetProcessName () << " from " << particle->GetParticleName ();
472- pmanager->RemoveProcess (proc);
473- }
474- }
475- // Ordered by G4PhysicsListHelper from the process subtype, as in the Geant4 monopole example
461+ // The only ionisation of the monopoles: they are defined as kPTUndefined, for which
462+ // Geant4-VMC attaches no process. Ordered by G4PhysicsListHelper from the process
463+ // subtype, as in the Geant4 monopole example
476464 G4PhysicsListHelper::GetPhysicsListHelper ()->RegisterProcess (new G4mplIonisation (mMagneticCharge ), particle);
477465
478466 // G4Transportation only looks up the field when the particle has a
0 commit comments