| 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 // This example is provided by the Geant4-DNA
27 // Any report or published results obtained us
28 // and the DNA geometry given in the Geom_DNA
29 // shall cite the following Geant4-DNA collabo
30 // [1] NIM B 298 (2013) 47-54
31 // [2] Med. Phys. 37 (2010) 4692-4708
32 // The Geant4-DNA web site is available at htt
33 //
34 /// \file DetectorConstruction.cc
35 /// \brief Implementation of the DetectorConst
36
37 #include "DetectorConstruction.hh"
38
39 // Geant4
40 #include "CLHEP/Units/SystemOfUnits.h"
41 #include "ChromosomeParameterisation.hh"
42
43 #include "G4Box.hh"
44 #include "G4Ellipsoid.hh"
45 #include "G4LogicalVolume.hh"
46 #include "G4NistManager.hh"
47 #include "G4Orb.hh"
48 #include "G4PVParameterised.hh"
49 #include "G4PVPlacement.hh"
50 #include "G4RotationMatrix.hh"
51 #include "G4Tubs.hh"
52 #include "G4UnionSolid.hh"
53 #include "G4VisAttributes.hh"
54 #include "globals.hh"
55
56 #define countof(x) (sizeof(x) / sizeof(x[0]))
57
58 using namespace std;
59 using CLHEP::degree;
60 using CLHEP::micrometer;
61 using CLHEP::mm;
62 using CLHEP::nanometer;
63
64 static G4VisAttributes visInvBlue(false, G4Col
65 static G4VisAttributes visInvWhite(false, G4Co
66 static G4VisAttributes visInvPink(false, G4Col
67 static G4VisAttributes visInvCyan(false, G4Col
68 static G4VisAttributes visInvRed(false, G4Colo
69 static G4VisAttributes visInvGreen(false, G4Co
70 static G4VisAttributes visBlue(true, G4Colour(
71 static G4VisAttributes visWhite(true, G4Colour
72 static G4VisAttributes visPink(true, G4Colour(
73 static G4VisAttributes visCyan(true, G4Colour(
74 static G4VisAttributes visRed(true, G4Colour(1
75 static G4VisAttributes visGreen(true, G4Colour
76
77 //....oooOO0OOooo........oooOO0OOooo........oo
78
79 DetectorConstruction::DetectorConstruction()
80 : G4VUserDetectorConstruction(), fBuildChrom
81 {}
82
83 //....oooOO0OOooo........oooOO0OOooo........oo
84
85 DetectorConstruction::~DetectorConstruction()
86
87 //....oooOO0OOooo........oooOO0OOooo........oo
88
89 G4VPhysicalVolume* DetectorConstruction::Const
90 {
91 DefineMaterials();
92 return ConstructDetector();
93 }
94
95 //....oooOO0OOooo........oooOO0OOooo........oo
96
97 void DetectorConstruction::DefineMaterials()
98 {
99 // Water is defined from NIST material datab
100 G4NistManager* man = G4NistManager::Instance
101 man->FindOrBuildMaterial("G4_WATER");
102 }
103
104 //....oooOO0OOooo........oooOO0OOooo........oo
105
106 void DetectorConstruction::LoadChromosome(cons
107 G4Lo
108 {
109 ChromosomeParameterisation* cp = new Chromos
110 new G4PVParameterised("box ros", logicBoxros
111
112 G4cout << filename << " done" << G4endl;
113 }
114
115 //....oooOO0OOooo........oooOO0OOooo........oo
116
117 G4VPhysicalVolume* DetectorConstruction::Const
118 {
119 if (fBuildBases == false && fBuildChromatine
120 G4cout << "===============================
121 G4cout << "WARNING from DetectorConstructi
122 G4cout << "As long as the flags fBuildBase
123 "false, the output root file wil
124 << G4endl;
125 G4cout << "This is intended for fast compu
126 G4cout << "===============================
127 }
128
129 G4String name;
130
131 /*******************************************
132 // World
133 /*******************************************
134
135 DefineMaterials();
136 G4Material* waterMaterial = G4Material::GetM
137
138 G4Box* solidWorld = new G4Box("world", 10.0
139 G4LogicalVolume* logicWorld = new G4LogicalV
140 G4PVPlacement* physiWorld =
141 new G4PVPlacement(0, G4ThreeVector(), "wor
142 logicWorld->SetVisAttributes(&visInvWhite);
143
144 /*******************************************
145 // Box nucleus
146 /*******************************************
147
148 G4Box* solidTin = new G4Box("tin", 13 * micr
149 G4LogicalVolume* logicTin = new G4LogicalVol
150 G4VPhysicalVolume* physiTin =
151 new G4PVPlacement(0, G4ThreeVector(), "tin
152 logicTin->SetVisAttributes(&visInvWhite);
153
154 /*******************************************
155 // Cell nucleus
156 /*******************************************
157
158 G4Ellipsoid* solidNucleus =
159 new G4Ellipsoid("nucleus", 11.82 * microme
160 G4LogicalVolume* logicNucleus = new G4Logica
161 G4VPhysicalVolume* physiNucleus =
162 new G4PVPlacement(0, G4ThreeVector(), "phy
163 logicNucleus->SetVisAttributes(&visPink);
164
165 /*******************************************
166 // Chromosomes territ
167 /*******************************************
168 // NOTE: The only supported values for the r
169 // 0 and 90 degrees on the Y axis.
170 G4double chromosomePositionSizeRotation[][7]
171 {4.467, 2.835, 0, 1.557, 1.557, 1.557, 90}
172 {-4.467, 2.835, 0, 1.557, 1.557, 1.557, 0}
173 {4.423, -2.831, 0, 1.553, 1.553, 1.553, 90
174 {-4.423, -2.831, 0, 1.553, 1.553, 1.553, 0
175 {1.455, 5.63, 0, 1.455, 1.455, 1.455, 0},
176 {-1.455, 5.63, 0, 1.455, 1.455, 1.455, 90}
177 {1.435, 0, 1.392, 1.435, 1.435, 1.435, 0},
178 {-1.435, 0, 1.392, 1.435, 1.435, 1.435, 90
179 {1.407, 0, -1.450, 1.407, 1.407, 1.407, 90
180 {-1.407, 0, -1.450, 1.407, 1.407, 1.407, 0
181 {1.380, -5.437, 0, 1.380, 1.380, 1.380, 0}
182 {-1.380, -5.437, 0, 1.380, 1.380, 1.380, 9
183 {1.347, 2.782, -1.150, 1.347, 1.347, 1.347
184 {-1.347, 2.782, -1.150, 1.347, 1.347, 1.34
185 {1.311, -2.746, -1.220, 1.311, 1.311, 1.31
186 {-1.311, -2.746, -1.220, 1.311, 1.311, 1.3
187 {7.251, -2.541, 0, 1.275, 1.275, 1.275, 0}
188 {-6.701, 0, -0.85, 1.275, 1.275, 1.275, 90
189 {4.148, 0, 1.278, 1.278, 1.278, 1.278, 90}
190 {-4.148, 0, 1.278, 1.278, 1.278, 1.278, 0}
191 {4.147, 0, -1.277, 1.277, 1.277, 1.277, 0}
192 {-4.147, 0, -1.277, 1.277, 1.277, 1.277, 9
193 {8.930, 0.006, 0, 1.272, 1.272, 1.272, 90}
194 {-7.296, 2.547, 0, 1.272, 1.272, 1.272, 90
195 {1.207, -2.642, 1.298, 1.207, 1.207, 1.207
196 {-1.207, -2.642, 1.298, 1.207, 1.207, 1.20
197 {1.176, 2.611, 1.368, 1.176, 1.176, 1.176,
198 {-1.176, 2.611, 1.368, 1.176, 1.176, 1.176
199 {4.065, 5.547, 0, 1.155, 1.155, 1.155, 90}
200 {-4.065, 5.547, 0, 1.155, 1.155, 1.155, 0}
201 {6.542, 0.159, 1.116, 1.116, 1.116, 1.116,
202 {-9.092, 0, 0, 1.116, 1.116, 1.116, 0},
203 {6.507, 0.159, -1.081, 1.081, 1.081, 1.081
204 {-7.057, -2.356, 0, 1.081, 1.081, 1.081, 9
205 {3.824, -5.448, 0, 1.064, 1.064, 1.064, 90
206 {-3.824, -5.448, 0, 1.064, 1.064, 1.064, 0
207 {5.883, -5.379, 0, 0.995, 0.995, 0.995, 0}
208 {-9.133, -2.111, 0, 0.995, 0.995, 0.995, 0
209 {6.215, 5.387, 0, 0.995, 0.995, 0.995, 0},
210 {-6.971, -4.432, 0, 0.995, 0.995, 0.995, 9
211 {9.583, 2.177, 0, 0.899, 0.899, 0.899, 90}
212 {-9.467, 2.03, 0, 0.899, 0.899, 0.899, 0},
213 {9.440, -2.180, 0, 0.914, 0.914, 0.914, 90
214 {-6.34, 0, 1.339, 0.914, 0.914, 0.914, 0},
215 {-6.947, 4.742, 0, 0.923, 0.923, 0.923, 90
216 {7.354, 2.605, 0, 1.330, 1.330, 1.330, 0}
217 };
218
219 G4RotationMatrix* rotch = new G4RotationMatr
220 rotch->rotateY(90 * degree);
221
222 vector<G4VPhysicalVolume*> physiBox(48);
223
224 for (unsigned int i = 0; i < countof(chromos
225 G4double* p = &chromosomePositionSizeRotat
226 G4double* size = &chromosomePositionSizeRo
227 G4double rotation = chromosomePositionSize
228 G4ThreeVector pos(p[0] * micrometer, p[1]
229 G4RotationMatrix* rot = rotation == 0 ? 0
230
231 ostringstream ss;
232 ss << "box" << (i / 2) + 1 << (i % 2 ? 'l'
233 name = ss.str();
234 ss.str("");
235 ss.clear();
236
237 /*
238 snprintf(name, countof(name), "box%d%c",
239 (i / 2) + 1, i % 2 ? 'l' : 'r');
240 */
241 G4Box* solidBox =
242 new G4Box(name, size[0] * micrometer, si
243 G4LogicalVolume* logicBox = new G4LogicalV
244 physiBox[i] = new G4PVPlacement(rot, pos,
245 logicBox->SetVisAttributes(&visBlue);
246 }
247
248 /*******************************************
249 // Box containing the chroma
250 /*******************************************
251
252 G4Tubs* solidBoxros = new G4Tubs("solid box
253 0 * degree,
254 G4LogicalVolume* logicBoxros = new G4Logical
255 logicBoxros->SetVisAttributes(&visInvBlue);
256
257 // Loading flower box position for each chro
258
259 for (int k = 0; k < 22; k++) {
260 ostringstream oss;
261 oss << "chromo" << k + 1 << ".dat";
262 name = oss.str();
263 oss.str("");
264 oss.clear();
265 // snprintf(name, countof(name), "chromo%d
266 LoadChromosome(name.c_str(), physiBox[k *
267 LoadChromosome(name.c_str(), physiBox[k *
268 }
269
270 LoadChromosome("chromoY.dat", physiBox[44],
271 LoadChromosome("chromoX.dat", physiBox[45],
272
273 /*******************************************
274 if (fBuildChromatineFiber) {
275 // chromatin fiber envelope
276 G4Tubs* solidEnv = new G4Tubs("chromatin f
277 0 * degree,
278 G4LogicalVolume* logicEnv = new G4LogicalV
279 logicEnv->SetVisAttributes(&visInvPink);
280
281 // Chromatin fiber position
282 for (G4int i = 0; i < 7; i++) {
283 G4RotationMatrix* rotFiber = new G4Rotat
284 rotFiber->rotateX(90 * degree);
285 rotFiber->rotateY(i * 25.72 * degree);
286 G4ThreeVector posFiber = G4ThreeVector(0
287 posFiber.rotateZ(i * 25.72 * degree);
288 new G4PVPlacement(rotFiber, posFiber, lo
289
290 rotFiber = new G4RotationMatrix;
291 rotFiber->rotateX(90 * degree);
292 rotFiber->rotateY((7 + i) * 25.72 * degr
293 posFiber = G4ThreeVector(0, 152 * nanome
294 posFiber.rotateZ((7 + i) * 25.72 * degre
295 new G4PVPlacement(rotFiber, posFiber, lo
296
297 rotFiber = new G4RotationMatrix;
298 rotFiber->rotateX(90 * degree);
299 rotFiber->rotateY((25.72 + (i - 14) * 51
300 posFiber = G4ThreeVector(-36.5 * nanomet
301 posFiber.rotateZ((i - 14) * 51.43 * degr
302 new G4PVPlacement(rotFiber, posFiber, lo
303
304 rotFiber = new G4RotationMatrix;
305 rotFiber->rotateX(90 * degree);
306 rotFiber->rotateY(180 * degree);
307 rotFiber->rotateY((i - 21) * 51.43 * deg
308 posFiber = G4ThreeVector(-103 * nanomete
309 posFiber.rotateZ((i - 21) * 51.43 * degr
310 new G4PVPlacement(rotFiber, posFiber, lo
311 }
312
313 if (fBuildBases) {
314 // Histones
315 G4Tubs* solidHistone = new G4Tubs("solid
316 0 * de
317 G4LogicalVolume* logicHistone =
318 new G4LogicalVolume(solidHistone, wate
319
320 // Base pair
321 G4Orb* solidBp1 = new G4Orb("blue sphere
322 G4LogicalVolume* logicBp1 = new G4Logica
323 G4Orb* solidBp2 = new G4Orb("pink sphere
324 G4LogicalVolume* logicBp2 = new G4Logica
325
326 // Phosphodiester group
327
328 G4Orb* solidSugar_48em1_nm = new G4Orb("
329
330 G4ThreeVector posi(0.180248 * nanometer,
331 G4UnionSolid* uniDNA =
332 new G4UnionSolid("move", solidSugar_48
333
334 G4ThreeVector posi2(-0.128248 * nanomete
335 G4UnionSolid* uniDNA2 =
336 new G4UnionSolid("move2", solidSugar_4
337
338 /***************************************
339 Phosphodiester group Position
340 ***************************************
341
342 for (G4int n = 2; n < 200; n++) {
343 G4double SP1[2][3] = {
344 {(-0.6 * nanometer) * cos(n * 0.26),
345 {(0.6 * nanometer) * cos(n * 0.26),
346 G4double matriceSP1[3][3] = {
347 {cos(n * 0.076), -sin(n * 0.076), 0}
348 G4double matriceSP2[2][3];
349
350 for (G4int i = 0; i < 3; i++) {
351 G4double sumSP1 = 0;
352 G4double sumSP2 = 0;
353 for (G4int j = 0; j < 3; j++) {
354 sumSP1 += matriceSP1[i][j] * SP1[0
355 sumSP2 += matriceSP1[i][j] * SP1[1
356 }
357 matriceSP2[0][i] = sumSP1;
358 matriceSP2[1][i] = sumSP2;
359 }
360
361 G4double heliceSP[3] = {(4.85 * nanome
362 (4.85 * nanome
363
364 for (G4int i = 0; i < 3; i++) {
365 matriceSP2[0][i] += heliceSP[i];
366 matriceSP2[1][i] += heliceSP[i];
367 }
368 G4ThreeVector posSugar1(matriceSP2[0][
369 (matriceSP2[0]
370 G4ThreeVector posSugar2(matriceSP2[1][
371 (matriceSP2[1]
372
373 ostringstream ss;
374 ss << "sugar_" << n;
375 name = ss.str().c_str();
376 ss.str("");
377 ss.clear();
378
379 // snprintf(name, countof(name), "sug
380 uniDNA = new G4UnionSolid(name, uniDNA
381
382 ss << "sugar_" << n;
383 name = ss.str().c_str();
384 ss.str("");
385 ss.clear();
386
387 // snprintf(name, countof(name), "sug
388 uniDNA2 = new G4UnionSolid(name, uniDN
389 }
390 G4LogicalVolume* logicSphere3 = new G4Lo
391 G4LogicalVolume* logicSphere4 = new G4Lo
392
393 /***************************************
394 Base pair Position
395 ***************************************
396 for (G4int n = 0; n < 200; n++) {
397 G4double bp1[2][3] = {
398 {(-0.34 * nanometer) * cos(n * 0.26)
399 {(0.34 * nanometer) * cos(n * 0.26),
400 G4double matriceBP1[3][3] = {
401 {cos(n * 0.076), -sin(n * 0.076), 0}
402 G4double matriceBP2[2][3];
403
404 for (G4int i = 0; i < 3; i++) {
405 G4double sumBP1 = 0;
406 G4double sumBP2 = 0;
407 for (G4int j = 0; j < 3; j++) {
408 sumBP1 += matriceBP1[i][j] * bp1[0
409 sumBP2 += matriceBP1[i][j] * bp1[1
410 }
411 matriceBP2[0][i] = sumBP1;
412 matriceBP2[1][i] = sumBP2;
413 }
414 G4double heliceBP[3] = {(4.8 * nanomet
415 (4.8 * nanomet
416
417 for (G4int i = 0; i < 3; i++) {
418 matriceBP2[0][i] += heliceBP[i];
419 matriceBP2[1][i] += heliceBP[i];
420 }
421 G4ThreeVector position1(matriceBP2[0][
422 matriceBP2[0][
423 G4ThreeVector position2(matriceBP2[1][
424 matriceBP2[1][
425
426 new G4PVPlacement(0, position1, logicB
427 new G4PVPlacement(0, position2, logicB
428 }
429
430 /***************************************
431 // Initial position of d
432 /***************************************
433 // DNA and histone positions
434 for (int j = 0; j < 90; j++) {
435 // DNA (bp-SP)
436 G4RotationMatrix* rotStrand1 = new G4R
437 rotStrand1->rotateZ(j * -51.43 * degre
438 G4ThreeVector posStrand1(-2.7 * nanome
439 (-69.9 * nano
440 posStrand1.rotateZ(j * 51.43 * degree)
441 new G4PVPlacement(rotStrand1, posStran
442 0);
443
444 G4RotationMatrix* rotStrand2 = new G4R
445 rotStrand2->rotateZ(j * -51.43 * degre
446 G4ThreeVector posStrand2(-2.7 * nanome
447 (-68.7 * nano
448 posStrand2.rotateZ(j * 51.43 * degree)
449 new G4PVPlacement(rotStrand2, posStran
450 0);
451
452 // histones
453 G4RotationMatrix* rotHistone = new G4R
454 rotHistone->rotateY(90 * degree);
455 rotHistone->rotateX(j * (-51.43 * degr
456 G4ThreeVector posHistone(0.0, 9.35 * n
457 posHistone.rotateZ(j * 51.43 * degree)
458 new G4PVPlacement(rotHistone, posHisto
459 }
460 /***************************************
461 // Visualisation
462 /***************************************
463
464 logicBp1->SetVisAttributes(&visInvCyan);
465 logicBp2->SetVisAttributes(&visInvPink);
466
467 logicSphere3->SetVisAttributes(&visInvWh
468 logicSphere4->SetVisAttributes(&visInvRe
469
470 logicHistone->SetVisAttributes(&visInvBl
471 }
472 }
473
474 G4cout << "Geometry has been loaded" << G4en
475 return physiWorld;
476 }
477