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 //*****************************************************************************
293 void FGeometryInit::createFlukaMatFile() {
294 // last modification Sara Vanini 1/III/99
295 // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case,
296 // according to the fluka standard. In addition,. they must be equal to the
297 // names of the fluka materials - see fluka manual - in order that the
298 // program load the right cross sections, and equal to the names included in
299 // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his
300 // own .pemf, in order to get the right cross sections loaded in memory.
302 #ifdef G4GEOMETRY_DEBUG
303 G4cout << "==> Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
304 G4cout << "================== FILEWR =================" << G4endl;
310 std::ofstream vos("Volumes_index.inp");
311 PrintRegionsMap(vos);
314 //Materials and compounds
315 BuildMaterialTables();
316 std::ofstream fos("flukaMat.inp");
317 PrintMaterialTables(fos);
319 PrintMagneticField(fos);
322 #ifdef G4GEOMETRY_DEBUG
323 G4cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
327 ////////////////////////////////////////////////////////////////////////
329 void FGeometryInit::BuildRegionsMap() {
330 #ifdef G4GEOMETRY_DEBUG
331 G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
334 //Find number of Volumes in physical volume store
335 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
336 unsigned int numVol = pVolStore->size();
338 G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore
339 << ") has " << numVol << " volumes. Iterating..."
342 for(unsigned int l=0; l < numVol; l++) {
343 //Get each of the physical volumes
344 G4VPhysicalVolume * physicalvolume = (*pVolStore)[l];
345 G4int iFlukaRegion = l+1;
346 fRegionVolumeMap[physicalvolume] = iFlukaRegion;
351 #ifdef G4GEOMETRY_DEBUG
352 G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
356 void FGeometryInit::PrintRegionsMap(std::ostream& os) {
357 #ifdef G4GEOMETRY_DEBUG
358 G4cout << "==> Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
362 PrintHeader(os, "GEANT4 VOLUMES");
364 //Iterate over all volumes in the map
365 for (RegionIterator i = fRegionVolumeMap.begin();
366 i != fRegionVolumeMap.end();
369 //Get info in the map
370 G4VPhysicalVolume* ptrVol = (*i).first;
371 int index = (*i).second;
373 //Print index and region name in some fixed format
374 os.setf(std::ios::left, std::ios::adjustfield);
375 os << setw10 << index;
376 os << std::setw(20) << ptrVol->GetName() << std::setw(20) << "";
378 //If volume is a replica... print some more stuff
379 if(ptrVol->IsReplicated()) {
383 G4double offset = -1;
384 G4bool consum = false;
385 ptrVol->GetReplicationData(axis, nRep, width, offset, consum);
386 os.setf(std::ios::left, std::ios::adjustfield);
387 os << setw10 << "Repetion Nb: " << std::setw(3) << nRep;
393 #ifdef G4GEOMETRY_DEBUG
394 G4cout << "<== Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
398 ////////////////////////////////////////////////////////////////////////
400 void FGeometryInit::BuildMediaMap()
402 fRegionMediumMap = new int[fNRegions+1];
403 for (RegionIterator i = fRegionVolumeMap.begin();
404 i != fRegionVolumeMap.end();
406 //Get info in the map
407 G4VPhysicalVolume* ptrVol = (*i).first;
408 int region = (*i).second;
409 G4int imed = fMediumVolumeMap[ptrVol];
410 fRegionMediumMap[region] = imed;
411 printf("BuildMediaMap %s %d %d\n",(ptrVol->GetName()).data(), region, imed);
416 G4int FGeometryInit::GetMedium(int region) const
418 return fRegionMediumMap[region];
422 void FGeometryInit::SetMediumFromName(const char* volName, int medium, int volid)
426 strncpy(tmp, volName, 4);
430 for (RegionIterator i = fRegionVolumeMap.begin();
431 i != fRegionVolumeMap.end();
434 //Get info in the map
435 G4VPhysicalVolume* ptrVol = (*i).first;
436 strncpy(name4, (ptrVol->GetName()).data(), 4);
438 for (int j = 0; j < 4; j++) {
439 if (name4[j] == '\0') {
440 for (int k = j; k < 4; k++) {
446 if (! strncmp(name4, tmp, 4)) {
447 fMediumVolumeMap[ptrVol] = medium;
448 fVolIdVolumeMap[ptrVol] = volid;
455 ////////////////////////////////////////////////////////////////////////
457 void FGeometryInit::BuildMaterialTables() {
458 #ifdef G4GEOMETRY_DEBUG
459 G4cout << "==> Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
462 //some terminal printout also
463 G4cout << "\t* Storing information..." << G4endl;
465 //The logic is the folloing:
466 //Get the Material Table and:
467 // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
468 // 2) For each single element material build a material equivalent
470 // 3.a) Build materials for each not already known element
471 // 3.b) Build the compound out of them
473 //Get the Material Table and iterate
474 const G4MaterialTable* matTable = G4Material::GetMaterialTable();
475 for (MatTableIterator i = matTable->begin(); i != matTable->end(); i++) {
477 //Get some basic material information
478 G4Material* material = (*i);
479 G4String matName = material->GetName();
480 const G4double matDensity = material->GetDensity();
481 const G4int nMatElements = material->GetNumberOfElements();
483 G4cout << "\t\t+ " << matName
484 << ": dens. = " << matDensity/(g/cm3) << "g/cm3"
485 << ", nElem = " << nMatElements << G4endl;
487 // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
488 // FlukaMaterial* is 0 in that case
489 if (matDensity <= 1.00e-10*g/cm3) {
490 G4FlukaMaterialMap[material] = 0;
491 G4cout << "\t\t Stored as vacuum" << G4endl;
493 // 2) For each single element material build a material equivalent
494 else if (nMatElements == 1) {
496 FlukaMaterial *flukamat =
497 BuildFlukaMaterialFromElement(material->GetElement(0),
500 G4FlukaMaterialMap[material] = flukamat;
501 G4cout << "\t\t Stored as " << flukamat->GetRealName() << G4endl;
503 } //else if (material->GetNumberOfElements() == 1)
506 // 3.a) Build materials for each not already known element
507 // 3.b) Build the compound out of them
509 FlukaCompound* flukacomp =
510 BuildFlukaCompoundFromMaterial(material);
511 G4FlukaCompoundMap[material] = flukacomp;
512 G4cout << "\t\t Stored as " << flukacomp->GetRealName() << G4endl;
516 #ifdef G4GEOMETRY_DEBUG
517 G4cout << "<== Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
522 FGeometryInit::BuildFlukaMaterialFromElement(const G4Element* element,
523 G4double matDensity) {
524 #ifdef G4GEOMETRY_DEBUG
525 G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromElement()"
529 //Get element and its properties
530 G4String elemName(ToFlukaString(element->GetName()));
532 FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(elemName);
533 if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
535 G4int nIsotopes = element->GetNumberOfIsotopes();
536 if (nIsotopes == 0) {
537 G4double elemA = element->GetA()/g;
538 G4double elemZ = element->GetZ();
540 if (elemA != G4int(elemA) && elemZ != G4int(elemZ)) {
541 G4cout << "WARNING: Element \'" << elemName
542 << "\' has non integer Z (" << elemZ << ") or A ("
547 flukamat = new FlukaMaterial(elemName,
552 else if (nIsotopes == 1) {
553 const G4Isotope* isotope = element->GetIsotope(0);
554 flukamat = BuildFlukaMaterialFromIsotope(isotope, matDensity);
557 FlukaCompound *flucomp = BuildFlukaCompoundFromElement(element,
559 flukamat = flucomp->GetFlukaMaterial();
562 #ifdef G4GEOMETRY_DEBUG
564 G4cout << "INFO: Element \'" << elemName
565 << "\' already exists in the DB. It will not be recreated."
572 #ifdef G4GEOMETRY_DEBUG
573 G4cout << "<== Flugg FGeometryInit::BuildFlukaMaterialFromElement()"
579 FGeometryInit::BuildFlukaMaterialFromIsotope(const G4Isotope* isotope,
580 G4double matDensity) {
581 #ifdef G4GEOMETRY_DEBUG
582 G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()"
585 G4String isoName(ToFlukaString(isotope->GetName()));
586 FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(isoName);
587 if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
588 G4int isoZ = isotope->GetZ();
589 G4double isoA = (isotope->GetA())/(g);
590 G4int isoN = isotope->GetN();
591 flukamat = new FlukaMaterial(isoName,
600 #ifdef G4GEOMETRY_DEBUG
601 G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()"
607 FGeometryInit::BuildFlukaCompoundFromMaterial(const G4Material* material) {
608 #ifdef G4GEOMETRY_DEBUG
609 G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()"
612 //Material properties
613 const G4double* elemFractions = material->GetFractionVector();
614 const G4int nMatElements = material->GetNumberOfElements();
615 const G4double matDensity = material->GetDensity();
616 G4String matName(ToFlukaString(material->GetName()));
617 FlukaCompound* flukacomp = new FlukaCompound(matName, matDensity/(g/cm3),
619 for (G4int i = 0; i < nMatElements; i++) {
620 FlukaMaterial *flukamat =
621 BuildFlukaMaterialFromElement(material->GetElement(i), 0.0);
623 flukacomp->AddElement(flukamat->GetIndex(), -elemFractions[i]);
629 #ifdef G4GEOMETRY_DEBUG
630 G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()"
636 FGeometryInit::BuildFlukaCompoundFromElement(const G4Element* element,
637 G4double matDensity) {
638 #ifdef G4GEOMETRY_DEBUG
639 G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromElement()"
642 G4int nIsotopes = element->GetNumberOfIsotopes();
643 //fraction of nb of atomes per volume (= volume fraction?)
644 const G4double* isoAbundance = element->GetRelativeAbundanceVector();
645 G4String elemName(ToFlukaString(element->GetName()));
647 //Material properties
648 FlukaCompound* flukacomp = new FlukaCompound(elemName, matDensity/(g/cm3),
650 for (G4int i = 0; i < nIsotopes; i++) {
651 FlukaMaterial *flukamat =
652 BuildFlukaMaterialFromIsotope(element->GetIsotope(i), 0.0);
654 flukacomp->AddElement(flukamat->GetIndex(), isoAbundance[i]);
660 #ifdef G4GEOMETRY_DEBUG
661 G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromElement()"
666 void FGeometryInit::PrintMaterialTables(std::ostream& os) {
667 #ifdef G4GEOMETRY_DEBUG
668 G4cout << "==> Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
671 PrintHeader(os, "GEANT4 MATERIALS AND COMPOUNDS");
673 //And some more stuff
674 size_t nIsotopes = G4Isotope::GetNumberOfIsotopes();
675 size_t nElements = G4Element::GetNumberOfElements();
676 size_t nMaterials = G4Material::GetNumberOfMaterials();
678 os << "* In Geant4 there are " << nMaterials << " materials" << G4endl;
679 os << "* In Geant4 there are " << nElements << " elements" << G4endl;
680 os << "* In Geant4 there are " << nIsotopes << " isotopes" << G4endl;
683 G4cout << "\t* Printing FLUKA materials..." << G4endl;
684 FlukaMaterial::PrintMaterialsByIndex(os);
685 //FlukaMaterial::PrintMaterialsByName(os);
688 G4cout << "\t* Printing FLUKA compounds..." << G4endl;
689 FlukaCompound::PrintCompounds(os);
691 #ifdef G4GEOMETRY_DEBUG
692 G4cout << "<== Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
696 ////////////////////////////////////////////////////////////////////////
698 void FGeometryInit::PrintAssignmat(std::ostream& os) {
699 #ifdef G4GEOMETRY_DEBUG
700 G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
703 //Find number of Volumes in physical volume store
704 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
705 unsigned int numVol = pVolStore->size();
707 G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore
708 << ") has " << numVol << " volumes. " << G4endl;
710 G4cout << "\t* Printing ASSIGNMAT..." << G4endl;
713 PrintHeader(os,"GEANT4 MATERIAL ASSIGNMENTS");
714 for(unsigned int l=0; l < numVol; l++) {
716 //Get each of the physical volumes
717 G4VPhysicalVolume * physicalvol = (*pVolStore)[l];
719 //Get index for that volume
720 G4int iFlukaRegion = fRegionVolumeMap[physicalvol];
722 //Find G4 material and navigate to its fluka compound/material
723 G4LogicalVolume * logicalVol = physicalvol->GetLogicalVolume();
724 G4Material* material = logicalVol->GetMaterial();
726 if (G4FlukaCompoundMap[material])
727 matIndex = G4FlukaCompoundMap[material]->GetIndex();
728 if (G4FlukaMaterialMap[material])
729 matIndex = G4FlukaMaterialMap[material]->GetIndex();
731 //Find if there is a magnetic field in the region
732 //check if Magnetic Field is present in the region
733 G4double flagField = 0.0;
734 G4FieldManager * pMagFieldMan = logicalVol->GetFieldManager();
735 if(pMagFieldMan && pMagFieldMan->GetDetectorField())
739 os << setw10 << "ASSIGNMAT ";
740 os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
741 os << setw10 << setfixed << G4double(matIndex);
742 os << setw10 << setfixed << G4double(iFlukaRegion);
743 os << setw10 << "0.0";
744 os << setw10 << setfixed << flagField;
750 #ifdef G4GEOMETRY_DEBUG
751 G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
756 void FGeometryInit::PrintMagneticField(std::ostream& os) {
757 #ifdef G4GEOMETRY_DEBUG
758 G4cout << "==> Flugg FGeometryInit::PrintMagneticField()" << G4endl;
761 G4cout << "\t* Printing Magnetic Field..." << G4endl;
763 if(fTransportationManager->GetFieldManager()->DoesFieldExist()) {
765 //get magnetic field pointer
766 const G4Field * pMagField =
767 fTransportationManager->GetFieldManager()->GetDetectorField();
771 //Check if it can be made a uniform magnetic field
772 const G4UniformMagField *pUnifMagField =
773 dynamic_cast<const G4UniformMagField*>(pMagField);
776 G4double point[4]; //it is not really used
777 pUnifMagField->GetFieldValue(point,B);
779 //write MGNFIELD card
780 PrintHeader(os,"GEANT4 MAGNETIC FIELD");
781 os << setw10 << "MGNFIELD ";
785 os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
786 os << setw10 << setfixed
787 << std::setprecision(4) << B[0]
793 G4cout << "WARNING: No Uniform Magnetic Field found." << G4endl;
794 G4cout << " Manual intervention might be needed." << G4endl;
798 G4cout << "\t No detector field found... " << G4endl;
799 } // end if magnetic field
801 G4cout << "\t No field found... " << G4endl;
803 #ifdef G4GEOMETRY_DEBUG
804 G4cout << "<== Flugg FGeometryInit::PrintMagneticField()" << G4endl;
808 int FGeometryInit::CurrentVolID(int ir, int& copyNo)
816 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
817 G4VPhysicalVolume * physicalvol = (*pVolStore)[ir- 1];
820 copyNo = physicalvol->GetCopyNo();
827 int id = fVolIdVolumeMap[physicalvol];
831 int FGeometryInit::CurrentVolOffID(int ir, int off, int& copyNo)
833 if (off == 0) return CurrentVolID(ir, copyNo);
835 G4PhysicalVolumeStore* pVolStore = G4PhysicalVolumeStore::GetInstance();
836 G4VPhysicalVolume* physicalvol = (*pVolStore)[ir- 1];
837 G4VPhysicalVolume* mother = physicalvol;
840 //============================================================================
842 // Check touchable depth
844 if (ptrTouchHist->GetHistoryDepth() < off) {
847 // Get the off-th mother
848 index = ptrTouchHist->GetHistoryDepth() - off;
849 // in the touchable history volumes are ordered
850 // from top volume up to mother volume;
851 // the touchable volume is not in the history
852 mother = ptrTouchHist->GetHistory()->GetVolume(index);
855 //============================================================================
860 G4cout << "Flugg FGeometryInit::CurrentVolOffID mother not found" << G4endl;
864 copyNo = mother ->GetCopyNo();
865 id = fVolIdVolumeMap[mother];
870 void FGeometryInit::Gmtod(double* xm, double* xd, int iflag)
872 // Transforms a position from the world reference frame
873 // to the current volume reference frame.
875 // Geant3 desription:
876 // ==================
877 // Computes coordinates XD (in DRS)
878 // from known coordinates XM in MRS
879 // The local reference system can be initialized by
880 // - the tracking routines and GMTOD used in GUSTEP
881 // - a call to GMEDIA(XM,NUMED)
882 // - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
883 // (inverse routine is GDTOM)
885 // If IFLAG=1 convert coordinates
886 // IFLAG=2 convert direction cosinus
889 FluggNavigator * ptrNavig = getNavigatorForTracking();
890 //setting variables (and dimension: Fluka uses cm.!)
891 G4ThreeVector pGlob(xm[0],xm[1],xm[2]);
895 pGlob *= 10.0; // in mm
896 // change because of geant 4 6.0
897 // pLoc = ptrNavig->ComputeLocalPoint(pGlob);
898 pLoc = ptrNavig->GetGlobalToLocalTransform().TransformPoint(pGlob);
900 pLoc /= 10.0; // in cm
901 } else if (iflag == 2) {
903 ptrNavig->ComputeLocalAxis(pGlob);
905 G4cout << "Flugg FGeometryInit::Gmtod called with undefined flag" << G4endl;
907 xd[0] = pLoc[0]; xd[1] = pLoc[1]; xd[2] = pLoc[2];
910 void FGeometryInit::Gdtom(double* xd, double* xm, int iflag)
912 // Transforms a position from the current volume reference frame
913 // to the world reference frame.
915 // Geant3 desription:
916 // ==================
917 // Computes coordinates XM (Master Reference System
918 // knowing the coordinates XD (Detector Ref System)
919 // The local reference system can be initialized by
920 // - the tracking routines and GDTOM used in GUSTEP
921 // - a call to GSCMED(NLEVEL,NAMES,NUMBER)
922 // (inverse routine is GMTOD)
924 // If IFLAG=1 convert coordinates
925 // IFLAG=2 convert direction cosinus
929 FluggNavigator * ptrNavig = getNavigatorForTracking();
930 G4ThreeVector pLoc(xd[0],xd[1],xd[2]);
934 pLoc *= 10.0; // in mm
935 pGlob = ptrNavig->GetLocalToGlobalTransform().
936 TransformPoint(pLoc);
937 pGlob /= 10.0; // in cm
938 } else if (iflag == 2) {
939 pGlob = ptrNavig->GetLocalToGlobalTransform().
942 G4cout << "Flugg FGeometryInit::Gdtom called with undefined flag" << G4endl;
945 xm[0] = pGlob[0]; xm[1] = pGlob[1]; xm[2] = pGlob[2];