First implementation. Needs cleanup.
authoriglez2 <iglez2@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Nov 2002 17:48:48 +0000 (17:48 +0000)
committeriglez2 <iglez2@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 7 Nov 2002 17:48:48 +0000 (17:48 +0000)
30 files changed:
Flugg/FGeometryInit.cxx [new file with mode: 0644]
Flugg/FGeometryInit.hh [new file with mode: 0644]
Flugg/FGeometryInit.icc [new file with mode: 0644]
Flugg/FluggNavigator.cxx [new file with mode: 0644]
Flugg/FluggNavigator.hh [new file with mode: 0644]
Flugg/G4FastSimulationManager.cxx [new file with mode: 0644]
Flugg/G4FastSimulationManager.hh [new file with mode: 0644]
Flugg/G4SDManager.cxx [new file with mode: 0644]
Flugg/G4SDManager.hh [new file with mode: 0644]
Flugg/G4VUserDetectorConstruction.cxx [new file with mode: 0644]
Flugg/G4VUserDetectorConstruction.hh [new file with mode: 0644]
Flugg/NavHistWithCount.hh [new file with mode: 0644]
Flugg/NavHistWithCount.icc [new file with mode: 0644]
Flugg/WrapDN.cxx [new file with mode: 0644]
Flugg/WrapG1.cxx [new file with mode: 0644]
Flugg/WrapG1RT.cxx [new file with mode: 0644]
Flugg/WrapIncrHist.cxx [new file with mode: 0644]
Flugg/WrapIniHist.cxx [new file with mode: 0644]
Flugg/WrapInit.cxx [new file with mode: 0644]
Flugg/WrapLookDB.cxx [new file with mode: 0644]
Flugg/WrapLookFX.cxx [new file with mode: 0644]
Flugg/WrapLookMG.cxx [new file with mode: 0644]
Flugg/WrapLookZ.cxx [new file with mode: 0644]
Flugg/WrapMag.cxx [new file with mode: 0644]
Flugg/WrapNorml.cxx [new file with mode: 0644]
Flugg/WrapReg.cxx [new file with mode: 0644]
Flugg/WrapSavHist.cxx [new file with mode: 0644]
Flugg/WrapUtils.cxx [new file with mode: 0644]
Flugg/WrapUtils.hh [new file with mode: 0644]
Flugg/Wrappers.hh [new file with mode: 0644]

diff --git a/Flugg/FGeometryInit.cxx b/Flugg/FGeometryInit.cxx
new file mode 100644 (file)
index 0000000..18b376c
--- /dev/null
@@ -0,0 +1,723 @@
+
+// Flugg tag 
+// modified 17/06/02: by I. Gonzalez. STL migration
+
+#include "FGeometryInit.hh"
+#include "FluggNavigator.hh"
+#include "WrapUtils.hh"
+#include <stdio.h>
+#include <iomanip.h>
+
+FGeometryInit * FGeometryInit::flagInstance=0;
+
+FGeometryInit* FGeometryInit::GetInstance()  {
+  cout << "==> Flugg::FGeometryInit::GetInstance(), instance=" << flagInstance
+       << endl;
+  if (!flagInstance) 
+    flagInstance = new FGeometryInit();
+  
+  cout << "<== Flugg::FGeometryInit::GetInstance(), instance=" << flagInstance
+       << endl;
+  return flagInstance;
+}  
+
+
+FGeometryInit::FGeometryInit():
+  fDetector(0),
+  fFieldManager(0),
+  myTopNode(0),
+  ptrGeoMan(0),
+  ptrArray(0),
+  ptrTouchHist(0),
+  ptrOldNavHist(0),
+  ptrTempNavHist(0),
+  ptrJrLtGeant(0),
+  fTransportationManager(G4TransportationManager::GetTransportationManager()){
+
+  cout << "==> Flugg FGeometryInit::FGeometryInit()" << endl;
+  cout << "\t+ Changing the G4Navigator for FluggNavigator..." << endl;
+  G4Navigator* actualnav = fTransportationManager->GetNavigatorForTracking();
+  if (actualnav) {
+    FluggNavigator* newnav = new FluggNavigator();
+    fTransportationManager->SetNavigatorForTracking(newnav);
+  }
+  else {
+    cerr << "ERROR: Could not find the actual G4Navigator" << endl;
+    abort();
+  }
+
+  
+  cout << "<== Flugg FGeometryInit::FGeometryInit()" << endl;
+
+}
+
+
+FGeometryInit::~FGeometryInit() {
+  cout << "==> Flugg FGeometryInit::~FGeometryInit()" << endl;
+  DeleteHistories();
+  ptrGeoMan->OpenGeometry();  
+  if (fTransportationManager)
+    delete fTransportationManager;
+  if (ptrJrLtGeant)
+    delete[] ptrJrLtGeant;
+  DelHistArray();
+  
+  //keep ATTENTION: never delete a pointer twice!
+  cout << "<== Flugg FGeometryInit::FGeometryInit()" << endl;
+}
+
+
+void FGeometryInit::closeGeometry() {
+  cout << "==> Flugg FGeometryInit::closeGeometry()" << endl;
+
+  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;
+  }
+
+  cout << "<== Flugg FGeometryInit::closeGeometry()" << endl;
+}
+
+//*************************************************************************
+
+void FGeometryInit::InitHistArray() {
+  cout << "==> Flugg FGeometryInit::InitHistArray()" << endl;
+  ptrArray = new G4int[1000000];
+  for(G4int i=0;i<1000000;i++) 
+    ptrArray[i]=0;
+  cout << "<== Flugg FGeometryInit::InitHistArray()" << endl;
+}
+
+
+
+//*************************************************************************
+//jrLtGeant stores all crossed lattice volume histories.
+
+void FGeometryInit::InitJrLtGeantArray() {
+  cout << "==> Flugg FGeometryInit::InitJrLtGeantArray()" << endl;
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "Initializing JrLtGeant array" << G4endl;
+#endif
+  ptrJrLtGeant = new G4int[10000];
+  for(G4int x=0;x<10000;x++) 
+    ptrJrLtGeant[x]=-1;
+  flagLttcGeant = -1;
+  cout << "<== Flugg FGeometryInit::InitJrLtGeantArray()" << endl;
+}
+
+
+void FGeometryInit::SetLttcFlagGeant(G4int newFlagLttc) {
+  cout << "==> Flugg FGeometryInit::SetLttcFlagGeant()" << endl;
+  // 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;
+  cout << "<== Flugg FGeometryInit::SetLttcFlagGeant()" << endl;
+}
+
+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() {  
+  cout << "==> Flugg FGeometryInit::InitHistories()" << endl;
+  //init utility histories with navigator history
+  ptrTouchHist = 
+    fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
+  ptrTempNavHist = 
+    fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();   
+  ptrOldNavHist = new G4TouchableHistory();
+  cout << "<== Flugg FGeometryInit::InitHistories()" << endl;
+}
+
+void FGeometryInit::DeleteHistories() {
+  cout << "==> Flugg FGeometryInit::DeleteHistories()" << endl;
+  delete ptrTouchHist;
+  delete ptrOldNavHist;
+  delete ptrTempNavHist;
+  
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "Deleting step-history objects at end of run!" << G4endl;
+#endif
+  cout << "<== Flugg FGeometryInit::DeleteHistories()" << endl;
+}
+
+
+void FGeometryInit::UpdateHistories(const G4NavigationHistory * history,
+                                   G4int flagHist) {
+  cout << "==> Flugg FGeometryInit::UpdateHistories()" << endl;
+  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
+  
+  cout << "<== Flugg FGeometryInit::UpdateHistories()" << endl;
+}
+
+//*****************************************************************************
+
+
+void FGeometryInit::createFlukaMatFile() {
+  cout << "==> Flugg FGeometryInit::createFlukaMatFile()" << endl;
+  // 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 << "================== FILEWR =================" << G4endl;
+#endif 
+  
+  //open file for output
+  ofstream fout("flukaMat.inp");  
+  
+  //PhysicalVolumeStore, Volume and MaterialTable pointers
+  G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
+  unsigned int  numVol = pVolStore->size();
+  
+  G4int* indexMatFluka = 0; 
+  static G4Material * ptrMat = 0;
+  unsigned int x = 0;
+
+  //Find a volume with a material associated
+  while(!ptrMat && x<numVol) {
+    G4VPhysicalVolume * ptrVol = (*pVolStore)[x];
+    G4LogicalVolume * ptrLogVol = ptrVol->GetLogicalVolume();
+    ptrMat = ptrLogVol->GetMaterial();
+    x++;
+  }
+  
+  if(ptrMat) {
+    static const G4MaterialTable * ptrMatTab = G4Material::GetMaterialTable();
+    
+    //number of materials, elements, variable initialisations
+    static size_t totNumMat = G4Material::GetNumberOfMaterials();
+#ifdef G4GEOMETRY_DEBUG
+    G4cout << "Number of materials: " << totNumMat << G4endl;
+#endif 
+
+    static size_t totNumElem = G4Element::GetNumberOfElements();
+#ifdef G4GEOMETRY_DEBUG
+    G4cout << "Number of elements: " << totNumMat << G4endl;
+#endif 
+
+
+    G4int * elemIndexInMATcard = new G4int[totNumElem];
+    for(unsigned int t=0; t<totNumElem; t++) 
+      elemIndexInMATcard[t] = 0;
+    static const G4IsotopeTable * ptrIsotTab = 0;
+    G4int initIsot = 0;
+    static size_t totNumIsot;
+    G4int* isotIndexInMATcard = 0;
+    G4int isotPresence = 0; 
+    
+    
+    // title
+    PrintHeader(fout,"GEANT4 ELEMENTS AND COMPOUNDS");
+    
+    // *** loop over G4Materials to assign Fluka index 
+    G4int indexCount = 3;
+    indexMatFluka = new G4int[totNumMat];
+    for(unsigned int i=0; i<totNumMat; i++) {
+      //pointer material, state 
+      ptrMat = (*ptrMatTab)[i];
+      G4double denMat = ptrMat->GetDensity();
+      G4String nameMat = ptrMat->GetName();
+      // Fluka index: bh=1, vacuum=2, others=3,4..
+      if(denMat<=1.00e-10*g/cm3)
+       //N.B. fluka density limit decided on XI-`98  
+       indexMatFluka[i]=2;     
+      else {  
+       indexMatFluka[i]=indexCount;
+       indexCount++;
+      }
+      
+      // *** write single-element material MATERIAL card
+      size_t numElem = ptrMat->GetNumberOfElements();
+      if(numElem==1) {
+       G4int index = indexMatFluka[i];
+       const G4Element * ptrElem = ptrMat->GetElement(0);
+       size_t indElemTab = ptrElem->GetIndex();
+       size_t numIsot = ptrElem->GetNumberOfIsotopes();
+       G4double A = (ptrElem->GetA())/(g);
+       if(!numIsot) {
+         if(index!=2 && !elemIndexInMATcard[indElemTab]) {
+           G4double Z = ptrElem->GetZ();
+           elemIndexInMATcard[indElemTab] = index;
+           G4String nameEl = ptrElem->GetName();
+           nameEl.toUpper();
+
+           //write on file MATERIAL card of element
+           PrintMaterial(fout, "MATERIAL  ",
+                         Z, A,
+                         denMat/(g/cm3),
+                         index,
+                         -1,
+                         nameEl);
+         }
+       }
+       if(numIsot==1) {
+         // G4Isotope pointer 
+         const G4Isotope* ptrIsot = ptrElem->GetIsotope(0);
+         size_t indIsotTab = ptrIsot->GetIndex();
+         //initialize variables
+         if(!initIsot) {
+           totNumIsot = ptrIsot->GetNumberOfIsotopes();
+           ptrIsotTab = ptrIsot->GetIsotopeTable(); 
+           isotIndexInMATcard = new G4int[totNumIsot];
+           for(unsigned int t=0; t<totNumIsot; t++) 
+             isotIndexInMATcard[t] = 0;
+           initIsot = 1;
+         }
+         if(!isotIndexInMATcard[indIsotTab]) {
+           //compute physical data and counters
+           G4int ZIs = ptrIsot->GetZ();
+           G4double AIs = (ptrIsot->GetA())/(g);
+           G4int NIs = ptrIsot->GetN();
+           G4String nameIsot = ptrIsot->GetName();
+           nameIsot.toUpper();
+           isotIndexInMATcard[indIsotTab] = index;
+           
+           //write on file MATERIAL card of isotope
+           PrintMaterial(fout, "MATERIAL  ",
+                         G4double(ZIs),
+                         AIs,
+                         denMat/(g/cm3),
+                         G4double(index),
+                         G4double(NIs),
+                         nameIsot);
+         }
+       }// end if(numIsot==1)
+      }// end if(numElem==1)
+    }// end for loop
+    
+    
+    // *** material definitions: elements, compound made of G4Elements
+    // or made of G4Materials
+    for(unsigned int j=0; j<totNumMat; j++) {
+      //pointer material, and material data 
+      ptrMat = (*ptrMatTab)[j];
+      size_t numElem = ptrMat->GetNumberOfElements();
+      G4double densityMat = (ptrMat->GetDensity())/(g/cm3);
+      G4String nameMat = ptrMat->GetName();
+      nameMat.toUpper();
+      isotPresence = 0;
+      
+      //fraction vector of compounds of the material
+      const G4double* ptrFracVect = ptrMat->GetFractionVector();
+      
+      //loop on elements of the material
+      for(unsigned int el=0; el<numElem; el++) {
+       //compute physical data, initialize variables
+       const G4Element * ptrElem = ptrMat->GetElement(el);
+       size_t indElemTab = ptrElem->GetIndex();
+       G4String nameElem = ptrElem->GetName();
+       nameElem.toUpper();
+       size_t numIsot = ptrElem->GetNumberOfIsotopes();
+       G4double A = (ptrElem->GetA())/(g);
+       
+       if(!numIsot) {
+         if(!elemIndexInMATcard[indElemTab]) {
+           G4double Z = ptrElem->GetZ();
+           G4double density = ptrFracVect[el]*densityMat;
+           elemIndexInMATcard[indElemTab] = indexCount;
+           
+           //write on file MATERIAL card of element
+           PrintMaterial(fout, "MATERIAL  ",
+                         Z, A,
+                         density,
+                         G4double(indexCount),
+                         -1,
+                         nameElem);
+         }
+       }
+       
+       else {
+         if(numIsot>=2) 
+           isotPresence = 1;
+                               //loop on isotopes
+         for(unsigned int nis=0; nis<numIsot; nis++) {
+           // G4Isotope pointer 
+           const G4Isotope* ptrIsot = ptrElem->GetIsotope(nis);
+           size_t indIsotTab = ptrIsot->GetIndex();
+           //initialize variables
+           if(!initIsot) {
+             totNumIsot = ptrIsot->GetNumberOfIsotopes();
+             ptrIsotTab = ptrIsot->GetIsotopeTable(); 
+             isotIndexInMATcard = new G4int[totNumIsot];
+             for(unsigned int t=0; t<totNumIsot; t++) 
+               isotIndexInMATcard[t] = 0;
+             initIsot = 1;
+           }
+           if(!isotIndexInMATcard[indIsotTab]) {
+             //compute physical data and counters
+             G4int ZIs = ptrIsot->GetZ();
+             G4double AIs = (ptrIsot->GetA())/(g);
+             G4int NIs = ptrIsot->GetN();
+             G4String nameIsot = ptrIsot->GetName();
+             nameIsot.toUpper();
+             G4double* ptrRelAbVect = ptrElem->GetRelativeAbundanceVector();
+             G4double density = 
+               ptrFracVect[el]*densityMat*ptrRelAbVect[nis]*AIs/A;
+             G4int index = indexCount;
+             isotIndexInMATcard[indIsotTab] = indexCount;
+             indexCount+=1;
+             
+             //write on file MATERIAL card of isotope
+             PrintMaterial(fout, "MATERIAL  ",
+                           G4double(ZIs), AIs,
+                           density,
+                           G4double(index),
+                           NIs,
+                           nameIsot);
+           }
+         }
+       }
+       
+      }
+      
+      if(numElem>1 || isotPresence==1) { 
+       // write MATERIAL+COMPOUND card specifing the compound
+       
+       //flags for writing COMPOUND card
+       G4int treCount=0;
+       
+       //make MATERIAL card for compound, start COMPOUND card
+       fout <<"*"<<G4endl;
+       fout <<"*   Define GEANT4 compound " << nameMat << G4endl;
+       PrintMaterial(fout, "MATERIAL  ",
+                     -1, -1,
+                     densityMat,
+                     G4double(indexMatFluka[j]),
+                     -1,
+                     nameMat);
+       fout << setw10 << "COMPOUND  ";
+       
+       
+       //write elements in COMPOUND card
+       for(unsigned int h=0; h<numElem; h++) {
+         const G4Element * ptrElemMat = ptrMat->GetElement(h);
+         size_t indexElemMat = ptrElemMat->GetIndex();
+         size_t numIsotElem = ptrElemMat->GetNumberOfIsotopes();
+         if(!numIsotElem) {
+           PrintCompound(fout, "COMPOUND  ",
+                         treCount,
+                         nameMat,
+                         -ptrFracVect[h],
+                         G4double(elemIndexInMATcard[indexElemMat]));
+
+           if(treCount==3)
+             treCount=0;
+           treCount+=1;
+         }
+         else {
+           G4double * ptrIsotAbbVect = 
+             ptrElemMat->GetRelativeAbundanceVector();
+           
+           for(unsigned int iso=0; iso<numIsotElem; iso++) {
+             const G4Isotope * ptrIsotElem =ptrElemMat->GetIsotope(iso);
+             size_t indexIsotMat = ptrIsotElem->GetIndex();
+             G4double isotAbundPerVol = 
+               ptrIsotAbbVect[iso]*Avogadro*densityMat*
+               ptrFracVect[h]/(ptrElemMat->GetA()/(g));
+             
+             PrintCompound(fout, "COMPOUND  ",
+                           treCount,
+                           nameMat,
+                           isotAbundPerVol,
+                           G4double(isotIndexInMATcard[indexIsotMat]));
+             if(treCount==3)
+               treCount=0;
+             treCount+=1;
+           }
+         }
+       }
+       
+       //end COMPOUND card
+       if(treCount==1) 
+         fout << setw10 << " " << setw10 << " "
+              << setw10 << " " << setw10 << " "
+              << nameMat << endl;
+       if(treCount==2) 
+         fout << setw10 << " " << setw10 << " "
+              << nameMat << endl;
+       if(treCount==3) 
+         fout << nameMat << endl;
+       fout << "*" << endl;
+      }
+      
+      
+    } // end for loop
+    delete elemIndexInMATcard;
+    if(initIsot) 
+      delete isotIndexInMATcard;
+
+  } // end if (ptrMat)
+  
+  // *** material-volume correspondence
+  PrintHeader(fout,"GEANT4 MATERIAL ASSIGNMENTS");
+  
+  //initializations
+  G4int indexMatOld = 0;
+  G4int indexRegFlukaFrom = 0;
+  G4int indexRegFlukaTo = 0;
+  G4int existTo = 0;
+  G4int flagField = 0;
+  G4int lastFlagField = 0;
+  
+  //open file for volume-index correspondence
+  ofstream vout("Volumes_index.inp");
+  
+  //... and write title
+  vout << "*" << endl;
+  vout << "********************  GEANT4 VOLUMES *******************\n";
+  vout << "*" << endl;
+  
+  //loop over all volumes...
+  for(unsigned int l=0;l<numVol;l++) {
+    //index of the region
+    G4VPhysicalVolume * ptrVol = (*pVolStore)[l];
+    G4LogicalVolume * ptrLogVol = ptrVol->GetLogicalVolume();
+    G4int indexRegFluka = l+1;
+    
+    
+    //write index volume and name on file Volumes_index.inp
+    vout.setf(G4std::ios::left,G4std::ios::adjustfield);
+    vout << setw10 << indexRegFluka;
+    vout << G4std::setw(20) << ptrVol->GetName() << G4std::setw(20) << "";
+    if(ptrVol->IsReplicated()) {
+      EAxis axis;
+      G4int nRep;
+      G4double width;
+      G4double offset;
+      G4bool consum;
+      ptrVol->GetReplicationData(axis,nRep,width,offset,consum);
+      vout.setf(G4std::ios::left,G4std::ios::adjustfield);
+      vout << setw10 << "Repetion Nb: " << G4std::setw(3) << nRep;
+    }
+    vout << endl;
+    
+    //check if Magnetic Field is present in the region
+    G4FieldManager * pMagFieldMan = ptrLogVol->GetFieldManager();
+    const G4Field * pMagField = 0;
+    if(pMagFieldMan) 
+      pMagField = pMagFieldMan->GetDetectorField();
+    if(pMagField)        
+      flagField = 1;
+    else         
+      flagField = 0;
+
+
+    //index of material in the region
+    G4Material * ptrMat = ptrLogVol->GetMaterial();
+    G4int indexMat; 
+    if(ptrMat) {
+      size_t indexMatGeant = ptrMat->GetIndex();
+      indexMat = indexMatFluka[indexMatGeant];
+    }
+    else 
+      indexMat = 2;
+       
+    //if materials are repeated
+    if(indexMat==indexMatOld && flagField==lastFlagField) {
+      indexRegFlukaTo=indexRegFluka;
+      existTo=1; 
+      if((l+1)==numVol) {      
+       fout << setw10 << G4double(indexRegFlukaTo);
+       fout << setw10 << "0.0";
+       fout << setw10 << G4double(flagField);
+       fout << setw10 << "0.0" << endl;
+      }
+      
+    }
+    else {
+      //write on file ASSIGNMAT card 
+      
+      //first complete last material card
+      if(!existTo) {
+       if(l) {
+         fout << setw10 << "0.0";
+         fout << setw10 << "0.0";
+         fout << setw10 << G4double(lastFlagField);
+         fout << setw10 << "0.0" << endl;
+       }
+      }
+      else {
+       fout << setw10 <<  G4double(indexRegFlukaTo);
+       fout << setw10 << "0.0";
+       fout << setw10 << G4double(lastFlagField);
+       fout << setw10 << "0.0" << endl;
+      }
+      
+      // begin material card           
+      fout << setw10 << "ASSIGNMAT ";
+      fout.setf(0,G4std::ios::floatfield);     
+      fout << setw10 << setfixed
+          << G4std::setprecision(1) << G4double(indexMat);
+      fout << setw10 << G4double(indexRegFluka);
+      
+      existTo=0;
+      indexRegFlukaFrom=indexRegFluka;
+      
+      if((l+1)==numVol) {      
+       fout << setw10 << "0.0";
+       fout << setw10 << "0.0";
+       fout << setw10 << G4double(flagField);
+       fout << setw10 << "0.0" << endl;
+      }
+    }
+    lastFlagField = flagField;
+    indexMatOld = indexMat; 
+  } // end of loop 
+  
+  //assign material 1 to black-hole=n.vol+1
+  fout << setw10 << "ASSIGNMAT ";
+  fout << setw10 << "1.0";
+  fout << setw10 << G4double(numVol+1);
+  fout << setw10 << "0.0";
+  fout << setw10 << "0.0";
+  fout << setw10 << "0.0";
+  fout << setw10 << "0.0" << endl;
+  
+  // *** magnetic field ***
+  if(fTransportationManager->GetFieldManager()->DoesFieldExist()) {
+    PrintHeader(fout,"GEANT4 MAGNETIC FIELD");
+    
+    //get magnetic field pointer
+    const G4Field * pMagField = 
+      fTransportationManager->GetFieldManager()->GetDetectorField();     
+    
+    //if uniform magnetic field, get value
+    G4double Bx=0.0;
+    G4double By=0.0;
+    G4double Bz=0.0;
+    
+    if(pMagField) {
+      //G4ThreeVector* Field = new G4ThreeVector(1.,2.,3.);
+      const G4UniformMagField *pUnifMagField = 
+       dynamic_cast<const G4UniformMagField*>(pMagField);
+      if(pUnifMagField) {
+       G4double *pFieldValue = 0;
+       G4double *point = new G4double[3];
+       point[0] = 0.;
+       point[1] = 0.;
+       point[2] = 0.;
+       pUnifMagField->GetFieldValue(point,pFieldValue);
+       //non capisco perche' l'instruzione seguente non fa linkare. Indaga!!
+       //per ora posso usare GetFieldValue: va bene lo stesso. 
+       //G4ThreeVector FieldValue = pUnifMagField->GetConstantFieldValue();
+       Bx = pFieldValue[0];
+       By = pFieldValue[1];
+       Bz = pFieldValue[2];
+      }
+      
+    }
+    
+    //write MGNFIELD card 
+    fout << setw10 << "MGNFIELD  ";
+    fout << setw10 << "";
+    fout << setw10 << "";
+    fout << setw10 << "";
+    fout.setf(0,G4std::ios::floatfield);
+    fout << setw10 << setfixed
+        << G4std::setprecision(4) << Bx
+        << setw10 << By
+        << setw10 << Bz 
+        << endl;
+  } // end if magnetic field
+  
+  vout.close();
+  fout.close();
+  delete [] indexMatFluka;
+  cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << endl;
+}
diff --git a/Flugg/FGeometryInit.hh b/Flugg/FGeometryInit.hh
new file mode 100644 (file)
index 0000000..6abeed0
--- /dev/null
@@ -0,0 +1,86 @@
+
+// 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"
+
+
+class FluggNavigator;
+
+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(); 
+
+   private:    
+       FGeometryInit();        //costructor
+       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;
+};
+
+
+//Include the file with the inline methods
+#include "FGeometryInit.icc"
+
+#endif  
diff --git a/Flugg/FGeometryInit.icc b/Flugg/FGeometryInit.icc
new file mode 100644 (file)
index 0000000..0ea504a
--- /dev/null
@@ -0,0 +1,56 @@
+#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() {
+  cout << "\t\t+ fTransportationManager = " << fTransportationManager << endl;
+  cout << "\t\t+ myTopNode = " << myTopNode << endl;
+  cout << "\t\t+ NavigatorForTracking = " 
+       << getNavigatorForTracking() << endl;
+  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;
+}
diff --git a/Flugg/FluggNavigator.cxx b/Flugg/FluggNavigator.cxx
new file mode 100644 (file)
index 0000000..6e4b26e
--- /dev/null
@@ -0,0 +1,55 @@
+//
+// ********************************************************************
+// * 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 "g4std/iomanip"
+
+#define G4DEBUG_NAVIGATION 1
+#define G4VERBOSE 1
+
+FluggNavigator::FluggNavigator() : 
+  G4Navigator()
+{
+  G4cout << "==> Flugg FluggNavigator constructor" << G4endl;
+  G4cout << "\t+fHistory=" << fHistory << ") ..." << G4endl;
+  ResetStackAndState();
+  G4cout << "<== Flugg FluggNavigator constructor" << G4endl;
+}
+    
+void FluggNavigator::UpdateNavigatorHistory(const G4NavigationHistory* newNavHistory)
+{
+  cout << "==> Flugg FluggNavigator::UpdateNavigatorHistory(" << newNavHistory 
+       << ")" << endl;
+  cout << "\t+fHistory=" << fHistory << ") ..." << G4endl;
+  ResetStackAndState();
+  fHistory = *newNavHistory;
+  SetupHierarchy();
+  cout << "<== Flugg FluggNavigator::UpdateNavigatorHistory(" << newNavHistory 
+       << ")" << endl;
+}
diff --git a/Flugg/FluggNavigator.hh b/Flugg/FluggNavigator.hh
new file mode 100644 (file)
index 0000000..b104aeb
--- /dev/null
@@ -0,0 +1,28 @@
+#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 G4std::ostream& operator << (G4std::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
diff --git a/Flugg/G4FastSimulationManager.cxx b/Flugg/G4FastSimulationManager.cxx
new file mode 100644 (file)
index 0000000..2eb4f01
--- /dev/null
@@ -0,0 +1,15 @@
+//---------------------------------------------------------------
+//
+//  dummy G4FastSimulationManager.cc
+//
+//---------------------------------------------------------------
+
+#include "G4FastSimulationManager.hh"
+#include "G4LogicalVolume.hh"
+
+
+G4FastSimulationManager::G4FastSimulationManager(G4LogicalVolume* logVol)
+{;}
+
+G4FastSimulationManager::~G4FastSimulationManager()
+{;}
diff --git a/Flugg/G4FastSimulationManager.hh b/Flugg/G4FastSimulationManager.hh
new file mode 100644 (file)
index 0000000..e47b7a1
--- /dev/null
@@ -0,0 +1,20 @@
+//---------------------------------------------------------------
+//
+//  dummy G4FastSimulationManager.hh
+//
+//---------------------------------------------------------------
+
+
+#ifndef G4FastSimulationManager_h
+#define G4FastSimulationManager_h 1
+
+#include "G4LogicalVolume.hh"
+
+class G4FastSimulationManager
+{
+public:  
+  G4FastSimulationManager(G4LogicalVolume* logVol);
+  ~G4FastSimulationManager();
+};
+
+#endif
diff --git a/Flugg/G4SDManager.cxx b/Flugg/G4SDManager.cxx
new file mode 100644 (file)
index 0000000..c8ea286
--- /dev/null
@@ -0,0 +1,21 @@
+// 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()
+{
+}
+
diff --git a/Flugg/G4SDManager.hh b/Flugg/G4SDManager.hh
new file mode 100644 (file)
index 0000000..84b1d2e
--- /dev/null
@@ -0,0 +1,29 @@
+// 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
+
diff --git a/Flugg/G4VUserDetectorConstruction.cxx b/Flugg/G4VUserDetectorConstruction.cxx
new file mode 100644 (file)
index 0000000..cb11c34
--- /dev/null
@@ -0,0 +1,20 @@
+// 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()
+{;}
+
diff --git a/Flugg/G4VUserDetectorConstruction.hh b/Flugg/G4VUserDetectorConstruction.hh
new file mode 100644 (file)
index 0000000..a3b4b1d
--- /dev/null
@@ -0,0 +1,40 @@
+// 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
+
+
+
+
diff --git a/Flugg/NavHistWithCount.hh b/Flugg/NavHistWithCount.hh
new file mode 100644 (file)
index 0000000..8fbb2c5
--- /dev/null
@@ -0,0 +1,43 @@
+
+// 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
+
+
+
diff --git a/Flugg/NavHistWithCount.icc b/Flugg/NavHistWithCount.icc
new file mode 100644 (file)
index 0000000..d45f4d5
--- /dev/null
@@ -0,0 +1,60 @@
+
+// Flugg tag 
+
+// 
+// NavHistWithCount.icc  - Sara Vanini
+// last modified: 2/II/1999
+//
+// NavHistWithCount.hh inline implementation
+//
+//
+
+inline NavHistWithCount::NavHistWithCount(const G4NavigationHistory &history)
+   {
+   G4cout << "NavHistWithCount::NavHistWithCount created" << G4endl;
+   fhistory = new G4NavigationHistory(history);
+   fDelate=0;
+   count=0;
+   }
+
+inline NavHistWithCount::~NavHistWithCount()
+   {
+   G4cout << "NavHistWithCount::NavHistWithCount deleted" << G4endl;
+   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;
+   }
+
+
+
diff --git a/Flugg/WrapDN.cxx b/Flugg/WrapDN.cxx
new file mode 100644 (file)
index 0000000..4d3021a
--- /dev/null
@@ -0,0 +1,42 @@
+
+// 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
diff --git a/Flugg/WrapG1.cxx b/Flugg/WrapG1.cxx
new file mode 100644 (file)
index 0000000..eb96394
--- /dev/null
@@ -0,0 +1,392 @@
+
+// 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;
+  static G4ThreeVector oldLocalPoint = 
+    ptrNavig->ComputeLocalPoint(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).
+  G4ThreeVector newLocalPoint = ptrNavig->ComputeLocalPoint(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;
+  
+  
+  //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;
+  oldLocalPoint = ptrNavig->ComputeLocalPoint(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
+}
+
+
+
+
diff --git a/Flugg/WrapG1RT.cxx b/Flugg/WrapG1RT.cxx
new file mode 100644 (file)
index 0000000..f3db603
--- /dev/null
@@ -0,0 +1,28 @@
+
+// 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
+}
diff --git a/Flugg/WrapIncrHist.cxx b/Flugg/WrapIncrHist.cxx
new file mode 100644 (file)
index 0000000..19d0d3d
--- /dev/null
@@ -0,0 +1,89 @@
+
+// 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 
+ }
+
+
+
+
+
+
+
+
+
+
diff --git a/Flugg/WrapIniHist.cxx b/Flugg/WrapIniHist.cxx
new file mode 100644 (file)
index 0000000..2ef173b
--- /dev/null
@@ -0,0 +1,85 @@
+
+// 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
+  }
+}
+
+
+
+
+
diff --git a/Flugg/WrapInit.cxx b/Flugg/WrapInit.cxx
new file mode 100644 (file)
index 0000000..f1427b0
--- /dev/null
@@ -0,0 +1,94 @@
+
+// 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
+}
+
+
+
diff --git a/Flugg/WrapLookDB.cxx b/Flugg/WrapLookDB.cxx
new file mode 100644 (file)
index 0000000..68070f2
--- /dev/null
@@ -0,0 +1,32 @@
+
+// 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
diff --git a/Flugg/WrapLookFX.cxx b/Flugg/WrapLookFX.cxx
new file mode 100644 (file)
index 0000000..5ffe449
--- /dev/null
@@ -0,0 +1,112 @@
+
+// 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"
+
+#define G4GEOMETRY_DEBUG 1
+
+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
+  }
+}
+
+
+
+
diff --git a/Flugg/WrapLookMG.cxx b/Flugg/WrapLookMG.cxx
new file mode 100644 (file)
index 0000000..01e3fc6
--- /dev/null
@@ -0,0 +1,325 @@
+
+// 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;
+}
+
+
+
diff --git a/Flugg/WrapLookZ.cxx b/Flugg/WrapLookZ.cxx
new file mode 100644 (file)
index 0000000..fd2302b
--- /dev/null
@@ -0,0 +1,115 @@
+
+// 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"
+
+#define G4GEOMETRY_DEBUG 1
+
+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
+}
+
+
+
+
diff --git a/Flugg/WrapMag.cxx b/Flugg/WrapMag.cxx
new file mode 100644 (file)
index 0000000..ae78fe9
--- /dev/null
@@ -0,0 +1,66 @@
+
+// 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;
+}
+
+
+
+
diff --git a/Flugg/WrapNorml.cxx b/Flugg/WrapNorml.cxx
new file mode 100644 (file)
index 0000000..ed65fc3
--- /dev/null
@@ -0,0 +1,210 @@
+
+// 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
+  if(ptrNavig->IsExitNormalValid()) {
+    normalLoc=ptrNavig->GetLocalExitNormal();
+    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
+}
+
+
diff --git a/Flugg/WrapReg.cxx b/Flugg/WrapReg.cxx
new file mode 100644 (file)
index 0000000..fe35764
--- /dev/null
@@ -0,0 +1,124 @@
+
+// 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;
+  }
+  
+}
+
+
+
+
+
diff --git a/Flugg/WrapSavHist.cxx b/Flugg/WrapSavHist.cxx
new file mode 100644 (file)
index 0000000..d78aab1
--- /dev/null
@@ -0,0 +1,147 @@
+
+// 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;
+    }
+  }
+}
+
+
+
+
+
+
diff --git a/Flugg/WrapUtils.cxx b/Flugg/WrapUtils.cxx
new file mode 100644 (file)
index 0000000..a00f99e
--- /dev/null
@@ -0,0 +1,266 @@
+#include "WrapUtils.hh"
+
+#include "FGeometryInit.hh"
+#include <iostream.h>
+#include <iomanip.h>
+
+////////////////////////////////////////////////////////////////////////
+// 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
+////////////////////////////////////////////////////////////////////////
+ostream& PrintHeader(ostream& os, const char* title) {
+  os << "*\n" << "*\n" << "*\n";
+  os << "*********************  " << title << " *********************\n"
+     << "*\n";
+  os << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..."
+     << endl;
+  os << "*" << endl;
+
+  return os;
+}
+
+////////////////////////////////////////////////////////////////////////
+// PrintMaterial
+////////////////////////////////////////////////////////////////////////
+ostream& PrintMaterial(ostream& os, const char* title,
+                      G4double Z, G4double A,
+                      G4double density,
+                      G4double index,
+                      G4double N,
+                      const char* name) {
+
+  os << setw10 << title;
+
+  os.setf(0,G4std::ios::floatfield);
+  if (Z < 0)
+    os << setw10 << " ";
+  else
+    os << setw10 
+       << setfixed
+       << G4std::setprecision(1) 
+       << Z;
+  
+  if (A < 0)
+    os << setw10 << " ";
+  else
+    os << setw10 << G4std::setprecision(3)
+       << A;
+
+  os.setf(0,G4std::ios::floatfield);
+  os << setw10 
+     << setscientific
+     << G4std::setprecision(3) 
+     << density;
+
+  os.setf(0,G4std::ios::floatfield);
+  os << setw10 
+     << setfixed
+     << G4std::setprecision(1) 
+     << index;
+
+  os << setw10 << " ";
+
+  if (N < 0)
+    os << setw10 << " ";
+  else
+    os << setw10 << N;
+
+  os << name << endl;
+
+  return os;
+}
+
+
+////////////////////////////////////////////////////////////////////////
+// PrintCompund
+////////////////////////////////////////////////////////////////////////
+ostream& PrintCompound(ostream& os, const char* title,
+                      G4int count,
+                      const char* name,
+                      G4double fraction,
+                      G4double index) {
+
+  
+  if(count==3) {
+    os << name << endl;
+    os << setw10 << "COMPOUND  ";
+  }
+  
+  os.setf(0,G4std::ios::floatfield);
+  os << setw10
+     << setscientific
+     << G4std::setprecision(2)
+     << fraction;
+  os.setf(0,G4std::ios::floatfield);
+  os << setw10
+     << setfixed
+     << G4std::setprecision(1)
+     << index;
+
+  return os;
+}
+
diff --git a/Flugg/WrapUtils.hh b/Flugg/WrapUtils.hh
new file mode 100644 (file)
index 0000000..fddb924
--- /dev/null
@@ -0,0 +1,50 @@
+
+// Flugg tag 
+
+///////////////////////////////////////////////////////////////////
+//
+// Some utility methods used for several wrapped functions
+//
+///////////////////////////////////////////////////////////////////
+
+#ifndef WrapUtils_hh
+#define WrapUtils_hh 1
+
+#include "G4NavigationHistory.hh"
+#include "G4ThreeVector.hh"
+#include <iostream.h>
+#include <iomanip.h>
+
+//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*);
+
+// Commonly printed things used in FGeometryInit.cc
+inline ostream& setw10(ostream& os) { return os << G4std::setw(10);}
+inline ostream& setscientific(ostream& os) { return os << G4std::setiosflags(G4std::ios::scientific);}
+inline ostream& setfixed(ostream& os) { return os << G4std::setiosflags(G4std::ios::fixed);}
+ostream& PrintHeader(ostream& os, const char* title);
+ostream& PrintMaterial(ostream& os, const char* title,
+                      G4double Z, G4double A,
+                      G4double density,
+                      G4double index,
+                      G4double N,
+                      const char* name);
+ostream& PrintCompound(ostream& os, const char* title,
+                      G4int count,
+                      const char* name,
+                      G4double fraction,
+                      G4double index);
+
+#endif
diff --git a/Flugg/Wrappers.hh b/Flugg/Wrappers.hh
new file mode 100644 (file)
index 0000000..707fad5
--- /dev/null
@@ -0,0 +1,109 @@
+
+// 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_
+
+
+#define G4GEOMETRY_DEBUG 1
+
+// 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
+