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