- FindChargedJets() added.
[u/mrichter/AliRoot.git] / Flugg / FGeometryInit.cxx
CommitLineData
26911512 1
2// Flugg tag
3// modified 17/06/02: by I. Gonzalez. STL migration
4
36c70081 5//#include <stdio.h>
6//#include <iomanip.h>
26911512 7#include "FGeometryInit.hh"
8#include "FluggNavigator.hh"
9#include "WrapUtils.hh"
26911512 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),
36c70081 32 fTransportationManager(G4TransportationManager::GetTransportationManager()),
26911512 33 myTopNode(0),
34 ptrGeoMan(0),
35 ptrArray(0),
36 ptrTouchHist(0),
37 ptrOldNavHist(0),
38 ptrTempNavHist(0),
36c70081 39 ptrJrLtGeant(0)
40{
26911512 41
e73e0522 42#ifdef G4GEOMETRY_DEBUG
43 G4cout << "==> Flugg FGeometryInit::FGeometryInit()" << G4endl;
44 G4cout << "\t+ Changing the G4Navigator for FluggNavigator..." << G4endl;
45#endif
26911512 46 G4Navigator* actualnav = fTransportationManager->GetNavigatorForTracking();
47 if (actualnav) {
48 FluggNavigator* newnav = new FluggNavigator();
49 fTransportationManager->SetNavigatorForTracking(newnav);
50 }
51 else {
36c70081 52 G4cerr << "ERROR: Could not find the actual G4Navigator" << G4endl;
26911512 53 abort();
54 }
55
56
e73e0522 57#ifdef G4GEOMETRY_DEBUG
58 G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
59#endif
26911512 60
61}
62
63
64FGeometryInit::~FGeometryInit() {
e73e0522 65 G4cout << "==> Flugg FGeometryInit::~FGeometryInit()" << G4endl;
26911512 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!
e73e0522 75 G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
26911512 76}
77
78
79void FGeometryInit::closeGeometry() {
e73e0522 80#ifdef G4GEOMETRY_DEBUG
81 G4cout << "==> Flugg FGeometryInit::closeGeometry()" << G4endl;
82#endif
26911512 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
e73e0522 99#ifdef G4GEOMETRY_DEBUG
100 G4cout << "<== Flugg FGeometryInit::closeGeometry()" << G4endl;
101#endif
26911512 102}
103
104//*************************************************************************
105
106void FGeometryInit::InitHistArray() {
e73e0522 107#ifdef G4GEOMETRY_DEBUG
108 G4cout << "==> Flugg FGeometryInit::InitHistArray()" << G4endl;
109#endif
26911512 110 ptrArray = new G4int[1000000];
111 for(G4int i=0;i<1000000;i++)
112 ptrArray[i]=0;
e73e0522 113#ifdef G4GEOMETRY_DEBUG
114 G4cout << "<== Flugg FGeometryInit::InitHistArray()" << G4endl;
115#endif
26911512 116}
117
118
119
120//*************************************************************************
121//jrLtGeant stores all crossed lattice volume histories.
122
123void FGeometryInit::InitJrLtGeantArray() {
26911512 124#ifdef G4GEOMETRY_DEBUG
e73e0522 125 G4cout << "==> Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
26911512 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;
e73e0522 132#ifdef G4GEOMETRY_DEBUG
133 G4cout << "<== Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
134#endif
26911512 135}
136
137
138void FGeometryInit::SetLttcFlagGeant(G4int newFlagLttc) {
e73e0522 139#ifdef G4GEOMETRY_DEBUG
140 G4cout << "==> Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
141#endif
26911512 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;
e73e0522 151#ifdef G4GEOMETRY_DEBUG
152 G4cout << "<== Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
153#endif
26911512 154}
155
156void 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
167void 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
182void FGeometryInit::InitHistories() {
e73e0522 183#ifdef G4GEOMETRY_DEBUG
184 G4cout << "==> Flugg FGeometryInit::InitHistories()" << G4endl;
185#endif
26911512 186 //init utility histories with navigator history
187 ptrTouchHist =
188 fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
189 ptrTempNavHist =
190 fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
191 ptrOldNavHist = new G4TouchableHistory();
e73e0522 192#ifdef G4GEOMETRY_DEBUG
193 G4cout << "<== Flugg FGeometryInit::InitHistories()" << G4endl;
194#endif
26911512 195}
196
197void FGeometryInit::DeleteHistories() {
e73e0522 198#ifdef G4GEOMETRY_DEBUG
199 G4cout << "==> Flugg FGeometryInit::DeleteHistories()" << G4endl;
200#endif
201
26911512 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;
e73e0522 208 G4cout << "<== Flugg FGeometryInit::DeleteHistories()" << G4endl;
26911512 209#endif
26911512 210}
211
212
213void FGeometryInit::UpdateHistories(const G4NavigationHistory * history,
214 G4int flagHist) {
e73e0522 215#ifdef G4GEOMETRY_DEBUG
216 G4cout << "==> Flugg FGeometryInit::UpdateHistories()" << G4endl;
217#endif
26911512 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
e73e0522 273#ifdef G4GEOMETRY_DEBUG
274 G4cout << "<== Flugg FGeometryInit::UpdateHistories()" << G4endl;
275#endif
26911512 276}
277
278//*****************************************************************************
279
280
281void FGeometryInit::createFlukaMatFile() {
26911512 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
e73e0522 291 G4cout << "==> Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
26911512 292 G4cout << "================== FILEWR =================" << G4endl;
293#endif
294
295 //open file for output
36c70081 296 G4std::ofstream fout("flukaMat.inp");
26911512 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 << " "
e73e0522 564 << nameMat << G4endl;
26911512 565 if(treCount==2)
566 fout << setw10 << " " << setw10 << " "
e73e0522 567 << nameMat << G4endl;
26911512 568 if(treCount==3)
e73e0522 569 fout << nameMat << G4endl;
570 fout << "*" << G4endl;
26911512 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
36c70081 593 G4std::ofstream vout("Volumes_index.inp");
26911512 594
595 //... and write title
e73e0522 596 vout << "*" << G4endl;
26911512 597 vout << "******************** GEANT4 VOLUMES *******************\n";
e73e0522 598 vout << "*" << G4endl;
26911512 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 }
e73e0522 622 vout << G4endl;
26911512 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);
e73e0522 653 fout << setw10 << "0.0" << G4endl;
26911512 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);
e73e0522 666 fout << setw10 << "0.0" << G4endl;
26911512 667 }
668 }
669 else {
670 fout << setw10 << G4double(indexRegFlukaTo);
671 fout << setw10 << "0.0";
672 fout << setw10 << G4double(lastFlagField);
e73e0522 673 fout << setw10 << "0.0" << G4endl;
26911512 674 }
675
676 // begin material card
677 fout << setw10 << "ASSIGNMAT ";
36c70081 678 fout.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
26911512 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);
e73e0522 690 fout << setw10 << "0.0" << G4endl;
26911512 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";
e73e0522 704 fout << setw10 << "0.0" << G4endl;
26911512 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 << "";
36c70081 745 fout.setf(static_cast<G4std::ios::fmtflags>(0),G4std::ios::floatfield);
26911512 746 fout << setw10 << setfixed
747 << G4std::setprecision(4) << Bx
748 << setw10 << By
749 << setw10 << Bz
e73e0522 750 << G4endl;
26911512 751 } // end if magnetic field
752
753 vout.close();
754 fout.close();
755 delete [] indexMatFluka;
e73e0522 756#ifdef G4GEOMETRY_DEBUG
757 G4cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
758#endif
26911512 759}