Back to home page

EIC code displayed by LXR
 
 

    

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