moleculardna/src/TimeStepAction.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 #include "TimeStepAction.hh"
0028
0029 #include "ChromosomeMapper.hh"
0030 #include "DNAGeometry.hh"
0031 #include "DNAHit.hh"
0032 #include "DetectorConstruction.hh"
0033 #include "EventAction.hh"
0034 #include "OctreeNode.hh"
0035
0036 #include "G4DNAMolecularReactionTable.hh"
0037 #include "G4ITTrackHolder.hh"
0038 #include "G4ITTrackingManager.hh"
0039 #include "G4Molecule.hh"
0040 #include "G4RunManager.hh"
0041
0042 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0043 TimeStepAction::TimeStepAction(EventAction* event)
0044   : G4UserTimeStepAction(),
0045     fEventAction(event),
0046     fRadicalKillDistance(4.5 * nm),
0047     fpChemistryTrackHolder(G4ITTrackHolder::Instance())
0048 {
0049   AddTimeStep(1 * picosecond, 0.5 * nanosecond);
0050   // ctor
0051 }
0052
0053 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0054
0055 TimeStepAction::~TimeStepAction() = default;
0056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0057
0058 void TimeStepAction::StartProcessing()
0059 {
0060   auto det = dynamic_cast<const DetectorConstruction*>(
0061     G4RunManager::GetRunManager()->GetUserDetectorConstruction());
0062   fDNAGeometry = det->GetDNAGeometry();
0063   if (fDNAGeometry == nullptr) {
0064     G4ExceptionDescription exceptionDescription;
0065     exceptionDescription << "fDNAGeometry is null";
0066     G4Exception(
0067       "TimeStepAction"
0068       "StartProcessing",
0069       "no fDNAGeometry", FatalException, exceptionDescription);
0070   }
0071
0072   fRadicalKillDistance = fDNAGeometry->GetRadicalKillDistance();
0073 }
0074
0075 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0076
0077 void TimeStepAction::UserPostTimeStepAction()
0078 {
0079   RadicalKillDistance();
0080 }
0081
0082 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0083
0084 void TimeStepAction::UserPreTimeStepAction() {}
0085
0086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0087
0088 void TimeStepAction::UserReactionAction(const G4Track&, const G4Track&,
0089                                         const std::vector<G4Track*>*)
0090 {}
0091
0092 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0093
0094 void TimeStepAction::RadicalKillDistance()
0095 {
0096   if (fpChemistryTrackHolder == nullptr) {
0097     G4ExceptionDescription exceptionDescription;
0098     exceptionDescription << "fpChemistryTrackHolder is null";
0099     G4Exception(
0100       "TimeStepAction"
0101       "RadicalKillDistance",
0102       "NO_fpChemistryTrackHolder", FatalException, exceptionDescription);
0103   }
0104   G4Track* trackToKill;
0105   G4TrackManyList::iterator it_begin = fpChemistryTrackHolder->GetMainList()->begin();
0106   G4TrackManyList::iterator it_end = fpChemistryTrackHolder->GetMainList()->end();
0107   while (it_begin != it_end) {
0108     trackToKill = nullptr;
0109
0110     const G4VTouchable* touchable = it_begin->GetTouchable();
0111     if (touchable == nullptr) {
0112       ++it_begin;
0113       continue;
0114     }
0115
0116     G4LogicalVolume* trackLV = touchable->GetVolume()->GetLogicalVolume();
0117     G4LogicalVolume* motherLV = touchable->GetVolume()->GetMotherLogical();
0118
0119     OctreeNode* octree_track = fDNAGeometry->GetTopOctreeNode(trackLV);
0120     OctreeNode* octree_mother = fDNAGeometry->GetTopOctreeNode(motherLV);
0121
0122     if ((octree_track == nullptr) && (octree_mother == nullptr)) {
0123       trackToKill = *it_begin;
0124     }
0125     else if (octree_track != nullptr) {
0126       const G4AffineTransform& trans = it_begin->GetTouchable()->GetHistory()->GetTopTransform();
0127       G4ThreeVector pos = trans.TransformPoint(it_begin->GetPosition());
0128       size_t n = octree_track->SearchOctree(pos, fRadicalKillDistance).size();
0129       if (n == 0) {
0130         trackToKill = *it_begin;
0131       }
0132       if (fDNAGeometry->IsInsideHistone(trackLV, pos)) {
0133         trackToKill = *it_begin;
0134       }
0135     }
0136     else {
0137       const G4AffineTransform& trans = it_begin->GetTouchable()->GetHistory()->GetTopTransform();
0138       G4ThreeVector pos = trans.TransformPoint(it_begin->GetPosition());
0139       if (fDNAGeometry->IsInsideHistone(trackLV, pos)) {
0140         trackToKill = *it_begin;
0141       }
0142     }
0143     ++it_begin;
0144     if (trackToKill != nullptr) {
0145       fpChemistryTrackHolder->PushToKill(trackToKill);
0146     }
0147   }
0148 }
0149
0150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......