| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 /// file:SteppingAction.cc
28 /// brief:
29 #include "SteppingAction.hh"
30
31 #include "ChromosomeMapper.hh"
32 #include "DNAGeometry.hh"
33 #include "DNAHit.hh"
34 #include "DetectorConstruction.hh"
35 #include "EventAction.hh"
36 #include "OctreeNode.hh"
37 #include "PlacementVolumeInfo.hh"
38 #include "UtilityFunctions.hh"
39
40 #include "G4AffineTransform.hh"
41 #include "G4LogicalVolume.hh"
42 #include "G4RunManager.hh"
43 #include "G4StepPoint.hh"
44 #include "G4ThreeVector.hh"
45 #include "G4VPhysicalVolume.hh"
46 #include "G4VTouchable.hh"
47 #include "Randomize.hh"
48
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50
51 SteppingAction::SteppingAction(EventAction* evt) : fEventAction(evt)
52 {
53 const auto* det = dynamic_cast<const DetectorConstruction*>(
54 G4RunManager::GetRunManager()->GetUserDetectorConstruction());
55
56 fDNAGeometry = det->GetDNAGeometry();
57 }
58
59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60
61 void SteppingAction::UserSteppingAction(const G4Step* aStep)
62 {
63 G4double edep = aStep->GetTotalEnergyDeposit();
64 if (edep == 0) {
65 return;
66 }
67
68 // Work out if we are in the Logical volume assigned to DNA
69 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
70
71 // remember the chem volume is the "real" volume.
72 G4LogicalVolume* dnaVol = fDNAGeometry->GetDNAChemVolumePointer();
73 G4int depth = touch->GetHistoryDepth();
74 G4bool isInDNARegion = false;
75
76 if (depth >= 0) {
77 fEventAction->AddCellEdep(edep); // dousatsu
78 }
79
80 while ((depth >= 0) && !isInDNARegion) {
81 if (touch->GetVolume(depth) != nullptr) {
82 isInDNARegion = isInDNARegion || (touch->GetVolume(depth)->GetLogicalVolume() == dnaVol);
83 }
84 else {
85 G4cout << "Null pointer volume in mother with depth of " << touch->GetHistoryDepth()
86 << " looking for volume at depth " << depth << G4endl;
87 assert(depth < 0);
88 }
89 depth--;
90 }
91
92 if (isInDNARegion) {
93 ProcessDNARegionHit(aStep);
94 }
95
96 G4ThreeVector prepos = aStep->GetPreStepPoint()->GetPosition();
97 if (fDNAGeometry->GetChromosomeMapper()->GetChromosome(prepos) != nullptr) {
98 DoChromosomeRegionHit(prepos, edep, false);
99 }
100
101 if (aStep->GetTrack()->GetTrackID() == 1) {
102 fEventAction->SetPrimStopPos(prepos);
103 G4double tl = aStep->GetStepLength();
104 fEventAction->AddTrackLengthCell(tl);
105 if (fDNAGeometry->GetChromosomeMapper()->GetChromosome(prepos) != nullptr) {
106 fEventAction->AddTrackLengthChro(tl);
107 }
108 }
109 }
110
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
112
113 void SteppingAction::ProcessDNARegionHit(const G4Step* aStep)
114 {
115 G4double edep = aStep->GetTotalEnergyDeposit();
116 const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
117 G4ThreeVector pos =
118 aStep->GetPreStepPoint()->GetPosition()
119 + G4UniformRand()
120 * (aStep->GetPostStepPoint()->GetPosition() - aStep->GetPreStepPoint()->GetPosition());
121 G4VPhysicalVolume* pv = aStep->GetPreStepPoint()->GetPhysicalVolume();
122 G4LogicalVolume* lv = pv->GetLogicalVolume();
123
124 std::vector<G4VPhysicalVolume*> molvec;
125 OctreeNode* parentNode;
126 G4AffineTransform transform = touch->GetHistory()->GetTopTransform();
127 parentNode = fDNAGeometry->GetTopOctreeNode(lv);
128 if (parentNode == nullptr) { // See if parent has an octree.
129 parentNode = fDNAGeometry->GetTopOctreeNode(pv->GetMotherLogical());
130 if (parentNode != nullptr) { // If mother is a node, update it to be the PV also
131 // These methods are unsafe if the pv is the top volume so they
132 // need to be in this if statement.
133 // lv = pv->GetMotherLogical();
134 pv = touch->GetVolume(1);
135 transform = touch->GetHistory()->GetTransform(touch->GetHistoryDepth() - 1);
136 }
137 }
138
139 G4ThreeVector localPos = transform.TransformPoint(pos);
140
141 if (parentNode != nullptr) { // Get Molecule Vector
142 molvec = parentNode->SearchOctree(localPos, fDNAGeometry->GetDirectInteractionRange());
143 }
144
145 // debug output for finding logical volumes.
146 if (fDNAGeometry->GetVerbosity() > 1) {
147 if (parentNode != nullptr) {
148 G4String lvname = pv->GetLogicalVolume()->GetName();
149 G4cout << "Found octree for logical volume: " << lvname << G4endl;
150 }
151 else {
152 // Don't bother doing error printouts if the particle
153 // is at the DNA mother volume level.
154 // Remember the Chemistry Volume is the "real" world
155 if (pv->GetLogicalVolume() != fDNAGeometry->GetDNAChemVolumePointer()) {
156 G4String lvname = pv->GetLogicalVolume()->GetName();
157 G4String motherlvname;
158 if (lv != nullptr) {
159 motherlvname = lv->GetName();
160 }
161 G4cout << "Could not find octree for logical volume." << G4endl
162 << "Particle in LV: " << lvname << ", with mother LV: " << motherlvname << G4endl;
163 }
164 }
165 }
166
167 if (!molvec.empty()) {
168 // NOTE: This loop can be optimised to avoid checking pairs twice
169 G4VPhysicalVolume* closest = molvec[0];
170 G4double dist = (closest->GetTranslation() - localPos).mag();
171 if (molvec.size() > 1) {
172 G4double newdist;
173 for (auto it = molvec.begin() + 1; it != molvec.end(); it++) {
174 newdist = ((*it)->GetTranslation() - localPos).mag();
175 if (newdist < dist) {
176 dist = newdist;
177 closest = (*it);
178 }
179 }
180 }
181 G4bool isDNAHit = (dist < fDNAGeometry->GetDirectInteractionRange());
182 if (isDNAHit) {
183 DoChromosomeDNAHit(pos, localPos, edep, dist, closest->GetName(), pv->GetName());
184 DoChromosomeRegionHit(pos, edep, true);
185 }
186 }
187 }
188
189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
190
191 void SteppingAction::DoChromosomeRegionHit(const G4ThreeVector& pos, const G4double& edep,
192 const G4bool& dnahit)
193 {
194 // TODO: chromosome should probably be based on position of center of
195 // parent volume, not of current pos
196 G4String chromosome = fDNAGeometry->GetChromosomeMapper()->GetCurrentChromosomeKey(pos);
197 if (!(chromosome.empty())) {
198 fEventAction->AddChromosomeEdep(chromosome, edep, dnahit);
199 }
200 }
201
202 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
203
204 void SteppingAction::DoChromosomeDNAHit(const G4ThreeVector& pos, const G4ThreeVector& localPos,
205 const G4double& edep, const G4double& dist,
206 const G4String& pvName, const G4String& motherName)
207 {
208 G4int placementIdx = std::stoi(utility::Split(motherName, '-').at(1));
209 std::vector<G4String> pvNameVec = utility::Split(pvName, '-');
210 molecule mol = utility::GetMoleculeEnum(pvNameVec.at(0));
211 G4int chainIdx = std::stoi(pvNameVec.at(1));
212 G4int strandIdx = std::stoi(pvNameVec.at(2));
213 int64_t baseIdx = std::stoll(pvNameVec.at(3));
214
215 // Convert to local chain and base index
216 chainIdx = fDNAGeometry->GetGlobalChain(placementIdx, chainIdx);
217 baseIdx += fDNAGeometry->GetStartIdx(placementIdx, chainIdx);
218
219 if (fDNAGeometry->GetStrandsFlipped(placementIdx)) {
220 strandIdx = (strandIdx + 1) % 2;
221 }
222 G4String chromosome = fDNAGeometry->GetChromosomeMapper()->GetCurrentChromosomeKey(pos);
223 const DNAHit* dnaHit = new DNAHit(mol, placementIdx, chainIdx, strandIdx, baseIdx, pos, localPos,
224 edep, dist, chromosome);
225
226 fEventAction->AddDNAHit(dnaHit);
227 }
228
229 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
230