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