Updated to new way of producing the input material/compound file
authoriglez2 <iglez2@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 21 Nov 2002 12:29:28 +0000 (12:29 +0000)
committeriglez2 <iglez2@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 21 Nov 2002 12:29:28 +0000 (12:29 +0000)
Flugg/FGeometryInit.cxx
Flugg/FGeometryInit.hh
Flugg/FlukaCompound.cxx [new file with mode: 0644]
Flugg/FlukaCompound.hh [new file with mode: 0644]
Flugg/FlukaLowMat.cxx [new file with mode: 0644]
Flugg/FlukaLowMat.hh [new file with mode: 0644]
Flugg/FlukaMaterial.cxx [new file with mode: 0644]
Flugg/FlukaMaterial.hh [new file with mode: 0644]
Flugg/WrapUtils.cxx
Flugg/WrapUtils.hh
Flugg/libFlugg.pkg

index 7c88b86..45b0494 100644 (file)
@@ -7,6 +7,8 @@
 #include "FGeometryInit.hh"
 #include "FluggNavigator.hh"
 #include "WrapUtils.hh"
+#include "FlukaMaterial.hh"
+#include "FlukaCompound.hh"
 
 FGeometryInit * FGeometryInit::flagInstance=0;
 
@@ -277,7 +279,6 @@ void FGeometryInit::UpdateHistories(const G4NavigationHistory * history,
 
 //*****************************************************************************
 
-
 void FGeometryInit::createFlukaMatFile() {
   // last modification Sara Vanini 1/III/99
   // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case,
@@ -291,469 +292,447 @@ void FGeometryInit::createFlukaMatFile() {
   G4cout << "==> Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
   G4cout << "================== FILEWR =================" << G4endl;
 #endif 
-  
-  //open file for output
-  G4std::ofstream fout("flukaMat.inp");  
-  
-  //PhysicalVolumeStore, Volume and MaterialTable pointers
+
+
+  //Regions map
+  BuildRegionsMap();
+  G4std::ofstream vos("Volumes_index.inp");
+  PrintRegionsMap(vos);
+  vos.close();
+
+  //Materials and compounds
+  BuildMaterialTables();
+  G4std::ofstream fos("flukaMat.inp");  
+  PrintMaterialTables(fos);
+  PrintAssignmat(fos);
+  PrintMagneticField(fos);
+  fos.close();
+
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
+#endif
+}
+
+////////////////////////////////////////////////////////////////////////
+// 
+void FGeometryInit::BuildRegionsMap() {
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
+#endif
+
+  //Find number of Volumes in physical volume store
   G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
-  unsigned int  numVol = pVolStore->size();
-  
-  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++;
+  unsigned int numVol = pVolStore->size();
+
+  G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore 
+        << ") has " << numVol << " volumes. Iterating..." 
+        << G4endl;
+
+  for(unsigned int l=0; l < numVol; l++) {
+    //Get each of the physical volumes
+    G4VPhysicalVolume * physicalvolume = (*pVolStore)[l];
+    G4int iFlukaRegion = l+1;
+    fRegionVolumeMap[physicalvolume] = iFlukaRegion;
   }
-  
-  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 
+  G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
+#endif
+}
 
-    static size_t totNumElem = G4Element::GetNumberOfElements();
+void FGeometryInit::PrintRegionsMap(G4std::ostream& os) {
 #ifdef G4GEOMETRY_DEBUG
-    G4cout << "Number of elements: " << totNumMat << G4endl;
-#endif 
+  G4cout << "==> Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
+#endif
 
+  //Print some header
+  PrintHeader(os, "GEANT4 VOLUMES");
 
-    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
-    
+  //Iterate over all volumes in the map
+  for (RegionIterator i = fRegionVolumeMap.begin(); 
+       i != fRegionVolumeMap.end(); 
+       i++) {
     
-    // *** 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 << G4endl;
-       if(treCount==2) 
-         fout << setw10 << " " << setw10 << " "
-              << nameMat << G4endl;
-       if(treCount==3) 
-         fout << nameMat << G4endl;
-       fout << "*" << G4endl;
-      }
-      
-      
-    } // end for loop
-    delete elemIndexInMATcard;
-    if(initIsot) 
-      delete isotIndexInMATcard;
+    //Get info in the map
+    G4VPhysicalVolume* ptrVol = (*i).first;
+    int index = (*i).second;
 
-  } // 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
-  G4std::ofstream vout("Volumes_index.inp");
-  
-  //... and write title
-  vout << "*" << G4endl;
-  vout << "********************  GEANT4 VOLUMES *******************\n";
-  vout << "*" << G4endl;
-  
-  //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) << "";
+    //Print index and region name in some fixed format
+    os.setf(G4std::ios::left, G4std::ios::adjustfield);
+    os << setw10 << index;
+    os << G4std::setw(20) << ptrVol->GetName() << G4std::setw(20) << "";
+
+    //If volume is a replica... print some more stuff
     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;
+      G4int nRep = -1;
+      G4double width = -1;
+      G4double offset = -1;
+      G4bool consum = false;
+      ptrVol->GetReplicationData(axis, nRep, width, offset, consum);
+      os.setf(G4std::ios::left, G4std::ios::adjustfield);
+      os << setw10 << "Repetion Nb: " << G4std::setw(3) << nRep;
     }
-    vout << G4endl;
+    os << G4endl;
     
-    //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" << G4endl;
-      }
-      
+  }
+  
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "<== Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
+#endif
+}
+
+////////////////////////////////////////////////////////////////////////
+// 
+void FGeometryInit::BuildMaterialTables() {
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "==> Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
+#endif
+
+  //some terminal printout also
+  G4cout << "\t* Storing information..." << G4endl;
+
+  //The logic is the folloing:
+  //Get the Material Table and:
+  // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
+  // 2) For each single element material build a material equivalent
+  // 3) For the rest:
+  //   3.a) Build materials for each not already known element
+  //   3.b) Build the compound out of them
+
+  //Get the Material Table and iterate
+  const G4MaterialTable* matTable = G4Material::GetMaterialTable();
+  for (MatTableIterator i = matTable->begin(); i != matTable->end(); i++) {
+
+    //Get some basic material information
+    G4Material* material = (*i);
+    G4String matName = material->GetName();
+    const G4double matDensity = material->GetDensity();
+    const G4int nMatElements  = material->GetNumberOfElements();
+
+    G4cout << "\t\t+ " << matName 
+          << ": dens. = " << matDensity/(g/cm3) << "g/cm3"
+          << ", nElem = " << nMatElements << G4endl;
+
+    // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
+    //    FlukaMaterial* is 0 in that case
+    if (matDensity <= 1.00e-10*g/cm3) {
+      G4FlukaMaterialMap[material] = 0;
+      G4cout << "\t\t  Stored as vacuum" << G4endl;
     }
-    else {
-      //write on file ASSIGNMAT card 
+    // 2) For each single element material build a material equivalent
+    else if (nMatElements == 1) {
       
-      //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" << G4endl;
-       }
-      }
-      else {
-       fout << setw10 <<  G4double(indexRegFlukaTo);
-       fout << setw10 << "0.0";
-       fout << setw10 << G4double(lastFlagField);
-       fout << setw10 << "0.0" << G4endl;
-      }
+      FlukaMaterial *flukamat = 
+       BuildFlukaMaterialFromElement(material->GetElement(0),
+                                     matDensity);
       
-      // begin material card           
-      fout << setw10 << "ASSIGNMAT ";
-      fout.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);  
-      fout << setw10 << setfixed
-          << G4std::setprecision(1) << G4double(indexMat);
-      fout << setw10 << G4double(indexRegFluka);
+      G4FlukaMaterialMap[material] = flukamat;
+      G4cout << "\t\t  Stored as " << flukamat->GetRealName() << G4endl;
       
-      existTo=0;
-      indexRegFlukaFrom=indexRegFluka;
+    } //else if (material->GetNumberOfElements() == 1)
+    
+    // 3) For the rest:
+    //   3.a) Build materials for each not already known element
+    //   3.b) Build the compound out of them
+    else {
+      FlukaCompound* flukacomp = 
+       BuildFlukaCompoundFromMaterial(material);
+      G4FlukaCompoundMap[material] = flukacomp;
+      G4cout << "\t\t  Stored as " << flukacomp->GetRealName() << G4endl;
+    } //else for case 3)
+  } //for (materials)
+  
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "<== Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
+#endif
+}
+
+FlukaMaterial* 
+FGeometryInit::BuildFlukaMaterialFromElement(const G4Element* element,
+                                            G4double matDensity) {
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromElement()" 
+        << G4endl;
+#endif
+
+  //Get element and its properties
+  G4String elemName(ToFlukaString(element->GetName()));
+
+  FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(elemName);
+  if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
+    //Check for isotopes
+    G4int nIsotopes = element->GetNumberOfIsotopes();
+    if (nIsotopes == 0) {
+      G4double elemA = element->GetA()/g;
+      G4double elemZ = element->GetZ();
       
-      if((l+1)==numVol) {      
-       fout << setw10 << "0.0";
-       fout << setw10 << "0.0";
-       fout << setw10 << G4double(flagField);
-       fout << setw10 << "0.0" << G4endl;
+      if (elemA != G4int(elemA) && elemZ != G4int(elemZ)) {
+       G4cout << "WARNING: Element \'" << elemName 
+              << "\' has non integer Z (" << elemZ << ") or A (" 
+              << elemA << ")"
+              << G4endl;
       }
+    
+      flukamat = new FlukaMaterial(elemName,
+                                  G4int(elemZ),
+                                  elemA,
+                                  matDensity/(g/cm3));
+    }
+    else if (nIsotopes == 1) {
+      const G4Isotope* isotope = element->GetIsotope(0);
+      flukamat = BuildFlukaMaterialFromIsotope(isotope, matDensity);
+    }
+    else {
+      FlukaCompound *flucomp = BuildFlukaCompoundFromElement(element,
+                                                            matDensity);
+      flukamat = flucomp->GetFlukaMaterial();      
     }
-    lastFlagField = flagField;
-    indexMatOld = indexMat; 
-  } // end of loop 
+  }
+#ifdef G4GEOMETRY_DEBUG
+  else {
+    G4cout << "INFO: Element \'" << elemName 
+          << "\' already exists in the DB. It will not be recreated."
+          << G4endl;
+  }
+#endif
+
+  return flukamat;
   
-  //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" << G4endl;
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "<== Flugg FGeometryInit::BuildFlukaMaterialFromElement()" 
+        << G4endl;
+#endif
+}
+
+FlukaMaterial* 
+FGeometryInit::BuildFlukaMaterialFromIsotope(const G4Isotope* isotope,
+                                            G4double matDensity) {
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()" 
+        << G4endl;
+#endif
+  G4String isoName(ToFlukaString(isotope->GetName()));
+  FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(isoName);
+  if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
+    G4int isoZ = isotope->GetZ();
+    G4double isoA = (isotope->GetA())/(g);
+    G4int isoN = isotope->GetN();
+    flukamat = new FlukaMaterial(isoName,
+                                isoZ,
+                                isoA,
+                                matDensity/(g/cm3),
+                                isoN);
+  }
+
+  return flukamat;
+
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()" 
+        << G4endl;
+#endif
+}
+
+FlukaCompound* 
+FGeometryInit::BuildFlukaCompoundFromMaterial(const G4Material* material) {
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()" 
+        << G4endl;
+#endif
+  //Material properties
+  const G4double* elemFractions = material->GetFractionVector();
+  const G4int nMatElements  = material->GetNumberOfElements();
+  const G4double matDensity = material->GetDensity();
+  G4String matName(ToFlukaString(material->GetName()));
+  FlukaCompound* flukacomp = new FlukaCompound(matName, matDensity/(g/cm3),
+                                              nMatElements);
+  for (G4int i = 0; i < nMatElements; i++) {
+    FlukaMaterial *flukamat = 
+      BuildFlukaMaterialFromElement(material->GetElement(i), 0.0);
+    
+    flukacomp->AddElement(flukamat->GetIndex(), -elemFractions[i]);
+    
+  } //for (elements)
+
+  return flukacomp;
+
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()" 
+        << G4endl;
+#endif
+}
+
+FlukaCompound* 
+FGeometryInit::BuildFlukaCompoundFromElement(const G4Element* element,
+                                            G4double matDensity) {
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromElement()" 
+        << G4endl;
+#endif
+  G4int nIsotopes = element->GetNumberOfIsotopes();
+  //fraction of nb of atomes per volume (= volume fraction?)
+  const G4double* isoAbundance =  element->GetRelativeAbundanceVector();
+  G4String elemName(ToFlukaString(element->GetName()));
+
+  //Material properties
+  FlukaCompound* flukacomp = new FlukaCompound(elemName, matDensity/(g/cm3),
+                                              nIsotopes);
+  for (G4int i = 0; i < nIsotopes; i++) {
+    FlukaMaterial *flukamat = 
+      BuildFlukaMaterialFromIsotope(element->GetIsotope(i), 0.0);
+    
+    flukacomp->AddElement(flukamat->GetIndex(), isoAbundance[i]);
+    
+  } //for (elements)
+
+  return flukacomp;
+
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromElement()" 
+        << G4endl;
+#endif
+}
+
+void FGeometryInit::PrintMaterialTables(G4std::ostream& os) {
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "==> Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
+#endif
+  //Print Header
+  PrintHeader(os, "GEANT4 MATERIALS AND COMPOUNDS");
   
-  // *** magnetic field ***
+  //And some more stuff  
+  size_t nIsotopes = G4Isotope::GetNumberOfIsotopes();
+  size_t nElements = G4Element::GetNumberOfElements();
+  size_t nMaterials = G4Material::GetNumberOfMaterials();
+
+  os << "* In Geant4 there are " << nMaterials << " materials" << endl;
+  os << "* In Geant4 there are " << nElements  << " elements"  << endl;
+  os << "* In Geant4 there are " << nIsotopes  << " isotopes"  << endl;
+
+  //Materials
+  G4cout << "\t* Printing FLUKA materials..." << G4endl;
+  FlukaMaterial::PrintMaterialsByIndex(os);
+  //FlukaMaterial::PrintMaterialsByName(os);
+
+  //Compounds
+  G4cout << "\t* Printing FLUKA compounds..." << G4endl;
+  FlukaCompound::PrintCompounds(os);
+
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "<== Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
+#endif
+}
+
+////////////////////////////////////////////////////////////////////////
+// 
+void FGeometryInit::PrintAssignmat(G4std::ostream& os) {
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
+#endif
+
+  //Find number of Volumes in physical volume store
+  G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
+  unsigned int numVol = pVolStore->size();
+
+  G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore 
+        << ") has " << numVol << " volumes. " << G4endl;
+
+  G4cout << "\t* Printing ASSIGNMAT..." << G4endl;
+
+
+  PrintHeader(os,"GEANT4 MATERIAL ASSIGNMENTS");
+  for(unsigned int l=0; l < numVol; l++) {
+
+    //Get each of the physical volumes
+    G4VPhysicalVolume * physicalvol = (*pVolStore)[l];
+
+    //Get index for that volume
+    G4int iFlukaRegion = fRegionVolumeMap[physicalvol];
+
+    //Find G4 material and navigate to its fluka compound/material
+    G4LogicalVolume * logicalVol = physicalvol->GetLogicalVolume();
+    G4Material* material = logicalVol->GetMaterial();
+    G4int matIndex = 2;
+    if (G4FlukaCompoundMap[material])
+      matIndex = G4FlukaCompoundMap[material]->GetIndex();
+    if (G4FlukaMaterialMap[material])
+      matIndex = G4FlukaMaterialMap[material]->GetIndex();
+
+    //Find if there is a magnetic field in the region
+    //check if Magnetic Field is present in the region
+    G4double flagField = 0.0;
+    G4FieldManager * pMagFieldMan = logicalVol->GetFieldManager();
+    if(pMagFieldMan && pMagFieldMan->GetDetectorField())
+      flagField = 1.0;
+    
+    //Print card
+    os << setw10 << "ASSIGNMAT ";
+    os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
+    os << setw10 << setfixed << G4double(matIndex);
+    os << setw10 << setfixed << G4double(iFlukaRegion);
+    os << setw10 << "0.0";
+    os << setw10 << setfixed << flagField;
+    os << G4endl;
+  }
+
+
+
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
+#endif
+}
+
+
+void FGeometryInit::PrintMagneticField(G4std::ostream& os) {
+#ifdef G4GEOMETRY_DEBUG
+  G4cout << "==> Flugg FGeometryInit::PrintMagneticField()" << G4endl;
+#endif
+
+  G4cout << "\t* Printing Magnetic Field..." << G4endl;
+
   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.);
+      //Check if it can be made a uniform magnetic field
       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];
+       G4double B[3];
+       G4double point[4]; //it is not really used
+       pUnifMagField->GetFieldValue(point,B);
+
+       //write MGNFIELD card 
+       PrintHeader(os,"GEANT4 MAGNETIC FIELD");
+       os << setw10 << "MGNFIELD  ";
+       os << setw10 << "";
+       os << setw10 << "";
+       os << setw10 << "";
+       os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
+       os << setw10 << setfixed
+          << G4std::setprecision(4) << B[0]
+          << setw10 << B[1]
+          << setw10 << B[2]
+          << G4endl;
+      }
+      else {
+       G4cout << "WARNING: No Uniform Magnetic Field found." << G4endl;
+       G4cout << "         Manual intervention might be needed." << G4endl;
       }
-      
     }
-    
-    //write MGNFIELD card 
-    fout << setw10 << "MGNFIELD  ";
-    fout << setw10 << "";
-    fout << setw10 << "";
-    fout << setw10 << "";
-    fout.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
-    fout << setw10 << setfixed
-        << G4std::setprecision(4) << Bx
-        << setw10 << By
-        << setw10 << Bz 
-        << G4endl;
+    else
+      G4cout << "\t  No detector field found... " << G4endl;
   } // end if magnetic field
-  
-  vout.close();
-  fout.close();
-  delete [] indexMatFluka;
+  else
+    G4cout << "\t  No field found... " << G4endl;
+
 #ifdef G4GEOMETRY_DEBUG
-  G4cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
+  G4cout << "<== Flugg FGeometryInit::PrintMagneticField()" << G4endl;
 #endif
 }
index 6abeed0..d50eac6 100644 (file)
 #include "G4TransportationManager.hh"
 
 
+#include <map>
+
 class FluggNavigator;
+class FlukaMaterial;
+class FlukaCompound;
+
+class FGeometryInit : public G4TransportationManager {
+
+public:
+  ~FGeometryInit();    //destructor
+  static FGeometryInit *GetInstance();
+  inline FluggNavigator *getNavigatorForTracking();
+  inline G4FieldManager * getFieldManager();
+  inline void setDetConstruction(G4VUserDetectorConstruction* detector);
+  inline void setDetector();
+  inline void setMotherVolume();
+  void createFlukaMatFile();
+  void closeGeometry();
+  
+  void PrintHistories();
+  void InitHistories();
+  void DeleteHistories();
+  void UpdateHistories(const G4NavigationHistory *, G4int);
+  inline G4TouchableHistory * GetTouchableHistory();
+  inline G4TouchableHistory * GetOldNavHist();
+  inline G4TouchableHistory * GetTempNavHist();
+  
+  void InitHistArray();
+  inline void DelHistArray();
+  inline G4int * GetHistArray();
+  
+  void InitJrLtGeantArray();
+  inline G4int * GetJrLtGeantArray();
+  inline G4int GetLttcFlagGeant();
+  void SetLttcFlagGeant(G4int);
+  void PrintJrLtGeant(); 
+
+protected:
+  void BuildRegionsMap();
+  void PrintRegionsMap(G4std::ostream& os);
+  void BuildMaterialTables();
+  FlukaMaterial* BuildFlukaMaterialFromElement(const G4Element* element, 
+                                              G4double matDensity);
+  FlukaMaterial* BuildFlukaMaterialFromIsotope(const G4Isotope* isotope,
+                                              G4double matDensity);
+  FlukaCompound* BuildFlukaCompoundFromMaterial(const G4Material* material); 
+  FlukaCompound* BuildFlukaCompoundFromElement(const G4Element* element,
+                                              G4double matDensity);
+  void PrintMaterialTables(G4std::ostream& os);
+  void PrintAssignmat(G4std::ostream& os);
+  void PrintMagneticField(G4std::ostream& os);
 
-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;
+private:    
+  FGeometryInit();     //costructor
+
+private:
+  G4VUserDetectorConstruction * fDetector;
+  G4FieldManager * fFieldManager; 
+  G4TransportationManager * fTransportationManager;
+  static FGeometryInit *flagInstance;
+  G4VPhysicalVolume * myTopNode;
+  G4GeometryManager * ptrGeoMan;
+  G4int * ptrArray;
+  G4TouchableHistory * ptrTouchHist;
+  G4TouchableHistory * ptrOldNavHist;
+  G4TouchableHistory * ptrTempNavHist;
+  G4int * ptrJrLtGeant;
+  G4int flagLttcGeant;
+
+  G4std::map<G4VPhysicalVolume*, int, G4std::less<G4VPhysicalVolume*> > fRegionVolumeMap;
+  G4std::map<G4Material*, FlukaMaterial*, G4std::less<G4Material*> > G4FlukaMaterialMap;
+  G4std::map<G4Material*, FlukaCompound*, G4std::less<G4Material*> > G4FlukaCompoundMap;
+  //G4int NOfMaterials;
 };
 
+typedef  G4std::map<G4VPhysicalVolume*, int, G4std::less<G4VPhysicalVolume*> >::iterator RegionIterator;
+typedef  G4std::vector<G4Material*>::const_iterator MatTableIterator;
+
 
 //Include the file with the inline methods
 #include "FGeometryInit.icc"
diff --git a/Flugg/FlukaCompound.cxx b/Flugg/FlukaCompound.cxx
new file mode 100644 (file)
index 0000000..a4eccb4
--- /dev/null
@@ -0,0 +1,109 @@
+#include "FlukaCompound.hh"
+#include "G4ios.hh"
+#include "WrapUtils.hh"
+
+FlukaCompoundsTable FlukaCompound::fFlukaCompounds;
+
+FlukaCompound::FlukaCompound(const G4String& name, 
+                            G4double density, 
+                            G4int nelem):
+  fNMaterials(nelem),
+  fNAdded(0) {
+  //Initialise arrays for element indices and fractions
+  fElIndex = new G4int[nelem];
+  fFraction = new G4double[nelem];
+  for (G4int i = 0; i < nelem; i++) {
+    fElIndex[i] = 0;
+    fFraction[i] = 0;
+  }
+
+  //Build the associated material
+  fFlukaMaterial = new FlukaMaterial(name, 0, 0, density);
+
+  //If it already exists it will have got a new name
+  G4String testname(fFlukaMaterial->GetRealName());
+  if (testname != name) {
+    G4cerr << "INFO: Found two materials with the same name! ("
+          << name << ")" << G4endl;
+    G4cerr << "      Renaming to " << testname << G4endl;
+  }
+  fFlukaCompounds[testname] = this;
+}
+
+FlukaCompound::~FlukaCompound() {
+  delete[] fElIndex;
+  delete[] fFraction;
+}
+
+
+void FlukaCompound::AddElement(G4int index, G4double fraction) {
+  if (fNAdded < fNMaterials) {
+    fElIndex[fNAdded] = index;
+    fFraction[fNAdded] = fraction;
+    fNAdded++;
+  }
+  else {
+    G4cerr << "ERROR: Trying to add to many elements to compound \'" 
+          << GetName() << "\'!" << G4endl;
+    G4cerr << "       Last element not added!!!" << G4endl;
+  }
+}
+
+G4std::ostream& FlukaCompound::PrintCompounds(G4std::ostream& os) {
+  PrintHeader(os, "COMPOUNDS");  
+  
+  for (FlukaCompoundsIterator i = fFlukaCompounds.begin(); 
+       i != fFlukaCompounds.end(); 
+       i++) {
+    FlukaCompound* flucomp     = (*i).second;
+    os << *flucomp;
+  }
+
+  return os;
+}
+
+G4std::ostream& operator<<(G4std::ostream& os, const FlukaCompound& flucomp) {
+  G4int nmats = flucomp.GetNMaterials();
+  G4String matName = flucomp.GetName();
+  G4String matRealName = flucomp.GetRealName().substr(0,8);
+  //Some comment
+  os << "* " << matName << " COMPOUND (" << nmats << ")" << G4endl;
+  
+  //Material card
+  os << *(flucomp.GetFlukaMaterial());
+
+  //The card
+  G4int counttothree = 0;
+  os << setw10 <<"COMPOUND  ";
+  for (G4int i = 0; i < nmats; i++) {
+    os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
+    os << setw10
+       << setfixed
+       << G4std::setprecision(6)
+       << flucomp.GetMaterialFraction(i);
+    os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
+    os << setw10
+       << setfixed
+       << G4std::setprecision(1)
+       << G4double(flucomp.GetMaterialIndex(i));
+    counttothree++;
+    if (counttothree == 3 ) {
+      os << matRealName;
+      os << G4endl;
+      if ( (i+1) != nmats)
+       os << setw10 <<"COMPOUND  ";
+      counttothree = 0;
+    }
+  }
+
+  if ((nmats % 3) != 0) {
+    //Unless we have 3, 6, 9... submaterials we need to put some empty
+    //space and the compound name
+    for (G4int i = 0; i < (3 - (nmats % 3)); i++)
+      os << setw10 << " " << setw10 << " ";
+    os << matRealName;
+    os << G4endl;
+  }
+  
+  return os;
+}
diff --git a/Flugg/FlukaCompound.hh b/Flugg/FlukaCompound.hh
new file mode 100644 (file)
index 0000000..ac5eebe
--- /dev/null
@@ -0,0 +1,75 @@
+#ifndef FLUKACOMPOUND_HH
+#define FLUKACOMPOUND_HH 1
+
+#include "G4String.hh"
+#include "FlukaMaterial.hh"
+
+#include <map>
+
+class FlukaCompound;
+typedef G4std::map<G4String, FlukaCompound*, G4std::less<G4String> > FlukaCompoundsTable;
+typedef G4std::map<G4String, FlukaCompound*, G4std::less<G4String> >::const_iterator FlukaCompoundsIterator;
+
+class FlukaCompound {
+public:
+
+  //Constructor
+  FlukaCompound(const G4String& name, G4double density, G4int nElem);
+  virtual ~FlukaCompound();
+
+  //Getters
+  // * Fluka name
+  G4String GetName() const {return fFlukaMaterial->GetName();}
+  // * Real name if the fluka name is duplicated
+  G4String GetRealName() const {return fFlukaMaterial->GetRealName();}
+  // * Density
+  G4double GetDensity() const {return fFlukaMaterial->GetDensity();}
+  // * Index (it comes from the material)
+  G4int    GetIndex() const {return fFlukaMaterial->GetIndex();}
+  // * Number of materials, and index and fraction for each material
+  G4int    GetNMaterials() const {return fNMaterials;}
+  G4int    GetMaterialIndex(G4int i) const {return fElIndex[i];}
+  G4double GetMaterialFraction(G4int i) const {return fFraction[i];}
+
+  // * Associated material
+  const FlukaMaterial* GetFlukaMaterial() const {return fFlukaMaterial;}
+  FlukaMaterial* GetFlukaMaterial() {return fFlukaMaterial;}
+
+  //Setters
+  void SetName(const G4String& n) {fFlukaMaterial->SetName(n);}
+  void SetDensity(G4double d) {fFlukaMaterial->SetDensity(d);}
+
+  //Other
+  void AddElement(G4int index, G4double fraction);
+
+  //Static
+  static inline const FlukaCompoundsTable* GetCompoundTable();
+  static inline const FlukaCompound* GetFlukaCompound(const G4String& name);
+  static G4std::ostream& PrintCompounds(G4std::ostream& os);
+
+public:
+  G4int     fNMaterials; //Number of elements in total
+  G4int     fNAdded;     //Number of elements added
+  G4int*    fElIndex;    //Array of element indices
+  G4double* fFraction;   //Array of element fracions
+
+  FlukaMaterial* fFlukaMaterial; //Each compound has a "dummy" mat associated
+
+  static   FlukaCompoundsTable fFlukaCompounds;
+
+};
+
+inline const FlukaCompoundsTable* FlukaCompound::GetCompoundTable() {
+  return &fFlukaCompounds;
+}
+
+inline const FlukaCompound* FlukaCompound::GetFlukaCompound(const G4String& name) { 
+  return fFlukaCompounds[name];
+}
+
+
+//Ostream operator
+G4std::ostream& operator<<(G4std::ostream& os, const FlukaCompound& flucomp);
+
+#endif
+
diff --git a/Flugg/FlukaLowMat.cxx b/Flugg/FlukaLowMat.cxx
new file mode 100644 (file)
index 0000000..52a060f
--- /dev/null
@@ -0,0 +1,22 @@
+#include "FlukaLowMat.hh"
+#include "WrapUtils.hh"
+
+FlukaLowMat::FlukaLowMat(const G4String& name, FlukaMaterial* fmat):
+  fName(name),
+  fFlukaMaterial(fmat){
+}
+
+G4std::ostream& operator<<(G4std::ostream& os, const FlukaLowMat& flowmat){
+  os << setw10 << "LOW-MAT   ";
+  os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
+  os << setw10 << setfixed << G4std::setprecision(1) 
+     << G4double(flowmat.GetIndex());
+  os << setw10 << " " 
+     << setw10 << " " 
+     << setw10 << " " 
+     << setw10 << " " 
+     << setw10 << " ";
+  os << flowmat.GetFlukaMaterial()->GetName() << G4endl;
+
+  return os;
+}
diff --git a/Flugg/FlukaLowMat.hh b/Flugg/FlukaLowMat.hh
new file mode 100644 (file)
index 0000000..88630f7
--- /dev/null
@@ -0,0 +1,33 @@
+#ifndef FLUKALOWMAT_HH
+#define FLUKALOWMAT_HH 1
+
+#include "G4String.hh"
+#include "FlukaMaterial.hh"
+
+class FlukaLowMat {
+public:
+
+  //Constructor
+  FlukaLowMat(const G4String& name, FlukaMaterial*);
+
+  //Getters
+  G4String GetName() const {return fName;}
+  inline G4int GetIndex() const;
+  const FlukaMaterial* GetFlukaMaterial() const {return fFlukaMaterial;}
+
+  //Setters
+  void SetName(const G4String& n) {fName = n;}
+
+public:
+  G4String fName;
+  FlukaMaterial* fFlukaMaterial;
+
+};
+
+
+inline G4int FlukaLowMat::GetIndex() const {return fFlukaMaterial->GetIndex();}
+
+G4std::ostream& operator<<(G4std::ostream& os, const FlukaLowMat& flowmat);
+
+#endif
+
diff --git a/Flugg/FlukaMaterial.cxx b/Flugg/FlukaMaterial.cxx
new file mode 100644 (file)
index 0000000..405a910
--- /dev/null
@@ -0,0 +1,142 @@
+#include "FlukaMaterial.hh"
+#include "FlukaLowMat.hh"
+#include "WrapUtils.hh"
+#include "G4ios.hh"
+
+FlukaMaterialsTable FlukaMaterial::fFlukaMaterials;
+FlukaMaterialsIndexTable FlukaMaterial::fFlukaIndexMaterials;
+
+FlukaMaterial::FlukaMaterial(const G4String& name, 
+                            G4int Z, G4double A, 
+                            G4double density,
+                            G4int N):
+  fName(name),
+  fZ(Z),
+  fA(A),
+  fDensity(density),
+  fN(N),
+  fFlukaLowMat(0) {
+
+  G4String testname(name);
+  G4int matrep = 1;
+  while (fFlukaMaterials[testname] && matrep < 100) {
+    matrep++;
+    char smatrep[3];
+    sprintf(smatrep,"%.2d",matrep);
+
+    testname = name;
+    if (testname.length() <= 6)
+      testname += smatrep;
+    else
+      testname.replace(6,testname.length()-6, smatrep, 2);
+
+#ifdef G4GEOMETRY_DEBUG
+    G4cout << "INFO: We found material \'" << name << " previously defined."
+          << G4endl;
+    G4cout << "         Checking if \'" << testname << "\' exists." << G4endl;
+#endif
+  }
+
+  if (matrep > 99) {
+    G4cerr << "ERROR: Too many materials with the same name. Exiting!"
+         << G4endl;
+    abort();
+  }
+  
+  fFlukaMaterials[testname] = this;
+  if (name != testname)
+    AddLowMat(testname);
+  fIndex = fFlukaMaterials.size() + 2;
+  fFlukaIndexMaterials[fIndex] = this;
+}
+
+FlukaMaterial::~FlukaMaterial() {
+  delete fFlukaLowMat;
+}
+
+void FlukaMaterial::AddLowMat(const G4String& name) {
+  fFlukaLowMat = new FlukaLowMat(name, this);
+}
+
+
+G4String FlukaMaterial::GetRealName() const {
+  if (fFlukaLowMat)
+    return fFlukaLowMat->GetName();
+  return GetName();
+}
+
+G4std::ostream& FlukaMaterial::PrintMaterialsByName(G4std::ostream& os) {
+  PrintHeader(os, "MATERIALS");
+  for (FlukaMaterialsIterator i = fFlukaMaterials.begin(); 
+       i != fFlukaMaterials.end(); 
+       i++) {
+    FlukaMaterial* flumat     = (*i).second;
+    
+    if (flumat->GetZ()) //Skip materials that describe only compounds
+      os << *flumat;
+  }
+  return os;
+}
+
+G4std::ostream& FlukaMaterial::PrintMaterialsByIndex(G4std::ostream& os) {
+  PrintHeader(os, "MATERIALS");
+  for (FlukaMaterialsIndexIterator i = fFlukaIndexMaterials.begin(); 
+       i != fFlukaIndexMaterials.end(); 
+       i++) {
+    FlukaMaterial* flumat     = (*i).second;
+    
+    if (flumat->GetZ()) //Skip materials that describe only compounds
+      os << *flumat;
+  }
+  return os;
+}
+
+G4std::ostream& operator<<(G4std::ostream& os, const FlukaMaterial& material){
+  os << setw10 << "MATERIAL  ";
+
+  os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
+  G4double Z = G4double(material.GetZ());
+  if (Z <= 0)
+    os << setw10 << " ";
+  else
+    os << setw10 
+       << setfixed
+       << G4std::setprecision(1) 
+       << Z;
+  
+  G4double A = material.GetA();
+  if (A <= 0)
+    os << setw10 << " ";
+  else
+    os << setw10 << G4std::setprecision(3)
+       << A;
+
+  G4double density = material.GetDensity();
+  if (density <=0)
+    density = 0.999;
+  os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
+  os << setw10 
+     << setscientific
+     << G4std::setprecision(3) 
+     << density;
+
+  os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
+  os << setw10 
+     << setfixed
+     << G4std::setprecision(1) 
+     << G4double(material.GetIndex());
+
+  os << setw10 << " ";
+  if (material.GetN())
+    os << setw10 << G4double(material.GetN());
+  else
+    os << setw10 << " ";
+
+
+  os << material.GetRealName().substr(0,8) << G4endl;
+
+  if (material.GetLowMat() && material.GetZ() != 0)
+    os << *(material.GetLowMat());
+
+  return os;
+}
diff --git a/Flugg/FlukaMaterial.hh b/Flugg/FlukaMaterial.hh
new file mode 100644 (file)
index 0000000..625e1b7
--- /dev/null
@@ -0,0 +1,74 @@
+#ifndef FLUKAMATERIAL_HH
+#define FLUKAMATERIAL_HH 1
+
+#include "G4String.hh"
+
+class FlukaLowMat;
+
+#include <map>
+
+class FlukaMaterial;
+typedef G4std::map<G4String, FlukaMaterial*, G4std::less<G4String> > FlukaMaterialsTable;
+typedef G4std::map<G4String, FlukaMaterial*, G4std::less<G4String> >::const_iterator FlukaMaterialsIterator;
+typedef G4std::map<G4int, FlukaMaterial*, G4std::less<G4int> > FlukaMaterialsIndexTable;
+typedef G4std::map<G4int, FlukaMaterial*, G4std::less<G4int> >::const_iterator FlukaMaterialsIndexIterator;
+
+class FlukaMaterial {
+public:
+
+  //Constructor
+  FlukaMaterial(const G4String& name, 
+               G4int Z, G4double A, 
+               G4double density,
+               G4int N=0);
+  virtual ~FlukaMaterial();
+
+  //Getters
+  G4int    GetZ() const {return fZ;}
+  G4double GetA() const {return fA;}
+  G4double GetDensity() const {return fDensity;}
+  G4int    GetN() const {return fN;}
+  G4String GetName() const {return fName;}
+  G4String GetRealName() const;
+  G4int    GetIndex() const {return fIndex;}
+  const FlukaLowMat* GetLowMat() const {return fFlukaLowMat;}
+
+
+  //Setters
+  void SetZ(G4int Z) {fZ = Z;}
+  void SetA(G4double A) {fA = A;}
+  void SetDensity(G4double d) {fDensity = d;}
+  void SetN(G4int N) {fN = N;}
+  void SetName(const G4String& n) {fName = n;}
+
+  //Static
+  static const FlukaMaterialsTable* GetMaterialTable() { 
+    return &fFlukaMaterials;}
+  static FlukaMaterial* GetFlukaMaterial(const G4String& name) { 
+    return fFlukaMaterials[name];}
+  static G4std::ostream& PrintMaterialsByName(G4std::ostream& os);
+  static G4std::ostream& PrintMaterialsByIndex(G4std::ostream& os);
+
+protected:
+  //Other
+  void AddLowMat(const G4String& name);
+
+public:
+  G4String fName;
+  G4int    fZ;
+  G4double fA;
+  G4double fDensity;
+  G4int    fN; //Isotope number
+  FlukaLowMat* fFlukaLowMat;
+
+  size_t   fIndex;
+
+  static   FlukaMaterialsTable fFlukaMaterials; //Map of name/material
+  static   FlukaMaterialsIndexTable fFlukaIndexMaterials; //Map index/material
+
+};
+
+G4std::ostream& operator<<(G4std::ostream& os, const FlukaMaterial& material);
+
+
+#endif
index ac0f8ca..9a0291e 100644 (file)
@@ -181,83 +181,17 @@ G4std::ostream& PrintHeader(G4std::ostream& os, const char* title) {
 }
 
 ////////////////////////////////////////////////////////////////////////
-// PrintMaterial
+// ToFlukaString converts an string to something usefull in FLUKA:
+// * Capital letters
+// * Only 8 letters
+// * Replace ' ' by '_'
 ////////////////////////////////////////////////////////////////////////
-G4std::ostream& PrintMaterial(G4std::ostream& os, const char* title,
-                      G4double Z, G4double A,
-                      G4double density,
-                      G4double index,
-                      G4double N,
-                      const char* name) {
-
-  os << setw10 << title;
-
-  os.setf(static_cast<G4std::ios::fmtflags>(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(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
-  os << setw10 
-     << setscientific
-     << G4std::setprecision(3) 
-     << density;
-
-  os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
-  os << setw10 
-     << setfixed
-     << G4std::setprecision(1) 
-     << index;
-
-  os << setw10 << " ";
-
-  if (N < 0)
-    os << setw10 << " ";
-  else
-    os << setw10 << N;
-
-  os << name << G4endl;
-
-  return os;
-}
-
-
-////////////////////////////////////////////////////////////////////////
-// PrintCompund
-////////////////////////////////////////////////////////////////////////
-G4std::ostream& PrintCompound(G4std::ostream& os, const char* title,
-                      G4int count,
-                      const char* name,
-                      G4double fraction,
-                      G4double index) {
-
-  
-  if(count==3) {
-    os << name << G4endl;
-    os << setw10 << "COMPOUND  ";
+G4String ToFlukaString(G4String str) {
+  str = str.substr(0,8);
+  str.toUpper();
+  for (unsigned int i = 0; i < str.length(); i++) {
+    if (str.at(i) == ' ')
+      str.replace(i,1,"_",1);
   }
-  
-  os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
-  os << setw10
-     << setscientific
-     << G4std::setprecision(2)
-     << fraction;
-  os.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
-  os << setw10
-     << setfixed
-     << G4std::setprecision(1)
-     << index;
-
-  return os;
+  return str;
 }
-
index 801581a..fe14ac6 100644 (file)
@@ -30,21 +30,11 @@ G4int StepAndLocation(G4ThreeVector &, const G4ThreeVector &,
 bool EqualHistories(const G4NavigationHistory*, 
                    const G4NavigationHistory*);
 
-// Commonly printed things used in FGeometryInit.cc
+// Common print utilities used in FGeometryInit.cc
 inline G4std::ostream& setw10(G4std::ostream& os) { return os << G4std::setw(10);}
 inline G4std::ostream& setscientific(G4std::ostream& os) { return os << G4std::setiosflags(G4std::ios::scientific);}
 inline G4std::ostream& setfixed(G4std::ostream& os) { return os << G4std::setiosflags(G4std::ios::fixed);}
 G4std::ostream& PrintHeader(G4std::ostream& os, const char* title);
-G4std::ostream& PrintMaterial(G4std::ostream& os, const char* title,
-                      G4double Z, G4double A,
-                      G4double density,
-                      G4double index,
-                      G4double N,
-                      const char* name);
-G4std::ostream& PrintCompound(G4std::ostream& os, const char* title,
-                      G4int count,
-                      const char* name,
-                      G4double fraction,
-                      G4double index);
+G4String ToFlukaString(G4String str);
 
 #endif
index 5767c12..dd70634 100644 (file)
@@ -1,13 +1,14 @@
 FLUGGDUMMYSRCS := G4SDManager.cxx G4VUserDetectorConstruction.cxx \
                  G4FastSimulationManager.cxx
-FLUGGDUMMYHDRS := $(filter-out source.h,$(FLUGGDUMMYSRCS:.cxx=.hh) )
+FLUGGDUMMYHDRS := $(FLUGGDUMMYSRCS:.cxx=.hh)
 FLUGGWRAPSRCS := FGeometryInit.cxx WrapDN.cxx WrapG1.cxx WrapG1RT.cxx \
                 WrapIncrHist.cxx WrapIniHist.cxx WrapInit.cxx WrapLookDB.cxx \
                 WrapLookFX.cxx WrapLookMG.cxx WrapLookZ.cxx WrapMag.cxx \
                 WrapNorml.cxx WrapReg.cxx WrapSavHist.cxx WrapUtils.cxx
 FLUGGWRAPHDRS := WrapUtils.hh Wrappers.hh NavHistWithCount.hh
-FLUGGOTHERSRCS := FluggNavigator.cxx
-FLUGGOTHERHDRS := $(filter-out source.h,$(FLUGGOTHERSRCS:.cxx=.hh) )
+FLUGGOTHERSRCS := FluggNavigator.cxx FlukaMaterial.cxx FlukaCompound.cxx \
+                 FlukaLowMat.cxx
+FLUGGOTHERHDRS := $(FLUGGOTHERSRCS:.cxx=.hh)
 
 
 # Sources