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