]>
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 | ||
0edf14e7 | 80 | void FGeometryInit::Init() { |
81 | // Build and initialize G4 geometry | |
82 | setDetector(); | |
83 | setMotherVolume(); | |
84 | closeGeometry(); | |
85 | InitHistories(); | |
86 | InitJrLtGeantArray(); | |
87 | InitHistArray(); | |
88 | createFlukaMatFile(); | |
89 | } | |
90 | ||
91 | ||
26911512 | 92 | void FGeometryInit::closeGeometry() { |
e73e0522 | 93 | #ifdef G4GEOMETRY_DEBUG |
94 | G4cout << "==> Flugg FGeometryInit::closeGeometry()" << G4endl; | |
95 | #endif | |
26911512 | 96 | |
97 | ptrGeoMan = G4GeometryManager::GetInstance(); | |
98 | if (ptrGeoMan) { | |
99 | ptrGeoMan->OpenGeometry(); | |
100 | ||
101 | //true argoment allows voxel construction; if false voxels are built | |
102 | //only for replicated volumes | |
103 | ptrGeoMan->CloseGeometry(true); | |
104 | } | |
105 | else { | |
106 | G4cerr << "ERROR in FLUGG: Could not get G4GeometryManager instance" | |
107 | << G4endl; | |
108 | G4cerr << " in FGeometryInit::closeGeometry. Exiting!!!" | |
109 | << G4endl; | |
110 | } | |
111 | ||
e73e0522 | 112 | #ifdef G4GEOMETRY_DEBUG |
113 | G4cout << "<== Flugg FGeometryInit::closeGeometry()" << G4endl; | |
114 | #endif | |
26911512 | 115 | } |
116 | ||
117 | //************************************************************************* | |
118 | ||
119 | void FGeometryInit::InitHistArray() { | |
e73e0522 | 120 | #ifdef G4GEOMETRY_DEBUG |
121 | G4cout << "==> Flugg FGeometryInit::InitHistArray()" << G4endl; | |
122 | #endif | |
26911512 | 123 | ptrArray = new G4int[1000000]; |
124 | for(G4int i=0;i<1000000;i++) | |
125 | ptrArray[i]=0; | |
e73e0522 | 126 | #ifdef G4GEOMETRY_DEBUG |
127 | G4cout << "<== Flugg FGeometryInit::InitHistArray()" << G4endl; | |
128 | #endif | |
26911512 | 129 | } |
130 | ||
131 | ||
132 | ||
133 | //************************************************************************* | |
134 | //jrLtGeant stores all crossed lattice volume histories. | |
135 | ||
136 | void FGeometryInit::InitJrLtGeantArray() { | |
26911512 | 137 | #ifdef G4GEOMETRY_DEBUG |
e73e0522 | 138 | G4cout << "==> Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl; |
26911512 | 139 | G4cout << "Initializing JrLtGeant array" << G4endl; |
140 | #endif | |
141 | ptrJrLtGeant = new G4int[10000]; | |
142 | for(G4int x=0;x<10000;x++) | |
143 | ptrJrLtGeant[x]=-1; | |
144 | flagLttcGeant = -1; | |
e73e0522 | 145 | #ifdef G4GEOMETRY_DEBUG |
146 | G4cout << "<== Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl; | |
147 | #endif | |
26911512 | 148 | } |
149 | ||
150 | ||
151 | void FGeometryInit::SetLttcFlagGeant(G4int newFlagLttc) { | |
e73e0522 | 152 | #ifdef G4GEOMETRY_DEBUG |
153 | G4cout << "==> Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl; | |
154 | #endif | |
26911512 | 155 | // Added by A.Solodkov |
156 | if (newFlagLttc >= 10000) { | |
157 | G4cout << "Problems in FGeometryInit::SetLttcFlagGeant" << G4endl; | |
158 | G4cout << "Index newFlagLttc=" << newFlagLttc << " is outside array bounds" | |
159 | << G4endl; | |
160 | G4cout << "Better to stop immediately !" << G4endl; | |
161 | exit(1); | |
162 | } | |
163 | flagLttcGeant = newFlagLttc; | |
e73e0522 | 164 | #ifdef G4GEOMETRY_DEBUG |
165 | G4cout << "<== Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl; | |
166 | #endif | |
26911512 | 167 | } |
168 | ||
169 | void FGeometryInit::PrintJrLtGeant() { | |
170 | #ifdef G4GEOMETRY_DEBUG | |
171 | //G4cout << "jrLtGeant:" << G4endl; | |
172 | //for(G4int y=0;y<=flagLttcGeant;y++) | |
173 | // | |
174 | // G4cout << "jrLtGeant[" << y << "]=" << ptrJrLtGeant[y] << G4endl; | |
175 | #endif | |
176 | } | |
177 | ||
178 | //************************************************************************** | |
179 | ||
180 | void FGeometryInit::PrintHistories() { | |
181 | /* | |
182 | #ifdef G4GEOMETRY_DEBUG | |
183 | G4cout << "Touch hist:" << G4endl; | |
184 | G4cout << *(ptrTouchHist->GetHistory()) << G4endl; | |
185 | G4cout << "Tmp hist:" << G4endl; | |
186 | G4cout << *(ptrTempNavHist->GetHistory()) << G4endl; | |
187 | G4cout << "Old hist:" << G4endl; | |
188 | G4cout << *(ptrOldNavHist->GetHistory()) << G4endl; | |
189 | #endif | |
190 | */ | |
191 | } | |
192 | ||
193 | ||
194 | ||
195 | void FGeometryInit::InitHistories() { | |
e73e0522 | 196 | #ifdef G4GEOMETRY_DEBUG |
197 | G4cout << "==> Flugg FGeometryInit::InitHistories()" << G4endl; | |
198 | #endif | |
26911512 | 199 | //init utility histories with navigator history |
200 | ptrTouchHist = | |
201 | fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory(); | |
202 | ptrTempNavHist = | |
203 | fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory(); | |
204 | ptrOldNavHist = new G4TouchableHistory(); | |
e73e0522 | 205 | #ifdef G4GEOMETRY_DEBUG |
206 | G4cout << "<== Flugg FGeometryInit::InitHistories()" << G4endl; | |
207 | #endif | |
26911512 | 208 | } |
209 | ||
210 | void FGeometryInit::DeleteHistories() { | |
e73e0522 | 211 | #ifdef G4GEOMETRY_DEBUG |
212 | G4cout << "==> Flugg FGeometryInit::DeleteHistories()" << G4endl; | |
213 | #endif | |
214 | ||
26911512 | 215 | delete ptrTouchHist; |
216 | delete ptrOldNavHist; | |
217 | delete ptrTempNavHist; | |
218 | ||
219 | #ifdef G4GEOMETRY_DEBUG | |
220 | G4cout << "Deleting step-history objects at end of run!" << G4endl; | |
e73e0522 | 221 | G4cout << "<== Flugg FGeometryInit::DeleteHistories()" << G4endl; |
26911512 | 222 | #endif |
26911512 | 223 | } |
224 | ||
225 | ||
226 | void FGeometryInit::UpdateHistories(const G4NavigationHistory * history, | |
227 | G4int flagHist) { | |
e73e0522 | 228 | #ifdef G4GEOMETRY_DEBUG |
229 | G4cout << "==> Flugg FGeometryInit::UpdateHistories()" << G4endl; | |
230 | #endif | |
26911512 | 231 | PrintHistories(); |
232 | ||
233 | #ifdef G4GEOMETRY_DEBUG | |
234 | G4cout << "...updating histories!" << G4endl; | |
235 | #endif | |
236 | ||
237 | G4VPhysicalVolume * pPhysVol = history->GetTopVolume(); | |
238 | ||
239 | switch (flagHist) { | |
240 | case 0: { | |
241 | //this is the case when a new history is given to the | |
242 | //navigator and old history has to be resetted | |
243 | //touchable history has not been updated jet, so: | |
244 | ||
245 | ptrTouchHist->UpdateYourself(pPhysVol,history); | |
246 | ptrTempNavHist->UpdateYourself(pPhysVol,history); | |
247 | G4NavigationHistory * ptrOldNavHistNotConst = | |
248 | const_cast<G4NavigationHistory * >(ptrOldNavHist->GetHistory()); | |
249 | ptrOldNavHistNotConst->Reset(); | |
250 | ptrOldNavHistNotConst->Clear(); | |
251 | PrintHistories(); | |
252 | break; | |
253 | } //case 0 | |
254 | ||
255 | case 1: { | |
256 | //this is the case when a new history is given to the | |
257 | //navigator but old history has to be kept (e.g. LOOKZ | |
258 | //is call during an event); | |
259 | //touchable history has not been updated jet, so: | |
260 | ||
261 | ptrTouchHist->UpdateYourself(pPhysVol,history); | |
262 | ptrTempNavHist->UpdateYourself(pPhysVol,history); | |
263 | PrintHistories(); | |
264 | break; | |
265 | } //case 1 | |
266 | ||
267 | case 2: { | |
268 | //this is the case when the touchable history has been | |
269 | //updated by a LocateGlobalPointAndUpdateTouchable call | |
270 | ||
271 | G4VPhysicalVolume * pPhysVolTemp = ptrTempNavHist->GetVolume(); | |
272 | ptrOldNavHist->UpdateYourself(pPhysVolTemp, | |
273 | ptrTempNavHist->GetHistory()); | |
274 | ||
275 | ptrTempNavHist->UpdateYourself(pPhysVol,history); | |
276 | PrintHistories(); | |
277 | break; | |
278 | } //case 2 | |
279 | ||
280 | default: { | |
281 | G4cout <<" ERROR in updating step-histories!" << G4endl; | |
282 | break; | |
283 | } //default | |
284 | } //switch | |
285 | ||
e73e0522 | 286 | #ifdef G4GEOMETRY_DEBUG |
287 | G4cout << "<== Flugg FGeometryInit::UpdateHistories()" << G4endl; | |
288 | #endif | |
26911512 | 289 | } |
290 | ||
291 | //***************************************************************************** | |
b9d7c32a | 292 | int FGeometryInit::GetLastMaterialIndex() const |
293 | { | |
294 | // Get last material index as known by FLUKA | |
295 | const FlukaMaterialsTable *matTable = FlukaMaterial::GetMaterialTable(); | |
296 | int matsize = matTable->size(); | |
297 | return matsize+2; | |
298 | } | |
26911512 | 299 | |
26911512 | 300 | void FGeometryInit::createFlukaMatFile() { |
26911512 | 301 | // last modification Sara Vanini 1/III/99 |
302 | // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case, | |
303 | // according to the fluka standard. In addition,. they must be equal to the | |
304 | // names of the fluka materials - see fluka manual - in order that the | |
305 | // program load the right cross sections, and equal to the names included in | |
306 | // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his | |
307 | // own .pemf, in order to get the right cross sections loaded in memory. | |
308 | ||
309 | #ifdef G4GEOMETRY_DEBUG | |
e73e0522 | 310 | G4cout << "==> Flugg FGeometryInit::createFlukaMatFile()" << G4endl; |
26911512 | 311 | G4cout << "================== FILEWR =================" << G4endl; |
312 | #endif | |
26d97e06 | 313 | |
314 | ||
315 | //Regions map | |
316 | BuildRegionsMap(); | |
0edf14e7 | 317 | std::ofstream vos("Volumes_index.inp"); |
26d97e06 | 318 | PrintRegionsMap(vos); |
319 | vos.close(); | |
320 | ||
321 | //Materials and compounds | |
322 | BuildMaterialTables(); | |
0edf14e7 | 323 | std::ofstream fos("flukaMat.inp"); |
26d97e06 | 324 | PrintMaterialTables(fos); |
325 | PrintAssignmat(fos); | |
326 | PrintMagneticField(fos); | |
327 | fos.close(); | |
328 | ||
329 | #ifdef G4GEOMETRY_DEBUG | |
330 | G4cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << G4endl; | |
331 | #endif | |
332 | } | |
333 | ||
334 | //////////////////////////////////////////////////////////////////////// | |
335 | // | |
336 | void FGeometryInit::BuildRegionsMap() { | |
337 | #ifdef G4GEOMETRY_DEBUG | |
338 | G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl; | |
339 | #endif | |
340 | ||
341 | //Find number of Volumes in physical volume store | |
26911512 | 342 | G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance(); |
26d97e06 | 343 | unsigned int numVol = pVolStore->size(); |
344 | ||
345 | G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore | |
346 | << ") has " << numVol << " volumes. Iterating..." | |
347 | << G4endl; | |
348 | ||
349 | for(unsigned int l=0; l < numVol; l++) { | |
350 | //Get each of the physical volumes | |
351 | G4VPhysicalVolume * physicalvolume = (*pVolStore)[l]; | |
352 | G4int iFlukaRegion = l+1; | |
353 | fRegionVolumeMap[physicalvolume] = iFlukaRegion; | |
26911512 | 354 | } |
26d97e06 | 355 | |
356 | ||
357 | ||
26911512 | 358 | #ifdef G4GEOMETRY_DEBUG |
26d97e06 | 359 | G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl; |
360 | #endif | |
361 | } | |
26911512 | 362 | |
0edf14e7 | 363 | void FGeometryInit::PrintRegionsMap(std::ostream& os) { |
26911512 | 364 | #ifdef G4GEOMETRY_DEBUG |
26d97e06 | 365 | G4cout << "==> Flugg FGeometryInit::PrintRegionsMap()" << G4endl; |
366 | #endif | |
26911512 | 367 | |
26d97e06 | 368 | //Print some header |
369 | PrintHeader(os, "GEANT4 VOLUMES"); | |
26911512 | 370 | |
26d97e06 | 371 | //Iterate over all volumes in the map |
372 | for (RegionIterator i = fRegionVolumeMap.begin(); | |
373 | i != fRegionVolumeMap.end(); | |
374 | i++) { | |
26911512 | 375 | |
26d97e06 | 376 | //Get info in the map |
377 | G4VPhysicalVolume* ptrVol = (*i).first; | |
378 | int index = (*i).second; | |
26911512 | 379 | |
26d97e06 | 380 | //Print index and region name in some fixed format |
0edf14e7 | 381 | os.setf(std::ios::left, std::ios::adjustfield); |
26d97e06 | 382 | os << setw10 << index; |
0edf14e7 | 383 | os << std::setw(20) << ptrVol->GetName() << std::setw(20) << ""; |
26d97e06 | 384 | |
385 | //If volume is a replica... print some more stuff | |
26911512 | 386 | if(ptrVol->IsReplicated()) { |
387 | EAxis axis; | |
26d97e06 | 388 | G4int nRep = -1; |
389 | G4double width = -1; | |
390 | G4double offset = -1; | |
391 | G4bool consum = false; | |
392 | ptrVol->GetReplicationData(axis, nRep, width, offset, consum); | |
0edf14e7 | 393 | os.setf(std::ios::left, std::ios::adjustfield); |
394 | os << setw10 << "Repetion Nb: " << std::setw(3) << nRep; | |
26911512 | 395 | } |
26d97e06 | 396 | os << G4endl; |
26911512 | 397 | |
26d97e06 | 398 | } |
399 | ||
400 | #ifdef G4GEOMETRY_DEBUG | |
401 | G4cout << "<== Flugg FGeometryInit::PrintRegionsMap()" << G4endl; | |
402 | #endif | |
403 | } | |
404 | ||
bf547b2f | 405 | //////////////////////////////////////////////////////////////////////// |
1617d9fa | 406 | // |
407 | void FGeometryInit::BuildMediaMap() | |
408 | { | |
409 | fRegionMediumMap = new int[fNRegions+1]; | |
410 | for (RegionIterator i = fRegionVolumeMap.begin(); | |
411 | i != fRegionVolumeMap.end(); | |
412 | i++) { | |
413 | //Get info in the map | |
414 | G4VPhysicalVolume* ptrVol = (*i).first; | |
415 | int region = (*i).second; | |
416 | G4int imed = fMediumVolumeMap[ptrVol]; | |
417 | fRegionMediumMap[region] = imed; | |
418 | printf("BuildMediaMap %s %d %d\n",(ptrVol->GetName()).data(), region, imed); | |
419 | ||
420 | } | |
421 | } | |
422 | ||
423 | G4int FGeometryInit::GetMedium(int region) const | |
424 | { | |
425 | return fRegionMediumMap[region]; | |
426 | } | |
427 | ||
428 | ||
6a53de92 | 429 | void FGeometryInit::SetMediumFromName(const char* volName, int medium, int volid) |
1617d9fa | 430 | { |
431 | char name4[5]; | |
432 | char tmp[5]; | |
433 | strncpy(tmp, volName, 4); | |
434 | tmp[4]='\0'; | |
435 | fNRegions = 0; | |
436 | ||
bf547b2f | 437 | for (RegionIterator i = fRegionVolumeMap.begin(); |
438 | i != fRegionVolumeMap.end(); | |
439 | i++) { | |
1617d9fa | 440 | fNRegions++; |
bf547b2f | 441 | //Get info in the map |
442 | G4VPhysicalVolume* ptrVol = (*i).first; | |
1617d9fa | 443 | strncpy(name4, (ptrVol->GetName()).data(), 4); |
444 | name4[4]='\0'; | |
445 | for (int j = 0; j < 4; j++) { | |
446 | if (name4[j] == '\0') { | |
447 | for (int k = j; k < 4; k++) { | |
448 | name4[k] = ' '; | |
449 | } | |
450 | break; | |
451 | } | |
452 | } | |
6a53de92 | 453 | if (! strncmp(name4, tmp, 4)) { |
454 | fMediumVolumeMap[ptrVol] = medium; | |
455 | fVolIdVolumeMap[ptrVol] = volid; | |
456 | } | |
bf547b2f | 457 | } |
bf547b2f | 458 | } |
459 | ||
460 | ||
461 | ||
26d97e06 | 462 | //////////////////////////////////////////////////////////////////////// |
463 | // | |
464 | void FGeometryInit::BuildMaterialTables() { | |
465 | #ifdef G4GEOMETRY_DEBUG | |
466 | G4cout << "==> Flugg FGeometryInit::BuildMaterialTables()" << G4endl; | |
467 | #endif | |
468 | ||
469 | //some terminal printout also | |
470 | G4cout << "\t* Storing information..." << G4endl; | |
471 | ||
472 | //The logic is the folloing: | |
473 | //Get the Material Table and: | |
474 | // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum | |
475 | // 2) For each single element material build a material equivalent | |
476 | // 3) For the rest: | |
477 | // 3.a) Build materials for each not already known element | |
478 | // 3.b) Build the compound out of them | |
479 | ||
480 | //Get the Material Table and iterate | |
481 | const G4MaterialTable* matTable = G4Material::GetMaterialTable(); | |
482 | for (MatTableIterator i = matTable->begin(); i != matTable->end(); i++) { | |
483 | ||
484 | //Get some basic material information | |
485 | G4Material* material = (*i); | |
486 | G4String matName = material->GetName(); | |
487 | const G4double matDensity = material->GetDensity(); | |
488 | const G4int nMatElements = material->GetNumberOfElements(); | |
489 | ||
490 | G4cout << "\t\t+ " << matName | |
491 | << ": dens. = " << matDensity/(g/cm3) << "g/cm3" | |
492 | << ", nElem = " << nMatElements << G4endl; | |
493 | ||
494 | // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum | |
495 | // FlukaMaterial* is 0 in that case | |
496 | if (matDensity <= 1.00e-10*g/cm3) { | |
497 | G4FlukaMaterialMap[material] = 0; | |
498 | G4cout << "\t\t Stored as vacuum" << G4endl; | |
26911512 | 499 | } |
26d97e06 | 500 | // 2) For each single element material build a material equivalent |
501 | else if (nMatElements == 1) { | |
26911512 | 502 | |
26d97e06 | 503 | FlukaMaterial *flukamat = |
504 | BuildFlukaMaterialFromElement(material->GetElement(0), | |
505 | matDensity); | |
26911512 | 506 | |
26d97e06 | 507 | G4FlukaMaterialMap[material] = flukamat; |
508 | G4cout << "\t\t Stored as " << flukamat->GetRealName() << G4endl; | |
26911512 | 509 | |
26d97e06 | 510 | } //else if (material->GetNumberOfElements() == 1) |
511 | ||
512 | // 3) For the rest: | |
513 | // 3.a) Build materials for each not already known element | |
514 | // 3.b) Build the compound out of them | |
515 | else { | |
516 | FlukaCompound* flukacomp = | |
517 | BuildFlukaCompoundFromMaterial(material); | |
518 | G4FlukaCompoundMap[material] = flukacomp; | |
519 | G4cout << "\t\t Stored as " << flukacomp->GetRealName() << G4endl; | |
520 | } //else for case 3) | |
521 | } //for (materials) | |
522 | ||
523 | #ifdef G4GEOMETRY_DEBUG | |
524 | G4cout << "<== Flugg FGeometryInit::BuildMaterialTables()" << G4endl; | |
525 | #endif | |
526 | } | |
527 | ||
528 | FlukaMaterial* | |
529 | FGeometryInit::BuildFlukaMaterialFromElement(const G4Element* element, | |
530 | G4double matDensity) { | |
531 | #ifdef G4GEOMETRY_DEBUG | |
532 | G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromElement()" | |
533 | << G4endl; | |
534 | #endif | |
535 | ||
536 | //Get element and its properties | |
537 | G4String elemName(ToFlukaString(element->GetName())); | |
538 | ||
539 | FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(elemName); | |
540 | if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) { | |
541 | //Check for isotopes | |
542 | G4int nIsotopes = element->GetNumberOfIsotopes(); | |
543 | if (nIsotopes == 0) { | |
544 | G4double elemA = element->GetA()/g; | |
545 | G4double elemZ = element->GetZ(); | |
26911512 | 546 | |
26d97e06 | 547 | if (elemA != G4int(elemA) && elemZ != G4int(elemZ)) { |
548 | G4cout << "WARNING: Element \'" << elemName | |
549 | << "\' has non integer Z (" << elemZ << ") or A (" | |
550 | << elemA << ")" | |
551 | << G4endl; | |
26911512 | 552 | } |
26d97e06 | 553 | |
554 | flukamat = new FlukaMaterial(elemName, | |
555 | G4int(elemZ), | |
556 | elemA, | |
557 | matDensity/(g/cm3)); | |
558 | } | |
559 | else if (nIsotopes == 1) { | |
560 | const G4Isotope* isotope = element->GetIsotope(0); | |
561 | flukamat = BuildFlukaMaterialFromIsotope(isotope, matDensity); | |
562 | } | |
563 | else { | |
564 | FlukaCompound *flucomp = BuildFlukaCompoundFromElement(element, | |
565 | matDensity); | |
566 | flukamat = flucomp->GetFlukaMaterial(); | |
26911512 | 567 | } |
26d97e06 | 568 | } |
569 | #ifdef G4GEOMETRY_DEBUG | |
570 | else { | |
571 | G4cout << "INFO: Element \'" << elemName | |
572 | << "\' already exists in the DB. It will not be recreated." | |
573 | << G4endl; | |
574 | } | |
575 | #endif | |
576 | ||
577 | return flukamat; | |
26911512 | 578 | |
26d97e06 | 579 | #ifdef G4GEOMETRY_DEBUG |
580 | G4cout << "<== Flugg FGeometryInit::BuildFlukaMaterialFromElement()" | |
581 | << G4endl; | |
582 | #endif | |
583 | } | |
584 | ||
585 | FlukaMaterial* | |
586 | FGeometryInit::BuildFlukaMaterialFromIsotope(const G4Isotope* isotope, | |
587 | G4double matDensity) { | |
588 | #ifdef G4GEOMETRY_DEBUG | |
589 | G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()" | |
590 | << G4endl; | |
591 | #endif | |
592 | G4String isoName(ToFlukaString(isotope->GetName())); | |
593 | FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(isoName); | |
594 | if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) { | |
595 | G4int isoZ = isotope->GetZ(); | |
596 | G4double isoA = (isotope->GetA())/(g); | |
597 | G4int isoN = isotope->GetN(); | |
598 | flukamat = new FlukaMaterial(isoName, | |
599 | isoZ, | |
600 | isoA, | |
601 | matDensity/(g/cm3), | |
602 | isoN); | |
603 | } | |
604 | ||
605 | return flukamat; | |
606 | ||
607 | #ifdef G4GEOMETRY_DEBUG | |
608 | G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()" | |
609 | << G4endl; | |
610 | #endif | |
611 | } | |
612 | ||
613 | FlukaCompound* | |
614 | FGeometryInit::BuildFlukaCompoundFromMaterial(const G4Material* material) { | |
615 | #ifdef G4GEOMETRY_DEBUG | |
616 | G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()" | |
617 | << G4endl; | |
618 | #endif | |
619 | //Material properties | |
620 | const G4double* elemFractions = material->GetFractionVector(); | |
621 | const G4int nMatElements = material->GetNumberOfElements(); | |
622 | const G4double matDensity = material->GetDensity(); | |
623 | G4String matName(ToFlukaString(material->GetName())); | |
624 | FlukaCompound* flukacomp = new FlukaCompound(matName, matDensity/(g/cm3), | |
625 | nMatElements); | |
626 | for (G4int i = 0; i < nMatElements; i++) { | |
627 | FlukaMaterial *flukamat = | |
628 | BuildFlukaMaterialFromElement(material->GetElement(i), 0.0); | |
629 | ||
630 | flukacomp->AddElement(flukamat->GetIndex(), -elemFractions[i]); | |
631 | ||
632 | } //for (elements) | |
633 | ||
634 | return flukacomp; | |
635 | ||
636 | #ifdef G4GEOMETRY_DEBUG | |
637 | G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()" | |
638 | << G4endl; | |
639 | #endif | |
640 | } | |
641 | ||
642 | FlukaCompound* | |
643 | FGeometryInit::BuildFlukaCompoundFromElement(const G4Element* element, | |
644 | G4double matDensity) { | |
645 | #ifdef G4GEOMETRY_DEBUG | |
646 | G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromElement()" | |
647 | << G4endl; | |
648 | #endif | |
649 | G4int nIsotopes = element->GetNumberOfIsotopes(); | |
650 | //fraction of nb of atomes per volume (= volume fraction?) | |
651 | const G4double* isoAbundance = element->GetRelativeAbundanceVector(); | |
652 | G4String elemName(ToFlukaString(element->GetName())); | |
653 | ||
654 | //Material properties | |
655 | FlukaCompound* flukacomp = new FlukaCompound(elemName, matDensity/(g/cm3), | |
656 | nIsotopes); | |
657 | for (G4int i = 0; i < nIsotopes; i++) { | |
658 | FlukaMaterial *flukamat = | |
659 | BuildFlukaMaterialFromIsotope(element->GetIsotope(i), 0.0); | |
660 | ||
661 | flukacomp->AddElement(flukamat->GetIndex(), isoAbundance[i]); | |
662 | ||
663 | } //for (elements) | |
664 | ||
665 | return flukacomp; | |
666 | ||
667 | #ifdef G4GEOMETRY_DEBUG | |
668 | G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromElement()" | |
669 | << G4endl; | |
670 | #endif | |
671 | } | |
672 | ||
0edf14e7 | 673 | void FGeometryInit::PrintMaterialTables(std::ostream& os) { |
26d97e06 | 674 | #ifdef G4GEOMETRY_DEBUG |
675 | G4cout << "==> Flugg FGeometryInit::PrintMaterialTables()" << G4endl; | |
676 | #endif | |
677 | //Print Header | |
678 | PrintHeader(os, "GEANT4 MATERIALS AND COMPOUNDS"); | |
26911512 | 679 | |
26d97e06 | 680 | //And some more stuff |
681 | size_t nIsotopes = G4Isotope::GetNumberOfIsotopes(); | |
682 | size_t nElements = G4Element::GetNumberOfElements(); | |
683 | size_t nMaterials = G4Material::GetNumberOfMaterials(); | |
684 | ||
0e22711e | 685 | os << "* In Geant4 there are " << nMaterials << " materials" << G4endl; |
686 | os << "* In Geant4 there are " << nElements << " elements" << G4endl; | |
687 | os << "* In Geant4 there are " << nIsotopes << " isotopes" << G4endl; | |
26d97e06 | 688 | |
689 | //Materials | |
690 | G4cout << "\t* Printing FLUKA materials..." << G4endl; | |
691 | FlukaMaterial::PrintMaterialsByIndex(os); | |
692 | //FlukaMaterial::PrintMaterialsByName(os); | |
693 | ||
694 | //Compounds | |
695 | G4cout << "\t* Printing FLUKA compounds..." << G4endl; | |
696 | FlukaCompound::PrintCompounds(os); | |
697 | ||
698 | #ifdef G4GEOMETRY_DEBUG | |
699 | G4cout << "<== Flugg FGeometryInit::PrintMaterialTables()" << G4endl; | |
700 | #endif | |
701 | } | |
702 | ||
703 | //////////////////////////////////////////////////////////////////////// | |
704 | // | |
0edf14e7 | 705 | void FGeometryInit::PrintAssignmat(std::ostream& os) { |
26d97e06 | 706 | #ifdef G4GEOMETRY_DEBUG |
707 | G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl; | |
708 | #endif | |
709 | ||
710 | //Find number of Volumes in physical volume store | |
711 | G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance(); | |
712 | unsigned int numVol = pVolStore->size(); | |
713 | ||
714 | G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore | |
715 | << ") has " << numVol << " volumes. " << G4endl; | |
716 | ||
717 | G4cout << "\t* Printing ASSIGNMAT..." << G4endl; | |
718 | ||
719 | ||
720 | PrintHeader(os,"GEANT4 MATERIAL ASSIGNMENTS"); | |
721 | for(unsigned int l=0; l < numVol; l++) { | |
722 | ||
723 | //Get each of the physical volumes | |
724 | G4VPhysicalVolume * physicalvol = (*pVolStore)[l]; | |
725 | ||
726 | //Get index for that volume | |
727 | G4int iFlukaRegion = fRegionVolumeMap[physicalvol]; | |
728 | ||
729 | //Find G4 material and navigate to its fluka compound/material | |
730 | G4LogicalVolume * logicalVol = physicalvol->GetLogicalVolume(); | |
731 | G4Material* material = logicalVol->GetMaterial(); | |
732 | G4int matIndex = 2; | |
733 | if (G4FlukaCompoundMap[material]) | |
734 | matIndex = G4FlukaCompoundMap[material]->GetIndex(); | |
735 | if (G4FlukaMaterialMap[material]) | |
736 | matIndex = G4FlukaMaterialMap[material]->GetIndex(); | |
737 | ||
738 | //Find if there is a magnetic field in the region | |
739 | //check if Magnetic Field is present in the region | |
740 | G4double flagField = 0.0; | |
741 | G4FieldManager * pMagFieldMan = logicalVol->GetFieldManager(); | |
742 | if(pMagFieldMan && pMagFieldMan->GetDetectorField()) | |
743 | flagField = 1.0; | |
744 | ||
745 | //Print card | |
746 | os << setw10 << "ASSIGNMAT "; | |
0edf14e7 | 747 | os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); |
26d97e06 | 748 | os << setw10 << setfixed << G4double(matIndex); |
749 | os << setw10 << setfixed << G4double(iFlukaRegion); | |
750 | os << setw10 << "0.0"; | |
751 | os << setw10 << setfixed << flagField; | |
752 | os << G4endl; | |
753 | } | |
754 | ||
755 | ||
756 | ||
757 | #ifdef G4GEOMETRY_DEBUG | |
758 | G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl; | |
759 | #endif | |
760 | } | |
761 | ||
762 | ||
0edf14e7 | 763 | void FGeometryInit::PrintMagneticField(std::ostream& os) { |
26d97e06 | 764 | #ifdef G4GEOMETRY_DEBUG |
765 | G4cout << "==> Flugg FGeometryInit::PrintMagneticField()" << G4endl; | |
766 | #endif | |
767 | ||
768 | G4cout << "\t* Printing Magnetic Field..." << G4endl; | |
769 | ||
26911512 | 770 | if(fTransportationManager->GetFieldManager()->DoesFieldExist()) { |
26911512 | 771 | |
772 | //get magnetic field pointer | |
773 | const G4Field * pMagField = | |
774 | fTransportationManager->GetFieldManager()->GetDetectorField(); | |
775 | ||
26911512 | 776 | |
777 | if(pMagField) { | |
26d97e06 | 778 | //Check if it can be made a uniform magnetic field |
26911512 | 779 | const G4UniformMagField *pUnifMagField = |
780 | dynamic_cast<const G4UniformMagField*>(pMagField); | |
781 | if(pUnifMagField) { | |
26d97e06 | 782 | G4double B[3]; |
783 | G4double point[4]; //it is not really used | |
784 | pUnifMagField->GetFieldValue(point,B); | |
785 | ||
786 | //write MGNFIELD card | |
787 | PrintHeader(os,"GEANT4 MAGNETIC FIELD"); | |
788 | os << setw10 << "MGNFIELD "; | |
789 | os << setw10 << ""; | |
790 | os << setw10 << ""; | |
791 | os << setw10 << ""; | |
0edf14e7 | 792 | os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield); |
26d97e06 | 793 | os << setw10 << setfixed |
0edf14e7 | 794 | << std::setprecision(4) << B[0] |
26d97e06 | 795 | << setw10 << B[1] |
796 | << setw10 << B[2] | |
797 | << G4endl; | |
798 | } | |
799 | else { | |
800 | G4cout << "WARNING: No Uniform Magnetic Field found." << G4endl; | |
801 | G4cout << " Manual intervention might be needed." << G4endl; | |
26911512 | 802 | } |
26911512 | 803 | } |
26d97e06 | 804 | else |
805 | G4cout << "\t No detector field found... " << G4endl; | |
26911512 | 806 | } // end if magnetic field |
26d97e06 | 807 | else |
808 | G4cout << "\t No field found... " << G4endl; | |
809 | ||
e73e0522 | 810 | #ifdef G4GEOMETRY_DEBUG |
26d97e06 | 811 | G4cout << "<== Flugg FGeometryInit::PrintMagneticField()" << G4endl; |
e73e0522 | 812 | #endif |
26911512 | 813 | } |
6a53de92 | 814 | |
815 | int FGeometryInit::CurrentVolID(int ir, int& copyNo) | |
816 | { | |
086d0acf | 817 | if (ir == 0) |
818 | { | |
819 | copyNo = -1; | |
820 | return -1; | |
821 | } | |
822 | ||
6a53de92 | 823 | G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance(); |
086d0acf | 824 | G4VPhysicalVolume * physicalvol = (*pVolStore)[ir- 1]; |
825 | ||
826 | if (physicalvol) { | |
827 | copyNo = physicalvol->GetCopyNo(); | |
828 | } else { | |
829 | copyNo = -1; | |
830 | return -1; | |
831 | } | |
832 | ||
833 | ||
6a53de92 | 834 | int id = fVolIdVolumeMap[physicalvol]; |
835 | return id; | |
836 | } | |
837 | ||
838 | int FGeometryInit::CurrentVolOffID(int ir, int off, int& copyNo) | |
839 | { | |
840 | if (off == 0) return CurrentVolID(ir, copyNo); | |
841 | ||
842 | G4PhysicalVolumeStore* pVolStore = G4PhysicalVolumeStore::GetInstance(); | |
843 | G4VPhysicalVolume* physicalvol = (*pVolStore)[ir- 1]; | |
844 | G4VPhysicalVolume* mother = physicalvol; | |
845 | ||
0edf14e7 | 846 | int index; |
847 | //============================================================================ | |
848 | if (mother) { | |
849 | // Check touchable depth | |
850 | // | |
851 | if (ptrTouchHist->GetHistoryDepth() < off) { | |
852 | mother = 0; | |
853 | } else { | |
854 | // Get the off-th mother | |
855 | index = ptrTouchHist->GetHistoryDepth() - off; | |
856 | // in the touchable history volumes are ordered | |
857 | // from top volume up to mother volume; | |
858 | // the touchable volume is not in the history | |
859 | mother = ptrTouchHist->GetHistory()->GetVolume(index); | |
860 | } | |
861 | } | |
862 | //============================================================================ | |
6a53de92 | 863 | |
864 | int id; | |
865 | ||
866 | if (!mother) { | |
867 | G4cout << "Flugg FGeometryInit::CurrentVolOffID mother not found" << G4endl; | |
868 | id = -1; | |
869 | copyNo = -1; | |
870 | } else { | |
871 | copyNo = mother ->GetCopyNo(); | |
872 | id = fVolIdVolumeMap[mother]; | |
873 | } | |
874 | return id; | |
875 | } | |
dc37cac6 | 876 | |
877 | void FGeometryInit::Gmtod(double* xm, double* xd, int iflag) | |
878 | { | |
879 | // Transforms a position from the world reference frame | |
880 | // to the current volume reference frame. | |
881 | // | |
882 | // Geant3 desription: | |
883 | // ================== | |
884 | // Computes coordinates XD (in DRS) | |
885 | // from known coordinates XM in MRS | |
886 | // The local reference system can be initialized by | |
887 | // - the tracking routines and GMTOD used in GUSTEP | |
888 | // - a call to GMEDIA(XM,NUMED) | |
889 | // - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER) | |
890 | // (inverse routine is GDTOM) | |
891 | // | |
892 | // If IFLAG=1 convert coordinates | |
893 | // IFLAG=2 convert direction cosinus | |
894 | // | |
895 | // --- | |
896 | FluggNavigator * ptrNavig = getNavigatorForTracking(); | |
897 | //setting variables (and dimension: Fluka uses cm.!) | |
898 | G4ThreeVector pGlob(xm[0],xm[1],xm[2]); | |
1b4955bb | 899 | G4ThreeVector pLoc; |
dc37cac6 | 900 | |
901 | if (iflag == 1) { | |
1b4955bb | 902 | pGlob *= 10.0; // in mm |
0edf14e7 | 903 | // change because of geant 4 6.0 |
904 | // pLoc = ptrNavig->ComputeLocalPoint(pGlob); | |
905 | pLoc = ptrNavig->GetGlobalToLocalTransform().TransformPoint(pGlob); | |
906 | ||
1b4955bb | 907 | pLoc /= 10.0; // in cm |
dc37cac6 | 908 | } else if (iflag == 2) { |
909 | pLoc = | |
910 | ptrNavig->ComputeLocalAxis(pGlob); | |
911 | } else { | |
912 | G4cout << "Flugg FGeometryInit::Gmtod called with undefined flag" << G4endl; | |
913 | } | |
dc37cac6 | 914 | xd[0] = pLoc[0]; xd[1] = pLoc[1]; xd[2] = pLoc[2]; |
915 | } | |
916 | ||
917 | void FGeometryInit::Gdtom(double* xd, double* xm, int iflag) | |
918 | { | |
919 | // Transforms a position from the current volume reference frame | |
920 | // to the world reference frame. | |
921 | // | |
922 | // Geant3 desription: | |
923 | // ================== | |
924 | // Computes coordinates XM (Master Reference System | |
925 | // knowing the coordinates XD (Detector Ref System) | |
926 | // The local reference system can be initialized by | |
927 | // - the tracking routines and GDTOM used in GUSTEP | |
928 | // - a call to GSCMED(NLEVEL,NAMES,NUMBER) | |
929 | // (inverse routine is GMTOD) | |
930 | // | |
931 | // If IFLAG=1 convert coordinates | |
932 | // IFLAG=2 convert direction cosinus | |
933 | // | |
934 | // --- | |
935 | ||
936 | FluggNavigator * ptrNavig = getNavigatorForTracking(); | |
937 | G4ThreeVector pLoc(xd[0],xd[1],xd[2]); | |
1b4955bb | 938 | |
dc37cac6 | 939 | G4ThreeVector pGlob; |
940 | if (iflag == 1) { | |
1b4955bb | 941 | pLoc *= 10.0; // in mm |
dc37cac6 | 942 | pGlob = ptrNavig->GetLocalToGlobalTransform(). |
943 | TransformPoint(pLoc); | |
1b4955bb | 944 | pGlob /= 10.0; // in cm |
dc37cac6 | 945 | } else if (iflag == 2) { |
946 | pGlob = ptrNavig->GetLocalToGlobalTransform(). | |
947 | TransformAxis(pLoc); | |
948 | } else { | |
949 | G4cout << "Flugg FGeometryInit::Gdtom called with undefined flag" << G4endl; | |
950 | } | |
951 | ||
952 | xm[0] = pGlob[0]; xm[1] = pGlob[1]; xm[2] = pGlob[2]; | |
953 | } |