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