6acbbfc5032317a007621f7befc0e9036dc626c7
[u/mrichter/AliRoot.git] / TFluka / TFluka.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 //
19 // Realisation of the TVirtualMC interface for the FLUKA code
20 // (See official web side http://www.fluka.org/).
21 //
22 // This implementation makes use of the TGeo geometry modeller.
23 // User configuration is via automatic generation of FLUKA input cards.
24 //
25 // Authors:
26 // A. Fasso
27 // E. Futo
28 // A. Gheata
29 // A. Morsch
30 //
31
32 #include <Riostream.h>
33
34 #include "TFluka.h"
35 #include "TFlukaCodes.h"
36 #include "TCallf77.h"      //For the fortran calls
37 #include "Fdblprc.h"       //(DBLPRC) fluka common
38 #include "Fsourcm.h"       //(SOURCM) fluka common
39 #include "Fgenstk.h"       //(GENSTK)  fluka common
40 #include "Fiounit.h"       //(IOUNIT) fluka common
41 #include "Fpaprop.h"       //(PAPROP) fluka common
42 #include "Fpart.h"         //(PART)   fluka common
43 #include "Ftrackr.h"       //(TRACKR) fluka common
44 #include "Fpaprop.h"       //(PAPROP) fluka common
45 #include "Ffheavy.h"       //(FHEAVY) fluka common
46 #include "Fopphst.h"       //(OPPHST) fluka common
47 #include "Fflkstk.h"       //(FLKSTK) fluka common
48 #include "Fstepsz.h"       //(STEPSZ) fluka common
49 #include "Fopphst.h"       //(OPPHST) fluka common
50 #include "Fltclcm.h"       //(LTCLCM) fluka common
51
52 #include "TVirtualMC.h"
53 #include "TMCProcess.h"
54 #include "TGeoManager.h"
55 #include "TGeoMaterial.h"
56 #include "TGeoMedium.h"
57 #include "TFlukaMCGeometry.h"
58 #include "TGeoMCGeometry.h"
59 #include "TFlukaCerenkov.h"
60 #include "TFlukaConfigOption.h"
61 #include "TFlukaScoringOption.h"
62 #include "TLorentzVector.h"
63 #include "TArrayI.h"
64 #include "TArrayD.h"
65
66 // Fluka methods that may be needed.
67 #ifndef WIN32 
68 # define flukam  flukam_
69 # define fluka_openinp fluka_openinp_
70 # define fluka_openout fluka_openout_
71 # define fluka_closeinp fluka_closeinp_
72 # define mcihad mcihad_
73 # define mpdgha mpdgha_
74 # define newplo newplo_
75 #else 
76 # define flukam  FLUKAM
77 # define fluka_openinp FLUKA_OPENINP
78 # define fluka_openout FLUKA_OPENOUT
79 # define fluka_closeinp FLUKA_CLOSEINP
80 # define mcihad MCIHAD
81 # define mpdgha MPDGHA
82 # define newplo NEWPLO
83 #endif
84
85 extern "C" 
86 {
87   //
88   // Prototypes for FLUKA functions
89   //
90   void type_of_call flukam(const int&);
91   void type_of_call newplo();
92   void type_of_call fluka_openinp(const int&, DEFCHARA);
93   void type_of_call fluka_openout(const int&, DEFCHARA);
94   void type_of_call fluka_closeinp(const int&);
95   int  type_of_call mcihad(const int&);
96   int  type_of_call mpdgha(const int&);
97 }
98
99 //
100 // Class implementation for ROOT
101 //
102 ClassImp(TFluka)
103
104 //
105 //----------------------------------------------------------------------------
106 // TFluka constructors and destructors.
107 //______________________________________________________________________________
108 TFluka::TFluka()
109   :TVirtualMC(),
110    fVerbosityLevel(0),
111    fInputFileName(""),
112    fUserConfig(0), 
113    fUserScore(0)
114
115   //
116   // Default constructor
117   //
118    fGeneratePemf = kFALSE;
119    fNVolumes = 0;
120    fCurrentFlukaRegion = -1;
121    fNewReg = -1;
122    fGeom = 0;
123    fMCGeo = 0;
124    fMaterials = 0;
125    fDummyBoundary = 0;
126    fFieldFlag = 1;
127    fStopped   = 0;
128    fStopEvent = 0;
129    fStopRun   = 0;
130    fNEvent    = 0;
131
132  
133 //______________________________________________________________________________ 
134 TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
135   :TVirtualMC("TFluka",title, isRootGeometrySupported),
136    fVerbosityLevel(verbosity),
137    fInputFileName(""),
138    fTrackIsEntering(0),
139    fTrackIsExiting(0),
140    fTrackIsNew(0),
141    fUserConfig(new TObjArray(100)),
142    fUserScore(new TObjArray(100)) 
143 {
144   // create geometry interface
145    if (fVerbosityLevel >=3)
146        cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
147    SetCoreInputFileName();
148    SetInputFileName();
149    SetGeneratePemf(kFALSE);
150    fNVolumes      = 0;
151    fCurrentFlukaRegion = -1;
152    fNewReg = -1;
153    fDummyBoundary = 0;
154    fFieldFlag = 1;
155    fGeneratePemf = kFALSE;
156    fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kTRUE);
157    fGeom  = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
158    if (verbosity > 2) fGeom->SetDebugMode(kTRUE);
159    fMaterials = 0;
160    fStopped   = 0;
161    fStopEvent = 0;
162    fStopRun   = 0;
163    fNEvent    = 0;
164    PrintHeader();
165 }
166
167 //______________________________________________________________________________ 
168 TFluka::~TFluka() {
169 // Destructor
170     if (fVerbosityLevel >=3)
171         cout << "<== TFluka::~TFluka() destructor called." << endl;
172     
173     delete fGeom;
174     delete fMCGeo;
175     
176     if (fUserConfig) {
177         fUserConfig->Delete();
178         delete fUserConfig;
179     }
180     
181     if (fUserScore) {
182         fUserScore->Delete();
183         delete fUserScore;
184     }
185 }
186
187 //
188 //______________________________________________________________________________
189 // TFluka control methods
190 //______________________________________________________________________________ 
191 void TFluka::Init() {
192 //
193 //  Geometry initialisation
194 //
195     if (fVerbosityLevel >=3) cout << "==> TFluka::Init() called." << endl;
196     
197     if (!gGeoManager) new TGeoManager("geom", "FLUKA geometry");
198     fApplication->ConstructGeometry();
199     if (!gGeoManager->IsClosed()) {
200        TGeoVolume *top = (TGeoVolume*)gGeoManager->GetListOfVolumes()->First();
201        gGeoManager->SetTopVolume(top);
202        gGeoManager->CloseGeometry("di");
203     } else {
204        TGeoNodeCache *cache = gGeoManager->GetCache();
205        if (!cache->HasIdArray()) {
206           Warning("Init", "Node ID tracking must be enabled with TFluka: enabling...\n");
207           cache->BuildIdArray();
208        }   
209     }           
210     fNVolumes = fGeom->NofVolumes();
211     fGeom->CreateFlukaMatFile("flukaMat.inp");   
212     if (fVerbosityLevel >=3) {
213        printf("== Number of volumes: %i\n ==", fNVolumes);
214        cout << "\t* InitPhysics() - Prepare input file to be called" << endl; 
215     }
216
217     fApplication->InitGeometry();
218     
219 }
220
221
222 //______________________________________________________________________________ 
223 void TFluka::FinishGeometry() {
224 //
225 // Build-up table with region to medium correspondance
226 //
227   if (fVerbosityLevel >=3) {
228     cout << "==> TFluka::FinishGeometry() called." << endl;
229     printf("----FinishGeometry - nothing to do with TGeo\n");
230     cout << "<== TFluka::FinishGeometry() called." << endl;
231   }  
232
233
234 //______________________________________________________________________________ 
235 void TFluka::BuildPhysics() {
236 //
237 //  Prepare FLUKA input files and call FLUKA physics initialisation
238 //
239     
240     if (fVerbosityLevel >=3)
241         cout << "==> TFluka::BuildPhysics() called." << endl;
242
243     
244     if (fVerbosityLevel >=3) {
245         TList *medlist = gGeoManager->GetListOfMedia();
246         TIter next(medlist);
247         TGeoMedium*   med = 0x0;
248         TGeoMaterial* mat = 0x0;
249         Int_t ic = 0;
250         
251         while((med = (TGeoMedium*)next()))
252         {
253             mat = med->GetMaterial();
254             printf("Medium %5d %12s %5d %5d\n", ic, (med->GetName()), med->GetId(), mat->GetIndex());
255             ic++;
256         }
257     }
258     
259     //
260     // At this stage we have the information on materials and cuts available.
261     // Now create the pemf file
262     
263     if (fGeneratePemf) fGeom->CreatePemfFile();
264     
265     //
266     // Prepare input file with the current physics settings
267     
268     InitPhysics(); 
269 //  Open fortran files    
270     const char* fname = fInputFileName;
271     fluka_openinp(lunin, PASSCHARA(fname));
272     fluka_openout(11, PASSCHARA("fluka.out"));
273 //  Read input cards    
274     GLOBAL.lfdrtr = true;
275     flukam(1);
276 //  Close input file
277     fluka_closeinp(lunin);
278 //  Finish geometry    
279     FinishGeometry();
280 }  
281
282 //______________________________________________________________________________ 
283 void TFluka::ProcessEvent() {
284 //
285 // Process one event
286 //
287     if (fStopRun) {
288         Warning("ProcessEvent", "User Run Abortion: No more events handled !\n");
289         fNEvent += 1;
290         return;
291     }
292
293     if (fVerbosityLevel >=3)
294         cout << "==> TFluka::ProcessEvent() called." << endl;
295     fApplication->GeneratePrimaries();
296     SOURCM.lsouit = true;
297     flukam(1);
298     if (fVerbosityLevel >=3)
299         cout << "<== TFluka::ProcessEvent() called." << endl;
300     //
301     // Increase event number
302     //
303     fNEvent += 1;
304 }
305
306 //______________________________________________________________________________ 
307 Bool_t TFluka::ProcessRun(Int_t nevent) {
308 //
309 // Run steering
310 //
311
312   if (fVerbosityLevel >=3)
313     cout << "==> TFluka::ProcessRun(" << nevent << ") called." 
314          << endl;
315
316   if (fVerbosityLevel >=2) {
317     cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
318     cout << "\t* Calling flukam again..." << endl;
319   }
320
321   Int_t todo = TMath::Abs(nevent);
322   for (Int_t ev = 0; ev < todo; ev++) {
323       fApplication->BeginEvent();
324       ProcessEvent();
325       fApplication->FinishEvent();
326   }
327
328   if (fVerbosityLevel >=3)
329     cout << "<== TFluka::ProcessRun(" << nevent << ") called." 
330          << endl;
331   // Write fluka specific scoring output
332   newplo();
333   
334   return kTRUE;
335 }
336
337 //_____________________________________________________________________________
338 // methods for building/management of geometry
339
340 // functions from GCONS 
341 //____________________________________________________________________________ 
342 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,  
343                     Float_t &dens, Float_t &radl, Float_t &absl,
344                     Float_t* /*ubuf*/, Int_t& /*nbuf*/) {
345 //
346    TGeoMaterial *mat;
347    TIter next (gGeoManager->GetListOfMaterials());
348    while ((mat = (TGeoMaterial*)next())) {
349      if (mat->GetUniqueID() == (UInt_t)imat) break;
350    }
351    if (!mat) {
352       Error("Gfmate", "no material with index %i found", imat);
353       return;
354    }
355    sprintf(name, "%s", mat->GetName());
356    a = mat->GetA();
357    z = mat->GetZ();
358    dens = mat->GetDensity();
359    radl = mat->GetRadLen();
360    absl = mat->GetIntLen();
361
362
363 //______________________________________________________________________________ 
364 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,  
365                     Double_t &dens, Double_t &radl, Double_t &absl,
366                     Double_t* /*ubuf*/, Int_t& /*nbuf*/) {
367 //
368    TGeoMaterial *mat;
369    TIter next (gGeoManager->GetListOfMaterials());
370    while ((mat = (TGeoMaterial*)next())) {
371      if (mat->GetUniqueID() == (UInt_t)imat) break;
372    }
373    if (!mat) {
374       Error("Gfmate", "no material with index %i found", imat);
375       return;
376    }
377    sprintf(name, "%s", mat->GetName());
378    a = mat->GetA();
379    z = mat->GetZ();
380    dens = mat->GetDensity();
381    radl = mat->GetRadLen();
382    absl = mat->GetIntLen();
383
384
385 // detector composition
386 //______________________________________________________________________________ 
387 void TFluka::Material(Int_t& kmat, const char* name, Double_t a, 
388                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
389                       Float_t* buf, Int_t nwbuf) {
390 //
391    Double_t* dbuf = fGeom->CreateDoubleArray(buf, nwbuf);  
392    Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf);
393    delete [] dbuf;
394
395
396 //______________________________________________________________________________ 
397 void TFluka::Material(Int_t& kmat, const char* name, Double_t a, 
398                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
399                       Double_t* /*buf*/, Int_t /*nwbuf*/) {
400 //
401 // Define a material
402   TGeoMaterial *mat;
403   kmat = gGeoManager->GetListOfMaterials()->GetSize();
404   if ((z-Int_t(z)) > 1E-3) {
405      mat = fGeom->GetMakeWrongMaterial(z);
406      if (mat) {
407         mat->SetRadLen(radl,absl);
408         mat->SetUniqueID(kmat);
409         return;
410      }
411   }      
412   gGeoManager->Material(name, a, z, dens, kmat, radl, absl);
413
414
415 //______________________________________________________________________________ 
416 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a, 
417                      Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
418 //
419 // Define a material mixture
420 //
421   Double_t* da = fGeom->CreateDoubleArray(a, TMath::Abs(nlmat));  
422   Double_t* dz = fGeom->CreateDoubleArray(z, TMath::Abs(nlmat));  
423   Double_t* dwmat = fGeom->CreateDoubleArray(wmat, TMath::Abs(nlmat));  
424
425   Mixture(kmat, name, da, dz, dens, nlmat, dwmat);
426   for (Int_t i=0; i<nlmat; i++) {
427     a[i] = da[i]; z[i] = dz[i]; wmat[i] = dwmat[i];
428   }  
429
430   delete [] da;
431   delete [] dz;
432   delete [] dwmat;
433
434
435 //______________________________________________________________________________ 
436 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a, 
437                      Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
438 //
439   // Defines mixture OR COMPOUND IMAT as composed by 
440   // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
441   // 
442   // If NLMAT > 0 then wmat contains the proportion by
443   // weights of each basic material in the mixture. 
444   // 
445   // If nlmat < 0 then WMAT contains the number of atoms 
446   // of a given kind into the molecule of the COMPOUND
447   // In this case, WMAT in output is changed to relative
448   // weigths.
449   //
450   Int_t i,j;
451   if (nlmat < 0) {
452      nlmat = - nlmat;
453      Double_t amol = 0;
454      for (i=0;i<nlmat;i++) {
455         amol += a[i]*wmat[i];
456      }
457      for (i=0;i<nlmat;i++) {
458         wmat[i] *= a[i]/amol;
459      }
460   }
461   kmat = gGeoManager->GetListOfMaterials()->GetSize();
462   // Check if we have elements with fractional Z
463   TGeoMaterial *mat = 0;
464   TGeoMixture *mix = 0;
465   Bool_t mixnew = kFALSE;
466   for (i=0; i<nlmat; i++) {
467      if (z[i]-Int_t(z[i]) < 1E-3) continue;
468      // We have found an element with fractional Z -> loop mixtures to look for it
469      for (j=0; j<kmat; j++) {
470         mat = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->At(j);
471         if (!mat) break;
472         if (!mat->IsMixture()) continue;
473         mix = (TGeoMixture*)mat;
474         if (TMath::Abs(z[i]-mix->GetZ()) >1E-3) continue;
475         mixnew = kTRUE;
476         break;
477      }
478      if (!mixnew) Warning("Mixture","%s : cannot find component %i with fractional Z=%f\n", name, i, z[i]);
479      break;
480   }   
481   if (mixnew) {
482      Int_t nlmatnew = nlmat+mix->GetNelements()-1;
483      Double_t *anew = new Double_t[nlmatnew];
484      Double_t *znew = new Double_t[nlmatnew];
485      Double_t *wmatnew = new Double_t[nlmatnew];
486      Int_t ind=0;
487      for (j=0; j<nlmat; j++) {
488         if (j==i) continue;
489         anew[ind] = a[j];
490         znew[ind] = z[j];
491         wmatnew[ind] = wmat[j];
492         ind++;
493      }
494      for (j=0; j<mix->GetNelements(); j++) {
495         anew[ind] = mix->GetAmixt()[j];
496         znew[ind] = mix->GetZmixt()[j];
497         wmatnew[ind] = wmat[i]*mix->GetWmixt()[j];
498         ind++;
499      }
500      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
501      delete [] anew;
502      delete [] znew;
503      delete [] wmatnew;
504      return;
505   }   
506   // Now we need to compact identical elements within the mixture
507   // First check if this happens   
508   mixnew = kFALSE;  
509   for (i=0; i<nlmat-1; i++) {
510      for (j=i+1; j<nlmat; j++) {
511         if (z[i] == z[j]) {
512            mixnew = kTRUE;
513            break;
514         }
515      }   
516      if (mixnew) break;
517   }   
518   if (mixnew) {
519      Int_t nlmatnew = 0;
520      Double_t *anew = new Double_t[nlmat];
521      Double_t *znew = new Double_t[nlmat];
522      memset(znew, 0, nlmat*sizeof(Double_t));
523      Double_t *wmatnew = new Double_t[nlmat];
524      Bool_t skipi;
525      for (i=0; i<nlmat; i++) {
526         skipi = kFALSE;
527         for (j=0; j<nlmatnew; j++) {
528            if (z[i] == z[j]) {
529               wmatnew[j] += wmat[i];
530               skipi = kTRUE;
531               break;
532            }
533         }   
534         if (skipi) continue;    
535         anew[nlmatnew] = a[i];
536         znew[nlmatnew] = z[i];
537         wmatnew[nlmatnew] = wmat[i];
538         nlmatnew++;
539      }
540      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
541      delete [] anew;
542      delete [] znew;
543      delete [] wmatnew;
544      return;     
545    }
546    gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
547
548
549 //______________________________________________________________________________ 
550 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, 
551                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, 
552                     Double_t stemax, Double_t deemax, Double_t epsil, 
553                     Double_t stmin, Float_t* ubuf, Int_t nbuf) { 
554   // Define a medium
555   // 
556   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
557   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, 
558              epsil, stmin, ubuf, nbuf);
559
560
561 //______________________________________________________________________________ 
562 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, 
563                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, 
564                     Double_t stemax, Double_t deemax, Double_t epsil, 
565                     Double_t stmin, Double_t* ubuf, Int_t nbuf) { 
566   // Define a medium
567   // 
568   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
569   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, 
570              epsil, stmin, ubuf, nbuf);
571
572
573 //______________________________________________________________________________ 
574 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX, 
575                     Double_t thetaY, Double_t phiY, Double_t thetaZ, 
576                     Double_t phiZ) {
577 //                   
578   krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
579   fMCGeo->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ); 
580
581
582 //______________________________________________________________________________ 
583 void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
584 //
585 //
586 //
587    Bool_t process = kFALSE;
588    if (strncmp(param, "DCAY",  4) == 0 ||
589        strncmp(param, "PAIR",  4) == 0 ||
590        strncmp(param, "COMP",  4) == 0 ||
591        strncmp(param, "PHOT",  4) == 0 ||
592        strncmp(param, "PFIS",  4) == 0 ||
593        strncmp(param, "DRAY",  4) == 0 ||
594        strncmp(param, "ANNI",  4) == 0 ||
595        strncmp(param, "BREM",  4) == 0 ||
596        strncmp(param, "MUNU",  4) == 0 ||
597        strncmp(param, "CKOV",  4) == 0 ||
598        strncmp(param, "HADR",  4) == 0 ||
599        strncmp(param, "LOSS",  4) == 0 ||
600        strncmp(param, "MULS",  4) == 0 ||
601        strncmp(param, "RAYL",  4) == 0) 
602    {
603        process = kTRUE;
604    } 
605    
606    if (process) {
607        SetProcess(param, Int_t (parval), itmed);
608    } else {
609        SetCut(param, parval, itmed);
610    }
611 }    
612
613 // functions from GGEOM 
614 //_____________________________________________________________________________
615 void TFluka::Gsatt(const char *name, const char *att, Int_t val)
616
617   // Set visualisation attributes for one volume
618   char vname[5];
619   fGeom->Vname(name,vname);
620   char vatt[5];
621   fGeom->Vname(att,vatt);
622   gGeoManager->SetVolumeAttribute(vname, vatt, val);
623 }
624
625 //______________________________________________________________________________ 
626 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,  
627                      Float_t *upar, Int_t np)  {
628 //
629     return fMCGeo->Gsvolu(name, shape, nmed, upar, np); 
630 }
631
632 //______________________________________________________________________________ 
633 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,  
634                      Double_t *upar, Int_t np)  {
635 //
636     return fMCGeo->Gsvolu(name, shape, nmed, upar, np); 
637 }
638  
639 //______________________________________________________________________________ 
640 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv, 
641                    Int_t iaxis) {
642 //
643     fMCGeo->Gsdvn(name, mother, ndiv, iaxis); 
644
645
646 //______________________________________________________________________________ 
647 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv, 
648                     Int_t iaxis, Double_t c0i, Int_t numed) {
649 //
650     fMCGeo->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed); 
651
652
653 //______________________________________________________________________________ 
654 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step, 
655                    Int_t iaxis, Int_t numed, Int_t ndvmx) {
656 //      
657     fMCGeo->Gsdvt(name, mother, step, iaxis, numed, ndvmx); 
658
659
660 //______________________________________________________________________________ 
661 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step, 
662                     Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) { 
663 //
664     fMCGeo->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx); 
665
666
667 //______________________________________________________________________________ 
668 void TFluka::Gsord(const char * /*name*/, Int_t /*iax*/) {
669 //
670 // Nothing to do with TGeo
671
672
673 //______________________________________________________________________________ 
674 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,  
675                    Double_t x, Double_t y, Double_t z, Int_t irot, 
676                    const char *konly) {
677 //
678   fMCGeo->Gspos(name, nr, mother, x, y, z, irot, konly); 
679
680
681 //______________________________________________________________________________ 
682 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,  
683                     Double_t x, Double_t y, Double_t z, Int_t irot,
684                     const char *konly, Float_t *upar, Int_t np)  {
685   //
686   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); 
687
688
689 //______________________________________________________________________________ 
690 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,  
691                     Double_t x, Double_t y, Double_t z, Int_t irot,
692                     const char *konly, Double_t *upar, Int_t np)  {
693   //
694   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); 
695
696
697 //______________________________________________________________________________ 
698 void TFluka::Gsbool(const char* /*onlyVolName*/, const char* /*manyVolName*/) {
699 //
700 // Nothing to do with TGeo
701 }
702
703 //______________________________________________________________________
704 Bool_t TFluka::GetTransformation(const TString &volumePath,TGeoHMatrix &mat)
705 {
706     // Returns the Transformation matrix between the volume specified
707     // by the path volumePath and the Top or mater volume. The format
708     // of the path volumePath is as follows (assuming ALIC is the Top volume)
709     // "/ALIC_1/DDIP_1/S05I_2/S05H_1/S05G_3". Here ALIC is the top most
710     // or master volume which has only 1 instance of. Of all of the daughter
711     // volumes of ALICE, DDIP volume copy #1 is indicated. Similarly for
712     // the daughter volume of DDIP is S05I copy #2 and so on.
713     // Inputs:
714     //   TString& volumePath  The volume path to the specific volume
715     //                        for which you want the matrix. Volume name
716     //                        hierarchy is separated by "/" while the
717     //                        copy number is appended using a "_".
718     // Outputs:
719     //  TGeoHMatrix &mat      A matrix with its values set to those
720     //                        appropriate to the Local to Master transformation
721     // Return:
722     //   A logical value if kFALSE then an error occurred and no change to
723     //   mat was made.
724
725    // We have to preserve the modeler state
726    return fMCGeo->GetTransformation(volumePath, mat);
727 }   
728    
729 //______________________________________________________________________
730 Bool_t TFluka::GetShape(const TString &volumePath,TString &shapeType,
731                         TArrayD &par)
732 {
733     // Returns the shape and its parameters for the volume specified
734     // by volumeName.
735     // Inputs:
736     //   TString& volumeName  The volume name
737     // Outputs:
738     //   TString &shapeType   Shape type
739     //   TArrayD &par         A TArrayD of parameters with all of the
740     //                        parameters of the specified shape.
741     // Return:
742     //   A logical indicating whether there was an error in getting this
743     //   information
744    return fMCGeo->GetShape(volumePath, shapeType, par);
745 }
746    
747 //______________________________________________________________________
748 Bool_t TFluka::GetMaterial(const TString &volumeName,
749                             TString &name,Int_t &imat,
750                             Double_t &a,Double_t &z,Double_t &dens,
751                             Double_t &radl,Double_t &inter,TArrayD &par)
752 {
753     // Returns the Material and its parameters for the volume specified
754     // by volumeName.
755     // Note, Geant3 stores and uses mixtures as an element with an effective
756     // Z and A. Consequently, if the parameter Z is not integer, then
757     // this material represents some sort of mixture.
758     // Inputs:
759     //   TString& volumeName  The volume name
760     // Outputs:
761     //   TSrting   &name       Material name
762     //   Int_t     &imat       Material index number
763     //   Double_t  &a          Average Atomic mass of material
764     //   Double_t  &z          Average Atomic number of material
765     //   Double_t  &dens       Density of material [g/cm^3]
766     //   Double_t  &radl       Average radiation length of material [cm]
767     //   Double_t  &inter      Average interaction length of material [cm]
768     //   TArrayD   &par        A TArrayD of user defined parameters.
769     // Return:
770     //   kTRUE if no errors
771    return fMCGeo->GetMaterial(volumeName,name,imat,a,z,dens,radl,inter,par);
772 }
773
774 //______________________________________________________________________
775 Bool_t TFluka::GetMedium(const TString &volumeName,TString &name,
776                          Int_t &imed,Int_t &nmat,Int_t &isvol,Int_t &ifield,
777                          Double_t &fieldm,Double_t &tmaxfd,Double_t &stemax,
778                          Double_t &deemax,Double_t &epsil, Double_t &stmin,
779                          TArrayD &par)
780 {
781     // Returns the Medium and its parameters for the volume specified
782     // by volumeName.
783     // Inputs:
784     //   TString& volumeName  The volume name.
785     // Outputs:
786     //   TString  &name       Medium name
787     //   Int_t    &nmat       Material number defined for this medium
788     //   Int_t    &imed       The medium index number
789     //   Int_t    &isvol      volume number defined for this medium
790     //   Int_t    &iflield    Magnetic field flag
791     //   Double_t &fieldm     Magnetic field strength
792     //   Double_t &tmaxfd     Maximum angle of deflection per step
793     //   Double_t &stemax     Maximum step size
794     //   Double_t &deemax     Maximum fraction of energy allowed to be lost
795     //                        to continuous process.
796     //   Double_t &epsil      Boundary crossing precision
797     //   Double_t &stmin      Minimum step size allowed
798     //   TArrayD  &par        A TArrayD of user parameters with all of the
799     //                        parameters of the specified medium.
800     // Return:
801     //   kTRUE if there where no errors
802    return fMCGeo->GetMedium(volumeName,name,imed,nmat,isvol,ifield,fieldm,tmaxfd,stemax,deemax,epsil,stmin,par);
803 }         
804
805 //______________________________________________________________________________ 
806 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
807                          Float_t* absco, Float_t* effic, Float_t* rindex) {
808 //
809 // Set Cerenkov properties for medium itmed
810 //
811 // npckov: number of sampling points
812 // ppckov: energy values
813 // absco:  absorption length
814 // effic:  quantum efficiency
815 // rindex: refraction index
816 //
817 //
818 //  
819 //  Create object holding Cerenkov properties
820 //  
821     TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex);
822 //
823 //  Pass object to medium
824     TGeoMedium* medium = gGeoManager->GetMedium(itmed);
825     medium->SetCerenkovProperties(cerenkovProperties);
826 }  
827
828 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
829                          Float_t* absco, Float_t* effic, Float_t* rindex, Float_t* rfl) {
830 //
831 // Set Cerenkov properties for medium itmed
832 //
833 // npckov: number of sampling points
834 // ppckov: energy values
835 // absco:  absorption length
836 // effic:  quantum efficiency
837 // rindex: refraction index
838 // rfl:    reflectivity for boundary to medium itmed
839 //
840 //  
841 //  Create object holding Cerenkov properties
842 //  
843     TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex, rfl);
844 //
845 //  Pass object to medium
846     TGeoMedium* medium = gGeoManager->GetMedium(itmed);
847     medium->SetCerenkovProperties(cerenkovProperties);
848 }  
849
850
851 //______________________________________________________________________________ 
852 void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t * /*ppckov*/,
853                          Double_t * /*absco*/, Double_t * /*effic*/, Double_t * /*rindex*/) {
854 //
855 //  Double_t version not implemented
856 }  
857
858 void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t* /*ppckov*/,
859                          Double_t* /*absco*/, Double_t* /*effic*/, Double_t* /*rindex*/, Double_t* /*rfl*/) { 
860 //
861 // //  Double_t version not implemented
862 }
863
864 // Euclid
865 //______________________________________________________________________________ 
866 void TFluka::WriteEuclid(const char* /*fileName*/, const char* /*topVol*/, 
867                           Int_t /*number*/, Int_t /*nlevel*/) {
868 //
869 // Not with TGeo
870    Warning("WriteEuclid", "Not implemented !");
871
872
873
874
875 //_____________________________________________________________________________
876 // methods needed by the stepping
877 //____________________________________________________________________________ 
878
879 Int_t TFluka::GetMedium() const {
880 //
881 //  Get the medium number for the current fluka region
882 //
883     return fGeom->GetMedium(); // this I need to check due to remapping !!!
884 }
885
886 //____________________________________________________________________________ 
887 Int_t TFluka::GetDummyRegion() const
888 {
889 // Returns index of the dummy region.
890    return fGeom->GetDummyRegion();
891 }   
892
893 //____________________________________________________________________________ 
894 Int_t TFluka::GetDummyLattice() const
895 {
896 // Returns index of the dummy lattice.
897    return fGeom->GetDummyLattice();
898 }   
899
900 //____________________________________________________________________________ 
901 // particle table usage
902 // ID <--> PDG transformations
903 //_____________________________________________________________________________
904 Int_t TFluka::IdFromPDG(Int_t pdg) const 
905 {
906     //
907     // Return Fluka code from PDG and pseudo ENDF code
908     
909     // Catch the feedback photons
910     if (pdg == 50000051) return (kFLUKAoptical);
911     // MCIHAD() goes from pdg to fluka internal.
912     Int_t intfluka = mcihad(pdg);
913     // KPTOIP array goes from internal to official
914     return GetFlukaKPTOIP(intfluka);
915 }
916
917 //______________________________________________________________________________ 
918 Int_t TFluka::PDGFromId(Int_t id) const 
919 {
920   //
921   // Return PDG code and pseudo ENDF code from Fluka code
922   //                      Alpha     He3       Triton    Deuteron  gen. ion  opt. photon   
923     Int_t idSpecial[6] = {10020040, 10020030, 10010030, 10010020, 10000000, 50000050};
924   // IPTOKP array goes from official to internal
925
926     if (id == kFLUKAoptical) {
927 // Cerenkov photon
928         if (fVerbosityLevel >= 3)
929             printf("\n PDGFromId: Cerenkov Photon \n");
930         return  50000050;
931     }
932 // Error id    
933     if (id == 0 || id < kFLUKAcodemin || id > kFLUKAcodemax) {
934         if (fVerbosityLevel >= 3)
935             printf("PDGFromId: Error id = 0\n");
936         return -1;
937     }
938 // Good id    
939     if (id > 0) {
940         Int_t intfluka = GetFlukaIPTOKP(id);
941         if (intfluka == 0) {
942             if (fVerbosityLevel >= 3)
943                 printf("PDGFromId: Error intfluka = 0: %d\n", id);
944             return -1;
945         } else if (intfluka < 0) {
946             if (fVerbosityLevel >= 3)
947                 printf("PDGFromId: Error intfluka < 0: %d\n", id);
948             return -1;
949         }
950         if (fVerbosityLevel >= 3)
951             printf("mpdgha called with %d %d \n", id, intfluka);
952         return mpdgha(intfluka);
953     } else {
954         // ions and optical photons
955         return idSpecial[id - kFLUKAcodemin];
956     }
957 }
958
959 void TFluka::StopTrack()
960 {
961     // Set stopping conditions
962     // Works for photons and charged particles
963     fStopped = kTRUE;
964 }
965   
966 //_____________________________________________________________________________
967 // methods for physics management
968 //____________________________________________________________________________ 
969 //
970 // set methods
971 //
972
973 void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imed)
974 {
975 //  Set process user flag for material imat
976 //
977 //    
978 //  Update if already in the list
979 //
980     TIter next(fUserConfig);
981     TFlukaConfigOption* proc;
982     while((proc = (TFlukaConfigOption*)next()))
983     { 
984         if (proc->Medium() == imed) {
985             proc->SetProcess(flagName, flagValue);
986             return;
987         }
988     }
989     proc = new TFlukaConfigOption(imed);
990     proc->SetProcess(flagName, flagValue);
991     fUserConfig->Add(proc);
992 }
993
994 //______________________________________________________________________________ 
995 Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue)
996 {
997 //  Set process user flag 
998 //
999 //    
1000     SetProcess(flagName, flagValue, -1);
1001     return kTRUE;  
1002 }
1003
1004 //______________________________________________________________________________ 
1005 void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
1006 {
1007 // Set user cut value for material imed
1008 //
1009     TIter next(fUserConfig);
1010     TFlukaConfigOption* proc;
1011     while((proc = (TFlukaConfigOption*)next()))
1012     { 
1013         if (proc->Medium() == imed) {
1014             proc->SetCut(cutName, cutValue);
1015             return;
1016         }
1017     }
1018
1019     proc = new TFlukaConfigOption(imed);
1020     proc->SetCut(cutName, cutValue);
1021     fUserConfig->Add(proc);
1022 }
1023
1024 //______________________________________________________________________________ 
1025 Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue)
1026 {
1027 // Set user cut value 
1028 //
1029 //    
1030     SetCut(cutName, cutValue, -1);
1031     return kTRUE;
1032 }
1033
1034
1035 void TFluka::SetUserScoring(const char* option, Int_t npr, char* outfile, Float_t* what)
1036 {
1037 //
1038 // Adds a user scoring option to the list
1039 //
1040     TFlukaScoringOption* opt = new TFlukaScoringOption(option, "User Scoring", npr,outfile,what);
1041     fUserScore->Add(opt);
1042 }
1043 //______________________________________________________________________________
1044 void TFluka::SetUserScoring(const char* option, Int_t npr, char* outfile, Float_t* what, const char* det1, const char* det2, const char* det3)
1045 {
1046 //
1047 // Adds a user scoring option to the list
1048 //
1049     TFlukaScoringOption* opt = new TFlukaScoringOption(option, "User Scoring", npr, outfile, what, det1, det2, det3);
1050     fUserScore->Add(opt);
1051 }
1052
1053 //______________________________________________________________________________ 
1054 Double_t TFluka::Xsec(char*, Double_t, Int_t, Int_t)
1055 {
1056   Warning("Xsec", "Not yet implemented.!\n"); return -1.;
1057 }
1058
1059
1060 //______________________________________________________________________________ 
1061 void TFluka::InitPhysics()
1062 {
1063 //
1064 // Physics initialisation with preparation of FLUKA input cards
1065 //
1066 // Construct file names
1067     FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
1068     TString sFlukaVmcCoreInp = getenv("ALICE_ROOT");
1069     sFlukaVmcCoreInp +="/TFluka/input/";
1070     TString sFlukaVmcTmp = "flukaMat.inp";
1071     TString sFlukaVmcInp = GetInputFileName();
1072     sFlukaVmcCoreInp += GetCoreInputFileName();
1073     
1074 // Open files 
1075     if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
1076         Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
1077         exit(1);
1078     }
1079     if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
1080         Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcTmp.Data());
1081         exit(1);
1082     }
1083     if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
1084         Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcInp.Data());
1085         exit(1);
1086     }
1087
1088 // Copy core input file 
1089     Char_t sLine[255];
1090     Float_t fEventsPerRun;
1091     
1092     while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
1093         if (strncmp(sLine,"GEOEND",6) != 0)
1094             fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
1095         else {
1096             fprintf(pFlukaVmcInp,"GEOEND\n");   // add GEOEND card
1097             goto flukamat;
1098         }
1099     } // end of while until GEOEND card
1100     
1101
1102  flukamat:
1103     while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
1104         fprintf(pFlukaVmcInp,"%s\n",sLine);
1105     }
1106     
1107     while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) { 
1108         if (strncmp(sLine,"START",5) != 0)
1109             fprintf(pFlukaVmcInp,"%s\n",sLine);
1110         else {
1111             sscanf(sLine+10,"%10f",&fEventsPerRun);
1112             goto fin;
1113         }
1114     } //end of while until START card
1115     
1116  fin:
1117
1118     
1119 // Pass information to configuration objects
1120     
1121     Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
1122     TFlukaConfigOption::SetStaticInfo(pFlukaVmcInp, 3, fLastMaterial, fGeom);
1123     
1124     TIter next(fUserConfig);
1125     TFlukaConfigOption* proc;
1126     while((proc = dynamic_cast<TFlukaConfigOption*> (next()))) proc->WriteFlukaInputCards();
1127 //
1128 // Process Fluka specific scoring options
1129 //
1130     TFlukaScoringOption::SetStaticInfo(pFlukaVmcInp, fGeom);
1131     Float_t loginp        = 49.0;
1132     Int_t inp             = 0;
1133     Int_t nscore          = fUserScore->GetEntries();
1134     
1135     TFlukaScoringOption *mopo = 0;
1136     TFlukaScoringOption *mopi = 0;
1137
1138     for (Int_t isc = 0; isc < nscore; isc++) 
1139     {
1140         mopo = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isc));
1141         char*    fileName = mopo->GetFileName();
1142         Int_t    size     = strlen(fileName);
1143         Float_t  lun      = -1.;
1144 //
1145 // Check if new output file has to be opened
1146         for (Int_t isci = 0; isci < isc; isci++) {
1147
1148             
1149             mopi = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isci));
1150             if(strncmp(mopi->GetFileName(), fileName, size)==0) {
1151                 // 
1152                 // No, the file already exists
1153                 lun = mopi->GetLun();
1154                 mopo->SetLun(lun);
1155                 break;
1156             }
1157         } // inner loop
1158
1159         if (lun == -1.) {
1160             // Open new output file
1161             inp++;
1162             mopo->SetLun(loginp + inp);
1163             mopo->WriteOpenFlukaFile();
1164         }
1165         mopo->WriteFlukaInputCards();
1166     }
1167
1168 // Add RANDOMIZ card
1169     fprintf(pFlukaVmcInp,"RANDOMIZ  %10.1f%10.0f\n", 1., Float_t(gRandom->GetSeed()));
1170 // Add START and STOP card
1171     fprintf(pFlukaVmcInp,"START     %10.1f\n",fEventsPerRun);
1172     fprintf(pFlukaVmcInp,"STOP      \n");
1173    
1174   
1175 // Close files
1176    fclose(pFlukaVmcCoreInp);
1177    fclose(pFlukaVmcFlukaMat);
1178    fclose(pFlukaVmcInp);
1179
1180
1181 //
1182 // Initialisation needed for Cerenkov photon production and transport
1183     TObjArray *matList = GetFlukaMaterials();
1184     Int_t nmaterial =  matList->GetEntriesFast();
1185     fMaterials = new Int_t[nmaterial+3];
1186     
1187     for (Int_t im = 0; im < nmaterial; im++)
1188     {
1189         TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
1190         Int_t idmat = material->GetIndex();
1191         fMaterials[idmat] = im;
1192     }
1193 } // end of InitPhysics
1194
1195
1196 //______________________________________________________________________________ 
1197 void TFluka::SetMaxStep(Double_t step)
1198 {
1199 // Set the maximum step size
1200     if (step > 1.e4) return;
1201     
1202     Int_t mreg, latt;
1203     fGeom->GetCurrentRegion(mreg, latt);
1204     STEPSZ.stepmx[mreg - 1] = step;
1205 }
1206
1207
1208 Double_t TFluka::MaxStep() const
1209 {
1210 // Return the maximum for current medium
1211     Int_t mreg, latt;
1212     fGeom->GetCurrentRegion(mreg, latt);
1213     return (STEPSZ.stepmx[mreg - 1]);
1214 }
1215
1216 //______________________________________________________________________________ 
1217 void TFluka::SetMaxNStep(Int_t)
1218 {
1219 // SetMaxNStep is dummy procedure in TFluka !
1220   if (fVerbosityLevel >=3)
1221   cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
1222 }
1223
1224 //______________________________________________________________________________ 
1225 void TFluka::SetUserDecay(Int_t)
1226 {
1227 // SetUserDecay is dummy procedure in TFluka !
1228   if (fVerbosityLevel >=3)
1229   cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
1230 }
1231
1232 //
1233 // dynamic properties
1234 //
1235 //______________________________________________________________________________ 
1236 void TFluka::TrackPosition(TLorentzVector& position) const
1237 {
1238 // Return the current position in the master reference frame of the
1239 // track being transported
1240 // TRACKR.atrack = age of the particle
1241 // TRACKR.xtrack = x-position of the last point
1242 // TRACKR.ytrack = y-position of the last point
1243 // TRACKR.ztrack = z-position of the last point
1244   FlukaCallerCode_t caller = GetCaller();
1245   if (caller == kENDRAW    || caller == kUSDRAW || 
1246       caller == kBXExiting || caller == kBXEntering || 
1247       caller == kUSTCKV) { 
1248     position.SetX(GetXsco());
1249     position.SetY(GetYsco());
1250     position.SetZ(GetZsco());
1251     position.SetT(TRACKR.atrack);
1252   }
1253   else if (caller == kMGDRAW) { 
1254     position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1255     position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1256     position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1257     position.SetT(TRACKR.atrack);
1258   }
1259   else if (caller == kSODRAW) { 
1260     position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1261     position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1262     position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1263     position.SetT(0);
1264   } else if (caller == kMGResumedTrack) { 
1265     position.SetX(TRACKR.spausr[0]);
1266     position.SetY(TRACKR.spausr[1]);
1267     position.SetZ(TRACKR.spausr[2]);
1268     position.SetT(TRACKR.spausr[3]);
1269   }
1270   else
1271     Warning("TrackPosition","position not available");
1272 }
1273
1274 //______________________________________________________________________________ 
1275 void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
1276 {
1277 // Return the current position in the master reference frame of the
1278 // track being transported
1279 // TRACKR.atrack = age of the particle
1280 // TRACKR.xtrack = x-position of the last point
1281 // TRACKR.ytrack = y-position of the last point
1282 // TRACKR.ztrack = z-position of the last point
1283   FlukaCallerCode_t caller = GetCaller();
1284   if (caller == kENDRAW    || caller == kUSDRAW || 
1285       caller == kBXExiting || caller == kBXEntering || 
1286       caller == kUSTCKV) { 
1287     x = GetXsco();
1288     y = GetYsco();
1289     z = GetZsco();
1290   }
1291   else if (caller == kMGDRAW || caller == kSODRAW) { 
1292     x = TRACKR.xtrack[TRACKR.ntrack];
1293     y = TRACKR.ytrack[TRACKR.ntrack];
1294     z = TRACKR.ztrack[TRACKR.ntrack];
1295   }
1296   else if (caller == kMGResumedTrack) {
1297     x = TRACKR.spausr[0];
1298     y = TRACKR.spausr[1];
1299     z = TRACKR.spausr[2];
1300   }
1301   else
1302     Warning("TrackPosition","position not available");
1303 }
1304
1305 //______________________________________________________________________________ 
1306 void TFluka::TrackMomentum(TLorentzVector& momentum) const
1307 {
1308 // Return the direction and the momentum (GeV/c) of the track
1309 // currently being transported
1310 // TRACKR.ptrack = momentum of the particle (not always defined, if
1311 //               < 0 must be obtained from etrack) 
1312 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1313 // TRACKR.etrack = total energy of the particle
1314 // TRACKR.jtrack = identity number of the particle
1315 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1316   FlukaCallerCode_t  caller = GetCaller();
1317   FlukaProcessCode_t icode  = GetIcode();
1318   
1319   if (caller != kEEDRAW && caller != kMGResumedTrack && 
1320       (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
1321     if (TRACKR.ptrack >= 0) {
1322       momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
1323       momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
1324       momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
1325       momentum.SetE(TRACKR.etrack);
1326       return;
1327     }
1328     else {
1329       Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
1330       momentum.SetPx(p*TRACKR.cxtrck);
1331       momentum.SetPy(p*TRACKR.cytrck);
1332       momentum.SetPz(p*TRACKR.cztrck);
1333       momentum.SetE(TRACKR.etrack);
1334       return;
1335     }
1336   } else if  (caller == kMGResumedTrack) {
1337     momentum.SetPx(TRACKR.spausr[4]);
1338     momentum.SetPy(TRACKR.spausr[5]);
1339     momentum.SetPz(TRACKR.spausr[6]);
1340     momentum.SetE (TRACKR.spausr[7]);
1341     return;
1342   } else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
1343       momentum.SetPx(0.);
1344       momentum.SetPy(0.);
1345       momentum.SetPz(0.);
1346       momentum.SetE(TrackMass());
1347   }
1348   else
1349     Warning("TrackMomentum","momentum not available");
1350 }
1351
1352 //______________________________________________________________________________ 
1353 void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e) const
1354 {
1355 // Return the direction and the momentum (GeV/c) of the track
1356 // currently being transported
1357 // TRACKR.ptrack = momentum of the particle (not always defined, if
1358 //               < 0 must be obtained from etrack) 
1359 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1360 // TRACKR.etrack = total energy of the particle
1361 // TRACKR.jtrack = identity number of the particle
1362 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1363   FlukaCallerCode_t   caller = GetCaller();
1364   FlukaProcessCode_t  icode  = GetIcode();
1365   if (caller != kEEDRAW && caller != kMGResumedTrack && 
1366       (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
1367     if (TRACKR.ptrack >= 0) {
1368       px = TRACKR.ptrack*TRACKR.cxtrck;
1369       py = TRACKR.ptrack*TRACKR.cytrck;
1370       pz = TRACKR.ptrack*TRACKR.cztrck;
1371       e  = TRACKR.etrack;
1372       return;
1373     }
1374     else {
1375       Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) *  ParticleMassFPC(TRACKR.jtrack));
1376       px = p*TRACKR.cxtrck;
1377       py = p*TRACKR.cytrck;
1378       pz = p*TRACKR.cztrck;
1379       e  = TRACKR.etrack;
1380       return;
1381     }
1382   } else if (caller == kMGResumedTrack) {
1383       px = TRACKR.spausr[4];
1384       py = TRACKR.spausr[5];
1385       pz = TRACKR.spausr[6];
1386       e  = TRACKR.spausr[7];
1387       return;
1388   } else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
1389       px = 0.;
1390       py = 0.;
1391       pz = 0.;
1392       e  = TrackMass();
1393   }
1394   else
1395     Warning("TrackMomentum","momentum not available");
1396 }
1397
1398 //______________________________________________________________________________ 
1399 Double_t TFluka::TrackStep() const
1400 {
1401 // Return the length in centimeters of the current step
1402 // TRACKR.ctrack = total curved path
1403   FlukaCallerCode_t caller = GetCaller();
1404   if (caller == kBXEntering || caller == kBXExiting || 
1405       caller == kENDRAW     || caller == kUSDRAW || 
1406       caller == kUSTCKV     || caller == kMGResumedTrack)
1407     return 0.0;
1408   else if (caller == kMGDRAW)
1409     return TRACKR.ctrack;
1410   else {
1411     Warning("TrackStep", "track step not available");
1412     return 0.0;
1413   }  
1414 }
1415
1416 //______________________________________________________________________________ 
1417 Double_t TFluka::TrackLength() const
1418 {
1419 // TRACKR.cmtrck = cumulative curved path since particle birth
1420   FlukaCallerCode_t caller = GetCaller();
1421   if (caller == kBXEntering || caller == kBXExiting || 
1422       caller == kENDRAW || caller == kUSDRAW || caller == kMGDRAW || 
1423       caller == kUSTCKV) 
1424     return TRACKR.cmtrck;
1425   else if (caller == kMGResumedTrack) 
1426     return TRACKR.spausr[8];
1427   else {
1428     Warning("TrackLength", "track length not available");
1429     return 0.0;
1430   } 
1431 }
1432
1433 //______________________________________________________________________________ 
1434 Double_t TFluka::TrackTime() const
1435 {
1436 // Return the current time of flight of the track being transported
1437 // TRACKR.atrack = age of the particle
1438   FlukaCallerCode_t caller = GetCaller();
1439   if (caller == kBXEntering || caller == kBXExiting || 
1440       caller == kENDRAW     || caller == kUSDRAW    || caller == kMGDRAW || 
1441       caller == kUSTCKV)
1442     return TRACKR.atrack;
1443   else if (caller == kMGResumedTrack)
1444     return TRACKR.spausr[3];
1445   else {
1446     Warning("TrackTime", "track time not available");
1447     return 0.0;
1448   }   
1449 }
1450
1451 //______________________________________________________________________________ 
1452 Double_t TFluka::Edep() const
1453 {
1454 // Energy deposition
1455 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
1456 // -->local energy deposition (the value and the point are not recorded in TRACKR)
1457 //    but in the variable "rull" of the procedure "endraw.cxx"
1458 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
1459 // -->no energy loss along the track
1460 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
1461 // -->energy loss distributed along the track
1462 // TRACKR.dtrack = energy deposition of the jth deposition event
1463
1464   // If coming from bxdraw we have 2 steps of 0 length and 0 edep
1465   // If coming from usdraw we just signal particle production - no edep
1466   // If just first time after resuming, no edep for the primary
1467   FlukaCallerCode_t caller = GetCaller();
1468   if (caller == kBXExiting || caller == kBXEntering || 
1469       caller == kUSDRAW    || caller == kMGResumedTrack) return 0.0;
1470   Double_t sum = 0;
1471   for ( Int_t j=0;j<TRACKR.mtrack;j++) {
1472       sum +=TRACKR.dtrack[j];  
1473   }
1474   if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
1475       return fRull + sum;
1476   else {
1477       return sum;
1478   }
1479 }
1480
1481 //______________________________________________________________________________ 
1482 Int_t TFluka::TrackPid() const
1483 {
1484 // Return the id of the particle transported
1485 // TRACKR.jtrack = identity number of the particle
1486   FlukaCallerCode_t caller = GetCaller();
1487   if (caller != kEEDRAW) {
1488       return PDGFromId(TRACKR.jtrack);
1489   }
1490   else
1491     return -1000;
1492 }
1493
1494 //______________________________________________________________________________ 
1495 Double_t TFluka::TrackCharge() const
1496 {
1497 // Return charge of the track currently transported
1498 // PAPROP.ichrge = electric charge of the particle
1499 // TRACKR.jtrack = identity number of the particle
1500   FlukaCallerCode_t caller = GetCaller();
1501   if (caller != kEEDRAW) 
1502     return PAPROP.ichrge[TRACKR.jtrack+6];
1503   else
1504     return -1000.0;
1505 }
1506
1507 //______________________________________________________________________________ 
1508 Double_t TFluka::TrackMass() const
1509 {
1510 // PAPROP.am = particle mass in GeV
1511 // TRACKR.jtrack = identity number of the particle
1512   FlukaCallerCode_t caller = GetCaller();
1513   if (caller != kEEDRAW)
1514     return PAPROP.am[TRACKR.jtrack+6];
1515   else
1516     return -1000.0;
1517 }
1518
1519 //______________________________________________________________________________ 
1520 Double_t TFluka::Etot() const
1521 {
1522 // TRACKR.etrack = total energy of the particle
1523   FlukaCallerCode_t caller = GetCaller();
1524   if (caller != kEEDRAW)
1525     return TRACKR.etrack;
1526   else
1527     return -1000.0;
1528 }
1529
1530 //
1531 // track status
1532 //
1533 //______________________________________________________________________________ 
1534 Bool_t   TFluka::IsNewTrack() const
1535 {
1536 // Return true for the first call of Stepping()
1537    return fTrackIsNew;
1538 }
1539
1540 void     TFluka::SetTrackIsNew(Bool_t flag)
1541 {
1542 // Return true for the first call of Stepping()
1543    fTrackIsNew = flag;
1544
1545 }
1546
1547
1548 //______________________________________________________________________________ 
1549 Bool_t   TFluka::IsTrackInside() const
1550 {
1551 // True if the track is not at the boundary of the current volume
1552 // In Fluka a step is always inside one kind of material
1553 // If the step would go behind the region of one material,
1554 // it will be shortened to reach only the boundary.
1555 // Therefore IsTrackInside() is always true.
1556   FlukaCallerCode_t caller = GetCaller();
1557   if (caller == kBXEntering || caller == kBXExiting)
1558     return 0;
1559   else
1560     return 1;
1561 }
1562
1563 //______________________________________________________________________________ 
1564 Bool_t   TFluka::IsTrackEntering() const
1565 {
1566 // True if this is the first step of the track in the current volume
1567
1568   FlukaCallerCode_t caller = GetCaller();
1569   if (caller == kBXEntering)
1570     return 1;
1571   else return 0;
1572 }
1573
1574 //______________________________________________________________________________ 
1575 Bool_t   TFluka::IsTrackExiting() const
1576 {
1577 // True if track is exiting volume
1578 //
1579   FlukaCallerCode_t caller = GetCaller();
1580   if (caller == kBXExiting)
1581     return 1;
1582   else return 0;
1583 }
1584
1585 //______________________________________________________________________________ 
1586 Bool_t   TFluka::IsTrackOut() const
1587 {
1588 // True if the track is out of the setup
1589 // means escape
1590   FlukaProcessCode_t icode = GetIcode();
1591     
1592   if (icode == kKASKADescape ||
1593       icode == kEMFSCOescape ||
1594       icode == kKASNEUescape ||
1595       icode == kKASHEAescape ||
1596       icode == kKASOPHescape) 
1597        return 1;
1598   else return 0;
1599 }
1600
1601 //______________________________________________________________________________ 
1602 Bool_t   TFluka::IsTrackDisappeared() const
1603 {
1604 // All inelastic interactions and decays
1605 // fIcode from usdraw
1606   FlukaProcessCode_t icode = GetIcode();
1607   if (icode == kKASKADinelint    || // inelastic interaction
1608       icode == kKASKADdecay      || // particle decay
1609       icode == kKASKADdray       || // delta ray generation by hadron
1610       icode == kKASKADpair       || // direct pair production
1611       icode == kKASKADbrems      || // bremsstrahlung (muon)
1612       icode == kEMFSCObrems      || // bremsstrahlung (electron)
1613       icode == kEMFSCOmoller     || // Moller scattering
1614       icode == kEMFSCObhabha     || // Bhaba scattering
1615       icode == kEMFSCOanniflight || // in-flight annihilation
1616       icode == kEMFSCOannirest   || // annihilation at rest
1617       icode == kEMFSCOpair       || // pair production
1618       icode == kEMFSCOcompton    || // Compton scattering
1619       icode == kEMFSCOphotoel    || // Photoelectric effect
1620       icode == kKASNEUhadronic   || // hadronic interaction
1621       icode == kKASHEAdray          // delta-ray
1622       ) return 1;
1623   else return 0;
1624 }
1625
1626 //______________________________________________________________________________ 
1627 Bool_t   TFluka::IsTrackStop() const
1628 {
1629 // True if the track energy has fallen below the threshold
1630 // means stopped by signal or below energy threshold
1631   FlukaProcessCode_t icode = GetIcode();
1632   if (icode == kKASKADstopping  ||
1633       icode == kKASKADtimekill  ||
1634       icode == kEMFSCOstopping1 ||
1635       icode == kEMFSCOstopping2 ||
1636       icode == kEMFSCOtimekill  ||
1637       icode == kKASNEUstopping  ||
1638       icode == kKASNEUtimekill  ||
1639       icode == kKASHEAtimekill  ||
1640       icode == kKASOPHtimekill) return 1;
1641   else return 0;
1642 }
1643
1644 //______________________________________________________________________________ 
1645 Bool_t   TFluka::IsTrackAlive() const
1646 {
1647 // means not disappeared or not out
1648   if (IsTrackDisappeared() || IsTrackOut() ) return 0;
1649   else return 1;
1650 }
1651
1652 //
1653 // secondaries
1654 //
1655
1656 //______________________________________________________________________________ 
1657 Int_t TFluka::NSecondaries() const
1658
1659 {
1660 // Number of secondary particles generated in the current step
1661 // GENSTK.np = number of secondaries except light and heavy ions
1662 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
1663     FlukaCallerCode_t caller = GetCaller();
1664     if (caller == kUSDRAW)  // valid only after usdraw
1665         return GENSTK.np + FHEAVY.npheav;
1666     else if (caller == kUSTCKV) {
1667         // Cerenkov Photon production
1668         return fNCerenkov;
1669     }
1670     return 0;
1671 } // end of NSecondaries
1672
1673 //______________________________________________________________________________ 
1674 void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
1675                 TLorentzVector& position, TLorentzVector& momentum)
1676 {
1677 // Copy particles from secondary stack to vmc stack
1678 //
1679
1680     FlukaCallerCode_t caller = GetCaller();
1681     if (caller == kUSDRAW) {  // valid only after usdraw
1682         if (GENSTK.np > 0) {
1683             // Hadronic interaction
1684             if (isec >= 0 && isec < GENSTK.np) {
1685                 particleId = PDGFromId(GENSTK.kpart[isec]);
1686                 position.SetX(fXsco);
1687                 position.SetY(fYsco);
1688                 position.SetZ(fZsco);
1689                 position.SetT(TRACKR.atrack);
1690                 momentum.SetPx(GENSTK.plr[isec]*GENSTK.cxr[isec]);
1691                 momentum.SetPy(GENSTK.plr[isec]*GENSTK.cyr[isec]);
1692                 momentum.SetPz(GENSTK.plr[isec]*GENSTK.czr[isec]);
1693                 momentum.SetE(GENSTK.tki[isec] + PAPROP.am[GENSTK.kpart[isec]+6]);
1694             }
1695             else if (isec >= GENSTK.np && isec < GENSTK.np + FHEAVY.npheav) {
1696                 Int_t jsec = isec - GENSTK.np;
1697                 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
1698                 position.SetX(fXsco);
1699                 position.SetY(fYsco);
1700                 position.SetZ(fZsco);
1701                 position.SetT(TRACKR.atrack);
1702                 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
1703                 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
1704                 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
1705                 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6) 
1706                     momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
1707                 else if (FHEAVY.tkheav[jsec] > 6)
1708                     momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
1709             }
1710             else
1711                 Warning("GetSecondary","isec out of range");
1712         } 
1713     } else if (caller == kUSTCKV) {
1714         Int_t index = OPPHST.lstopp - isec;
1715         position.SetX(OPPHST.xoptph[index]);
1716         position.SetY(OPPHST.yoptph[index]);
1717         position.SetZ(OPPHST.zoptph[index]);
1718         position.SetT(OPPHST.agopph[index]);
1719         Double_t p = OPPHST.poptph[index];
1720         
1721         momentum.SetPx(p * OPPHST.txopph[index]);
1722         momentum.SetPy(p * OPPHST.tyopph[index]);
1723         momentum.SetPz(p * OPPHST.tzopph[index]);
1724         momentum.SetE(p);
1725     }
1726     else
1727         Warning("GetSecondary","no secondaries available");
1728     
1729 } // end of GetSecondary
1730
1731
1732 //______________________________________________________________________________ 
1733 TMCProcess TFluka::ProdProcess(Int_t) const
1734
1735 {
1736 // Name of the process that has produced the secondary particles
1737 // in the current step
1738
1739     Int_t mugamma = (TRACKR.jtrack == kFLUKAphoton || 
1740                      TRACKR.jtrack == kFLUKAmuplus || 
1741                      TRACKR.jtrack == kFLUKAmuminus);
1742     FlukaProcessCode_t icode = GetIcode();
1743
1744     if      (icode == kKASKADdecay)                                   return kPDecay;
1745     else if (icode == kKASKADpair || icode == kEMFSCOpair)            return kPPair;
1746     else if (icode == kEMFSCOcompton)                                 return kPCompton;
1747     else if (icode == kEMFSCOphotoel)                                 return kPPhotoelectric;
1748     else if (icode == kKASKADbrems      || icode == kEMFSCObrems)     return kPBrem;
1749     else if (icode == kKASKADdray       || icode == kKASHEAdray)      return kPDeltaRay;
1750     else if (icode == kEMFSCOmoller     || icode == kEMFSCObhabha)    return kPDeltaRay;
1751     else if (icode == kEMFSCOanniflight || icode == kEMFSCOannirest)  return kPAnnihilation;
1752     else if (icode == kKASKADinelint) {
1753         if (!mugamma)                                                 return kPHadronic;
1754         else if (TRACKR.jtrack == kFLUKAphoton)                       return kPPhotoFission;
1755         else                                                          return kPMuonNuclear;
1756     }
1757     else if (icode == kEMFSCOrayleigh)                                return kPRayleigh;
1758 // Fluka codes 100, 300 and 400 still to be investigasted
1759     else                                                              return kPNoProcess;
1760 }
1761
1762
1763 Int_t TFluka::StepProcesses(TArrayI &proc) const
1764 {
1765   //
1766   // Return processes active in the current step
1767   //
1768     FlukaProcessCode_t icode = GetIcode();
1769     proc.Set(1);
1770     TMCProcess iproc;
1771     switch (icode) {
1772     case kKASKADtimekill:
1773     case kEMFSCOtimekill:
1774     case kKASNEUtimekill:
1775     case kKASHEAtimekill:
1776     case kKASOPHtimekill:
1777         iproc =  kPTOFlimit;
1778         break;
1779     case kKASKADstopping:
1780     case kKASKADescape:
1781     case kEMFSCOstopping1:
1782     case kEMFSCOstopping2:
1783     case kEMFSCOescape:
1784     case kKASNEUstopping:
1785     case kKASNEUescape:
1786     case kKASHEAescape:
1787     case kKASOPHescape:
1788         iproc = kPStop;
1789         break;
1790     case kKASOPHabsorption:
1791         iproc = kPLightAbsorption;
1792         break;
1793     case kKASOPHrefraction:
1794         iproc = kPLightRefraction;
1795     case kEMSCOlocaledep : 
1796         iproc = kPPhotoelectric;
1797         break;
1798     default:
1799         iproc = ProdProcess(0);
1800     }
1801     proc[0] = iproc;
1802     return 1;
1803 }
1804 //______________________________________________________________________________ 
1805 Int_t TFluka::VolId2Mate(Int_t id) const
1806 {
1807 //
1808 // Returns the material number for a given volume ID
1809 //
1810    return fMCGeo->VolId2Mate(id);
1811 }
1812
1813 //______________________________________________________________________________ 
1814 const char* TFluka::VolName(Int_t id) const
1815 {
1816 //
1817 // Returns the volume name for a given volume ID
1818 //
1819    return fMCGeo->VolName(id);
1820 }
1821
1822 //______________________________________________________________________________ 
1823 Int_t TFluka::VolId(const Text_t* volName) const
1824 {
1825 //
1826 // Converts from volume name to volume ID.
1827 // Time consuming. (Only used during set-up)
1828 // Could be replaced by hash-table
1829 //
1830     char sname[20];
1831     Int_t len;
1832     strncpy(sname, volName, len = strlen(volName));
1833     sname[len] = 0;
1834     while (sname[len - 1] == ' ') sname[--len] = 0;
1835     return fMCGeo->VolId(sname);
1836 }
1837
1838 //______________________________________________________________________________ 
1839 Int_t TFluka::CurrentVolID(Int_t& copyNo) const
1840 {
1841 //
1842 // Return the logical id and copy number corresponding to the current fluka region
1843 //
1844   if (gGeoManager->IsOutside()) return 0;
1845   TGeoNode *node = gGeoManager->GetCurrentNode();
1846   copyNo = node->GetNumber();
1847   Int_t id = node->GetVolume()->GetNumber();
1848   return id;
1849
1850
1851 //______________________________________________________________________________ 
1852 Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
1853 {
1854 //
1855 // Return the logical id and copy number of off'th mother 
1856 // corresponding to the current fluka region
1857 //
1858   if (off<0 || off>gGeoManager->GetLevel()) return 0;
1859   if (off==0) return CurrentVolID(copyNo);
1860   TGeoNode *node = gGeoManager->GetMother(off);
1861   if (!node) return 0;
1862   copyNo = node->GetNumber();
1863   return node->GetVolume()->GetNumber();
1864 }
1865
1866 //______________________________________________________________________________ 
1867 const char* TFluka::CurrentVolName() const
1868 {
1869 //
1870 // Return the current volume name
1871 //
1872   if (gGeoManager->IsOutside()) return 0;
1873   return gGeoManager->GetCurrentVolume()->GetName();
1874 }
1875
1876 //______________________________________________________________________________ 
1877 const char* TFluka::CurrentVolOffName(Int_t off) const
1878 {
1879 //
1880 // Return the volume name of the off'th mother of the current volume
1881 //
1882   if (off<0 || off>gGeoManager->GetLevel()) return 0;
1883   if (off==0) return CurrentVolName();
1884   TGeoNode *node = gGeoManager->GetMother(off);
1885   if (!node) return 0;
1886   return node->GetVolume()->GetName();
1887 }
1888
1889 const char* TFluka::CurrentVolPath() {
1890   // Return the current volume path
1891   return gGeoManager->GetPath(); 
1892 }
1893 //______________________________________________________________________________ 
1894 Int_t TFluka::CurrentMaterial(Float_t & a, Float_t & z, 
1895                       Float_t & dens, Float_t & radl, Float_t & absl) const
1896 {
1897 //
1898 //  Return the current medium number and material properties
1899 //
1900   Int_t copy;
1901   Int_t id  =  TFluka::CurrentVolID(copy);
1902   Int_t med =  TFluka::VolId2Mate(id);
1903   TGeoVolume*     vol = gGeoManager->GetCurrentVolume();
1904   TGeoMaterial*   mat = vol->GetMaterial();
1905   a    = mat->GetA();
1906   z    = mat->GetZ();
1907   dens = mat->GetDensity();
1908   radl = mat->GetRadLen();
1909   absl = mat->GetIntLen();
1910   
1911   return med;
1912 }
1913
1914 //______________________________________________________________________________ 
1915 void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
1916 {
1917 // Transforms a position from the world reference frame
1918 // to the current volume reference frame.
1919 //
1920 //  Geant3 desription:
1921 //  ==================
1922 //       Computes coordinates XD (in DRS) 
1923 //       from known coordinates XM in MRS 
1924 //       The local reference system can be initialized by
1925 //         - the tracking routines and GMTOD used in GUSTEP
1926 //         - a call to GMEDIA(XM,NUMED)
1927 //         - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER) 
1928 //             (inverse routine is GDTOM) 
1929 //
1930 //        If IFLAG=1  convert coordinates 
1931 //           IFLAG=2  convert direction cosinus
1932 //
1933 // ---
1934    Double_t xmL[3], xdL[3];
1935    Int_t i;
1936    for (i=0;i<3;i++) xmL[i]=xm[i];
1937    if (iflag == 1) gGeoManager->MasterToLocal(xmL,xdL);
1938    else            gGeoManager->MasterToLocalVect(xmL,xdL);
1939    for (i=0;i<3;i++) xd[i] = xdL[i];
1940 }
1941   
1942 //______________________________________________________________________________ 
1943 void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
1944 {
1945 //
1946 // See Gmtod(Float_t*, Float_t*, Int_t)
1947 //
1948    if (iflag == 1) gGeoManager->MasterToLocal(xm,xd);
1949    else            gGeoManager->MasterToLocalVect(xm,xd);
1950 }
1951
1952 //______________________________________________________________________________ 
1953 void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
1954 {
1955 // Transforms a position from the current volume reference frame
1956 // to the world reference frame.
1957 //
1958 //  Geant3 desription:
1959 //  ==================
1960 //  Computes coordinates XM (Master Reference System
1961 //  knowing the coordinates XD (Detector Ref System)
1962 //  The local reference system can be initialized by
1963 //    - the tracking routines and GDTOM used in GUSTEP
1964 //    - a call to GSCMED(NLEVEL,NAMES,NUMBER)
1965 //        (inverse routine is GMTOD)
1966 // 
1967 //   If IFLAG=1  convert coordinates
1968 //      IFLAG=2  convert direction cosinus
1969 //
1970 // ---
1971    Double_t xmL[3], xdL[3];
1972    Int_t i;
1973    for (i=0;i<3;i++) xdL[i] = xd[i];
1974    if (iflag == 1) gGeoManager->LocalToMaster(xdL,xmL);
1975    else            gGeoManager->LocalToMasterVect(xdL,xmL);
1976    for (i=0;i<3;i++) xm[i]=xmL[i];
1977 }
1978
1979 //______________________________________________________________________________ 
1980 void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
1981 {
1982 //
1983 // See Gdtom(Float_t*, Float_t*, Int_t)
1984 //
1985    if (iflag == 1) gGeoManager->LocalToMaster(xd,xm);
1986    else            gGeoManager->LocalToMasterVect(xd,xm);
1987 }
1988
1989 //______________________________________________________________________________
1990 TObjArray *TFluka::GetFlukaMaterials()
1991 {
1992 //
1993 // Get array of Fluka materials
1994    return fGeom->GetMatList();
1995 }   
1996
1997 //______________________________________________________________________________
1998 void TFluka::SetMreg(Int_t l, Int_t lttc) 
1999 {
2000 // Set current fluka region
2001    fCurrentFlukaRegion = l;
2002    fGeom->SetMreg(l,lttc);
2003 }
2004
2005
2006
2007
2008 TString TFluka::ParticleName(Int_t pdg) const
2009 {
2010     // Return particle name for particle with pdg code pdg.
2011     Int_t ifluka = IdFromPDG(pdg);
2012     return TString((CHPPRP.btype[ifluka - kFLUKAcodemin]), 8);
2013 }
2014  
2015
2016 Double_t TFluka::ParticleMass(Int_t pdg) const
2017 {
2018     // Return particle mass for particle with pdg code pdg.
2019     Int_t ifluka = IdFromPDG(pdg);
2020     return (PAPROP.am[ifluka - kFLUKAcodemin]);
2021 }
2022
2023 Double_t TFluka::ParticleMassFPC(Int_t fpc) const
2024 {
2025     // Return particle mass for particle with Fluka particle code fpc
2026     return (PAPROP.am[fpc - kFLUKAcodemin]);
2027 }
2028
2029 Double_t TFluka::ParticleCharge(Int_t pdg) const
2030 {
2031     // Return particle charge for particle with pdg code pdg.
2032     Int_t ifluka = IdFromPDG(pdg);
2033     return Double_t(PAPROP.ichrge[ifluka - kFLUKAcodemin]);
2034 }
2035
2036 Double_t TFluka::ParticleLifeTime(Int_t pdg) const
2037 {
2038     // Return particle lifetime for particle with pdg code pdg.
2039     Int_t ifluka = IdFromPDG(pdg);
2040     return (PAPROP.tmnlf[ifluka - kFLUKAcodemin]);
2041 }
2042
2043 void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t& charge, Float_t& tlife)
2044 {
2045     // Retrieve particle properties for particle with pdg code pdg.
2046     
2047     strcpy(name, ParticleName(pdg).Data());
2048     type   = ParticleMCType(pdg);
2049     mass   = ParticleMass(pdg);
2050     charge = ParticleCharge(pdg);
2051     tlife  = ParticleLifeTime(pdg);
2052 }
2053
2054 void TFluka::PrintHeader()
2055 {
2056     //
2057     // Print a header
2058     printf("\n");
2059     printf("\n");    
2060     printf("------------------------------------------------------------------------------\n");
2061     printf("- You are using the TFluka Virtual Monte Carlo Interface to FLUKA.           -\n");    
2062     printf("- Please see the file fluka.out for FLUKA output and licensing information.  -\n");    
2063     printf("------------------------------------------------------------------------------\n");
2064     printf("\n");
2065     printf("\n");    
2066 }
2067
2068
2069 #define pshckp pshckp_
2070 #define ustckv ustckv_
2071
2072
2073 extern "C" {
2074   void pshckp(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
2075               Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
2076               Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
2077   {
2078     //
2079     // Pushes one cerenkov photon to the stack
2080     //
2081     
2082     TFluka* fluka =  (TFluka*) gMC;
2083     TVirtualMCStack* cppstack = fluka->GetStack();
2084     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
2085     cppstack->PushTrack(0, parent, 50000050,
2086                         px, py, pz, e,
2087                         vx, vy, vz, tof,
2088                         polx, poly, polz,
2089                         kPCerenkov, ntr, wgt, 0);
2090   }
2091     
2092     void ustckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
2093     {
2094         //
2095         // Calls stepping in order to signal cerenkov production
2096         //
2097         TFluka *fluka = (TFluka*)gMC;
2098         fluka->SetMreg(mreg,LTCLCM.mlatm1);
2099         fluka->SetXsco(x);
2100         fluka->SetYsco(y);
2101         fluka->SetZsco(z);
2102         fluka->SetNCerenkov(nphot);
2103         fluka->SetCaller(kUSTCKV);
2104         if (fluka->GetVerbosityLevel() >= 3) 
2105         (TVirtualMCApplication::Instance())->Stepping();
2106         
2107     }
2108 }
2109