]> 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
track counter implemented in ESD reader
[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//*****************************************************************************
b9d7c32a 292int FGeometryInit::GetLastMaterialIndex() const
293{
294// Get last material index as known by FLUKA
295 const FlukaMaterialsTable *matTable = FlukaMaterial::GetMaterialTable();
296 int matsize = matTable->size();
297 return matsize+2;
298}
26911512 299
26911512 300void FGeometryInit::createFlukaMatFile() {
26911512 301 // last modification Sara Vanini 1/III/99
302 // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case,
303 // according to the fluka standard. In addition,. they must be equal to the
304 // names of the fluka materials - see fluka manual - in order that the
305 // program load the right cross sections, and equal to the names included in
306 // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his
307 // own .pemf, in order to get the right cross sections loaded in memory.
308
309#ifdef G4GEOMETRY_DEBUG
e73e0522 310 G4cout << "==> Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
26911512 311 G4cout << "================== FILEWR =================" << G4endl;
312#endif
26d97e06 313
314
315 //Regions map
316 BuildRegionsMap();
0edf14e7 317 std::ofstream vos("Volumes_index.inp");
26d97e06 318 PrintRegionsMap(vos);
319 vos.close();
320
321 //Materials and compounds
322 BuildMaterialTables();
0edf14e7 323 std::ofstream fos("flukaMat.inp");
26d97e06 324 PrintMaterialTables(fos);
325 PrintAssignmat(fos);
326 PrintMagneticField(fos);
327 fos.close();
328
329#ifdef G4GEOMETRY_DEBUG
330 G4cout << "<== Flugg FGeometryInit::createFlukaMatFile()" << G4endl;
331#endif
332}
333
334////////////////////////////////////////////////////////////////////////
335//
336void FGeometryInit::BuildRegionsMap() {
337#ifdef G4GEOMETRY_DEBUG
338 G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
339#endif
340
341 //Find number of Volumes in physical volume store
26911512 342 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
26d97e06 343 unsigned int numVol = pVolStore->size();
344
345 G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore
346 << ") has " << numVol << " volumes. Iterating..."
347 << G4endl;
348
349 for(unsigned int l=0; l < numVol; l++) {
350 //Get each of the physical volumes
351 G4VPhysicalVolume * physicalvolume = (*pVolStore)[l];
352 G4int iFlukaRegion = l+1;
353 fRegionVolumeMap[physicalvolume] = iFlukaRegion;
26911512 354 }
26d97e06 355
356
357
26911512 358#ifdef G4GEOMETRY_DEBUG
26d97e06 359 G4cout << "==> Flugg FGeometryInit::BuildRegionsMap()" << G4endl;
360#endif
361}
26911512 362
0edf14e7 363void FGeometryInit::PrintRegionsMap(std::ostream& os) {
26911512 364#ifdef G4GEOMETRY_DEBUG
26d97e06 365 G4cout << "==> Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
366#endif
26911512 367
26d97e06 368 //Print some header
369 PrintHeader(os, "GEANT4 VOLUMES");
26911512 370
26d97e06 371 //Iterate over all volumes in the map
372 for (RegionIterator i = fRegionVolumeMap.begin();
373 i != fRegionVolumeMap.end();
374 i++) {
26911512 375
26d97e06 376 //Get info in the map
377 G4VPhysicalVolume* ptrVol = (*i).first;
378 int index = (*i).second;
26911512 379
26d97e06 380 //Print index and region name in some fixed format
0edf14e7 381 os.setf(std::ios::left, std::ios::adjustfield);
26d97e06 382 os << setw10 << index;
0edf14e7 383 os << std::setw(20) << ptrVol->GetName() << std::setw(20) << "";
26d97e06 384
385 //If volume is a replica... print some more stuff
26911512 386 if(ptrVol->IsReplicated()) {
387 EAxis axis;
26d97e06 388 G4int nRep = -1;
389 G4double width = -1;
390 G4double offset = -1;
391 G4bool consum = false;
392 ptrVol->GetReplicationData(axis, nRep, width, offset, consum);
0edf14e7 393 os.setf(std::ios::left, std::ios::adjustfield);
394 os << setw10 << "Repetion Nb: " << std::setw(3) << nRep;
26911512 395 }
26d97e06 396 os << G4endl;
26911512 397
26d97e06 398 }
399
400#ifdef G4GEOMETRY_DEBUG
401 G4cout << "<== Flugg FGeometryInit::PrintRegionsMap()" << G4endl;
402#endif
403}
404
bf547b2f 405////////////////////////////////////////////////////////////////////////
1617d9fa 406//
407void FGeometryInit::BuildMediaMap()
408{
409 fRegionMediumMap = new int[fNRegions+1];
410 for (RegionIterator i = fRegionVolumeMap.begin();
411 i != fRegionVolumeMap.end();
412 i++) {
413 //Get info in the map
414 G4VPhysicalVolume* ptrVol = (*i).first;
415 int region = (*i).second;
416 G4int imed = fMediumVolumeMap[ptrVol];
417 fRegionMediumMap[region] = imed;
418 printf("BuildMediaMap %s %d %d\n",(ptrVol->GetName()).data(), region, imed);
419
420 }
421}
422
423G4int FGeometryInit::GetMedium(int region) const
424{
425 return fRegionMediumMap[region];
426}
427
428
6a53de92 429void FGeometryInit::SetMediumFromName(const char* volName, int medium, int volid)
1617d9fa 430 {
431 char name4[5];
432 char tmp[5];
433 strncpy(tmp, volName, 4);
434 tmp[4]='\0';
435 fNRegions = 0;
436
bf547b2f 437 for (RegionIterator i = fRegionVolumeMap.begin();
438 i != fRegionVolumeMap.end();
439 i++) {
1617d9fa 440 fNRegions++;
bf547b2f 441 //Get info in the map
442 G4VPhysicalVolume* ptrVol = (*i).first;
1617d9fa 443 strncpy(name4, (ptrVol->GetName()).data(), 4);
444 name4[4]='\0';
445 for (int j = 0; j < 4; j++) {
446 if (name4[j] == '\0') {
447 for (int k = j; k < 4; k++) {
448 name4[k] = ' ';
449 }
450 break;
451 }
452 }
6a53de92 453 if (! strncmp(name4, tmp, 4)) {
454 fMediumVolumeMap[ptrVol] = medium;
455 fVolIdVolumeMap[ptrVol] = volid;
456 }
bf547b2f 457 }
bf547b2f 458}
459
460
461
26d97e06 462////////////////////////////////////////////////////////////////////////
463//
464void FGeometryInit::BuildMaterialTables() {
465#ifdef G4GEOMETRY_DEBUG
466 G4cout << "==> Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
467#endif
468
469 //some terminal printout also
470 G4cout << "\t* Storing information..." << G4endl;
471
472 //The logic is the folloing:
473 //Get the Material Table and:
474 // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
475 // 2) For each single element material build a material equivalent
476 // 3) For the rest:
477 // 3.a) Build materials for each not already known element
478 // 3.b) Build the compound out of them
479
480 //Get the Material Table and iterate
481 const G4MaterialTable* matTable = G4Material::GetMaterialTable();
482 for (MatTableIterator i = matTable->begin(); i != matTable->end(); i++) {
483
484 //Get some basic material information
485 G4Material* material = (*i);
486 G4String matName = material->GetName();
487 const G4double matDensity = material->GetDensity();
488 const G4int nMatElements = material->GetNumberOfElements();
489
490 G4cout << "\t\t+ " << matName
491 << ": dens. = " << matDensity/(g/cm3) << "g/cm3"
492 << ", nElem = " << nMatElements << G4endl;
493
494 // 1) For materials with density <= 1.00e-10*g/cm3 assign vacuum
495 // FlukaMaterial* is 0 in that case
496 if (matDensity <= 1.00e-10*g/cm3) {
497 G4FlukaMaterialMap[material] = 0;
498 G4cout << "\t\t Stored as vacuum" << G4endl;
26911512 499 }
26d97e06 500 // 2) For each single element material build a material equivalent
501 else if (nMatElements == 1) {
26911512 502
26d97e06 503 FlukaMaterial *flukamat =
504 BuildFlukaMaterialFromElement(material->GetElement(0),
505 matDensity);
26911512 506
26d97e06 507 G4FlukaMaterialMap[material] = flukamat;
508 G4cout << "\t\t Stored as " << flukamat->GetRealName() << G4endl;
26911512 509
26d97e06 510 } //else if (material->GetNumberOfElements() == 1)
511
512 // 3) For the rest:
513 // 3.a) Build materials for each not already known element
514 // 3.b) Build the compound out of them
515 else {
516 FlukaCompound* flukacomp =
517 BuildFlukaCompoundFromMaterial(material);
518 G4FlukaCompoundMap[material] = flukacomp;
519 G4cout << "\t\t Stored as " << flukacomp->GetRealName() << G4endl;
520 } //else for case 3)
521 } //for (materials)
522
523#ifdef G4GEOMETRY_DEBUG
524 G4cout << "<== Flugg FGeometryInit::BuildMaterialTables()" << G4endl;
525#endif
526}
527
528FlukaMaterial*
529FGeometryInit::BuildFlukaMaterialFromElement(const G4Element* element,
530 G4double matDensity) {
531#ifdef G4GEOMETRY_DEBUG
532 G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromElement()"
533 << G4endl;
534#endif
535
536 //Get element and its properties
537 G4String elemName(ToFlukaString(element->GetName()));
538
539 FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(elemName);
540 if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
541 //Check for isotopes
542 G4int nIsotopes = element->GetNumberOfIsotopes();
543 if (nIsotopes == 0) {
544 G4double elemA = element->GetA()/g;
545 G4double elemZ = element->GetZ();
26911512 546
26d97e06 547 if (elemA != G4int(elemA) && elemZ != G4int(elemZ)) {
548 G4cout << "WARNING: Element \'" << elemName
549 << "\' has non integer Z (" << elemZ << ") or A ("
550 << elemA << ")"
551 << G4endl;
26911512 552 }
26d97e06 553
554 flukamat = new FlukaMaterial(elemName,
555 G4int(elemZ),
556 elemA,
557 matDensity/(g/cm3));
558 }
559 else if (nIsotopes == 1) {
560 const G4Isotope* isotope = element->GetIsotope(0);
561 flukamat = BuildFlukaMaterialFromIsotope(isotope, matDensity);
562 }
563 else {
564 FlukaCompound *flucomp = BuildFlukaCompoundFromElement(element,
565 matDensity);
566 flukamat = flucomp->GetFlukaMaterial();
26911512 567 }
26d97e06 568 }
569#ifdef G4GEOMETRY_DEBUG
570 else {
571 G4cout << "INFO: Element \'" << elemName
572 << "\' already exists in the DB. It will not be recreated."
573 << G4endl;
574 }
575#endif
576
577 return flukamat;
26911512 578
26d97e06 579#ifdef G4GEOMETRY_DEBUG
580 G4cout << "<== Flugg FGeometryInit::BuildFlukaMaterialFromElement()"
581 << G4endl;
582#endif
583}
584
585FlukaMaterial*
586FGeometryInit::BuildFlukaMaterialFromIsotope(const G4Isotope* isotope,
587 G4double matDensity) {
588#ifdef G4GEOMETRY_DEBUG
589 G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()"
590 << G4endl;
591#endif
592 G4String isoName(ToFlukaString(isotope->GetName()));
593 FlukaMaterial* flukamat = FlukaMaterial::GetFlukaMaterial(isoName);
594 if (matDensity != 0 || (matDensity == 0 && flukamat == 0)) {
595 G4int isoZ = isotope->GetZ();
596 G4double isoA = (isotope->GetA())/(g);
597 G4int isoN = isotope->GetN();
598 flukamat = new FlukaMaterial(isoName,
599 isoZ,
600 isoA,
601 matDensity/(g/cm3),
602 isoN);
603 }
604
605 return flukamat;
606
607#ifdef G4GEOMETRY_DEBUG
608 G4cout << "==> Flugg FGeometryInit::BuildFlukaMaterialFromIsotope()"
609 << G4endl;
610#endif
611}
612
613FlukaCompound*
614FGeometryInit::BuildFlukaCompoundFromMaterial(const G4Material* material) {
615#ifdef G4GEOMETRY_DEBUG
616 G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()"
617 << G4endl;
618#endif
619 //Material properties
620 const G4double* elemFractions = material->GetFractionVector();
621 const G4int nMatElements = material->GetNumberOfElements();
622 const G4double matDensity = material->GetDensity();
623 G4String matName(ToFlukaString(material->GetName()));
624 FlukaCompound* flukacomp = new FlukaCompound(matName, matDensity/(g/cm3),
625 nMatElements);
626 for (G4int i = 0; i < nMatElements; i++) {
627 FlukaMaterial *flukamat =
628 BuildFlukaMaterialFromElement(material->GetElement(i), 0.0);
629
630 flukacomp->AddElement(flukamat->GetIndex(), -elemFractions[i]);
631
632 } //for (elements)
633
634 return flukacomp;
635
636#ifdef G4GEOMETRY_DEBUG
637 G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromMaterial()"
638 << G4endl;
639#endif
640}
641
642FlukaCompound*
643FGeometryInit::BuildFlukaCompoundFromElement(const G4Element* element,
644 G4double matDensity) {
645#ifdef G4GEOMETRY_DEBUG
646 G4cout << "==> Flugg FGeometryInit::BuildFlukaCompoundFromElement()"
647 << G4endl;
648#endif
649 G4int nIsotopes = element->GetNumberOfIsotopes();
650 //fraction of nb of atomes per volume (= volume fraction?)
651 const G4double* isoAbundance = element->GetRelativeAbundanceVector();
652 G4String elemName(ToFlukaString(element->GetName()));
653
654 //Material properties
655 FlukaCompound* flukacomp = new FlukaCompound(elemName, matDensity/(g/cm3),
656 nIsotopes);
657 for (G4int i = 0; i < nIsotopes; i++) {
658 FlukaMaterial *flukamat =
659 BuildFlukaMaterialFromIsotope(element->GetIsotope(i), 0.0);
660
661 flukacomp->AddElement(flukamat->GetIndex(), isoAbundance[i]);
662
663 } //for (elements)
664
665 return flukacomp;
666
667#ifdef G4GEOMETRY_DEBUG
668 G4cout << "<== Flugg FGeometryInit::BuildFlukaCompoundFromElement()"
669 << G4endl;
670#endif
671}
672
0edf14e7 673void FGeometryInit::PrintMaterialTables(std::ostream& os) {
26d97e06 674#ifdef G4GEOMETRY_DEBUG
675 G4cout << "==> Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
676#endif
677 //Print Header
678 PrintHeader(os, "GEANT4 MATERIALS AND COMPOUNDS");
26911512 679
26d97e06 680 //And some more stuff
681 size_t nIsotopes = G4Isotope::GetNumberOfIsotopes();
682 size_t nElements = G4Element::GetNumberOfElements();
683 size_t nMaterials = G4Material::GetNumberOfMaterials();
684
0e22711e 685 os << "* In Geant4 there are " << nMaterials << " materials" << G4endl;
686 os << "* In Geant4 there are " << nElements << " elements" << G4endl;
687 os << "* In Geant4 there are " << nIsotopes << " isotopes" << G4endl;
26d97e06 688
689 //Materials
690 G4cout << "\t* Printing FLUKA materials..." << G4endl;
691 FlukaMaterial::PrintMaterialsByIndex(os);
692 //FlukaMaterial::PrintMaterialsByName(os);
693
694 //Compounds
695 G4cout << "\t* Printing FLUKA compounds..." << G4endl;
696 FlukaCompound::PrintCompounds(os);
697
698#ifdef G4GEOMETRY_DEBUG
699 G4cout << "<== Flugg FGeometryInit::PrintMaterialTables()" << G4endl;
700#endif
701}
702
703////////////////////////////////////////////////////////////////////////
704//
0edf14e7 705void FGeometryInit::PrintAssignmat(std::ostream& os) {
26d97e06 706#ifdef G4GEOMETRY_DEBUG
707 G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
708#endif
709
710 //Find number of Volumes in physical volume store
711 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
712 unsigned int numVol = pVolStore->size();
713
714 G4cout << "\t* G4PhysicalVolumeStore (" << pVolStore
715 << ") has " << numVol << " volumes. " << G4endl;
716
717 G4cout << "\t* Printing ASSIGNMAT..." << G4endl;
718
719
720 PrintHeader(os,"GEANT4 MATERIAL ASSIGNMENTS");
721 for(unsigned int l=0; l < numVol; l++) {
722
723 //Get each of the physical volumes
724 G4VPhysicalVolume * physicalvol = (*pVolStore)[l];
725
726 //Get index for that volume
727 G4int iFlukaRegion = fRegionVolumeMap[physicalvol];
728
729 //Find G4 material and navigate to its fluka compound/material
730 G4LogicalVolume * logicalVol = physicalvol->GetLogicalVolume();
731 G4Material* material = logicalVol->GetMaterial();
732 G4int matIndex = 2;
733 if (G4FlukaCompoundMap[material])
734 matIndex = G4FlukaCompoundMap[material]->GetIndex();
735 if (G4FlukaMaterialMap[material])
736 matIndex = G4FlukaMaterialMap[material]->GetIndex();
737
738 //Find if there is a magnetic field in the region
739 //check if Magnetic Field is present in the region
740 G4double flagField = 0.0;
741 G4FieldManager * pMagFieldMan = logicalVol->GetFieldManager();
742 if(pMagFieldMan && pMagFieldMan->GetDetectorField())
743 flagField = 1.0;
744
745 //Print card
746 os << setw10 << "ASSIGNMAT ";
0edf14e7 747 os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
26d97e06 748 os << setw10 << setfixed << G4double(matIndex);
749 os << setw10 << setfixed << G4double(iFlukaRegion);
750 os << setw10 << "0.0";
751 os << setw10 << setfixed << flagField;
752 os << G4endl;
753 }
754
755
756
757#ifdef G4GEOMETRY_DEBUG
758 G4cout << "==> Flugg FGeometryInit::PrintAssignmat()" << G4endl;
759#endif
760}
761
762
0edf14e7 763void FGeometryInit::PrintMagneticField(std::ostream& os) {
26d97e06 764#ifdef G4GEOMETRY_DEBUG
765 G4cout << "==> Flugg FGeometryInit::PrintMagneticField()" << G4endl;
766#endif
767
768 G4cout << "\t* Printing Magnetic Field..." << G4endl;
769
26911512 770 if(fTransportationManager->GetFieldManager()->DoesFieldExist()) {
26911512 771
772 //get magnetic field pointer
773 const G4Field * pMagField =
774 fTransportationManager->GetFieldManager()->GetDetectorField();
775
26911512 776
777 if(pMagField) {
26d97e06 778 //Check if it can be made a uniform magnetic field
26911512 779 const G4UniformMagField *pUnifMagField =
780 dynamic_cast<const G4UniformMagField*>(pMagField);
781 if(pUnifMagField) {
26d97e06 782 G4double B[3];
783 G4double point[4]; //it is not really used
784 pUnifMagField->GetFieldValue(point,B);
785
786 //write MGNFIELD card
787 PrintHeader(os,"GEANT4 MAGNETIC FIELD");
788 os << setw10 << "MGNFIELD ";
789 os << setw10 << "";
790 os << setw10 << "";
791 os << setw10 << "";
0edf14e7 792 os.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
26d97e06 793 os << setw10 << setfixed
0edf14e7 794 << std::setprecision(4) << B[0]
26d97e06 795 << setw10 << B[1]
796 << setw10 << B[2]
797 << G4endl;
798 }
799 else {
800 G4cout << "WARNING: No Uniform Magnetic Field found." << G4endl;
801 G4cout << " Manual intervention might be needed." << G4endl;
26911512 802 }
26911512 803 }
26d97e06 804 else
805 G4cout << "\t No detector field found... " << G4endl;
26911512 806 } // end if magnetic field
26d97e06 807 else
808 G4cout << "\t No field found... " << G4endl;
809
e73e0522 810#ifdef G4GEOMETRY_DEBUG
26d97e06 811 G4cout << "<== Flugg FGeometryInit::PrintMagneticField()" << G4endl;
e73e0522 812#endif
26911512 813}
6a53de92 814
815int FGeometryInit::CurrentVolID(int ir, int& copyNo)
816{
086d0acf 817 if (ir == 0)
818 {
819 copyNo = -1;
820 return -1;
821 }
822
6a53de92 823 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
086d0acf 824 G4VPhysicalVolume * physicalvol = (*pVolStore)[ir- 1];
825
826 if (physicalvol) {
827 copyNo = physicalvol->GetCopyNo();
828 } else {
829 copyNo = -1;
830 return -1;
831 }
832
833
6a53de92 834 int id = fVolIdVolumeMap[physicalvol];
835 return id;
836}
837
838int FGeometryInit::CurrentVolOffID(int ir, int off, int& copyNo)
839{
840 if (off == 0) return CurrentVolID(ir, copyNo);
841
842 G4PhysicalVolumeStore* pVolStore = G4PhysicalVolumeStore::GetInstance();
843 G4VPhysicalVolume* physicalvol = (*pVolStore)[ir- 1];
844 G4VPhysicalVolume* mother = physicalvol;
845
0edf14e7 846 int index;
847//============================================================================
848 if (mother) {
849 // Check touchable depth
850 //
851 if (ptrTouchHist->GetHistoryDepth() < off) {
852 mother = 0;
853 } else {
854 // Get the off-th mother
855 index = ptrTouchHist->GetHistoryDepth() - off;
856 // in the touchable history volumes are ordered
857 // from top volume up to mother volume;
858 // the touchable volume is not in the history
859 mother = ptrTouchHist->GetHistory()->GetVolume(index);
860 }
861 }
862//============================================================================
6a53de92 863
864 int id;
865
866 if (!mother) {
867 G4cout << "Flugg FGeometryInit::CurrentVolOffID mother not found" << G4endl;
868 id = -1;
869 copyNo = -1;
870 } else {
871 copyNo = mother ->GetCopyNo();
872 id = fVolIdVolumeMap[mother];
873 }
874 return id;
875}
dc37cac6 876
877void FGeometryInit::Gmtod(double* xm, double* xd, int iflag)
878{
879// Transforms a position from the world reference frame
880// to the current volume reference frame.
881//
882// Geant3 desription:
883// ==================
884// Computes coordinates XD (in DRS)
885// from known coordinates XM in MRS
886// The local reference system can be initialized by
887// - the tracking routines and GMTOD used in GUSTEP
888// - a call to GMEDIA(XM,NUMED)
889// - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
890// (inverse routine is GDTOM)
891//
892// If IFLAG=1 convert coordinates
893// IFLAG=2 convert direction cosinus
894//
895// ---
896 FluggNavigator * ptrNavig = getNavigatorForTracking();
897 //setting variables (and dimension: Fluka uses cm.!)
898 G4ThreeVector pGlob(xm[0],xm[1],xm[2]);
1b4955bb 899 G4ThreeVector pLoc;
dc37cac6 900
901 if (iflag == 1) {
1b4955bb 902 pGlob *= 10.0; // in mm
0edf14e7 903// change because of geant 4 6.0
904// pLoc = ptrNavig->ComputeLocalPoint(pGlob);
905 pLoc = ptrNavig->GetGlobalToLocalTransform().TransformPoint(pGlob);
906
1b4955bb 907 pLoc /= 10.0; // in cm
dc37cac6 908 } else if (iflag == 2) {
909 pLoc =
910 ptrNavig->ComputeLocalAxis(pGlob);
911 } else {
912 G4cout << "Flugg FGeometryInit::Gmtod called with undefined flag" << G4endl;
913 }
dc37cac6 914 xd[0] = pLoc[0]; xd[1] = pLoc[1]; xd[2] = pLoc[2];
915}
916
917void FGeometryInit::Gdtom(double* xd, double* xm, int iflag)
918{
919// Transforms a position from the current volume reference frame
920// to the world reference frame.
921//
922// Geant3 desription:
923// ==================
924// Computes coordinates XM (Master Reference System
925// knowing the coordinates XD (Detector Ref System)
926// The local reference system can be initialized by
927// - the tracking routines and GDTOM used in GUSTEP
928// - a call to GSCMED(NLEVEL,NAMES,NUMBER)
929// (inverse routine is GMTOD)
930//
931// If IFLAG=1 convert coordinates
932// IFLAG=2 convert direction cosinus
933//
934// ---
935
936 FluggNavigator * ptrNavig = getNavigatorForTracking();
937 G4ThreeVector pLoc(xd[0],xd[1],xd[2]);
1b4955bb 938
dc37cac6 939 G4ThreeVector pGlob;
940 if (iflag == 1) {
1b4955bb 941 pLoc *= 10.0; // in mm
dc37cac6 942 pGlob = ptrNavig->GetLocalToGlobalTransform().
943 TransformPoint(pLoc);
1b4955bb 944 pGlob /= 10.0; // in cm
dc37cac6 945 } else if (iflag == 2) {
946 pGlob = ptrNavig->GetLocalToGlobalTransform().
947 TransformAxis(pLoc);
948 } else {
949 G4cout << "Flugg FGeometryInit::Gdtom called with undefined flag" << G4endl;
950 }
951
952 xm[0] = pGlob[0]; xm[1] = pGlob[1]; xm[2] = pGlob[2];
953}