Some corrections. (A. Gheata)
[u/mrichter/AliRoot.git] / TFluka / TFlukaMCGeometry.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 // $Id$
17
18 // Class TFlukaMCGeometry
19 // --------------------
20 // Implementation of the TVirtualMCGeometry interface
21 // for defining and using TGeo geometry with FLUKA.
22 // This allows the FLUKA MonteCarlo to run with the TGeo 
23 // geometrical modeller
24 // Author: Andrei Gheata 10/07/2003
25
26 #include "Riostream.h"
27 #include "TCallf77.h"
28 #include "TFluka.h"
29 #include "TFlukaMCGeometry.h"
30 #include "TFlukaConfigOption.h"
31 #include "TGeoManager.h" 
32 #include "TGeoVolume.h" 
33 #include "TObjString.h"
34 #include "Fepisor.h"
35
36 #ifndef WIN32 
37 # define idnrwr idnrwr_
38 # define g1wr   g1wr_
39 # define g1rtwr g1rtwr_
40 # define conhwr conhwr_
41 # define inihwr inihwr_
42 # define jomiwr jomiwr_
43 # define lkdbwr lkdbwr_
44 # define lkfxwr lkfxwr_
45 # define lkmgwr lkmgwr_
46 # define lkwr lkwr_
47 # define magfld magfld_
48 # define nrmlwr nrmlwr_
49 # define rgrpwr rgrpwr_
50 # define isvhwr isvhwr_
51
52 #else
53
54 # define idnrwr IDNRWR
55 # define g1wr   G1WR
56 # define g1rtwr G1RTWR
57 # define conhwr CONHWR
58 # define inihwr INIHWR
59 # define jomiwr JOMIWR
60 # define lkdbwr LKDBWR
61 # define lkfxwr LKFXWR
62 # define lkmgwr LKMGWR
63 # define lkwr   LKWR
64 # define magfld MAGFLD
65 # define nrmlwr NRMLWR
66 # define rgrpwr RGRPWR
67 # define isvhwr ISVHWR
68
69 #endif
70
71 //____________________________________________________________________________ 
72 extern "C" 
73 {
74    //
75    // Prototypes for FLUKA navigation methods
76    //
77    Int_t type_of_call idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/);
78    void  type_of_call   g1wr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, 
79                              Double_t * /*pV*/,  Int_t & /*oldReg*/ , const Int_t & /*oldLttc*/, Double_t & /*propStep*/,
80                              Int_t & /*nascFlag*/, Double_t & /*retStep*/, Int_t & /*newReg*/,
81                                   Double_t & /*saf*/, Int_t & /*newLttc*/, Int_t & /*LttcFlag*/,
82                              Double_t *s /*Lt*/, Int_t * /*jrLt*/);
83    
84    void  type_of_call g1rtwr();
85    void  type_of_call conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/); 
86    void  type_of_call inihwr(Int_t & /*intHist*/);
87    void  type_of_call jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
88                              Int_t & /*flukaReg*/);
89    void  type_of_call lkdbwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
90                              Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
91                              Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/);
92    void  type_of_call lkfxwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
93                              Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
94                              Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/);
95    void  type_of_call lkmgwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
96                              Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
97                                        Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
98    void  type_of_call   lkwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
99                              Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
100                                   Int_t & /*newReg*/, Int_t & /*flagErr*/, Int_t & /*newLttc*/);
101 //   void  type_of_call magfld(const Double_t & /*pX*/, const Double_t & /*pY*/, const Double_t & /*pZ*/,
102 //                             Double_t & /*cosBx*/, Double_t & /*cosBy*/, Double_t & /*cosBz*/, 
103 //                             Double_t & /*Bmag*/, Int_t & /*reg*/, Int_t & /*idiscflag*/);        
104    void  type_of_call nrmlwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
105                              Double_t & /*pVx*/, Double_t & /*pVy*/, Double_t & /*pVz*/,
106                                   Double_t * /*norml*/, const Int_t & /*oldReg*/, 
107                                   const Int_t & /*newReg*/, Int_t & /*flagErr*/);
108    void  type_of_call rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
109                              Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/);
110    Int_t type_of_call isvhwr(const Int_t & /*fCheck*/, const Int_t & /*intHist*/);
111 };
112    
113 // TFluka global pointer
114 TFluka *gFluka = 0;
115 TFlukaMCGeometry *gMCGeom = 0;
116 Int_t gNstep = 0;
117
118 ClassImp(TFlukaMCGeometry)
119
120 TFlukaMCGeometry* TFlukaMCGeometry::fgInstance=0;
121
122 //_____________________________________________________________________________
123 TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title) 
124                  : TNamed(name, title)
125 {
126   //
127   // Standard constructor
128   //
129   fDebug        = kFALSE;
130   fLastMaterial = 0;
131   fNextRegion   = 0;
132   fNextLattice  = 0;
133   fRegionList   = 0;
134   gFluka = (TFluka*)gMC;
135   gMCGeom = this;
136   gNstep = 0;
137   fMatList = new TObjArray(256);
138   fMatNames = new TObjArray(256);
139 }
140
141 //_____________________________________________________________________________
142 TFlukaMCGeometry::TFlukaMCGeometry()
143 {    
144   //
145   // Default constructor
146   //
147   fDebug        = kFALSE;
148   fLastMaterial = 0;
149   fNextRegion   = 0;
150   fNextLattice  = 0;
151   fRegionList   = 0;
152   gFluka = (TFluka*)gMC;
153   gMCGeom = this;
154   gNstep = 0;
155   fMatList = 0;
156   fMatNames = 0;
157 }
158
159 //_____________________________________________________________________________
160 TFlukaMCGeometry::~TFlukaMCGeometry() 
161 {
162   //
163   // Destructor
164   //
165   fgInstance=0;
166   if (fRegionList) delete [] fRegionList;
167   if (fMatList) delete fMatList;
168   if (fMatNames) {fMatNames->Delete(); delete fMatNames;}
169   if (gGeoManager) delete gGeoManager;
170 }
171
172 //
173 // private methods
174 //
175 //_____________________________________________________________________________
176 TFlukaMCGeometry::TFlukaMCGeometry(const TFlukaMCGeometry &)
177   : TNamed()
178 {    
179   //
180   // Copy constructor
181   //
182 }
183
184 //_____________________________________________________________________________
185 Double_t* TFlukaMCGeometry::CreateDoubleArray(Float_t* array, Int_t size) const
186 {
187 // Converts Float_t* array to Double_t*,
188 // !! The new array has to be deleted by user.
189 // ---
190
191   Double_t* doubleArray;
192   if (size>0) {
193     doubleArray = new Double_t[size]; 
194     for (Int_t i=0; i<size; i++) doubleArray[i] = array[i];
195   }
196   else {
197     //doubleArray = 0; 
198     doubleArray = new Double_t[1]; 
199   }  
200   return doubleArray;
201 }
202 //
203 // public methods
204
205
206 //_____________________________________________________________________________
207 Int_t TFlukaMCGeometry::GetMedium() const
208 {
209 // Get current medium number
210    Int_t imed = 0;
211    TGeoNode *node = gGeoManager->GetCurrentNode();
212    if (!node) imed = gGeoManager->GetTopNode()->GetVolume()->GetMedium()->GetId();
213    else       imed = node->GetVolume()->GetMedium()->GetId();
214    if (fDebug) printf("GetMedium=%i\n", imed);
215    return imed;
216 }
217
218 //_____________________________________________________________________________
219 Int_t TFlukaMCGeometry::GetFlukaMaterial(Int_t imed) const
220 {
221 // Returns FLUKA material index for medium IMED
222    TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1);
223    if (!med) {
224       Error("GetFlukaMaterial", "MEDIUM %i nor found", imed);
225       return -1;
226    }
227    Int_t imatfl = med->GetMaterial()->GetIndex();
228    return imatfl;   
229 }
230
231 //_____________________________________________________________________________
232 Int_t *TFlukaMCGeometry::GetRegionList(Int_t imed, Int_t &nreg)
233 {
234 // Get an ordered list of regions matching a given medium number
235    nreg = 0;
236    if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
237    TIter next(gGeoManager->GetListOfUVolumes());
238    TGeoVolume *vol;
239    Int_t imedium, ireg;
240    while ((vol = (TGeoVolume*)next())) {
241       imedium = vol->GetMedium()->GetId();
242       if (imedium == imed) {
243          ireg = vol->GetNumber();
244          fRegionList[nreg++] = ireg;
245       }
246    }
247    return fRegionList;
248 }         
249
250 //_____________________________________________________________________________
251 Int_t *TFlukaMCGeometry::GetMaterialList(Int_t imat, Int_t &nreg)
252 {
253 // Get an ordered list of regions matching a given medium number
254    nreg = 0;
255    if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
256    TIter next(gGeoManager->GetListOfUVolumes());
257    TGeoVolume *vol;
258    Int_t imaterial, ireg;
259    while ((vol = (TGeoVolume*)next())) {
260       imaterial = vol->GetMedium()->GetMaterial()->GetIndex();
261       if (imaterial == imat) {
262          ireg = vol->GetNumber();
263          fRegionList[nreg++] = ireg;
264       }
265    }
266    return fRegionList;
267 }         
268  
269 //_____________________________________________________________________________
270 Int_t TFlukaMCGeometry::NofVolumes() const 
271 {
272   //
273   // Return total number of volumes in the geometry
274   //
275
276   return gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1;
277 }
278   
279 //_____________________________________________________________________________
280 TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z)
281 {
282 // Try to replace a wrongly-defined material
283    static Double_t kz[23] = {7.3, 17.8184, 7.2167, 10.856, 8.875, 8.9, 7.177,
284       25.72, 6.2363, 7.1315, 47.7056, 10.6467, 7.8598, 2.10853, 10.6001, 9.1193, 
285       15.3383, 4.55,   9.6502, 6.4561, 21.7963, 29.8246, 15.4021};
286
287    Int_t ind;
288    Double_t dz;
289    for (ind=0; ind<23; ind++) {
290       dz = TMath::Abs(z-kz[ind]);
291       if (dz<1E-4) break;
292    }
293    if (ind>22) {
294       printf("Cannot patch material with Z=%g\n", z);
295       return 0;
296    }
297    TGeoMixture *mix = 0;
298    TGeoElement *element;
299    TGeoElementTable *table = TGeoElementTable::Instance();
300    switch (ind) {
301       case 0: // AIR
302          mix = new TGeoMixture("AIR", 4, 0.001205);
303          element = table->GetElement(6); // C
304          mix->DefineElement(0, element, 0.000124);
305          element = table->GetElement(7); // N
306          mix->DefineElement(1, element, 0.755267);
307          element = table->GetElement(8); // O
308          mix->DefineElement(2, element, 0.231781);
309          element = table->GetElement(18); // AR
310          mix->DefineElement(3, element, 0.012827);
311          break;
312       case 1: //SDD SI CHIP
313          mix = new TGeoMixture("SDD_SI", 6, 2.4485);
314          element = table->GetElement(1);
315          mix->DefineElement(0, element, 0.004367771);
316          element = table->GetElement(6);
317          mix->DefineElement(1, element, 0.039730642);
318          element = table->GetElement(7);
319          mix->DefineElement(2, element, 0.001396798);
320          element = table->GetElement(8);
321          mix->DefineElement(3, element, 0.01169634);
322          element = table->GetElement(14);
323          mix->DefineElement(4, element, 0.844665);
324          element = table->GetElement(47);
325          mix->DefineElement(5, element, 0.09814344903);
326          break;
327       case 2:  // WATER
328          mix = new TGeoMixture("WATER", 2, 1.0);
329          element = table->GetElement(1);
330          mix->DefineElement(0, element, 0.111898344);
331          element = table->GetElement(8);
332          mix->DefineElement(1, element, 0.888101656);
333          break;
334       case 3: // CERAMICS
335          mix = new TGeoMixture("CERAMICS", 5, 3.6);
336          element = table->GetElement(8);
337          mix->DefineElement(0, element, 0.59956);
338          element = table->GetElement(13);
339          mix->DefineElement(1, element, 0.3776);
340          element = table->GetElement(14);
341          mix->DefineElement(2, element, 0.00933);
342          element = table->GetElement(24);
343          mix->DefineElement(3, element, 0.002);
344          element = table->GetElement(25);
345          mix->DefineElement(4, element, 0.0115);
346          break;
347       case 4: // EPOXY
348          mix = new TGeoMixture("G10FR4", 4, 1.8);
349          element = table->GetElement(1);
350          mix->DefineElement(0, element, 0.19);
351          element = table->GetElement(6);
352          mix->DefineElement(1, element, 0.18);
353          element = table->GetElement(8);
354          mix->DefineElement(2, element, 0.35);
355          element = table->GetElement(14);
356          mix->DefineElement(3, element, 0.28);
357          break;
358       case 5: // EPOXY
359          mix = new TGeoMixture("G10FR4", 4, 1.8);
360          element = table->GetElement(1);
361          mix->DefineElement(0, element, 0.19);
362          element = table->GetElement(6);
363          mix->DefineElement(1, element, 0.18);
364          element = table->GetElement(8);
365          mix->DefineElement(2, element, 0.35);
366          element = table->GetElement(14);
367          mix->DefineElement(3, element, 0.28);
368          break;
369       case 6: // KAPTON
370          mix = new TGeoMixture("KAPTON", 4, 1.3);
371          element = table->GetElement(1);
372          mix->DefineElement(0, element, 0.026363415);
373          element = table->GetElement(6);
374          mix->DefineElement(1, element, 0.6911272);
375          element = table->GetElement(7);
376          mix->DefineElement(2, element, 0.073271325);
377          element = table->GetElement(8);
378          mix->DefineElement(3, element, 0.209238060);
379          break;
380       case 7: // INOX
381          mix = new TGeoMixture("INOX", 9, 7.9);
382          element = table->GetElement(6);
383          mix->DefineElement(0, element, 0.0003);
384          element = table->GetElement(14);
385          mix->DefineElement(1, element, 0.01);          
386          element = table->GetElement(15);
387          mix->DefineElement(2, element, 0.00045);
388          element = table->GetElement(16);
389          mix->DefineElement(3, element, 0.0003);
390          element = table->GetElement(24);
391          mix->DefineElement(4, element, 0.17);
392          element = table->GetElement(25);
393          mix->DefineElement(5, element, 0.02);
394          element = table->GetElement(26);
395          mix->DefineElement(6, element, 0.654);
396          element = table->GetElement(28);
397          mix->DefineElement(7, element, 0.12);
398          element = table->GetElement(42);
399          mix->DefineElement(8, element, 0.025);
400          break;
401       case 8: // ROHACELL
402          mix = new TGeoMixture("ROHACELL", 4, 0.05);
403          element = table->GetElement(1);
404          mix->DefineElement(0, element, 0.07836617);
405          element = table->GetElement(6);
406          mix->DefineElement(1, element, 0.64648941);
407          element = table->GetElement(7);
408          mix->DefineElement(2, element, 0.08376983);
409          element = table->GetElement(8);
410          mix->DefineElement(3, element, 0.19137459);
411          break;
412       case 9: // SDD-C-AL
413          mix = new TGeoMixture("SDD-C-AL", 5, 1.9837);
414          element = table->GetElement(1);
415          mix->DefineElement(0, element, 0.022632);
416          element = table->GetElement(6);
417          mix->DefineElement(1, element, 0.8176579);
418          element = table->GetElement(7);
419          mix->DefineElement(2, element, 0.0093488);
420          element = table->GetElement(8);
421          mix->DefineElement(3, element, 0.0503618);
422          element = table->GetElement(13);
423          mix->DefineElement(4, element, 0.1);
424          break;
425       case 10: // X7R-CAP
426          mix = new TGeoMixture("X7R-CAP", 7, 6.72);
427          element = table->GetElement(8);
428          mix->DefineElement(0, element, 0.085975822);
429          element = table->GetElement(22);
430          mix->DefineElement(1, element, 0.084755042);
431          element = table->GetElement(28);
432          mix->DefineElement(2, element, 0.038244751);
433          element = table->GetElement(29);
434          mix->DefineElement(3, element, 0.009471271);
435          element = table->GetElement(50);
436          mix->DefineElement(4, element, 0.321736471);
437          element = table->GetElement(56);
438          mix->DefineElement(5, element, 0.251639432);
439          element = table->GetElement(82);
440          mix->DefineElement(6, element, 0.2081768);
441          break;
442       case 11: // SDD ruby sph. Al2O3
443          mix = new TGeoMixture("AL2O3", 2, 3.97);
444          element = table->GetElement(8);
445          mix->DefineElement(0, element, 0.5293);
446          element = table->GetElement(13);
447          mix->DefineElement(1, element, 0.4707);
448          break;
449       case 12: // SDD HV microcable
450          mix = new TGeoMixture("HV-CABLE", 5, 1.6087);
451          element = table->GetElement(1);
452          mix->DefineElement(0, element, 0.01983871336);
453          element = table->GetElement(6);
454          mix->DefineElement(1, element, 0.520088819984);
455          element = table->GetElement(7);
456          mix->DefineElement(2, element, 0.0551367996);
457          element = table->GetElement(8);
458          mix->DefineElement(3, element, 0.157399667056);
459          element = table->GetElement(13);
460          mix->DefineElement(4, element, 0.247536);
461          break;
462       case 13: //SDD LV+signal cable
463          mix = new TGeoMixture("LV-CABLE", 5, 2.1035);
464          element = table->GetElement(1);
465          mix->DefineElement(0, element, 0.0082859922);
466          element = table->GetElement(6);
467          mix->DefineElement(1, element, 0.21722436468);
468          element = table->GetElement(7);
469          mix->DefineElement(2, element, 0.023028867);
470          element = table->GetElement(8);
471          mix->DefineElement(3, element, 0.06574077612);
472          element = table->GetElement(13);
473          mix->DefineElement(4, element, 0.68572);
474          break;
475       case 14: //SDD hybrid microcab
476          mix = new TGeoMixture("HYB-CAB", 5, 2.0502);
477          element = table->GetElement(1);
478          mix->DefineElement(0, element, 0.00926228815);
479          element = table->GetElement(6);
480          mix->DefineElement(1, element, 0.24281879711);
481          element = table->GetElement(7);
482          mix->DefineElement(2, element, 0.02574224025);
483          element = table->GetElement(8);
484          mix->DefineElement(3, element, 0.07348667449);
485          element = table->GetElement(13);
486          mix->DefineElement(4, element, 0.64869);
487          break;
488       case 15: //SDD anode microcab
489          mix = new TGeoMixture("ANOD-CAB", 5, 1.7854);
490          element = table->GetElement(1);
491          mix->DefineElement(0, element, 0.0128595919215);
492          element = table->GetElement(6);
493          mix->DefineElement(1, element, 0.392653705471);
494          element = table->GetElement(7);
495          mix->DefineElement(2, element, 0.041626868025);
496          element = table->GetElement(8);
497          mix->DefineElement(3, element, 0.118832707289);
498          element = table->GetElement(13);
499          mix->DefineElement(4, element, 0.431909);
500          break;
501       case 16: // inox/alum
502          mix = new TGeoMixture("INOX-AL", 5, 3.0705);
503          element = table->GetElement(13);
504          mix->DefineElement(0, element, 0.816164);
505          element = table->GetElement(14);
506          mix->DefineElement(1, element, 0.000919182);
507          element = table->GetElement(24);
508          mix->DefineElement(2, element, 0.0330906);
509          element = table->GetElement(26);
510          mix->DefineElement(3, element, 0.131443);
511          element = table->GetElement(28);
512          mix->DefineElement(4, element, 0.0183836);
513       case 17: // MYLAR
514          mix = new TGeoMixture("MYLAR", 3, 1.39);
515          element = table->GetElement(1);
516          mix->DefineElement(0, element, 0.0416667);
517          element = table->GetElement(6);
518          mix->DefineElement(1, element, 0.625);
519          element = table->GetElement(8);
520          mix->DefineElement(2, element, 0.333333);
521          break;
522       case 18: // SPDBUS(AL+KPT+EPOX)   - unknown composition
523          mix = new TGeoMixture("SPDBUS", 1, 1.906);
524          element = table->GetElement(9);
525          mix->DefineElement(0, element, 1.);
526          z = element->Z();
527          break;
528       case 19: // SDD/SSD rings   - unknown composition
529          mix = new TGeoMixture("SDDRINGS", 1, 1.8097);
530          element = table->GetElement(6);
531          mix->DefineElement(0, element, 1.);
532          z = element->Z();
533          break;
534       case 20: // SPD end ladder   - unknown composition
535          mix = new TGeoMixture("SPDEL", 1, 3.6374);
536          element = table->GetElement(22);
537          mix->DefineElement(0, element, 1.);
538          z = element->Z();
539          break;
540       case 21: // SDD end ladder   - unknown composition
541          mix = new TGeoMixture("SDDEL", 1, 0.3824);
542          element = table->GetElement(30);
543          mix->DefineElement(0, element, 1.);
544          z = element->Z();
545          break;
546       case 22: // SSD end ladder   - unknown composition
547          mix = new TGeoMixture("SSDEL", 1, 0.68);
548          element = table->GetElement(16);
549          mix->DefineElement(0, element, 1.);
550          z = element->Z();
551          break;
552    }
553    mix->SetZ(z);      
554    printf("Patched with mixture %s\n", mix->GetName());
555    return mix;
556 }   
557
558 //_____________________________________________________________________________
559 void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname)
560 {
561   // ==== from FLUGG ====
562   // NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case,
563   // according to the fluka standard. In addition,. they must be equal to the
564   // names of the fluka materials - see fluka manual - in order that the 
565   // program load the right cross sections, and equal to the names included in
566   // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his
567   // own .pemf, in order to get the right cross sections loaded in memory.
568
569
570    TString sname;
571    gGeoManager->Export("flgeom.root");
572    if (fname) sname = fname;
573    else       sname = "flukaMat.inp";
574    ofstream out;
575    out.open(sname.Data(), ios::out);
576    if (!out.good()) {
577       Fatal("CreateFlukaMatFile", "could not open file %s for writing", sname.Data());
578       return;
579    }
580    PrintHeader(out, "MATERIALS AND COMPOUNDS");
581    PrintHeader(out, "MATERIALS");   
582    Int_t i,j,idmat;
583    Int_t counttothree, nelem;
584    Double_t a,z,rho, w;
585    TGeoElementTable *table = TGeoElementTable::Instance();
586    TGeoElement *element;
587    element = table->GetElement(13);
588    element->SetTitle("ALUMINUM"); // this is how FLUKA likes it ...
589    element = table->GetElement(15);
590    element->SetTitle("PHOSPHO");  // same story ...
591 //   element = table->GetElement(10);
592 //   element->SetTitle("ARGON");  // NEON not in neutron xsec table
593    Int_t nelements = table->GetNelements();
594    TList *matlist = gGeoManager->GetListOfMaterials();
595 //   TList *medlist = gGeoManager->GetListOfMedia();
596 //   Int_t nmed = medlist->GetSize();
597    TIter next(matlist);
598    Int_t nmater = matlist->GetSize();
599    Int_t nfmater = 0;
600    TGeoMaterial *mat;
601    TGeoMixture *mix = 0;
602    TString matname;
603    TObjString *objstr;
604    // Create all needed elements
605    for (Int_t i=1; i<nelements; i++) {
606       element = table->GetElement(i);
607       // skip elements which are not defined
608       if (!element->IsUsed() && !element->IsDefined()) continue;
609       matname = element->GetTitle();
610       ToFlukaString(matname);
611       rho = 0.999;
612
613       mat = new TGeoMaterial(matname, element->A(), element->Z(), rho);
614       mat->SetIndex(nfmater+3);
615       mat->SetUsed(kTRUE);
616       fMatList->Add(mat);
617       objstr = new TObjString(matname.Data());
618       fMatNames->Add(objstr);
619       nfmater++;
620    }
621    
622    fIndmat = nfmater;
623 //   TGeoMedium *med;
624    // Adjust material names and add them to FLUKA list
625    for (i=0; i<nmater; i++) {
626       mat = (TGeoMaterial*)matlist->At(i);
627       if (!mat->IsUsed()) continue;      
628       z = mat->GetZ();
629       a = mat->GetA();
630       rho = mat->GetDensity();
631       if (mat->GetZ()<0.001) {
632          mat->SetIndex(2); // vacuum, built-in inside FLUKA
633          continue;
634       } 
635       matname = mat->GetName();
636       FlukaMatName(matname);
637 /*
638       // material with one element: create it as mixture since it can be duplicated    
639       if (!mat->IsMixture()) {
640          // normal materials
641          mix = new TGeoMixture(matname.Data(), 1, rho);
642          mix->DefineElement(0, mat->GetElement(), 1.);
643          mat->SetIndex(nfmater+3);
644          for (j=0; j<nmed; j++) {
645             med = (TGeoMedium*)medlist->At(j);
646             if (med->GetMaterial() == mat) {
647                med->SetMaterial(mix);
648                if (mat->GetCerenkovProperties()) {
649                   mix->SetCerenkovProperties(mat->GetCerenkovProperties());
650                   mat->SetCerenkovProperties(0);
651                }
652                break;
653             }
654          }                              
655          mat = (TGeoMaterial*)mix;
656       }
657 */    
658       mat->SetIndex(nfmater+3);
659       objstr = new TObjString(matname.Data());
660       fMatList->Add(mat);
661       fMatNames->Add(objstr);
662       nfmater++;
663    }   
664
665    // Dump all elements with MATERIAL cards         
666    for (i=0; i<nfmater; i++) {
667       mat = (TGeoMaterial*)fMatList->At(i);
668 //      mat->SetUsed(kFALSE);
669       mix = 0;
670       out << setw(10) << "MATERIAL  ";
671       out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
672       objstr = (TObjString*)fMatNames->At(i);
673       matname = objstr->GetString();
674       z = mat->GetZ();
675       a = mat->GetA();
676       rho = mat->GetDensity();
677       if (mat->IsMixture()) {
678          out << setw(10) << " ";
679          out << setw(10) << " ";
680          mix = (TGeoMixture*)mat;
681       } else {   
682          out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << z;
683          out << setw(10) << setprecision(3) << a;
684       }
685       out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
686       out << setw(10) << setiosflags(ios::scientific) << setprecision(3) << rho;
687       out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
688       out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex());   
689       out << setw(10) << " ";
690       out << setw(10) << " ";
691       out << setw(8) << matname.Data() << endl;
692       if (!mix) {
693          // add LOW-MAT card for NEON to associate with ARGON neutron xsec
694          if (z==10) {
695             out << setw(10) << "LOW-MAT   ";
696             out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
697             out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex());
698             out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << 18.;
699             out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << -2.;
700             out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << 293.;
701             out << setw(10) << " ";
702             out << setw(10) << " ";
703 //            out << setw(8) << matname.Data() << endl;
704             out << setw(8) << " " << endl;
705          } 
706          else { 
707             element = table->GetElement((int)z);
708             TString elename = element->GetTitle();
709             ToFlukaString(elename);
710             if( matname.CompareTo( elename ) != 0 ) {
711                out << setw(10) << "LOW-MAT   ";
712                out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
713                out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex());
714                out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << z;
715                out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << " ";
716                out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << " ";
717                out << setw(10) << " ";
718                out << setw(10) << " ";
719                // missing material at Low Energy Cross Section Table
720                if( (int)z==10 || (int)z==21 || (int)z==34 || (int)z==37 || (int)z==39 || (int)z==44 ||
721                    (int)z==45 || (int)z==46 || (int)z==52 || (int)z==57 || (int)z==59 || (int)z==60 ||
722                    (int)z==61 || (int)z==65 || (int)z==66 || (int)z==67 || (int)z==68 || (int)z==69 ||
723                    (int)z==70 || (int)z==71 || (int)z==72 || (int)z==76 || (int)z==77 || (int)z==78 ||
724                    (int)z==81 || (int)z==84 || (int)z==85 || (int)z==86 || (int)z==87 || (int)z==88 ||
725                    (int)z==89 || (int)z==91 )
726                   out << setw(8) << "UNKNOWN " << endl;
727                else
728                   out << setw(8) << elename.Data() << endl;
729    //               out << setw(8) << " " << endl;
730             }
731          }
732          continue;
733       }   
734       counttothree = 0;
735       out << setw(10) << "COMPOUND  ";
736       nelem = mix->GetNelements();
737       objstr = (TObjString*)fMatNames->At(i);
738       matname = objstr->GetString();
739       for (j=0; j<nelem; j++) {
740          w = (mix->GetWmixt())[j];
741          if (w<0.00001) w=0.00001;
742          z = (mix->GetZmixt())[j];       
743          a = (mix->GetAmixt())[j];
744          idmat = GetElementIndex(Int_t(z));
745          if (!idmat) Error("CreateFlukaMatFile", "element with Z=%f not found", z);
746          out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
747          out << setw(10) << setiosflags(ios::fixed) << setprecision(6) << -w;   
748          out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
749          out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(idmat);
750          counttothree++;
751          if (counttothree == 3) {
752             out << matname.Data();
753             out << endl;
754             if ( (j+1) != nelem) out << setw(10) << "COMPOUND  ";
755             counttothree = 0;
756          } 
757       }               
758       if (nelem%3) {
759          for (j=0; j<(3-(nelem%3)); j++)
760             out << setw(10) << " " << setw(10) << " ";
761          out << matname.Data();
762          out << endl;
763       } 
764    }     
765    Int_t nvols = gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1;
766    TGeoVolume *vol;
767    // Now print the material assignments
768    Double_t flagfield = 0.;
769    printf("#############################################################\n");
770    if (gFluka->IsFieldEnabled()) {
771       flagfield = 1.;
772       printf("Magnetic field enabled\n");
773    } else printf("Magnetic field disabled\n");   
774    printf("#############################################################\n");
775    
776    PrintHeader(out, "TGEO MATERIAL ASSIGNMENTS");   
777    for (i=1; i<=nvols; i++) {
778        
779       vol = gGeoManager->GetVolume(i);
780       mat = vol->GetMedium()->GetMaterial();
781       idmat = mat->GetIndex();
782       for (Int_t j=0; j<nfmater; j++) {
783          mat = (TGeoMaterial*)fMatList->At(j);
784          if (mat->GetIndex() == idmat) mat->SetUsed(kTRUE);
785       }   
786
787       Float_t hasfield  = (vol->GetMedium()->GetParam(1) > 0) ? flagfield : 0.;
788       out << "* Assigning material:   " << vol->GetMedium()->GetMaterial()->GetName() << "   to Volume: " << vol->GetName();
789       out << endl;
790
791       out << setw(10) << "ASSIGNMAT ";
792       out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
793       out << setw(10) << setiosflags(ios::fixed) << Double_t(idmat);
794       out << setw(10) << setiosflags(ios::fixed) << Double_t(i);
795       out << setw(10) << "0.0";
796       out << setw(10) << "0.0";
797       out << setw(10) << setiosflags(ios::fixed) <<  hasfield;
798       out << setw(10) << "0.0";
799       out << endl;
800    }
801    out.close();
802    fLastMaterial = nfmater+2;
803    
804    if (!gFluka->IsGeneratePemf()) {
805        if (gSystem->AccessPathName("FlukaVmc.pemf")) Fatal("CreateFlukaMatFile", "No pemf file in working directory");
806        return;
807    }
808 }
809
810 void TFlukaMCGeometry::CreatePemfFile()
811 {
812     //
813     // Steering routine to write and process peg files producing the pemf input 
814     //
815     char number[20];
816     Int_t countMatOK     = 0;
817     Int_t countElemError = 0;
818     Int_t countNoStern   = 0;
819     Int_t countMixError  = 0;
820     Int_t countGas       = 0;
821     Int_t countPemfError = 0;
822     Int_t i;
823     TGeoMaterial* mat = 0x0;
824     TString sname;
825     
826     for (i = fIndmat; i < fLastMaterial - 2; i++) {
827         printf("Write Peg Files %d\n", i);
828         
829         mat = (TGeoMaterial*)fMatList->At(i);
830         if (!mat->IsUsed()) continue;
831         sname = "mat";
832         sprintf(number, "%d", i);
833         sname.Append(number);
834         cout << endl;
835         cout << endl;
836         cout << "******************************************************************************" << endl;
837         cout << "******************************************************************************" << endl;
838         cout << endl;
839         WritePegFile(i, &countNoStern, &countElemError, &countMixError, &countGas);
840         sname.Prepend("$FLUPRO/pemf/rpemf peg/");
841         gSystem->Exec(sname.Data());
842         
843         // check if the pemf file was created
844         TString sname = Form("peg/mat%d.pemf", i);
845         ifstream in( sname.Data() );
846         if ( in ) {
847             countMatOK++;
848             in.close();
849         } else {
850             cout << "ERROR Fail to create the pemf file " << sname << endl;
851             countPemfError++; 
852         }
853     }
854     cout << "Materials (pemf created)   " << countMatOK         << endl;
855     cout << "Not Sternheimer par. found  " << countNoStern   << endl;
856     cout << "Elements with error definitions (Z not integer)  " << countElemError      << endl;
857     cout << "Mixtures with error definitions (Z not integer) " << countMixError  << endl;
858     cout << "Posible Gas (rho < 0.01) " << countGas           << endl;
859     // cout << "Posible Gas (without pressure information) " << countGasError           << endl;
860     cout << "Pemf files Error    " << countPemfError     << endl;
861     cout << endl << endl;
862     
863     sname = "cat peg/*.pemf > peg/FlukaVmc.pemf";         
864     gSystem->Exec(sname.Data());
865     sname = "mv peg/FlukaVmc.pemf FlukaVmc.pemf";
866     gSystem->Exec(sname.Data());
867 }
868
869 //_____________________________________________________________________________
870 void TFlukaMCGeometry::WritePegFile(Int_t imat, Int_t *NoStern, Int_t *ElemError,
871                        Int_t *MixError, Int_t *countGas) const
872 {
873    // Write the .peg file for one material
874    
875    TGeoMaterial *mat = (TGeoMaterial*)fMatList->At(imat);
876    TString name = ((TObjString*)fMatNames->At(imat))->GetString();
877    TString line;
878    char number[20];
879    TGeoElement *elem = mat->GetElement();
880    name = name.Strip();
881    TString sname = "mat";
882    sprintf(number, "%d", imat);
883    sname.Append(number);
884    sname.Append(".peg");
885    sname.Prepend("peg/");
886    ofstream out;
887    out.open(sname.Data(), ios::out);
888    if (!out.good()) return;
889    Double_t dens = mat->GetDensity();
890    TGeoMixture *mix = 0;
891    Int_t nel = 1;
892    Int_t i;
893    if (mat->IsMixture()) {
894       mix = (TGeoMixture*)mat;
895       nel = mix->GetNelements();
896    } 
897      
898    if (nel==1) {
899       cout  << "( Element ) " << name << "  Z=" << mat->GetZ() << " Rho " << mat->GetDensity() << endl;
900
901       Double_t zel = mat->GetZ();
902       if( (zel-Int_t(zel))>0.001 || zel < 1 ) {
903          cout << " ERROR: A Element with not integer Z=" << zel << endl;
904          cout << endl;
905          (*ElemError)++;
906          return;
907       }
908       
909       out << "ELEM" << endl;
910       out << " &INP IRAYL=1, RHO=" << dens << ", " << endl;
911       
912       // check for the Sternheimer parameters
913       Double_t *issb_parm = GetISSB( mat->GetDensity(), 1, &zel, 0 );
914       if( issb_parm[0] > 0 && issb_parm[1] > 0 ) {
915          cout << "Sternheimer parameters found" << endl;
916          out << ", ISSB=1, IEV=" << issb_parm[0] << ", CBAR=" << issb_parm[1]
917              << ", X0=" << issb_parm[2] << "," << endl;
918          out << "X1=" <<issb_parm[3] <<", AFACT="<<issb_parm[4] <<", SK="
919              << issb_parm[5] << ", DELTA0=" << issb_parm[6];
920       } 
921       else {
922          cout << "WARNING: Strange element, Sternheimer parameters  not found" << endl;
923         (*NoStern)++;
924       }
925
926       if (dens<0.01) {
927         (*countGas)++;
928         out << " GASP=1." << endl;
929       }
930       
931       out << " &END" <<  endl;
932       out << name.Data() << endl;
933       out << elem->GetName() << endl;
934       
935    } 
936    else {
937    
938       cout  << "( Mixture )  " << name << "  Rho " << dens << " nElem " << nel << endl;
939     
940       Double_t *zt = new Double_t[nel];
941       Double_t *wt = new Double_t[nel];
942       for (int j=0; j<nel; j++) {
943          zt[j] = (mix->GetZmixt())[j];
944          wt[j] = (mix->GetWmixt())[j];
945          if( (zt[j]-Int_t(zt[j])) > 0.001 || zt[j] < 1 ) {
946             cout << "ERROR Mixture " << name << " with an element with not integer Z=" << zt[j] << endl;
947             cout << endl;
948             (*MixError)++;
949             // just continue since the mixtures are not patch, 
950             // but the final release should include the return   
951             //  return;         
952          }
953       }
954       Double_t *issb_parm = GetISSB( mat->GetDensity(), nel, zt, wt );
955       out << "MIXT" << endl;
956       out << " &INP IRAYL=1, NE=" << nel << ", RHOZ=" << wt[0] << ",";
957       line = Form(" &INP IRAYL=1, NE=%d RHOZ=%g", nel, wt[0]);
958       for(int j=1; j<nel; j++) {
959          out << " " << wt[j] << ",";
960          line += Form(" %g,", wt[j] );
961          if( line.Length() > 60 ) { out << endl; line = ""; }
962       }
963       out << " RHO=" << mat->GetDensity() << ", ";
964       line += Form(" RHO=%g, ", mat->GetDensity());
965       if( line.Length() > 60 ) { out << endl; line = ""; }
966       
967       if( issb_parm[0] > 0 && issb_parm[1] > 0 ) {
968          cout << "Sternheimer parameters found" << endl;
969          out << " ISSB=1, IEV=" << issb_parm[0] << ",";
970          line += Form(" ISSB=1, IEV=%g,", issb_parm[0]);
971          if( line.Length() > 60 ) { out << endl; line = "";  }
972          out << " CBAR=" << issb_parm[1] << ",";
973          line += Form(" CBAR=%g,",issb_parm[1]);
974          if( line.Length() > 60 ) { out << endl; line = "";  }
975          out << " X0=" << issb_parm[2] << ",";
976          line += Form(" X0=%g,", issb_parm[2]);
977          if( line.Length() > 60 ) { out << endl; line = "";  }
978          out << " X1=" << issb_parm[3] << ",";
979          line += Form(" X1=%g,", issb_parm[3]);
980          if( line.Length() > 60 ) { out << endl; line = "";  }
981          out << " AFACT="<< issb_parm[4] << ",";
982          line += Form(" AFACT=%g,", issb_parm[4]);
983          if( line.Length() > 60 ) { out << endl; line = "";  }
984          out << " SK=" << issb_parm[5] << ",";
985          line += Form(" SK=%g,", issb_parm[5]);
986          if( line.Length() > 60 ) { out << endl; line = "";  }
987       }
988       else {
989          cout << "Sternheimer parameters  not found" << endl;
990          (*NoStern)++;
991       }
992       
993       if (dens<0.01){
994          (*countGas)++;
995          out << " GASP=1." << endl;
996       }
997       
998       out << " &END" <<  endl;
999       out << name.Data() << endl;
1000       for (i=0; i<nel; i++) {
1001          elem = mix->GetElement(i);
1002          line = elem->GetName();
1003          if (line.Length()==1) line.Append(" ");
1004          out << line.Data() << " ";
1005       }
1006       out << endl;
1007       
1008       delete [] zt;
1009       delete [] wt;
1010    }
1011    
1012    Double_t ue = 3000000.; // [MeV]
1013    Double_t up = 3000000.; // [MeV]
1014    Double_t ae = -1.;
1015    Double_t ap = -1.;
1016    
1017    
1018    TObjArray* cutList = ((TFluka*) gMC)->GetListOfUserConfigs();
1019    TIter next(cutList);
1020    TFlukaConfigOption* proc;
1021    
1022    while((proc = (TFlukaConfigOption*)next()))
1023    { 
1024        if (proc->Medium() == mat->GetIndex()) {
1025            ap = proc->Cut(kCUTGAM);
1026            ae = proc->Cut(kCUTELE);
1027            if (ap == -1.) ap =  TFlukaConfigOption::DefaultCut(kCUTGAM);
1028            if (ae == -1.) ae =  TFlukaConfigOption::DefaultCut(kCUTELE);
1029            break;
1030        }
1031    }
1032
1033    if (ap == -1.) ap =  TFlukaConfigOption::DefaultCut(kCUTGAM);
1034    if (ae == -1.) ae =  TFlukaConfigOption::DefaultCut(kCUTELE);
1035
1036    ap *= 1000.;                         // [MeV]
1037    ae  = (ae + 0.00051099906) * 1000.;  // [MeV]
1038    
1039    out << "ENER" << endl;
1040    out << " $INP AE=" << ae << ", UE=" << ue <<", AP=" << ap << ", UP=" << up << " $END" << endl;
1041    out << "PWLF" << endl;
1042    out << " $INP NALE=300, NALG=400, NALR=100 $END" << endl;
1043    out << "DECK" << endl;
1044    out << " $INP $END" << endl;
1045    out << "TEST" << endl;
1046    out << " $INP $END" << endl;
1047    out.close();
1048 }
1049
1050 Double_t * TFlukaMCGeometry::GetISSB(Double_t rho, Int_t nElem, Double_t *zelem, Double_t *welem ) const
1051 {
1052    // Read the density effect parameters
1053    // from R.M. Sternheimer et al. Atomic Data
1054    // and Nuclear Data Tables, Vol. 30 No. 2
1055    //
1056    // return the parameters if the element/mixture match with one of the list
1057    // otherwise returns the parameters set to 0
1058    
1059    struct sternheimerData {
1060       TString     longname;           // element/mixture name
1061       Int_t       nelems;             // number of constituents N
1062       Int_t       Z[20];              //[nelems] Z
1063       Double_t    wt[20];             //[nelems] weight fraction
1064       Double_t    density;            // g/cm3
1065       Double_t    iev;                // Average Ion potential (eV)
1066                                       // ****   Sternheimer parameters  ****
1067       Double_t    cbar;               // CBAR
1068       Double_t    x0;                 // X0
1069       Double_t    x1;                 // X1
1070       Double_t    afact;              // AFACT
1071       Double_t    sk;                 // SK
1072       Double_t    delta0;             // DELTA0
1073    };
1074    
1075    TString     shortname;
1076    TString     formula;
1077    Int_t       num;
1078    char        state;
1079    
1080    static Double_t parameters[7];
1081    memset( parameters, 0, sizeof(Double_t) );
1082
1083    static sternheimerData sternDataArray[300];
1084    static Bool_t isFileRead = kFALSE;
1085    
1086    // Read the data file if is needed
1087    if( isFileRead == kFALSE ) {
1088       TString sSternheimerInp = getenv("ALICE_ROOT");
1089       sSternheimerInp +="/TFluka/input/Sternheimer.data";
1090    
1091       ifstream in(sSternheimerInp);
1092       char line[100];
1093       in.getline(line, 100);   
1094       in.getline(line, 100);   
1095       in.getline(line, 100);   
1096       in.getline(line, 100);   
1097       in.getline(line, 100);   
1098       in.getline(line, 100);   
1099       
1100       
1101       Int_t is = 0;
1102       while( !in.eof() ) {
1103          in >> shortname >> num     >> sternDataArray[is].nelems 
1104             >> sternDataArray[is].longname  >> formula >> state;
1105          if( in.eof() ) break;
1106          for(int i=0; i<sternDataArray[is].nelems; i++) {
1107             in >> sternDataArray[is].Z[i] >> sternDataArray[is].wt[i]; 
1108          }
1109          in >> sternDataArray[is].density; 
1110          in >> sternDataArray[is].iev; 
1111          in >> sternDataArray[is].cbar; 
1112          in >> sternDataArray[is].x0; 
1113          in >> sternDataArray[is].x1; 
1114          in >> sternDataArray[is].afact; 
1115          in >> sternDataArray[is].sk;
1116          if( sternDataArray[is].nelems == 1 ) in >> sternDataArray[is].delta0;
1117          is++;
1118       }
1119       isFileRead = kTRUE;
1120       in.close();
1121    }   
1122    
1123    Int_t is = 0;
1124    while( is < 280 ) {
1125    
1126       // check for elements
1127       if( sternDataArray[is].nelems == 1 && nElem == 1
1128           && sternDataArray[is].Z[0] == Int_t(*zelem)
1129           && TMath::Abs( (sternDataArray[is].density - rho)/sternDataArray[is].density ) < 0.1 ) {
1130          cout << sternDataArray[is].longname << "   #elems:" <<  sternDataArray[is].nelems << "  Rho:" 
1131               << sternDataArray[is].density << endl;
1132          cout << sternDataArray[is].iev   << " " 
1133               << sternDataArray[is].cbar  << " " 
1134               << sternDataArray[is].x0    << " " 
1135               << sternDataArray[is].x1    << " " 
1136               << sternDataArray[is].afact << " " 
1137               << sternDataArray[is].sk    << " " 
1138               << sternDataArray[is].delta0 << endl;
1139          
1140          parameters[0] = sternDataArray[is].iev;
1141          parameters[1] = sternDataArray[is].cbar;
1142          parameters[2] = sternDataArray[is].x0;
1143          parameters[3] = sternDataArray[is].x1;
1144          parameters[4] = sternDataArray[is].afact;
1145          parameters[5] = sternDataArray[is].sk;
1146          parameters[6] = sternDataArray[is].delta0;
1147          return parameters;        
1148       }
1149       
1150       // check for mixture
1151       int nmatch = 0;
1152       if( sternDataArray[is].nelems > 1 && sternDataArray[is].nelems == nElem ) {
1153          for(int j=0; j<sternDataArray[is].nelems; j++) {
1154             if( sternDataArray[is].Z[j] == Int_t(zelem[j]) && 
1155                TMath::Abs( (sternDataArray[is].wt[j] - welem[j])/sternDataArray[is].wt[j] ) < 0.1 )
1156             nmatch++;            
1157          }
1158       }
1159
1160       if( sternDataArray[is].nelems > 1 && 
1161           TMath::Abs( (sternDataArray[is].density - rho)/sternDataArray[is].density ) < 0.1 
1162           && nmatch == sternDataArray[is].nelems ) {
1163          cout << sternDataArray[is].longname << "   #elem:" <<  sternDataArray[is].nelems << "  Rho:" 
1164               << sternDataArray[is].density << endl;
1165          cout << sternDataArray[is].iev   << " " 
1166               << sternDataArray[is].cbar  << " " 
1167               << sternDataArray[is].x0    << " " 
1168               << sternDataArray[is].x1    << " " 
1169               << sternDataArray[is].afact << " " 
1170               << sternDataArray[is].sk    << " " 
1171               << sternDataArray[is].delta0 << endl;
1172
1173          parameters[0] = sternDataArray[is].iev;
1174          parameters[1] = sternDataArray[is].cbar;
1175          parameters[2] = sternDataArray[is].x0;
1176          parameters[3] = sternDataArray[is].x1;
1177          parameters[4] = sternDataArray[is].afact;
1178          parameters[5] = sternDataArray[is].sk;
1179          parameters[6] = 0;
1180          return parameters;        
1181       }
1182       is++; 
1183    }   
1184    return parameters;        
1185 }
1186
1187 //_____________________________________________________________________________
1188 void TFlukaMCGeometry::PrintHeader(ofstream &out, const char *text) const
1189 {
1190 // Print a FLUKA header.
1191   out << "*\n" << "*\n" << "*\n";
1192   out << "*********************  " << text << " *********************\n"
1193      << "*\n";
1194   out << "*...+....1....+....2....+....3....+....4....+....5....+....6....+....7..."
1195      << endl;
1196   out << "*" << endl;
1197 }
1198
1199 //_____________________________________________________________________________
1200 Int_t TFlukaMCGeometry::RegionId() const
1201 {
1202 // Returns current region id <-> TGeo node id
1203    if (gGeoManager->IsOutside()) return 0;
1204    return gGeoManager->GetCurrentNode()->GetUniqueID();
1205 }
1206
1207 //_____________________________________________________________________________
1208 Int_t TFlukaMCGeometry::GetElementIndex(Int_t z) const
1209 {
1210 // Get index of a material having a given Z element.
1211    TIter next(fMatList);
1212    TGeoMaterial *mat;
1213    Int_t index = 0;
1214    while ((mat=(TGeoMaterial*)next())) {
1215       if (mat->IsMixture()) continue;
1216       if (mat->GetElement()->Z() == z) return mat->GetIndex();
1217    }
1218    return index;   
1219 }
1220
1221 //_____________________________________________________________________________
1222 void TFlukaMCGeometry::SetMreg(Int_t mreg)
1223 {
1224 // Update if needed next history;
1225    if (gFluka->GetDummyBoundary()==2) {
1226       gGeoManager->CdNode(fNextLattice-1);
1227       return;
1228    }   
1229    Int_t curreg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
1230    if (mreg==curreg) return;
1231    if (mreg==fNextRegion) {
1232       if (fNextLattice!=999999999) gGeoManager->CdNode(fNextLattice-1);
1233       return;
1234    } else {
1235       if (mreg == fCurrentRegion) {
1236          if (fCurrentLattice!=999999999) gGeoManager->CdNode(fCurrentLattice-1);
1237          return;
1238       }   
1239    }     
1240    if (fDebug) printf("ERROR: mreg=%i neither current nor next region\n", mreg);
1241 }
1242
1243 //_____________________________________________________________________________
1244 void TFlukaMCGeometry::SetCurrentRegion(Int_t mreg, Int_t latt)
1245 {
1246 // Set index/history for next entered region
1247    fCurrentRegion = mreg;
1248    fCurrentLattice = latt;
1249 }   
1250
1251 //_____________________________________________________________________________
1252 void TFlukaMCGeometry::SetNextRegion(Int_t mreg, Int_t latt)
1253 {
1254 // Set index/history for next entered region
1255    fNextRegion = mreg;
1256    fNextLattice = latt;
1257 }   
1258
1259 //_____________________________________________________________________________
1260 void TFlukaMCGeometry::ToFlukaString(TString &str) const
1261 {
1262 // ToFlukaString converts an string to something usefull in FLUKA:
1263 // * Capital letters
1264 // * Only 8 letters
1265 // * Replace ' ' by '_'
1266    if (str.Length()<8) {
1267       str += "        ";
1268    }   
1269    str.Remove(8);
1270    Int_t ilast;
1271    for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break;
1272    str.ToUpper();
1273    for (Int_t pos=0; pos<ilast; pos++)
1274       if (str(pos)==' ') str.Replace(pos,1,"_",1);
1275    return;
1276 }   
1277
1278 //_____________________________________________________________________________
1279 void TFlukaMCGeometry::FlukaMatName(TString &str) const
1280 {
1281 // Convert a name to upper case 8 chars.
1282    ToFlukaString(str);
1283    Int_t ilast;
1284    for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break;
1285    if (ilast>5) ilast = 5;
1286    char number[3];
1287    TIter next(fMatNames);
1288    TObjString *objstr;
1289    TString matname;
1290    Int_t index = 0;
1291    while ((objstr=(TObjString*)next())) {
1292       matname = objstr->GetString();
1293       if (matname == str) {
1294          index++;
1295          if (index<10) {
1296             number[0] = '0';
1297             sprintf(&number[1], "%d", index);
1298          } else if (index<100) {
1299             sprintf(number, "%d", index);            
1300          } else {
1301             Error("FlukaMatName", "Too many materials %s", str.Data());
1302             return;
1303          }
1304          str.Replace(ilast+1, 2, number);
1305          str.Remove(8);
1306       }   
1307    }   
1308 }   
1309          
1310 //______________________________________________________________________________
1311 void TFlukaMCGeometry::Vname(const char *name, char *vname) const
1312 {
1313   //
1314   //  convert name to upper case. Make vname at least 4 chars
1315   //
1316   Int_t l = strlen(name);
1317   Int_t i;
1318   l = l < 4 ? l : 4;
1319   for (i=0;i<l;i++) vname[i] = toupper(name[i]);
1320   for (i=l;i<4;i++) vname[i] = ' ';
1321   vname[4] = 0;      
1322 }
1323
1324
1325 // FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers
1326
1327 //_____________________________________________________________________________
1328 Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/)
1329 {
1330 //   from FLUGG:
1331 // Wrapper for setting DNEAR option on fluka side. Must return 0 
1332 // if user doesn't want Fluka to use DNEAR to compute the 
1333 // step (the same effect is obtained with the GLOBAL (WHAT(3)=-1)
1334 // card in fluka input), returns 1 if user wants Fluka always to 
1335 // use DNEAR (in this case, be sure that GEANT4 DNEAR is unique, 
1336 // coming from all directions!!!)
1337    if (gMCGeom->IsDebugging()) printf("========== Dummy IDNRWR\n");
1338    return 0;
1339 }
1340
1341 //_____________________________________________________________________________
1342 void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz, 
1343           Double_t *pV,  Int_t &oldReg , const Int_t &oldLttc, Double_t &propStep,
1344           Int_t &/*nascFlag*/, Double_t &retStep, Int_t &newReg,
1345                Double_t &saf, Int_t &newLttc, Int_t &lttcFlag,
1346           Double_t *sLt, Int_t *jrLt)
1347
1348 {
1349    // Initialize FLUKa point and direction;
1350    gNstep++;
1351    NORLAT.xn[0] = pSx;
1352    NORLAT.xn[1] = pSy;
1353    NORLAT.xn[2] = pSz;
1354    NORLAT.wn[0] = pV[0];
1355    NORLAT.wn[1] = pV[1];
1356    NORLAT.wn[2] = pV[2];
1357
1358    if (gMCGeom->IsDebugging()) {
1359       printf("========== Inside G1WR\n");
1360       printf("   point/dir:(%14.9f, %14.9f, %14.9f, %g, %g, %g)\n", pSx,pSy,pSz,pV[0],pV[1],pV[2]);
1361       printf("   oldReg=%i  oldLttc=%i  pstep=%f\n",oldReg, oldLttc, propStep);
1362    }
1363    Int_t olttc = oldLttc;
1364    if (oldLttc<0) {
1365       gGeoManager->FindNode(pSx,pSy,pSz);
1366       olttc = gGeoManager->GetCurrentNodeId()+1;
1367    }   
1368    gMCGeom->SetCurrentRegion(oldReg, olttc);
1369    // Initialize default return values
1370    lttcFlag = 0;
1371    jrLt[lttcFlag] = olttc;
1372    sLt[lttcFlag] = propStep;
1373    jrLt[lttcFlag+1] = -1;
1374    sLt[lttcFlag+1] = 0.;
1375    newReg = oldReg;
1376    newLttc = olttc;
1377    // check if dummy boundary flag is set
1378    Int_t curLttc, curReg;
1379    if (gFluka->IsDummyBoundary()) {
1380 //      printf("Dummy boundary intercepted. Point is: %f, %f, %f\n", pSx, pSy, pSz);
1381       Bool_t crossedDummy = (olttc == TFlukaMCGeometry::kLttcVirtual)?kTRUE:kFALSE;
1382       if (crossedDummy) {
1383       // FLUKA crossed the dummy boundary - update new region/history
1384          retStep = 0.;
1385          saf = 0.;
1386          gMCGeom->GetNextRegion(newReg, newLttc);
1387          gMCGeom->SetMreg(newReg);
1388          if (gMCGeom->IsDebugging()) printf("   virtual newReg=%i newLttc=%i\n", newReg, newLttc);
1389          sLt[lttcFlag] = 0.; // null step in current region
1390          lttcFlag++;
1391          jrLt[lttcFlag] = newLttc;
1392          sLt[lttcFlag] = 0.; // null step in next region
1393          jrLt[lttcFlag+1] = -1;
1394          sLt[lttcFlag+1] = 0.;
1395          gFluka->SetDummyBoundary(0);
1396          return;
1397       }   
1398    }   
1399       
1400    // Reset outside flag
1401    gGeoManager->SetOutside(kFALSE);
1402    
1403    // Reset dummy boundary flag
1404    gFluka->SetDummyBoundary(0); 
1405     
1406    curLttc = gGeoManager->GetCurrentNodeId()+1;
1407    curReg = gGeoManager->GetCurrentVolume()->GetNumber();
1408    if (olttc != curLttc) {
1409       // FLUKA crossed the boundary : we trust that the given point is really there,
1410       // so we just update TGeo state
1411       gGeoManager->CdNode(olttc-1);
1412       curLttc = gGeoManager->GetCurrentNodeId()+1;
1413       curReg  = gGeoManager->GetCurrentVolume()->GetNumber();
1414       if (gMCGeom->IsDebugging()) printf("   re-initialized point: curReg=%i  curLttc=%i\n", curReg, curLttc);
1415    }  
1416    // Now the current TGeo state reflects the FLUKA state 
1417    Double_t extra = 1.E-10;
1418
1419      
1420 //      printf("ERROR: (%f, %f, %f)\n",pSx,pSy,pSz);
1421    if (gMCGeom->IsDebugging()) printf("   current path: %s\n", gGeoManager->GetPath());
1422    gGeoManager->SetCurrentPoint(pSx+extra*pV[0], pSy+extra*pV[1], pSz+extra*pV[2]);
1423    gGeoManager->SetCurrentDirection(pV);
1424    gGeoManager->FindNextBoundary(-propStep);
1425    Double_t snext = gGeoManager->GetStep();
1426
1427    if (snext<=0.0) {
1428       saf = 0.0;
1429       newReg = -3;
1430       sLt[lttcFlag] = 0.0;
1431       if (gMCGeom->IsDebugging()) printf("BACK SCATTERING\n");
1432       return;
1433    }  
1434
1435    snext += extra;
1436    saf = gGeoManager->GetSafeDistance();
1437    saf -= extra;
1438    if (saf<0) saf=0.0;
1439    else       saf -= saf*3.0e-09;
1440    NORLAT.distn = snext;
1441    NORLAT.xn[0] += snext*pV[0];
1442    NORLAT.xn[1] += snext*pV[1];
1443    NORLAT.xn[2] += snext*pV[2];
1444 //   saf = 0.0; // !!! TEMPORARY FOR TESTING MAGSPHF - TO BE REMOVED 
1445    if (snext>propStep) {
1446    // Next boundary further than proposed step, which is approved
1447       retStep = propStep;
1448       sLt[lttcFlag] = propStep;
1449       return;
1450    }
1451    // The next boundary is closer. We try to cross it.
1452    gGeoManager->SetCurrentPoint(pSx,pSy,pSz);
1453    Double_t *point = gGeoManager->GetCurrentPoint();
1454    Double_t *dir = gGeoManager->GetCurrentDirection();
1455    Double_t pt[3];
1456    memcpy(pt, point, 3*sizeof(Double_t));
1457    
1458    Int_t i;
1459    for (i=0;i<3;i++) point[i] += (snext+1E-6)*dir[i];
1460    // locate next region
1461    gGeoManager->FindNode();
1462    newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
1463    gGeoManager->SetCurrentPoint(pt);
1464    newReg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
1465    if (gMCGeom->IsDebugging()) printf("   newReg=%i newLttc=%i\n", newReg, newLttc);
1466
1467    // We really crossed the boundary, but is it the same region ?
1468    gMCGeom->SetNextRegion(newReg, newLttc);
1469    if (newReg==oldReg && newLttc!=olttc) {
1470       // Virtual boundary between replicants
1471       if (gMCGeom->IsDebugging()) printf("   DUMMY boundary\n");
1472       newReg = 1;  // cheat FLUKA telling it it crossed the TOP region
1473       newLttc = TFlukaMCGeometry::kLttcVirtual;
1474       // mark that next boundary is virtual
1475       gFluka->SetDummyBoundary(1);
1476    } 
1477    retStep = snext;
1478    sLt[lttcFlag] = snext;
1479    lttcFlag++;
1480    jrLt[lttcFlag] = newLttc;
1481    sLt[lttcFlag] = snext;
1482    jrLt[lttcFlag+1] = -1;
1483    sLt[lttcFlag+1] = 0.;      
1484
1485    if (newLttc!=olttc) {
1486       if (gGeoManager->IsOutside()) gGeoManager->SetOutside(kFALSE);
1487       gGeoManager->CdNode(olttc-1);
1488    }   
1489    if (gMCGeom->IsDebugging()) {
1490       printf("=> snext=%g safe=%g\n", snext, saf);
1491       for (Int_t i=0; i<lttcFlag+1; i++) printf("   jrLt[%i]=%i  sLt[%i]=%g\n", i,jrLt[i],i,sLt[i]);
1492    }   
1493    if (gMCGeom->IsDebugging()) printf("<= G1WR (in: %s)\n", gGeoManager->GetPath());
1494 }
1495
1496 //_____________________________________________________________________________
1497 void g1rtwr()
1498 {
1499    if (gMCGeom->IsDebugging()) printf("========== Dummy G1RTWR\n");
1500
1501
1502 //_____________________________________________________________________________
1503 void conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/)
1504 {
1505    if (gMCGeom->IsDebugging()) printf("========== Dummy CONHWR\n");
1506 }
1507
1508 //_____________________________________________________________________________
1509 void inihwr(Int_t &intHist)
1510 {
1511    if (gMCGeom->IsDebugging()) printf("========== Inside INIHWR -> reinitializing history: %i\n", intHist);
1512    if (gGeoManager->IsOutside()) gGeoManager->CdTop();
1513    if (intHist<=0) {
1514 //      printf("=== wrong history number\n");
1515       return;
1516    }
1517    if (intHist==0) gGeoManager->CdTop();
1518    else gGeoManager->CdNode(intHist-1);
1519    if (gMCGeom->IsDebugging()) {
1520       printf(" --- current path: %s\n", gGeoManager->GetPath());
1521       printf("<= INIHWR\n");
1522    }   
1523 }
1524
1525 //_____________________________________________________________________________
1526 void  jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
1527              Int_t &flukaReg)
1528 {
1529 // Geometry initialization wrapper called by FLUKAM. Provides to FLUKA the
1530 // number of regions (volumes in TGeo)
1531    // build application geometry
1532    if (gMCGeom->IsDebugging()) printf("========== Inside JOMIWR\n");
1533    flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast();
1534    if (gMCGeom->IsDebugging()) printf("<= JOMIWR: last region=%i\n", flukaReg);
1535 }   
1536
1537 //_____________________________________________________________________________
1538 void lkdbwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1539             Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1540             Int_t &newReg, Int_t &flagErr, Int_t &newLttc)             
1541 {
1542    if (gMCGeom->IsDebugging()) {
1543       printf("========== Inside LKDBWR (%f, %f, %f)\n",pSx, pSy, pSz);
1544       printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1545       printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1546    }   
1547    lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,newReg,flagErr,newLttc);
1548 }
1549
1550 //_____________________________________________________________________________
1551 void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1552             Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1553             Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
1554 {
1555    if (gMCGeom->IsDebugging()) {
1556       printf("========== Inside LKFXWR (%f, %f, %f)\n",pSx, pSy, pSz);
1557       printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1558       printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1559    }   
1560    lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,newReg,flagErr,newLttc);
1561 }
1562
1563 //_____________________________________________________________________________
1564 void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1565             Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1566                       Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
1567 {
1568    if (gMCGeom->IsDebugging()) {
1569       printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz);
1570       printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1571       printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1572    }   
1573    lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,newReg,flagErr,newLttc);
1574 }
1575
1576 //_____________________________________________________________________________
1577 void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1578           Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
1579                Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
1580 {
1581    if (gMCGeom->IsDebugging()) {
1582       printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz);
1583       printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
1584       printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
1585    }   
1586    NORLAT.xn[0] = pSx;
1587    NORLAT.xn[1] = pSy;
1588    NORLAT.xn[2] = pSz;
1589    NORLAT.wn[0] = pV[0];
1590    NORLAT.wn[1] = pV[1];
1591    NORLAT.wn[2] = pV[2];
1592    TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
1593    if (gGeoManager->IsOutside()) {
1594       newReg = gMCGeom->NofVolumes()+1;
1595 //      newLttc = gGeoManager->GetCurrentNodeId();
1596       newLttc = 999999999;
1597       if (gMCGeom->IsDebugging()) {
1598          printf("OUTSIDE\n");
1599          printf("  out: newReg=%i newLttc=%i\n", newReg, newLttc);
1600          printf("<= LKMGWR\n");
1601       }   
1602       flagErr = newReg;
1603       return;
1604    } 
1605    newReg = node->GetVolume()->GetNumber();
1606    newLttc = gGeoManager->GetCurrentNodeId()+1;
1607    gMCGeom->SetNextRegion(newReg, newLttc);
1608    flagErr = newReg;
1609    if (gMCGeom->IsDebugging()) {
1610       printf("  out: newReg=%i newLttc=%i in %s\n", newReg, newLttc, gGeoManager->GetPath());
1611       printf("<= LKWR\n");
1612    }   
1613 }
1614
1615 //_____________________________________________________________________________
1616 void nrmlwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
1617             Double_t &pVx, Double_t &pVy, Double_t &pVz,
1618                  Double_t *norml, const Int_t &oldReg, 
1619                  const Int_t &newReg, Int_t &flagErr)
1620 {
1621    if (gMCGeom->IsDebugging()) {
1622       printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz);
1623       printf("   oldReg=%i, newReg=%i\n", oldReg,newReg);
1624    }   
1625    gGeoManager->SetCurrentPoint(NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]);
1626    gGeoManager->SetCurrentDirection(NORLAT.wn[0], NORLAT.wn[1], NORLAT.wn[2]);
1627    Double_t *dnorm = gGeoManager->FindNormalFast();
1628    flagErr = 0;
1629    if (!dnorm) {
1630       printf("   ERROR: Cannot compute fast normal\n");
1631       flagErr = 1;
1632       norml[0] = -pVx;   
1633       norml[1] = -pVy;   
1634       norml[2] = -pVz; 
1635    }
1636    norml[0] = -dnorm[0];   
1637    norml[1] = -dnorm[1];   
1638    norml[2] = -dnorm[2]; 
1639    if (gMCGeom->IsDebugging()) {
1640       printf("   normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]);  
1641       printf("<= NRMLWR\n");
1642    }   
1643 }
1644
1645 //_____________________________________________________________________________
1646 void rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
1647             Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/)
1648 {
1649    if (gMCGeom->IsDebugging()) printf("=> Dummy RGRPWR\n");
1650 }
1651
1652 //_____________________________________________________________________________
1653 Int_t isvhwr(const Int_t &check, const Int_t & intHist)
1654 {
1655 //   from FLUGG:
1656 // Wrapper for saving current navigation history (fCheck=default) 
1657 // and returning its pointer. If fCheck=-1 copy of history pointed 
1658 // by intHist is made in NavHistWithCount object, and its pointer 
1659 // is returned. fCheck=1 and fCheck=2 cases are only in debugging 
1660 // version: an array is created by means of FGeometryInit functions
1661 // (but could be a static int * ptrArray = new int[10000] with 
1662 // file scope as well) that stores a flag for deleted/undeleted 
1663 // histories and at the end of event is checked to verify that 
1664 // all saved history objects have been deleted.
1665
1666 // For TGeo, just return the current node ID. No copy need to be made.
1667
1668    if (gMCGeom->IsDebugging()) printf("=> Inside ISVHWR\n");
1669    if (check<0) return intHist;
1670    Int_t histInt = gGeoManager->GetCurrentNodeId()+1;
1671    if (gMCGeom->IsDebugging()) printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath());
1672    return histInt;
1673 }
1674
1675
1676
1677