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