First implementation. Needs cleanup.
[u/mrichter/AliRoot.git] / Flugg / FGeometryInit.cxx
CommitLineData
26911512 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
11FGeometryInit * FGeometryInit::flagInstance=0;
12
13FGeometryInit* 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
25FGeometryInit::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
55FGeometryInit::~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
70void 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
93void 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
106void 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
119void 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
133void 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
144void 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
159void 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
170void 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
183void 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
247void 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}