26911512 |
1 | |
2 | // Flugg tag |
3 | // modified 17/06/02: by I. Gonzalez. STL migration |
4 | |
36c70081 |
5 | //#include <stdio.h> |
6 | //#include <iomanip.h> |
26911512 |
7 | #include "FGeometryInit.hh" |
8 | #include "FluggNavigator.hh" |
9 | #include "WrapUtils.hh" |
26d97e06 |
10 | #include "FlukaMaterial.hh" |
11 | #include "FlukaCompound.hh" |
26911512 |
12 | |
13 | FGeometryInit * FGeometryInit::flagInstance=0; |
14 | |
15 | FGeometryInit* FGeometryInit::GetInstance() { |
e73e0522 |
16 | #ifdef G4GEOMETRY_DEBUG |
17 | G4cout << "==> Flugg::FGeometryInit::GetInstance(), instance=" |
18 | << flagInstance << G4endl; |
19 | #endif |
26911512 |
20 | if (!flagInstance) |
21 | flagInstance = new FGeometryInit(); |
22 | |
e73e0522 |
23 | #ifdef G4GEOMETRY_DEBUG |
24 | G4cout << "<== Flugg::FGeometryInit::GetInstance(), instance=" |
25 | << flagInstance << G4endl; |
26 | #endif |
26911512 |
27 | return flagInstance; |
28 | } |
29 | |
30 | |
31 | FGeometryInit::FGeometryInit(): |
32 | fDetector(0), |
33 | fFieldManager(0), |
36c70081 |
34 | fTransportationManager(G4TransportationManager::GetTransportationManager()), |
26911512 |
35 | myTopNode(0), |
36 | ptrGeoMan(0), |
37 | ptrArray(0), |
38 | ptrTouchHist(0), |
39 | ptrOldNavHist(0), |
40 | ptrTempNavHist(0), |
36c70081 |
41 | ptrJrLtGeant(0) |
42 | { |
26911512 |
43 | |
e73e0522 |
44 | #ifdef G4GEOMETRY_DEBUG |
45 | G4cout << "==> Flugg FGeometryInit::FGeometryInit()" << G4endl; |
46 | G4cout << "\t+ Changing the G4Navigator for FluggNavigator..." << G4endl; |
47 | #endif |
26911512 |
48 | G4Navigator* actualnav = fTransportationManager->GetNavigatorForTracking(); |
49 | if (actualnav) { |
50 | FluggNavigator* newnav = new FluggNavigator(); |
51 | fTransportationManager->SetNavigatorForTracking(newnav); |
52 | } |
53 | else { |
36c70081 |
54 | G4cerr << "ERROR: Could not find the actual G4Navigator" << G4endl; |
26911512 |
55 | abort(); |
56 | } |
57 | |
58 | |
e73e0522 |
59 | #ifdef G4GEOMETRY_DEBUG |
60 | G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl; |
61 | #endif |
26911512 |
62 | |
63 | } |
64 | |
65 | |
66 | FGeometryInit::~FGeometryInit() { |
e73e0522 |
67 | G4cout << "==> Flugg FGeometryInit::~FGeometryInit()" << G4endl; |
26911512 |
68 | DeleteHistories(); |
69 | ptrGeoMan->OpenGeometry(); |
70 | if (fTransportationManager) |
71 | delete fTransportationManager; |
72 | if (ptrJrLtGeant) |
73 | delete[] ptrJrLtGeant; |
74 | DelHistArray(); |
75 | |
76 | //keep ATTENTION: never delete a pointer twice! |
e73e0522 |
77 | G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl; |
26911512 |
78 | } |
79 | |
80 | |
81 | void FGeometryInit::closeGeometry() { |
e73e0522 |
82 | #ifdef G4GEOMETRY_DEBUG |
83 | G4cout << "==> Flugg FGeometryInit::closeGeometry()" << G4endl; |
84 | #endif |
26911512 |
85 | |
86 | ptrGeoMan = G4GeometryManager::GetInstance(); |
87 | if (ptrGeoMan) { |
88 | ptrGeoMan->OpenGeometry(); |
89 | |
90 | //true argoment allows voxel construction; if false voxels are built |
91 | //only for replicated volumes |
92 | ptrGeoMan->CloseGeometry(true); |
93 | } |
94 | else { |
95 | G4cerr << "ERROR in FLUGG: Could not get G4GeometryManager instance" |
96 | << G4endl; |
97 | G4cerr << " in FGeometryInit::closeGeometry. Exiting!!!" |
98 | << G4endl; |
99 | } |
100 | |
e73e0522 |
101 | #ifdef G4GEOMETRY_DEBUG |
102 | G4cout << "<== Flugg FGeometryInit::closeGeometry()" << G4endl; |
103 | #endif |
26911512 |
104 | } |
105 | |
106 | //************************************************************************* |
107 | |
108 | void FGeometryInit::InitHistArray() { |
e73e0522 |
109 | #ifdef G4GEOMETRY_DEBUG |
110 | G4cout << "==> Flugg FGeometryInit::InitHistArray()" << G4endl; |
111 | #endif |
26911512 |
112 | ptrArray = new G4int[1000000]; |
113 | for(G4int i=0;i<1000000;i++) |
114 | ptrArray[i]=0; |
e73e0522 |
115 | #ifdef G4GEOMETRY_DEBUG |
116 | G4cout << "<== Flugg FGeometryInit::InitHistArray()" << G4endl; |
117 | #endif |
26911512 |
118 | } |
119 | |
120 | |
121 | |
122 | //************************************************************************* |
123 | //jrLtGeant stores all crossed lattice volume histories. |
124 | |
125 | void FGeometryInit::InitJrLtGeantArray() { |
26911512 |
126 | #ifdef G4GEOMETRY_DEBUG |
e73e0522 |
127 | G4cout << "==> Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl; |
26911512 |
128 | G4cout << "Initializing JrLtGeant array" << G4endl; |
129 | #endif |
130 | ptrJrLtGeant = new G4int[10000]; |
131 | for(G4int x=0;x<10000;x++) |
132 | ptrJrLtGeant[x]=-1; |
133 | flagLttcGeant = -1; |
e73e0522 |
134 | #ifdef G4GEOMETRY_DEBUG |
135 | G4cout << "<== Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl; |
136 | #endif |
26911512 |
137 | } |
138 | |
139 | |
140 | void FGeometryInit::SetLttcFlagGeant(G4int newFlagLttc) { |
e73e0522 |
141 | #ifdef G4GEOMETRY_DEBUG |
142 | G4cout << "==> Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl; |
143 | #endif |
26911512 |
144 | // Added by A.Solodkov |
145 | if (newFlagLttc >= 10000) { |
146 | G4cout << "Problems in FGeometryInit::SetLttcFlagGeant" << G4endl; |
147 | G4cout << "Index newFlagLttc=" << newFlagLttc << " is outside array bounds" |
148 | << G4endl; |
149 | G4cout << "Better to stop immediately !" << G4endl; |
150 | exit(1); |
151 | } |
152 | flagLttcGeant = newFlagLttc; |
e73e0522 |
153 | #ifdef G4GEOMETRY_DEBUG |
154 | G4cout << "<== Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl; |
155 | #endif |
26911512 |
156 | } |
157 | |
158 | void FGeometryInit::PrintJrLtGeant() { |
159 | #ifdef G4GEOMETRY_DEBUG |
160 | //G4cout << "jrLtGeant:" << G4endl; |
161 | //for(G4int y=0;y<=flagLttcGeant;y++) |
162 | // |
163 | // G4cout << "jrLtGeant[" << y << "]=" << ptrJrLtGeant[y] << G4endl; |
164 | #endif |
165 | } |
166 | |
167 | //************************************************************************** |
168 | |
169 | void FGeometryInit::PrintHistories() { |
170 | /* |
171 | #ifdef G4GEOMETRY_DEBUG |
172 | G4cout << "Touch hist:" << G4endl; |
173 | G4cout << *(ptrTouchHist->GetHistory()) << G4endl; |
174 | G4cout << "Tmp hist:" << G4endl; |
175 | G4cout << *(ptrTempNavHist->GetHistory()) << G4endl; |
176 | G4cout << "Old hist:" << G4endl; |
177 | G4cout << *(ptrOldNavHist->GetHistory()) << G4endl; |
178 | #endif |
179 | */ |
180 | } |
181 | |
182 | |
183 | |
184 | void FGeometryInit::InitHistories() { |
e73e0522 |
185 | #ifdef G4GEOMETRY_DEBUG |
186 | G4cout << "==> Flugg FGeometryInit::InitHistories()" << G4endl; |
187 | #endif |
26911512 |
188 | //init utility histories with navigator history |
189 | ptrTouchHist = |
190 | fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory(); |
191 | ptrTempNavHist = |
192 | fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory(); |
193 | ptrOldNavHist = new G4TouchableHistory(); |
e73e0522 |
194 | #ifdef G4GEOMETRY_DEBUG |
195 | G4cout << "<== Flugg FGeometryInit::InitHistories()" << G4endl; |
196 | #endif |
26911512 |
197 | } |
198 | |
199 | void FGeometryInit::DeleteHistories() { |
e73e0522 |
200 | #ifdef G4GEOMETRY_DEBUG |
201 | G4cout << "==> Flugg FGeometryInit::DeleteHistories()" << G4endl; |
202 | #endif |
203 | |
26911512 |
204 | delete ptrTouchHist; |
205 | delete ptrOldNavHist; |
206 | delete ptrTempNavHist; |
207 | |
208 | #ifdef G4GEOMETRY_DEBUG |
209 | G4cout << "Deleting step-history objects at end of run!" << G4endl; |
e73e0522 |
210 | G4cout << "<== Flugg FGeometryInit::DeleteHistories()" << G4endl; |
26911512 |
211 | #endif |
26911512 |
212 | } |
213 | |
214 | |
215 | void FGeometryInit::UpdateHistories(const G4NavigationHistory * history, |
216 | G4int flagHist) { |
e73e0522 |
217 | #ifdef G4GEOMETRY_DEBUG |
218 | G4cout << "==> Flugg FGeometryInit::UpdateHistories()" << G4endl; |
219 | #endif |
26911512 |
220 | PrintHistories(); |
221 | |
222 | #ifdef G4GEOMETRY_DEBUG |
223 | G4cout << "...updating histories!" << G4endl; |
224 | #endif |
225 | |
226 | G4VPhysicalVolume * pPhysVol = history->GetTopVolume(); |
227 | |
228 | switch (flagHist) { |
229 | case 0: { |
230 | //this is the case when a new history is given to the |
231 | //navigator and old history has to be resetted |
232 | //touchable history has not been updated jet, so: |
233 | |
234 | ptrTouchHist->UpdateYourself(pPhysVol,history); |
235 | ptrTempNavHist->UpdateYourself(pPhysVol,history); |
236 | G4NavigationHistory * ptrOldNavHistNotConst = |
237 | const_cast<G4NavigationHistory * >(ptrOldNavHist->GetHistory()); |
238 | ptrOldNavHistNotConst->Reset(); |
239 | ptrOldNavHistNotConst->Clear(); |
240 | PrintHistories(); |
241 | break; |
242 | } //case 0 |
243 | |
244 | case 1: { |
245 | //this is the case when a new history is given to the |
246 | //navigator but old history has to be kept (e.g. LOOKZ |
247 | //is call during an event); |
248 | //touchable history has not been updated jet, so: |
249 | |
250 | ptrTouchHist->UpdateYourself(pPhysVol,history); |
251 | ptrTempNavHist->UpdateYourself(pPhysVol,history); |
252 | PrintHistories(); |
253 | break; |
254 | } //case 1 |
255 | |
256 | case 2: { |
257 | //this is the case when the touchable history has been |
258 | //updated by a LocateGlobalPointAndUpdateTouchable call |
259 | |
260 | G4VPhysicalVolume * pPhysVolTemp = ptrTempNavHist->GetVolume(); |
261 | ptrOldNavHist->UpdateYourself(pPhysVolTemp, |
262 | ptrTempNavHist->GetHistory()); |
263 | |
264 | ptrTempNavHist->UpdateYourself(pPhysVol,history); |
265 | PrintHistories(); |
266 | break; |
267 | } //case 2 |
268 | |
269 | default: { |
270 | G4cout <<" ERROR in updating step-histories!" << G4endl; |
271 | break; |
272 | } //default |
273 | } //switch |
274 | |
e73e0522 |
275 | #ifdef G4GEOMETRY_DEBUG |
276 | G4cout << "<== Flugg FGeometryInit::UpdateHistories()" << G4endl; |
277 | #endif |
26911512 |
278 | } |
279 | |
280 | //***************************************************************************** |
281 | |
26911512 |
282 | void FGeometryInit::createFlukaMatFile() { |
26911512 |
283 | // last modification Sara Vanini 1/III/99 |
284 | // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case, |
285 | // according to the fluka standard. In addition,. they must be equal to the |
286 | // names of the fluka materials - see fluka manual - in order that the |
287 | // program load the right cross sections, and equal to the names included in |
288 | // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his |
289 | // own .pemf, in order to get the right cross sections loaded in memory. |
290 | |
291 | #ifdef G4GEOMETRY_DEBUG |
e73e0522 |
292 | G4cout << "==> Flugg FGeometryInit::createFlukaMatFile()" << G4endl; |
26911512 |
293 | G4cout << "================== FILEWR =================" << G4endl; |
294 | #endif |
26d97e06 |
295 | |
296 | |
297 | //Regions map |
298 | BuildRegionsMap(); |
299 | G4std::ofstream vos("Volumes_index.inp"); |
300 | PrintRegionsMap(vos); |
301 | vos.close(); |
302 | |
303 | //Materials and compounds |
304 | BuildMaterialTables(); |
305 | G4std::ofstream fos("flukaMat.inp"); |
306 | PrintMaterialTables(fos); |
307 | PrintAssignmat(fos); |
308 | PrintMagneticField(fos); |
309 | fos.close(); |
310 | |
311 | #ifdef G4GEOMETRY_DEBUG |
312 | G4cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << G4endl; |
313 | #endif |
314 | } |
315 | |
316 | //////////////////////////////////////////////////////////////////////// |
317 | // |
318 | void FGeometryInit::BuildRegionsMap() { |
319 | #ifdef G4GEOMETRY_DEBUG |
320 | G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl; |
321 | #endif |
322 | |
323 | //Find number of Volumes in physical volume store |
26911512 |
324 | G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance(); |
26d97e06 |
325 | unsigned int numVol = pVolStore->size(); |
326 | |
327 | G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore |
328 | << ") has " << numVol << " volumes. Iterating..." |
329 | << G4endl; |
330 | |
331 | for(unsigned int l=0; l < numVol; l++) { |
332 | //Get each of the physical volumes |
333 | G4VPhysicalVolume * physicalvolume = (*pVolStore)[l]; |
334 | G4int iFlukaRegion = l+1; |
335 | fRegionVolumeMap[physicalvolume] = iFlukaRegion; |
26911512 |
336 | } |
26d97e06 |
337 | |
338 | |
339 | |
26911512 |
340 | #ifdef G4GEOMETRY_DEBUG |
26d97e06 |
341 | G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl; |
342 | #endif |
343 | } |
26911512 |
344 | |
26d97e06 |
345 | void FGeometryInit::PrintRegionsMap(G4std::ostream& os) { |
26911512 |
346 | #ifdef G4GEOMETRY_DEBUG |
26d97e06 |
347 | G4cout << "==> Flugg FGeometryInit::PrintRegionsMap()" << G4endl; |
348 | #endif |
26911512 |
349 | |
26d97e06 |
350 | //Print some header |
351 | PrintHeader(os, "GEANT4 VOLUMES"); |
26911512 |
352 | |
26d97e06 |
353 | //Iterate over all volumes in the map |
354 | for (RegionIterator i = fRegionVolumeMap.begin(); |
355 | i != fRegionVolumeMap.end(); |
356 | i++) { |
26911512 |
357 | |
26d97e06 |
358 | //Get info in the map |
359 | G4VPhysicalVolume* ptrVol = (*i).first; |
360 | int index = (*i).second; |
26911512 |
361 | |
26d97e06 |
362 | //Print index and region name in some fixed format |
363 | os.setf(G4std::ios::left, G4std::ios::adjustfield); |
364 | os << setw10 << index; |
365 | os << G4std::setw(20) << ptrVol->GetName() << G4std::setw(20) << ""; |
366 | |
367 | //If volume is a replica... print some more stuff |
26911512 |
368 | if(ptrVol->IsReplicated()) { |
369 | EAxis axis; |
26d97e06 |
370 | G4int nRep = -1; |
371 | G4double width = -1; |
372 | G4double offset = -1; |
373 | G4bool consum = false; |
374 | ptrVol->GetReplicationData(axis, nRep, width, offset, consum); |
375 | os.setf(G4std::ios::left, G4std::ios::adjustfield); |
376 | os << setw10 << "Repetion Nb: " << G4std::setw(3) << nRep; |
26911512 |
377 | } |
26d97e06 |
378 | os << G4endl; |
26911512 |
379 | |
26d97e06 |
380 | } |
381 | |
382 | #ifdef G4GEOMETRY_DEBUG |
383 | G4cout << "<== Flugg FGeometryInit::PrintRegionsMap()" << G4endl; |
384 | #endif |
385 | } |
386 | |
bf547b2f |
387 | //////////////////////////////////////////////////////////////////////// |
388 | // |
389 | G4int FGeometryInit::GetRegionFromName(const char* volName) const { |
390 | for (RegionIterator i = fRegionVolumeMap.begin(); |
391 | i != fRegionVolumeMap.end(); |
392 | i++) { |
393 | |
394 | //Get info in the map |
395 | G4VPhysicalVolume* ptrVol = (*i).first; |
396 | if (ptrVol->GetName() == volName) |
397 | return ((*i).second); |
398 | } |
399 | return -1; |
400 | } |
401 | |
402 | |
403 | |
26d97e06 |
404 | //////////////////////////////////////////////////////////////////////// |
405 | // |
406 | void FGeometryInit::BuildMaterialTables() { |
407 | #ifdef G4GEOMETRY_DEBUG |
408 | G4cout << "==> Flugg FGeometryInit::BuildMaterialTables()" << G4endl; |
409 | #endif |
410 | |
411 | //some terminal printout also |
412 | G4cout << "\t* Storing information..." << G4endl; |
413 | |
414 | //The logic is the folloing: |
415 | //Get the Material Table and: |
416 | // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum |
417 | // 2) For each single element material build a material equivalent |
418 | // 3) For the rest: |
419 | // 3.a) Build materials for each not already known element |
420 | // 3.b) Build the compound out of them |
421 | |
422 | //Get the Material Table and iterate |
423 | const G4MaterialTable* matTable = G4Material::GetMaterialTable(); |
424 | for (MatTableIterator i = matTable->begin(); i != matTable->end(); i++) { |
425 | |
426 | //Get some basic material information |
427 | G4Material* material = (*i); |
428 | G4String matName = material->GetName(); |
429 | const G4double matDensity = material->GetDensity(); |
430 | const G4int nMatElements = material->GetNumberOfElements(); |
431 | |
432 | G4cout << "\t\t+ " << matName |
433 | << ": dens. = " << matDensity/(g/cm3) << "g/cm3" |
434 | << ", nElem = " << nMatElements << G4endl; |
435 | |
436 | // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum |
437 | // FlukaMaterial* is 0 in that case |
438 | if (matDensity <= 1.00e-10*g/cm3) { |
439 | G4FlukaMaterialMap[material] = 0; |
440 | G4cout << "\t\t Stored as vacuum" << G4endl; |
26911512 |
441 | } |
26d97e06 |
442 | // 2) For each single element material build a material equivalent |
443 | else if (nMatElements == 1) { |
26911512 |
444 | |
26d97e06 |
445 | FlukaMaterial *flukamat = |
446 | BuildFlukaMaterialFromElement(material->GetElement(0), |
447 | matDensity); |
26911512 |
448 | |
26d97e06 |
449 | G4FlukaMaterialMap[material] = flukamat; |
450 | G4cout << "\t\t Stored as " << flukamat->GetRealName() << G4endl; |
26911512 |
451 | |
26d97e06 |
452 | } //else if (material->GetNumberOfElements() == 1) |
453 | |
454 | // 3) For the rest: |
455 | // 3.a) Build materials for each not already known element |
456 | // 3.b) Build the compound out of them |
457 | else { |
458 | FlukaCompound* flukacomp = |
459 | BuildFlukaCompoundFromMaterial(material); |
460 | G4FlukaCompoundMap[material] = flukacomp; |
461 | G4cout << "\t\t Stored as " << flukacomp->GetRealName() << G4endl; |
462 | } //else for case 3) |
463 | } //for (materials) |
464 | |
465 | #ifdef G4GEOMETRY_DEBUG |
466 | G4cout << "<== Flugg FGeometryInit::BuildMaterialTables()" << G4endl; |
467 | #endif |
468 | } |
469 | |
470 | FlukaMaterial* |
471 | FGeometryInit::BuildFlukaMaterialFromElement(const G4Element* element, |
472 | G4double matDensity) { |
473 | #ifdef G4GEOMETRY_DEBUG |
474 | G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromElement()" |
475 | << G4endl; |
476 | #endif |
477 | |
478 | //Get element and its properties |
479 | G4String elemName(ToFlukaString(element->GetName())); |
480 | |
481 | FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(elemName); |
482 | if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) { |
483 | //Check for isotopes |
484 | G4int nIsotopes = element->GetNumberOfIsotopes(); |
485 | if (nIsotopes == 0) { |
486 | G4double elemA = element->GetA()/g; |
487 | G4double elemZ = element->GetZ(); |
26911512 |
488 | |
26d97e06 |
489 | if (elemA != G4int(elemA) && elemZ != G4int(elemZ)) { |
490 | G4cout << "WARNING: Element \'" << elemName |
491 | << "\' has non integer Z (" << elemZ << ") or A (" |
492 | << elemA << ")" |
493 | << G4endl; |
26911512 |
494 | } |
26d97e06 |
495 | |
496 | flukamat = new FlukaMaterial(elemName, |
497 | G4int(elemZ), |
498 | elemA, |
499 | matDensity/(g/cm3)); |
500 | } |
501 | else if (nIsotopes == 1) { |
502 | const G4Isotope* isotope = element->GetIsotope(0); |
503 | flukamat = BuildFlukaMaterialFromIsotope(isotope, matDensity); |
504 | } |
505 | else { |
506 | FlukaCompound *flucomp = BuildFlukaCompoundFromElement(element, |
507 | matDensity); |
508 | flukamat = flucomp->GetFlukaMaterial(); |
26911512 |
509 | } |
26d97e06 |
510 | } |
511 | #ifdef G4GEOMETRY_DEBUG |
512 | else { |
513 | G4cout << "INFO: Element \'" << elemName |
514 | << "\' already exists in the DB. It will not be recreated." |
515 | << G4endl; |
516 | } |
517 | #endif |
518 | |
519 | return flukamat; |
26911512 |
520 | |
26d97e06 |
521 | #ifdef G4GEOMETRY_DEBUG |
522 | G4cout << "<== Flugg FGeometryInit::BuildFlukaMaterialFromElement()" |
523 | << G4endl; |
524 | #endif |
525 | } |
526 | |
527 | FlukaMaterial* |
528 | FGeometryInit::BuildFlukaMaterialFromIsotope(const G4Isotope* isotope, |
529 | G4double matDensity) { |
530 | #ifdef G4GEOMETRY_DEBUG |
531 | G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()" |
532 | << G4endl; |
533 | #endif |
534 | G4String isoName(ToFlukaString(isotope->GetName())); |
535 | FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(isoName); |
536 | if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) { |
537 | G4int isoZ = isotope->GetZ(); |
538 | G4double isoA = (isotope->GetA())/(g); |
539 | G4int isoN = isotope->GetN(); |
540 | flukamat = new FlukaMaterial(isoName, |
541 | isoZ, |
542 | isoA, |
543 | matDensity/(g/cm3), |
544 | isoN); |
545 | } |
546 | |
547 | return flukamat; |
548 | |
549 | #ifdef G4GEOMETRY_DEBUG |
550 | G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()" |
551 | << G4endl; |
552 | #endif |
553 | } |
554 | |
555 | FlukaCompound* |
556 | FGeometryInit::BuildFlukaCompoundFromMaterial(const G4Material* material) { |
557 | #ifdef G4GEOMETRY_DEBUG |
558 | G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()" |
559 | << G4endl; |
560 | #endif |
561 | //Material properties |
562 | const G4double* elemFractions = material->GetFractionVector(); |
563 | const G4int nMatElements = material->GetNumberOfElements(); |
564 | const G4double matDensity = material->GetDensity(); |
565 | G4String matName(ToFlukaString(material->GetName())); |
566 | FlukaCompound* flukacomp = new FlukaCompound(matName, matDensity/(g/cm3), |
567 | nMatElements); |
568 | for (G4int i = 0; i < nMatElements; i++) { |
569 | FlukaMaterial *flukamat = |
570 | BuildFlukaMaterialFromElement(material->GetElement(i), 0.0); |
571 | |
572 | flukacomp->AddElement(flukamat->GetIndex(), -elemFractions[i]); |
573 | |
574 | } //for (elements) |
575 | |
576 | return flukacomp; |
577 | |
578 | #ifdef G4GEOMETRY_DEBUG |
579 | G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()" |
580 | << G4endl; |
581 | #endif |
582 | } |
583 | |
584 | FlukaCompound* |
585 | FGeometryInit::BuildFlukaCompoundFromElement(const G4Element* element, |
586 | G4double matDensity) { |
587 | #ifdef G4GEOMETRY_DEBUG |
588 | G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromElement()" |
589 | << G4endl; |
590 | #endif |
591 | G4int nIsotopes = element->GetNumberOfIsotopes(); |
592 | //fraction of nb of atomes per volume (= volume fraction?) |
593 | const G4double* isoAbundance = element->GetRelativeAbundanceVector(); |
594 | G4String elemName(ToFlukaString(element->GetName())); |
595 | |
596 | //Material properties |
597 | FlukaCompound* flukacomp = new FlukaCompound(elemName, matDensity/(g/cm3), |
598 | nIsotopes); |
599 | for (G4int i = 0; i < nIsotopes; i++) { |
600 | FlukaMaterial *flukamat = |
601 | BuildFlukaMaterialFromIsotope(element->GetIsotope(i), 0.0); |
602 | |
603 | flukacomp->AddElement(flukamat->GetIndex(), isoAbundance[i]); |
604 | |
605 | } //for (elements) |
606 | |
607 | return flukacomp; |
608 | |
609 | #ifdef G4GEOMETRY_DEBUG |
610 | G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromElement()" |
611 | << G4endl; |
612 | #endif |
613 | } |
614 | |
615 | void FGeometryInit::PrintMaterialTables(G4std::ostream& os) { |
616 | #ifdef G4GEOMETRY_DEBUG |
617 | G4cout << "==> Flugg FGeometryInit::PrintMaterialTables()" << G4endl; |
618 | #endif |
619 | //Print Header |
620 | PrintHeader(os, "GEANT4 MATERIALS AND COMPOUNDS"); |
26911512 |
621 | |
26d97e06 |
622 | //And some more stuff |
623 | size_t nIsotopes = G4Isotope::GetNumberOfIsotopes(); |
624 | size_t nElements = G4Element::GetNumberOfElements(); |
625 | size_t nMaterials = G4Material::GetNumberOfMaterials(); |
626 | |
0e22711e |
627 | os << "* In Geant4 there are " << nMaterials << " materials" << G4endl; |
628 | os << "* In Geant4 there are " << nElements << " elements" << G4endl; |
629 | os << "* In Geant4 there are " << nIsotopes << " isotopes" << G4endl; |
26d97e06 |
630 | |
631 | //Materials |
632 | G4cout << "\t* Printing FLUKA materials..." << G4endl; |
633 | FlukaMaterial::PrintMaterialsByIndex(os); |
634 | //FlukaMaterial::PrintMaterialsByName(os); |
635 | |
636 | //Compounds |
637 | G4cout << "\t* Printing FLUKA compounds..." << G4endl; |
638 | FlukaCompound::PrintCompounds(os); |
639 | |
640 | #ifdef G4GEOMETRY_DEBUG |
641 | G4cout << "<== Flugg FGeometryInit::PrintMaterialTables()" << G4endl; |
642 | #endif |
643 | } |
644 | |
645 | //////////////////////////////////////////////////////////////////////// |
646 | // |
647 | void FGeometryInit::PrintAssignmat(G4std::ostream& os) { |
648 | #ifdef G4GEOMETRY_DEBUG |
649 | G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl; |
650 | #endif |
651 | |
652 | //Find number of Volumes in physical volume store |
653 | G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance(); |
654 | unsigned int numVol = pVolStore->size(); |
655 | |
656 | G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore |
657 | << ") has " << numVol << " volumes. " << G4endl; |
658 | |
659 | G4cout << "\t* Printing ASSIGNMAT..." << G4endl; |
660 | |
661 | |
662 | PrintHeader(os,"GEANT4 MATERIAL ASSIGNMENTS"); |
663 | for(unsigned int l=0; l < numVol; l++) { |
664 | |
665 | //Get each of the physical volumes |
666 | G4VPhysicalVolume * physicalvol = (*pVolStore)[l]; |
667 | |
668 | //Get index for that volume |
669 | G4int iFlukaRegion = fRegionVolumeMap[physicalvol]; |
670 | |
671 | //Find G4 material and navigate to its fluka compound/material |
672 | G4LogicalVolume * logicalVol = physicalvol->GetLogicalVolume(); |
673 | G4Material* material = logicalVol->GetMaterial(); |
674 | G4int matIndex = 2; |
675 | if (G4FlukaCompoundMap[material]) |
676 | matIndex = G4FlukaCompoundMap[material]->GetIndex(); |
677 | if (G4FlukaMaterialMap[material]) |
678 | matIndex = G4FlukaMaterialMap[material]->GetIndex(); |
679 | |
680 | //Find if there is a magnetic field in the region |
681 | //check if Magnetic Field is present in the region |
682 | G4double flagField = 0.0; |
683 | G4FieldManager * pMagFieldMan = logicalVol->GetFieldManager(); |
684 | if(pMagFieldMan && pMagFieldMan->GetDetectorField()) |
685 | flagField = 1.0; |
686 | |
687 | //Print card |
688 | os << setw10 << "ASSIGNMAT "; |
689 | os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield); |
690 | os << setw10 << setfixed << G4double(matIndex); |
691 | os << setw10 << setfixed << G4double(iFlukaRegion); |
692 | os << setw10 << "0.0"; |
693 | os << setw10 << setfixed << flagField; |
694 | os << G4endl; |
695 | } |
696 | |
697 | |
698 | |
699 | #ifdef G4GEOMETRY_DEBUG |
700 | G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl; |
701 | #endif |
702 | } |
703 | |
704 | |
705 | void FGeometryInit::PrintMagneticField(G4std::ostream& os) { |
706 | #ifdef G4GEOMETRY_DEBUG |
707 | G4cout << "==> Flugg FGeometryInit::PrintMagneticField()" << G4endl; |
708 | #endif |
709 | |
710 | G4cout << "\t* Printing Magnetic Field..." << G4endl; |
711 | |
26911512 |
712 | if(fTransportationManager->GetFieldManager()->DoesFieldExist()) { |
26911512 |
713 | |
714 | //get magnetic field pointer |
715 | const G4Field * pMagField = |
716 | fTransportationManager->GetFieldManager()->GetDetectorField(); |
717 | |
26911512 |
718 | |
719 | if(pMagField) { |
26d97e06 |
720 | //Check if it can be made a uniform magnetic field |
26911512 |
721 | const G4UniformMagField *pUnifMagField = |
722 | dynamic_cast<const G4UniformMagField*>(pMagField); |
723 | if(pUnifMagField) { |
26d97e06 |
724 | G4double B[3]; |
725 | G4double point[4]; //it is not really used |
726 | pUnifMagField->GetFieldValue(point,B); |
727 | |
728 | //write MGNFIELD card |
729 | PrintHeader(os,"GEANT4 MAGNETIC FIELD"); |
730 | os << setw10 << "MGNFIELD "; |
731 | os << setw10 << ""; |
732 | os << setw10 << ""; |
733 | os << setw10 << ""; |
734 | os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield); |
735 | os << setw10 << setfixed |
736 | << G4std::setprecision(4) << B[0] |
737 | << setw10 << B[1] |
738 | << setw10 << B[2] |
739 | << G4endl; |
740 | } |
741 | else { |
742 | G4cout << "WARNING: No Uniform Magnetic Field found." << G4endl; |
743 | G4cout << " Manual intervention might be needed." << G4endl; |
26911512 |
744 | } |
26911512 |
745 | } |
26d97e06 |
746 | else |
747 | G4cout << "\t No detector field found... " << G4endl; |
26911512 |
748 | } // end if magnetic field |
26d97e06 |
749 | else |
750 | G4cout << "\t No field found... " << G4endl; |
751 | |
e73e0522 |
752 | #ifdef G4GEOMETRY_DEBUG |
26d97e06 |
753 | G4cout << "<== Flugg FGeometryInit::PrintMagneticField()" << G4endl; |
e73e0522 |
754 | #endif |
26911512 |
755 | } |