]>
Commit | Line | Data |
---|---|---|
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 | ||
387 | //////////////////////////////////////////////////////////////////////// | |
388 | // | |
389 | void FGeometryInit::BuildMaterialTables() { | |
390 | #ifdef G4GEOMETRY_DEBUG | |
391 | G4cout << "==> Flugg FGeometryInit::BuildMaterialTables()" << G4endl; | |
392 | #endif | |
393 | ||
394 | //some terminal printout also | |
395 | G4cout << "\t* Storing information..." << G4endl; | |
396 | ||
397 | //The logic is the folloing: | |
398 | //Get the Material Table and: | |
399 | // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum | |
400 | // 2) For each single element material build a material equivalent | |
401 | // 3) For the rest: | |
402 | // 3.a) Build materials for each not already known element | |
403 | // 3.b) Build the compound out of them | |
404 | ||
405 | //Get the Material Table and iterate | |
406 | const G4MaterialTable* matTable = G4Material::GetMaterialTable(); | |
407 | for (MatTableIterator i = matTable->begin(); i != matTable->end(); i++) { | |
408 | ||
409 | //Get some basic material information | |
410 | G4Material* material = (*i); | |
411 | G4String matName = material->GetName(); | |
412 | const G4double matDensity = material->GetDensity(); | |
413 | const G4int nMatElements = material->GetNumberOfElements(); | |
414 | ||
415 | G4cout << "\t\t+ " << matName | |
416 | << ": dens. = " << matDensity/(g/cm3) << "g/cm3" | |
417 | << ", nElem = " << nMatElements << G4endl; | |
418 | ||
419 | // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum | |
420 | // FlukaMaterial* is 0 in that case | |
421 | if (matDensity <= 1.00e-10*g/cm3) { | |
422 | G4FlukaMaterialMap[material] = 0; | |
423 | G4cout << "\t\t Stored as vacuum" << G4endl; | |
26911512 | 424 | } |
26d97e06 | 425 | // 2) For each single element material build a material equivalent |
426 | else if (nMatElements == 1) { | |
26911512 | 427 | |
26d97e06 | 428 | FlukaMaterial *flukamat = |
429 | BuildFlukaMaterialFromElement(material->GetElement(0), | |
430 | matDensity); | |
26911512 | 431 | |
26d97e06 | 432 | G4FlukaMaterialMap[material] = flukamat; |
433 | G4cout << "\t\t Stored as " << flukamat->GetRealName() << G4endl; | |
26911512 | 434 | |
26d97e06 | 435 | } //else if (material->GetNumberOfElements() == 1) |
436 | ||
437 | // 3) For the rest: | |
438 | // 3.a) Build materials for each not already known element | |
439 | // 3.b) Build the compound out of them | |
440 | else { | |
441 | FlukaCompound* flukacomp = | |
442 | BuildFlukaCompoundFromMaterial(material); | |
443 | G4FlukaCompoundMap[material] = flukacomp; | |
444 | G4cout << "\t\t Stored as " << flukacomp->GetRealName() << G4endl; | |
445 | } //else for case 3) | |
446 | } //for (materials) | |
447 | ||
448 | #ifdef G4GEOMETRY_DEBUG | |
449 | G4cout << "<== Flugg FGeometryInit::BuildMaterialTables()" << G4endl; | |
450 | #endif | |
451 | } | |
452 | ||
453 | FlukaMaterial* | |
454 | FGeometryInit::BuildFlukaMaterialFromElement(const G4Element* element, | |
455 | G4double matDensity) { | |
456 | #ifdef G4GEOMETRY_DEBUG | |
457 | G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromElement()" | |
458 | << G4endl; | |
459 | #endif | |
460 | ||
461 | //Get element and its properties | |
462 | G4String elemName(ToFlukaString(element->GetName())); | |
463 | ||
464 | FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(elemName); | |
465 | if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) { | |
466 | //Check for isotopes | |
467 | G4int nIsotopes = element->GetNumberOfIsotopes(); | |
468 | if (nIsotopes == 0) { | |
469 | G4double elemA = element->GetA()/g; | |
470 | G4double elemZ = element->GetZ(); | |
26911512 | 471 | |
26d97e06 | 472 | if (elemA != G4int(elemA) && elemZ != G4int(elemZ)) { |
473 | G4cout << "WARNING: Element \'" << elemName | |
474 | << "\' has non integer Z (" << elemZ << ") or A (" | |
475 | << elemA << ")" | |
476 | << G4endl; | |
26911512 | 477 | } |
26d97e06 | 478 | |
479 | flukamat = new FlukaMaterial(elemName, | |
480 | G4int(elemZ), | |
481 | elemA, | |
482 | matDensity/(g/cm3)); | |
483 | } | |
484 | else if (nIsotopes == 1) { | |
485 | const G4Isotope* isotope = element->GetIsotope(0); | |
486 | flukamat = BuildFlukaMaterialFromIsotope(isotope, matDensity); | |
487 | } | |
488 | else { | |
489 | FlukaCompound *flucomp = BuildFlukaCompoundFromElement(element, | |
490 | matDensity); | |
491 | flukamat = flucomp->GetFlukaMaterial(); | |
26911512 | 492 | } |
26d97e06 | 493 | } |
494 | #ifdef G4GEOMETRY_DEBUG | |
495 | else { | |
496 | G4cout << "INFO: Element \'" << elemName | |
497 | << "\' already exists in the DB. It will not be recreated." | |
498 | << G4endl; | |
499 | } | |
500 | #endif | |
501 | ||
502 | return flukamat; | |
26911512 | 503 | |
26d97e06 | 504 | #ifdef G4GEOMETRY_DEBUG |
505 | G4cout << "<== Flugg FGeometryInit::BuildFlukaMaterialFromElement()" | |
506 | << G4endl; | |
507 | #endif | |
508 | } | |
509 | ||
510 | FlukaMaterial* | |
511 | FGeometryInit::BuildFlukaMaterialFromIsotope(const G4Isotope* isotope, | |
512 | G4double matDensity) { | |
513 | #ifdef G4GEOMETRY_DEBUG | |
514 | G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()" | |
515 | << G4endl; | |
516 | #endif | |
517 | G4String isoName(ToFlukaString(isotope->GetName())); | |
518 | FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(isoName); | |
519 | if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) { | |
520 | G4int isoZ = isotope->GetZ(); | |
521 | G4double isoA = (isotope->GetA())/(g); | |
522 | G4int isoN = isotope->GetN(); | |
523 | flukamat = new FlukaMaterial(isoName, | |
524 | isoZ, | |
525 | isoA, | |
526 | matDensity/(g/cm3), | |
527 | isoN); | |
528 | } | |
529 | ||
530 | return flukamat; | |
531 | ||
532 | #ifdef G4GEOMETRY_DEBUG | |
533 | G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()" | |
534 | << G4endl; | |
535 | #endif | |
536 | } | |
537 | ||
538 | FlukaCompound* | |
539 | FGeometryInit::BuildFlukaCompoundFromMaterial(const G4Material* material) { | |
540 | #ifdef G4GEOMETRY_DEBUG | |
541 | G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()" | |
542 | << G4endl; | |
543 | #endif | |
544 | //Material properties | |
545 | const G4double* elemFractions = material->GetFractionVector(); | |
546 | const G4int nMatElements = material->GetNumberOfElements(); | |
547 | const G4double matDensity = material->GetDensity(); | |
548 | G4String matName(ToFlukaString(material->GetName())); | |
549 | FlukaCompound* flukacomp = new FlukaCompound(matName, matDensity/(g/cm3), | |
550 | nMatElements); | |
551 | for (G4int i = 0; i < nMatElements; i++) { | |
552 | FlukaMaterial *flukamat = | |
553 | BuildFlukaMaterialFromElement(material->GetElement(i), 0.0); | |
554 | ||
555 | flukacomp->AddElement(flukamat->GetIndex(), -elemFractions[i]); | |
556 | ||
557 | } //for (elements) | |
558 | ||
559 | return flukacomp; | |
560 | ||
561 | #ifdef G4GEOMETRY_DEBUG | |
562 | G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()" | |
563 | << G4endl; | |
564 | #endif | |
565 | } | |
566 | ||
567 | FlukaCompound* | |
568 | FGeometryInit::BuildFlukaCompoundFromElement(const G4Element* element, | |
569 | G4double matDensity) { | |
570 | #ifdef G4GEOMETRY_DEBUG | |
571 | G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromElement()" | |
572 | << G4endl; | |
573 | #endif | |
574 | G4int nIsotopes = element->GetNumberOfIsotopes(); | |
575 | //fraction of nb of atomes per volume (= volume fraction?) | |
576 | const G4double* isoAbundance = element->GetRelativeAbundanceVector(); | |
577 | G4String elemName(ToFlukaString(element->GetName())); | |
578 | ||
579 | //Material properties | |
580 | FlukaCompound* flukacomp = new FlukaCompound(elemName, matDensity/(g/cm3), | |
581 | nIsotopes); | |
582 | for (G4int i = 0; i < nIsotopes; i++) { | |
583 | FlukaMaterial *flukamat = | |
584 | BuildFlukaMaterialFromIsotope(element->GetIsotope(i), 0.0); | |
585 | ||
586 | flukacomp->AddElement(flukamat->GetIndex(), isoAbundance[i]); | |
587 | ||
588 | } //for (elements) | |
589 | ||
590 | return flukacomp; | |
591 | ||
592 | #ifdef G4GEOMETRY_DEBUG | |
593 | G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromElement()" | |
594 | << G4endl; | |
595 | #endif | |
596 | } | |
597 | ||
598 | void FGeometryInit::PrintMaterialTables(G4std::ostream& os) { | |
599 | #ifdef G4GEOMETRY_DEBUG | |
600 | G4cout << "==> Flugg FGeometryInit::PrintMaterialTables()" << G4endl; | |
601 | #endif | |
602 | //Print Header | |
603 | PrintHeader(os, "GEANT4 MATERIALS AND COMPOUNDS"); | |
26911512 | 604 | |
26d97e06 | 605 | //And some more stuff |
606 | size_t nIsotopes = G4Isotope::GetNumberOfIsotopes(); | |
607 | size_t nElements = G4Element::GetNumberOfElements(); | |
608 | size_t nMaterials = G4Material::GetNumberOfMaterials(); | |
609 | ||
610 | os << "* In Geant4 there are " << nMaterials << " materials" << endl; | |
611 | os << "* In Geant4 there are " << nElements << " elements" << endl; | |
612 | os << "* In Geant4 there are " << nIsotopes << " isotopes" << endl; | |
613 | ||
614 | //Materials | |
615 | G4cout << "\t* Printing FLUKA materials..." << G4endl; | |
616 | FlukaMaterial::PrintMaterialsByIndex(os); | |
617 | //FlukaMaterial::PrintMaterialsByName(os); | |
618 | ||
619 | //Compounds | |
620 | G4cout << "\t* Printing FLUKA compounds..." << G4endl; | |
621 | FlukaCompound::PrintCompounds(os); | |
622 | ||
623 | #ifdef G4GEOMETRY_DEBUG | |
624 | G4cout << "<== Flugg FGeometryInit::PrintMaterialTables()" << G4endl; | |
625 | #endif | |
626 | } | |
627 | ||
628 | //////////////////////////////////////////////////////////////////////// | |
629 | // | |
630 | void FGeometryInit::PrintAssignmat(G4std::ostream& os) { | |
631 | #ifdef G4GEOMETRY_DEBUG | |
632 | G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl; | |
633 | #endif | |
634 | ||
635 | //Find number of Volumes in physical volume store | |
636 | G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance(); | |
637 | unsigned int numVol = pVolStore->size(); | |
638 | ||
639 | G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore | |
640 | << ") has " << numVol << " volumes. " << G4endl; | |
641 | ||
642 | G4cout << "\t* Printing ASSIGNMAT..." << G4endl; | |
643 | ||
644 | ||
645 | PrintHeader(os,"GEANT4 MATERIAL ASSIGNMENTS"); | |
646 | for(unsigned int l=0; l < numVol; l++) { | |
647 | ||
648 | //Get each of the physical volumes | |
649 | G4VPhysicalVolume * physicalvol = (*pVolStore)[l]; | |
650 | ||
651 | //Get index for that volume | |
652 | G4int iFlukaRegion = fRegionVolumeMap[physicalvol]; | |
653 | ||
654 | //Find G4 material and navigate to its fluka compound/material | |
655 | G4LogicalVolume * logicalVol = physicalvol->GetLogicalVolume(); | |
656 | G4Material* material = logicalVol->GetMaterial(); | |
657 | G4int matIndex = 2; | |
658 | if (G4FlukaCompoundMap[material]) | |
659 | matIndex = G4FlukaCompoundMap[material]->GetIndex(); | |
660 | if (G4FlukaMaterialMap[material]) | |
661 | matIndex = G4FlukaMaterialMap[material]->GetIndex(); | |
662 | ||
663 | //Find if there is a magnetic field in the region | |
664 | //check if Magnetic Field is present in the region | |
665 | G4double flagField = 0.0; | |
666 | G4FieldManager * pMagFieldMan = logicalVol->GetFieldManager(); | |
667 | if(pMagFieldMan && pMagFieldMan->GetDetectorField()) | |
668 | flagField = 1.0; | |
669 | ||
670 | //Print card | |
671 | os << setw10 << "ASSIGNMAT "; | |
672 | os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield); | |
673 | os << setw10 << setfixed << G4double(matIndex); | |
674 | os << setw10 << setfixed << G4double(iFlukaRegion); | |
675 | os << setw10 << "0.0"; | |
676 | os << setw10 << setfixed << flagField; | |
677 | os << G4endl; | |
678 | } | |
679 | ||
680 | ||
681 | ||
682 | #ifdef G4GEOMETRY_DEBUG | |
683 | G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl; | |
684 | #endif | |
685 | } | |
686 | ||
687 | ||
688 | void FGeometryInit::PrintMagneticField(G4std::ostream& os) { | |
689 | #ifdef G4GEOMETRY_DEBUG | |
690 | G4cout << "==> Flugg FGeometryInit::PrintMagneticField()" << G4endl; | |
691 | #endif | |
692 | ||
693 | G4cout << "\t* Printing Magnetic Field..." << G4endl; | |
694 | ||
26911512 | 695 | if(fTransportationManager->GetFieldManager()->DoesFieldExist()) { |
26911512 | 696 | |
697 | //get magnetic field pointer | |
698 | const G4Field * pMagField = | |
699 | fTransportationManager->GetFieldManager()->GetDetectorField(); | |
700 | ||
26911512 | 701 | |
702 | if(pMagField) { | |
26d97e06 | 703 | //Check if it can be made a uniform magnetic field |
26911512 | 704 | const G4UniformMagField *pUnifMagField = |
705 | dynamic_cast<const G4UniformMagField*>(pMagField); | |
706 | if(pUnifMagField) { | |
26d97e06 | 707 | G4double B[3]; |
708 | G4double point[4]; //it is not really used | |
709 | pUnifMagField->GetFieldValue(point,B); | |
710 | ||
711 | //write MGNFIELD card | |
712 | PrintHeader(os,"GEANT4 MAGNETIC FIELD"); | |
713 | os << setw10 << "MGNFIELD "; | |
714 | os << setw10 << ""; | |
715 | os << setw10 << ""; | |
716 | os << setw10 << ""; | |
717 | os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield); | |
718 | os << setw10 << setfixed | |
719 | << G4std::setprecision(4) << B[0] | |
720 | << setw10 << B[1] | |
721 | << setw10 << B[2] | |
722 | << G4endl; | |
723 | } | |
724 | else { | |
725 | G4cout << "WARNING: No Uniform Magnetic Field found." << G4endl; | |
726 | G4cout << " Manual intervention might be needed." << G4endl; | |
26911512 | 727 | } |
26911512 | 728 | } |
26d97e06 | 729 | else |
730 | G4cout << "\t No detector field found... " << G4endl; | |
26911512 | 731 | } // end if magnetic field |
26d97e06 | 732 | else |
733 | G4cout << "\t No field found... " << G4endl; | |
734 | ||
e73e0522 | 735 | #ifdef G4GEOMETRY_DEBUG |
26d97e06 | 736 | G4cout << "<== Flugg FGeometryInit::PrintMagneticField()" << G4endl; |
e73e0522 | 737 | #endif |
26911512 | 738 | } |