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