| 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 // Authors: S. Meylan and C. Villagrasa (IRSN, 26 // Authors: S. Meylan and C. Villagrasa (IRSN, France)
27 // Updated: H. Tran (IRSN), France: 20/12/2018 27 // Updated: H. Tran (IRSN), France: 20/12/2018
28 // 28 //
29 29
30 #include "DNAParser.hh" 30 #include "DNAParser.hh"
31 << 31 #include <fstream>
32 #include "G4Box.hh" << 32 #include "G4SystemOfUnits.hh"
33 #include "G4DNAChemistryManager.hh" 33 #include "G4DNAChemistryManager.hh"
>> 34 #include "G4VPhysicalVolume.hh"
34 #include "G4LogicalVolume.hh" 35 #include "G4LogicalVolume.hh"
35 #include "G4Material.hh" << 36 #include "G4VPhysicalVolume.hh"
>> 37 #include "G4VSolid.hh"
>> 38 #include "G4PVPlacement.hh"
36 #include "G4NistManager.hh" 39 #include "G4NistManager.hh"
>> 40 #include "G4Material.hh"
>> 41 #include "G4Box.hh"
37 #include "G4Orb.hh" 42 #include "G4Orb.hh"
38 #include "G4PVPlacement.hh" <<
39 #include "G4SubtractionSolid.hh" 43 #include "G4SubtractionSolid.hh"
40 #include "G4SystemOfUnits.hh" <<
41 #include "G4VPhysicalVolume.hh" <<
42 #include "G4VSolid.hh" <<
43 <<
44 #include <fstream> <<
45 44
46 //....oooOO0OOooo........oooOO0OOooo........oo 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
47 // Molecule struct def 46 // Molecule struct def
48 struct DNAParser::Molecule 47 struct DNAParser::Molecule
49 { 48 {
50 Molecule(std::string name, int copyNumber, << 49 Molecule(std::string name,
51 double waterRadius, const std::st << 50 int copyNumber,
52 : fName(name), << 51 const G4ThreeVector& position,
53 fMaterial(material), << 52 double radius,
54 fCopyNumber(copyNumber), << 53 double waterRadius,
55 fPosition(position), << 54 const std::string& material,
56 fRadius(radius), << 55 int strand)
57 fRadiusWater(waterRadius), << 56 : fName(name),
58 fStrand(strand) << 57 fMaterial(material),
>> 58 fCopyNumber(copyNumber),
>> 59 fPosition(position),
>> 60 fRadius(radius),
>> 61 fRadiusWater(waterRadius),
>> 62 fStrand(strand)
59 {} 63 {}
60 64
61 G4String fName; 65 G4String fName;
62 G4String fMaterial; 66 G4String fMaterial;
63 G4int fCopyNumber; 67 G4int fCopyNumber;
64 G4ThreeVector fPosition; 68 G4ThreeVector fPosition;
65 G4double fRadius; 69 G4double fRadius;
66 G4double fRadiusWater; 70 G4double fRadiusWater;
67 G4int fStrand; 71 G4int fStrand;
68 72
69 G4bool operator<(const Molecule& rhs) cons << 73 G4bool operator<(const Molecule& rhs) const
>> 74 {
>> 75 return (fPosition.z() < rhs.fPosition.z());
>> 76 }
70 }; 77 };
71 78
72 //....oooOO0OOooo........oooOO0OOooo........oo 79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73 DNAParser::DNAParser() : fSize(0), fGeoName("" << 80 DNAParser::DNAParser()
>> 81 : fSize(0)
>> 82 , fGeoName("")
>> 83 , fpWater(nullptr)
>> 84 , fpGun(new G4MoleculeGun())
74 { 85 {
75 EnumParser(); << 86 EnumParser();
76 } 87 }
77 88
78 //....oooOO0OOooo........oooOO0OOooo........oo 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
79 90
80 DNAParser::~DNAParser() = default; 91 DNAParser::~DNAParser() = default;
81 92
82 //....oooOO0OOooo........oooOO0OOooo........oo 93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83 94
84 G4LogicalVolume* DNAParser::CreateLogicVolume( 95 G4LogicalVolume* DNAParser::CreateLogicVolume()
85 { 96 {
86 G4NistManager* pMan = G4NistManager::Instanc << 97 G4NistManager * pMan = G4NistManager::Instance();
87 fpWater = pMan->FindOrBuildMaterial("G4_WATE << 98 fpWater = pMan->FindOrBuildMaterial("G4_WATER");
88 <<
89 G4String boxNameSolid = fGeoName + "_solid"; <<
90 <<
91 G4Box* pBoxSolid = new G4Box(boxNameSolid, f <<
92 <<
93 G4String boxNameLogic = fGeoName + "_logic"; <<
94 99
95 G4LogicalVolume* pBoxLogic = new G4LogicalVo << 100 G4String boxNameSolid = fGeoName + "_solid";
96 << 101
97 std::sort(fMolecules.begin(), fMolecules.end << 102 G4Box* pBoxSolid = new G4Box(boxNameSolid,
98 << 103 fSize/2,
99 for (int i = 0, ie = fMolecules.size(); i < << 104 fSize/2,
100 G4String name = fMolecules[i].fName; << 105 fSize/2);
101 G4double radius = fMolecules[i].fRadius; << 106
102 G4double waterRadius = fMolecules[i].fRadi << 107 G4String boxNameLogic = fGeoName + "_logic";
103 G4ThreeVector moleculePosition = fMolecule << 108
104 << 109 G4LogicalVolume* pBoxLogic = new G4LogicalVolume(pBoxSolid,
105 int copyNum = fMolecules[i].fCopyNumber; << 110 fpWater,
106 << 111 boxNameLogic);
107 G4Orb* pMoleculeWaterSolid = nullptr; << 112
108 G4VSolid* pMoleculeWaterCutSolid = nullptr << 113 std::sort( fMolecules.begin(), fMolecules.end() );
109 G4LogicalVolume* pMoleculeWaterLogic = nul << 114
110 G4VPhysicalVolume* pPhysicalVolumeWater = << 115 for(int i = 0, ie = fMolecules.size(); i < ie; ++i)
111 << 116 {
112 G4double tolerence = 1e-4 * nm; << 117 G4String name = fMolecules[i].fName;
113 << 118 G4double radius = fMolecules[i].fRadius;
114 if (waterRadius > tolerence) { << 119 G4double waterRadius = fMolecules[i].fRadiusWater;
115 G4String nameWaterSolid = name + "_water << 120 G4ThreeVector moleculePosition = fMolecules[i].fPosition;
116 pMoleculeWaterSolid = new G4Orb(nameWate << 121
117 pMoleculeWaterCutSolid = << 122 int copyNum = fMolecules[i].fCopyNumber;
118 CreateCutSolid(pMoleculeWaterSolid, fM << 123
119 << 124 G4Orb* pMoleculeWaterSolid = nullptr;
120 G4String nameWaterLogic = name + "_water << 125 G4VSolid* pMoleculeWaterCutSolid = nullptr;
121 << 126 G4LogicalVolume* pMoleculeWaterLogic = nullptr;
122 pMoleculeWaterLogic = new G4LogicalVolum << 127 G4VPhysicalVolume* pPhysicalVolumeWater = nullptr;
123 G4String nameWaterPhys = name + "_water" << 128
124 pPhysicalVolumeWater = new G4PVPlacement << 129 G4double tolerence = 1e-4 * nm;
125 << 130
126 << 131 if(waterRadius > tolerence)
127 auto it = fEnumMap.find(nameWaterPhys); << 132 {
128 if (it != fEnumMap.end()) { << 133 G4String nameWaterSolid = name + "_water_solid";
129 fGeometryMap.insert(std::make_pair(pPh << 134 pMoleculeWaterSolid = new G4Orb(nameWaterSolid, waterRadius);
130 } << 135 pMoleculeWaterCutSolid = CreateCutSolid(pMoleculeWaterSolid,
131 else { << 136 fMolecules[i],
132 G4ExceptionDescription exceptionDescri << 137 fMolecules,
133 exceptionDescription << nameWaterPhys << 138 false);
134 G4Exception("DNAParser::CreateLogicVol << 139
135 exceptionDescription); << 140 G4String nameWaterLogic = name + "_water_logic";
136 } << 141
137 } << 142 pMoleculeWaterLogic = new G4LogicalVolume(pMoleculeWaterCutSolid,
138 // Dna volume part << 143 fpWater,
>> 144 nameWaterLogic);
>> 145 G4String nameWaterPhys = name + "_water";
>> 146 pPhysicalVolumeWater = new G4PVPlacement(0,
>> 147 moleculePosition,
>> 148 pMoleculeWaterLogic,
>> 149 nameWaterPhys,
>> 150 pBoxLogic,
>> 151 false,
>> 152 copyNum);
>> 153
>> 154 auto it = fEnumMap.find(nameWaterPhys);
>> 155 if(it != fEnumMap.end())
>> 156 {
>> 157 fGeometryMap.insert(std::make_pair(pPhysicalVolumeWater,
>> 158 it->second));
>> 159 }
>> 160 else
>> 161 {
>> 162 G4ExceptionDescription exceptionDescription;
>> 163 exceptionDescription <<nameWaterPhys
>> 164 <<" could not be exsited";
>> 165 G4Exception("DNAParser::CreateLogicVolume()",
>> 166 "DNAParser001", FatalException,
>> 167 exceptionDescription);
>> 168 }
>> 169 }
>> 170 // Dna volume part
139 171
140 G4Orb* pMoleculeSolid = nullptr; << 172 G4Orb* pMoleculeSolid = nullptr;
141 G4VSolid* pMoleculeCutSolid = nullptr; << 173 G4VSolid* pMoleculeCutSolid = nullptr;
142 G4LogicalVolume* pMoleculeLogic = nullptr; << 174 G4LogicalVolume* pMoleculeLogic = nullptr;
143 G4VPhysicalVolume* pPhysicalVolumeMolecule << 175 G4VPhysicalVolume* pPhysicalVolumeMolecule = nullptr;
144 << 176
145 G4String nameSolid = fMolecules[i].fName + << 177 G4String nameSolid = fMolecules[i].fName + "_solid";
146 pMoleculeSolid = new G4Orb(nameSolid, radi << 178 pMoleculeSolid = new G4Orb(nameSolid, radius);
147 << 179
148 pMoleculeCutSolid = CreateCutSolid(pMolecu << 180 pMoleculeCutSolid = CreateCutSolid(pMoleculeSolid,
149 << 181 fMolecules[i],
150 G4String nameLogic = name + "_logic"; << 182 fMolecules,
151 << 183 true);
152 pMoleculeLogic = new G4LogicalVolume(pMole << 184
153 if (waterRadius > tolerence) { << 185 G4String nameLogic = name + "_logic";
154 G4ThreeVector position(0, 0, 0); << 186
155 std::string namePhys = name; << 187 pMoleculeLogic = new G4LogicalVolume(pMoleculeCutSolid,
156 pPhysicalVolumeMolecule = new G4PVPlacem << 188 fpWater,
157 << 189 nameLogic);
158 << 190 if(waterRadius > tolerence)
159 auto it = fEnumMap.find(namePhys); << 191 {
160 if (it != fEnumMap.end()) { << 192 G4ThreeVector position(0,0,0);
161 fGeometryMap.insert(std::make_pair(pPh << 193 std::string namePhys = name;
162 } << 194 pPhysicalVolumeMolecule = new G4PVPlacement(0,
163 else { << 195 position,
164 G4ExceptionDescription exceptionDescri << 196 pMoleculeLogic,
165 exceptionDescription << namePhys << " << 197 namePhys,
166 G4Exception("DNAParser::CreateLogicVol << 198 pMoleculeWaterLogic,
167 exceptionDescription); << 199 false,
168 } << 200 copyNum);
169 } << 201
170 else { << 202
171 G4ThreeVector position = moleculePositio << 203 auto it = fEnumMap.find(namePhys);
172 G4String namePhys = name; << 204 if(it != fEnumMap.end())
173 pPhysicalVolumeMolecule = << 205 {
174 new G4PVPlacement(0, position, pMolecu << 206 fGeometryMap.insert(std::make_pair(pPhysicalVolumeMolecule,
175 << 207 it->second));
176 auto it = fEnumMap.find(namePhys); << 208 }
177 << 209 else
178 if (it != fEnumMap.end()) { << 210 {
179 fGeometryMap.insert(std::make_pair(pPh << 211 G4ExceptionDescription exceptionDescription;
180 } << 212 exceptionDescription <<namePhys
181 else { << 213 <<" could not be exsited";
182 G4ExceptionDescription exceptionDescri << 214 G4Exception("DNAParser::CreateLogicVolume()",
183 exceptionDescription << namePhys << " << 215 "DNAParser002", FatalException,
184 G4Exception("DNAParser::CreateLogicVol << 216 exceptionDescription);
185 exceptionDescription); << 217 }
186 } << 218 }
>> 219 else
>> 220 {
>> 221 G4ThreeVector position = moleculePosition;
>> 222 G4String namePhys = name;
>> 223 pPhysicalVolumeMolecule = new G4PVPlacement(0,
>> 224 position,
>> 225 pMoleculeLogic,
>> 226 namePhys,
>> 227 pBoxLogic,
>> 228 false,
>> 229 copyNum);
>> 230
>> 231 auto it = fEnumMap.find(namePhys);
>> 232
>> 233 if(it != fEnumMap.end())
>> 234 {
>> 235 fGeometryMap.insert(std::make_pair(pPhysicalVolumeMolecule,
>> 236 it->second));
>> 237 }
>> 238 else
>> 239 {
>> 240 G4ExceptionDescription exceptionDescription;
>> 241 exceptionDescription <<namePhys
>> 242 <<" could not be exsited";
>> 243 G4Exception("DNAParser::CreateLogicVolume()",
>> 244 "DNAParser003", FatalException,
>> 245 exceptionDescription);
>> 246 }
>> 247 }
187 } 248 }
188 } << 249 fMolecules.clear();
189 fMolecules.clear(); << 250 fRadiusMap.clear();
190 fRadiusMap.clear(); << 251 fWaterRadiusMap.clear();
191 fWaterRadiusMap.clear(); << 252 return pBoxLogic;
192 return pBoxLogic; <<
193 } 253 }
194 254
195 //....oooOO0OOooo........oooOO0OOooo........oo 255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
196 256
197 void DNAParser::ParseFile(const std::string& f 257 void DNAParser::ParseFile(const std::string& fileName)
198 { 258 {
199 fMolecules.clear(); << 259 fMolecules.clear();
200 fRadiusMap.clear(); << 260 fRadiusMap.clear();
201 fWaterRadiusMap.clear(); << 261 fWaterRadiusMap.clear();
202 <<
203 std::ifstream file(fileName.c_str()); <<
204 <<
205 if (!file.is_open()) { <<
206 G4ExceptionDescription exceptionDescriptio <<
207 exceptionDescription << fileName << "could <<
208 G4Exception("DNAParser::ParseFile()", "DNA <<
209 } <<
210 std::string line; <<
211 while (std::getline(file, line)) { <<
212 if (line.empty()) continue; <<
213 std::istringstream issLine(line); <<
214 std::string firstItem; <<
215 issLine >> firstItem; <<
216 <<
217 if (firstItem == "_Name") { <<
218 G4String name = ""; <<
219 issLine >> name; <<
220 fGeoName = name; <<
221 } <<
222 if (firstItem == "_Size") { <<
223 G4double size; <<
224 issLine >> size; <<
225 size *= nm; <<
226 fSize = size; <<
227 } <<
228 if (firstItem == "_Radius") { <<
229 G4String name; <<
230 issLine >> name; <<
231 G4double radius; <<
232 issLine >> radius; <<
233 radius *= nm; <<
234 G4double waterRadius; <<
235 issLine >> waterRadius; <<
236 waterRadius *= nm; <<
237 fRadiusMap[name] = radius; <<
238 fWaterRadiusMap[name] = waterRadius; <<
239 } <<
240 if (firstItem == "_pl") { <<
241 std::string name; <<
242 issLine >> name; <<
243 G4String material; <<
244 issLine >> material; <<
245 <<
246 G4int strand; <<
247 issLine >> strand; <<
248 <<
249 G4int copyNumber; <<
250 issLine >> copyNumber; <<
251 <<
252 G4double x; <<
253 issLine >> x; <<
254 x *= nm; <<
255 <<
256 G4double y; <<
257 issLine >> y; <<
258 y *= nm; <<
259 <<
260 G4double z; <<
261 issLine >> z; <<
262 z *= nm; <<
263 <<
264 Molecule molecule(name, copyNumber, G4Th <<
265 fWaterRadiusMap[name], <<
266 fMolecules.push_back(molecule); <<
267 <<
268 auto itt = fEnumMap.find(name); <<
269 <<
270 if (itt != fEnumMap.end()) { <<
271 if (itt->second != DNAVolumeType::phos <<
272 if (itt->second == DNAVolumeType::de <<
273 || itt->second == DNAVolumeType: <<
274 { <<
275 name = "Deoxyribose"; <<
276 } <<
277 if (itt->second == DNAVolumeType::ba <<
278 name = "Adenine"; <<
279 } <<
280 if (itt->second == DNAVolumeType::ba <<
281 name = "Thymine"; <<
282 } <<
283 if (itt->second == DNAVolumeType::ba <<
284 name = "Cytosine"; <<
285 } <<
286 if (itt->second == DNAVolumeType::ba <<
287 name = "Guanine"; <<
288 } <<
289 if (itt->second == DNAVolumeType::hi <<
290 name = "Histone"; <<
291 } <<
292 262
293 fpGun->AddMolecule(name, G4ThreeVect << 263 std::ifstream file(fileName.c_str());
294 } << 264
295 } << 265 if(!file.is_open())
296 else { << 266 {
297 G4ExceptionDescription exceptionDescri 267 G4ExceptionDescription exceptionDescription;
298 exceptionDescription << name << " coul << 268 exceptionDescription <<fileName
299 G4Exception("DNAParser::ParseFile()", << 269 <<"could not be opened";
300 } << 270 G4Exception("DNAParser::ParseFile()",
>> 271 "DNAParser002", FatalException,
>> 272 exceptionDescription);
301 } 273 }
302 } << 274 std::string line;
303 file.close(); << 275 while(std::getline(file, line) )
>> 276 {
>> 277 if( line.empty() ) continue;
>> 278 std::istringstream issLine(line);
>> 279 std::string firstItem;
>> 280 issLine >> firstItem;
>> 281
>> 282 if(firstItem == "_Name")
>> 283 {
>> 284 G4String name = "";
>> 285 issLine >> name;
>> 286 fGeoName = name;
>> 287 }
>> 288 if(firstItem == "_Size")
>> 289 {
>> 290 G4double size;
>> 291 issLine >> size;
>> 292 size *= nm;
>> 293 fSize = size;
>> 294 }
>> 295 if(firstItem == "_Radius")
>> 296 {
>> 297 G4String name;
>> 298 issLine >> name;
>> 299 G4double radius;
>> 300 issLine >> radius;
>> 301 radius *= nm;
>> 302 G4double waterRadius;
>> 303 issLine >> waterRadius;
>> 304 waterRadius *= nm;
>> 305 fRadiusMap[name] = radius;
>> 306 fWaterRadiusMap[name] = waterRadius;
>> 307 }
>> 308 if(firstItem == "_pl")
>> 309 {
>> 310 std::string name;
>> 311 issLine >> name;
>> 312 G4String material;
>> 313 issLine >> material;
>> 314
>> 315 G4int strand;
>> 316 issLine >> strand;
>> 317
>> 318 G4int copyNumber;
>> 319 issLine >> copyNumber;
>> 320
>> 321 G4double x;
>> 322 issLine >> x;
>> 323 x *= nm;
>> 324
>> 325 G4double y;
>> 326 issLine >> y;
>> 327 y *= nm;
>> 328
>> 329 G4double z;
>> 330 issLine >> z;
>> 331 z *= nm;
>> 332
>> 333 Molecule molecule(name,
>> 334 copyNumber,
>> 335 G4ThreeVector(x, y, z),
>> 336 fRadiusMap[name],
>> 337 fWaterRadiusMap[name],
>> 338 material,
>> 339 strand);
>> 340 fMolecules.push_back(molecule);
>> 341
>> 342 auto itt = fEnumMap.find(name);
>> 343
>> 344 if(itt != fEnumMap.end())
>> 345 {
>> 346
>> 347 if(itt->second != DNAVolumeType::phosphate1 &&
>> 348 itt->second != DNAVolumeType::phosphate2)
>> 349 {
>> 350 if(itt->second == DNAVolumeType::deoxyribose1 ||
>> 351 itt->second == DNAVolumeType::deoxyribose2)
>> 352 {
>> 353 name = "Deoxyribose";
>> 354 }
>> 355 if(itt->second == DNAVolumeType::base_adenine)
>> 356 {
>> 357 name = "Adenine";
>> 358 }
>> 359 if(itt->second == DNAVolumeType::base_thymine)
>> 360 {
>> 361 name = "Thymine";
>> 362 }
>> 363 if(itt->second == DNAVolumeType::base_cytosine)
>> 364 {
>> 365 name = "Cytosine";
>> 366 }
>> 367 if(itt->second == DNAVolumeType::base_guanine)
>> 368 {
>> 369 name = "Guanine";
>> 370 }
>> 371 if(itt->second == DNAVolumeType::histone)
>> 372 {
>> 373 name = "Histone";
>> 374 }
>> 375
>> 376 fpGun->AddMolecule(name,
>> 377 G4ThreeVector(x, y, z),
>> 378 1*picosecond);
>> 379 }
>> 380 }
>> 381 else
>> 382 {
>> 383 G4ExceptionDescription exceptionDescription;
>> 384 exceptionDescription <<name
>> 385 <<" could not be exsited";
>> 386 G4Exception("DNAParser::ParseFile()",
>> 387 "DNAParser005", FatalException,
>> 388 exceptionDescription);
>> 389 }
>> 390 }
>> 391 }
>> 392 file.close();
304 } 393 }
305 394
306 //....oooOO0OOooo........oooOO0OOooo........oo 395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
307 396
308 G4VSolid* DNAParser::CreateCutSolid(G4Orb* pSo << 397 G4VSolid* DNAParser::CreateCutSolid(G4Orb *pSolidOrbRef,
309 std::vecto << 398 Molecule &molRef,
>> 399 std::vector<Molecule> &molList,
>> 400 G4bool in)
310 { 401 {
311 // The idea behing this method is to cut ove << 402 // The idea behing this method is to cut overlap volumes by selecting
312 // one of them (the reference) and checking << 403 //one of them (the reference) and checking all the other volumes (the targets).
313 // If a reference and a target volumes are c << 404 // If a reference and a target volumes are close enough to overlap they will be cut.
314 // The reference is already selected when we << 405 // The reference is already selected when we enter this method.
315 // Use the tiny space to differentiate the f << 406 // Use the tiny space to differentiate the frontiers (may not be necessary)
316 G4double tinySpace = 0.001 * nm; << 407 G4double tinySpace = 0.001*nm;
317 << 408
318 G4SubtractionSolid* pSolidCut(NULL); << 409 G4SubtractionSolid* pSolidCut(NULL);
319 G4bool isCutted = false; << 410 G4bool isCutted = false;
320 G4bool isOurVol = false; << 411 G4bool isOurVol = false;
321 G4double radiusRef; << 412 G4double radiusRef;
322 << 413
323 if (molRef.fRadiusWater == 0) { << 414 if(molRef.fRadiusWater == 0)
324 radiusRef = molRef.fRadius; << 415 {
325 } << 416 radiusRef = molRef.fRadius;
326 else { << 417 }
327 radiusRef = molRef.fRadiusWater; << 418 else
328 } << 419 {
329 << 420 radiusRef = molRef.fRadiusWater;
330 G4ThreeVector posRef = molRef.fPosition; << 421 }
331 << 422
332 if (std::abs(posRef.x()) + radiusRef > fSize << 423 G4ThreeVector posRef = molRef.fPosition;
333 || std::abs(posRef.y()) + radiusRef > fS << 424
334 || std::abs(posRef.z()) + radiusRef > fS << 425 if(std::abs(posRef.x() ) + radiusRef > fSize/2 // along x
335 { << 426 || std::abs(posRef.y() ) + radiusRef > fSize/2 // along y
336 G4Box* solidBox = new G4Box("solid_box_for << 427 || std::abs(posRef.z() ) + radiusRef > fSize/2) // along z
337 G4ThreeVector posBox; << 428 {
338 G4RotationMatrix rotMat; << 429 G4Box* solidBox = new G4Box("solid_box_for_cut",
339 rotMat.rotate(0, G4ThreeVector(0, 0, 1)); << 430 fSize/2,
340 << 431 fSize/2,
341 if (std::abs(posRef.y() + radiusRef) > fSi << 432 fSize/2);
342 posBox = -posRef + G4ThreeVector(0, fSiz << 433 G4ThreeVector posBox;
343 if (!isCutted) { << 434 G4RotationMatrix rotMat;
344 pSolidCut = new G4SubtractionSolid("so << 435 rotMat.rotate(0, G4ThreeVector(0,0,1) );
345 isCutted = true; << 436
346 } << 437 if(std::abs( posRef.y() + radiusRef ) > fSize/2 )
347 else { << 438 {
348 pSolidCut = new G4SubtractionSolid("so << 439 posBox = -posRef + G4ThreeVector(0,fSize,0);
349 } << 440 if(!isCutted)
350 } << 441 {
351 // Down << 442 pSolidCut =
352 if (std::abs(posRef.y() - radiusRef) > fSi << 443 new G4SubtractionSolid("solidCut",
353 posBox = -posRef + G4ThreeVector(0, -fSi << 444 pSolidOrbRef,
354 << 445 solidBox,
355 if (!isCutted) { << 446 &rotMat,
356 pSolidCut = new G4SubtractionSolid("so << 447 posBox);
357 isCutted = true; << 448 isCutted = true;
358 } << 449 }
359 else { << 450 else
360 pSolidCut = new G4SubtractionSolid("so << 451 {
361 } << 452 pSolidCut =
362 } << 453 new G4SubtractionSolid("solidCut",
363 // Left << 454 pSolidCut,
364 if (std::abs(posRef.x() + radiusRef) > fSi << 455 solidBox,
365 posBox = -posRef + G4ThreeVector(fSize, << 456 &rotMat,
366 if (!isCutted) { << 457 posBox);
367 pSolidCut = new G4SubtractionSolid("so << 458 }
368 isCutted = true; << 459 }
369 } << 460 // Down
370 else { << 461 if(std::abs( posRef.y() - radiusRef ) > fSize/2 )
371 pSolidCut = new G4SubtractionSolid("so << 462 {
372 } << 463 posBox = -posRef + G4ThreeVector(0,-fSize,0);
373 } << 464
>> 465 if(!isCutted)
>> 466 {
>> 467 pSolidCut = new G4SubtractionSolid("solidCut",
>> 468 pSolidOrbRef,
>> 469 solidBox,
>> 470 &rotMat,
>> 471 posBox);
>> 472 isCutted = true;
>> 473 }
>> 474 else
>> 475 {
>> 476 pSolidCut = new G4SubtractionSolid("solidCut",
>> 477 pSolidCut,
>> 478 solidBox,
>> 479 &rotMat,
>> 480 posBox);
>> 481 }
>> 482 }
>> 483 // Left
>> 484 if(std::abs( posRef.x() + radiusRef ) > fSize/2 )
>> 485 {
>> 486 posBox = -posRef + G4ThreeVector(fSize,0,0);
>> 487 if(!isCutted)
>> 488 {
>> 489 pSolidCut = new G4SubtractionSolid("solidCut",
>> 490 pSolidOrbRef,
>> 491 solidBox,
>> 492 &rotMat,
>> 493 posBox);
>> 494 isCutted = true;
>> 495 }
>> 496 else
>> 497 {
>> 498 pSolidCut = new G4SubtractionSolid("solidCut",
>> 499 pSolidCut,
>> 500 solidBox,
>> 501 &rotMat,
>> 502 posBox);
>> 503 }
>> 504 }
374 505
375 // Right << 506 // Right
376 if (std::abs(posRef.x() - radiusRef) > fSi << 507 if(std::abs( posRef.x() - radiusRef ) > fSize/2 )
377 posBox = -posRef + G4ThreeVector(-fSize, << 508 {
378 if (!isCutted) { << 509 posBox = -posRef + G4ThreeVector(-fSize,0,0);
379 pSolidCut = new G4SubtractionSolid("so << 510 if(!isCutted)
380 isCutted = true; << 511 {
381 } << 512 pSolidCut = new G4SubtractionSolid("solidCut",
382 else { << 513 pSolidOrbRef,
383 pSolidCut = new G4SubtractionSolid("so << 514 solidBox,
384 } << 515 &rotMat,
385 } << 516 posBox);
386 // Forward << 517 isCutted = true;
387 if (std::abs(posRef.z() + radiusRef) > fSi << 518 }
388 posBox = -posRef + G4ThreeVector(0, 0, f << 519 else
389 if (!isCutted) { << 520 {
390 pSolidCut = new G4SubtractionSolid("so << 521 pSolidCut = new G4SubtractionSolid("solidCut",
391 isCutted = true; << 522 pSolidCut,
392 } << 523 solidBox,
393 else { << 524 &rotMat,
394 pSolidCut = new G4SubtractionSolid("so << 525 posBox);
395 } << 526 }
396 } << 527 }
>> 528 // Forward
>> 529 if(std::abs( posRef.z() + radiusRef ) > fSize/2 )
>> 530 {
>> 531 posBox = -posRef + G4ThreeVector(0,0,fSize);
>> 532 if(!isCutted)
>> 533 {
>> 534 pSolidCut = new G4SubtractionSolid("solidCut",
>> 535 pSolidOrbRef,
>> 536 solidBox,
>> 537 &rotMat,
>> 538 posBox);
>> 539 isCutted = true;
>> 540 }
>> 541 else
>> 542 {
>> 543 pSolidCut = new G4SubtractionSolid("solidCut",
>> 544 pSolidCut,
>> 545 solidBox,
>> 546 &rotMat,
>> 547 posBox);
>> 548 }
>> 549 }
397 550
398 // Backward << 551 // Backward
399 if (std::abs(posRef.z() - radiusRef) > fSi << 552 if(std::abs( posRef.z() - radiusRef ) > fSize/2 )
400 posBox = -posRef + G4ThreeVector(0, 0, - << 553 {
401 if (!isCutted) { << 554 posBox = -posRef + G4ThreeVector(0,0,-fSize);
402 pSolidCut = new G4SubtractionSolid("so << 555 if(!isCutted)
403 isCutted = true; << 556 {
404 } << 557 pSolidCut = new G4SubtractionSolid("solidCut",
405 else { << 558 pSolidOrbRef,
406 pSolidCut = new G4SubtractionSolid("so << 559 solidBox,
407 } << 560 &rotMat,
>> 561 posBox);
>> 562 isCutted = true;
>> 563 }
>> 564 else
>> 565 {
>> 566 pSolidCut = new G4SubtractionSolid("solidCut",
>> 567 pSolidCut,
>> 568 solidBox,
>> 569 &rotMat,
>> 570 posBox);
>> 571 }
>> 572 }
408 } 573 }
409 } <<
410 <<
411 for (int i = 0, ie = molList.size(); i < ie; <<
412 G4ThreeVector posTar = molList[i].fPositio <<
413 G4double rTar = posRef.z(); <<
414 G4double zTar = posTar.z(); <<
415 574
416 if (zTar > rTar + 20 * nm) { << 575 for(int i=0, ie=molList.size(); i<ie; ++i)
417 break; << 576 {
418 } << 577 G4ThreeVector posTar = molList[i].fPosition;
419 else if (zTar < rTar - 20 * nm) { << 578 G4double rTar = posRef.z();
420 continue; << 579 G4double zTar = posTar.z();
421 } << 580
>> 581 if(zTar>rTar+20*nm)
>> 582 {
>> 583 break;
>> 584 }
>> 585 else if(zTar<rTar-20*nm)
>> 586 {
>> 587 continue;
>> 588 }
422 589
423 G4double radiusTar; << 590 G4double radiusTar;
424 if (molList[i].fRadiusWater == 0) { << 591 if(molList[i].fRadiusWater == 0)
425 radiusTar = molList[i].fRadius; << 592 {
426 } << 593 radiusTar = molList[i].fRadius;
427 else { << 594 }
428 radiusTar = molList[i].fRadiusWater; << 595 else
>> 596 {
>> 597 radiusTar = molList[i].fRadiusWater;
>> 598 }
>> 599
>> 600 G4double distance = std::abs( (posRef - posTar).getR() );
>> 601
>> 602 if(distance==0 && !isOurVol)
>> 603 {
>> 604 isOurVol = true;
>> 605 continue;
>> 606 }
>> 607 else if(distance == 0 && isOurVol)
>> 608 {
>> 609 G4ExceptionDescription exceptionDescription;
>> 610 exceptionDescription <<"CreateCutSolid: Two volumes "
>> 611 <<"are placed at the same position.";
>> 612 G4Exception("DNAParser::CreateCutSolid()",
>> 613 "DNAParser004", FatalException,
>> 614 exceptionDescription);
>> 615 }
>> 616 else if(distance <= radiusRef+radiusTar)
>> 617 {
>> 618 G4Box* solidBox = new G4Box("solid_box_for_cut",
>> 619 2*radiusTar,
>> 620 2*radiusTar,
>> 621 2*radiusTar);
>> 622 G4ThreeVector diff = posTar - posRef;
>> 623 G4double d = (std::pow(radiusRef,2) -
>> 624 std::pow(radiusTar,2) +
>> 625 std::pow(distance,2) )/(2*distance) +
>> 626 solidBox->GetZHalfLength() -
>> 627 tinySpace;
>> 628 if(in)
>> 629 {
>> 630 d -= 2*tinySpace;
>> 631 }
>> 632 G4ThreeVector pos = d *( diff/diff.getR() );
>> 633 G4double phi = std::acos(pos.getZ()/pos.getR());
>> 634 G4double theta = std::acos( pos.getX() /
>> 635 ( pos.getR() *
>> 636 std::cos(CLHEP::pi/2. - phi) ) );
>> 637
>> 638 if(pos.getY()<0)
>> 639 {
>> 640 theta = -theta;
>> 641 }
>> 642
>> 643 G4ThreeVector rotAxisForPhi(1*nm,0.,0.);
>> 644 rotAxisForPhi.rotateZ(theta + CLHEP::pi/2.);
>> 645
>> 646 G4RotationMatrix *rotMat = new G4RotationMatrix;
>> 647 rotMat->rotate(-phi, rotAxisForPhi);
>> 648
>> 649 G4ThreeVector rotZAxis(0.,0.,1*nm);
>> 650 rotMat->rotate(theta, rotZAxis);
>> 651
>> 652 if(!isCutted)
>> 653 {
>> 654 pSolidCut = new G4SubtractionSolid("solidCut",
>> 655 pSolidOrbRef,
>> 656 solidBox,
>> 657 rotMat,
>> 658 pos);
>> 659 }
>> 660 else
>> 661 {
>> 662 pSolidCut = new G4SubtractionSolid("solidCut",
>> 663 pSolidCut,
>> 664 solidBox,
>> 665 rotMat,
>> 666 pos);
>> 667 }
>> 668 isCutted = true;
>> 669 }
429 } 670 }
430 671
431 G4double distance = std::abs((posRef - pos << 672 if(isCutted)
432 << 673 {
433 if (distance == 0 && !isOurVol) { << 674 return pSolidCut;
434 isOurVol = true; << 675 }
435 continue; << 676 else
436 } << 677 {
437 else if (distance == 0 && isOurVol) { << 678 return pSolidOrbRef;
438 G4ExceptionDescription exceptionDescript <<
439 exceptionDescription << "CreateCutSolid: <<
440 << "are placed at t <<
441 G4Exception("DNAParser::CreateCutSolid() <<
442 exceptionDescription); <<
443 } 679 }
444 else if (distance <= radiusRef + radiusTar <<
445 G4Box* solidBox = new G4Box("solid_box_f <<
446 G4ThreeVector diff = posTar - posRef; <<
447 G4double d = <<
448 (std::pow(radiusRef, 2) - std::pow(rad <<
449 + solidBox->GetZHalfLength() - tinySpa <<
450 if (in) { <<
451 d -= 2 * tinySpace; <<
452 } <<
453 G4ThreeVector pos = d * (diff / diff.get <<
454 G4double phi = std::acos(pos.getZ() / po <<
455 G4double theta = std::acos(pos.getX() / <<
456 <<
457 if (pos.getY() < 0) { <<
458 theta = -theta; <<
459 } <<
460 <<
461 G4ThreeVector rotAxisForPhi(1 * nm, 0., <<
462 rotAxisForPhi.rotateZ(theta + CLHEP::pi <<
463 <<
464 G4RotationMatrix* rotMat = new G4Rotatio <<
465 rotMat->rotate(-phi, rotAxisForPhi); <<
466 <<
467 G4ThreeVector rotZAxis(0., 0., 1 * nm); <<
468 rotMat->rotate(theta, rotZAxis); <<
469 <<
470 if (!isCutted) { <<
471 pSolidCut = new G4SubtractionSolid("so <<
472 } <<
473 else { <<
474 pSolidCut = new G4SubtractionSolid("so <<
475 } <<
476 isCutted = true; <<
477 } <<
478 } <<
479 <<
480 if (isCutted) { <<
481 return pSolidCut; <<
482 } <<
483 else { <<
484 return pSolidOrbRef; <<
485 } <<
486 } 680 }
487 681
488 //....oooOO0OOooo........oooOO0OOooo........oo 682 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
489 683
490 void DNAParser::EnumParser() 684 void DNAParser::EnumParser()
491 { << 685 {
492 fEnumMap["deoxyribose1"] = deoxyribose1; << 686 fEnumMap["deoxyribose1"] = deoxyribose1;
493 fEnumMap["phosphate1"] = phosphate1; << 687 fEnumMap["phosphate1"] = phosphate1;
494 fEnumMap["deoxyribose2"] = deoxyribose2; << 688 fEnumMap["deoxyribose2"] = deoxyribose2;
495 fEnumMap["phosphate2"] = phosphate2; << 689 fEnumMap["phosphate2"] = phosphate2;
496 << 690
497 fEnumMap["base_cytosine"] = base_cytosine; << 691 fEnumMap["base_cytosine"] = base_cytosine;
498 fEnumMap["base_guanine"] = base_guanine; << 692 fEnumMap["base_guanine"] = base_guanine;
499 fEnumMap["base_thymine"] = base_thymine; << 693 fEnumMap["base_thymine"] = base_thymine;
500 fEnumMap["base_adenine"] = base_adenine; << 694 fEnumMap["base_adenine"] = base_adenine;
501 << 695
502 fEnumMap["deoxyribose1_water"] = deoxyribose << 696 fEnumMap["deoxyribose1_water"] = deoxyribose1_water;
503 fEnumMap["phosphate1_water"] = phosphate1_wa << 697 fEnumMap["phosphate1_water"] = phosphate1_water;
504 fEnumMap["deoxyribose2_water"] = deoxyribose << 698 fEnumMap["deoxyribose2_water"] = deoxyribose2_water;
505 fEnumMap["phosphate2_water"] = phosphate2_wa << 699 fEnumMap["phosphate2_water"] = phosphate2_water;
506 << 700
507 fEnumMap["base_adenine_water"] = base_adenin << 701
508 fEnumMap["base_guanine_water"] = base_guanin << 702 fEnumMap["base_adenine_water"] = base_adenine_water;
509 fEnumMap["base_cytosine_water"] = base_cytos << 703 fEnumMap["base_guanine_water"] = base_guanine_water;
510 fEnumMap["base_thymine_water"] = base_thymin << 704 fEnumMap["base_cytosine_water"] = base_cytosine_water;
511 << 705 fEnumMap["base_thymine_water"] = base_thymine_water;
512 fEnumMap["histone"] = histone; << 706
513 fEnumMap["physWorld"] = physWorld; << 707 fEnumMap["histone"] = histone;
514 fEnumMap["VoxelStraight"] = VoxelStraight; << 708 fEnumMap["physWorld"] = physWorld;
>> 709 fEnumMap["VoxelStraight"] = VoxelStraight;
515 } 710 }
>> 711
516 712