Interface TGeo -> FLUKA (A. Gheata)
[u/mrichter/AliRoot.git] / TFluka / TFlukaMCGeometry.cxx
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 //____________________________________________________________________________ 
49 extern "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
91 TFluka *fluka = 0;
92 TFlukaMCGeometry *mcgeom = 0;
93
94 ClassImp(TFlukaMCGeometry)
95
96 TFlukaMCGeometry* TFlukaMCGeometry::fgInstance=0;
97
98 //_____________________________________________________________________________
99 TFlukaMCGeometry::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 //_____________________________________________________________________________
110 TFlukaMCGeometry::TFlukaMCGeometry()
111   : TVirtualMCGeometry()
112 {    
113   //
114   // Default constructor
115   //
116   fluka = (TFluka*)gMC;
117   mcgeom = this;
118 }
119
120 //_____________________________________________________________________________
121 TFlukaMCGeometry::~TFlukaMCGeometry() 
122 {
123   //
124   // Destructor
125   //
126   fgInstance=0;
127   if (gGeoManager) delete gGeoManager;
128 }
129
130 //
131 // private methods
132 //
133 //_____________________________________________________________________________
134 TFlukaMCGeometry::TFlukaMCGeometry(const TFlukaMCGeometry &)
135   : TVirtualMCGeometry()
136 {    
137   //
138   // Copy constructor
139   //
140 }
141
142 //_____________________________________________________________________________
143 Double_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 //_____________________________________________________________________________
163 void 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 //_____________________________________________________________________________
187 void 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 //_____________________________________________________________________________
211 void 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 //_____________________________________________________________________________
239 void 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 //_____________________________________________________________________________
267 void 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 //_____________________________________________________________________________
298 void 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 //_____________________________________________________________________________
331 Int_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 //_____________________________________________________________________________
343 void 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 //_____________________________________________________________________________
374 void 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 //_____________________________________________________________________________
403 void 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 //_____________________________________________________________________________
424 Int_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 //_____________________________________________________________________________
444 Int_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 //_____________________________________________________________________________
467 void  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 //_____________________________________________________________________________
491 void  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 //_____________________________________________________________________________
509 void  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 //_____________________________________________________________________________
532 void  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 //_____________________________________________________________________________
555 void  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 //_____________________________________________________________________________
576 void  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 //_____________________________________________________________________________
609 void  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 //_____________________________________________________________________________
624 void  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 //_____________________________________________________________________________
647 Int_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 //_____________________________________________________________________________
663 const 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 //_____________________________________________________________________________
678 Int_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 //_____________________________________________________________________________
688 Int_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 //_____________________________________________________________________________
705 Int_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 //_____________________________________________________________________________
717 Int_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 //_____________________________________________________________________________
732 const 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 //_____________________________________________________________________________
742 const 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 //_____________________________________________________________________________
758 void 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 //_____________________________________________________________________________
773 void 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 //_____________________________________________________________________________
795 void 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 //_____________________________________________________________________________
802 void 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 //_____________________________________________________________________________
825 void 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 //_____________________________________________________________________________
832 void 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 //_____________________________________________________________________________
1088 void 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 //_____________________________________________________________________________
1100 Int_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 //_____________________________________________________________________________
1108 void 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 //______________________________________________________________________________
1126 void 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 //_____________________________________________________________________________
1143 void  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 //_____________________________________________________________________________
1154 Int_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 //_____________________________________________________________________________
1167 Int_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 //_____________________________________________________________________________
1188 void 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