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