]> git.uio.no Git - u/mrichter/AliRoot.git/blame - Flugg/source/Wrappers/src/FGeometryInit.cc
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Removing unnecessary files for new the new Flugg
[u/mrichter/AliRoot.git] / Flugg / source / Wrappers / src / FGeometryInit.cc
CommitLineData
1fcb989d 1// $Id$
2// Flugg tag $Name$
3
4#include "FGeometryInit.hh"
5#include <stdio.h>
6
7FGeometryInit * FGeometryInit::flagInstance=0;
8
9FGeometryInit* FGeometryInit::GetInstance()
10{
11 if (!flagInstance) new FGeometryInit();
12
13 return flagInstance;
14}
15
16
17FGeometryInit::FGeometryInit()
18 {
19 flagInstance = this;
20 fTransportationManager = G4TransportationManager::GetTransportationManager();
21 }
22
23
24FGeometryInit::~FGeometryInit()
25 {
26 DeleteHistories();
27 ptrGeoMan->OpenGeometry();
28 delete fTransportationManager;
29 delete ptrJrLtGeant;
30 DelHistArray();
31
32 //keep ATTENTION: never delete a pointer twice!
33 }
34
35
36G4Navigator* FGeometryInit::getNavigatorForTracking()
37{
38 return fTransportationManager->GetNavigatorForTracking();
39}
40
41void FGeometryInit::setDetConstruction(G4VUserDetectorConstruction* detector)
42{
43 fDetector = detector;;
44}
45
46void FGeometryInit::setDetector()
47{
48 myTopNode = fDetector->Construct();
49}
50
51void FGeometryInit::setMotherVolume()
52{
53 fTransportationManager->GetNavigatorForTracking()->SetWorldVolume(myTopNode);
54}
55
56void FGeometryInit::closeGeometry()
57{
58 ptrGeoMan = G4GeometryManager::GetInstance();
59 ptrGeoMan->OpenGeometry();
60
61 //true argoment allows voxel construction; if false voxels are built
62 //only for replicated volumes
63 ptrGeoMan->CloseGeometry(true);
64}
65
66G4FieldManager * FGeometryInit::getFieldManager()
67{
68 return fTransportationManager->GetFieldManager();
69}
70
71//*************************************************************************
72
73void FGeometryInit::InitHistArray()
74{
75 ptrArray = new G4int[1000000];
76 for(G4int i=0;i<1000000;i++) ptrArray[i]=0;
77}
78
79void FGeometryInit::DelHistArray()
80{
81 delete ptrArray;
82}
83
84G4int * FGeometryInit::GetHistArray()
85{
86 return ptrArray;
87}
88
89
90
91//*************************************************************************
92//jrLtGeant stores all crossed lattice volume histories.
93
94void FGeometryInit::InitJrLtGeantArray()
95{
96#ifdef G4GEOMETRY_DEBUG
97 G4cout<<"Initializing JrLtGeant array"<<G4endl;
98#endif
99 ptrJrLtGeant = new G4int[10000];
100 for(G4int x=0;x<10000;x++) ptrJrLtGeant[x]=-1;
101 flagLttcGeant = -1;
102}
103
104
105G4int * FGeometryInit::GetJrLtGeantArray()
106{
107 return ptrJrLtGeant;
108}
109
110
111G4int FGeometryInit::GetLttcFlagGeant()
112{
113 return flagLttcGeant;
114}
115
116void FGeometryInit::SetLttcFlagGeant(G4int newFlagLttc)
117{
118 // Added by A.Solodkov
119 if (newFlagLttc >= 10000) {
120 G4cout<<"Problems in FGeometryInit::SetLttcFlagGeant"<<G4endl;
121 G4cout<<"Index newFlagLttc="<<newFlagLttc<<" is outside array bounds"<<G4endl;
122 G4cout<<"Better to stop immediately !"<<G4endl;
123 exit(1);
124 }
125 flagLttcGeant = newFlagLttc;
126}
127
128void FGeometryInit::PrintJrLtGeant()
129{
130#ifdef G4GEOMETRY_DEBUG
131 //G4cout<<"jrLtGeant:"<<G4endl;
132 //for(G4int y=0;y<=flagLttcGeant;y++)
133 //
134 // G4cout<<"jrLtGeant["<<y<<"]="<<ptrJrLtGeant[y]<<G4endl;
135#endif
136}
137
138//**************************************************************************
139
140void FGeometryInit::PrintHistories()
141{
142/*
143#ifdef G4GEOMETRY_DEBUG
144 G4cout<<"Touch hist:"<<G4endl;
145 G4cout<<*(ptrTouchHist->GetHistory())<<G4endl;
146 G4cout<<"Tmp hist:"<<G4endl;
147 G4cout<<*(ptrTempNavHist->GetHistory())<<G4endl;
148 G4cout<<"Old hist:"<<G4endl;
149 G4cout<<*(ptrOldNavHist->GetHistory())<<G4endl;
150#endif
151*/
152}
153
154
155
156
157void FGeometryInit::InitHistories()
158{
159#ifdef G4GEOMETRY_DEBUG
160 G4cout <<" InitHistories start" << G4endl;
161#endif
162 //init utility histories with navigator history
163
164 G4cout << fTransportationManager<< G4endl;
165 G4cout << fTransportationManager->GetNavigatorForTracking()<< G4endl;
166
167 ptrTouchHist = fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
168 G4cout << "Touchable history " << G4endl;
169 ptrTempNavHist = fTransportationManager->GetNavigatorForTracking()->CreateTouchableHistory();
170 ptrOldNavHist = new G4TouchableHistory();
171#ifdef G4GEOMETRY_DEBUG
172 G4cout <<" InitHistories end" << G4endl;
173#endif
174}
175
176void FGeometryInit::DeleteHistories()
177 {
178 delete ptrTouchHist;
179 delete ptrOldNavHist;
180 delete ptrTempNavHist;
181
182#ifdef G4GEOMETRY_DEBUG
183 G4cout<<"Deleting step-history objects at end of run!"<<G4endl;
184#endif
185 }
186
187G4TouchableHistory * FGeometryInit::GetTouchableHistory()
188{
189 return ptrTouchHist;
190}
191
192G4TouchableHistory * FGeometryInit::GetOldNavHist()
193{
194 return ptrOldNavHist;
195}
196
197G4TouchableHistory * FGeometryInit::GetTempNavHist()
198{
199 return ptrTempNavHist;
200}
201
202
203void FGeometryInit::UpdateHistories(const G4NavigationHistory * history,
204 const G4int & flagHist)
205{
206 PrintHistories();
207
208#ifdef G4GEOMETRY_DEBUG
209 G4cout<<"...updating histories!"<<G4endl;
210#endif
211
212 G4VPhysicalVolume * pPhysVol = history->GetTopVolume();
213
214 switch (flagHist)
215 {
216 case 0:
217 {
218 //this is the case when a new history is given to the
219 //navigator and old history has to be resetted
220 //touchable history has not been updated jet, so:
221
222 ptrTouchHist->UpdateYourself(pPhysVol,history);
223 ptrTempNavHist->UpdateYourself(pPhysVol,history);
224 G4NavigationHistory * ptrOldNavHistNotConst =
225 const_cast<G4NavigationHistory * >(ptrOldNavHist->GetHistory());
226 ptrOldNavHistNotConst->Reset();
227 ptrOldNavHistNotConst->Clear();
228 PrintHistories();
229 break;
230 }
231
232 case 1:
233 {
234 //this is the case when a new history is given to the
235 //navigator but old history has to be kept (e.g. LOOKZ
236 //is call during an event);
237 //touchable history has not been updated jet, so:
238
239 ptrTouchHist->UpdateYourself(pPhysVol,history);
240 ptrTempNavHist->UpdateYourself(pPhysVol,history);
241 PrintHistories();
242 break;
243 }
244
245 case 2:
246 {
247 //this is the case when the touchable history has been
248 //updated by a LocateGlobalPointAndUpdateTouchable call
249
250 G4VPhysicalVolume * pPhysVolTemp = ptrTempNavHist->GetVolume();
251 ptrOldNavHist->UpdateYourself(pPhysVolTemp,
252 ptrTempNavHist->GetHistory());
253
254 ptrTempNavHist->UpdateYourself(pPhysVol,history);
255 PrintHistories();
256 break;
257 }
258 default:
259 {
260 G4cout<<"ERROR in updating step-histories!"<<G4endl;
261 break;
262 }
263 }
264
265}
266
267//*****************************************************************************
268
269
270void FGeometryInit::createFlukaMatFile()
271{
272// ultima modifica Sara Vanini 1/III/99
273// NOMI DI ELEMENTI E COMPOSTI: i nomi devono essere scritti maiuscolo,
274// secondo lo standard di fluka. Devono inoltre essere uguali ai nomi dei
275// materiali fluka - vedere il manuale di fluka - perche` il programma
276// carichi le giuste sezioni d`urto, e uguali ai nomi inclusi nel .pemf.
277// Altrimenti l`utente deve definirsi le CARDS LOW-MAT, e costruirsi il
278// .pemf, per avere le giuste sezioni d`urto caricate in memoria.
279
280 //flag
281#ifdef G4GEOMETRY_DEBUG
282 G4cout<<"================== FILEWR ================="<<G4endl;
283#endif
284
285 //open file for output
286 ofstream fout("flukaMat.inp");
287
288 //PhysicalVolumeStore, Volume and MaterialTable pointers
289 G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
290 G4int numVol = G4int(pVolStore->size());
291
292 G4int* indexMatFluka = 0;
293 static G4Material * ptrMat = 0;
294 G4int x = 0;
295 while(!ptrMat && x<numVol)
296 {
297 G4VPhysicalVolume * ptrVol = (*pVolStore)[x];
298 G4LogicalVolume * ptrLogVol = ptrVol->GetLogicalVolume();
299 ptrMat = ptrLogVol->GetMaterial();
300 x+=1;
301 }
302
303 if(ptrMat)
304 {
305 static const G4MaterialTable * ptrMatTab = G4Material::GetMaterialTable();
306
307 //number of materials, elements, variable initialisations
308 static size_t totNumMat = G4Material::GetNumberOfMaterials();
309#ifdef G4GEOMETRY_DEBUG
310 G4cout<<"Number of materials: "<<totNumMat<<G4endl;
311#endif
312 const G4Element * ptrElem = ptrMat->GetElement(0);
313 static size_t totNumElem = ptrElem->GetNumberOfElements();
314 static const G4ElementTable * ptrElemTab = ptrElem->GetElementTable();
315 G4int * elemIndexInMATcard = new G4int[totNumElem];
316 for(G4int t=0; t<totNumElem; t++)
317 elemIndexInMATcard[t] = 0;
318 static const G4IsotopeTable * ptrIsotTab;
319 G4int initIsot = 0;
320 static size_t totNumIsot;
321 G4int* isotIndexInMATcard = 0;
322 G4int isotPresence = 0;
323
324
325 // title
326 fout<<"*\n"<<"*\n"<<"*\n";
327 fout<<"********************* GEANT4 ELEMENTS AND COMPOUNDS *********************\n"<<"*\n";
328 fout<<("*...+....1....+....2....+....3....+....4....+....5....+....6....+....7...")<<G4endl;
329 fout<<("*")<<G4endl;
330
331 // *** loop over G4Materials to assign Fluka index
332 G4int indexCount=3;
333 indexMatFluka = new G4int[totNumMat];
334 for(G4int i=0; i<totNumMat; i++)
335 {
336 //pointer material, state
337 ptrMat = (*ptrMatTab)[i];
338 G4double denMat = ptrMat->GetDensity();
339 G4String nameMat = ptrMat->GetName();
340 // Fluka index: bh=1, vacuum=2, others=3,4..
341 if(denMat<=1.00e-10*g/cm3)
342 //N.B. fluka density limit decided on XI-`98
343 {
344 indexMatFluka[i]=2;
345 }
346 else
347 {
348 indexMatFluka[i]=indexCount;
349 indexCount+=1;
350 }
351
352 // *** write single-element material MATERIAL card
353 size_t numElem = ptrMat->GetNumberOfElements();
354 if(numElem==1)
355 {
356 G4int index = indexMatFluka[i];
357 const G4Element * ptrElem = ptrMat->GetElement(0);
358 size_t indElemTab = ptrElem->GetIndex();
359 size_t numIsot = ptrElem->GetNumberOfIsotopes();
360 G4double A = (ptrElem->GetA())/(g);
361 if(!numIsot)
362 {
363 if(index!=2 && !elemIndexInMATcard[indElemTab])
364 {
365 G4double Z = ptrElem->GetZ();
366 elemIndexInMATcard[indElemTab] = index;
367 G4String nameEl = ptrElem->GetName();
368 nameEl.toUpper();
369
370 //write on file MATERIAL card of element
371 fout<<G4std::setw(10)<<"MATERIAL ";
372 fout.setf(0,G4std::ios::floatfield);
373 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::fixed)
374 <<G4std::setprecision(1)<<Z;
375 fout<<G4std::setw(10)<<G4std::setprecision(3)<<A;
376 fout.setf(0,G4std::ios::floatfield);
377 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::scientific)
378 <<G4std::setprecision(3)<<denMat/(g/cm3);
379 fout.setf(0,G4std::ios::floatfield);
380 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::fixed)<<
381 G4std::setprecision(1)<<G4double(index);
382 fout<<G4std::setw(10)<<" ";
383 fout<<G4std::setw(10)<<" ";
384 fout<<nameEl<<G4endl;
385 }
386 }
387 if(numIsot==1)
388 {
389 // G4Isotope pointer
390 const G4Isotope* ptrIsot = ptrElem->GetIsotope(0);
391 size_t indIsotTab = ptrIsot->GetIndex();
392 //initialize variables
393 if(!initIsot)
394 {
395 totNumIsot = ptrIsot->GetNumberOfIsotopes();
396 ptrIsotTab = ptrIsot->GetIsotopeTable();
397 isotIndexInMATcard = new G4int[totNumIsot];
398 for(G4int t=0; t<totNumIsot; t++) isotIndexInMATcard[t] = 0;
399 initIsot = 1;
400 }
401 if(!isotIndexInMATcard[indIsotTab])
402 {//compute physical data and counters
403 G4int ZIs = ptrIsot->GetZ();
404 G4double AIs = (ptrIsot->GetA())/(g);
405 G4int NIs = ptrIsot->GetN();
406 G4String nameIsot = ptrIsot->GetName();
407 nameIsot.toUpper();
408 G4double* ptrRelAbVect = ptrElem->GetRelativeAbundanceVector();
409 isotIndexInMATcard[indIsotTab] = index;
410
411 //write on file MATERIAL card of isotope
412 fout<<G4std::setw(10)<<"MATERIAL ";
413 fout.setf(0,G4std::ios::floatfield);
414 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::fixed)
415 <<G4std::setprecision(1)<<G4double(ZIs);
416 fout<<G4std::setw(10)<<G4std::setprecision(3)<<AIs;
417 fout.setf(0,ios::floatfield);
418 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::scientific)
419 <<G4std::setprecision(3)<<denMat/(g/cm3);
420 fout.setf(0,G4std::ios::floatfield);
421 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::fixed)
422 <<G4std::setprecision(1)<<G4double(index);
423 fout<<G4std::setw(10)<<" ";
424 fout<<G4std::setw(10)<<G4double(NIs);
425 fout<<nameIsot<<G4endl;
426 }
427 }// end if(numIsot==1)
428 }// end if(numElem==1)
429 }// end for loop
430
431
432 // *** material definitions: elements, compound made of G4Elements
433 // or made of G4Materials
434 for(G4int j=0; j<totNumMat; j++)
435 {
436 //pointer material, and material data
437 ptrMat = (*ptrMatTab)[j];
438 size_t numElem = ptrMat->GetNumberOfElements();
439 G4double densityMat = (ptrMat->GetDensity())/(g/cm3);
440 G4String nameMat = ptrMat->GetName();
441 nameMat.toUpper();
442 isotPresence = 0;
443
444 //fraction vector of compounds of the material
445 const G4double* ptrFracVect = ptrMat->GetFractionVector();
446
447 //loop on elements of the material
448 for(G4int el=0; el<numElem; el++)
449 {
450 //compute physical data, initialize variables
451 const G4Element * ptrElem = ptrMat->GetElement(el);
452 size_t indElemTab = ptrElem->GetIndex();
453 G4String nameElem = ptrElem->GetName();
454 nameElem.toUpper();
455 size_t numIsot = ptrElem->GetNumberOfIsotopes();
456 G4double A = (ptrElem->GetA())/(g);
457
458 if(!numIsot)
459 {
460 if(!elemIndexInMATcard[indElemTab])
461 {
462 G4double Z = ptrElem->GetZ();
463 G4double density = ptrFracVect[el]*densityMat;
464 elemIndexInMATcard[indElemTab] = indexCount;
465
466 //write on file MATERIAL card of element
467 fout<<G4std::setw(10)<<"MATERIAL ";
468 fout.setf(0,G4std::ios::floatfield);
469 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::fixed)
470 <<G4std::setprecision(1)<<Z;
471 fout<<G4std::setw(10)<<G4std::setprecision(3)<<A;
472 fout.setf(0,G4std::ios::floatfield);
473 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::scientific)
474 <<G4std::setprecision(3)<<density;
475 fout.setf(0,G4std::ios::floatfield);
476 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::fixed)<<
477 G4std::setprecision(1)<<G4double(indexCount);
478 fout<<G4std::setw(10)<<" ";
479 fout<<G4std::setw(10)<<" ";
480 fout<<nameElem<<G4endl;
481 indexCount+=1;
482 }
483 }
484
485 else
486 {
487 if(numIsot>=2) isotPresence = 1;
488 //loop on isotopes
489 for(G4int nis=0; nis<numIsot; nis++)
490 {
491 // G4Isotope pointer
492 const G4Isotope* ptrIsot = ptrElem->GetIsotope(nis);
493 size_t indIsotTab = ptrIsot->GetIndex();
494 //initialize variables
495 if(!initIsot)
496 {
497 totNumIsot = ptrIsot->GetNumberOfIsotopes();
498 ptrIsotTab = ptrIsot->GetIsotopeTable();
499 isotIndexInMATcard = new G4int[totNumIsot];
500 for(G4int t=0; t<totNumIsot; t++) isotIndexInMATcard[t] = 0;
501 initIsot = 1;
502 }
503 if(!isotIndexInMATcard[indIsotTab])
504 {//compute physical data and counters
505 G4int ZIs = ptrIsot->GetZ();
506 G4double AIs = (ptrIsot->GetA())/(g);
507 G4int NIs = ptrIsot->GetN();
508 G4String nameIsot = ptrIsot->GetName();
509 nameIsot.toUpper();
510 G4double* ptrRelAbVect = ptrElem->GetRelativeAbundanceVector();
511 G4double density = ptrFracVect[el]*densityMat*
512 ptrRelAbVect[nis]*AIs/A;
513 G4int index = indexCount;
514 isotIndexInMATcard[indIsotTab] = indexCount;
515 indexCount+=1;
516
517 //write on file MATERIAL card of isotope
518 fout<<G4std::setw(10)<<"MATERIAL ";
519 fout.setf(0,G4std::ios::floatfield);
520 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::fixed)
521 <<G4std::setprecision(1)<<G4double(ZIs);
522 fout<<G4std::setw(10)<<G4std::setprecision(3)<<AIs;
523 fout.setf(0,G4std::ios::floatfield);
524 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::scientific)
525 <<G4std::setprecision(3)<<density;
526 fout.setf(0,G4std::ios::floatfield);
527 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::fixed)
528 <<G4std::setprecision(1)<<G4double(index);
529 fout<<G4std::setw(10)<<" ";
530 fout<<G4std::setw(10)<<G4double(NIs);
531 fout<<nameIsot<<G4endl;
532 }
533 }
534 }
535
536 }
537
538 if(numElem>1 || isotPresence==1)
539 {
540 // write MATERIAL+COMPOUND card specifing the compound
541
542 //flags for writing COMPOUND card
543 G4int treCount=0;
544
545 //make MATERIAL card for compound, start COMPOUND card
546 fout<<"*"<<G4endl;
547 fout<<"* Define GEANT4 compound "<<nameMat<<G4endl;
548 fout<<G4std::setw(10)<<"MATERIAL ";
549 fout.setf(0,G4std::ios::floatfield);
550 fout<<G4std::setw(10)<<" "<<G4std::setw(10)<<" ";
551 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::scientific)
552 <<G4std::setprecision(3)<<densityMat;
553 fout.setf(0,G4std::ios::floatfield);
554 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::fixed)<<G4std::setprecision(1)
555 <<G4double(indexMatFluka[j]);
556 fout<<G4std::setw(10)<<" ";
557 fout<<G4std::setw(10)<<" ";
558 fout<<nameMat<<G4endl;
559 fout<<G4std::setw(10)<<"COMPOUND ";
560
561
562 //write elements in COMPOUND card
563 for(G4int h=0; h<numElem; h++)
564 {
565 const G4Element * ptrElemMat = ptrMat->GetElement(h);
566 size_t indexElemMat = ptrElemMat->GetIndex();
567 size_t numIsotElem = ptrElemMat->GetNumberOfIsotopes();
568 if(!numIsotElem)
569 {
570 if(treCount==3)
571 {
572 treCount=0;
573 fout<<nameMat<<G4endl;
574 fout<<G4std::setw(10)<<"COMPOUND ";
575 }
576
577 fout.setf(0,G4std::ios::floatfield);
578 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::scientific)
579 <<G4std::setprecision(2)<<-ptrFracVect[h];
580 fout.setf(0,G4std::ios::floatfield);
581 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::fixed)<<
582 G4std::setprecision(1)<<G4double(elemIndexInMATcard[indexElemMat]);
583 treCount+=1;
584 }
585 else
586 {
587 G4double * ptrIsotAbbVect = ptrElemMat->GetRelativeAbundanceVector();
588
589 for(G4int iso=0; iso<numIsotElem; iso++)
590 {
591 const G4Isotope * ptrIsotElem =ptrElemMat->GetIsotope(iso);
592 size_t indexIsotMat = ptrIsotElem->GetIndex();
593 G4double isotAbundPerVol =
594 ptrIsotAbbVect[iso]*Avogadro*densityMat*
595 ptrFracVect[h]/(ptrElemMat->GetA()/(g));
596
597 if(treCount==3)
598 {
599 treCount=0;
600 fout<<nameMat<<G4endl;
601 fout<<G4std::setw(10)<<"COMPOUND ";
602 }
603 fout.setf(0,G4std::ios::floatfield);
604 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::scientific)
605 <<G4std::setprecision(2)<<isotAbundPerVol;
606 fout.setf(0,G4std::ios::floatfield);
607 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::fixed)<<
608 G4std::setprecision(1)<<G4double(isotIndexInMATcard[indexIsotMat]);
609 treCount+=1;
610 }
611 }
612 }
613
614 //end COMPOUND card
615 if(treCount==1) fout<<G4std::setw(10)<<" "<<G4std::setw(10)<<" "<<
616 G4std::setw(10)<<" "<<G4std::setw(10)<<" "<<nameMat<<G4endl;
617 if(treCount==2) fout<<G4std::setw(10)<<" "<<G4std::setw(10)<<
618 " "<<nameMat<<G4endl;
619 if(treCount==3) fout<<nameMat<<G4endl;
620 fout<<"*"<<G4endl;
621 }
622
623
624 } // end for loop
625 delete elemIndexInMATcard;
626 if(initIsot) delete isotIndexInMATcard;
627
628 } // end if (ptrMat)
629
630 // *** material-volume correspondence
631 fout<<"*\n"<<"*\n"<<"*\n";
632 fout<<"******************** GEANT4 MATERIAL ASSIGNMENTS *******************\n"<<"*\n";
633 fout<<("*...+....1....+....2....+....3....+....4....+....5....+....6....+....7...")<<G4endl;
634 fout<<("*")<<G4endl;
635
636 //initializations
637 G4int indexMatOld = 0;
638 G4int indexRegFlukaFrom = 0;
639 G4int indexRegFlukaTo = 0;
640 G4int existTo = 0;
641 G4int flagField = 0;
642 G4int lastFlagField = 0;
643
644 //open file for volume-index correspondence
645 ofstream vout("Volumes_index.inp");
646
647 //... and write title
648 vout<<"*"<<G4endl;
649 vout<<"******************** GEANT4 VOLUMES *******************\n";
650 vout<<("*")<<G4endl;
651
652 //loop su tutti i volumi
653 for(G4int l=0;l<numVol;l++)
654 {
655 //index of the region
656 G4VPhysicalVolume * ptrVol = (*pVolStore)[l];
657 G4LogicalVolume * ptrLogVol = ptrVol->GetLogicalVolume();
658 G4int indexRegFluka = l+1;
659
660
661 //write index volume and name on file Volumes_index.inp
662 vout.setf(G4std::ios::left,G4std::ios::adjustfield);
663 vout<<G4std::setw(10)<<indexRegFluka;
664 vout<<G4std::setw(20)<<ptrVol->GetName()<<G4std::setw(20)<<"";
665 if(ptrVol->IsReplicated())
666 {
667 EAxis axis;
668 G4int nRep;
669 G4double width;
670 G4double offset;
671 G4bool consum;
672 ptrVol->GetReplicationData(axis,nRep,width,offset,consum);
673 vout.setf(G4std::ios::left,G4std::ios::adjustfield);
674 vout<<G4std::setw(10)<<"Repetion Nb: "<<G4std::setw(3)<<nRep;
675 }
676 vout<<G4endl;
677
678 //check if Magnetic Field is present in the region
679 G4FieldManager * pMagFieldMan = ptrLogVol->GetFieldManager();
680 const G4Field * pMagField = 0;
681 if(pMagFieldMan) pMagField = pMagFieldMan->GetDetectorField();
682 if(pMagField) flagField = 1;
683 else flagField = 0;
684
685
686 //index of material in the region
687 G4Material * ptrMat = ptrLogVol->GetMaterial();
688 G4int indexMat;
689 if(ptrMat)
690 {
691 size_t indexMatGeant = ptrMat->GetIndex();
692 indexMat = indexMatFluka[indexMatGeant];
693 }
694 else indexMat = 2;
695
696 //if materials are repeated
697 if(indexMat==indexMatOld && flagField==lastFlagField)
698 {
699 indexRegFlukaTo=indexRegFluka;
700 existTo=1;
701 if(l==(numVol-1))
702 {
703 fout<<G4std::setw(10)<<G4double(indexRegFlukaTo);
704 fout<<G4std::setw(10)<<"0.0";
705 fout<<G4std::setw(10)<<G4double(flagField);
706 fout<<G4std::setw(10)<<"0.0"<<G4endl;
707 }
708
709 }
710 else
711 {
712 //write on file ASSIGNMAT card
713
714 //first complete last material card
715 if(!existTo)
716 {
717 if(l)
718 {
719 fout<<G4std::setw(10)<<"0.0";
720 fout<<G4std::setw(10)<<"0.0";
721 fout<<G4std::setw(10)<<G4double(lastFlagField);
722 fout<<G4std::setw(10)<<"0.0"<<G4endl;
723 }
724 }
725 else
726 {
727 fout<<G4std::setw(10)<<G4double(indexRegFlukaTo);
728 fout<<G4std::setw(10)<<"0.0";
729 fout<<G4std::setw(10)<<G4double(lastFlagField);
730 fout<<G4std::setw(10)<<"0.0"<<G4endl;
731 }
732
733 // begin material card
734 fout<<G4std::setw(10)<<"ASSIGNMAT ";
735 fout.setf(0,G4std::ios::floatfield);
736 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::fixed)<<
737 G4std::setprecision(1)<<G4double(indexMat);
738 fout<<G4std::setw(10)<<G4double(indexRegFluka);
739
740 existTo=0;
741 indexRegFlukaFrom=indexRegFluka;
742
743 if(l==(numVol-1))
744 {
745 fout<<G4std::setw(10)<<"0.0";
746 fout<<G4std::setw(10)<<"0.0";
747 fout<<G4std::setw(10)<<G4double(flagField);
748 fout<<G4std::setw(10)<<"0.0"<<G4endl;
749 }
750 }
751 lastFlagField = flagField;
752 indexMatOld = indexMat;
753 } // end of loop ??
754
755 //assign material 1 to black-hole=n.vol+1
756 fout<<G4std::setw(10)<<"ASSIGNMAT ";
757 fout<<G4std::setw(10)<<"1.0";
758 fout<<G4std::setw(10)<<G4double(numVol+1);
759 fout<<G4std::setw(10)<<"0.0";
760 fout<<G4std::setw(10)<<"0.0";
761 fout<<G4std::setw(10)<<"0.0";
762 fout<<G4std::setw(10)<<"0.0"<<G4endl;
763
764 // *** magnetic field ***
765 if(fTransportationManager->GetFieldManager()->DoesFieldExist())
766 {
767 fout<<"*\n"<<"*\n"<<"*\n";
768 fout<<"*********************** GEANT4 MAGNETIC FIELD ************************\n"<<"*\n";
769 fout<<("*...+....1....+....2....+....3....+....4....+....5....+....6....+....7...")<<G4endl;
770 fout<<("*")<<G4endl;
771
772 //get magnetic field pointer
773 const G4Field * pMagField = fTransportationManager->GetFieldManager()->GetDetectorField();
774
775 //if uniform magnetic field, get value
776 G4double Bx=0.0;
777 G4double By=0.0;
778 G4double Bz=0.0;
779
780 if(pMagField)
781 {
782 G4ThreeVector*Field = new G4ThreeVector(1.,2.,3.);
783 const G4UniformMagField *pUnifMagField =
784 dynamic_cast<const G4UniformMagField*>(pMagField);
785 if(pUnifMagField)
786 {
787 G4double *pFieldValue = 0;
788 G4double *point = new G4double[3];
789 point[0] = 0.;
790 point[1] = 0.;
791 point[2] = 0.;
792 pUnifMagField->GetFieldValue(point,pFieldValue);
793 //non capisco perche' l'instruzione seguente non fa linkare. Indaga!!
794 //per ora posso usare GetFieldValue: va bene lo stesso.
795 //G4ThreeVector FieldValue = pUnifMagField->GetConstantFieldValue();
796 Bx = pFieldValue[0];
797 By = pFieldValue[1];
798 Bz = pFieldValue[2];
799 }
800
801 }
802
803 //write MGNFIELD card
804 fout<<G4std::setw(10)<<"MGNFIELD ";
805 fout<<G4std::setw(10)<<"";
806 fout<<G4std::setw(10)<<"";
807 fout<<G4std::setw(10)<<"";
808 fout.setf(0,G4std::ios::floatfield);
809 fout<<G4std::setw(10)<<G4std::setiosflags(G4std::ios::fixed)<<G4std::setprecision(4)<<Bx;
810 fout<<G4std::setw(10)<<By;
811 fout<<G4std::setw(10)<<Bz<<G4endl;
812 } // end if magnetic field
813
814 vout.close();
815 fout.close();
816 delete [] indexMatFluka;
817}