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