Cleanup of output messages...
[u/mrichter/AliRoot.git] / 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 }