1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 #include <Riostream.h>
20 #include "TClonesArray.h"
21 #include "TFlukaGeo.h"
22 #include "TCallf77.h" //For the fortran calls
23 #include "Fdblprc.h" //(DBLPRC) fluka common
24 #include "Fepisor.h" //(EPISOR) fluka common
25 #include "Ffinuc.h" //(FINUC) fluka common
26 #include "Fiounit.h" //(IOUNIT) fluka common
27 #include "Fpaprop.h" //(PAPROP) fluka common
28 #include "Fpart.h" //(PART) fluka common
29 #include "Ftrackr.h" //(TRACKR) fluka common
30 #include "Fpaprop.h" //(PAPROP) fluka common
31 #include "Ffheavy.h" //(FHEAVY) fluka common
33 #include "TVirtualMC.h"
34 #include "TGeoManager.h"
35 #include "TFlukaMCGeometry.h"
37 #include "TLorentzVector.h"
39 // Fluka methods that may be needed.
41 # define flukam flukam_
42 # define fluka_openinp fluka_openinp_
43 # define fluka_closeinp fluka_closeinp_
44 # define mcihad mcihad_
45 # define mpdgha mpdgha_
47 # define flukam FLUKAM
48 # define fluka_openinp FLUKA_OPENINP
49 # define fluka_closeinp FLUKA_CLOSEINP
50 # define mcihad MCIHAD
51 # define mpdgha MPDGHA
57 // Prototypes for FLUKA functions
59 void type_of_call flukam(const int&);
60 void type_of_call fluka_openinp(const int&, DEFCHARA);
61 void type_of_call fluka_closeinp(const int&);
62 int type_of_call mcihad(const int&);
63 int type_of_call mpdgha(const int&);
67 // Class implementation for ROOT
72 //----------------------------------------------------------------------------
73 // TFluka constructors and destructors.
74 //______________________________________________________________________________
81 // Default constructor
84 fCurrentFlukaRegion = -1;
88 //______________________________________________________________________________
89 TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
90 :TVirtualMC("TFluka",title, isRootGeometrySupported),
91 fVerbosityLevel(verbosity),
97 // create geometry interface
98 if (fVerbosityLevel >=3)
99 cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
102 fCurrentFlukaRegion = -1;
103 fGeom = new TFlukaMCGeometry("geom", "ALICE geometry");
106 //______________________________________________________________________________
108 if (fVerbosityLevel >=3)
109 cout << "==> TFluka::~TFluka() destructor called." << endl;
113 if (fVerbosityLevel >=3)
114 cout << "<== TFluka::~TFluka() destructor called." << endl;
118 //______________________________________________________________________________
119 // TFluka control methods
120 //______________________________________________________________________________
121 void TFluka::Init() {
123 if (fVerbosityLevel >=3)
124 cout << "==> TFluka::Init() called." << endl;
126 if (!gGeoManager) new TGeoManager("geom", "FLUKA geometry");
127 fApplication->ConstructGeometry();
128 TGeoVolume *top = (TGeoVolume*)gGeoManager->GetListOfVolumes()->First();
129 gGeoManager->SetTopVolume(top);
130 gGeoManager->CloseGeometry("di");
131 gGeoManager->DefaultColors(); // to be removed
132 fNVolumes = fGeom->NofVolumes();
133 printf("== Number of volumes: %i\n ==", fNVolumes);
134 fGeom->CreateFlukaMatFile("flukaMat.inp");
135 cout << "\t* InitPhysics() - Prepare input file to be called" << endl;
136 // now we have TGeo geometry created and we have to patch alice.inp
137 // with the material mapping file FlukaMat.inp
138 InitPhysics(); // prepare input file with the current physics settings
139 cout << "\t* InitPhysics() - Prepare input file was called" << endl;
141 if (fVerbosityLevel >=2)
142 cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
143 << ") in fluka..." << endl;
144 GLOBAL.lfdrtr = true;
146 if (fVerbosityLevel >=2)
147 cout << "\t* Opening file " << sInputFileName << endl;
148 const char* fname = sInputFileName;
149 fluka_openinp(lunin, PASSCHARA(fname));
151 if (fVerbosityLevel >=2)
152 cout << "\t* Calling flukam..." << endl;
155 if (fVerbosityLevel >=2)
156 cout << "\t* Closing file " << sInputFileName << endl;
157 fluka_closeinp(lunin);
161 if (fVerbosityLevel >=3)
162 cout << "<== TFluka::Init() called." << endl;
165 //______________________________________________________________________________
166 void TFluka::FinishGeometry() {
168 // Build-up table with region to medium correspondance
170 if (fVerbosityLevel >=3)
171 cout << "==> TFluka::FinishGeometry() called." << endl;
173 printf("----FinishGeometry - nothing to do with TGeo\n");
175 if (fVerbosityLevel >=3)
176 cout << "<== TFluka::FinishGeometry() called." << endl;
179 //______________________________________________________________________________
180 void TFluka::BuildPhysics() {
181 if (fVerbosityLevel >=3)
182 cout << "==> TFluka::BuildPhysics() called." << endl;
185 if (fVerbosityLevel >=3)
186 cout << "<== TFluka::BuildPhysics() called." << endl;
189 //______________________________________________________________________________
190 void TFluka::ProcessEvent() {
191 if (fVerbosityLevel >=3)
192 cout << "==> TFluka::ProcessEvent() called." << endl;
193 fApplication->GeneratePrimaries();
194 EPISOR.lsouit = true;
196 if (fVerbosityLevel >=3)
197 cout << "<== TFluka::ProcessEvent() called." << endl;
200 //______________________________________________________________________________
201 void TFluka::ProcessRun(Int_t nevent) {
202 if (fVerbosityLevel >=3)
203 cout << "==> TFluka::ProcessRun(" << nevent << ") called."
206 if (fVerbosityLevel >=2) {
207 cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
208 cout << "\t* Calling flukam again..." << endl;
210 fApplication->InitGeometry();
211 fApplication->BeginEvent();
213 fApplication->FinishEvent();
214 if (fVerbosityLevel >=3)
215 cout << "<== TFluka::ProcessRun(" << nevent << ") called."
220 //_____________________________________________________________________________
221 // methods for building/management of geometry
223 // functions from GCONS
224 //____________________________________________________________________________
225 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
226 Float_t &dens, Float_t &radl, Float_t &absl,
227 Float_t* ubuf, Int_t& nbuf) {
229 fGeom->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
232 //______________________________________________________________________________
233 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
234 Double_t &dens, Double_t &radl, Double_t &absl,
235 Double_t* ubuf, Int_t& nbuf) {
237 fGeom->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
240 // detector composition
241 //______________________________________________________________________________
242 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
243 Double_t z, Double_t dens, Double_t radl, Double_t absl,
244 Float_t* buf, Int_t nwbuf) {
246 fGeom->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf);
249 //______________________________________________________________________________
250 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
251 Double_t z, Double_t dens, Double_t radl, Double_t absl,
252 Double_t* buf, Int_t nwbuf) {
254 fGeom->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf);
257 //______________________________________________________________________________
258 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
259 Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
261 fGeom->Mixture(kmat, name, a, z, dens, nlmat, wmat);
264 //______________________________________________________________________________
265 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
266 Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
268 fGeom->Mixture(kmat, name, a, z, dens, nlmat, wmat);
271 //______________________________________________________________________________
272 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
273 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
274 Double_t stemax, Double_t deemax, Double_t epsil,
275 Double_t stmin, Float_t* ubuf, Int_t nbuf) {
277 fGeom->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
278 epsil, stmin, ubuf, nbuf);
281 //______________________________________________________________________________
282 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
283 Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
284 Double_t stemax, Double_t deemax, Double_t epsil,
285 Double_t stmin, Double_t* ubuf, Int_t nbuf) {
287 fGeom->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
288 epsil, stmin, ubuf, nbuf);
291 //______________________________________________________________________________
292 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
293 Double_t thetaY, Double_t phiY, Double_t thetaZ,
296 fGeom->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
299 //______________________________________________________________________________
300 void TFluka::Gstpar(Int_t /*itmed*/, const char */*param*/, Double_t /*parval*/) {
302 // Is it needed with TGeo ??? - to clear-up
303 Warning("Gstpar", "Not implemented with TGeo");
306 // functions from GGEOM
307 //_____________________________________________________________________________
308 void TFluka::Gsatt(const char *name, const char *att, Int_t val)
310 fGeom->Gsatt(name,att, val);
313 //______________________________________________________________________________
314 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
315 Float_t *upar, Int_t np) {
317 return fGeom->Gsvolu(name, shape, nmed, upar, np);
320 //______________________________________________________________________________
321 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
322 Double_t *upar, Int_t np) {
324 return fGeom->Gsvolu(name, shape, nmed, upar, np);
327 //______________________________________________________________________________
328 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
331 fGeom->Gsdvn(name, mother, ndiv, iaxis);
334 //______________________________________________________________________________
335 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
336 Int_t iaxis, Double_t c0i, Int_t numed) {
338 fGeom->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
341 //______________________________________________________________________________
342 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
343 Int_t iaxis, Int_t numed, Int_t ndvmx) {
345 fGeom->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
348 //______________________________________________________________________________
349 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
350 Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
352 fGeom->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
355 //______________________________________________________________________________
356 void TFluka::Gsord(const char * /*name*/, Int_t /*iax*/) {
358 // Nothing to do with TGeo
361 //______________________________________________________________________________
362 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
363 Double_t x, Double_t y, Double_t z, Int_t irot,
366 fGeom->Gspos(name, nr, mother, x, y, z, irot, konly);
369 //______________________________________________________________________________
370 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
371 Double_t x, Double_t y, Double_t z, Int_t irot,
372 const char *konly, Float_t *upar, Int_t np) {
374 fGeom->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
377 //______________________________________________________________________________
378 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
379 Double_t x, Double_t y, Double_t z, Int_t irot,
380 const char *konly, Double_t *upar, Int_t np) {
382 fGeom->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
385 //______________________________________________________________________________
386 void TFluka::Gsbool(const char* /*onlyVolName*/, const char* /*manyVolName*/) {
388 // Nothing to do with TGeo
389 Warning("Gsbool", "Not implemented with TGeo");
392 //______________________________________________________________________________
393 void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Float_t */*ppckov*/,
394 Float_t * /*absco*/, Float_t * /*effic*/, Float_t * /*rindex*/) {
396 // Not implemented with TGeo - what G4 did ? Any FLUKA card generated?
397 Warning("SetCerenkov", "Not implemented with TGeo");
400 //______________________________________________________________________________
401 void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t * /*ppckov*/,
402 Double_t * /*absco*/, Double_t * /*effic*/, Double_t * /*rindex*/) {
404 // Not implemented with TGeo - what G4 did ? Any FLUKA card generated?
405 Warning("SetCerenkov", "Not implemented with TGeo");
409 //______________________________________________________________________________
410 void TFluka::WriteEuclid(const char* /*fileName*/, const char* /*topVol*/,
411 Int_t /*number*/, Int_t /*nlevel*/) {
414 Warning("WriteEuclid", "Not implemented with TGeo");
419 //_____________________________________________________________________________
420 // methods needed by the stepping
421 //____________________________________________________________________________
423 Int_t TFluka::GetMedium() const {
425 // Get the medium number for the current fluka region
427 return fGeom->GetMedium(); // this I need to check due to remapping !!!
432 //____________________________________________________________________________
433 // particle table usage
434 // ID <--> PDG transformations
435 //_____________________________________________________________________________
436 Int_t TFluka::IdFromPDG(Int_t pdg) const
439 // Return Fluka code from PDG and pseudo ENDF code
441 // Catch the feedback photons
442 if (pdg == 50000051) return (-1);
443 // MCIHAD() goes from pdg to fluka internal.
444 Int_t intfluka = mcihad(pdg);
445 // KPTOIP array goes from internal to official
446 return GetFlukaKPTOIP(intfluka);
449 //______________________________________________________________________________
450 Int_t TFluka::PDGFromId(Int_t id) const
453 // Return PDG code and pseudo ENDF code from Fluka code
455 // IPTOKP array goes from official to internal
459 if (fVerbosityLevel >= 1)
460 printf("\n PDGFromId: Cerenkov Photon \n");
465 if (fVerbosityLevel >= 1)
466 printf("PDGFromId: Error id = 0\n");
470 Int_t intfluka = GetFlukaIPTOKP(id);
472 if (fVerbosityLevel >= 1)
473 printf("PDGFromId: Error intfluka = 0: %d\n", id);
475 } else if (intfluka < 0) {
476 if (fVerbosityLevel >= 1)
477 printf("PDGFromId: Error intfluka < 0: %d\n", id);
480 if (fVerbosityLevel >= 3)
481 printf("mpdgha called with %d %d \n", id, intfluka);
482 // MPDGHA() goes from fluka internal to pdg.
483 return mpdgha(intfluka);
486 //_____________________________________________________________________________
487 // methods for physics management
488 //____________________________________________________________________________
493 //______________________________________________________________________________
494 void TFluka::SetProcess(const char* flagName, Int_t flagValue)
497 if (iNbOfProc < 100) {
498 for (i=0; i<iNbOfProc; i++) {
499 if (strcmp(&sProcessFlag[i][0],flagName) == 0) {
500 iProcessValue[iNbOfProc] = flagValue;
504 strcpy(&sProcessFlag[iNbOfProc][0],flagName);
505 iProcessValue[iNbOfProc++] = flagValue;
508 cout << "Nb of SetProcess calls exceeds 100 - ignored" << endl;
510 iNbOfProc = iNbOfProc;
513 //______________________________________________________________________________
514 void TFluka::SetCut(const char* cutName, Double_t cutValue)
517 if (iNbOfCut < 100) {
518 for (i=0; i<iNbOfCut; i++) {
519 if (strcmp(&sCutFlag[i][0],cutName) == 0) {
520 fCutValue[iNbOfCut] = cutValue;
524 strcpy(&sCutFlag[iNbOfCut][0],cutName);
525 fCutValue[iNbOfCut++] = cutValue;
528 cout << "Nb of SetCut calls exceeds 100 - ignored" << endl;
533 //______________________________________________________________________________
534 Double_t TFluka::Xsec(char*, Double_t, Int_t, Int_t)
536 printf("WARNING: Xsec not yet implemented !\n"); return -1.;
540 //______________________________________________________________________________
541 void TFluka::InitPhysics()
545 Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
546 printf(" last FLUKA material is %g\n", fLastMaterial);
548 // construct file names
549 TString sAliceCoreInp = getenv("ALICE_ROOT");
550 sAliceCoreInp +="/TFluka/input/";
551 TString sAliceTmp = "flukaMat.inp";
552 TString sAliceInp = GetInputFileName();
553 sAliceCoreInp += GetCoreInputFileName();
554 ifstream AliceCoreInp(sAliceCoreInp.Data());
555 ifstream AliceFlukaMat(sAliceTmp.Data());
556 ofstream AliceInp(sAliceInp.Data());
558 // copy core input file
560 Float_t fEventsPerRun;
562 while (AliceCoreInp.getline(sLine,255)) {
563 if (strncmp(sLine,"GEOEND",6) != 0)
564 AliceInp << sLine << endl; // copy until GEOEND card
566 AliceInp << "GEOEND" << endl; // add GEOEND card
569 } // end of while until GEOEND card
572 while (AliceFlukaMat.getline(sLine,255)) { // copy flukaMat.inp file
573 AliceInp << sLine << endl;
576 while (AliceCoreInp.getline(sLine,255)) {
577 if (strncmp(sLine,"START",5) != 0)
578 AliceInp << sLine << endl;
580 sscanf(sLine+10,"%10f",&fEventsPerRun);
583 } //end of while until START card
586 // in G3 the process control values meaning can be different for
587 // different processes, but for most of them is:
588 // 0 process is not activated
589 // 1 process is activated WITH generation of secondaries
590 // 2 process is activated WITHOUT generation of secondaries
591 // if process does not generate secondaries => 1 same as 2
600 // Loop over number of SetProcess calls
601 AliceInp << "*----------------------------------------------------------------------------- ";
603 AliceInp << "*----- The following data are generated from SetProcess and SetCut calls ----- ";
605 AliceInp << "*----------------------------------------------------------------------------- ";
607 for (i=0; i<iNbOfProc; i++) {
610 // G3 default value: 1
611 // G4 processes: G4eplusAnnihilation/G4IeplusAnnihilation
614 // flag = 0 no annihilation
615 // flag = 1 annihilation, decays processed
616 // flag = 2 annihilation, no decay product stored
617 // gMC ->SetProcess("ANNI",1); // EMFCUT -1. 0. 0. 3. lastmat 0. ANNH-THR
618 if (strncmp(&sProcessFlag[i][0],"ANNI",4) == 0) {
619 if (iProcessValue[i] == 1 || iProcessValue[i] == 2) {
622 AliceInp << "*Kinetic energy threshold (GeV) for e+ annihilation - resets to default=0.";
624 AliceInp << "*Generated from call: SetProcess('ANNI',1) or SetProcess('ANNI',2)";
626 AliceInp << setw(10) << "EMFCUT ";
627 AliceInp << setiosflags(ios::scientific) << setprecision(5);
628 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
629 AliceInp << setw(10) << -1.0; // kinetic energy threshold (GeV) for e+ annihilation (resets to default=0)
630 AliceInp << setw(10) << 0.0; // not used
631 AliceInp << setw(10) << 0.0; // not used
632 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
633 AliceInp << setw(10) << setprecision(2);
634 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
635 AliceInp << setprecision(1);
636 AliceInp << setw(10) << 1.0; // step length in assigning indices
637 AliceInp << setw(8) << "ANNH-THR";
640 else if (iProcessValue[i] == 0) {
643 AliceInp << "*No annihilation - no FLUKA card generated";
645 AliceInp << "*Generated from call: SetProcess('ANNI',0)";
651 AliceInp << "*Illegal flag value in SetProcess('ANNI',?) call.";
653 AliceInp << "*No FLUKA card generated";
658 // bremsstrahlung and pair production are both activated
659 // G3 default value: 1
660 // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung,
661 // G4MuBremsstrahlung/G4IMuBremsstrahlung,
662 // G4LowEnergyBremstrahlung
663 // Particles: e-/e+; mu+/mu-
665 // flag = 0 no bremsstrahlung
666 // flag = 1 bremsstrahlung, photon processed
667 // flag = 2 bremsstrahlung, no photon stored
668 // gMC ->SetProcess("BREM",1); // PAIRBREM 2. 0. 0. 3. lastmat
669 // EMFCUT -1. 0. 0. 3. lastmat 0. ELPO-THR
670 // G3 default value: 1
671 // G4 processes: G4GammaConversion,
672 // G4MuPairProduction/G4IMuPairProduction
673 // G4LowEnergyGammaConversion
674 // Particles: gamma, mu
676 // flag = 0 no delta rays
677 // flag = 1 delta rays, secondaries processed
678 // flag = 2 delta rays, no secondaries stored
679 // gMC ->SetProcess("PAIR",1); // PAIRBREM 1. 0. 0. 3. lastmat
680 // EMFCUT 0. 0. -1. 3. lastmat 0. PHOT-THR
681 else if ((strncmp(&sProcessFlag[i][0],"PAIR",4) == 0) && (iProcessValue[i] == 1 || iProcessValue[i] == 2)) {
682 for (j=0; j<iNbOfProc; j++) {
683 if ((strncmp(&sProcessFlag[j][0],"BREM",4) == 0) && (iProcessValue[j] == 1 || iProcessValue[j] == 2)) {
686 AliceInp << "*Bremsstrahlung and pair production by muons and charged hadrons both activated";
688 AliceInp << "*Generated from call: SetProcess('BREM',1) and SetProcess('PAIR',1)";
690 AliceInp << "*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.";
692 AliceInp << "*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.";
694 AliceInp << setw(10) << "PAIRBREM ";
695 AliceInp << setiosflags(ios::scientific) << setprecision(5);
696 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
697 AliceInp << setw(10) << 3.0; // bremsstrahlung and pair production by muons and charged hadrons both are activated
698 // direct pair production by muons
699 // G4 particles: "e-", "e+"
700 // G3 default value: 0.01 GeV
701 //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons
703 for (k=0; k<iNbOfCut; k++) {
704 if (strncmp(&sCutFlag[k][0],"PPCUTM",6) == 0) fCut = fCutValue[k];
706 AliceInp << setiosflags(ios::scientific) << setprecision(5);
707 AliceInp << setw(10) << fCut; // e+, e- kinetic energy threshold (in GeV) for explicit pair production.
708 // muon and hadron bremsstrahlung
709 // G4 particles: "gamma"
710 // G3 default value: CUTGAM=0.001 GeV
711 //gMC ->SetCut("BCUTM",cut); // cut for muon and hadron bremsstrahlung
713 for (k=0; k<iNbOfCut; k++) {
714 if (strncmp(&sCutFlag[k][0],"BCUTM",5) == 0) fCut = fCutValue[k];
716 AliceInp << setiosflags(ios::scientific) << setprecision(5);
717 AliceInp << setw(10) << fCut; // photon energy threshold (GeV) for explicit bremsstrahlung production
718 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
719 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
720 AliceInp << setw(10) << setprecision(2);
721 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
727 AliceInp << "*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.";
729 AliceInp << "*Generated from call: SetProcess('BREM',1);";
731 AliceInp << setw(10) << "EMFCUT ";
733 for (k=0; k<iNbOfCut; k++) {
734 if (strncmp(&sCutFlag[k][0],"BCUTE",5) == 0) fCut = fCutValue[k];
736 AliceInp << setiosflags(ios::scientific) << setprecision(5);
737 AliceInp << setw(10) << fCut; // kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
738 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
739 AliceInp << setw(10) << 0.0; // not used
740 AliceInp << setw(10) << 0.0; // not used
741 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
742 AliceInp << setw(10) << setprecision(2);
743 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
744 AliceInp << setprecision(1);
745 AliceInp << setw(10) << 1.0; // step length in assigning indices
746 AliceInp << setw(8) << "ELPO-THR";
752 AliceInp << "*Pair production by electrons is activated";
754 AliceInp << "*Generated from call: SetProcess('PAIR',1);";
756 AliceInp << setw(10) << "EMFCUT ";
757 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
758 AliceInp << setw(10) << 0.0; // energy threshold (GeV) for Compton scattering (= 0.0 : ignored)
759 AliceInp << setw(10) << 0.0; // energy threshold (GeV) for Photoelectric (= 0.0 : ignored)
761 for (j=0; j<iNbOfCut; j++) {
762 if (strncmp(&sCutFlag[j][0],"CUTGAM",6) == 0) fCut = fCutValue[j];
764 AliceInp << setiosflags(ios::scientific) << setprecision(5);
765 AliceInp << setw(10) << fCut; // energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored)
766 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
767 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
768 AliceInp << setprecision(2);
769 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
770 AliceInp << setprecision(1);
771 AliceInp << setw(10) << 1.0; // step length in assigning indices
772 AliceInp << setw(8) << "PHOT-THR";
775 } // end of if for BREM
776 } // end of loop for BREM
778 // only pair production by muons and charged hadrons is activated
781 AliceInp << "*Pair production by muons and charged hadrons is activated";
783 AliceInp << "*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)";
785 AliceInp << "*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.";
787 AliceInp << setw(10) << "PAIRBREM ";
788 AliceInp << setiosflags(ios::scientific) << setprecision(5);
789 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
790 AliceInp << setw(10) << 1.0; // pair production by muons and charged hadrons is activated
791 // direct pair production by muons
792 // G4 particles: "e-", "e+"
793 // G3 default value: 0.01 GeV
794 //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons
795 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
796 AliceInp << setw(10) << 0.0; // e+, e- kinetic energy threshold (in GeV) for explicit pair production.
797 AliceInp << setw(10) << 0.0; // no explicit bremsstrahlung production is simulated
798 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
799 AliceInp << setprecision(2);
800 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
806 AliceInp << "*Pair production by electrons is activated";
808 AliceInp << "*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)";
810 AliceInp << setw(10) << "EMFCUT ";
811 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
812 AliceInp << setw(10) << 0.0; // energy threshold (GeV) for Compton scattering (= 0.0 : ignored)
813 AliceInp << setw(10) << 0.0; // energy threshold (GeV) for Photoelectric (= 0.0 : ignored)
816 for (j=0; j<iNbOfCut; j++) {
817 if (strncmp(&sCutFlag[j][0],"CUTGAM",6) == 0) fCut = fCutValue[j];
819 AliceInp << setiosflags(ios::scientific) << setprecision(5);
820 AliceInp << setw(10) << fCut; // energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored)
821 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
822 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
823 AliceInp << setprecision(2);
824 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
825 AliceInp << setprecision(1);
826 AliceInp << setw(10) << 1.0; // step length in assigning indices
827 AliceInp << setw(8) << "PHOT-THR";
832 } // end of if for PAIR
837 // G3 default value: 1
838 // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung,
839 // G4MuBremsstrahlung/G4IMuBremsstrahlung,
840 // G4LowEnergyBremstrahlung
841 // Particles: e-/e+; mu+/mu-
843 // flag = 0 no bremsstrahlung
844 // flag = 1 bremsstrahlung, photon processed
845 // flag = 2 bremsstrahlung, no photon stored
846 // gMC ->SetProcess("BREM",1); // PAIRBREM 2. 0. 0. 3. lastmat
847 // EMFCUT -1. 0. 0. 3. lastmat 0. ELPO-THR
848 else if (strncmp(&sProcessFlag[i][0],"BREM",4) == 0) {
849 for (j=0; j<iNbOfProc; j++) {
850 if ((strncmp(&sProcessFlag[j][0],"PAIR",4) == 0) && iProcessValue[j] == 1) goto NOBREM;
852 if (iProcessValue[i] == 1 || iProcessValue[i] == 2) {
855 AliceInp << "*Bremsstrahlung by muons and charged hadrons is activated";
857 AliceInp << "*Generated from call: SetProcess('BREM',1) or SetProcess('BREM',2)";
859 AliceInp << "*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.";
861 AliceInp << setw(10) << "PAIRBREM ";
862 AliceInp << setiosflags(ios::scientific) << setprecision(5);
863 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
864 AliceInp << setw(10) << 2.0; // bremsstrahlung by muons and charged hadrons is activated
865 AliceInp << setw(10) << 0.0; // no meaning
866 // muon and hadron bremsstrahlung
867 // G4 particles: "gamma"
868 // G3 default value: CUTGAM=0.001 GeV
869 //gMC ->SetCut("BCUTM",cut); // cut for muon and hadron bremsstrahlung
871 for (j=0; j<iNbOfCut; j++) {
872 if (strncmp(&sCutFlag[j][0],"BCUTM",5) == 0) fCut = fCutValue[j];
874 AliceInp << setw(10) << fCut; // photon energy threshold (GeV) for explicit bremsstrahlung production
875 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
876 AliceInp << setw(10) << setprecision(2);
877 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
883 AliceInp << "*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.";
885 AliceInp << "*Generated from call: SetProcess('BREM',1);";
887 AliceInp << setw(10) << "EMFCUT ";
888 AliceInp << setiosflags(ios::scientific) << setprecision(5);
889 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
890 AliceInp << setw(10) << -1.0; // kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
891 AliceInp << setw(10) << 0.0; // not used
892 AliceInp << setw(10) << 0.0; // not used
893 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
894 AliceInp << setw(10) << setprecision(2);
895 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
896 AliceInp << setprecision(1);
897 AliceInp << setw(10) << 1.0; // step length in assigning indices
898 AliceInp << setw(8) << "ELPO-THR";
901 else if (iProcessValue[i] == 0) {
904 AliceInp << "*No bremsstrahlung - no FLUKA card generated";
906 AliceInp << "*Generated from call: SetProcess('BREM',0)";
912 AliceInp << "*Illegal flag value in SetProcess('BREM',?) call.";
914 AliceInp << "*No FLUKA card generated";
919 } // end of else if (strncmp(&sProcessFlag[i][0],"BREM",4) == 0)
922 // Cerenkov photon generation
923 // G3 default value: 0
924 // G4 process: G4Cerenkov
926 // Particles: charged
928 // flag = 0 no Cerenkov photon generation
929 // flag = 1 Cerenkov photon generation
930 // flag = 2 Cerenkov photon generation with primary stopped at each step
931 //xx gMC ->SetProcess("CKOV",1); // ??? Cerenkov photon generation
932 else if (strncmp(&sProcessFlag[i][0],"CKOV",4) == 0) {
933 if (iProcessValue[i] == 1 || iProcessValue[i] == 2) {
936 AliceInp << "*Cerenkov photon generation";
938 AliceInp << "*Generated from call: SetProcess('CKOV',1) or SetProcess('CKOV',2)";
940 AliceInp << setw(10) << "OPT-PROD ";
941 AliceInp << setiosflags(ios::scientific) << setprecision(5);
942 AliceInp << setw(10) << 2.07e-9 ; // minimum Cerenkov photon emission energy (in GeV!). Default: 2.07E-9 GeV (corresponding to 600 nm)
943 AliceInp << setw(10) << 4.96e-9; // maximum Cerenkov photon emission energy (in GeV!). Default: 4.96E-9 GeV (corresponding to 250 nm)
944 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
945 AliceInp << setw(10) << 0.0; // not used
946 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
947 AliceInp << setprecision(2);
948 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
949 AliceInp << setprecision(1);
950 AliceInp << setw(10) << 1.0; // step length in assigning indices
951 AliceInp << setw(8) << "CERENKOV";
954 else if (iProcessValue[i] == 0) {
957 AliceInp << "*No Cerenkov photon generation";
959 AliceInp << "*Generated from call: SetProcess('CKOV',0)";
961 AliceInp << setw(10) << "OPT-PROD ";
962 AliceInp << setiosflags(ios::scientific) << setprecision(5);
963 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
964 AliceInp << setw(10) << 0.0; // not used
965 AliceInp << setw(10) << 0.0; // not used
966 AliceInp << setw(10) << 0.0; // not used
967 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
968 AliceInp << setprecision(2);
969 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
970 AliceInp << setprecision(1);
971 AliceInp << setw(10) << 1.0; // step length in assigning indices
972 AliceInp << setw(8) << "CERE-OFF";
978 AliceInp << "*Illegal flag value in SetProcess('CKOV',?) call.";
980 AliceInp << "*No FLUKA card generated";
983 } // end of else if (strncmp(&sProcessFlag[i][0],"CKOV",4) == 0)
986 // Compton scattering
987 // G3 default value: 1
988 // G4 processes: G4ComptonScattering,
989 // G4LowEnergyCompton,
990 // G4PolarizedComptonScattering
993 // flag = 0 no Compton scattering
994 // flag = 1 Compton scattering, electron processed
995 // flag = 2 Compton scattering, no electron stored
996 // gMC ->SetProcess("COMP",1); // EMFCUT -1. 0. 0. 3. lastmat 0. PHOT-THR
997 else if (strncmp(&sProcessFlag[i][0],"COMP",4) == 0) {
998 if (iProcessValue[i] == 1 || iProcessValue[i] == 2) {
1001 AliceInp << "*Energy threshold (GeV) for Compton scattering - resets to default=0.";
1003 AliceInp << "*Generated from call: SetProcess('COMP',1);";
1005 AliceInp << setw(10) << "EMFCUT ";
1006 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1007 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1008 AliceInp << setw(10) << -1.0; // energy threshold (GeV) for Compton scattering - resets to default=0.
1009 AliceInp << setw(10) << 0.0; // not used
1010 AliceInp << setw(10) << 0.0; // not used
1011 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
1012 AliceInp << setprecision(2);
1013 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
1014 AliceInp << setprecision(1);
1015 AliceInp << setw(10) << 1.0; // step length in assigning indices
1016 AliceInp << setw(8) << "PHOT-THR";
1019 else if (iProcessValue[i] == 0) {
1022 AliceInp << "*No Compton scattering - no FLUKA card generated";
1024 AliceInp << "*Generated from call: SetProcess('COMP',0)";
1030 AliceInp << "*Illegal flag value in SetProcess('COMP',?) call.";
1032 AliceInp << "*No FLUKA card generated";
1035 } // end of else if (strncmp(&sProcessFlag[i][0],"COMP",4) == 0)
1038 // G3 default value: 1
1039 // G4 process: G4Decay
1041 // Particles: all which decay is applicable for
1043 // flag = 0 no decays
1044 // flag = 1 decays, secondaries processed
1045 // flag = 2 decays, no secondaries stored
1046 //gMC ->SetProcess("DCAY",1); // not available
1047 else if ((strncmp(&sProcessFlag[i][0],"DCAY",4) == 0) && iProcessValue[i] == 1)
1048 cout << "SetProcess for flag=" << &sProcessFlag[i][0] << " value=" << iProcessValue[i] << " not avaliable!" << endl;
1051 // G3 default value: 2
1052 // !! G4 treats delta rays in different way
1053 // G4 processes: G4eIonisation/G4IeIonization,
1054 // G4MuIonisation/G4IMuIonization,
1055 // G4hIonisation/G4IhIonisation
1056 // Particles: charged
1058 // flag = 0 no energy loss
1059 // flag = 1 restricted energy loss fluctuations
1060 // flag = 2 complete energy loss fluctuations
1061 // flag = 3 same as 1
1062 // flag = 4 no energy loss fluctuations
1063 // gMC ->SetProcess("DRAY",0); // DELTARAY 1.E+6 0. 0. 3. lastmat 0.
1064 else if (strncmp(&sProcessFlag[i][0],"DRAY",4) == 0) {
1065 if (iProcessValue[i] == 0 || iProcessValue[i] == 4) {
1068 AliceInp << "*Kinetic energy threshold (GeV) for delta ray production";
1070 AliceInp << "*Generated from call: SetProcess('DRAY',0) or SetProcess('DRAY',4)";
1072 AliceInp << "*No delta ray production by muons - threshold set artificially high";
1074 AliceInp << setw(10) << "DELTARAY ";
1075 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1076 AliceInp << setw(10) << 1.0e+6; // kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
1077 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1078 AliceInp << setw(10) << 0.0; // ignored
1079 AliceInp << setw(10) << 0.0; // ignored
1080 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
1081 AliceInp << setw(10) << setprecision(2);
1082 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
1083 AliceInp << setprecision(1);
1084 AliceInp << setw(10) << 1.0; // step length in assigning indices
1087 else if (iProcessValue[i] == 1 || iProcessValue[i] == 2 || iProcessValue[i] == 3) {
1090 AliceInp << "*Kinetic energy threshold (GeV) for delta ray production";
1092 AliceInp << "*Generated from call: SetProcess('DRAY',flag), flag=1,2,3";
1094 AliceInp << "*Delta ray production by muons switched on";
1096 AliceInp << "*Energy threshold set by call SetCut('DCUTM',cut) or set to 0.";
1098 AliceInp << setw(10) << "DELTARAY ";
1099 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1101 for (j=0; j<iNbOfCut; j++) {
1102 if (strncmp(&sCutFlag[j][0],"DCUTM",5) == 0) fCut = fCutValue[j];
1104 AliceInp << setw(10) << fCut; // kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
1105 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1106 AliceInp << setw(10) << 0.0; // ignored
1107 AliceInp << setw(10) << 0.0; // ignored
1108 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
1109 AliceInp << setw(10) << setprecision(2);
1110 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
1111 AliceInp << setprecision(1);
1112 AliceInp << setw(10) << 1.0; // step length in assigning indices
1118 AliceInp << "*Illegal flag value in SetProcess('DRAY',?) call.";
1120 AliceInp << "*No FLUKA card generated";
1123 } // end of else if (strncmp(&sProcessFlag[i][0],"DRAY",4) == 0)
1126 // G3 default value: 1
1127 // G4 processes: all defined by TG4PhysicsConstructorHadron
1129 // Particles: hadrons
1131 // flag = 0 no multiple scattering
1132 // flag = 1 hadronic interactions, secondaries processed
1133 // flag = 2 hadronic interactions, no secondaries stored
1134 // gMC ->SetProcess("HADR",1); // ??? hadronic process
1135 //Select pure GEANH (HADR 1) or GEANH/NUCRIN (HADR 3) ?????
1136 else if (strncmp(&sProcessFlag[i][0],"HADR",4) == 0) {
1137 if (iProcessValue[i] == 1 || iProcessValue[i] == 2) {
1140 AliceInp << "*Hadronic interaction is ON by default in FLUKA";
1142 AliceInp << "*No FLUKA card generated";
1145 else if (iProcessValue[i] == 0) {
1148 AliceInp << "*Hadronic interaction is set OFF";
1150 AliceInp << "*Generated from call: SetProcess('HADR',0);";
1152 AliceInp << setw(10) << "MULSOPT ";
1153 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1154 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1155 AliceInp << setw(10) << 0.0; // ignored
1156 AliceInp << setw(10) << 3.0; // multiple scattering for hadrons and muons is completely suppressed
1157 AliceInp << setw(10) << 0.0; // no spin-relativistic corrections
1158 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
1159 AliceInp << setprecision(2);
1160 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
1167 AliceInp << "*Illegal flag value in SetProcess('HADR',?) call.";
1169 AliceInp << "*No FLUKA card generated";
1172 } // end of else if (strncmp(&sProcessFlag[i][0],"HADR",4) == 0)
1176 // G3 default value: 2
1177 // G4 processes: G4eIonisation/G4IeIonization,
1178 // G4MuIonisation/G4IMuIonization,
1179 // G4hIonisation/G4IhIonisation
1181 // Particles: charged
1183 // flag=0 no energy loss
1184 // flag=1 restricted energy loss fluctuations
1185 // flag=2 complete energy loss fluctuations
1187 // flag=4 no energy loss fluctuations
1188 // If the value ILOSS is changed, then (in G3) cross-sections and energy
1189 // loss tables must be recomputed via the command 'PHYSI'
1190 // gMC ->SetProcess("LOSS",2); // ??? IONFLUCT ? energy loss
1191 else if (strncmp(&sProcessFlag[i][0],"LOSS",4) == 0) {
1192 if (iProcessValue[i] == 2) { // complete energy loss fluctuations
1195 AliceInp << "*Complete energy loss fluctuations do not exist in FLUKA";
1197 AliceInp << "*Generated from call: SetProcess('LOSS',2);";
1199 AliceInp << "*flag=2=complete energy loss fluctuations";
1201 AliceInp << "*No input card generated";
1204 else if (iProcessValue[i] == 1 || iProcessValue[i] == 3) { // restricted energy loss fluctuations
1207 AliceInp << "*Restricted energy loss fluctuations";
1209 AliceInp << "*Generated from call: SetProcess('LOSS',1) or SetProcess('LOSS',3)";
1211 AliceInp << setw(10) << "IONFLUCT ";
1212 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1213 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1214 AliceInp << setw(10) << 1.0; // restricted energy loss fluctuations (for hadrons and muons) switched on
1215 AliceInp << setw(10) << 1.0; // restricted energy loss fluctuations (for e+ and e-) switched on
1216 AliceInp << setw(10) << 1.0; // minimal accuracy
1217 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
1218 AliceInp << setprecision(2);
1219 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
1222 else if (iProcessValue[i] == 4) { // no energy loss fluctuations
1225 AliceInp << "*No energy loss fluctuations";
1227 AliceInp << "*Generated from call: SetProcess('LOSS',4)";
1229 AliceInp << setw(10) << -1.0; // restricted energy loss fluctuations (for hadrons and muons) switched off
1230 AliceInp << setw(10) << -1.0; // restricted energy loss fluctuations (for e+ and e-) switched off
1231 AliceInp << setw(10) << 1.0; // minimal accuracy
1232 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
1233 AliceInp << setprecision(2);
1234 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
1240 AliceInp << "*Illegal flag value in SetProcess('LOSS',?) call.";
1242 AliceInp << "*No FLUKA card generated";
1245 } // end of else if (strncmp(&sProcessFlag[i][0],"LOSS",4) == 0)
1248 // multiple scattering
1249 // G3 default value: 1
1250 // G4 process: G4MultipleScattering/G4IMultipleScattering
1252 // Particles: charged
1254 // flag = 0 no multiple scattering
1255 // flag = 1 Moliere or Coulomb scattering
1256 // flag = 2 Moliere or Coulomb scattering
1257 // flag = 3 Gaussian scattering
1258 // gMC ->SetProcess("MULS",1); // MULSOPT multiple scattering
1259 else if (strncmp(&sProcessFlag[i][0],"MULS",4) == 0) {
1260 if (iProcessValue[i] == 1 || iProcessValue[i] == 2 || iProcessValue[i] == 3) {
1263 AliceInp << "*Multiple scattering is ON by default for e+e- and for hadrons/muons";
1265 AliceInp << "*No FLUKA card generated";
1268 else if (iProcessValue[i] == 0) {
1271 AliceInp << "*Multiple scattering is set OFF";
1273 AliceInp << "*Generated from call: SetProcess('MULS',0);";
1275 AliceInp << setw(10) << "MULSOPT ";
1276 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1277 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1278 AliceInp << setw(10) << 0.0; // ignored
1279 AliceInp << setw(10) << 3.0; // multiple scattering for hadrons and muons is completely suppressed
1280 AliceInp << setw(10) << 3.0; // multiple scattering for e+ and e- is completely suppressed
1281 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
1282 AliceInp << setprecision(2);
1283 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
1289 AliceInp << "*Illegal flag value in SetProcess('MULS',?) call.";
1291 AliceInp << "*No FLUKA card generated";
1294 } // end of else if (strncmp(&sProcessFlag[i][0],"MULS",4) == 0)
1297 // muon nuclear interaction
1298 // G3 default value: 0
1299 // G4 processes: G4MuNuclearInteraction,
1300 // G4MuonMinusCaptureAtRest
1304 // flag = 0 no muon-nuclear interaction
1305 // flag = 1 nuclear interaction, secondaries processed
1306 // flag = 2 nuclear interaction, secondaries not processed
1307 // gMC ->SetProcess("MUNU",1); // MUPHOTON 1. 0. 0. 3. lastmat
1308 else if (strncmp(&sProcessFlag[i][0],"MUNU",4) == 0) {
1309 if (iProcessValue[i] == 1) {
1312 AliceInp << "*Muon nuclear interactions with production of secondary hadrons";
1314 AliceInp << "*Generated from call: SetProcess('MUNU',1);";
1316 AliceInp << setw(10) << "MUPHOTON ";
1317 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1318 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1319 AliceInp << setw(10) << 1.0; // full simulation of muon nuclear interactions and production of secondary hadrons
1320 AliceInp << setw(10) << 0.0; // ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
1321 AliceInp << setw(10) << 0.0; // fraction of rho-like interactions ( must be < 1) - Default = 0.75.
1322 AliceInp << setprecision(1);
1323 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
1324 AliceInp << setprecision(2);
1325 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
1328 else if (iProcessValue[i] == 2) {
1331 AliceInp << "*Muon nuclear interactions without production of secondary hadrons";
1333 AliceInp << "*Generated from call: SetProcess('MUNU',2);";
1335 AliceInp << setw(10) << "MUPHOTON ";
1336 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1337 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1338 AliceInp << setw(10) << 2.0; // full simulation of muon nuclear interactions and production of secondary hadrons
1339 AliceInp << setw(10) << 0.0; // ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
1340 AliceInp << setw(10) << 0.0; // fraction of rho-like interactions ( must be < 1) - Default = 0.75.
1341 AliceInp << setprecision(1);
1342 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
1343 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
1346 else if (iProcessValue[i] == 0) {
1349 AliceInp << "*No muon nuclear interaction - no FLUKA card generated";
1351 AliceInp << "*Generated from call: SetProcess('MUNU',0)";
1357 AliceInp << "*Illegal flag value in SetProcess('MUNU',?) call.";
1359 AliceInp << "*No FLUKA card generated";
1362 } // end of else if (strncmp(&sProcessFlag[i][0],"MUNU",4) == 0)
1366 // G3 default value: 0
1371 // gMC ->SetProcess("PFIS",0); // PHOTONUC -1. 0. 0. 3. lastmat 0.
1372 // flag = 0 no photon fission
1373 // flag = 1 photon fission, secondaries processed
1374 // flag = 2 photon fission, no secondaries stored
1375 else if (strncmp(&sProcessFlag[i][0],"PFIS",4) == 0) {
1376 if (iProcessValue[i] == 0) {
1379 AliceInp << "*No photonuclear interactions";
1381 AliceInp << "*Generated from call: SetProcess('PFIS',0);";
1383 AliceInp << setw(10) << "PHOTONUC ";
1384 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1385 AliceInp << setw(10) << -1.0; // no photonuclear interactions
1386 AliceInp << setw(10) << 0.0; // not used
1387 AliceInp << setw(10) << 0.0; // not used
1388 AliceInp << setw(10) << 3.0; // upper bound of the material indices in which the respective thresholds apply
1389 AliceInp << setprecision(2);
1390 AliceInp << setw(10) << fLastMaterial;
1391 AliceInp << setprecision(1); // upper bound of the material indices in which the respective thresholds apply
1392 AliceInp << setprecision(1);
1393 AliceInp << setw(10) << 1.0; // step length in assigning indices
1396 else if (iProcessValue[i] == 1) {
1399 AliceInp << "*Photon nuclear interactions are activated at all energies";
1401 AliceInp << "*Generated from call: SetProcess('PFIS',1);";
1403 AliceInp << setw(10) << "PHOTONUC ";
1404 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1405 AliceInp << setw(10) << 1.0; // photonuclear interactions are activated at all energies
1406 AliceInp << setw(10) << 0.0; // not used
1407 AliceInp << setw(10) << 0.0; // not used
1408 AliceInp << setprecision(2);
1409 AliceInp << setw(10) << 3.0; // upper bound of the material indices in which the respective thresholds apply
1410 AliceInp << setw(10) << fLastMaterial;
1411 AliceInp << setprecision(1); // upper bound of the material indices in which the respective thresholds apply
1412 AliceInp << setprecision(1);
1413 AliceInp << setw(10) << 1.0; // step length in assigning indices
1416 else if (iProcessValue[i] == 0) {
1419 AliceInp << "*No photofission - no FLUKA card generated";
1421 AliceInp << "*Generated from call: SetProcess('PFIS',0)";
1427 AliceInp << "*Illegal flag value in SetProcess('PFIS',?) call.";
1429 AliceInp << "*No FLUKA card generated";
1435 // photo electric effect
1436 // G3 default value: 1
1437 // G4 processes: G4PhotoElectricEffect
1438 // G4LowEnergyPhotoElectric
1441 // flag = 0 no photo electric effect
1442 // flag = 1 photo electric effect, electron processed
1443 // flag = 2 photo electric effect, no electron stored
1444 // gMC ->SetProcess("PHOT",1); // EMFCUT 0. -1. 0. 3. lastmat 0. PHOT-THR
1445 else if (strncmp(&sProcessFlag[i][0],"PHOT",4) == 0) {
1446 if (iProcessValue[i] == 1 || iProcessValue[i] == 2) {
1449 AliceInp << "*Photo electric effect is activated";
1451 AliceInp << "*Generated from call: SetProcess('PHOT',1);";
1453 AliceInp << setw(10) << "EMFCUT ";
1454 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1455 AliceInp << setw(10) << 0.0; // ignored
1456 AliceInp << setw(10) << -1.0; // resets to default=0.
1457 AliceInp << setw(10) << 0.0; // ignored
1458 AliceInp << setw(10) << 3.0; // upper bound of the material indices in which the respective thresholds apply
1459 AliceInp << setprecision(2);
1460 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
1461 AliceInp << setprecision(1);
1462 AliceInp << setw(10) << 1.0; // step length in assigning indices
1463 AliceInp << setw(8) << "PHOT-THR";
1466 else if (iProcessValue[i] == 0) {
1469 AliceInp << "*No photo electric effect - no FLUKA card generated";
1471 AliceInp << "*Generated from call: SetProcess('PHOT',0)";
1477 AliceInp << "*Illegal flag value in SetProcess('PHOT',?) call.";
1479 AliceInp << "*No FLUKA card generated";
1482 } // else if (strncmp(&sProcessFlag[i][0],"PHOT",4) == 0)
1485 // Rayleigh scattering
1486 // G3 default value: 0
1487 // G4 process: G4OpRayleigh
1489 // Particles: optical photon
1491 // flag = 0 Rayleigh scattering off
1492 // flag = 1 Rayleigh scattering on
1493 //xx gMC ->SetProcess("RAYL",1);
1494 else if (strncmp(&sProcessFlag[i][0],"RAYL",4) == 0) {
1495 if (iProcessValue[i] == 1) {
1498 AliceInp << "*Rayleigh scattering is ON by default in FLUKA";
1500 AliceInp << "*No FLUKA card generated";
1503 else if (iProcessValue[i] == 0) {
1506 AliceInp << "*Rayleigh scattering is set OFF";
1508 AliceInp << "*Generated from call: SetProcess('RAYL',0);";
1510 AliceInp << setw(10) << "EMFRAY ";
1511 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1512 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1513 AliceInp << setw(10) << -1.0; // no Rayleigh scattering and no binding corrections for Compton
1514 AliceInp << setw(10) << 3.0; // lower bound of the material indices in which the respective thresholds apply
1515 AliceInp << setprecision(2);
1516 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
1522 AliceInp << "*Illegal flag value in SetProcess('RAYL',?) call.";
1524 AliceInp << "*No FLUKA card generated";
1527 } // end of else if (strncmp(&sProcessFlag[i][0],"RAYL",4) == 0)
1530 // synchrotron radiation in magnetic field
1531 // G3 default value: 0
1532 // G4 process: G4SynchrotronRadiation
1536 // flag = 0 no synchrotron radiation
1537 // flag = 1 synchrotron radiation
1538 //xx gMC ->SetProcess("SYNC",1); // synchrotron radiation generation
1539 else if (strncmp(&sProcessFlag[i][0],"SYNC",4) == 0) {
1542 AliceInp << "*Synchrotron radiation generation is NOT implemented in FLUKA";
1544 AliceInp << "*No FLUKA card generated";
1549 // Automatic calculation of tracking medium parameters
1550 // flag = 0 no automatic calculation
1551 // flag = 1 automatic calculation
1552 //xx gMC ->SetProcess("AUTO",1); // ??? automatic computation of the tracking medium parameters
1553 else if (strncmp(&sProcessFlag[i][0],"AUTO",4) == 0) {
1556 AliceInp << "*Automatic calculation of tracking medium parameters is always ON in FLUKA";
1558 AliceInp << "*No FLUKA card generated";
1563 // To control energy loss fluctuation model
1564 // flag = 0 Urban model
1565 // flag = 1 PAI model
1566 // flag = 2 PAI+ASHO model (not active at the moment)
1567 //xx gMC ->SetProcess("STRA",1); // ??? energy fluctuation model
1568 else if (strncmp(&sProcessFlag[i][0],"STRA",4) == 0) {
1569 if (iProcessValue[i] == 0 || iProcessValue[i] == 2 || iProcessValue[i] == 3) {
1572 AliceInp << "*Ionization energy losses calculation is activated";
1574 AliceInp << "*Generated from call: SetProcess('STRA',n);, n=0,1,2";
1576 AliceInp << setw(10) << "IONFLUCT ";
1577 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1578 AliceInp << setw(10) << 1.0; // restricted energy loss fluctuations
1579 // (for hadrons and muons) switched on
1580 AliceInp << setw(10) << 1.0; // restricted energy loss fluctuations
1581 // (for e+ and e-) switched on
1582 AliceInp << setw(10) << 1.0; // minimal accuracy
1583 AliceInp << setw(10) << 3.0; // upper bound of the material indices in
1584 // which the respective thresholds apply
1585 AliceInp << setprecision(2);
1586 AliceInp << setw(10) << fLastMaterial; // upper bound of the material indices in which the respective thresholds apply
1587 AliceInp << setprecision(1);
1588 AliceInp << setw(10) << 1.0; // step length in assigning indices
1594 AliceInp << "*Illegal flag value in SetProcess('STRA',?) call.";
1596 AliceInp << "*No FLUKA card generated";
1599 } // else if (strncmp(&sProcessFlag[i][0],"STRA",4) == 0)
1604 else { // processes not yet treated
1606 // light photon absorption (Cerenkov photons)
1607 // it is turned on when Cerenkov process is turned on
1608 // G3 default value: 0
1609 // G4 process: G4OpAbsorption, G4OpBoundaryProcess
1611 // Particles: optical photon
1613 // flag = 0 no absorption of Cerenkov photons
1614 // flag = 1 absorption of Cerenkov photons
1615 // gMC ->SetProcess("LABS",2); // ??? Cerenkov light absorption
1619 cout << "SetProcess for flag=" << &sProcessFlag[i][0] << " value=" << iProcessValue[i] << " not yet implemented!" << endl;
1621 } //end of loop number of SetProcess calls
1624 // Loop over number of SetCut calls
1625 for (Int_t i=0; i<iNbOfCut; i++) {
1627 // cuts used in SetProcess calls
1628 if (strncmp(&sCutFlag[i][0],"BCUTM",5) == 0) continue;
1629 else if (strncmp(&sCutFlag[i][0],"BCUTE",5) == 0) continue;
1630 else if (strncmp(&sCutFlag[i][0],"DCUTM",5) == 0) continue;
1631 else if (strncmp(&sCutFlag[i][0],"PPCUTM",6) == 0) continue;
1634 // G4 particles: "gamma"
1635 // G3 default value: 0.001 GeV
1636 //gMC ->SetCut("CUTGAM",cut); // cut for gammas
1637 else if (strncmp(&sCutFlag[i][0],"CUTGAM",6) == 0) {
1640 AliceInp << "*Cut for gamma";
1642 AliceInp << "*Generated from call: SetCut('CUTGAM',cut);";
1644 AliceInp << setw(10) << "PART-THR ";
1645 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1646 AliceInp << setw(10) << -fCutValue[i];
1647 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1648 AliceInp << setw(10) << 7.0;
1653 // G4 particles: "e-"
1655 // G3 default value: 0.001 GeV
1656 //gMC ->SetCut("CUTELE",cut); // cut for e+,e-
1657 else if (strncmp(&sCutFlag[i][0],"CUTELE",6) == 0) {
1660 AliceInp << "*Cut for electrons";
1662 AliceInp << "*Generated from call: SetCut('CUTELE',cut);";
1664 AliceInp << setw(10) << "PART-THR ";
1665 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1666 AliceInp << setw(10) << -fCutValue[i];
1667 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1668 AliceInp << setw(10) << 3.0;
1669 AliceInp << setw(10) << 4.0;
1670 AliceInp << setw(10) << 1.0;
1675 // G4 particles: of type "baryon", "meson", "nucleus" with zero charge
1676 // G3 default value: 0.01 GeV
1677 //gMC ->SetCut("CUTNEU",cut); // cut for neutral hadrons
1678 else if (strncmp(&sCutFlag[i][0],"CUTNEU",6) == 0) {
1681 AliceInp << "*Cut for neutral hadrons";
1683 AliceInp << "*Generated from call: SetCut('CUTNEU',cut);";
1685 AliceInp << setw(10) << "PART-THR ";
1686 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1687 AliceInp << setw(10) << -fCutValue[i];
1688 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1689 AliceInp << setw(10) << 8.0; // Neutron
1690 AliceInp << setw(10) << 9.0; // Antineutron
1692 AliceInp << setw(10) << "PART-THR ";
1693 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1694 AliceInp << setw(10) << -fCutValue[i];
1695 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1696 AliceInp << setw(10) << 12.0; // Kaon zero long
1697 AliceInp << setw(10) << 12.0; // Kaon zero long
1699 AliceInp << setw(10) << "PART-THR ";
1700 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1701 AliceInp << setw(10) << -fCutValue[i];
1702 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1703 AliceInp << setw(10) << 17.0; // Lambda, 18=Antilambda
1704 AliceInp << setw(10) << 19.0; // Kaon zero short
1706 AliceInp << setw(10) << "PART-THR ";
1707 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1708 AliceInp << setw(10) << -fCutValue[i];
1709 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1710 AliceInp << setw(10) << 22.0; // Sigma zero, Pion zero, Kaon zero
1711 AliceInp << setw(10) << 25.0; // Antikaon zero
1713 AliceInp << setw(10) << "PART-THR ";
1714 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1715 AliceInp << setw(10) << -fCutValue[i];
1716 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1717 AliceInp << setw(10) << 32.0; // Antisigma zero
1718 AliceInp << setw(10) << 32.0; // Antisigma zero
1720 AliceInp << setw(10) << "PART-THR ";
1721 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1722 AliceInp << setw(10) << -fCutValue[i];
1723 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1724 AliceInp << setw(10) << 34.0; // Xi zero
1725 AliceInp << setw(10) << 35.0; // AntiXi zero
1727 AliceInp << setw(10) << "PART-THR ";
1728 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1729 AliceInp << setw(10) << -fCutValue[i];
1730 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1731 AliceInp << setw(10) << 47.0; // D zero
1732 AliceInp << setw(10) << 48.0; // AntiD zero
1734 AliceInp << setw(10) << "PART-THR ";
1735 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1736 AliceInp << setw(10) << -fCutValue[i];
1737 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1738 AliceInp << setw(10) << 53.0; // Xi_c zero
1739 AliceInp << setw(10) << 53.0; // Xi_c zero
1741 AliceInp << setw(10) << "PART-THR ";
1742 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1743 AliceInp << setw(10) << -fCutValue[i];
1744 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1745 AliceInp << setw(10) << 55.0; // Xi'_c zero
1746 AliceInp << setw(10) << 56.0; // Omega_c zero
1748 AliceInp << setw(10) << "PART-THR ";
1749 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1750 AliceInp << setw(10) << -fCutValue[i];
1751 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1752 AliceInp << setw(10) << 59.0; // AntiXi_c zero
1753 AliceInp << setw(10) << 59.0; // AntiXi_c zero
1755 AliceInp << setw(10) << "PART-THR ";
1756 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1757 AliceInp << setw(10) << -fCutValue[i];
1758 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1759 AliceInp << setw(10) << 61.0; // AntiXi'_c zero
1760 AliceInp << setw(10) << 62.0; // AntiOmega_c zero
1765 // G4 particles: of type "baryon", "meson", "nucleus" with non-zero charge
1766 // G3 default value: 0.01 GeV
1767 //gMC ->SetCut("CUTHAD",cut); // cut for charged hadrons
1768 else if (strncmp(&sCutFlag[i][0],"CUTHAD",6) == 0) {
1771 AliceInp << "*Cut for charged hadrons";
1773 AliceInp << "*Generated from call: SetCut('CUTHAD',cut);";
1775 AliceInp << setw(10) << "PART-THR ";
1776 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1777 AliceInp << setw(10) << -fCutValue[i];
1778 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1779 AliceInp << setw(10) << 1.0; // Proton
1780 AliceInp << setw(10) << 2.0; // Antiproton
1782 AliceInp << setw(10) << "PART-THR ";
1783 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1784 AliceInp << setw(10) << -fCutValue[i];
1785 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1786 AliceInp << setw(10) << 13.0; // Positive Pion, Negative Pion, Positive Kaon
1787 AliceInp << setw(10) << 16.0; // Negative Kaon
1789 AliceInp << setw(10) << "PART-THR ";
1790 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1791 AliceInp << setw(10) << -fCutValue[i];
1792 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1793 AliceInp << setw(10) << 20.0; // Negative Sigma
1794 AliceInp << setw(10) << 16.0; // Positive Sigma
1796 AliceInp << setw(10) << "PART-THR ";
1797 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1798 AliceInp << setw(10) << -fCutValue[i];
1799 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1800 AliceInp << setw(10) << 31.0; // Antisigma minus
1801 AliceInp << setw(10) << 33.0; // Antisigma plus
1802 AliceInp << setprecision(1);
1803 AliceInp << setw(10) << 2.0; // step length
1805 AliceInp << setw(10) << "PART-THR ";
1806 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1807 AliceInp << setw(10) << -fCutValue[i];
1808 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1809 AliceInp << setw(10) << 36.0; // Negative Xi, Positive Xi, Omega minus
1810 AliceInp << setw(10) << 39.0; // Antiomega
1812 AliceInp << setw(10) << "PART-THR ";
1813 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1814 AliceInp << setw(10) << -fCutValue[i];
1815 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1816 AliceInp << setw(10) << 45.0; // D plus
1817 AliceInp << setw(10) << 46.0; // D minus
1819 AliceInp << setw(10) << "PART-THR ";
1820 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1821 AliceInp << setw(10) << -fCutValue[i];
1822 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1823 AliceInp << setw(10) << 49.0; // D_s plus, D_s minus, Lambda_c plus
1824 AliceInp << setw(10) << 52.0; // Xi_c plus
1826 AliceInp << setw(10) << "PART-THR ";
1827 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1828 AliceInp << setw(10) << -fCutValue[i];
1829 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1830 AliceInp << setw(10) << 54.0; // Xi'_c plus
1831 AliceInp << setw(10) << 60.0; // AntiXi'_c minus
1832 AliceInp << setprecision(1);
1833 AliceInp << setw(10) << 6.0; // step length
1835 AliceInp << setw(10) << "PART-THR ";
1836 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1837 AliceInp << setw(10) << -fCutValue[i];
1838 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(2);
1839 AliceInp << setw(10) << 57.0; // Antilambda_c minus
1840 AliceInp << setw(10) << 58.0; // AntiXi_c minus
1845 // G4 particles: "mu+", "mu-"
1846 // G3 default value: 0.01 GeV
1847 //gMC ->SetCut("CUTMUO",cut); // cut for mu+, mu-
1848 else if (strncmp(&sCutFlag[i][0],"CUTMUO",6) == 0) {
1851 AliceInp << "*Cut for muons";
1853 AliceInp << "*Generated from call: SetCut('CUTMUO',cut);";
1855 AliceInp << setw(10) << "PART-THR ";
1856 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1857 AliceInp << setw(10) << -fCutValue[i];
1858 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1859 AliceInp << setprecision(2);
1860 AliceInp << setw(10) << 10.0;
1861 AliceInp << setw(10) << 11.0;
1864 // delta-rays by electrons
1865 // G4 particles: "e-"
1866 // G3 default value: 10**4 GeV
1867 // gMC ->SetCut("DCUTE",cut); // cut for deltarays by electrons ???????????????
1868 else if (strncmp(&sCutFlag[i][0],"DCUTE",5) == 0) {
1871 AliceInp << "*Cut for delta rays by electrons ????????????";
1873 AliceInp << "*Generated from call: SetCut('DCUTE',cut);";
1875 AliceInp << setw(10) << "EMFCUT ";
1876 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1877 AliceInp << setw(10) << -fCutValue[i];
1878 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1879 AliceInp << setw(10) << 0.0;
1880 AliceInp << setw(10) << 0.0;
1881 AliceInp << setw(10) << 3.0;
1882 AliceInp << setprecision(2);
1883 AliceInp << setw(10) << fLastMaterial;
1884 AliceInp << setprecision(1);
1885 AliceInp << setw(10) << 1.0;
1890 // time of flight cut in seconds
1891 // G4 particles: all
1892 // G3 default value: 0.01 GeV
1893 //gMC ->SetCut("TOFMAX",tofmax); // time of flight cuts in seconds
1894 else if (strncmp(&sCutFlag[i][0],"TOFMAX",6) == 0) {
1897 AliceInp << "*Time of flight cuts in seconds";
1899 AliceInp << "*Generated from call: SetCut('TOFMAX',tofmax);";
1901 AliceInp << setw(10) << "TIME-CUT ";
1902 AliceInp << setiosflags(ios::scientific) << setprecision(5);
1903 AliceInp << setw(10) << fCutValue[i]*1.e9;
1904 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint) << setprecision(1);
1905 AliceInp << setw(10) << 0.0;
1906 AliceInp << setw(10) << 0.0;
1907 AliceInp << setw(10) << -6.0; // lower bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
1908 AliceInp << setprecision(2);
1909 AliceInp << setw(10) << 64.0; // upper bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
1910 AliceInp << setprecision(1);
1911 AliceInp << setw(10) << 1.0; // step length in assigning numbers
1916 cout << "SetCut for flag=" << &sCutFlag[i][0] << " value=" << fCutValue[i] << " not yet implemented!" << endl;
1918 } //end of loop over SeCut calls
1920 // Add START and STOP card
1921 AliceInp << setw(10) << "START ";
1922 AliceInp << setiosflags(ios::fixed) << setiosflags(ios::showpoint);
1923 AliceInp << setw(10) << fEventsPerRun;
1925 AliceInp << setw(10) << "STOP ";
1928 } // end of InitPhysics
1931 //______________________________________________________________________________
1932 void TFluka::SetMaxStep(Double_t)
1934 // SetMaxStep is dummy procedure in TFluka !
1935 if (fVerbosityLevel >=3)
1936 cout << "SetMaxStep is dummy procedure in TFluka !" << endl;
1939 //______________________________________________________________________________
1940 void TFluka::SetMaxNStep(Int_t)
1942 // SetMaxNStep is dummy procedure in TFluka !
1943 if (fVerbosityLevel >=3)
1944 cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
1947 //______________________________________________________________________________
1948 void TFluka::SetUserDecay(Int_t)
1950 // SetUserDecay is dummy procedure in TFluka !
1951 if (fVerbosityLevel >=3)
1952 cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
1956 // dynamic properties
1958 //______________________________________________________________________________
1959 void TFluka::TrackPosition(TLorentzVector& position) const
1961 // Return the current position in the master reference frame of the
1962 // track being transported
1963 // TRACKR.atrack = age of the particle
1964 // TRACKR.xtrack = x-position of the last point
1965 // TRACKR.ytrack = y-position of the last point
1966 // TRACKR.ztrack = z-position of the last point
1967 Int_t caller = GetCaller();
1968 if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
1969 position.SetX(GetXsco());
1970 position.SetY(GetYsco());
1971 position.SetZ(GetZsco());
1972 position.SetT(TRACKR.atrack);
1974 else if (caller == 4) { // mgdraw
1975 position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1976 position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1977 position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1978 position.SetT(TRACKR.atrack);
1980 else if (caller == 5) { // sodraw
1981 position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1982 position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1983 position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1987 Warning("TrackPosition","position not available");
1990 //______________________________________________________________________________
1991 void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
1993 // Return the current position in the master reference frame of the
1994 // track being transported
1995 // TRACKR.atrack = age of the particle
1996 // TRACKR.xtrack = x-position of the last point
1997 // TRACKR.ytrack = y-position of the last point
1998 // TRACKR.ztrack = z-position of the last point
1999 Int_t caller = GetCaller();
2000 if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
2005 else if (caller == 4 || caller == 5) { // mgdraw, sodraw
2006 x = TRACKR.xtrack[TRACKR.ntrack];
2007 y = TRACKR.ytrack[TRACKR.ntrack];
2008 z = TRACKR.ztrack[TRACKR.ntrack];
2011 Warning("TrackPosition","position not available");
2014 //______________________________________________________________________________
2015 void TFluka::TrackMomentum(TLorentzVector& momentum) const
2017 // Return the direction and the momentum (GeV/c) of the track
2018 // currently being transported
2019 // TRACKR.ptrack = momentum of the particle (not always defined, if
2020 // < 0 must be obtained from etrack)
2021 // TRACKR.cx,y,ztrck = direction cosines of the current particle
2022 // TRACKR.etrack = total energy of the particle
2023 // TRACKR.jtrack = identity number of the particle
2024 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
2025 Int_t caller = GetCaller();
2026 if (caller != 2) { // not eedraw
2027 if (TRACKR.ptrack >= 0) {
2028 momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
2029 momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
2030 momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
2031 momentum.SetE(TRACKR.etrack);
2035 Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
2036 momentum.SetPx(p*TRACKR.cxtrck);
2037 momentum.SetPy(p*TRACKR.cytrck);
2038 momentum.SetPz(p*TRACKR.cztrck);
2039 momentum.SetE(TRACKR.etrack);
2044 Warning("TrackMomentum","momentum not available");
2047 //______________________________________________________________________________
2048 void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e) const
2050 // Return the direction and the momentum (GeV/c) of the track
2051 // currently being transported
2052 // TRACKR.ptrack = momentum of the particle (not always defined, if
2053 // < 0 must be obtained from etrack)
2054 // TRACKR.cx,y,ztrck = direction cosines of the current particle
2055 // TRACKR.etrack = total energy of the particle
2056 // TRACKR.jtrack = identity number of the particle
2057 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
2058 Int_t caller = GetCaller();
2059 if (caller != 2) { // not eedraw
2060 if (TRACKR.ptrack >= 0) {
2061 px = TRACKR.ptrack*TRACKR.cxtrck;
2062 py = TRACKR.ptrack*TRACKR.cytrck;
2063 pz = TRACKR.ptrack*TRACKR.cztrck;
2068 Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
2069 px = p*TRACKR.cxtrck;
2070 py = p*TRACKR.cytrck;
2071 pz = p*TRACKR.cztrck;
2077 Warning("TrackMomentum","momentum not available");
2080 //______________________________________________________________________________
2081 Double_t TFluka::TrackStep() const
2083 // Return the length in centimeters of the current step
2084 // TRACKR.ctrack = total curved path
2085 Int_t caller = GetCaller();
2086 if (caller == 11 || caller==12 || caller == 3 || caller == 6) //bxdraw,endraw,usdraw
2088 else if (caller == 4) //mgdraw
2089 return TRACKR.ctrack;
2094 //______________________________________________________________________________
2095 Double_t TFluka::TrackLength() const
2097 // TRACKR.cmtrck = cumulative curved path since particle birth
2098 Int_t caller = GetCaller();
2099 if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
2100 return TRACKR.cmtrck;
2105 //______________________________________________________________________________
2106 Double_t TFluka::TrackTime() const
2108 // Return the current time of flight of the track being transported
2109 // TRACKR.atrack = age of the particle
2110 Int_t caller = GetCaller();
2111 if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
2112 return TRACKR.atrack;
2117 //______________________________________________________________________________
2118 Double_t TFluka::Edep() const
2120 // Energy deposition
2121 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
2122 // -->local energy deposition (the value and the point are not recorded in TRACKR)
2123 // but in the variable "rull" of the procedure "endraw.cxx"
2124 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
2125 // -->no energy loss along the track
2126 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
2127 // -->energy loss distributed along the track
2128 // TRACKR.dtrack = energy deposition of the jth deposition even
2130 // If coming from bxdraw we have 2 steps of 0 length and 0 edep
2131 Int_t caller = GetCaller();
2132 if (caller == 11 || caller==12) return 0.0;
2134 for ( Int_t j=0;j<TRACKR.mtrack;j++) {
2135 sum +=TRACKR.dtrack[j];
2137 if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
2144 //______________________________________________________________________________
2145 Int_t TFluka::TrackPid() const
2147 // Return the id of the particle transported
2148 // TRACKR.jtrack = identity number of the particle
2149 Int_t caller = GetCaller();
2150 if (caller != 2) // not eedraw
2151 return PDGFromId(TRACKR.jtrack);
2156 //______________________________________________________________________________
2157 Double_t TFluka::TrackCharge() const
2159 // Return charge of the track currently transported
2160 // PAPROP.ichrge = electric charge of the particle
2161 // TRACKR.jtrack = identity number of the particle
2162 Int_t caller = GetCaller();
2163 if (caller != 2) // not eedraw
2164 return PAPROP.ichrge[TRACKR.jtrack+6];
2169 //______________________________________________________________________________
2170 Double_t TFluka::TrackMass() const
2172 // PAPROP.am = particle mass in GeV
2173 // TRACKR.jtrack = identity number of the particle
2174 Int_t caller = GetCaller();
2175 if (caller != 2) // not eedraw
2176 return PAPROP.am[TRACKR.jtrack+6];
2181 //______________________________________________________________________________
2182 Double_t TFluka::Etot() const
2184 // TRACKR.etrack = total energy of the particle
2185 Int_t caller = GetCaller();
2186 if (caller != 2) // not eedraw
2187 return TRACKR.etrack;
2195 //______________________________________________________________________________
2196 Bool_t TFluka::IsNewTrack() const
2198 // Return true for the first call of Stepping()
2202 //______________________________________________________________________________
2203 Bool_t TFluka::IsTrackInside() const
2205 // True if the track is not at the boundary of the current volume
2206 // In Fluka a step is always inside one kind of material
2207 // If the step would go behind the region of one material,
2208 // it will be shortened to reach only the boundary.
2209 // Therefore IsTrackInside() is always true.
2210 Int_t caller = GetCaller();
2211 if (caller == 11 || caller==12) // bxdraw
2217 //______________________________________________________________________________
2218 Bool_t TFluka::IsTrackEntering() const
2220 // True if this is the first step of the track in the current volume
2222 Int_t caller = GetCaller();
2223 if (caller == 11) // bxdraw entering
2228 //______________________________________________________________________________
2229 Bool_t TFluka::IsTrackExiting() const
2231 Int_t caller = GetCaller();
2232 if (caller == 12) // bxdraw exiting
2237 //______________________________________________________________________________
2238 Bool_t TFluka::IsTrackOut() const
2240 // True if the track is out of the setup
2242 // Icode = 14: escape - call from Kaskad
2243 // Icode = 23: escape - call from Emfsco
2244 // Icode = 32: escape - call from Kasneu
2245 // Icode = 40: escape - call from Kashea
2246 // Icode = 51: escape - call from Kasoph
2251 fIcode == 51) return 1;
2255 //______________________________________________________________________________
2256 Bool_t TFluka::IsTrackDisappeared() const
2258 // means all inelastic interactions and decays
2259 // fIcode from usdraw
2260 if (fIcode == 101 || // inelastic interaction
2261 fIcode == 102 || // particle decay
2262 fIcode == 214 || // in-flight annihilation
2263 fIcode == 215 || // annihilation at rest
2264 fIcode == 217 || // pair production
2265 fIcode == 221) return 1;
2269 //______________________________________________________________________________
2270 Bool_t TFluka::IsTrackStop() const
2272 // True if the track energy has fallen below the threshold
2273 // means stopped by signal or below energy threshold
2274 // Icode = 12: stopping particle - call from Kaskad
2275 // Icode = 15: time kill - call from Kaskad
2276 // Icode = 21: below threshold, iarg=1 - call from Emfsco
2277 // Icode = 22: below threshold, iarg=2 - call from Emfsco
2278 // Icode = 24: time kill - call from Emfsco
2279 // Icode = 31: below threshold - call from Kasneu
2280 // Icode = 33: time kill - call from Kasneu
2281 // Icode = 41: time kill - call from Kashea
2282 // Icode = 52: time kill - call from Kasoph
2291 fIcode == 52) return 1;
2295 //______________________________________________________________________________
2296 Bool_t TFluka::IsTrackAlive() const
2298 // means not disappeared or not out
2299 if (IsTrackDisappeared() || IsTrackOut() ) return 0;
2307 //______________________________________________________________________________
2308 Int_t TFluka::NSecondaries() const
2309 // Number of secondary particles generated in the current step
2310 // FINUC.np = number of secondaries except light and heavy ions
2311 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
2313 Int_t caller = GetCaller();
2314 if (caller == 6) // valid only after usdraw
2315 return FINUC.np + FHEAVY.npheav;
2318 } // end of NSecondaries
2320 //______________________________________________________________________________
2321 void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
2322 TLorentzVector& position, TLorentzVector& momentum)
2324 Int_t caller = GetCaller();
2325 if (caller == 6) { // valid only after usdraw
2326 if (isec >= 0 && isec < FINUC.np) {
2327 particleId = PDGFromId(FINUC.kpart[isec]);
2328 position.SetX(fXsco);
2329 position.SetY(fYsco);
2330 position.SetZ(fZsco);
2331 position.SetT(TRACKR.atrack);
2332 momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
2333 momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
2334 momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
2335 momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
2337 else if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
2338 Int_t jsec = isec - FINUC.np;
2339 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
2340 position.SetX(fXsco);
2341 position.SetY(fYsco);
2342 position.SetZ(fZsco);
2343 position.SetT(TRACKR.atrack);
2344 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
2345 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
2346 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
2347 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
2348 momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
2349 else if (FHEAVY.tkheav[jsec] > 6)
2350 momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
2353 Warning("GetSecondary","isec out of range");
2356 Warning("GetSecondary","no secondaries available");
2357 } // end of GetSecondary
2359 //______________________________________________________________________________
2360 TMCProcess TFluka::ProdProcess(Int_t) const
2361 // Name of the process that has produced the secondary particles
2362 // in the current step
2364 const TMCProcess kIpNoProc = kPNoProcess;
2365 const TMCProcess kIpPDecay = kPDecay;
2366 const TMCProcess kIpPPair = kPPair;
2367 // const TMCProcess kIpPPairFromPhoton = kPPairFromPhoton;
2368 // const TMCProcess kIpPPairFromVirtualPhoton = kPPairFromVirtualPhoton;
2369 const TMCProcess kIpPCompton = kPCompton;
2370 const TMCProcess kIpPPhotoelectric = kPPhotoelectric;
2371 const TMCProcess kIpPBrem = kPBrem;
2372 // const TMCProcess kIpPBremFromHeavy = kPBremFromHeavy;
2373 // const TMCProcess kIpPBremFromElectronOrPositron = kPBremFromElectronOrPositron;
2374 const TMCProcess kIpPDeltaRay = kPDeltaRay;
2375 // const TMCProcess kIpPMoller = kPMoller;
2376 // const TMCProcess kIpPBhabha = kPBhabha;
2377 const TMCProcess kIpPAnnihilation = kPAnnihilation;
2378 // const TMCProcess kIpPAnnihilInFlight = kPAnnihilInFlight;
2379 // const TMCProcess kIpPAnnihilAtRest = kPAnnihilAtRest;
2380 const TMCProcess kIpPHadronic = kPHadronic;
2381 const TMCProcess kIpPMuonNuclear = kPMuonNuclear;
2382 const TMCProcess kIpPPhotoFission = kPPhotoFission;
2383 const TMCProcess kIpPRayleigh = kPRayleigh;
2384 // const TMCProcess kIpPCerenkov = kPCerenkov;
2385 // const TMCProcess kIpPSynchrotron = kPSynchrotron;
2387 Int_t mugamma = TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11;
2388 if (fIcode == 102) return kIpPDecay;
2389 else if (fIcode == 104 || fIcode == 217) return kIpPPair;
2390 // else if (fIcode == 104) return kIpPairFromPhoton;
2391 // else if (fIcode == 217) return kIpPPairFromVirtualPhoton;
2392 else if (fIcode == 219) return kIpPCompton;
2393 else if (fIcode == 221) return kIpPPhotoelectric;
2394 else if (fIcode == 105 || fIcode == 208) return kIpPBrem;
2395 // else if (fIcode == 105) return kIpPBremFromHeavy;
2396 // else if (fIcode == 208) return kPBremFromElectronOrPositron;
2397 else if (fIcode == 103 || fIcode == 400) return kIpPDeltaRay;
2398 else if (fIcode == 210 || fIcode == 212) return kIpPDeltaRay;
2399 // else if (fIcode == 210) return kIpPMoller;
2400 // else if (fIcode == 212) return kIpPBhabha;
2401 else if (fIcode == 214 || fIcode == 215) return kIpPAnnihilation;
2402 // else if (fIcode == 214) return kIpPAnnihilInFlight;
2403 // else if (fIcode == 215) return kIpPAnnihilAtRest;
2404 else if (fIcode == 101) return kIpPHadronic;
2405 else if (fIcode == 101) {
2406 if (!mugamma) return kIpPHadronic;
2407 else if (TRACKR.jtrack == 7) return kIpPPhotoFission;
2408 else return kIpPMuonNuclear;
2410 else if (fIcode == 225) return kIpPRayleigh;
2411 // Fluka codes 100, 300 and 400 still to be investigasted
2412 else return kIpNoProc;
2415 //Int_t StepProcesses(TArrayI &proc) const
2416 // Return processes active in the current step
2418 //ck = total energy of the particl ????????????????
2422 //______________________________________________________________________________
2423 Int_t TFluka::VolId2Mate(Int_t id) const
2426 // Returns the material number for a given volume ID
2428 return fGeom->VolId2Mate(id);
2431 //______________________________________________________________________________
2432 const char* TFluka::VolName(Int_t id) const
2435 // Returns the volume name for a given volume ID
2437 return fGeom->VolName(id);
2440 //______________________________________________________________________________
2441 Int_t TFluka::VolId(const Text_t* volName) const
2444 // Converts from volume name to volume ID.
2445 // Time consuming. (Only used during set-up)
2446 // Could be replaced by hash-table
2448 return fGeom->VolId(volName);
2451 //______________________________________________________________________________
2452 Int_t TFluka::CurrentVolID(Int_t& copyNo) const
2455 // Return the logical id and copy number corresponding to the current fluka region
2457 return fGeom->CurrentVolID(copyNo);
2460 //______________________________________________________________________________
2461 Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
2464 // Return the logical id and copy number of off'th mother
2465 // corresponding to the current fluka region
2467 return fGeom->CurrentVolOffID(off, copyNo);
2470 //______________________________________________________________________________
2471 const char* TFluka::CurrentVolName() const
2474 // Return the current volume name
2476 return fGeom->CurrentVolName();
2479 //______________________________________________________________________________
2480 const char* TFluka::CurrentVolOffName(Int_t off) const
2483 // Return the volume name of the off'th mother of the current volume
2485 return fGeom->CurrentVolOffName(off);
2488 //______________________________________________________________________________
2489 Int_t TFluka::CurrentMaterial(Float_t & /*a*/, Float_t & /*z*/,
2490 Float_t & /*dens*/, Float_t & /*radl*/, Float_t & /*absl*/) const
2493 // Return the current medium number ??? what about material properties
2496 Int_t id = TFluka::CurrentVolID(copy);
2497 Int_t med = TFluka::VolId2Mate(id);
2501 //______________________________________________________________________________
2502 void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
2504 // Transforms a position from the world reference frame
2505 // to the current volume reference frame.
2507 // Geant3 desription:
2508 // ==================
2509 // Computes coordinates XD (in DRS)
2510 // from known coordinates XM in MRS
2511 // The local reference system can be initialized by
2512 // - the tracking routines and GMTOD used in GUSTEP
2513 // - a call to GMEDIA(XM,NUMED)
2514 // - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
2515 // (inverse routine is GDTOM)
2517 // If IFLAG=1 convert coordinates
2518 // IFLAG=2 convert direction cosinus
2521 fGeom->Gmtod(xm,xd,iflag);
2524 //______________________________________________________________________________
2525 void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
2527 // Transforms a position from the world reference frame
2528 // to the current volume reference frame.
2530 // Geant3 desription:
2531 // ==================
2532 // Computes coordinates XD (in DRS)
2533 // from known coordinates XM in MRS
2534 // The local reference system can be initialized by
2535 // - the tracking routines and GMTOD used in GUSTEP
2536 // - a call to GMEDIA(XM,NUMED)
2537 // - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
2538 // (inverse routine is GDTOM)
2540 // If IFLAG=1 convert coordinates
2541 // IFLAG=2 convert direction cosinus
2544 fGeom->Gmtod(xm,xd,iflag);
2547 //______________________________________________________________________________
2548 void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
2550 // Transforms a position from the current volume reference frame
2551 // to the world reference frame.
2553 // Geant3 desription:
2554 // ==================
2555 // Computes coordinates XM (Master Reference System
2556 // knowing the coordinates XD (Detector Ref System)
2557 // The local reference system can be initialized by
2558 // - the tracking routines and GDTOM used in GUSTEP
2559 // - a call to GSCMED(NLEVEL,NAMES,NUMBER)
2560 // (inverse routine is GMTOD)
2562 // If IFLAG=1 convert coordinates
2563 // IFLAG=2 convert direction cosinus
2566 fGeom->Gdtom(xd,xm,iflag);
2569 //______________________________________________________________________________
2570 void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
2572 // Transforms a position from the current volume reference frame
2573 // to the world reference frame.
2575 // Geant3 desription:
2576 // ==================
2577 // Computes coordinates XM (Master Reference System
2578 // knowing the coordinates XD (Detector Ref System)
2579 // The local reference system can be initialized by
2580 // - the tracking routines and GDTOM used in GUSTEP
2581 // - a call to GSCMED(NLEVEL,NAMES,NUMBER)
2582 // (inverse routine is GMTOD)
2584 // If IFLAG=1 convert coordinates
2585 // IFLAG=2 convert direction cosinus
2588 fGeom->Gdtom(xd,xm,iflag);
2590 //______________________________________________________________________________
2591 void TFluka::SetMreg(Int_t l)
2593 // Set current fluka region
2594 fCurrentFlukaRegion = l;