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