UHDR/src/PeriodicBoundaryProcess.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file PeriodicBoundaryProcess.cc
0027 /// \brief Implementation of the PeriodicBoundaryProcess class
0028
0029 /*
0030  * Based on 'g4pbc'.
0031  * Copyright (c) 2020 Amentum Pty Ltd
0032  * team@amentum.space
0033  * The original open-source version of this code
0034  * may be found at https://github.com/amentumspace/g4pbc
0035  * Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
0036  * associated documentation files (the "Software"), to deal in the Software without restriction,
0037  * including without limitation the rights to use, copy, modify, merge, publish, distribute,
0038  * sublicense, and/or sell copies of the Software, and to permit persons to whom the Software
0039  * is furnished to do so, subject to the following conditions:
0040  * The above copyright notice and this permission notice shall be included in all copies
0041  * or substantial portions of the Software.
0042  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
0043  * NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
0044  * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
0045  * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
0046  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
0047  *
0048  *
0049  */
0050 #include "PeriodicBoundaryProcess.hh"
0051
0052 #include "LogicalVolumePeriodic.hh"
0053 #include "ParticleChangeForPeriodic.hh"
0054
0055 #include "G4EventManager.hh"
0056 #include "G4GeometryTolerance.hh"
0057 #include "G4Navigator.hh"
0058 #include "G4PathFinder.hh"
0059 #include "G4TransportationManager.hh"
0060 #include "G4UnitsTable.hh"
0061 #include "G4VTrajectory.hh"
0062 #include "G4ios.hh"
0063
0064 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0065
0066 PeriodicBoundaryProcess::PeriodicBoundaryProcess(const G4String& processName, G4ProcessType type,
0067                                                  G4bool per_x, G4bool per_y, G4bool per_z)
0068   : G4VDiscreteProcess(processName, type), fPeriodicX(per_x), fPeriodicY(per_y), fPeriodicZ(per_z)
0069 {
0070   fkCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
0071   pParticleChange = &fParticleChange;
0072   G4PathFinder::GetInstance()->SetVerboseLevel(0);
0073 }
0074
0075 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0076
0077 G4VParticleChange* PeriodicBoundaryProcess::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep)
0078 {
0079   fTheStatus = Undefined;
0080
0081   fParticleChange.InitializeForPostStep(aTrack);
0082
0083   const G4Step* pStep = &aStep;
0084
0085   G4bool isOnBoundary = (pStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary);
0086
0087   if (!isOnBoundary) {
0088     fTheStatus = NotAtBoundary;
0089     if (verboseLevel > 0) {
0090       BoundaryProcessVerbose();
0091     }
0092     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
0093   }
0094
0095   if (aTrack.GetParentID() == 0)  // not primary particle (electron ?)
0096   {
0097     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
0098   }
0099
0100   auto thePrePV = pStep->GetPreStepPoint()->GetPhysicalVolume();
0101   auto thePostPV = pStep->GetPostStepPoint()->GetPhysicalVolume();
0102
0103   if (verboseLevel > 0) {
0104     G4cout << " Particle at Boundary! " << G4endl;
0105     if (thePrePV) {
0106       G4cout << " thePrePV:  " << thePrePV->GetName() << G4endl;
0107     }
0108     if (thePostPV) {
0109       G4cout << " thePostPV: " << thePostPV->GetName() << G4endl;
0110     }
0111     G4cout << "step length " << aTrack.GetStepLength() << G4endl;
0112     G4cout << "ParentID : " << aTrack.GetParentID() << "  TrackID : " << aTrack.GetTrackID()
0113            << " Position : " << aTrack.GetPosition()
0114            << "  aTrack Energy : " << G4BestUnit(aTrack.GetKineticEnergy(), "Energy") << G4endl;
0115   }
0116
0117   // avoid trapped particles at boundaries by testing for minimum step length
0118   if (aTrack.GetStepLength() <= fkCarTolerance / 2) {
0119     fTheStatus = StepTooSmall;
0120     if (verboseLevel > 0) {
0121       BoundaryProcessVerbose();
0122     }
0123     return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
0124   }
0125
0126   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
0127
0128   // store the current values
0129   fOldMomentum = aParticle->GetMomentumDirection();
0130   fOldPolarization = aParticle->GetPolarization();
0131   fOldPosition = pStep->GetPostStepPoint()->GetPosition();
0132   fNewPosition = fOldPosition;
0133
0134   if (verboseLevel > 0) {
0135     G4cout << " Old Momentum Direction: " << fOldMomentum << G4endl;
0136     G4cout << " Old Position: " << fNewPosition << G4endl;
0137     G4cout << "Get aTrack.GetParentID() : " << aTrack.GetParentID() << G4endl;
0138   }
0139
0140   auto theGlobalPoint = pStep->GetPostStepPoint()->GetPosition();
0141
0142   G4bool valid = false;
0143   //  Use the new method for Exit Normal in global coordinates,
0144   //    which provides the normal more reliably.
0145   // get from the real-world navigator (process not applied to parallel worlds)
0146   fTheGlobalNormal = G4TransportationManager::GetTransportationManager()
0147                        ->GetNavigatorForTracking()
0148                        ->GetGlobalExitNormal(theGlobalPoint, &valid);
0149
0150   if (valid) {
0151     fTheGlobalNormal = -fTheGlobalNormal;
0152   }
0153   else {
0154     G4cout << "global normal " << fTheGlobalNormal << G4endl;
0155     G4ExceptionDescription ed;
0156     ed << " PeriodicBoundaryProcess::PostStepDoIt(): "
0157        << " The Navigator reports that it returned an invalid normal" << G4endl;
0158     G4Exception("G4PeriodicBoundaryProcess::PostStepDoIt", "PerBoun01", FatalException, ed,
0159                 "Invalid Surface Normal - Geometry must return valid surface normal");
0160   }
0161
0162   G4bool isWrongDirection = fOldMomentum * fTheGlobalNormal > 0.0;
0163   if (isWrongDirection) {
0164     if (verboseLevel > 0) {
0165       G4cout << "ftheGlobalNormal points in a wrong direction." << G4endl;
0166       G4cout << "Invalid Surface Normal - Geometry must return valid surface \
0167         normal pointing in the right direction"
0168              << G4endl;
0169     }
0170
0171     fTheGlobalNormal = -fTheGlobalNormal;
0172   }
0173
0174   if (thePostPV == nullptr) {
0175     G4ExceptionDescription ed;
0176     ed << " PeriodicBoundaryProcess::PostStepDoIt(): "
0177        << " thePostPV == nullptr" << G4endl;
0178     G4Exception("G4PeriodicBoundaryProcess::PostStepDoIt", "PB14", FatalException, ed,
0179                 "Invalid thePostPV");
0180   }
0181   auto lvol = thePostPV->GetLogicalVolume();
0182
0183   if (lvol == nullptr) {
0184     G4ExceptionDescription ed;
0185     ed << " PeriodicBoundaryProcess::PostStepDoIt(): "
0186        << " lvol == nullptr" << G4endl;
0187     G4Exception("G4PeriodicBoundaryProcess::PostStepDoIt", "PB12", FatalException, ed,
0188                 "Invalid lvol");
0189   }
0190
0191   if (verboseLevel > 0) {
0192     G4cout << "Post step logical " << lvol->GetName() << G4endl;
0193   }
0194
0195   G4LogicalVolume* dlvol = nullptr;
0196
0197   if (lvol->GetNoDaughters() > 0) {
0198     if (verboseLevel > 0) {
0199       G4cout << "eldest daughter " << lvol->GetDaughter(0)->GetName() << G4endl;
0200     }
0201
0202     dlvol = lvol->GetDaughter(0)->GetLogicalVolume();
0203   }
0204
0205   if (dlvol && dlvol->IsExtended()) {
0206     if (verboseLevel > 0) {
0207       G4cout << " Logical surface, periodic " << G4endl;
0208     }
0209
0210     G4bool on_x = fTheGlobalNormal.isParallel(G4ThreeVector(1, 0, 0));
0211     G4bool on_y = fTheGlobalNormal.isParallel(G4ThreeVector(0, 1, 0));
0212     G4bool on_z = fTheGlobalNormal.isParallel(G4ThreeVector(0, 0, 1));
0213
0214     // make sure that we are at a plane
0215     G4bool on_plane = (on_x || on_y || on_z);
0216
0217     if (!on_plane) {
0218       G4ExceptionDescription ed;
0219       ed << " G4PeriodicBoundaryProcess::ostStepDoIt(): "
0220          << " The particle is not on a surface of the cyclic world" << G4endl;
0221       G4Exception(
0222         "G4PeriodicBoundaryProcess::PostStepDoIt", "Periodic01", FatalException, ed,
0223         "Periodic boundary process must only occur for particle on periodic world surface");
0224     }
0225     else {
0226       G4bool on_x_and_periodic = (on_x && fPeriodicX);
0227       G4bool on_y_and_periodic = (on_y && fPeriodicY);
0228       G4bool on_z_and_periodic = (on_z && fPeriodicZ);
0229
0230       G4bool on_a_periodic_plane = (on_x_and_periodic || on_y_and_periodic || on_z_and_periodic);
0231
0232       if (on_a_periodic_plane) {
0233         if (verboseLevel > 0) {
0234           G4cout << " on periodic plane " << G4endl;
0235         }
0236
0237         fTheStatus = Cycling;
0238
0239         if (verboseLevel > 0) {
0240           G4cout << " periodic " << G4endl;
0241         }
0242
0243         if (verboseLevel > 0) {
0244           G4cout << "Global normal " << fTheGlobalNormal << G4endl;
0245         }
0246
0247         // translate a component of the position vector according to which plane we are on
0248         if (on_x_and_periodic) {
0249           fNewPosition.setX(-fNewPosition.x());
0250         }
0251         else if (on_y_and_periodic) {
0252           fNewPosition.setY(-fNewPosition.y());
0253         }
0254         else if (on_z_and_periodic) {
0255           fNewPosition.setZ(-fNewPosition.z());
0256         }
0257         else {
0258           G4cout << "global normal does not belong to periodic plane!!" << G4endl;
0259         }
0260
0261         fNewMomentum = fOldMomentum.unit();
0262         fNewPolarization = fOldPolarization.unit();
0263
0264         if (verboseLevel > 0) {
0265           G4cout << " New Position: " << fNewPosition << G4endl;
0266           G4cout << " New Momentum Direction: " << fNewMomentum << G4endl;
0267           G4cout << " New Polarization:       " << fNewPolarization << G4endl;
0268           BoundaryProcessVerbose();
0269         }
0270
0271         fParticleChange.ProposeMomentumDirection(fNewMomentum);
0272         fParticleChange.ProposePolarization(fNewPolarization);
0273         fParticleChange.ProposePosition(fNewPosition);
0274
0275         G4PathFinder::GetInstance()->ReLocate(fNewPosition);
0276         G4PathFinder::GetInstance()->ComputeSafety(fNewPosition);
0277
0278         // we must notify the navigator that we have moved the particle artificially
0279         G4Navigator* gNavigator =
0280           G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
0281
0282         // Inform the navigator that the previous Step calculated
0283         // by the geometry was taken in its entirety.
0284         gNavigator->SetGeometricallyLimitedStep();
0285
0286         // Search the geometrical hierarchy for the volumes deepest in the hierarchy
0287         // containing the point in the global coordinate space.
0288         gNavigator->LocateGlobalPointAndSetup(fNewPosition, &fNewMomentum, false,
0289                                               false);  // do not ignore direction
0290
0291         // Calculate the isotropic distance to the nearest boundary from the
0292         // specified point in the global coordinate system.
0293         gNavigator->ComputeSafety(fNewPosition);
0294
0295         // Locates the volume containing the specified global point.
0296         // force drawing of the step prior to periodic the particle
0297         auto evtm = G4EventManager::GetEventManager();
0298         auto tckm = evtm->GetTrackingManager();
0299         auto pTrajectory = tckm->GimmeTrajectory();
0300         if (pTrajectory) {
0301           pTrajectory->AppendStep(pStep);
0302         }
0303       }
0304     }
0305   }
0306   return &fParticleChange;
0307 }
0308
0309 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0310
0311 G4double PeriodicBoundaryProcess::GetMeanFreePath(const G4Track&, G4double,
0312                                                   G4ForceCondition* condition)
0313 {
0314   *condition = Forced;
0315   return DBL_MAX;
0316 }
0317
0318 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0319
0320 void PeriodicBoundaryProcess::BoundaryProcessVerbose()
0321 {
0322   if (fStatusMessages.count(fTheStatus) > 0) {
0323     G4cout << PeriodicBoundaryProcess::fStatusMessages[fTheStatus] << G4endl;
0324   }
0325 }
0326
0327 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......