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