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" |
26911512 |
10 | |
11 | FGeometryInit * FGeometryInit::flagInstance=0; |
12 | |
13 | FGeometryInit* FGeometryInit::GetInstance() { |
e73e0522 |
14 | #ifdef G4GEOMETRY_DEBUG |
15 | G4cout << "==> Flugg::FGeometryInit::GetInstance(), instance=" |
16 | << flagInstance << G4endl; |
17 | #endif |
26911512 |
18 | if (!flagInstance) |
19 | flagInstance = new FGeometryInit(); |
20 | |
e73e0522 |
21 | #ifdef G4GEOMETRY_DEBUG |
22 | G4cout << "<== Flugg::FGeometryInit::GetInstance(), instance=" |
23 | << flagInstance << G4endl; |
24 | #endif |
26911512 |
25 | return flagInstance; |
26 | } |
27 | |
28 | |
29 | FGeometryInit::FGeometryInit(): |
30 | fDetector(0), |
31 | fFieldManager(0), |
36c70081 |
32 | fTransportationManager(G4TransportationManager::GetTransportationManager()), |
26911512 |
33 | myTopNode(0), |
34 | ptrGeoMan(0), |
35 | ptrArray(0), |
36 | ptrTouchHist(0), |
37 | ptrOldNavHist(0), |
38 | ptrTempNavHist(0), |
36c70081 |
39 | ptrJrLtGeant(0) |
40 | { |
26911512 |
41 | |
e73e0522 |
42 | #ifdef G4GEOMETRY_DEBUG |
43 | G4cout << "==> Flugg FGeometryInit::FGeometryInit()" << G4endl; |
44 | G4cout << "\t+ Changing the G4Navigator for FluggNavigator..." << G4endl; |
45 | #endif |
26911512 |
46 | G4Navigator* actualnav = fTransportationManager->GetNavigatorForTracking(); |
47 | if (actualnav) { |
48 | FluggNavigator* newnav = new FluggNavigator(); |
49 | fTransportationManager->SetNavigatorForTracking(newnav); |
50 | } |
51 | else { |
36c70081 |
52 | G4cerr << "ERROR: Could not find the actual G4Navigator" << G4endl; |
26911512 |
53 | abort(); |
54 | } |
55 | |
56 | |
e73e0522 |
57 | #ifdef G4GEOMETRY_DEBUG |
58 | G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl; |
59 | #endif |
26911512 |
60 | |
61 | } |
62 | |
63 | |
64 | FGeometryInit::~FGeometryInit() { |
e73e0522 |
65 | G4cout << "==> Flugg FGeometryInit::~FGeometryInit()" << G4endl; |
26911512 |
66 | DeleteHistories(); |
67 | ptrGeoMan->OpenGeometry(); |
68 | if (fTransportationManager) |
69 | delete fTransportationManager; |
70 | if (ptrJrLtGeant) |
71 | delete[] ptrJrLtGeant; |
72 | DelHistArray(); |
73 | |
74 | //keep ATTENTION: never delete a pointer twice! |
e73e0522 |
75 | G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl; |
26911512 |
76 | } |
77 | |
78 | |
79 | void FGeometryInit::closeGeometry() { |
e73e0522 |
80 | #ifdef G4GEOMETRY_DEBUG |
81 | G4cout << "==> Flugg FGeometryInit::closeGeometry()" << G4endl; |
82 | #endif |
26911512 |
83 | |
84 | ptrGeoMan = G4GeometryManager::GetInstance(); |
85 | if (ptrGeoMan) { |
86 | ptrGeoMan->OpenGeometry(); |
87 | |
88 | //true argoment allows voxel construction; if false voxels are built |
89 | //only for replicated volumes |
90 | ptrGeoMan->CloseGeometry(true); |
91 | } |
92 | else { |
93 | G4cerr << "ERROR in FLUGG: Could not get G4GeometryManager instance" |
94 | << G4endl; |
95 | G4cerr << " in FGeometryInit::closeGeometry. Exiting!!!" |
96 | << G4endl; |
97 | } |
98 | |
e73e0522 |
99 | #ifdef G4GEOMETRY_DEBUG |
100 | G4cout << "<== Flugg FGeometryInit::closeGeometry()" << G4endl; |
101 | #endif |
26911512 |
102 | } |
103 | |
104 | //************************************************************************* |
105 | |
106 | void FGeometryInit::InitHistArray() { |
e73e0522 |
107 | #ifdef G4GEOMETRY_DEBUG |
108 | G4cout << "==> Flugg FGeometryInit::InitHistArray()" << G4endl; |
109 | #endif |
26911512 |
110 | ptrArray = new G4int[1000000]; |
111 | for(G4int i=0;i<1000000;i++) |
112 | ptrArray[i]=0; |
e73e0522 |
113 | #ifdef G4GEOMETRY_DEBUG |
114 | G4cout << "<== Flugg FGeometryInit::InitHistArray()" << G4endl; |
115 | #endif |
26911512 |
116 | } |
117 | |
118 | |
119 | |
120 | //************************************************************************* |
121 | //jrLtGeant stores all crossed lattice volume histories. |
122 | |
123 | void FGeometryInit::InitJrLtGeantArray() { |
26911512 |
124 | #ifdef G4GEOMETRY_DEBUG |
e73e0522 |
125 | G4cout << "==> Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl; |
26911512 |
126 | G4cout << "Initializing JrLtGeant array" << G4endl; |
127 | #endif |
128 | ptrJrLtGeant = new G4int[10000]; |
129 | for(G4int x=0;x<10000;x++) |
130 | ptrJrLtGeant[x]=-1; |
131 | flagLttcGeant = -1; |
e73e0522 |
132 | #ifdef G4GEOMETRY_DEBUG |
133 | G4cout << "<== Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl; |
134 | #endif |
26911512 |
135 | } |
136 | |
137 | |
138 | void FGeometryInit::SetLttcFlagGeant(G4int newFlagLttc) { |
e73e0522 |
139 | #ifdef G4GEOMETRY_DEBUG |
140 | G4cout << "==> Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl; |
141 | #endif |
26911512 |
142 | // Added by A.Solodkov |
143 | if (newFlagLttc >= 10000) { |
144 | G4cout << "Problems in FGeometryInit::SetLttcFlagGeant" << G4endl; |
145 | G4cout << "Index newFlagLttc=" << newFlagLttc << " is outside array bounds" |
146 | << G4endl; |
147 | G4cout << "Better to stop immediately !" << G4endl; |
148 | exit(1); |
149 | } |
150 | flagLttcGeant = newFlagLttc; |
e73e0522 |
151 | #ifdef G4GEOMETRY_DEBUG |
152 | G4cout << "<== Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl; |
153 | #endif |
26911512 |
154 | } |
155 | |
156 | void FGeometryInit::PrintJrLtGeant() { |
157 | #ifdef G4GEOMETRY_DEBUG |
158 | //G4cout << "jrLtGeant:" << G4endl; |
159 | //for(G4int y=0;y<=flagLttcGeant;y++) |
160 | // |
161 | // G4cout << "jrLtGeant[" << y << "]=" << ptrJrLtGeant[y] << G4endl; |
162 | #endif |
163 | } |
164 | |
165 | //************************************************************************** |
166 | |
167 | void FGeometryInit::PrintHistories() { |
168 | /* |
169 | #ifdef G4GEOMETRY_DEBUG |
170 | G4cout << "Touch hist:" << G4endl; |
171 | G4cout << *(ptrTouchHist->GetHistory()) << G4endl; |
172 | G4cout << "Tmp hist:" << G4endl; |
173 | G4cout << *(ptrTempNavHist->GetHistory()) << G4endl; |
174 | G4cout << "Old hist:" << G4endl; |
175 | G4cout << *(ptrOldNavHist->GetHistory()) << G4endl; |
176 | #endif |
177 | */ |
178 | } |
179 | |
180 | |
181 | |
182 | void FGeometryInit::InitHistories() { |
e73e0522 |
183 | #ifdef G4GEOMETRY_DEBUG |
184 | G4cout << "==> Flugg FGeometryInit::InitHistories()" << G4endl; |
185 | #endif |
26911512 |
186 | //init utility histories with navigator history |
187 | ptrTouchHist = |
188 | fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory(); |
189 | ptrTempNavHist = |
190 | fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory(); |
191 | ptrOldNavHist = new G4TouchableHistory(); |
e73e0522 |
192 | #ifdef G4GEOMETRY_DEBUG |
193 | G4cout << "<== Flugg FGeometryInit::InitHistories()" << G4endl; |
194 | #endif |
26911512 |
195 | } |
196 | |
197 | void FGeometryInit::DeleteHistories() { |
e73e0522 |
198 | #ifdef G4GEOMETRY_DEBUG |
199 | G4cout << "==> Flugg FGeometryInit::DeleteHistories()" << G4endl; |
200 | #endif |
201 | |
26911512 |
202 | delete ptrTouchHist; |
203 | delete ptrOldNavHist; |
204 | delete ptrTempNavHist; |
205 | |
206 | #ifdef G4GEOMETRY_DEBUG |
207 | G4cout << "Deleting step-history objects at end of run!" << G4endl; |
e73e0522 |
208 | G4cout << "<== Flugg FGeometryInit::DeleteHistories()" << G4endl; |
26911512 |
209 | #endif |
26911512 |
210 | } |
211 | |
212 | |
213 | void FGeometryInit::UpdateHistories(const G4NavigationHistory * history, |
214 | G4int flagHist) { |
e73e0522 |
215 | #ifdef G4GEOMETRY_DEBUG |
216 | G4cout << "==> Flugg FGeometryInit::UpdateHistories()" << G4endl; |
217 | #endif |
26911512 |
218 | PrintHistories(); |
219 | |
220 | #ifdef G4GEOMETRY_DEBUG |
221 | G4cout << "...updating histories!" << G4endl; |
222 | #endif |
223 | |
224 | G4VPhysicalVolume * pPhysVol = history->GetTopVolume(); |
225 | |
226 | switch (flagHist) { |
227 | case 0: { |
228 | //this is the case when a new history is given to the |
229 | //navigator and old history has to be resetted |
230 | //touchable history has not been updated jet, so: |
231 | |
232 | ptrTouchHist->UpdateYourself(pPhysVol,history); |
233 | ptrTempNavHist->UpdateYourself(pPhysVol,history); |
234 | G4NavigationHistory * ptrOldNavHistNotConst = |
235 | const_cast<G4NavigationHistory * >(ptrOldNavHist->GetHistory()); |
236 | ptrOldNavHistNotConst->Reset(); |
237 | ptrOldNavHistNotConst->Clear(); |
238 | PrintHistories(); |
239 | break; |
240 | } //case 0 |
241 | |
242 | case 1: { |
243 | //this is the case when a new history is given to the |
244 | //navigator but old history has to be kept (e.g. LOOKZ |
245 | //is call during an event); |
246 | //touchable history has not been updated jet, so: |
247 | |
248 | ptrTouchHist->UpdateYourself(pPhysVol,history); |
249 | ptrTempNavHist->UpdateYourself(pPhysVol,history); |
250 | PrintHistories(); |
251 | break; |
252 | } //case 1 |
253 | |
254 | case 2: { |
255 | //this is the case when the touchable history has been |
256 | //updated by a LocateGlobalPointAndUpdateTouchable call |
257 | |
258 | G4VPhysicalVolume * pPhysVolTemp = ptrTempNavHist->GetVolume(); |
259 | ptrOldNavHist->UpdateYourself(pPhysVolTemp, |
260 | ptrTempNavHist->GetHistory()); |
261 | |
262 | ptrTempNavHist->UpdateYourself(pPhysVol,history); |
263 | PrintHistories(); |
264 | break; |
265 | } //case 2 |
266 | |
267 | default: { |
268 | G4cout <<" ERROR in updating step-histories!" << G4endl; |
269 | break; |
270 | } //default |
271 | } //switch |
272 | |
e73e0522 |
273 | #ifdef G4GEOMETRY_DEBUG |
274 | G4cout << "<== Flugg FGeometryInit::UpdateHistories()" << G4endl; |
275 | #endif |
26911512 |
276 | } |
277 | |
278 | //***************************************************************************** |
279 | |
280 | |
281 | void FGeometryInit::createFlukaMatFile() { |
26911512 |
282 | // last modification Sara Vanini 1/III/99 |
283 | // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case, |
284 | // according to the fluka standard. In addition,. they must be equal to the |
285 | // names of the fluka materials - see fluka manual - in order that the |
286 | // program load the right cross sections, and equal to the names included in |
287 | // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his |
288 | // own .pemf, in order to get the right cross sections loaded in memory. |
289 | |
290 | #ifdef G4GEOMETRY_DEBUG |
e73e0522 |
291 | G4cout << "==> Flugg FGeometryInit::createFlukaMatFile()" << G4endl; |
26911512 |
292 | G4cout << "================== FILEWR =================" << G4endl; |
293 | #endif |
294 | |
295 | //open file for output |
36c70081 |
296 | G4std::ofstream fout("flukaMat.inp"); |
26911512 |
297 | |
298 | //PhysicalVolumeStore, Volume and MaterialTable pointers |
299 | G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance(); |
300 | unsigned int numVol = pVolStore->size(); |
301 | |
302 | G4int* indexMatFluka = 0; |
303 | static G4Material * ptrMat = 0; |
304 | unsigned int x = 0; |
305 | |
306 | //Find a volume with a material associated |
307 | while(!ptrMat && x<numVol) { |
308 | G4VPhysicalVolume * ptrVol = (*pVolStore)[x]; |
309 | G4LogicalVolume * ptrLogVol = ptrVol->GetLogicalVolume(); |
310 | ptrMat = ptrLogVol->GetMaterial(); |
311 | x++; |
312 | } |
313 | |
314 | if(ptrMat) { |
315 | static const G4MaterialTable * ptrMatTab = G4Material::GetMaterialTable(); |
316 | |
317 | //number of materials, elements, variable initialisations |
318 | static size_t totNumMat = G4Material::GetNumberOfMaterials(); |
319 | #ifdef G4GEOMETRY_DEBUG |
320 | G4cout << "Number of materials: " << totNumMat << G4endl; |
321 | #endif |
322 | |
323 | static size_t totNumElem = G4Element::GetNumberOfElements(); |
324 | #ifdef G4GEOMETRY_DEBUG |
325 | G4cout << "Number of elements: " << totNumMat << G4endl; |
326 | #endif |
327 | |
328 | |
329 | G4int * elemIndexInMATcard = new G4int[totNumElem]; |
330 | for(unsigned int t=0; t<totNumElem; t++) |
331 | elemIndexInMATcard[t] = 0; |
332 | static const G4IsotopeTable * ptrIsotTab = 0; |
333 | G4int initIsot = 0; |
334 | static size_t totNumIsot; |
335 | G4int* isotIndexInMATcard = 0; |
336 | G4int isotPresence = 0; |
337 | |
338 | |
339 | // title |
340 | PrintHeader(fout,"GEANT4 ELEMENTS AND COMPOUNDS"); |
341 | |
342 | // *** loop over G4Materials to assign Fluka index |
343 | G4int indexCount = 3; |
344 | indexMatFluka = new G4int[totNumMat]; |
345 | for(unsigned int i=0; i<totNumMat; i++) { |
346 | //pointer material, state |
347 | ptrMat = (*ptrMatTab)[i]; |
348 | G4double denMat = ptrMat->GetDensity(); |
349 | G4String nameMat = ptrMat->GetName(); |
350 | // Fluka index: bh=1, vacuum=2, others=3,4.. |
351 | if(denMat<=1.00e-10*g/cm3) |
352 | //N.B. fluka density limit decided on XI-`98 |
353 | indexMatFluka[i]=2; |
354 | else { |
355 | indexMatFluka[i]=indexCount; |
356 | indexCount++; |
357 | } |
358 | |
359 | // *** write single-element material MATERIAL card |
360 | size_t numElem = ptrMat->GetNumberOfElements(); |
361 | if(numElem==1) { |
362 | G4int index = indexMatFluka[i]; |
363 | const G4Element * ptrElem = ptrMat->GetElement(0); |
364 | size_t indElemTab = ptrElem->GetIndex(); |
365 | size_t numIsot = ptrElem->GetNumberOfIsotopes(); |
366 | G4double A = (ptrElem->GetA())/(g); |
367 | if(!numIsot) { |
368 | if(index!=2 && !elemIndexInMATcard[indElemTab]) { |
369 | G4double Z = ptrElem->GetZ(); |
370 | elemIndexInMATcard[indElemTab] = index; |
371 | G4String nameEl = ptrElem->GetName(); |
372 | nameEl.toUpper(); |
373 | |
374 | //write on file MATERIAL card of element |
375 | PrintMaterial(fout, "MATERIAL ", |
376 | Z, A, |
377 | denMat/(g/cm3), |
378 | index, |
379 | -1, |
380 | nameEl); |
381 | } |
382 | } |
383 | if(numIsot==1) { |
384 | // G4Isotope pointer |
385 | const G4Isotope* ptrIsot = ptrElem->GetIsotope(0); |
386 | size_t indIsotTab = ptrIsot->GetIndex(); |
387 | //initialize variables |
388 | if(!initIsot) { |
389 | totNumIsot = ptrIsot->GetNumberOfIsotopes(); |
390 | ptrIsotTab = ptrIsot->GetIsotopeTable(); |
391 | isotIndexInMATcard = new G4int[totNumIsot]; |
392 | for(unsigned int t=0; t<totNumIsot; t++) |
393 | isotIndexInMATcard[t] = 0; |
394 | initIsot = 1; |
395 | } |
396 | if(!isotIndexInMATcard[indIsotTab]) { |
397 | //compute physical data and counters |
398 | G4int ZIs = ptrIsot->GetZ(); |
399 | G4double AIs = (ptrIsot->GetA())/(g); |
400 | G4int NIs = ptrIsot->GetN(); |
401 | G4String nameIsot = ptrIsot->GetName(); |
402 | nameIsot.toUpper(); |
403 | isotIndexInMATcard[indIsotTab] = index; |
404 | |
405 | //write on file MATERIAL card of isotope |
406 | PrintMaterial(fout, "MATERIAL ", |
407 | G4double(ZIs), |
408 | AIs, |
409 | denMat/(g/cm3), |
410 | G4double(index), |
411 | G4double(NIs), |
412 | nameIsot); |
413 | } |
414 | }// end if(numIsot==1) |
415 | }// end if(numElem==1) |
416 | }// end for loop |
417 | |
418 | |
419 | // *** material definitions: elements, compound made of G4Elements |
420 | // or made of G4Materials |
421 | for(unsigned int j=0; j<totNumMat; j++) { |
422 | //pointer material, and material data |
423 | ptrMat = (*ptrMatTab)[j]; |
424 | size_t numElem = ptrMat->GetNumberOfElements(); |
425 | G4double densityMat = (ptrMat->GetDensity())/(g/cm3); |
426 | G4String nameMat = ptrMat->GetName(); |
427 | nameMat.toUpper(); |
428 | isotPresence = 0; |
429 | |
430 | //fraction vector of compounds of the material |
431 | const G4double* ptrFracVect = ptrMat->GetFractionVector(); |
432 | |
433 | //loop on elements of the material |
434 | for(unsigned int el=0; el<numElem; el++) { |
435 | //compute physical data, initialize variables |
436 | const G4Element * ptrElem = ptrMat->GetElement(el); |
437 | size_t indElemTab = ptrElem->GetIndex(); |
438 | G4String nameElem = ptrElem->GetName(); |
439 | nameElem.toUpper(); |
440 | size_t numIsot = ptrElem->GetNumberOfIsotopes(); |
441 | G4double A = (ptrElem->GetA())/(g); |
442 | |
443 | if(!numIsot) { |
444 | if(!elemIndexInMATcard[indElemTab]) { |
445 | G4double Z = ptrElem->GetZ(); |
446 | G4double density = ptrFracVect[el]*densityMat; |
447 | elemIndexInMATcard[indElemTab] = indexCount; |
448 | |
449 | //write on file MATERIAL card of element |
450 | PrintMaterial(fout, "MATERIAL ", |
451 | Z, A, |
452 | density, |
453 | G4double(indexCount), |
454 | -1, |
455 | nameElem); |
456 | } |
457 | } |
458 | |
459 | else { |
460 | if(numIsot>=2) |
461 | isotPresence = 1; |
462 | //loop on isotopes |
463 | for(unsigned int nis=0; nis<numIsot; nis++) { |
464 | // G4Isotope pointer |
465 | const G4Isotope* ptrIsot = ptrElem->GetIsotope(nis); |
466 | size_t indIsotTab = ptrIsot->GetIndex(); |
467 | //initialize variables |
468 | if(!initIsot) { |
469 | totNumIsot = ptrIsot->GetNumberOfIsotopes(); |
470 | ptrIsotTab = ptrIsot->GetIsotopeTable(); |
471 | isotIndexInMATcard = new G4int[totNumIsot]; |
472 | for(unsigned int t=0; t<totNumIsot; t++) |
473 | isotIndexInMATcard[t] = 0; |
474 | initIsot = 1; |
475 | } |
476 | if(!isotIndexInMATcard[indIsotTab]) { |
477 | //compute physical data and counters |
478 | G4int ZIs = ptrIsot->GetZ(); |
479 | G4double AIs = (ptrIsot->GetA())/(g); |
480 | G4int NIs = ptrIsot->GetN(); |
481 | G4String nameIsot = ptrIsot->GetName(); |
482 | nameIsot.toUpper(); |
483 | G4double* ptrRelAbVect = ptrElem->GetRelativeAbundanceVector(); |
484 | G4double density = |
485 | ptrFracVect[el]*densityMat*ptrRelAbVect[nis]*AIs/A; |
486 | G4int index = indexCount; |
487 | isotIndexInMATcard[indIsotTab] = indexCount; |
488 | indexCount+=1; |
489 | |
490 | //write on file MATERIAL card of isotope |
491 | PrintMaterial(fout, "MATERIAL ", |
492 | G4double(ZIs), AIs, |
493 | density, |
494 | G4double(index), |
495 | NIs, |
496 | nameIsot); |
497 | } |
498 | } |
499 | } |
500 | |
501 | } |
502 | |
503 | if(numElem>1 || isotPresence==1) { |
504 | // write MATERIAL+COMPOUND card specifing the compound |
505 | |
506 | //flags for writing COMPOUND card |
507 | G4int treCount=0; |
508 | |
509 | //make MATERIAL card for compound, start COMPOUND card |
510 | fout <<"*"<<G4endl; |
511 | fout <<"* Define GEANT4 compound " << nameMat << G4endl; |
512 | PrintMaterial(fout, "MATERIAL ", |
513 | -1, -1, |
514 | densityMat, |
515 | G4double(indexMatFluka[j]), |
516 | -1, |
517 | nameMat); |
518 | fout << setw10 << "COMPOUND "; |
519 | |
520 | |
521 | //write elements in COMPOUND card |
522 | for(unsigned int h=0; h<numElem; h++) { |
523 | const G4Element * ptrElemMat = ptrMat->GetElement(h); |
524 | size_t indexElemMat = ptrElemMat->GetIndex(); |
525 | size_t numIsotElem = ptrElemMat->GetNumberOfIsotopes(); |
526 | if(!numIsotElem) { |
527 | PrintCompound(fout, "COMPOUND ", |
528 | treCount, |
529 | nameMat, |
530 | -ptrFracVect[h], |
531 | G4double(elemIndexInMATcard[indexElemMat])); |
532 | |
533 | if(treCount==3) |
534 | treCount=0; |
535 | treCount+=1; |
536 | } |
537 | else { |
538 | G4double * ptrIsotAbbVect = |
539 | ptrElemMat->GetRelativeAbundanceVector(); |
540 | |
541 | for(unsigned int iso=0; iso<numIsotElem; iso++) { |
542 | const G4Isotope * ptrIsotElem =ptrElemMat->GetIsotope(iso); |
543 | size_t indexIsotMat = ptrIsotElem->GetIndex(); |
544 | G4double isotAbundPerVol = |
545 | ptrIsotAbbVect[iso]*Avogadro*densityMat* |
546 | ptrFracVect[h]/(ptrElemMat->GetA()/(g)); |
547 | |
548 | PrintCompound(fout, "COMPOUND ", |
549 | treCount, |
550 | nameMat, |
551 | isotAbundPerVol, |
552 | G4double(isotIndexInMATcard[indexIsotMat])); |
553 | if(treCount==3) |
554 | treCount=0; |
555 | treCount+=1; |
556 | } |
557 | } |
558 | } |
559 | |
560 | //end COMPOUND card |
561 | if(treCount==1) |
562 | fout << setw10 << " " << setw10 << " " |
563 | << setw10 << " " << setw10 << " " |
e73e0522 |
564 | << nameMat << G4endl; |
26911512 |
565 | if(treCount==2) |
566 | fout << setw10 << " " << setw10 << " " |
e73e0522 |
567 | << nameMat << G4endl; |
26911512 |
568 | if(treCount==3) |
e73e0522 |
569 | fout << nameMat << G4endl; |
570 | fout << "*" << G4endl; |
26911512 |
571 | } |
572 | |
573 | |
574 | } // end for loop |
575 | delete elemIndexInMATcard; |
576 | if(initIsot) |
577 | delete isotIndexInMATcard; |
578 | |
579 | } // end if (ptrMat) |
580 | |
581 | // *** material-volume correspondence |
582 | PrintHeader(fout,"GEANT4 MATERIAL ASSIGNMENTS"); |
583 | |
584 | //initializations |
585 | G4int indexMatOld = 0; |
586 | G4int indexRegFlukaFrom = 0; |
587 | G4int indexRegFlukaTo = 0; |
588 | G4int existTo = 0; |
589 | G4int flagField = 0; |
590 | G4int lastFlagField = 0; |
591 | |
592 | //open file for volume-index correspondence |
36c70081 |
593 | G4std::ofstream vout("Volumes_index.inp"); |
26911512 |
594 | |
595 | //... and write title |
e73e0522 |
596 | vout << "*" << G4endl; |
26911512 |
597 | vout << "******************** GEANT4 VOLUMES *******************\n"; |
e73e0522 |
598 | vout << "*" << G4endl; |
26911512 |
599 | |
600 | //loop over all volumes... |
601 | for(unsigned int l=0;l<numVol;l++) { |
602 | //index of the region |
603 | G4VPhysicalVolume * ptrVol = (*pVolStore)[l]; |
604 | G4LogicalVolume * ptrLogVol = ptrVol->GetLogicalVolume(); |
605 | G4int indexRegFluka = l+1; |
606 | |
607 | |
608 | //write index volume and name on file Volumes_index.inp |
609 | vout.setf(G4std::ios::left,G4std::ios::adjustfield); |
610 | vout << setw10 << indexRegFluka; |
611 | vout << G4std::setw(20) << ptrVol->GetName() << G4std::setw(20) << ""; |
612 | if(ptrVol->IsReplicated()) { |
613 | EAxis axis; |
614 | G4int nRep; |
615 | G4double width; |
616 | G4double offset; |
617 | G4bool consum; |
618 | ptrVol->GetReplicationData(axis,nRep,width,offset,consum); |
619 | vout.setf(G4std::ios::left,G4std::ios::adjustfield); |
620 | vout << setw10 << "Repetion Nb: " << G4std::setw(3) << nRep; |
621 | } |
e73e0522 |
622 | vout << G4endl; |
26911512 |
623 | |
624 | //check if Magnetic Field is present in the region |
625 | G4FieldManager * pMagFieldMan = ptrLogVol->GetFieldManager(); |
626 | const G4Field * pMagField = 0; |
627 | if(pMagFieldMan) |
628 | pMagField = pMagFieldMan->GetDetectorField(); |
629 | if(pMagField) |
630 | flagField = 1; |
631 | else |
632 | flagField = 0; |
633 | |
634 | |
635 | //index of material in the region |
636 | G4Material * ptrMat = ptrLogVol->GetMaterial(); |
637 | G4int indexMat; |
638 | if(ptrMat) { |
639 | size_t indexMatGeant = ptrMat->GetIndex(); |
640 | indexMat = indexMatFluka[indexMatGeant]; |
641 | } |
642 | else |
643 | indexMat = 2; |
644 | |
645 | //if materials are repeated |
646 | if(indexMat==indexMatOld && flagField==lastFlagField) { |
647 | indexRegFlukaTo=indexRegFluka; |
648 | existTo=1; |
649 | if((l+1)==numVol) { |
650 | fout << setw10 << G4double(indexRegFlukaTo); |
651 | fout << setw10 << "0.0"; |
652 | fout << setw10 << G4double(flagField); |
e73e0522 |
653 | fout << setw10 << "0.0" << G4endl; |
26911512 |
654 | } |
655 | |
656 | } |
657 | else { |
658 | //write on file ASSIGNMAT card |
659 | |
660 | //first complete last material card |
661 | if(!existTo) { |
662 | if(l) { |
663 | fout << setw10 << "0.0"; |
664 | fout << setw10 << "0.0"; |
665 | fout << setw10 << G4double(lastFlagField); |
e73e0522 |
666 | fout << setw10 << "0.0" << G4endl; |
26911512 |
667 | } |
668 | } |
669 | else { |
670 | fout << setw10 << G4double(indexRegFlukaTo); |
671 | fout << setw10 << "0.0"; |
672 | fout << setw10 << G4double(lastFlagField); |
e73e0522 |
673 | fout << setw10 << "0.0" << G4endl; |
26911512 |
674 | } |
675 | |
676 | // begin material card |
677 | fout << setw10 << "ASSIGNMAT "; |
36c70081 |
678 | fout.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield); |
26911512 |
679 | fout << setw10 << setfixed |
680 | << G4std::setprecision(1) << G4double(indexMat); |
681 | fout << setw10 << G4double(indexRegFluka); |
682 | |
683 | existTo=0; |
684 | indexRegFlukaFrom=indexRegFluka; |
685 | |
686 | if((l+1)==numVol) { |
687 | fout << setw10 << "0.0"; |
688 | fout << setw10 << "0.0"; |
689 | fout << setw10 << G4double(flagField); |
e73e0522 |
690 | fout << setw10 << "0.0" << G4endl; |
26911512 |
691 | } |
692 | } |
693 | lastFlagField = flagField; |
694 | indexMatOld = indexMat; |
695 | } // end of loop |
696 | |
697 | //assign material 1 to black-hole=n.vol+1 |
698 | fout << setw10 << "ASSIGNMAT "; |
699 | fout << setw10 << "1.0"; |
700 | fout << setw10 << G4double(numVol+1); |
701 | fout << setw10 << "0.0"; |
702 | fout << setw10 << "0.0"; |
703 | fout << setw10 << "0.0"; |
e73e0522 |
704 | fout << setw10 << "0.0" << G4endl; |
26911512 |
705 | |
706 | // *** magnetic field *** |
707 | if(fTransportationManager->GetFieldManager()->DoesFieldExist()) { |
708 | PrintHeader(fout,"GEANT4 MAGNETIC FIELD"); |
709 | |
710 | //get magnetic field pointer |
711 | const G4Field * pMagField = |
712 | fTransportationManager->GetFieldManager()->GetDetectorField(); |
713 | |
714 | //if uniform magnetic field, get value |
715 | G4double Bx=0.0; |
716 | G4double By=0.0; |
717 | G4double Bz=0.0; |
718 | |
719 | if(pMagField) { |
720 | //G4ThreeVector* Field = new G4ThreeVector(1.,2.,3.); |
721 | const G4UniformMagField *pUnifMagField = |
722 | dynamic_cast<const G4UniformMagField*>(pMagField); |
723 | if(pUnifMagField) { |
724 | G4double *pFieldValue = 0; |
725 | G4double *point = new G4double[3]; |
726 | point[0] = 0.; |
727 | point[1] = 0.; |
728 | point[2] = 0.; |
729 | pUnifMagField->GetFieldValue(point,pFieldValue); |
730 | //non capisco perche' l'instruzione seguente non fa linkare. Indaga!! |
731 | //per ora posso usare GetFieldValue: va bene lo stesso. |
732 | //G4ThreeVector FieldValue = pUnifMagField->GetConstantFieldValue(); |
733 | Bx = pFieldValue[0]; |
734 | By = pFieldValue[1]; |
735 | Bz = pFieldValue[2]; |
736 | } |
737 | |
738 | } |
739 | |
740 | //write MGNFIELD card |
741 | fout << setw10 << "MGNFIELD "; |
742 | fout << setw10 << ""; |
743 | fout << setw10 << ""; |
744 | fout << setw10 << ""; |
36c70081 |
745 | fout.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield); |
26911512 |
746 | fout << setw10 << setfixed |
747 | << G4std::setprecision(4) << Bx |
748 | << setw10 << By |
749 | << setw10 << Bz |
e73e0522 |
750 | << G4endl; |
26911512 |
751 | } // end if magnetic field |
752 | |
753 | vout.close(); |
754 | fout.close(); |
755 | delete [] indexMatFluka; |
e73e0522 |
756 | #ifdef G4GEOMETRY_DEBUG |
757 | G4cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << G4endl; |
758 | #endif |
26911512 |
759 | } |