git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Modifications for Mac
[u/mrichter/AliRoot.git] / TFluka / TFlukaMCGeometry.cxx
CommitLineData
8495a208 1// @(#):$Name$:$Id$
2// Author: Andrei Gheata 10/07/2003
3
4#include "TObjString.h"
5#include "TFluka.h"
6//#include "TVirtualMCApplication.h"
7#include "TFlukaMCGeometry.h"
8#include "TGeoManager.h"
9#include "TGeoVolume.h"
10
11#include "TCallf77.h"
12
13#ifndef WIN32
14# define idnrwr idnrwr_
15# define g1wr g1wr_
16# define g1rtwr g1rtwr_
17# define conhwr conhwr_
18# define inihwr inihwr_
19# define jomiwr jomiwr_
20# define lkdbwr lkdbwr_
21# define lkfxwr lkfxwr_
22# define lkmgwr lkmgwr_
23# define lkwr lkwr_
24# define magfld magfld_
25# define nrmlwr nrmlwr_
26# define rgrpwr rgrpwr_
27# define isvhwr isvhwr_
28
29#else
30
31# define idnrwr IDNRWR
32# define g1wr G1WR
33# define g1rtwr G1RTWR
34# define conhwr CONHWR
35# define inihwr INIHWR
36# define jomiwr JOMIWR
37# define lkdbwr LKDBWR
38# define lkfxwr LKFXWR
39# define lkmgwr LKMGWR
40# define lkwr LKWR
41# define magfld MAGFLD
42# define nrmlwr NRMLWR
43# define rgrpwr RGRPWR
44# define isvhwr ISVHWR
45
46#endif
47
48//____________________________________________________________________________
49extern "C"
50{
51 //
52 // Prototypes for FLUKA navigation methods
53 //
54 Int_t type_of_call idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/);
55 void type_of_call g1wr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
56 Double_t * /*pV*/, Int_t & /*oldReg*/ , const Int_t & /*oldLttc*/, Double_t & /*propStep*/,
57 Int_t & /*nascFlag*/, Double_t & /*retStep*/, Int_t & /*newReg*/,
58 Double_t & /*saf*/, Int_t & /*newLttc*/, Int_t & /*LttcFlag*/,
59 Double_t *s /*Lt*/, Int_t * /*jrLt*/);
60
61 void type_of_call g1rtwr();
62 void type_of_call conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/);
63 void type_of_call inihwr(Int_t & /*intHist*/);
64 void type_of_call jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
65 Int_t & /*flukaReg*/);
66 void type_of_call lkdbwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
67 Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
68 Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/);
69 void type_of_call lkfxwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
70 Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
71 Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/);
72 void type_of_call lkmgwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
73 Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
74 Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
75 void type_of_call lkwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
76 Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
77 Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/);
efde9b4d 78// void type_of_call magfld(const Double_t & /*pX*/, const Double_t & /*pY*/, const Double_t & /*pZ*/,
79// Double_t & /*cosBx*/, Double_t & /*cosBy*/, Double_t & /*cosBz*/,
80// Double_t & /*Bmag*/, Int_t & /*reg*/, Int_t & /*idiscflag*/);
8495a208 81 void type_of_call nrmlwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
82 Double_t & /*pVx*/, Double_t & /*pVy*/, Double_t & /*pVz*/,
83 Double_t * /*norml*/, const Int_t & /*oldReg*/,
84 const Int_t & /*newReg*/, Int_t & /*flagErr*/);
85 void type_of_call rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
86 Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/);
87 Int_t type_of_call isvhwr(const Int_t & /*fCheck*/, const Int_t & /*intHist*/);
88};
89
90// TFluka global pointer
91TFluka *fluka = 0;
92TFlukaMCGeometry *mcgeom = 0;
93
94ClassImp(TFlukaMCGeometry)
95
96TFlukaMCGeometry* TFlukaMCGeometry::fgInstance=0;
97
98//_____________________________________________________________________________
99TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title)
100 : TVirtualMCGeometry(name, title)
101{
102 //
103 // Standard constructor
104 //
efde9b4d 105 fLastMaterial = 0;
b1536e91 106 fNextRegion = 0;
107 fNextLattice = 0;
108 fRegionList = 0;
8495a208 109 fluka = (TFluka*)gMC;
110 mcgeom = this;
111}
112
113//_____________________________________________________________________________
114TFlukaMCGeometry::TFlukaMCGeometry()
115 : TVirtualMCGeometry()
116{
117 //
118 // Default constructor
119 //
efde9b4d 120 fLastMaterial = 0;
b1536e91 121 fNextRegion = 0;
122 fNextLattice = 0;
123 fRegionList = 0;
8495a208 124 fluka = (TFluka*)gMC;
125 mcgeom = this;
126}
127
128//_____________________________________________________________________________
129TFlukaMCGeometry::~TFlukaMCGeometry()
130{
131 //
132 // Destructor
133 //
134 fgInstance=0;
b1536e91 135 if (fRegionList) delete [] fRegionList;
8495a208 136 if (gGeoManager) delete gGeoManager;
137}
138
139//
140// private methods
141//
142//_____________________________________________________________________________
143TFlukaMCGeometry::TFlukaMCGeometry(const TFlukaMCGeometry &)
144 : TVirtualMCGeometry()
145{
146 //
147 // Copy constructor
148 //
149}
150
151//_____________________________________________________________________________
152Double_t* TFlukaMCGeometry::CreateDoubleArray(Float_t* array, Int_t size) const
153{
154// Converts Float_t* array to Double_t*,
155// !! The new array has to be deleted by user.
156// ---
157
158 Double_t* doubleArray;
159 if (size>0) {
160 doubleArray = new Double_t[size];
161 for (Int_t i=0; i<size; i++) doubleArray[i] = array[i];
162 }
163 else {
164 //doubleArray = 0;
165 doubleArray = new Double_t[1];
166 }
167 return doubleArray;
168}
169//
170// public methods
171//_____________________________________________________________________________
172void TFlukaMCGeometry::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
173 Float_t &dens, Float_t &radl, Float_t &absl,
174 Float_t* /*ubuf*/, Int_t& /*nbuf*/)
175{
176 printf("Gfmate %i\n", imat);
177 TGeoMaterial *mat;
178 TIter next (gGeoManager->GetListOfMaterials());
179 while ((mat = (TGeoMaterial*)next())) {
180 if (mat->GetUniqueID() == (UInt_t)imat) break;
181 }
182 if (!mat) {
183 Error("Gfmate", "no material with index %i found", imat);
184 return;
185 }
186 sprintf(name, "%s", mat->GetName());
187 a = mat->GetA();
188 z = mat->GetZ();
189 dens = mat->GetDensity();
190 radl = mat->GetRadLen();
191 absl = mat->GetIntLen();
192 printf(" ->material found : %s a=%g, z=%g, dens=%g, radl=%g, absl=%g\n", name, a,z,dens,radl,absl);
193}
194
195//_____________________________________________________________________________
196void TFlukaMCGeometry::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
197 Double_t &dens, Double_t &radl, Double_t &absl,
198 Double_t* /*ubuf*/, Int_t& /*nbuf*/)
199{
200 printf("Gfmate %i\n", imat);
201 TGeoMaterial *mat;
202 TIter next (gGeoManager->GetListOfMaterials());
203 while ((mat = (TGeoMaterial*)next())) {
204 if (mat->GetUniqueID() == (UInt_t)imat) break;
205 }
206 if (!mat) {
207 Error("Gfmate", "no material with index %i found", imat);
208 return;
209 }
210 sprintf(name, "%s", mat->GetName());
211 a = mat->GetA();
212 z = mat->GetZ();
213 dens = mat->GetDensity();
214 radl = mat->GetRadLen();
215 absl = mat->GetIntLen();
216 printf(" ->material found : %s a=%g, z=%g, dens=%g, radl=%g, absl=%g\n", name, a,z,dens,radl,absl);
217}
218
219//_____________________________________________________________________________
220void TFlukaMCGeometry::Material(Int_t& kmat, const char* name, Double_t a, Double_t z,
221 Double_t dens, Double_t radl, Double_t absl, Float_t* buf,
222 Int_t nwbuf)
223{
224 //
225 // Defines a Material
226 //
227 // kmat number assigned to the material
228 // name material name
229 // a atomic mass in au
230 // z atomic number
231 // dens density in g/cm3
232 // absl absorbtion length in cm
233 // if >=0 it is ignored and the program
234 // calculates it, if <0. -absl is taken
235 // radl radiation length in cm
236 // if >=0 it is ignored and the program
237 // calculates it, if <0. -radl is taken
238 // buf pointer to an array of user words
239 // nbuf number of user words
240 //
241
242 Double_t* dbuf = CreateDoubleArray(buf, nwbuf);
243 Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf);
244 delete [] dbuf;
245}
246
247//_____________________________________________________________________________
248void TFlukaMCGeometry::Material(Int_t& kmat, const char* name, Double_t a, Double_t z,
249 Double_t dens, Double_t radl, Double_t absl, Double_t* /*buf*/,
250 Int_t /*nwbuf*/)
251{
252 //
253 // Defines a Material
254 //
255 // kmat number assigned to the material
256 // name material name
257 // a atomic mass in au
258 // z atomic number
259 // dens density in g/cm3
260 // absl absorbtion length in cm
261 // if >=0 it is ignored and the program
262 // calculates it, if <0. -absl is taken
263 // radl radiation length in cm
264 // if >=0 it is ignored and the program
265 // calculates it, if <0. -radl is taken
266 // buf pointer to an array of user words
267 // nbuf number of user words
268 //
269
270 kmat = gGeoManager->GetListOfMaterials()->GetSize();
271 gGeoManager->Material(name, a, z, dens, kmat, radl, absl);
272 printf("Material %s: kmat=%i, a=%g, z=%g, dens=%g\n", name, kmat, a, z, dens);
273}
274
275//_____________________________________________________________________________
276void TFlukaMCGeometry::Mixture(Int_t& kmat, const char* name, Float_t* a, Float_t* z,
277 Double_t dens, Int_t nlmat, Float_t* wmat)
278{
279 //
280 // Defines mixture OR COMPOUND IMAT as composed by
281 // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
282 //
283 // If NLMAT > 0 then wmat contains the proportion by
284 // weights of each basic material in the mixture.
285 //
286 // If nlmat < 0 then WMAT contains the number of atoms
287 // of a given kind into the molecule of the COMPOUND
288 // In this case, WMAT in output is changed to relative
289 // weigths.
290 //
291
292 Double_t* da = CreateDoubleArray(a, TMath::Abs(nlmat));
293 Double_t* dz = CreateDoubleArray(z, TMath::Abs(nlmat));
294 Double_t* dwmat = CreateDoubleArray(wmat, TMath::Abs(nlmat));
295
296 Mixture(kmat, name, da, dz, dens, nlmat, dwmat);
297 for (Int_t i=0; i<nlmat; i++) {
298 a[i] = da[i]; z[i] = dz[i]; wmat[i] = dwmat[i];
299 }
300
301 delete [] da;
302 delete [] dz;
303 delete [] dwmat;
304}
305
306//_____________________________________________________________________________
307void TFlukaMCGeometry::Mixture(Int_t& kmat, const char* name, Double_t* a, Double_t* z,
308 Double_t dens, Int_t nlmat, Double_t* wmat)
309{
310 //
311 // Defines mixture OR COMPOUND IMAT as composed by
312 // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
313 //
314 // If NLMAT > 0 then wmat contains the proportion by
315 // weights of each basic material in the mixture.
316 //
317 // If nlmat < 0 then WMAT contains the number of atoms
318 // of a given kind into the molecule of the COMPOUND
319 // In this case, WMAT in output is changed to relative
320 // weigths.
321 //
322
323 if (nlmat < 0) {
324 nlmat = - nlmat;
325 Double_t amol = 0;
326 Int_t i;
327 for (i=0;i<nlmat;i++) {
328 amol += a[i]*wmat[i];
329 }
330 for (i=0;i<nlmat;i++) {
331 wmat[i] *= a[i]/amol;
332 }
333 }
334 kmat = gGeoManager->GetListOfMaterials()->GetSize();
335 printf("Mixture %s with %i elem: kmat=%i, dens=%g\n", name, nlmat, kmat, dens);
336 for (Int_t j=0; j<nlmat; j++) printf(" Elem %i: z=%g a=%g w=%g\n",j,z[j],a[j],wmat[j]);
337 gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
338}
339//_____________________________________________________________________________
340Int_t TFlukaMCGeometry::GetMedium() const
341{
342// Get current medium number
343 Int_t imed = 0;
344 TGeoNode *node = gGeoManager->GetCurrentNode();
345 if (!node) imed = gGeoManager->GetTopNode()->GetVolume()->GetMedium()->GetId();
346 else imed = node->GetVolume()->GetMedium()->GetId();
347 printf("GetMedium=%i\n", imed);
348 return imed;
349}
350
0bc73b6c 351//_____________________________________________________________________________
352Int_t TFlukaMCGeometry::GetFlukaMaterial(Int_t imed) const
353{
354// Returns FLUKA material index for medium IMED
355 TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1);
356 if (!med) {
357 Error("GetFlukaMaterial", "MEDIUM %i nor found", imed);
358 return -1;
359 }
360 Int_t imatfl = med->GetMaterial()->GetIndex();
361 return imatfl;
362}
363
b1536e91 364//_____________________________________________________________________________
365Int_t *TFlukaMCGeometry::GetRegionList(Int_t imed, Int_t &nreg)
366{
367// Get an ordered list of regions matching a given medium number
368 nreg = 0;
369 if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
370 TIter next(gGeoManager->GetListOfUVolumes());
371 TGeoVolume *vol;
372 Int_t imedium, ireg;
373 while ((vol = (TGeoVolume*)next())) {
374 imedium = vol->GetMedium()->GetId();
375 if (imedium == imed) {
376 ireg = vol->GetNumber();
377 fRegionList[nreg++] = ireg;
378 }
379 }
380 return fRegionList;
381}
382
383//_____________________________________________________________________________
384Int_t *TFlukaMCGeometry::GetMaterialList(Int_t imat, Int_t &nreg)
385{
386// Get an ordered list of regions matching a given medium number
387 nreg = 0;
388 if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
389 TIter next(gGeoManager->GetListOfUVolumes());
390 TGeoVolume *vol;
391 Int_t imaterial, ireg;
392 while ((vol = (TGeoVolume*)next())) {
393 imaterial = vol->GetMedium()->GetMaterial()->GetIndex();
394 if (imaterial == imat) {
395 ireg = vol->GetNumber();
396 fRegionList[nreg++] = ireg;
397 }
398 }
399 return fRegionList;
400}
8495a208 401//_____________________________________________________________________________
402void TFlukaMCGeometry::Medium(Int_t& kmed, const char* name, Int_t nmat, Int_t isvol,
403 Int_t ifield, Double_t fieldm, Double_t tmaxfd,
404 Double_t stemax, Double_t deemax, Double_t epsil,
405 Double_t stmin, Float_t* ubuf, Int_t nbuf)
406{
407 //
408 // kmed tracking medium number assigned
409 // name tracking medium name
410 // nmat material number
411 // isvol sensitive volume flag
412 // ifield magnetic field
413 // fieldm max. field value (kilogauss)
414 // tmaxfd max. angle due to field (deg/step)
415 // stemax max. step allowed
416 // deemax max. fraction of energy lost in a step
417 // epsil tracking precision (cm)
418 // stmin min. step due to continuous processes (cm)
419 //
420 // ifield = 0 if no magnetic field; ifield = -1 if user decision in guswim;
421 // ifield = 1 if tracking performed with g3rkuta; ifield = 2 if tracking
422 // performed with g3helix; ifield = 3 if tracking performed with g3helx3.
423 //
424
425 //printf("Creating mediuma: %s, numed=%d, nmat=%d\n",name,kmed,nmat);
426 Double_t* dubuf = CreateDoubleArray(ubuf, nbuf);
427 Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, epsil,
428 stmin, dubuf, nbuf);
429 delete [] dubuf;
430}
431
432//_____________________________________________________________________________
433void TFlukaMCGeometry::Medium(Int_t& kmed, const char* name, Int_t nmat, Int_t isvol,
434 Int_t ifield, Double_t fieldm, Double_t tmaxfd,
435 Double_t stemax, Double_t deemax, Double_t epsil,
436 Double_t stmin, Double_t* /*ubuf*/, Int_t /*nbuf*/)
437{
438 //
439 // kmed tracking medium number assigned
440 // name tracking medium name
441 // nmat material number
442 // isvol sensitive volume flag
443 // ifield magnetic field
444 // fieldm max. field value (kilogauss)
445 // tmaxfd max. angle due to field (deg/step)
446 // stemax max. step allowed
447 // deemax max. fraction of energy lost in a step
448 // epsil tracking precision (cm)
449 // stmin min. step due to continuos processes (cm)
450 //
451 // ifield = 0 if no magnetic field; ifield = -1 if user decision in guswim;
452 // ifield = 1 if tracking performed with g3rkuta; ifield = 2 if tracking
453 // performed with g3helix; ifield = 3 if tracking performed with g3helx3.
454 //
455
0bc73b6c 456 kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
8495a208 457 gGeoManager->Medium(name,kmed,nmat, isvol, ifield, fieldm, tmaxfd, stemax,deemax, epsil, stmin);
458 printf("Medium %s: kmed=%i, nmat=%i, isvol=%i\n", name, kmed, nmat,isvol);
459}
460
461//_____________________________________________________________________________
462void TFlukaMCGeometry::Matrix(Int_t& krot, Double_t thex, Double_t phix, Double_t they,
463 Double_t phiy, Double_t thez, Double_t phiz)
464{
465 //
466 // krot rotation matrix number assigned
467 // theta1 polar angle for axis i
468 // phi1 azimuthal angle for axis i
469 // theta2 polar angle for axis ii
470 // phi2 azimuthal angle for axis ii
471 // theta3 polar angle for axis iii
472 // phi3 azimuthal angle for axis iii
473 //
474 // it defines the rotation matrix number irot.
475 //
476
477 krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
478 gGeoManager->Matrix(krot, thex, phix, they, phiy, thez, phiz);
479 printf("Rotation %i defined\n", krot);
480}
481
482//_____________________________________________________________________________
483Int_t TFlukaMCGeometry::Gsvolu(const char *name, const char *shape, Int_t nmed,
484 Float_t *upar, Int_t npar)
485{
486 //
487 // NAME Volume name
488 // SHAPE Volume type
489 // NUMED Tracking medium number
490 // NPAR Number of shape parameters
491 // UPAR Vector containing shape parameters
492 //
493 // It creates a new volume in the JVOLUM data structure.
494 //
495
496 Double_t* dupar = CreateDoubleArray(upar, npar);
497 Int_t id = Gsvolu(name, shape, nmed, dupar, npar);
498 delete [] dupar;
499 return id;
500}
501
502//_____________________________________________________________________________
503Int_t TFlukaMCGeometry::Gsvolu(const char *name, const char *shape, Int_t nmed,
504 Double_t *upar, Int_t npar)
505{
506 //
507 // NAME Volume name
508 // SHAPE Volume type
509 // NUMED Tracking medium number
510 // NPAR Number of shape parameters
511 // UPAR Vector containing shape parameters
512 //
513 // It creates a new volume in the JVOLUM data structure.
514 //
515 char vname[5];
516 Vname(name,vname);
517 char vshape[5];
518 Vname(shape,vshape);
519
520 TGeoVolume* vol = gGeoManager->Volume(vname, shape, nmed, upar, npar);
521 printf("Volume %s: id=%i shape=%s, nmed=%i\n", vname, vol->GetNumber(), shape, nmed);
522 return vol->GetNumber();
523}
524
525//_____________________________________________________________________________
526void TFlukaMCGeometry::Gsdvn(const char *name, const char *mother, Int_t ndiv,
527 Int_t iaxis)
528{
529 //
530 // Create a new volume by dividing an existing one
531 //
532 // NAME Volume name
533 // MOTHER Mother volume name
534 // NDIV Number of divisions
535 // IAXIS Axis value
536 //
537 // X,Y,Z of CAXIS will be translated to 1,2,3 for IAXIS.
538 // It divides a previously defined volume.
539 //
540 char vname[5];
541 Vname(name,vname);
542 char vmother[5];
543 Vname(mother,vmother);
544
545 gGeoManager->Division(vname, vmother, iaxis, ndiv, 0, 0, 0, "n");
546 printf("Division %s: mother=%s iaxis=%i ndiv=%i\n", vname, vmother, iaxis, ndiv);
547}
548
549//_____________________________________________________________________________
550void TFlukaMCGeometry::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
551 Int_t iaxis, Double_t c0i, Int_t numed)
552{
553 //
554 // Create a new volume by dividing an existing one
555 //
556 // Divides mother into ndiv divisions called name
557 // along axis iaxis starting at coordinate value c0.
558 // the new volume created will be medium number numed.
559 //
560 char vname[5];
561 Vname(name,vname);
562 char vmother[5];
563 Vname(mother,vmother);
564
565 gGeoManager->Division(vname, vmother, iaxis, ndiv, c0i, 0, numed, "nx");
566}
567//_____________________________________________________________________________
568void TFlukaMCGeometry::Gsdvt(const char *name, const char *mother, Double_t step,
569 Int_t iaxis, Int_t numed, Int_t /*ndvmx*/)
570{
571 //
572 // Create a new volume by dividing an existing one
573 //
574 // Divides MOTHER into divisions called NAME along
575 // axis IAXIS in steps of STEP. If not exactly divisible
576 // will make as many as possible and will centre them
577 // with respect to the mother. Divisions will have medium
578 // number NUMED. If NUMED is 0, NUMED of MOTHER is taken.
579 // NDVMX is the expected maximum number of divisions
580 // (If 0, no protection tests are performed)
581 //
582 char vname[5];
583 Vname(name,vname);
584 char vmother[5];
585 Vname(mother,vmother);
586
587 gGeoManager->Division(vname, vmother, iaxis, 0, 0, step, numed, "s");
588}
589
590//_____________________________________________________________________________
591void TFlukaMCGeometry::Gsdvt2(const char *name, const char *mother, Double_t step,
592 Int_t iaxis, Double_t c0, Int_t numed, Int_t /*ndvmx*/)
593{
594 //
595 // Create a new volume by dividing an existing one
596 //
597 // Divides MOTHER into divisions called NAME along
598 // axis IAXIS starting at coordinate value C0 with step
599 // size STEP.
600 // The new volume created will have medium number NUMED.
601 // If NUMED is 0, NUMED of mother is taken.
602 // NDVMX is the expected maximum number of divisions
603 // (If 0, no protection tests are performed)
604 //
605 char vname[5];
606 Vname(name,vname);
607 char vmother[5];
608 Vname(mother,vmother);
609
610 gGeoManager->Division(vname, vmother, iaxis, 0, c0, step, numed, "sx");
611}
612
613//_____________________________________________________________________________
614void TFlukaMCGeometry::Gsord(const char * /*name*/, Int_t /*iax*/)
615{
616 //
617 // Flags volume CHNAME whose contents will have to be ordered
618 // along axis IAX, by setting the search flag to -IAX
619 // IAX = 1 X axis
620 // IAX = 2 Y axis
621 // IAX = 3 Z axis
622 // IAX = 4 Rxy (static ordering only -> GTMEDI)
623 // IAX = 14 Rxy (also dynamic ordering -> GTNEXT)
624 // IAX = 5 Rxyz (static ordering only -> GTMEDI)
625 // IAX = 15 Rxyz (also dynamic ordering -> GTNEXT)
626 // IAX = 6 PHI (PHI=0 => X axis)
627 // IAX = 7 THETA (THETA=0 => Z axis)
628 //
629
630 // TBC - keep this function
631 // nothing to be done for TGeo //xx
632}
633
634//_____________________________________________________________________________
635void TFlukaMCGeometry::Gspos(const char *name, Int_t nr, const char *mother, Double_t x,
636 Double_t y, Double_t z, Int_t irot, const char *konly)
637{
638 //
639 // Position a volume into an existing one
640 //
641 // NAME Volume name
642 // NUMBER Copy number of the volume
643 // MOTHER Mother volume name
644 // X X coord. of the volume in mother ref. sys.
645 // Y Y coord. of the volume in mother ref. sys.
646 // Z Z coord. of the volume in mother ref. sys.
647 // IROT Rotation matrix number w.r.t. mother ref. sys.
648 // ONLY ONLY/MANY flag
649 //
650 // It positions a previously defined volume in the mother.
651 //
652
653 TString only = konly;
654 only.ToLower();
655 Bool_t isOnly = kFALSE;
656 if (only.Contains("only")) isOnly = kTRUE;
657 char vname[5];
658 Vname(name,vname);
659 char vmother[5];
660 Vname(mother,vmother);
661
662 Double_t *upar=0;
663 gGeoManager->Node(vname, nr, vmother, x, y, z, irot, isOnly, upar);
664 printf("Adding daughter %s to %s: cpy=%i irot=%i only=%s\n", vname,vmother,nr,irot,only.Data());
665}
666
667//_____________________________________________________________________________
668void TFlukaMCGeometry::Gsposp(const char *name, Int_t nr, const char *mother,
669 Double_t x, Double_t y, Double_t z, Int_t irot,
670 const char *konly, Float_t *upar, Int_t np )
671{
672 //
673 // Place a copy of generic volume NAME with user number
674 // NR inside MOTHER, with its parameters UPAR(1..NP)
675 //
676
677 Double_t* dupar = CreateDoubleArray(upar, np);
678 Gsposp(name, nr, mother, x, y, z, irot, konly, dupar, np);
679 delete [] dupar;
680}
681
682//_____________________________________________________________________________
683void TFlukaMCGeometry::Gsposp(const char *name, Int_t nr, const char *mother,
684 Double_t x, Double_t y, Double_t z, Int_t irot,
685 const char *konly, Double_t *upar, Int_t np )
686{
687 //
688 // Place a copy of generic volume NAME with user number
689 // NR inside MOTHER, with its parameters UPAR(1..NP)
690 //
691
692 TString only = konly;
693 only.ToLower();
694 Bool_t isOnly = kFALSE;
695 if (only.Contains("only")) isOnly = kTRUE;
696 char vname[5];
697 Vname(name,vname);
698 char vmother[5];
699 Vname(mother,vmother);
700
701 gGeoManager->Node(vname,nr,vmother, x,y,z,irot,isOnly,upar,np);
702 printf("Adding daughter(s) %s to %s: cpy=%i irot=%i only=%s\n", vname,vmother,nr,irot,only.Data());
703}
704
705//_____________________________________________________________________________
706Int_t TFlukaMCGeometry::VolId(const Text_t *name) const
707{
708 //
709 // Return the unique numeric identifier for volume name
710 //
711
712 Int_t uid = gGeoManager->GetUID(name);
713 if (uid<0) {
714 printf("VolId: Volume %s not found\n",name);
715 return 0;
716 }
717 printf("VolId for %s: %i\n", name, uid);
718 return uid;
719}
720
721//_____________________________________________________________________________
722const char* TFlukaMCGeometry::VolName(Int_t id) const
723{
724 //
725 // Return the volume name given the volume identifier
726 //
727 TGeoVolume *volume = gGeoManager->GetVolume(id);
728 if (!volume) {
729 Error("VolName","volume with id=%d does not exist",id);
730 return "NULL";
731 }
732 printf("VolName for id=%i: %s\n", id, volume->GetName());
733 return volume->GetName();
734}
735
736//_____________________________________________________________________________
737Int_t TFlukaMCGeometry::NofVolumes() const
738{
739 //
740 // Return total number of volumes in the geometry
741 //
742
743 return gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1;
744}
745
746//_____________________________________________________________________________
747Int_t TFlukaMCGeometry::VolId2Mate(Int_t id) const
748{
749 //
750 // Return material number for a given volume id
751 //
752 TGeoVolume *volume = gGeoManager->GetVolume(id);
753 if (!volume) {
754 Error("VolId2Mate","volume with id=%d does not exist",id);
755 return 0;
756 }
757 TGeoMedium *med = volume->GetMedium();
758 if (!med) return 0;
759 printf("VolId2Mate id=%i: idmed=%i\n", id, med->GetId());
760 return med->GetId();
761}
762
763//_____________________________________________________________________________
764Int_t TFlukaMCGeometry::CurrentVolID(Int_t& copyNo) const
765{
766 // Returns the current volume ID and copy number
767 if (gGeoManager->IsOutside()) return 0;
768 TGeoNode *node = gGeoManager->GetCurrentNode();
769 copyNo = node->GetNumber();
770 Int_t id = node->GetVolume()->GetNumber();
771 printf("CurrentVolId(cpy=%i) = %i\n", copyNo, id);
772 return id;
773}
774
775//_____________________________________________________________________________
776Int_t TFlukaMCGeometry::CurrentVolOffID(Int_t off, Int_t& copyNo) const
777{
778 // Return the current volume "off" upward in the geometrical tree
779 // ID and copy number
780 if (off<0 || off>gGeoManager->GetLevel()) return 0;
781 if (off==0) return CurrentVolID(copyNo);
782 TGeoNode *node = gGeoManager->GetMother(off);
783 if (!node) return 0;
784 copyNo = node->GetNumber();
785 printf("CurrentVolOffId(off=%i,cpy=%i) = %i\n", off,copyNo,node->GetVolume()->GetNumber() );
786 return node->GetVolume()->GetNumber();
787}
788// FLUKA specific
789
790//_____________________________________________________________________________
791const char* TFlukaMCGeometry::CurrentVolName() const
792{
793 //
794 // Returns the current volume name
795 //
796 if (gGeoManager->IsOutside()) return 0;
797 printf("CurrentVolName : %s\n", gGeoManager->GetCurrentVolume()->GetName());
798 return gGeoManager->GetCurrentVolume()->GetName();
799}
800//_____________________________________________________________________________
801const char* TFlukaMCGeometry::CurrentVolOffName(Int_t off) const
802{
803 //
804 // Return the current volume "off" upward in the geometrical tree
805 // ID, name and copy number
806 // if name=0 no name is returned
807 //
808 if (off<0 || off>gGeoManager->GetLevel()) return 0;
809 if (off==0) return CurrentVolName();
810 TGeoNode *node = gGeoManager->GetMother(off);
811 if (!node) return 0;
812 printf("CurrentVolOffName(off=%i) : %s\n", off,node->GetVolume()->GetName());
813 return node->GetVolume()->GetName();
814}
815
816//_____________________________________________________________________________
817void TFlukaMCGeometry::Gsatt(const char *name, const char *att, Int_t val)
818{
819 //
820 // NAME Volume name
821 // IOPT Name of the attribute to be set
822 // IVAL Value to which the attribute is to be set
823 // see: TFluka::Gsatt
824 char vname[5];
825 Vname(name,vname);
826 char vatt[5];
827 Vname(att,vatt);
828 gGeoManager->SetVolumeAttribute(vname, vatt, val);
829}
830
831//_____________________________________________________________________________
832void TFlukaMCGeometry::Gdtom(Float_t *xd, Float_t *xm, Int_t iflag)
833{
834 //
835 // Computes coordinates XM (Master Reference System
836 // knowing the coordinates XD (Detector Ref System)
837 // The local reference system can be initialized by
838 // - the tracking routines and GDTOM used in GUSTEP
839 // - a call to GSCMED(NLEVEL,NAMES,NUMBER)
840 // (inverse routine is GMTOD)
841 //
842 // If IFLAG=1 convert coordinates
843 // IFLAG=2 convert direction cosinus
844 //
845 Double_t XM[3], XD[3];
846 Int_t i;
847 for (i=0;i<3;i++) XD[i] = xd[i];
848 if (iflag == 1) gGeoManager->LocalToMaster(XD,XM);
849 else gGeoManager->LocalToMasterVect(XD,XM);
850 for (i=0;i<3;i++) xm[i]=XM[i];
851}
852
853//_____________________________________________________________________________
854void TFlukaMCGeometry::Gdtom(Double_t *xd, Double_t *xm, Int_t iflag)
855{
856 if (iflag == 1) gGeoManager->LocalToMaster(xd,xm);
857 else gGeoManager->LocalToMasterVect(xd,xm);
858}
859
860//_____________________________________________________________________________
861void TFlukaMCGeometry::Gmtod(Float_t *xm, Float_t *xd, Int_t iflag)
862{
863 //
864 // Computes coordinates XD (in DRS)
865 // from known coordinates XM in MRS
866 // The local reference system can be initialized by
867 // - the tracking routines and GMTOD used in GUSTEP
868 // - a call to GMEDIA(XM,NUMED,CHECK)
869 // - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
870 // (inverse routine is GDTOM)
871 //
872 // If IFLAG=1 convert coordinates
873 // IFLAG=2 convert direction cosinus
874 //
875 Double_t XM[3], XD[3];
876 Int_t i;
877 for (i=0;i<3;i++) XM[i]=xm[i];
878 if (iflag == 1) gGeoManager->MasterToLocal(XM,XD);
879 else gGeoManager->MasterToLocalVect(XM,XD);
880 for (i=0;i<3;i++) xd[i] = XD[i];
881}
882
883//_____________________________________________________________________________
884void TFlukaMCGeometry::Gmtod(Double_t *xm, Double_t *xd, Int_t iflag)
885{
886 if (iflag == 1) gGeoManager->MasterToLocal(xm,xd);
887 else gGeoManager->MasterToLocalVect(xm,xd);
888}
889
890//_____________________________________________________________________________
efde9b4d 891void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname)
8495a208 892{
893 // ==== from FLUGG ====
894 // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case,
895 // according to the fluka standard. In addition,. they must be equal to the
896 // names of the fluka materials - see fluka manual - in order that the
897 // program load the right cross sections, and equal to the names included in
898 // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his
899 // own .pemf, in order to get the right cross sections loaded in memory.
900
901 Int_t zelem[128];
902 static char elNames[220] = {
903 // 1 ============================= 5 ==================================== 10 ===================================== 15 ===
904 'H','_','H','E','L','I','B','E','B','_','C','_','N','_','O','_','F','_','N','E','N','A','M','G','A','L','S','I','P','_',
905 'S','_','C','L','A','R','K','_','C','A','S','C','T','I','V','_','C','R','M','N','F','E','C','O','N','I','C','U','Z','N',
906 'G','A','G','E','A','S','S','E','B','R','K','R','R','B','S','R','Y','_','Z','R','N','B','M','O','T','C','R','U','R','H',
907 'P','D','A','G','C','D','I','N','S','N','S','B','T','E','I','_','X','E','C','S','B','A','L','A','C','E','P','R','N','D',
908 'P','M','S','M','E','U','G','D','T','B','D','Y','H','O','E','R','T','M','Y','B','L','U','H','F','T','A','W','_','R','E',
909 'O','S','I','R','P','T','A','U','H','G','T','L','P','B','B','I','P','O','A','T','R','N','F','R','R','A','A','C','T','H',
910 'P','A','U','_','N','P','P','U','A','M','C','M','B','K','C','F','E','S','F','M','M','D','N','O','L','R','R','F','D','B',
911 'S','G','B','H','H','S','M','T','D','S'};
912 memset(zelem, 0, 128*sizeof(Int_t));
913 TString sname;
914 gGeoManager->Export("flgeom.root");
915 if (fname) sname = fname;
916 else sname = "flukaMat.inp";
917 ofstream out;
918 out.open(sname.Data(), ios::out);
919 if (!out.good()) {
920 Fatal("CreateFlukaMatFile", "could not open file %s for writing", sname.Data());
921 return;
922 }
923 PrintHeader(out, "MATERIALS AND COMPOUNDS");
924 PrintHeader(out, "MATERIALS");
925 TList *matlist = gGeoManager->GetListOfMaterials();
926 TIter next(matlist);
927 Int_t nmater = matlist->GetSize();
928 Int_t nfmater = 0;
929 TObjArray *listfluka = new TObjArray(nmater+50);
930 TObjArray *listflukanames = new TObjArray(nmater+50);
931 TGeoMaterial *mat, *matorig;
932 TGeoMixture *mix = 0;
933 TString matname;
934 TObjString *objstr, *objstrother;
935 Int_t i,j,k,idmat;
936 Bool_t done;
937 Int_t nelem, nidmat;
938 Double_t amat,zmat,rhomat;
939 Double_t zel, ael, wel, rho;
940 char elname[8] = {' ',' ','_', 'E','L','E','M','\0'};
941 char digit[3];
efde9b4d 942 Bool_t found = kFALSE;
8495a208 943
944 printf("Creating materials and compounds\n");
945 for (i=0; i<nmater; i++) {
946 mat = (TGeoMaterial*)matlist->At(i);
efde9b4d 947 if (mat->GetZ()<1E-1) {
948 mat->SetIndex(2); // vacuum, built-in inside FLUKA
949 continue;
950 }
951// printf("material: %s index=%i: Z=%f A=%f rho=%f\n", mat->GetName(), mat->GetIndex(),mat->GetZ(),mat->GetA(),mat->GetDensity());
8495a208 952 matorig = gGeoManager->FindDuplicateMaterial(mat);
953 if (matorig) {
954 idmat = matorig->GetIndex();
955 mat->SetIndex(idmat);
efde9b4d 956// printf(" -> found a duplicate: %s with index %i\n", matorig->GetName(), idmat);
8495a208 957 matorig = 0;
958 } else {
efde9b4d 959// printf(" Adding to temp list with index %i\n", nfmater+3);
8495a208 960 listfluka->Add(mat);
961 mat->SetIndex(nfmater+3);
962 matorig = mat;
963 objstr = new TObjString(mat->GetName());
964 listflukanames->Add(objstr);
965 nfmater++;
966 // look if name is existing
967 nidmat = 0;
968 matname = objstr->GetString();
969 ToFlukaString(matname);
970 objstr->SetString(matname.Data());
971 done = kFALSE;
972 while (!done) {
973 if (nfmater == 1) break;
974 for (j=0; j<nfmater-1; j++) {
975 objstrother = (TObjString*)listflukanames->At(j);
976 if (objstr->IsEqual(objstrother)) {
977 // we have to change the name
978 if (nidmat>98) {
979 Error("CreateFlukaMatFile", "too many materials having same name");
980 return;
981 }
982 nidmat++;
983 k = matname.Index(" ");
984 if (k<0 || k>6) k=6;
985 if (nidmat>9) {
986 sprintf(digit, "%d", nidmat);
987 } else {
988 digit[0] = '0';
989 sprintf(&digit[1], "%d", nidmat);
990 }
991 matname.Insert(k,digit);
992 matname.Remove(8);
993 objstr->SetString(matname.Data());
994 break;
995 }
996 if (j == nfmater-2) {
997 done = kTRUE;
998 break;
999 }
1000 }
1001 }
efde9b4d 1002// printf(" newmat name: %s\n", matname.Data());
8495a208 1003 }
1004 // now we have unique materials with unique names in the lists
1005
8495a208 1006 if (matorig && matorig->IsMixture()) {
1007 // create dummy materials for elements
1008 rho = 0.999;
1009 mix = (TGeoMixture*)matorig;
1010 nelem = mix->GetNelements();
efde9b4d 1011// printf(" material is a MIXTURE with %i elements:\n", nelem);
8495a208 1012 for (j=0; j<nelem; j++) {
efde9b4d 1013 found = kFALSE;
8495a208 1014 zel = (mix->GetZmixt())[j];
efde9b4d 1015 ael = (mix->GetAmixt())[j];
1016// printf(" Zelem[%i] = %g\n",j,zel);
8495a208 1017 if ((zel-Int_t(zel))>0.01) {
efde9b4d 1018 TGeoMaterial *mat1;
1019 for (Int_t imat=0; imat<nfmater; imat++) {
1020 mat1 = (TGeoMaterial*)listfluka->At(imat);
1021 if (TMath::Abs(mat1->GetZ()-zel)>1E-4) continue;
1022 if (TMath::Abs(mat1->GetA()-ael)>1E-4) continue;
1023 found = kTRUE;
1024 break;
1025 }
1026 if (!found) Warning("CreateFlukaMatFile", "element with Z=%f\n", zel);
8495a208 1027 }
efde9b4d 1028 if (!zelem[Int_t(zel)] && !found) {
8495a208 1029 // write fluka element
1030 memcpy(elname, &elNames[2*Int_t(zel-1)], 2);
1031 zelem[Int_t(zel)] = 1;
8495a208 1032 mat = new TGeoMaterial(elname, ael, zel, rho);
1033 mat->SetIndex(nfmater+3);
efde9b4d 1034// printf(" element not in list: new material %s at index=%i, Z=%g, A=%g, dummyrho=%g\n",
1035// elname,nfmater+3,zel,ael,rho);
8495a208 1036 listfluka->Add(mat);
1037 objstr = new TObjString(elname);
1038 listflukanames->Add(objstr);
1039 nfmater++;
1040 }
1041 }
1042 }
1043 }
1044 // now dump materials in the file
efde9b4d 1045// printf("DUMPING %i materials\n", nfmater);
8495a208 1046 for (i=0; i<nfmater; i++) {
1047 mat = (TGeoMaterial*)listfluka->At(i);
1048 out << setw(10) << "MATERIAL ";
1049 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
b0853588 1050// matname = mat->GetName();
1051 objstr = (TObjString*)listflukanames->At(i);
1052 matname = objstr->GetString();
8495a208 1053 ToFlukaString(matname);
1054 zmat = mat->GetZ();
efde9b4d 1055 if (zmat-Int_t(zmat)>0.01) {
1056 if (zmat-Int_t(zmat)>0.5) zmat = Int_t(zmat)+1.;
1057 else zmat = Int_t(zmat);
1058 }
8495a208 1059 amat = mat->GetA();
1060 rhomat = mat->GetDensity();
1061 // write material card
1062 if (mat->IsMixture()) {
1063 out << setw(10) << " ";
1064 out << setw(10) << " ";
1065 mix = (TGeoMixture*)mat;
1066 } else {
1067 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << zmat;
1068 out << setw(10) << setprecision(3) << amat;
1069 }
1070 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
1071 out << setw(10) << setiosflags(ios::scientific) << setprecision(3) << rhomat;
1072 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
1073 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(i+3);
1074 out << setw(10) << " ";
efde9b4d 1075 out << setw(10) << " ";
1076 out << setw(8) << matname.Data() << endl;
8495a208 1077 }
1078 // write mixture header
1079 PrintHeader(out, "COMPOUNDS");
1080 Int_t counttothree;
1081 TGeoMaterial *element;
1082 for (i=0; i<nfmater; i++) {
1083 mat = (TGeoMaterial*)listfluka->At(i);
1084 if (!mat->IsMixture()) continue;
1085 mix = (TGeoMixture*)mat;
1086 counttothree = 0;
1087 out << setw(10) << "COMPOUND ";
1088 nelem = mix->GetNelements();
1089 objstr = (TObjString*)listflukanames->At(i);
1090 matname = objstr->GetString();
efde9b4d 1091// printf("MIXTURE %s with index %i having %i elements\n", matname.Data(), mat->GetIndex(),nelem);
8495a208 1092 for (j=0; j<nelem; j++) {
1093 // dump mixture cards
efde9b4d 1094// printf(" #elem %i: Z=%g, A=%g, W=%g\n", j, (mix->GetZmixt())[j],
1095// (mix->GetAmixt())[j],(mix->GetWmixt())[j]);
8495a208 1096 wel = (mix->GetWmixt())[j];
1097 zel = (mix->GetZmixt())[j];
efde9b4d 1098 ael = (mix->GetAmixt())[j];
1099 if (zel-Int_t(zel)>0.01) {
1100 // loop the temporary list
1101 element = 0;
1102 TGeoMaterial *mat1;
1103 for (Int_t imat=0; imat<i; imat++) {
1104 mat1 = (TGeoMaterial*)listfluka->At(imat);
1105 if (TMath::Abs(mat1->GetZ()-zel)>1E-4) continue;
1106 if (TMath::Abs(mat1->GetA()-ael)>1E-4) continue;
1107 element = mat1;
1108 break;
1109 }
1110 } else {
1111 memcpy(elname, &elNames[2*Int_t(zel-1)], 2);
1112 element = (TGeoMaterial*)listfluka->FindObject(elname);
1113 }
8495a208 1114 if (!element) {
1115 Error("CreateFlukaMatFile", "Element Z=%g %s not found", zel, elname);
1116 return;
1117 }
1118 idmat = element->GetIndex();
efde9b4d 1119// printf("element %s , index=%i\n", element->GetName(), idmat);
8495a208 1120 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
1121 out << setw(10) << setiosflags(ios::fixed) << setprecision(6) << -wel;
1122 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
1123 out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(idmat);
1124 counttothree++;
1125 if (counttothree == 3) {
1126 out << matname.Data();
1127 out << endl;
1128 if ( (j+1) != nelem) out << setw(10) << "COMPOUND ";
1129 counttothree = 0;
1130 }
1131 }
1132 //Unless we have 3, 6, 9... submaterials we need to put some empty
1133 //space and the compound name
1134 if (nelem%3) {
1135 for (j=0; j<(3-(nelem%3)); j++)
1136 out << setw(10) << " " << setw(10) << " ";
1137 out << matname.Data();
1138 out << endl;
1139 }
1140 }
1141
1142 // Now print the list of regions (volumes in TGeo)
1143 Int_t nvols = gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1;
1144 TGeoVolume *vol;
1145/*
1146 PrintHeader(out, "TGEO VOLUMES");
1147 for (i=1; i<=nvols; i++) {
1148 vol = gGeoManager->GetVolume(i);
1149 out.setf(std::ios::left, std::ios::adjustfield);
1150 out << setw(10) << i;
1151 out << setw(20) << vol->GetName() << endl;
1152 }
1153*/
1154 // Now print the material assignments
1155 Double_t flagfield;
1156 PrintHeader(out, "TGEO MATERIAL ASSIGNMENTS");
1157 for (i=1; i<=nvols; i++) {
1158 vol = gGeoManager->GetVolume(i);
1159 mat = vol->GetMedium()->GetMaterial();
1160 idmat = mat->GetIndex();
b0853588 1161// flagfield = (vol->GetField())?1.:0.;
05265ca9 1162 flagfield = 1.;
8495a208 1163 out << setw(10) << "ASSIGNMAT ";
1164 out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
1165 out << setw(10) << setiosflags(ios::fixed) << Double_t(idmat);
1166 out << setw(10) << setiosflags(ios::fixed) << Double_t(i);
1167 out << setw(10) << "0.0";
b0853588 1168 out << setw(10) << "0.0";
8495a208 1169 out << setw(10) << setiosflags(ios::fixed) << flagfield;
b0853588 1170 out << setw(10) << "0.0";
8495a208 1171 out << endl;
1172 }
1173 delete listfluka;
1174 listflukanames->Delete();
1175 delete listflukanames;
1176 out.close();
efde9b4d 1177 fLastMaterial = nfmater+2;
8495a208 1178}
1179
1180//_____________________________________________________________________________
1181void TFlukaMCGeometry::PrintHeader(ofstream &out, const char *text) const
1182{
1183// Print a FLUKA header.
1184 out << "*\n" << "*\n" << "*\n";
1185 out << "********************* " << text << " *********************\n"
1186 << "*\n";
1187 out << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..."
1188 << endl;
1189 out << "*" << endl;
1190}
1191
1192//_____________________________________________________________________________
1193Int_t TFlukaMCGeometry::RegionId() const
1194{
1195// Returns current region id <-> TGeo node id
1196 if (gGeoManager->IsOutside()) return 0;
1197 return gGeoManager->GetCurrentNode()->GetUniqueID();
1198}
05265ca9 1199//_____________________________________________________________________________
1200void TFlukaMCGeometry::SetMreg(Int_t mreg)
1201{
1202// Update if needed next history;
1203 Int_t curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
1204 if (mreg==curreg) return;
1205 if (mreg==fNextRegion) {
1206 if (fNextLattice!=999999999) gGeoManager->CdNode(fNextLattice-1);
1207 return;
1208 }
1209 printf("ERROR: mreg=%i neither current nor next region\n", mreg);
1210}
1211
1212//_____________________________________________________________________________
1213void TFlukaMCGeometry::SetNextRegion(Int_t mreg, Int_t latt)
1214{
1215// Set index/history for next entered region
1216 fNextRegion = mreg;
1217 fNextLattice = latt;
1218}
8495a208 1219
1220//_____________________________________________________________________________
1221void TFlukaMCGeometry::ToFlukaString(TString &str) const
1222{
1223// ToFlukaString converts an string to something usefull in FLUKA:
1224// * Capital letters
1225// * Only 8 letters
1226// * Replace ' ' by '_'
1227 if (str.Length()<8) {
1228 str += " ";
1229 }
1230 str.Remove(8);
1231 Int_t ilast;
1232 for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break;
1233 str.ToUpper();
1234 for (Int_t pos=0; pos<ilast; pos++)
1235 if (str(pos)==' ') str.Replace(pos,1,"_",1);
1236 return;
1237}
1238//______________________________________________________________________________
1239void TFlukaMCGeometry::Vname(const char *name, char *vname) const
1240{
1241 //
1242 // convert name to upper case. Make vname at least 4 chars
1243 //
1244 Int_t l = strlen(name);
1245 Int_t i;
1246 l = l < 4 ? l : 4;
1247 for (i=0;i<l;i++) vname[i] = toupper(name[i]);
1248 for (i=l;i<4;i++) vname[i] = ' ';
1249 vname[4] = 0;
1250}
1251
1252
1253// FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers
1254
8495a208 1255//_____________________________________________________________________________
1256Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/)
1257{
1258// from FLUGG:
1259// Wrapper for setting DNEAR option on fluka side. Must return 0
1260// if user doesn't want Fluka to use DNEAR to compute the
1261// step (the same effect is obtained with the GLOBAL (WHAT(3)=-1)
1262// card in fluka input), returns 1 if user wants Fluka always to
1263// use DNEAR (in this case, be sure that GEANT4 DNEAR is unique,
1264// coming from all directions!!!)
2573ac89 1265 printf("========== Dummy IDNRWR\n");
8495a208 1266 return 0;
1267}
1268
8495a208 1269//_____________________________________________________________________________
1270void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1271 Double_t *pV, Int_t &oldReg , const Int_t &oldLttc, Double_t & propStep,
1272 Int_t & /*nascFlag*/, Double_t &retStep, Int_t &newReg,
2573ac89 1273 Double_t &saf, Int_t &newLttc, Int_t &lttcFlag,
8495a208 1274 Double_t *sLt, Int_t *jrLt)
1275{
1276// from FLUGG:
1277// Wrapper for geometry tracking: returns approved step of
1278// particle and all variables that fluka G1 computes.
1279
1280 // Initialize current point/direction
2573ac89 1281 printf("========== Inside G1WR\n");
05265ca9 1282 printf(" point/dir:(%14.9f, %14.9f, %14.9f, %g, %g, %g)\n", pSx,pSy,pSz,pV[0],pV[1],pV[2]);
8495a208 1283 gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
1284 gGeoManager->SetCurrentDirection(pV);
2573ac89 1285 printf(" oldReg=%i oldLttc=%i pstep=%f\n",oldReg, oldLttc, propStep);
05265ca9 1286 if (oldLttc==999999999) printf("WOOPS - wrong old lattice\n");
1287 if (gGeoManager->IsOutside()) {
1288 gGeoManager->SetOutside(kFALSE);
1289 gGeoManager->CdTop();
1290 }
2573ac89 1291 Int_t curLttc = gGeoManager->GetCurrentNodeId()+1;
05265ca9 1292 Int_t curreg = gGeoManager->GetCurrentVolume()->GetNumber();
2573ac89 1293 printf(" curReg=%i curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath());
1294 Bool_t regsame = (curreg==oldReg)?kTRUE:kFALSE;
1295 if (!regsame) printf(" REGIONS DOES NOT MATCH\n");
1296 if (oldLttc != curLttc) {
1297 printf(" HISTORIES DOES NOT MATCH\n");
1298 gGeoManager->CdNode(oldLttc-1);
1299 curLttc = gGeoManager->GetCurrentNodeId()+1;
05265ca9 1300 curreg = gGeoManager->GetCurrentVolume()->GetNumber();
2573ac89 1301 printf(" re-initialized point: curReg=%i curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath());
1302 }
1303 lttcFlag = 0;
8495a208 1304 sLt[lttcFlag] = 0.;
2573ac89 1305 jrLt[lttcFlag] = curLttc;
8495a208 1306 // now 'oldregion' contains the real region, matching or not the old history
1307
1308 // Compute geometry step/safety within physical step limit
2573ac89 1309// newReg = oldregion;
1310 Double_t *point = gGeoManager->GetCurrentPoint();
1311 Double_t *dir = gGeoManager->GetCurrentDirection();
8495a208 1312 Double_t steptot = 0.;
1313 Double_t snext = 0.;
1314 Int_t istep = 0;
1315 Bool_t done = kFALSE;
05265ca9 1316 Double_t pst;
2573ac89 1317 Int_t i;
8495a208 1318 while (!done) {
2573ac89 1319 gGeoManager->FindNextBoundary(-propStep);
8495a208 1320 snext = gGeoManager->GetStep();
2573ac89 1321 printf(" FindNextBoundary(%g) snext=%g\n", propStep, snext);
1322 if (steptot == 0) {
1323 saf = gGeoManager->GetSafeDistance();
1324 printf(" Safety: %g\n", saf);
1325 }
1326 sLt[lttcFlag] = propStep;
1327 jrLt[lttcFlag] = gGeoManager->GetCurrentNodeId()+1;
1328 lttcFlag++; //1
05265ca9 1329 sLt[lttcFlag] = 0.;
1330 jrLt[lttcFlag] = -1;
2573ac89 1331 newReg = curreg;
1332 newLttc = oldLttc;
8495a208 1333 if (snext<propStep) {
1334 // There is a boundary on the way.
1335 // Make a step=snext+1E-6 to force boundary crossing
2573ac89 1336 lttcFlag--; // 0
8495a208 1337 steptot += snext;
1338 sLt[lttcFlag] = snext;
1339 retStep = snext;
2573ac89 1340// lttcFlag++;
1341 // make the step to get into the next region
1342 for (i=0;i<3;i++) point[i]+=(snext+1E-6)*dir[i];
1343 gGeoManager->FindNode();
1344 istep = 0;
1345 printf(" boundary: step made %g\n", snext);
1346 while (gGeoManager->IsSameLocation() && steptot<propStep) {
1347 if (istep>1E3) {
1348 printf("Geometry error: could not cross boundary after extra 10 microns\n");
1349 return;
1350 }
1351 for (i=0;i<3;i++) point[i]+=1E-6*dir[i];
1352 gGeoManager->FindNode();
8495a208 1353 sLt[lttcFlag] += 1E-6;
1354 retStep = sLt[lttcFlag];
1355 steptot += 1E-6;
1356 istep++;
2573ac89 1357 }
2573ac89 1358 if (steptot>propStep) {printf("Error\n");return;}
8495a208 1359 // we managed to cross the boundary -> in which region
1360 newReg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
05265ca9 1361 lttcFlag++; //1
2573ac89 1362 newLttc = (gGeoManager->IsOutside())?999999999:gGeoManager->GetCurrentNodeId()+1;
05265ca9 1363 sLt[lttcFlag] = snext; // at 1
2573ac89 1364 jrLt[lttcFlag] = newLttc;
05265ca9 1365 sLt[lttcFlag+1] = 0.;
1366 jrLt[lttcFlag+1] = -1;
1367 // !!!!!!!!!!
1368
1369 while (newReg==oldReg && steptot<propStep) {
1370 printf(" Entered SAME region... continue\n");
1371 pst = propStep-steptot;
1372 gGeoManager->FindNextBoundary(-pst);
1373 snext = gGeoManager->GetStep();
1374 steptot += snext;
1375 if (snext<pst) {
1376 printf("Found new boundary\n");
1377 sLt[lttcFlag] = snext;
1378 retStep = steptot; // ???
1379 for (i=0;i<3;i++) point[i]+=(snext+1E-6)*dir[i];
1380 steptot += 1E-6;
1381 gGeoManager->FindNode();
1382 if (gGeoManager->IsSameLocation()) {
1383 printf("Cannot cross boundary\n");
1384 break;
1385 }
1386 newReg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
1387 newLttc = (gGeoManager->IsOutside())?999999999:gGeoManager->GetCurrentNodeId()+1;
1388 printf("Found newreg=%i, newLttc=%i, lttFlag is: %i\n", newReg, newLttc, lttcFlag);
1389 sLt[lttcFlag-1] += snext; // correct step in old region
1390 sLt[lttcFlag] = propStep-snext;
1391 jrLt[lttcFlag] = newLttc;
1392 sLt[lttcFlag+1] = 0.;
1393 jrLt[lttcFlag+1] = -1;
1394 if (newReg != oldReg) break; // lttcFlag=1
1395 lttcFlag++;
1396 } else {
1397 printf("Not crossing next\n");
1398 lttcFlag--; //0
1399 retStep=steptot;
1400 sLt[lttcFlag] = retStep;
1401 sLt[lttcFlag+1] = 0.;
1402 jrLt[lttcFlag+1] = -1;
1403 done = kTRUE;
1404 }
1405 }
1406
1407 lttcFlag++; //2
2573ac89 1408 if (!gGeoManager->IsOutside()) {
2573ac89 1409 printf(" ENTERED region %i, newLttc=%i in: %s\n", newReg,newLttc,gGeoManager->GetPath());
1410 } else printf(" EXIT GEOMETRY: BLKHOLE reg=%i\n", newReg);
1411 }
1412 // no boundary within proposed step
1413 lttcFlag--;
1414 done = kTRUE;
1415 }
1416 printf("=> newReg=%i newLttc=%i lttcFlag=%i\n", newReg, newLttc, lttcFlag);
05265ca9 1417 mcgeom->SetNextRegion(newReg, newLttc);
2573ac89 1418 printf("=> snext=%g safe=%g\n", snext, saf);
1419 for (Int_t i=0; i<lttcFlag+1; i++) printf(" jrLt[%i]=%i sLt[%i]=%g\n", i,jrLt[i],i,sLt[i]);
1420 if (newLttc!=oldLttc) {
1421 if (gGeoManager->IsOutside()) {
1422 gGeoManager->SetOutside(kFALSE);
1423 gGeoManager->CdTop();
1424 }
1425 gGeoManager->CdNode(oldLttc-1);
8495a208 1426 }
2573ac89 1427 printf("<= G1WR (in: %s)\n", gGeoManager->GetPath());
8495a208 1428}
1429
efde9b4d 1430//_____________________________________________________________________________
1431void g1rtwr()
1432{
2573ac89 1433 printf("========== Dummy G1RTWR\n");
efde9b4d 1434}
1435
1436//_____________________________________________________________________________
1437void conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/)
1438{
2573ac89 1439 printf("========== Dummy CONHWR\n");
efde9b4d 1440}
1441
1442//_____________________________________________________________________________
2573ac89 1443void inihwr(Int_t &intHist)
efde9b4d 1444{
2573ac89 1445 printf("========== Inside INIHWR -> reinitializing history: %i\n", intHist);
1446 if (gGeoManager->IsOutside()) gGeoManager->CdTop();
1447 if (intHist<=0) {
1448// printf("=== wrong history number\n");
1449 return;
1450 }
1451 if (intHist==0) gGeoManager->CdTop();
1452 else gGeoManager->CdNode(intHist-1);
1453 printf(" --- current path: %s\n", gGeoManager->GetPath());
1454 printf("<= INIHWR\n");
efde9b4d 1455}
1456
1457//_____________________________________________________________________________
1458void jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
1459 Int_t &flukaReg)
1460{
1461// Geometry initialization wrapper called by FLUKAM. Provides to FLUKA the
1462// number of regions (volumes in TGeo)
1463 // build application geometry
2573ac89 1464 printf("========== Inside JOMIWR\n");
1465 flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast();
1466 printf("<= JOMIWR: last region=%i\n", flukaReg);
efde9b4d 1467}
1468
1469//_____________________________________________________________________________
1470void lkdbwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1471 Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
1472 Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
1473{
2573ac89 1474 printf("========== Inside LKDBWR (%f, %f, %f)\n",pSx, pSy, pSz);
efde9b4d 1475// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1476 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1477 TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
1478 if (gGeoManager->IsOutside()) {
05265ca9 1479 printf("OUTSIDE\n");
efde9b4d 1480 newReg = mcgeom->NofVolumes()+1;
05265ca9 1481// newLttc = gGeoManager->GetCurrentNodeId();
1482 newLttc = 999999999;
1483 printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
1484 printf("<= LKMGWR\n");
1485 flagErr = newReg;
1486 return;
efde9b4d 1487 }
1488 newReg = node->GetVolume()->GetNumber();
1489 newLttc = gGeoManager->GetCurrentNodeId()+1;
1490 flagErr = newReg;
1491 printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
2573ac89 1492 printf("<= LKDBWR\n");
efde9b4d 1493}
1494
1495//_____________________________________________________________________________
1496void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1497 Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
1498 Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
1499{
2573ac89 1500 printf("========== Inside LKFXWR (%f, %f, %f)\n",pSx, pSy, pSz);
efde9b4d 1501// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1502 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1503 TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
1504 if (gGeoManager->IsOutside()) {
05265ca9 1505 printf("OUTSIDE\n");
efde9b4d 1506 newReg = mcgeom->NofVolumes()+1;
05265ca9 1507// newLttc = gGeoManager->GetCurrentNodeId();
1508 newLttc = 999999999;
1509 printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
1510 printf("<= LKMGWR\n");
1511 flagErr = newReg;
1512 return;
efde9b4d 1513 }
1514 newReg = node->GetVolume()->GetNumber();
1515 newLttc = gGeoManager->GetCurrentNodeId()+1;
1516 flagErr = newReg;
1517 printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
2573ac89 1518 printf("<= LKFXWR\n");
efde9b4d 1519}
1520
1521//_____________________________________________________________________________
1522void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1523 Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
1524 Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
1525{
2573ac89 1526 printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz);
efde9b4d 1527// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1528 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1529 TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
1530 if (gGeoManager->IsOutside()) {
05265ca9 1531 printf("OUTSIDE\n");
efde9b4d 1532 newReg = mcgeom->NofVolumes()+1;
05265ca9 1533// newLttc = gGeoManager->GetCurrentNodeId();
1534 newLttc = 999999999;
1535 printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
1536 printf("<= LKMGWR\n");
1537 flagErr = newReg;
1538 return;
efde9b4d 1539 }
1540 newReg = node->GetVolume()->GetNumber();
1541 newLttc = gGeoManager->GetCurrentNodeId()+1;
1542 flagErr = newReg;
1543 printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
2573ac89 1544 printf("<= LKMGWR\n");
efde9b4d 1545}
1546
1547//_____________________________________________________________________________
1548void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1549 Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
1550 Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
1551{
2573ac89 1552 printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz);
efde9b4d 1553// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1554 printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1555 TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
1556 if (gGeoManager->IsOutside()) {
05265ca9 1557 printf("OUTSIDE\n");
efde9b4d 1558 newReg = mcgeom->NofVolumes()+1;
05265ca9 1559// newLttc = gGeoManager->GetCurrentNodeId();
1560 newLttc = 999999999;
1561 printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
1562 printf("<= LKMGWR\n");
1563 flagErr = newReg;
1564 return;
efde9b4d 1565 }
1566 newReg = node->GetVolume()->GetNumber();
1567 newLttc = gGeoManager->GetCurrentNodeId()+1;
1568 flagErr = newReg;
2573ac89 1569 printf(" out: newReg=%i newLttc=%i in %s\n", newReg, newLttc, gGeoManager->GetPath());
1570 printf("<= LKWR\n");
efde9b4d 1571}
1572
1573//_____________________________________________________________________________
2573ac89 1574void nrmlwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1575 Double_t &pVx, Double_t &pVy, Double_t &pVz,
1576 Double_t *norml, const Int_t &oldReg,
1577 const Int_t &newReg, Int_t &flagErr)
efde9b4d 1578{
2573ac89 1579 printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz);
1580 printf(" oldReg=%i, newReg=%i\n", oldReg,newReg);
1581 Int_t curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
1582 Int_t curLttc = gGeoManager->GetCurrentNodeId()+1;
1583 printf(" curReg=%i, curLttc=%i in: %s\n", curreg, curLttc, gGeoManager->GetPath());
1584 Bool_t regsame = (curreg==oldReg)?kTRUE:kFALSE;
1585 gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
1586 gGeoManager->SetCurrentDirection(pVx,pVy,pVz);
1587 if (!regsame) {
1588 printf(" REGIONS DOEN NOT MATCH\n");
1589 gGeoManager->FindNode();
1590 curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
1591 curLttc = gGeoManager->GetCurrentNodeId()+1;
1592 printf(" re-initialized point: curReg=%i curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath());
1593 }
1594 Double_t *dnorm = gGeoManager->FindNormalFast();
1595 flagErr = 0;
1596 if (!dnorm) {
1597 printf(" ERROR: Cannot compute fast normal\n");
1598 flagErr = 1;
1599 norml[0] = -pVx;
1600 norml[1] = -pVy;
1601 norml[2] = -pVz;
1602 }
1603 norml[0] = -dnorm[0];
1604 norml[1] = -dnorm[1];
1605 norml[2] = -dnorm[2];
1606 printf(" normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]);
1607 curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
1608 curLttc = gGeoManager->GetCurrentNodeId()+1;
1609 printf(" final location: curReg=%i, curLttc=%i in %s\n", curreg,curLttc,gGeoManager->GetPath());
1610 printf("<= NRMLWR\n");
efde9b4d 1611}
1612
1613//_____________________________________________________________________________
1614void rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
1615 Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/)
1616{
2573ac89 1617 printf("=> Dummy RGRPWR\n");
efde9b4d 1618}
1619
1620//_____________________________________________________________________________
1621Int_t isvhwr(const Int_t &check, const Int_t & intHist)
1622{
1623// from FLUGG:
1624// Wrapper for saving current navigation history (fCheck=default)
1625// and returning its pointer. If fCheck=-1 copy of history pointed
1626// by intHist is made in NavHistWithCount object, and its pointer
1627// is returned. fCheck=1 and fCheck=2 cases are only in debugging
1628// version: an array is created by means of FGeometryInit functions
1629// (but could be a static int * ptrArray = new int[10000] with
1630// file scope as well) that stores a flag for deleted/undeleted
1631// histories and at the end of event is checked to verify that
1632// all saved history objects have been deleted.
1633
1634// For TGeo, just return the current node ID. No copy need to be made.
1635
1636 printf("=> Inside ISVHWR\n");
1637 if (check<0) return intHist;
2573ac89 1638 Int_t histInt = gGeoManager->GetCurrentNodeId()+1;
1639 printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath());
efde9b4d 1640 return histInt;
1641}
1642
1643
8495a208 1644
1645