| 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 // This example is provided by the Geant4-DNA collaboration
27 // Any report or published results obtained using the Geant4-DNA software
28 // shall cite the following Geant4-DNA collaboration publication:
29 // Med. Phys. 37 (2010) 4692-4708
30 // Delage et al. PDB4DNA: implementation of DNA geometry from the Protein Data
31 // Bank (PDB) description for Geant4-DNA Monte-Carlo
32 // simulations (submitted to Comput. Phys. Commun.)
33 // The Geant4-DNA web site is available at http://geant4-dna.org
34 //
35 //
36 /// \file DetectorConstruction.cc
37 /// \brief Implementation of the DetectorConstruction class
38
39 #include "DetectorConstruction.hh"
40
41 #include "DetectorMessenger.hh"
42
43 #include "G4Box.hh"
44 #include "G4Colour.hh"
45 #include "G4GeometryManager.hh"
46 #include "G4LogicalVolume.hh"
47 #include "G4LogicalVolumeStore.hh"
48 #include "G4Material.hh"
49 #include "G4NistManager.hh"
50 #include "G4Orb.hh"
51 #include "G4PVPlacement.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4PhysicalVolumeStore.hh"
54 #include "G4SolidStore.hh"
55 #include "G4SystemOfUnits.hh"
56 #include "G4Tubs.hh"
57 #include "G4VisAttributes.hh"
58
59 #ifdef G4MULTITHREADED
60 # include "G4MTRunManager.hh"
61 #else
62 # include "G4RunManager.hh"
63 #endif
64
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
66
67 DetectorConstruction::DetectorConstruction()
68 : G4VUserDetectorConstruction(), fpMessenger(0), fCheckOverlaps(false)
69 {
70 // Select a PDB file name by default
71 // otherwise modified by the LoadPDBfile messenger
72 //
73 fPdbFileName = G4String("1ZBB.pdb");
74 fPdbFileStatus = 0;
75 fChosenOption = 11;
76
77 fpDefaultMaterial = 0;
78 fpWaterMaterial = 0;
79 fpMessenger = new DetectorMessenger(this);
80 }
81
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83
84 DetectorConstruction::~DetectorConstruction() {}
85
86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
87
88 G4VPhysicalVolume* DetectorConstruction::Construct()
89 {
90 fChosenOption = 11; // Draw atomic with bounding box (visualisation by default)
91 fPdbFileStatus = 0; // There is no PDB file loaded at this stage
92
93 // Define materials and geometry
94 G4VPhysicalVolume* worldPV;
95 worldPV = DefineVolumes(fPdbFileName, fChosenOption);
96 return worldPV;
97 }
98
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
100
101 void DetectorConstruction::ConstructMaterials()
102 {
103 //[G4_WATER]
104 //
105 G4NistManager* nistManager = G4NistManager::Instance();
106 G4bool fromIsotopes = false;
107 nistManager->FindOrBuildMaterial("G4_WATER", fromIsotopes);
108 fpWaterMaterial = G4Material::GetMaterial("G4_WATER");
109
110 //[Vacuum]
111 G4double a; // mass of a mole;
112 G4double z; // z=mean number of protons;
113 G4double density;
114 new G4Material("Galactic", z = 1., a = 1.01 * g / mole, density = universe_mean_density,
115 kStateGas, 2.73 * kelvin, 3.e-18 * pascal);
116 fpDefaultMaterial = G4Material::GetMaterial("Galactic");
117 }
118
119 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
120
121 void DetectorConstruction::CheckMaterials()
122 {
123 if (!fpDefaultMaterial)
124 G4Exception("DetectorConstruction::CheckMaterials", "DEFAULT_MATERIAL_NOT_INIT_1",
125 FatalException, "Default material not initialized.");
126
127 if (!fpWaterMaterial)
128 G4Exception("DetectorConstruction::CheckMaterials", "WATER_MATERIAL_NOT_INIT", FatalException,
129 "Water material not initialized.");
130 }
131
132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
133
134 G4VPhysicalVolume* DetectorConstruction::ConstructWorld()
135 {
136 // Geometry parameters
137 G4double worldSize = 1000 * 1 * angstrom;
138
139 if (!fpDefaultMaterial) {
140 G4Exception("DetectorConstruction::ConstructWorld", "DEFAULT_MATERIAL_NOT_INIT_2",
141 FatalException, "Default material not initialized.");
142 }
143
144 //
145 // World
146 //
147 G4VSolid* worldS = new G4Box("World", // its name
148 worldSize / 2, worldSize / 2, worldSize / 2); // its size
149
150 G4LogicalVolume* worldLV = new G4LogicalVolume(worldS, // its solid
151 fpDefaultMaterial, // its material
152 "World"); // its name
153
154 G4VisAttributes* MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Gray()));
155 MyVisAtt_ZZ->SetVisibility(false);
156 worldLV->SetVisAttributes(MyVisAtt_ZZ);
157
158 G4VPhysicalVolume* worldPV = new G4PVPlacement(0, // no rotation
159 G4ThreeVector(), // at (0,0,0)
160 worldLV, // its logical volume
161 "World", // its name
162 0, // its mother volume
163 false, // no boolean operation
164 0, // copy number
165 true); // checking overlaps forced to YES
166
167 return worldPV;
168 }
169
170 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
171
172 G4VPhysicalVolume* DetectorConstruction::DefineVolumes(G4String filename, unsigned short int option)
173 {
174 // Clean old geometry
175 //
176 G4GeometryManager::GetInstance()->OpenGeometry();
177 G4PhysicalVolumeStore::GetInstance()->Clean();
178 G4LogicalVolumeStore::GetInstance()->Clean();
179 G4SolidStore::GetInstance()->Clean();
180
181 // Define Materials
182 //
183 ConstructMaterials();
184
185 // Reconstruct world not to superimpose on geometries
186 G4VPhysicalVolume* worldPV;
187 G4LogicalVolume* worldLV;
188 worldPV = ConstructWorld();
189 worldLV = worldPV->GetLogicalVolume();
190
191 // List of molecules inside pdb file separated by TER keyword
192 fpMoleculeList = NULL;
193 fpBarycenterList = NULL;
194
195 //'fPDBlib.load' is currently done for each 'DefineVolumes' call.
196 int verbosity = 0; // To be implemented
197 unsigned short int isDNA;
198 fpMoleculeList = fPDBlib.Load(filename, isDNA, verbosity);
199 if (fpMoleculeList != NULL && isDNA == 1) {
200 fPDBlib.ComputeNbNucleotidsPerStrand(fpMoleculeList);
201 fpBarycenterList = fPDBlib.ComputeNucleotideBarycenters(fpMoleculeList);
202 G4cout << "This PDB file is DNA." << G4endl;
203 }
204
205 if (fpMoleculeList != NULL) {
206 fPdbFileStatus = 1;
207 }
208
209 if (option == 1) {
210 AtomisticView(worldLV, fpMoleculeList, 1.0);
211 }
212 else if (option == 2) {
213 BarycenterView(worldLV, fpBarycenterList);
214 }
215 else if (option == 3) {
216 ResiduesView(worldLV, fpBarycenterList);
217 }
218 else if (option == 10) {
219 DrawBoundingVolume(worldLV, fpMoleculeList);
220 }
221 else if (option == 11) {
222 AtomisticView(worldLV, fpMoleculeList, 1.0);
223 DrawBoundingVolume(worldLV, fpMoleculeList);
224 }
225 else if (option == 12) {
226 BarycenterView(worldLV, fpBarycenterList);
227 DrawBoundingVolume(worldLV, fpMoleculeList);
228 }
229 else if (option == 13) {
230 ResiduesView(worldLV, fpBarycenterList);
231 DrawBoundingVolume(worldLV, fpMoleculeList);
232 }
233 // Always return the physical World
234 //
235 return worldPV;
236 }
237
238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
239
240 PDBlib DetectorConstruction::GetPDBlib()
241 {
242 return fPDBlib;
243 }
244
245 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
246
247 Barycenter* DetectorConstruction::GetBarycenterList()
248 {
249 return fpBarycenterList;
250 }
251
252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
253
254 Molecule* DetectorConstruction::GetMoleculeList()
255 {
256 return fpMoleculeList;
257 }
258
259 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
260
261 //////////////////////////////////////////////////
262 ///////////////// BEGIN atomistic representation
263 //
264 void DetectorConstruction::AtomisticView(G4LogicalVolume* worldLV, Molecule* moleculeListTemp,
265 double atomSizeFactor)
266 {
267 CheckMaterials();
268
269 // All solids are defined for different color attributes :
270 G4double sphereSize = atomSizeFactor * 1 * angstrom;
271 G4VSolid* atomS_H = new G4Orb("Sphere", sphereSize * 1.2);
272 G4VSolid* atomS_C = new G4Orb("Sphere", sphereSize * 1.7);
273 G4VSolid* atomS_O = new G4Orb("Sphere", sphereSize * 1.52);
274 G4VSolid* atomS_N = new G4Orb("Sphere", sphereSize * 1.55);
275 G4VSolid* atomS_S = new G4Orb("Sphere", sphereSize * 1.8);
276 G4VSolid* atomS_P = new G4Orb("Sphere", sphereSize * 1.8);
277 G4VSolid* atomS_X = new G4Orb("Sphere", sphereSize); // Undefined
278
279 // Logical volumes and color table for atoms :
280 G4LogicalVolume* atomLV_H =
281 new G4LogicalVolume(atomS_H, fpWaterMaterial, "atomLV_H"); // its solid, material, name
282 G4VisAttributes* MyVisAtt_H = new G4VisAttributes(G4Colour(G4Colour::White()));
283 MyVisAtt_H->SetForceSolid(true);
284 atomLV_H->SetVisAttributes(MyVisAtt_H);
285
286 G4LogicalVolume* atomLV_C =
287 new G4LogicalVolume(atomS_C, fpWaterMaterial, "atomLV_C"); // its solid, material, name
288 G4VisAttributes* MyVisAtt_C =
289 new G4VisAttributes(G4Colour(G4Colour::Gray())); // Black in CPK, but bad rendering
290 MyVisAtt_C->SetForceSolid(true);
291 atomLV_C->SetVisAttributes(MyVisAtt_C);
292
293 G4LogicalVolume* atomLV_O =
294 new G4LogicalVolume(atomS_O, fpWaterMaterial, "atomLV_O"); // its solid, material, name
295 G4VisAttributes* MyVisAtt_O = new G4VisAttributes(G4Colour(G4Colour::Red()));
296 MyVisAtt_O->SetForceSolid(true);
297 atomLV_O->SetVisAttributes(MyVisAtt_O);
298
299 G4LogicalVolume* atomLV_N =
300 new G4LogicalVolume(atomS_N, fpWaterMaterial, "atomLV_N"); // its solid, material, name
301 G4VisAttributes* MyVisAtt_N = new G4VisAttributes(G4Colour(G4Colour(0., 0., 0.5))); // Dark blue
302 MyVisAtt_N->SetForceSolid(true);
303 atomLV_N->SetVisAttributes(MyVisAtt_N);
304
305 G4LogicalVolume* atomLV_S =
306 new G4LogicalVolume(atomS_S, fpWaterMaterial, "atomLV_S"); // its solid, material, name
307 G4VisAttributes* MyVisAtt_S = new G4VisAttributes(G4Colour(G4Colour::Yellow()));
308 MyVisAtt_S->SetForceSolid(true);
309 atomLV_S->SetVisAttributes(MyVisAtt_S);
310
311 G4LogicalVolume* atomLV_P =
312 new G4LogicalVolume(atomS_P, fpWaterMaterial, "atomLV_P"); // its solid, material, name
313 G4VisAttributes* MyVisAtt_P = new G4VisAttributes(G4Colour(G4Colour(1.0, 0.5, 0.))); // Orange
314 MyVisAtt_P->SetForceSolid(true);
315 atomLV_P->SetVisAttributes(MyVisAtt_P);
316
317 G4LogicalVolume* atomLV_X =
318 new G4LogicalVolume(atomS_X, fpWaterMaterial, "atomLV_X"); // its solid, material, name
319 G4VisAttributes* MyVisAtt_X =
320 new G4VisAttributes(G4Colour(G4Colour(1.0, 0.75, 0.8))); // Pink, other elements in CKP
321 MyVisAtt_X->SetForceSolid(true);
322 atomLV_X->SetVisAttributes(MyVisAtt_X);
323
324 // Placement and physical volume construction from memory
325 Residue* residueListTemp;
326 Atom* AtomTemp;
327
328 int nbAtomTot = 0; // Number total of atoms
329 int nbAtomH = 0, nbAtomC = 0, nbAtomO = 0, nbAtomN = 0, nbAtomS = 0, nbAtomP = 0;
330 int nbAtomX = 0;
331 int k = 0;
332
333 while (moleculeListTemp) {
334 residueListTemp = moleculeListTemp->GetFirst();
335
336 k++;
337 while (residueListTemp) {
338 AtomTemp = residueListTemp->GetFirst();
339
340 int startFrom = 0;
341 int upTo = residueListTemp->fNbAtom; // case Base+sugar+phosphat (all atoms)
342 for (int i = 0; i < (upTo + startFrom); i++) // Phosphat or Sugar or Base
343 {
344 if (AtomTemp->fElement.compare("H") == 0) {
345 nbAtomH++;
346 new G4PVPlacement(0,
347 G4ThreeVector(AtomTemp->fX * 1 * angstrom, AtomTemp->fY * 1 * angstrom,
348 AtomTemp->fZ * 1 * angstrom),
349 atomLV_H, "atomP", worldLV, false, 0, fCheckOverlaps);
350 }
351 else if (AtomTemp->fElement.compare("C") == 0) {
352 nbAtomC++;
353 new G4PVPlacement(0,
354 G4ThreeVector(AtomTemp->fX * 1 * angstrom, AtomTemp->fY * 1 * angstrom,
355 AtomTemp->fZ * 1 * angstrom),
356 atomLV_C, "atomP", worldLV, false, 0, fCheckOverlaps);
357 }
358 else if (AtomTemp->fElement.compare("O") == 0) {
359 nbAtomO++;
360 new G4PVPlacement(0,
361 G4ThreeVector(AtomTemp->fX * 1 * angstrom, AtomTemp->fY * 1 * angstrom,
362 AtomTemp->fZ * 1 * angstrom),
363 atomLV_O, "atomP", worldLV, false, 0, fCheckOverlaps);
364 }
365 else if (AtomTemp->fElement.compare("N") == 0) {
366 nbAtomN++;
367 new G4PVPlacement(0,
368 G4ThreeVector(AtomTemp->fX * 1 * angstrom, AtomTemp->fY * 1 * angstrom,
369 AtomTemp->fZ * 1 * angstrom),
370 atomLV_N, "atomP", worldLV, false, 0, fCheckOverlaps);
371 }
372 else if (AtomTemp->fElement.compare("S") == 0) {
373 nbAtomS++;
374 new G4PVPlacement(0,
375 G4ThreeVector(AtomTemp->fX * 1 * angstrom, AtomTemp->fY * 1 * angstrom,
376 AtomTemp->fZ * 1 * angstrom),
377 atomLV_S, "atomP", worldLV, false, 0, fCheckOverlaps);
378 }
379 else if (AtomTemp->fElement.compare("P") == 0) {
380 nbAtomP++;
381 new G4PVPlacement(0,
382 G4ThreeVector(AtomTemp->fX * 1 * angstrom, AtomTemp->fY * 1 * angstrom,
383 AtomTemp->fZ * 1 * angstrom),
384 atomLV_P, "atomP", worldLV, false, 0, fCheckOverlaps);
385 }
386 else {
387 nbAtomX++;
388 new G4PVPlacement(0,
389 G4ThreeVector(AtomTemp->fX * 1 * angstrom, AtomTemp->fY * 1 * angstrom,
390 AtomTemp->fZ * 1 * angstrom),
391 atomLV_X, "atomP", worldLV, false, 0, fCheckOverlaps);
392 }
393
394 nbAtomTot++;
395 //}//End if (i>=startFrom)
396 AtomTemp = AtomTemp->GetNext();
397 } // end of for ( i=0 ; i < residueListTemp->nbAtom ; i++)
398
399 residueListTemp = residueListTemp->GetNext();
400 } // end of while (residueListTemp)
401
402 moleculeListTemp = moleculeListTemp->GetNext();
403 } // end of while (moleculeListTemp)
404
405 G4cout << "**************** atomisticView(...) ****************" << G4endl;
406 G4cout << "Number of loaded chains = " << k << G4endl;
407 G4cout << "Number of Atoms = " << nbAtomTot << G4endl;
408 G4cout << "Number of Hydrogens = " << nbAtomH << G4endl;
409 G4cout << "Number of Carbons = " << nbAtomC << G4endl;
410 G4cout << "Number of Oxygens = " << nbAtomO << G4endl;
411 G4cout << "Number of Nitrogens = " << nbAtomN << G4endl;
412 G4cout << "Number of Sulfurs = " << nbAtomS << G4endl;
413 G4cout << "Number of Phosphorus = " << nbAtomP << G4endl;
414 G4cout << "Number of undifined atoms =" << nbAtomX << G4endl << G4endl;
415 }
416
417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
418
419 //////////////////////////////////////////////////
420 ///////////////// BEGIN representation for nucleotide barycenters
421 //
422 void DetectorConstruction::BarycenterView(G4LogicalVolume* worldLV, Barycenter* barycenterListTemp)
423 {
424 CheckMaterials();
425
426 G4VSolid* atomS_ZZ;
427 G4LogicalVolume* atomLV_ZZ;
428 G4VisAttributes* MyVisAtt_ZZ;
429
430 while (barycenterListTemp) {
431 atomS_ZZ = new G4Orb("Sphere", (barycenterListTemp->GetRadius()) * angstrom);
432 atomLV_ZZ = new G4LogicalVolume(atomS_ZZ, fpWaterMaterial, "atomLV_ZZ");
433 MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Magenta()));
434 MyVisAtt_ZZ->SetForceSolid(true);
435 atomLV_ZZ->SetVisAttributes(MyVisAtt_ZZ);
436
437 new G4PVPlacement(0,
438 G4ThreeVector(barycenterListTemp->fCenterX * 1 * angstrom,
439 barycenterListTemp->fCenterY * 1 * angstrom,
440 barycenterListTemp->fCenterZ * 1 * angstrom),
441 atomLV_ZZ, "atomZZ", worldLV, false, 0, fCheckOverlaps);
442
443 barycenterListTemp = barycenterListTemp->GetNext();
444 } // end of while (moleculeListTemp)
445 }
446
447 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
448
449 //////////////////////////////////////////////////
450 ///////////////// BEGIN representation for Base Sugar Phosphate barycenters
451 //
452 void DetectorConstruction::ResiduesView(G4LogicalVolume* worldLV, Barycenter* barycenterListTemp)
453 {
454 CheckMaterials();
455 G4VisAttributes* MyVisAtt_ZZ;
456
457 G4VSolid* tubS1_ZZ;
458 G4LogicalVolume* tubLV1_ZZ;
459 G4VSolid* tubS2_ZZ;
460 G4LogicalVolume* tubLV2_ZZ;
461
462 G4VSolid* AS_ZZ;
463 G4LogicalVolume* ALV_ZZ;
464 G4VSolid* BS_ZZ;
465 G4LogicalVolume* BLV_ZZ;
466 G4VSolid* CS_ZZ;
467 G4LogicalVolume* CLV_ZZ;
468
469 while (barycenterListTemp) {
470 // 3 spheres to Base, Sugar, Phosphate)
471 AS_ZZ = new G4Orb("Sphere", 1. * angstrom);
472 ALV_ZZ = new G4LogicalVolume(AS_ZZ, fpWaterMaterial, "ALV_ZZ");
473 MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Blue()));
474 MyVisAtt_ZZ->SetForceSolid(true);
475 ALV_ZZ->SetVisAttributes(MyVisAtt_ZZ);
476 new G4PVPlacement(0,
477 G4ThreeVector(barycenterListTemp->fCenterBaseX * angstrom,
478 barycenterListTemp->fCenterBaseY * angstrom,
479 barycenterListTemp->fCenterBaseZ * angstrom),
480 ALV_ZZ, "AZZ", worldLV, false, 0, fCheckOverlaps);
481
482 BS_ZZ = new G4Orb("Sphere", 1. * angstrom);
483 BLV_ZZ = new G4LogicalVolume(BS_ZZ, fpWaterMaterial, "BLV_ZZ");
484 MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Red()));
485 MyVisAtt_ZZ->SetForceSolid(true);
486 BLV_ZZ->SetVisAttributes(MyVisAtt_ZZ);
487 new G4PVPlacement(0,
488 G4ThreeVector((barycenterListTemp->fCenterPhosphateX) * angstrom,
489 (barycenterListTemp->fCenterPhosphateY) * angstrom,
490 (barycenterListTemp->fCenterPhosphateZ) * angstrom),
491 BLV_ZZ, "BZZ", worldLV, false, 0, fCheckOverlaps);
492
493 CS_ZZ = new G4Orb("Sphere", 1. * angstrom);
494 CLV_ZZ = new G4LogicalVolume(CS_ZZ, fpWaterMaterial, "CLV_ZZ");
495 MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Yellow()));
496 MyVisAtt_ZZ->SetForceSolid(true);
497 CLV_ZZ->SetVisAttributes(MyVisAtt_ZZ);
498 new G4PVPlacement(0,
499 G4ThreeVector(barycenterListTemp->fCenterSugarX * angstrom,
500 barycenterListTemp->fCenterSugarY * angstrom,
501 barycenterListTemp->fCenterSugarZ * angstrom),
502 CLV_ZZ, "CZZ", worldLV, false, 0, fCheckOverlaps);
503
504 // 2 cylinders (Base<=>Sugar and Sugar<=>Phosphat)
505 // Cylinder Base<=>Sugar
506 tubS1_ZZ = new G4Tubs(
507 "Cylinder", 0., 0.5 * angstrom,
508 std::sqrt((barycenterListTemp->fCenterBaseX - barycenterListTemp->fCenterSugarX)
509 * (barycenterListTemp->fCenterBaseX - barycenterListTemp->fCenterSugarX)
510 + (barycenterListTemp->fCenterBaseY - barycenterListTemp->fCenterSugarY)
511 * (barycenterListTemp->fCenterBaseY - barycenterListTemp->fCenterSugarY)
512 + (barycenterListTemp->fCenterBaseZ - barycenterListTemp->fCenterSugarZ)
513 * (barycenterListTemp->fCenterBaseZ - barycenterListTemp->fCenterSugarZ))
514 / 2 * angstrom,
515 0., 2. * pi);
516
517 tubLV1_ZZ = new G4LogicalVolume(tubS1_ZZ, fpWaterMaterial, "tubLV_ZZ");
518 MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Green()));
519 MyVisAtt_ZZ->SetForceSolid(true);
520 tubLV1_ZZ->SetVisAttributes(MyVisAtt_ZZ);
521
522 G4double Ux = barycenterListTemp->fCenterBaseX - barycenterListTemp->fCenterSugarX;
523 G4double Uy = barycenterListTemp->fCenterBaseY - barycenterListTemp->fCenterSugarY;
524 G4double Uz = barycenterListTemp->fCenterBaseZ - barycenterListTemp->fCenterSugarZ;
525 G4double llUll = std::sqrt(Ux * Ux + Uy * Uy + Uz * Uz);
526
527 Ux = Ux / llUll;
528 Uy = Uy / llUll;
529 Uz = Uz / llUll;
530
531 G4ThreeVector direction = G4ThreeVector(Ux, Uy, Uz);
532 G4double theta_euler = direction.theta();
533 G4double phi_euler = direction.phi();
534 G4double psi_euler = 0;
535
536 // Warning : clhep Euler constructor build inverse matrix !
537 G4RotationMatrix rotm1Inv = G4RotationMatrix(phi_euler + pi / 2, theta_euler, psi_euler);
538 G4RotationMatrix rotm1 = rotm1Inv.inverse();
539 G4ThreeVector translm1 = G4ThreeVector(
540 (barycenterListTemp->fCenterBaseX + barycenterListTemp->fCenterSugarX) / 2. * angstrom,
541 (barycenterListTemp->fCenterBaseY + barycenterListTemp->fCenterSugarY) / 2. * angstrom,
542 (barycenterListTemp->fCenterBaseZ + barycenterListTemp->fCenterSugarZ) / 2. * angstrom);
543 G4Transform3D transform1 = G4Transform3D(rotm1, translm1);
544 new G4PVPlacement(transform1, // rotation translation
545 tubLV1_ZZ, "atomZZ", worldLV, false, 0, fCheckOverlaps);
546
547 // Cylinder Sugar<=>Phosphat
548 tubS2_ZZ = new G4Tubs(
549 "Cylinder2", 0., 0.5 * angstrom,
550 std::sqrt((barycenterListTemp->fCenterSugarX - barycenterListTemp->fCenterPhosphateX)
551 * (barycenterListTemp->fCenterSugarX - barycenterListTemp->fCenterPhosphateX)
552 + (barycenterListTemp->fCenterSugarY - barycenterListTemp->fCenterPhosphateY)
553 * (barycenterListTemp->fCenterSugarY - barycenterListTemp->fCenterPhosphateY)
554 + (barycenterListTemp->fCenterSugarZ - barycenterListTemp->fCenterPhosphateZ)
555 * (barycenterListTemp->fCenterSugarZ - barycenterListTemp->fCenterPhosphateZ))
556 / 2 * angstrom,
557 0., 2. * pi);
558
559 tubLV2_ZZ = new G4LogicalVolume(tubS2_ZZ, fpWaterMaterial, "tubLV2_ZZ");
560 MyVisAtt_ZZ = new G4VisAttributes(G4Colour(1, 0.5, 0));
561 MyVisAtt_ZZ->SetForceSolid(true);
562 tubLV2_ZZ->SetVisAttributes(MyVisAtt_ZZ);
563
564 Ux = barycenterListTemp->fCenterSugarX - barycenterListTemp->fCenterPhosphateX;
565 Uy = barycenterListTemp->fCenterSugarY - barycenterListTemp->fCenterPhosphateY;
566 Uz = barycenterListTemp->fCenterSugarZ - barycenterListTemp->fCenterPhosphateZ;
567 llUll = std::sqrt(Ux * Ux + Uy * Uy + Uz * Uz);
568
569 Ux = Ux / llUll;
570 Uy = Uy / llUll;
571 Uz = Uz / llUll;
572
573 direction = G4ThreeVector(Ux, Uy, Uz);
574 theta_euler = direction.theta();
575 phi_euler = direction.phi();
576 psi_euler = 0;
577
578 // Warning : clhep Euler constructor build inverse matrix !
579 rotm1Inv = G4RotationMatrix(phi_euler + pi / 2, theta_euler, psi_euler);
580 rotm1 = rotm1Inv.inverse();
581 translm1 = G4ThreeVector(
582 (barycenterListTemp->fCenterSugarX + barycenterListTemp->fCenterPhosphateX) / 2. * angstrom,
583 (barycenterListTemp->fCenterSugarY + barycenterListTemp->fCenterPhosphateY) / 2. * angstrom,
584 (barycenterListTemp->fCenterSugarZ + barycenterListTemp->fCenterPhosphateZ) / 2. * angstrom);
585 transform1 = G4Transform3D(rotm1, translm1);
586 new G4PVPlacement(transform1, // rotation translation
587 tubLV2_ZZ, "atomZZ", worldLV, false, 0, fCheckOverlaps);
588
589 barycenterListTemp = barycenterListTemp->GetNext();
590 } // end of while sur moleculeListTemp
591 }
592
593 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
594
595 //////////////////////////////////////////////////
596 ///////////////// BEGIN draw a bounding volume
597 //
598 void DetectorConstruction::DrawBoundingVolume(G4LogicalVolume* worldLV, Molecule* moleculeListTemp)
599 {
600 CheckMaterials();
601
602 double dX, dY, dZ; // Dimensions for bounding volume
603 double tX, tY, tZ; // Translation for bounding volume
604 fPDBlib.ComputeBoundingVolumeParams(moleculeListTemp, dX, dY, dZ, tX, tY, tZ);
605
606 //[Geometry] Build a box :
607 G4VSolid* boundingS =
608 new G4Box("Bounding", dX * 1 * angstrom, dY * 1 * angstrom, dZ * 1 * angstrom);
609
610 G4LogicalVolume* boundingLV = new G4LogicalVolume(boundingS, fpWaterMaterial, "BoundingLV");
611
612 G4RotationMatrix* pRot = new G4RotationMatrix();
613
614 G4VisAttributes* MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Gray()));
615 boundingLV->SetVisAttributes(MyVisAtt_ZZ);
616
617 new G4PVPlacement(pRot,
618 G4ThreeVector(tX * 1 * angstrom, tY * 1 * angstrom,
619 tZ * 1 * angstrom), // at (x,y,z)
620 boundingLV, "boundingPV", worldLV, false, 0, fCheckOverlaps);
621 }
622
623 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
624
625 void DetectorConstruction::LoadPDBfile(G4String fileName)
626 {
627 G4cout << "Load PDB file : " << fileName << "." << G4endl << G4endl;
628 fPdbFileName = fileName;
629 #ifdef G4MULTITHREADED
630 G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, fChosenOption));
631 G4MTRunManager::GetRunManager()->ReinitializeGeometry();
632 #else
633 G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, fChosenOption));
634 #endif
635 }
636
637 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
638
639 void DetectorConstruction::BuildBoundingVolume()
640 {
641 if (fPdbFileStatus > 0) // a PDB file has been loaded
642 {
643 G4cout << "Build only world volume and bounding volume"
644 " for computation."
645 << G4endl << G4endl;
646 #ifdef G4MULTITHREADED
647 G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 10));
648 G4MTRunManager::GetRunManager()->ReinitializeGeometry();
649 #else
650 G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 10));
651 #endif
652 }
653 else
654 G4cout << "PDB file not found!" << G4endl << G4endl;
655 }
656
657 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
658
659 void DetectorConstruction::DrawAtoms_()
660 {
661 if (fPdbFileStatus > 0) // a PDB file has been loaded
662 {
663 #ifdef G4MULTITHREADED
664 G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 1));
665 G4MTRunManager::GetRunManager()->ReinitializeGeometry();
666 #else
667 G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 1));
668 #endif
669 }
670 else
671 G4cout << "PDB file not found!" << G4endl << G4endl;
672 }
673
674 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
675
676 void DetectorConstruction::DrawNucleotides_()
677 {
678 if (fPdbFileStatus > 0) // a PDB file has been loaded
679 {
680 #ifdef G4MULTITHREADED
681 G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 2));
682 G4MTRunManager::GetRunManager()->ReinitializeGeometry();
683 #else
684 G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 2));
685 #endif
686 }
687 else
688 G4cout << "PDB file not found!" << G4endl << G4endl;
689 }
690
691 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
692
693 void DetectorConstruction::DrawResidues_()
694 {
695 if (fPdbFileStatus > 0) // a PDB file has been loaded
696 {
697 #ifdef G4MULTITHREADED
698 G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 3));
699 G4MTRunManager::GetRunManager()->ReinitializeGeometry();
700 #else
701 G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 3));
702 #endif
703 }
704 else
705 G4cout << "PDB file not found!" << G4endl << G4endl;
706 }
707
708 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
709
710 void DetectorConstruction::DrawAtomsWithBounding_()
711 {
712 if (fPdbFileStatus > 0) // a PDB file has been loaded
713 {
714 #ifdef G4MULTITHREADED
715 G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 11));
716 G4MTRunManager::GetRunManager()->ReinitializeGeometry();
717 #else
718 G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 11));
719 #endif
720 }
721 else
722 G4cout << "PDB file not found!" << G4endl << G4endl;
723 }
724
725 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
726
727 void DetectorConstruction::DrawNucleotidesWithBounding_()
728 {
729 if (fPdbFileStatus > 0) // a PDB file has been loaded
730 {
731 #ifdef G4MULTITHREADED
732 G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 12));
733 G4MTRunManager::GetRunManager()->ReinitializeGeometry();
734 #else
735 G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 12));
736 #endif
737 }
738 else
739 G4cout << "PDB file not found!" << G4endl << G4endl;
740 }
741
742 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
743
744 void DetectorConstruction::DrawResiduesWithBounding_()
745 {
746 if (fPdbFileStatus > 0) // a PDB file has been loaded
747 {
748 #ifdef G4MULTITHREADED
749 G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 13));
750 G4MTRunManager::GetRunManager()->ReinitializeGeometry();
751 #else
752 G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 13));
753 #endif
754 }
755 else
756 G4cout << "PDB file not found!" << G4endl << G4endl;
757 }
758