| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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........oo
50
51 SteppingAction::SteppingAction(EventAction* ev
52 {
53 const auto* det = dynamic_cast<const Detecto
54 G4RunManager::GetRunManager()->GetUserDete
55
56 fDNAGeometry = det->GetDNAGeometry();
57 }
58
59 //....oooOO0OOooo........oooOO0OOooo........oo
60
61 void SteppingAction::UserSteppingAction(const
62 {
63 G4double edep = aStep->GetTotalEnergyDeposit
64 if (edep == 0) {
65 return;
66 }
67
68 // Work out if we are in the Logical volume
69 const G4VTouchable* touch = aStep->GetPreSte
70
71 // remember the chem volume is the "real" vo
72 G4LogicalVolume* dnaVol = fDNAGeometry->GetD
73 G4int depth = touch->GetHistoryDepth();
74 G4bool isInDNARegion = false;
75
76 if (depth >= 0) {
77 fEventAction->AddCellEdep(edep); // dousa
78 }
79
80 while ((depth >= 0) && !isInDNARegion) {
81 if (touch->GetVolume(depth) != nullptr) {
82 isInDNARegion = isInDNARegion || (touch-
83 }
84 else {
85 G4cout << "Null pointer volume in mother
86 << " looking for volume at depth
87 assert(depth < 0);
88 }
89 depth--;
90 }
91
92 if (isInDNARegion) {
93 ProcessDNARegionHit(aStep);
94 }
95
96 G4ThreeVector prepos = aStep->GetPreStepPoin
97 if (fDNAGeometry->GetChromosomeMapper()->Get
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()->G
106 fEventAction->AddTrackLengthChro(tl);
107 }
108 }
109 }
110
111 //....oooOO0OOooo........oooOO0OOooo........oo
112
113 void SteppingAction::ProcessDNARegionHit(const
114 {
115 G4double edep = aStep->GetTotalEnergyDeposit
116 const G4VTouchable* touch = aStep->GetPreSte
117 G4ThreeVector pos =
118 aStep->GetPreStepPoint()->GetPosition()
119 + G4UniformRand()
120 * (aStep->GetPostStepPoint()->GetPosit
121 G4VPhysicalVolume* pv = aStep->GetPreStepPoi
122 G4LogicalVolume* lv = pv->GetLogicalVolume()
123
124 std::vector<G4VPhysicalVolume*> molvec;
125 OctreeNode* parentNode;
126 G4AffineTransform transform = touch->GetHist
127 parentNode = fDNAGeometry->GetTopOctreeNode(
128 if (parentNode == nullptr) { // See if pare
129 parentNode = fDNAGeometry->GetTopOctreeNod
130 if (parentNode != nullptr) { // If mother
131 // These methods are unsafe if the pv is
132 // need to be in this if statement.
133 // lv = pv->GetMotherLogical();
134 pv = touch->GetVolume(1);
135 transform = touch->GetHistory()->GetTran
136 }
137 }
138
139 G4ThreeVector localPos = transform.Transform
140
141 if (parentNode != nullptr) { // Get Molecul
142 molvec = parentNode->SearchOctree(localPos
143 }
144
145 // debug output for finding logical volumes.
146 if (fDNAGeometry->GetVerbosity() > 1) {
147 if (parentNode != nullptr) {
148 G4String lvname = pv->GetLogicalVolume()
149 G4cout << "Found octree for logical volu
150 }
151 else {
152 // Don't bother doing error printouts if
153 // is at the DNA mother volume level.
154 // Remember the Chemistry Volume is the
155 if (pv->GetLogicalVolume() != fDNAGeomet
156 G4String lvname = pv->GetLogicalVolume
157 G4String motherlvname;
158 if (lv != nullptr) {
159 motherlvname = lv->GetName();
160 }
161 G4cout << "Could not find octree for l
162 << "Particle in LV: " << lvname
163 }
164 }
165 }
166
167 if (!molvec.empty()) {
168 // NOTE: This loop can be optimised to avo
169 G4VPhysicalVolume* closest = molvec[0];
170 G4double dist = (closest->GetTranslation()
171 if (molvec.size() > 1) {
172 G4double newdist;
173 for (auto it = molvec.begin() + 1; it !=
174 newdist = ((*it)->GetTranslation() - l
175 if (newdist < dist) {
176 dist = newdist;
177 closest = (*it);
178 }
179 }
180 }
181 G4bool isDNAHit = (dist < fDNAGeometry->Ge
182 if (isDNAHit) {
183 DoChromosomeDNAHit(pos, localPos, edep,
184 DoChromosomeRegionHit(pos, edep, true);
185 }
186 }
187 }
188
189 //....oooOO0OOooo........oooOO0OOooo........oo
190
191 void SteppingAction::DoChromosomeRegionHit(con
192 con
193 {
194 // TODO: chromosome should probably be based
195 // parent volume, not of current pos
196 G4String chromosome = fDNAGeometry->GetChrom
197 if (!(chromosome.empty())) {
198 fEventAction->AddChromosomeEdep(chromosome
199 }
200 }
201
202 //....oooOO0OOooo........oooOO0OOooo........oo
203
204 void SteppingAction::DoChromosomeDNAHit(const
205 const
206 const
207 {
208 G4int placementIdx = std::stoi(utility::Spli
209 std::vector<G4String> pvNameVec = utility::S
210 molecule mol = utility::GetMoleculeEnum(pvNa
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(plac
217 baseIdx += fDNAGeometry->GetStartIdx(placeme
218
219 if (fDNAGeometry->GetStrandsFlipped(placemen
220 strandIdx = (strandIdx + 1) % 2;
221 }
222 G4String chromosome = fDNAGeometry->GetChrom
223 const DNAHit* dnaHit = new DNAHit(mol, place
224 edep, dist
225
226 fEventAction->AddDNAHit(dnaHit);
227 }
228
229 //....oooOO0OOooo........oooOO0OOooo........oo
230