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