| 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: DNAGeometry.hh
28 /// brief: This file handls MolecularDNA Geometry
29
30 #ifndef MOLECULAR_DNA_GEOMETRY_HH
31 #define MOLECULAR_DNA_GEOMETRY_HH
32
33 #include "ChromosomeMapper.hh"
34 #include "DNAWorld.hh"
35 #include "DamageModel.hh"
36 #include "MoleculeList.hh"
37 #include "UtilityFunctions.hh"
38
39 #include "G4LogicalVolume.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4ThreeVector.hh"
42 #include "G4VDNAMolecularGeometry.hh"
43 #include "G4VPhysicalVolume.hh"
44 #include "G4VTouchable.hh"
45 #include "globals.hh"
46
47 #include <array>
48 #include <map>
49 #include <unordered_map>
50 #include <vector>
51
52 class G4Material;
53
54 class G4MolecularConfiguration;
55
56 class G4VisAttributes;
57
58 class OctreeNode;
59
60 class DNAGeometryMessenger;
61
62 class PlacementVolumeInfo;
63
64 using placement = std::tuple<G4LogicalVolume*, int64_t, int64_t, G4ThreeVector,
65 G4ThreeVector>; // dousatsu
66 // G4int, G4int, G4ThreeVector, G4ThreeVector> placement;//ORG
67
68 // Holder to describe the global co-ordinates of each DNA chain
69 // in a given placement volume
70 // Format << array of global chain indices for local indices 0, 1, 2, 3>
71 // < array of start indices for the base pairs along each chain,
72 // given by GLOBAL chain index>
73 // < array of end indices for the base pairs along each chain,
74 // given by GLOBAL chain index>
75 // < Boolean to indicate if the strand ID should be flipped due
76 // to a 180 degree rotation of the DNA >
77 // < Boolean to indicate if volume has actually been placed >
78 using placement_transform =
79 std::tuple<std::array<int, 8>, std::array<int64_t, 8>, std::array<int64_t, 8>, G4bool, G4bool>;
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
81
82 struct molecule_t
83 {
84 G4String fname;
85 G4String fshape;
86 int64_t fchain_id;
87 int64_t fstrand_id;
88 int64_t fbase_idx;
89 G4ThreeVector fposition;
90 G4ThreeVector frotation;
91 G4ThreeVector fsize;
92 };
93
94 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
95
96 struct compareLVByName
97 {
98 bool operator()(const G4LogicalVolume* lhs, const G4LogicalVolume* rhs) const;
99 };
100
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
102
103 class DNAGeometry : public G4VDNAMolecularGeometry
104 {
105 public:
106 DNAGeometry();
107
108 ~DNAGeometry() override;
109
110 void BuildDNA(G4LogicalVolume*);
111
112 inline DNAWorld* GetDNAWorld();
113
114 void AddVoxelFile(const G4String& vname, const G4String& fname, const G4bool& twist)
115 {
116 fVoxelNames[vname] = fname;
117 fVoxelTwist[vname] = twist;
118 }
119
120 void SetFractalFilename(const G4String& fname)
121 {
122 fFractalCurveFile = fname;
123 FillVoxelVectors();
124 }
125
126 void SetFractalAnglesAsPi(const G4bool& value) { fAnglesAsPi = value; }
127
128 void EnableCustomMoleculeSizes(const G4bool& value) { fEnableCustomMoleculeSizes = value; }
129
130 void AddChangeMoleculeSize(const G4String& name, const G4ThreeVector& size);
131
132 inline void SetVoxelSideLength(const G4ThreeVector& length) { fVoxelSideLength = length; }
133
134 inline void SetFractalScaling(const G4ThreeVector& length) { fFractalScaling = length; }
135
136 OctreeNode* GetTopOctreeNode(G4LogicalVolume*) const;
137
138 const PlacementVolumeInfo* GetPVInfo(const G4LogicalVolume*) const;
139
140 inline auto GetChromosomeMapper() const { return fpChromosomeMapper; };
141
142 inline void SetOverlaps(G4bool overlap) { fCheckOverlaps = overlap; };
143
144 inline G4bool GetOverlaps() const { return fCheckOverlaps; };
145
146 inline void SetDrawCellVolumes(G4bool b) { fDrawCellVolumes = b; };
147
148 inline G4bool GetDrawCellVolumes() const { return fDrawCellVolumes; };
149
150 inline void SetVerbosity(G4int i) { fVerbosity = i; };
151
152 inline G4int GetVerbosity() const { return fVerbosity; };
153
154 inline void SetSmartless(G4int i) { fSmartless = i; };
155
156 inline G4int GetSmartless() const { return fSmartless; };
157
158 inline void SetMediumMaterial(G4Material* mat) { fpMediumMaterial = mat; };
159
160 inline G4LogicalVolume* GetDNAChemVolumePointer() const;
161
162 inline G4LogicalVolume* GetDNAPhysicsVolumePointer() const;
163
164 G4LogicalVolume* GetMatchingChemVolume(G4LogicalVolume*) const;
165
166 G4LogicalVolume* GetMatchingPhysVolume(G4LogicalVolume*) const;
167
168 void PrintOctreeStats();
169
170 inline void SetDirectInteractionRange(G4double r) { fDirectInteractionRange = r; };
171
172 void SetRadicalKillDistance(G4double r) { fRadicalKillDistance = r; };
173
174 void SetHistoneScav(G4bool hf) { fUseHistoneScav = hf; };
175
176 inline G4double GetDirectInteractionRange() const { return fDirectInteractionRange; };
177
178 G4double GetRadicalKillDistance() const { return fRadicalKillDistance; };
179
180 const DamageModel* GetDamageModel() { return fpDamageModel; };
181
182 // Tests
183 void ChromosomeTest();
184
185 void BasePairIndexTest();
186
187 void UniqueIDTest();
188
189 int64_t GetGlobalUniqueIDTest( // dousatsu
190 int64_t, int64_t, int64_t, int64_t, int64_t) const; // dousatsu
191
192 private:
193 DNAGeometryMessenger* fpMessenger;
194 G4String fFractalCurveFile;
195 G4bool fAnglesAsPi = false;
196 G4bool fEnableCustomMoleculeSizes = false;
197 G4bool fCheckOverlaps = false;
198 G4bool fDrawCellVolumes = false;
199 G4int fVerbosity = 0;
200 G4bool fUseHistoneScav = true;
201 // vis attrs to draw cells
202 G4VisAttributes* fpDrawCellAttrs;
203 // Materials
204 G4Material* fpMediumMaterial;
205 G4Material *fpSugarMaterial, *fpPhosphateMaterial;
206 G4Material *fpGuanineMaterial, *fpAdenineMaterial, *fpThymineMaterial;
207 G4Material *fpCytosineMaterial, *fpHistoneMaterial;
208
209 // Placement of molecules
210 void FillParameterisedSpace();
211
212 G4VPhysicalVolume* PlacePhosphate(G4LogicalVolume*, const molecule_t&, const molecule_t&);
213
214 G4VPhysicalVolume* PlaceSugar(G4LogicalVolume*, const molecule_t&, const molecule_t&);
215
216 G4VPhysicalVolume* PlaceBase(G4LogicalVolume*, const molecule_t&, const molecule_t&,
217 const molecule_t&, const molecule_t&);
218
219 G4VPhysicalVolume* PlaceHistone(G4LogicalVolume*, // dousatsu
220 const molecule_t&);
221
222 // defining voxels
223 std::pair<G4LogicalVolume*, PlacementVolumeInfo*> LoadVoxelVolume(const G4String&,
224 const G4String&);
225
226 void FillVoxelVectors();
227
228 G4LogicalVolume *fpDNAVolumePhys{}, *fpDNAVolumeChem{};
229 std::map<const G4LogicalVolume*, G4LogicalVolume*> fPhysToChem;
230 std::map<const G4LogicalVolume*, G4LogicalVolume*> fChemToPhys;
231 std::vector<G4ThreeVector> fVoxelPositions;
232 std::vector<G4ThreeVector> fVoxelRotations;
233 std::vector<G4String> fVoxelTypes;
234 std::vector<G4int> fVoxelIndices;
235 // std::vector<G4LogicalVolume *> fVoxelLogicalVolumes;
236 G4ThreeVector fVoxelSideLength = G4ThreeVector(75 * nm, 75 * nm, 75 * nm);
237 G4ThreeVector fFractalScaling = G4ThreeVector(1, 1, 1);
238 std::map<G4String, G4String> fVoxelNames; // For these, the first el.
239 std::map<G4String, G4bool> fVoxelTwist; // is also lv->GetName()
240
241 std::map<const G4LogicalVolume*, PlacementVolumeInfo*> fInfoMap;
242 std::map<G4String, G4ThreeVector> fMoleculeSizes;
243 G4double fSmartless = 2.0;
244 G4double fRadicalKillDistance = 10 * nm, fDirectInteractionRange = 6 * angstrom;
245 G4double fHistoneInteractionRadius = 25 * angstrom;
246
247 ChromosomeMapper* fpChromosomeMapper;
248 DamageModel* fpDamageModel;
249 DNAWorld* fpDNAPhysicsWorld;
250
251 // Members for placement transformations
252 public:
253 inline int64_t GetNumberOfPlacements() const
254 {
255 return (int64_t)fPlacementTransformations.size();
256 }
257
258 G4int GetGlobalChain(G4int vol_idx, G4int local_chain) const
259 {
260 return std::get<0>(fPlacementTransformations[vol_idx])[local_chain];
261 };
262
263 G4int GetLocalChain(G4int vol_idx, G4int global_chain) const;
264
265 inline int64_t GetStartIdx(int64_t vol_idx,
266 int64_t global_chn) const // dousatsu
267 {
268 return std::get<1>(fPlacementTransformations[vol_idx])[global_chn];
269 };
270
271 inline int64_t GetEndIdx(int64_t vol_idx,
272 int64_t global_chn) const // dousatsu
273 {
274 return std::get<2>(fPlacementTransformations[vol_idx])[global_chn];
275 };
276
277 int64_t GetMaxBPIdx() const;
278
279 inline G4int GetNumberOfChains() const
280 {
281 return std::get<0>(fPlacementTransformations[0]).size();
282 };
283
284 inline G4bool GetStrandsFlipped(G4int vol_idx) const
285 {
286 return std::get<3>(fPlacementTransformations[vol_idx]);
287 };
288
289 private:
290 std::vector<placement_transform> fPlacementTransformations;
291
292 void AddNewPlacement(const G4LogicalVolume*, std::array<int, 8>, G4bool, G4bool);
293
294 void AddFourChainPlacement(std::vector<placement>::iterator, std::vector<placement>::iterator,
295 G4bool);
296
297 void AddSingleChainPlacement(std::vector<placement>::iterator, std::vector<placement>::iterator,
298 G4bool);
299
300 // Lookup methods for chemistry handling
301 public:
302 // Unique ID generators
303 inline int64_t GetGlobalPairID(G4int, G4int, int64_t) const;
304
305 int64_t GetGlobalUniqueID(G4VPhysicalVolume*, // dousatsu
306 const G4VTouchable*) const;
307
308 // And reconstructors from unique ID
309 int64_t GetBasePairFromUniqueID(int64_t) const;
310
311 static inline molecule GetMoleculeFromUniqueID(int64_t);
312
313 G4Material* GetMaterialFromUniqueID(int64_t) const;
314
315 G4int GetChainIDFromUniqueID(int64_t) const;
316
317 G4int GetStrandIDFromUniqueID(int64_t) const;
318
319 G4int GetPlacementIndexFromUniqueID(int64_t) const;
320
321 void FindNearbyMolecules(const G4LogicalVolume* motherLogical,
322 const G4ThreeVector& localPosition,
323 std::vector<G4VPhysicalVolume*>& out, double searchRange) override;
324
325 G4bool IsInsideHistone(const G4LogicalVolume* motherLogical,
326 const G4ThreeVector& localPosition) const;
327 };
328
329 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
330
331 inline int64_t DNAGeometry::GetGlobalPairID(G4int place_idx, G4int chain_idx,
332 int64_t base_idx) const
333 {
334 return base_idx + this->GetStartIdx(place_idx, chain_idx);
335 }
336
337 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
338
339 inline molecule DNAGeometry::GetMoleculeFromUniqueID(int64_t idx)
340 {
341 return (molecule)(idx % ::molecule::Count);
342 }
343
344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
345
346 inline G4LogicalVolume* DNAGeometry::GetDNAPhysicsVolumePointer() const
347 {
348 return fpDNAVolumePhys;
349 }
350 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
351
352 inline DNAWorld* DNAGeometry::GetDNAWorld()
353 {
354 return fpDNAPhysicsWorld;
355 }
356
357 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
358
359 inline G4LogicalVolume* DNAGeometry::GetDNAChemVolumePointer() const
360 {
361 return fpDNAVolumeChem;
362 }
363 #endif // MOLECULAR_DNA_GEOMETRY_HH
364