Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-05-07 08:06:39 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file DetectorConstruction.cc
0027 /// \brief Implementation of the DetectorConstruction class
0028 
0029 // This example is provided by the Geant4-DNA collaboration
0030 // Any report or published results obtained using the Geant4-DNA software
0031 // shall cite the following Geant4-DNA collaboration publication:
0032 // Med. Phys. 37 (2010) 4692-4708
0033 // Delage et al. PDB4DNA: implementation of DNA geometry from the Protein Data
0034 //                  Bank (PDB) description for Geant4-DNA Monte-Carlo
0035 //                  simulations (submitted to Comput. Phys. Commun.)
0036 // The Geant4-DNA web site is available at http://geant4-dna.org
0037 //
0038 //
0039 
0040 #include "DetectorConstruction.hh"
0041 
0042 #include "DetectorMessenger.hh"
0043 
0044 #include "G4Box.hh"
0045 #include "G4Colour.hh"
0046 #include "G4GeometryManager.hh"
0047 #include "G4LogicalVolume.hh"
0048 #include "G4LogicalVolumeStore.hh"
0049 #include "G4Material.hh"
0050 #include "G4NistManager.hh"
0051 #include "G4Orb.hh"
0052 #include "G4PVPlacement.hh"
0053 #include "G4PhysicalConstants.hh"
0054 #include "G4PhysicalVolumeStore.hh"
0055 #include "G4SolidStore.hh"
0056 #include "G4SystemOfUnits.hh"
0057 #include "G4Tubs.hh"
0058 #include "G4VisAttributes.hh"
0059 
0060 #ifdef G4MULTITHREADED
0061 #  include "G4MTRunManager.hh"
0062 #else
0063 #  include "G4RunManager.hh"
0064 #endif
0065 
0066 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0067 
0068 DetectorConstruction::DetectorConstruction()
0069   : G4VUserDetectorConstruction(), fpMessenger(0), fCheckOverlaps(false)
0070 {
0071   // Select a PDB file name by default
0072   // otherwise modified by the LoadPDBfile messenger
0073   //
0074   fPdbFileName = G4String("1ZBB.pdb");
0075   fPdbFileStatus = 0;
0076   fChosenOption = 11;
0077 
0078   fpDefaultMaterial = 0;
0079   fpWaterMaterial = 0;
0080   fpMessenger = new DetectorMessenger(this);
0081 }
0082 
0083 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0084 
0085 DetectorConstruction::~DetectorConstruction() {}
0086 
0087 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0088 
0089 G4VPhysicalVolume* DetectorConstruction::Construct()
0090 {
0091   fChosenOption = 11;  // Draw atomic with bounding box (visualisation by default)
0092   fPdbFileStatus = 0;  // There is no PDB file loaded at this stage
0093 
0094   // Define materials and geometry
0095   G4VPhysicalVolume* worldPV;
0096   worldPV = DefineVolumes(fPdbFileName, fChosenOption);
0097   return worldPV;
0098 }
0099 
0100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0101 
0102 void DetectorConstruction::ConstructMaterials()
0103 {
0104   //[G4_WATER]
0105   //
0106   G4NistManager* nistManager = G4NistManager::Instance();
0107   G4bool fromIsotopes = false;
0108   nistManager->FindOrBuildMaterial("G4_WATER", fromIsotopes);
0109   fpWaterMaterial = G4Material::GetMaterial("G4_WATER");
0110 
0111   //[Vacuum]
0112   G4double a;  // mass of a mole;
0113   G4double z;  // z=mean number of protons;
0114   G4double density;
0115   new G4Material("Galactic", z = 1., a = 1.01 * g / mole, density = universe_mean_density,
0116                  kStateGas, 2.73 * kelvin, 3.e-18 * pascal);
0117   fpDefaultMaterial = G4Material::GetMaterial("Galactic");
0118 }
0119 
0120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0121 
0122 void DetectorConstruction::CheckMaterials()
0123 {
0124   if (!fpDefaultMaterial)
0125     G4Exception("DetectorConstruction::CheckMaterials", "DEFAULT_MATERIAL_NOT_INIT_1",
0126                 FatalException, "Default material not initialized.");
0127 
0128   if (!fpWaterMaterial)
0129     G4Exception("DetectorConstruction::CheckMaterials", "WATER_MATERIAL_NOT_INIT", FatalException,
0130                 "Water material not initialized.");
0131 }
0132 
0133 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0134 
0135 G4VPhysicalVolume* DetectorConstruction::ConstructWorld()
0136 {
0137   // Geometry parameters
0138   G4double worldSize = 1000 * 1 * angstrom;
0139 
0140   if (!fpDefaultMaterial) {
0141     G4Exception("DetectorConstruction::ConstructWorld", "DEFAULT_MATERIAL_NOT_INIT_2",
0142                 FatalException, "Default material not initialized.");
0143   }
0144 
0145   //
0146   // World
0147   //
0148   G4VSolid* worldS = new G4Box("World",  // its name
0149                                worldSize / 2, worldSize / 2, worldSize / 2);  // its size
0150 
0151   G4LogicalVolume* worldLV = new G4LogicalVolume(worldS,  // its solid
0152                                                  fpDefaultMaterial,  // its material
0153                                                  "World");  // its name
0154 
0155   G4VisAttributes* MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Gray()));
0156   MyVisAtt_ZZ->SetVisibility(false);
0157   worldLV->SetVisAttributes(MyVisAtt_ZZ);
0158 
0159   G4VPhysicalVolume* worldPV = new G4PVPlacement(0,  // no rotation
0160                                                  G4ThreeVector(),  // at (0,0,0)
0161                                                  worldLV,  // its logical volume
0162                                                  "World",  // its name
0163                                                  0,  // its mother  volume
0164                                                  false,  // no boolean operation
0165                                                  0,  // copy number
0166                                                  true);  // checking overlaps forced to YES
0167 
0168   return worldPV;
0169 }
0170 
0171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0172 
0173 G4VPhysicalVolume* DetectorConstruction::DefineVolumes(G4String filename, unsigned short int option)
0174 {
0175   // Clean old geometry
0176   //
0177   G4GeometryManager::GetInstance()->OpenGeometry();
0178   G4PhysicalVolumeStore::GetInstance()->Clean();
0179   G4LogicalVolumeStore::GetInstance()->Clean();
0180   G4SolidStore::GetInstance()->Clean();
0181 
0182   // Define Materials
0183   //
0184   ConstructMaterials();
0185 
0186   // Reconstruct world not to superimpose on geometries
0187   G4VPhysicalVolume* worldPV;
0188   G4LogicalVolume* worldLV;
0189   worldPV = ConstructWorld();
0190   worldLV = worldPV->GetLogicalVolume();
0191 
0192   // List of molecules inside pdb file separated by TER keyword
0193   fpMoleculeList = NULL;
0194   fpBarycenterList = NULL;
0195 
0196   //'fPDBlib.load' is currently done for each 'DefineVolumes' call.
0197   int verbosity = 0;  // To be implemented
0198   unsigned short int isDNA;
0199   fpMoleculeList = fPDBlib.Load(filename, isDNA, verbosity);
0200   if (fpMoleculeList != NULL && isDNA == 1) {
0201     fPDBlib.ComputeNbNucleotidsPerStrand(fpMoleculeList);
0202     fpBarycenterList = fPDBlib.ComputeNucleotideBarycenters(fpMoleculeList);
0203     G4cout << "This PDB file is DNA." << G4endl;
0204   }
0205 
0206   if (fpMoleculeList != NULL) {
0207     fPdbFileStatus = 1;
0208   }
0209 
0210   if (option == 1) {
0211     AtomisticView(worldLV, fpMoleculeList, 1.0);
0212   }
0213   else if (option == 2) {
0214     BarycenterView(worldLV, fpBarycenterList);
0215   }
0216   else if (option == 3) {
0217     ResiduesView(worldLV, fpBarycenterList);
0218   }
0219   else if (option == 10) {
0220     DrawBoundingVolume(worldLV, fpMoleculeList);
0221   }
0222   else if (option == 11) {
0223     AtomisticView(worldLV, fpMoleculeList, 1.0);
0224     DrawBoundingVolume(worldLV, fpMoleculeList);
0225   }
0226   else if (option == 12) {
0227     BarycenterView(worldLV, fpBarycenterList);
0228     DrawBoundingVolume(worldLV, fpMoleculeList);
0229   }
0230   else if (option == 13) {
0231     ResiduesView(worldLV, fpBarycenterList);
0232     DrawBoundingVolume(worldLV, fpMoleculeList);
0233   }
0234   // Always return the physical World
0235   //
0236   return worldPV;
0237 }
0238 
0239 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0240 
0241 PDBlib DetectorConstruction::GetPDBlib()
0242 {
0243   return fPDBlib;
0244 }
0245 
0246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0247 
0248 Barycenter* DetectorConstruction::GetBarycenterList()
0249 {
0250   return fpBarycenterList;
0251 }
0252 
0253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0254 
0255 Molecule* DetectorConstruction::GetMoleculeList()
0256 {
0257   return fpMoleculeList;
0258 }
0259 
0260 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0261 
0262 //////////////////////////////////////////////////
0263 ///////////////// BEGIN atomistic representation
0264 //
0265 void DetectorConstruction::AtomisticView(G4LogicalVolume* worldLV, Molecule* moleculeListTemp,
0266                                          double atomSizeFactor)
0267 {
0268   CheckMaterials();
0269 
0270   // All solids are defined for different color attributes :
0271   G4double sphereSize = atomSizeFactor * 1 * angstrom;
0272   G4VSolid* atomS_H = new G4Orb("Sphere", sphereSize * 1.2);
0273   G4VSolid* atomS_C = new G4Orb("Sphere", sphereSize * 1.7);
0274   G4VSolid* atomS_O = new G4Orb("Sphere", sphereSize * 1.52);
0275   G4VSolid* atomS_N = new G4Orb("Sphere", sphereSize * 1.55);
0276   G4VSolid* atomS_S = new G4Orb("Sphere", sphereSize * 1.8);
0277   G4VSolid* atomS_P = new G4Orb("Sphere", sphereSize * 1.8);
0278   G4VSolid* atomS_X = new G4Orb("Sphere", sphereSize);  // Undefined
0279 
0280   // Logical volumes and color table for atoms :
0281   G4LogicalVolume* atomLV_H =
0282     new G4LogicalVolume(atomS_H, fpWaterMaterial, "atomLV_H");  // its solid, material, name
0283   G4VisAttributes* MyVisAtt_H = new G4VisAttributes(G4Colour(G4Colour::White()));
0284   MyVisAtt_H->SetForceSolid(true);
0285   atomLV_H->SetVisAttributes(MyVisAtt_H);
0286 
0287   G4LogicalVolume* atomLV_C =
0288     new G4LogicalVolume(atomS_C, fpWaterMaterial, "atomLV_C");  // its solid, material, name
0289   G4VisAttributes* MyVisAtt_C =
0290     new G4VisAttributes(G4Colour(G4Colour::Gray()));  // Black in CPK, but bad rendering
0291   MyVisAtt_C->SetForceSolid(true);
0292   atomLV_C->SetVisAttributes(MyVisAtt_C);
0293 
0294   G4LogicalVolume* atomLV_O =
0295     new G4LogicalVolume(atomS_O, fpWaterMaterial, "atomLV_O");  // its solid, material, name
0296   G4VisAttributes* MyVisAtt_O = new G4VisAttributes(G4Colour(G4Colour::Red()));
0297   MyVisAtt_O->SetForceSolid(true);
0298   atomLV_O->SetVisAttributes(MyVisAtt_O);
0299 
0300   G4LogicalVolume* atomLV_N =
0301     new G4LogicalVolume(atomS_N, fpWaterMaterial, "atomLV_N");  // its solid, material, name
0302   G4VisAttributes* MyVisAtt_N = new G4VisAttributes(G4Colour(G4Colour(0., 0., 0.5)));  // Dark blue
0303   MyVisAtt_N->SetForceSolid(true);
0304   atomLV_N->SetVisAttributes(MyVisAtt_N);
0305 
0306   G4LogicalVolume* atomLV_S =
0307     new G4LogicalVolume(atomS_S, fpWaterMaterial, "atomLV_S");  // its solid, material, name
0308   G4VisAttributes* MyVisAtt_S = new G4VisAttributes(G4Colour(G4Colour::Yellow()));
0309   MyVisAtt_S->SetForceSolid(true);
0310   atomLV_S->SetVisAttributes(MyVisAtt_S);
0311 
0312   G4LogicalVolume* atomLV_P =
0313     new G4LogicalVolume(atomS_P, fpWaterMaterial, "atomLV_P");  // its solid, material, name
0314   G4VisAttributes* MyVisAtt_P = new G4VisAttributes(G4Colour(G4Colour(1.0, 0.5, 0.)));  // Orange
0315   MyVisAtt_P->SetForceSolid(true);
0316   atomLV_P->SetVisAttributes(MyVisAtt_P);
0317 
0318   G4LogicalVolume* atomLV_X =
0319     new G4LogicalVolume(atomS_X, fpWaterMaterial, "atomLV_X");  // its solid, material, name
0320   G4VisAttributes* MyVisAtt_X =
0321     new G4VisAttributes(G4Colour(G4Colour(1.0, 0.75, 0.8)));  // Pink, other elements in CKP
0322   MyVisAtt_X->SetForceSolid(true);
0323   atomLV_X->SetVisAttributes(MyVisAtt_X);
0324 
0325   // Placement and physical volume construction from memory
0326   Residue* residueListTemp;
0327   Atom* AtomTemp;
0328 
0329   int nbAtomTot = 0;  // Number total of atoms
0330   int nbAtomH = 0, nbAtomC = 0, nbAtomO = 0, nbAtomN = 0, nbAtomS = 0, nbAtomP = 0;
0331   int nbAtomX = 0;
0332   int k = 0;
0333 
0334   while (moleculeListTemp) {
0335     residueListTemp = moleculeListTemp->GetFirst();
0336 
0337     k++;
0338     while (residueListTemp) {
0339       AtomTemp = residueListTemp->GetFirst();
0340 
0341       int startFrom = 0;
0342       int upTo = residueListTemp->fNbAtom;  // case Base+sugar+phosphat (all atoms)
0343       for (int i = 0; i < (upTo + startFrom); i++)  // Phosphat or Sugar or Base
0344       {
0345         if (AtomTemp->fElement.compare("H") == 0) {
0346           nbAtomH++;
0347           new G4PVPlacement(0,
0348                             G4ThreeVector(AtomTemp->fX * 1 * angstrom, AtomTemp->fY * 1 * angstrom,
0349                                           AtomTemp->fZ * 1 * angstrom),
0350                             atomLV_H, "atomP", worldLV, false, 0, fCheckOverlaps);
0351         }
0352         else if (AtomTemp->fElement.compare("C") == 0) {
0353           nbAtomC++;
0354           new G4PVPlacement(0,
0355                             G4ThreeVector(AtomTemp->fX * 1 * angstrom, AtomTemp->fY * 1 * angstrom,
0356                                           AtomTemp->fZ * 1 * angstrom),
0357                             atomLV_C, "atomP", worldLV, false, 0, fCheckOverlaps);
0358         }
0359         else if (AtomTemp->fElement.compare("O") == 0) {
0360           nbAtomO++;
0361           new G4PVPlacement(0,
0362                             G4ThreeVector(AtomTemp->fX * 1 * angstrom, AtomTemp->fY * 1 * angstrom,
0363                                           AtomTemp->fZ * 1 * angstrom),
0364                             atomLV_O, "atomP", worldLV, false, 0, fCheckOverlaps);
0365         }
0366         else if (AtomTemp->fElement.compare("N") == 0) {
0367           nbAtomN++;
0368           new G4PVPlacement(0,
0369                             G4ThreeVector(AtomTemp->fX * 1 * angstrom, AtomTemp->fY * 1 * angstrom,
0370                                           AtomTemp->fZ * 1 * angstrom),
0371                             atomLV_N, "atomP", worldLV, false, 0, fCheckOverlaps);
0372         }
0373         else if (AtomTemp->fElement.compare("S") == 0) {
0374           nbAtomS++;
0375           new G4PVPlacement(0,
0376                             G4ThreeVector(AtomTemp->fX * 1 * angstrom, AtomTemp->fY * 1 * angstrom,
0377                                           AtomTemp->fZ * 1 * angstrom),
0378                             atomLV_S, "atomP", worldLV, false, 0, fCheckOverlaps);
0379         }
0380         else if (AtomTemp->fElement.compare("P") == 0) {
0381           nbAtomP++;
0382           new G4PVPlacement(0,
0383                             G4ThreeVector(AtomTemp->fX * 1 * angstrom, AtomTemp->fY * 1 * angstrom,
0384                                           AtomTemp->fZ * 1 * angstrom),
0385                             atomLV_P, "atomP", worldLV, false, 0, fCheckOverlaps);
0386         }
0387         else {
0388           nbAtomX++;
0389           new G4PVPlacement(0,
0390                             G4ThreeVector(AtomTemp->fX * 1 * angstrom, AtomTemp->fY * 1 * angstrom,
0391                                           AtomTemp->fZ * 1 * angstrom),
0392                             atomLV_X, "atomP", worldLV, false, 0, fCheckOverlaps);
0393         }
0394 
0395         nbAtomTot++;
0396         //}//End if (i>=startFrom)
0397         AtomTemp = AtomTemp->GetNext();
0398       }  // end of for (  i=0 ; i < residueListTemp->nbAtom ; i++)
0399 
0400       residueListTemp = residueListTemp->GetNext();
0401     }  // end of while (residueListTemp)
0402 
0403     moleculeListTemp = moleculeListTemp->GetNext();
0404   }  // end of while (moleculeListTemp)
0405 
0406   G4cout << "**************** atomisticView(...) ****************" << G4endl;
0407   G4cout << "Number of loaded chains = " << k << G4endl;
0408   G4cout << "Number of Atoms      = " << nbAtomTot << G4endl;
0409   G4cout << "Number of Hydrogens  = " << nbAtomH << G4endl;
0410   G4cout << "Number of Carbons    = " << nbAtomC << G4endl;
0411   G4cout << "Number of Oxygens    = " << nbAtomO << G4endl;
0412   G4cout << "Number of Nitrogens  = " << nbAtomN << G4endl;
0413   G4cout << "Number of Sulfurs    = " << nbAtomS << G4endl;
0414   G4cout << "Number of Phosphorus = " << nbAtomP << G4endl;
0415   G4cout << "Number of undifined atoms =" << nbAtomX << G4endl << G4endl;
0416 }
0417 
0418 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0419 
0420 //////////////////////////////////////////////////
0421 ///////////////// BEGIN representation for nucleotide barycenters
0422 //
0423 void DetectorConstruction::BarycenterView(G4LogicalVolume* worldLV, Barycenter* barycenterListTemp)
0424 {
0425   CheckMaterials();
0426 
0427   G4VSolid* atomS_ZZ;
0428   G4LogicalVolume* atomLV_ZZ;
0429   G4VisAttributes* MyVisAtt_ZZ;
0430 
0431   while (barycenterListTemp) {
0432     atomS_ZZ = new G4Orb("Sphere", (barycenterListTemp->GetRadius()) * angstrom);
0433     atomLV_ZZ = new G4LogicalVolume(atomS_ZZ, fpWaterMaterial, "atomLV_ZZ");
0434     MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Magenta()));
0435     MyVisAtt_ZZ->SetForceSolid(true);
0436     atomLV_ZZ->SetVisAttributes(MyVisAtt_ZZ);
0437 
0438     new G4PVPlacement(0,
0439                       G4ThreeVector(barycenterListTemp->fCenterX * 1 * angstrom,
0440                                     barycenterListTemp->fCenterY * 1 * angstrom,
0441                                     barycenterListTemp->fCenterZ * 1 * angstrom),
0442                       atomLV_ZZ, "atomZZ", worldLV, false, 0, fCheckOverlaps);
0443 
0444     barycenterListTemp = barycenterListTemp->GetNext();
0445   }  // end of while (moleculeListTemp)
0446 }
0447 
0448 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0449 
0450 //////////////////////////////////////////////////
0451 ///////////////// BEGIN representation for Base Sugar Phosphate barycenters
0452 //
0453 void DetectorConstruction::ResiduesView(G4LogicalVolume* worldLV, Barycenter* barycenterListTemp)
0454 {
0455   CheckMaterials();
0456   G4VisAttributes* MyVisAtt_ZZ;
0457 
0458   G4VSolid* tubS1_ZZ;
0459   G4LogicalVolume* tubLV1_ZZ;
0460   G4VSolid* tubS2_ZZ;
0461   G4LogicalVolume* tubLV2_ZZ;
0462 
0463   G4VSolid* AS_ZZ;
0464   G4LogicalVolume* ALV_ZZ;
0465   G4VSolid* BS_ZZ;
0466   G4LogicalVolume* BLV_ZZ;
0467   G4VSolid* CS_ZZ;
0468   G4LogicalVolume* CLV_ZZ;
0469 
0470   while (barycenterListTemp) {
0471     // 3 spheres to Base, Sugar, Phosphate)
0472     AS_ZZ = new G4Orb("Sphere", 1. * angstrom);
0473     ALV_ZZ = new G4LogicalVolume(AS_ZZ, fpWaterMaterial, "ALV_ZZ");
0474     MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Blue()));
0475     MyVisAtt_ZZ->SetForceSolid(true);
0476     ALV_ZZ->SetVisAttributes(MyVisAtt_ZZ);
0477     new G4PVPlacement(0,
0478                       G4ThreeVector(barycenterListTemp->fCenterBaseX * angstrom,
0479                                     barycenterListTemp->fCenterBaseY * angstrom,
0480                                     barycenterListTemp->fCenterBaseZ * angstrom),
0481                       ALV_ZZ, "AZZ", worldLV, false, 0, fCheckOverlaps);
0482 
0483     BS_ZZ = new G4Orb("Sphere", 1. * angstrom);
0484     BLV_ZZ = new G4LogicalVolume(BS_ZZ, fpWaterMaterial, "BLV_ZZ");
0485     MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Red()));
0486     MyVisAtt_ZZ->SetForceSolid(true);
0487     BLV_ZZ->SetVisAttributes(MyVisAtt_ZZ);
0488     new G4PVPlacement(0,
0489                       G4ThreeVector((barycenterListTemp->fCenterPhosphateX) * angstrom,
0490                                     (barycenterListTemp->fCenterPhosphateY) * angstrom,
0491                                     (barycenterListTemp->fCenterPhosphateZ) * angstrom),
0492                       BLV_ZZ, "BZZ", worldLV, false, 0, fCheckOverlaps);
0493 
0494     CS_ZZ = new G4Orb("Sphere", 1. * angstrom);
0495     CLV_ZZ = new G4LogicalVolume(CS_ZZ, fpWaterMaterial, "CLV_ZZ");
0496     MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Yellow()));
0497     MyVisAtt_ZZ->SetForceSolid(true);
0498     CLV_ZZ->SetVisAttributes(MyVisAtt_ZZ);
0499     new G4PVPlacement(0,
0500                       G4ThreeVector(barycenterListTemp->fCenterSugarX * angstrom,
0501                                     barycenterListTemp->fCenterSugarY * angstrom,
0502                                     barycenterListTemp->fCenterSugarZ * angstrom),
0503                       CLV_ZZ, "CZZ", worldLV, false, 0, fCheckOverlaps);
0504 
0505     // 2 cylinders (Base<=>Sugar and Sugar<=>Phosphat)
0506     //  Cylinder Base<=>Sugar
0507     tubS1_ZZ = new G4Tubs(
0508       "Cylinder", 0., 0.5 * angstrom,
0509       std::sqrt((barycenterListTemp->fCenterBaseX - barycenterListTemp->fCenterSugarX)
0510                   * (barycenterListTemp->fCenterBaseX - barycenterListTemp->fCenterSugarX)
0511                 + (barycenterListTemp->fCenterBaseY - barycenterListTemp->fCenterSugarY)
0512                     * (barycenterListTemp->fCenterBaseY - barycenterListTemp->fCenterSugarY)
0513                 + (barycenterListTemp->fCenterBaseZ - barycenterListTemp->fCenterSugarZ)
0514                     * (barycenterListTemp->fCenterBaseZ - barycenterListTemp->fCenterSugarZ))
0515         / 2 * angstrom,
0516       0., 2. * pi);
0517 
0518     tubLV1_ZZ = new G4LogicalVolume(tubS1_ZZ, fpWaterMaterial, "tubLV_ZZ");
0519     MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Green()));
0520     MyVisAtt_ZZ->SetForceSolid(true);
0521     tubLV1_ZZ->SetVisAttributes(MyVisAtt_ZZ);
0522 
0523     G4double Ux = barycenterListTemp->fCenterBaseX - barycenterListTemp->fCenterSugarX;
0524     G4double Uy = barycenterListTemp->fCenterBaseY - barycenterListTemp->fCenterSugarY;
0525     G4double Uz = barycenterListTemp->fCenterBaseZ - barycenterListTemp->fCenterSugarZ;
0526     G4double llUll = std::sqrt(Ux * Ux + Uy * Uy + Uz * Uz);
0527 
0528     Ux = Ux / llUll;
0529     Uy = Uy / llUll;
0530     Uz = Uz / llUll;
0531 
0532     G4ThreeVector direction = G4ThreeVector(Ux, Uy, Uz);
0533     G4double theta_euler = direction.theta();
0534     G4double phi_euler = direction.phi();
0535     G4double psi_euler = 0;
0536 
0537     // Warning : clhep Euler constructor build inverse matrix !
0538     G4RotationMatrix rotm1Inv = G4RotationMatrix(phi_euler + pi / 2, theta_euler, psi_euler);
0539     G4RotationMatrix rotm1 = rotm1Inv.inverse();
0540     G4ThreeVector translm1 = G4ThreeVector(
0541       (barycenterListTemp->fCenterBaseX + barycenterListTemp->fCenterSugarX) / 2. * angstrom,
0542       (barycenterListTemp->fCenterBaseY + barycenterListTemp->fCenterSugarY) / 2. * angstrom,
0543       (barycenterListTemp->fCenterBaseZ + barycenterListTemp->fCenterSugarZ) / 2. * angstrom);
0544     G4Transform3D transform1 = G4Transform3D(rotm1, translm1);
0545     new G4PVPlacement(transform1,  // rotation translation
0546                       tubLV1_ZZ, "atomZZ", worldLV, false, 0, fCheckOverlaps);
0547 
0548     // Cylinder Sugar<=>Phosphat
0549     tubS2_ZZ = new G4Tubs(
0550       "Cylinder2", 0., 0.5 * angstrom,
0551       std::sqrt((barycenterListTemp->fCenterSugarX - barycenterListTemp->fCenterPhosphateX)
0552                   * (barycenterListTemp->fCenterSugarX - barycenterListTemp->fCenterPhosphateX)
0553                 + (barycenterListTemp->fCenterSugarY - barycenterListTemp->fCenterPhosphateY)
0554                     * (barycenterListTemp->fCenterSugarY - barycenterListTemp->fCenterPhosphateY)
0555                 + (barycenterListTemp->fCenterSugarZ - barycenterListTemp->fCenterPhosphateZ)
0556                     * (barycenterListTemp->fCenterSugarZ - barycenterListTemp->fCenterPhosphateZ))
0557         / 2 * angstrom,
0558       0., 2. * pi);
0559 
0560     tubLV2_ZZ = new G4LogicalVolume(tubS2_ZZ, fpWaterMaterial, "tubLV2_ZZ");
0561     MyVisAtt_ZZ = new G4VisAttributes(G4Colour(1, 0.5, 0));
0562     MyVisAtt_ZZ->SetForceSolid(true);
0563     tubLV2_ZZ->SetVisAttributes(MyVisAtt_ZZ);
0564 
0565     Ux = barycenterListTemp->fCenterSugarX - barycenterListTemp->fCenterPhosphateX;
0566     Uy = barycenterListTemp->fCenterSugarY - barycenterListTemp->fCenterPhosphateY;
0567     Uz = barycenterListTemp->fCenterSugarZ - barycenterListTemp->fCenterPhosphateZ;
0568     llUll = std::sqrt(Ux * Ux + Uy * Uy + Uz * Uz);
0569 
0570     Ux = Ux / llUll;
0571     Uy = Uy / llUll;
0572     Uz = Uz / llUll;
0573 
0574     direction = G4ThreeVector(Ux, Uy, Uz);
0575     theta_euler = direction.theta();
0576     phi_euler = direction.phi();
0577     psi_euler = 0;
0578 
0579     // Warning : clhep Euler constructor build inverse matrix !
0580     rotm1Inv = G4RotationMatrix(phi_euler + pi / 2, theta_euler, psi_euler);
0581     rotm1 = rotm1Inv.inverse();
0582     translm1 = G4ThreeVector(
0583       (barycenterListTemp->fCenterSugarX + barycenterListTemp->fCenterPhosphateX) / 2. * angstrom,
0584       (barycenterListTemp->fCenterSugarY + barycenterListTemp->fCenterPhosphateY) / 2. * angstrom,
0585       (barycenterListTemp->fCenterSugarZ + barycenterListTemp->fCenterPhosphateZ) / 2. * angstrom);
0586     transform1 = G4Transform3D(rotm1, translm1);
0587     new G4PVPlacement(transform1,  // rotation translation
0588                       tubLV2_ZZ, "atomZZ", worldLV, false, 0, fCheckOverlaps);
0589 
0590     barycenterListTemp = barycenterListTemp->GetNext();
0591   }  // end of while sur moleculeListTemp
0592 }
0593 
0594 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0595 
0596 //////////////////////////////////////////////////
0597 ///////////////// BEGIN draw a bounding volume
0598 //
0599 void DetectorConstruction::DrawBoundingVolume(G4LogicalVolume* worldLV, Molecule* moleculeListTemp)
0600 {
0601   CheckMaterials();
0602 
0603   double dX, dY, dZ;  // Dimensions for bounding volume
0604   double tX, tY, tZ;  // Translation for bounding volume
0605   fPDBlib.ComputeBoundingVolumeParams(moleculeListTemp, dX, dY, dZ, tX, tY, tZ);
0606 
0607   //[Geometry] Build a box :
0608   G4VSolid* boundingS =
0609     new G4Box("Bounding", dX * 1 * angstrom, dY * 1 * angstrom, dZ * 1 * angstrom);
0610 
0611   G4LogicalVolume* boundingLV = new G4LogicalVolume(boundingS, fpWaterMaterial, "BoundingLV");
0612 
0613   G4RotationMatrix* pRot = new G4RotationMatrix();
0614 
0615   G4VisAttributes* MyVisAtt_ZZ = new G4VisAttributes(G4Colour(G4Colour::Gray()));
0616   boundingLV->SetVisAttributes(MyVisAtt_ZZ);
0617 
0618   new G4PVPlacement(pRot,
0619                     G4ThreeVector(tX * 1 * angstrom, tY * 1 * angstrom,
0620                                   tZ * 1 * angstrom),  // at (x,y,z)
0621                     boundingLV, "boundingPV", worldLV, false, 0, fCheckOverlaps);
0622 }
0623 
0624 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0625 
0626 void DetectorConstruction::LoadPDBfile(G4String fileName)
0627 {
0628   G4cout << "Load PDB file : " << fileName << "." << G4endl << G4endl;
0629   fPdbFileName = fileName;
0630 #ifdef G4MULTITHREADED
0631   G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, fChosenOption));
0632   G4MTRunManager::GetRunManager()->ReinitializeGeometry();
0633 #else
0634   G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, fChosenOption));
0635 #endif
0636 }
0637 
0638 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0639 
0640 void DetectorConstruction::BuildBoundingVolume()
0641 {
0642   if (fPdbFileStatus > 0)  // a PDB file has been loaded
0643   {
0644     G4cout << "Build only world volume and bounding volume"
0645               " for computation."
0646            << G4endl << G4endl;
0647 #ifdef G4MULTITHREADED
0648     G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 10));
0649     G4MTRunManager::GetRunManager()->ReinitializeGeometry();
0650 #else
0651     G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 10));
0652 #endif
0653   }
0654   else
0655     G4cout << "PDB file not found!" << G4endl << G4endl;
0656 }
0657 
0658 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0659 
0660 void DetectorConstruction::DrawAtoms_()
0661 {
0662   if (fPdbFileStatus > 0)  // a PDB file has been loaded
0663   {
0664 #ifdef G4MULTITHREADED
0665     G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 1));
0666     G4MTRunManager::GetRunManager()->ReinitializeGeometry();
0667 #else
0668     G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 1));
0669 #endif
0670   }
0671   else
0672     G4cout << "PDB file not found!" << G4endl << G4endl;
0673 }
0674 
0675 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0676 
0677 void DetectorConstruction::DrawNucleotides_()
0678 {
0679   if (fPdbFileStatus > 0)  // a PDB file has been loaded
0680   {
0681 #ifdef G4MULTITHREADED
0682     G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 2));
0683     G4MTRunManager::GetRunManager()->ReinitializeGeometry();
0684 #else
0685     G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 2));
0686 #endif
0687   }
0688   else
0689     G4cout << "PDB file not found!" << G4endl << G4endl;
0690 }
0691 
0692 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0693 
0694 void DetectorConstruction::DrawResidues_()
0695 {
0696   if (fPdbFileStatus > 0)  // a PDB file has been loaded
0697   {
0698 #ifdef G4MULTITHREADED
0699     G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 3));
0700     G4MTRunManager::GetRunManager()->ReinitializeGeometry();
0701 #else
0702     G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 3));
0703 #endif
0704   }
0705   else
0706     G4cout << "PDB file not found!" << G4endl << G4endl;
0707 }
0708 
0709 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0710 
0711 void DetectorConstruction::DrawAtomsWithBounding_()
0712 {
0713   if (fPdbFileStatus > 0)  // a PDB file has been loaded
0714   {
0715 #ifdef G4MULTITHREADED
0716     G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 11));
0717     G4MTRunManager::GetRunManager()->ReinitializeGeometry();
0718 #else
0719     G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 11));
0720 #endif
0721   }
0722   else
0723     G4cout << "PDB file not found!" << G4endl << G4endl;
0724 }
0725 
0726 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0727 
0728 void DetectorConstruction::DrawNucleotidesWithBounding_()
0729 {
0730   if (fPdbFileStatus > 0)  // a PDB file has been loaded
0731   {
0732 #ifdef G4MULTITHREADED
0733     G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 12));
0734     G4MTRunManager::GetRunManager()->ReinitializeGeometry();
0735 #else
0736     G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 12));
0737 #endif
0738   }
0739   else
0740     G4cout << "PDB file not found!" << G4endl << G4endl;
0741 }
0742 
0743 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0744 
0745 void DetectorConstruction::DrawResiduesWithBounding_()
0746 {
0747   if (fPdbFileStatus > 0)  // a PDB file has been loaded
0748   {
0749 #ifdef G4MULTITHREADED
0750     G4MTRunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 13));
0751     G4MTRunManager::GetRunManager()->ReinitializeGeometry();
0752 #else
0753     G4RunManager::GetRunManager()->DefineWorldVolume(DefineVolumes(fPdbFileName, 13));
0754 #endif
0755   }
0756   else
0757     G4cout << "PDB file not found!" << G4endl << G4endl;
0758 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!