1ab0ab9ab9e15b9aca111ec559ad5fb9328b4f99
[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
49 #include "TVirtualMC.h"
50 #include "TMCProcess.h"
51 #include "TGeoManager.h"
52 #include "TGeoMaterial.h"
53 #include "TGeoMedium.h"
54 #include "TFlukaMCGeometry.h"
55 #include "TGeoMCGeometry.h"
56 #include "TFlukaCerenkov.h"
57 #include "TFlukaConfigOption.h"
58 #include "TFlukaScoringOption.h"
59 #include "TLorentzVector.h"
60 #include "TArrayI.h"
61
62 // Fluka methods that may be needed.
63 #ifndef WIN32 
64 # define flukam  flukam_
65 # define fluka_openinp fluka_openinp_
66 # define fluka_closeinp fluka_closeinp_
67 # define mcihad mcihad_
68 # define mpdgha mpdgha_
69 #else 
70 # define flukam  FLUKAM
71 # define fluka_openinp FLUKA_OPENINP
72 # define fluka_closeinp FLUKA_CLOSEINP
73 # define mcihad MCIHAD
74 # define mpdgha MPDGHA
75 #endif
76
77 extern "C" 
78 {
79   //
80   // Prototypes for FLUKA functions
81   //
82   void type_of_call flukam(const int&);
83   void type_of_call fluka_openinp(const int&, DEFCHARA);
84   void type_of_call fluka_closeinp(const int&);
85   int  type_of_call mcihad(const int&);
86   int  type_of_call mpdgha(const int&);
87 }
88
89 //
90 // Class implementation for ROOT
91 //
92 ClassImp(TFluka)
93
94 //
95 //----------------------------------------------------------------------------
96 // TFluka constructors and destructors.
97 //______________________________________________________________________________
98 TFluka::TFluka()
99   :TVirtualMC(),
100    fVerbosityLevel(0),
101    fInputFileName(""),
102    fProcesses(0), 
103    fCuts(0),
104    fUserScore(0)
105
106   //
107   // Default constructor
108   //
109    fGeneratePemf = kFALSE;
110    fNVolumes = 0;
111    fCurrentFlukaRegion = -1;
112    fGeom = 0;
113    fMCGeo = 0;
114    fMaterials = 0;
115    fDummyBoundary = 0;
116    fFieldFlag = 1;
117    fStopped   = 0;
118    fStopEvent = 0;
119    fStopRun   = 0;
120    fNEvent    = 0;
121
122  
123 //______________________________________________________________________________ 
124 TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
125   :TVirtualMC("TFluka",title, isRootGeometrySupported),
126    fVerbosityLevel(verbosity),
127    fInputFileName(""),
128    fTrackIsEntering(0),
129    fTrackIsExiting(0),
130    fTrackIsNew(0),
131    fProcesses(new TObjArray(100)),
132    fCuts(new TObjArray(100)), 
133    fUserScore(new TObjArray(100)) 
134 {
135   // create geometry interface
136    if (fVerbosityLevel >=3)
137        cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
138    SetCoreInputFileName();
139    SetInputFileName();
140    SetGeneratePemf(kFALSE);
141    fNVolumes      = 0;
142    fCurrentFlukaRegion = -1;
143    fDummyBoundary = 0;
144    fFieldFlag = 1;
145    fGeneratePemf = kFALSE;
146    fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kTRUE);
147    fGeom = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
148    if (verbosity > 2) fGeom->SetDebugMode(kTRUE);
149    fMaterials = 0;
150    fStopped   = 0;
151    fStopEvent = 0;
152    fStopRun   = 0;
153    fNEvent    = 0;
154 }
155
156 //______________________________________________________________________________ 
157 TFluka::~TFluka() {
158 // Destructor
159     if (fVerbosityLevel >=3)
160         cout << "<== TFluka::~TFluka() destructor called." << endl;
161     
162     delete fGeom;
163     delete fMCGeo;
164     
165     if (fCuts) {
166         fCuts->Delete();
167         delete fCuts;
168     }
169
170     if (fProcesses) {
171         fProcesses->Delete();
172         delete fProcesses;
173     }
174
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     fNVolumes = fGeom->NofVolumes();
195     fGeom->CreateFlukaMatFile("flukaMat.inp");   
196     if (fVerbosityLevel >=3) {
197        printf("== Number of volumes: %i\n ==", fNVolumes);
198        cout << "\t* InitPhysics() - Prepare input file to be called" << endl; 
199     }   
200     // now we have TGeo geometry created and we have to patch FlukaVmc.inp
201     // with the material mapping file FlukaMat.inp
202 }
203
204
205 //______________________________________________________________________________ 
206 void TFluka::FinishGeometry() {
207 //
208 // Build-up table with region to medium correspondance
209 //
210   if (fVerbosityLevel >=3) {
211     cout << "==> TFluka::FinishGeometry() called." << endl;
212     printf("----FinishGeometry - nothing to do with TGeo\n");
213     cout << "<== TFluka::FinishGeometry() called." << endl;
214   }  
215
216
217 //______________________________________________________________________________ 
218 void TFluka::BuildPhysics() {
219 //
220 //  Prepare FLUKA input files and call FLUKA physics initialisation
221 //
222     
223     if (fVerbosityLevel >=3)
224         cout << "==> TFluka::BuildPhysics() called." << endl;
225 // Prepare input file with the current physics settings
226     InitPhysics(); 
227     cout << "\t* InitPhysics() - Prepare input file was called" << endl; 
228     
229     if (fVerbosityLevel >=2)
230         cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
231              << ") in fluka..." << endl;
232     GLOBAL.lfdrtr = true;
233     
234     if (fVerbosityLevel >=2)
235         cout << "\t* Opening file " << fInputFileName << endl;
236     const char* fname = fInputFileName;
237     fluka_openinp(lunin, PASSCHARA(fname));
238     
239     if (fVerbosityLevel >=2)
240         cout << "\t* Calling flukam..." << endl;
241     flukam(1);
242     
243     if (fVerbosityLevel >=2)
244         cout << "\t* Closing file " << fInputFileName << endl;
245     fluka_closeinp(lunin);
246     
247     FinishGeometry();
248     
249     if (fVerbosityLevel >=3)
250         cout << "<== TFluka::Init() called." << endl;
251     
252     
253     if (fVerbosityLevel >=3)
254         cout << "<== TFluka::BuildPhysics() called." << endl;
255 }  
256
257 //______________________________________________________________________________ 
258 void TFluka::ProcessEvent() {
259 //
260 // Process one event
261 //
262     if (fStopRun) {
263         printf("User Run Abortion: No more events handled !\n");
264         fNEvent += 1;
265         return;
266     }
267
268     if (fVerbosityLevel >=3)
269         cout << "==> TFluka::ProcessEvent() called." << endl;
270     fApplication->GeneratePrimaries();
271     EPISOR.lsouit = true;
272     flukam(1);
273     if (fVerbosityLevel >=3)
274         cout << "<== TFluka::ProcessEvent() called." << endl;
275     //
276     // Increase event number
277     //
278     fNEvent += 1;
279 }
280
281 //______________________________________________________________________________ 
282 Bool_t TFluka::ProcessRun(Int_t nevent) {
283 //
284 // Run steering
285 //
286
287   if (fVerbosityLevel >=3)
288     cout << "==> TFluka::ProcessRun(" << nevent << ") called." 
289          << endl;
290
291   if (fVerbosityLevel >=2) {
292     cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
293     cout << "\t* Calling flukam again..." << endl;
294   }
295
296   fApplication->InitGeometry();
297   Int_t todo = TMath::Abs(nevent);
298   for (Int_t ev = 0; ev < todo; ev++) {
299       fApplication->BeginEvent();
300       ProcessEvent();
301       fApplication->FinishEvent();
302   }
303
304   if (fVerbosityLevel >=3)
305     cout << "<== TFluka::ProcessRun(" << nevent << ") called." 
306          << endl;
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   TGeoMaterial *mat;
375   kmat = gGeoManager->GetListOfMaterials()->GetSize();
376   if ((z-Int_t(z)) > 1E-3) {
377      mat = fGeom->GetMakeWrongMaterial(z);
378      if (mat) {
379         mat->SetRadLen(radl,absl);
380         mat->SetUniqueID(kmat);
381         return;
382      }
383   }      
384   gGeoManager->Material(name, a, z, dens, kmat, radl, absl);
385
386
387 //______________________________________________________________________________ 
388 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a, 
389                      Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
390 //
391   Double_t* da = fGeom->CreateDoubleArray(a, TMath::Abs(nlmat));  
392   Double_t* dz = fGeom->CreateDoubleArray(z, TMath::Abs(nlmat));  
393   Double_t* dwmat = fGeom->CreateDoubleArray(wmat, TMath::Abs(nlmat));  
394
395   Mixture(kmat, name, da, dz, dens, nlmat, dwmat);
396   for (Int_t i=0; i<nlmat; i++) {
397     a[i] = da[i]; z[i] = dz[i]; wmat[i] = dwmat[i];
398   }  
399
400   delete [] da;
401   delete [] dz;
402   delete [] dwmat;
403
404
405 //______________________________________________________________________________ 
406 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a, 
407                      Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
408 //
409   // Defines mixture OR COMPOUND IMAT as composed by 
410   // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
411   // 
412   // If NLMAT > 0 then wmat contains the proportion by
413   // weights of each basic material in the mixture. 
414   // 
415   // If nlmat < 0 then WMAT contains the number of atoms 
416   // of a given kind into the molecule of the COMPOUND
417   // In this case, WMAT in output is changed to relative
418   // weigths.
419   //
420   Int_t i,j;
421   if (nlmat < 0) {
422      nlmat = - nlmat;
423      Double_t amol = 0;
424      for (i=0;i<nlmat;i++) {
425         amol += a[i]*wmat[i];
426      }
427      for (i=0;i<nlmat;i++) {
428         wmat[i] *= a[i]/amol;
429      }
430   }
431   kmat = gGeoManager->GetListOfMaterials()->GetSize();
432   // Check if we have elements with fractional Z
433   TGeoMaterial *mat = 0;
434   TGeoMixture *mix = 0;
435   Bool_t mixnew = kFALSE;
436   for (i=0; i<nlmat; i++) {
437      if (z[i]-Int_t(z[i]) < 1E-3) continue;
438      // We have found an element with fractional Z -> loop mixtures to look for it
439      for (j=0; j<kmat; j++) {
440         mat = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->At(j);
441         if (!mat) break;
442         if (!mat->IsMixture()) continue;
443         mix = (TGeoMixture*)mat;
444         if (TMath::Abs(z[i]-mix->GetZ()) >1E-3) continue;
445 //        printf(" FOUND component %i as mixture %s\n", i, mat->GetName());
446         mixnew = kTRUE;
447         break;
448      }
449      if (!mixnew) Warning("Mixture","%s : cannot find component %i with fractional Z=%f\n", name, i, z[i]);
450      break;
451   }   
452   if (mixnew) {
453      Int_t nlmatnew = nlmat+mix->GetNelements()-1;
454      Double_t *anew = new Double_t[nlmatnew];
455      Double_t *znew = new Double_t[nlmatnew];
456      Double_t *wmatnew = new Double_t[nlmatnew];
457      Int_t ind=0;
458      for (j=0; j<nlmat; j++) {
459         if (j==i) continue;
460         anew[ind] = a[j];
461         znew[ind] = z[j];
462         wmatnew[ind] = wmat[j];
463         ind++;
464      }
465      for (j=0; j<mix->GetNelements(); j++) {
466         anew[ind] = mix->GetAmixt()[j];
467         znew[ind] = mix->GetZmixt()[j];
468         wmatnew[ind] = wmat[i]*mix->GetWmixt()[j];
469         ind++;
470      }
471      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
472      delete [] anew;
473      delete [] znew;
474      delete [] wmatnew;
475      return;
476   }   
477   // Now we need to compact identical elements within the mixture
478   // First check if this happens   
479   mixnew = kFALSE;  
480   for (i=0; i<nlmat-1; i++) {
481      for (j=i+1; j<nlmat; j++) {
482         if (z[i] == z[j]) {
483            mixnew = kTRUE;
484            break;
485         }
486      }   
487      if (mixnew) break;
488   }   
489   if (mixnew) {
490      Int_t nlmatnew = 0;
491      Double_t *anew = new Double_t[nlmat];
492      Double_t *znew = new Double_t[nlmat];
493      memset(znew, 0, nlmat*sizeof(Double_t));
494      Double_t *wmatnew = new Double_t[nlmat];
495      Bool_t skipi;
496      for (i=0; i<nlmat; i++) {
497         skipi = kFALSE;
498         for (j=0; j<nlmatnew; j++) {
499            if (z[i] == z[j]) {
500               wmatnew[j] += wmat[i];
501               skipi = kTRUE;
502               break;
503            }
504         }   
505         if (skipi) continue;    
506         anew[nlmatnew] = a[i];
507         znew[nlmatnew] = z[i];
508         wmatnew[nlmatnew] = wmat[i];
509         nlmatnew++;
510      }
511      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
512      delete [] anew;
513      delete [] znew;
514      delete [] wmatnew;
515      return;     
516    }
517    gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
518
519
520 //______________________________________________________________________________ 
521 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, 
522                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, 
523                     Double_t stemax, Double_t deemax, Double_t epsil, 
524                     Double_t stmin, Float_t* ubuf, Int_t nbuf) { 
525   //
526   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
527   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, 
528              epsil, stmin, ubuf, nbuf);
529
530
531 //______________________________________________________________________________ 
532 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, 
533                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, 
534                     Double_t stemax, Double_t deemax, Double_t epsil, 
535                     Double_t stmin, Double_t* ubuf, Int_t nbuf) { 
536   //
537   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
538   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, 
539              epsil, stmin, ubuf, nbuf);
540
541
542 //______________________________________________________________________________ 
543 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX, 
544                     Double_t thetaY, Double_t phiY, Double_t thetaZ, 
545                     Double_t phiZ) {
546 //                   
547   krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
548   fMCGeo->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ); 
549
550
551 //______________________________________________________________________________ 
552 void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
553 //
554 //
555     
556    if (fVerbosityLevel >=3) printf("Gstpar called with %6d %5s %12.4e %6d\n", itmed, param, parval, fGeom->GetFlukaMaterial(itmed));
557  
558    Bool_t process = kFALSE;
559    if (strncmp(param, "DCAY",  4) == 0 ||
560        strncmp(param, "PAIR",  4) == 0 ||
561        strncmp(param, "COMP",  4) == 0 ||
562        strncmp(param, "PHOT",  4) == 0 ||
563        strncmp(param, "PFIS",  4) == 0 ||
564        strncmp(param, "DRAY",  4) == 0 ||
565        strncmp(param, "ANNI",  4) == 0 ||
566        strncmp(param, "BREM",  4) == 0 ||
567        strncmp(param, "MUNU",  4) == 0 ||
568        strncmp(param, "CKOV",  4) == 0 ||
569        strncmp(param, "HADR",  4) == 0 ||
570        strncmp(param, "LOSS",  4) == 0 ||
571        strncmp(param, "MULS",  4) == 0 ||
572        strncmp(param, "RAYL",  4) == 0) 
573    {
574        process = kTRUE;
575    } 
576    if (process) {
577        SetProcess(param, Int_t (parval), fGeom->GetFlukaMaterial(itmed));
578    } else {
579        SetCut(param, parval, fGeom->GetFlukaMaterial(itmed));
580    }
581 }    
582
583 // functions from GGEOM 
584 //_____________________________________________________________________________
585 void TFluka::Gsatt(const char *name, const char *att, Int_t val)
586
587   // Set visualisation attributes for one volume
588   char vname[5];
589   fGeom->Vname(name,vname);
590   char vatt[5];
591   fGeom->Vname(att,vatt);
592   gGeoManager->SetVolumeAttribute(vname, vatt, val);
593 }
594
595 //______________________________________________________________________________ 
596 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,  
597                      Float_t *upar, Int_t np)  {
598 //
599     return fMCGeo->Gsvolu(name, shape, nmed, upar, np); 
600 }
601
602 //______________________________________________________________________________ 
603 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,  
604                      Double_t *upar, Int_t np)  {
605 //
606     return fMCGeo->Gsvolu(name, shape, nmed, upar, np); 
607 }
608  
609 //______________________________________________________________________________ 
610 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv, 
611                    Int_t iaxis) {
612 //
613     fMCGeo->Gsdvn(name, mother, ndiv, iaxis); 
614
615
616 //______________________________________________________________________________ 
617 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv, 
618                     Int_t iaxis, Double_t c0i, Int_t numed) {
619 //
620     fMCGeo->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed); 
621
622
623 //______________________________________________________________________________ 
624 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step, 
625                    Int_t iaxis, Int_t numed, Int_t ndvmx) {
626 //      
627     fMCGeo->Gsdvt(name, mother, step, iaxis, numed, ndvmx); 
628
629
630 //______________________________________________________________________________ 
631 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step, 
632                     Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) { 
633 //
634     fMCGeo->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx); 
635
636
637 //______________________________________________________________________________ 
638 void TFluka::Gsord(const char * /*name*/, Int_t /*iax*/) {
639 //
640 // Nothing to do with TGeo
641
642
643 //______________________________________________________________________________ 
644 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,  
645                    Double_t x, Double_t y, Double_t z, Int_t irot, 
646                    const char *konly) {
647 //
648   fMCGeo->Gspos(name, nr, mother, x, y, z, irot, konly); 
649
650
651 //______________________________________________________________________________ 
652 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,  
653                     Double_t x, Double_t y, Double_t z, Int_t irot,
654                     const char *konly, Float_t *upar, Int_t np)  {
655   //
656   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); 
657
658
659 //______________________________________________________________________________ 
660 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,  
661                     Double_t x, Double_t y, Double_t z, Int_t irot,
662                     const char *konly, Double_t *upar, Int_t np)  {
663   //
664   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); 
665
666
667 //______________________________________________________________________________ 
668 void TFluka::Gsbool(const char* /*onlyVolName*/, const char* /*manyVolName*/) {
669 //
670 // Nothing to do with TGeo
671 }
672
673 //______________________________________________________________________________ 
674 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
675                          Float_t* absco, Float_t* effic, Float_t* rindex) {
676 //
677 // Set Cerenkov properties for medium itmed
678 //
679 // npckov: number of sampling points
680 // ppckov: energy values
681 // absco:  absorption length
682 // effic:  quantum efficiency
683 // rindex: refraction index
684 //
685 //
686 //  
687 //  Create object holding Cerenkov properties
688 //  
689     TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex);
690 //
691 //  Pass object to medium
692     TGeoMedium* medium = gGeoManager->GetMedium(itmed);
693     medium->SetCerenkovProperties(cerenkovProperties);
694 }  
695
696 //______________________________________________________________________________ 
697 void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t * /*ppckov*/,
698                          Double_t * /*absco*/, Double_t * /*effic*/, Double_t * /*rindex*/) {
699 //
700 // Not implemented with TGeo - what G4 did ? Any FLUKA card generated?
701    Warning("SetCerenkov", "Not implemented with TGeo");
702 }  
703     
704 // Euclid
705 //______________________________________________________________________________ 
706 void TFluka::WriteEuclid(const char* /*fileName*/, const char* /*topVol*/, 
707                           Int_t /*number*/, Int_t /*nlevel*/) {
708 //
709 // Not with TGeo
710    Warning("WriteEuclid", "Not implemented with TGeo");
711
712
713
714
715 //_____________________________________________________________________________
716 // methods needed by the stepping
717 //____________________________________________________________________________ 
718
719 Int_t TFluka::GetMedium() const {
720 //
721 //  Get the medium number for the current fluka region
722 //
723     return fGeom->GetMedium(); // this I need to check due to remapping !!!
724 }
725
726
727
728 //____________________________________________________________________________ 
729 // particle table usage
730 // ID <--> PDG transformations
731 //_____________________________________________________________________________
732 Int_t TFluka::IdFromPDG(Int_t pdg) const 
733 {
734     //
735     // Return Fluka code from PDG and pseudo ENDF code
736     
737     // Catch the feedback photons
738     if (pdg == 50000051) return (-1);
739     // MCIHAD() goes from pdg to fluka internal.
740     Int_t intfluka = mcihad(pdg);
741     // KPTOIP array goes from internal to official
742     return GetFlukaKPTOIP(intfluka);
743 }
744
745 //______________________________________________________________________________ 
746 Int_t TFluka::PDGFromId(Int_t id) const 
747 {
748   //
749   // Return PDG code and pseudo ENDF code from Fluka code
750   //                      Alpha     He3       Triton    Deuteron  gen. ion  opt. photon   
751     Int_t idSpecial[6] = {10020040, 10020030, 10010030, 10010020, 10000000, 50000050};
752   // IPTOKP array goes from official to internal
753
754     if (id == -1) {
755 // Cerenkov photon
756         if (fVerbosityLevel >= 3)
757             printf("\n PDGFromId: Cerenkov Photon \n");
758         return  50000050;
759     }
760 // Error id    
761     if (id == 0 || id < -6 || id > 250) {
762         if (fVerbosityLevel >= 3)
763             printf("PDGFromId: Error id = 0\n");
764         return -1;
765     }
766 // Good id    
767     if (id > 0) {
768         Int_t intfluka = GetFlukaIPTOKP(id);
769         if (intfluka == 0) {
770             if (fVerbosityLevel >= 3)
771                 printf("PDGFromId: Error intfluka = 0: %d\n", id);
772             return -1;
773         } else if (intfluka < 0) {
774             if (fVerbosityLevel >= 3)
775                 printf("PDGFromId: Error intfluka < 0: %d\n", id);
776             return -1;
777         }
778         if (fVerbosityLevel >= 3)
779             printf("mpdgha called with %d %d \n", id, intfluka);
780         // MPDGHA() goes from fluka internal to pdg.
781         return mpdgha(intfluka);
782     } else {
783         // ions and optical photons
784         return idSpecial[id + 6];
785     }
786 }
787
788 void TFluka::StopTrack()
789 {
790     // Set stopping conditions
791     // Works for photons and charged particles
792     fStopped = kTRUE;
793 }
794   
795 //_____________________________________________________________________________
796 // methods for physics management
797 //____________________________________________________________________________ 
798 //
799 // set methods
800 //
801
802 void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imed)
803 {
804 //  Set process user flag for material imat
805 //
806     TFlukaConfigOption* proc = new TFlukaConfigOption(flagName, flagValue, imed);
807     fProcesses->Add(proc);
808 }
809
810 //______________________________________________________________________________ 
811 Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue)
812 {
813 //  Set process user flag 
814 //
815 //    
816 //  Update if already in the list
817 //
818
819     TIter next(fProcesses);
820     TFlukaConfigOption* proc;
821     while((proc = (TFlukaConfigOption*)next()))
822     { 
823         if (strcmp(proc->GetName(), flagName) == 0) {
824             proc->SetFlag(flagValue);
825             proc->SetMedium(-1);
826             return kTRUE;
827          }
828     }
829 //
830 // If not create a new process
831 //    
832
833     proc = new TFlukaConfigOption(flagName, flagValue);
834     fProcesses->Add(proc);
835     
836     return kTRUE;  
837 }
838
839 //______________________________________________________________________________ 
840 void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
841 {
842 // Set user cut value for material imed
843 //
844     TFlukaConfigOption* cut = new TFlukaConfigOption(cutName, cutValue, imed);
845     fCuts->Add(cut);
846 }
847
848 //______________________________________________________________________________ 
849 Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue)
850 {
851 // Set user cut value 
852 //
853 //    
854 //  Update if already in the list
855 //
856
857     TIter next(fCuts);
858     TFlukaConfigOption* cut;
859     while((cut = (TFlukaConfigOption*)next()))
860     { 
861         if (strcmp(cut->GetName(), cutName) == 0) {
862             cut->SetCut(cutValue);
863             return kTRUE;
864          }
865     }
866 //
867 // If not create a new process
868 //    
869
870     cut = new TFlukaConfigOption(cutName, cutValue);
871     fCuts->Add(cut);
872     
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   Int_t j, k;
901   Double_t fCut;
902
903   FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
904
905   Double_t zero  = 0.0;
906   Double_t one   = 1.0;
907   Double_t two   = 2.0;
908   Double_t three = 3.0;
909
910   Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
911   if (fVerbosityLevel >= 3) printf("   last FLUKA material is %g\n", fLastMaterial);
912
913   // Prepare  Cerenkov
914   TObjArray *matList = GetFlukaMaterials();
915   Int_t nmaterial =  matList->GetEntriesFast();
916   fMaterials = new Int_t[nmaterial+3];
917               
918 // construct file names
919
920   TString sFlukaVmcCoreInp = getenv("ALICE_ROOT");
921   sFlukaVmcCoreInp +="/TFluka/input/";
922   TString sFlukaVmcTmp = "flukaMat.inp";
923   TString sFlukaVmcInp = GetInputFileName();
924   sFlukaVmcCoreInp += GetCoreInputFileName();
925
926 // open files 
927
928   if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
929       printf("\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
930       exit(1);
931   }
932   if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
933       printf("\nCannot open file %s\n",sFlukaVmcTmp.Data());
934       exit(1);
935   }
936   if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
937       printf("\nCannot open file %s\n",sFlukaVmcInp.Data());
938       exit(1);
939   }
940
941 // copy core input file 
942   Char_t sLine[255];
943   Float_t fEventsPerRun;
944   
945   while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
946       if (strncmp(sLine,"GEOEND",6) != 0)
947           fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
948       else {
949           fprintf(pFlukaVmcInp,"GEOEND\n");   // add GEOEND card
950           goto flukamat;
951       }
952   } // end of while until GEOEND card
953   
954
955  flukamat:
956   while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
957       fprintf(pFlukaVmcInp,"%s\n",sLine);
958   }
959   
960   while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) { 
961       if (strncmp(sLine,"START",5) != 0)
962           fprintf(pFlukaVmcInp,"%s\n",sLine);
963       else {
964           sscanf(sLine+10,"%10f",&fEventsPerRun);
965       goto fin;
966       }
967   } //end of while until START card
968   
969 fin:
970 // in G3 the process control values meaning can be different for
971 // different processes, but for most of them is:
972 //   0  process is not activated
973 //   1  process is activated WITH generation of secondaries
974 //   2  process is activated WITHOUT generation of secondaries
975 // if process does not generate secondaries => 1 same as 2
976 //
977 // Exceptions:
978 //   MULS:  also 3
979 //   LOSS:  also 3, 4
980 //   RAYL:  only 0,1
981 //   HADR:  may be > 2
982 //
983  
984 // Loop over number of SetProcess calls 
985   fprintf(pFlukaVmcInp,"*----------------------------------------------------------------------------- \n");
986   fprintf(pFlukaVmcInp,"*----- The following data are generated from SetProcess and SetCut calls ----- \n");
987   fprintf(pFlukaVmcInp,"*----------------------------------------------------------------------------- \n");
988
989 // Outer loop over processes
990   TIter next(fProcesses);
991   TFlukaConfigOption *proc;
992 // Inner loop over processes
993   TIter nextp(fProcesses);
994   TFlukaConfigOption *procp;
995 // Loop over cuts
996   TIter nextc(fCuts);
997   TFlukaConfigOption *cut = 0x0;
998
999   while((proc = (TFlukaConfigOption*)next())) {
1000       Float_t matMin = three;
1001       Float_t matMax = fLastMaterial;
1002       Bool_t  global = kTRUE;
1003       if (proc->Medium() != -1) {
1004           matMin = Float_t(proc->Medium());
1005           matMax = matMin;
1006           global = kFALSE;
1007       }
1008       
1009     // annihilation
1010     // G3 default value: 1
1011     // G4 processes: G4eplusAnnihilation/G4IeplusAnnihilation
1012     // Particles: e+
1013     // Physics:   EM
1014     // flag = 0 no annihilation
1015     // flag = 1 annihilation, decays processed
1016     // flag = 2 annihilation, no decay product stored
1017     // gMC ->SetProcess("ANNI",1); // EMFCUT   -1.   0.  0. 3. lastmat 0. ANNH-THR
1018       if (strncmp(proc->GetName(),"ANNI",4) == 0) {
1019           if (proc->Flag() == 1 || proc->Flag() == 2) {
1020               fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for e+ annihilation - resets to default=0.\n");
1021               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('ANNI',1) or SetProcess('ANNI',2)\n");
1022               // -one = kinetic energy threshold (GeV) for e+ annihilation (resets to default=0)
1023               // zero = not used
1024               // zero = not used
1025               // matMin = lower bound of the material indices in which the respective thresholds apply
1026               // matMax = upper bound of the material indices in which the respective thresholds apply
1027               // one = step length in assigning indices
1028               // "ANNH-THR"; 
1029               fprintf(pFlukaVmcInp,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",-one,zero,zero,matMin,matMax,one);
1030           }
1031           else if (proc->Flag() == 0) {
1032               fprintf(pFlukaVmcInp,"*\n*No annihilation - no FLUKA card generated\n");
1033               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('ANNI',0)\n");
1034           }
1035           else  {
1036               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('ANNI',?) call.\n");
1037               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1038           }
1039       }
1040     
1041     // bremsstrahlung and pair production are both activated
1042     // G3 default value: 1
1043     // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung,
1044     //               G4MuBremsstrahlung/G4IMuBremsstrahlung,
1045     //               G4LowEnergyBremstrahlung
1046     // Particles: e-/e+; mu+/mu-
1047     // Physics:   EM
1048     // flag = 0 no bremsstrahlung
1049     // flag = 1 bremsstrahlung, photon processed
1050     // flag = 2 bremsstrahlung, no photon stored
1051     // gMC ->SetProcess("BREM",1); // PAIRBREM  2.   0.  0. 3. lastmat
1052                                  // EMFCUT   -1.   0.  0. 3. lastmat 0. ELPO-THR
1053     // G3 default value: 1
1054     // G4 processes: G4GammaConversion,
1055     //               G4MuPairProduction/G4IMuPairProduction
1056     //               G4LowEnergyGammaConversion
1057     // Particles: gamma, mu
1058     // Physics:   EM
1059     // flag = 0 no delta rays
1060     // flag = 1 delta rays, secondaries processed
1061     // flag = 2 delta rays, no secondaries stored
1062     // gMC ->SetProcess("PAIR",1); // PAIRBREM  1.   0.  0. 3. lastmat
1063                                  // EMFCUT    0.   0. -1. 3. lastmat 0. PHOT-THR
1064     else if ((strncmp(proc->GetName(),"PAIR",4) == 0) && (proc->Flag() == 1 || proc->Flag() == 2)) {
1065
1066         nextp.Reset();
1067         
1068         while ((procp = (TFlukaConfigOption*)nextp())) {
1069             if ((strncmp(procp->GetName(),"BREM",4) == 0) && 
1070                 (proc->Flag() == 1 || procp->Flag() == 2) &&
1071                 (procp->Medium() == proc->Medium())) {
1072                 fprintf(pFlukaVmcInp,"*\n*Bremsstrahlung and pair production by muons and charged hadrons both activated\n");
1073                 fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1) and SetProcess('PAIR',1)\n");
1074                 fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.\n");
1075                 fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.\n");
1076                 // three = bremsstrahlung and pair production by muons and charged hadrons both are activated
1077                 fprintf(pFlukaVmcInp,"PAIRBREM  %10.1f",three);
1078                 // direct pair production by muons
1079                 // G4 particles: "e-", "e+"
1080                 // G3 default value: 0.01 GeV
1081                 //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons
1082                 fCut = 0.0;
1083                 nextc.Reset();
1084                 while ((cut = (TFlukaConfigOption*)nextc())) {
1085                     if (strncmp(cut->GetName(), "PPCUTM", 6) == 0 &&
1086                         (cut->Medium() == proc->Medium())) fCut = cut->Cut();
1087                 }
1088                 fprintf(pFlukaVmcInp,"%10.4g",fCut);
1089                 // fCut; = e+, e- kinetic energy threshold (in GeV) for explicit pair production.
1090                 // muon and hadron bremsstrahlung
1091                 // G4 particles: "gamma"
1092                 // G3 default value: CUTGAM=0.001 GeV
1093                 //gMC ->SetCut("BCUTM",cut);  // cut for muon and hadron bremsstrahlung
1094                 fCut = 0.0;
1095                 nextc.Reset();
1096                 while ((cut = (TFlukaConfigOption*)nextc())) {
1097                     if (strncmp(cut->GetName(), "BCUTM", 5) == 0 &&
1098                         (cut->Medium() == proc->Medium())) fCut = cut->Cut();
1099                 }
1100                 fprintf(pFlukaVmcInp,"%10.4g%10.1f%10.1f\n",fCut,matMin,matMax);
1101                 // fCut = photon energy threshold (GeV) for explicit bremsstrahlung production
1102                 // matMin = lower bound of the material indices in which the respective thresholds apply
1103                 // matMax = upper bound of the material indices in which the respective thresholds apply
1104                 
1105                 // for e+ and e-
1106                 fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n");
1107                 fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1);\n");
1108                 fCut = -1.0;
1109                 nextc.Reset();
1110                 while ((cut = (TFlukaConfigOption*)nextc())) {
1111                     if (strncmp(cut->GetName(), "BCUTE", 5) == 0 &&
1112                         (cut->Medium() == proc->Medium())) fCut = cut->Cut();
1113                 }
1114                 //fCut = kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
1115                 // zero = not used
1116                 // zero = not used
1117                 // matMin = lower bound of the material indices in which the respective thresholds apply
1118                 // matMax = upper bound of the material indices in which the respective thresholds apply
1119                 // one = step length in assigning indices
1120                 // "ELPO-THR"; 
1121                 fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",fCut,zero,zero,matMin,matMax,one);
1122                 
1123           // for e+ and e-
1124                 fprintf(pFlukaVmcInp,"*\n*Pair production by electrons is activated\n");
1125                 fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PAIR',1);\n");
1126                 fCut = -1.0;
1127                 nextc.Reset();
1128                 while ((cut = (TFlukaConfigOption*)nextc())) {
1129                     if (strncmp(cut->GetName(), "CUTGAM", 6) == 0 &&
1130                         (cut->Medium() == proc->Medium())) fCut = cut->Cut();
1131                 }
1132                 // fCut = energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored)
1133                 // matMin = lower bound of the material indices in which the respective thresholds apply
1134                 // matMax =  upper bound of the material indices in which the respective thresholds apply
1135                 // one = step length in assigning indices
1136                 fprintf(pFlukaVmcInp,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,fCut,matMin,matMax,one);
1137                 goto BOTH;
1138             } // end of if for BREM
1139         } // end of loop for BREM
1140         
1141         // only pair production by muons and charged hadrons is activated
1142         fprintf(pFlukaVmcInp,"*\n*Pair production by muons and charged hadrons is activated\n");
1143         fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n");
1144         fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.\n");
1145         // direct pair production by muons
1146         // G4 particles: "e-", "e+"
1147         // G3 default value: 0.01 GeV
1148         //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons
1149         // one = pair production by muons and charged hadrons is activated
1150         // zero = e+, e- kinetic energy threshold (in GeV) for explicit pair production.
1151         // zero = no explicit bremsstrahlung production is simulated
1152         // matMin = lower bound of the material indices in which the respective thresholds apply
1153         // matMax = upper bound of the material indices in which the respective thresholds apply
1154         fprintf(pFlukaVmcInp,"PAIRBREM  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax);
1155         
1156         // for e+ and e-
1157         fprintf(pFlukaVmcInp,"*\n*Pair production by electrons is activated\n");
1158         fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n");
1159         fCut = -1.0;
1160         nextc.Reset();
1161         while ((cut = (TFlukaConfigOption*)nextc())) {
1162             if (strncmp(cut->GetName(), "CUTGAM", 6) == 0 &&
1163                 (cut->Medium() == proc->Medium())) fCut = cut->Cut();
1164         }
1165         // zero = energy threshold (GeV) for Compton scattering (= 0.0 : ignored)
1166         // zero = energy threshold (GeV) for Photoelectric (= 0.0 : ignored)
1167         // fCut = energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored)
1168         // matMin = lower bound of the material indices in which the respective thresholds apply
1169         // matMax = upper bound of the material indices in which the respective thresholds apply
1170         // one = step length in assigning indices
1171         fprintf(pFlukaVmcInp,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,fCut,matMin,matMax,one);
1172       
1173     BOTH:
1174         k = 0;
1175     } // end of if for PAIR
1176       
1177       
1178       
1179       // bremsstrahlung
1180       // G3 default value: 1
1181       // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung,
1182       //               G4MuBremsstrahlung/G4IMuBremsstrahlung,
1183       //               G4LowEnergyBremstrahlung
1184       // Particles: e-/e+; mu+/mu-
1185       // Physics:   EM
1186       // flag = 0 no bremsstrahlung
1187       // flag = 1 bremsstrahlung, photon processed
1188       // flag = 2 bremsstrahlung, no photon stored
1189       // gMC ->SetProcess("BREM",1); // PAIRBREM  2.   0.  0. 3. lastmat
1190       // EMFCUT   -1.   0.  0. 3. lastmat 0. ELPO-THR
1191       else if (strncmp(proc->GetName(),"BREM",4) == 0) {
1192           nextp.Reset();
1193           while((procp = (TFlukaConfigOption*)nextp())) {
1194               if ((strncmp(procp->GetName(),"PAIR",4) == 0) && 
1195                   procp->Flag() == 1 &&
1196                   (procp->Medium() == proc->Medium())) goto NOBREM;
1197           }
1198           if (proc->Flag() == 1 || proc->Flag() == 2) { 
1199               fprintf(pFlukaVmcInp,"*\n*Bremsstrahlung by muons and charged hadrons is activated\n");
1200               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1) or SetProcess('BREM',2)\n");
1201               fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.\n");
1202               // two = bremsstrahlung by muons and charged hadrons is activated
1203               // zero = no meaning
1204               // muon and hadron bremsstrahlung
1205               // G4 particles: "gamma"
1206               // G3 default value: CUTGAM=0.001 GeV
1207               //gMC ->SetCut("BCUTM",cut);  // cut for muon and hadron bremsstrahlung
1208               fCut = 0.0;
1209               nextc.Reset();
1210               while ((cut = (TFlukaConfigOption*)nextc())) {
1211                   if (strncmp(cut->GetName(), "BCUTM", 5) == 0 &&
1212                       (cut->Medium() == proc->Medium())) fCut = cut->Cut();
1213               }
1214               // fCut = photon energy threshold (GeV) for explicit bremsstrahlung production
1215               // matMin = lower bound of the material indices in which the respective thresholds apply
1216               // matMax = upper bound of the material indices in which the respective thresholds apply
1217               fprintf(pFlukaVmcInp,"PAIRBREM  %10.1f%10.1f%10.4g%10.1f%10.1f\n",two,zero,fCut,matMin,matMax);
1218               
1219               // for e+ and e-
1220               fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n");
1221               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1);");
1222               // - one = kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
1223               // zero = not used
1224               // zero = not used
1225               // matMin = lower bound of the material indices in which the respective thresholds apply
1226               // matMax = upper bound of the material indices in which the respective thresholds apply
1227               // one = step length in assigning indices
1228               //"ELPO-THR"; 
1229               fprintf(pFlukaVmcInp,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",-one,zero,zero,matMin,matMax,one);
1230           }
1231           else if (proc->Flag() == 0) {
1232               fprintf(pFlukaVmcInp,"*\n*No bremsstrahlung - no FLUKA card generated\n");
1233               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',0)\n");
1234           }
1235           else  {
1236               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('BREM',?) call.\n");
1237               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1238           }
1239       NOBREM:
1240           j = 0;
1241       } // end of else if (strncmp(proc->GetName(),"BREM",4) == 0)
1242       
1243       // Cerenkov photon generation
1244       // G3 default value: 0
1245       // G4 process: G4Cerenkov
1246       // 
1247       // Particles: charged
1248       // Physics:   Optical
1249       // flag = 0 no Cerenkov photon generation
1250       // flag = 1 Cerenkov photon generation
1251       // flag = 2 Cerenkov photon generation with primary stopped at each step
1252       //xx gMC ->SetProcess("CKOV",1); // ??? Cerenkov photon generation
1253       
1254       else if (strncmp(proc->GetName(),"CKOV",4) == 0) {
1255           if ((proc->Flag() == 1 || proc->Flag() == 2) && global) {
1256               // Write comments
1257               fprintf(pFlukaVmcInp, "* \n"); 
1258               fprintf(pFlukaVmcInp, "*Cerenkov photon generation\n"); 
1259               fprintf(pFlukaVmcInp, "*Generated from call: SetProcess('CKOV',1) or SetProcess('CKOV',2)\n"); 
1260               // Loop over media 
1261               for (Int_t im = 0; im < nmaterial; im++)
1262               {
1263                   TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
1264                   Int_t idmat = material->GetIndex();
1265
1266                   if (!global && idmat != proc->Medium()) continue;
1267                   
1268                   fMaterials[idmat] = im;
1269                   // Skip media with no Cerenkov properties
1270                   TFlukaCerenkov* cerenkovProp;
1271                   if (!(cerenkovProp = dynamic_cast<TFlukaCerenkov*>(material->GetCerenkovProperties()))) continue;
1272                   //
1273                   // This medium has Cerenkov properties 
1274                   //
1275                   //
1276                   // Write OPT-PROD card for each medium 
1277                   Float_t  emin  = cerenkovProp->GetMinimumEnergy();
1278                   Float_t  emax  = cerenkovProp->GetMaximumEnergy();          
1279                   fprintf(pFlukaVmcInp, "OPT-PROD  %10.4g%10.4g%10.4g%10.4g%10.4g%10.4gCERENKOV\n", emin, emax, 0., 
1280                           Float_t(idmat), Float_t(idmat), 0.); 
1281                   //
1282                   // Write OPT-PROP card for each medium 
1283                   // Forcing FLUKA to call user routines (queffc.cxx, rflctv.cxx, rfrndx.cxx)
1284                   //
1285                   fprintf(pFlukaVmcInp, "OPT-PROP  %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fWV-LIMIT\n",  
1286                           cerenkovProp->GetMinimumWavelength(),
1287                           cerenkovProp->GetMaximumWavelength(), 
1288                           cerenkovProp->GetMaximumWavelength(), 
1289                           Float_t(idmat), Float_t(idmat), 0.0);
1290                   
1291                   if (cerenkovProp->IsMetal()) {
1292                       fprintf(pFlukaVmcInp, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fMETAL\n",  
1293                               -100., -100., -100., 
1294                               Float_t(idmat), Float_t(idmat), 0.0);
1295                   } else {
1296                       fprintf(pFlukaVmcInp, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",  
1297                               -100., -100., -100., 
1298                               Float_t(idmat), Float_t(idmat), 0.0);
1299                   }
1300                   
1301                   
1302                   for (Int_t j = 0; j < 3; j++) {
1303                       fprintf(pFlukaVmcInp, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f&\n",  
1304                               -100., -100., -100., 
1305                               Float_t(idmat), Float_t(idmat), 0.0);
1306                   }
1307                   // Photon detection efficiency user defined
1308                   
1309                   if (cerenkovProp->IsSensitive())
1310                       fprintf(pFlukaVmcInp, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fSENSITIV\n",  
1311                               -100., -100., -100., 
1312                               Float_t(idmat), Float_t(idmat), 0.0);
1313                   
1314               } // materials
1315           } else if (proc->Flag() == 0) {
1316               fprintf(pFlukaVmcInp,"*\n*No Cerenkov photon generation\n");
1317               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('CKOV',0)\n");
1318               // zero = not used
1319               // zero = not used
1320               // zero = not used
1321               // matMin = lower bound of the material indices in which the respective thresholds apply
1322               // matMax = upper bound of the material indices in which the respective thresholds apply
1323               // one = step length in assigning indices
1324               //"CERE-OFF"; 
1325               fprintf(pFlukaVmcInp,"OPT-PROD  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",zero,zero,zero,matMin,matMax,one);
1326           }
1327           else  {
1328               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('CKOV',?) call.\n");
1329               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1330           }
1331       } // end of else if (strncmp(proc->GetName(),"CKOV",4) == 0)
1332       
1333       // Compton scattering
1334       // G3 default value: 1
1335       // G4 processes: G4ComptonScattering,
1336       //               G4LowEnergyCompton,
1337       //               G4PolarizedComptonScattering
1338       // Particles: gamma
1339       // Physics:   EM
1340       // flag = 0 no Compton scattering
1341       // flag = 1 Compton scattering, electron processed
1342       // flag = 2 Compton scattering, no electron stored
1343       // gMC ->SetProcess("COMP",1); // EMFCUT   -1.   0.  0. 3. lastmat 0. PHOT-THR
1344       else if (strncmp(proc->GetName(),"COMP",4) == 0) {
1345           if (proc->Flag() == 1 || proc->Flag() == 2) { 
1346               fprintf(pFlukaVmcInp,"*\n*Energy threshold (GeV) for Compton scattering - resets to default=0.\n");
1347               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('COMP',1);\n");
1348               // - one = energy threshold (GeV) for Compton scattering - resets to default=0.
1349               // zero = not used
1350               // zero = not used
1351               // matMin = lower bound of the material indices in which the respective thresholds apply
1352               // matMax = upper bound of the material indices in which the respective thresholds apply
1353               // one = step length in assigning indices
1354               //"PHOT-THR"; 
1355               fprintf(pFlukaVmcInp,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",-one,zero,zero,matMin,matMax,one);
1356           }
1357           else if (proc->Flag() == 0) {
1358               fprintf(pFlukaVmcInp,"*\n*No Compton scattering - no FLUKA card generated\n");
1359               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('COMP',0)\n");
1360           }
1361           else  {
1362               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('COMP',?) call.\n");
1363               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1364           }
1365       } // end of else if (strncmp(proc->GetName(),"COMP",4) == 0)
1366       
1367       // decay
1368       // G3 default value: 1
1369       // G4 process: G4Decay
1370       // 
1371       // Particles: all which decay is applicable for
1372       // Physics:   General
1373       // flag = 0 no decays
1374       // flag = 1 decays, secondaries processed
1375       // flag = 2 decays, no secondaries stored
1376       //gMC ->SetProcess("DCAY",0); // not available
1377       else if ((strncmp(proc->GetName(),"DCAY",4) == 0) && proc->Flag() == 0) 
1378           cout << "SetProcess for flag =" << proc->GetName() << " value=" << proc->Flag() << " not avaliable!" << endl;
1379       else if ((strncmp(proc->GetName(),"DCAY",4) == 0) && proc->Flag() == 1) {
1380           // Nothing to do decays are switched on by default
1381       }
1382       
1383       
1384       // delta-ray
1385       // G3 default value: 2
1386       // !! G4 treats delta rays in different way
1387       // G4 processes: G4eIonisation/G4IeIonization,
1388       //               G4MuIonisation/G4IMuIonization,
1389       //               G4hIonisation/G4IhIonisation
1390       // Particles: charged
1391       // Physics:   EM
1392       // flag = 0 no energy loss
1393       // flag = 1 restricted energy loss fluctuations
1394       // flag = 2 complete energy loss fluctuations
1395       // flag = 3 same as 1
1396       // flag = 4 no energy loss fluctuations
1397       // gMC ->SetProcess("DRAY",0); // DELTARAY 1.E+6 0.  0. 3. lastmat 0.
1398       else if (strncmp(proc->GetName(),"DRAY",4) == 0) {
1399           if (proc->Flag() == 0 || proc->Flag() == 4) {
1400               fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n");
1401               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('DRAY',0) or SetProcess('DRAY',4)\n");
1402               fprintf(pFlukaVmcInp,"*No delta ray production by muons - threshold set artificially high\n");
1403               Double_t emin = 1.0e+6; // kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
1404               // zero = ignored
1405               // zero = ignored
1406               // matMin = lower bound of the material indices in which the respective thresholds apply
1407               // matMax = upper bound of the material indices in which the respective thresholds apply
1408               // one = step length in assigning indices
1409               fprintf(pFlukaVmcInp,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n",emin,zero,zero,matMin,matMax,one);
1410           }
1411           else if (proc->Flag() == 1 || proc->Flag() == 2 || proc->Flag() == 3) {
1412               fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n");
1413               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('DRAY',flag), flag=1,2,3\n");
1414               fprintf(pFlukaVmcInp,"*Delta ray production by muons switched on\n");
1415               fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('DCUTM',cut) or set to 1.0e+6.\n");
1416               fCut = 1.0e+6;
1417               nextc.Reset();
1418               while ((cut = (TFlukaConfigOption*)nextc())) {
1419                   if (strncmp(cut->GetName(), "DCUTM", 5) == 0 &&
1420                       cut->Medium() == proc->Medium()) fCut = cut->Cut();
1421               }
1422               // fCut = kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
1423               // zero = ignored
1424               // zero = ignored
1425               // matMin = lower bound of the material indices in which the respective thresholds apply
1426               // matMax =  upper bound of the material indices in which the respective thresholds apply
1427               // one = step length in assigning indices
1428               fprintf(pFlukaVmcInp,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n",fCut,zero,zero,matMin,matMax,one);
1429           }
1430           else  {
1431               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('DRAY',?) call.\n");
1432               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1433           }
1434       } // end of else if (strncmp(proc->GetName(),"DRAY",4) == 0)
1435       
1436       // hadronic process
1437       // G3 default value: 1
1438       // G4 processes: all defined by TG4PhysicsConstructorHadron
1439       //  
1440       // Particles: hadrons
1441       // Physics:   Hadron
1442       // flag = 0 no multiple scattering
1443       // flag = 1 hadronic interactions, secondaries processed
1444       // flag = 2 hadronic interactions, no secondaries stored
1445       // gMC ->SetProcess("HADR",1); // ??? hadronic process
1446       //Select pure GEANH (HADR 1) or GEANH/NUCRIN (HADR 3) ?????
1447       else if (strncmp(proc->GetName(),"HADR",4) == 0) {
1448           if (proc->Flag() == 1 || proc->Flag() == 2) {
1449               fprintf(pFlukaVmcInp,"*\n*Hadronic interaction is ON by default in FLUKA\n");
1450               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1451           }
1452           else if (proc->Flag() == 0) {
1453               fprintf(pFlukaVmcInp,"*\n*Hadronic interaction is set OFF\n");
1454               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('HADR',0);\n");
1455               fprintf(pFlukaVmcInp,"*Switching off hadronic interactions not foreseen in FLUKA\n");
1456               fprintf(pFlukaVmcInp,"THRESHOL  %10.1f%10.1f%10.1f%10.1e%10.1f\n",zero, zero, zero, 1.e10, zero);
1457           }
1458           else  {
1459               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('HADR',?) call.\n");
1460               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1461           }
1462       } // end of else if (strncmp(proc->GetName(),"HADR",4) == 0)
1463       
1464       
1465       // energy loss
1466       // G3 default value: 2
1467       // G4 processes: G4eIonisation/G4IeIonization,
1468       //               G4MuIonisation/G4IMuIonization,
1469       //               G4hIonisation/G4IhIonisation
1470       // 
1471       // Particles: charged
1472       // Physics:   EM
1473       // flag=0 no energy loss
1474       // flag=1 restricted energy loss fluctuations
1475       // flag=2 complete energy loss fluctuations
1476       // flag=3 same as 1
1477       // flag=4 no energy loss fluctuations
1478       // If the value ILOSS is changed, then (in G3) cross-sections and energy
1479       // loss tables must be recomputed via the command 'PHYSI'
1480       // gMC ->SetProcess("LOSS",2); // ??? IONFLUCT ? energy loss
1481       else if (strncmp(proc->GetName(),"LOSS",4) == 0) {
1482           if (proc->Flag() == 2) { // complete energy loss fluctuations
1483               fprintf(pFlukaVmcInp,"*\n*Complete energy loss fluctuations do not exist in FLUKA\n");
1484               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('LOSS',2);\n");
1485               fprintf(pFlukaVmcInp,"*flag=2=complete energy loss fluctuations\n");
1486               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1487           }
1488           else if (proc->Flag() == 1 || proc->Flag() == 3) { // restricted energy loss fluctuations
1489               fprintf(pFlukaVmcInp,"*\n*Restricted energy loss fluctuations\n");
1490               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('LOSS',1) or SetProcess('LOSS',3)\n");
1491               // one = restricted energy loss fluctuations (for hadrons and muons) switched on
1492               // one = restricted energy loss fluctuations (for e+ and e-) switched on
1493               // one = minimal accuracy
1494               // matMin = lower bound of the material indices in which the respective thresholds apply
1495               // upper bound of the material indices in which the respective thresholds apply
1496               fprintf(pFlukaVmcInp,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,one,one,matMin,matMax);
1497           }
1498           else if (proc->Flag() == 4) { // no energy loss fluctuations
1499               fprintf(pFlukaVmcInp,"*\n*No energy loss fluctuations\n");
1500               fprintf(pFlukaVmcInp,"*\n*Generated from call: SetProcess('LOSS',4)\n");
1501               // - one = restricted energy loss fluctuations (for hadrons and muons) switched off
1502               // - one = restricted energy loss fluctuations (for e+ and e-) switched off
1503               // one = minimal accuracy
1504               // matMin = lower bound of the material indices in which the respective thresholds apply
1505               // matMax = upper bound of the material indices in which the respective thresholds apply
1506               fprintf(pFlukaVmcInp,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,-one,one,matMin,matMax);
1507           }
1508           else  {
1509               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('LOSS',?) call.\n");
1510               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1511           }
1512       } // end of else if (strncmp(proc->GetName(),"LOSS",4) == 0)
1513       
1514       
1515       // multiple scattering
1516       // G3 default value: 1
1517       // G4 process: G4MultipleScattering/G4IMultipleScattering
1518       // 
1519       // Particles: charged
1520       // Physics:   EM
1521       // flag = 0 no multiple scattering
1522       // flag = 1 Moliere or Coulomb scattering
1523       // flag = 2 Moliere or Coulomb scattering
1524       // flag = 3 Gaussian scattering
1525       // gMC ->SetProcess("MULS",1); // MULSOPT multiple scattering
1526       else if (strncmp(proc->GetName(),"MULS",4) == 0) {
1527           if (proc->Flag() == 1 || proc->Flag() == 2 || proc->Flag() == 3) {
1528               fprintf(pFlukaVmcInp,"*\n*Multiple scattering is ON by default for e+e- and for hadrons/muons\n");
1529               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1530           }
1531           else if (proc->Flag() == 0) {
1532               fprintf(pFlukaVmcInp,"*\n*Multiple scattering is set OFF\n");
1533               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('MULS',0);\n");
1534               // zero = ignored
1535               // three = multiple scattering for hadrons and muons is completely suppressed
1536               // three = multiple scattering for e+ and e- is completely suppressed
1537               // matMin = lower bound of the material indices in which the respective thresholds apply
1538               // matMax = upper bound of the material indices in which the respective thresholds apply
1539               fprintf(pFlukaVmcInp,"MULSOPT   %10.1f%10.1f%10.1f%10.1f%10.1f\n",zero,three,three,matMin,matMax);
1540           }
1541           else  {
1542               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('MULS',?) call.\n");
1543               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1544           }
1545       } // end of else if (strncmp(proc->GetName(),"MULS",4) == 0)
1546       
1547
1548       // muon nuclear interaction
1549       // G3 default value: 0
1550       // G4 processes: G4MuNuclearInteraction,
1551       // G4MuonMinusCaptureAtRest
1552       // 
1553       // Particles: mu
1554       // Physics:   Not set
1555       // flag = 0 no muon-nuclear interaction
1556       // flag = 1 nuclear interaction, secondaries processed
1557       // flag = 2 nuclear interaction, secondaries not processed
1558       // gMC ->SetProcess("MUNU",1); // MUPHOTON  1.   0.  0. 3. lastmat
1559       else if (strncmp(proc->GetName(),"MUNU",4) == 0) {
1560           if (proc->Flag() == 1) {
1561               fprintf(pFlukaVmcInp,"*\n*Muon nuclear interactions with production of secondary hadrons\n");
1562               fprintf(pFlukaVmcInp,"*\n*Generated from call: SetProcess('MUNU',1);\n");
1563               // one = full simulation of muon nuclear interactions and production of secondary hadrons
1564               // zero = ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
1565               // zero = fraction of rho-like interactions ( must be < 1) - Default = 0.75.
1566               // matMin = lower bound of the material indices in which the respective thresholds apply
1567               // matMax = upper bound of the material indices in which the respective thresholds apply
1568               fprintf(pFlukaVmcInp,"MUPHOTON  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax);
1569           }
1570           else if (proc->Flag() == 2) {
1571               fprintf(pFlukaVmcInp,"*\n*Muon nuclear interactions without production of secondary hadrons\n");
1572               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('MUNU',2);\n");
1573               // two = full simulation of muon nuclear interactions and production of secondary hadrons
1574               // zero = ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
1575               // zero = fraction of rho-like interactions ( must be < 1) - Default = 0.75.
1576               // matMin = lower bound of the material indices in which the respective thresholds apply
1577               // matMax = upper bound of the material indices in which the respective thresholds apply
1578               fprintf(pFlukaVmcInp,"MUPHOTON  %10.1f%10.1f%10.1f%10.1f%10.1f\n",two,zero,zero,matMin,matMax);
1579           }
1580           else if (proc->Flag() == 0) {
1581               fprintf(pFlukaVmcInp,"*\n*No muon nuclear interaction - no FLUKA card generated\n");
1582               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('MUNU',0)\n");
1583           }
1584           else  {
1585               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('MUNU',?) call.\n");
1586               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1587           }
1588       } // end of else if (strncmp(proc->GetName(),"MUNU",4) == 0)
1589       
1590       
1591       // photofission
1592       // G3 default value: 0
1593       // G4 process: ??
1594       //
1595       // Particles: gamma
1596       // Physics:   ??
1597       // gMC ->SetProcess("PFIS",0); // PHOTONUC -1.   0.  0. 3. lastmat 0.
1598       // flag = 0 no photon fission
1599       // flag = 1 photon fission, secondaries processed
1600       // flag = 2 photon fission, no secondaries stored
1601       else if (strncmp(proc->GetName(),"PFIS",4) == 0) {
1602           if (proc->Flag() == 0) {
1603               fprintf(pFlukaVmcInp,"*\n*No photonuclear interactions\n");
1604               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PFIS',0);\n");
1605               // - one = no photonuclear interactions
1606               // zero = not used
1607               // zero = not used
1608               // matMin = lower bound of the material indices in which the respective thresholds apply
1609               // matMax = upper bound of the material indices in which the respective thresholds apply
1610               fprintf(pFlukaVmcInp,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,zero,zero,matMin,matMax);
1611           }
1612           else if (proc->Flag() == 1) {
1613               fprintf(pFlukaVmcInp,"*\n*Photon nuclear interactions are activated at all energies\n");
1614               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PFIS',1);\n");
1615               // one = photonuclear interactions are activated at all energies
1616               // zero = not used
1617               // zero = not used
1618               // matMin = lower bound of the material indices in which the respective thresholds apply
1619               // matMax = upper bound of the material indices in which the respective thresholds apply
1620               fprintf(pFlukaVmcInp,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax);
1621           }
1622           else if (proc->Flag() == 0) {
1623               fprintf(pFlukaVmcInp,"*\n*No photofission - no FLUKA card generated\n");
1624               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PFIS',0)\n");
1625           }
1626           else {
1627               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('PFIS',?) call.\n");
1628               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1629           }
1630       }
1631
1632  
1633       // photo electric effect
1634       // G3 default value: 1
1635       // G4 processes: G4PhotoElectricEffect
1636       //               G4LowEnergyPhotoElectric
1637       // Particles: gamma
1638       // Physics:   EM
1639       // flag = 0 no photo electric effect
1640       // flag = 1 photo electric effect, electron processed
1641       // flag = 2 photo electric effect, no electron stored
1642       // gMC ->SetProcess("PHOT",1); // EMFCUT    0.  -1.  0. 3. lastmat 0. PHOT-THR
1643       else if (strncmp(proc->GetName(),"PHOT",4) == 0) {
1644           if (proc->Flag() == 1 || proc->Flag() == 2) {
1645               fprintf(pFlukaVmcInp,"*\n*Photo electric effect is activated\n");
1646               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PHOT',1);\n");
1647               // zero = ignored
1648               // - one = resets to default=0.
1649               // zero = ignored
1650               // matMin = lower bound of the material indices in which the respective thresholds apply
1651               // matMax = upper bound of the material indices in which the respective thresholds apply
1652               // one = step length in assigning indices
1653               //"PHOT-THR"; 
1654               fprintf(pFlukaVmcInp,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",zero,-one,zero,matMin,matMax,one);
1655           }
1656           else if (proc->Flag() == 0) {
1657               fprintf(pFlukaVmcInp,"*\n*No photo electric effect - no FLUKA card generated\n");
1658               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PHOT',0)\n");
1659           }
1660           else {
1661               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('PHOT',?) call.\n");
1662               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1663           }
1664       } // else if (strncmp(proc->GetName(),"PHOT",4) == 0)
1665       
1666       
1667       // Rayleigh scattering
1668       // G3 default value: 0
1669       // G4 process: G4OpRayleigh
1670       // 
1671       // Particles: optical photon
1672       // Physics:   Optical
1673       // flag = 0 Rayleigh scattering off
1674       // flag = 1 Rayleigh scattering on
1675       //xx gMC ->SetProcess("RAYL",1);
1676       else if (strncmp(proc->GetName(),"RAYL",4) == 0) {
1677           if (proc->Flag() == 1) {
1678               fprintf(pFlukaVmcInp,"*\n*Rayleigh scattering is ON by default in FLUKA\n");
1679               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1680           }
1681           else if (proc->Flag() == 0) {
1682               fprintf(pFlukaVmcInp,"*\n*Rayleigh scattering is set OFF\n");
1683               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('RAYL',0);\n");
1684               // - one = no Rayleigh scattering and no binding corrections for Compton
1685               // matMin = lower bound of the material indices in which the respective thresholds apply
1686               // matMax = upper bound of the material indices in which the respective thresholds apply
1687               fprintf(pFlukaVmcInp,"EMFRAY    %10.1f%10.1f%10.1f%10.1f\n",-one,three,matMin,matMax);
1688           }
1689           else  {
1690               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('RAYL',?) call.\n");
1691               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1692           }
1693       } // end of else if (strncmp(proc->GetName(),"RAYL",4) == 0)
1694       
1695       
1696       // synchrotron radiation in magnetic field
1697       // G3 default value: 0
1698       // G4 process: G4SynchrotronRadiation
1699       // 
1700       // Particles: ??
1701       // Physics:   Not set
1702       // flag = 0 no synchrotron radiation
1703       // flag = 1 synchrotron radiation
1704       //xx gMC ->SetProcess("SYNC",1); // synchrotron radiation generation
1705       else if (strncmp(proc->GetName(),"SYNC",4) == 0) {
1706           fprintf(pFlukaVmcInp,"*\n*Synchrotron radiation generation is NOT implemented in FLUKA\n");
1707           fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1708       }
1709       
1710       
1711       // Automatic calculation of tracking medium parameters
1712       // flag = 0 no automatic calculation
1713       // flag = 1 automatic calculation
1714       //xx gMC ->SetProcess("AUTO",1); // ??? automatic computation of the tracking medium parameters
1715       else if (strncmp(proc->GetName(),"AUTO",4) == 0) {
1716           fprintf(pFlukaVmcInp,"*\n*Automatic calculation of tracking medium parameters is always ON in FLUKA\n");
1717           fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1718       }
1719       
1720       
1721       // To control energy loss fluctuation model
1722       // flag = 0 Urban model
1723       // flag = 1 PAI model
1724       // flag = 2 PAI+ASHO model (not active at the moment)
1725       //xx gMC ->SetProcess("STRA",1); // ??? energy fluctuation model
1726       else if (strncmp(proc->GetName(),"STRA",4) == 0) {
1727           if (proc->Flag() == 0 || proc->Flag() == 2 || proc->Flag() == 3) {
1728               fprintf(pFlukaVmcInp,"*\n*Ionization energy losses calculation is activated\n");
1729               fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('STRA',n);, n=0,1,2\n");
1730               // one = restricted energy loss fluctuations (for hadrons and muons) switched on
1731               // one = restricted energy loss fluctuations (for e+ and e-) switched on
1732               // one = minimal accuracy
1733               // matMin = lower bound of the material indices in which the respective thresholds apply
1734               // matMax = upper bound of the material indices in which the respective thresholds apply
1735               fprintf(pFlukaVmcInp,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,one,one,matMin,matMax);
1736           }
1737           else {
1738               fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('STRA',?) call.\n");
1739               fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
1740           }
1741       } // else if (strncmp(proc->GetName(),"STRA",4) == 0)
1742       
1743
1744
1745
1746       else { // processes not yet treated
1747           
1748           // light photon absorption (Cerenkov photons)
1749           // it is turned on when Cerenkov process is turned on
1750           // G3 default value: 0
1751           // G4 process: G4OpAbsorption, G4OpBoundaryProcess
1752           // 
1753           // Particles: optical photon
1754           // Physics:   Optical
1755           // flag = 0 no absorption of Cerenkov photons
1756           // flag = 1 absorption of Cerenkov photons
1757           // gMC ->SetProcess("LABS",2); // ??? Cerenkov light absorption
1758           
1759
1760
1761           cout << "SetProcess for flag=" << proc->GetName() << " value=" << proc->Flag() << " not yet implemented!" << endl;
1762       }
1763   } //end of loop number of SetProcess calls
1764
1765  
1766 // Loop over number of SetCut calls  
1767             
1768   nextc.Reset();
1769   while ((cut = (TFlukaConfigOption*)nextc())) {
1770       Float_t matMin = three;
1771       Float_t matMax = fLastMaterial;
1772       Bool_t global  = kTRUE;
1773       if (cut->Medium() != -1) {
1774         matMin = Float_t(cut->Medium());
1775         matMax = matMin;
1776         global = kFALSE;
1777       }
1778
1779       // cuts handled in SetProcess calls
1780       if (strncmp(cut->GetName(),"BCUTM",5) == 0) continue;
1781       else if (strncmp(cut->GetName(),"BCUTE",5) == 0) continue;
1782       else if (strncmp(cut->GetName(),"DCUTM",5) == 0) continue;
1783       else if (strncmp(cut->GetName(),"PPCUTM",6) == 0) continue;
1784       
1785       // delta-rays by electrons
1786       // G4 particles: "e-"
1787       // G3 default value: 10**4 GeV
1788       // gMC ->SetCut("DCUTE",cut);  // cut for deltarays by electrons 
1789       else if (strncmp(cut->GetName(),"DCUTE",5) == 0) {
1790         fprintf(pFlukaVmcInp,"*\n*Cut for delta rays by electrons\n");
1791         fprintf(pFlukaVmcInp,"*Generated from call: SetCut('DCUTE',cut);\n");
1792         // -cut->Cut();
1793         // zero = ignored
1794         // zero = ignored
1795         // matMin = lower bound of the material indices in which the respective thresholds apply
1796         // matMax = upper bound of the material indices in which the respective thresholds apply
1797         // loop over materials for EMFCUT FLUKA cards
1798         for (j=0; j < matMax-matMin+1; j++) {
1799           Int_t nreg, imat, *reglist;
1800           Float_t ireg;
1801           imat = (Int_t) matMin + j;
1802           reglist = fGeom->GetMaterialList(imat, nreg);
1803           // loop over regions of a given material
1804           for (k=0; k<nreg; k++) {
1805             ireg = reglist[k];
1806             fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f\n",-cut->Cut(),zero,zero,ireg,ireg);
1807           }
1808         }
1809         fprintf(pFlukaVmcInp,"DELTARAY  %10.4g%10.3f%10.3f%10.1f%10.1f%10.1f\n",cut->Cut(), 100., 1.03, matMin, matMax, 1.0);
1810       } // end of if for delta-rays by electrons
1811     
1812
1813       // gammas
1814       // G4 particles: "gamma"
1815       // G3 default value: 0.001 GeV
1816       // gMC ->SetCut("CUTGAM",cut); // cut for gammas
1817       
1818       else if (strncmp(cut->GetName(),"CUTGAM",6) == 0 && global) {
1819         fprintf(pFlukaVmcInp,"*\n*Cut for gamma\n");
1820         fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTGAM',cut);\n");
1821         // -cut->Cut();
1822         // 7.0 = lower bound of the particle id-numbers to which the cut-off
1823         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f\n",-cut->Cut(),7.0);
1824       }
1825       else if (strncmp(cut->GetName(),"CUTGAM",6) == 0 && !global) {
1826         fprintf(pFlukaVmcInp,"*\n*Cut specific to  material for gamma\n");
1827         fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTGAM',cut);\n");
1828         // cut->Cut();
1829         // loop over materials for EMFCUT FLUKA cards
1830         for (j=0; j < matMax-matMin+1; j++) {
1831           Int_t nreg, imat, *reglist;
1832           Float_t ireg;
1833           imat = (Int_t) matMin + j;
1834           reglist = fGeom->GetMaterialList(imat, nreg);
1835           // loop over regions of a given material
1836           for (Int_t k=0; k<nreg; k++) {
1837             ireg = reglist[k];
1838             fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", zero, cut->Cut(), zero, ireg, ireg, one);
1839           }
1840         }
1841       } // end of else if for gamma
1842
1843
1844       // electrons
1845       // G4 particles: "e-"
1846       // ?? positrons
1847       // G3 default value: 0.001 GeV
1848       //gMC ->SetCut("CUTELE",cut); // cut for e+,e-
1849       else if (strncmp(cut->GetName(),"CUTELE",6) == 0 && global) {
1850         fprintf(pFlukaVmcInp,"*\n*Cut for electrons\n");
1851         fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTELE',cut);\n");
1852         // -cut->Cut();
1853         // three = lower bound of the particle id-numbers to which the cut-off
1854         // 4.0 = upper bound of the particle id-numbers to which the cut-off
1855         // one = step length in assigning numbers
1856         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f%10.1f\n",-cut->Cut(),three,4.0,one);
1857       }
1858       else if (strncmp(cut->GetName(),"CUTELE",6) == 0 && !global) {
1859         fprintf(pFlukaVmcInp,"*\n*Cut specific to material for electrons\n");
1860         fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTELE',cut);\n");
1861         // -cut->Cut();
1862         // loop over materials for EMFCUT FLUKA cards
1863         for (j=0; j < matMax-matMin+1; j++) {
1864           Int_t nreg, imat, *reglist;
1865           Float_t ireg;
1866           imat = (Int_t) matMin + j;
1867           reglist = fGeom->GetMaterialList(imat, nreg);
1868           // loop over regions of a given material
1869           for (k=0; k<nreg; k++) {
1870             ireg = reglist[k];
1871             fprintf(pFlukaVmcInp,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -cut->Cut(), zero, zero, ireg, ireg, one);
1872           }
1873         }
1874       } // end of else if for electrons
1875
1876     
1877       // neutral hadrons
1878       // G4 particles: of type "baryon", "meson", "nucleus" with zero charge
1879       // G3 default value: 0.01 GeV
1880       //gMC ->SetCut("CUTNEU",cut); // cut for neutral hadrons
1881       else if (strncmp(cut->GetName(),"CUTNEU",6) == 0 && global) {
1882         fprintf(pFlukaVmcInp,"*\n*Cut for neutral hadrons\n");
1883         fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTNEU',cut);\n");
1884           
1885         // 8.0 = Neutron
1886         // 9.0 = Antineutron
1887         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),8.0,9.0);
1888           
1889         // 12.0 = Kaon zero long
1890         // 12.0 = Kaon zero long
1891         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),12.0,12.0);
1892           
1893         // 17.0 = Lambda, 18.0 = Antilambda
1894         // 19.0 = Kaon zero short
1895         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),17.0,19.0);
1896           
1897         // 22.0 = Sigma zero, Pion zero, Kaon zero
1898         // 25.0 = Antikaon zero
1899         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),22.0,25.0);
1900           
1901         // 32.0 = Antisigma zero
1902         // 32.0 = Antisigma zero
1903         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),32.0,32.0);
1904           
1905         // 34.0 = Xi zero
1906         // 35.0 = AntiXi zero
1907         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),34.0,35.0);
1908           
1909         // 47.0 = D zero
1910         // 48.0 = AntiD zero
1911         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),47.0,48.0);
1912           
1913         // 53.0 = Xi_c zero
1914         // 53.0 = Xi_c zero
1915         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),53.0,53.0);
1916           
1917         // 55.0 = Xi'_c zero
1918         // 56.0 = Omega_c zero
1919         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),55.0,56.0);
1920           
1921         // 59.0 = AntiXi_c zero
1922         // 59.0 = AntiXi_c zero
1923         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),59.0,59.0);
1924           
1925         // 61.0 = AntiXi'_c zero
1926         // 62.0 = AntiOmega_c zero
1927         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),61.0,62.0);
1928       }
1929       
1930       // charged hadrons
1931       // G4 particles: of type "baryon", "meson", "nucleus" with non-zero charge
1932       // G3 default value: 0.01 GeV
1933       //gMC ->SetCut("CUTHAD",cut); // cut for charged hadrons
1934       else if (strncmp(cut->GetName(),"CUTHAD",6) == 0 && global) {
1935         fprintf(pFlukaVmcInp,"*\n*Cut for charged hadrons\n");
1936         fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTHAD',cut);\n");
1937           
1938         // 1.0 = Proton
1939         // 2.0 = Antiproton
1940         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),1.0,2.0);
1941           
1942         // 13.0 = Positive Pion, Negative Pion, Positive Kaon
1943         // 16.0 = Negative Kaon
1944         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),13.0,16.0);
1945           
1946         // 20.0 = Negative Sigma
1947         // 21.0 = Positive Sigma
1948         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),20.0,21.0);
1949           
1950         // 31.0 = Antisigma minus
1951         // 33.0 = Antisigma plus
1952         // 2.0 = step length
1953         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f%10.1f\n",-cut->Cut(),31.0,33.0,2.0);
1954           
1955         // 36.0 = Negative Xi, Positive Xi, Omega minus
1956         // 39.0 = Antiomega
1957         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),36.0,39.0);
1958           
1959         // 45.0 = D plus
1960         // 46.0 = D minus
1961         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),45.0,46.0);
1962           
1963         // 49.0 = D_s plus, D_s minus, Lambda_c plus
1964         // 52.0 = Xi_c plus
1965         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),49.0,52.0);
1966           
1967         // 54.0 = Xi'_c plus
1968         // 60.0 = AntiXi'_c minus
1969         // 6.0 = step length
1970         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f%10.1f\n",-cut->Cut(),54.0,60.0,6.0);
1971           
1972         // 57.0 = Antilambda_c minus
1973         // 58.0 = AntiXi_c minus
1974         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),57.0,58.0);
1975       }
1976
1977       // muons
1978       // G4 particles: "mu+", "mu-"
1979       // G3 default value: 0.01 GeV
1980       //gMC ->SetCut("CUTMUO",cut); // cut for mu+, mu-
1981       else if (strncmp(cut->GetName(),"CUTMUO",6)== 0 && global) {
1982         fprintf(pFlukaVmcInp,"*\n*Cut for muons\n");
1983         fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTMUO',cut);\n");
1984         // 10.0 = Muon+
1985         // 11.0 = Muon-
1986         fprintf(pFlukaVmcInp,"PART-THR  %10.4g%10.1f%10.1f\n",-cut->Cut(),10.0,11.0);
1987       }
1988       
1989       //
1990       // time of flight cut in seconds
1991       // G4 particles: all
1992       // G3 default value: 0.01 GeV
1993       //gMC ->SetCut("TOFMAX",tofmax); // time of flight cuts in seconds
1994       else if (strncmp(cut->GetName(),"TOFMAX",6) == 0) {
1995         fprintf(pFlukaVmcInp,"*\n*Time of flight cuts in seconds\n");
1996         fprintf(pFlukaVmcInp,"*Generated from call: SetCut('TOFMAX',tofmax);\n");
1997         // zero = ignored
1998         // zero = ignored
1999         // -6.0 = lower bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
2000         // 64.0 = upper bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
2001         fprintf(pFlukaVmcInp,"TIME-CUT  %10.4g%10.1f%10.1f%10.1f%10.1f\n",cut->Cut()*1.e9,zero,zero,-6.0,64.0);
2002       }
2003       
2004       else if (global){
2005         cout << "SetCut for flag=" << cut->GetName() << " value=" << cut->Cut() << " not yet implemented!" << endl;
2006       }
2007       else {
2008         cout << "SetCut for flag=" << cut->GetName() << " value=" << cut->Cut() << " (material specific) not yet implemented!" << endl;
2009       }
2010       
2011   } //end of loop over SetCut calls
2012   
2013 // Add START and STOP card
2014   fprintf(pFlukaVmcInp,"START     %10.1f\n",fEventsPerRun);
2015   fprintf(pFlukaVmcInp,"STOP      \n");
2016    
2017   
2018 // Close files
2019   
2020    fclose(pFlukaVmcCoreInp);
2021    fclose(pFlukaVmcFlukaMat);
2022    fclose(pFlukaVmcInp);
2023    
2024 } // end of InitPhysics
2025
2026
2027 //______________________________________________________________________________ 
2028 void TFluka::SetMaxStep(Double_t step)
2029 {
2030 // Set the maximum step size
2031     if (step > 1.e4) return;
2032     
2033     Int_t mreg, latt;
2034     fGeom->GetCurrentRegion(mreg, latt);
2035     STEPSZ.stepmx[mreg - 1] = step;
2036 }
2037
2038 //______________________________________________________________________________ 
2039 void TFluka::SetMaxNStep(Int_t)
2040 {
2041 // SetMaxNStep is dummy procedure in TFluka !
2042   if (fVerbosityLevel >=3)
2043   cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
2044 }
2045
2046 //______________________________________________________________________________ 
2047 void TFluka::SetUserDecay(Int_t)
2048 {
2049 // SetUserDecay is dummy procedure in TFluka !
2050   if (fVerbosityLevel >=3)
2051   cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
2052 }
2053
2054 //
2055 // dynamic properties
2056 //
2057 //______________________________________________________________________________ 
2058 void TFluka::TrackPosition(TLorentzVector& position) const
2059 {
2060 // Return the current position in the master reference frame of the
2061 // track being transported
2062 // TRACKR.atrack = age of the particle
2063 // TRACKR.xtrack = x-position of the last point
2064 // TRACKR.ytrack = y-position of the last point
2065 // TRACKR.ztrack = z-position of the last point
2066   Int_t caller = GetCaller();
2067   if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
2068     position.SetX(GetXsco());
2069     position.SetY(GetYsco());
2070     position.SetZ(GetZsco());
2071     position.SetT(TRACKR.atrack);
2072   }
2073   else if (caller == 4) { // mgdraw
2074     position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
2075     position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
2076     position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
2077     position.SetT(TRACKR.atrack);
2078   }
2079   else if (caller == 5) { // sodraw
2080     position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
2081     position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
2082     position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
2083     position.SetT(0);
2084   }
2085   else
2086     Warning("TrackPosition","position not available");
2087 }
2088
2089 //______________________________________________________________________________ 
2090 void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
2091 {
2092 // Return the current position in the master reference frame of the
2093 // track being transported
2094 // TRACKR.atrack = age of the particle
2095 // TRACKR.xtrack = x-position of the last point
2096 // TRACKR.ytrack = y-position of the last point
2097 // TRACKR.ztrack = z-position of the last point
2098   Int_t caller = GetCaller();
2099   if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
2100     x = GetXsco();
2101     y = GetYsco();
2102     z = GetZsco();
2103   }
2104   else if (caller == 4 || caller == 5) { // mgdraw, sodraw
2105     x = TRACKR.xtrack[TRACKR.ntrack];
2106     y = TRACKR.ytrack[TRACKR.ntrack];
2107     z = TRACKR.ztrack[TRACKR.ntrack];
2108   }
2109   else
2110     Warning("TrackPosition","position not available");
2111 }
2112
2113 //______________________________________________________________________________ 
2114 void TFluka::TrackMomentum(TLorentzVector& momentum) const
2115 {
2116 // Return the direction and the momentum (GeV/c) of the track
2117 // currently being transported
2118 // TRACKR.ptrack = momentum of the particle (not always defined, if
2119 //               < 0 must be obtained from etrack) 
2120 // TRACKR.cx,y,ztrck = direction cosines of the current particle
2121 // TRACKR.etrack = total energy of the particle
2122 // TRACKR.jtrack = identity number of the particle
2123 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
2124   Int_t caller = GetCaller();
2125   if (caller != 2) { // not eedraw 
2126     if (TRACKR.ptrack >= 0) {
2127       momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
2128       momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
2129       momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
2130       momentum.SetE(TRACKR.etrack);
2131       return;
2132     }
2133     else {
2134       Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
2135       momentum.SetPx(p*TRACKR.cxtrck);
2136       momentum.SetPy(p*TRACKR.cytrck);
2137       momentum.SetPz(p*TRACKR.cztrck);
2138       momentum.SetE(TRACKR.etrack);
2139       return;
2140     }
2141   }
2142   else
2143     Warning("TrackMomentum","momentum not available");
2144 }
2145
2146 //______________________________________________________________________________ 
2147 void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e) const
2148 {
2149 // Return the direction and the momentum (GeV/c) of the track
2150 // currently being transported
2151 // TRACKR.ptrack = momentum of the particle (not always defined, if
2152 //               < 0 must be obtained from etrack) 
2153 // TRACKR.cx,y,ztrck = direction cosines of the current particle
2154 // TRACKR.etrack = total energy of the particle
2155 // TRACKR.jtrack = identity number of the particle
2156 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
2157   Int_t caller = GetCaller();
2158   if (caller != 2) { // not eedraw 
2159     if (TRACKR.ptrack >= 0) {
2160       px = TRACKR.ptrack*TRACKR.cxtrck;
2161       py = TRACKR.ptrack*TRACKR.cytrck;
2162       pz = TRACKR.ptrack*TRACKR.cztrck;
2163       e = TRACKR.etrack;
2164       return;
2165     }
2166     else {
2167       Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
2168       px = p*TRACKR.cxtrck;
2169       py = p*TRACKR.cytrck;
2170       pz = p*TRACKR.cztrck;
2171       e = TRACKR.etrack;
2172       return;
2173     }
2174   }
2175   else
2176     Warning("TrackMomentum","momentum not available");
2177 }
2178
2179 //______________________________________________________________________________ 
2180 Double_t TFluka::TrackStep() const
2181 {
2182 // Return the length in centimeters of the current step
2183 // TRACKR.ctrack = total curved path
2184   Int_t caller = GetCaller();
2185   if (caller == 11 || caller==12 || caller == 3 || caller == 6) //bxdraw,endraw,usdraw
2186     return 0.0;
2187   else if (caller == 4) //mgdraw
2188     return TRACKR.ctrack;
2189   else
2190     return -1.0;
2191 }
2192
2193 //______________________________________________________________________________ 
2194 Double_t TFluka::TrackLength() const
2195 {
2196 // TRACKR.cmtrck = cumulative curved path since particle birth
2197   Int_t caller = GetCaller();
2198   if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
2199     return TRACKR.cmtrck;
2200   else 
2201     return -1.0;
2202 }
2203
2204 //______________________________________________________________________________ 
2205 Double_t TFluka::TrackTime() const
2206 {
2207 // Return the current time of flight of the track being transported
2208 // TRACKR.atrack = age of the particle
2209   Int_t caller = GetCaller();
2210   if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
2211     return TRACKR.atrack;
2212   else 
2213     return -1;
2214 }
2215
2216 //______________________________________________________________________________ 
2217 Double_t TFluka::Edep() const
2218 {
2219 // Energy deposition
2220 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
2221 // -->local energy deposition (the value and the point are not recorded in TRACKR)
2222 //    but in the variable "rull" of the procedure "endraw.cxx"
2223 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
2224 // -->no energy loss along the track
2225 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
2226 // -->energy loss distributed along the track
2227 // TRACKR.dtrack = energy deposition of the jth deposition event
2228
2229   // If coming from bxdraw we have 2 steps of 0 length and 0 edep
2230   Int_t caller = GetCaller();
2231   if (caller == 11 || caller==12) return 0.0;
2232   Double_t sum = 0;
2233   for ( Int_t j=0;j<TRACKR.mtrack;j++) {
2234     sum +=TRACKR.dtrack[j];  
2235   }
2236   if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
2237     return fRull + sum;
2238   else {
2239     return sum;
2240   }
2241 }
2242
2243 //______________________________________________________________________________ 
2244 Int_t TFluka::TrackPid() const
2245 {
2246 // Return the id of the particle transported
2247 // TRACKR.jtrack = identity number of the particle
2248   Int_t caller = GetCaller();
2249   if (caller != 2) { // not eedraw 
2250       return PDGFromId(TRACKR.jtrack);
2251   }
2252   else
2253     return -1000;
2254 }
2255
2256 //______________________________________________________________________________ 
2257 Double_t TFluka::TrackCharge() const
2258 {
2259 // Return charge of the track currently transported
2260 // PAPROP.ichrge = electric charge of the particle
2261 // TRACKR.jtrack = identity number of the particle
2262   Int_t caller = GetCaller();
2263   if (caller != 2)  // not eedraw 
2264     return PAPROP.ichrge[TRACKR.jtrack+6];
2265   else
2266     return -1000.0;
2267 }
2268
2269 //______________________________________________________________________________ 
2270 Double_t TFluka::TrackMass() const
2271 {
2272 // PAPROP.am = particle mass in GeV
2273 // TRACKR.jtrack = identity number of the particle
2274   Int_t caller = GetCaller();
2275   if (caller != 2)  // not eedraw 
2276     return PAPROP.am[TRACKR.jtrack+6];
2277   else
2278     return -1000.0;
2279 }
2280
2281 //______________________________________________________________________________ 
2282 Double_t TFluka::Etot() const
2283 {
2284 // TRACKR.etrack = total energy of the particle
2285   Int_t caller = GetCaller();
2286   if (caller != 2)  // not eedraw
2287     return TRACKR.etrack;
2288   else
2289     return -1000.0;
2290 }
2291
2292 //
2293 // track status
2294 //
2295 //______________________________________________________________________________ 
2296 Bool_t   TFluka::IsNewTrack() const
2297 {
2298 // Return true for the first call of Stepping()
2299    return fTrackIsNew;
2300 }
2301
2302 void     TFluka::SetTrackIsNew(Bool_t flag)
2303 {
2304 // Return true for the first call of Stepping()
2305    fTrackIsNew = flag;
2306
2307 }
2308
2309
2310 //______________________________________________________________________________ 
2311 Bool_t   TFluka::IsTrackInside() const
2312 {
2313 // True if the track is not at the boundary of the current volume
2314 // In Fluka a step is always inside one kind of material
2315 // If the step would go behind the region of one material,
2316 // it will be shortened to reach only the boundary.
2317 // Therefore IsTrackInside() is always true.
2318   Int_t caller = GetCaller();
2319   if (caller == 11 || caller==12)  // bxdraw
2320     return 0;
2321   else
2322     return 1;
2323 }
2324
2325 //______________________________________________________________________________ 
2326 Bool_t   TFluka::IsTrackEntering() const
2327 {
2328 // True if this is the first step of the track in the current volume
2329
2330   Int_t caller = GetCaller();
2331   if (caller == 11)  // bxdraw entering
2332     return 1;
2333   else return 0;
2334 }
2335
2336 //______________________________________________________________________________ 
2337 Bool_t   TFluka::IsTrackExiting() const
2338 {
2339 // True if track is exiting volume
2340 //
2341   Int_t caller = GetCaller();
2342   if (caller == 12)  // bxdraw exiting
2343     return 1;
2344   else return 0;
2345 }
2346
2347 //______________________________________________________________________________ 
2348 Bool_t   TFluka::IsTrackOut() const
2349 {
2350 // True if the track is out of the setup
2351 // means escape
2352 // Icode = 14: escape - call from Kaskad
2353 // Icode = 23: escape - call from Emfsco
2354 // Icode = 32: escape - call from Kasneu
2355 // Icode = 40: escape - call from Kashea
2356 // Icode = 51: escape - call from Kasoph
2357   if (fIcode == 14 ||
2358       fIcode == 23 ||
2359       fIcode == 32 ||
2360       fIcode == 40 ||
2361       fIcode == 51) return 1;
2362   else return 0;
2363 }
2364
2365 //______________________________________________________________________________ 
2366 Bool_t   TFluka::IsTrackDisappeared() const
2367 {
2368 // means all inelastic interactions and decays
2369 // fIcode from usdraw
2370   if (fIcode == 101 || // inelastic interaction
2371       fIcode == 102 || // particle decay
2372       fIcode == 103 || // delta ray generation by hadron
2373       fIcode == 104 || // direct pair production
2374       fIcode == 105 || // bremsstrahlung (muon)
2375       fIcode == 208 || // bremsstrahlung (electron)
2376       fIcode == 214 || // in-flight annihilation
2377       fIcode == 215 || // annihilation at rest
2378       fIcode == 217 || // pair production
2379       fIcode == 219 || // Compton scattering
2380       fIcode == 221 || // Photoelectric effect
2381       fIcode == 300 || // hadronic interaction
2382       fIcode == 400    // delta-ray
2383       ) return 1;
2384   else return 0;
2385 }
2386
2387 //______________________________________________________________________________ 
2388 Bool_t   TFluka::IsTrackStop() const
2389 {
2390 // True if the track energy has fallen below the threshold
2391 // means stopped by signal or below energy threshold
2392 // Icode = 12: stopping particle       - call from Kaskad
2393 // Icode = 15: time kill               - call from Kaskad
2394 // Icode = 21: below threshold, iarg=1 - call from Emfsco
2395 // Icode = 22: below threshold, iarg=2 - call from Emfsco
2396 // Icode = 24: time kill               - call from Emfsco
2397 // Icode = 31: below threshold         - call from Kasneu
2398 // Icode = 33: time kill               - call from Kasneu
2399 // Icode = 41: time kill               - call from Kashea
2400 // Icode = 52: time kill               - call from Kasoph
2401   if (fIcode == 12 ||
2402       fIcode == 15 ||
2403       fIcode == 21 ||
2404       fIcode == 22 ||
2405       fIcode == 24 ||
2406       fIcode == 31 ||
2407       fIcode == 33 ||
2408       fIcode == 41 ||
2409       fIcode == 52) return 1;
2410   else return 0;
2411 }
2412
2413 //______________________________________________________________________________ 
2414 Bool_t   TFluka::IsTrackAlive() const
2415 {
2416 // means not disappeared or not out
2417   if (IsTrackDisappeared() || IsTrackOut() ) return 0;
2418   else return 1;
2419 }
2420
2421 //
2422 // secondaries
2423 //
2424
2425 //______________________________________________________________________________ 
2426 Int_t TFluka::NSecondaries() const
2427
2428 {
2429 // Number of secondary particles generated in the current step
2430 // FINUC.np = number of secondaries except light and heavy ions
2431 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
2432   Int_t caller = GetCaller();
2433   if (caller == 6)  // valid only after usdraw
2434     return FINUC.np + FHEAVY.npheav;
2435   else
2436     return 0;
2437 } // end of NSecondaries
2438
2439 //______________________________________________________________________________ 
2440 void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
2441                 TLorentzVector& position, TLorentzVector& momentum)
2442 {
2443 // Copy particles from secondary stack to vmc stack
2444 //
2445
2446   Int_t caller = GetCaller();
2447   if (caller == 6) {  // valid only after usdraw
2448     if (isec >= 0 && isec < FINUC.np) {
2449       particleId = PDGFromId(FINUC.kpart[isec]);
2450       position.SetX(fXsco);
2451       position.SetY(fYsco);
2452       position.SetZ(fZsco);
2453       position.SetT(TRACKR.atrack);
2454       momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
2455       momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
2456       momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
2457       momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
2458     }
2459     else if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
2460       Int_t jsec = isec - FINUC.np;
2461       particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
2462       position.SetX(fXsco);
2463       position.SetY(fYsco);
2464       position.SetZ(fZsco);
2465       position.SetT(TRACKR.atrack);
2466       momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
2467       momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
2468       momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
2469       if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6) 
2470         momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
2471       else if (FHEAVY.tkheav[jsec] > 6)
2472         momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
2473     }
2474     else
2475       Warning("GetSecondary","isec out of range");
2476   }
2477   else
2478     Warning("GetSecondary","no secondaries available");
2479 } // end of GetSecondary
2480
2481 //______________________________________________________________________________ 
2482 TMCProcess TFluka::ProdProcess(Int_t) const
2483
2484 {
2485 // Name of the process that has produced the secondary particles
2486 // in the current step
2487
2488     Int_t mugamma = (TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11);
2489
2490     if      (fIcode == 102)                  return kPDecay;
2491     else if (fIcode == 104 || fIcode == 217) return kPPair;
2492     else if (fIcode == 219)                  return kPCompton;
2493     else if (fIcode == 221)                  return kPPhotoelectric;
2494     else if (fIcode == 105 || fIcode == 208) return kPBrem;
2495     else if (fIcode == 103 || fIcode == 400) return kPDeltaRay;
2496     else if (fIcode == 210 || fIcode == 212) return kPDeltaRay;
2497     else if (fIcode == 214 || fIcode == 215) return kPAnnihilation;
2498     else if (fIcode == 101)                  return kPHadronic;
2499     else if (fIcode == 101) {
2500         if (!mugamma)                        return kPHadronic;
2501         else if (TRACKR.jtrack == 7)         return kPPhotoFission;
2502         else return kPMuonNuclear;
2503     }
2504     else if (fIcode == 225)                  return kPRayleigh;
2505 // Fluka codes 100, 300 and 400 still to be investigasted
2506     else                                     return kPNoProcess;
2507 }
2508
2509
2510 Int_t TFluka::StepProcesses(TArrayI &proc) const
2511 {
2512   //
2513   // Return processes active in the current step
2514   //
2515     proc.Set(1);
2516     TMCProcess iproc;
2517     switch (fIcode) {
2518     case 15:
2519     case 24:
2520     case 33:
2521     case 41:
2522     case 52:
2523         iproc =  kPTOFlimit;
2524         break;
2525     case 12:
2526     case 14:
2527     case 21:
2528     case 22:
2529     case 23:
2530     case 31:
2531     case 32:
2532     case 40:
2533     case 51:
2534         iproc =  kPStop;
2535         break;
2536     case 50:
2537         iproc = kPLightAbsorption;
2538         break;
2539     case 20: 
2540         iproc = kPPhotoelectric;
2541         break;
2542     default:
2543         iproc = ProdProcess(0);
2544     }
2545     proc[0] = iproc;
2546     return 1;
2547 }
2548 //______________________________________________________________________________ 
2549 Int_t TFluka::VolId2Mate(Int_t id) const
2550 {
2551 //
2552 // Returns the material number for a given volume ID
2553 //
2554    return fMCGeo->VolId2Mate(id);
2555 }
2556
2557 //______________________________________________________________________________ 
2558 const char* TFluka::VolName(Int_t id) const
2559 {
2560 //
2561 // Returns the volume name for a given volume ID
2562 //
2563    return fMCGeo->VolName(id);
2564 }
2565
2566 //______________________________________________________________________________ 
2567 Int_t TFluka::VolId(const Text_t* volName) const
2568 {
2569 //
2570 // Converts from volume name to volume ID.
2571 // Time consuming. (Only used during set-up)
2572 // Could be replaced by hash-table
2573 //
2574    return fMCGeo->VolId(volName);
2575 }
2576
2577 //______________________________________________________________________________ 
2578 Int_t TFluka::CurrentVolID(Int_t& copyNo) const
2579 {
2580 //
2581 // Return the logical id and copy number corresponding to the current fluka region
2582 //
2583   if (gGeoManager->IsOutside()) return 0;
2584   TGeoNode *node = gGeoManager->GetCurrentNode();
2585   copyNo = node->GetNumber();
2586   Int_t id = node->GetVolume()->GetNumber();
2587   return id;
2588
2589
2590 //______________________________________________________________________________ 
2591 Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
2592 {
2593 //
2594 // Return the logical id and copy number of off'th mother 
2595 // corresponding to the current fluka region
2596 //
2597   if (off<0 || off>gGeoManager->GetLevel()) return 0;
2598   if (off==0) return CurrentVolID(copyNo);
2599   TGeoNode *node = gGeoManager->GetMother(off);
2600   if (!node) return 0;
2601   copyNo = node->GetNumber();
2602   return node->GetVolume()->GetNumber();
2603 }
2604
2605 //______________________________________________________________________________ 
2606 const char* TFluka::CurrentVolName() const
2607 {
2608 //
2609 // Return the current volume name
2610 //
2611   if (gGeoManager->IsOutside()) return 0;
2612   return gGeoManager->GetCurrentVolume()->GetName();
2613 }
2614
2615 //______________________________________________________________________________ 
2616 const char* TFluka::CurrentVolOffName(Int_t off) const
2617 {
2618 //
2619 // Return the volume name of the off'th mother of the current volume
2620 //
2621   if (off<0 || off>gGeoManager->GetLevel()) return 0;
2622   if (off==0) return CurrentVolName();
2623   TGeoNode *node = gGeoManager->GetMother(off);
2624   if (!node) return 0;
2625   return node->GetVolume()->GetName();
2626 }
2627
2628 //______________________________________________________________________________ 
2629 Int_t TFluka::CurrentMaterial(Float_t & /*a*/, Float_t & /*z*/, 
2630                       Float_t & /*dens*/, Float_t & /*radl*/, Float_t & /*absl*/) const
2631 {
2632 //
2633 //  Return the current medium number  ??? what about material properties
2634 //
2635   Int_t copy;
2636   Int_t id  =  TFluka::CurrentVolID(copy);
2637   Int_t med =  TFluka::VolId2Mate(id);
2638   return med;
2639 }
2640
2641 //______________________________________________________________________________ 
2642 void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
2643 {
2644 // Transforms a position from the world reference frame
2645 // to the current volume reference frame.
2646 //
2647 //  Geant3 desription:
2648 //  ==================
2649 //       Computes coordinates XD (in DRS) 
2650 //       from known coordinates XM in MRS 
2651 //       The local reference system can be initialized by
2652 //         - the tracking routines and GMTOD used in GUSTEP
2653 //         - a call to GMEDIA(XM,NUMED)
2654 //         - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER) 
2655 //             (inverse routine is GDTOM) 
2656 //
2657 //        If IFLAG=1  convert coordinates 
2658 //           IFLAG=2  convert direction cosinus
2659 //
2660 // ---
2661    Double_t xmL[3], xdL[3];
2662    Int_t i;
2663    for (i=0;i<3;i++) xmL[i]=xm[i];
2664    if (iflag == 1) gGeoManager->MasterToLocal(xmL,xdL);
2665    else            gGeoManager->MasterToLocalVect(xmL,xdL);
2666    for (i=0;i<3;i++) xd[i] = xdL[i];
2667 }
2668   
2669 //______________________________________________________________________________ 
2670 void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
2671 {
2672    if (iflag == 1) gGeoManager->MasterToLocal(xm,xd);
2673    else            gGeoManager->MasterToLocalVect(xm,xd);
2674 }
2675
2676 //______________________________________________________________________________ 
2677 void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
2678 {
2679 // Transforms a position from the current volume reference frame
2680 // to the world reference frame.
2681 //
2682 //  Geant3 desription:
2683 //  ==================
2684 //  Computes coordinates XM (Master Reference System
2685 //  knowing the coordinates XD (Detector Ref System)
2686 //  The local reference system can be initialized by
2687 //    - the tracking routines and GDTOM used in GUSTEP
2688 //    - a call to GSCMED(NLEVEL,NAMES,NUMBER)
2689 //        (inverse routine is GMTOD)
2690 // 
2691 //   If IFLAG=1  convert coordinates
2692 //      IFLAG=2  convert direction cosinus
2693 //
2694 // ---
2695    Double_t xmL[3], xdL[3];
2696    Int_t i;
2697    for (i=0;i<3;i++) xdL[i] = xd[i];
2698    if (iflag == 1) gGeoManager->LocalToMaster(xdL,xmL);
2699    else            gGeoManager->LocalToMasterVect(xdL,xmL);
2700    for (i=0;i<3;i++) xm[i]=xmL[i];
2701 }
2702
2703 //______________________________________________________________________________ 
2704 void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
2705 {
2706    if (iflag == 1) gGeoManager->LocalToMaster(xd,xm);
2707    else            gGeoManager->LocalToMasterVect(xd,xm);
2708 }
2709
2710 //______________________________________________________________________________
2711 TObjArray *TFluka::GetFlukaMaterials()
2712 {
2713    return fGeom->GetMatList();
2714 }   
2715
2716 //______________________________________________________________________________
2717 void TFluka::SetMreg(Int_t l) 
2718 {
2719 // Set current fluka region
2720    fCurrentFlukaRegion = l;
2721    fGeom->SetMreg(l);
2722 }
2723
2724
2725
2726
2727 TString TFluka::ParticleName(Int_t pdg) const
2728 {
2729     // Return particle name for particle with pdg code pdg.
2730     Int_t ifluka = IdFromPDG(pdg);
2731     return TString((CHPPRP.btype[ifluka+6]), 8);
2732 }
2733  
2734
2735 Double_t TFluka::ParticleMass(Int_t pdg) const
2736 {
2737     // Return particle mass for particle with pdg code pdg.
2738     Int_t ifluka = IdFromPDG(pdg);
2739     return (PAPROP.am[ifluka+6]);
2740 }
2741
2742 Double_t TFluka::ParticleCharge(Int_t pdg) const
2743 {
2744     // Return particle charge for particle with pdg code pdg.
2745     Int_t ifluka = IdFromPDG(pdg);
2746     return Double_t(PAPROP.ichrge[ifluka+6]);
2747 }
2748
2749 Double_t TFluka::ParticleLifeTime(Int_t pdg) const
2750 {
2751     // Return particle lifetime for particle with pdg code pdg.
2752     Int_t ifluka = IdFromPDG(pdg);
2753     return (PAPROP.thalf[ifluka+6]);
2754 }
2755
2756 void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t& charge, Float_t& tlife)
2757 {
2758     // Retrieve particle properties for particle with pdg code pdg.
2759     
2760     strcpy(name, ParticleName(pdg).Data());
2761     type   = ParticleMCType(pdg);
2762     mass   = ParticleMass(pdg);
2763     charge = ParticleCharge(pdg);
2764     tlife  = ParticleLifeTime(pdg);
2765 }
2766
2767
2768
2769 #define pushcerenkovphoton pushcerenkovphoton_
2770
2771
2772 extern "C" {
2773     void pushcerenkovphoton(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
2774                             Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
2775                             Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
2776     {
2777         //
2778         // Pushes one cerenkov photon to the stack
2779         //
2780         
2781         TFluka* fluka =  (TFluka*) gMC;
2782         TVirtualMCStack* cppstack = fluka->GetStack();
2783         Int_t parent =  TRACKR.ispusr[mkbmx2-1];
2784         cppstack->PushTrack(0, parent, 50000050,
2785                             px, py, pz, e,
2786                             vx, vy, vz, tof,
2787                             polx, poly, polz,
2788                             kPCerenkov, ntr, wgt, 0); 
2789     }
2790 }
2791