Flugg/source/Wrappers/src/FGeometryInit.cc

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