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