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