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