f384afd791c14c00d0360f1674d56e56f840e865
[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 #include <TList.h>
34
35 #include "TFluka.h"
36 #include "TFlukaIon.h"
37 #include "TFlukaCodes.h"
38 #include "TCallf77.h"      //For the fortran calls
39 #include "Fdblprc.h"       //(DBLPRC) fluka common
40 #include "Fsourcm.h"       //(SOURCM) fluka common
41 #include "Fgenstk.h"       //(GENSTK)  fluka common
42 #include "Fiounit.h"       //(IOUNIT) fluka common
43 #include "Fpaprop.h"       //(PAPROP) fluka common
44 #include "Fpart.h"         //(PART)   fluka common
45 #include "Ftrackr.h"       //(TRACKR) fluka common
46 #include "Fpaprop.h"       //(PAPROP) fluka common
47 #include "Ffheavy.h"       //(FHEAVY) fluka common
48 #include "Fopphst.h"       //(OPPHST) fluka common
49 #include "Fflkstk.h"       //(FLKSTK) fluka common
50 #include "Fstepsz.h"       //(STEPSZ) fluka common
51 #include "Fopphst.h"       //(OPPHST) fluka common
52 #include "Fltclcm.h"       //(LTCLCM) fluka common
53 #include "Falldlt.h"       //(ALLDLT) fluka common
54
55 #include "TVirtualMC.h"
56 #include "TMCProcess.h"
57 #include "TGeoManager.h"
58 #include "TGeoMaterial.h"
59 #include "TGeoMedium.h"
60 #include "TFlukaMCGeometry.h"
61 #include "TGeoMCGeometry.h"
62 #include "TFlukaCerenkov.h"
63 #include "TFlukaConfigOption.h"
64 #include "TFlukaScoringOption.h"
65 #include "TLorentzVector.h"
66 #include "TArrayI.h"
67 #include "TArrayD.h"
68 #include "TDatabasePDG.h"
69 #include "TStopwatch.h"
70
71
72 // Fluka methods that may be needed.
73 #ifndef WIN32 
74 # define flukam  flukam_
75 # define fluka_openinp fluka_openinp_
76 # define fluka_openout fluka_openout_
77 # define fluka_closeinp fluka_closeinp_
78 # define mcihad mcihad_
79 # define mpdgha mpdgha_
80 # define newplo newplo_
81 # define genout genout_
82 # define flkend flkend_
83 #else 
84 # define flukam  FLUKAM
85 # define fluka_openinp FLUKA_OPENINP
86 # define fluka_openout FLUKA_OPENOUT
87 # define fluka_closeinp FLUKA_CLOSEINP
88 # define mcihad MCIHAD
89 # define mpdgha MPDGHA
90 # define newplo NEWPLO
91 # define genout GENOUT
92 # define flkend FLKEND
93 #endif
94
95 extern "C" 
96 {
97   //
98   // Prototypes for FLUKA functions
99   //
100   void type_of_call flukam(const int&);
101   void type_of_call newplo();
102   void type_of_call genout();
103   void type_of_call flkend();
104   void type_of_call fluka_openinp(const int&, DEFCHARA);
105   void type_of_call fluka_openout(const int&, DEFCHARA);
106   void type_of_call fluka_closeinp(const int&);
107   int  type_of_call mcihad(const int&);
108   int  type_of_call mpdgha(const int&);
109 }
110
111 //
112 // Class implementation for ROOT
113 //
114 ClassImp(TFluka)
115
116 //
117 //----------------------------------------------------------------------------
118 // TFluka constructors and destructors.
119 //______________________________________________________________________________
120 TFluka::TFluka()
121   :TVirtualMC(),
122    fVerbosityLevel(0),
123    fNEvent(0),
124    fInputFileName(""),
125    fCoreInputFileName(""),
126    fCaller(kNoCaller),
127    fIcode(kNoProcess),
128    fNewReg(-1),
129    fRull(0),
130    fXsco(0),
131    fYsco(0),
132    fZsco(0),
133    fPItime(0),
134    fPIlength(0),
135    fNPI(0),
136    fTrackIsEntering(kFALSE),
137    fTrackIsExiting(kFALSE),
138    fTrackIsNew(kFALSE),
139    fFieldFlag(kTRUE),
140    fDummyBoundary(kFALSE),
141    fStopped(kFALSE),
142    fStopEvent(kFALSE),
143    fStopRun(kFALSE),
144    fPrimaryElectronIndex(-1),
145    fLowEnergyNeutronTransport(kFALSE),
146    fMaterials(0),
147    fNVolumes(0),
148    fCurrentFlukaRegion(-1),
149    fNCerenkov(0),
150    fGeom(0),
151    fMCGeo(0),
152    fUserConfig(0), 
153    fUserScore(0),
154    fUserIons(0)
155
156   //
157   // Default constructor
158   //
159     for (Int_t i = 0; i < 4; i++) fPint[i] = 0.;
160
161  
162 //______________________________________________________________________________ 
163 TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
164   :TVirtualMC("TFluka",title, isRootGeometrySupported),
165    fVerbosityLevel(verbosity),
166    fNEvent(0),
167    fInputFileName(""),
168    fCoreInputFileName(""),
169    fCaller(kNoCaller),
170    fIcode(kNoProcess),
171    fNewReg(-1),
172    fRull(0),
173    fXsco(0),
174    fYsco(0),
175    fZsco(0),
176    fPItime(0),
177    fPIlength(0),
178    fNPI(0),
179    fTrackIsEntering(kFALSE),
180    fTrackIsExiting(kFALSE),
181    fTrackIsNew(kFALSE),
182    fFieldFlag(kTRUE),
183    fDummyBoundary(kFALSE),
184    fStopped(kFALSE),
185    fStopEvent(kFALSE),
186    fStopRun(kFALSE),
187    fPrimaryElectronIndex(-1),
188    fLowEnergyNeutronTransport(kFALSE),
189    fMaterials(0),
190    fNVolumes(0),
191    fCurrentFlukaRegion(-1),
192    fNCerenkov(0),
193    fGeom(0),
194    fMCGeo(0),
195    fUserConfig(new TObjArray(100)),
196    fUserScore(new TObjArray(100)),
197    fUserIons(0)
198 {
199   // create geometry interface
200     for (Int_t i = 0; i < 4; i++) fPint[i] = 0.;
201     
202    if (fVerbosityLevel >=3)
203        cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
204    SetCoreInputFileName();
205    SetInputFileName();
206    fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kFALSE);
207    fGeom  = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
208    if (verbosity > 2) fGeom->SetDebugMode(kTRUE);
209    PrintHeader();
210 }
211
212 //______________________________________________________________________________ 
213 TFluka::~TFluka()
214 {
215     // Destructor
216     if (fVerbosityLevel >=3)
217         cout << "<== TFluka::~TFluka() destructor called." << endl;
218     if (fMaterials) delete [] fMaterials;
219     
220 //    delete fGeom;
221 //    delete fMCGeo;
222     
223     if (fUserConfig) {
224         fUserConfig->Delete();
225         delete fUserConfig;
226     }
227     
228     if (fUserScore) {
229         fUserScore->Delete();
230         delete fUserScore;
231     }
232 }
233
234 //
235 //______________________________________________________________________________
236 // TFluka control methods
237 //______________________________________________________________________________ 
238 void TFluka::Init() {
239 //
240 //  Geometry initialisation
241 //
242     if (fVerbosityLevel >=3) cout << "==> TFluka::Init() called." << endl;
243     
244     if (!gGeoManager) new TGeoManager("geom", "FLUKA geometry");
245     fApplication->ConstructGeometry();
246     if (!gGeoManager->IsClosed()) {
247        TGeoVolume *top = (TGeoVolume*)gGeoManager->GetListOfVolumes()->First();
248        gGeoManager->SetTopVolume(top);
249        gGeoManager->CloseGeometry("di");
250     } else {
251        TGeoNodeCache *cache = gGeoManager->GetCache();
252        if (!cache->HasIdArray()) {
253           Warning("Init", "Node ID tracking must be enabled with TFluka: enabling...\n");
254           cache->BuildIdArray();
255        }   
256     }           
257     fNVolumes = fGeom->NofVolumes();
258     fGeom->CreateFlukaMatFile("flukaMat.inp");   
259     if (fVerbosityLevel >=3) {
260        printf("== Number of volumes: %i\n ==", fNVolumes);
261        cout << "\t* InitPhysics() - Prepare input file to be called" << endl; 
262     }
263
264     fApplication->InitGeometry();
265     fApplication->ConstructOpGeometry();
266     //
267     // Add ions to PDG Data base
268     //
269      AddParticlesToPdgDataBase();
270     //
271 }
272
273
274 //______________________________________________________________________________ 
275 void TFluka::FinishGeometry() {
276 //
277 // Build-up table with region to medium correspondance
278 //
279   if (fVerbosityLevel >=3) {
280     cout << "==> TFluka::FinishGeometry() called." << endl;
281     printf("----FinishGeometry - applying misalignment if any\n");
282     cout << "<== TFluka::FinishGeometry() called." << endl;
283   }  
284   TVirtualMCApplication::Instance()->MisalignGeometry();
285
286
287 //______________________________________________________________________________ 
288 void TFluka::BuildPhysics() {
289 //
290 //  Prepare FLUKA input files and call FLUKA physics initialisation
291 //
292     
293     if (fVerbosityLevel >=3)
294         cout << "==> TFluka::BuildPhysics() called." << endl;
295
296     
297     if (fVerbosityLevel >=3) {
298         TList *medlist = gGeoManager->GetListOfMedia();
299         TIter next(medlist);
300         TGeoMedium*   med = 0x0;
301         TGeoMaterial* mat = 0x0;
302         Int_t ic = 0;
303         
304         while((med = (TGeoMedium*)next()))
305         {
306             mat = med->GetMaterial();
307             printf("Medium %5d %12s %5d %5d\n", ic, (med->GetName()), med->GetId(), mat->GetIndex());
308             ic++;
309         }
310     }
311     
312
313     // Prepare input file with the current physics settings
314     
315     InitPhysics(); 
316 //  Open fortran files    
317     const char* fname = fInputFileName;
318     fluka_openinp(lunin, PASSCHARA(fname));
319     fluka_openout(11, PASSCHARA("fluka.out"));
320 //  Read input cards    
321     cout << "==> TFluka::BuildPhysics() Read input cards." << endl;
322     TStopwatch timer;
323     timer.Start();
324     GLOBAL.lfdrtr = true;
325     flukam(1);
326     cout << "<== TFluka::BuildPhysics() Read input cards End"
327          << Form(" R:%.2fs C:%.2fs", timer.RealTime(),timer.CpuTime()) << endl;
328 //  Close input file
329     fluka_closeinp(lunin);
330 //  Finish geometry    
331     FinishGeometry();
332 }  
333
334 //______________________________________________________________________________ 
335 void TFluka::ProcessEvent() {
336 //
337 // Process one event
338 //
339     if (fStopRun) {
340         Warning("ProcessEvent", "User Run Abortion: No more events handled !\n");
341         fNEvent += 1;
342         return;
343     }
344
345     if (fVerbosityLevel >=3)
346         cout << "==> TFluka::ProcessEvent() called." << endl;
347     fApplication->GeneratePrimaries();
348     SOURCM.lsouit = true;
349     flukam(1);
350     if (fVerbosityLevel >=3)
351         cout << "<== TFluka::ProcessEvent() called." << endl;
352     //
353     // Increase event number
354     //
355     fNEvent += 1;
356 }
357
358 //______________________________________________________________________________ 
359 Bool_t TFluka::ProcessRun(Int_t nevent) {
360 //
361 // Run steering
362 //
363     
364   if (fVerbosityLevel >=3)
365     cout << "==> TFluka::ProcessRun(" << nevent << ") called." 
366          << endl;
367
368   if (fVerbosityLevel >=2) {
369     cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
370     cout << "\t* Calling flukam again..." << endl;
371   }
372
373   Int_t todo = TMath::Abs(nevent);
374   for (Int_t ev = 0; ev < todo; ev++) {
375       TStopwatch timer;
376       timer.Start();
377       fApplication->BeginEvent();
378       ProcessEvent();
379       fApplication->FinishEvent();
380       cout << "Event: "<< ev
381            << Form(" R:%.2fs C:%.2fs", timer.RealTime(),timer.CpuTime()) << endl;
382   }
383
384   if (fVerbosityLevel >=3)
385     cout << "<== TFluka::ProcessRun(" << nevent << ") called." 
386          << endl;
387   
388   // Write fluka specific scoring output
389   genout();
390   newplo();
391   flkend();
392   
393   return kTRUE;
394 }
395
396 //_____________________________________________________________________________
397 // methods for building/management of geometry
398
399 // functions from GCONS 
400 //____________________________________________________________________________ 
401 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,  
402                     Float_t &dens, Float_t &radl, Float_t &absl,
403                     Float_t* /*ubuf*/, Int_t& /*nbuf*/) {
404 //
405    TGeoMaterial *mat;
406    TIter next (gGeoManager->GetListOfMaterials());
407    while ((mat = (TGeoMaterial*)next())) {
408      if (mat->GetUniqueID() == (UInt_t)imat) break;
409    }
410    if (!mat) {
411       Error("Gfmate", "no material with index %i found", imat);
412       return;
413    }
414    sprintf(name, "%s", mat->GetName());
415    a = mat->GetA();
416    z = mat->GetZ();
417    dens = mat->GetDensity();
418    radl = mat->GetRadLen();
419    absl = mat->GetIntLen();
420
421
422 //______________________________________________________________________________ 
423 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,  
424                     Double_t &dens, Double_t &radl, Double_t &absl,
425                     Double_t* /*ubuf*/, Int_t& /*nbuf*/) {
426 //
427    TGeoMaterial *mat;
428    TIter next (gGeoManager->GetListOfMaterials());
429    while ((mat = (TGeoMaterial*)next())) {
430      if (mat->GetUniqueID() == (UInt_t)imat) break;
431    }
432    if (!mat) {
433       Error("Gfmate", "no material with index %i found", imat);
434       return;
435    }
436    sprintf(name, "%s", mat->GetName());
437    a = mat->GetA();
438    z = mat->GetZ();
439    dens = mat->GetDensity();
440    radl = mat->GetRadLen();
441    absl = mat->GetIntLen();
442
443
444 // detector composition
445 //______________________________________________________________________________ 
446 void TFluka::Material(Int_t& kmat, const char* name, Double_t a, 
447                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
448                       Float_t* buf, Int_t nwbuf) {
449 //
450    Double_t* dbuf = fGeom->CreateDoubleArray(buf, nwbuf);  
451    Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf);
452    delete [] dbuf;
453
454
455 //______________________________________________________________________________ 
456 void TFluka::Material(Int_t& kmat, const char* name, Double_t a, 
457                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
458                       Double_t* /*buf*/, Int_t /*nwbuf*/) {
459 //
460 // Define a material
461   TGeoMaterial *mat;
462   kmat = gGeoManager->GetListOfMaterials()->GetSize();
463   if ((z-Int_t(z)) > 1E-3) {
464      mat = fGeom->GetMakeWrongMaterial(z);
465      if (mat) {
466         mat->SetRadLen(radl,absl);
467         mat->SetUniqueID(kmat);
468         return;
469      }
470   }      
471   gGeoManager->Material(name, a, z, dens, kmat, radl, absl);
472
473
474 //______________________________________________________________________________ 
475 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a, 
476                      Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
477 //
478 // Define a material mixture
479 //
480   Double_t* da = fGeom->CreateDoubleArray(a, TMath::Abs(nlmat));  
481   Double_t* dz = fGeom->CreateDoubleArray(z, TMath::Abs(nlmat));  
482   Double_t* dwmat = fGeom->CreateDoubleArray(wmat, TMath::Abs(nlmat));  
483
484   Mixture(kmat, name, da, dz, dens, nlmat, dwmat);
485   for (Int_t i=0; i<nlmat; i++) {
486     a[i] = da[i]; z[i] = dz[i]; wmat[i] = dwmat[i];
487   }  
488
489   delete [] da;
490   delete [] dz;
491   delete [] dwmat;
492
493
494 //______________________________________________________________________________ 
495 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a, 
496                      Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
497 //
498   // Defines mixture OR COMPOUND IMAT as composed by 
499   // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
500   // 
501   // If NLMAT > 0 then wmat contains the proportion by
502   // weights of each basic material in the mixture. 
503   // 
504   // If nlmat < 0 then WMAT contains the number of atoms 
505   // of a given kind into the molecule of the COMPOUND
506   // In this case, WMAT in output is changed to relative
507   // weigths.
508   //
509   Int_t i,j;
510   if (nlmat < 0) {
511      nlmat = - nlmat;
512      Double_t amol = 0;
513      for (i=0;i<nlmat;i++) {
514         amol += a[i]*wmat[i];
515      }
516      for (i=0;i<nlmat;i++) {
517         wmat[i] *= a[i]/amol;
518      }
519   }
520   kmat = gGeoManager->GetListOfMaterials()->GetSize();
521   // Check if we have elements with fractional Z
522   TGeoMaterial *mat = 0;
523   TGeoMixture *mix = 0;
524   Bool_t mixnew = kFALSE;
525   for (i=0; i<nlmat; i++) {
526      if (z[i]-Int_t(z[i]) < 1E-3) continue;
527      // We have found an element with fractional Z -> loop mixtures to look for it
528      for (j=0; j<kmat; j++) {
529         mat = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->At(j);
530         if (!mat) break;
531         if (!mat->IsMixture()) continue;
532         mix = (TGeoMixture*)mat;
533         if (TMath::Abs(z[i]-mix->GetZ()) >1E-3) continue;
534         mixnew = kTRUE;
535         break;
536      }
537      if (!mixnew) Warning("Mixture","%s : cannot find component %i with fractional Z=%f\n", name, i, z[i]);
538      break;
539   }   
540   if (mixnew) {
541      Int_t nlmatnew = nlmat+mix->GetNelements()-1;
542      Double_t *anew = new Double_t[nlmatnew];
543      Double_t *znew = new Double_t[nlmatnew];
544      Double_t *wmatnew = new Double_t[nlmatnew];
545      Int_t ind=0;
546      for (j=0; j<nlmat; j++) {
547         if (j==i) continue;
548         anew[ind] = a[j];
549         znew[ind] = z[j];
550         wmatnew[ind] = wmat[j];
551         ind++;
552      }
553      for (j=0; j<mix->GetNelements(); j++) {
554         anew[ind] = mix->GetAmixt()[j];
555         znew[ind] = mix->GetZmixt()[j];
556         wmatnew[ind] = wmat[i]*mix->GetWmixt()[j];
557         ind++;
558      }
559      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
560      delete [] anew;
561      delete [] znew;
562      delete [] wmatnew;
563      return;
564   }   
565   // Now we need to compact identical elements within the mixture
566   // First check if this happens   
567   mixnew = kFALSE;  
568   for (i=0; i<nlmat-1; i++) {
569      for (j=i+1; j<nlmat; j++) {
570         if (z[i] == z[j]) {
571            mixnew = kTRUE;
572            break;
573         }
574      }   
575      if (mixnew) break;
576   }   
577   if (mixnew) {
578      Int_t nlmatnew = 0;
579      Double_t *anew = new Double_t[nlmat];
580      Double_t *znew = new Double_t[nlmat];
581      memset(znew, 0, nlmat*sizeof(Double_t));
582      Double_t *wmatnew = new Double_t[nlmat];
583      Bool_t skipi;
584      for (i=0; i<nlmat; i++) {
585         skipi = kFALSE;
586         for (j=0; j<nlmatnew; j++) {
587            if (z[i] == z[j]) {
588               wmatnew[j] += wmat[i];
589               skipi = kTRUE;
590               break;
591            }
592         }   
593         if (skipi) continue;    
594         anew[nlmatnew] = a[i];
595         znew[nlmatnew] = z[i];
596         wmatnew[nlmatnew] = wmat[i];
597         nlmatnew++;
598      }
599      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
600      delete [] anew;
601      delete [] znew;
602      delete [] wmatnew;
603      return;     
604    }
605   gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
606
607
608 //______________________________________________________________________________ 
609 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, 
610                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
611                     Double_t stemax, Double_t deemax, Double_t epsil,
612                     Double_t stmin, Float_t* ubuf, Int_t nbuf) {
613   // Define a medium
614   // 
615   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
616   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, 
617              epsil, stmin, ubuf, nbuf);
618
619
620 //______________________________________________________________________________ 
621 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, 
622                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
623                     Double_t stemax, Double_t deemax, Double_t epsil,
624                     Double_t stmin, Double_t* ubuf, Int_t nbuf) {
625   // Define a medium
626   // 
627   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
628   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, 
629              epsil, stmin, ubuf, nbuf);
630
631
632 //______________________________________________________________________________ 
633 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX, 
634                     Double_t thetaY, Double_t phiY, Double_t thetaZ,
635                     Double_t phiZ) {
636 //        
637   krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
638   fMCGeo->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ); 
639
640
641 //______________________________________________________________________________ 
642 void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
643 //
644 //
645 //
646    Bool_t process = kFALSE;
647    Bool_t modelp  = kFALSE;
648    
649    if (strncmp(param, "DCAY",  4) == 0 ||
650        strncmp(param, "PAIR",  4) == 0 ||
651        strncmp(param, "COMP",  4) == 0 ||
652        strncmp(param, "PHOT",  4) == 0 ||
653        strncmp(param, "PFIS",  4) == 0 ||
654        strncmp(param, "DRAY",  4) == 0 ||
655        strncmp(param, "ANNI",  4) == 0 ||
656        strncmp(param, "BREM",  4) == 0 ||
657        strncmp(param, "MUNU",  4) == 0 ||
658        strncmp(param, "CKOV",  4) == 0 ||
659        strncmp(param, "HADR",  4) == 0 ||
660        strncmp(param, "LOSS",  4) == 0 ||
661        strncmp(param, "MULS",  4) == 0 ||
662        strncmp(param, "RAYL",  4) == 0 ||
663        strncmp(param, "STRA",  4) == 0) 
664    {
665        process = kTRUE;
666    } 
667    
668    if (strncmp(param, "PRIMIO_N",  8) == 0 ||
669        strncmp(param, "PRIMIO_E",  8) == 0)
670    {
671        modelp = kTRUE;
672    }
673    
674    if (process) {
675        // Process switch
676        SetProcess(param, Int_t (parval), itmed);
677    } else if (modelp) {
678        // Model parameters
679        SetModelParameter(param, parval, itmed);
680    } else {
681        // Cuts
682        SetCut(param, parval, itmed);
683    }
684    
685    
686 }    
687
688 // functions from GGEOM 
689 //_____________________________________________________________________________
690 void TFluka::Gsatt(const char *name, const char *att, Int_t val)
691
692   // Set visualisation attributes for one volume
693   char vname[5];
694   fGeom->Vname(name,vname);
695   char vatt[5];
696   fGeom->Vname(att,vatt);
697   gGeoManager->SetVolumeAttribute(vname, vatt, val);
698 }
699
700 //______________________________________________________________________________ 
701 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,  
702                      Float_t *upar, Int_t np)  {
703 //
704     return fMCGeo->Gsvolu(name, shape, nmed, upar, np); 
705 }
706
707 //______________________________________________________________________________ 
708 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,  
709                      Double_t *upar, Int_t np)  {
710 //
711     return fMCGeo->Gsvolu(name, shape, nmed, upar, np); 
712 }
713  
714 //______________________________________________________________________________ 
715 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv, 
716                    Int_t iaxis) {
717 //
718     fMCGeo->Gsdvn(name, mother, ndiv, iaxis); 
719
720
721 //______________________________________________________________________________ 
722 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv, 
723                     Int_t iaxis, Double_t c0i, Int_t numed) {
724 //
725     fMCGeo->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed); 
726
727
728 //______________________________________________________________________________ 
729 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step, 
730                    Int_t iaxis, Int_t numed, Int_t ndvmx) {
731 //        
732     fMCGeo->Gsdvt(name, mother, step, iaxis, numed, ndvmx); 
733
734
735 //______________________________________________________________________________ 
736 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step, 
737                     Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
738 //
739     fMCGeo->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx); 
740
741
742 //______________________________________________________________________________ 
743 void TFluka::Gsord(const char * /*name*/, Int_t /*iax*/) {
744 //
745 // Nothing to do with TGeo
746
747
748 //______________________________________________________________________________ 
749 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,  
750                    Double_t x, Double_t y, Double_t z, Int_t irot,
751                    const char *konly) {
752 //
753   fMCGeo->Gspos(name, nr, mother, x, y, z, irot, konly); 
754
755
756 //______________________________________________________________________________ 
757 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,  
758                     Double_t x, Double_t y, Double_t z, Int_t irot,
759                     const char *konly, Float_t *upar, Int_t np)  {
760   //
761   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); 
762
763
764 //______________________________________________________________________________ 
765 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,  
766                     Double_t x, Double_t y, Double_t z, Int_t irot,
767                     const char *konly, Double_t *upar, Int_t np)  {
768   //
769   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); 
770
771
772 //______________________________________________________________________________ 
773 void TFluka::Gsbool(const char* /*onlyVolName*/, const char* /*manyVolName*/) {
774 //
775 // Nothing to do with TGeo
776 }
777
778 //______________________________________________________________________
779 Bool_t TFluka::GetTransformation(const TString &volumePath,TGeoHMatrix &mat)
780 {
781     // Returns the Transformation matrix between the volume specified
782     // by the path volumePath and the Top or mater volume. The format
783     // of the path volumePath is as follows (assuming ALIC is the Top volume)
784     // "/ALIC_1/DDIP_1/S05I_2/S05H_1/S05G_3". Here ALIC is the top most
785     // or master volume which has only 1 instance of. Of all of the daughter
786     // volumes of ALICE, DDIP volume copy #1 is indicated. Similarly for
787     // the daughter volume of DDIP is S05I copy #2 and so on.
788     // Inputs:
789     //   TString& volumePath  The volume path to the specific volume
790     //                        for which you want the matrix. Volume name
791     //                        hierarchy is separated by "/" while the
792     //                        copy number is appended using a "_".
793     // Outputs:
794     //  TGeoHMatrix &mat      A matrix with its values set to those
795     //                        appropriate to the Local to Master transformation
796     // Return:
797     //   A logical value if kFALSE then an error occurred and no change to
798     //   mat was made.
799
800    // We have to preserve the modeler state
801    return fMCGeo->GetTransformation(volumePath, mat);
802 }   
803    
804 //______________________________________________________________________
805 Bool_t TFluka::GetShape(const TString &volumePath,TString &shapeType,
806                         TArrayD &par)
807 {
808     // Returns the shape and its parameters for the volume specified
809     // by volumeName.
810     // Inputs:
811     //   TString& volumeName  The volume name
812     // Outputs:
813     //   TString &shapeType   Shape type
814     //   TArrayD &par         A TArrayD of parameters with all of the
815     //                        parameters of the specified shape.
816     // Return:
817     //   A logical indicating whether there was an error in getting this
818     //   information
819    return fMCGeo->GetShape(volumePath, shapeType, par);
820 }
821    
822 //______________________________________________________________________
823 Bool_t TFluka::GetMaterial(const TString &volumeName,
824                             TString &name,Int_t &imat,
825                             Double_t &a,Double_t &z,Double_t &dens,
826                             Double_t &radl,Double_t &inter,TArrayD &par)
827 {
828     // Returns the Material and its parameters for the volume specified
829     // by volumeName.
830     // Note, Geant3 stores and uses mixtures as an element with an effective
831     // Z and A. Consequently, if the parameter Z is not integer, then
832     // this material represents some sort of mixture.
833     // Inputs:
834     //   TString& volumeName  The volume name
835     // Outputs:
836     //   TSrting   &name       Material name
837     //   Int_t     &imat       Material index number
838     //   Double_t  &a          Average Atomic mass of material
839     //   Double_t  &z          Average Atomic number of material
840     //   Double_t  &dens       Density of material [g/cm^3]
841     //   Double_t  &radl       Average radiation length of material [cm]
842     //   Double_t  &inter      Average interaction length of material [cm]
843     //   TArrayD   &par        A TArrayD of user defined parameters.
844     // Return:
845     //   kTRUE if no errors
846    return fMCGeo->GetMaterial(volumeName,name,imat,a,z,dens,radl,inter,par);
847 }
848
849 //______________________________________________________________________
850 Bool_t TFluka::GetMedium(const TString &volumeName,TString &name,
851                          Int_t &imed,Int_t &nmat,Int_t &isvol,Int_t &ifield,
852                          Double_t &fieldm,Double_t &tmaxfd,Double_t &stemax,
853                          Double_t &deemax,Double_t &epsil, Double_t &stmin,
854                          TArrayD &par)
855 {
856     // Returns the Medium and its parameters for the volume specified
857     // by volumeName.
858     // Inputs:
859     //   TString& volumeName  The volume name.
860     // Outputs:
861     //   TString  &name       Medium name
862     //   Int_t    &nmat       Material number defined for this medium
863     //   Int_t    &imed       The medium index number
864     //   Int_t    &isvol      volume number defined for this medium
865     //   Int_t    &iflield    Magnetic field flag
866     //   Double_t &fieldm     Magnetic field strength
867     //   Double_t &tmaxfd     Maximum angle of deflection per step
868     //   Double_t &stemax     Maximum step size
869     //   Double_t &deemax     Maximum fraction of energy allowed to be lost
870     //                        to continuous process.
871     //   Double_t &epsil      Boundary crossing precision
872     //   Double_t &stmin      Minimum step size allowed
873     //   TArrayD  &par        A TArrayD of user parameters with all of the
874     //                        parameters of the specified medium.
875     // Return:
876     //   kTRUE if there where no errors
877    return fMCGeo->GetMedium(volumeName,name,imed,nmat,isvol,ifield,fieldm,tmaxfd,stemax,deemax,epsil,stmin,par);
878 }         
879
880 //______________________________________________________________________________ 
881 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
882                          Float_t* absco, Float_t* effic, Float_t* rindex) {
883 //
884 // Set Cerenkov properties for medium itmed
885 //
886 // npckov: number of sampling points
887 // ppckov: energy values
888 // absco:  absorption length
889 // effic:  quantum efficiency
890 // rindex: refraction index
891 //
892 //
893 //  
894 //  Create object holding Cerenkov properties
895 // 
896     
897     TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex);
898 //
899 //  Pass object to medium
900     TGeoMedium* medium = gGeoManager->GetMedium(itmed);
901     medium->SetCerenkovProperties(cerenkovProperties);
902 }  
903
904 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
905                          Float_t* absco, Float_t* effic, Float_t* rindex, Float_t* rfl) {
906 //
907 // Set Cerenkov properties for medium itmed
908 //
909 // npckov: number of sampling points
910 // ppckov: energy values
911 // absco:  absorption length
912 // effic:  quantum efficiency
913 // rindex: refraction index
914 // rfl:    reflectivity for boundary to medium itmed
915 //
916 //  
917 //  Create object holding Cerenkov properties
918 //  
919     TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex, rfl);
920 //
921 //  Pass object to medium
922     TGeoMedium* medium = gGeoManager->GetMedium(itmed);
923     medium->SetCerenkovProperties(cerenkovProperties);
924 }  
925
926
927 //______________________________________________________________________________ 
928 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
929                          Double_t *absco, Double_t *effic, Double_t *rindex) {
930 //
931 // Set Cerenkov properties for medium itmed
932 //
933 // npckov: number of sampling points
934 // ppckov: energy values
935 // absco:  absorption length
936 // effic:  quantum efficiency
937 // rindex: refraction index
938 //
939
940 //
941 //  Double_t version 
942   Float_t* fppckov = CreateFloatArray(ppckov, npckov);
943   Float_t* fabsco  = CreateFloatArray(absco,  npckov);
944   Float_t* feffic  = CreateFloatArray(effic,  npckov);
945   Float_t* frindex = CreateFloatArray(rindex, npckov);
946
947   SetCerenkov(itmed, npckov, fppckov, fabsco, feffic, frindex);
948
949   delete [] fppckov;
950   delete [] fabsco;
951   delete [] feffic;
952   delete [] frindex;
953 }  
954
955 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t* ppckov,
956                          Double_t* absco, Double_t* effic, Double_t* rindex, Double_t* rfl) {
957 //
958 // Set Cerenkov properties for medium itmed
959 //
960 // npckov: number of sampling points
961 // ppckov: energy values
962 // absco:  absorption length
963 // effic:  quantum efficiency
964 // rindex: refraction index
965 // rfl:    reflectivity for boundary to medium itmed
966 //
967
968 //
969 // //  Double_t version 
970   Float_t* fppckov = CreateFloatArray(ppckov, npckov);
971   Float_t* fabsco  = CreateFloatArray(absco,  npckov);
972   Float_t* feffic  = CreateFloatArray(effic,  npckov);
973   Float_t* frindex = CreateFloatArray(rindex, npckov);
974   Float_t* frfl    = CreateFloatArray(rfl,    npckov);
975
976   SetCerenkov(itmed, npckov, fppckov, fabsco, feffic, frindex, frfl);
977
978   delete [] fppckov;
979   delete [] fabsco;
980   delete [] feffic;
981   delete [] frindex;
982   delete [] frfl;
983 }
984
985 // Euclid
986 //______________________________________________________________________________ 
987 void TFluka::WriteEuclid(const char* /*fileName*/, const char* /*topVol*/, 
988                           Int_t /*number*/, Int_t /*nlevel*/) {
989 //
990 // Not with TGeo
991    Warning("WriteEuclid", "Not implemented !");
992
993
994
995
996 //_____________________________________________________________________________
997 // methods needed by the stepping
998 //____________________________________________________________________________ 
999
1000 Int_t TFluka::GetMedium() const {
1001 //
1002 //  Get the medium number for the current fluka region
1003 //
1004     if (gGeoManager->IsOutside()) {
1005         return (-1);
1006     } else {
1007         return (fGeom->GetMedium()); // this I need to check due to remapping !!!
1008     }
1009 }
1010
1011 //____________________________________________________________________________ 
1012 Int_t TFluka::GetDummyRegion() const
1013 {
1014 // Returns index of the dummy region.
1015    return fGeom->GetDummyRegion();
1016 }   
1017
1018 //____________________________________________________________________________ 
1019 Int_t TFluka::GetDummyLattice() const
1020 {
1021 // Returns index of the dummy lattice.
1022    return fGeom->GetDummyLattice();
1023 }   
1024
1025 //____________________________________________________________________________ 
1026 // particle table usage
1027 // ID <--> PDG transformations
1028 //_____________________________________________________________________________
1029 Int_t TFluka::IdFromPDG(Int_t pdg) const 
1030 {
1031
1032     //
1033     // Return Fluka code from PDG and pseudo ENDF code
1034     Int_t idSpecial[4] = {TFlukaIon::GetIonPdg(2,4), 
1035                           TFlukaIon::GetIonPdg(2,3), 
1036                           TFlukaIon::GetIonPdg(1,3), 
1037                           TFlukaIon::GetIonPdg(1,2)};    
1038     // Catch the feedback photons
1039     if (pdg == 50000051) return (kFLUKAoptical);
1040
1041     // Light ions (d,t,h3,alpha)
1042     for (Int_t i = 0; i < 4; i++) {
1043         if (pdg == idSpecial[i]) return (i + kFLUKAcodemin);
1044     }
1045     
1046     // Heavy ions
1047     if (pdg > TFlukaIon::GetIonPdg(1,1)) return (-2);
1048
1049     // MCIHAD() goes from pdg to fluka internal.
1050     Int_t intfluka = mcihad(pdg);
1051     // KPTOIP array goes from internal to official
1052     return GetFlukaKPTOIP(intfluka);
1053 }
1054
1055 //______________________________________________________________________________ 
1056 Int_t TFluka::PDGFromId(Int_t id) const 
1057 {
1058   //
1059   // Return PDG code and pseudo ENDF code from Fluka code
1060     Int_t idSpecial[6] = {TFlukaIon::GetIonPdg(2,4), // alpha 
1061                           TFlukaIon::GetIonPdg(2,3), // He3
1062                           TFlukaIon::GetIonPdg(1,3), // triton
1063                           TFlukaIon::GetIonPdg(1,2), // deuteron 
1064                           TFlukaIon::GetIonPdg(0,0), // gen. ion 
1065                           50000050};
1066   // IPTOKP array goes from official to internal
1067
1068     if (id == kFLUKAoptical) {
1069 // Cerenkov photon
1070 //        if (fVerbosityLevel >= 3)
1071 //            printf("\n PDGFromId: Cerenkov Photon \n");
1072         return  50000050;
1073     }
1074 // Error id    
1075     if (id == 0 || id < kFLUKAcodemin || id > kFLUKAcodemax) {
1076         if (fVerbosityLevel >= 3)
1077             printf("PDGFromId: Error id = 0 %5d %5d\n", id, fCaller);
1078         return -1;
1079     }
1080 // Good id    
1081     if (id > 0) {
1082         Int_t intfluka = GetFlukaIPTOKP(id);
1083         if (intfluka == 0) {
1084             if (fVerbosityLevel >= 3)
1085                 printf("PDGFromId: Error intfluka = 0: %d\n", id);
1086             return -1;
1087         } else if (intfluka < 0) {
1088             if (fVerbosityLevel >= 3)
1089                 printf("PDGFromId: Error intfluka < 0: %d\n", id);
1090             return -1;
1091         }
1092 //        if (fVerbosityLevel >= 3)
1093 //            printf("mpdgha called with %d %d \n", id, intfluka);
1094         return mpdgha(intfluka);
1095     } else {
1096         // ions and optical photons
1097         return idSpecial[id - kFLUKAcodemin];
1098     }
1099 }
1100
1101 void TFluka::StopTrack()
1102 {
1103     // Set stopping conditions
1104     // Works for photons and charged particles
1105     fStopped = kTRUE;
1106 }
1107   
1108 //_____________________________________________________________________________
1109 // methods for physics management
1110 //____________________________________________________________________________ 
1111 //
1112 // set methods
1113 //
1114
1115 void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imed)
1116 {
1117 //  Set process user flag for material imat
1118 //
1119 //    
1120 //  Update if already in the list
1121 //
1122     TIter next(fUserConfig);
1123     TFlukaConfigOption* proc;
1124     while((proc = (TFlukaConfigOption*)next()))
1125     { 
1126         if (proc->Medium() == imed) {
1127             proc->SetProcess(flagName, flagValue);
1128             return;
1129         }
1130     }
1131     proc = new TFlukaConfigOption(imed);
1132     proc->SetProcess(flagName, flagValue);
1133     fUserConfig->Add(proc);
1134 }
1135
1136 //______________________________________________________________________________ 
1137 Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue)
1138 {
1139 //  Set process user flag 
1140 //
1141 //    
1142     SetProcess(flagName, flagValue, -1);
1143     return kTRUE;  
1144 }
1145
1146 //______________________________________________________________________________ 
1147 void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
1148 {
1149 // Set user cut value for material imed
1150 //
1151     TIter next(fUserConfig);
1152     TFlukaConfigOption* proc;
1153     while((proc = (TFlukaConfigOption*)next()))
1154     { 
1155         if (proc->Medium() == imed) {
1156             proc->SetCut(cutName, cutValue);
1157             return;
1158         }
1159     }
1160
1161     proc = new TFlukaConfigOption(imed);
1162     proc->SetCut(cutName, cutValue);
1163     fUserConfig->Add(proc);
1164 }
1165
1166
1167 //______________________________________________________________________________ 
1168 void TFluka::SetModelParameter(const char* parName, Double_t parValue, Int_t imed)
1169 {
1170 // Set model parameter for material imed
1171 //
1172     TIter next(fUserConfig);
1173     TFlukaConfigOption* proc;
1174     while((proc = (TFlukaConfigOption*)next()))
1175     { 
1176         if (proc->Medium() == imed) {
1177             proc->SetModelParameter(parName, parValue);
1178             return;
1179         }
1180     }
1181
1182     proc = new TFlukaConfigOption(imed);
1183     proc->SetModelParameter(parName, parValue);
1184     fUserConfig->Add(proc);
1185 }
1186
1187 //______________________________________________________________________________ 
1188 Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue)
1189 {
1190 // Set user cut value 
1191 //
1192 //    
1193     SetCut(cutName, cutValue, -1);
1194     return kTRUE;
1195 }
1196
1197
1198 void TFluka::SetUserScoring(const char* option, const char* sdum, Int_t npr, char* outfile, Float_t* what)
1199 {
1200 //
1201 // Adds a user scoring option to the list
1202 //
1203     TFlukaScoringOption* opt = new TFlukaScoringOption(option, sdum, npr,outfile,what);
1204     fUserScore->Add(opt);
1205 }
1206 //______________________________________________________________________________
1207 void TFluka::SetUserScoring(const char* option, const char* sdum, Int_t npr, char* outfile, Float_t* what, 
1208                             const char* det1, const char* det2, const char* det3)
1209 {
1210 //
1211 // Adds a user scoring option to the list
1212 //
1213     TFlukaScoringOption* opt = new TFlukaScoringOption(option, sdum, npr, outfile, what, det1, det2, det3);
1214     fUserScore->Add(opt);
1215 }
1216
1217 //______________________________________________________________________________ 
1218 Double_t TFluka::Xsec(char*, Double_t, Int_t, Int_t)
1219 {
1220   Warning("Xsec", "Not yet implemented.!\n"); return -1.;
1221 }
1222
1223
1224 //______________________________________________________________________________ 
1225 void TFluka::InitPhysics()
1226 {
1227 //
1228 // Physics initialisation with preparation of FLUKA input cards
1229 //
1230 // Construct file names
1231     FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
1232     TString sFlukaVmcTmp = "flukaMat.inp";
1233     TString sFlukaVmcInp = GetInputFileName();
1234     TString sFlukaVmcCoreInp = GetCoreInputFileName();
1235     
1236 // Open files 
1237     if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
1238         Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
1239         exit(1);
1240     }
1241     if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
1242         Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcTmp.Data());
1243         exit(1);
1244     }
1245     if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
1246         Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcInp.Data());
1247         exit(1);
1248     }
1249
1250 // Copy core input file 
1251     Char_t sLine[255];
1252     Float_t fEventsPerRun;
1253     
1254     while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
1255         if (strncmp(sLine,"GEOEND",6) != 0)
1256             fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
1257         else {
1258             fprintf(pFlukaVmcInp,"GEOEND\n");   // add GEOEND card
1259             goto flukamat;
1260         }
1261     } // end of while until GEOEND card
1262     
1263
1264  flukamat:
1265     while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
1266         fprintf(pFlukaVmcInp,"%s\n",sLine);
1267     }
1268     
1269     while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) { 
1270         if (strncmp(sLine,"START",5) != 0)
1271             fprintf(pFlukaVmcInp,"%s\n",sLine);
1272         else {
1273             sscanf(sLine+10,"%10f",&fEventsPerRun);
1274             goto fin;
1275         }
1276     } //end of while until START card
1277     
1278  fin:
1279
1280     
1281 // Pass information to configuration objects
1282     
1283     Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
1284     TFlukaConfigOption::SetStaticInfo(pFlukaVmcInp, 3, fLastMaterial, fGeom);
1285     
1286     TIter next(fUserConfig);
1287     TFlukaConfigOption* proc;
1288     while((proc = dynamic_cast<TFlukaConfigOption*> (next()))) proc->WriteFlukaInputCards();
1289 //
1290 // Process Fluka specific scoring options
1291 //
1292     TFlukaScoringOption::SetStaticInfo(pFlukaVmcInp, fGeom);
1293     Float_t loginp        = -49.0;
1294     Int_t inp             = 0;
1295     Int_t nscore          = fUserScore->GetEntries();
1296     
1297     TFlukaScoringOption *mopo = 0;
1298     TFlukaScoringOption *mopi = 0;
1299
1300     for (Int_t isc = 0; isc < nscore; isc++) 
1301     {
1302         mopo = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isc));
1303         char*    fileName = mopo->GetFileName();
1304         Int_t    size     = strlen(fileName);
1305         Float_t  lun      = -1.;
1306 //
1307 // Check if new output file has to be opened
1308         for (Int_t isci = 0; isci < isc; isci++) {
1309
1310         
1311             mopi = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isci));
1312             if(strncmp(mopi->GetFileName(), fileName, size)==0) {
1313                 //
1314                 // No, the file already exists
1315                 lun = mopi->GetLun();
1316                 mopo->SetLun(lun);
1317                 break;
1318             }
1319         } // inner loop
1320
1321         if (lun == -1.) {
1322             // Open new output file
1323             inp++;
1324             mopo->SetLun(loginp + inp);
1325             mopo->WriteOpenFlukaFile();
1326         }
1327         mopo->WriteFlukaInputCards();
1328     }
1329
1330 // Add RANDOMIZ card
1331     fprintf(pFlukaVmcInp,"RANDOMIZ  %10.1f%10.0f\n", 1., Float_t(gRandom->GetSeed()));
1332 // User defined ion
1333     if (fUserIons) TFlukaIon::WriteUserInputCard(pFlukaVmcInp);
1334 // Add START and STOP card
1335     fprintf(pFlukaVmcInp,"START     %10.1f\n",fEventsPerRun);
1336     fprintf(pFlukaVmcInp,"STOP      \n");
1337    
1338   
1339 // Close files
1340    fclose(pFlukaVmcCoreInp);
1341    fclose(pFlukaVmcFlukaMat);
1342    fclose(pFlukaVmcInp);
1343
1344
1345 //
1346 // Initialisation needed for Cerenkov photon production and transport
1347     TObjArray *matList = GetFlukaMaterials();
1348     Int_t nmaterial =  matList->GetEntriesFast();
1349     fMaterials = new Int_t[nmaterial+25];
1350     
1351     for (Int_t im = 0; im < nmaterial; im++)
1352     {
1353         TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
1354         Int_t idmat = material->GetIndex();
1355         fMaterials[idmat] = im;
1356     }
1357 } // end of InitPhysics
1358
1359
1360 //______________________________________________________________________________ 
1361 void TFluka::SetMaxStep(Double_t step)
1362 {
1363 // Set the maximum step size
1364 //    if (step > 1.e4) return;
1365     
1366 //    Int_t mreg=0, latt=0;
1367 //    fGeom->GetCurrentRegion(mreg, latt);
1368
1369     
1370     Int_t mreg = fGeom->GetCurrentRegion();
1371     STEPSZ.stepmx[mreg - 1] = step;
1372 }
1373
1374
1375 Double_t TFluka::MaxStep() const
1376 {
1377 // Return the maximum for current medium
1378     Int_t mreg, latt;
1379     fGeom->GetCurrentRegion(mreg, latt);
1380     return (STEPSZ.stepmx[mreg - 1]);
1381 }
1382
1383 //______________________________________________________________________________ 
1384 void TFluka::SetMaxNStep(Int_t)
1385 {
1386 // SetMaxNStep is dummy procedure in TFluka !
1387   if (fVerbosityLevel >=3)
1388   cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
1389 }
1390
1391 //______________________________________________________________________________ 
1392 void TFluka::SetUserDecay(Int_t)
1393 {
1394 // SetUserDecay is dummy procedure in TFluka !
1395   if (fVerbosityLevel >=3)
1396   cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
1397 }
1398
1399 //
1400 // dynamic properties
1401 //
1402 //______________________________________________________________________________ 
1403 void TFluka::TrackPosition(TLorentzVector& position) const
1404 {
1405 // Return the current position in the master reference frame of the
1406 // track being transported
1407 // TRACKR.atrack = age of the particle
1408 // TRACKR.xtrack = x-position of the last point
1409 // TRACKR.ytrack = y-position of the last point
1410 // TRACKR.ztrack = z-position of the last point
1411   FlukaCallerCode_t caller = GetCaller();
1412   if (caller == kENDRAW    || caller == kUSDRAW || 
1413       caller == kBXExiting || caller == kBXEntering || 
1414       caller == kUSTCKV) { 
1415       position.SetX(GetXsco());
1416       position.SetY(GetYsco());
1417       position.SetZ(GetZsco());
1418       position.SetT(TRACKR.atrack);
1419   }
1420   else if (caller == kMGDRAW) {
1421       Int_t i = -1;
1422       if ((i = fPrimaryElectronIndex) > -1) {
1423           // Primary Electron Ionisation
1424           Double_t x, y, z, t;
1425           GetPrimaryElectronPosition(i, x, y, z, t);
1426           position.SetX(x);
1427           position.SetY(y);
1428           position.SetZ(z);
1429           position.SetT(t);
1430       } else {
1431           position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1432           position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1433           position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1434           position.SetT(TRACKR.atrack);
1435       }
1436   }
1437   else if (caller == kSODRAW) { 
1438       Int_t ist = FLKSTK.npflka;
1439       position.SetX(FLKSTK.xflk[ist]);
1440       position.SetY(FLKSTK.yflk[ist]);
1441       position.SetZ(FLKSTK.zflk[ist]);
1442       position.SetT(FLKSTK.agestk[ist]);
1443   } else if (caller == kMGResumedTrack) { 
1444       position.SetX(TRACKR.spausr[0]);
1445       position.SetY(TRACKR.spausr[1]);
1446       position.SetZ(TRACKR.spausr[2]);
1447       position.SetT(TRACKR.spausr[3]);
1448   }
1449   else
1450       Warning("TrackPosition","position not available");
1451 }
1452
1453 //______________________________________________________________________________ 
1454 void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
1455 {
1456 // Return the current position in the master reference frame of the
1457 // track being transported
1458 // TRACKR.atrack = age of the particle
1459 // TRACKR.xtrack = x-position of the last point
1460 // TRACKR.ytrack = y-position of the last point
1461 // TRACKR.ztrack = z-position of the last point
1462   FlukaCallerCode_t caller = GetCaller();
1463   if (caller == kENDRAW    || caller == kUSDRAW || 
1464       caller == kBXExiting || caller == kBXEntering || 
1465       caller == kUSTCKV) { 
1466       x = GetXsco();
1467       y = GetYsco();
1468       z = GetZsco();
1469   }
1470   else if (caller == kMGDRAW) { 
1471       Int_t i = -1;
1472       if ((i = fPrimaryElectronIndex) > -1) {
1473           Double_t t;
1474           GetPrimaryElectronPosition(i, x, y, z, t);
1475       } else {
1476           x = TRACKR.xtrack[TRACKR.ntrack];
1477           y = TRACKR.ytrack[TRACKR.ntrack];
1478           z = TRACKR.ztrack[TRACKR.ntrack];
1479       }
1480   }
1481   else if (caller == kSODRAW) { 
1482       Int_t ist = FLKSTK.npflka;
1483       x = FLKSTK.xflk[ist];
1484       y = FLKSTK.yflk[ist];
1485       z = FLKSTK.zflk[ist];
1486   }
1487   else if (caller == kMGResumedTrack) {
1488       x = TRACKR.spausr[0];
1489       y = TRACKR.spausr[1];
1490       z = TRACKR.spausr[2];
1491   }
1492   else
1493       Warning("TrackPosition","position not available");
1494 }
1495
1496 //______________________________________________________________________________ 
1497 void TFluka::TrackMomentum(TLorentzVector& momentum) const
1498 {
1499 // Return the direction and the momentum (GeV/c) of the track
1500 // currently being transported
1501 // TRACKR.ptrack = momentum of the particle (not always defined, if
1502 //               < 0 must be obtained from etrack) 
1503 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1504 // TRACKR.etrack = total energy of the particle
1505 // TRACKR.jtrack = identity number of the particle
1506 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1507   FlukaCallerCode_t  caller = GetCaller();
1508   FlukaProcessCode_t icode  = GetIcode();
1509   
1510   if (caller  != kEEDRAW         && 
1511       caller  != kMGResumedTrack && 
1512       caller  != kSODRAW         &&
1513       caller  != kUSDRAW         &&
1514       (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
1515       if (TRACKR.ptrack >= 0) {
1516           momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
1517           momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
1518           momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
1519           momentum.SetE(TRACKR.etrack);
1520           return;
1521       }
1522       else {
1523           Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
1524           momentum.SetPx(p*TRACKR.cxtrck);
1525           momentum.SetPy(p*TRACKR.cytrck);
1526           momentum.SetPz(p*TRACKR.cztrck);
1527           momentum.SetE(TRACKR.etrack);
1528           return;
1529       }
1530   } else if  (caller == kMGResumedTrack) {
1531       momentum.SetPx(TRACKR.spausr[4]);
1532       momentum.SetPy(TRACKR.spausr[5]);
1533       momentum.SetPz(TRACKR.spausr[6]);
1534       momentum.SetE (TRACKR.spausr[7]);
1535       return;
1536   } else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
1537       momentum.SetPx(0.);
1538       momentum.SetPy(0.);
1539       momentum.SetPz(0.);
1540       momentum.SetE(TrackMass());
1541       
1542   } else if (caller == kSODRAW) {
1543       Int_t ist  = FLKSTK.npflka;
1544       Double_t p = FLKSTK.pmoflk[ist];
1545       Int_t ifl  = FLKSTK.iloflk[ist];
1546       Double_t m = PAPROP.am[ifl + 6];
1547       Double_t e = TMath::Sqrt(p * p + m * m);
1548       momentum.SetPx(p * FLKSTK.txflk[ist]);
1549       momentum.SetPy(p * FLKSTK.tyflk[ist]);
1550       momentum.SetPz(p * FLKSTK.tzflk[ist]);
1551       momentum.SetE(e);
1552   } else if (caller == kUSDRAW) {
1553       if (icode == kEMFSCObrems  || 
1554           icode == kEMFSCOmoller || 
1555           icode == kEMFSCObhabha || 
1556           icode == kEMFSCOcompton ) 
1557       {
1558           momentum.SetPx(fPint[0]);
1559           momentum.SetPy(fPint[1]);
1560           momentum.SetPz(fPint[2]);
1561           momentum.SetE(fPint[3]);
1562       } else if (icode == kKASKADdray  || 
1563                  icode == kKASKADbrems || 
1564                  icode == kKASKADpair) {
1565           momentum.SetPx(GENSTK.plr[0] * GENSTK.cxr[0]);
1566           momentum.SetPy(GENSTK.plr[0] * GENSTK.cyr[0]);
1567           momentum.SetPz(GENSTK.plr[0] * GENSTK.czr[0]);
1568           momentum.SetE (GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6]);
1569       } else {
1570           Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack 
1571                             - ParticleMassFPC(TRACKR.jtrack) 
1572                             * ParticleMassFPC(TRACKR.jtrack));
1573           momentum.SetPx(p*TRACKR.cxtrck);
1574           momentum.SetPy(p*TRACKR.cytrck);
1575           momentum.SetPz(p*TRACKR.cztrck);
1576           momentum.SetE(TRACKR.etrack);
1577       }
1578   }
1579   else
1580     Warning("TrackMomentum","momentum not available");
1581 }
1582
1583 //______________________________________________________________________________ 
1584 void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e) const
1585 {
1586 // Return the direction and the momentum (GeV/c) of the track
1587 // currently being transported
1588 // TRACKR.ptrack = momentum of the particle (not always defined, if
1589 //               < 0 must be obtained from etrack) 
1590 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1591 // TRACKR.etrack = total energy of the particle
1592 // TRACKR.jtrack = identity number of the particle
1593 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1594   FlukaCallerCode_t   caller = GetCaller();
1595   FlukaProcessCode_t  icode  = GetIcode();
1596   if (caller != kEEDRAW         && 
1597       caller != kMGResumedTrack && 
1598       caller != kSODRAW         &&
1599       caller != kUSDRAW         &&
1600       (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
1601     if (TRACKR.ptrack >= 0) {
1602       px = TRACKR.ptrack*TRACKR.cxtrck;
1603       py = TRACKR.ptrack*TRACKR.cytrck;
1604       pz = TRACKR.ptrack*TRACKR.cztrck;
1605       e  = TRACKR.etrack;
1606       return;
1607     }
1608     else {
1609       Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) *  ParticleMassFPC(TRACKR.jtrack));
1610       px = p*TRACKR.cxtrck;
1611       py = p*TRACKR.cytrck;
1612       pz = p*TRACKR.cztrck;
1613       e  = TRACKR.etrack;
1614       return;
1615     }
1616   } else if (caller == kMGResumedTrack) {
1617       px = TRACKR.spausr[4];
1618       py = TRACKR.spausr[5];
1619       pz = TRACKR.spausr[6];
1620       e  = TRACKR.spausr[7];
1621       return;
1622   } else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
1623       px = 0.;
1624       py = 0.;
1625       pz = 0.;
1626       e  = TrackMass();
1627   } else if (caller == kSODRAW) {
1628       Int_t ist  = FLKSTK.npflka;
1629       Double_t p = FLKSTK.pmoflk[ist];
1630       Int_t ifl  = FLKSTK.iloflk[ist];
1631       Double_t m = PAPROP.am[ifl + 6];
1632                e = TMath::Sqrt(p * p + m * m);
1633       px = p * FLKSTK.txflk[ist];
1634       py = p * FLKSTK.tyflk[ist];
1635       pz = p * FLKSTK.tzflk[ist];
1636   } else if (caller == kUSDRAW) {
1637       if (icode == kEMFSCObrems  || 
1638           icode == kEMFSCOmoller || 
1639           icode == kEMFSCObhabha || 
1640           icode == kEMFSCOcompton ) 
1641       {
1642           px = fPint[0];
1643           py = fPint[1];
1644           pz = fPint[2];
1645           e  = fPint[3];
1646       } else if (icode == kKASKADdray  || 
1647                  icode == kKASKADbrems || 
1648                  icode == kKASKADpair) {
1649           px = GENSTK.plr[0] * GENSTK.cxr[0];
1650           py = GENSTK.plr[0] * GENSTK.cyr[0];
1651           pz = GENSTK.plr[0] * GENSTK.czr[0];
1652           e  = GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6];
1653       } else {
1654           Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) *  ParticleMassFPC(TRACKR.jtrack));
1655           px = p*TRACKR.cxtrck;
1656           py = p*TRACKR.cytrck;
1657           pz = p*TRACKR.cztrck;
1658           e  = TRACKR.etrack;
1659       }
1660   }
1661   else
1662       Warning("TrackMomentum","momentum not available");
1663 }
1664
1665 //______________________________________________________________________________ 
1666 Double_t TFluka::TrackStep() const
1667 {
1668 // Return the length in centimeters of the current step
1669 // TRACKR.ctrack = total curved path
1670     FlukaCallerCode_t caller = GetCaller();
1671     if (caller == kMGDRAW) {
1672         Int_t i;
1673         if ((i = fPrimaryElectronIndex) > -1) {
1674             if (i > 0) {
1675                 return (fPIlength[i] - fPIlength[i-1]); 
1676             } else {
1677                 Double_t s (TRACKR.ctrack - (fPIlength[fNPI - 1] - fPIlength[0]));
1678                 return s;
1679             }
1680         } else {
1681             return TRACKR.ctrack;
1682         }
1683     } else if (caller == kBXEntering || caller == kBXExiting || 
1684                caller == kENDRAW     || caller == kUSDRAW || 
1685                caller == kUSTCKV     || caller == kMGResumedTrack ||
1686                caller == kSODRAW)
1687     {
1688         return 0.0;
1689     } else {
1690         Warning("TrackStep", "track step not available");
1691         return 0.0;
1692     }  
1693 }
1694
1695 //______________________________________________________________________________ 
1696 Double_t TFluka::TrackLength() const
1697 {
1698 // TRACKR.cmtrck = cumulative curved path since particle birth
1699   FlukaCallerCode_t caller = GetCaller();
1700   if (caller == kMGDRAW) {
1701       Int_t i;
1702       if ((i = fPrimaryElectronIndex) > -1) {
1703           return fPIlength[i];
1704       } else {
1705           return TRACKR.cmtrck;
1706       }
1707       
1708   } else if (caller == kBXEntering || caller == kBXExiting || 
1709              caller == kENDRAW || caller == kUSDRAW || caller == kUSTCKV) 
1710       return TRACKR.cmtrck;
1711   else if (caller == kMGResumedTrack) 
1712       return TRACKR.spausr[8];
1713   else if (caller == kSODRAW)
1714       return 0.0;
1715   else {
1716       Warning("TrackLength", "track length not available for caller %5d \n", caller);
1717       return 0.0;
1718   }
1719 }
1720
1721
1722 //______________________________________________________________________________ 
1723 Double_t TFluka::TrackTime() const
1724 {
1725 // Return the current time of flight of the track being transported
1726 // TRACKR.atrack = age of the particle
1727   FlukaCallerCode_t caller = GetCaller();
1728   if (caller == kMGDRAW) {
1729       Int_t i;
1730       if ((i = fPrimaryElectronIndex) > -1) {
1731           Double_t t = fPItime[i];
1732           return t;
1733       } else {
1734           return TRACKR.atrack;
1735       }
1736   } else if (caller == kBXEntering || caller == kBXExiting || 
1737              caller == kENDRAW     || caller == kUSDRAW    || 
1738              caller == kUSTCKV)
1739     return TRACKR.atrack;
1740   else if (caller == kMGResumedTrack)
1741     return TRACKR.spausr[3];
1742   else if (caller == kSODRAW) {
1743       return (FLKSTK.agestk[FLKSTK.npflka]);
1744   }
1745   else {
1746     Warning("TrackTime", "track time not available");
1747     return 0.0;
1748   }   
1749 }
1750
1751 //______________________________________________________________________________ 
1752 Double_t TFluka::Edep() const
1753 {
1754 // Energy deposition
1755 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
1756 // -->local energy deposition (the value and the point are not recorded in TRACKR)
1757 //    but in the variable "rull" of the procedure "endraw.cxx"
1758 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
1759 // -->no energy loss along the track
1760 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
1761 // -->energy loss distributed along the track
1762 // TRACKR.dtrack = energy deposition of the jth deposition event
1763
1764   // If coming from bxdraw we have 2 steps of 0 length and 0 edep
1765   // If coming from usdraw we just signal particle production - no edep
1766   // If just first time after resuming, no edep for the primary
1767   FlukaCallerCode_t caller = GetCaller();
1768     
1769   if (caller == kBXExiting || caller == kBXEntering || 
1770       caller == kUSDRAW    || caller == kMGResumedTrack ||
1771       caller == kSODRAW) 
1772       return 0.0;
1773   Double_t sum = 0;
1774   Int_t i = -1;
1775   
1776   // Material with primary ionisation activated but number of primary electrons nprim = 0
1777   if (fPrimaryElectronIndex == -2) return 0.0;
1778   // nprim > 0
1779   if ((i = fPrimaryElectronIndex) > -1) {
1780       // Primary ionisation
1781       sum = GetPrimaryElectronKineticEnergy(i);
1782       if (sum > 100.) {
1783           printf("edep > 100. %d %d %f \n", i, ALLDLT.nalldl, sum);
1784       }
1785       return sum;
1786   } else {
1787       // Normal ionisation
1788       if (TRACKR.mtrack > 1) printf("Edep: %6d\n", TRACKR.mtrack);
1789       
1790       for ( Int_t j=0;j<TRACKR.mtrack;j++) {
1791           sum +=TRACKR.dtrack[j];  
1792       }
1793       if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
1794           return fRull + sum;
1795       else {
1796           return sum;
1797       }
1798   }
1799 }
1800
1801 //______________________________________________________________________________ 
1802 Int_t TFluka::CorrectFlukaId() const
1803 {
1804    // since we don't put photons and e- created bellow transport cut on the vmc stack
1805    // and there is a call to endraw for energy deposition for each of them
1806    // and they have the track number of their parent, but different identity (pdg)
1807    // so we want to assign also their parent identity.
1808
1809    if( (IsTrackStop())
1810         && TRACKR.ispusr[mkbmx2 - 4] == TRACKR.ispusr[mkbmx2 - 1]
1811         && TRACKR.jtrack != TRACKR.ispusr[mkbmx2 - 3] ) {
1812       if (fVerbosityLevel >=3)
1813          cout << "CorrectFlukaId() for icode=" << GetIcode()
1814                << " track=" << TRACKR.ispusr[mkbmx2 - 1]
1815                << " current PDG=" << PDGFromId(TRACKR.jtrack)
1816                << " assign parent PDG=" << PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
1817       return TRACKR.ispusr[mkbmx2 - 3]; // assign parent identity
1818    }
1819    if (TRACKR.jtrack <= 64){
1820        return TRACKR.jtrack;
1821    } else {
1822        return TRACKR.j0trck;
1823    }
1824 }
1825
1826
1827 //______________________________________________________________________________ 
1828 Int_t TFluka::TrackPid() const
1829 {
1830 // Return the id of the particle transported
1831 // TRACKR.jtrack = identity number of the particle
1832   FlukaCallerCode_t caller = GetCaller();
1833   if (caller != kEEDRAW && caller != kSODRAW) {
1834      return PDGFromId( CorrectFlukaId() );
1835   }
1836   else if (caller == kSODRAW) {
1837       return PDGFromId(FLKSTK.iloflk[FLKSTK.npflka]);
1838   }
1839   else
1840     return -1000;
1841 }
1842
1843 //______________________________________________________________________________ 
1844 Double_t TFluka::TrackCharge() const
1845 {
1846 // Return charge of the track currently transported
1847 // PAPROP.ichrge = electric charge of the particle
1848 // TRACKR.jtrack = identity number of the particle
1849     
1850   FlukaCallerCode_t caller = GetCaller();
1851   if (caller != kEEDRAW && caller != kSODRAW) 
1852      return PAPROP.ichrge[CorrectFlukaId() + 6];
1853   else if (caller == kSODRAW) {
1854       Int_t ifl =  PDGFromId(FLKSTK.iloflk[FLKSTK.npflka]);
1855       return PAPROP.ichrge[ifl + 6];
1856   }
1857   else
1858     return -1000.0;
1859 }
1860
1861 //______________________________________________________________________________ 
1862 Double_t TFluka::TrackMass() const
1863 {
1864 // PAPROP.am = particle mass in GeV
1865 // TRACKR.jtrack = identity number of the particle
1866   FlukaCallerCode_t caller = GetCaller();
1867   if (caller != kEEDRAW && caller != kSODRAW)
1868      return PAPROP.am[CorrectFlukaId()+6];
1869   else if (caller == kSODRAW) {
1870       Int_t ifl =  FLKSTK.iloflk[FLKSTK.npflka];
1871       return PAPROP.am[ifl + 6];
1872   }
1873   else
1874     return -1000.0;
1875 }
1876
1877 //______________________________________________________________________________ 
1878 Double_t TFluka::Etot() const
1879 {
1880 // TRACKR.etrack = total energy of the particle
1881   FlukaCallerCode_t  caller = GetCaller();
1882   FlukaProcessCode_t icode  = GetIcode();
1883   if (caller != kEEDRAW && caller != kSODRAW && caller != kUSDRAW)
1884   {
1885       return TRACKR.etrack;
1886   } else if (caller == kUSDRAW) {
1887       if (icode == kEMFSCObrems  || 
1888           icode == kEMFSCOmoller || 
1889           icode == kEMFSCObhabha || 
1890           icode == kEMFSCOcompton ) {
1891           return  fPint[3];
1892       }
1893       else if (icode == kKASKADdray  || 
1894                icode == kKASKADbrems || 
1895                icode == kKASKADpair) {
1896           return (GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6]);      
1897       } else {
1898           return TRACKR.etrack;
1899       }
1900       
1901   }
1902   else if (caller == kSODRAW) {
1903       Int_t ist  = FLKSTK.npflka;
1904       Double_t p = FLKSTK.pmoflk[ist];
1905       Int_t ifl  = FLKSTK.iloflk[ist];
1906       Double_t m = PAPROP.am[ifl + 6];
1907       Double_t e = TMath::Sqrt(p * p + m * m);
1908       return e;
1909   }
1910   printf("Etot %5d %5d \n", caller, icode);
1911   
1912   return -1000.0;
1913 }
1914
1915 //
1916 // track status
1917 //
1918 //______________________________________________________________________________ 
1919 Bool_t   TFluka::IsNewTrack() const
1920 {
1921 // Return true for the first call of Stepping()
1922    return fTrackIsNew;
1923 }
1924
1925 void     TFluka::SetTrackIsNew(Bool_t flag)
1926 {
1927 // Return true for the first call of Stepping()
1928    fTrackIsNew = flag;
1929
1930 }
1931
1932
1933 //______________________________________________________________________________ 
1934 Bool_t   TFluka::IsTrackInside() const
1935 {
1936 // True if the track is not at the boundary of the current volume
1937 // In Fluka a step is always inside one kind of material
1938 // If the step would go behind the region of one material,
1939 // it will be shortened to reach only the boundary.
1940 // Therefore IsTrackInside() is always true.
1941   FlukaCallerCode_t caller = GetCaller();
1942   if (caller == kBXEntering || caller == kBXExiting)
1943     return 0;
1944   else
1945     return 1;
1946 }
1947
1948 //______________________________________________________________________________ 
1949 Bool_t   TFluka::IsTrackEntering() const
1950 {
1951 // True if this is the first step of the track in the current volume
1952
1953   FlukaCallerCode_t caller = GetCaller();
1954   if (caller == kBXEntering)
1955     return 1;
1956   else return 0;
1957 }
1958
1959 //______________________________________________________________________________ 
1960 Bool_t   TFluka::IsTrackExiting() const
1961 {
1962 // True if track is exiting volume
1963 //
1964   FlukaCallerCode_t caller = GetCaller();
1965   if (caller == kBXExiting)
1966     return 1;
1967   else return 0;
1968 }
1969
1970 //______________________________________________________________________________ 
1971 Bool_t   TFluka::IsTrackOut() const
1972 {
1973 // True if the track is out of the setup
1974 // means escape
1975   FlukaProcessCode_t icode = GetIcode();
1976     
1977   if (icode == kKASKADescape ||
1978       icode == kEMFSCOescape ||
1979       icode == kKASNEUescape ||
1980       icode == kKASHEAescape ||
1981       icode == kKASOPHescape) 
1982        return 1;
1983   else return 0;
1984 }
1985
1986 //______________________________________________________________________________ 
1987 Bool_t   TFluka::IsTrackDisappeared() const
1988 {
1989 // All inelastic interactions and decays
1990 // fIcode from usdraw
1991   FlukaProcessCode_t icode = GetIcode();
1992   if (icode == kKASKADinelint    || // inelastic interaction
1993       icode == kKASKADdecay      || // particle decay
1994       icode == kKASKADdray       || // delta ray generation by hadron
1995       icode == kKASKADpair       || // direct pair production
1996       icode == kKASKADbrems      || // bremsstrahlung (muon)
1997       icode == kEMFSCObrems      || // bremsstrahlung (electron)
1998       icode == kEMFSCOmoller     || // Moller scattering
1999       icode == kEMFSCObhabha     || // Bhaba scattering
2000       icode == kEMFSCOanniflight || // in-flight annihilation
2001       icode == kEMFSCOannirest   || // annihilation at rest
2002       icode == kEMFSCOpair       || // pair production
2003       icode == kEMFSCOcompton    || // Compton scattering
2004       icode == kEMFSCOphotoel    || // Photoelectric effect
2005       icode == kKASNEUhadronic   || // hadronic interaction
2006       icode == kKASHEAdray          // delta-ray
2007       ) return 1;
2008   else return 0;
2009 }
2010
2011 //______________________________________________________________________________ 
2012 Bool_t   TFluka::IsTrackStop() const
2013 {
2014 // True if the track energy has fallen below the threshold
2015 // means stopped by signal or below energy threshold
2016   FlukaProcessCode_t icode = GetIcode();
2017   if (icode == kKASKADstopping  || // stopping particle
2018       icode == kKASKADtimekill  || // time kill 
2019       icode == kEMFSCOstopping1 || // below user-defined cut-off
2020       icode == kEMFSCOstopping2 || // below user cut-off
2021       icode == kEMFSCOtimekill  || // time kill
2022       icode == kKASNEUstopping  || // neutron below threshold
2023       icode == kKASNEUtimekill  || // time kill
2024       icode == kKASHEAtimekill  || // time kill
2025       icode == kKASOPHtimekill) return 1; // time kill
2026   else return 0;
2027 }
2028
2029 //______________________________________________________________________________ 
2030 Bool_t   TFluka::IsTrackAlive() const
2031 {
2032 // Means not disappeared or not out
2033     FlukaProcessCode_t icode = GetIcode();
2034     
2035     if (IsTrackOut()               || 
2036         IsTrackStop()              ||
2037         icode == kKASKADinelint    || // inelastic interaction
2038         icode == kKASKADdecay      || // particle decay
2039         icode == kEMFSCOanniflight || // in-flight annihilation
2040         icode == kEMFSCOannirest   || // annihilation at rest
2041         icode == kEMFSCOpair       || // pair production
2042         icode == kEMFSCOphotoel    || // Photoelectric effect
2043         icode == kKASNEUhadronic      // hadronic interaction
2044         ) 
2045     {
2046         // Exclude the cases for which the particle has disappeared (paused) but will reappear later (= alive).
2047         return 0;
2048     } else {
2049         return 1;
2050     }
2051 }
2052
2053 //
2054 // secondaries
2055 //
2056
2057 //______________________________________________________________________________ 
2058 Int_t TFluka::NSecondaries() const
2059
2060 {
2061 // Number of secondary particles generated in the current step
2062 // GENSTK.np = number of secondaries except light and heavy ions
2063 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
2064     FlukaCallerCode_t caller = GetCaller();
2065     if (caller == kUSDRAW)  // valid only after usdraw
2066         return GENSTK.np + FHEAVY.npheav;
2067     else if (caller == kUSTCKV) {
2068         // Cerenkov Photon production
2069         return fNCerenkov;
2070     }
2071     return 0;
2072 } // end of NSecondaries
2073
2074 //______________________________________________________________________________ 
2075 void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
2076                 TLorentzVector& position, TLorentzVector& momentum)
2077 {
2078 // Copy particles from secondary stack to vmc stack
2079 //
2080
2081     FlukaCallerCode_t caller = GetCaller();
2082     if (caller == kUSDRAW) {  // valid only after usdraw
2083         if (GENSTK.np > 0) {
2084             // Hadronic interaction
2085             if (isec >= 0 && isec < GENSTK.np) {
2086                 particleId = PDGFromId(GENSTK.kpart[isec]);
2087                 position.SetX(fXsco);
2088                 position.SetY(fYsco);
2089                 position.SetZ(fZsco);
2090                 position.SetT(TRACKR.atrack);
2091                 momentum.SetPx(GENSTK.plr[isec]*GENSTK.cxr[isec]);
2092                 momentum.SetPy(GENSTK.plr[isec]*GENSTK.cyr[isec]);
2093                 momentum.SetPz(GENSTK.plr[isec]*GENSTK.czr[isec]);
2094                 momentum.SetE(GENSTK.tki[isec] + PAPROP.am[GENSTK.kpart[isec]+6]);
2095             }
2096             else if (isec >= GENSTK.np && isec < GENSTK.np + FHEAVY.npheav) {
2097                 Int_t jsec = isec - GENSTK.np;
2098                 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
2099                 position.SetX(fXsco);
2100                 position.SetY(fYsco);
2101                 position.SetZ(fZsco);
2102                 position.SetT(TRACKR.atrack);
2103                 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
2104                 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
2105                 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
2106                 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
2107                     momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
2108                 else if (FHEAVY.tkheav[jsec] > 6)
2109                     momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
2110             }
2111             else
2112                 Warning("GetSecondary","isec out of range");
2113         }
2114     } else if (caller == kUSTCKV) {
2115         Int_t index = OPPHST.lstopp - isec;
2116         position.SetX(OPPHST.xoptph[index]);
2117         position.SetY(OPPHST.yoptph[index]);
2118         position.SetZ(OPPHST.zoptph[index]);
2119         position.SetT(OPPHST.agopph[index]);
2120         Double_t p = OPPHST.poptph[index];
2121         
2122         momentum.SetPx(p * OPPHST.txopph[index]);
2123         momentum.SetPy(p * OPPHST.tyopph[index]);
2124         momentum.SetPz(p * OPPHST.tzopph[index]);
2125         momentum.SetE(p);
2126     }
2127     else
2128         Warning("GetSecondary","no secondaries available");
2129     
2130 } // end of GetSecondary
2131
2132
2133 //______________________________________________________________________________ 
2134 TMCProcess TFluka::ProdProcess(Int_t) const
2135
2136 {
2137 // Name of the process that has produced the secondary particles
2138 // in the current step
2139
2140     Int_t mugamma = (TRACKR.jtrack == kFLUKAphoton || 
2141                      TRACKR.jtrack == kFLUKAmuplus ||
2142                      TRACKR.jtrack == kFLUKAmuminus);
2143     FlukaProcessCode_t icode = GetIcode();
2144
2145     if      (icode == kKASKADdecay)                                   return kPDecay;
2146     else if (icode == kKASKADpair || icode == kEMFSCOpair)            return kPPair;
2147     else if (icode == kEMFSCOcompton)                                 return kPCompton;
2148     else if (icode == kEMFSCOphotoel)                                 return kPPhotoelectric;
2149     else if (icode == kKASKADbrems      || icode == kEMFSCObrems)     return kPBrem;
2150     else if (icode == kKASKADdray       || icode == kKASHEAdray)      return kPDeltaRay;
2151     else if (icode == kEMFSCOmoller     || icode == kEMFSCObhabha)    return kPDeltaRay;
2152     else if (icode == kEMFSCOanniflight || icode == kEMFSCOannirest)  return kPAnnihilation;
2153     else if (icode == kKASKADinelint) {
2154         if (!mugamma)                                                 return kPHadronic;
2155         else if (TRACKR.jtrack == kFLUKAphoton)                       return kPPhotoFission;
2156         else                                                          return kPMuonNuclear;
2157     }
2158     else if (icode == kEMFSCOrayleigh)                                return kPRayleigh;
2159 // Fluka codes 100, 300 and 400 still to be investigasted
2160     else                                                              return kPNoProcess;
2161 }
2162
2163
2164 Int_t TFluka::StepProcesses(TArrayI &proc) const
2165 {
2166   //
2167   // Return processes active in the current step
2168   //
2169     FlukaProcessCode_t icode   = GetIcode();
2170     FlukaCallerCode_t  caller  = GetCaller();
2171     proc.Set(1);
2172     TMCProcess iproc;
2173     if (caller == kBXEntering || caller == kBXExiting || caller == kEEDRAW || caller == kSODRAW) {
2174         iproc = kPTransportation;
2175     }
2176     else if (caller == kUSTCKV) {
2177         iproc = kPCerenkov;
2178     } else {
2179         switch (icode) {
2180         case kEMFSCO:
2181             if (Edep() > 0.) {
2182                 iproc = kPEnergyLoss;
2183             } else {
2184                 iproc = kPTransportation;
2185             }
2186             break;
2187         case kKASKAD:
2188             if (Edep() > 0.) {
2189                 iproc = kPEnergyLoss;
2190             } else {
2191                 iproc = kPTransportation;
2192             }
2193             break;
2194         case kKASHEA:
2195         case kKASNEU:
2196         case kKASOPH:
2197         case kKASKADescape:
2198         case kEMFSCOescape:
2199         case kKASNEUescape:
2200         case kKASHEAescape:
2201         case kKASOPHescape:
2202             iproc = kPTransportation;
2203             break;
2204         case kKASKADtimekill:
2205         case kEMFSCOtimekill:
2206         case kKASNEUtimekill:
2207         case kKASHEAtimekill:
2208         case kKASOPHtimekill:
2209             iproc =  kPTOFlimit;
2210             break;
2211         case kKASKADstopping:
2212         case kEMFSCOstopping1:
2213         case kEMFSCOstopping2:
2214         case kKASNEUstopping:
2215             iproc = kPStop;
2216             break; 
2217         case kKASKADinelint:
2218         case kKASNEUhadronic:
2219             iproc = kPHadronic;
2220             break;
2221         case kKASKADinelarecoil:
2222             iproc = kPHadronic;
2223             break;
2224         case kKASKADnelint:
2225             iproc = kPHElastic;
2226             break;
2227         case kKASOPHabsorption:
2228             iproc = kPLightAbsorption;
2229             break;
2230         case kKASOPHrefraction:
2231             iproc = kPLightRefraction;
2232             break;
2233         case kEMFSCOlocaldep : 
2234             iproc = kPPhotoelectric;
2235             break;
2236         default:
2237             iproc = ProdProcess(0);
2238         }
2239     }
2240     
2241     proc[0] = iproc;
2242     return 1;
2243 }
2244 //______________________________________________________________________________ 
2245 Int_t TFluka::VolId2Mate(Int_t id) const
2246 {
2247 //
2248 // Returns the material number for a given volume ID
2249 //
2250    return fMCGeo->VolId2Mate(id);
2251 }
2252
2253 //______________________________________________________________________________ 
2254 const char* TFluka::VolName(Int_t id) const
2255 {
2256 //
2257 // Returns the volume name for a given volume ID
2258 //
2259    return fMCGeo->VolName(id);
2260 }
2261
2262 Int_t TFluka::MediumId(const Text_t* mediumName) const
2263 {
2264     //
2265     // Return the unique medium id for medium with name mediumName
2266     TList *medlist = gGeoManager->GetListOfMedia();
2267     TGeoMedium* med = (TGeoMedium*) medlist->FindObject(mediumName);
2268     if (med) {
2269         return (med->GetId());
2270     } else {
2271         return (-1);
2272     }
2273 }
2274
2275 //______________________________________________________________________________ 
2276 Int_t TFluka::VolId(const Text_t* volName) const
2277 {
2278 //
2279 // Converts from volume name to volume ID.
2280 // Time consuming. (Only used during set-up)
2281 // Could be replaced by hash-table
2282 //
2283     char sname[20];
2284     Int_t len;
2285     strncpy(sname, volName, len = strlen(volName));
2286     sname[len] = 0;
2287     while (sname[len - 1] == ' ') sname[--len] = 0;
2288     return fMCGeo->VolId(sname);
2289 }
2290
2291 //______________________________________________________________________________ 
2292 Int_t TFluka::CurrentVolID(Int_t& copyNo) const
2293 {
2294 //
2295 // Return the logical id and copy number corresponding to the current fluka region
2296 //
2297   if (gGeoManager->IsOutside()) return 0;
2298   TGeoNode *node = gGeoManager->GetCurrentNode();
2299   copyNo = node->GetNumber();
2300   Int_t id = node->GetVolume()->GetNumber();
2301   return id;
2302
2303
2304 //______________________________________________________________________________ 
2305 Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
2306 {
2307 //
2308 // Return the logical id and copy number of off'th mother 
2309 // corresponding to the current fluka region
2310 //
2311   if (off<0 || off>gGeoManager->GetLevel()) return 0;
2312   if (off==0) return CurrentVolID(copyNo);
2313   TGeoNode *node = gGeoManager->GetMother(off);
2314   if (!node) return 0;
2315   copyNo = node->GetNumber();
2316   return node->GetVolume()->GetNumber();
2317 }
2318
2319 //______________________________________________________________________________ 
2320 const char* TFluka::CurrentVolName() const
2321 {
2322 //
2323 // Return the current volume name
2324 //
2325   if (gGeoManager->IsOutside()) return "OutOfWorld";
2326   return gGeoManager->GetCurrentVolume()->GetName();
2327 }
2328
2329 //______________________________________________________________________________ 
2330 const char* TFluka::CurrentVolOffName(Int_t off) const
2331 {
2332 //
2333 // Return the volume name of the off'th mother of the current volume
2334 //
2335   if (off<0 || off>gGeoManager->GetLevel()) return 0;
2336   if (off==0) return CurrentVolName();
2337   TGeoNode *node = gGeoManager->GetMother(off);
2338   if (!node) return 0;
2339   return node->GetVolume()->GetName();
2340 }
2341
2342 const char* TFluka::CurrentVolPath() {
2343   // Return the current volume path
2344   return gGeoManager->GetPath(); 
2345 }
2346 //______________________________________________________________________________ 
2347 Int_t TFluka::CurrentMaterial(Float_t & a, Float_t & z, 
2348                       Float_t & dens, Float_t & radl, Float_t & absl) const
2349 {
2350 //
2351 //  Return the current medium number and material properties
2352 //
2353   Int_t copy;
2354   Int_t id  =  TFluka::CurrentVolID(copy);
2355   Int_t med =  TFluka::VolId2Mate(id);
2356   TGeoVolume*     vol = gGeoManager->GetCurrentVolume();
2357   TGeoMaterial*   mat = vol->GetMaterial();
2358   a    = mat->GetA();
2359   z    = mat->GetZ();
2360   dens = mat->GetDensity();
2361   radl = mat->GetRadLen();
2362   absl = mat->GetIntLen();
2363   
2364   return med;
2365 }
2366
2367 //______________________________________________________________________________ 
2368 void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
2369 {
2370 // Transforms a position from the world reference frame
2371 // to the current volume reference frame.
2372 //
2373 //  Geant3 desription:
2374 //  ==================
2375 //       Computes coordinates XD (in DRS) 
2376 //       from known coordinates XM in MRS 
2377 //       The local reference system can be initialized by
2378 //         - the tracking routines and GMTOD used in GUSTEP
2379 //         - a call to GMEDIA(XM,NUMED)
2380 //         - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER) 
2381 //             (inverse routine is GDTOM) 
2382 //
2383 //        If IFLAG=1  convert coordinates 
2384 //           IFLAG=2  convert direction cosinus
2385 //
2386 // ---
2387    Double_t xmL[3], xdL[3];
2388    Int_t i;
2389    for (i=0;i<3;i++) xmL[i]=xm[i];
2390    if (iflag == 1) gGeoManager->MasterToLocal(xmL,xdL);
2391    else            gGeoManager->MasterToLocalVect(xmL,xdL);
2392    for (i=0;i<3;i++) xd[i] = xdL[i];
2393 }
2394   
2395 //______________________________________________________________________________ 
2396 void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
2397 {
2398 //
2399 // See Gmtod(Float_t*, Float_t*, Int_t)
2400 //
2401    if (iflag == 1) gGeoManager->MasterToLocal(xm,xd);
2402    else            gGeoManager->MasterToLocalVect(xm,xd);
2403 }
2404
2405 //______________________________________________________________________________ 
2406 void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
2407 {
2408 // Transforms a position from the current volume reference frame
2409 // to the world reference frame.
2410 //
2411 //  Geant3 desription:
2412 //  ==================
2413 //  Computes coordinates XM (Master Reference System
2414 //  knowing the coordinates XD (Detector Ref System)
2415 //  The local reference system can be initialized by
2416 //    - the tracking routines and GDTOM used in GUSTEP
2417 //    - a call to GSCMED(NLEVEL,NAMES,NUMBER)
2418 //        (inverse routine is GMTOD)
2419 // 
2420 //   If IFLAG=1  convert coordinates
2421 //      IFLAG=2  convert direction cosinus
2422 //
2423 // ---
2424    Double_t xmL[3], xdL[3];
2425    Int_t i;
2426    for (i=0;i<3;i++) xdL[i] = xd[i];
2427    if (iflag == 1) gGeoManager->LocalToMaster(xdL,xmL);
2428    else            gGeoManager->LocalToMasterVect(xdL,xmL);
2429    for (i=0;i<3;i++) xm[i]=xmL[i];
2430 }
2431
2432 //______________________________________________________________________________ 
2433 void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
2434 {
2435 //
2436 // See Gdtom(Float_t*, Float_t*, Int_t)
2437 //
2438    if (iflag == 1) gGeoManager->LocalToMaster(xd,xm);
2439    else            gGeoManager->LocalToMasterVect(xd,xm);
2440 }
2441
2442 //______________________________________________________________________________
2443 TObjArray *TFluka::GetFlukaMaterials()
2444 {
2445 //
2446 // Get array of Fluka materials
2447    return fGeom->GetMatList();
2448 }   
2449
2450 //______________________________________________________________________________
2451 void TFluka::SetMreg(Int_t l, Int_t lttc) 
2452 {
2453 // Set current fluka region
2454    fCurrentFlukaRegion = l;
2455    fGeom->SetMreg(l,lttc);
2456 }
2457
2458
2459
2460
2461 //______________________________________________________________________________
2462 TString TFluka::ParticleName(Int_t pdg) const
2463 {
2464     // Return particle name for particle with pdg code pdg.
2465     Int_t ifluka = IdFromPDG(pdg);
2466     return TString((CHPPRP.btype[ifluka - kFLUKAcodemin]), 8);
2467 }
2468  
2469
2470 //______________________________________________________________________________
2471 Double_t TFluka::ParticleMass(Int_t pdg) const
2472 {
2473     // Return particle mass for particle with pdg code pdg.
2474     Int_t ifluka = IdFromPDG(pdg);
2475     return (PAPROP.am[ifluka - kFLUKAcodemin]);
2476 }
2477
2478 //______________________________________________________________________________
2479 Double_t TFluka::ParticleMassFPC(Int_t fpc) const
2480 {
2481     // Return particle mass for particle with Fluka particle code fpc
2482     return (PAPROP.am[fpc - kFLUKAcodemin]);
2483 }
2484
2485 //______________________________________________________________________________
2486 Double_t TFluka::ParticleCharge(Int_t pdg) const
2487 {
2488     // Return particle charge for particle with pdg code pdg.
2489     Int_t ifluka = IdFromPDG(pdg);
2490     return Double_t(PAPROP.ichrge[ifluka - kFLUKAcodemin]);
2491 }
2492
2493 //______________________________________________________________________________
2494 Double_t TFluka::ParticleLifeTime(Int_t pdg) const
2495 {
2496     // Return particle lifetime for particle with pdg code pdg.
2497     Int_t ifluka = IdFromPDG(pdg);
2498     return (PAPROP.tmnlf[ifluka - kFLUKAcodemin]);
2499 }
2500
2501 //______________________________________________________________________________
2502 void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t& charge, Float_t& tlife)
2503 {
2504     // Retrieve particle properties for particle with pdg code pdg.
2505     
2506     strcpy(name, ParticleName(pdg).Data());
2507     type   = ParticleMCType(pdg);
2508     mass   = ParticleMass(pdg);
2509     charge = ParticleCharge(pdg);
2510     tlife  = ParticleLifeTime(pdg);
2511 }
2512
2513 //______________________________________________________________________________
2514 void TFluka::PrintHeader()
2515 {
2516     //
2517     // Print a header
2518     printf("\n");
2519     printf("\n");    
2520     printf("------------------------------------------------------------------------------\n");
2521     printf("- You are using the TFluka Virtual Monte Carlo Interface to FLUKA.           -\n");    
2522     printf("- Please see the file fluka.out for FLUKA output and licensing information.  -\n");    
2523     printf("------------------------------------------------------------------------------\n");
2524     printf("\n");
2525     printf("\n");    
2526 }
2527
2528
2529 #define pshckp pshckp_
2530 #define ustckv ustckv_
2531
2532
2533 extern "C" {
2534   void pshckp(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
2535               Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
2536               Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
2537   {
2538     //
2539     // Pushes one cerenkov photon to the stack
2540     //
2541     
2542     TFluka* fluka =  (TFluka*) gMC;
2543     TVirtualMCStack* cppstack = fluka->GetStack();
2544     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
2545     cppstack->PushTrack(0, parent, 50000050,
2546                         px, py, pz, e,
2547                         vx, vy, vz, tof,
2548                         polx, poly, polz,
2549                         kPCerenkov, ntr, wgt, 0);
2550     if (fluka->GetVerbosityLevel() >= 3)
2551             printf("pshckp: track=%d parent=%d lattc=%d %s\n", ntr, parent, TRACKR.lt1trk, fluka->CurrentVolName());
2552   }
2553     
2554     void ustckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
2555     {
2556         //
2557         // Calls stepping in order to signal cerenkov production
2558         //
2559         TFluka *fluka = (TFluka*)gMC;
2560         fluka->SetMreg(mreg, TRACKR.lt1trk); //LTCLCM.mlatm1);
2561         fluka->SetXsco(x);
2562         fluka->SetYsco(y);
2563         fluka->SetZsco(z);
2564         fluka->SetNCerenkov(nphot);
2565         fluka->SetCaller(kUSTCKV);
2566         if (fluka->GetVerbosityLevel() >= 3)
2567             printf("ustckv: %10d mreg=%d lattc=%d  newlat=%d (%f, %f, %f) edep=%f vol=%s\n",
2568                     nphot, mreg, TRACKR.lt1trk, LTCLCM.newlat, x, y, z, fluka->Edep(), fluka->CurrentVolName());
2569    
2570     // check region lattice consistency (debug Ernesto)
2571     // *****************************************************
2572    Int_t nodeId;
2573    Int_t volId = fluka->CurrentVolID(nodeId);
2574    Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
2575
2576    if( mreg != volId  && !gGeoManager->IsOutside() ) {
2577        cout << "  ustckv:   track=" << TRACKR.ispusr[mkbmx2-1] << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
2578             << " icode=" << fluka->GetIcode() << " gNstep=" << fluka->GetNstep() << endl
2579             << "               fluka   mreg=" << mreg << " mlttc=" << TRACKR.lt1trk << endl
2580             << "               TGeo   volId=" << volId << " crtlttc=" << crtlttc << endl
2581             << "     common TRACKR   lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
2582             << "     common LTCLCM   newlat=" << LTCLCM.newlat << " mlatld=" <<  LTCLCM.mlatld << endl
2583             << "                     mlatm1=" << LTCLCM.mlatm1 << " mltsen=" <<  LTCLCM.mltsen << endl
2584             << "                     mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
2585         if( TRACKR.lt1trk == crtlttc ) cout << "   *************************************************************" << endl;
2586     }
2587     // *****************************************************
2588
2589
2590
2591         (TVirtualMCApplication::Instance())->Stepping();
2592     }
2593 }
2594
2595 //______________________________________________________________________________
2596 void TFluka::AddParticlesToPdgDataBase() const
2597 {
2598
2599 //
2600 // Add particles to the PDG data base
2601
2602     TDatabasePDG *pdgDB = TDatabasePDG::Instance();
2603
2604     const Double_t kAu2Gev   = 0.9314943228;
2605     const Double_t khSlash   = 1.0545726663e-27;
2606     const Double_t kErg2Gev  = 1/1.6021773349e-3;
2607     const Double_t khShGev   = khSlash*kErg2Gev;
2608     const Double_t kYear2Sec = 3600*24*365.25;
2609 //
2610 // Ions
2611 //
2612   pdgDB->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE,
2613                      0,3,"Ion",TFlukaIon::GetIonPdg(1,2));
2614   pdgDB->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE,
2615                      khShGev/(12.33*kYear2Sec),3,"Ion",TFlukaIon::GetIonPdg(1,3));
2616   pdgDB->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,
2617                      khShGev/(12.33*kYear2Sec),6,"Ion",TFlukaIon::GetIonPdg(2,4));
2618   pdgDB->AddParticle("HE3","HE3",3*kAu2Gev+14.931e-3,kFALSE,
2619                      0,6,"Ion",TFlukaIon::GetIonPdg(2,3));
2620 //
2621 // Default user ion
2622   TFlukaIon::AddIon(12, 6);
2623   
2624 //
2625 //
2626 //
2627 // Special particles
2628 //
2629   pdgDB->AddParticle("Cherenkov","Cherenkov",0,kFALSE,
2630                      0,0,"Special",GetSpecialPdg(50));
2631   pdgDB->AddParticle("FeedbackPhoton","FeedbackPhoton",0,kFALSE,
2632                      0,0,"Special",GetSpecialPdg(51));
2633 }
2634
2635
2636 //
2637 // Info about primary ionization electrons
2638 //
2639
2640 //______________________________________________________________________________
2641 Int_t TFluka::GetNPrimaryElectrons()
2642 {
2643     // Get number of primary electrons
2644     return ALLDLT.nalldl;
2645 }
2646
2647 //______________________________________________________________________________
2648 Double_t TFluka::GetPrimaryElectronKineticEnergy(Int_t i) const
2649 {
2650     // Returns kinetic energy of primary electron i
2651
2652     Double_t ekin = -1.;
2653     
2654     if (i >= 0 && i < ALLDLT.nalldl) {
2655         ekin =  ALLDLT.talldl[i];
2656     } else {
2657         Warning("GetPrimaryElectronKineticEnergy",
2658                 "Primary electron index out of range %d %d \n",
2659                 i, ALLDLT.nalldl);
2660     }
2661     return ekin;
2662 }
2663
2664 void TFluka::GetPrimaryElectronPosition(Int_t i, Double_t& x, Double_t& y, Double_t& z, Double_t& t) const
2665 {
2666     // Returns position  of primary electron i
2667         if (i >= 0 && i < ALLDLT.nalldl) {
2668             x = ALLDLT.xalldl[i];
2669             y = ALLDLT.yalldl[i];
2670             z = ALLDLT.zalldl[i];
2671             t = ALLDLT.talldl[i];
2672             return;
2673         } else {
2674             Warning("GetPrimaryElectronPosition",
2675                     "Primary electron index out of range %d %d \n",
2676                     i, ALLDLT.nalldl);
2677             return;
2678         }
2679         return;
2680 }
2681
2682
2683 //__________________________________________________________________
2684 Int_t TFluka::GetSpecialPdg(Int_t number) const
2685 {
2686 // Numbering for special particles
2687
2688   return 50000000 + number;
2689 }                
2690
2691      
2692 void  TFluka::PrimaryIonisationStepping(Int_t nprim)
2693 {
2694 // Call Stepping for primary ionisation electrons
2695 // Protection against nprim > mxalld
2696 // Multiple steps for nprim > 0
2697     Int_t i;
2698     fNPI = nprim;
2699     if (nprim > 0) {
2700         CalcPrimaryIonisationTime();
2701         for (i = 0; i < nprim; i++) {
2702             SetCurrentPrimaryElectronIndex(i);
2703             (TVirtualMCApplication::Instance())->Stepping();
2704             if (i == 0) SetTrackIsNew(kFALSE);
2705         }       
2706     } else {
2707         // No primary electron ionisation
2708         // Call Stepping anyway but flag nprim = 0 as index = -2
2709         SetCurrentPrimaryElectronIndex(-2);
2710         (TVirtualMCApplication::Instance())->Stepping();
2711     }
2712     // Reset the index
2713     SetCurrentPrimaryElectronIndex(-1);
2714 }
2715
2716 //______________________________________________________________________
2717 Float_t* TFluka::CreateFloatArray(Double_t* array, Int_t size) const
2718 {
2719 // Converts Double_t* array to Float_t*,
2720 // !! The new array has to be deleted by user.
2721 // ---
2722     
2723   Float_t* floatArray;
2724   if (size>0) {
2725     floatArray = new Float_t[size];
2726     for (Int_t i=0; i<size; i++)
2727       if (array[i] >= FLT_MAX ) 
2728         floatArray[i] = FLT_MAX/100.;
2729       else      
2730         floatArray[i] = array[i];
2731   }
2732   else {
2733     //floatArray = 0;
2734     floatArray = new Float_t[1];
2735   }
2736   return floatArray;
2737 }
2738
2739 void TFluka::CalcPrimaryIonisationTime()
2740 {
2741     // Calculates the primary ionisation time
2742     if (fPItime) delete [] fPItime; 
2743     fPItime = new Double_t[fNPI];
2744     if (fPIlength) delete [] fPIlength;
2745     fPIlength = new Double_t[fNPI];
2746     //
2747     Double_t px, py, pz, e, t;
2748     TrackMomentum(px, py, pz, e);
2749     Double_t p    = TMath::Sqrt(px * px + py * py + pz * pz);
2750     Double_t beta = p / e;
2751     Double_t x0, y0, z0;
2752     fPItime[fNPI -1]   = TRACKR.atrack;
2753     fPIlength[fNPI -1] = TRACKR.cmtrck;
2754     GetPrimaryElectronPosition(fNPI - 1, x0, y0, z0, t);        
2755     if (fNPI > 1) {
2756         for (Int_t i = fNPI - 2; i > -1; i--) {
2757             Double_t x, y, z, t;
2758             GetPrimaryElectronPosition(i, x, y, z, t);  
2759             Double_t ds = TMath::Sqrt((x-x0) * (x-x0) + (y-y0) * (y-y0) + (z-z0) * (z-z0));
2760             fPItime[i]   = fPItime[i+1]   - ds / (beta * 2.99792458e10);
2761             fPIlength[i] = fPIlength[i+1] - ds;
2762             x0 = x; y0 = y; z0 = z;
2763         }
2764     }
2765     
2766 }
2767
2768 Bool_t TFluka::DefineIon(const char* name , Int_t z, Int_t a, Int_t q, Double_t exE, Double_t mass)
2769 {
2770     // User defined ion that can be used as a primary
2771     fUserIons = kTRUE;
2772     TFlukaIon::AddIon(name, z, a, q,exE, mass);
2773     return kTRUE;
2774 }