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