]> git.uio.no Git - u/mrichter/AliRoot.git/blame - Flugg/FGeometryInit.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Updates needed for move to geant4.6
[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"
26d97e06 10#include "FlukaMaterial.hh"
11#include "FlukaCompound.hh"
26911512 12
13FGeometryInit * FGeometryInit::flagInstance=0;
14
15FGeometryInit* FGeometryInit::GetInstance() {
e73e0522 16#ifdef G4GEOMETRY_DEBUG
17 G4cout << "==> Flugg::FGeometryInit::GetInstance(), instance="
18 << flagInstance << G4endl;
19#endif
26911512 20 if (!flagInstance)
21 flagInstance = new FGeometryInit();
22
e73e0522 23#ifdef G4GEOMETRY_DEBUG
24 G4cout << "<== Flugg::FGeometryInit::GetInstance(), instance="
25 << flagInstance << G4endl;
26#endif
26911512 27 return flagInstance;
28}
29
30
31FGeometryInit::FGeometryInit():
32 fDetector(0),
33 fFieldManager(0),
36c70081 34 fTransportationManager(G4TransportationManager::GetTransportationManager()),
26911512 35 myTopNode(0),
36 ptrGeoMan(0),
37 ptrArray(0),
38 ptrTouchHist(0),
39 ptrOldNavHist(0),
40 ptrTempNavHist(0),
36c70081 41 ptrJrLtGeant(0)
42{
26911512 43
e73e0522 44#ifdef G4GEOMETRY_DEBUG
45 G4cout << "==> Flugg FGeometryInit::FGeometryInit()" << G4endl;
46 G4cout << "\t+ Changing the G4Navigator for FluggNavigator..." << G4endl;
47#endif
26911512 48 G4Navigator* actualnav = fTransportationManager->GetNavigatorForTracking();
49 if (actualnav) {
50 FluggNavigator* newnav = new FluggNavigator();
51 fTransportationManager->SetNavigatorForTracking(newnav);
52 }
53 else {
36c70081 54 G4cerr << "ERROR: Could not find the actual G4Navigator" << G4endl;
26911512 55 abort();
56 }
57
58
e73e0522 59#ifdef G4GEOMETRY_DEBUG
60 G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
61#endif
26911512 62
63}
64
65
66FGeometryInit::~FGeometryInit() {
e73e0522 67 G4cout << "==> Flugg FGeometryInit::~FGeometryInit()" << G4endl;
26911512 68 DeleteHistories();
69 ptrGeoMan->OpenGeometry();
70 if (fTransportationManager)
71 delete fTransportationManager;
72 if (ptrJrLtGeant)
73 delete[] ptrJrLtGeant;
74 DelHistArray();
75
76 //keep ATTENTION: never delete a pointer twice!
e73e0522 77 G4cout << "<== Flugg FGeometryInit::FGeometryInit()" << G4endl;
26911512 78}
79
0edf14e7 80void FGeometryInit::Init() {
81// Build and initialize G4 geometry
82 setDetector();
83 setMotherVolume();
84 closeGeometry();
85 InitHistories();
86 InitJrLtGeantArray();
87 InitHistArray();
88 createFlukaMatFile();
89}
90
91
26911512 92void FGeometryInit::closeGeometry() {
e73e0522 93#ifdef G4GEOMETRY_DEBUG
94 G4cout << "==> Flugg FGeometryInit::closeGeometry()" << G4endl;
95#endif
26911512 96
97 ptrGeoMan = G4GeometryManager::GetInstance();
98 if (ptrGeoMan) {
99 ptrGeoMan->OpenGeometry();
100
101 //true argoment allows voxel construction; if false voxels are built
102 //only for replicated volumes
103 ptrGeoMan->CloseGeometry(true);
104 }
105 else {
106 G4cerr << "ERROR in FLUGG: Could not get G4GeometryManager instance"
107 << G4endl;
108 G4cerr << " in FGeometryInit::closeGeometry. Exiting!!!"
109 << G4endl;
110 }
111
e73e0522 112#ifdef G4GEOMETRY_DEBUG
113 G4cout << "<== Flugg FGeometryInit::closeGeometry()" << G4endl;
114#endif
26911512 115}
116
117//*************************************************************************
118
119void FGeometryInit::InitHistArray() {
e73e0522 120#ifdef G4GEOMETRY_DEBUG
121 G4cout << "==> Flugg FGeometryInit::InitHistArray()" << G4endl;
122#endif
26911512 123 ptrArray = new G4int[1000000];
124 for(G4int i=0;i<1000000;i++)
125 ptrArray[i]=0;
e73e0522 126#ifdef G4GEOMETRY_DEBUG
127 G4cout << "<== Flugg FGeometryInit::InitHistArray()" << G4endl;
128#endif
26911512 129}
130
131
132
133//*************************************************************************
134//jrLtGeant stores all crossed lattice volume histories.
135
136void FGeometryInit::InitJrLtGeantArray() {
26911512 137#ifdef G4GEOMETRY_DEBUG
e73e0522 138 G4cout << "==> Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
26911512 139 G4cout << "Initializing JrLtGeant array" << G4endl;
140#endif
141 ptrJrLtGeant = new G4int[10000];
142 for(G4int x=0;x<10000;x++)
143 ptrJrLtGeant[x]=-1;
144 flagLttcGeant = -1;
e73e0522 145#ifdef G4GEOMETRY_DEBUG
146 G4cout << "<== Flugg FGeometryInit::InitJrLtGeantArray()" << G4endl;
147#endif
26911512 148}
149
150
151void FGeometryInit::SetLttcFlagGeant(G4int newFlagLttc) {
e73e0522 152#ifdef G4GEOMETRY_DEBUG
153 G4cout << "==> Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
154#endif
26911512 155 // Added by A.Solodkov
156 if (newFlagLttc >= 10000) {
157 G4cout << "Problems in FGeometryInit::SetLttcFlagGeant" << G4endl;
158 G4cout << "Index newFlagLttc=" << newFlagLttc << " is outside array bounds"
159 << G4endl;
160 G4cout << "Better to stop immediately !" << G4endl;
161 exit(1);
162 }
163 flagLttcGeant = newFlagLttc;
e73e0522 164#ifdef G4GEOMETRY_DEBUG
165 G4cout << "<== Flugg FGeometryInit::SetLttcFlagGeant()" << G4endl;
166#endif
26911512 167}
168
169void FGeometryInit::PrintJrLtGeant() {
170#ifdef G4GEOMETRY_DEBUG
171 //G4cout << "jrLtGeant:" << G4endl;
172 //for(G4int y=0;y<=flagLttcGeant;y++)
173 //
174 // G4cout << "jrLtGeant[" << y << "]=" << ptrJrLtGeant[y] << G4endl;
175#endif
176}
177
178//**************************************************************************
179
180void FGeometryInit::PrintHistories() {
181 /*
182 #ifdef G4GEOMETRY_DEBUG
183 G4cout << "Touch hist:" << G4endl;
184 G4cout << *(ptrTouchHist->GetHistory()) << G4endl;
185 G4cout << "Tmp hist:" << G4endl;
186 G4cout << *(ptrTempNavHist->GetHistory()) << G4endl;
187 G4cout << "Old hist:" << G4endl;
188 G4cout << *(ptrOldNavHist->GetHistory()) << G4endl;
189 #endif
190 */
191}
192
193
194
195void FGeometryInit::InitHistories() {
e73e0522 196#ifdef G4GEOMETRY_DEBUG
197 G4cout << "==> Flugg FGeometryInit::InitHistories()" << G4endl;
198#endif
26911512 199 //init utility histories with navigator history
200 ptrTouchHist =
201 fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
202 ptrTempNavHist =
203 fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
204 ptrOldNavHist = new G4TouchableHistory();
e73e0522 205#ifdef G4GEOMETRY_DEBUG
206 G4cout << "<== Flugg FGeometryInit::InitHistories()" << G4endl;
207#endif
26911512 208}
209
210void FGeometryInit::DeleteHistories() {
e73e0522 211#ifdef G4GEOMETRY_DEBUG
212 G4cout << "==> Flugg FGeometryInit::DeleteHistories()" << G4endl;
213#endif
214
26911512 215 delete ptrTouchHist;
216 delete ptrOldNavHist;
217 delete ptrTempNavHist;
218
219#ifdef G4GEOMETRY_DEBUG
220 G4cout << "Deleting step-history objects at end of run!" << G4endl;
e73e0522 221 G4cout << "<== Flugg FGeometryInit::DeleteHistories()" << G4endl;
26911512 222#endif
26911512 223}
224
225
226void FGeometryInit::UpdateHistories(const G4NavigationHistory * history,
227 G4int flagHist) {
e73e0522 228#ifdef G4GEOMETRY_DEBUG
229 G4cout << "==> Flugg FGeometryInit::UpdateHistories()" << G4endl;
230#endif
26911512 231 PrintHistories();
232
233#ifdef G4GEOMETRY_DEBUG
234 G4cout << "...updating histories!" << G4endl;
235#endif
236
237 G4VPhysicalVolume * pPhysVol = history->GetTopVolume();
238
239 switch (flagHist) {
240 case 0: {
241 //this is the case when a new history is given to the
242 //navigator and old history has to be resetted
243 //touchable history has not been updated jet, so:
244
245 ptrTouchHist->UpdateYourself(pPhysVol,history);
246 ptrTempNavHist->UpdateYourself(pPhysVol,history);
247 G4NavigationHistory * ptrOldNavHistNotConst =
248 const_cast<G4NavigationHistory * >(ptrOldNavHist->GetHistory());
249 ptrOldNavHistNotConst->Reset();
250 ptrOldNavHistNotConst->Clear();
251 PrintHistories();
252 break;
253 } //case 0
254
255 case 1: {
256 //this is the case when a new history is given to the
257 //navigator but old history has to be kept (e.g. LOOKZ
258 //is call during an event);
259 //touchable history has not been updated jet, so:
260
261 ptrTouchHist->UpdateYourself(pPhysVol,history);
262 ptrTempNavHist->UpdateYourself(pPhysVol,history);
263 PrintHistories();
264 break;
265 } //case 1
266
267 case 2: {
268 //this is the case when the touchable history has been
269 //updated by a LocateGlobalPointAndUpdateTouchable call
270
271 G4VPhysicalVolume * pPhysVolTemp = ptrTempNavHist->GetVolume();
272 ptrOldNavHist->UpdateYourself(pPhysVolTemp,
273 ptrTempNavHist->GetHistory());
274
275 ptrTempNavHist->UpdateYourself(pPhysVol,history);
276 PrintHistories();
277 break;
278 } //case 2
279
280 default: {
281 G4cout <<" ERROR in updating step-histories!" << G4endl;
282 break;
283 } //default
284 } //switch
285
e73e0522 286#ifdef G4GEOMETRY_DEBUG
287 G4cout << "<== Flugg FGeometryInit::UpdateHistories()" << G4endl;
288#endif
26911512 289}
290
291//*****************************************************************************
292
26911512 293void FGeometryInit::createFlukaMatFile() {
26911512 294 // last modification Sara Vanini 1/III/99
295 // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case,
296 // according to the fluka standard. In addition,. they must be equal to the
297 // names of the fluka materials - see fluka manual - in order that the
298 // program load the right cross sections, and equal to the names included in
299 // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his
300 // own .pemf, in order to get the right cross sections loaded in memory.
301
302#ifdef G4GEOMETRY_DEBUG
e73e0522 303 G4cout << "==> Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
26911512 304 G4cout << "================== FILEWR =================" << G4endl;
305#endif
26d97e06 306
307
308 //Regions map
309 BuildRegionsMap();
0edf14e7 310 std::ofstream vos("Volumes_index.inp");
26d97e06 311 PrintRegionsMap(vos);
312 vos.close();
313
314 //Materials and compounds
315 BuildMaterialTables();
0edf14e7 316 std::ofstream fos("flukaMat.inp");
26d97e06 317 PrintMaterialTables(fos);
318 PrintAssignmat(fos);
319 PrintMagneticField(fos);
320 fos.close();
321
322#ifdef G4GEOMETRY_DEBUG
323 G4cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
324#endif
325}
326
327////////////////////////////////////////////////////////////////////////
328//
329void FGeometryInit::BuildRegionsMap() {
330#ifdef G4GEOMETRY_DEBUG
331 G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
332#endif
333
334 //Find number of Volumes in physical volume store
26911512 335 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
26d97e06 336 unsigned int numVol = pVolStore->size();
337
338 G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore
339 << ") has " << numVol << " volumes. Iterating..."
340 << G4endl;
341
342 for(unsigned int l=0; l < numVol; l++) {
343 //Get each of the physical volumes
344 G4VPhysicalVolume * physicalvolume = (*pVolStore)[l];
345 G4int iFlukaRegion = l+1;
346 fRegionVolumeMap[physicalvolume] = iFlukaRegion;
26911512 347 }
26d97e06 348
349
350
26911512 351#ifdef G4GEOMETRY_DEBUG
26d97e06 352 G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
353#endif
354}
26911512 355
0edf14e7 356void FGeometryInit::PrintRegionsMap(std::ostream& os) {
26911512 357#ifdef G4GEOMETRY_DEBUG
26d97e06 358 G4cout << "==> Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
359#endif
26911512 360
26d97e06 361 //Print some header
362 PrintHeader(os, "GEANT4 VOLUMES");
26911512 363
26d97e06 364 //Iterate over all volumes in the map
365 for (RegionIterator i = fRegionVolumeMap.begin();
366 i != fRegionVolumeMap.end();
367 i++) {
26911512 368
26d97e06 369 //Get info in the map
370 G4VPhysicalVolume* ptrVol = (*i).first;
371 int index = (*i).second;
26911512 372
26d97e06 373 //Print index and region name in some fixed format
0edf14e7 374 os.setf(std::ios::left, std::ios::adjustfield);
26d97e06 375 os << setw10 << index;
0edf14e7 376 os << std::setw(20) << ptrVol->GetName() << std::setw(20) << "";
26d97e06 377
378 //If volume is a replica... print some more stuff
26911512 379 if(ptrVol->IsReplicated()) {
380 EAxis axis;
26d97e06 381 G4int nRep = -1;
382 G4double width = -1;
383 G4double offset = -1;
384 G4bool consum = false;
385 ptrVol->GetReplicationData(axis, nRep, width, offset, consum);
0edf14e7 386 os.setf(std::ios::left, std::ios::adjustfield);
387 os << setw10 << "Repetion Nb: " << std::setw(3) << nRep;
26911512 388 }
26d97e06 389 os << G4endl;
26911512 390
26d97e06 391 }
392
393#ifdef G4GEOMETRY_DEBUG
394 G4cout << "<== Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
395#endif
396}
397
bf547b2f 398////////////////////////////////////////////////////////////////////////
1617d9fa 399//
400void FGeometryInit::BuildMediaMap()
401{
402 fRegionMediumMap = new int[fNRegions+1];
403 for (RegionIterator i = fRegionVolumeMap.begin();
404 i != fRegionVolumeMap.end();
405 i++) {
406 //Get info in the map
407 G4VPhysicalVolume* ptrVol = (*i).first;
408 int region = (*i).second;
409 G4int imed = fMediumVolumeMap[ptrVol];
410 fRegionMediumMap[region] = imed;
411 printf("BuildMediaMap %s %d %d\n",(ptrVol->GetName()).data(), region, imed);
412
413 }
414}
415
416G4int FGeometryInit::GetMedium(int region) const
417{
418 return fRegionMediumMap[region];
419}
420
421
6a53de92 422void FGeometryInit::SetMediumFromName(const char* volName, int medium, int volid)
1617d9fa 423 {
424 char name4[5];
425 char tmp[5];
426 strncpy(tmp, volName, 4);
427 tmp[4]='\0';
428 fNRegions = 0;
429
bf547b2f 430 for (RegionIterator i = fRegionVolumeMap.begin();
431 i != fRegionVolumeMap.end();
432 i++) {
1617d9fa 433 fNRegions++;
bf547b2f 434 //Get info in the map
435 G4VPhysicalVolume* ptrVol = (*i).first;
1617d9fa 436 strncpy(name4, (ptrVol->GetName()).data(), 4);
437 name4[4]='\0';
438 for (int j = 0; j < 4; j++) {
439 if (name4[j] == '\0') {
440 for (int k = j; k < 4; k++) {
441 name4[k] = ' ';
442 }
443 break;
444 }
445 }
6a53de92 446 if (! strncmp(name4, tmp, 4)) {
447 fMediumVolumeMap[ptrVol] = medium;
448 fVolIdVolumeMap[ptrVol] = volid;
449 }
bf547b2f 450 }
bf547b2f 451}
452
453
454
26d97e06 455////////////////////////////////////////////////////////////////////////
456//
457void FGeometryInit::BuildMaterialTables() {
458#ifdef G4GEOMETRY_DEBUG
459 G4cout << "==> Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
460#endif
461
462 //some terminal printout also
463 G4cout << "\t* Storing information..." << G4endl;
464
465 //The logic is the folloing:
466 //Get the Material Table and:
467 // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
468 // 2) For each single element material build a material equivalent
469 // 3) For the rest:
470 // 3.a) Build materials for each not already known element
471 // 3.b) Build the compound out of them
472
473 //Get the Material Table and iterate
474 const G4MaterialTable* matTable = G4Material::GetMaterialTable();
475 for (MatTableIterator i = matTable->begin(); i != matTable->end(); i++) {
476
477 //Get some basic material information
478 G4Material* material = (*i);
479 G4String matName = material->GetName();
480 const G4double matDensity = material->GetDensity();
481 const G4int nMatElements = material->GetNumberOfElements();
482
483 G4cout << "\t\t+ " << matName
484 << ": dens. = " << matDensity/(g/cm3) << "g/cm3"
485 << ", nElem = " << nMatElements << G4endl;
486
487 // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
488 // FlukaMaterial* is 0 in that case
489 if (matDensity <= 1.00e-10*g/cm3) {
490 G4FlukaMaterialMap[material] = 0;
491 G4cout << "\t\t Stored as vacuum" << G4endl;
26911512 492 }
26d97e06 493 // 2) For each single element material build a material equivalent
494 else if (nMatElements == 1) {
26911512 495
26d97e06 496 FlukaMaterial *flukamat =
497 BuildFlukaMaterialFromElement(material->GetElement(0),
498 matDensity);
26911512 499
26d97e06 500 G4FlukaMaterialMap[material] = flukamat;
501 G4cout << "\t\t Stored as " << flukamat->GetRealName() << G4endl;
26911512 502
26d97e06 503 } //else if (material->GetNumberOfElements() == 1)
504
505 // 3) For the rest:
506 // 3.a) Build materials for each not already known element
507 // 3.b) Build the compound out of them
508 else {
509 FlukaCompound* flukacomp =
510 BuildFlukaCompoundFromMaterial(material);
511 G4FlukaCompoundMap[material] = flukacomp;
512 G4cout << "\t\t Stored as " << flukacomp->GetRealName() << G4endl;
513 } //else for case 3)
514 } //for (materials)
515
516#ifdef G4GEOMETRY_DEBUG
517 G4cout << "<== Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
518#endif
519}
520
521FlukaMaterial*
522FGeometryInit::BuildFlukaMaterialFromElement(const G4Element* element,
523 G4double matDensity) {
524#ifdef G4GEOMETRY_DEBUG
525 G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromElement()"
526 << G4endl;
527#endif
528
529 //Get element and its properties
530 G4String elemName(ToFlukaString(element->GetName()));
531
532 FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(elemName);
533 if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
534 //Check for isotopes
535 G4int nIsotopes = element->GetNumberOfIsotopes();
536 if (nIsotopes == 0) {
537 G4double elemA = element->GetA()/g;
538 G4double elemZ = element->GetZ();
26911512 539
26d97e06 540 if (elemA != G4int(elemA) && elemZ != G4int(elemZ)) {
541 G4cout << "WARNING: Element \'" << elemName
542 << "\' has non integer Z (" << elemZ << ") or A ("
543 << elemA << ")"
544 << G4endl;
26911512 545 }
26d97e06 546
547 flukamat = new FlukaMaterial(elemName,
548 G4int(elemZ),
549 elemA,
550 matDensity/(g/cm3));
551 }
552 else if (nIsotopes == 1) {
553 const G4Isotope* isotope = element->GetIsotope(0);
554 flukamat = BuildFlukaMaterialFromIsotope(isotope, matDensity);
555 }
556 else {
557 FlukaCompound *flucomp = BuildFlukaCompoundFromElement(element,
558 matDensity);
559 flukamat = flucomp->GetFlukaMaterial();
26911512 560 }
26d97e06 561 }
562#ifdef G4GEOMETRY_DEBUG
563 else {
564 G4cout << "INFO: Element \'" << elemName
565 << "\' already exists in the DB. It will not be recreated."
566 << G4endl;
567 }
568#endif
569
570 return flukamat;
26911512 571
26d97e06 572#ifdef G4GEOMETRY_DEBUG
573 G4cout << "<== Flugg FGeometryInit::BuildFlukaMaterialFromElement()"
574 << G4endl;
575#endif
576}
577
578FlukaMaterial*
579FGeometryInit::BuildFlukaMaterialFromIsotope(const G4Isotope* isotope,
580 G4double matDensity) {
581#ifdef G4GEOMETRY_DEBUG
582 G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()"
583 << G4endl;
584#endif
585 G4String isoName(ToFlukaString(isotope->GetName()));
586 FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(isoName);
587 if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
588 G4int isoZ = isotope->GetZ();
589 G4double isoA = (isotope->GetA())/(g);
590 G4int isoN = isotope->GetN();
591 flukamat = new FlukaMaterial(isoName,
592 isoZ,
593 isoA,
594 matDensity/(g/cm3),
595 isoN);
596 }
597
598 return flukamat;
599
600#ifdef G4GEOMETRY_DEBUG
601 G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()"
602 << G4endl;
603#endif
604}
605
606FlukaCompound*
607FGeometryInit::BuildFlukaCompoundFromMaterial(const G4Material* material) {
608#ifdef G4GEOMETRY_DEBUG
609 G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()"
610 << G4endl;
611#endif
612 //Material properties
613 const G4double* elemFractions = material->GetFractionVector();
614 const G4int nMatElements = material->GetNumberOfElements();
615 const G4double matDensity = material->GetDensity();
616 G4String matName(ToFlukaString(material->GetName()));
617 FlukaCompound* flukacomp = new FlukaCompound(matName, matDensity/(g/cm3),
618 nMatElements);
619 for (G4int i = 0; i < nMatElements; i++) {
620 FlukaMaterial *flukamat =
621 BuildFlukaMaterialFromElement(material->GetElement(i), 0.0);
622
623 flukacomp->AddElement(flukamat->GetIndex(), -elemFractions[i]);
624
625 } //for (elements)
626
627 return flukacomp;
628
629#ifdef G4GEOMETRY_DEBUG
630 G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()"
631 << G4endl;
632#endif
633}
634
635FlukaCompound*
636FGeometryInit::BuildFlukaCompoundFromElement(const G4Element* element,
637 G4double matDensity) {
638#ifdef G4GEOMETRY_DEBUG
639 G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromElement()"
640 << G4endl;
641#endif
642 G4int nIsotopes = element->GetNumberOfIsotopes();
643 //fraction of nb of atomes per volume (= volume fraction?)
644 const G4double* isoAbundance = element->GetRelativeAbundanceVector();
645 G4String elemName(ToFlukaString(element->GetName()));
646
647 //Material properties
648 FlukaCompound* flukacomp = new FlukaCompound(elemName, matDensity/(g/cm3),
649 nIsotopes);
650 for (G4int i = 0; i < nIsotopes; i++) {
651 FlukaMaterial *flukamat =
652 BuildFlukaMaterialFromIsotope(element->GetIsotope(i), 0.0);
653
654 flukacomp->AddElement(flukamat->GetIndex(), isoAbundance[i]);
655
656 } //for (elements)
657
658 return flukacomp;
659
660#ifdef G4GEOMETRY_DEBUG
661 G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromElement()"
662 << G4endl;
663#endif
664}
665
0edf14e7 666void FGeometryInit::PrintMaterialTables(std::ostream& os) {
26d97e06 667#ifdef G4GEOMETRY_DEBUG
668 G4cout << "==> Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
669#endif
670 //Print Header
671 PrintHeader(os, "GEANT4 MATERIALS AND COMPOUNDS");
26911512 672
26d97e06 673 //And some more stuff
674 size_t nIsotopes = G4Isotope::GetNumberOfIsotopes();
675 size_t nElements = G4Element::GetNumberOfElements();
676 size_t nMaterials = G4Material::GetNumberOfMaterials();
677
0e22711e 678 os << "* In Geant4 there are " << nMaterials << " materials" << G4endl;
679 os << "* In Geant4 there are " << nElements << " elements" << G4endl;
680 os << "* In Geant4 there are " << nIsotopes << " isotopes" << G4endl;
26d97e06 681
682 //Materials
683 G4cout << "\t* Printing FLUKA materials..." << G4endl;
684 FlukaMaterial::PrintMaterialsByIndex(os);
685 //FlukaMaterial::PrintMaterialsByName(os);
686
687 //Compounds
688 G4cout << "\t* Printing FLUKA compounds..." << G4endl;
689 FlukaCompound::PrintCompounds(os);
690
691#ifdef G4GEOMETRY_DEBUG
692 G4cout << "<== Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
693#endif
694}
695
696////////////////////////////////////////////////////////////////////////
697//
0edf14e7 698void FGeometryInit::PrintAssignmat(std::ostream& os) {
26d97e06 699#ifdef G4GEOMETRY_DEBUG
700 G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
701#endif
702
703 //Find number of Volumes in physical volume store
704 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
705 unsigned int numVol = pVolStore->size();
706
707 G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore
708 << ") has " << numVol << " volumes. " << G4endl;
709
710 G4cout << "\t* Printing ASSIGNMAT..." << G4endl;
711
712
713 PrintHeader(os,"GEANT4 MATERIAL ASSIGNMENTS");
714 for(unsigned int l=0; l < numVol; l++) {
715
716 //Get each of the physical volumes
717 G4VPhysicalVolume * physicalvol = (*pVolStore)[l];
718
719 //Get index for that volume
720 G4int iFlukaRegion = fRegionVolumeMap[physicalvol];
721
722 //Find G4 material and navigate to its fluka compound/material
723 G4LogicalVolume * logicalVol = physicalvol->GetLogicalVolume();
724 G4Material* material = logicalVol->GetMaterial();
725 G4int matIndex = 2;
726 if (G4FlukaCompoundMap[material])
727 matIndex = G4FlukaCompoundMap[material]->GetIndex();
728 if (G4FlukaMaterialMap[material])
729 matIndex = G4FlukaMaterialMap[material]->GetIndex();
730
731 //Find if there is a magnetic field in the region
732 //check if Magnetic Field is present in the region
733 G4double flagField = 0.0;
734 G4FieldManager * pMagFieldMan = logicalVol->GetFieldManager();
735 if(pMagFieldMan && pMagFieldMan->GetDetectorField())
736 flagField = 1.0;
737
738 //Print card
739 os << setw10 << "ASSIGNMAT ";
0edf14e7 740 os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
26d97e06 741 os << setw10 << setfixed << G4double(matIndex);
742 os << setw10 << setfixed << G4double(iFlukaRegion);
743 os << setw10 << "0.0";
744 os << setw10 << setfixed << flagField;
745 os << G4endl;
746 }
747
748
749
750#ifdef G4GEOMETRY_DEBUG
751 G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
752#endif
753}
754
755
0edf14e7 756void FGeometryInit::PrintMagneticField(std::ostream& os) {
26d97e06 757#ifdef G4GEOMETRY_DEBUG
758 G4cout << "==> Flugg FGeometryInit::PrintMagneticField()" << G4endl;
759#endif
760
761 G4cout << "\t* Printing Magnetic Field..." << G4endl;
762
26911512 763 if(fTransportationManager->GetFieldManager()->DoesFieldExist()) {
26911512 764
765 //get magnetic field pointer
766 const G4Field * pMagField =
767 fTransportationManager->GetFieldManager()->GetDetectorField();
768
26911512 769
770 if(pMagField) {
26d97e06 771 //Check if it can be made a uniform magnetic field
26911512 772 const G4UniformMagField *pUnifMagField =
773 dynamic_cast<const G4UniformMagField*>(pMagField);
774 if(pUnifMagField) {
26d97e06 775 G4double B[3];
776 G4double point[4]; //it is not really used
777 pUnifMagField->GetFieldValue(point,B);
778
779 //write MGNFIELD card
780 PrintHeader(os,"GEANT4 MAGNETIC FIELD");
781 os << setw10 << "MGNFIELD ";
782 os << setw10 << "";
783 os << setw10 << "";
784 os << setw10 << "";
0edf14e7 785 os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
26d97e06 786 os << setw10 << setfixed
0edf14e7 787 << std::setprecision(4) << B[0]
26d97e06 788 << setw10 << B[1]
789 << setw10 << B[2]
790 << G4endl;
791 }
792 else {
793 G4cout << "WARNING: No Uniform Magnetic Field found." << G4endl;
794 G4cout << " Manual intervention might be needed." << G4endl;
26911512 795 }
26911512 796 }
26d97e06 797 else
798 G4cout << "\t No detector field found... " << G4endl;
26911512 799 } // end if magnetic field
26d97e06 800 else
801 G4cout << "\t No field found... " << G4endl;
802
e73e0522 803#ifdef G4GEOMETRY_DEBUG
26d97e06 804 G4cout << "<== Flugg FGeometryInit::PrintMagneticField()" << G4endl;
e73e0522 805#endif
26911512 806}
6a53de92 807
808int FGeometryInit::CurrentVolID(int ir, int& copyNo)
809{
810 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
811 G4VPhysicalVolume * physicalvol = (*pVolStore)[ir- 1];
812 copyNo = physicalvol->GetCopyNo();
813 int id = fVolIdVolumeMap[physicalvol];
814 return id;
815}
816
817int FGeometryInit::CurrentVolOffID(int ir, int off, int& copyNo)
818{
819 if (off == 0) return CurrentVolID(ir, copyNo);
820
821 G4PhysicalVolumeStore* pVolStore = G4PhysicalVolumeStore::GetInstance();
822 G4VPhysicalVolume* physicalvol = (*pVolStore)[ir- 1];
823 G4VPhysicalVolume* mother = physicalvol;
824
0edf14e7 825 int index;
826//============================================================================
827 if (mother) {
828 // Check touchable depth
829 //
830 if (ptrTouchHist->GetHistoryDepth() < off) {
831 mother = 0;
832 } else {
833 // Get the off-th mother
834 index = ptrTouchHist->GetHistoryDepth() - off;
835 // in the touchable history volumes are ordered
836 // from top volume up to mother volume;
837 // the touchable volume is not in the history
838 mother = ptrTouchHist->GetHistory()->GetVolume(index);
839 }
840 }
841//============================================================================
6a53de92 842
843 int id;
844
845 if (!mother) {
846 G4cout << "Flugg FGeometryInit::CurrentVolOffID mother not found" << G4endl;
847 id = -1;
848 copyNo = -1;
849 } else {
850 copyNo = mother ->GetCopyNo();
851 id = fVolIdVolumeMap[mother];
852 }
853 return id;
854}
dc37cac6 855
856void FGeometryInit::Gmtod(double* xm, double* xd, int iflag)
857{
858// Transforms a position from the world reference frame
859// to the current volume reference frame.
860//
861// Geant3 desription:
862// ==================
863// Computes coordinates XD (in DRS)
864// from known coordinates XM in MRS
865// The local reference system can be initialized by
866// - the tracking routines and GMTOD used in GUSTEP
867// - a call to GMEDIA(XM,NUMED)
868// - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
869// (inverse routine is GDTOM)
870//
871// If IFLAG=1 convert coordinates
872// IFLAG=2 convert direction cosinus
873//
874// ---
875 FluggNavigator * ptrNavig = getNavigatorForTracking();
876 //setting variables (and dimension: Fluka uses cm.!)
877 G4ThreeVector pGlob(xm[0],xm[1],xm[2]);
1b4955bb 878 G4ThreeVector pLoc;
dc37cac6 879
880 if (iflag == 1) {
1b4955bb 881 pGlob *= 10.0; // in mm
0edf14e7 882// change because of geant 4 6.0
883// pLoc = ptrNavig->ComputeLocalPoint(pGlob);
884 pLoc = ptrNavig->GetGlobalToLocalTransform().TransformPoint(pGlob);
885
1b4955bb 886 pLoc /= 10.0; // in cm
dc37cac6 887 } else if (iflag == 2) {
888 pLoc =
889 ptrNavig->ComputeLocalAxis(pGlob);
890 } else {
891 G4cout << "Flugg FGeometryInit::Gmtod called with undefined flag" << G4endl;
892 }
dc37cac6 893 xd[0] = pLoc[0]; xd[1] = pLoc[1]; xd[2] = pLoc[2];
894}
895
896void FGeometryInit::Gdtom(double* xd, double* xm, int iflag)
897{
898// Transforms a position from the current volume reference frame
899// to the world reference frame.
900//
901// Geant3 desription:
902// ==================
903// Computes coordinates XM (Master Reference System
904// knowing the coordinates XD (Detector Ref System)
905// The local reference system can be initialized by
906// - the tracking routines and GDTOM used in GUSTEP
907// - a call to GSCMED(NLEVEL,NAMES,NUMBER)
908// (inverse routine is GMTOD)
909//
910// If IFLAG=1 convert coordinates
911// IFLAG=2 convert direction cosinus
912//
913// ---
914
915 FluggNavigator * ptrNavig = getNavigatorForTracking();
916 G4ThreeVector pLoc(xd[0],xd[1],xd[2]);
1b4955bb 917
dc37cac6 918 G4ThreeVector pGlob;
919 if (iflag == 1) {
1b4955bb 920 pLoc *= 10.0; // in mm
dc37cac6 921 pGlob = ptrNavig->GetLocalToGlobalTransform().
922 TransformPoint(pLoc);
1b4955bb 923 pGlob /= 10.0; // in cm
dc37cac6 924 } else if (iflag == 2) {
925 pGlob = ptrNavig->GetLocalToGlobalTransform().
926 TransformAxis(pLoc);
927 } else {
928 G4cout << "Flugg FGeometryInit::Gdtom called with undefined flag" << G4endl;
929 }
930
931 xm[0] = pGlob[0]; xm[1] = pGlob[1]; xm[2] = pGlob[2];
932}