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