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