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