moleculardna/src/SteppingAction.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 /// file:SteppingAction.cc
0028 /// brief:
0029 #include "SteppingAction.hh"
0030
0031 #include "ChromosomeMapper.hh"
0032 #include "DNAGeometry.hh"
0033 #include "DNAHit.hh"
0034 #include "DetectorConstruction.hh"
0035 #include "EventAction.hh"
0036 #include "OctreeNode.hh"
0037 #include "PlacementVolumeInfo.hh"
0038 #include "UtilityFunctions.hh"
0039
0040 #include "G4AffineTransform.hh"
0041 #include "G4LogicalVolume.hh"
0042 #include "G4RunManager.hh"
0043 #include "G4StepPoint.hh"
0044 #include "G4ThreeVector.hh"
0045 #include "G4VPhysicalVolume.hh"
0046 #include "G4VTouchable.hh"
0047 #include "Randomize.hh"
0048
0049 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0050
0051 SteppingAction::SteppingAction(EventAction* evt) : fEventAction(evt)
0052 {
0053   const auto* det = dynamic_cast<const DetectorConstruction*>(
0054     G4RunManager::GetRunManager()->GetUserDetectorConstruction());
0055
0056   fDNAGeometry = det->GetDNAGeometry();
0057 }
0058
0059 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0060
0061 void SteppingAction::UserSteppingAction(const G4Step* aStep)
0062 {
0063   G4double edep = aStep->GetTotalEnergyDeposit();
0064   if (edep == 0) {
0065     return;
0066   }
0067
0068   // Work out if we are in the Logical volume assigned to DNA
0069   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0070
0071   // remember the chem volume is the "real" volume.
0072   G4LogicalVolume* dnaVol = fDNAGeometry->GetDNAChemVolumePointer();
0073   G4int depth = touch->GetHistoryDepth();
0074   G4bool isInDNARegion = false;
0075
0076   if (depth >= 0) {
0077     fEventAction->AddCellEdep(edep);  // dousatsu
0078   }
0079
0080   while ((depth >= 0) && !isInDNARegion) {
0081     if (touch->GetVolume(depth) != nullptr) {
0082       isInDNARegion = isInDNARegion || (touch->GetVolume(depth)->GetLogicalVolume() == dnaVol);
0083     }
0084     else {
0085       G4cout << "Null pointer volume in mother with depth of " << touch->GetHistoryDepth()
0086              << " looking for volume at depth " << depth << G4endl;
0087       assert(depth < 0);
0088     }
0089     depth--;
0090   }
0091
0092   if (isInDNARegion) {
0093     ProcessDNARegionHit(aStep);
0094   }
0095
0096   G4ThreeVector prepos = aStep->GetPreStepPoint()->GetPosition();
0097   if (fDNAGeometry->GetChromosomeMapper()->GetChromosome(prepos) != nullptr) {
0098     DoChromosomeRegionHit(prepos, edep, false);
0099   }
0100
0101   if (aStep->GetTrack()->GetTrackID() == 1) {
0102     fEventAction->SetPrimStopPos(prepos);
0103     G4double tl = aStep->GetStepLength();
0104     fEventAction->AddTrackLengthCell(tl);
0105     if (fDNAGeometry->GetChromosomeMapper()->GetChromosome(prepos) != nullptr) {
0106       fEventAction->AddTrackLengthChro(tl);
0107     }
0108   }
0109 }
0110
0111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0112
0113 void SteppingAction::ProcessDNARegionHit(const G4Step* aStep)
0114 {
0115   G4double edep = aStep->GetTotalEnergyDeposit();
0116   const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
0117   G4ThreeVector pos =
0118     aStep->GetPreStepPoint()->GetPosition()
0119     + G4UniformRand()
0120         * (aStep->GetPostStepPoint()->GetPosition() - aStep->GetPreStepPoint()->GetPosition());
0121   G4VPhysicalVolume* pv = aStep->GetPreStepPoint()->GetPhysicalVolume();
0122   G4LogicalVolume* lv = pv->GetLogicalVolume();
0123
0124   std::vector<G4VPhysicalVolume*> molvec;
0125   OctreeNode* parentNode;
0126   G4AffineTransform transform = touch->GetHistory()->GetTopTransform();
0127   parentNode = fDNAGeometry->GetTopOctreeNode(lv);
0128   if (parentNode == nullptr) {  // See if parent has an octree.
0129     parentNode = fDNAGeometry->GetTopOctreeNode(pv->GetMotherLogical());
0130     if (parentNode != nullptr) {  // If mother is a node, update it to be the PV also
0131       // These methods are unsafe if the pv is the top volume so they
0132       // need to be in this if statement.
0133       // lv = pv->GetMotherLogical();
0134       pv = touch->GetVolume(1);
0135       transform = touch->GetHistory()->GetTransform(touch->GetHistoryDepth() - 1);
0136     }
0137   }
0138
0139   G4ThreeVector localPos = transform.TransformPoint(pos);
0140
0141   if (parentNode != nullptr) {  // Get Molecule Vector
0142     molvec = parentNode->SearchOctree(localPos, fDNAGeometry->GetDirectInteractionRange());
0143   }
0144
0145   // debug output for finding logical volumes.
0146   if (fDNAGeometry->GetVerbosity() > 1) {
0147     if (parentNode != nullptr) {
0148       G4String lvname = pv->GetLogicalVolume()->GetName();
0149       G4cout << "Found octree for logical volume: " << lvname << G4endl;
0150     }
0151     else {
0152       // Don't bother doing error printouts if the particle
0153       // is at the DNA mother volume level.
0154       // Remember the Chemistry Volume is the "real" world
0155       if (pv->GetLogicalVolume() != fDNAGeometry->GetDNAChemVolumePointer()) {
0156         G4String lvname = pv->GetLogicalVolume()->GetName();
0157         G4String motherlvname;
0158         if (lv != nullptr) {
0159           motherlvname = lv->GetName();
0160         }
0161         G4cout << "Could not find octree for logical volume." << G4endl
0162                << "Particle in LV: " << lvname << ", with mother LV: " << motherlvname << G4endl;
0163       }
0164     }
0165   }
0166
0167   if (!molvec.empty()) {
0168     // NOTE: This loop can be optimised to avoid checking pairs twice
0169     G4VPhysicalVolume* closest = molvec[0];
0170     G4double dist = (closest->GetTranslation() - localPos).mag();
0171     if (molvec.size() > 1) {
0172       G4double newdist;
0173       for (auto it = molvec.begin() + 1; it != molvec.end(); it++) {
0174         newdist = ((*it)->GetTranslation() - localPos).mag();
0175         if (newdist < dist) {
0176           dist = newdist;
0177           closest = (*it);
0178         }
0179       }
0180     }
0181     G4bool isDNAHit = (dist < fDNAGeometry->GetDirectInteractionRange());
0182     if (isDNAHit) {
0183       DoChromosomeDNAHit(pos, localPos, edep, dist, closest->GetName(), pv->GetName());
0184       DoChromosomeRegionHit(pos, edep, true);
0185     }
0186   }
0187 }
0188
0189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0190
0191 void SteppingAction::DoChromosomeRegionHit(const G4ThreeVector& pos, const G4double& edep,
0192                                            const G4bool& dnahit)
0193 {
0194   // TODO: chromosome should probably be based on position of center of
0195   //       parent volume, not of current pos
0196   G4String chromosome = fDNAGeometry->GetChromosomeMapper()->GetCurrentChromosomeKey(pos);
0197   if (!(chromosome.empty())) {
0198     fEventAction->AddChromosomeEdep(chromosome, edep, dnahit);
0199   }
0200 }
0201
0202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0203
0204 void SteppingAction::DoChromosomeDNAHit(const G4ThreeVector& pos, const G4ThreeVector& localPos,
0205                                         const G4double& edep, const G4double& dist,
0206                                         const G4String& pvName, const G4String& motherName)
0207 {
0208   G4int placementIdx = std::stoi(utility::Split(motherName, '-').at(1));
0209   std::vector<G4String> pvNameVec = utility::Split(pvName, '-');
0210   molecule mol = utility::GetMoleculeEnum(pvNameVec.at(0));
0211   G4int chainIdx = std::stoi(pvNameVec.at(1));
0212   G4int strandIdx = std::stoi(pvNameVec.at(2));
0213   int64_t baseIdx = std::stoll(pvNameVec.at(3));
0214
0215   // Convert to local chain and base index
0216   chainIdx = fDNAGeometry->GetGlobalChain(placementIdx, chainIdx);
0217   baseIdx += fDNAGeometry->GetStartIdx(placementIdx, chainIdx);
0218
0219   if (fDNAGeometry->GetStrandsFlipped(placementIdx)) {
0220     strandIdx = (strandIdx + 1) % 2;
0221   }
0222   G4String chromosome = fDNAGeometry->GetChromosomeMapper()->GetCurrentChromosomeKey(pos);
0223   const DNAHit* dnaHit = new DNAHit(mol, placementIdx, chainIdx, strandIdx, baseIdx, pos, localPos,
0224                                     edep, dist, chromosome);
0225
0226   fEventAction->AddDNAHit(dnaHit);
0227 }
0228
0229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......