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