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