| 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