Methods needed for TFluka added.
[u/mrichter/AliRoot.git] / Flugg / FGeometryInit.cxx
1
2 // Flugg tag 
3 // modified 17/06/02: by I. Gonzalez. STL migration
4
5 //#include <stdio.h>
6 //#include <iomanip.h>
7 #include "FGeometryInit.hh"
8 #include "FluggNavigator.hh"
9 #include "WrapUtils.hh"
10 #include "FlukaMaterial.hh"
11 #include "FlukaCompound.hh"
12
13 FGeometryInit * FGeometryInit::flagInstance=0;
14
15 FGeometryInit* FGeometryInit::GetInstance()  {
16 #ifdef G4GEOMETRY_DEBUG
17   G4cout << "==> Flugg::FGeometryInit::GetInstance(), instance=" 
18          << flagInstance << G4endl;
19 #endif
20   if (!flagInstance) 
21     flagInstance = new FGeometryInit();
22   
23 #ifdef G4GEOMETRY_DEBUG
24   G4cout << "<== Flugg::FGeometryInit::GetInstance(), instance=" 
25          << flagInstance << G4endl;
26 #endif
27   return flagInstance;
28 }  
29
30
31 FGeometryInit::FGeometryInit():
32   fDetector(0),
33   fFieldManager(0),
34   fTransportationManager(G4TransportationManager::GetTransportationManager()),
35   myTopNode(0),
36   ptrGeoMan(0),
37   ptrArray(0),
38   ptrTouchHist(0),
39   ptrOldNavHist(0),
40   ptrTempNavHist(0),
41   ptrJrLtGeant(0)
42 {
43
44 #ifdef G4GEOMETRY_DEBUG
45   G4cout << "==> Flugg FGeometryInit::FGeometryInit()" << G4endl;
46   G4cout << "\t+ Changing the G4Navigator for FluggNavigator..." << G4endl;
47 #endif
48   G4Navigator* actualnav = fTransportationManager->GetNavigatorForTracking();
49   if (actualnav) {
50     FluggNavigator* newnav = new FluggNavigator();
51     fTransportationManager->SetNavigatorForTracking(newnav);
52   }
53   else {
54     G4cerr << "ERROR: Could not find the actual G4Navigator" << G4endl;
55     abort();
56   }
57
58   
59 #ifdef G4GEOMETRY_DEBUG
60   G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
61 #endif
62
63 }
64
65
66 FGeometryInit::~FGeometryInit() {
67   G4cout << "==> Flugg FGeometryInit::~FGeometryInit()" << G4endl;
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!
77   G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
78 }
79
80
81 void FGeometryInit::closeGeometry() {
82 #ifdef G4GEOMETRY_DEBUG
83   G4cout << "==> Flugg FGeometryInit::closeGeometry()" << G4endl;
84 #endif
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
101 #ifdef G4GEOMETRY_DEBUG
102   G4cout << "<== Flugg FGeometryInit::closeGeometry()" << G4endl;
103 #endif
104 }
105
106 //*************************************************************************
107
108 void FGeometryInit::InitHistArray() {
109 #ifdef G4GEOMETRY_DEBUG
110   G4cout << "==> Flugg FGeometryInit::InitHistArray()" << G4endl;
111 #endif
112   ptrArray = new G4int[1000000];
113   for(G4int i=0;i<1000000;i++) 
114     ptrArray[i]=0;
115 #ifdef G4GEOMETRY_DEBUG
116   G4cout << "<== Flugg FGeometryInit::InitHistArray()" << G4endl;
117 #endif
118 }
119
120
121
122 //*************************************************************************
123 //jrLtGeant stores all crossed lattice volume histories.
124
125 void FGeometryInit::InitJrLtGeantArray() {
126 #ifdef G4GEOMETRY_DEBUG
127   G4cout << "==> Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
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;
134 #ifdef G4GEOMETRY_DEBUG
135   G4cout << "<== Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
136 #endif
137 }
138
139
140 void FGeometryInit::SetLttcFlagGeant(G4int newFlagLttc) {
141 #ifdef G4GEOMETRY_DEBUG
142   G4cout << "==> Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
143 #endif
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;
153 #ifdef G4GEOMETRY_DEBUG
154   G4cout << "<== Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
155 #endif
156 }
157
158 void FGeometryInit::PrintJrLtGeant() {
159 #ifdef G4GEOMETRY_DEBUG
160   //G4cout << "jrLtGeant:" << G4endl;
161   //for(G4int y=0;y<=flagLttcGeant;y++)
162   //
163   //     G4cout << "jrLtGeant[" << y << "]=" << ptrJrLtGeant[y] << G4endl;
164 #endif
165 }
166
167 //**************************************************************************
168
169 void FGeometryInit::PrintHistories() {
170   /*
171     #ifdef G4GEOMETRY_DEBUG
172     G4cout << "Touch hist:" << G4endl;
173     G4cout << *(ptrTouchHist->GetHistory()) << G4endl;
174     G4cout << "Tmp hist:" << G4endl;
175     G4cout << *(ptrTempNavHist->GetHistory()) << G4endl;
176     G4cout << "Old hist:" << G4endl;
177     G4cout << *(ptrOldNavHist->GetHistory()) << G4endl;
178     #endif
179  */
180 }
181
182
183
184 void FGeometryInit::InitHistories() {  
185 #ifdef G4GEOMETRY_DEBUG
186   G4cout << "==> Flugg FGeometryInit::InitHistories()" << G4endl;
187 #endif
188   //init utility histories with navigator history
189   ptrTouchHist = 
190     fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
191   ptrTempNavHist = 
192     fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();   
193   ptrOldNavHist = new G4TouchableHistory();
194 #ifdef G4GEOMETRY_DEBUG
195   G4cout << "<== Flugg FGeometryInit::InitHistories()" << G4endl;
196 #endif
197 }
198
199 void FGeometryInit::DeleteHistories() {
200 #ifdef G4GEOMETRY_DEBUG
201   G4cout << "==> Flugg FGeometryInit::DeleteHistories()" << G4endl;
202 #endif
203
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;
210   G4cout << "<== Flugg FGeometryInit::DeleteHistories()" << G4endl;
211 #endif
212 }
213
214
215 void FGeometryInit::UpdateHistories(const G4NavigationHistory * history,
216                                     G4int flagHist) {
217 #ifdef G4GEOMETRY_DEBUG
218   G4cout << "==> Flugg FGeometryInit::UpdateHistories()" << G4endl;
219 #endif
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   
275 #ifdef G4GEOMETRY_DEBUG
276   G4cout << "<== Flugg FGeometryInit::UpdateHistories()" << G4endl;
277 #endif
278 }
279
280 //*****************************************************************************
281
282 void FGeometryInit::createFlukaMatFile() {
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
292   G4cout << "==> Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
293   G4cout << "================== FILEWR =================" << G4endl;
294 #endif 
295
296
297   //Regions map
298   BuildRegionsMap();
299   G4std::ofstream vos("Volumes_index.inp");
300   PrintRegionsMap(vos);
301   vos.close();
302
303   //Materials and compounds
304   BuildMaterialTables();
305   G4std::ofstream fos("flukaMat.inp");  
306   PrintMaterialTables(fos);
307   PrintAssignmat(fos);
308   PrintMagneticField(fos);
309   fos.close();
310
311 #ifdef G4GEOMETRY_DEBUG
312   G4cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
313 #endif
314 }
315
316 ////////////////////////////////////////////////////////////////////////
317 // 
318 void FGeometryInit::BuildRegionsMap() {
319 #ifdef G4GEOMETRY_DEBUG
320   G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
321 #endif
322
323   //Find number of Volumes in physical volume store
324   G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
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;
336   }
337
338
339
340 #ifdef G4GEOMETRY_DEBUG
341   G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
342 #endif
343 }
344
345 void FGeometryInit::PrintRegionsMap(G4std::ostream& os) {
346 #ifdef G4GEOMETRY_DEBUG
347   G4cout << "==> Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
348 #endif
349
350   //Print some header
351   PrintHeader(os, "GEANT4 VOLUMES");
352
353   //Iterate over all volumes in the map
354   for (RegionIterator i = fRegionVolumeMap.begin(); 
355        i != fRegionVolumeMap.end(); 
356        i++) {
357     
358     //Get info in the map
359     G4VPhysicalVolume* ptrVol = (*i).first;
360     int index = (*i).second;
361
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
368     if(ptrVol->IsReplicated()) {
369       EAxis axis;
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;
377     }
378     os << G4endl;
379     
380   }
381   
382 #ifdef G4GEOMETRY_DEBUG
383   G4cout << "<== Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
384 #endif
385 }
386
387 ////////////////////////////////////////////////////////////////////////
388 //
389 void    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
405 G4int FGeometryInit::GetMedium(int region) const
406 {
407     return fRegionMediumMap[region];
408 }
409
410
411 void  FGeometryInit::SetMediumFromName(const char* volName, int medium, int volid) 
412  {
413   char name4[5];
414   char tmp[5];
415   strncpy(tmp, volName, 4);
416   tmp[4]='\0';
417   fNRegions = 0;
418   
419   for (RegionIterator i = fRegionVolumeMap.begin(); 
420        i != fRegionVolumeMap.end(); 
421        i++) {
422     fNRegions++;
423     //Get info in the map
424     G4VPhysicalVolume* ptrVol = (*i).first;
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)) {
436         fMediumVolumeMap[ptrVol] = medium;
437         fVolIdVolumeMap[ptrVol]  = volid;
438     }
439   }
440 }
441
442
443
444 ////////////////////////////////////////////////////////////////////////
445 // 
446 void FGeometryInit::BuildMaterialTables() {
447 #ifdef G4GEOMETRY_DEBUG
448   G4cout << "==> Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
449 #endif
450
451   //some terminal printout also
452   G4cout << "\t* Storing information..." << G4endl;
453
454   //The logic is the folloing:
455   //Get the Material Table and:
456   // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
457   // 2) For each single element material build a material equivalent
458   // 3) For the rest:
459   //   3.a) Build materials for each not already known element
460   //   3.b) Build the compound out of them
461
462   //Get the Material Table and iterate
463   const G4MaterialTable* matTable = G4Material::GetMaterialTable();
464   for (MatTableIterator i = matTable->begin(); i != matTable->end(); i++) {
465
466     //Get some basic material information
467     G4Material* material = (*i);
468     G4String matName = material->GetName();
469     const G4double matDensity = material->GetDensity();
470     const G4int nMatElements  = material->GetNumberOfElements();
471
472     G4cout << "\t\t+ " << matName 
473            << ": dens. = " << matDensity/(g/cm3) << "g/cm3"
474            << ", nElem = " << nMatElements << G4endl;
475
476     // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
477     //    FlukaMaterial* is 0 in that case
478     if (matDensity <= 1.00e-10*g/cm3) {
479       G4FlukaMaterialMap[material] = 0;
480       G4cout << "\t\t  Stored as vacuum" << G4endl;
481     }
482     // 2) For each single element material build a material equivalent
483     else if (nMatElements == 1) {
484       
485       FlukaMaterial *flukamat = 
486         BuildFlukaMaterialFromElement(material->GetElement(0),
487                                       matDensity);
488       
489       G4FlukaMaterialMap[material] = flukamat;
490       G4cout << "\t\t  Stored as " << flukamat->GetRealName() << G4endl;
491       
492     } //else if (material->GetNumberOfElements() == 1)
493     
494     // 3) For the rest:
495     //   3.a) Build materials for each not already known element
496     //   3.b) Build the compound out of them
497     else {
498       FlukaCompound* flukacomp = 
499         BuildFlukaCompoundFromMaterial(material);
500       G4FlukaCompoundMap[material] = flukacomp;
501       G4cout << "\t\t  Stored as " << flukacomp->GetRealName() << G4endl;
502     } //else for case 3)
503   } //for (materials)
504   
505 #ifdef G4GEOMETRY_DEBUG
506   G4cout << "<== Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
507 #endif
508 }
509
510 FlukaMaterial* 
511 FGeometryInit::BuildFlukaMaterialFromElement(const G4Element* element,
512                                              G4double matDensity) {
513 #ifdef G4GEOMETRY_DEBUG
514   G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromElement()" 
515          << G4endl;
516 #endif
517
518   //Get element and its properties
519   G4String elemName(ToFlukaString(element->GetName()));
520
521   FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(elemName);
522   if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
523     //Check for isotopes
524     G4int nIsotopes = element->GetNumberOfIsotopes();
525     if (nIsotopes == 0) {
526       G4double elemA = element->GetA()/g;
527       G4double elemZ = element->GetZ();
528       
529       if (elemA != G4int(elemA) && elemZ != G4int(elemZ)) {
530         G4cout << "WARNING: Element \'" << elemName 
531                << "\' has non integer Z (" << elemZ << ") or A (" 
532                << elemA << ")"
533                << G4endl;
534       }
535     
536       flukamat = new FlukaMaterial(elemName,
537                                    G4int(elemZ),
538                                    elemA,
539                                    matDensity/(g/cm3));
540     }
541     else if (nIsotopes == 1) {
542       const G4Isotope* isotope = element->GetIsotope(0);
543       flukamat = BuildFlukaMaterialFromIsotope(isotope, matDensity);
544     }
545     else {
546       FlukaCompound *flucomp = BuildFlukaCompoundFromElement(element,
547                                                              matDensity);
548       flukamat = flucomp->GetFlukaMaterial();      
549     }
550   }
551 #ifdef G4GEOMETRY_DEBUG
552   else {
553     G4cout << "INFO: Element \'" << elemName 
554            << "\' already exists in the DB. It will not be recreated."
555            << G4endl;
556   }
557 #endif
558
559   return flukamat;
560   
561 #ifdef G4GEOMETRY_DEBUG
562   G4cout << "<== Flugg FGeometryInit::BuildFlukaMaterialFromElement()" 
563          << G4endl;
564 #endif
565 }
566
567 FlukaMaterial* 
568 FGeometryInit::BuildFlukaMaterialFromIsotope(const G4Isotope* isotope,
569                                              G4double matDensity) {
570 #ifdef G4GEOMETRY_DEBUG
571   G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()" 
572          << G4endl;
573 #endif
574   G4String isoName(ToFlukaString(isotope->GetName()));
575   FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(isoName);
576   if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
577     G4int isoZ = isotope->GetZ();
578     G4double isoA = (isotope->GetA())/(g);
579     G4int isoN = isotope->GetN();
580     flukamat = new FlukaMaterial(isoName,
581                                  isoZ,
582                                  isoA,
583                                  matDensity/(g/cm3),
584                                  isoN);
585   }
586
587   return flukamat;
588
589 #ifdef G4GEOMETRY_DEBUG
590   G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()" 
591          << G4endl;
592 #endif
593 }
594
595 FlukaCompound* 
596 FGeometryInit::BuildFlukaCompoundFromMaterial(const G4Material* material) {
597 #ifdef G4GEOMETRY_DEBUG
598   G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()" 
599          << G4endl;
600 #endif
601   //Material properties
602   const G4double* elemFractions = material->GetFractionVector();
603   const G4int nMatElements  = material->GetNumberOfElements();
604   const G4double matDensity = material->GetDensity();
605   G4String matName(ToFlukaString(material->GetName()));
606   FlukaCompound* flukacomp = new FlukaCompound(matName, matDensity/(g/cm3),
607                                                nMatElements);
608   for (G4int i = 0; i < nMatElements; i++) {
609     FlukaMaterial *flukamat = 
610       BuildFlukaMaterialFromElement(material->GetElement(i), 0.0);
611     
612     flukacomp->AddElement(flukamat->GetIndex(), -elemFractions[i]);
613     
614   } //for (elements)
615
616   return flukacomp;
617
618 #ifdef G4GEOMETRY_DEBUG
619   G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()" 
620          << G4endl;
621 #endif
622 }
623
624 FlukaCompound* 
625 FGeometryInit::BuildFlukaCompoundFromElement(const G4Element* element,
626                                              G4double matDensity) {
627 #ifdef G4GEOMETRY_DEBUG
628   G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromElement()" 
629          << G4endl;
630 #endif
631   G4int nIsotopes = element->GetNumberOfIsotopes();
632   //fraction of nb of atomes per volume (= volume fraction?)
633   const G4double* isoAbundance =  element->GetRelativeAbundanceVector();
634   G4String elemName(ToFlukaString(element->GetName()));
635
636   //Material properties
637   FlukaCompound* flukacomp = new FlukaCompound(elemName, matDensity/(g/cm3),
638                                                nIsotopes);
639   for (G4int i = 0; i < nIsotopes; i++) {
640     FlukaMaterial *flukamat = 
641       BuildFlukaMaterialFromIsotope(element->GetIsotope(i), 0.0);
642     
643     flukacomp->AddElement(flukamat->GetIndex(), isoAbundance[i]);
644     
645   } //for (elements)
646
647   return flukacomp;
648
649 #ifdef G4GEOMETRY_DEBUG
650   G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromElement()" 
651          << G4endl;
652 #endif
653 }
654
655 void FGeometryInit::PrintMaterialTables(G4std::ostream& os) {
656 #ifdef G4GEOMETRY_DEBUG
657   G4cout << "==> Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
658 #endif
659   //Print Header
660   PrintHeader(os, "GEANT4 MATERIALS AND COMPOUNDS");
661   
662   //And some more stuff  
663   size_t nIsotopes = G4Isotope::GetNumberOfIsotopes();
664   size_t nElements = G4Element::GetNumberOfElements();
665   size_t nMaterials = G4Material::GetNumberOfMaterials();
666
667   os << "* In Geant4 there are " << nMaterials << " materials" << G4endl;
668   os << "* In Geant4 there are " << nElements  << " elements"  << G4endl;
669   os << "* In Geant4 there are " << nIsotopes  << " isotopes"  << G4endl;
670
671   //Materials
672   G4cout << "\t* Printing FLUKA materials..." << G4endl;
673   FlukaMaterial::PrintMaterialsByIndex(os);
674   //FlukaMaterial::PrintMaterialsByName(os);
675
676   //Compounds
677   G4cout << "\t* Printing FLUKA compounds..." << G4endl;
678   FlukaCompound::PrintCompounds(os);
679
680 #ifdef G4GEOMETRY_DEBUG
681   G4cout << "<== Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
682 #endif
683 }
684
685 ////////////////////////////////////////////////////////////////////////
686 // 
687 void FGeometryInit::PrintAssignmat(G4std::ostream& os) {
688 #ifdef G4GEOMETRY_DEBUG
689   G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
690 #endif
691
692   //Find number of Volumes in physical volume store
693   G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
694   unsigned int numVol = pVolStore->size();
695
696   G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore 
697          << ") has " << numVol << " volumes. " << G4endl;
698
699   G4cout << "\t* Printing ASSIGNMAT..." << G4endl;
700
701
702   PrintHeader(os,"GEANT4 MATERIAL ASSIGNMENTS");
703   for(unsigned int l=0; l < numVol; l++) {
704
705     //Get each of the physical volumes
706     G4VPhysicalVolume * physicalvol = (*pVolStore)[l];
707
708     //Get index for that volume
709     G4int iFlukaRegion = fRegionVolumeMap[physicalvol];
710
711     //Find G4 material and navigate to its fluka compound/material
712     G4LogicalVolume * logicalVol = physicalvol->GetLogicalVolume();
713     G4Material* material = logicalVol->GetMaterial();
714     G4int matIndex = 2;
715     if (G4FlukaCompoundMap[material])
716       matIndex = G4FlukaCompoundMap[material]->GetIndex();
717     if (G4FlukaMaterialMap[material])
718       matIndex = G4FlukaMaterialMap[material]->GetIndex();
719
720     //Find if there is a magnetic field in the region
721     //check if Magnetic Field is present in the region
722     G4double flagField = 0.0;
723     G4FieldManager * pMagFieldMan = logicalVol->GetFieldManager();
724     if(pMagFieldMan && pMagFieldMan->GetDetectorField())
725       flagField = 1.0;
726     
727     //Print card
728     os << setw10 << "ASSIGNMAT ";
729     os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
730     os << setw10 << setfixed << G4double(matIndex);
731     os << setw10 << setfixed << G4double(iFlukaRegion);
732     os << setw10 << "0.0";
733     os << setw10 << setfixed << flagField;
734     os << G4endl;
735   }
736
737
738
739 #ifdef G4GEOMETRY_DEBUG
740   G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
741 #endif
742 }
743
744
745 void FGeometryInit::PrintMagneticField(G4std::ostream& os) {
746 #ifdef G4GEOMETRY_DEBUG
747   G4cout << "==> Flugg FGeometryInit::PrintMagneticField()" << G4endl;
748 #endif
749
750   G4cout << "\t* Printing Magnetic Field..." << G4endl;
751
752   if(fTransportationManager->GetFieldManager()->DoesFieldExist()) {
753     
754     //get magnetic field pointer
755     const G4Field * pMagField = 
756       fTransportationManager->GetFieldManager()->GetDetectorField();     
757     
758     
759     if(pMagField) {
760       //Check if it can be made a uniform magnetic field
761       const G4UniformMagField *pUnifMagField = 
762         dynamic_cast<const G4UniformMagField*>(pMagField);
763       if(pUnifMagField) {
764         G4double B[3];
765         G4double point[4]; //it is not really used
766         pUnifMagField->GetFieldValue(point,B);
767
768         //write MGNFIELD card 
769         PrintHeader(os,"GEANT4 MAGNETIC FIELD");
770         os << setw10 << "MGNFIELD  ";
771         os << setw10 << "";
772         os << setw10 << "";
773         os << setw10 << "";
774         os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
775         os << setw10 << setfixed
776            << G4std::setprecision(4) << B[0]
777            << setw10 << B[1]
778            << setw10 << B[2]
779            << G4endl;
780       }
781       else {
782         G4cout << "WARNING: No Uniform Magnetic Field found." << G4endl;
783         G4cout << "         Manual intervention might be needed." << G4endl;
784       }
785     }
786     else
787       G4cout << "\t  No detector field found... " << G4endl;
788   } // end if magnetic field
789   else
790     G4cout << "\t  No field found... " << G4endl;
791
792 #ifdef G4GEOMETRY_DEBUG
793   G4cout << "<== Flugg FGeometryInit::PrintMagneticField()" << G4endl;
794 #endif
795 }
796
797 int FGeometryInit::CurrentVolID(int ir, int& copyNo)
798 {
799     G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
800     G4VPhysicalVolume * physicalvol = (*pVolStore)[ir- 1];
801     copyNo =  physicalvol->GetCopyNo();
802     int id = fVolIdVolumeMap[physicalvol];
803     return id;
804 }
805
806 int FGeometryInit::CurrentVolOffID(int ir, int off, int& copyNo)
807 {
808     if (off == 0) return CurrentVolID(ir, copyNo);
809     
810     G4PhysicalVolumeStore* pVolStore = G4PhysicalVolumeStore::GetInstance();
811     G4VPhysicalVolume*     physicalvol = (*pVolStore)[ir- 1];
812     G4VPhysicalVolume*     mother = physicalvol; 
813
814     int level = off;
815     while (level > 0) { 
816         if (mother) mother = mother->GetMother();
817         level--;
818     }
819     
820     int id;
821     
822     if (!mother) {
823         G4cout << "Flugg FGeometryInit::CurrentVolOffID mother not found" << G4endl;
824         id = -1;
825         copyNo = -1;
826     } else {
827         copyNo =  mother ->GetCopyNo();
828         id =  fVolIdVolumeMap[mother];
829     }
830     return id;
831 }