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