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