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
  11 FGeometryInit * FGeometryInit::flagInstance=0;
  12
  13 FGeometryInit* FGeometryInit::GetInstance()  {
  14 #ifdef G4GEOMETRY_DEBUG
  15   G4cout << "==> Flugg::FGeometryInit::GetInstance(), instance="
  16          << flagInstance << G4endl;
  17 #endif
  18   if (!flagInstance)
  19     flagInstance = new FGeometryInit();
  20
  21 #ifdef G4GEOMETRY_DEBUG
  22   G4cout << "<== Flugg::FGeometryInit::GetInstance(), instance="
  23          << flagInstance << G4endl;
  24 #endif
  25   return flagInstance;
  26 }
  27
  28
  29 FGeometryInit::FGeometryInit():
  30   fDetector(0),
  31   fFieldManager(0),
  32   fTransportationManager(G4TransportationManager::GetTransportationManager()),
  33   myTopNode(0),
  34   ptrGeoMan(0),
  35   ptrArray(0),
  36   ptrTouchHist(0),
  37   ptrOldNavHist(0),
  38   ptrTempNavHist(0),
  39   ptrJrLtGeant(0)
  40 {
  41
  42 #ifdef G4GEOMETRY_DEBUG
  43   G4cout << "==> Flugg FGeometryInit::FGeometryInit()" << G4endl;
  44   G4cout << "\t+ Changing the G4Navigator for FluggNavigator..." << G4endl;
  45 #endif
  46   G4Navigator* actualnav = fTransportationManager->GetNavigatorForTracking();
  47   if (actualnav) {
  48     FluggNavigator* newnav = new FluggNavigator();
  49     fTransportationManager->SetNavigatorForTracking(newnav);
  50   }
  51   else {
  52     G4cerr << "ERROR: Could not find the actual G4Navigator" << G4endl;
  53     abort();
  54   }
  55
  56
  57 #ifdef G4GEOMETRY_DEBUG
  58   G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
  59 #endif
  60
  61 }
  62
  63
  64 FGeometryInit::~FGeometryInit() {
  65   G4cout << "==> Flugg FGeometryInit::~FGeometryInit()" << G4endl;
  66   DeleteHistories();
  67   ptrGeoMan->OpenGeometry();
  68   if (fTransportationManager)
  69     delete fTransportationManager;
  70   if (ptrJrLtGeant)
  71     delete[] ptrJrLtGeant;
  72   DelHistArray();
  73
  74   //keep ATTENTION: never delete a pointer twice!
  75   G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
  76 }
  77
  78
  79 void FGeometryInit::closeGeometry() {
  80 #ifdef G4GEOMETRY_DEBUG
  81   G4cout << "==> Flugg FGeometryInit::closeGeometry()" << G4endl;
  82 #endif
  83
  84   ptrGeoMan = G4GeometryManager::GetInstance();
  85   if (ptrGeoMan) {
  86     ptrGeoMan->OpenGeometry();
  87
  88     //true argoment allows voxel construction; if false voxels are built
  89     //only for replicated volumes
  90     ptrGeoMan->CloseGeometry(true);
  91   }
  92   else {
  93     G4cerr << "ERROR in FLUGG: Could not get G4GeometryManager instance"
  94            << G4endl;
  95     G4cerr << "                in FGeometryInit::closeGeometry. Exiting!!!"
  96            << G4endl;
  97   }
  98
  99 #ifdef G4GEOMETRY_DEBUG
 100   G4cout << "<== Flugg FGeometryInit::closeGeometry()" << G4endl;
 101 #endif
 102 }
 103
 104 //*************************************************************************
 105
 106 void FGeometryInit::InitHistArray() {
 107 #ifdef G4GEOMETRY_DEBUG
 108   G4cout << "==> Flugg FGeometryInit::InitHistArray()" << G4endl;
 109 #endif
 110   ptrArray = new G4int[1000000];
 111   for(G4int i=0;i<1000000;i++)
 112     ptrArray[i]=0;
 113 #ifdef G4GEOMETRY_DEBUG
 114   G4cout << "<== Flugg FGeometryInit::InitHistArray()" << G4endl;
 115 #endif
 116 }
 117
 118
 119
 120 //*************************************************************************
 121 //jrLtGeant stores all crossed lattice volume histories.
 122
 123 void FGeometryInit::InitJrLtGeantArray() {
 124 #ifdef G4GEOMETRY_DEBUG
 125   G4cout << "==> Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
 126   G4cout << "Initializing JrLtGeant array" << G4endl;
 127 #endif
 128   ptrJrLtGeant = new G4int[10000];
 129   for(G4int x=0;x<10000;x++)
 130     ptrJrLtGeant[x]=-1;
 131   flagLttcGeant = -1;
 132 #ifdef G4GEOMETRY_DEBUG
 133   G4cout << "<== Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
 134 #endif
 135 }
 136
 137
 138 void FGeometryInit::SetLttcFlagGeant(G4int newFlagLttc) {
 139 #ifdef G4GEOMETRY_DEBUG
 140   G4cout << "==> Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
 141 #endif
 142   // Added by A.Solodkov
 143   if (newFlagLttc >= 10000) {
 144     G4cout << "Problems in FGeometryInit::SetLttcFlagGeant" << G4endl;
 145     G4cout << "Index newFlagLttc=" << newFlagLttc << " is outside array bounds"
 146            << G4endl;
 147     G4cout << "Better to stop immediately !" << G4endl;
 148     exit(1);
 149   }
 150   flagLttcGeant = newFlagLttc;
 151 #ifdef G4GEOMETRY_DEBUG
 152   G4cout << "<== Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
 153 #endif
 154 }
 155
 156 void FGeometryInit::PrintJrLtGeant() {
 157 #ifdef G4GEOMETRY_DEBUG
 158   //G4cout << "jrLtGeant:" << G4endl;
 159   //for(G4int y=0;y<=flagLttcGeant;y++)
 160   //
 161   //     G4cout << "jrLtGeant[" << y << "]=" << ptrJrLtGeant[y] << G4endl;
 162 #endif
 163 }
 164
 165 //**************************************************************************
 166
 167 void FGeometryInit::PrintHistories() {
 168   /*
 169     #ifdef G4GEOMETRY_DEBUG
 170     G4cout << "Touch hist:" << G4endl;
 171     G4cout << *(ptrTouchHist->GetHistory()) << G4endl;
 172     G4cout << "Tmp hist:" << G4endl;
 173     G4cout << *(ptrTempNavHist->GetHistory()) << G4endl;
 174     G4cout << "Old hist:" << G4endl;
 175     G4cout << *(ptrOldNavHist->GetHistory()) << G4endl;
 176     #endif
 177  */
 178 }
 179
 180
 181
 182 void FGeometryInit::InitHistories() {
 183 #ifdef G4GEOMETRY_DEBUG
 184   G4cout << "==> Flugg FGeometryInit::InitHistories()" << G4endl;
 185 #endif
 186   //init utility histories with navigator history
 187   ptrTouchHist =
 188     fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
 189   ptrTempNavHist =
 190     fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
 191   ptrOldNavHist = new G4TouchableHistory();
 192 #ifdef G4GEOMETRY_DEBUG
 193   G4cout << "<== Flugg FGeometryInit::InitHistories()" << G4endl;
 194 #endif
 195 }
 196
 197 void FGeometryInit::DeleteHistories() {
 198 #ifdef G4GEOMETRY_DEBUG
 199   G4cout << "==> Flugg FGeometryInit::DeleteHistories()" << G4endl;
 200 #endif
 201
 202   delete ptrTouchHist;
 203   delete ptrOldNavHist;
 204   delete ptrTempNavHist;
 205
 206 #ifdef G4GEOMETRY_DEBUG
 207   G4cout << "Deleting step-history objects at end of run!" << G4endl;
 208   G4cout << "<== Flugg FGeometryInit::DeleteHistories()" << G4endl;
 209 #endif
 210 }
 211
 212
 213 void FGeometryInit::UpdateHistories(const G4NavigationHistory * history,
 214                                     G4int flagHist) {
 215 #ifdef G4GEOMETRY_DEBUG
 216   G4cout << "==> Flugg FGeometryInit::UpdateHistories()" << G4endl;
 217 #endif
 218   PrintHistories();
 219
 220 #ifdef G4GEOMETRY_DEBUG
 221   G4cout << "...updating histories!" << G4endl;
 222 #endif
 223
 224   G4VPhysicalVolume * pPhysVol = history->GetTopVolume();
 225
 226   switch (flagHist) {
 227   case 0: {
 228     //this is the case when a new history is given to the
 229     //navigator and old history has to be resetted
 230     //touchable history has not been updated jet, so:
 231
 232     ptrTouchHist->UpdateYourself(pPhysVol,history);
 233     ptrTempNavHist->UpdateYourself(pPhysVol,history);
 234     G4NavigationHistory * ptrOldNavHistNotConst =
 235       const_cast<G4NavigationHistory * >(ptrOldNavHist->GetHistory());
 236     ptrOldNavHistNotConst->Reset();
 237     ptrOldNavHistNotConst->Clear();
 238     PrintHistories();
 239     break;
 240   } //case 0
 241
 242   case 1: {
 243     //this is the case when a new history is given to the
 244     //navigator but old history has to be kept (e.g. LOOKZ
 245     //is call during an event);
 246     //touchable history has not been updated jet, so:
 247
 248     ptrTouchHist->UpdateYourself(pPhysVol,history);
 249     ptrTempNavHist->UpdateYourself(pPhysVol,history);
 250     PrintHistories();
 251     break;
 252   } //case 1
 253
 254   case 2: {
 255     //this is the case when the touchable history has been
 256     //updated by a LocateGlobalPointAndUpdateTouchable call
 257
 258     G4VPhysicalVolume * pPhysVolTemp = ptrTempNavHist->GetVolume();
 259     ptrOldNavHist->UpdateYourself(pPhysVolTemp,
 260                                   ptrTempNavHist->GetHistory());
 261
 262     ptrTempNavHist->UpdateYourself(pPhysVol,history);
 263     PrintHistories();
 264     break;
 265   } //case 2
 266
 267   default: {
 268     G4cout <<" ERROR in updating step-histories!" << G4endl;
 269     break;
 270   } //default
 271   } //switch
 272
 273 #ifdef G4GEOMETRY_DEBUG
 274   G4cout << "<== Flugg FGeometryInit::UpdateHistories()" << G4endl;
 275 #endif
 276 }
 277
 278 //*****************************************************************************
 279
 280
 281 void FGeometryInit::createFlukaMatFile() {
 282   // last modification Sara Vanini 1/III/99
 283   // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case,
 284   // according to the fluka standard. In addition,. they must be equal to the
 285   // names of the fluka materials - see fluka manual - in order that the
 286   // program load the right cross sections, and equal to the names included in
 287   // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his
 288   // own .pemf, in order to get the right cross sections loaded in memory.
 289
 290 #ifdef G4GEOMETRY_DEBUG
 291   G4cout << "==> Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
 292   G4cout << "================== FILEWR =================" << G4endl;
 293 #endif
 294
 295   //open file for output
 296   G4std::ofstream fout("flukaMat.inp");
 297
 298   //PhysicalVolumeStore, Volume and MaterialTable pointers
 299   G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
 300   unsigned int  numVol = pVolStore->size();
 301
 302   G4int* indexMatFluka = 0;
 303   static G4Material * ptrMat = 0;
 304   unsigned int x = 0;
 305
 306   //Find a volume with a material associated
 307   while(!ptrMat && x<numVol) {
 308     G4VPhysicalVolume * ptrVol = (*pVolStore)[x];
 309     G4LogicalVolume * ptrLogVol = ptrVol->GetLogicalVolume();
 310     ptrMat = ptrLogVol->GetMaterial();
 311     x++;
 312   }
 313
 314   if(ptrMat) {
 315     static const G4MaterialTable * ptrMatTab = G4Material::GetMaterialTable();
 316
 317     //number of materials, elements, variable initialisations
 318     static size_t totNumMat = G4Material::GetNumberOfMaterials();
 319 #ifdef G4GEOMETRY_DEBUG
 320     G4cout << "Number of materials: " << totNumMat << G4endl;
 321 #endif
 322
 323     static size_t totNumElem = G4Element::GetNumberOfElements();
 324 #ifdef G4GEOMETRY_DEBUG
 325     G4cout << "Number of elements: " << totNumMat << G4endl;
 326 #endif
 327
 328
 329     G4int * elemIndexInMATcard = new G4int[totNumElem];
 330     for(unsigned int t=0; t<totNumElem; t++)
 331       elemIndexInMATcard[t] = 0;
 332     static const G4IsotopeTable * ptrIsotTab = 0;
 333     G4int initIsot = 0;
 334     static size_t totNumIsot;
 335     G4int* isotIndexInMATcard = 0;
 336     G4int isotPresence = 0;
 337
 338
 339     // title
 340     PrintHeader(fout,"GEANT4 ELEMENTS AND COMPOUNDS");
 341
 342     // *** loop over G4Materials to assign Fluka index
 343     G4int indexCount = 3;
 344     indexMatFluka = new G4int[totNumMat];
 345     for(unsigned int i=0; i<totNumMat; i++) {
 346       //pointer material, state
 347       ptrMat = (*ptrMatTab)[i];
 348       G4double denMat = ptrMat->GetDensity();
 349       G4String nameMat = ptrMat->GetName();
 350       // Fluka index: bh=1, vacuum=2, others=3,4..
 351       if(denMat<=1.00e-10*g/cm3)
 352         //N.B. fluka density limit decided on XI-`98
 353         indexMatFluka[i]=2;
 354       else {
 355         indexMatFluka[i]=indexCount;
 356         indexCount++;
 357       }
 358
 359       // *** write single-element material MATERIAL card
 360       size_t numElem = ptrMat->GetNumberOfElements();
 361       if(numElem==1) {
 362         G4int index = indexMatFluka[i];
 363         const G4Element * ptrElem = ptrMat->GetElement(0);
 364         size_t indElemTab = ptrElem->GetIndex();
 365         size_t numIsot = ptrElem->GetNumberOfIsotopes();
 366         G4double A = (ptrElem->GetA())/(g);
 367         if(!numIsot) {
 368           if(index!=2 && !elemIndexInMATcard[indElemTab]) {
 369             G4double Z = ptrElem->GetZ();
 370             elemIndexInMATcard[indElemTab] = index;
 371             G4String nameEl = ptrElem->GetName();
 372             nameEl.toUpper();
 373
 374             //write on file MATERIAL card of element
 375             PrintMaterial(fout, "MATERIAL  ",
 376                           Z, A,
 377                           denMat/(g/cm3),
 378                           index,
 379                           -1,
 380                           nameEl);
 381           }
 382         }
 383         if(numIsot==1) {
 384           // G4Isotope pointer
 385           const G4Isotope* ptrIsot = ptrElem->GetIsotope(0);
 386           size_t indIsotTab = ptrIsot->GetIndex();
 387           //initialize variables
 388           if(!initIsot) {
 389             totNumIsot = ptrIsot->GetNumberOfIsotopes();
 390             ptrIsotTab = ptrIsot->GetIsotopeTable();
 391             isotIndexInMATcard = new G4int[totNumIsot];
 392             for(unsigned int t=0; t<totNumIsot; t++)
 393               isotIndexInMATcard[t] = 0;
 394             initIsot = 1;
 395           }
 396           if(!isotIndexInMATcard[indIsotTab]) {
 397             //compute physical data and counters
 398             G4int ZIs = ptrIsot->GetZ();
 399             G4double AIs = (ptrIsot->GetA())/(g);
 400             G4int NIs = ptrIsot->GetN();
 401             G4String nameIsot = ptrIsot->GetName();
 402             nameIsot.toUpper();
 403             isotIndexInMATcard[indIsotTab] = index;
 404
 405             //write on file MATERIAL card of isotope
 406             PrintMaterial(fout, "MATERIAL  ",
 407                           G4double(ZIs),
 408                           AIs,
 409                           denMat/(g/cm3),
 410                           G4double(index),
 411                           G4double(NIs),
 412                           nameIsot);
 413           }
 414         }// end if(numIsot==1)
 415       }// end if(numElem==1)
 416     }// end for loop
 417
 418
 419     // *** material definitions: elements, compound made of G4Elements
 420     // or made of G4Materials
 421     for(unsigned int j=0; j<totNumMat; j++) {
 422       //pointer material, and material data
 423       ptrMat = (*ptrMatTab)[j];
 424       size_t numElem = ptrMat->GetNumberOfElements();
 425       G4double densityMat = (ptrMat->GetDensity())/(g/cm3);
 426       G4String nameMat = ptrMat->GetName();
 427       nameMat.toUpper();
 428       isotPresence = 0;
 429
 430       //fraction vector of compounds of the material
 431       const G4double* ptrFracVect = ptrMat->GetFractionVector();
 432
 433       //loop on elements of the material
 434       for(unsigned int el=0; el<numElem; el++) {
 435         //compute physical data, initialize variables
 436         const G4Element * ptrElem = ptrMat->GetElement(el);
 437         size_t indElemTab = ptrElem->GetIndex();
 438         G4String nameElem = ptrElem->GetName();
 439         nameElem.toUpper();
 440         size_t numIsot = ptrElem->GetNumberOfIsotopes();
 441         G4double A = (ptrElem->GetA())/(g);
 442
 443         if(!numIsot) {
 444           if(!elemIndexInMATcard[indElemTab]) {
 445             G4double Z = ptrElem->GetZ();
 446             G4double density = ptrFracVect[el]*densityMat;
 447             elemIndexInMATcard[indElemTab] = indexCount;
 448
 449             //write on file MATERIAL card of element
 450             PrintMaterial(fout, "MATERIAL  ",
 451                           Z, A,
 452                           density,
 453                           G4double(indexCount),
 454                           -1,
 455                           nameElem);
 456           }
 457         }
 458
 459         else {
 460           if(numIsot>=2)
 461             isotPresence = 1;
 462                                 //loop on isotopes
 463           for(unsigned int nis=0; nis<numIsot; nis++) {
 464             // G4Isotope pointer
 465             const G4Isotope* ptrIsot = ptrElem->GetIsotope(nis);
 466             size_t indIsotTab = ptrIsot->GetIndex();
 467             //initialize variables
 468             if(!initIsot) {
 469               totNumIsot = ptrIsot->GetNumberOfIsotopes();
 470               ptrIsotTab = ptrIsot->GetIsotopeTable();
 471               isotIndexInMATcard = new G4int[totNumIsot];
 472               for(unsigned int t=0; t<totNumIsot; t++)
 473                 isotIndexInMATcard[t] = 0;
 474               initIsot = 1;
 475             }
 476             if(!isotIndexInMATcard[indIsotTab]) {
 477               //compute physical data and counters
 478               G4int ZIs = ptrIsot->GetZ();
 479               G4double AIs = (ptrIsot->GetA())/(g);
 480               G4int NIs = ptrIsot->GetN();
 481               G4String nameIsot = ptrIsot->GetName();
 482               nameIsot.toUpper();
 483               G4double* ptrRelAbVect = ptrElem->GetRelativeAbundanceVector();
 484               G4double density =
 485                 ptrFracVect[el]*densityMat*ptrRelAbVect[nis]*AIs/A;
 486               G4int index = indexCount;
 487               isotIndexInMATcard[indIsotTab] = indexCount;
 488               indexCount+=1;
 489
 490               //write on file MATERIAL card of isotope
 491               PrintMaterial(fout, "MATERIAL  ",
 492                             G4double(ZIs), AIs,
 493                             density,
 494                             G4double(index),
 495                             NIs,
 496                             nameIsot);
 497             }
 498           }
 499         }
 500
 501       }
 502
 503       if(numElem>1 || isotPresence==1) {
 504         // write MATERIAL+COMPOUND card specifing the compound
 505
 506         //flags for writing COMPOUND card
 507         G4int treCount=0;
 508
 509         //make MATERIAL card for compound, start COMPOUND card
 510         fout <<"*"<<G4endl;
 511         fout <<"*   Define GEANT4 compound " << nameMat << G4endl;
 512         PrintMaterial(fout, "MATERIAL  ",
 513                       -1, -1,
 514                       densityMat,
 515                       G4double(indexMatFluka[j]),
 516                       -1,
 517                       nameMat);
 518         fout << setw10 << "COMPOUND  ";
 519
 520
 521         //write elements in COMPOUND card
 522         for(unsigned int h=0; h<numElem; h++) {
 523           const G4Element * ptrElemMat = ptrMat->GetElement(h);
 524           size_t indexElemMat = ptrElemMat->GetIndex();
 525           size_t numIsotElem = ptrElemMat->GetNumberOfIsotopes();
 526           if(!numIsotElem) {
 527             PrintCompound(fout, "COMPOUND  ",
 528                           treCount,
 529                           nameMat,
 530                           -ptrFracVect[h],
 531                           G4double(elemIndexInMATcard[indexElemMat]));
 532
 533             if(treCount==3)
 534               treCount=0;
 535             treCount+=1;
 536           }
 537           else {
 538             G4double * ptrIsotAbbVect =
 539               ptrElemMat->GetRelativeAbundanceVector();
 540
 541             for(unsigned int iso=0; iso<numIsotElem; iso++) {
 542               const G4Isotope * ptrIsotElem =ptrElemMat->GetIsotope(iso);
 543               size_t indexIsotMat = ptrIsotElem->GetIndex();
 544               G4double isotAbundPerVol =
 545                 ptrIsotAbbVect[iso]*Avogadro*densityMat*
 546                 ptrFracVect[h]/(ptrElemMat->GetA()/(g));
 547
 548               PrintCompound(fout, "COMPOUND  ",
 549                             treCount,
 550                             nameMat,
 551                             isotAbundPerVol,
 552                             G4double(isotIndexInMATcard[indexIsotMat]));
 553               if(treCount==3)
 554                 treCount=0;
 555               treCount+=1;
 556             }
 557           }
 558         }
 559
 560         //end COMPOUND card
 561         if(treCount==1)
 562           fout << setw10 << " " << setw10 << " "
 563                << setw10 << " " << setw10 << " "
 564                << nameMat << G4endl;
 565         if(treCount==2)
 566           fout << setw10 << " " << setw10 << " "
 567                << nameMat << G4endl;
 568         if(treCount==3)
 569           fout << nameMat << G4endl;
 570         fout << "*" << G4endl;
 571       }
 572
 573
 574     } // end for loop
 575     delete elemIndexInMATcard;
 576     if(initIsot)
 577       delete isotIndexInMATcard;
 578
 579   } // end if (ptrMat)
 580
 581   // *** material-volume correspondence
 582   PrintHeader(fout,"GEANT4 MATERIAL ASSIGNMENTS");
 583
 584   //initializations
 585   G4int indexMatOld = 0;
 586   G4int indexRegFlukaFrom = 0;
 587   G4int indexRegFlukaTo = 0;
 588   G4int existTo = 0;
 589   G4int flagField = 0;
 590   G4int lastFlagField = 0;
 591
 592   //open file for volume-index correspondence
 593   G4std::ofstream vout("Volumes_index.inp");
 594
 595   //... and write title
 596   vout << "*" << G4endl;
 597   vout << "********************  GEANT4 VOLUMES *******************\n";
 598   vout << "*" << G4endl;
 599
 600   //loop over all volumes...
 601   for(unsigned int l=0;l<numVol;l++) {
 602     //index of the region
 603     G4VPhysicalVolume * ptrVol = (*pVolStore)[l];
 604     G4LogicalVolume * ptrLogVol = ptrVol->GetLogicalVolume();
 605     G4int indexRegFluka = l+1;
 606
 607
 608     //write index volume and name on file Volumes_index.inp
 609     vout.setf(G4std::ios::left,G4std::ios::adjustfield);
 610     vout << setw10 << indexRegFluka;
 611     vout << G4std::setw(20) << ptrVol->GetName() << G4std::setw(20) << "";
 612     if(ptrVol->IsReplicated()) {
 613       EAxis axis;
 614       G4int nRep;
 615       G4double width;
 616       G4double offset;
 617       G4bool consum;
 618       ptrVol->GetReplicationData(axis,nRep,width,offset,consum);
 619       vout.setf(G4std::ios::left,G4std::ios::adjustfield);
 620       vout << setw10 << "Repetion Nb: " << G4std::setw(3) << nRep;
 621     }
 622     vout << G4endl;
 623
 624     //check if Magnetic Field is present in the region
 625     G4FieldManager * pMagFieldMan = ptrLogVol->GetFieldManager();
 626     const G4Field * pMagField = 0;
 627     if(pMagFieldMan)
 628       pMagField = pMagFieldMan->GetDetectorField();
 629     if(pMagField)
 630       flagField = 1;
 631     else
 632       flagField = 0;
 633
 634
 635     //index of material in the region
 636     G4Material * ptrMat = ptrLogVol->GetMaterial();
 637     G4int indexMat;
 638     if(ptrMat) {
 639       size_t indexMatGeant = ptrMat->GetIndex();
 640       indexMat = indexMatFluka[indexMatGeant];
 641     }
 642     else
 643       indexMat = 2;
 644
 645     //if materials are repeated
 646     if(indexMat==indexMatOld && flagField==lastFlagField) {
 647       indexRegFlukaTo=indexRegFluka;
 648       existTo=1;
 649       if((l+1)==numVol) {
 650         fout << setw10 << G4double(indexRegFlukaTo);
 651         fout << setw10 << "0.0";
 652         fout << setw10 << G4double(flagField);
 653         fout << setw10 << "0.0" << G4endl;
 654       }
 655
 656     }
 657     else {
 658       //write on file ASSIGNMAT card
 659
 660       //first complete last material card
 661       if(!existTo) {
 662         if(l) {
 663           fout << setw10 << "0.0";
 664           fout << setw10 << "0.0";
 665           fout << setw10 << G4double(lastFlagField);
 666           fout << setw10 << "0.0" << G4endl;
 667         }
 668       }
 669       else {
 670         fout << setw10 <<  G4double(indexRegFlukaTo);
 671         fout << setw10 << "0.0";
 672         fout << setw10 << G4double(lastFlagField);
 673         fout << setw10 << "0.0" << G4endl;
 674       }
 675
 676       // begin material card
 677       fout << setw10 << "ASSIGNMAT ";
 678       fout.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
 679       fout << setw10 << setfixed
 680            << G4std::setprecision(1) << G4double(indexMat);
 681       fout << setw10 << G4double(indexRegFluka);
 682
 683       existTo=0;
 684       indexRegFlukaFrom=indexRegFluka;
 685
 686       if((l+1)==numVol) {
 687         fout << setw10 << "0.0";
 688         fout << setw10 << "0.0";
 689         fout << setw10 << G4double(flagField);
 690         fout << setw10 << "0.0" << G4endl;
 691       }
 692     }
 693     lastFlagField = flagField;
 694     indexMatOld = indexMat;
 695   } // end of loop
 696
 697   //assign material 1 to black-hole=n.vol+1
 698   fout << setw10 << "ASSIGNMAT ";
 699   fout << setw10 << "1.0";
 700   fout << setw10 << G4double(numVol+1);
 701   fout << setw10 << "0.0";
 702   fout << setw10 << "0.0";
 703   fout << setw10 << "0.0";
 704   fout << setw10 << "0.0" << G4endl;
 705
 706   // *** magnetic field ***
 707   if(fTransportationManager->GetFieldManager()->DoesFieldExist()) {
 708     PrintHeader(fout,"GEANT4 MAGNETIC FIELD");
 709
 710     //get magnetic field pointer
 711     const G4Field * pMagField =
 712       fTransportationManager->GetFieldManager()->GetDetectorField();
 713
 714     //if uniform magnetic field, get value
 715     G4double Bx=0.0;
 716     G4double By=0.0;
 717     G4double Bz=0.0;
 718
 719     if(pMagField) {
 720       //G4ThreeVector* Field = new G4ThreeVector(1.,2.,3.);
 721       const G4UniformMagField *pUnifMagField =
 722         dynamic_cast<const G4UniformMagField*>(pMagField);
 723       if(pUnifMagField) {
 724         G4double *pFieldValue = 0;
 725         G4double *point = new G4double[3];
 726         point[0] = 0.;
 727         point[1] = 0.;
 728         point[2] = 0.;
 729         pUnifMagField->GetFieldValue(point,pFieldValue);
 730         //non capisco perche' l'instruzione seguente non fa linkare. Indaga!!
 731         //per ora posso usare GetFieldValue: va bene lo stesso.
 732         //G4ThreeVector FieldValue = pUnifMagField->GetConstantFieldValue();
 733         Bx = pFieldValue[0];
 734         By = pFieldValue[1];
 735         Bz = pFieldValue[2];
 736       }
 737
 738     }
 739
 740     //write MGNFIELD card
 741     fout << setw10 << "MGNFIELD  ";
 742     fout << setw10 << "";
 743     fout << setw10 << "";
 744     fout << setw10 << "";
 745     fout.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
 746     fout << setw10 << setfixed
 747          << G4std::setprecision(4) << Bx
 748          << setw10 << By
 749          << setw10 << Bz
 750          << G4endl;
 751   } // end if magnetic field
 752
 753   vout.close();
 754   fout.close();
 755   delete [] indexMatFluka;
 756 #ifdef G4GEOMETRY_DEBUG
 757   G4cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
 758 #endif
 759 }