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