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