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