Flugg/FGeometryInit.cxx

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