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