Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 // 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 // $Id$
36 //
39 
40 #include "DetectorConstruction.hh"
41 
42 #include "DetectorMessenger.hh"
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 "G4PhysicalConstants.hh"
52 #include "G4PhysicalVolumeStore.hh"
53 #include "G4PVPlacement.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(),
69  fpMessenger(0),
70  fCheckOverlaps(false)
71 {
72  //Select a PDB file name by default
73  //otherwise modified by the LoadPDBfile messenger
74  //
75  fPdbFileName=G4String("1ZBB.pdb");
76  fPdbFileStatus=0;
77  fChosenOption=11;
78 
79  fpDefaultMaterial=0;
80  fpWaterMaterial=0;
81  fpMessenger = new DetectorMessenger(this);
82 }
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
85 
86 DetectorConstruction::~DetectorConstruction()
87 {
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91 
92 G4VPhysicalVolume* DetectorConstruction::Construct()
93 {
94  fChosenOption=11;//Draw atomic with bounding box (visualisation by default)
95  fPdbFileStatus=0;//There is no PDB file loaded at this stage
96 
97  //Define materials and geometry
98  G4VPhysicalVolume* worldPV;
99  worldPV=DefineVolumes(fPdbFileName,fChosenOption);
100  return worldPV;
101 }
102 
103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
104 
105 void DetectorConstruction::ConstructMaterials()
106 {
107  //[G4_WATER]
108  //
109  G4NistManager* nistManager = G4NistManager::Instance();
110  G4bool fromIsotopes = false;
111  nistManager->FindOrBuildMaterial("G4_WATER", fromIsotopes);
112  fpWaterMaterial = G4Material::GetMaterial("G4_WATER");
113 
114  //[Vacuum]
115  G4double a; // mass of a mole;
116  G4double z; // z=mean number of protons;
117  G4double density;
118  new G4Material("Galactic", z=1., a=1.01*g/mole,density= universe_mean_density,
119  kStateGas, 2.73*kelvin, 3.e-18*pascal);
120  fpDefaultMaterial = G4Material::GetMaterial("Galactic");
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
124 
125 void DetectorConstruction::CheckMaterials()
126 {
127  if(!fpDefaultMaterial)
128  G4Exception("DetectorConstruction::CheckMaterials",
129  "DEFAULT_MATERIAL_NOT_INIT_1",
130  FatalException,
131  "Default material not initialized.");
132 
133  if(!fpWaterMaterial)
134  G4Exception("DetectorConstruction::CheckMaterials",
135  "WATER_MATERIAL_NOT_INIT",
136  FatalException,
137  "Water material not initialized.");
138 }
139 
140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
141 
142 G4VPhysicalVolume* DetectorConstruction::ConstructWorld()
143 {
144  // Geometry parameters
145  G4double worldSize = 1000*1*angstrom;
146 
147  if ( !fpDefaultMaterial )
148  {
149  G4Exception("DetectorConstruction::ConstructWorld",
150  "DEFAULT_MATERIAL_NOT_INIT_2",
151  FatalException,
152  "Default material not initialized.");
153  }
154 
155  //
156  // World
157  //
158  G4VSolid* worldS
159  = new G4Box("World", // its name
160  worldSize/2, worldSize/2, worldSize/2); // its size
161 
162  G4LogicalVolume*
163  worldLV
164  = new G4LogicalVolume(
165  worldS, // its solid
166  fpDefaultMaterial, // its material
167  "World"); // its name
168 
169  G4VisAttributes *MyVisAtt_ZZ = new G4VisAttributes(
170  G4Colour(G4Colour::Gray()));
171  MyVisAtt_ZZ ->SetVisibility (false);
172  worldLV->SetVisAttributes(MyVisAtt_ZZ);
173 
174  G4VPhysicalVolume*
175  worldPV
176  = new G4PVPlacement(
177  0, // no rotation
178  G4ThreeVector(), // at (0,0,0)
179  worldLV, // its logical volume
180  "World", // its name
181  0, // its mother volume
182  false, // no boolean operation
183  0, // copy number
184  true); // checking overlaps forced to YES
185 
186  return worldPV;
187 }
188 
189 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
190 
191 G4VPhysicalVolume* DetectorConstruction::DefineVolumes(G4String filename,
192  unsigned short int option)
193 {
194  //Clean old geometry
195  //
196  G4GeometryManager::GetInstance()->OpenGeometry();
197  G4PhysicalVolumeStore::GetInstance()->Clean();
198  G4LogicalVolumeStore::GetInstance()->Clean();
199  G4SolidStore::GetInstance()->Clean();
200 
201  // Define Materials
202  //
203  ConstructMaterials();
204 
205  //Reconstruct world not to superimpose on geometries
206  G4VPhysicalVolume* worldPV;
207  G4LogicalVolume* worldLV;
208  worldPV=ConstructWorld();
209  worldLV=worldPV->GetLogicalVolume();
210 
211  //List of molecules inside pdb file separated by TER keyword
212  fpMoleculeList=NULL;
213  fpBarycenterList=NULL;
214 
215  //'fPDBlib.load' is currently done for each 'DefineVolumes' call.
216  int verbosity=0; //To be implemented
217  unsigned short int isDNA;
218  fpMoleculeList = fPDBlib.Load(filename,isDNA,verbosity);
219  if (fpMoleculeList!=NULL && isDNA==1)
220  {
221  fPDBlib.ComputeNbNucleotidsPerStrand(fpMoleculeList );
222  fpBarycenterList=fPDBlib.ComputeNucleotideBarycenters(fpMoleculeList );
223  G4cout<<"This PDB file is DNA."<<G4endl;
224  }
225 
226  if (fpMoleculeList!=NULL)
227  {
228  fPdbFileStatus=1;
229  }
230 
231  if (option==1)
232  {
233  AtomisticView( worldLV,fpMoleculeList,1.0);
234  }
235  else
236  if (option==2)
237  {
238  BarycenterView(worldLV,fpBarycenterList);
239  }
240  else
241  if (option==3)
242  {
243  ResiduesView(worldLV,fpBarycenterList);
244  }
245  else
246  if (option==10)
247  {
248  DrawBoundingVolume( worldLV,fpMoleculeList);
249  }
250  else
251  if (option==11)
252  {
253  AtomisticView( worldLV,fpMoleculeList,1.0);
254  DrawBoundingVolume( worldLV,fpMoleculeList);
255  }
256  else
257  if (option==12)
258  {
259  BarycenterView(worldLV,fpBarycenterList);
260  DrawBoundingVolume( worldLV,fpMoleculeList);
261  }
262  else
263  if (option==13)
264  {
265  ResiduesView(worldLV,fpBarycenterList);
266  DrawBoundingVolume( worldLV,fpMoleculeList);
267  }
268  // Always return the physical World
269  //
270  return worldPV;
271 }
272 
273 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
274 
275 PDBlib DetectorConstruction::GetPDBlib()
276 {
277  return fPDBlib;
278 }
279 
280 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
281 
282 Barycenter *DetectorConstruction::GetBarycenterList()
283 {
284  return fpBarycenterList;
285 }
286 
287 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
288 
289 Molecule *DetectorConstruction::GetMoleculeList()
290 {
291  return fpMoleculeList;
292 }
293 
294 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
295 
298 //
299 void DetectorConstruction::AtomisticView(G4LogicalVolume* worldLV,
300  Molecule *moleculeListTemp, double atomSizeFactor)
301 {
302  CheckMaterials();
303 
304  // All solids are defined for different color attributes :
305  G4double sphereSize = atomSizeFactor*1*angstrom;
306  G4VSolid* atomS_H = new G4Orb("Sphere", sphereSize*1.2);
307  G4VSolid* atomS_C = new G4Orb("Sphere", sphereSize*1.7);
308  G4VSolid* atomS_O = new G4Orb("Sphere", sphereSize*1.52);
309  G4VSolid* atomS_N = new G4Orb("Sphere", sphereSize*1.55);
310  G4VSolid* atomS_S = new G4Orb("Sphere", sphereSize*1.8);
311  G4VSolid* atomS_P = new G4Orb("Sphere", sphereSize*1.8);
312  G4VSolid* atomS_X = new G4Orb("Sphere", sphereSize); //Undefined
313 
314  //Logical volumes and color table for atoms :
315  G4LogicalVolume* atomLV_H = new G4LogicalVolume(
316  atomS_H, fpWaterMaterial, "atomLV_H"); // its solid, material, name
317  G4VisAttributes * MyVisAtt_H = new G4VisAttributes(
318  G4Colour(G4Colour::White()));
319  MyVisAtt_H->SetForceSolid(true);
320  atomLV_H->SetVisAttributes(MyVisAtt_H);
321 
322  G4LogicalVolume* atomLV_C = new G4LogicalVolume(
323  atomS_C, fpWaterMaterial, "atomLV_C"); // its solid, material, name
324  G4VisAttributes * MyVisAtt_C = new G4VisAttributes(
325  G4Colour(G4Colour::Gray()));//Black in CPK, but bad rendering
326  MyVisAtt_C->SetForceSolid(true);
327  atomLV_C->SetVisAttributes(MyVisAtt_C);
328 
329  G4LogicalVolume* atomLV_O = new G4LogicalVolume(
330  atomS_O, fpWaterMaterial, "atomLV_O"); // its solid, material, name
331  G4VisAttributes * MyVisAtt_O = new G4VisAttributes(
332  G4Colour(G4Colour::Red()));
333  MyVisAtt_O->SetForceSolid(true);
334  atomLV_O->SetVisAttributes(MyVisAtt_O);
335 
336  G4LogicalVolume* atomLV_N = new G4LogicalVolume(
337  atomS_N, fpWaterMaterial, "atomLV_N"); // its solid, material, name
338  G4VisAttributes * MyVisAtt_N = new G4VisAttributes(
339  G4Colour(G4Colour(0.,0.,0.5)));//Dark blue
340  MyVisAtt_N->SetForceSolid(true);
341  atomLV_N->SetVisAttributes(MyVisAtt_N);
342 
343  G4LogicalVolume* atomLV_S = new G4LogicalVolume(
344  atomS_S, fpWaterMaterial, "atomLV_S"); // its solid, material, name
345  G4VisAttributes * MyVisAtt_S = new G4VisAttributes(G4Colour(
346  G4Colour::Yellow()));
347  MyVisAtt_S->SetForceSolid(true);
348  atomLV_S->SetVisAttributes(MyVisAtt_S);
349 
350  G4LogicalVolume* atomLV_P = new G4LogicalVolume(
351  atomS_P, fpWaterMaterial, "atomLV_P"); // its solid, material, name
352  G4VisAttributes * MyVisAtt_P = new G4VisAttributes(
353  G4Colour(G4Colour(1.0,0.5,0.)));//Orange
354  MyVisAtt_P->SetForceSolid(true);
355  atomLV_P->SetVisAttributes(MyVisAtt_P);
356 
357  G4LogicalVolume* atomLV_X = new G4LogicalVolume(
358  atomS_X, fpWaterMaterial, "atomLV_X"); // its solid, material, name
359  G4VisAttributes * MyVisAtt_X = new G4VisAttributes(
360  G4Colour(G4Colour(1.0,0.75,0.8)));//Pink, other elements in CKP
361  MyVisAtt_X->SetForceSolid(true);
362  atomLV_X->SetVisAttributes(MyVisAtt_X);
363 
364  //Placement and physical volume construction from memory
365  Residue *residueListTemp;
366  Atom *AtomTemp;
367 
368  int nbAtomTot=0; //Number total of atoms
369  int nbAtomH=0, nbAtomC=0, nbAtomO=0, nbAtomN=0, nbAtomS=0, nbAtomP=0;
370  int nbAtomX=0;
371  int k=0;
372 
373  while (moleculeListTemp)
374  {
375  residueListTemp = moleculeListTemp->GetFirst();
376 
377  k++;
378  int j=0;
379  while (residueListTemp)
380  {
381  AtomTemp=residueListTemp->GetFirst();
382  j++;
383 
384  int startFrom=0;
385  int upTo=residueListTemp->fNbAtom; //case Base+sugar+phosphat (all atoms)
386  for (int i=0 ; i < (upTo + startFrom) ; i++) //Phosphat or Sugar or Base
387  {
388 
389  if (AtomTemp->fElement.compare("H") == 0)
390  {
391  nbAtomH++;
392  new G4PVPlacement(0,
393  G4ThreeVector(AtomTemp->fX*1*angstrom,
394  AtomTemp->fY*1*angstrom,
395  AtomTemp->fZ*1*angstrom),
396  atomLV_H,
397  "atomP",
398  worldLV,
399  false,
400  0,
401  fCheckOverlaps);
402  }
403  else if (AtomTemp->fElement.compare("C") == 0)
404  {
405  nbAtomC++;
406  new G4PVPlacement(0,
407  G4ThreeVector(AtomTemp->fX*1*angstrom,
408  AtomTemp->fY*1*angstrom,
409  AtomTemp->fZ*1*angstrom),
410  atomLV_C,
411  "atomP",
412  worldLV,
413  false,
414  0,
415  fCheckOverlaps);
416  }
417  else if (AtomTemp->fElement.compare("O") == 0)
418  {
419  nbAtomO++;
420  new G4PVPlacement(0,
421  G4ThreeVector(AtomTemp->fX*1*angstrom,
422  AtomTemp->fY*1*angstrom,
423  AtomTemp->fZ*1*angstrom),
424  atomLV_O,
425  "atomP",
426  worldLV,
427  false,
428  0,
429  fCheckOverlaps);
430  }
431  else if (AtomTemp->fElement.compare("N") == 0)
432  {
433  nbAtomN++;
434  new G4PVPlacement(0,
435  G4ThreeVector(AtomTemp->fX*1*angstrom,
436  AtomTemp->fY*1*angstrom,
437  AtomTemp->fZ*1*angstrom),
438  atomLV_N,
439  "atomP",
440  worldLV,
441  false,
442  0,
443  fCheckOverlaps);
444  }
445  else if (AtomTemp->fElement.compare("S") == 0)
446  {
447  nbAtomS++;
448  new G4PVPlacement(0,
449  G4ThreeVector(AtomTemp->fX*1*angstrom,
450  AtomTemp->fY*1*angstrom,
451  AtomTemp->fZ*1*angstrom),
452  atomLV_S,
453  "atomP",
454  worldLV,
455  false,
456  0,
457  fCheckOverlaps);
458  }
459  else if (AtomTemp->fElement.compare("P") == 0)
460  {
461  nbAtomP++;
462  new G4PVPlacement(0,
463  G4ThreeVector(AtomTemp->fX*1*angstrom,
464  AtomTemp->fY*1*angstrom,
465  AtomTemp->fZ*1*angstrom),
466  atomLV_P,
467  "atomP",
468  worldLV,
469  false,
470  0,
471  fCheckOverlaps);
472  }
473  else
474  {
475  nbAtomX++;
476  new G4PVPlacement(0,
477  G4ThreeVector(AtomTemp->fX*1*angstrom,
478  AtomTemp->fY*1*angstrom,
479  AtomTemp->fZ*1*angstrom),
480  atomLV_X,
481  "atomP",
482  worldLV,
483  false,
484  0,
485  fCheckOverlaps);
486  }
487 
488  nbAtomTot++;
489  //}//End if (i>=startFrom)
490  AtomTemp=AtomTemp->GetNext();
491  }//end of for ( i=0 ; i < residueListTemp->nbAtom ; i++)
492 
493  residueListTemp=residueListTemp->GetNext();
494  }//end of while (residueListTemp)
495 
496  moleculeListTemp=moleculeListTemp->GetNext();
497  }//end of while (moleculeListTemp)
498 
499  G4cout << "**************** atomisticView(...) ****************" <<G4endl;
500  G4cout << "Number of loaded chains = " << k <<G4endl;
501  G4cout << "Number of Atoms = " << nbAtomTot <<G4endl;
502  G4cout << "Number of Hydrogens = " << nbAtomH <<G4endl;
503  G4cout << "Number of Carbons = " << nbAtomC <<G4endl;
504  G4cout << "Number of Oxygens = " << nbAtomO <<G4endl;
505  G4cout << "Number of Nitrogens = " << nbAtomN <<G4endl;
506  G4cout << "Number of Sulfurs = " << nbAtomS <<G4endl;
507  G4cout << "Number of Phosphorus = " << nbAtomP <<G4endl;
508  G4cout << "Number of undifined atoms =" << nbAtomX <<G4endl<<G4endl;
509 }
510 
511 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
512 
515 //
516 void DetectorConstruction::BarycenterView(G4LogicalVolume* worldLV,
517  Barycenter *barycenterListTemp)
518 {
519  CheckMaterials();
520 
521  G4VSolid* atomS_ZZ;
522  G4LogicalVolume* atomLV_ZZ;
523  G4VisAttributes* MyVisAtt_ZZ;
524 
525  int k=0;
526 
527  while (barycenterListTemp)
528  {
529  k++;
530 
531  atomS_ZZ = new G4Orb("Sphere", (barycenterListTemp->GetRadius())*angstrom);
532  atomLV_ZZ = new G4LogicalVolume(atomS_ZZ,fpWaterMaterial,"atomLV_ZZ");
533  MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Magenta()));
534  MyVisAtt_ZZ->SetForceSolid(true);
535  atomLV_ZZ->SetVisAttributes(MyVisAtt_ZZ);
536 
537  new G4PVPlacement(0, G4ThreeVector(
538  barycenterListTemp->fCenterX*1*angstrom,
539  barycenterListTemp->fCenterY*1*angstrom,
540  barycenterListTemp->fCenterZ*1*angstrom),
541  atomLV_ZZ,
542  "atomZZ",
543  worldLV,
544  false,
545  0,
546  fCheckOverlaps);
547 
548  barycenterListTemp=barycenterListTemp->GetNext();
549  }//end of while (moleculeListTemp)
550 }
551 
552 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
553 
556 //
557 void DetectorConstruction::ResiduesView(G4LogicalVolume* worldLV,
558  Barycenter *barycenterListTemp)
559 {
560  CheckMaterials();
561  G4VisAttributes* MyVisAtt_ZZ;
562 
563  G4VSolid* tubS1_ZZ;
564  G4LogicalVolume* tubLV1_ZZ;
565  G4VSolid* tubS2_ZZ;
566  G4LogicalVolume* tubLV2_ZZ;
567 
568  G4VSolid* AS_ZZ;
569  G4LogicalVolume* ALV_ZZ;
570  G4VSolid* BS_ZZ;
571  G4LogicalVolume* BLV_ZZ;
572  G4VSolid* CS_ZZ;
573  G4LogicalVolume* CLV_ZZ;
574 
575  int k=0;
576 
577  while (barycenterListTemp)
578  {
579  k++;
580 
581  //3 spheres to Base, Sugar, Phosphate)
582  AS_ZZ = new G4Orb("Sphere", 1.*angstrom);
583  ALV_ZZ = new G4LogicalVolume(AS_ZZ,fpWaterMaterial, "ALV_ZZ");
584  MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Blue()));
585  MyVisAtt_ZZ->SetForceSolid(true);
586  ALV_ZZ->SetVisAttributes(MyVisAtt_ZZ);
587  new G4PVPlacement(0,
588  G4ThreeVector(barycenterListTemp->fCenterBaseX*angstrom,
589  barycenterListTemp->fCenterBaseY*angstrom,
590  barycenterListTemp->fCenterBaseZ*angstrom),
591  ALV_ZZ,
592  "AZZ",
593  worldLV,
594  false,
595  0,
596  fCheckOverlaps);
597 
598  BS_ZZ = new G4Orb("Sphere", 1.*angstrom);
599  BLV_ZZ = new G4LogicalVolume(BS_ZZ,fpWaterMaterial, "BLV_ZZ");
600  MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Red()));
601  MyVisAtt_ZZ->SetForceSolid(true);
602  BLV_ZZ->SetVisAttributes(MyVisAtt_ZZ);
603  new G4PVPlacement(0,
604  G4ThreeVector(
605  (barycenterListTemp->fCenterPhosphateX)*angstrom,
606  (barycenterListTemp->fCenterPhosphateY)*angstrom,
607  (barycenterListTemp->fCenterPhosphateZ)*angstrom),
608  BLV_ZZ,
609  "BZZ",
610  worldLV,
611  false,
612  0,
613  fCheckOverlaps);
614 
615  CS_ZZ = new G4Orb("Sphere", 1.*angstrom);
616  CLV_ZZ = new G4LogicalVolume(CS_ZZ,fpWaterMaterial, "CLV_ZZ");
617  MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Yellow()));
618  MyVisAtt_ZZ->SetForceSolid(true);
619  CLV_ZZ->SetVisAttributes(MyVisAtt_ZZ);
620  new G4PVPlacement(0,
621  G4ThreeVector(
622  barycenterListTemp->fCenterSugarX*angstrom,
623  barycenterListTemp->fCenterSugarY*angstrom,
624  barycenterListTemp->fCenterSugarZ*angstrom),
625  CLV_ZZ,
626  "CZZ",
627  worldLV,
628  false,
629  0,
630  fCheckOverlaps);
631 
632  //2 cylinders (Base<=>Sugar and Sugar<=>Phosphat)
633  // Cylinder Base<=>Sugar
634  tubS1_ZZ = new G4Tubs( "Cylinder",
635  0.,
636  0.5*angstrom,
637  std::sqrt (
638  (barycenterListTemp->fCenterBaseX-barycenterListTemp->fCenterSugarX)
639  * (barycenterListTemp->fCenterBaseX-barycenterListTemp->fCenterSugarX)
640  + (barycenterListTemp->fCenterBaseY-barycenterListTemp->fCenterSugarY)
641  * (barycenterListTemp->fCenterBaseY-barycenterListTemp->fCenterSugarY)
642  + (barycenterListTemp->fCenterBaseZ-barycenterListTemp->fCenterSugarZ)
643  * (barycenterListTemp->fCenterBaseZ-barycenterListTemp->fCenterSugarZ)
644  ) /2 *angstrom,
645  0.,
646  2.*pi );
647 
648  tubLV1_ZZ = new G4LogicalVolume(tubS1_ZZ,fpWaterMaterial, "tubLV_ZZ");
649  MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Green()));
650  MyVisAtt_ZZ->SetForceSolid(true);
651  tubLV1_ZZ->SetVisAttributes(MyVisAtt_ZZ);
652 
653  G4double Ux= barycenterListTemp->fCenterBaseX-
654  barycenterListTemp->fCenterSugarX;
655  G4double Uy= barycenterListTemp->fCenterBaseY-
656  barycenterListTemp->fCenterSugarY;
657  G4double Uz= barycenterListTemp->fCenterBaseZ-
658  barycenterListTemp->fCenterSugarZ;
659  G4double llUll=std::sqrt(Ux*Ux+Uy*Uy+Uz*Uz);
660 
661  Ux=Ux/llUll;
662  Uy=Uy/llUll;
663  Uz=Uz/llUll;
664 
665  G4ThreeVector direction = G4ThreeVector(Ux,Uy,Uz);
666  G4double theta_euler = direction.theta();
667  G4double phi_euler = direction.phi();
668  G4double psi_euler = 0;
669 
670  //Warning : clhep Euler constructor build inverse matrix !
671  G4RotationMatrix rotm1Inv = G4RotationMatrix(phi_euler+pi/2,
672  theta_euler,
673  psi_euler);
674  G4RotationMatrix rotm1 = rotm1Inv.inverse();
675  G4ThreeVector translm1 = G4ThreeVector(
676  (barycenterListTemp->fCenterBaseX+barycenterListTemp->fCenterSugarX)/2.
677  *angstrom,
678  (barycenterListTemp->fCenterBaseY+barycenterListTemp->fCenterSugarY)/2.
679  *angstrom,
680  (barycenterListTemp->fCenterBaseZ+barycenterListTemp->fCenterSugarZ)/2.
681  *angstrom);
682  G4Transform3D transform1 = G4Transform3D(rotm1,translm1);
683  new G4PVPlacement(transform1, // rotation translation
684  tubLV1_ZZ,"atomZZ",worldLV,false, 0, fCheckOverlaps);
685 
686  //Cylinder Sugar<=>Phosphat
687  tubS2_ZZ = new G4Tubs( "Cylinder2",
688  0.,
689  0.5*angstrom,
690  std::sqrt (
691  (barycenterListTemp->fCenterSugarX-barycenterListTemp->fCenterPhosphateX)
692  * (barycenterListTemp->fCenterSugarX-barycenterListTemp->fCenterPhosphateX)
693  + (barycenterListTemp->fCenterSugarY-barycenterListTemp->fCenterPhosphateY)
694  * (barycenterListTemp->fCenterSugarY-barycenterListTemp->fCenterPhosphateY)
695  + (barycenterListTemp->fCenterSugarZ-barycenterListTemp->fCenterPhosphateZ)
696  * (barycenterListTemp->fCenterSugarZ-barycenterListTemp->fCenterPhosphateZ)
697  ) /2 *angstrom,
698  0.,
699  2.*pi );
700 
701  tubLV2_ZZ = new G4LogicalVolume(tubS2_ZZ, fpWaterMaterial, "tubLV2_ZZ");
702  MyVisAtt_ZZ = new G4VisAttributes(G4Colour(1,0.5,0));
703  MyVisAtt_ZZ->SetForceSolid(true);
704  tubLV2_ZZ->SetVisAttributes(MyVisAtt_ZZ);
705 
706  Ux= barycenterListTemp->fCenterSugarX-
707  barycenterListTemp->fCenterPhosphateX;
708  Uy= barycenterListTemp->fCenterSugarY-
709  barycenterListTemp->fCenterPhosphateY;
710  Uz= barycenterListTemp->fCenterSugarZ-
711  barycenterListTemp->fCenterPhosphateZ;
712  llUll=std::sqrt(Ux*Ux+Uy*Uy+Uz*Uz);
713 
714  Ux=Ux/llUll;
715  Uy=Uy/llUll;
716  Uz=Uz/llUll;
717 
718  direction = G4ThreeVector(Ux,Uy,Uz);
719  theta_euler = direction.theta();
720  phi_euler = direction.phi();
721  psi_euler = 0;
722 
723  //Warning : clhep Euler constructor build inverse matrix !
724  rotm1Inv = G4RotationMatrix(phi_euler+pi/2,theta_euler,psi_euler);
725  rotm1 = rotm1Inv.inverse();
726  translm1 = G4ThreeVector(
727  (barycenterListTemp->fCenterSugarX+
728  barycenterListTemp->fCenterPhosphateX)/2.*angstrom,
729  (barycenterListTemp->fCenterSugarY+
730  barycenterListTemp->fCenterPhosphateY)/2.*angstrom,
731  (barycenterListTemp->fCenterSugarZ+
732  barycenterListTemp->fCenterPhosphateZ)/2.*angstrom);
733  transform1 = G4Transform3D(rotm1,translm1);
734  new G4PVPlacement(transform1, // rotation translation
735  tubLV2_ZZ,
736  "atomZZ",
737  worldLV,
738  false,
739  0,
740  fCheckOverlaps);
741 
742  barycenterListTemp=barycenterListTemp->GetNext();
743  }//end of while sur moleculeListTemp
744 }
745 
746 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
747 
750 //
751 void DetectorConstruction::DrawBoundingVolume(G4LogicalVolume* worldLV,
752  Molecule *moleculeListTemp)
753 {
754  CheckMaterials();
755 
756  double dX,dY,dZ;//Dimensions for bounding volume
757  double tX,tY,tZ;//Translation for bounding volume
758  fPDBlib.ComputeBoundingVolumeParams(moleculeListTemp,
759  dX, dY, dZ,
760  tX, tY, tZ);
761 
762  //[Geometry] Build a box :
763  G4VSolid* boundingS
764  = new G4Box("Bounding", dX*1*angstrom, dY*1*angstrom, dZ*1*angstrom);
765 
766  G4LogicalVolume* boundingLV
767  = new G4LogicalVolume(boundingS, fpWaterMaterial, "BoundingLV");
768 
769  G4RotationMatrix *pRot = new G4RotationMatrix();
770 
771  G4VisAttributes *MyVisAtt_ZZ = new G4VisAttributes(
772  G4Colour(G4Colour::Gray()));
773  boundingLV->SetVisAttributes(MyVisAtt_ZZ);
774 
775  new G4PVPlacement(pRot,
776  G4ThreeVector(
777  tX*1*angstrom,
778  tY*1*angstrom,
779  tZ*1*angstrom), // at (x,y,z)
780  boundingLV,
781  "boundingPV",
782  worldLV
783  ,false,
784  0,
785  fCheckOverlaps);
786 }
787 
788 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
789 
790 void DetectorConstruction::LoadPDBfile(G4String fileName)
791 {
792  G4cout<<"Load PDB file : "<<fileName<<"."<<G4endl<<G4endl;
793  fPdbFileName=fileName;
794 #ifdef G4MULTITHREADED
795  G4MTRunManager::GetRunManager()->DefineWorldVolume(
796  DefineVolumes(fPdbFileName,fChosenOption)
797  );
798  G4MTRunManager::GetRunManager()->ReinitializeGeometry();
799 #else
800  G4RunManager::GetRunManager()->DefineWorldVolume(
801  DefineVolumes(fPdbFileName,fChosenOption)
802  );
803 #endif
804 }
805 
806 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
807 
808 void DetectorConstruction::BuildBoundingVolume()
809 {
810  if (fPdbFileStatus>0) //a PDB file has been loaded
811  {
812  G4cout<<"Build only world volume and bounding volume"
813  " for computation."<<G4endl<<G4endl;
814 #ifdef G4MULTITHREADED
815  G4MTRunManager::GetRunManager()->DefineWorldVolume(
816  DefineVolumes(fPdbFileName,10)
817  );
818  G4MTRunManager::GetRunManager()->ReinitializeGeometry();
819 #else
820  G4RunManager::GetRunManager()->DefineWorldVolume(
821  DefineVolumes(fPdbFileName,10)
822  );
823 #endif
824  }
825  else G4cout<<"PDB file not found!"<<G4endl<<G4endl;
826 }
827 
828 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
829 
830 void DetectorConstruction::DrawAtoms_()
831 {
832  if (fPdbFileStatus>0) //a PDB file has been loaded
833  {
834 #ifdef G4MULTITHREADED
835  G4MTRunManager::GetRunManager()->DefineWorldVolume(
836  DefineVolumes(fPdbFileName,1)
837  );
838  G4MTRunManager::GetRunManager()->ReinitializeGeometry();
839 #else
840  G4RunManager::GetRunManager()->DefineWorldVolume(
841  DefineVolumes(fPdbFileName,1)
842  );
843 #endif
844  }
845  else G4cout<<"PDB file not found!"<<G4endl<<G4endl;
846 }
847 
848 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
849 
850 void DetectorConstruction::DrawNucleotides_()
851 {
852  if (fPdbFileStatus>0) //a PDB file has been loaded
853  {
854 #ifdef G4MULTITHREADED
855  G4MTRunManager::GetRunManager()->DefineWorldVolume(
856  DefineVolumes(fPdbFileName,2)
857  );
858  G4MTRunManager::GetRunManager()->ReinitializeGeometry();
859 #else
860  G4RunManager::GetRunManager()->DefineWorldVolume(
861  DefineVolumes(fPdbFileName,2)
862  );
863 #endif
864  }
865  else G4cout<<"PDB file not found!"<<G4endl<<G4endl;
866 }
867 
868 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
869 
870 void DetectorConstruction::DrawResidues_()
871 {
872  if (fPdbFileStatus>0) //a PDB file has been loaded
873  {
874 #ifdef G4MULTITHREADED
875  G4MTRunManager::GetRunManager()->DefineWorldVolume(
876  DefineVolumes(fPdbFileName,3)
877  );
878  G4MTRunManager::GetRunManager()->ReinitializeGeometry();
879 #else
880  G4RunManager::GetRunManager()->DefineWorldVolume(
881  DefineVolumes(fPdbFileName,3)
882  );
883 #endif
884  }
885  else G4cout<<"PDB file not found!"<<G4endl<<G4endl;
886 }
887 
888 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
889 
890 void DetectorConstruction::DrawAtomsWithBounding_()
891 {
892  if (fPdbFileStatus>0) //a PDB file has been loaded
893  {
894 #ifdef G4MULTITHREADED
895  G4MTRunManager::GetRunManager()->DefineWorldVolume(
896  DefineVolumes(fPdbFileName,11)
897  );
898  G4MTRunManager::GetRunManager()->ReinitializeGeometry();
899 #else
900  G4RunManager::GetRunManager()->DefineWorldVolume(
901  DefineVolumes(fPdbFileName,11)
902  );
903 #endif
904  }
905  else G4cout<<"PDB file not found!"<<G4endl<<G4endl;
906 }
907 
908 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
909 
910 void DetectorConstruction::DrawNucleotidesWithBounding_()
911 {
912  if (fPdbFileStatus>0) //a PDB file has been loaded
913  {
914 #ifdef G4MULTITHREADED
915  G4MTRunManager::GetRunManager()->DefineWorldVolume(
916  DefineVolumes(fPdbFileName,12)
917  );
918  G4MTRunManager::GetRunManager()->ReinitializeGeometry();
919 #else
920  G4RunManager::GetRunManager()->DefineWorldVolume(
921  DefineVolumes(fPdbFileName,12)
922  );
923 #endif
924  }
925  else G4cout<<"PDB file not found!"<<G4endl<<G4endl;
926 }
927 
928 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
929 
930 void DetectorConstruction::DrawResiduesWithBounding_()
931 {
932  if (fPdbFileStatus>0) //a PDB file has been loaded
933  {
934 #ifdef G4MULTITHREADED
935  G4MTRunManager::GetRunManager()->DefineWorldVolume(
936  DefineVolumes(fPdbFileName,13)
937  );
938  G4MTRunManager::GetRunManager()->ReinitializeGeometry();
939 #else
940  G4RunManager::GetRunManager()->DefineWorldVolume(
941  DefineVolumes(fPdbFileName,13)
942  );
943 #endif
944  }
945  else G4cout<<"PDB file not found!"<<G4endl<<G4endl;
946 }
Barycenter::fCenterBaseX
double fCenterBaseX
&quot;X coordinate&quot; of this Base Barycenter
Definition: PDBbarycenter.hh:92
G4Colour::Gray
static G4Colour Gray()
Definition: G4Colour.hh:144
Barycenter::fCenterPhosphateZ
double fCenterPhosphateZ
&quot;Z coordinate&quot; of this Phosphate Barycenter
Definition: PDBbarycenter.hh:102
Barycenter
Molecule Class.
Definition: PDBbarycenter.hh:54
Barycenter::fCenterPhosphateX
double fCenterPhosphateX
&quot;X coordinate&quot; of this Phosphate Barycenter
Definition: PDBbarycenter.hh:100
Barycenter::fCenterPhosphateY
double fCenterPhosphateY
&quot;Y coordinate&quot; of this Phosphate Barycenter
Definition: PDBbarycenter.hh:101
G4NistManager::FindOrBuildMaterial
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
Definition: G4NistManager.hh:474
PDBlib
PDBlib Class.
Definition: PDBlib.hh:57
PDBlib::Load
Molecule * Load(const std::string &filename, unsigned short int &isDNA, unsigned short int verbose)
Load PDB file into memory.
Definition: PDBlib.cc:83
Molecule
Molecule Class.
Definition: PDBmolecule.hh:58
Barycenter::fCenterY
double fCenterY
&quot;Y coordinate&quot; of this nucelotide Barycenter
Definition: PDBbarycenter.hh:89
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4RunManager.hh
G4RotationMatrix
CLHEP::HepRotation G4RotationMatrix
Definition: G4RotationMatrix.hh:43
G4Colour::Green
static G4Colour Green()
Definition: G4Colour.hh:149
Barycenter::fCenterX
double fCenterX
&quot;X coordinate&quot; of this nucelotide Barycenter
Definition: PDBbarycenter.hh:88
Barycenter::fCenterSugarX
double fCenterSugarX
&quot;X coordinate&quot; of this Sugar Barycenter
Definition: PDBbarycenter.hh:96
a
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
G4Material::GetMaterial
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:602
Residue::fNbAtom
int fNbAtom
Number of atoms into a residue (usage before vector)
Definition: PDBresidue.hh:88
python.hepunit.universe_mean_density
int universe_mean_density
Definition: hepunit.py:307
G4Box
Definition: G4Box.hh:64
G4NistManager.hh
PDBlib::ComputeNucleotideBarycenters
Barycenter * ComputeNucleotideBarycenters(Molecule *moleculeListTemp)
Compute nucleotide barycenter from memory.
Definition: PDBlib.cc:452
DetectorConstruction::Construct
G4VPhysicalVolume * Construct()
Definition: DetectorConstruction.cc:70
G4Tubs
Definition: G4Tubs.hh:85
Molecule::GetFirst
Residue * GetFirst()
Get the first Residue.
Definition: PDBmolecule.cc:82
Residue::GetNext
Residue * GetNext()
Get the next residue.
Definition: PDBresidue.cc:67
G4Colour::Magenta
static G4Colour Magenta()
Definition: G4Colour.hh:152
G4SolidStore::Clean
static void Clean()
Definition: G4SolidStore.cc:79
G4Box.hh
Atom::fZ
double fZ
Z orthogonal coordinates in Angstroms.
Definition: PDBatom.hh:104
Residue
Residue Class.
Definition: PDBresidue.hh:58
Barycenter::fCenterSugarY
double fCenterSugarY
&quot;Y coordinate&quot; of this Sugar Barycenter
Definition: PDBbarycenter.hh:97
kStateGas
Definition: G4Material.hh:114
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:68
Barycenter::fCenterBaseZ
double fCenterBaseZ
&quot;Z coordinate&quot; of this Base Barycenter
Definition: PDBbarycenter.hh:94
G4VisAttributes::SetForceSolid
void SetForceSolid(G4bool=true)
Definition: G4VisAttributes.cc:173
CLHEP::HepRotation::inverse
HepRotation inverse() const
DetectorConstruction::GetBarycenterList
Barycenter * GetBarycenterList()
Definition: DetectorConstruction.cc:282
G4RunManager::DefineWorldVolume
virtual void DefineWorldVolume(G4VPhysicalVolume *worldVol, G4bool topologyIsChanged=true)
Definition: G4RunManager.cc:653
G4PhysicalVolumeStore::GetInstance
static G4PhysicalVolumeStore * GetInstance()
Definition: G4PhysicalVolumeStore.cc:187
G4SolidStore.hh
g
function g(Y1, Y2, PT2)
Definition: hijing1.383.f:5205
G4VisAttributes.hh
G4cout
G4GLOB_DLL std::ostream G4cout
G4Tubs.hh
Atom
Atom Class.
Definition: PDBatom.hh:56
G4VisAttributes::SetVisibility
void SetVisibility(G4bool=true)
G4bool
bool G4bool
Definition: G4Types.hh:79
PDBlib::ComputeNbNucleotidsPerStrand
void ComputeNbNucleotidsPerStrand(Molecule *moleculeListTemp)
Compute number of nucleotide per strand.
Definition: PDBlib.cc:726
G4Material.hh
G4Colour.hh
Barycenter::fCenterSugarZ
double fCenterSugarZ
&quot;Z coordinate&quot; of this Sugar Barycenter
Definition: PDBbarycenter.hh:98
pascal
#define pascal
Definition: SystemOfUnits.h:206
G4LogicalVolumeStore::GetInstance
static G4LogicalVolumeStore * GetInstance()
Definition: G4LogicalVolumeStore.cc:181
G4SolidStore::GetInstance
static G4SolidStore * GetInstance()
Definition: G4SolidStore.cc:185
G4Transform3D
HepGeom::Transform3D G4Transform3D
Definition: G4Transform3D.hh:35
G4Colour::Blue
static G4Colour Blue()
Definition: G4Colour.hh:150
Molecule::GetNext
Molecule * GetNext()
information about molecule (not implemented)
Definition: PDBmolecule.cc:75
CLHEP::Hep3Vector::phi
double phi() const
G4GeometryManager::GetInstance
static G4GeometryManager * GetInstance()
Definition: G4GeometryManager.cc:135
kelvin
static constexpr double kelvin
Definition: G4SIunits.hh:281
G4Orb
Definition: G4Orb.hh:61
G4RunManager::ReinitializeGeometry
void ReinitializeGeometry(G4bool destroyFirst=false, G4bool prop=true)
Definition: G4RunManager.cc:923
Barycenter::GetRadius
double GetRadius()
Get the distance between the farther atom and nucleotide barycenter.
Definition: PDBbarycenter.cc:127
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4PhysicalConstants.hh
CLHEP::Hep3Vector::theta
double theta() const
Barycenter::GetNext
Barycenter * GetNext()
Get the next Barycenter.
Definition: PDBbarycenter.cc:87
G4RunManager::GetRunManager
static G4RunManager * GetRunManager()
Definition: G4RunManager.cc:79
G4LogicalVolume.hh
G4PhysicalVolumeStore.hh
G4VPhysicalVolume::GetLogicalVolume
G4LogicalVolume * GetLogicalVolume() const
Residue::GetFirst
Atom * GetFirst()
Get the first atom.
Definition: PDBresidue.cc:74
FatalException
Definition: G4ExceptionSeverity.hh:60
PDBlib::ComputeBoundingVolumeParams
void ComputeBoundingVolumeParams(Molecule *moleculeListTemp, double &dX, double &dY, double &dZ, double &tX, double &tY, double &tZ)
Compute the corresponding bounding volume parameters.
Definition: PDBlib.cc:662
test.z
tuple z
Definition: test.py:28
Atom::fX
double fX
X orthogonal coordinates in Angstroms.
Definition: PDBatom.hh:102
DetectorConstruction::LoadPDBfile
void LoadPDBfile(G4String fileName)
Definition: DetectorConstruction.cc:790
Barycenter::fCenterBaseY
double fCenterBaseY
&quot;Y coordinate&quot; of this Base Barycenter
Definition: PDBbarycenter.hh:93
Atom::GetNext
Atom * GetNext()
Returns the next Atom.
Definition: PDBatom.cc:67
G4endl
#define G4endl
Definition: G4ios.hh:61
angstrom
static constexpr double angstrom
Definition: G4SIunits.hh:102
Atom::fY
double fY
Y orthogonal coordinates in Angstroms.
Definition: PDBatom.hh:103
G4GeometryManager::OpenGeometry
void OpenGeometry(G4VPhysicalVolume *vol=0)
Definition: G4GeometryManager.cc:105
pi
static constexpr double pi
Definition: G4SIunits.hh:75
Atom::fElement
std::string fElement
Element symbol extracted from &#39;atom name&#39;.
Definition: PDBatom.hh:107
G4Colour::Red
static G4Colour Red()
Definition: G4Colour.hh:148
Barycenter::fCenterZ
double fCenterZ
&quot;Z coordinate&quot; of this nucelotide Barycenter
Definition: PDBbarycenter.hh:90
G4MTRunManager.hh
G4Orb.hh
G4LogicalVolumeStore.hh
G4double
double G4double
Definition: G4Types.hh:76
G4SystemOfUnits.hh
G4PVPlacement.hh
G4Colour::Yellow
static G4Colour Yellow()
Definition: G4Colour.hh:153
mole
static constexpr double mole
Definition: G4SIunits.hh:286
G4LogicalVolume::SetVisAttributes
void SetVisAttributes(const G4VisAttributes *pVA)
G4GeometryManager.hh
G4Colour::White
static G4Colour White()
Definition: G4Colour.hh:143