Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
examples/advanced/dnageometry/src/DetectorConstruction.cc
Go to the documentation of this file.
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 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained using the Geant4-DNA software
28 // and the DNA geometry given in the Geom_DNA example
29 // shall cite the following Geant4-DNA collaboration publications:
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 http://geant4-dna.org
33 //
35 #include "DetectorConstruction.hh"
36 #include "Randomize.hh"
37 #include "time.h"
38 #include "G4Region.hh"
39 #include "G4ProductionCuts.hh"
40 #include "CLHEP/Random/Randomize.h"
41 #include "G4PVParameterised.hh"
42 #include <fstream>
43 #include <vector>
44 #include "G4SystemOfUnits.hh"
45 
46 #define countof(x) (sizeof(x) / sizeof(x[0]))
47 
48 using namespace std;
49 
50 static G4VisAttributes visInvBlue(false, G4Colour(0.0, 0.0, 1.0));
51 static G4VisAttributes visInvWhite(false, G4Colour(1.0, 1.0, 1.0));
52 static G4VisAttributes visInvPink(false, G4Colour(1.0, 0.0, 1.0));
53 static G4VisAttributes visInvCyan(false, G4Colour(0.0, 1.0, 1.0));
54 static G4VisAttributes visInvRed(false, G4Colour(1.0, 0.0, 0.0));
55 static G4VisAttributes visInvGreen(false, G4Colour(0.0, 1.0, 0.0));
56 static G4VisAttributes visBlue(true, G4Colour(0.0, 0.0, 1.0));
57 static G4VisAttributes visWhite(true, G4Colour(1.0, 1.0, 1.0));
58 static G4VisAttributes visPink(true, G4Colour(1.0, 0.0, 1.0));
59 static G4VisAttributes visCyan(true, G4Colour(0.0, 1.0, 1.0));
60 static G4VisAttributes visRed(true, G4Colour(1.0, 0.0, 0.0));
61 static G4VisAttributes visGreen(true, G4Colour(0.0, 1.0, 0.0));
62 
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64 
66 :physiWorld(NULL), logicWorld(NULL), solidWorld(NULL)
67 {}
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70 
72 {}
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 
77 {
78  DefineMaterials();
79  return ConstructDetector();
80 }
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 
84 void DetectorConstruction::DefineMaterials()
85 {
86  // Water is defined from NIST material database
88  waterMaterial = man->FindOrBuildMaterial("G4_WATER");
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 
93 void DetectorConstruction::LoadChromosome(const char* filename, G4VPhysicalVolume* chromBox)
94 {
96  new G4PVParameterised("box ros", logicBoxros, chromBox, kUndefined, cp->getNumRosettes(), cp);
97 
98  G4cout << filename << " done" << G4endl;
99 }
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102 
103 G4VPhysicalVolume* DetectorConstruction::ConstructDetector()
104 {
105  char name[256];
106 
107  /**************************************************************************************************************/
108  // World
109  /*************************************************************************************************************/
110 
111  solidWorld = new G4Box("world", 10.0*mm, 10.0*mm, 10.0*mm);
112  logicWorld = new G4LogicalVolume(solidWorld, waterMaterial, "world");
113  physiWorld = new G4PVPlacement(0, G4ThreeVector(), "world", logicWorld, NULL, false, 0);
114  logicWorld->SetVisAttributes(&visInvWhite);
115 
116  /**************************************************************************************************************/
117  // Box nucleus
118  /*************************************************************************************************************/
119 
120  solidTin = new G4Box("tin", 13*micrometer, 10*micrometer, 5*micrometer);
121  logicTin = new G4LogicalVolume(solidTin, waterMaterial, "tin");
122  physiTin = new G4PVPlacement(0, G4ThreeVector(), "tin", logicTin, physiWorld, false, 0);
123  logicTin->SetVisAttributes(&visInvWhite);
124 
125  /*************************************************************************************************************/
126  // Cell nucleus
127  /*************************************************************************************************************/
128 
129  solidNucleus = new G4Ellipsoid("nucleus", 11.82*micrometer, 8.52*micrometer, 3*micrometer, 0, 0);
130  logicNucleus = new G4LogicalVolume(solidNucleus, waterMaterial, "logic nucleus");
131  physiNucleus = new G4PVPlacement(0,G4ThreeVector(), "physi nucleus", logicNucleus, physiTin, false, 0);
132  logicNucleus->SetVisAttributes(&visPink);
133 
134  /*************************************************************************************************************/
135  // Chromosomes territories
136  /*************************************************************************************************************/
137  // NOTE: The only supported values for the rotation are 0 and 90 degrees on the Y axis.
138  float chromosomePositionSizeRotation[][7] = {
139  {4.467, 2.835, 0, 1.557, 1.557, 1.557, 90},
140  {-4.467, 2.835, 0, 1.557, 1.557, 1.557, 0},
141  {4.423, -2.831, 0, 1.553, 1.553, 1.553, 90},
142  {-4.423, -2.831, 0, 1.553, 1.553, 1.553, 0},
143  {1.455, 5.63, 0, 1.455, 1.455, 1.455, 0},
144  {-1.455, 5.63, 0, 1.455, 1.455, 1.455, 90},
145  {1.435, 0, 1.392, 1.435, 1.435, 1.435, 0},
146  {-1.435, 0, 1.392, 1.435, 1.435, 1.435, 90},
147  {1.407, 0, -1.450, 1.407, 1.407, 1.407, 90}, // 5 right
148  {-1.407, 0, -1.450, 1.407, 1.407, 1.407, 0}, // 5 left
149  {1.380, -5.437, 0, 1.380, 1.380, 1.380, 0},
150  {-1.380, -5.437, 0, 1.380, 1.380, 1.380, 90},
151  {1.347, 2.782, -1.150, 1.347, 1.347, 1.347, 90},
152  {-1.347, 2.782, -1.150, 1.347, 1.347, 1.347, 0},
153  {1.311, -2.746, -1.220, 1.311, 1.311, 1.311, 90},
154  {-1.311, -2.746, -1.220, 1.311, 1.311, 1.311, 0},
155  {7.251, -2.541, 0, 1.275, 1.275, 1.275, 0},
156  {-6.701, 0, -0.85, 1.275, 1.275, 1.275, 90},
157  {4.148, 0, 1.278, 1.278, 1.278, 1.278, 90}, // 10 right
158  {-4.148, 0, 1.278, 1.278, 1.278, 1.278, 0}, // 10 left
159  {4.147, 0, -1.277, 1.277, 1.277, 1.277, 0},
160  {-4.147, 0, -1.277, 1.277, 1.277, 1.277, 90},
161  {8.930, 0.006, 0, 1.272, 1.272, 1.272, 90},
162  {-7.296, 2.547, 0, 1.272, 1.272, 1.272, 90},
163  {1.207, -2.642, 1.298, 1.207, 1.207, 1.207, 0},
164  {-1.207, -2.642, 1.298, 1.207, 1.207, 1.207, 90},
165  {1.176, 2.611, 1.368, 1.176, 1.176, 1.176, 0},
166  {-1.176, 2.611, 1.368, 1.176, 1.176, 1.176, 90},
167  {4.065, 5.547, 0, 1.155, 1.155, 1.155, 90}, // 15 right
168  {-4.065, 5.547, 0, 1.155, 1.155, 1.155, 0}, // 15 left
169  {6.542, 0.159, 1.116, 1.116, 1.116, 1.116, 0},
170  {-9.092, 0, 0, 1.116, 1.116, 1.116, 0},
171  {6.507, 0.159, -1.081, 1.081, 1.081, 1.081, 90},
172  {-7.057, -2.356, 0, 1.081, 1.081, 1.081, 90},
173  {3.824, -5.448, 0, 1.064, 1.064, 1.064, 90},
174  {-3.824, -5.448, 0, 1.064, 1.064, 1.064, 0},
175  {5.883, -5.379, 0, 0.995, 0.995, 0.995, 0},
176  {-9.133, -2.111, 0, 0.995, 0.995, 0.995, 0},
177  {6.215, 5.387, 0, 0.995, 0.995, 0.995, 0}, // 20 right
178  {-6.971, -4.432, 0, 0.995, 0.995, 0.995, 90}, // 20 left
179  {9.583, 2.177, 0, 0.899, 0.899, 0.899, 90},
180  {-9.467, 2.03, 0, 0.899, 0.899, 0.899, 0},
181  {9.440, -2.180, 0, 0.914, 0.914, 0.914, 90},
182  {-6.34, 0, 1.339, 0.914, 0.914, 0.914, 0},
183  {-6.947, 4.742, 0, 0.923, 0.923, 0.923, 90}, // Y
184  {7.354, 2.605, 0, 1.330, 1.330, 1.330, 0} // X
185  };
186  G4RotationMatrix* rotch = new G4RotationMatrix;
187  rotch->rotateY(90 * degree);
188 
189  for (unsigned int i = 0; i < countof(chromosomePositionSizeRotation); i++)
190  {
191  float* p = &chromosomePositionSizeRotation[i][0];
192  float* size = &chromosomePositionSizeRotation[i][3];
193  float rotation = chromosomePositionSizeRotation[i][6];
194  G4ThreeVector pos(p[0] * micrometer, p[1] * micrometer, p[2] * micrometer);
195  G4RotationMatrix* rot = rotation == 0 ? 0 : rotch;
196 
197  snprintf(name, countof(name), "box%d%c", (i / 2) + 1, i % 2 ? 'l' : 'r');
198  G4Box* solidBox = new G4Box(name, size[0] * micrometer, size[1] * micrometer, size[2] * micrometer);
199  G4LogicalVolume* logicBox = new G4LogicalVolume(solidBox, waterMaterial, name);
200  physiBox[i] = new G4PVPlacement(rot, pos, "chromo", logicBox, physiNucleus, false, 0);
201  logicBox->SetVisAttributes(&visBlue);
202  }
203 
204  /***************************************************************************************************/
205 
206 // chromatin fiber envelope
207  G4Tubs* solidEnv = new G4Tubs("chromatin fiber", 0, 15.4 * nanometer, 80.5 * nanometer, 0 * degree, 360 * degree);
208  logicEnv = new G4LogicalVolume(solidEnv, waterMaterial, "LV chromatin fiber");
209  logicEnv->SetVisAttributes(&visInvPink);
210 
211 
212 // Histones
213  solidHistone = new G4Tubs("solid histone", 0, 3.25 * nanometer, 2.85 * nanometer, 0*degree,360*degree);
214  logicHistone = new G4LogicalVolume(solidHistone,waterMaterial,"logic histone");
215 
216 
217 //Base pair
218  solidBp1 = new G4Orb("blue sphere",0.17*nanometer);
219  logicBp1 = new G4LogicalVolume(solidBp1,waterMaterial,"logic blue sphere");
220  solidBp2 = new G4Orb("pink sphere",0.17*nanometer);
221  logicBp2 = new G4LogicalVolume(solidBp2,waterMaterial,"logic pink sphere");
222 
223 
224 //Phosphodiester group
225  solidSugar1 = new G4Orb("sugar 1", 0.48*nanometer);
226  solidSugar2 = new G4Orb("sugar 2", 0.48*nanometer);
227  G4ThreeVector posi(0.180248* nanometer,0.32422* nanometer,0.00784 * nanometer);
228  uniDNA = new G4UnionSolid("move", solidSugar1, solidSugar2, 0, posi);
229 
230  solidSugar3 = new G4Orb("sugar 3", 0.48*nanometer);
231  solidSugar4 = new G4Orb("sugar 4", 0.48*nanometer);
232  G4ThreeVector posi2(-0.128248* nanometer,0.41227* nanometer,0.03584 * nanometer);
233  uniDNA2 = new G4UnionSolid("move2", solidSugar3, solidSugar4, 0, posi2);
234 
235 
236 
237 /**********************************************************************************************************************
238  Phosphodiester group Position
239 **********************************************************************************************************************/
240  for(G4int n=2;n<200;n++)
241  {
242  G4float SP1[2][3]={{(-0.6*nanometer)*cos(n*0.26),0,(0.6*nanometer)*sin(n*0.26)},{(0.6*nanometer)*cos(n*0.26),0,(-0.6*nanometer)*sin(0.26*n)}};
243  G4float matriceSP1[3][3]={{cos(n*0.076),-sin(n*0.076),0},{sin(n*0.076),cos(n*0.076),0},{0,0,1}};
244  G4float matriceSP2[2][3];
245 
246  for(G4int i=0;i<3;i++)
247  {
248  G4float sumSP1=0;
249  G4float sumSP2=0;
250  for(G4int j=0;j<3;j++)
251  {
252  sumSP1+=matriceSP1[i][j]*SP1[0][j];
253  sumSP2+=matriceSP1[i][j]*SP1[1][j];
254  }
255  matriceSP2[0][i] = sumSP1;
256  matriceSP2[1][i] = sumSP2;
257  }
258  G4float heliceSP[3]={(4.85*nanometer)*cos(n*0.076),(4.85*nanometer)*sin(n*0.076),(n*0.026*nanometer)};
259  for(G4int i=0;i<3;i++)
260  {
261  matriceSP2[0][i]+=heliceSP[i];
262  matriceSP2[1][i]+=heliceSP[i];
263  }
264  G4ThreeVector posSugar1(matriceSP2[0][2],matriceSP2[0][1],(matriceSP2[0][0])-(4.25*nanometer));
265  G4ThreeVector posSugar2(matriceSP2[1][2],matriceSP2[1][1],(matriceSP2[1][0])-(5.45*nanometer));
266 
267  snprintf(name, countof(name), "sugar %d", n);
268  solidSphere3 = new G4Orb(name, 0.48*nanometer);
269  uniDNA = new G4UnionSolid("move", uniDNA, solidSphere3, 0, posSugar1);
270 
271  snprintf(name, countof(name), "sugar %d", n);
272  solidSphere4 = new G4Orb(name, 0.48*nanometer);
273  uniDNA2 = new G4UnionSolid("move2", uniDNA2, solidSphere4, 0, posSugar2);
274 }
275  logicSphere3 = new G4LogicalVolume(uniDNA,waterMaterial,"logic sugar 2");
276  logicSphere4 = new G4LogicalVolume(uniDNA2,waterMaterial,"logic sugar 4");
277 
278 /**********************************************************************************************************************
279  Base pair Position
280 **********************************************************************************************************************/
281  for(G4int n=0;n<200;n++)
282  {
283  G4float bp1[2][3]={{(-0.34*nanometer)*cos(n*0.26),0,(0.34*nanometer)*sin(n*0.26)},{(0.34*nanometer)*cos(n*0.26),0,(-0.34*nanometer)*sin(0.26*n)}};
284  G4float matriceBP1[3][3]={{cos(n*0.076),-sin(n*0.076),0},{sin(n*0.076),cos(n*0.076),0},{0,0,1}};
285  G4float matriceBP2[2][3];
286 
287  for(G4int i=0;i<3;i++)
288  {
289  G4float sumBP1=0;
290  G4float sumBP2=0;
291  for(G4int j=0;j<3;j++)
292  {
293  sumBP1+=matriceBP1[i][j]*bp1[0][j];
294  sumBP2+=matriceBP1[i][j]*bp1[1][j];
295  }
296  matriceBP2[0][i] = sumBP1;
297  matriceBP2[1][i] = sumBP2;
298  }
299  G4float heliceBP[3]={(4.8*nanometer)*cos(n*0.076),(4.8*nanometer)*sin(n*0.076),n*0.026*nanometer};
300 
301  for(G4int i=0;i<3;i++)
302  {
303  matriceBP2[0][i]+=heliceBP[i];
304  matriceBP2[1][i]+=heliceBP[i];
305  }
306  G4ThreeVector position1(matriceBP2[0][2],matriceBP2[0][1],matriceBP2[0][0]-(4.25*nanometer));
307  G4ThreeVector position2(matriceBP2[1][2],matriceBP2[1][1],matriceBP2[1][0]-(5.45*nanometer));
308 
309  physiBp1[n] = new G4PVPlacement(0,position1,logicBp1,"physi blue sphere",logicSphere3,false,0);
310  physiBp2[n] = new G4PVPlacement(0,position2,logicBp2,"physi pink sphere",logicSphere4,false,0);
311 }
312 
313 
314  /**************************************************************************************************************/
315  // Box continaing the chromatin flowers
316  /*************************************************************************************************************/
317 
318  solidBoxros = new G4Tubs("solid box ros", 0*nanometer, 399*nanometer, 20*nanometer, 0*degree, 360*degree);
319  logicBoxros = new G4LogicalVolume(solidBoxros, waterMaterial, "box ros");
320  logicBoxros->SetVisAttributes(&visInvBlue);
321 
322  /**************************************************************************************************************/
323  // Initial position of different elements
324  /*************************************************************************************************************/
325 
326  // Chromatin fiber position
327  for (G4int i = 0; i < 7; i++) {
328  G4RotationMatrix* rotFiber = new G4RotationMatrix;
329  rotFiber->rotateX(90*degree);
330  rotFiber->rotateY(i *25.72*degree);
331  G4ThreeVector posFiber = G4ThreeVector(0, 152*nanometer, 0);
332  posFiber.rotateZ(i*25.72*degree);
333  new G4PVPlacement(rotFiber, posFiber, logicEnv, "physi env", logicBoxros, false, 0);
334 
335  rotFiber = new G4RotationMatrix;
336  rotFiber->rotateX(90 * degree);
337  rotFiber->rotateY((7 + i) * 25.72 *degree);
338  posFiber = G4ThreeVector(0, 152 *nanometer, 0);
339  posFiber.rotateZ((7 + i) * 25.72*degree);
340  new G4PVPlacement(rotFiber, posFiber, logicEnv, "physi env", logicBoxros, false, 0);
341 
342  rotFiber = new G4RotationMatrix;
343  rotFiber->rotateX(90 * degree);
344  rotFiber->rotateY((25.72 + (i - 14) * 51.43) * degree);
345  posFiber = G4ThreeVector(-36.5 * nanometer, 312 * nanometer, 0);
346  posFiber.rotateZ((i - 14) * 51.43 * degree);
347  new G4PVPlacement(rotFiber, posFiber, logicEnv, "physi env", logicBoxros, false, 0);
348 
349  rotFiber = new G4RotationMatrix;
350  rotFiber->rotateX(90 * degree);
351  rotFiber->rotateY(180 * degree);
352  rotFiber->rotateY((i - 21) * 51.43 * degree);
353  posFiber = G4ThreeVector(-103 * nanometer, 297 * nanometer, 0);
354  posFiber.rotateZ((i - 21) * 51.43 * degree);
355  new G4PVPlacement(rotFiber, posFiber, logicEnv, "physi env", logicBoxros, false, 0);
356 
357 
358  }
359 
360 
361 // DNA and histone positions
362  for (int j = 0; j < 90; j++) {
363  // DNA (bp-SP)
364  G4RotationMatrix* rotStrand1 = new G4RotationMatrix;
365  rotStrand1->rotateZ(j*-51.43*degree);
366  G4ThreeVector posStrand1(-2.7*nanometer,9.35*nanometer,(-69.9*nanometer)+(j*1.67*nanometer));
367  posStrand1.rotateZ(j*51.43*degree);
368  new G4PVPlacement(rotStrand1,posStrand1,logicSphere3,"physi sugar 2",logicEnv,false,0);
369 
370  G4RotationMatrix* rotStrand2 = new G4RotationMatrix;
371  rotStrand2->rotateZ(j* -51.43*degree);
372  G4ThreeVector posStrand2(-2.7*nanometer,9.35*nanometer,(-68.7*nanometer)+(j*1.67*nanometer));
373  posStrand2.rotateZ(j*51.43*degree);
374  new G4PVPlacement(rotStrand2,posStrand2,logicSphere4,"physi sugar 4",logicEnv,false,0);
375 
376  // histones
377  G4RotationMatrix* rotHistone = new G4RotationMatrix;
378  rotHistone->rotateY(90*degree);
379  rotHistone->rotateX(j*(-51.43*degree));
380  G4ThreeVector posHistone(0.0,9.35 *nanometer,(-74.15+ j*1.67)*nanometer);
381  posHistone.rotateZ(j*51.43*degree);
382  new G4PVPlacement(rotHistone, posHistone, logicHistone, "PV histone", logicEnv, false, 0);
383  }
384 
385 //Loading flower box position for each chromosome territory
386 
387  for (int k = 0; k < 22; k++)
388  {
389  snprintf(name, countof(name), "chromo%d.dat", k + 1);
390  LoadChromosome(name, physiBox[k * 2]);
391  LoadChromosome(name, physiBox[k * 2 + 1]);
392  }
393 
394  LoadChromosome("chromoY.dat", physiBox[44]);
395  LoadChromosome("chromoX.dat", physiBox[45]);
396 
397  /**************************************************************************************************************/
398  // Visualisation colors
399  /*************************************************************************************************************/
400 
401  logicBp1->SetVisAttributes(&visInvCyan);
402  logicBp2->SetVisAttributes(&visInvPink);
403 
404  logicSphere3->SetVisAttributes(&visInvWhite);
405  logicSphere4->SetVisAttributes(&visInvRed);
406 
407  logicHistone->SetVisAttributes(&visInvBlue);
408 
409 
410  G4cout << "Geometry has been loaded" << G4endl;
411  return physiWorld;
412 }
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
CLHEP::Hep3Vector G4ThreeVector
HepRotation & rotateX(double delta)
Definition: Rotation.cc:66
CLHEP::HepRotation G4RotationMatrix
Definition: G4Box.hh:63
const char * p
Definition: xmltok.h:285
float G4float
Definition: G4Types.hh:77
Definition: G4Tubs.hh:84
const XML_Char * name
HepRotation & rotateY(double delta)
Definition: Rotation.cc:79
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
int nanometer
Definition: hepunit.py:35
G4GLOB_DLL std::ostream G4cout
tuple degree
Definition: hepunit.py:69
Hep3Vector & rotateZ(double)
Definition: ThreeVector.cc:144
const G4int n
Definition: G4Orb.hh:60
HepRotation & rotateZ(double delta)
Definition: Rotation.cc:92
#define G4endl
Definition: G4ios.hh:61
int micrometer
Definition: hepunit.py:34
void SetVisAttributes(const G4VisAttributes *pVA)