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