3 // modified 17/06/02: by I. Gonzalez. STL migration
7 #include "FGeometryInit.hh"
8 #include "FluggNavigator.hh"
9 #include "WrapUtils.hh"
10 #include "FlukaMaterial.hh"
11 #include "FlukaCompound.hh"
13 FGeometryInit * FGeometryInit::flagInstance=0;
15 FGeometryInit* FGeometryInit::GetInstance() {
16 #ifdef G4GEOMETRY_DEBUG
17 G4cout << "==> Flugg::FGeometryInit::GetInstance(), instance="
18 << flagInstance << G4endl;
21 flagInstance = new FGeometryInit();
23 #ifdef G4GEOMETRY_DEBUG
24 G4cout << "<== Flugg::FGeometryInit::GetInstance(), instance="
25 << flagInstance << G4endl;
31 FGeometryInit::FGeometryInit():
34 fTransportationManager(G4TransportationManager::GetTransportationManager()),
44 #ifdef G4GEOMETRY_DEBUG
45 G4cout << "==> Flugg FGeometryInit::FGeometryInit()" << G4endl;
46 G4cout << "\t+ Changing the G4Navigator for FluggNavigator..." << G4endl;
48 G4Navigator* actualnav = fTransportationManager->GetNavigatorForTracking();
50 FluggNavigator* newnav = new FluggNavigator();
51 fTransportationManager->SetNavigatorForTracking(newnav);
54 G4cerr << "ERROR: Could not find the actual G4Navigator" << G4endl;
59 #ifdef G4GEOMETRY_DEBUG
60 G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
66 FGeometryInit::~FGeometryInit() {
67 G4cout << "==> Flugg FGeometryInit::~FGeometryInit()" << G4endl;
69 ptrGeoMan->OpenGeometry();
70 if (fTransportationManager)
71 delete fTransportationManager;
73 delete[] ptrJrLtGeant;
76 //keep ATTENTION: never delete a pointer twice!
77 G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
80 void FGeometryInit::Init() {
81 // Build and initialize G4 geometry
92 void FGeometryInit::closeGeometry() {
93 #ifdef G4GEOMETRY_DEBUG
94 G4cout << "==> Flugg FGeometryInit::closeGeometry()" << G4endl;
97 ptrGeoMan = G4GeometryManager::GetInstance();
99 ptrGeoMan->OpenGeometry();
101 //true argoment allows voxel construction; if false voxels are built
102 //only for replicated volumes
103 ptrGeoMan->CloseGeometry(true);
106 G4cerr << "ERROR in FLUGG: Could not get G4GeometryManager instance"
108 G4cerr << " in FGeometryInit::closeGeometry. Exiting!!!"
112 #ifdef G4GEOMETRY_DEBUG
113 G4cout << "<== Flugg FGeometryInit::closeGeometry()" << G4endl;
117 //*************************************************************************
119 void FGeometryInit::InitHistArray() {
120 #ifdef G4GEOMETRY_DEBUG
121 G4cout << "==> Flugg FGeometryInit::InitHistArray()" << G4endl;
123 ptrArray = new G4int[1000000];
124 for(G4int i=0;i<1000000;i++)
126 #ifdef G4GEOMETRY_DEBUG
127 G4cout << "<== Flugg FGeometryInit::InitHistArray()" << G4endl;
133 //*************************************************************************
134 //jrLtGeant stores all crossed lattice volume histories.
136 void FGeometryInit::InitJrLtGeantArray() {
137 #ifdef G4GEOMETRY_DEBUG
138 G4cout << "==> Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
139 G4cout << "Initializing JrLtGeant array" << G4endl;
141 ptrJrLtGeant = new G4int[10000];
142 for(G4int x=0;x<10000;x++)
145 #ifdef G4GEOMETRY_DEBUG
146 G4cout << "<== Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
151 void FGeometryInit::SetLttcFlagGeant(G4int newFlagLttc) {
152 #ifdef G4GEOMETRY_DEBUG
153 G4cout << "==> Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
155 // Added by A.Solodkov
156 if (newFlagLttc >= 10000) {
157 G4cout << "Problems in FGeometryInit::SetLttcFlagGeant" << G4endl;
158 G4cout << "Index newFlagLttc=" << newFlagLttc << " is outside array bounds"
160 G4cout << "Better to stop immediately !" << G4endl;
163 flagLttcGeant = newFlagLttc;
164 #ifdef G4GEOMETRY_DEBUG
165 G4cout << "<== Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
169 void FGeometryInit::PrintJrLtGeant() {
170 #ifdef G4GEOMETRY_DEBUG
171 //G4cout << "jrLtGeant:" << G4endl;
172 //for(G4int y=0;y<=flagLttcGeant;y++)
174 // G4cout << "jrLtGeant[" << y << "]=" << ptrJrLtGeant[y] << G4endl;
178 //**************************************************************************
180 void FGeometryInit::PrintHistories() {
182 #ifdef G4GEOMETRY_DEBUG
183 G4cout << "Touch hist:" << G4endl;
184 G4cout << *(ptrTouchHist->GetHistory()) << G4endl;
185 G4cout << "Tmp hist:" << G4endl;
186 G4cout << *(ptrTempNavHist->GetHistory()) << G4endl;
187 G4cout << "Old hist:" << G4endl;
188 G4cout << *(ptrOldNavHist->GetHistory()) << G4endl;
195 void FGeometryInit::InitHistories() {
196 #ifdef G4GEOMETRY_DEBUG
197 G4cout << "==> Flugg FGeometryInit::InitHistories()" << G4endl;
199 //init utility histories with navigator history
201 fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
203 fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
204 ptrOldNavHist = new G4TouchableHistory();
205 #ifdef G4GEOMETRY_DEBUG
206 G4cout << "<== Flugg FGeometryInit::InitHistories()" << G4endl;
210 void FGeometryInit::DeleteHistories() {
211 #ifdef G4GEOMETRY_DEBUG
212 G4cout << "==> Flugg FGeometryInit::DeleteHistories()" << G4endl;
216 delete ptrOldNavHist;
217 delete ptrTempNavHist;
219 #ifdef G4GEOMETRY_DEBUG
220 G4cout << "Deleting step-history objects at end of run!" << G4endl;
221 G4cout << "<== Flugg FGeometryInit::DeleteHistories()" << G4endl;
226 void FGeometryInit::UpdateHistories(const G4NavigationHistory * history,
228 #ifdef G4GEOMETRY_DEBUG
229 G4cout << "==> Flugg FGeometryInit::UpdateHistories()" << G4endl;
233 #ifdef G4GEOMETRY_DEBUG
234 G4cout << "...updating histories!" << G4endl;
237 G4VPhysicalVolume * pPhysVol = history->GetTopVolume();
241 //this is the case when a new history is given to the
242 //navigator and old history has to be resetted
243 //touchable history has not been updated jet, so:
245 ptrTouchHist->UpdateYourself(pPhysVol,history);
246 ptrTempNavHist->UpdateYourself(pPhysVol,history);
247 G4NavigationHistory * ptrOldNavHistNotConst =
248 const_cast<G4NavigationHistory * >(ptrOldNavHist->GetHistory());
249 ptrOldNavHistNotConst->Reset();
250 ptrOldNavHistNotConst->Clear();
256 //this is the case when a new history is given to the
257 //navigator but old history has to be kept (e.g. LOOKZ
258 //is call during an event);
259 //touchable history has not been updated jet, so:
261 ptrTouchHist->UpdateYourself(pPhysVol,history);
262 ptrTempNavHist->UpdateYourself(pPhysVol,history);
268 //this is the case when the touchable history has been
269 //updated by a LocateGlobalPointAndUpdateTouchable call
271 G4VPhysicalVolume * pPhysVolTemp = ptrTempNavHist->GetVolume();
272 ptrOldNavHist->UpdateYourself(pPhysVolTemp,
273 ptrTempNavHist->GetHistory());
275 ptrTempNavHist->UpdateYourself(pPhysVol,history);
281 G4cout <<" ERROR in updating step-histories!" << G4endl;
286 #ifdef G4GEOMETRY_DEBUG
287 G4cout << "<== Flugg FGeometryInit::UpdateHistories()" << G4endl;
291 //*****************************************************************************
292 int FGeometryInit::GetLastMaterialIndex() const
294 // Get last material index as known by FLUKA
295 const FlukaMaterialsTable *matTable = FlukaMaterial::GetMaterialTable();
296 int matsize = matTable->size();
300 void FGeometryInit::createFlukaMatFile() {
301 // last modification Sara Vanini 1/III/99
302 // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case,
303 // according to the fluka standard. In addition,. they must be equal to the
304 // names of the fluka materials - see fluka manual - in order that the
305 // program load the right cross sections, and equal to the names included in
306 // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his
307 // own .pemf, in order to get the right cross sections loaded in memory.
309 #ifdef G4GEOMETRY_DEBUG
310 G4cout << "==> Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
311 G4cout << "================== FILEWR =================" << G4endl;
317 std::ofstream vos("Volumes_index.inp");
318 PrintRegionsMap(vos);
321 //Materials and compounds
322 BuildMaterialTables();
323 std::ofstream fos("flukaMat.inp");
324 PrintMaterialTables(fos);
326 PrintMagneticField(fos);
329 #ifdef G4GEOMETRY_DEBUG
330 G4cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
334 ////////////////////////////////////////////////////////////////////////
336 void FGeometryInit::BuildRegionsMap() {
337 #ifdef G4GEOMETRY_DEBUG
338 G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
341 //Find number of Volumes in physical volume store
342 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
343 unsigned int numVol = pVolStore->size();
345 G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore
346 << ") has " << numVol << " volumes. Iterating..."
349 for(unsigned int l=0; l < numVol; l++) {
350 //Get each of the physical volumes
351 G4VPhysicalVolume * physicalvolume = (*pVolStore)[l];
352 G4int iFlukaRegion = l+1;
353 fRegionVolumeMap[physicalvolume] = iFlukaRegion;
358 #ifdef G4GEOMETRY_DEBUG
359 G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
363 void FGeometryInit::PrintRegionsMap(std::ostream& os) {
364 #ifdef G4GEOMETRY_DEBUG
365 G4cout << "==> Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
369 PrintHeader(os, "GEANT4 VOLUMES");
371 //Iterate over all volumes in the map
372 for (RegionIterator i = fRegionVolumeMap.begin();
373 i != fRegionVolumeMap.end();
376 //Get info in the map
377 G4VPhysicalVolume* ptrVol = (*i).first;
378 int index = (*i).second;
380 //Print index and region name in some fixed format
381 os.setf(std::ios::left, std::ios::adjustfield);
382 os << setw10 << index;
383 os << std::setw(20) << ptrVol->GetName() << std::setw(20) << "";
385 //If volume is a replica... print some more stuff
386 if(ptrVol->IsReplicated()) {
390 G4double offset = -1;
391 G4bool consum = false;
392 ptrVol->GetReplicationData(axis, nRep, width, offset, consum);
393 os.setf(std::ios::left, std::ios::adjustfield);
394 os << setw10 << "Repetion Nb: " << std::setw(3) << nRep;
400 #ifdef G4GEOMETRY_DEBUG
401 G4cout << "<== Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
405 ////////////////////////////////////////////////////////////////////////
407 void FGeometryInit::BuildMediaMap()
409 fRegionMediumMap = new int[fNRegions+1];
410 for (RegionIterator i = fRegionVolumeMap.begin();
411 i != fRegionVolumeMap.end();
413 //Get info in the map
414 G4VPhysicalVolume* ptrVol = (*i).first;
415 int region = (*i).second;
416 G4int imed = fMediumVolumeMap[ptrVol];
417 fRegionMediumMap[region] = imed;
418 printf("BuildMediaMap %s %d %d\n",(ptrVol->GetName()).data(), region, imed);
423 G4int FGeometryInit::GetMedium(int region) const
425 return fRegionMediumMap[region];
429 void FGeometryInit::SetMediumFromName(const char* volName, int medium, int volid)
433 strncpy(tmp, volName, 4);
437 for (RegionIterator i = fRegionVolumeMap.begin();
438 i != fRegionVolumeMap.end();
441 //Get info in the map
442 G4VPhysicalVolume* ptrVol = (*i).first;
443 strncpy(name4, (ptrVol->GetName()).data(), 4);
445 for (int j = 0; j < 4; j++) {
446 if (name4[j] == '\0') {
447 for (int k = j; k < 4; k++) {
453 if (! strncmp(name4, tmp, 4)) {
454 fMediumVolumeMap[ptrVol] = medium;
455 fVolIdVolumeMap[ptrVol] = volid;
462 ////////////////////////////////////////////////////////////////////////
464 void FGeometryInit::BuildMaterialTables() {
465 #ifdef G4GEOMETRY_DEBUG
466 G4cout << "==> Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
469 //some terminal printout also
470 G4cout << "\t* Storing information..." << G4endl;
472 //The logic is the folloing:
473 //Get the Material Table and:
474 // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
475 // 2) For each single element material build a material equivalent
477 // 3.a) Build materials for each not already known element
478 // 3.b) Build the compound out of them
480 //Get the Material Table and iterate
481 const G4MaterialTable* matTable = G4Material::GetMaterialTable();
482 for (MatTableIterator i = matTable->begin(); i != matTable->end(); i++) {
484 //Get some basic material information
485 G4Material* material = (*i);
486 G4String matName = material->GetName();
487 const G4double matDensity = material->GetDensity();
488 const G4int nMatElements = material->GetNumberOfElements();
490 G4cout << "\t\t+ " << matName
491 << ": dens. = " << matDensity/(g/cm3) << "g/cm3"
492 << ", nElem = " << nMatElements << G4endl;
494 // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
495 // FlukaMaterial* is 0 in that case
496 if (matDensity <= 1.00e-10*g/cm3) {
497 G4FlukaMaterialMap[material] = 0;
498 G4cout << "\t\t Stored as vacuum" << G4endl;
500 // 2) For each single element material build a material equivalent
501 else if (nMatElements == 1) {
503 FlukaMaterial *flukamat =
504 BuildFlukaMaterialFromElement(material->GetElement(0),
507 G4FlukaMaterialMap[material] = flukamat;
508 G4cout << "\t\t Stored as " << flukamat->GetRealName() << G4endl;
510 } //else if (material->GetNumberOfElements() == 1)
513 // 3.a) Build materials for each not already known element
514 // 3.b) Build the compound out of them
516 FlukaCompound* flukacomp =
517 BuildFlukaCompoundFromMaterial(material);
518 G4FlukaCompoundMap[material] = flukacomp;
519 G4cout << "\t\t Stored as " << flukacomp->GetRealName() << G4endl;
523 #ifdef G4GEOMETRY_DEBUG
524 G4cout << "<== Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
529 FGeometryInit::BuildFlukaMaterialFromElement(const G4Element* element,
530 G4double matDensity) {
531 #ifdef G4GEOMETRY_DEBUG
532 G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromElement()"
536 //Get element and its properties
537 G4String elemName(ToFlukaString(element->GetName()));
539 FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(elemName);
540 if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
542 G4int nIsotopes = element->GetNumberOfIsotopes();
543 if (nIsotopes == 0) {
544 G4double elemA = element->GetA()/g;
545 G4double elemZ = element->GetZ();
547 if (elemA != G4int(elemA) && elemZ != G4int(elemZ)) {
548 G4cout << "WARNING: Element \'" << elemName
549 << "\' has non integer Z (" << elemZ << ") or A ("
554 flukamat = new FlukaMaterial(elemName,
559 else if (nIsotopes == 1) {
560 const G4Isotope* isotope = element->GetIsotope(0);
561 flukamat = BuildFlukaMaterialFromIsotope(isotope, matDensity);
564 FlukaCompound *flucomp = BuildFlukaCompoundFromElement(element,
566 flukamat = flucomp->GetFlukaMaterial();
569 #ifdef G4GEOMETRY_DEBUG
571 G4cout << "INFO: Element \'" << elemName
572 << "\' already exists in the DB. It will not be recreated."
579 #ifdef G4GEOMETRY_DEBUG
580 G4cout << "<== Flugg FGeometryInit::BuildFlukaMaterialFromElement()"
586 FGeometryInit::BuildFlukaMaterialFromIsotope(const G4Isotope* isotope,
587 G4double matDensity) {
588 #ifdef G4GEOMETRY_DEBUG
589 G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()"
592 G4String isoName(ToFlukaString(isotope->GetName()));
593 FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(isoName);
594 if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
595 G4int isoZ = isotope->GetZ();
596 G4double isoA = (isotope->GetA())/(g);
597 G4int isoN = isotope->GetN();
598 flukamat = new FlukaMaterial(isoName,
607 #ifdef G4GEOMETRY_DEBUG
608 G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()"
614 FGeometryInit::BuildFlukaCompoundFromMaterial(const G4Material* material) {
615 #ifdef G4GEOMETRY_DEBUG
616 G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()"
619 //Material properties
620 const G4double* elemFractions = material->GetFractionVector();
621 const G4int nMatElements = material->GetNumberOfElements();
622 const G4double matDensity = material->GetDensity();
623 G4String matName(ToFlukaString(material->GetName()));
624 FlukaCompound* flukacomp = new FlukaCompound(matName, matDensity/(g/cm3),
626 for (G4int i = 0; i < nMatElements; i++) {
627 FlukaMaterial *flukamat =
628 BuildFlukaMaterialFromElement(material->GetElement(i), 0.0);
630 flukacomp->AddElement(flukamat->GetIndex(), -elemFractions[i]);
636 #ifdef G4GEOMETRY_DEBUG
637 G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()"
643 FGeometryInit::BuildFlukaCompoundFromElement(const G4Element* element,
644 G4double matDensity) {
645 #ifdef G4GEOMETRY_DEBUG
646 G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromElement()"
649 G4int nIsotopes = element->GetNumberOfIsotopes();
650 //fraction of nb of atomes per volume (= volume fraction?)
651 const G4double* isoAbundance = element->GetRelativeAbundanceVector();
652 G4String elemName(ToFlukaString(element->GetName()));
654 //Material properties
655 FlukaCompound* flukacomp = new FlukaCompound(elemName, matDensity/(g/cm3),
657 for (G4int i = 0; i < nIsotopes; i++) {
658 FlukaMaterial *flukamat =
659 BuildFlukaMaterialFromIsotope(element->GetIsotope(i), 0.0);
661 flukacomp->AddElement(flukamat->GetIndex(), isoAbundance[i]);
667 #ifdef G4GEOMETRY_DEBUG
668 G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromElement()"
673 void FGeometryInit::PrintMaterialTables(std::ostream& os) {
674 #ifdef G4GEOMETRY_DEBUG
675 G4cout << "==> Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
678 PrintHeader(os, "GEANT4 MATERIALS AND COMPOUNDS");
680 //And some more stuff
681 size_t nIsotopes = G4Isotope::GetNumberOfIsotopes();
682 size_t nElements = G4Element::GetNumberOfElements();
683 size_t nMaterials = G4Material::GetNumberOfMaterials();
685 os << "* In Geant4 there are " << nMaterials << " materials" << G4endl;
686 os << "* In Geant4 there are " << nElements << " elements" << G4endl;
687 os << "* In Geant4 there are " << nIsotopes << " isotopes" << G4endl;
690 G4cout << "\t* Printing FLUKA materials..." << G4endl;
691 FlukaMaterial::PrintMaterialsByIndex(os);
692 //FlukaMaterial::PrintMaterialsByName(os);
695 G4cout << "\t* Printing FLUKA compounds..." << G4endl;
696 FlukaCompound::PrintCompounds(os);
698 #ifdef G4GEOMETRY_DEBUG
699 G4cout << "<== Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
703 ////////////////////////////////////////////////////////////////////////
705 void FGeometryInit::PrintAssignmat(std::ostream& os) {
706 #ifdef G4GEOMETRY_DEBUG
707 G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
710 //Find number of Volumes in physical volume store
711 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
712 unsigned int numVol = pVolStore->size();
714 G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore
715 << ") has " << numVol << " volumes. " << G4endl;
717 G4cout << "\t* Printing ASSIGNMAT..." << G4endl;
720 PrintHeader(os,"GEANT4 MATERIAL ASSIGNMENTS");
721 for(unsigned int l=0; l < numVol; l++) {
723 //Get each of the physical volumes
724 G4VPhysicalVolume * physicalvol = (*pVolStore)[l];
726 //Get index for that volume
727 G4int iFlukaRegion = fRegionVolumeMap[physicalvol];
729 //Find G4 material and navigate to its fluka compound/material
730 G4LogicalVolume * logicalVol = physicalvol->GetLogicalVolume();
731 G4Material* material = logicalVol->GetMaterial();
733 if (G4FlukaCompoundMap[material])
734 matIndex = G4FlukaCompoundMap[material]->GetIndex();
735 if (G4FlukaMaterialMap[material])
736 matIndex = G4FlukaMaterialMap[material]->GetIndex();
738 //Find if there is a magnetic field in the region
739 //check if Magnetic Field is present in the region
740 G4double flagField = 0.0;
741 G4FieldManager * pMagFieldMan = logicalVol->GetFieldManager();
742 if(pMagFieldMan && pMagFieldMan->GetDetectorField())
746 os << setw10 << "ASSIGNMAT ";
747 os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
748 os << setw10 << setfixed << G4double(matIndex);
749 os << setw10 << setfixed << G4double(iFlukaRegion);
750 os << setw10 << "0.0";
751 os << setw10 << setfixed << flagField;
757 #ifdef G4GEOMETRY_DEBUG
758 G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
763 void FGeometryInit::PrintMagneticField(std::ostream& os) {
764 #ifdef G4GEOMETRY_DEBUG
765 G4cout << "==> Flugg FGeometryInit::PrintMagneticField()" << G4endl;
768 G4cout << "\t* Printing Magnetic Field..." << G4endl;
770 if(fTransportationManager->GetFieldManager()->DoesFieldExist()) {
772 //get magnetic field pointer
773 const G4Field * pMagField =
774 fTransportationManager->GetFieldManager()->GetDetectorField();
778 //Check if it can be made a uniform magnetic field
779 const G4UniformMagField *pUnifMagField =
780 dynamic_cast<const G4UniformMagField*>(pMagField);
783 G4double point[4]; //it is not really used
784 pUnifMagField->GetFieldValue(point,B);
786 //write MGNFIELD card
787 PrintHeader(os,"GEANT4 MAGNETIC FIELD");
788 os << setw10 << "MGNFIELD ";
792 os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
793 os << setw10 << setfixed
794 << std::setprecision(4) << B[0]
800 G4cout << "WARNING: No Uniform Magnetic Field found." << G4endl;
801 G4cout << " Manual intervention might be needed." << G4endl;
805 G4cout << "\t No detector field found... " << G4endl;
806 } // end if magnetic field
808 G4cout << "\t No field found... " << G4endl;
810 #ifdef G4GEOMETRY_DEBUG
811 G4cout << "<== Flugg FGeometryInit::PrintMagneticField()" << G4endl;
815 int FGeometryInit::CurrentVolID(int ir, int& copyNo)
823 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
824 G4VPhysicalVolume * physicalvol = (*pVolStore)[ir- 1];
827 copyNo = physicalvol->GetCopyNo();
834 int id = fVolIdVolumeMap[physicalvol];
838 int FGeometryInit::CurrentVolOffID(int ir, int off, int& copyNo)
840 if (off == 0) return CurrentVolID(ir, copyNo);
842 G4PhysicalVolumeStore* pVolStore = G4PhysicalVolumeStore::GetInstance();
843 G4VPhysicalVolume* physicalvol = (*pVolStore)[ir- 1];
844 G4VPhysicalVolume* mother = physicalvol;
847 //============================================================================
849 // Check touchable depth
851 if (ptrTouchHist->GetHistoryDepth() < off) {
854 // Get the off-th mother
855 index = ptrTouchHist->GetHistoryDepth() - off;
856 // in the touchable history volumes are ordered
857 // from top volume up to mother volume;
858 // the touchable volume is not in the history
859 mother = ptrTouchHist->GetHistory()->GetVolume(index);
862 //============================================================================
867 G4cout << "Flugg FGeometryInit::CurrentVolOffID mother not found" << G4endl;
871 copyNo = mother ->GetCopyNo();
872 id = fVolIdVolumeMap[mother];
877 void FGeometryInit::Gmtod(double* xm, double* xd, int iflag)
879 // Transforms a position from the world reference frame
880 // to the current volume reference frame.
882 // Geant3 desription:
883 // ==================
884 // Computes coordinates XD (in DRS)
885 // from known coordinates XM in MRS
886 // The local reference system can be initialized by
887 // - the tracking routines and GMTOD used in GUSTEP
888 // - a call to GMEDIA(XM,NUMED)
889 // - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
890 // (inverse routine is GDTOM)
892 // If IFLAG=1 convert coordinates
893 // IFLAG=2 convert direction cosinus
896 FluggNavigator * ptrNavig = getNavigatorForTracking();
897 //setting variables (and dimension: Fluka uses cm.!)
898 G4ThreeVector pGlob(xm[0],xm[1],xm[2]);
902 pGlob *= 10.0; // in mm
903 // change because of geant 4 6.0
904 // pLoc = ptrNavig->ComputeLocalPoint(pGlob);
905 pLoc = ptrNavig->GetGlobalToLocalTransform().TransformPoint(pGlob);
907 pLoc /= 10.0; // in cm
908 } else if (iflag == 2) {
910 ptrNavig->ComputeLocalAxis(pGlob);
912 G4cout << "Flugg FGeometryInit::Gmtod called with undefined flag" << G4endl;
914 xd[0] = pLoc[0]; xd[1] = pLoc[1]; xd[2] = pLoc[2];
917 void FGeometryInit::Gdtom(double* xd, double* xm, int iflag)
919 // Transforms a position from the current volume reference frame
920 // to the world reference frame.
922 // Geant3 desription:
923 // ==================
924 // Computes coordinates XM (Master Reference System
925 // knowing the coordinates XD (Detector Ref System)
926 // The local reference system can be initialized by
927 // - the tracking routines and GDTOM used in GUSTEP
928 // - a call to GSCMED(NLEVEL,NAMES,NUMBER)
929 // (inverse routine is GMTOD)
931 // If IFLAG=1 convert coordinates
932 // IFLAG=2 convert direction cosinus
936 FluggNavigator * ptrNavig = getNavigatorForTracking();
937 G4ThreeVector pLoc(xd[0],xd[1],xd[2]);
941 pLoc *= 10.0; // in mm
942 pGlob = ptrNavig->GetLocalToGlobalTransform().
943 TransformPoint(pLoc);
944 pGlob /= 10.0; // in cm
945 } else if (iflag == 2) {
946 pGlob = ptrNavig->GetLocalToGlobalTransform().
949 G4cout << "Flugg FGeometryInit::Gdtom called with undefined flag" << G4endl;
952 xm[0] = pGlob[0]; xm[1] = pGlob[1]; xm[2] = pGlob[2];