Using AliRICHv3
[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() {
e73e0522 14#ifdef G4GEOMETRY_DEBUG
15 G4cout << "==> Flugg::FGeometryInit::GetInstance(), instance="
16 << flagInstance << G4endl;
17#endif
26911512 18 if (!flagInstance)
19 flagInstance = new FGeometryInit();
20
e73e0522 21#ifdef G4GEOMETRY_DEBUG
22 G4cout << "<== Flugg::FGeometryInit::GetInstance(), instance="
23 << flagInstance << G4endl;
24#endif
26911512 25 return flagInstance;
26}
27
28
29FGeometryInit::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
e73e0522 41#ifdef G4GEOMETRY_DEBUG
42 G4cout << "==> Flugg FGeometryInit::FGeometryInit()" << G4endl;
43 G4cout << "\t+ Changing the G4Navigator for FluggNavigator..." << G4endl;
44#endif
26911512 45 G4Navigator* actualnav = fTransportationManager->GetNavigatorForTracking();
46 if (actualnav) {
47 FluggNavigator* newnav = new FluggNavigator();
48 fTransportationManager->SetNavigatorForTracking(newnav);
49 }
50 else {
e73e0522 51 cerr << "ERROR: Could not find the actual G4Navigator" << G4endl;
26911512 52 abort();
53 }
54
55
e73e0522 56#ifdef G4GEOMETRY_DEBUG
57 G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
58#endif
26911512 59
60}
61
62
63FGeometryInit::~FGeometryInit() {
e73e0522 64 G4cout << "==> Flugg FGeometryInit::~FGeometryInit()" << G4endl;
26911512 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!
e73e0522 74 G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
26911512 75}
76
77
78void FGeometryInit::closeGeometry() {
e73e0522 79#ifdef G4GEOMETRY_DEBUG
80 G4cout << "==> Flugg FGeometryInit::closeGeometry()" << G4endl;
81#endif
26911512 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
e73e0522 98#ifdef G4GEOMETRY_DEBUG
99 G4cout << "<== Flugg FGeometryInit::closeGeometry()" << G4endl;
100#endif
26911512 101}
102
103//*************************************************************************
104
105void FGeometryInit::InitHistArray() {
e73e0522 106#ifdef G4GEOMETRY_DEBUG
107 G4cout << "==> Flugg FGeometryInit::InitHistArray()" << G4endl;
108#endif
26911512 109 ptrArray = new G4int[1000000];
110 for(G4int i=0;i<1000000;i++)
111 ptrArray[i]=0;
e73e0522 112#ifdef G4GEOMETRY_DEBUG
113 G4cout << "<== Flugg FGeometryInit::InitHistArray()" << G4endl;
114#endif
26911512 115}
116
117
118
119//*************************************************************************
120//jrLtGeant stores all crossed lattice volume histories.
121
122void FGeometryInit::InitJrLtGeantArray() {
26911512 123#ifdef G4GEOMETRY_DEBUG
e73e0522 124 G4cout << "==> Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
26911512 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;
e73e0522 131#ifdef G4GEOMETRY_DEBUG
132 G4cout << "<== Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
133#endif
26911512 134}
135
136
137void FGeometryInit::SetLttcFlagGeant(G4int newFlagLttc) {
e73e0522 138#ifdef G4GEOMETRY_DEBUG
139 G4cout << "==> Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
140#endif
26911512 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;
e73e0522 150#ifdef G4GEOMETRY_DEBUG
151 G4cout << "<== Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
152#endif
26911512 153}
154
155void 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
166void 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
181void FGeometryInit::InitHistories() {
e73e0522 182#ifdef G4GEOMETRY_DEBUG
183 G4cout << "==> Flugg FGeometryInit::InitHistories()" << G4endl;
184#endif
26911512 185 //init utility histories with navigator history
186 ptrTouchHist =
187 fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
188 ptrTempNavHist =
189 fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
190 ptrOldNavHist = new G4TouchableHistory();
e73e0522 191#ifdef G4GEOMETRY_DEBUG
192 G4cout << "<== Flugg FGeometryInit::InitHistories()" << G4endl;
193#endif
26911512 194}
195
196void FGeometryInit::DeleteHistories() {
e73e0522 197#ifdef G4GEOMETRY_DEBUG
198 G4cout << "==> Flugg FGeometryInit::DeleteHistories()" << G4endl;
199#endif
200
26911512 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;
e73e0522 207 G4cout << "<== Flugg FGeometryInit::DeleteHistories()" << G4endl;
26911512 208#endif
26911512 209}
210
211
212void FGeometryInit::UpdateHistories(const G4NavigationHistory * history,
213 G4int flagHist) {
e73e0522 214#ifdef G4GEOMETRY_DEBUG
215 G4cout << "==> Flugg FGeometryInit::UpdateHistories()" << G4endl;
216#endif
26911512 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
e73e0522 272#ifdef G4GEOMETRY_DEBUG
273 G4cout << "<== Flugg FGeometryInit::UpdateHistories()" << G4endl;
274#endif
26911512 275}
276
277//*****************************************************************************
278
279
280void FGeometryInit::createFlukaMatFile() {
26911512 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
e73e0522 290 G4cout << "==> Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
26911512 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 << " "
e73e0522 563 << nameMat << G4endl;
26911512 564 if(treCount==2)
565 fout << setw10 << " " << setw10 << " "
e73e0522 566 << nameMat << G4endl;
26911512 567 if(treCount==3)
e73e0522 568 fout << nameMat << G4endl;
569 fout << "*" << G4endl;
26911512 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
e73e0522 595 vout << "*" << G4endl;
26911512 596 vout << "******************** GEANT4 VOLUMES *******************\n";
e73e0522 597 vout << "*" << G4endl;
26911512 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 }
e73e0522 621 vout << G4endl;
26911512 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);
e73e0522 652 fout << setw10 << "0.0" << G4endl;
26911512 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);
e73e0522 665 fout << setw10 << "0.0" << G4endl;
26911512 666 }
667 }
668 else {
669 fout << setw10 << G4double(indexRegFlukaTo);
670 fout << setw10 << "0.0";
671 fout << setw10 << G4double(lastFlagField);
e73e0522 672 fout << setw10 << "0.0" << G4endl;
26911512 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);
e73e0522 689 fout << setw10 << "0.0" << G4endl;
26911512 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";
e73e0522 703 fout << setw10 << "0.0" << G4endl;
26911512 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
e73e0522 749 << G4endl;
26911512 750 } // end if magnetic field
751
752 vout.close();
753 fout.close();
754 delete [] indexMatFluka;
e73e0522 755#ifdef G4GEOMETRY_DEBUG
756 G4cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
757#endif
26911512 758}