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