+++ /dev/null
-
-// Flugg tag
-// modified 17/06/02: by I. Gonzalez. STL migration
-
-//#include <stdio.h>
-//#include <iomanip.h>
-#include "FGeometryInit.hh"
-#include "FluggNavigator.hh"
-#include "WrapUtils.hh"
-#include "FlukaMaterial.hh"
-#include "FlukaCompound.hh"
-
-FGeometryInit * FGeometryInit::flagInstance=0;
-
-FGeometryInit* FGeometryInit::GetInstance() {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg::FGeometryInit::GetInstance(), instance="
- << flagInstance << G4endl;
-#endif
- if (!flagInstance)
- flagInstance = new FGeometryInit();
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg::FGeometryInit::GetInstance(), instance="
- << flagInstance << G4endl;
-#endif
- return flagInstance;
-}
-
-
-FGeometryInit::FGeometryInit():
- fDetector(0),
- fFieldManager(0),
- fTransportationManager(G4TransportationManager::GetTransportationManager()),
- myTopNode(0),
- ptrGeoMan(0),
- ptrArray(0),
- ptrTouchHist(0),
- ptrOldNavHist(0),
- ptrTempNavHist(0),
- ptrJrLtGeant(0)
-{
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::FGeometryInit()" << G4endl;
- G4cout << "\t+ Changing the G4Navigator for FluggNavigator..." << G4endl;
-#endif
- G4Navigator* actualnav = fTransportationManager->GetNavigatorForTracking();
- if (actualnav) {
- FluggNavigator* newnav = new FluggNavigator();
- fTransportationManager->SetNavigatorForTracking(newnav);
- }
- else {
- G4cerr << "ERROR: Could not find the actual G4Navigator" << G4endl;
- abort();
- }
-
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
-#endif
-
-}
-
-
-FGeometryInit::~FGeometryInit() {
- G4cout << "==> Flugg FGeometryInit::~FGeometryInit()" << G4endl;
- DeleteHistories();
- ptrGeoMan->OpenGeometry();
- if (fTransportationManager)
- delete fTransportationManager;
- if (ptrJrLtGeant)
- delete[] ptrJrLtGeant;
- DelHistArray();
-
- //keep ATTENTION: never delete a pointer twice!
- G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
-}
-
-void FGeometryInit::Init() {
-// Build and initialize G4 geometry
- setDetector();
- setMotherVolume();
- closeGeometry();
- InitHistories();
- InitJrLtGeantArray();
- InitHistArray();
- createFlukaMatFile();
-}
-
-
-void FGeometryInit::closeGeometry() {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::closeGeometry()" << G4endl;
-#endif
-
- ptrGeoMan = G4GeometryManager::GetInstance();
- if (ptrGeoMan) {
- ptrGeoMan->OpenGeometry();
-
- //true argoment allows voxel construction; if false voxels are built
- //only for replicated volumes
- ptrGeoMan->CloseGeometry(true);
- }
- else {
- G4cerr << "ERROR in FLUGG: Could not get G4GeometryManager instance"
- << G4endl;
- G4cerr << " in FGeometryInit::closeGeometry. Exiting!!!"
- << G4endl;
- }
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FGeometryInit::closeGeometry()" << G4endl;
-#endif
-}
-
-//*************************************************************************
-
-void FGeometryInit::InitHistArray() {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::InitHistArray()" << G4endl;
-#endif
- ptrArray = new G4int[1000000];
- for(G4int i=0;i<1000000;i++)
- ptrArray[i]=0;
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FGeometryInit::InitHistArray()" << G4endl;
-#endif
-}
-
-
-
-//*************************************************************************
-//jrLtGeant stores all crossed lattice volume histories.
-
-void FGeometryInit::InitJrLtGeantArray() {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
- G4cout << "Initializing JrLtGeant array" << G4endl;
-#endif
- ptrJrLtGeant = new G4int[10000];
- for(G4int x=0;x<10000;x++)
- ptrJrLtGeant[x]=-1;
- flagLttcGeant = -1;
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
-#endif
-}
-
-
-void FGeometryInit::SetLttcFlagGeant(G4int newFlagLttc) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
-#endif
- // Added by A.Solodkov
- if (newFlagLttc >= 10000) {
- G4cout << "Problems in FGeometryInit::SetLttcFlagGeant" << G4endl;
- G4cout << "Index newFlagLttc=" << newFlagLttc << " is outside array bounds"
- << G4endl;
- G4cout << "Better to stop immediately !" << G4endl;
- exit(1);
- }
- flagLttcGeant = newFlagLttc;
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
-#endif
-}
-
-void FGeometryInit::PrintJrLtGeant() {
-#ifdef G4GEOMETRY_DEBUG
- //G4cout << "jrLtGeant:" << G4endl;
- //for(G4int y=0;y<=flagLttcGeant;y++)
- //
- // G4cout << "jrLtGeant[" << y << "]=" << ptrJrLtGeant[y] << G4endl;
-#endif
-}
-
-//**************************************************************************
-
-void FGeometryInit::PrintHistories() {
- /*
- #ifdef G4GEOMETRY_DEBUG
- G4cout << "Touch hist:" << G4endl;
- G4cout << *(ptrTouchHist->GetHistory()) << G4endl;
- G4cout << "Tmp hist:" << G4endl;
- G4cout << *(ptrTempNavHist->GetHistory()) << G4endl;
- G4cout << "Old hist:" << G4endl;
- G4cout << *(ptrOldNavHist->GetHistory()) << G4endl;
- #endif
- */
-}
-
-
-
-void FGeometryInit::InitHistories() {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::InitHistories()" << G4endl;
-#endif
- //init utility histories with navigator history
- ptrTouchHist =
- fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
- ptrTempNavHist =
- fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
- ptrOldNavHist = new G4TouchableHistory();
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FGeometryInit::InitHistories()" << G4endl;
-#endif
-}
-
-void FGeometryInit::DeleteHistories() {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::DeleteHistories()" << G4endl;
-#endif
-
- delete ptrTouchHist;
- delete ptrOldNavHist;
- delete ptrTempNavHist;
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "Deleting step-history objects at end of run!" << G4endl;
- G4cout << "<== Flugg FGeometryInit::DeleteHistories()" << G4endl;
-#endif
-}
-
-
-void FGeometryInit::UpdateHistories(const G4NavigationHistory * history,
- G4int flagHist) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::UpdateHistories()" << G4endl;
-#endif
- PrintHistories();
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "...updating histories!" << G4endl;
-#endif
-
- G4VPhysicalVolume * pPhysVol = history->GetTopVolume();
-
- switch (flagHist) {
- case 0: {
- //this is the case when a new history is given to the
- //navigator and old history has to be resetted
- //touchable history has not been updated jet, so:
-
- ptrTouchHist->UpdateYourself(pPhysVol,history);
- ptrTempNavHist->UpdateYourself(pPhysVol,history);
- G4NavigationHistory * ptrOldNavHistNotConst =
- const_cast<G4NavigationHistory * >(ptrOldNavHist->GetHistory());
- ptrOldNavHistNotConst->Reset();
- ptrOldNavHistNotConst->Clear();
- PrintHistories();
- break;
- } //case 0
-
- case 1: {
- //this is the case when a new history is given to the
- //navigator but old history has to be kept (e.g. LOOKZ
- //is call during an event);
- //touchable history has not been updated jet, so:
-
- ptrTouchHist->UpdateYourself(pPhysVol,history);
- ptrTempNavHist->UpdateYourself(pPhysVol,history);
- PrintHistories();
- break;
- } //case 1
-
- case 2: {
- //this is the case when the touchable history has been
- //updated by a LocateGlobalPointAndUpdateTouchable call
-
- G4VPhysicalVolume * pPhysVolTemp = ptrTempNavHist->GetVolume();
- ptrOldNavHist->UpdateYourself(pPhysVolTemp,
- ptrTempNavHist->GetHistory());
-
- ptrTempNavHist->UpdateYourself(pPhysVol,history);
- PrintHistories();
- break;
- } //case 2
-
- default: {
- G4cout <<" ERROR in updating step-histories!" << G4endl;
- break;
- } //default
- } //switch
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FGeometryInit::UpdateHistories()" << G4endl;
-#endif
-}
-
-//*****************************************************************************
-int FGeometryInit::GetLastMaterialIndex() const
-{
-// Get last material index as known by FLUKA
- const FlukaMaterialsTable *matTable = FlukaMaterial::GetMaterialTable();
- int matsize = matTable->size();
- return matsize+2;
-}
-
-void FGeometryInit::createFlukaMatFile() {
- // last modification Sara Vanini 1/III/99
- // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case,
- // according to the fluka standard. In addition,. they must be equal to the
- // names of the fluka materials - see fluka manual - in order that the
- // program load the right cross sections, and equal to the names included in
- // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his
- // own .pemf, in order to get the right cross sections loaded in memory.
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
- G4cout << "================== FILEWR =================" << G4endl;
-#endif
-
-
- //Regions map
- BuildRegionsMap();
- std::ofstream vos("Volumes_index.inp");
- PrintRegionsMap(vos);
- vos.close();
-
- //Materials and compounds
- BuildMaterialTables();
- std::ofstream fos("flukaMat.inp");
- PrintMaterialTables(fos);
- PrintAssignmat(fos);
- PrintMagneticField(fos);
- fos.close();
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
-#endif
-}
-
-////////////////////////////////////////////////////////////////////////
-//
-void FGeometryInit::BuildRegionsMap() {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
-#endif
-
- //Find number of Volumes in physical volume store
- G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
- unsigned int numVol = pVolStore->size();
-
- G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore
- << ") has " << numVol << " volumes. Iterating..."
- << G4endl;
-
- for(unsigned int l=0; l < numVol; l++) {
- //Get each of the physical volumes
- G4VPhysicalVolume * physicalvolume = (*pVolStore)[l];
- G4int iFlukaRegion = l+1;
- fRegionVolumeMap[physicalvolume] = iFlukaRegion;
- }
-
-
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
-#endif
-}
-
-void FGeometryInit::PrintRegionsMap(std::ostream& os) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
-#endif
-
- //Print some header
- PrintHeader(os, "GEANT4 VOLUMES");
-
- //Iterate over all volumes in the map
- for (RegionIterator i = fRegionVolumeMap.begin();
- i != fRegionVolumeMap.end();
- i++) {
-
- //Get info in the map
- G4VPhysicalVolume* ptrVol = (*i).first;
- int index = (*i).second;
-
- //Print index and region name in some fixed format
- os.setf(std::ios::left, std::ios::adjustfield);
- os << setw10 << index;
- os << std::setw(20) << ptrVol->GetName() << std::setw(20) << "";
-
- //If volume is a replica... print some more stuff
- if(ptrVol->IsReplicated()) {
- EAxis axis;
- G4int nRep = -1;
- G4double width = -1;
- G4double offset = -1;
- G4bool consum = false;
- ptrVol->GetReplicationData(axis, nRep, width, offset, consum);
- os.setf(std::ios::left, std::ios::adjustfield);
- os << setw10 << "Repetion Nb: " << std::setw(3) << nRep;
- }
- os << G4endl;
-
- }
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
-#endif
-}
-
-////////////////////////////////////////////////////////////////////////
-//
-void FGeometryInit::BuildMediaMap()
-{
- fRegionMediumMap = new int[fNRegions+1];
- for (RegionIterator i = fRegionVolumeMap.begin();
- i != fRegionVolumeMap.end();
- i++) {
- //Get info in the map
- G4VPhysicalVolume* ptrVol = (*i).first;
- int region = (*i).second;
- G4int imed = fMediumVolumeMap[ptrVol];
- fRegionMediumMap[region] = imed;
- printf("BuildMediaMap %s %d %d\n",(ptrVol->GetName()).data(), region, imed);
-
- }
-}
-
-G4int FGeometryInit::GetMedium(int region) const
-{
- return fRegionMediumMap[region];
-}
-
-
-void FGeometryInit::SetMediumFromName(const char* volName, int medium, int volid)
- {
- char name4[5];
- char tmp[5];
- strncpy(tmp, volName, 4);
- tmp[4]='\0';
- fNRegions = 0;
-
- for (RegionIterator i = fRegionVolumeMap.begin();
- i != fRegionVolumeMap.end();
- i++) {
- fNRegions++;
- //Get info in the map
- G4VPhysicalVolume* ptrVol = (*i).first;
- strncpy(name4, (ptrVol->GetName()).data(), 4);
- name4[4]='\0';
- for (int j = 0; j < 4; j++) {
- if (name4[j] == '\0') {
- for (int k = j; k < 4; k++) {
- name4[k] = ' ';
- }
- break;
- }
- }
- if (! strncmp(name4, tmp, 4)) {
- fMediumVolumeMap[ptrVol] = medium;
- fVolIdVolumeMap[ptrVol] = volid;
- }
- }
-}
-
-
-
-////////////////////////////////////////////////////////////////////////
-//
-void FGeometryInit::BuildMaterialTables() {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
-#endif
-
- //some terminal printout also
- G4cout << "\t* Storing information..." << G4endl;
-
- //The logic is the folloing:
- //Get the Material Table and:
- // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
- // 2) For each single element material build a material equivalent
- // 3) For the rest:
- // 3.a) Build materials for each not already known element
- // 3.b) Build the compound out of them
-
- //Get the Material Table and iterate
- const G4MaterialTable* matTable = G4Material::GetMaterialTable();
- for (MatTableIterator i = matTable->begin(); i != matTable->end(); i++) {
-
- //Get some basic material information
- G4Material* material = (*i);
- G4String matName = material->GetName();
- const G4double matDensity = material->GetDensity();
- const G4int nMatElements = material->GetNumberOfElements();
-
- G4cout << "\t\t+ " << matName
- << ": dens. = " << matDensity/(g/cm3) << "g/cm3"
- << ", nElem = " << nMatElements << G4endl;
-
- // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
- // FlukaMaterial* is 0 in that case
- if (matDensity <= 1.00e-10*g/cm3) {
- G4FlukaMaterialMap[material] = 0;
- G4cout << "\t\t Stored as vacuum" << G4endl;
- }
- // 2) For each single element material build a material equivalent
- else if (nMatElements == 1) {
-
- FlukaMaterial *flukamat =
- BuildFlukaMaterialFromElement(material->GetElement(0),
- matDensity);
-
- G4FlukaMaterialMap[material] = flukamat;
- G4cout << "\t\t Stored as " << flukamat->GetRealName() << G4endl;
-
- } //else if (material->GetNumberOfElements() == 1)
-
- // 3) For the rest:
- // 3.a) Build materials for each not already known element
- // 3.b) Build the compound out of them
- else {
- FlukaCompound* flukacomp =
- BuildFlukaCompoundFromMaterial(material);
- G4FlukaCompoundMap[material] = flukacomp;
- G4cout << "\t\t Stored as " << flukacomp->GetRealName() << G4endl;
- } //else for case 3)
- } //for (materials)
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
-#endif
-}
-
-FlukaMaterial*
-FGeometryInit::BuildFlukaMaterialFromElement(const G4Element* element,
- G4double matDensity) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromElement()"
- << G4endl;
-#endif
-
- //Get element and its properties
- G4String elemName(ToFlukaString(element->GetName()));
-
- FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(elemName);
- if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
- //Check for isotopes
- G4int nIsotopes = element->GetNumberOfIsotopes();
- if (nIsotopes == 0) {
- G4double elemA = element->GetA()/g;
- G4double elemZ = element->GetZ();
-
- if (elemA != G4int(elemA) && elemZ != G4int(elemZ)) {
- G4cout << "WARNING: Element \'" << elemName
- << "\' has non integer Z (" << elemZ << ") or A ("
- << elemA << ")"
- << G4endl;
- }
-
- flukamat = new FlukaMaterial(elemName,
- G4int(elemZ),
- elemA,
- matDensity/(g/cm3));
- }
- else if (nIsotopes == 1) {
- const G4Isotope* isotope = element->GetIsotope(0);
- flukamat = BuildFlukaMaterialFromIsotope(isotope, matDensity);
- }
- else {
- FlukaCompound *flucomp = BuildFlukaCompoundFromElement(element,
- matDensity);
- flukamat = flucomp->GetFlukaMaterial();
- }
- }
-#ifdef G4GEOMETRY_DEBUG
- else {
- G4cout << "INFO: Element \'" << elemName
- << "\' already exists in the DB. It will not be recreated."
- << G4endl;
- }
-#endif
-
- return flukamat;
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FGeometryInit::BuildFlukaMaterialFromElement()"
- << G4endl;
-#endif
-}
-
-FlukaMaterial*
-FGeometryInit::BuildFlukaMaterialFromIsotope(const G4Isotope* isotope,
- G4double matDensity) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()"
- << G4endl;
-#endif
- G4String isoName(ToFlukaString(isotope->GetName()));
- FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(isoName);
- if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
- G4int isoZ = isotope->GetZ();
- G4double isoA = (isotope->GetA())/(g);
- G4int isoN = isotope->GetN();
- flukamat = new FlukaMaterial(isoName,
- isoZ,
- isoA,
- matDensity/(g/cm3),
- isoN);
- }
-
- return flukamat;
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()"
- << G4endl;
-#endif
-}
-
-FlukaCompound*
-FGeometryInit::BuildFlukaCompoundFromMaterial(const G4Material* material) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()"
- << G4endl;
-#endif
- //Material properties
- const G4double* elemFractions = material->GetFractionVector();
- const G4int nMatElements = material->GetNumberOfElements();
- const G4double matDensity = material->GetDensity();
- G4String matName(ToFlukaString(material->GetName()));
- FlukaCompound* flukacomp = new FlukaCompound(matName, matDensity/(g/cm3),
- nMatElements);
- for (G4int i = 0; i < nMatElements; i++) {
- FlukaMaterial *flukamat =
- BuildFlukaMaterialFromElement(material->GetElement(i), 0.0);
-
- flukacomp->AddElement(flukamat->GetIndex(), -elemFractions[i]);
-
- } //for (elements)
-
- return flukacomp;
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()"
- << G4endl;
-#endif
-}
-
-FlukaCompound*
-FGeometryInit::BuildFlukaCompoundFromElement(const G4Element* element,
- G4double matDensity) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromElement()"
- << G4endl;
-#endif
- G4int nIsotopes = element->GetNumberOfIsotopes();
- //fraction of nb of atomes per volume (= volume fraction?)
- const G4double* isoAbundance = element->GetRelativeAbundanceVector();
- G4String elemName(ToFlukaString(element->GetName()));
-
- //Material properties
- FlukaCompound* flukacomp = new FlukaCompound(elemName, matDensity/(g/cm3),
- nIsotopes);
- for (G4int i = 0; i < nIsotopes; i++) {
- FlukaMaterial *flukamat =
- BuildFlukaMaterialFromIsotope(element->GetIsotope(i), 0.0);
-
- flukacomp->AddElement(flukamat->GetIndex(), isoAbundance[i]);
-
- } //for (elements)
-
- return flukacomp;
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromElement()"
- << G4endl;
-#endif
-}
-
-void FGeometryInit::PrintMaterialTables(std::ostream& os) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
-#endif
- //Print Header
- PrintHeader(os, "GEANT4 MATERIALS AND COMPOUNDS");
-
- //And some more stuff
- size_t nIsotopes = G4Isotope::GetNumberOfIsotopes();
- size_t nElements = G4Element::GetNumberOfElements();
- size_t nMaterials = G4Material::GetNumberOfMaterials();
-
- os << "* In Geant4 there are " << nMaterials << " materials" << G4endl;
- os << "* In Geant4 there are " << nElements << " elements" << G4endl;
- os << "* In Geant4 there are " << nIsotopes << " isotopes" << G4endl;
-
- //Materials
- G4cout << "\t* Printing FLUKA materials..." << G4endl;
- FlukaMaterial::PrintMaterialsByIndex(os);
- //FlukaMaterial::PrintMaterialsByName(os);
-
- //Compounds
- G4cout << "\t* Printing FLUKA compounds..." << G4endl;
- FlukaCompound::PrintCompounds(os);
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
-#endif
-}
-
-////////////////////////////////////////////////////////////////////////
-//
-void FGeometryInit::PrintAssignmat(std::ostream& os) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
-#endif
-
- //Find number of Volumes in physical volume store
- G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
- unsigned int numVol = pVolStore->size();
-
- G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore
- << ") has " << numVol << " volumes. " << G4endl;
-
- G4cout << "\t* Printing ASSIGNMAT..." << G4endl;
-
-
- PrintHeader(os,"GEANT4 MATERIAL ASSIGNMENTS");
- for(unsigned int l=0; l < numVol; l++) {
-
- //Get each of the physical volumes
- G4VPhysicalVolume * physicalvol = (*pVolStore)[l];
-
- //Get index for that volume
- G4int iFlukaRegion = fRegionVolumeMap[physicalvol];
-
- //Find G4 material and navigate to its fluka compound/material
- G4LogicalVolume * logicalVol = physicalvol->GetLogicalVolume();
- G4Material* material = logicalVol->GetMaterial();
- G4int matIndex = 2;
- if (G4FlukaCompoundMap[material])
- matIndex = G4FlukaCompoundMap[material]->GetIndex();
- if (G4FlukaMaterialMap[material])
- matIndex = G4FlukaMaterialMap[material]->GetIndex();
-
- //Find if there is a magnetic field in the region
- //check if Magnetic Field is present in the region
- G4double flagField = 0.0;
- G4FieldManager * pMagFieldMan = logicalVol->GetFieldManager();
- if(pMagFieldMan && pMagFieldMan->GetDetectorField())
- flagField = 1.0;
-
- //Print card
- os << setw10 << "ASSIGNMAT ";
- os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
- os << setw10 << setfixed << G4double(matIndex);
- os << setw10 << setfixed << G4double(iFlukaRegion);
- os << setw10 << "0.0";
- os << setw10 << setfixed << flagField;
- os << G4endl;
- }
-
-
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
-#endif
-}
-
-
-void FGeometryInit::PrintMagneticField(std::ostream& os) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FGeometryInit::PrintMagneticField()" << G4endl;
-#endif
-
- G4cout << "\t* Printing Magnetic Field..." << G4endl;
-
- if(fTransportationManager->GetFieldManager()->DoesFieldExist()) {
-
- //get magnetic field pointer
- const G4Field * pMagField =
- fTransportationManager->GetFieldManager()->GetDetectorField();
-
-
- if(pMagField) {
- //Check if it can be made a uniform magnetic field
- const G4UniformMagField *pUnifMagField =
- dynamic_cast<const G4UniformMagField*>(pMagField);
- if(pUnifMagField) {
- G4double B[3];
- G4double point[4]; //it is not really used
- pUnifMagField->GetFieldValue(point,B);
-
- //write MGNFIELD card
- PrintHeader(os,"GEANT4 MAGNETIC FIELD");
- os << setw10 << "MGNFIELD ";
- os << setw10 << "";
- os << setw10 << "";
- os << setw10 << "";
- os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
- os << setw10 << setfixed
- << std::setprecision(4) << B[0]
- << setw10 << B[1]
- << setw10 << B[2]
- << G4endl;
- }
- else {
- G4cout << "WARNING: No Uniform Magnetic Field found." << G4endl;
- G4cout << " Manual intervention might be needed." << G4endl;
- }
- }
- else
- G4cout << "\t No detector field found... " << G4endl;
- } // end if magnetic field
- else
- G4cout << "\t No field found... " << G4endl;
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FGeometryInit::PrintMagneticField()" << G4endl;
-#endif
-}
-
-int FGeometryInit::CurrentVolID(int ir, int& copyNo)
-{
- if (ir == 0)
- {
- copyNo = -1;
- return -1;
- }
-
- G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
- G4VPhysicalVolume * physicalvol = (*pVolStore)[ir- 1];
-
- if (physicalvol) {
- copyNo = physicalvol->GetCopyNo();
- } else {
- copyNo = -1;
- return -1;
- }
-
-
- int id = fVolIdVolumeMap[physicalvol];
- return id;
-}
-
-int FGeometryInit::CurrentVolOffID(int ir, int off, int& copyNo)
-{
- if (off == 0) return CurrentVolID(ir, copyNo);
-
- G4PhysicalVolumeStore* pVolStore = G4PhysicalVolumeStore::GetInstance();
- G4VPhysicalVolume* physicalvol = (*pVolStore)[ir- 1];
- G4VPhysicalVolume* mother = physicalvol;
-
- int index;
-//============================================================================
- if (mother) {
- // Check touchable depth
- //
- if (ptrTouchHist->GetHistoryDepth() < off) {
- mother = 0;
- } else {
- // Get the off-th mother
- index = ptrTouchHist->GetHistoryDepth() - off;
- // in the touchable history volumes are ordered
- // from top volume up to mother volume;
- // the touchable volume is not in the history
- mother = ptrTouchHist->GetHistory()->GetVolume(index);
- }
- }
-//============================================================================
-
- int id;
-
- if (!mother) {
- G4cout << "Flugg FGeometryInit::CurrentVolOffID mother not found" << G4endl;
- id = -1;
- copyNo = -1;
- } else {
- copyNo = mother ->GetCopyNo();
- id = fVolIdVolumeMap[mother];
- }
- return id;
-}
-
-void FGeometryInit::Gmtod(double* xm, double* xd, int iflag)
-{
-// Transforms a position from the world reference frame
-// to the current volume reference frame.
-//
-// Geant3 desription:
-// ==================
-// Computes coordinates XD (in DRS)
-// from known coordinates XM in MRS
-// The local reference system can be initialized by
-// - the tracking routines and GMTOD used in GUSTEP
-// - a call to GMEDIA(XM,NUMED)
-// - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
-// (inverse routine is GDTOM)
-//
-// If IFLAG=1 convert coordinates
-// IFLAG=2 convert direction cosinus
-//
-// ---
- FluggNavigator * ptrNavig = getNavigatorForTracking();
- //setting variables (and dimension: Fluka uses cm.!)
- G4ThreeVector pGlob(xm[0],xm[1],xm[2]);
- G4ThreeVector pLoc;
-
- if (iflag == 1) {
- pGlob *= 10.0; // in mm
-// change because of geant 4 6.0
-// pLoc = ptrNavig->ComputeLocalPoint(pGlob);
- pLoc = ptrNavig->GetGlobalToLocalTransform().TransformPoint(pGlob);
-
- pLoc /= 10.0; // in cm
- } else if (iflag == 2) {
- pLoc =
- ptrNavig->ComputeLocalAxis(pGlob);
- } else {
- G4cout << "Flugg FGeometryInit::Gmtod called with undefined flag" << G4endl;
- }
- xd[0] = pLoc[0]; xd[1] = pLoc[1]; xd[2] = pLoc[2];
-}
-
-void FGeometryInit::Gdtom(double* xd, double* xm, int iflag)
-{
-// Transforms a position from the current volume reference frame
-// to the world reference frame.
-//
-// Geant3 desription:
-// ==================
-// Computes coordinates XM (Master Reference System
-// knowing the coordinates XD (Detector Ref System)
-// The local reference system can be initialized by
-// - the tracking routines and GDTOM used in GUSTEP
-// - a call to GSCMED(NLEVEL,NAMES,NUMBER)
-// (inverse routine is GMTOD)
-//
-// If IFLAG=1 convert coordinates
-// IFLAG=2 convert direction cosinus
-//
-// ---
-
- FluggNavigator * ptrNavig = getNavigatorForTracking();
- G4ThreeVector pLoc(xd[0],xd[1],xd[2]);
-
- G4ThreeVector pGlob;
- if (iflag == 1) {
- pLoc *= 10.0; // in mm
- pGlob = ptrNavig->GetLocalToGlobalTransform().
- TransformPoint(pLoc);
- pGlob /= 10.0; // in cm
- } else if (iflag == 2) {
- pGlob = ptrNavig->GetLocalToGlobalTransform().
- TransformAxis(pLoc);
- } else {
- G4cout << "Flugg FGeometryInit::Gdtom called with undefined flag" << G4endl;
- }
-
- xm[0] = pGlob[0]; xm[1] = pGlob[1]; xm[2] = pGlob[2];
-}
+++ /dev/null
-
-// Flugg tag
-
-// modified 10/IX/99 for including preStepPoint
-// modified 28/IX/99 for delating allocated memory
-// modified 4/X/99 function FreeMemory
-// modified 2/III/00 base class G4TransportationManager included
-// modified 20/III/00 PrintHistories() included
-
-
-#ifndef FGeometryInit_h
-#define FGeometryInit_h 1
-
-//#include "g4std/fstream"
-//#include "g4std/iomanip"
-//#include "g4rw/cstring.h"
-#include "globals.hh"
-
-#include "G4LogicalVolume.hh"
-#include "G4PhysicalVolumeStore.hh"
-#include "G4VPhysicalVolume.hh"
-#include "G4Material.hh"
-#include "G4MaterialTable.hh"
-#include "G4Isotope.hh"
-#include "G4VUserDetectorConstruction.hh"
-#include "G4TouchableHistory.hh"
-#include "G4GeometryManager.hh"
-#include "G4FieldManager.hh"
-#include "G4UniformMagField.hh"
-#include "G4TransportationManager.hh"
-
-
-#include <map>
-
-class FluggNavigator;
-class FlukaMaterial;
-class FlukaCompound;
-
-class FGeometryInit : public G4TransportationManager {
-
-public:
- ~FGeometryInit(); //destructor
- static FGeometryInit *GetInstance();
- inline FluggNavigator *getNavigatorForTracking();
- inline G4FieldManager * getFieldManager();
- inline void setDetConstruction(G4VUserDetectorConstruction* detector);
- inline void setDetector();
- inline void setMotherVolume();
- void createFlukaMatFile();
- void closeGeometry();
-
- void PrintHistories();
- void InitHistories();
- void DeleteHistories();
- void UpdateHistories(const G4NavigationHistory *, G4int);
- inline G4TouchableHistory * GetTouchableHistory();
- inline G4TouchableHistory * GetOldNavHist();
- inline G4TouchableHistory * GetTempNavHist();
-
- void InitHistArray();
- inline void DelHistArray();
- inline G4int * GetHistArray();
-
- void InitJrLtGeantArray();
- inline G4int * GetJrLtGeantArray();
- inline G4int GetLttcFlagGeant();
- void SetLttcFlagGeant(G4int);
- void PrintJrLtGeant();
- void Init();
-
- //Map access methods
- void BuildMediaMap();
- void SetMediumFromName(const char* volName, int med, int volid);
- //G4int GetRegionFromName(const char* volName) const;
- int GetLastMaterialIndex() const;
- G4int GetMedium(int) const;
- int CurrentVolID(int ir, int& copyNo);
- int CurrentVolOffID(int ir, int off, int& copyNo);
- void Gmtod(double* xm, double* xd, int iflag);
- void Gdtom(double* xd, double* xm, int iflag);
-
-protected:
- void BuildRegionsMap();
- void PrintRegionsMap(std::ostream& os);
- void BuildMaterialTables();
- FlukaMaterial* BuildFlukaMaterialFromElement(const G4Element* element,
- G4double matDensity);
- FlukaMaterial* BuildFlukaMaterialFromIsotope(const G4Isotope* isotope,
- G4double matDensity);
- FlukaCompound* BuildFlukaCompoundFromMaterial(const G4Material* material);
- FlukaCompound* BuildFlukaCompoundFromElement(const G4Element* element,
- G4double matDensity);
- void PrintMaterialTables(std::ostream& os);
- void PrintAssignmat(std::ostream& os);
- void PrintMagneticField(std::ostream& os);
-
-private:
- FGeometryInit(); //costructor
-
-private:
- G4VUserDetectorConstruction * fDetector;
- G4FieldManager * fFieldManager;
- G4TransportationManager * fTransportationManager;
- static FGeometryInit *flagInstance;
- G4VPhysicalVolume * myTopNode;
- G4GeometryManager * ptrGeoMan;
- G4int * ptrArray;
- G4TouchableHistory * ptrTouchHist;
- G4TouchableHistory * ptrOldNavHist;
- G4TouchableHistory * ptrTempNavHist;
- G4int * ptrJrLtGeant;
- G4int flagLttcGeant;
- G4int fNRegions;
- int* fRegionMediumMap;
-
- std::map<G4VPhysicalVolume*, int, std::less<G4VPhysicalVolume*> > fRegionVolumeMap;
- std::map<G4VPhysicalVolume*, int, std::less<G4VPhysicalVolume*> > fMediumVolumeMap;
- std::map<G4VPhysicalVolume*, int, std::less<G4VPhysicalVolume*> > fVolIdVolumeMap;
-
- std::map<G4Material*, FlukaMaterial*, std::less<G4Material*> > G4FlukaMaterialMap;
- std::map<G4Material*, FlukaCompound*, std::less<G4Material*> > G4FlukaCompoundMap;
- //G4int NOfMaterials;
-};
-
-typedef std::map<G4VPhysicalVolume*, int, std::less<G4VPhysicalVolume*> >::const_iterator RegionIterator;
-typedef std::vector<G4Material*>::const_iterator MatTableIterator;
-
-
-//Include the file with the inline methods
-#include "FGeometryInit.icc"
-
-#endif
+++ /dev/null
-//#include <iostream.h>
-#include <FluggNavigator.hh>
-
-FluggNavigator* FGeometryInit::getNavigatorForTracking() {
- G4Navigator* g4nav = fTransportationManager->GetNavigatorForTracking();
- return ((FluggNavigator*) g4nav);
-}
-
-void FGeometryInit::setDetConstruction(G4VUserDetectorConstruction* detector) {
- fDetector = detector;;
-}
-
-void FGeometryInit::setDetector() {
- myTopNode = fDetector->Construct();
-}
-
-void FGeometryInit::setMotherVolume() {
- getNavigatorForTracking()->SetWorldVolume(myTopNode);
-}
-
-G4FieldManager * FGeometryInit::getFieldManager() {
- return fTransportationManager->GetFieldManager();
-}
-
-void FGeometryInit::DelHistArray() {
- delete[] ptrArray;
-}
-
-G4int * FGeometryInit::GetHistArray() {
- return ptrArray;
-}
-
-G4int * FGeometryInit::GetJrLtGeantArray() {
- return ptrJrLtGeant;
-}
-
-
-G4int FGeometryInit::GetLttcFlagGeant() {
- return flagLttcGeant;
-}
-
-G4TouchableHistory * FGeometryInit::GetTouchableHistory() {
- return ptrTouchHist;
-}
-
-G4TouchableHistory * FGeometryInit::GetOldNavHist() {
- return ptrOldNavHist;
-}
-
-G4TouchableHistory * FGeometryInit::GetTempNavHist() {
- return ptrTempNavHist;
-}
+++ /dev/null
-//
-// ********************************************************************
-// * DISCLAIMER *
-// * *
-// * The following disclaimer summarizes all the specific disclaimers *
-// * of contributors to this software. The specific disclaimers,which *
-// * govern, are listed with their locations in: *
-// * http://cern.ch/geant4/license *
-// * *
-// * Neither the authors of this software system, nor their employing *
-// * institutes,nor the agencies providing financial support for this *
-// * work make any representation or warranty, express or implied, *
-// * regarding this software system or assume any liability for its *
-// * use. *
-// * *
-// * This code implementation is the intellectual property of the *
-// * GEANT4 collaboration. *
-// * By copying, distributing or modifying the Program (or any work *
-// * based on the Program) you indicate your acceptance of this *
-// * statement, and all its terms. *
-// ********************************************************************
-//
-//
-
-// GEANT4 tag $ Name: $
-//
-// class FluggNavigator Implementation Paul Kent July 95/96
-
-#include "FluggNavigator.hh"
-#include "G4ios.hh"
-#include <iomanip.h>
-
-#ifdef G4GEOMETRY_DEBUG
-# define G4DEBUG_NAVIGATION 1
-# define G4VERBOSE 1
-#endif
-
-FluggNavigator::FluggNavigator() :
- G4Navigator()
-{
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "==> Flugg FluggNavigator constructor" << G4endl;
-#endif
-
- ResetStackAndState();
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "<== Flugg FluggNavigator constructor" << G4endl;
-#endif
-}
-
-void FluggNavigator::UpdateNavigatorHistory(const G4NavigationHistory* newNavHistory)
-{
-#ifdef G4GEOMETRY_DEBUG
- cout << "==> Flugg FluggNavigator::UpdateNavigatorHistory(" << newNavHistory
- << ")" << endl;
- cout << "\t+fHistory=" << fHistory << ") ..." << G4endl;
-#endif
-
- ResetStackAndState();
- fHistory = *newNavHistory;
- SetupHierarchy();
-
-#ifdef G4GEOMETRY_DEBUG
- cout << "<== Flugg FluggNavigator::UpdateNavigatorHistory(" << newNavHistory
- << ")" << endl;
-#endif
-}
+++ /dev/null
-#ifndef FLUGGNAVIGATOR_H
-#define FLUGGNAVIGATOR_H 1
-
-//To have access to the private data members, let's play a dirty trick
-#define private protected
-#include "G4Navigator.hh"
-#undef private
-
-
-class FluggNavigator : public G4Navigator
-{
- public:
-
- friend std::ostream& operator << (std::ostream &os, const FluggNavigator &n);
-
- FluggNavigator();
- // Constructor - initialisers and setup.
-
- virtual ~FluggNavigator() {}
- // Destructor. No actions.
-
- // flugg member function: reinitialization of navigator history with
- // secondary particle history banked on fluka side
- void UpdateNavigatorHistory(const G4NavigationHistory* newNavHistory);
-
-};
-
-#endif
+++ /dev/null
-#include "FlukaCompound.hh"
-#include "G4ios.hh"
-#include "WrapUtils.hh"
-
-FlukaCompoundsTable FlukaCompound::fFlukaCompounds;
-
-FlukaCompound::FlukaCompound(const G4String& name,
- G4double density,
- G4int nelem):
- fNMaterials(nelem),
- fNAdded(0) {
- //Initialise arrays for element indices and fractions
- fElIndex = new G4int[nelem];
- fFraction = new G4double[nelem];
- for (G4int i = 0; i < nelem; i++) {
- fElIndex[i] = 0;
- fFraction[i] = 0;
- }
-
- //Build the associated material
- fFlukaMaterial = new FlukaMaterial(name, 0, 0, density);
-
- //If it already exists it will have got a new name
- G4String testname(fFlukaMaterial->GetRealName());
- if (testname != name) {
- G4cerr << "INFO: Found two materials with the same name! ("
- << name << ")" << G4endl;
- G4cerr << " Renaming to " << testname << G4endl;
- }
- fFlukaCompounds[testname] = this;
-}
-
-FlukaCompound::~FlukaCompound() {
- delete[] fElIndex;
- delete[] fFraction;
-}
-
-
-void FlukaCompound::AddElement(G4int index, G4double fraction) {
- if (fNAdded < fNMaterials) {
- fElIndex[fNAdded] = index;
- fFraction[fNAdded] = fraction;
- fNAdded++;
- }
- else {
- G4cerr << "ERROR: Trying to add to many elements to compound \'"
- << GetName() << "\'!" << G4endl;
- G4cerr << " Last element not added!!!" << G4endl;
- }
-}
-
-std::ostream& FlukaCompound::PrintCompounds(std::ostream& os) {
- PrintHeader(os, "COMPOUNDS");
-
- for (FlukaCompoundsIterator i = fFlukaCompounds.begin();
- i != fFlukaCompounds.end();
- i++) {
- FlukaCompound* flucomp = (*i).second;
- os << *flucomp;
- }
-
- return os;
-}
-
-std::ostream& operator<<(std::ostream& os, const FlukaCompound& flucomp) {
- G4int nmats = flucomp.GetNMaterials();
- G4String matName = flucomp.GetName();
- G4String matRealName = flucomp.GetRealName().substr(0,8);
- //Some comment
- os << "* " << matName << " COMPOUND (" << nmats << ")" << G4endl;
-
- //Material card
- //os << *(flucomp.GetFlukaMaterial());
-
- //The card
- G4int counttothree = 0;
- os << setw10 <<"COMPOUND ";
- for (G4int i = 0; i < nmats; i++) {
- os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
- os << setw10
- << setfixed
- << std::setprecision(6)
- << flucomp.GetMaterialFraction(i);
- os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
- os << setw10
- << setfixed
- << std::setprecision(1)
- << G4double(flucomp.GetMaterialIndex(i));
- counttothree++;
- if (counttothree == 3 ) {
- os << matRealName;
- os << G4endl;
- if ( (i+1) != nmats)
- os << setw10 <<"COMPOUND ";
- counttothree = 0;
- }
- }
-
- if ((nmats % 3) != 0) {
- //Unless we have 3, 6, 9... submaterials we need to put some empty
- //space and the compound name
- for (G4int i = 0; i < (3 - (nmats % 3)); i++)
- os << setw10 << " " << setw10 << " ";
- os << matRealName;
- os << G4endl;
- }
-
- return os;
-}
+++ /dev/null
-#ifndef FLUKACOMPOUND_HH
-#define FLUKACOMPOUND_HH 1
-
-#include "G4String.hh"
-#include "FlukaMaterial.hh"
-
-#include <map>
-
-class FlukaCompound;
-typedef std::map<G4String, FlukaCompound*, std::less<G4String> > FlukaCompoundsTable;
-typedef std::map<G4String, FlukaCompound*, std::less<G4String> >::const_iterator FlukaCompoundsIterator;
-
-class FlukaCompound {
-public:
-
- //Constructor
- FlukaCompound(const G4String& name, G4double density, G4int nElem);
- virtual ~FlukaCompound();
-
- //Getters
- // * Fluka name
- G4String GetName() const {return fFlukaMaterial->GetName();}
- // * Real name if the fluka name is duplicated
- G4String GetRealName() const {return fFlukaMaterial->GetRealName();}
- // * Density
- G4double GetDensity() const {return fFlukaMaterial->GetDensity();}
- // * Index (it comes from the material)
- G4int GetIndex() const {return fFlukaMaterial->GetIndex();}
- // * Number of materials, and index and fraction for each material
- G4int GetNMaterials() const {return fNMaterials;}
- G4int GetMaterialIndex(G4int i) const {return fElIndex[i];}
- G4double GetMaterialFraction(G4int i) const {return fFraction[i];}
-
- // * Associated material
- const FlukaMaterial* GetFlukaMaterial() const {return fFlukaMaterial;}
- FlukaMaterial* GetFlukaMaterial() {return fFlukaMaterial;}
-
- //Setters
- void SetName(const G4String& n) {fFlukaMaterial->SetName(n);}
- void SetDensity(G4double d) {fFlukaMaterial->SetDensity(d);}
-
- //Other
- void AddElement(G4int index, G4double fraction);
-
- //Static
- static inline const FlukaCompoundsTable* GetCompoundTable();
- static inline const FlukaCompound* GetFlukaCompound(const G4String& name);
- static std::ostream& PrintCompounds(std::ostream& os);
-
-public:
- G4int fNMaterials; //Number of elements in total
- G4int fNAdded; //Number of elements added
- G4int* fElIndex; //Array of element indices
- G4double* fFraction; //Array of element fracions
-
- FlukaMaterial* fFlukaMaterial; //Each compound has a "dummy" mat associated
-
- static FlukaCompoundsTable fFlukaCompounds;
-
-};
-
-inline const FlukaCompoundsTable* FlukaCompound::GetCompoundTable() {
- return &fFlukaCompounds;
-}
-
-inline const FlukaCompound* FlukaCompound::GetFlukaCompound(const G4String& name) {
- return fFlukaCompounds[name];
-}
-
-
-//Ostream operator
-std::ostream& operator<<(std::ostream& os, const FlukaCompound& flucomp);
-
-#endif
-
+++ /dev/null
-#include "FlukaLowMat.hh"
-#include "WrapUtils.hh"
-
-FlukaLowMat::FlukaLowMat(const G4String& name, FlukaMaterial* fmat):
- fName(name),
- fFlukaMaterial(fmat){
-}
-
-std::ostream& operator<<(std::ostream& os, const FlukaLowMat& flowmat){
- os << setw10 << "LOW-MAT ";
- os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
- os << setw10 << setfixed << std::setprecision(1)
- << G4double(flowmat.GetIndex());
- os << setw10 << " "
- << setw10 << " "
- << setw10 << " "
- << setw10 << " "
- << setw10 << " ";
- os << flowmat.GetFlukaMaterial()->GetName() << G4endl;
-
- return os;
-}
+++ /dev/null
-#ifndef FLUKALOWMAT_HH
-#define FLUKALOWMAT_HH 1
-
-#include "G4String.hh"
-#include "FlukaMaterial.hh"
-
-class FlukaLowMat {
-public:
-
- //Constructor
- FlukaLowMat(const G4String& name, FlukaMaterial*);
-
- //Getters
- G4String GetName() const {return fName;}
- inline G4int GetIndex() const;
- const FlukaMaterial* GetFlukaMaterial() const {return fFlukaMaterial;}
-
- //Setters
- void SetName(const G4String& n) {fName = n;}
-
-public:
- G4String fName;
- FlukaMaterial* fFlukaMaterial;
-
-};
-
-
-inline G4int FlukaLowMat::GetIndex() const {return fFlukaMaterial->GetIndex();}
-
-std::ostream& operator<<(std::ostream& os, const FlukaLowMat& flowmat);
-
-#endif
-
+++ /dev/null
-#include "FlukaMaterial.hh"
-#include "FlukaLowMat.hh"
-#include "WrapUtils.hh"
-#include "G4ios.hh"
-
-FlukaMaterialsTable FlukaMaterial::fFlukaMaterials;
-FlukaMaterialsIndexTable FlukaMaterial::fFlukaIndexMaterials;
-
-FlukaMaterial::FlukaMaterial(const G4String& name,
- G4int Z, G4double A,
- G4double density,
- G4int N):
- fName(name),
- fZ(Z),
- fA(A),
- fDensity(density),
- fN(N),
- fFlukaLowMat(0) {
-
- G4String testname(name);
- G4int matrep = 1;
- while (fFlukaMaterials[testname] && matrep < 100) {
- matrep++;
- char smatrep[3];
- sprintf(smatrep,"%.2d",matrep);
-
- testname = name;
- if (testname.length() <= 6)
- testname += smatrep;
- else
- testname.replace(6,testname.length()-6, smatrep, 2);
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "INFO: We found material \'" << name << " previously defined."
- << G4endl;
- G4cout << " Checking if \'" << testname << "\' exists." << G4endl;
-#endif
- }
-
- if (matrep > 99) {
- G4cerr << "ERROR: Too many materials with the same name. Exiting!"
- << G4endl;
- abort();
- }
-
- fFlukaMaterials[testname] = this;
- if (name != testname)
- AddLowMat(testname);
- fIndex = fFlukaMaterials.size() + 2;
- fFlukaIndexMaterials[fIndex] = this;
-}
-
-FlukaMaterial::~FlukaMaterial() {
- delete fFlukaLowMat;
-}
-
-void FlukaMaterial::AddLowMat(const G4String& name) {
- fFlukaLowMat = new FlukaLowMat(name, this);
-}
-
-
-G4String FlukaMaterial::GetRealName() const {
- if (fFlukaLowMat)
- return fFlukaLowMat->GetName();
- return GetName();
-}
-
-std::ostream& FlukaMaterial::PrintMaterialsByName(std::ostream& os) {
- PrintHeader(os, "MATERIALS");
- for (FlukaMaterialsIterator i = fFlukaMaterials.begin();
- i != fFlukaMaterials.end();
- i++) {
- FlukaMaterial* flumat = (*i).second;
-
- //if (flumat->GetZ()) //Skip materials that describe only compounds
- os << *flumat;
- }
- return os;
-}
-
-std::ostream& FlukaMaterial::PrintMaterialsByIndex(std::ostream& os) {
- PrintHeader(os, "MATERIALS");
- for (FlukaMaterialsIndexIterator i = fFlukaIndexMaterials.begin();
- i != fFlukaIndexMaterials.end();
- i++) {
- FlukaMaterial* flumat = (*i).second;
-
- //if (flumat->GetZ()) //Skip materials that describe only compounds
- os << *flumat;
- }
- return os;
-}
-
-std::ostream& operator<<(std::ostream& os, const FlukaMaterial& material){
- os << setw10 << "MATERIAL ";
-
- os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
- G4double Z = G4double(material.GetZ());
- if (Z <= 0)
- os << setw10 << " ";
- else
- os << setw10
- << setfixed
- << std::setprecision(1)
- << Z;
-
- G4double A = material.GetA();
- if (A <= 0)
- os << setw10 << " ";
- else
- os << setw10 << std::setprecision(3)
- << A;
-
- G4double density = material.GetDensity();
- if (density <=0)
- density = 0.999;
- os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
- os << setw10
- << setscientific
- << std::setprecision(3)
- << density;
-
- os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
- os << setw10
- << setfixed
- << std::setprecision(1)
- << G4double(material.GetIndex());
-
- os << setw10 << " ";
- if (material.GetN())
- os << setw10 << G4double(material.GetN());
- else
- os << setw10 << " ";
-
-
- os << material.GetRealName().substr(0,8) << G4endl;
-
- if (material.GetLowMat() && material.GetZ() != 0)
- os << *(material.GetLowMat());
-
- return os;
-}
+++ /dev/null
-#ifndef FLUKAMATERIAL_HH
-#define FLUKAMATERIAL_HH 1
-
-#include "G4String.hh"
-
-class FlukaLowMat;
-
-#include <map>
-
-class FlukaMaterial;
-typedef std::map<G4String, FlukaMaterial*, std::less<G4String> > FlukaMaterialsTable;
-typedef std::map<G4String, FlukaMaterial*, std::less<G4String> >::const_iterator FlukaMaterialsIterator;
-typedef std::map<G4int, FlukaMaterial*, std::less<G4int> > FlukaMaterialsIndexTable;
-typedef std::map<G4int, FlukaMaterial*, std::less<G4int> >::const_iterator FlukaMaterialsIndexIterator;
-
-class FlukaMaterial {
-public:
-
- //Constructor
- FlukaMaterial(const G4String& name,
- G4int Z, G4double A,
- G4double density,
- G4int N=0);
- virtual ~FlukaMaterial();
-
- //Getters
- G4int GetZ() const {return fZ;}
- G4double GetA() const {return fA;}
- G4double GetDensity() const {return fDensity;}
- G4int GetN() const {return fN;}
- G4String GetName() const {return fName;}
- G4String GetRealName() const;
- G4int GetIndex() const {return fIndex;}
- const FlukaLowMat* GetLowMat() const {return fFlukaLowMat;}
-
-
- //Setters
- void SetZ(G4int Z) {fZ = Z;}
- void SetA(G4double A) {fA = A;}
- void SetDensity(G4double d) {fDensity = d;}
- void SetN(G4int N) {fN = N;}
- void SetName(const G4String& n) {fName = n;}
-
- //Static
- static const FlukaMaterialsTable* GetMaterialTable() {
- return &fFlukaMaterials;}
- static FlukaMaterial* GetFlukaMaterial(const G4String& name) {
- return fFlukaMaterials[name];}
- static std::ostream& PrintMaterialsByName(std::ostream& os);
- static std::ostream& PrintMaterialsByIndex(std::ostream& os);
-
-protected:
- //Other
- void AddLowMat(const G4String& name);
-
-public:
- G4String fName;
- G4int fZ;
- G4double fA;
- G4double fDensity;
- G4int fN; //Isotope number
- FlukaLowMat* fFlukaLowMat;
-
- size_t fIndex;
-
- static FlukaMaterialsTable fFlukaMaterials; //Map of name/material
- static FlukaMaterialsIndexTable fFlukaIndexMaterials; //Map index/material
-
-};
-
-std::ostream& operator<<(std::ostream& os, const FlukaMaterial& material);
-
-
-#endif
+++ /dev/null
-//---------------------------------------------------------------
-//
-// dummy G4FastSimulationManager.cc
-//
-//---------------------------------------------------------------
-
-#include "G4FastSimulationManager.hh"
-#include "G4LogicalVolume.hh"
-
-
-G4FastSimulationManager::G4FastSimulationManager(G4LogicalVolume* logVol)
-{;}
-
-G4FastSimulationManager::~G4FastSimulationManager()
-{;}
+++ /dev/null
-//---------------------------------------------------------------
-//
-// dummy G4FastSimulationManager.hh
-//
-//---------------------------------------------------------------
-
-
-#ifndef G4FastSimulationManager_h
-#define G4FastSimulationManager_h 1
-
-#include "G4LogicalVolume.hh"
-
-class G4FastSimulationManager
-{
-public:
- G4FastSimulationManager(G4LogicalVolume* logVol);
- ~G4FastSimulationManager();
-};
-
-#endif
+++ /dev/null
-// This code implementation is the intellectual property of
-// the GEANT4 collaboration.
-//
-// By copying, distributing or modifying the Program (or any work
-// based on the Program) you indicate your acceptance of this statement,
-// and all its terms.
-//
-
-// GEANT4 tag
-//
-
-#include "G4SDManager.hh"
-
-G4SDManager::G4SDManager()
-{
-}
-
-G4SDManager::~G4SDManager()
-{
-}
-
+++ /dev/null
-// This code implementation is the intellectual property of
-// the GEANT4 collaboration.
-//
-// By copying, distributing or modifying the Program (or any work
-// based on the Program) you indicate your acceptance of this statement,
-// and all its terms.
-//
-
-// GEANT4 tag
-//
-
-#ifndef G4SDManager_h
-#define G4SDManager_h 1
-
-class G4SDManager
-{
- protected:
- G4SDManager();
-
- public:
- ~G4SDManager();
-
-};
-
-
-
-
-#endif
-
+++ /dev/null
-// This code implementation is the intellectual property of
-// the GEANT4 collaboration.
-//
-// By copying, distributing or modifying the Program (or any work
-// based on the Program) you indicate your acceptance of this statement,
-// and all its terms.
-//
-
-// GEANT4 tag
-//
-
-#include "G4VUserDetectorConstruction.hh"
-#include "G4VPhysicalVolume.hh"
-
-G4VUserDetectorConstruction::G4VUserDetectorConstruction()
-{;}
-
-G4VUserDetectorConstruction::~G4VUserDetectorConstruction()
-{;}
-
+++ /dev/null
-// This code implementation is the intellectual property of
-// the GEANT4 collaboration.
-//
-// By copying, distributing or modifying the Program (or any work
-// based on the Program) you indicate your acceptance of this statement,
-// and all its terms.
-//
-
-// GEANT4 tag
-//
-
-#ifndef G4VUserDetectorConstruction_h
-#define G4VUserDetectorConstruction_h 1
-
-class G4VPhysicalVolume;
-
-// class description:
-//
-// This is the abstract base class of the user's mandatory initialization class
-// for detector setup. It has only one pure virtual method Construct() which is
-// invoked by G4RunManager when it's Initialize() method is invoked.
-// The Construct() method must return the G4VPhysicalVolume pointer which represents
-// the world volume.
-//
-
-class G4VUserDetectorConstruction
-{
- public:
- G4VUserDetectorConstruction();
- virtual ~G4VUserDetectorConstruction();
-
- public:
- virtual G4VPhysicalVolume* Construct() = 0;
-};
-
-#endif
-
-
-
-
+++ /dev/null
-
-// Flugg tag
-
-//
-// NavHistWithCount.hh - Sara Vanini
-// last modified 2/II/'99
-// This class stores a G4NavigationHistory and it adds to it a
-// counter for secondary particles.
-//
-//
-
-
-#ifndef NAVHISTWITHCOUNT_HH
-#define NAVHISTWITHCOUNT_HH
-
-#include "G4NavigationHistory.hh"
-
-
-class NavHistWithCount
-{
-public:
- NavHistWithCount(const G4NavigationHistory &history);
- ~NavHistWithCount();
- void UpdateCount(G4int incrCount);
- G4NavigationHistory* GetNavHistPtr();
- G4int GetCount();
- G4int GetDelateFlag();
- G4int GetCheckInd();
- void SaveCheckInd(G4int);
-
-private:
- G4int count;
- G4NavigationHistory * fhistory;
- G4int fDelate;
- G4int fCheck;
-};
-
-#include "NavHistWithCount.icc"
-
-#endif
-
-
-
+++ /dev/null
-
-// Flugg tag
-
-//
-// NavHistWithCount.icc - Sara Vanini
-// last modified: 2/II/1999
-//
-// NavHistWithCount.hh inline implementation
-//
-//
-
-inline NavHistWithCount::NavHistWithCount(const G4NavigationHistory &history)
- {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "NavHistWithCount::NavHistWithCount created" << G4endl;
-#endif
- fhistory = new G4NavigationHistory(history);
- fDelate=0;
- count=0;
- }
-
-inline NavHistWithCount::~NavHistWithCount()
- {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "NavHistWithCount::NavHistWithCount deleted" << G4endl;
-#endif
- delete fhistory;
- fDelate=1;
- fhistory=0;
- count=0;
- }
-
-inline void NavHistWithCount::UpdateCount(G4int incrCount)
- {
- count=count+incrCount;
- }
-
-inline G4NavigationHistory * NavHistWithCount::GetNavHistPtr()
- {
- return fhistory;
- }
-
-inline G4int NavHistWithCount::GetCount()
- {
- return count;
- }
-
-inline G4int NavHistWithCount::GetDelateFlag()
- {
- return fDelate;
- }
-
-inline G4int NavHistWithCount::GetCheckInd()
- {
- return fCheck;
- }
-
-inline void NavHistWithCount::SaveCheckInd(G4int index)
- {
- fCheck=index;
- }
-
-
-
+++ /dev/null
-
-// Flugg tag
-
-///////////////////////////////////////////////////////////////////
-//
-// WrapDN.hh - Sara Vanini - 17/XI/98
-//
-// Wrapper for setting DNEAR option on fluka side. Must return 0
-// if user doesn't want Fluka to use DNEAR to compute the
-// step (the same effect is obtained with the GLOBAL (WHAT(3)=-1)
-// card in fluka input), returns 1 if user wants Fluka always to
-// use DNEAR (in this case, be sure that GEANT4 DNEAR is unique,
-// coming from all directions!!!).
-//
-// modified 24.10.01: by I. Hrivnacova
-// functions declarations separated from implementation
-// (moved to Wrappers.hh);
-//
-///////////////////////////////////////////////////////////////////
-
-//#ifndef idnrwr
-//#define idnrwr idnrwr_
-
-#include "Wrappers.hh"
-#include "globals.hh"
-
-G4int idnrwr(const G4int & nreg, const G4int & mlat)
-
-{
-//flag
-#ifdef G4GEOMETRY_DEBUG
- G4cout<<"================== IDNRWR ================="<<G4endl;
-#endif
-
-// returns 0 if user doesn't want Fluka to use DNEAR to compute the
-// step (the same effect is obtained with the GLOBAL (WHAT(3)=-1)
-// card in fluka input), returns 1 if (be sure that GEANT4 DNEAR is unique,
-// coming from all directions!!!) user wants Fluka always to use DNEAR.
-
-return 0;
-}
-//#endif
+++ /dev/null
-
-// Flugg tag
-
-///////////////////////////////////////////////////////////////////
-//
-// WrapG1.hh - Sara Vanini
-//
-// Wrapper for geometry tracking: returns approved step of
-// particle and alla variables that fluka G1 computes.
-//
-// modified 11/III/99 : included jrLtGeant array for storing
-// lattice histories; fixed jump-on-boundaries and back-scattering
-// modified 20/IV/99 : history counters of jrLt array are
-// incremented when created and decremented when check is required:
-// this for not deleting histories twice.
-// modified 18/X/99 : LocateGlobalPointWithinVolume used when position
-// is changed from last located point.
-// modified 1/III/00 : update utilities histories when crossing
-// identical volume boundaries.
-//
-// modified 22/III/00 : fixed LttcFlag and jrLt return values.
-// modified 12/VI/00 : end-step on Boundary bug fixed.
-// modified 5/VII/00 : boundary not seen by G4 geometry bug fixed.
-// modified 24.10.01: by I. Hrivnacova
-// functions declarations separated from implementation
-// (moved to Wrappers.hh);
-//
-/////////////////////////////////////////////////////////////////////
-
-#include "Wrappers.hh"
-#include "WrapUtils.hh"
-#include "FGeometryInit.hh"
-#include "NavHistWithCount.hh"
-#include "G4VPhysicalVolume.hh"
-#include "G4NavigationHistory.hh"
-#include "FluggNavigator.hh"
-#include "G4PhysicalVolumeStore.hh"
-#include "globals.hh"
-
-
-void g1wr(G4double& pSx, G4double& pSy, G4double& pSz, G4double* pV,
- G4int& oldReg, const G4int& oldLttc, G4double& propStep,
- G4int& nascFlag, G4double& retStep, G4int& newReg,
- G4double& saf, G4int& newLttc, G4int& LttcFlag,
- G4double* sLt, G4int* jrLt)
-{
- //flag
-#ifdef G4GEOMETRY_DEBUG
- G4cout<<"============= G1WR =============="<<G4endl;
-#endif
-
- //static G4int count=0;
- //count+=1;
- //G4cout<<"contatore G1="<<count<<G4endl;
-
- ///////////////////////// INPUT ///////////////////////////
- //Geoinit, Navigator, TouchableHistory pointers
- static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
- FluggNavigator * ptrNavig = ptrGeoInit->getNavigatorForTracking();
- G4TouchableHistory * ptrTouchHist = ptrGeoInit->GetTouchableHistory();
- G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
-
- //setting variables (and dimension: Fluka uses cm.!)
- G4ThreeVector partLoc(pSx,pSy,pSz);
- partLoc *= 10.0; // in millimeters!
- static G4ThreeVector partLocOld = partLoc;
-
-// change because of geant4 6.0
-//static G4ThreeVector oldLocalPoint =
-// ptrNavig->ComputeLocalPoint(partLocOld);
- static G4ThreeVector oldLocalPoint =
- ptrNavig->GetGlobalToLocalTransform().TransformPoint(partLocOld);
-
- G4ThreeVector pVec(pV[0],pV[1],pV[2]);
- const G4double physStep=G4double(propStep*10.);
-
- G4int LttcFlagGeant = ptrGeoInit->GetLttcFlagGeant();
-
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout.precision(10);
- G4cout << "Position (cm):" << pSx << "," << pSy << "," << pSz << G4endl;
- G4cout << "Direction: " << pVec << G4endl;
- G4cout << "Proposed step :"<< propStep << G4endl;
-#endif
-
-
-
- ///////////////////// FLUKA/G4 REGIONS COMPARISON //////////////////////
- //get oldReg pointer and G4volume of previous step end-point pointer
- G4VPhysicalVolume * ptrOldReg = (*pVolStore)[oldReg-1];
- G4VPhysicalVolume * lastVolume = ptrTouchHist->GetVolume();
- G4int touVolCpNb = 0;
- if (lastVolume)
- touVolCpNb = lastVolume->GetCopyNo();
-
-#ifdef G4GEOMETRY_DEBUG
- G4int fluVolCpNb = ptrOldReg->GetCopyNo();
- G4cout << "Fluka volume before step: " << ptrOldReg->GetName()
- << "," << fluVolCpNb<< G4endl;
- if(lastVolume)
- G4cout << "G4 Touch Hist volume: "
- << lastVolume->GetName() << "," << touVolCpNb << G4endl;
- G4cout << "------------------------------------------------" << G4endl;
-#endif
-
-
- //if volume is changed, this is a new particle tracking, or fluka tries
- //to reach a boundary more softly with the same particle. In this case
- //fluka restart tracking from old history, in general. For tracking in
- //lattice-volume fluka goes back to one of the previous lattice volumes.
- //Check if ptrOldReg is equal to old history, or to one of the N lattice
- //volumes stored in jrLt. Then reinitialise step-histories! Otherwise
- //must relocate.
- //NB. jrLtGeant stores lattice volume histories until LttcFlag==-1,
- //then all histories are checked and deleted, and array is reinitialised
-
-
- G4int haveHistNb = -1;
- G4int newRegErr=0;
- G4int indHist = LttcFlagGeant;
- G4int * jrLtGeant = ptrGeoInit->GetJrLtGeantArray();
-
-
- G4NavigationHistory* ptrLttcHist=0;
- if(oldLttc)
- ptrLttcHist = reinterpret_cast<NavHistWithCount*>(oldLttc)->GetNavHistPtr();
-
-
- while(indHist>=0 && haveHistNb==-1) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout<<"Searching in jrLtc...."<<G4endl;
-#endif
- if(oldLttc==jrLtGeant[indHist])
- haveHistNb=indHist;
- indHist-=1;
- }
-
- if(haveHistNb!=-1) {
- //fluka history found in jrLtGeant
- if(haveHistNb<LttcFlagGeant) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout<<"* Fluka reaches boundary more softly..."<<G4endl;
- G4cout<<"* Re-initializing FluggNavigator history"<<G4endl;
- G4cout<<"* and updating step-histories"<<G4endl;
-#endif
-
- ptrNavig->UpdateNavigatorHistory(ptrLttcHist);
- ptrGeoInit->UpdateHistories(ptrLttcHist,0);
- }
-#ifdef G4GEOMETRY_DEBUG
- if(haveHistNb==LttcFlagGeant)
- G4cout << "Continuing step...." << G4endl;
-#endif
- jrLt[0]=oldLttc;
- }
- else {
- //not found fluka history in jrLttGeant!
- G4cout << "* ERROR! Geometry not correctly initialised in fluka history!"
- << G4endl;
-
- //relocation!
- ptrNavig->LocateGlobalPointAndUpdateTouchable(partLoc,0,ptrTouchHist,true);
-
- G4cout << "* ATTENTION: point relocation in: "
- << ptrTouchHist->GetVolume()->GetName() << G4endl;
-
- ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),1);
-
- if(ptrTouchHist->GetVolume() != ptrOldReg) {
- G4cout << "* ERROR! Point not in fluka volume!" << G4endl;
- newRegErr=-3;
- }
-
- //save new history in jrLt[0] and increment its counter
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* ISVHWR call to store new NavHistWithCount in jrLt[0]"
- << G4endl;
-#endif
- jrLt[0]=isvhwr(0,0);
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* CONHWR call to increment counter" << G4endl;
-#endif
- G4int incrCount2=1;
- conhwr(jrLt[0],&incrCount2);
- }
-
-
- //jrLtGeant - history check: decrement counter and delete
- //histories, if LttcFlag=-1, then reinitialise array with -1.
- if(LttcFlag==-1 && LttcFlagGeant>=0) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* CONHWR call to check and delete histories in jrLtGeant[]:"
- << G4endl;
-#endif
-
- for(G4int ind=0; ind<=LttcFlagGeant; ind++) {
- G4int incrCount1=-1;
- if(jrLtGeant[ind]!=jrLt[0]) conhwr(jrLtGeant[ind],&incrCount1);
- jrLtGeant[ind]=-1;
- }
- LttcFlagGeant=-1;
- }
-
-
- //update jrLt and sLt arrays
- G4int end = 99;
- if(LttcFlag>=0)
- end=LttcFlag;
-
- // Added by A.Solodkov
- if (end>=100) {
- G4cout << "Problems in WrapG1 routine" << G4endl;
- G4cout << "Index LttcFlag=" << end << " is outside array bounds" << G4endl;
- G4cout << "Better to stop immediately !" << G4endl;
- exit(1);
- }
- //jrLt re-initialization with -1 (jrLt[0] is already set)
- for(G4int vv=1;vv<=end;vv++)
- jrLt[vv]=-1;
- //sLt re-initialization
- for(G4int vs=0;vs<=end;vs++)
- sLt[vs]=0;
-
- LttcFlag=0;
-
-
- ///////////////////// COMPUTE STEP ////////////////////////
- //update Navigator private flags and voxel stack if point is
- //changed from last located point (otherwise troubles come
- //when fluka changes location or particle because G4 computes
- //from last located point).
-
-// change because of geant4 6.0
-//G4ThreeVector newLocalPoint = ptrNavig->ComputeLocalPoint(partLoc);
- G4ThreeVector newLocalPoint =
- ptrNavig->GetGlobalToLocalTransform().TransformPoint(partLoc);
-
-
-
- G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2();
- if(moveLenSq>=kCarTolerance*kCarTolerance)
- ptrNavig->LocateGlobalPointWithinVolume(partLoc);
-
-
- //compute step and new location
- newReg = oldReg;
- G4double appStep = 0;
- G4double safety = 0;
- G4bool onBoundaryRound = false;
- G4bool crossBound = false;
- G4double physStepTmp = physStep;
- G4bool flagError = false;
-
- while( (oldReg==newReg && appStep<physStepTmp) || onBoundaryRound ) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* Computing step..." << G4endl;
-#endif
-
- //update variables
- oldReg = newReg;
-
- if(onBoundaryRound) {
- physStepTmp=10.e-10;
- //compute step and location: returns newReg
- newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
- safety,onBoundaryRound,flagError,oldReg);
- physStepTmp=0.;
- crossBound=true;
- }
- else {
- physStepTmp -= appStep;
- //compute step and location: returns newReg
- newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
- safety,onBoundaryRound,flagError,oldReg);
- }
-
- G4bool EqualHist;
- EqualHist = EqualHistories(ptrTouchHist->GetHistory(),
- ptrGeoInit->GetTempNavHist()->GetHistory());
-
- if(!EqualHist && flagError) {
- pVec=-pVec;
- newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
- safety,onBoundaryRound,flagError,oldReg);
- pVec=-pVec;
- physStepTmp+=1;
- newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
- safety,onBoundaryRound,flagError,oldReg);
-
- EqualHist = EqualHistories(ptrTouchHist->GetHistory(),
- ptrGeoInit->GetTempNavHist()->GetHistory());
- }
-
- //update sLt
- G4double pas = G4double(appStep);
- sLt[LttcFlag] += pas/cm;
- safety = (oldReg!=newReg)?0.:safety;
-
- if(!EqualHist) {
- //end-step point is on boundary between volumes;
- //==> update step-histories, save new NavHistWithCount
- //and save its pointer in jrLt; save step in sLt.
-
- //set onBoundaryRound=false to avoid re-compute step!
- onBoundaryRound=false;
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* History is changed!" << G4endl;
- G4cout << "* updating step-histories, jrLt, LttcFlag" << G4endl;
-#endif
-
- ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
-
- LttcFlag += 1;
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* ISVHWR call to store new NavHistWithCount in jrLt"
- << G4endl;
-#endif
-
- jrLt[LttcFlag] = isvhwr(0,0);
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* CONHWR call to increment counter" << G4endl;
-#endif
- G4int incrCount3=1;
- conhwr(jrLt[LttcFlag],&incrCount3);
-
- sLt[LttcFlag] = sLt[LttcFlag-1];
- }
- }
-
- ////////////////////// OUTPUT //////////////////////////
- //If back-scattering occured, and fluka is in the wrong region, return -3.
- //(N. B. Boundary between replicans are not seen when physStep=distance
- //particle-boundary: in this case step=kInfinity and history is unchanged.
- //Following step is =0 and then history changes.)
- if(nascFlag<0 && !appStep && physStep && newReg!=oldReg && !crossBound) {
- //don't need to compare histories because boundary between
- //identical volumes in different replicans are not seen
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* Back-scattering!" << G4endl;
-#endif
-
- newReg=-3;
- }
- else if(newRegErr<0)
- newReg=newRegErr;
-//======================= Endre print ===================
-printf("\nWrapG1 mreg=%d",newReg);
-//======================= Endre print ===================
-
- //compute output variables (in cm.!)
- //final step
- retStep = sLt[LttcFlag];
-
- //safety (Fluka sottracts a bit to safety to be sure
- //not to jump on a boundary)
- G4double s = G4double(safety);
- s -= s*3.0e-09;
- saf=s/cm;
-
- //update wrapper utility variables
- //copy jrLt in jrLtGeant
- G4int start=0;
- if(haveHistNb!=-1 && LttcFlagGeant!=-1)
- start=1;
- for(G4int lt=start;lt<=LttcFlag;lt++)
- jrLtGeant[LttcFlagGeant+1+lt-start]=jrLt[lt];
- LttcFlagGeant+=(1+LttcFlag-start);
- newLttc = jrLt[LttcFlag];
- ptrGeoInit->SetLttcFlagGeant(LttcFlagGeant);
-
- partLocOld=partLoc;
-
-// change because of geant4 6.0
-//oldLocalPoint = ptrNavig->ComputeLocalPoint(partLocOld);
- oldLocalPoint =
- ptrNavig->GetGlobalToLocalTransform().TransformPoint(partLocOld);
-
- //compute new position
- G4ThreeVector oldPos = G4ThreeVector(pSx,pSy,pSz);
- G4ThreeVector newPos = oldPos + retStep*pVec;
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "New position: " << newPos << G4endl;
- G4cout << "Output region: " << newReg << G4endl;
- G4cout << "G4 safety (cm): " << (safety*0.1) << G4endl;
- G4cout << "Fluka safety (cm): " << saf << G4endl;
- G4cout << "Approved step: " << retStep << G4endl;
- G4cout << "LttcFlag = " << LttcFlag << G4endl;
- for(G4int i=0;i<=LttcFlag+1;i++) {
- G4cout << "jrLt[" << i <<"]=" << jrLt[i] << G4endl;
- G4cout << "sLt[" << i <<"]=" << sLt[i] << G4endl;
- }
-
- G4cout << "LttcFlagGeant =" << LttcFlagGeant << G4endl;
- for(G4int ib=0;ib<=LttcFlagGeant+1;ib++) {
- G4cout << "jrLtGeant[" << ib <<"]=" << jrLtGeant[ib] << G4endl;
- }
- G4cout << "newLttc=" << newLttc << G4endl;
-#endif
-}
-
-
-
-
+++ /dev/null
-
-// Flugg tag
-
-///////////////////////////////////////////////////////////////////
-//
-// WrapG1.hh - Sara Vanini
-//
-// Dummy wrapper for geometry tracking.
-//
-// modified 24.10.01: by I. Hrivnacova
-// functions declarations separated from implementation
-// (moved to Wrappers.hh);
-//
-///////////////////////////////////////////////////////////////////
-
-
-#include "Wrappers.hh"
-#include "globals.hh"
-
-void g1rtwr(void)
-{
-//flag
-#ifdef G4GEOMETRY_DEBUG
- G4cout<<"============ G1RTWR ============="<<G4endl;
-#endif
-
-//dummy wrapper
-}
+++ /dev/null
-
-// Flugg tag
-
-///////////////////////////////////////////////////////////////////
-//
-// WrapIncrHist.hh - Sara Vanini
-//
-// Wrapper for updating secondary particles history counter.
-// If counter=0 the history is deleted.
-//
-// modified 14/I/9999
-// modified 24.10.01: by I. Hrivnacova
-// functions declarations separated from implementation
-// (moved to Wrappers.hh);
-//
-//////////////////////////////////////////////////////////////////
-
-
-#include "Wrappers.hh"
-#include "FGeometryInit.hh"
-#include "NavHistWithCount.hh"
-#include "globals.hh"
-
-
-void conhwr(G4int& intHist, G4int* incrCount)
-{
-//flag
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "============= CONHWR ==============" << G4endl;
- G4cout << "Ptr History = " << intHist << G4endl;
-#endif
-
- //get NavHistWithCount pointer
- if(intHist!=-1) {
- NavHistWithCount* ptrNavHistCount=
- reinterpret_cast<NavHistWithCount*>(intHist);
-
- //for debugging...
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "Secondary counter=" << ptrNavHistCount->GetCount();
- if(*incrCount>0)
- G4cout << "+" << *incrCount << G4endl;
- if(*incrCount<0)
- G4cout<< *incrCount << G4endl;
- if(*incrCount==0)
- G4cout << G4endl;
-#endif
-
- //update secondary particles counter
- ptrNavHistCount->UpdateCount(*incrCount);
-
- //delete history if counter=0 or if counter=-1
- G4int counter = ptrNavHistCount->GetCount();
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "Counter = " << counter << G4endl;
-#endif
- if(!counter || counter==-1) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "Deleting Nav Hist object..." << G4endl;
-#endif
- /*
- //for history checking....
- G4int index = ptrNavHistCount->GetCheckInd();
- static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
- G4int * ptrArray = ptrGeoInit->GetHistArray();
- ptrArray[index]=0;
- */
-
- delete ptrNavHistCount;
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "end delete" << G4endl;
-#endif
- intHist=-1;
- }
- }
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "============= Out of CONHWR ==============" << G4endl;
-#endif
- }
-
-
-
-
-
-
-
-
-
-
+++ /dev/null
-
-// Flugg tag
-
-///////////////////////////////////////////////////////////////////
-//
-// WrapIniHist.hh - Sara Vanini
-//
-// Wrapper for reinitialization of FluggNavigator history.
-//
-// modified 14/I/99
-// modified 24.10.01: by I. Hrivnacova
-// functions declarations separated from implementation
-// (moved to Wrappers.hh);
-//
-///////////////////////////////////////////////////////////////////
-
-
-#include "Wrappers.hh"
-#include "FGeometryInit.hh"
-#include "NavHistWithCount.hh"
-#include "G4NavigationHistory.hh"
-#include "FluggNavigator.hh"
-#include "globals.hh"
-
-void inihwr(G4int& intHist)
-{
-//flag
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "============= INIHWR ==============" << G4endl;
- G4cout << "Ptr History=" <<intHist<< G4endl;
-#endif
- if(intHist==-1) {
- G4cout << "ERROR! This history has been deleted!" << G4endl;
- return;
- }
- else {
- //get NavHistWithCount,G4NavigationHistory,FGeometryInit,
- //FluggNavigator pointers
- NavHistWithCount* ptrNavHistCount=
- reinterpret_cast<NavHistWithCount*>(intHist);
- G4NavigationHistory* ptrNavHist=ptrNavHistCount->GetNavHistPtr();
- static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
- FluggNavigator* ptrNavig = ptrGeoInit->getNavigatorForTracking();
-
- //reinitialize navigator history
- ptrNavig->UpdateNavigatorHistory(ptrNavHist);
-
- //update utility histories: touch,temp, and reset old history
- ptrGeoInit->UpdateHistories(ptrNavHist,0);
-
- //save new history in jrLtGeant if not present
- G4int LttcFlagGeant = ptrGeoInit->GetLttcFlagGeant();
- G4int * jrLtGeant = ptrGeoInit->GetJrLtGeantArray();
-
- G4bool intHistInJrLtGeant = false;
- for(G4int h=0; h<=LttcFlagGeant; h++)
- if(jrLtGeant[h]==intHist)
- intHistInJrLtGeant = true;
-
- if(!intHistInJrLtGeant) {
- LttcFlagGeant += 1;
- ptrGeoInit->SetLttcFlagGeant(LttcFlagGeant);
-
- jrLtGeant[LttcFlagGeant]=intHist;
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* CONHWR call to increment counter" << G4endl;
-#endif
- G4int incrCount=1;
- conhwr(jrLtGeant[LttcFlagGeant],&incrCount);
- }
-
- //print history....
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "History reinitialized in:" << G4endl;
- G4cout << *ptrNavHist << G4endl;
- ptrGeoInit->PrintJrLtGeant();
-#endif
- }
-}
-
-
-
-
-
+++ /dev/null
-
-// Flugg tag
-
-///////////////////////////////////////////////////////////////////
-//
-// WrapInit.hh - Sara Vanini
-//
-// Wrapper for geometry initialisation.
-//
-// modified 12-IV-2000
-// modified 24.10.01: by I. Hrivnacova
-// functions declarations separated from implementation
-// (moved to Wrappers.hh);
-//
-//////////////////////////////////////////////////////////////////
-
-
-#include "Wrappers.hh"
-#include "FGeometryInit.hh"
-#include "G4PhysicalVolumeStore.hh"
-#include "G4VPhysicalVolume.hh"
-#include "globals.hh"
-
-void jomiwr(const G4int & nge, const G4int& lin, const G4int& lou,
- G4int& flukaReg)
-{
-//flag
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "================== JOMIWR =================" << G4endl;
-#endif
-
-
- //Geoinit Pointer
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "\t *==> JOMIWR: Getting FGeometry..." << G4endl;
-#endif
- FGeometryInit * ptrGeoInit=FGeometryInit::GetInstance();
-
- //initialize geometry:construct detector and set world volume
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "\t *==> JOMIWR: Setting the detector..." << G4endl;
-#endif
-// ptrGeoInit->setDetector();
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "\t *==> JOMIWR: Setting mother volume..." << G4endl;
-#endif
-// ptrGeoInit->setMotherVolume();
-
- //close geometry for optimization
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "\t *==> JOMIWR: Closing geometry..." << G4endl;
-#endif
-// ptrGeoInit->closeGeometry();
-
- //initialize wrappers utility histories at the beginning of run and set flag
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "\t *==> JOMIWR: InitHistories..." << G4endl;
-#endif
-// ptrGeoInit->InitHistories();
-
- //initialize lattice array
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "\t *==> JOMIWR: Init lattice array..." << G4endl;
-#endif
-// ptrGeoInit->InitJrLtGeantArray();
-
- //initialize debug-array
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "\t *==> JOMIWR: Init debug array..." << G4endl;
-#endif
-// ptrGeoInit->InitHistArray();
-
- //create Fluka material cards in flukaMat.inp file
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "\t *==> JOMIWR: Init fluka materials..." << G4endl;
-#endif
-// ptrGeoInit->createFlukaMatFile();
-
- //returns number of volumes + 1
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "\t *==> JOMIWR: Returning..." << G4endl;
-#endif
- G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
- G4int numVol = G4int(pVolStore->size());
- flukaReg = numVol + 1;
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "Number of volumes + 1: " << flukaReg << G4endl;
- G4cout << "================== Out of JOMIWR =================" << G4endl;
-#endif
-}
-
-
-
+++ /dev/null
-
-// Flugg tag
-
-///////////////////////////////////////////////////////////////////
-//
-// WrapG1.hh - Sara Vanini
-//
-// Dummy wrapper (in fluka: for geometry debugging)
-//
-//////////////////////////////////////////////////////////////////
-
-//#ifndef lkdbwr
-//#define lkdbwr lkdbwr_
-
-#include "Wrappers.hh"
-#include "globals.hh"
-
-void lkdbwr(G4double& pSx, G4double& pSy, G4double& pSz,
- G4double* pV, const G4int& oldReg, const G4int& oldLttc,
- G4int& newReg, G4int& flagErr, G4int& newLttc)
-{
- //flag
-#ifdef G4GEOMETRY_DEBUG
- G4cout<<"============= LKDBWR =============="<<G4endl;
-#endif
-
- //return region number and dummy variables
- newReg=0;
- newLttc=0;
- flagErr=-1;
-}
-//#endif
+++ /dev/null
-
-// Flugg tag
-
-///////////////////////////////////////////////////////////////////
-//
-// WrapLookFX.hh - Sara Vanini - 24/III/00
-//
-// Wrapper for localisation of particle to fix particular conditions.
-// At the moment is the same as WrapLookZ.hh.
-//
-// modified 24.10.01: by I. Hrivnacova
-// functions declarations separated from implementation
-// (moved to Wrappers.hh);
-// modified 17/06/02: by I. Gonzalez. STL migration
-//
-//////////////////////////////////////////////////////////////////
-
-
-#include "Wrappers.hh"
-#include "FGeometryInit.hh"
-#include "G4VPhysicalVolume.hh"
-#include "FluggNavigator.hh"
-#include "G4ThreeVector.hh"
-#include "G4PhysicalVolumeStore.hh"
-#include "globals.hh"
-
-void lkfxwr(G4double& pSx, G4double& pSy, G4double& pSz,
- G4double* pV, const G4int& oldReg, const G4int& oldLttc,
- G4int& newReg, G4int& flagErr, G4int& newLttc)
-{
- //flag
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "======= LKFXWR =======" << G4endl;
-#endif
-
- //FGeometryInit, navigator, volumeStore pointers
- static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
- FluggNavigator * ptrNavig = ptrGeoInit->getNavigatorForTracking();
- G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
-
- //coordinates in mm.
- G4ThreeVector pSource(pSx,pSy,pSz);
- pSource *= 10.0; //in millimeters!
-
- //locate point and update histories
- G4TouchableHistory * ptrTouchableHistory =
- ptrGeoInit->GetTouchableHistory();
- ptrNavig->LocateGlobalPointAndUpdateTouchable(pSource,0,
- ptrTouchableHistory,true);
- //updating tmp but not old histories, they are useful in
- //case of RGRPWR call, or when fluka, after a LOOKZ call,
- //descards step for multiple scattering and returns to old history
- //NO, after lattice-fix we don't need old history anymore!
-
- ptrGeoInit->UpdateHistories(ptrTouchableHistory->GetHistory(),0);
- G4VPhysicalVolume * located = ptrTouchableHistory->GetVolume();
-
- //if volume not found, out of mother volume: returns "number of volumes"+1
- if(!located) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "Out of mother volume!";
-#endif
- G4int numVol = G4int(pVolStore->size());
- newReg = numVol + 1;
- }
- else {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* ISVHWR call to store current NavHistWithCount in jrLtGeant"
- << G4endl;
-#endif
-
- //save history in jrLtGeant and increment counter
- G4int * jrLtGeant = ptrGeoInit->GetJrLtGeantArray();
- G4int LttcFlagGeant = ptrGeoInit->GetLttcFlagGeant();
- LttcFlagGeant += 1;
- jrLtGeant[LttcFlagGeant] = isvhwr(0,0);
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* CONHWR call to increment counter" << G4endl;
-#endif
- G4int incrCount=1;
- conhwr(jrLtGeant[LttcFlagGeant],&incrCount);
-
- //update LttcFlagGeant
- ptrGeoInit->SetLttcFlagGeant(LttcFlagGeant);
-
- //return region number and dummy variables
- G4int volIndex = ~0;
- for (unsigned int i=0; i<pVolStore->size(); i++)
- if ((*pVolStore)[i] == located)
- volIndex = i;
- // G4int volIndex=G4int(pVolStore->index(located));
- if (volIndex==(~0)) {
- G4cerr << "FLUGG: Problem in file WrapLookFX trying to find volume after step" << G4endl;
- exit(-999);
- }
- //G4int volIndex=G4int(pVolStore->index(located));
- newReg=volIndex+1;
- newLttc=jrLtGeant[LttcFlagGeant];
- flagErr=newReg;
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "LKFXWR Located Physical volume = ";
- G4cout << located->GetName() << G4endl;
-#endif
- }
-}
-
-
-
-
+++ /dev/null
-
-// Flugg tag
-
-///////////////////////////////////////////////////////////////////
-//
-// WrapLookMG.hh - Sara Vanini 26/X/99
-//
-// Wrapper for localisation of particle for magnetic field tracking
-//
-// modified 13/IV/00: check if point belongs to one of the lattice
-// histories stored in jrLtGeant
-//
-// modified 24.10.01: by I. Hrivnacova
-// functions declarations separated from implementation
-// (moved to Wrappers.hh);
-// modified 17/06/02: by I. Gonzalez. STL migration
-//
-///////////////////////////////////////////////////////////////////
-
-
-#include "Wrappers.hh"
-#include "FGeometryInit.hh"
-#include "NavHistWithCount.hh"
-#include "G4PhysicalVolumeStore.hh"
-#include "G4NormalNavigation.hh"
-#include "G4VoxelNavigation.hh"
-#include "G4ParameterisedNavigation.hh"
-#include "G4ReplicaNavigation.hh"
-#include "globals.hh"
-
-
-//auxiliary function declarations
-G4bool PointLocate(const G4NavigationHistory &,
- const G4ThreeVector &,
- const G4ThreeVector *,
- G4int &);
-EVolume CharacteriseDaughters(const G4LogicalVolume *);
-
-
-void lkmgwr(G4double& pSx, G4double& pSy, G4double& pSz,
- G4double* pV, const G4int& oldReg, const G4int& oldLttc,
- G4int& flagErr, G4int& newReg, G4int& newLttc)
-{
-//flag
-#ifdef G4GEOMETRY_DEBUG
- G4cout<<"============= LKMGWR =============="<<G4endl;
-#endif
-
- //Geoinit, Navigator, etc. pointers
- static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
- G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
-
- //setting variables (and dimension: Fluka uses cm.!)
- G4ThreeVector globalPointcm(pSx,pSy,pSz);
- G4ThreeVector globalPoint = globalPointcm * 10.; //in mm.
- G4ThreeVector globalDirection(pV[0],pV[1],pV[2]);
-
- //get jrLtGeant array and initialize variables
- G4int * jrLtGeant = ptrGeoInit->GetJrLtGeantArray();
- G4int LttcFlagGeant = ptrGeoInit->GetLttcFlagGeant();
- G4bool belongToVolume = false;
- G4int i = LttcFlagGeant;
- G4int regionIn = 0;
- newReg = -1;
- newLttc = -1;
-
- while(!belongToVolume && i>=0) {
- //get history from jrLtGeant
- NavHistWithCount * ptrNavHistCount =
- reinterpret_cast<NavHistWithCount*>(jrLtGeant[i]);
- const G4NavigationHistory * ptrNavHist =
- ptrNavHistCount->GetNavHistPtr();
-
- //check if globalPoint belongs to volume (can't call
- //LocateGlobalPoint... because of flag settings)
- belongToVolume = PointLocate(*ptrNavHist,globalPoint,
- &globalDirection,regionIn);
-
-
- //if point belongs to surface, check scalar product direction*normal
- if(regionIn==-100) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout<<"On surface!"<<G4endl;
-#endif
-
- //if entering, then point belongs to volume,
- //if exiting, point doesn't belong to volume.
- G4int oldReg,flagErr,reg;
- G4double x,y,z,px,py,pz;
- G4double * norml = new G4double[3];
- x = globalPoint.x();
- y = globalPoint.y();
- z = globalPoint.z();
- px = globalDirection.x();
- py = globalDirection.y();
- pz = globalDirection.z();
-
- nrmlwr(x,y,z,px,py,pz,norml,oldReg,reg,flagErr);
-
- G4ThreeVector normal(norml[0],norml[1],norml[2]);
-#ifdef G4GEOMETRY_DEBUG
- // G4cout<<"Scalar product="<<globalDirection.dot(normal)<<G4endl;
-#endif
- if(globalDirection.dot(normal)>0) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout<<"entering volume!"<<G4endl;
-#endif
- G4int volIndex = ~0;
- for (unsigned int i=0; i<pVolStore->size(); i++)
- if ((*pVolStore)[i] == ptrNavHist->GetTopVolume())
- volIndex = i;
- if (volIndex==(~0)) {
- G4cerr << "FLUGG: Problem in routine WrapG1 tryingto find volume after step" << G4endl;
- exit(-999);
- }
- //regionIn = G4int(pVolStore->index(ptrNavHist->GetTopVolume()))+1;
- regionIn = volIndex+1;
- belongToVolume = true;
- }
- }
-
- i -= 1;
- }
-
- //output variables
- if(belongToVolume) {
- newReg = regionIn;
- newLttc = jrLtGeant[i+1];
- }
-
- flagErr = pVolStore->size() + 2; //flagErr=fluka region number + 1
-
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "Global point (cm): " << globalPointcm << G4endl;
- G4cout << "Direction: " << globalDirection << G4endl;
- if(newReg!=-1)
- G4cout << "Point belongs to region " << newReg
- << " ptr history=" << newLttc << G4endl;
- else
- G4cout << "No containing volume found!" << G4endl;
-#endif
-}
-
-
-
-
-//returns false if the point doesn't belong to history top volume (can either
-//belong to one of its daughter volumes or none), otherwise returns true.
-
-G4bool PointLocate(const G4NavigationHistory & blockedHistConst,
- const G4ThreeVector & globalPoint,
- const G4ThreeVector * pGlobalDirection,
- G4int & reg)
-{
-
- //variables
- G4NavigationHistory * fHistory = new G4NavigationHistory(blockedHistConst);
- G4bool belongVolume;
-
- //G4 flags resetted (see: ResetStackAndState)
- G4bool fExiting = false;
- G4VPhysicalVolume * fBlockedPhysicalVolume=0;
- G4int fBlockedReplicaNo=-1;
- G4bool fLocatedOnEdge=false;
-
- G4bool notKnownContained=true,noResult;
- G4VPhysicalVolume *targetPhysical;
- G4LogicalVolume *targetLogical;
- G4VSolid *targetSolid;
- G4ThreeVector localPoint,localDirection;
- EInside insideCode;
-
-
- //Helpers/Utility classes
- G4NormalNavigation fnormalNav;
- G4VoxelNavigation fvoxelNav;
- G4ParameterisedNavigation fparamNav;
- G4ReplicaNavigation freplicaNav;
- G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
-
- //local variables
- localPoint=fHistory->GetTopTransform().TransformPoint(globalPoint);
- localDirection=fHistory->GetTopTransform().TransformPoint(*pGlobalDirection);
-
- //search if volume contains point
- if (fHistory->GetTopVolumeType()!=kReplica) {
- targetSolid=fHistory->GetTopVolume()->GetLogicalVolume()->GetSolid();
- insideCode=targetSolid->Inside(localPoint);
- }
- else {
- insideCode=freplicaNav.BackLocate(*fHistory,globalPoint,
- localPoint,fExiting,notKnownContained);
- // !CARE! if notKnownContained returns false then the point is within
- // the containing placement volume of the replica(s). If insidecode
- // will result in the history being backed up one level, then the
- // local point returned is the point in the system of this new level
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "Replica: fExiting=" << fExiting << G4endl;
- G4cout << " notKnownContained=" << notKnownContained << G4endl;
-#endif
- }
-
- if (insideCode==kOutside)
- belongVolume = false;
-
- else if (insideCode==kSurface) {
- belongVolume = false;
- reg = -100;
- }
- else {
- // Search downwards in daughter volumes
- // until deepest containing volume found
- //
- // 3 Cases:
- //
- // o Parameterised daughters
- // =>Must be one G4PVParameterised daughter & voxels
- // o Positioned daughters & voxels
- // o Positioned daughters & no voxels
-
- belongVolume = true;
- noResult=true;
-
- do {
- // Determine `type' of current mother volume
- targetPhysical=fHistory->GetTopVolume();
- targetLogical=targetPhysical->GetLogicalVolume();
- switch(CharacteriseDaughters(targetLogical)) {
- case kNormal:
- if (targetLogical->GetVoxelHeader()) {
- noResult=fvoxelNav.LevelLocate(*fHistory,
- fBlockedPhysicalVolume,
- fBlockedReplicaNo,
- globalPoint,
- pGlobalDirection,
- fLocatedOnEdge,
- localPoint);
- }
- else {
- noResult=fnormalNav.LevelLocate(*fHistory,
- fBlockedPhysicalVolume,
- fBlockedReplicaNo,
- globalPoint,
- pGlobalDirection,
- fLocatedOnEdge,
- localPoint);
- }
- break;
- case kReplica:
- noResult=freplicaNav.LevelLocate(*fHistory,
- fBlockedPhysicalVolume,
- fBlockedReplicaNo,
- globalPoint,
- pGlobalDirection,
- fLocatedOnEdge,
- localPoint);
- break;
- case kParameterised:
- noResult=fparamNav.LevelLocate(*fHistory,
- fBlockedPhysicalVolume,
- fBlockedReplicaNo,
- globalPoint,
- pGlobalDirection,
- fLocatedOnEdge,
- localPoint);
- break;
- } //switch
-
- // LevelLocate search in the first daughter level.
- // LevelLocate returns noResult=true if it finds a daughter volume
- // in which globalPoint is inside (or on the surface). So point
- // doesn't belong only to mother volume ==> belongVolume=false.
-
- if (noResult) {
- // The blocked volume no longer valid - it was for another level
- fBlockedPhysicalVolume= 0;
- fBlockedReplicaNo= -1;
- belongVolume=false;
- }
- } while (noResult);
-
- //on exit targetPhysical is the volume globalPoint belongs to;
- G4int volIndex = ~0;
- for (unsigned int i=0; i<pVolStore->size(); i++)
- if ((*pVolStore)[i] == targetPhysical)
- volIndex = i;
- if (volIndex==(~0)) {
- G4cerr << "FLUGG: Problem in routine WrapG1 tryingto find volume after step" << G4endl;
- exit(-999);
- }
- //reg = G4int(pVolStore->index(targetPhysical))+1;
- reg = volIndex+1;
-
- }
-
- delete fHistory;
- return belongVolume;
-}
-
-EVolume CharacteriseDaughters(const G4LogicalVolume *pLog)
-{
- EVolume type;
- EAxis axis;
- G4int nReplicas;
- G4double width,offset;
- G4bool consuming;
- G4VPhysicalVolume *pVol;
-
- if (pLog->GetNoDaughters()==1) {
- pVol=pLog->GetDaughter(0);
- if (pVol->IsReplicated()) {
- pVol->GetReplicationData(axis,nReplicas,width,offset,consuming);
- type=(consuming) ? kReplica : kParameterised;
- }
- else
- type=kNormal;
- }
- else
- type=kNormal;
- return type;
-}
-
-
-
+++ /dev/null
-
-// Flugg tag
-
-///////////////////////////////////////////////////////////////////
-//
-// WrapLookZ.hh - Sara Vanini
-//
-// Wrapper for localisation of starting point of particle.
-//
-// modified 20/III/00: history initialization moved to ISVHWR
-// modified 13/IV/00: located history saved in jrLtcGeant and
-// incremented
-// modified 24.10.01: by I. Hrivnacova
-// functions declarations separated from implementation
-// (moved to Wrappers.hh);
-// modified 17/06/02: by I. Gonzalez. STL migration
-//
-//////////////////////////////////////////////////////////////////
-
-
-#include "Wrappers.hh"
-#include "FGeometryInit.hh"
-#include "G4VPhysicalVolume.hh"
-#include "FluggNavigator.hh"
-#include "G4ThreeVector.hh"
-#include "G4PhysicalVolumeStore.hh"
-#include "globals.hh"
-
-
-void lkwr(G4double& pSx, G4double& pSy, G4double& pSz,
- G4double* pV, const G4int& oldReg, const G4int& oldLttc,
- G4int& newReg, G4int& flagErr, G4int& newLttc)
-{
- //flag
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "======= LKWR =======" << G4endl;
-#endif
-
- //FGeometryInit, navigator, volumeStore pointers
- static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
- FluggNavigator * ptrNavig = ptrGeoInit->getNavigatorForTracking();
- G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
-
- //coordinates in mm.
- G4ThreeVector pSource(pSx,pSy,pSz);
- pSource *= 10.0; //in millimeters!
-
- //locate point and update histories
- G4TouchableHistory * ptrTouchableHistory =
- ptrGeoInit->GetTouchableHistory();
- ptrNavig->LocateGlobalPointAndUpdateTouchable(pSource,0,
- ptrTouchableHistory,true);
- //updating tmp but not old histories, they are useful in
- //case of RGRPWR call, or when fluka, after a LOOKZ call,
- //descards step for multiple scattering and returns to old history
- //NO, after lattice-fix we don't need old history anymore!
-
- ptrGeoInit->UpdateHistories(ptrTouchableHistory->GetHistory(),0);
- G4VPhysicalVolume * located = ptrTouchableHistory->GetVolume();
-
- //if volume not found, out of mother volume: returns "number of volumes"+1
- if(!located) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "Out of mother volume!";
-#endif
- G4int numVol = G4int(pVolStore->size());
- newReg = numVol + 1;
- }
- else {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* ISVHWR call to store current NavHistWithCount in jrLtGeant"
- << G4endl;
-#endif
-
- //save history in jrLtGeant and increment counter
- G4int * jrLtGeant = ptrGeoInit->GetJrLtGeantArray();
- G4int LttcFlagGeant = ptrGeoInit->GetLttcFlagGeant();
- LttcFlagGeant += 1;
- jrLtGeant[LttcFlagGeant] = isvhwr(0,0);
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* CONHWR call to increment counter" << G4endl;
-#endif
- G4int incrCount=1;
- conhwr(jrLtGeant[LttcFlagGeant],&incrCount);
-
- //update LttcFlagGeant
- ptrGeoInit->SetLttcFlagGeant(LttcFlagGeant);
-
- //return region number and dummy variables
- G4int volIndex = ~0;
- for (unsigned int i=0; i<pVolStore->size(); i++)
- if ((*pVolStore)[i] == located) volIndex = i;
- //G4int volIndex=G4int(pVolStore->index(located));
- if (volIndex==(~0)) {
- G4cerr << "FLUGG: Problem in routine WrapG1 tryingto find volume after step" << G4endl;
- exit(-999);
- }
- newReg=volIndex+1;
- newLttc=jrLtGeant[LttcFlagGeant];
- flagErr=newReg;
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "LKWR Located Physical volume = ";
- G4cout << located->GetName() << G4endl;
-#endif
- }
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "======= Out of LKWR =======" << G4endl;
-#endif
-}
-
-
-
-
+++ /dev/null
-
-// Flugg tag
-
-///////////////////////////////////////////////////////////////////
-//
-// WrapMag.hh - Sara Vanini
-//
-// Wrapper for geometry tracking in magnetic field: returns magnetic
-// field values in a given position.
-//
-// modified 26/X/1998
-// modified 18/XI/1999
-// modified 24.10.01: by I. Hrivnacova
-// functions declarations separated from implementation
-// (moved to Wrappers.hh);
-//
-/////////////////////////////////////////////////////////////////
-
-
-#include "Wrappers.hh"
-#include "FGeometryInit.hh"
-#include "globals.hh"
-
-void magfld(const G4double& pX, const G4double& pY, const G4double& pZ,
- G4double& cosBx, G4double& cosBy, G4double& cosBz,
- G4double& Bmag, G4int& reg, G4int& idiscflag)
-
-{
- //flag
-#ifdef G4GEOMETRY_DEBUG
- G4cout<<"================== MAGFLD ================="<<G4endl;
-#endif
-
- //Geoinit Pointer
- FGeometryInit * ptrGeoInit=FGeometryInit::GetInstance();
-
- //get FieldManager, Field pointers for magnetic field handling
- G4FieldManager * pFieldMgr = ptrGeoInit->getFieldManager();
- const G4Field * ptrField = pFieldMgr->GetDetectorField();
-
- //compute field
- G4double point[3];
- point[0] = pX*10.;
- point[1] = pY*10.;
- point[2] = pZ*10.;
- G4double B[3];
- ptrField->GetFieldValue(point,B);
- G4double Bnor = sqrt(sqr(B[0])+sqr(B[1])+sqr(B[2]));
- if(Bnor) {
- cosBx = B[0]/Bnor;
- cosBy = B[1]/Bnor;
- cosBz = B[2]/Bnor;
- }
- else {
- cosBx = 0;
- cosBy = 0;
- cosBz = 1;
- }
-
- Bmag = Bnor/tesla;
- idiscflag = 0;
-}
-
-
-
-
+++ /dev/null
-
-// Flugg tag
-
-////////////////////////////////////////////////////////////////////
-//
-// WrapNorml.hh - Sara Vanini
-//
-// Wrapper for computing normal unit-vector in global coordinates.
-//
-// Fluka requires normal vector exiting from final position (of
-// particle) volume, that is: entering in volume of initial position.
-// Geant4 always computes the normal vector exiting from the volume.
-// In GetLocalExitNormal() call the volume is the pre-step volume
-// (so G4 normal vector sign is opposite of fluka-required normal).
-// If IsLocalExitNormalValid=false, normal value is computed from
-// init-step volume (in this case sign must change), or from
-// end-step volume (the sign is the same). Normal vector is computed
-// always on boundary between volumes, in global coordinates (to take
-// rotation of parameterised volumes in hierarchy in consideration).
-// So: nrmlwr returns inwards pointing unit normal of the shape for
-// surface closest to the point returned by the navigator (last step
-// end-point).
-//
-// modified 10/III/99
-// modified 25/V/00
-// modified 7/VI/00 for boundary-crossing in case of relocation
-// modified 5/VII/00 geometry error on boundary fixed
-// modified 24.10.01: by I. Hrivnacova
-// functions declarations separated from implementation
-// (moved to Wrappers.hh);
-//
-////////////////////////////////////////////////////////////////////
-
-#include "Wrappers.hh"
-#include "WrapUtils.hh"
-#include "FGeometryInit.hh"
-#include "G4VPhysicalVolume.hh"
-#include "FluggNavigator.hh"
-#include "G4ThreeVector.hh"
-#include "globals.hh"
-
-
-void nrmlwr(G4double& /*pSx*/, G4double& /*pSy*/, G4double& /*pSz*/,
- G4double& pVx, G4double& pVy, G4double& pVz,
- G4double* norml, const G4int& oldReg,
- const G4int& /*newReg*/, G4int& flagErr)
-{
- //flag
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "============ NRMLWR-DBG =============" << G4endl;
-#endif
-
- //dummy variables
- flagErr=0;
-
- //navigator pointer
- static FGeometryInit * ptrGeoInit;
- ptrGeoInit = FGeometryInit::GetInstance();
- FluggNavigator * ptrNavig = ptrGeoInit->getNavigatorForTracking();
-
- //variables
- G4ThreeVector normalLoc;
- G4ThreeVector normalGlob;
-
- //normal computing
- G4bool valid;
- normalLoc=ptrNavig->GetLocalExitNormal(&valid);
- if(valid) {
- normalLoc *= -1;
-
- //global cooordinates normal
- normalGlob = ptrNavig->GetLocalToGlobalTransform().
- TransformAxis(normalLoc);
- }
- else {
- G4VPhysicalVolume *touchvolume;
- G4ThreeVector theLocalPoint;
-
- if(ptrNavig->EnteredDaughterVolume()) {
- //volume from touchable history
- G4TouchableHistory *ptrTouchableHistory=ptrGeoInit->
- GetTouchableHistory();
- touchvolume = ptrTouchableHistory->GetVolume();
-
- //local point from navigator and normal
- theLocalPoint = ptrNavig->GetCurrentLocalCoordinate();
- normalLoc = touchvolume->GetLogicalVolume()->GetSolid()->
- SurfaceNormal(theLocalPoint);
-
- //global cooordinates normal
- normalGlob = ptrNavig->GetLocalToGlobalTransform().
- TransformAxis(normalLoc);
- }
- else {
- //volume from old history
- const G4NavigationHistory * ptrOldNavHist =
- ptrGeoInit->GetOldNavHist()->GetHistory();
- touchvolume = ptrOldNavHist->GetTopVolume();
-
- if(!touchvolume) {
- // Old history has been reseted by LOOKZ relocation,
- // so is necessary to track back and forward to find
- // the right histories.
-
-
- //////////// COMPUTE STEP BACKWARD //////////////////
- G4ThreeVector theGlobalPoint = ptrNavig->GetLocalToGlobalTransform().
- TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
-
- //compute step and new location
- G4ThreeVector pVec(pVx,pVy,pVz);
- pVec=-pVec;
- G4double appStep = 0;
- G4double safety = 0;
- G4bool onBoundary = false;
- G4double physStep = 1000000000;
- G4int newRegStep;
- G4ThreeVector partLoc = theGlobalPoint;
- G4bool fErr=false;
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "Old history not found" << G4endl;
- G4cout << "* NRML needs boundary-crossing: computing step backward..."
- << G4endl;
-#endif
-
- //compute step and location
- newRegStep=StepAndLocation(partLoc,pVec,physStep,
- appStep,safety,onBoundary,
- fErr,oldReg);
-
- G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::
- GetInstance();
-
- if(appStep<physStep && newRegStep!=G4int(pVolStore->size())+1) {
- //end-step point is on boundary between volumes;
- //==> update step-histories
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* updating step-histories" << G4endl;
-#endif
-
- G4TouchableHistory * ptrTouchHist =
- ptrGeoInit->GetTouchableHistory();
- ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
- }
- else {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "ERROR! Boundary not found" << G4endl;
-#endif
- }
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* computing step forward..." << G4endl;
-#endif
- pVec=-pVec;
- safety = 0;
- onBoundary = false;
-
- //compute step and location for boundary crossing
- newRegStep=StepAndLocation(partLoc,pVec,physStep,
- appStep,safety,onBoundary,fErr,oldReg);
- if(appStep<physStep) {
- //end-step point is on boundary between volumes;
- //==> update step-histories
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* updating step-histories" << G4endl;
-#endif
-
- G4TouchableHistory * ptrTouchHist =
- ptrGeoInit->GetTouchableHistory();
- ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
- }
- }
-
- // now touchvolume exist.
- // local point from navigator and global point
- // N.B. if particle has exited world volume,
- // FluggNavigator doesn't update lastLocatedPoint.
- // So be carefull in building geometry always to have a
- // big world volume that fluka won't exit.
-
- touchvolume = ptrOldNavHist->GetTopVolume();
-
- G4ThreeVector theGlobalPoint = ptrNavig->GetLocalToGlobalTransform().
- TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
- theLocalPoint = ptrOldNavHist->GetTopTransform().
- TransformPoint(theGlobalPoint);
- normalLoc = (touchvolume->GetLogicalVolume()->GetSolid()->
- SurfaceNormal(theLocalPoint));
- normalLoc *= -1;
-
- //global cooordinates normal
- normalGlob = ptrOldNavHist->GetTopTransform().
- Inverse().TransformAxis(normalLoc);
- }
- }
-
- //return normal:
- norml[0]=G4double(normalGlob.x());
- norml[1]=G4double(normalGlob.y());
- norml[2]=G4double(normalGlob.z());
-
- //for debugging
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "Normal: " << normalGlob << G4endl;
-#endif
-}
-
-
+++ /dev/null
-
-// Flugg tag
-
-///////////////////////////////////////////////////////////////////
-//
-// WrapReg.hh - Sara Vanini
-//
-// Wrapper for scoring hits: previous step end-point is taken from
-// history (and compared with fluka region index, flukaReg),
-// then the wrapper returns all the information regarding the
-// volume tree, i.e. returns indMother[] array with all the
-// mother volumes index and repMother[] array with all the
-// mother volumes repetition number.
-//
-// modified: 16/III/99
-// modified: 14/IV/00 ptrLttc included
-// modified: 24.10.00: by I. Hrivnacova
-// functions declarations separated from implementation
-// (moved to Wrappers.hh);
-// modified: 17/06/02 by I. Gonzalez. STL migration.
-//
-///////////////////////////////////////////////////////////////////
-
-#include "Wrappers.hh"
-#include "FGeometryInit.hh"
-#include "NavHistWithCount.hh"
-#include "G4VPhysicalVolume.hh"
-#include "G4ThreeVector.hh"
-#include "G4PhysicalVolumeStore.hh"
-#include "globals.hh"
-
-
-void rgrpwr(const G4int& flukaReg, const G4int& ptrLttc, G4int& g4Reg,
- G4int* indMother, G4int* repMother, G4int& depthFluka)
-{
- //flag
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "============= RGRPWR ==============" << G4endl;
- G4cout << "ptrLttc=" << ptrLttc << G4endl;
-#endif
-
- //Geoinit, Navigator, VolStore pointers
- static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
- G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
-
- //get jrLtGeant array and flag
- G4int * jrLtGeant = ptrGeoInit->GetJrLtGeantArray();
- G4int LttcFlagGeant = ptrGeoInit->GetLttcFlagGeant();
-
- G4bool foundHistory = false;
- G4int i = LttcFlagGeant;
-
- while(!foundHistory && i>=0) {
- if(jrLtGeant[i]==ptrLttc) foundHistory = true;
- i -= 1;
- }
-
- if(!foundHistory) {
- G4cout << "* ERROR! History not in jrLtGeant!" << G4endl;
- //only in debugging version....
- assert(foundHistory);
- }
- else {
- //get history pointer from ptrLttc
- NavHistWithCount* ptrNavHistCount =
- reinterpret_cast<NavHistWithCount*>(ptrLttc);
- G4NavigationHistory* ptrNavHist = ptrNavHistCount->GetNavHistPtr();
-
- G4VPhysicalVolume* ptrVolHistory = 0;
- G4int volHistIndex = 0;
- G4int depth = ptrNavHist->GetDepth();
- for(G4int h=0; h<=depth; h++) {
- ptrVolHistory = ptrNavHist->GetVolume(h);
- //
- volHistIndex = ~0;
- for (unsigned int i=0; i<pVolStore->size(); i++)
- if ((*pVolStore)[i] == ptrVolHistory)
- volHistIndex = i;
- if (volHistIndex==(~0)) {
- G4cerr << "FLUGG: Problem in routine WrapReg tryingto find volume after step" << G4endl;
- exit(-999);
- }
- //volHistIndex = G4int(pVolStore->index(ptrVolHistory));
- //
- indMother[h] = volHistIndex+1;
- if(ptrVolHistory->IsReplicated()) {
- //true if volume is replica or parameterized,
- //false for placements; repetition numbers
- //are set: 1,2,3,etc.; 0 for placed volumes.
- repMother[h] = 1+G4int(ptrNavHist->GetReplicaNo(h));
- }
- else {
- repMother[h] = 0;
- }
-
-#ifdef G4GEOMETRY_DEBUG
- // G4cout<<"Level="<<h<<" : ";
- // G4cout<<"Region="<<indMother[h]<<" (repetition="<<
- // repMother[h]<<")"<<G4endl;
-#endif
- }
-
- //compute new region index
- G4int volIndex = ~0;
- for (unsigned int i=0; i<pVolStore->size(); i++)
- if ((*pVolStore)[i] == ptrVolHistory)
- volIndex = i;
- //G4int volIndex=G4int(pVolStore->index(ptrVolHistory));
- if (volIndex==(~0)) {
- G4cerr << "FLUGG: Problem in routine WrapReg tryingto find volume after step" << G4endl;
- exit(-999);
- }
-
- g4Reg=volIndex+1;
-
- depthFluka = depth;
- }
-
-}
-
-
-
-
-
+++ /dev/null
-
-// Flugg tag
-
-///////////////////////////////////////////////////////////////////
-//
-// WrapSavHist.hh - Sara Vanini
-//
-// Wrapper for saving current navigation history (fCheck=default)
-// and returning its pointer. If fCheck=-1 copy of history pointed
-// by intHist is made in NavHistWithCount object, and its pointer
-// is returned. fCheck=1 and fCheck=2 cases are only in debugging
-// version: an array is created by means of FGeometryInit functions
-// (but could be a static int * ptrArray = new int[10000] with
-// file scope as well) that stores a flag for deleted/undeleted
-// histories and at the end of event is checked to verify that
-// all saved history objects have been deleted.
-//
-// modified 6/III/99: history check array implemented
-// modified 14/IV/00: fCheck=-1 case modified
-// modified 24.10.01: by I. Hrivnacova
-// functions declarations separated from implementation
-// (moved to Wrappers.hh);
-//
-//////////////////////////////////////////////////////////////////
-
-
-#include "Wrappers.hh"
-#include "FGeometryInit.hh"
-#include "NavHistWithCount.hh"
-#include "G4TouchableHistory.hh"
-#include "globals.hh"
-
-G4int isvhwr(const G4int& fCheck, const G4int& intHist)
-{
- //flag
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "============= ISVHWR ==============" << G4endl;
- G4cout << "fCheck=" << fCheck << G4endl;
- if(fCheck==-1)
- G4cout << "intHist=" << intHist << G4endl;
-#endif
-
- //Geoinit pointer
- static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
- static G4int j=0;
-
- switch (fCheck) {
- case 1:
- {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "Start event." << G4endl;
-#endif
-
- return 0;
- }
-
- case 2:
- {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "End event. Check histories in debug-array..." << G4endl;
-#endif
-
- //check that all fluka-banked histories have been delated
- // commented by A.Solodkov
- /*
- G4int * ptrArray = ptrGeoInit->GetHistArray();
-
- for(G4int k=0;k<j;k++) {
- NavHistWithCount* ptrNavHistCount=reinterpret_cast
- <NavHistWithCount*>(ptrArray[k]);
-
- if(ptrArray[k] && !ptrNavHistCount->GetDelateFlag())
- G4cout << "WARNING! History pointed by " <<ptrArray[k]<<
- " has not been deleted at end of event" << G4endl;
- }
-
- //reinitialise debug histories array
- for(G4int i=0;i<1000000;i++) ptrArray[i]=0;
- */
- j=0;
-
- return 0;
- }
-
- case -1:
- {
- //get history object from intHist
- NavHistWithCount* ptrNavHistCountCopy =
- reinterpret_cast<NavHistWithCount*>(intHist);
- G4NavigationHistory* ptrNavHistCopy =
- ptrNavHistCountCopy->GetNavHistPtr();
-
- //copy history in another NavHistWithCount object
- //and save index of check-array
- NavHistWithCount * ptrNavHistCount =
- new NavHistWithCount(*(ptrNavHistCopy));
- ptrNavHistCount->SaveCheckInd(j);
- G4int intHistCopy = G4int(ptrNavHistCount);
-
- //store ptr in array
- // commented by PS
- //G4int * ptrArray = ptrGeoInit->GetHistArray();
- //ptrArray[j]=intHistCopy;
- j+=1;
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "Copying history..." << G4endl;
- G4cout << "Ptr History-copy =" <<intHistCopy<< G4endl;
- G4cout<<*ptrNavHistCopy<< G4endl;
-#endif
-
- return intHistCopy;
- }
-
- default:
- {
- //save G4NavigationHistory and index of check-array
- NavHistWithCount *ptrNavHistCount =
- new NavHistWithCount(*(ptrGeoInit->GetTempNavHist()->GetHistory()));
-
- ptrNavHistCount->SaveCheckInd(j);
- G4int histInt=G4int(ptrNavHistCount);
-
- //store ptr in array
- // comented by PS
- //G4int * ptrArray = ptrGeoInit->GetHistArray();
- // ptrArray[j]=histInt;
- j+=1;
-
-#ifdef G4GEOMETRY_DEBUG
- //TouchableHistory
- G4TouchableHistory * ptrTouchHist = ptrGeoInit->GetTouchableHistory();
- G4cout << "Saving history..." << G4endl;
- G4cout << "Ptr saved History=" << histInt << G4endl;
- G4cout << *(ptrTouchHist->GetHistory()) << G4endl;
-#endif
-
- return histInt;
- }
- }
-}
-
-
-
-
-
-
+++ /dev/null
-#include "WrapUtils.hh"
-#include "FGeometryInit.hh"
-
-////////////////////////////////////////////////////////////////////////
-// StepAndLocation
-////////////////////////////////////////////////////////////////////////
-
-G4int StepAndLocation(G4ThreeVector &partLoc, const G4ThreeVector &pVec,
- const G4double &physStep, G4double &appStep,
- G4double &safety, G4bool & onBound,
- G4bool & flagErr, const G4int & oldReg)
-{
- //NB onBound=true in the particular case when particle is in boundary
- //rounding BUT geometry hasn't crossed it.
- //flagErr=true when particle is on boundary but step=infinity,
- //newSafety>10.e-10 and lLocate function locates end step point in
- // a new region. This is absurb!
-
-
- //Geoinit, Navigator, TouchableHistory, etc. pointers
- static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
- FluggNavigator * ptrNavig = ptrGeoInit->getNavigatorForTracking();
- G4TouchableHistory * ptrTouchableHistory = ptrGeoInit->GetTouchableHistory();
- G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
-
- //compute step
- G4double step = ptrNavig->ComputeStep(partLoc,pVec,physStep,safety);
-
- //compute approved step
- if(step<physStep) {
- //NB:to check if point is really on boundary:
- //see G4AuxiliaryNavServices.hh
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* Boundary crossing!" << G4endl;
-#endif
-
- appStep=step;
- ptrNavig->SetGeometricallyLimitedStep();
-
- //if SetGeometricallyLimitedStep() is called, next
- //LocateGlobalPointAndUpdateTouchable(...,true) will
- //start searching from the end of the previous step,
- //otherwise will start from last located point.
- //(e.g.PropagatorInField and SteppingManger.
- }
- else {
- //step=kInfinity (the nearest boundary is >physStep)
- //Fluka always wants step lenght approved, so:
- appStep=physStep;
- }
-
- //new position
- partLoc+=(pVec*appStep);
-
- //locate point and update touchable history
- ptrNavig->LocateGlobalPointAndUpdateTouchable(partLoc,0,
- ptrTouchableHistory,true);
-
- G4VPhysicalVolume * touchvolume=ptrTouchableHistory->GetVolume();
-
- //if volume not found, out of mother volume:
- //returns [number of volumes]+1
- G4int newReg;
- if(!touchvolume) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "Out of mother volume! (G1WR)" << G4endl;
-#endif
- G4int numVol = G4int(pVolStore->size());
- newReg = numVol + 1;
- }
- else {
-#ifdef G4GEOMETRY_DEBUG
- G4cout <<"Volume after step: " <<touchvolume->GetName() << G4endl;
-#endif
-
- //compute new volume index
- G4int volIndex = ~0;
- for (unsigned int i=0; i<pVolStore->size(); i++)
- if ((*pVolStore)[i] == touchvolume)
- volIndex = i;
- if (volIndex==(~0)) {
- G4cerr << "ERROR in FLUGG: Problem in routine WrapG1 trying to find volume after step" << G4endl;
- exit(-999);
- }
- newReg=volIndex+1;
- }
-
- onBound=false;
- flagErr=false;
- //check if position is in a boundary rounding while volume isn't changed
- if(step==kInfinity || step==physStep) {
- G4ThreeVector globalpoint = ptrNavig->GetLocalToGlobalTransform().
- TransformPoint(ptrNavig->GetCurrentLocalCoordinate());
-
- //compute new safety
- G4double newSafetydbg = ptrNavig->ComputeSafety(globalpoint,DBL_MAX );
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "|From StepAndLocation:" << G4endl;
- G4cout << "|step=" << step << G4endl;
- G4cout << "|phystep=" << physStep << G4endl;
- G4cout << "|safety=" << safety << G4endl;
- G4cout << "|newSafetydbg =" << newSafetydbg << G4endl;
-#endif
-
- if(newSafetydbg<1.0e-10)
- onBound=true;
- else if(newReg != oldReg) {
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "New Volume but ComputeStep didn't notice!" << G4endl;
-#endif
- flagErr=true;
- }
- }
-
- //return
- return newReg;
-}
-
-
-
-
-////////////////////////////////////////////////////////////////////////
-// EqualHistories
-////////////////////////////////////////////////////////////////////////
-
-bool EqualHistories(const G4NavigationHistory* ptrFirstHist,
- const G4NavigationHistory* ptrSecHist)
-{
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "* Testing if histories are equal..." << G4endl;
-#endif
-
- G4int depth1 = ptrFirstHist->GetDepth();
- G4int depth2 = ptrSecHist->GetDepth();
-
- if(depth1!=depth2)
- return false;
-
- for(G4int w=0;w<=depth1;w++) {
- if (*(ptrFirstHist->GetVolume(w))==*(ptrSecHist->GetVolume(w))) {
- //kNormal volume
- if(ptrFirstHist->GetVolumeType(w) == kNormal) {
- if ( ptrFirstHist->GetVolume(w)->GetCopyNo() !=
- ptrSecHist->GetVolume(w)->GetCopyNo() )
- return false;
- }
-
- //Replica or Parametric volume
- else {
- if ( ptrFirstHist->GetReplicaNo(w) !=
- ptrSecHist->GetReplicaNo(w) )
- return false;
- }
- }
- else
- return false;
- }
-
-#ifdef G4GEOMETRY_DEBUG
- G4cout << "Histories are equal!" << G4endl;
-#endif
-
- return true;
-}
-
-
-
-////////////////////////////////////////////////////////////////////////
-// PrintHeader
-////////////////////////////////////////////////////////////////////////
-std::ostream& PrintHeader(std::ostream& os, const char* title) {
- os << "*\n" << "*\n" << "*\n";
- os << "********************* " << title << " *********************\n"
- << "*\n";
- os << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..."
- << G4endl;
- os << "*" << G4endl;
-
- return os;
-}
-
-////////////////////////////////////////////////////////////////////////
-// ToFlukaString converts an string to something usefull in FLUKA:
-// * Capital letters
-// * Only 8 letters
-// * Replace ' ' by '_'
-////////////////////////////////////////////////////////////////////////
-G4String ToFlukaString(G4String str) {
- str = str.substr(0,8);
- str.toUpper();
- for (unsigned int i = 0; i < str.length(); i++) {
- if (str.at(i) == ' ')
- str.replace(i,1,"_",1);
- }
- return str;
-}
+++ /dev/null
-
-// Flugg tag
-
-///////////////////////////////////////////////////////////////////
-//
-// Some utility methods used for several wrapped functions
-//
-///////////////////////////////////////////////////////////////////
-
-#ifndef WrapUtils_hh
-#define WrapUtils_hh 1
-
-//#include <iostream.h>
-#include <iomanip>
-#include "G4NavigationHistory.hh"
-#include "G4ThreeVector.hh"
-
-//Forward declarations
-
-//StepAndLocation declaration. Used in:
-// - WrapG1.cc
-// - WrapNorml.cc
-G4int StepAndLocation(G4ThreeVector &, const G4ThreeVector &,
- const G4double &, G4double &, G4double &, G4bool &,
- G4bool &, const G4int &);
-
-//EqualHistories declaration: true if histories are identical, otherwise false.
-//Used in:
-// - WrapG1.cc
-bool EqualHistories(const G4NavigationHistory*,
- const G4NavigationHistory*);
-
-// Common print utilities used in FGeometryInit.cc
-inline std::ostream& setw10(std::ostream& os) { return os << std::setw(10);}
-inline std::ostream& setscientific(std::ostream& os) { return os << std::setiosflags(std::ios::scientific);}
-inline std::ostream& setfixed(std::ostream& os) { return os << std::setiosflags(std::ios::fixed);}
-std::ostream& PrintHeader(std::ostream& os, const char* title);
-G4String ToFlukaString(G4String str);
-
-#endif
+++ /dev/null
-
-// Flugg tag
-
-///////////////////////////////////////////////////////////////////
-//
-// Interface for Flugg Wrappers
-//
-///////////////////////////////////////////////////////////////////
-
-#ifndef WRAPPERS_HH
-#define WRAPPERS_HH
-
-#include "globals.hh"
-
-#define idnrwr idnrwr_
-#define g1wr g1wr_
-#define g1rtwr g1rtwr_
-#define conhwr conhwr_
-#define inihwr inihwr_
-#define jomiwr jomiwr_
-#define lkdbwr lkdbwr_
-#define lkfxwr lkfxwr_
-#define lkmgwr lkmgwr_
-#define lkwr lkwr_
-#define magfld magfld_
-#define nrmlwr nrmlwr_
-#define rgrpwr rgrpwr_
-#define isvhwr isvhwr_
-
-
-// WrapDN
-
-extern "C" G4int idnrwr(const G4int & nreg, const G4int & mlat);
-
-// WrapG1
-
-extern "C" void g1wr(G4double& pSx, G4double& pSy, G4double& pSz, G4double* pV,
- G4int& oldReg, const G4int& oldLttc, G4double& propStep,
- G4int& nascFlag, G4double& retStep, G4int& newReg,
- G4double& saf, G4int& newLttc, G4int& LttcFlag,
- G4double* sLt, G4int* jrLt);
-
-// WrapG1RT
-
-extern "C" void g1rtwr(void);
-
-// WrapIncrHist
-
-extern "C" void conhwr(G4int& intHist, G4int* incrCount);
-
-// WrapIniHist
-
-extern "C" void inihwr(G4int& intHist);
-
-// WrapInit
-
-extern "C" void jomiwr(const G4int & nge, const G4int& lin, const G4int& lou,
- G4int& flukaReg);
-
-// WrapLookDB
-
-extern "C" void lkdbwr(G4double& pSx, G4double& pSy, G4double& pSz,
- G4double* pV, const G4int& oldReg, const G4int& oldLttc,
- G4int& newReg, G4int& flagErr, G4int& newLttc);
-
-// WrapLookFX
-
-extern "C" void lkfxwr(G4double& pSx, G4double& pSy, G4double& pSz,
- G4double* pV, const G4int& oldReg, const G4int& oldLttc,
- G4int& newReg, G4int& flagErr, G4int& newLttc);
-
-// WrapLookMG
-
-extern "C" void lkmgwr(G4double& pSx, G4double& pSy, G4double& pSz,
- G4double* pV, const G4int& oldReg, const G4int& oldLttc,
- G4int& flagErr, G4int& newReg, G4int& newLttc);
-
-// WrapLookZ
-
-extern "C" void lkwr(G4double& pSx, G4double& pSy, G4double& pSz,
- G4double* pV, const G4int& oldReg, const G4int& oldLttc,
- G4int& newReg, G4int& flagErr, G4int& newLttc);
-
-// WrapMag
-
-extern "C" void magfld(const G4double& pX, const G4double& pY, const G4double& pZ,
- G4double& cosBx, G4double& cosBy, G4double& cosBz,
- G4double& Bmag, G4int& reg, G4int& idiscflag);
-
-// WrapNorml
-
-extern "C" void nrmlwr(G4double& pSx, G4double& pSy, G4double& pSz,
- G4double& pVx, G4double& pVy, G4double& pVz,
- G4double* norml, const G4int& oldReg,
- const G4int& newReg, G4int& flagErr);
-
-// WrapReg
-
-extern "C" void rgrpwr(const G4int& flukaReg, const G4int& ptrLttc, G4int& g4Reg,
- G4int* indMother, G4int* repMother, G4int& depthFluka);
-
-// WrapSavHist
-
-extern "C" G4int isvhwr(const G4int& fCheck, const G4int& intHist);
-
-#endif //WRAPPERS_HH
-
+++ /dev/null
-$Id$
------------------------------------------------------------------
-
-Tags (history):
-===============
-
- 25.10.2001:
- flugg-03_geant4-3-1:
- Wrappers directory restructured to enable creating a library;
- extra fortran source code handled by makefiles (binmake.gmk);
- extended unsufficient arrays size in FGeometryInit.
-
- 31.5.2001:
- flugg-02_geant4-3-1:
- Upgrade to Geant4 3.1 release;
- added scripts for Geant4 upgrades.
-
- 16.5.2001:
- flugg-01_geant4-1-0:
- The cvs keywords (Id, Name) added in the top of files;
- all built-in types changed to G4 typedefs;
- G4std namespace introduced;
- corrections in Fluka input files and in the makefiles.
-
- 10.5.2001:
- flugg-00_geant4-1-0:
- Import of flugg0.4 based on Geant4 1.0
-
+++ /dev/null
-# $Id$
-# Flugg tag $Name$
-
-====================================================================
-====================================================================
-READ-ME-FLUGG - Sara Vanini - June 2000
-Upgrade to Geant4 3.1, I. Hrivnacova, 30 May 2001
-
-(FLUGG documentation: ATL-SOFT-98-039 and ATL-SOFT-99-004.)
-====================================================================
-====================================================================
-
-
-====================================================================
-= 1. Flugg Installation =
-====================================================================
-
-1.1 Supported computers and operating systems
---------------------------------------------------------------------
-FLUGG is supported under the following operating systems:
-
- - Flavors of Unix (from vendors HP, DEC)
- - Linux on PC with g++ (egcs compiler)
-
-Currently, this is the set of flavors which can be associated with
-the environment variable $G4SYSTEM to
-identify the system architecture and compiler used:
-
- UNIX - HP-UX v.10.20, aCC v.1.23 G4SYSTEM: HP-aCC
- DEC-OSF/1 v.4.0, cxx 6.1 G4SYSTEM: DEC-cxx
- Linux - Linux RedHat 6.1, g++ egcs 1.1.2 G4SYSTEM: Linux-g++
-
-1.2 Required software
--------------------------------------------------------------------
-To run FLUGG, the following software must be properly installed in
-your computing environment:
-
- - C++ compiler (compiler from Unix vendor, g++ or Visual C++
- for Windows systems);
- - CLHEP library (see CLHEP reference guide http://wwwinfo.
- cern.ch/asd/lhc++/clhep/manual/RefGuide/index.html);
- - Native STL or ObjectSpace STL
- (see http://www.objectspace.com/)
- - GNU Make (note: g++ preprocessing is used to build file
- dependencies) is also used and a UNIX shell;
- - FLUGG toolkit;
- - FLUKA libraries;
-
-1.3 FLUGG environment
---------------------------------------------------------------------
-Before proceeding with the installation, you need to define some key
-environment variables in your user environment, in order to specify
-where all software components are placed and to set some compilation
-options:
-
-FLUKA path to fluka binary files;
-FORLIB path to the fortran library;
-FLUPRO path to fluka program
-FLUKAOBJ_PATH path to fluka objects;
-FLUKAOBJ name of fluka objects (user's routines, ecc.);
-
-FLUGGINSTALL path where the FLUGG toolkit tree is installed;
-G4SYSTEM set to one of the flavors listed in the above section to
- specify the kind of architecture and compiler used;
-G4GEOMETRY_DEBUG if set the additional debug information is printed
- in the log file (default = unset)
-G4LIB_BUILD_SHARED if set shared libraries are built, static ones otherwise
- (default = unset)
-
-CLHEP_BASE_DIR: path to the CLHEP installation
-OSPACE_BASE_DIR: path to the ObjectSpace STL installation (in case
- ObjectSpace STL implementation is used in place of the system's
- native STL).
-
-
-1.4 How to make G4 libraries
---------------------------------------------------------------------
-At this point, do the following to start building the compilation
-and installation of the kernel libraries.
-You can choose to build libraries in one of two ways, according to
-the needs and system resources. From $FLUGGINSTALL/source:
-
- 1.gmake global
- This will make global libraries, one for each major category.
- 2.gmake
- This will make one library for each "leaf" category (maximum
- library granularity) and produce automatically a map
- of library use and dependencies.
-
-The standard build procedure assumes global libraries if they exist.
-Advantages of using approach 2. can be noticed mainly concerning
-library and program build speed, which in some cases can be improved
-also of a factor 2 or 3 compared to the "global library" approach.
-Using the "granular library" approach a fairly big number of "leaf"
-libraries is produced, dependencies and linking list are however
-evaluated and generated automatically on the fly. The top-level
-GNUmakefile in $FLUGGINSTALL/source parses the dependency files of
-FLUGG and produces libname.map in $G4LIB. libname.map is produced
-by the tool liblist, whose source code is in $FLUGGINSTALL/config.
-When building a binary application the script binmake.gmk in
-$FLUGGINSTALL/config will parse the user's dependency files and use
-libname.map to determine through liblist the required libraries
-to add to the linking list. Only the required libraries will be
-loaded at link time.
-
-
-
-1.5 FLUKA libraries
---------------------------------------------------------------------
-For FLUKA libraries and FLUKA data files, pemf file, etc, you could
-take a look at FLUKA manual.
-
-At this point, you'll be ready to start building your first
-FLUGG application.
-
-
-====================================================================
-= 2. FluGG - Fluka + Geant4 Geometry for Simulation in HEP =
-====================================================================
-
-2.1 How to build a FLUGG example
---------------------------------------------------------------------
-You can create a FLUGG example in the directory:
-$(FLUGGINSTALL)/examples/fluggEx/emptyEx
-
-a)GEOMETRY
-The detector definition requires the representation of its geometrical
-elements, their materials and electronics properties, together with
-user defined properties. FLUGG geometrical representation of
-detector elements requires implementing the following G4 classes:
- - Detector Construction (named "MyDetectorConstruction");
- - Detector Parameterization (optional);
- - Magnetic Field Construction (optional);
-
-You must define these classes in include directory and their
-implementation in src of the above directory. (see http://wwwinfo.cern.ch/
-asd/geant4/G4UsersDocuments/UsersGuides/ForApplicationDeveloper/html/
-index.html for details on how to write a G4 geometry).
-
-As regards geometry construction, remember that replicas and
-parameterized volumes are handled as fluka lattice-volumes, so
-boundaries between them are not seen during tracking time.
-Different placements of identical physical volumes, instead,
-correspond to different fluka-regions (boundaries are seen!).
-
-As regards magnetic field construction:
-1) Fluka needs flagged region for magnetic field. So it is
-imperative to flag logical volumes (in Detector Construction
-file) where magnetic field is present for FLUGG simulation
-(in GEANT4 this isn't necessary, you can define only
-magnetic field too).
-
-How volumes can be flagged:
- 1.a) with flag in G4LogicalVolume constructor:
- logicBox = new G4LogicalVolume(solidBox,BoxMaterial,
- "Box",fieldMgr,0,0);
- 1.b) with member function:
- logicBigBox->SetFieldManager(fieldMgr,true);
- *true flag implies that the field is extended to all
- volume daughters (for field in all detector set this
- flag for world volume!).
-
-2) If field is defined with G4UniformMgnField (so it's uniform)
-the field value is written in fluka card MGNFIELD (in flukaMat.inp
-file). If field is defined with a user-implemented function,
-MAGFLD wrapper gives the local value as fluka calls it run time.
-
-
-b)FLUGG ENVIRONMENT VARIABLES
-Define environment variables described in section 1.3.
-
-c)FLUGG EXECUTABLE
-In the directory: $(FLUGGINSTALL)/examples/fluggEx/emptyEx,
-run "gmake" (or "gmake clean" and then "gmake" to recompile all files).
-Gmake creates the executable ($(FLUGGINSTALL)/bin/Linux-g++/mainFLUGG)
-compiling FLUGG geometry classes, linking G4 geometry libraries
-and fluka libraries (libflukahp.a).
-
-d)FLUKA INPUT
-Now, to run FLUGG, you need fluka-stile input file. You can see FLUKA
-manual for details, and input and pemf file examples in /flugg/
-examples/fluka. Remember that GEOBEGIN-GEOEND card is dummy (you
-have to put no lines between GEOBEGIN and GEOEND, geometry input comes
-from G4 classes!). As regards material specification: in GEANT4
-toolkit, materials and material-volume assignments are specified in
-geometry input classes, in the detector constructor file. When FLUGG
-initializes the geometry, material information is read from
-GEANT4 detector description, and translated into FLUKA-formatted
-input cards.
-So a first dummy run must be executed; the newly created file containing
-GEANT4 material specifications and volume-material assignments
-(flukaMat.inp) must be included into the FLUKA input file where
-additional properties can still be defined. Other useful information,
-like FLUGG geometry volume names and indexes (together with repetition
-numbers for replicated or parameterized volumes!) are dump in
-Volumes_index.inp.
-
-Now you are ready to run mainFLUGG (see rfluka script in /flugg/
-examples/fluka).
-
-
-2.2 Examples module
---------------------------------------------------------------------
-This module collects a set of user examples aimed to demonstrate to
-a user how to make correct use of the FLUGG toolkit by implementing
-those user-classes which a Geant4 user is supposed to customize
-in order to define his/her own detector geometry setup.
-This set of examples covers some possible general use-cases
-for actual detector simulation for HEP.
-
- AlAuAl
- - series of slab of different elements (Al, Au, Al), which total
- thickness is of the order of the electron range;
- - multiple scattering at boundaries (handle delicate situations such
- as grazing angles, backscattering, deflections at boundaries);
- Fluka input files: alaual.inp, wa_50m.pemf
-
- BiasEx
- - 500 MeV protons on a thick Cu target: the generated neutrons
- propagate in a concrete shield;
- - fluka biasing techniques are applied to achieve variance
- reduction;
- Fluka input files: Bias.inp, wa_50m.pemf
-
- MagSphereNotRep
- - test-geometry with spheres and tubs;
- - tracking in magnetic field;
- Fluka input files: MagSphereNotRep.inp, t36.pemf
-
- T36flugg
- - Test-36 em-hadronic calorimeter;
- - Full "ordinary" processes;
- - no magnetic field;
- Fluka input files: T36.inp, t36.pemf
-
-
-
-Appendix A. FLUKA input --> FLUGG input
-********************************************************
-When transposing a fluka input in a flugg input you need
-to:
-
-a) Update the following fluka cards substituting fluka
-region indexes with G4 geometry volume indexes: (in
-parenthesis the WHAT() with region numbers):
-
-ASSIGNMAT (2,3,4)
-BEAMPOS (6)
-BIASING (4,5,6)
-DETECT (6; contin.card:2-6)
-EMF-BIAS (4,5,6)
-EMFCUT (4,5,6)
-EMFRAY (2,3,4)
-EXPTRANS (3,4,5)
-GEOBEGIN/END no WHAT exept COMBINAT; throw
- away everything between GEOBEGIN/END
-GLOBAL (1)
-LOW-BIAS (4,5,6)
-LOW-DOWN (4,5,6)
-RESNUCLEI (5,6)
-STEPSIZE (3,4,5)
-TIMECUT (5,6)
-USRBDX (4,5)
-USRBIN (4,5,6; contin.card:1,2,3,4,5,6)
-USRCOLL (4)
-USRTRACK (4)
-USRYIELD (4,5)
-WW-FACTO (4,5,6)
-
-b)integrate the following material cards with information:
-DELTARAY, EMF-BIAS, EMF-CUT, EMF-FIX, EMF-FLUO, EVXTEST,
-FLUKAFIX, LAM-BIAS, LANDAU, LOW-MAT, MAT-PROP, MULSOPT,
-OPT-PROP, PAIRBREM, PHOTONUC; and in MAPA and PEMF file,
-which contain detailed information about materials.
-
-When comparing fluka and FLUGG runs, REMEMBER to:
-1) check that coordinate system is the same in fluka geometry
-and in GEANT4 geometry; if it is not, update position and beam
-direction.
-2) check how events are the scored!
-
-
-Appendix B. GEANT4 example --> FLUGG example:
-****************************************************************
-0) Copy /src and /include geant4 examples directories in emptyEx/
-1) Delete all classes exept:
- - Detector Construction;
- - Detector Parameterization;
- - Magnetic Field Construction;
-2) Delete corresponding "#include ..." lines in Detector
-Construction source file and variables definitions in the constructor,
-and corresponding class declarations and variable declaration
-in Detector Construction include file.
-2b) Delete #include "G4RunManager.hh" and lines where RunManager is
-invoked.
-3) Delete Detector messenger, initialization and destruction, delete
-code lines where Sensitive Detector is set (be aware of memory leaks!).
-4) Delete code lines where parameterization of processes for fast
-simulation is set (and visualization of ghost volumes).
-Delete commands for interactive definition of the calorimeter
-
-
-
-Appendix C. How to update FLUGG with new Geant4 releases
-****************************************************************
-For updating FLUGG with new Geant4 releases the following
-operations are required:
-
-1) from flugg/source: gmake clean_all to remove the previous
-installation;
-
-2) unzip and untar Geant4 release;
-
-3) set environment variable for paths to Geant4 (G4INSTALL)
- and Flugg (FLUGGINSTALL)
-
-4) run update_source.sh script
-
-This will replace files placed in include and src sub-directories
-of the following directories:
-geant4/source/geometry ---> flugg/source/geometry
-geant4/source/graphics_reps ---> flugg/source/graphics_reps
-geant4/source/material ---> flugg/source/material
-geant4/source/global ---> flugg/source/global
-geant4/source/intercoms ---> flugg/source/intercoms
-The Geant4 makefiles are modified for Flugg; old Flugg
-makefiles are kept and renamed to GNUmakefile.flugg and kept.
-
-5) now you must add "UpdateNavigatorHistory" member function to
-G4Navigator class; insert the following declaration in
-"flugg/source/geometry/volumes/include/G4Navigator.hh" file:
-
-// flugg member function: reinitialization of navigator history with
-// secondary particle history banked on fluka side
- void UpdateNavigatorHistory(const G4NavigationHistory* newNavHistory);
-
-6) And append the following definition to
-"flugg/source/geometry/volumes/src/G4Navigator.cc" file:
-
-void G4Navigator::UpdateNavigatorHistory(const G4NavigationHistory* newNavHistory)
-{
- ResetStackAndState();
- fHistory = *newNavHistory;
- SetupHierarchy();
-}
-
-7) run update_config.sh
-
-This will replace makefiles placed in config directory
-with new ones from Geant4 modified for Flugg;
-old makefiles are kept and renamed (by adding .flugg extension).
-!!! binmake.gmk has to be updated by hands.
-In case no major modifications were done in config
-in new Geant4 release the old binmake.gmk may be used without
-modofications.
-
+++ /dev/null
-FLUGGDUMMYSRCS := G4SDManager.cxx G4VUserDetectorConstruction.cxx \
- G4FastSimulationManager.cxx
-FLUGGDUMMYHDRS := $(FLUGGDUMMYSRCS:.cxx=.hh)
-FLUGGWRAPSRCS := WrapDN.cxx WrapG1.cxx WrapG1RT.cxx \
- WrapIncrHist.cxx WrapIniHist.cxx WrapInit.cxx WrapLookDB.cxx \
- WrapLookFX.cxx WrapLookMG.cxx WrapLookZ.cxx \
- WrapNorml.cxx WrapReg.cxx WrapSavHist.cxx WrapUtils.cxx
-FLUGGWRAPHDRS := WrapUtils.hh Wrappers.hh
-FLUGGOTHERSRCS := FluggNavigator.cxx FGeometryInit.cxx \
- FlukaMaterial.cxx FlukaCompound.cxx \
- FlukaLowMat.cxx
-FLUGGOTHERHDRS := $(FLUGGOTHERSRCS:.cxx=.hh) NavHistWithCount.hh
-
-
-# Sources
-SRCS:= $(FLUGGDUMMYSRCS) $(FLUGGWRAPSRCS) $(FLUGGOTHERSRCS)
-
-
-# Headers
-HDRS:= $(FLUGGDUMMYHDRS) $(FLUGGWRAPHDRS) $(FLUGGOTHERHDRS)
-
-# ROOT Dictionary
-# DHDR:= FluggLinkDef.h
-
-# Extra includes and libraries
-EINCLUDE:= Flugg $(G4INSTALL)/include $(CLHEP_BASE_DIR)/include
-ELIBSDIR:= $(G4INSTALL)/lib/$(G4SYSTEM)
-ELIBS := G4brep G4csg G4geomBoolean G4geometrymng G4globman \
- G4graphics_reps G4hepnumerics G4intercoms G4magneticfield \
- G4materials G4specsolids G4volumes
-
-######################
-#Handle Geant4 flags
-######################
-
-# If G4DEBUG or G4NO_OPTIMISE are not specified,
-# the default compilation is optimised ...
-#
-ifdef G4DEBUG
- CXXFLAGS += -DG4DEBUG
-else
- ifndef G4NO_OPTIMISE
- CXXFLAGS += -DG4OPTIMISE
- endif
-endif
-# Verbosity code can be left out (for better performance)
-# by defining G4_NO_VERBOSE.
-#
-ifndef G4_NO_VERBOSE
- CXXFLAGS += -DG4VERBOSE
-endif
-# Trajectory related classes can be left out (for better performance)
-# by defining G4_NO_STORE_TRAJECTORY.
-#
-ifndef G4_NO_STORE_TRAJECTORY
- CPPFLAGS += -DG4_STORE_TRAJECTORY
-endif
-
-# If GEOMETRY_DEBUG is defined, a lot of information is printed out
-# from FLUGG
-#
-ifdef G4GEOMETRY_DEBUG
- CXXFLAGS += -DG4GEOMETRY_DEBUG
-endif
-
-CXXFLAGS += -DGNU_GCC -DG4USE_STL -DG4USE_STD_NAMESPACE
-
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff