Flugg/FGeometryInit.cxx

   1
   2 // Flugg tag
   3 // modified 17/06/02: by I. Gonzalez. STL migration
   4
   5 //#include <stdio.h>
   6 //#include <iomanip.h>
   7 #include "FGeometryInit.hh"
   8 #include "FluggNavigator.hh"
   9 #include "WrapUtils.hh"
  10 #include "FlukaMaterial.hh"
  11 #include "FlukaCompound.hh"
  12
  13 FGeometryInit * FGeometryInit::flagInstance=0;
  14
  15 FGeometryInit* FGeometryInit::GetInstance()  {
  16 #ifdef G4GEOMETRY_DEBUG
  17   G4cout << "==> Flugg::FGeometryInit::GetInstance(), instance="
  18          << flagInstance << G4endl;
  19 #endif
  20   if (!flagInstance)
  21     flagInstance = new FGeometryInit();
  22
  23 #ifdef G4GEOMETRY_DEBUG
  24   G4cout << "<== Flugg::FGeometryInit::GetInstance(), instance="
  25          << flagInstance << G4endl;
  26 #endif
  27   return flagInstance;
  28 }
  29
  30
  31 FGeometryInit::FGeometryInit():
  32   fDetector(0),
  33   fFieldManager(0),
  34   fTransportationManager(G4TransportationManager::GetTransportationManager()),
  35   myTopNode(0),
  36   ptrGeoMan(0),
  37   ptrArray(0),
  38   ptrTouchHist(0),
  39   ptrOldNavHist(0),
  40   ptrTempNavHist(0),
  41   ptrJrLtGeant(0)
  42 {
  43
  44 #ifdef G4GEOMETRY_DEBUG
  45   G4cout << "==> Flugg FGeometryInit::FGeometryInit()" << G4endl;
  46   G4cout << "\t+ Changing the G4Navigator for FluggNavigator..." << G4endl;
  47 #endif
  48   G4Navigator* actualnav = fTransportationManager->GetNavigatorForTracking();
  49   if (actualnav) {
  50     FluggNavigator* newnav = new FluggNavigator();
  51     fTransportationManager->SetNavigatorForTracking(newnav);
  52   }
  53   else {
  54     G4cerr << "ERROR: Could not find the actual G4Navigator" << G4endl;
  55     abort();
  56   }
  57
  58
  59 #ifdef G4GEOMETRY_DEBUG
  60   G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
  61 #endif
  62
  63 }
  64
  65
  66 FGeometryInit::~FGeometryInit() {
  67   G4cout << "==> Flugg FGeometryInit::~FGeometryInit()" << G4endl;
  68   DeleteHistories();
  69   ptrGeoMan->OpenGeometry();
  70   if (fTransportationManager)
  71     delete fTransportationManager;
  72   if (ptrJrLtGeant)
  73     delete[] ptrJrLtGeant;
  74   DelHistArray();
  75
  76   //keep ATTENTION: never delete a pointer twice!
  77   G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
  78 }
  79
  80 void FGeometryInit::Init() {
  81 // Build and initialize G4 geometry
  82    setDetector();
  83    setMotherVolume();
  84    closeGeometry();
  85    InitHistories();
  86    InitJrLtGeantArray();
  87    InitHistArray();
  88    createFlukaMatFile();
  89 }
  90
  91
  92 void FGeometryInit::closeGeometry() {
  93 #ifdef G4GEOMETRY_DEBUG
  94   G4cout << "==> Flugg FGeometryInit::closeGeometry()" << G4endl;
  95 #endif
  96
  97   ptrGeoMan = G4GeometryManager::GetInstance();
  98   if (ptrGeoMan) {
  99     ptrGeoMan->OpenGeometry();
 100
 101     //true argoment allows voxel construction; if false voxels are built
 102     //only for replicated volumes
 103     ptrGeoMan->CloseGeometry(true);
 104   }
 105   else {
 106     G4cerr << "ERROR in FLUGG: Could not get G4GeometryManager instance"
 107            << G4endl;
 108     G4cerr << "                in FGeometryInit::closeGeometry. Exiting!!!"
 109            << G4endl;
 110   }
 111
 112 #ifdef G4GEOMETRY_DEBUG
 113   G4cout << "<== Flugg FGeometryInit::closeGeometry()" << G4endl;
 114 #endif
 115 }
 116
 117 //*************************************************************************
 118
 119 void FGeometryInit::InitHistArray() {
 120 #ifdef G4GEOMETRY_DEBUG
 121   G4cout << "==> Flugg FGeometryInit::InitHistArray()" << G4endl;
 122 #endif
 123   ptrArray = new G4int[1000000];
 124   for(G4int i=0;i<1000000;i++)
 125     ptrArray[i]=0;
 126 #ifdef G4GEOMETRY_DEBUG
 127   G4cout << "<== Flugg FGeometryInit::InitHistArray()" << G4endl;
 128 #endif
 129 }
 130
 131
 132
 133 //*************************************************************************
 134 //jrLtGeant stores all crossed lattice volume histories.
 135
 136 void FGeometryInit::InitJrLtGeantArray() {
 137 #ifdef G4GEOMETRY_DEBUG
 138   G4cout << "==> Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
 139   G4cout << "Initializing JrLtGeant array" << G4endl;
 140 #endif
 141   ptrJrLtGeant = new G4int[10000];
 142   for(G4int x=0;x<10000;x++)
 143     ptrJrLtGeant[x]=-1;
 144   flagLttcGeant = -1;
 145 #ifdef G4GEOMETRY_DEBUG
 146   G4cout << "<== Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
 147 #endif
 148 }
 149
 150
 151 void FGeometryInit::SetLttcFlagGeant(G4int newFlagLttc) {
 152 #ifdef G4GEOMETRY_DEBUG
 153   G4cout << "==> Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
 154 #endif
 155   // Added by A.Solodkov
 156   if (newFlagLttc >= 10000) {
 157     G4cout << "Problems in FGeometryInit::SetLttcFlagGeant" << G4endl;
 158     G4cout << "Index newFlagLttc=" << newFlagLttc << " is outside array bounds"
 159            << G4endl;
 160     G4cout << "Better to stop immediately !" << G4endl;
 161     exit(1);
 162   }
 163   flagLttcGeant = newFlagLttc;
 164 #ifdef G4GEOMETRY_DEBUG
 165   G4cout << "<== Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
 166 #endif
 167 }
 168
 169 void FGeometryInit::PrintJrLtGeant() {
 170 #ifdef G4GEOMETRY_DEBUG
 171   //G4cout << "jrLtGeant:" << G4endl;
 172   //for(G4int y=0;y<=flagLttcGeant;y++)
 173   //
 174   //     G4cout << "jrLtGeant[" << y << "]=" << ptrJrLtGeant[y] << G4endl;
 175 #endif
 176 }
 177
 178 //**************************************************************************
 179
 180 void FGeometryInit::PrintHistories() {
 181   /*
 182     #ifdef G4GEOMETRY_DEBUG
 183     G4cout << "Touch hist:" << G4endl;
 184     G4cout << *(ptrTouchHist->GetHistory()) << G4endl;
 185     G4cout << "Tmp hist:" << G4endl;
 186     G4cout << *(ptrTempNavHist->GetHistory()) << G4endl;
 187     G4cout << "Old hist:" << G4endl;
 188     G4cout << *(ptrOldNavHist->GetHistory()) << G4endl;
 189     #endif
 190  */
 191 }
 192
 193
 194
 195 void FGeometryInit::InitHistories() {
 196 #ifdef G4GEOMETRY_DEBUG
 197   G4cout << "==> Flugg FGeometryInit::InitHistories()" << G4endl;
 198 #endif
 199   //init utility histories with navigator history
 200   ptrTouchHist =
 201     fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
 202   ptrTempNavHist =
 203     fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
 204   ptrOldNavHist = new G4TouchableHistory();
 205 #ifdef G4GEOMETRY_DEBUG
 206   G4cout << "<== Flugg FGeometryInit::InitHistories()" << G4endl;
 207 #endif
 208 }
 209
 210 void FGeometryInit::DeleteHistories() {
 211 #ifdef G4GEOMETRY_DEBUG
 212   G4cout << "==> Flugg FGeometryInit::DeleteHistories()" << G4endl;
 213 #endif
 214
 215   delete ptrTouchHist;
 216   delete ptrOldNavHist;
 217   delete ptrTempNavHist;
 218
 219 #ifdef G4GEOMETRY_DEBUG
 220   G4cout << "Deleting step-history objects at end of run!" << G4endl;
 221   G4cout << "<== Flugg FGeometryInit::DeleteHistories()" << G4endl;
 222 #endif
 223 }
 224
 225
 226 void FGeometryInit::UpdateHistories(const G4NavigationHistory * history,
 227                                     G4int flagHist) {
 228 #ifdef G4GEOMETRY_DEBUG
 229   G4cout << "==> Flugg FGeometryInit::UpdateHistories()" << G4endl;
 230 #endif
 231   PrintHistories();
 232
 233 #ifdef G4GEOMETRY_DEBUG
 234   G4cout << "...updating histories!" << G4endl;
 235 #endif
 236
 237   G4VPhysicalVolume * pPhysVol = history->GetTopVolume();
 238
 239   switch (flagHist) {
 240   case 0: {
 241     //this is the case when a new history is given to the
 242     //navigator and old history has to be resetted
 243     //touchable history has not been updated jet, so:
 244
 245     ptrTouchHist->UpdateYourself(pPhysVol,history);
 246     ptrTempNavHist->UpdateYourself(pPhysVol,history);
 247     G4NavigationHistory * ptrOldNavHistNotConst =
 248       const_cast<G4NavigationHistory * >(ptrOldNavHist->GetHistory());
 249     ptrOldNavHistNotConst->Reset();
 250     ptrOldNavHistNotConst->Clear();
 251     PrintHistories();
 252     break;
 253   } //case 0
 254
 255   case 1: {
 256     //this is the case when a new history is given to the
 257     //navigator but old history has to be kept (e.g. LOOKZ
 258     //is call during an event);
 259     //touchable history has not been updated jet, so:
 260
 261     ptrTouchHist->UpdateYourself(pPhysVol,history);
 262     ptrTempNavHist->UpdateYourself(pPhysVol,history);
 263     PrintHistories();
 264     break;
 265   } //case 1
 266
 267   case 2: {
 268     //this is the case when the touchable history has been
 269     //updated by a LocateGlobalPointAndUpdateTouchable call
 270
 271     G4VPhysicalVolume * pPhysVolTemp = ptrTempNavHist->GetVolume();
 272     ptrOldNavHist->UpdateYourself(pPhysVolTemp,
 273                                   ptrTempNavHist->GetHistory());
 274
 275     ptrTempNavHist->UpdateYourself(pPhysVol,history);
 276     PrintHistories();
 277     break;
 278   } //case 2
 279
 280   default: {
 281     G4cout <<" ERROR in updating step-histories!" << G4endl;
 282     break;
 283   } //default
 284   } //switch
 285
 286 #ifdef G4GEOMETRY_DEBUG
 287   G4cout << "<== Flugg FGeometryInit::UpdateHistories()" << G4endl;
 288 #endif
 289 }
 290
 291 //*****************************************************************************
 292 int FGeometryInit::GetLastMaterialIndex() const
 293 {
 294 // Get last material index as known by FLUKA
 295    const FlukaMaterialsTable *matTable = FlukaMaterial::GetMaterialTable();
 296    int matsize = matTable->size();
 297    return matsize+2;
 298 }
 299
 300 void FGeometryInit::createFlukaMatFile() {
 301   // last modification Sara Vanini 1/III/99
 302   // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case,
 303   // according to the fluka standard. In addition,. they must be equal to the
 304   // names of the fluka materials - see fluka manual - in order that the
 305   // program load the right cross sections, and equal to the names included in
 306   // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his
 307   // own .pemf, in order to get the right cross sections loaded in memory.
 308
 309 #ifdef G4GEOMETRY_DEBUG
 310   G4cout << "==> Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
 311   G4cout << "================== FILEWR =================" << G4endl;
 312 #endif
 313
 314
 315   //Regions map
 316   BuildRegionsMap();
 317   std::ofstream vos("Volumes_index.inp");
 318   PrintRegionsMap(vos);
 319   vos.close();
 320
 321   //Materials and compounds
 322   BuildMaterialTables();
 323   std::ofstream fos("flukaMat.inp");
 324   PrintMaterialTables(fos);
 325   PrintAssignmat(fos);
 326   PrintMagneticField(fos);
 327   fos.close();
 328
 329 #ifdef G4GEOMETRY_DEBUG
 330   G4cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
 331 #endif
 332 }
 333
 334 ////////////////////////////////////////////////////////////////////////
 335 //
 336 void FGeometryInit::BuildRegionsMap() {
 337 #ifdef G4GEOMETRY_DEBUG
 338   G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
 339 #endif
 340
 341   //Find number of Volumes in physical volume store
 342   G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
 343   unsigned int numVol = pVolStore->size();
 344
 345   G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore
 346          << ") has " << numVol << " volumes. Iterating..."
 347          << G4endl;
 348
 349   for(unsigned int l=0; l < numVol; l++) {
 350     //Get each of the physical volumes
 351     G4VPhysicalVolume * physicalvolume = (*pVolStore)[l];
 352     G4int iFlukaRegion = l+1;
 353     fRegionVolumeMap[physicalvolume] = iFlukaRegion;
 354   }
 355
 356
 357
 358 #ifdef G4GEOMETRY_DEBUG
 359   G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
 360 #endif
 361 }
 362
 363 void FGeometryInit::PrintRegionsMap(std::ostream& os) {
 364 #ifdef G4GEOMETRY_DEBUG
 365   G4cout << "==> Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
 366 #endif
 367
 368   //Print some header
 369   PrintHeader(os, "GEANT4 VOLUMES");
 370
 371   //Iterate over all volumes in the map
 372   for (RegionIterator i = fRegionVolumeMap.begin();
 373        i != fRegionVolumeMap.end();
 374        i++) {
 375
 376     //Get info in the map
 377     G4VPhysicalVolume* ptrVol = (*i).first;
 378     int index = (*i).second;
 379
 380     //Print index and region name in some fixed format
 381     os.setf(std::ios::left, std::ios::adjustfield);
 382     os << setw10 << index;
 383     os << std::setw(20) << ptrVol->GetName() << std::setw(20) << "";
 384
 385     //If volume is a replica... print some more stuff
 386     if(ptrVol->IsReplicated()) {
 387       EAxis axis;
 388       G4int nRep = -1;
 389       G4double width = -1;
 390       G4double offset = -1;
 391       G4bool consum = false;
 392       ptrVol->GetReplicationData(axis, nRep, width, offset, consum);
 393       os.setf(std::ios::left, std::ios::adjustfield);
 394       os << setw10 << "Repetion Nb: " << std::setw(3) << nRep;
 395     }
 396     os << G4endl;
 397
 398   }
 399
 400 #ifdef G4GEOMETRY_DEBUG
 401   G4cout << "<== Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
 402 #endif
 403 }
 404
 405 ////////////////////////////////////////////////////////////////////////
 406 //
 407 void    FGeometryInit::BuildMediaMap()
 408 {
 409     fRegionMediumMap = new int[fNRegions+1];
 410     for (RegionIterator i = fRegionVolumeMap.begin();
 411          i != fRegionVolumeMap.end();
 412          i++) {
 413         //Get info in the map
 414         G4VPhysicalVolume* ptrVol = (*i).first;
 415         int region = (*i).second;
 416         G4int imed = fMediumVolumeMap[ptrVol];
 417         fRegionMediumMap[region] = imed;
 418         printf("BuildMediaMap %s %d %d\n",(ptrVol->GetName()).data(), region, imed);
 419
 420     }
 421 }
 422
 423 G4int FGeometryInit::GetMedium(int region) const
 424 {
 425     return fRegionMediumMap[region];
 426 }
 427
 428
 429 void  FGeometryInit::SetMediumFromName(const char* volName, int medium, int volid)
 430  {
 431   char name4[5];
 432   char tmp[5];
 433   strncpy(tmp, volName, 4);
 434   tmp[4]='\0';
 435   fNRegions = 0;
 436
 437   for (RegionIterator i = fRegionVolumeMap.begin();
 438        i != fRegionVolumeMap.end();
 439        i++) {
 440     fNRegions++;
 441     //Get info in the map
 442     G4VPhysicalVolume* ptrVol = (*i).first;
 443     strncpy(name4, (ptrVol->GetName()).data(), 4);
 444     name4[4]='\0';
 445     for (int j = 0; j < 4; j++) {
 446         if (name4[j] == '\0') {
 447             for (int k = j; k < 4; k++) {
 448                 name4[k] = ' ';
 449             }
 450             break;
 451         }
 452     }
 453     if (! strncmp(name4, tmp, 4)) {
 454         fMediumVolumeMap[ptrVol] = medium;
 455         fVolIdVolumeMap[ptrVol]  = volid;
 456     }
 457   }
 458 }
 459
 460
 461
 462 ////////////////////////////////////////////////////////////////////////
 463 //
 464 void FGeometryInit::BuildMaterialTables() {
 465 #ifdef G4GEOMETRY_DEBUG
 466   G4cout << "==> Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
 467 #endif
 468
 469   //some terminal printout also
 470   G4cout << "\t* Storing information..." << G4endl;
 471
 472   //The logic is the folloing:
 473   //Get the Material Table and:
 474   // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
 475   // 2) For each single element material build a material equivalent
 476   // 3) For the rest:
 477   //   3.a) Build materials for each not already known element
 478   //   3.b) Build the compound out of them
 479
 480   //Get the Material Table and iterate
 481   const G4MaterialTable* matTable = G4Material::GetMaterialTable();
 482   for (MatTableIterator i = matTable->begin(); i != matTable->end(); i++) {
 483
 484     //Get some basic material information
 485     G4Material* material = (*i);
 486     G4String matName = material->GetName();
 487     const G4double matDensity = material->GetDensity();
 488     const G4int nMatElements  = material->GetNumberOfElements();
 489
 490     G4cout << "\t\t+ " << matName
 491            << ": dens. = " << matDensity/(g/cm3) << "g/cm3"
 492            << ", nElem = " << nMatElements << G4endl;
 493
 494     // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
 495     //    FlukaMaterial* is 0 in that case
 496     if (matDensity <= 1.00e-10*g/cm3) {
 497       G4FlukaMaterialMap[material] = 0;
 498       G4cout << "\t\t  Stored as vacuum" << G4endl;
 499     }
 500     // 2) For each single element material build a material equivalent
 501     else if (nMatElements == 1) {
 502
 503       FlukaMaterial *flukamat =
 504         BuildFlukaMaterialFromElement(material->GetElement(0),
 505                                       matDensity);
 506
 507       G4FlukaMaterialMap[material] = flukamat;
 508       G4cout << "\t\t  Stored as " << flukamat->GetRealName() << G4endl;
 509
 510     } //else if (material->GetNumberOfElements() == 1)
 511
 512     // 3) For the rest:
 513     //   3.a) Build materials for each not already known element
 514     //   3.b) Build the compound out of them
 515     else {
 516       FlukaCompound* flukacomp =
 517         BuildFlukaCompoundFromMaterial(material);
 518       G4FlukaCompoundMap[material] = flukacomp;
 519       G4cout << "\t\t  Stored as " << flukacomp->GetRealName() << G4endl;
 520     } //else for case 3)
 521   } //for (materials)
 522
 523 #ifdef G4GEOMETRY_DEBUG
 524   G4cout << "<== Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
 525 #endif
 526 }
 527
 528 FlukaMaterial*
 529 FGeometryInit::BuildFlukaMaterialFromElement(const G4Element* element,
 530                                              G4double matDensity) {
 531 #ifdef G4GEOMETRY_DEBUG
 532   G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromElement()"
 533          << G4endl;
 534 #endif
 535
 536   //Get element and its properties
 537   G4String elemName(ToFlukaString(element->GetName()));
 538
 539   FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(elemName);
 540   if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
 541     //Check for isotopes
 542     G4int nIsotopes = element->GetNumberOfIsotopes();
 543     if (nIsotopes == 0) {
 544       G4double elemA = element->GetA()/g;
 545       G4double elemZ = element->GetZ();
 546
 547       if (elemA != G4int(elemA) && elemZ != G4int(elemZ)) {
 548         G4cout << "WARNING: Element \'" << elemName
 549                << "\' has non integer Z (" << elemZ << ") or A ("
 550                << elemA << ")"
 551                << G4endl;
 552       }
 553
 554       flukamat = new FlukaMaterial(elemName,
 555                                    G4int(elemZ),
 556                                    elemA,
 557                                    matDensity/(g/cm3));
 558     }
 559     else if (nIsotopes == 1) {
 560       const G4Isotope* isotope = element->GetIsotope(0);
 561       flukamat = BuildFlukaMaterialFromIsotope(isotope, matDensity);
 562     }
 563     else {
 564       FlukaCompound *flucomp = BuildFlukaCompoundFromElement(element,
 565                                                              matDensity);
 566       flukamat = flucomp->GetFlukaMaterial();
 567     }
 568   }
 569 #ifdef G4GEOMETRY_DEBUG
 570   else {
 571     G4cout << "INFO: Element \'" << elemName
 572            << "\' already exists in the DB. It will not be recreated."
 573            << G4endl;
 574   }
 575 #endif
 576
 577   return flukamat;
 578
 579 #ifdef G4GEOMETRY_DEBUG
 580   G4cout << "<== Flugg FGeometryInit::BuildFlukaMaterialFromElement()"
 581          << G4endl;
 582 #endif
 583 }
 584
 585 FlukaMaterial*
 586 FGeometryInit::BuildFlukaMaterialFromIsotope(const G4Isotope* isotope,
 587                                              G4double matDensity) {
 588 #ifdef G4GEOMETRY_DEBUG
 589   G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()"
 590          << G4endl;
 591 #endif
 592   G4String isoName(ToFlukaString(isotope->GetName()));
 593   FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(isoName);
 594   if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
 595     G4int isoZ = isotope->GetZ();
 596     G4double isoA = (isotope->GetA())/(g);
 597     G4int isoN = isotope->GetN();
 598     flukamat = new FlukaMaterial(isoName,
 599                                  isoZ,
 600                                  isoA,
 601                                  matDensity/(g/cm3),
 602                                  isoN);
 603   }
 604
 605   return flukamat;
 606
 607 #ifdef G4GEOMETRY_DEBUG
 608   G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()"
 609          << G4endl;
 610 #endif
 611 }
 612
 613 FlukaCompound*
 614 FGeometryInit::BuildFlukaCompoundFromMaterial(const G4Material* material) {
 615 #ifdef G4GEOMETRY_DEBUG
 616   G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()"
 617          << G4endl;
 618 #endif
 619   //Material properties
 620   const G4double* elemFractions = material->GetFractionVector();
 621   const G4int nMatElements  = material->GetNumberOfElements();
 622   const G4double matDensity = material->GetDensity();
 623   G4String matName(ToFlukaString(material->GetName()));
 624   FlukaCompound* flukacomp = new FlukaCompound(matName, matDensity/(g/cm3),
 625                                                nMatElements);
 626   for (G4int i = 0; i < nMatElements; i++) {
 627     FlukaMaterial *flukamat =
 628       BuildFlukaMaterialFromElement(material->GetElement(i), 0.0);
 629
 630     flukacomp->AddElement(flukamat->GetIndex(), -elemFractions[i]);
 631
 632   } //for (elements)
 633
 634   return flukacomp;
 635
 636 #ifdef G4GEOMETRY_DEBUG
 637   G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()"
 638          << G4endl;
 639 #endif
 640 }
 641
 642 FlukaCompound*
 643 FGeometryInit::BuildFlukaCompoundFromElement(const G4Element* element,
 644                                              G4double matDensity) {
 645 #ifdef G4GEOMETRY_DEBUG
 646   G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromElement()"
 647          << G4endl;
 648 #endif
 649   G4int nIsotopes = element->GetNumberOfIsotopes();
 650   //fraction of nb of atomes per volume (= volume fraction?)
 651   const G4double* isoAbundance =  element->GetRelativeAbundanceVector();
 652   G4String elemName(ToFlukaString(element->GetName()));
 653
 654   //Material properties
 655   FlukaCompound* flukacomp = new FlukaCompound(elemName, matDensity/(g/cm3),
 656                                                nIsotopes);
 657   for (G4int i = 0; i < nIsotopes; i++) {
 658     FlukaMaterial *flukamat =
 659       BuildFlukaMaterialFromIsotope(element->GetIsotope(i), 0.0);
 660
 661     flukacomp->AddElement(flukamat->GetIndex(), isoAbundance[i]);
 662
 663   } //for (elements)
 664
 665   return flukacomp;
 666
 667 #ifdef G4GEOMETRY_DEBUG
 668   G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromElement()"
 669          << G4endl;
 670 #endif
 671 }
 672
 673 void FGeometryInit::PrintMaterialTables(std::ostream& os) {
 674 #ifdef G4GEOMETRY_DEBUG
 675   G4cout << "==> Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
 676 #endif
 677   //Print Header
 678   PrintHeader(os, "GEANT4 MATERIALS AND COMPOUNDS");
 679
 680   //And some more stuff
 681   size_t nIsotopes = G4Isotope::GetNumberOfIsotopes();
 682   size_t nElements = G4Element::GetNumberOfElements();
 683   size_t nMaterials = G4Material::GetNumberOfMaterials();
 684
 685   os << "* In Geant4 there are " << nMaterials << " materials" << G4endl;
 686   os << "* In Geant4 there are " << nElements  << " elements"  << G4endl;
 687   os << "* In Geant4 there are " << nIsotopes  << " isotopes"  << G4endl;
 688
 689   //Materials
 690   G4cout << "\t* Printing FLUKA materials..." << G4endl;
 691   FlukaMaterial::PrintMaterialsByIndex(os);
 692   //FlukaMaterial::PrintMaterialsByName(os);
 693
 694   //Compounds
 695   G4cout << "\t* Printing FLUKA compounds..." << G4endl;
 696   FlukaCompound::PrintCompounds(os);
 697
 698 #ifdef G4GEOMETRY_DEBUG
 699   G4cout << "<== Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
 700 #endif
 701 }
 702
 703 ////////////////////////////////////////////////////////////////////////
 704 //
 705 void FGeometryInit::PrintAssignmat(std::ostream& os) {
 706 #ifdef G4GEOMETRY_DEBUG
 707   G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
 708 #endif
 709
 710   //Find number of Volumes in physical volume store
 711   G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
 712   unsigned int numVol = pVolStore->size();
 713
 714   G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore
 715          << ") has " << numVol << " volumes. " << G4endl;
 716
 717   G4cout << "\t* Printing ASSIGNMAT..." << G4endl;
 718
 719
 720   PrintHeader(os,"GEANT4 MATERIAL ASSIGNMENTS");
 721   for(unsigned int l=0; l < numVol; l++) {
 722
 723     //Get each of the physical volumes
 724     G4VPhysicalVolume * physicalvol = (*pVolStore)[l];
 725
 726     //Get index for that volume
 727     G4int iFlukaRegion = fRegionVolumeMap[physicalvol];
 728
 729     //Find G4 material and navigate to its fluka compound/material
 730     G4LogicalVolume * logicalVol = physicalvol->GetLogicalVolume();
 731     G4Material* material = logicalVol->GetMaterial();
 732     G4int matIndex = 2;
 733     if (G4FlukaCompoundMap[material])
 734       matIndex = G4FlukaCompoundMap[material]->GetIndex();
 735     if (G4FlukaMaterialMap[material])
 736       matIndex = G4FlukaMaterialMap[material]->GetIndex();
 737
 738     //Find if there is a magnetic field in the region
 739     //check if Magnetic Field is present in the region
 740     G4double flagField = 0.0;
 741     G4FieldManager * pMagFieldMan = logicalVol->GetFieldManager();
 742     if(pMagFieldMan && pMagFieldMan->GetDetectorField())
 743       flagField = 1.0;
 744
 745     //Print card
 746     os << setw10 << "ASSIGNMAT ";
 747     os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
 748     os << setw10 << setfixed << G4double(matIndex);
 749     os << setw10 << setfixed << G4double(iFlukaRegion);
 750     os << setw10 << "0.0";
 751     os << setw10 << setfixed << flagField;
 752     os << G4endl;
 753   }
 754
 755
 756
 757 #ifdef G4GEOMETRY_DEBUG
 758   G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
 759 #endif
 760 }
 761
 762
 763 void FGeometryInit::PrintMagneticField(std::ostream& os) {
 764 #ifdef G4GEOMETRY_DEBUG
 765   G4cout << "==> Flugg FGeometryInit::PrintMagneticField()" << G4endl;
 766 #endif
 767
 768   G4cout << "\t* Printing Magnetic Field..." << G4endl;
 769
 770   if(fTransportationManager->GetFieldManager()->DoesFieldExist()) {
 771
 772     //get magnetic field pointer
 773     const G4Field * pMagField =
 774       fTransportationManager->GetFieldManager()->GetDetectorField();
 775
 776
 777     if(pMagField) {
 778       //Check if it can be made a uniform magnetic field
 779       const G4UniformMagField *pUnifMagField =
 780         dynamic_cast<const G4UniformMagField*>(pMagField);
 781       if(pUnifMagField) {
 782         G4double B[3];
 783         G4double point[4]; //it is not really used
 784         pUnifMagField->GetFieldValue(point,B);
 785
 786         //write MGNFIELD card
 787         PrintHeader(os,"GEANT4 MAGNETIC FIELD");
 788         os << setw10 << "MGNFIELD  ";
 789         os << setw10 << "";
 790         os << setw10 << "";
 791         os << setw10 << "";
 792         os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
 793         os << setw10 << setfixed
 794            << std::setprecision(4) << B[0]
 795            << setw10 << B[1]
 796            << setw10 << B[2]
 797            << G4endl;
 798       }
 799       else {
 800         G4cout << "WARNING: No Uniform Magnetic Field found." << G4endl;
 801         G4cout << "         Manual intervention might be needed." << G4endl;
 802       }
 803     }
 804     else
 805       G4cout << "\t  No detector field found... " << G4endl;
 806   } // end if magnetic field
 807   else
 808     G4cout << "\t  No field found... " << G4endl;
 809
 810 #ifdef G4GEOMETRY_DEBUG
 811   G4cout << "<== Flugg FGeometryInit::PrintMagneticField()" << G4endl;
 812 #endif
 813 }
 814
 815 int FGeometryInit::CurrentVolID(int ir, int& copyNo)
 816 {
 817     if (ir == 0)
 818     {
 819         copyNo = -1;
 820         return -1;
 821     }
 822
 823     G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
 824     G4VPhysicalVolume   * physicalvol = (*pVolStore)[ir- 1];
 825
 826     if (physicalvol) {
 827         copyNo =  physicalvol->GetCopyNo();
 828     } else {
 829         copyNo = -1;
 830         return -1;
 831     }
 832
 833
 834     int id = fVolIdVolumeMap[physicalvol];
 835     return id;
 836 }
 837
 838 int FGeometryInit::CurrentVolOffID(int ir, int off, int& copyNo)
 839 {
 840     if (off == 0) return CurrentVolID(ir, copyNo);
 841
 842     G4PhysicalVolumeStore* pVolStore = G4PhysicalVolumeStore::GetInstance();
 843     G4VPhysicalVolume*     physicalvol = (*pVolStore)[ir- 1];
 844     G4VPhysicalVolume*     mother = physicalvol;
 845
 846     int index;
 847 //============================================================================
 848     if (mother) {
 849   // Check touchable depth
 850   //
 851        if (ptrTouchHist->GetHistoryDepth() < off) {
 852           mother = 0;
 853        } else {
 854           // Get the off-th mother
 855           index = ptrTouchHist->GetHistoryDepth() - off;
 856           // in the touchable history volumes are ordered
 857           // from top volume up to mother volume;
 858           // the touchable volume is not in the history
 859           mother = ptrTouchHist->GetHistory()->GetVolume(index);
 860        }
 861     }
 862 //============================================================================
 863
 864     int id;
 865
 866     if (!mother) {
 867         G4cout << "Flugg FGeometryInit::CurrentVolOffID mother not found" << G4endl;
 868         id = -1;
 869         copyNo = -1;
 870     } else {
 871         copyNo =  mother ->GetCopyNo();
 872         id =  fVolIdVolumeMap[mother];
 873     }
 874     return id;
 875 }
 876
 877 void  FGeometryInit::Gmtod(double* xm, double* xd, int iflag)
 878 {
 879 // Transforms a position from the world reference frame
 880 // to the current volume reference frame.
 881 //
 882 //  Geant3 desription:
 883 //  ==================
 884 //       Computes coordinates XD (in DRS)
 885 //       from known coordinates XM in MRS
 886 //       The local reference system can be initialized by
 887 //         - the tracking routines and GMTOD used in GUSTEP
 888 //         - a call to GMEDIA(XM,NUMED)
 889 //         - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
 890 //             (inverse routine is GDTOM)
 891 //
 892 //        If IFLAG=1  convert coordinates
 893 //           IFLAG=2  convert direction cosinus
 894 //
 895 // ---
 896     FluggNavigator        * ptrNavig     = getNavigatorForTracking();
 897     //setting variables (and dimension: Fluka uses cm.!)
 898     G4ThreeVector pGlob(xm[0],xm[1],xm[2]);
 899         G4ThreeVector pLoc;
 900
 901     if (iflag == 1) {
 902         pGlob *= 10.0; // in mm
 903 // change because of geant 4 6.0
 904 //      pLoc = ptrNavig->ComputeLocalPoint(pGlob);
 905         pLoc = ptrNavig->GetGlobalToLocalTransform().TransformPoint(pGlob);
 906
 907         pLoc /= 10.0;  // in cm
 908     } else if (iflag == 2) {
 909         pLoc =
 910             ptrNavig->ComputeLocalAxis(pGlob);
 911     } else {
 912         G4cout << "Flugg FGeometryInit::Gmtod called with undefined flag" << G4endl;
 913     }
 914     xd[0] = pLoc[0]; xd[1] = pLoc[1]; xd[2] = pLoc[2];
 915 }
 916
 917 void  FGeometryInit::Gdtom(double* xd, double* xm, int iflag)
 918 {
 919 // Transforms a position from the current volume reference frame
 920 // to the world reference frame.
 921 //
 922 //  Geant3 desription:
 923 //  ==================
 924 //  Computes coordinates XM (Master Reference System
 925 //  knowing the coordinates XD (Detector Ref System)
 926 //  The local reference system can be initialized by
 927 //    - the tracking routines and GDTOM used in GUSTEP
 928 //    - a call to GSCMED(NLEVEL,NAMES,NUMBER)
 929 //        (inverse routine is GMTOD)
 930 //
 931 //   If IFLAG=1  convert coordinates
 932 //      IFLAG=2  convert direction cosinus
 933 //
 934 // ---
 935
 936     FluggNavigator        * ptrNavig     = getNavigatorForTracking();
 937     G4ThreeVector pLoc(xd[0],xd[1],xd[2]);
 938
 939     G4ThreeVector pGlob;
 940      if (iflag == 1) {
 941          pLoc  *= 10.0; // in mm
 942          pGlob = ptrNavig->GetLocalToGlobalTransform().
 943              TransformPoint(pLoc);
 944          pGlob /= 10.0; // in cm
 945      } else if (iflag == 2) {
 946          pGlob = ptrNavig->GetLocalToGlobalTransform().
 947              TransformAxis(pLoc);
 948      } else {
 949          G4cout << "Flugg FGeometryInit::Gdtom called with undefined flag" << G4endl;
 950      }
 951
 952      xm[0] = pGlob[0]; xm[1] = pGlob[1]; xm[2] = pGlob[2];
 953 }