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