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