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