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