//
#include <Riostream.h>
+#include <TList.h>
#include "TFluka.h"
+#include "TFlukaCodes.h"
#include "TCallf77.h" //For the fortran calls
#include "Fdblprc.h" //(DBLPRC) fluka common
-#include "Fepisor.h" //(EPISOR) fluka common
-#include "Ffinuc.h" //(FINUC) fluka common
+#include "Fsourcm.h" //(SOURCM) fluka common
+#include "Fgenstk.h" //(GENSTK) fluka common
#include "Fiounit.h" //(IOUNIT) fluka common
#include "Fpaprop.h" //(PAPROP) fluka common
#include "Fpart.h" //(PART) fluka common
#include "Fpaprop.h" //(PAPROP) fluka common
#include "Ffheavy.h" //(FHEAVY) fluka common
#include "Fopphst.h" //(OPPHST) fluka common
-#include "Fstack.h" //(STACK) fluka common
+#include "Fflkstk.h" //(FLKSTK) fluka common
#include "Fstepsz.h" //(STEPSZ) fluka common
#include "Fopphst.h" //(OPPHST) fluka common
+#include "Fltclcm.h" //(LTCLCM) fluka common
+#include "Falldlt.h" //(ALLDLT) fluka common
#include "TVirtualMC.h"
#include "TMCProcess.h"
#include "TFlukaScoringOption.h"
#include "TLorentzVector.h"
#include "TArrayI.h"
+#include "TArrayD.h"
+#include "TDatabasePDG.h"
+#include "TStopwatch.h"
+
// Fluka methods that may be needed.
#ifndef WIN32
# define flukam flukam_
# define fluka_openinp fluka_openinp_
+# define fluka_openout fluka_openout_
# define fluka_closeinp fluka_closeinp_
# define mcihad mcihad_
# define mpdgha mpdgha_
+# define newplo newplo_
+# define genout genout_
+# define flkend flkend_
#else
# define flukam FLUKAM
# define fluka_openinp FLUKA_OPENINP
+# define fluka_openout FLUKA_OPENOUT
# define fluka_closeinp FLUKA_CLOSEINP
# define mcihad MCIHAD
# define mpdgha MPDGHA
+# define newplo NEWPLO
+# define genout GENOUT
+# define flkend FLKEND
#endif
extern "C"
// Prototypes for FLUKA functions
//
void type_of_call flukam(const int&);
+ void type_of_call newplo();
+ void type_of_call genout();
+ void type_of_call flkend();
void type_of_call fluka_openinp(const int&, DEFCHARA);
+ void type_of_call fluka_openout(const int&, DEFCHARA);
void type_of_call fluka_closeinp(const int&);
int type_of_call mcihad(const int&);
int type_of_call mpdgha(const int&);
TFluka::TFluka()
:TVirtualMC(),
fVerbosityLevel(0),
+ fNEvent(0),
fInputFileName(""),
- fProcesses(0),
- fCuts(0),
+ fCoreInputFileName(""),
+ fCaller(kNoCaller),
+ fIcode(kNoProcess),
+ fNewReg(-1),
+ fRull(0),
+ fXsco(0),
+ fYsco(0),
+ fZsco(0),
+ fTrackIsEntering(kFALSE),
+ fTrackIsExiting(kFALSE),
+ fTrackIsNew(kFALSE),
+ fFieldFlag(kTRUE),
+ fDummyBoundary(kFALSE),
+ fStopped(kFALSE),
+ fStopEvent(kFALSE),
+ fStopRun(kFALSE),
+ fPrimaryElectronIndex(-1),
+ fMaterials(0),
+ fNVolumes(0),
+ fCurrentFlukaRegion(-1),
+ fNCerenkov(0),
+ fGeom(0),
+ fMCGeo(0),
+ fUserConfig(0),
fUserScore(0)
{
//
// Default constructor
//
- fGeneratePemf = kFALSE;
- fNVolumes = 0;
- fCurrentFlukaRegion = -1;
- fGeom = 0;
- fMCGeo = 0;
- fMaterials = 0;
- fDummyBoundary = 0;
- fFieldFlag = 1;
- fStopped = 0;
- fStopEvent = 0;
- fStopRun = 0;
- fNEvent = 0;
+ for (Int_t i = 0; i < 4; i++) fPint[i] = 0.;
}
//______________________________________________________________________________
TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
:TVirtualMC("TFluka",title, isRootGeometrySupported),
fVerbosityLevel(verbosity),
+ fNEvent(0),
fInputFileName(""),
- fTrackIsEntering(0),
- fTrackIsExiting(0),
- fTrackIsNew(0),
- fProcesses(new TObjArray(100)),
- fCuts(new TObjArray(100)),
+ fCoreInputFileName(""),
+ fCaller(kNoCaller),
+ fIcode(kNoProcess),
+ fNewReg(-1),
+ fRull(0),
+ fXsco(0),
+ fYsco(0),
+ fZsco(0),
+ fTrackIsEntering(kFALSE),
+ fTrackIsExiting(kFALSE),
+ fTrackIsNew(kFALSE),
+ fFieldFlag(kTRUE),
+ fDummyBoundary(kFALSE),
+ fStopped(kFALSE),
+ fStopEvent(kFALSE),
+ fStopRun(kFALSE),
+ fPrimaryElectronIndex(-1),
+ fMaterials(0),
+ fNVolumes(0),
+ fCurrentFlukaRegion(-1),
+ fNCerenkov(0),
+ fGeom(0),
+ fMCGeo(0),
+ fUserConfig(new TObjArray(100)),
fUserScore(new TObjArray(100))
{
// create geometry interface
+ for (Int_t i = 0; i < 4; i++) fPint[i] = 0.;
+
if (fVerbosityLevel >=3)
cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
SetCoreInputFileName();
SetInputFileName();
- SetGeneratePemf(kFALSE);
- fNVolumes = 0;
- fCurrentFlukaRegion = -1;
- fDummyBoundary = 0;
- fFieldFlag = 1;
- fGeneratePemf = kFALSE;
- fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kTRUE);
- fGeom = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
+ fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kFALSE);
+ fGeom = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
if (verbosity > 2) fGeom->SetDebugMode(kTRUE);
- fMaterials = 0;
- fStopped = 0;
- fStopEvent = 0;
- fStopRun = 0;
- fNEvent = 0;
+ PrintHeader();
}
//______________________________________________________________________________
-TFluka::~TFluka() {
-// Destructor
+TFluka::~TFluka()
+{
+ // Destructor
if (fVerbosityLevel >=3)
- cout << "<== TFluka::~TFluka() destructor called." << endl;
+ cout << "<== TFluka::~TFluka() destructor called." << endl;
+ if (fMaterials) delete [] fMaterials;
- delete fGeom;
- delete fMCGeo;
+// delete fGeom;
+// delete fMCGeo;
- if (fCuts) {
- fCuts->Delete();
- delete fCuts;
+ if (fUserConfig) {
+ fUserConfig->Delete();
+ delete fUserConfig;
}
-
- if (fProcesses) {
- fProcesses->Delete();
- delete fProcesses;
+
+ if (fUserScore) {
+ fUserScore->Delete();
+ delete fUserScore;
}
-
-
}
//
if (!gGeoManager) new TGeoManager("geom", "FLUKA geometry");
fApplication->ConstructGeometry();
- TGeoVolume *top = (TGeoVolume*)gGeoManager->GetListOfVolumes()->First();
- gGeoManager->SetTopVolume(top);
- gGeoManager->CloseGeometry("di");
- gGeoManager->DefaultColors(); // to be removed
+ if (!gGeoManager->IsClosed()) {
+ TGeoVolume *top = (TGeoVolume*)gGeoManager->GetListOfVolumes()->First();
+ gGeoManager->SetTopVolume(top);
+ gGeoManager->CloseGeometry("di");
+ } else {
+ TGeoNodeCache *cache = gGeoManager->GetCache();
+ if (!cache->HasIdArray()) {
+ Warning("Init", "Node ID tracking must be enabled with TFluka: enabling...\n");
+ cache->BuildIdArray();
+ }
+ }
fNVolumes = fGeom->NofVolumes();
fGeom->CreateFlukaMatFile("flukaMat.inp");
if (fVerbosityLevel >=3) {
printf("== Number of volumes: %i\n ==", fNVolumes);
cout << "\t* InitPhysics() - Prepare input file to be called" << endl;
- }
- // now we have TGeo geometry created and we have to patch FlukaVmc.inp
- // with the material mapping file FlukaMat.inp
+ }
+
+ fApplication->InitGeometry();
+ fApplication->ConstructOpGeometry();
+ //
+ // Add ions to PDG Data base
+ //
+ AddParticlesToPdgDataBase();
+ //
}
//
if (fVerbosityLevel >=3) {
cout << "==> TFluka::FinishGeometry() called." << endl;
- printf("----FinishGeometry - nothing to do with TGeo\n");
+ printf("----FinishGeometry - applying misalignment if any\n");
cout << "<== TFluka::FinishGeometry() called." << endl;
}
+ TVirtualMCApplication::Instance()->MisalignGeometry();
}
//______________________________________________________________________________
//
if (fVerbosityLevel >=3)
- cout << "==> TFluka::BuildPhysics() called." << endl;
-// Prepare input file with the current physics settings
- InitPhysics();
- cout << "\t* InitPhysics() - Prepare input file was called" << endl;
+ cout << "==> TFluka::BuildPhysics() called." << endl;
+
- if (fVerbosityLevel >=2)
- cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
- << ") in fluka..." << endl;
- GLOBAL.lfdrtr = true;
+ if (fVerbosityLevel >=3) {
+ TList *medlist = gGeoManager->GetListOfMedia();
+ TIter next(medlist);
+ TGeoMedium* med = 0x0;
+ TGeoMaterial* mat = 0x0;
+ Int_t ic = 0;
+
+ while((med = (TGeoMedium*)next()))
+ {
+ mat = med->GetMaterial();
+ printf("Medium %5d %12s %5d %5d\n", ic, (med->GetName()), med->GetId(), mat->GetIndex());
+ ic++;
+ }
+ }
+
+
+ // Prepare input file with the current physics settings
- if (fVerbosityLevel >=2)
- cout << "\t* Opening file " << fInputFileName << endl;
+ InitPhysics();
+// Open fortran files
const char* fname = fInputFileName;
fluka_openinp(lunin, PASSCHARA(fname));
-
- if (fVerbosityLevel >=2)
- cout << "\t* Calling flukam..." << endl;
+ fluka_openout(11, PASSCHARA("fluka.out"));
+// Read input cards
+ cout << "==> TFluka::BuildPhysics() Read input cards." << endl;
+ TStopwatch timer;
+ timer.Start();
+ GLOBAL.lfdrtr = true;
flukam(1);
-
- if (fVerbosityLevel >=2)
- cout << "\t* Closing file " << fInputFileName << endl;
+ cout << "<== TFluka::BuildPhysics() Read input cards End"
+ << Form(" R:%.2fs C:%.2fs", timer.RealTime(),timer.CpuTime()) << endl;
+// Close input file
fluka_closeinp(lunin);
-
+// Finish geometry
FinishGeometry();
-
- if (fVerbosityLevel >=3)
- cout << "<== TFluka::Init() called." << endl;
-
-
- if (fVerbosityLevel >=3)
- cout << "<== TFluka::BuildPhysics() called." << endl;
}
//______________________________________________________________________________
// Process one event
//
if (fStopRun) {
- printf("User Run Abortion: No more events handled !\n");
- fNEvent += 1;
- return;
+ Warning("ProcessEvent", "User Run Abortion: No more events handled !\n");
+ fNEvent += 1;
+ return;
}
if (fVerbosityLevel >=3)
- cout << "==> TFluka::ProcessEvent() called." << endl;
+ cout << "==> TFluka::ProcessEvent() called." << endl;
fApplication->GeneratePrimaries();
- EPISOR.lsouit = true;
+ SOURCM.lsouit = true;
flukam(1);
if (fVerbosityLevel >=3)
- cout << "<== TFluka::ProcessEvent() called." << endl;
+ cout << "<== TFluka::ProcessEvent() called." << endl;
//
// Increase event number
//
//
// Run steering
//
-
+
if (fVerbosityLevel >=3)
cout << "==> TFluka::ProcessRun(" << nevent << ") called."
- << endl;
+ << endl;
if (fVerbosityLevel >=2) {
cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
cout << "\t* Calling flukam again..." << endl;
}
- fApplication->InitGeometry();
Int_t todo = TMath::Abs(nevent);
for (Int_t ev = 0; ev < todo; ev++) {
+ TStopwatch timer;
+ timer.Start();
fApplication->BeginEvent();
ProcessEvent();
fApplication->FinishEvent();
+ cout << "Event: "<< ev
+ << Form(" R:%.2fs C:%.2fs", timer.RealTime(),timer.CpuTime()) << endl;
}
if (fVerbosityLevel >=3)
cout << "<== TFluka::ProcessRun(" << nevent << ") called."
- << endl;
+ << endl;
+
+ // Write fluka specific scoring output
+ genout();
+ newplo();
+ flkend();
+
return kTRUE;
}
// functions from GCONS
//____________________________________________________________________________
void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
- Float_t &dens, Float_t &radl, Float_t &absl,
- Float_t* /*ubuf*/, Int_t& /*nbuf*/) {
+ Float_t &dens, Float_t &radl, Float_t &absl,
+ Float_t* /*ubuf*/, Int_t& /*nbuf*/) {
//
TGeoMaterial *mat;
TIter next (gGeoManager->GetListOfMaterials());
//______________________________________________________________________________
void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
- Double_t &dens, Double_t &radl, Double_t &absl,
- Double_t* /*ubuf*/, Int_t& /*nbuf*/) {
+ Double_t &dens, Double_t &radl, Double_t &absl,
+ Double_t* /*ubuf*/, Int_t& /*nbuf*/) {
//
TGeoMaterial *mat;
TIter next (gGeoManager->GetListOfMaterials());
// detector composition
//______________________________________________________________________________
void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
- Double_t z, Double_t dens, Double_t radl, Double_t absl,
- Float_t* buf, Int_t nwbuf) {
+ Double_t z, Double_t dens, Double_t radl, Double_t absl,
+ Float_t* buf, Int_t nwbuf) {
//
Double_t* dbuf = fGeom->CreateDoubleArray(buf, nwbuf);
Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf);
//______________________________________________________________________________
void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
- Double_t z, Double_t dens, Double_t radl, Double_t absl,
- Double_t* /*buf*/, Int_t /*nwbuf*/) {
+ Double_t z, Double_t dens, Double_t radl, Double_t absl,
+ Double_t* /*buf*/, Int_t /*nwbuf*/) {
//
+// Define a material
TGeoMaterial *mat;
kmat = gGeoManager->GetListOfMaterials()->GetSize();
if ((z-Int_t(z)) > 1E-3) {
//______________________________________________________________________________
void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
- Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
+ Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
+//
+// Define a material mixture
//
Double_t* da = fGeom->CreateDoubleArray(a, TMath::Abs(nlmat));
Double_t* dz = fGeom->CreateDoubleArray(z, TMath::Abs(nlmat));
//______________________________________________________________________________
void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
- Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
+ Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
//
// Defines mixture OR COMPOUND IMAT as composed by
// THE BASIC NLMAT materials defined by arrays A,Z and WMAT
// In this case, WMAT in output is changed to relative
// weigths.
//
+ printf("Mixture %5d %10s %5d \n", kmat, name, nlmat);
+
Int_t i,j;
if (nlmat < 0) {
nlmat = - nlmat;
if (!mat->IsMixture()) continue;
mix = (TGeoMixture*)mat;
if (TMath::Abs(z[i]-mix->GetZ()) >1E-3) continue;
-// printf(" FOUND component %i as mixture %s\n", i, mat->GetName());
mixnew = kTRUE;
break;
}
delete [] wmatnew;
return;
}
- gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
+ printf("Mixture (2) %5d %10s %5d \n", kmat, name, nlmat);
+ gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
}
//______________________________________________________________________________
void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
- Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
- Double_t stemax, Double_t deemax, Double_t epsil,
- Double_t stmin, Float_t* ubuf, Int_t nbuf) {
- //
+ Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
+ Double_t stemax, Double_t deemax, Double_t epsil,
+ Double_t stmin, Float_t* ubuf, Int_t nbuf) {
+ // Define a medium
+ //
kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
- epsil, stmin, ubuf, nbuf);
+ epsil, stmin, ubuf, nbuf);
}
//______________________________________________________________________________
void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
- Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
- Double_t stemax, Double_t deemax, Double_t epsil,
- Double_t stmin, Double_t* ubuf, Int_t nbuf) {
- //
+ Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
+ Double_t stemax, Double_t deemax, Double_t epsil,
+ Double_t stmin, Double_t* ubuf, Int_t nbuf) {
+ // Define a medium
+ //
kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
- epsil, stmin, ubuf, nbuf);
+ epsil, stmin, ubuf, nbuf);
}
//______________________________________________________________________________
void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
- Double_t thetaY, Double_t phiY, Double_t thetaZ,
- Double_t phiZ) {
-//
+ Double_t thetaY, Double_t phiY, Double_t thetaZ,
+ Double_t phiZ) {
+//
krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
fMCGeo->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
}
void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
//
//
-// Check if material is used
- if (fVerbosityLevel >=3)
- printf("Gstpar called with %6d %5s %12.4e %6d\n", itmed, param, parval, fGeom->GetFlukaMaterial(itmed));
- Int_t* reglist;
- Int_t nreg;
- reglist = fGeom->GetMaterialList(itmed, nreg);
- if (nreg == 0) return;
//
Bool_t process = kFALSE;
+ Bool_t modelp = kFALSE;
+
if (strncmp(param, "DCAY", 4) == 0 ||
strncmp(param, "PAIR", 4) == 0 ||
strncmp(param, "COMP", 4) == 0 ||
{
process = kTRUE;
}
+
+ if (strncmp(param, "PRIMIO_N", 8) == 0 ||
+ strncmp(param, "PRIMIO_E", 8) == 0)
+ {
+ modelp = kTRUE;
+ }
+
if (process) {
- SetProcess(param, Int_t (parval), fGeom->GetFlukaMaterial(itmed));
+ // Process switch
+ SetProcess(param, Int_t (parval), itmed);
+ } else if (modelp) {
+ // Model parameters
+ SetModelParameter(param, parval, itmed);
} else {
- SetCut(param, parval, fGeom->GetFlukaMaterial(itmed));
+ // Cuts
+ SetCut(param, parval, itmed);
}
+
+
}
// functions from GGEOM
//______________________________________________________________________________
Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
- Float_t *upar, Int_t np) {
+ Float_t *upar, Int_t np) {
//
return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
}
//______________________________________________________________________________
Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
- Double_t *upar, Int_t np) {
+ Double_t *upar, Int_t np) {
//
return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
}
//______________________________________________________________________________
void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
- Int_t iaxis) {
+ Int_t iaxis) {
//
fMCGeo->Gsdvn(name, mother, ndiv, iaxis);
}
//______________________________________________________________________________
void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
- Int_t iaxis, Double_t c0i, Int_t numed) {
+ Int_t iaxis, Double_t c0i, Int_t numed) {
//
fMCGeo->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
}
//______________________________________________________________________________
void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
- Int_t iaxis, Int_t numed, Int_t ndvmx) {
-//
+ Int_t iaxis, Int_t numed, Int_t ndvmx) {
+//
fMCGeo->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
}
//______________________________________________________________________________
void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
- Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
+ Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
//
fMCGeo->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
}
//______________________________________________________________________________
void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
- Double_t x, Double_t y, Double_t z, Int_t irot,
- const char *konly) {
+ Double_t x, Double_t y, Double_t z, Int_t irot,
+ const char *konly) {
//
fMCGeo->Gspos(name, nr, mother, x, y, z, irot, konly);
}
//______________________________________________________________________________
void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
- Double_t x, Double_t y, Double_t z, Int_t irot,
- const char *konly, Float_t *upar, Int_t np) {
+ Double_t x, Double_t y, Double_t z, Int_t irot,
+ const char *konly, Float_t *upar, Int_t np) {
//
fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
}
//______________________________________________________________________________
void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
- Double_t x, Double_t y, Double_t z, Int_t irot,
- const char *konly, Double_t *upar, Int_t np) {
+ Double_t x, Double_t y, Double_t z, Int_t irot,
+ const char *konly, Double_t *upar, Int_t np) {
//
fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
}
// Nothing to do with TGeo
}
+//______________________________________________________________________
+Bool_t TFluka::GetTransformation(const TString &volumePath,TGeoHMatrix &mat)
+{
+ // Returns the Transformation matrix between the volume specified
+ // by the path volumePath and the Top or mater volume. The format
+ // of the path volumePath is as follows (assuming ALIC is the Top volume)
+ // "/ALIC_1/DDIP_1/S05I_2/S05H_1/S05G_3". Here ALIC is the top most
+ // or master volume which has only 1 instance of. Of all of the daughter
+ // volumes of ALICE, DDIP volume copy #1 is indicated. Similarly for
+ // the daughter volume of DDIP is S05I copy #2 and so on.
+ // Inputs:
+ // TString& volumePath The volume path to the specific volume
+ // for which you want the matrix. Volume name
+ // hierarchy is separated by "/" while the
+ // copy number is appended using a "_".
+ // Outputs:
+ // TGeoHMatrix &mat A matrix with its values set to those
+ // appropriate to the Local to Master transformation
+ // Return:
+ // A logical value if kFALSE then an error occurred and no change to
+ // mat was made.
+
+ // We have to preserve the modeler state
+ return fMCGeo->GetTransformation(volumePath, mat);
+}
+
+//______________________________________________________________________
+Bool_t TFluka::GetShape(const TString &volumePath,TString &shapeType,
+ TArrayD &par)
+{
+ // Returns the shape and its parameters for the volume specified
+ // by volumeName.
+ // Inputs:
+ // TString& volumeName The volume name
+ // Outputs:
+ // TString &shapeType Shape type
+ // TArrayD &par A TArrayD of parameters with all of the
+ // parameters of the specified shape.
+ // Return:
+ // A logical indicating whether there was an error in getting this
+ // information
+ return fMCGeo->GetShape(volumePath, shapeType, par);
+}
+
+//______________________________________________________________________
+Bool_t TFluka::GetMaterial(const TString &volumeName,
+ TString &name,Int_t &imat,
+ Double_t &a,Double_t &z,Double_t &dens,
+ Double_t &radl,Double_t &inter,TArrayD &par)
+{
+ // Returns the Material and its parameters for the volume specified
+ // by volumeName.
+ // Note, Geant3 stores and uses mixtures as an element with an effective
+ // Z and A. Consequently, if the parameter Z is not integer, then
+ // this material represents some sort of mixture.
+ // Inputs:
+ // TString& volumeName The volume name
+ // Outputs:
+ // TSrting &name Material name
+ // Int_t &imat Material index number
+ // Double_t &a Average Atomic mass of material
+ // Double_t &z Average Atomic number of material
+ // Double_t &dens Density of material [g/cm^3]
+ // Double_t &radl Average radiation length of material [cm]
+ // Double_t &inter Average interaction length of material [cm]
+ // TArrayD &par A TArrayD of user defined parameters.
+ // Return:
+ // kTRUE if no errors
+ return fMCGeo->GetMaterial(volumeName,name,imat,a,z,dens,radl,inter,par);
+}
+
+//______________________________________________________________________
+Bool_t TFluka::GetMedium(const TString &volumeName,TString &name,
+ Int_t &imed,Int_t &nmat,Int_t &isvol,Int_t &ifield,
+ Double_t &fieldm,Double_t &tmaxfd,Double_t &stemax,
+ Double_t &deemax,Double_t &epsil, Double_t &stmin,
+ TArrayD &par)
+{
+ // Returns the Medium and its parameters for the volume specified
+ // by volumeName.
+ // Inputs:
+ // TString& volumeName The volume name.
+ // Outputs:
+ // TString &name Medium name
+ // Int_t &nmat Material number defined for this medium
+ // Int_t &imed The medium index number
+ // Int_t &isvol volume number defined for this medium
+ // Int_t &iflield Magnetic field flag
+ // Double_t &fieldm Magnetic field strength
+ // Double_t &tmaxfd Maximum angle of deflection per step
+ // Double_t &stemax Maximum step size
+ // Double_t &deemax Maximum fraction of energy allowed to be lost
+ // to continuous process.
+ // Double_t &epsil Boundary crossing precision
+ // Double_t &stmin Minimum step size allowed
+ // TArrayD &par A TArrayD of user parameters with all of the
+ // parameters of the specified medium.
+ // Return:
+ // kTRUE if there where no errors
+ return fMCGeo->GetMedium(volumeName,name,imed,nmat,isvol,ifield,fieldm,tmaxfd,stemax,deemax,epsil,stmin,par);
+}
+
//______________________________________________________________________________
void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
- Float_t* absco, Float_t* effic, Float_t* rindex) {
+ Float_t* absco, Float_t* effic, Float_t* rindex) {
//
// Set Cerenkov properties for medium itmed
//
//
//
// Create object holding Cerenkov properties
-//
+//
+
TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex);
//
// Pass object to medium
medium->SetCerenkovProperties(cerenkovProperties);
}
+void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
+ Float_t* absco, Float_t* effic, Float_t* rindex, Float_t* rfl) {
+//
+// Set Cerenkov properties for medium itmed
+//
+// npckov: number of sampling points
+// ppckov: energy values
+// absco: absorption length
+// effic: quantum efficiency
+// rindex: refraction index
+// rfl: reflectivity for boundary to medium itmed
+//
+//
+// Create object holding Cerenkov properties
+//
+ TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex, rfl);
+//
+// Pass object to medium
+ TGeoMedium* medium = gGeoManager->GetMedium(itmed);
+ medium->SetCerenkovProperties(cerenkovProperties);
+}
+
+
//______________________________________________________________________________
-void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t * /*ppckov*/,
- Double_t * /*absco*/, Double_t * /*effic*/, Double_t * /*rindex*/) {
+void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
+ Double_t *absco, Double_t *effic, Double_t *rindex) {
+//
+// Set Cerenkov properties for medium itmed
//
-// Not implemented with TGeo - what G4 did ? Any FLUKA card generated?
- Warning("SetCerenkov", "Not implemented with TGeo");
+// npckov: number of sampling points
+// ppckov: energy values
+// absco: absorption length
+// effic: quantum efficiency
+// rindex: refraction index
+//
+
+//
+// Double_t version
+ Float_t* fppckov = CreateFloatArray(ppckov, npckov);
+ Float_t* fabsco = CreateFloatArray(absco, npckov);
+ Float_t* feffic = CreateFloatArray(effic, npckov);
+ Float_t* frindex = CreateFloatArray(rindex, npckov);
+
+ SetCerenkov(itmed, npckov, fppckov, fabsco, feffic, frindex);
+
+ delete [] fppckov;
+ delete [] fabsco;
+ delete [] feffic;
+ delete [] frindex;
}
-
+
+void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t* ppckov,
+ Double_t* absco, Double_t* effic, Double_t* rindex, Double_t* rfl) {
+//
+// Set Cerenkov properties for medium itmed
+//
+// npckov: number of sampling points
+// ppckov: energy values
+// absco: absorption length
+// effic: quantum efficiency
+// rindex: refraction index
+// rfl: reflectivity for boundary to medium itmed
+//
+
+//
+// // Double_t version
+ Float_t* fppckov = CreateFloatArray(ppckov, npckov);
+ Float_t* fabsco = CreateFloatArray(absco, npckov);
+ Float_t* feffic = CreateFloatArray(effic, npckov);
+ Float_t* frindex = CreateFloatArray(rindex, npckov);
+ Float_t* frfl = CreateFloatArray(rfl, npckov);
+
+ SetCerenkov(itmed, npckov, fppckov, fabsco, feffic, frindex, frfl);
+
+ delete [] fppckov;
+ delete [] fabsco;
+ delete [] feffic;
+ delete [] frindex;
+ delete [] frfl;
+}
+
// Euclid
//______________________________________________________________________________
void TFluka::WriteEuclid(const char* /*fileName*/, const char* /*topVol*/,
Int_t /*number*/, Int_t /*nlevel*/) {
//
// Not with TGeo
- Warning("WriteEuclid", "Not implemented with TGeo");
+ Warning("WriteEuclid", "Not implemented !");
}
//
// Get the medium number for the current fluka region
//
- return fGeom->GetMedium(); // this I need to check due to remapping !!!
+ if (gGeoManager->IsOutside()) {
+ return (-1);
+ } else {
+ return (fGeom->GetMedium()); // this I need to check due to remapping !!!
+ }
}
+//____________________________________________________________________________
+Int_t TFluka::GetDummyRegion() const
+{
+// Returns index of the dummy region.
+ return fGeom->GetDummyRegion();
+}
+//____________________________________________________________________________
+Int_t TFluka::GetDummyLattice() const
+{
+// Returns index of the dummy lattice.
+ return fGeom->GetDummyLattice();
+}
//____________________________________________________________________________
// particle table usage
// Return Fluka code from PDG and pseudo ENDF code
// Catch the feedback photons
- if (pdg == 50000051) return (-1);
+ if (pdg == 50000051) return (kFLUKAoptical);
// MCIHAD() goes from pdg to fluka internal.
Int_t intfluka = mcihad(pdg);
// KPTOIP array goes from internal to official
//
// Return PDG code and pseudo ENDF code from Fluka code
// Alpha He3 Triton Deuteron gen. ion opt. photon
- Int_t idSpecial[6] = {10020040, 10020030, 10010030, 10010020, 10000000, 50000050};
+ Int_t idSpecial[6] = {GetIonPdg(2,4), GetIonPdg(2, 3), GetIonPdg(1,3), GetIonPdg(1,2), GetIonPdg(0,0), 50000050};
// IPTOKP array goes from official to internal
- if (id == -1) {
+ if (id == kFLUKAoptical) {
// Cerenkov photon
- if (fVerbosityLevel >= 3)
- printf("\n PDGFromId: Cerenkov Photon \n");
- return 50000050;
+// if (fVerbosityLevel >= 3)
+// printf("\n PDGFromId: Cerenkov Photon \n");
+ return 50000050;
}
// Error id
- if (id == 0 || id < -6 || id > 250) {
- if (fVerbosityLevel >= 3)
- printf("PDGFromId: Error id = 0\n");
- return -1;
+ if (id == 0 || id < kFLUKAcodemin || id > kFLUKAcodemax) {
+ if (fVerbosityLevel >= 3)
+ printf("PDGFromId: Error id = 0 %5d %5d\n", id, fCaller);
+ return -1;
}
// Good id
if (id > 0) {
- Int_t intfluka = GetFlukaIPTOKP(id);
- if (intfluka == 0) {
- if (fVerbosityLevel >= 3)
- printf("PDGFromId: Error intfluka = 0: %d\n", id);
- return -1;
- } else if (intfluka < 0) {
- if (fVerbosityLevel >= 3)
- printf("PDGFromId: Error intfluka < 0: %d\n", id);
- return -1;
- }
- if (fVerbosityLevel >= 3)
- printf("mpdgha called with %d %d \n", id, intfluka);
- // MPDGHA() goes from fluka internal to pdg.
- return mpdgha(intfluka);
+ Int_t intfluka = GetFlukaIPTOKP(id);
+ if (intfluka == 0) {
+ if (fVerbosityLevel >= 3)
+ printf("PDGFromId: Error intfluka = 0: %d\n", id);
+ return -1;
+ } else if (intfluka < 0) {
+ if (fVerbosityLevel >= 3)
+ printf("PDGFromId: Error intfluka < 0: %d\n", id);
+ return -1;
+ }
+// if (fVerbosityLevel >= 3)
+// printf("mpdgha called with %d %d \n", id, intfluka);
+ return mpdgha(intfluka);
} else {
- // ions and optical photons
- return idSpecial[id + 6];
+ // ions and optical photons
+ return idSpecial[id - kFLUKAcodemin];
}
}
{
// Set process user flag for material imat
//
- TFlukaConfigOption* proc = new TFlukaConfigOption(flagName, flagValue, imed);
- fProcesses->Add(proc);
+//
+// Update if already in the list
+//
+ TIter next(fUserConfig);
+ TFlukaConfigOption* proc;
+ while((proc = (TFlukaConfigOption*)next()))
+ {
+ if (proc->Medium() == imed) {
+ proc->SetProcess(flagName, flagValue);
+ return;
+ }
+ }
+ proc = new TFlukaConfigOption(imed);
+ proc->SetProcess(flagName, flagValue);
+ fUserConfig->Add(proc);
}
//______________________________________________________________________________
// Set process user flag
//
//
-// Update if already in the list
-//
+ SetProcess(flagName, flagValue, -1);
+ return kTRUE;
+}
- TIter next(fProcesses);
+//______________________________________________________________________________
+void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
+{
+// Set user cut value for material imed
+//
+ TIter next(fUserConfig);
TFlukaConfigOption* proc;
while((proc = (TFlukaConfigOption*)next()))
{
- if (strcmp(proc->GetName(), flagName) == 0) {
- proc->SetFlag(flagValue);
- proc->SetMedium(-1);
- return kTRUE;
- }
+ if (proc->Medium() == imed) {
+ proc->SetCut(cutName, cutValue);
+ return;
+ }
}
-//
-// If not create a new process
-//
- proc = new TFlukaConfigOption(flagName, flagValue);
- fProcesses->Add(proc);
-
- return kTRUE;
+ proc = new TFlukaConfigOption(imed);
+ proc->SetCut(cutName, cutValue);
+ fUserConfig->Add(proc);
}
+
//______________________________________________________________________________
-void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
+void TFluka::SetModelParameter(const char* parName, Double_t parValue, Int_t imed)
{
-// Set user cut value for material imed
+// Set model parameter for material imed
//
- printf("TFluka::SetCut %s %e %d \n", cutName, cutValue, imed);
-
- TFlukaConfigOption* cut = new TFlukaConfigOption(cutName, cutValue, imed);
- fCuts->Add(cut);
+ TIter next(fUserConfig);
+ TFlukaConfigOption* proc;
+ while((proc = (TFlukaConfigOption*)next()))
+ {
+ if (proc->Medium() == imed) {
+ proc->SetModelParameter(parName, parValue);
+ return;
+ }
+ }
+
+ proc = new TFlukaConfigOption(imed);
+ proc->SetModelParameter(parName, parValue);
+ fUserConfig->Add(proc);
}
//______________________________________________________________________________
// Set user cut value
//
//
-// Update if already in the list
-//
+ SetCut(cutName, cutValue, -1);
+ return kTRUE;
+}
- TIter next(fCuts);
- TFlukaConfigOption* cut;
- while((cut = (TFlukaConfigOption*)next()))
- {
- if (strcmp(cut->GetName(), cutName) == 0) {
- cut->SetCut(cutValue);
- return kTRUE;
- }
- }
-//
-// If not create a new process
-//
- cut = new TFlukaConfigOption(cutName, cutValue);
- fCuts->Add(cut);
-
- return kTRUE;
+void TFluka::SetUserScoring(const char* option, const char* sdum, Int_t npr, char* outfile, Float_t* what)
+{
+//
+// Adds a user scoring option to the list
+//
+ TFlukaScoringOption* opt = new TFlukaScoringOption(option, sdum, npr,outfile,what);
+ fUserScore->Add(opt);
}
-
-void TFluka::SetUserScoring(const char* option, Int_t npar, Float_t what[12])
+//______________________________________________________________________________
+void TFluka::SetUserScoring(const char* option, const char* sdum, Int_t npr, char* outfile, Float_t* what,
+ const char* det1, const char* det2, const char* det3)
{
//
-// Ads a user scoring option to th list
+// Adds a user scoring option to the list
//
- TFlukaScoringOption* opt = new TFlukaScoringOption(option, "User Scoring", npar, what);
+ TFlukaScoringOption* opt = new TFlukaScoringOption(option, sdum, npr, outfile, what, det1, det2, det3);
fUserScore->Add(opt);
}
-
//______________________________________________________________________________
Double_t TFluka::Xsec(char*, Double_t, Int_t, Int_t)
{
- printf("WARNING: Xsec not yet implemented !\n"); return -1.;
+ Warning("Xsec", "Not yet implemented.!\n"); return -1.;
}
//
// Physics initialisation with preparation of FLUKA input cards
//
- printf("=>InitPhysics\n");
- Int_t j, k;
- Double_t fCut;
-
- FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
-
- Double_t zero = 0.0;
- Double_t one = 1.0;
- Double_t two = 2.0;
- Double_t three = 3.0;
-
- Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
- if (fVerbosityLevel >= 3) printf(" last FLUKA material is %g\n", fLastMaterial);
-
- // Prepare Cerenkov
- TObjArray *matList = GetFlukaMaterials();
- Int_t nmaterial = matList->GetEntriesFast();
- fMaterials = new Int_t[nmaterial+3];
-
-// construct file names
-
- TString sFlukaVmcCoreInp = getenv("ALICE_ROOT");
- sFlukaVmcCoreInp +="/TFluka/input/";
- TString sFlukaVmcTmp = "flukaMat.inp";
- TString sFlukaVmcInp = GetInputFileName();
- sFlukaVmcCoreInp += GetCoreInputFileName();
-
-// open files
-
- if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
- printf("\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
- exit(1);
- }
- if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
- printf("\nCannot open file %s\n",sFlukaVmcTmp.Data());
- exit(1);
- }
- if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
- printf("\nCannot open file %s\n",sFlukaVmcInp.Data());
- exit(1);
- }
-
-// copy core input file
- Char_t sLine[255];
- Float_t fEventsPerRun;
-
- while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
- if (strncmp(sLine,"GEOEND",6) != 0)
- fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
- else {
- fprintf(pFlukaVmcInp,"GEOEND\n"); // add GEOEND card
- goto flukamat;
- }
- } // end of while until GEOEND card
-
-
- flukamat:
- while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
- fprintf(pFlukaVmcInp,"%s\n",sLine);
- }
-
- while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
- if (strncmp(sLine,"START",5) != 0)
- fprintf(pFlukaVmcInp,"%s\n",sLine);
- else {
- sscanf(sLine+10,"%10f",&fEventsPerRun);
- goto fin;
- }
- } //end of while until START card
-
-fin:
-// in G3 the process control values meaning can be different for
-// different processes, but for most of them is:
-// 0 process is not activated
-// 1 process is activated WITH generation of secondaries
-// 2 process is activated WITHOUT generation of secondaries
-// if process does not generate secondaries => 1 same as 2
-//
-// Exceptions:
-// MULS: also 3
-// LOSS: also 3, 4
-// RAYL: only 0,1
-// HADR: may be > 2
-//
-
-// Loop over number of SetProcess calls
- fprintf(pFlukaVmcInp,"*----------------------------------------------------------------------------- \n");
- fprintf(pFlukaVmcInp,"*----- The following data are generated from SetProcess and SetCut calls ----- \n");
- fprintf(pFlukaVmcInp,"*----------------------------------------------------------------------------- \n");
-
-// Outer loop over processes
- TIter next(fProcesses);
- TFlukaConfigOption *proc;
-// Inner loop over processes
- TIter nextp(fProcesses);
- TFlukaConfigOption *procp;
-// Loop over cuts
- TIter nextc(fCuts);
- TFlukaConfigOption *cut = 0x0;
-
- while((proc = (TFlukaConfigOption*)next())) {
- Float_t matMin = three;
- Float_t matMax = fLastMaterial;
- Bool_t global = kTRUE;
- if (proc->Medium() != -1) {
- Int_t mat;
- if ((mat = proc->Medium()) >= GetFlukaMaterials()->GetEntries()) continue;
- matMin = Float_t(mat);
- matMax = matMin;
- global = kFALSE;
-
- fprintf(pFlukaVmcInp,"*\n*Material specific process setting for #%8d \n", mat);
- printf("Process for %d \n", mat);
- TGeoMaterial* material = (TGeoMaterial*) (GetFlukaMaterials())->At(GetMaterialIndex(mat));
- printf("Process for %d %s\n", mat, material->GetName());
-
- }
-
- // annihilation
- // G3 default value: 1
- // G4 processes: G4eplusAnnihilation/G4IeplusAnnihilation
- // Particles: e+
- // Physics: EM
- // flag = 0 no annihilation
- // flag = 1 annihilation, decays processed
- // flag = 2 annihilation, no decay product stored
- // gMC ->SetProcess("ANNI",1); // EMFCUT -1. 0. 0. 3. lastmat 0. ANNH-THR
- if (strncmp(proc->GetName(),"ANNI",4) == 0) {
- if (proc->Flag() == 1 || proc->Flag() == 2) {
- fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for e+ annihilation - resets to default=0.\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('ANNI',1) or SetProcess('ANNI',2)\n");
- // -one = kinetic energy threshold (GeV) for e+ annihilation (resets to default=0)
- // zero = not used
- // zero = not used
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- // one = step length in assigning indices
- // "ANNH-THR";
- fprintf(pFlukaVmcInp,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",-one,zero,zero,matMin,matMax,one);
- }
- else if (proc->Flag() == 0) {
- fprintf(pFlukaVmcInp,"*\n*No annihilation - no FLUKA card generated\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('ANNI',0)\n");
- }
- else {
- fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('ANNI',?) call.\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- }
+// Construct file names
+ FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
+ TString sFlukaVmcCoreInp = getenv("ALICE_ROOT");
+ sFlukaVmcCoreInp +="/TFluka/input/";
+ TString sFlukaVmcTmp = "flukaMat.inp";
+ TString sFlukaVmcInp = GetInputFileName();
+ sFlukaVmcCoreInp += GetCoreInputFileName();
- // bremsstrahlung and pair production are both activated
- // G3 default value: 1
- // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung,
- // G4MuBremsstrahlung/G4IMuBremsstrahlung,
- // G4LowEnergyBremstrahlung
- // Particles: e-/e+; mu+/mu-
- // Physics: EM
- // flag = 0 no bremsstrahlung
- // flag = 1 bremsstrahlung, photon processed
- // flag = 2 bremsstrahlung, no photon stored
- // gMC ->SetProcess("BREM",1); // PAIRBREM 2. 0. 0. 3. lastmat
- // EMFCUT -1. 0. 0. 3. lastmat 0. ELPO-THR
- // G3 default value: 1
- // G4 processes: G4GammaConversion,
- // G4MuPairProduction/G4IMuPairProduction
- // G4LowEnergyGammaConversion
- // Particles: gamma, mu
- // Physics: EM
- // flag = 0 no delta rays
- // flag = 1 delta rays, secondaries processed
- // flag = 2 delta rays, no secondaries stored
- // gMC ->SetProcess("PAIR",1); // PAIRBREM 1. 0. 0. 3. lastmat
- // EMFCUT 0. 0. -1. 3. lastmat 0. PHOT-THR
- else if ((strncmp(proc->GetName(),"PAIR",4) == 0) && (proc->Flag() == 1 || proc->Flag() == 2)) {
-
- nextp.Reset();
-
- while ((procp = (TFlukaConfigOption*)nextp())) {
- if ((strncmp(procp->GetName(),"BREM",4) == 0) &&
- (proc->Flag() == 1 || procp->Flag() == 2) &&
- (procp->Medium() == proc->Medium())) {
- fprintf(pFlukaVmcInp,"*\n*Bremsstrahlung and pair production by muons and charged hadrons both activated\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1) and SetProcess('PAIR',1)\n");
- fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.\n");
- fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.\n");
- // three = bremsstrahlung and pair production by muons and charged hadrons both are activated
- fprintf(pFlukaVmcInp,"PAIRBREM %10.1f",three);
- // direct pair production by muons
- // G4 particles: "e-", "e+"
- // G3 default value: 0.01 GeV
- //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons
- fCut = 0.0;
- nextc.Reset();
- while ((cut = (TFlukaConfigOption*)nextc())) {
- if (strncmp(cut->GetName(), "PPCUTM", 6) == 0 &&
- (cut->Medium() == proc->Medium())) fCut = cut->Cut();
- }
- fprintf(pFlukaVmcInp,"%10.4g",fCut);
- // fCut; = e+, e- kinetic energy threshold (in GeV) for explicit pair production.
- // muon and hadron bremsstrahlung
- // G4 particles: "gamma"
- // G3 default value: CUTGAM=0.001 GeV
- //gMC ->SetCut("BCUTM",cut); // cut for muon and hadron bremsstrahlung
- fCut = 0.0;
- nextc.Reset();
- while ((cut = (TFlukaConfigOption*)nextc())) {
- if (strncmp(cut->GetName(), "BCUTM", 5) == 0 &&
- (cut->Medium() == proc->Medium())) fCut = cut->Cut();
- }
- fprintf(pFlukaVmcInp,"%10.4g%10.1f%10.1f\n",fCut,matMin,matMax);
- // fCut = photon energy threshold (GeV) for explicit bremsstrahlung production
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
-
- // for e+ and e-
- fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1);\n");
- fCut = -1.0;
- nextc.Reset();
- while ((cut = (TFlukaConfigOption*)nextc())) {
- if (strncmp(cut->GetName(), "BCUTE", 5) == 0 &&
- (cut->Medium() == proc->Medium())) fCut = cut->Cut();
- }
- //fCut = kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
- // zero = not used
- // zero = not used
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- // one = step length in assigning indices
- // "ELPO-THR";
- fprintf(pFlukaVmcInp,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",fCut,zero,zero,matMin,matMax,one);
-
- // for e+ and e-
- fprintf(pFlukaVmcInp,"*\n*Pair production by electrons is activated\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PAIR',1);\n");
- fCut = -1.0;
- nextc.Reset();
- while ((cut = (TFlukaConfigOption*)nextc())) {
- if (strncmp(cut->GetName(), "CUTGAM", 6) == 0 &&
- (cut->Medium() == proc->Medium())) fCut = cut->Cut();
- }
- // fCut = energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored)
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- // one = step length in assigning indices
- fprintf(pFlukaVmcInp,"EMFCUT %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,fCut,matMin,matMax,one);
- goto BOTH;
- } // end of if for BREM
- } // end of loop for BREM
-
- // only pair production by muons and charged hadrons is activated
- fprintf(pFlukaVmcInp,"*\n*Pair production by muons and charged hadrons is activated\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n");
- fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('PPCUTM',cut) or set to 0.\n");
- // direct pair production by muons
- // G4 particles: "e-", "e+"
- // G3 default value: 0.01 GeV
- //gMC ->SetCut("PPCUTM",cut); // total energy cut for direct pair prod. by muons
- // one = pair production by muons and charged hadrons is activated
- // zero = e+, e- kinetic energy threshold (in GeV) for explicit pair production.
- // zero = no explicit bremsstrahlung production is simulated
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- fprintf(pFlukaVmcInp,"PAIRBREM %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax);
-
- // for e+ and e-
- fprintf(pFlukaVmcInp,"*\n*Pair production by electrons is activated\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PAIR',1) or SetProcess('PAIR',2)\n");
- fCut = -1.0;
- nextc.Reset();
- while ((cut = (TFlukaConfigOption*)nextc())) {
- if (strncmp(cut->GetName(), "CUTGAM", 6) == 0 &&
- (cut->Medium() == proc->Medium())) fCut = cut->Cut();
- }
- // zero = energy threshold (GeV) for Compton scattering (= 0.0 : ignored)
- // zero = energy threshold (GeV) for Photoelectric (= 0.0 : ignored)
- // fCut = energy threshold (GeV) for gamma pair production (< 0.0 : resets to default, = 0.0 : ignored)
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- // one = step length in assigning indices
- fprintf(pFlukaVmcInp,"EMFCUT %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",zero,zero,fCut,matMin,matMax,one);
-
- BOTH:
- k = 0;
- } // end of if for PAIR
-
-
-
- // bremsstrahlung
- // G3 default value: 1
- // G4 processes: G4eBremsstrahlung/G4IeBremsstrahlung,
- // G4MuBremsstrahlung/G4IMuBremsstrahlung,
- // G4LowEnergyBremstrahlung
- // Particles: e-/e+; mu+/mu-
- // Physics: EM
- // flag = 0 no bremsstrahlung
- // flag = 1 bremsstrahlung, photon processed
- // flag = 2 bremsstrahlung, no photon stored
- // gMC ->SetProcess("BREM",1); // PAIRBREM 2. 0. 0. 3. lastmat
- // EMFCUT -1. 0. 0. 3. lastmat 0. ELPO-THR
- else if (strncmp(proc->GetName(),"BREM",4) == 0) {
- nextp.Reset();
- while((procp = (TFlukaConfigOption*)nextp())) {
- if ((strncmp(procp->GetName(),"PAIR",4) == 0) &&
- procp->Flag() == 1 &&
- (procp->Medium() == proc->Medium())) goto NOBREM;
- }
- if (proc->Flag() == 1 || proc->Flag() == 2) {
- fprintf(pFlukaVmcInp,"*\n*Bremsstrahlung by muons and charged hadrons is activated\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1) or SetProcess('BREM',2)\n");
- fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('BCUTM',cut) or set to 0.\n");
- // two = bremsstrahlung by muons and charged hadrons is activated
- // zero = no meaning
- // muon and hadron bremsstrahlung
- // G4 particles: "gamma"
- // G3 default value: CUTGAM=0.001 GeV
- //gMC ->SetCut("BCUTM",cut); // cut for muon and hadron bremsstrahlung
- fCut = 0.0;
- nextc.Reset();
- while ((cut = (TFlukaConfigOption*)nextc())) {
- if (strncmp(cut->GetName(), "BCUTM", 5) == 0 &&
- (cut->Medium() == proc->Medium())) fCut = cut->Cut();
- }
- // fCut = photon energy threshold (GeV) for explicit bremsstrahlung production
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- fprintf(pFlukaVmcInp,"PAIRBREM %10.1f%10.1f%10.4g%10.1f%10.1f\n",two,zero,fCut,matMin,matMax);
-
- // for e+ and e-
- fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for e+/e- bremsstrahlung - resets to default=0.\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',1);");
- // - one = kinetic energy threshold (GeV) for e+/e- bremsstrahlung (resets to default=0)
- // zero = not used
- // zero = not used
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- // one = step length in assigning indices
- //"ELPO-THR";
- fprintf(pFlukaVmcInp,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",-one,zero,zero,matMin,matMax,one);
- }
- else if (proc->Flag() == 0) {
- fprintf(pFlukaVmcInp,"*\n*No bremsstrahlung - no FLUKA card generated\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('BREM',0)\n");
- }
- else {
- fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('BREM',?) call.\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- NOBREM:
- j = 0;
- } // end of else if (strncmp(proc->GetName(),"BREM",4) == 0)
-
- // Cerenkov photon generation
- // G3 default value: 0
- // G4 process: G4Cerenkov
- //
- // Particles: charged
- // Physics: Optical
- // flag = 0 no Cerenkov photon generation
- // flag = 1 Cerenkov photon generation
- // flag = 2 Cerenkov photon generation with primary stopped at each step
- //xx gMC ->SetProcess("CKOV",1); // ??? Cerenkov photon generation
-
- else if (strncmp(proc->GetName(),"CKOV",4) == 0) {
- if ((proc->Flag() == 1 || proc->Flag() == 2) && global) {
- // Write comments
- fprintf(pFlukaVmcInp, "* \n");
- fprintf(pFlukaVmcInp, "*Cerenkov photon generation\n");
- fprintf(pFlukaVmcInp, "*Generated from call: SetProcess('CKOV',1) or SetProcess('CKOV',2)\n");
- // Loop over media
- for (Int_t im = 0; im < nmaterial; im++)
- {
- TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
- Int_t idmat = material->GetIndex();
-
- if (!global && idmat != proc->Medium()) continue;
-
- fMaterials[idmat] = im;
- // Skip media with no Cerenkov properties
- TFlukaCerenkov* cerenkovProp;
- if (!(cerenkovProp = dynamic_cast<TFlukaCerenkov*>(material->GetCerenkovProperties()))) continue;
- //
- // This medium has Cerenkov properties
- //
- //
- // Write OPT-PROD card for each medium
- Float_t emin = cerenkovProp->GetMinimumEnergy();
- Float_t emax = cerenkovProp->GetMaximumEnergy();
- fprintf(pFlukaVmcInp, "OPT-PROD %10.4g%10.4g%10.4g%10.4g%10.4g%10.4gCERENKOV\n", emin, emax, 0.,
- Float_t(idmat), Float_t(idmat), 0.);
- //
- // Write OPT-PROP card for each medium
- // Forcing FLUKA to call user routines (queffc.cxx, rflctv.cxx, rfrndx.cxx)
- //
- fprintf(pFlukaVmcInp, "OPT-PROP %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fWV-LIMIT\n",
- cerenkovProp->GetMinimumWavelength(),
- cerenkovProp->GetMaximumWavelength(),
- cerenkovProp->GetMaximumWavelength(),
- Float_t(idmat), Float_t(idmat), 0.0);
-
- if (cerenkovProp->IsMetal()) {
- fprintf(pFlukaVmcInp, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fMETAL\n",
- -100., -100., -100.,
- Float_t(idmat), Float_t(idmat), 0.0);
- } else {
- fprintf(pFlukaVmcInp, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",
- -100., -100., -100.,
- Float_t(idmat), Float_t(idmat), 0.0);
- }
-
-
- for (Int_t j = 0; j < 3; j++) {
- fprintf(pFlukaVmcInp, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f&\n",
- -100., -100., -100.,
- Float_t(idmat), Float_t(idmat), 0.0);
- }
- // Photon detection efficiency user defined
-
- if (cerenkovProp->IsSensitive())
- fprintf(pFlukaVmcInp, "OPT-PROP %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fSENSITIV\n",
- -100., -100., -100.,
- Float_t(idmat), Float_t(idmat), 0.0);
-
- } // materials
- } else if (proc->Flag() == 0) {
- fprintf(pFlukaVmcInp,"*\n*No Cerenkov photon generation\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('CKOV',0)\n");
- // zero = not used
- // zero = not used
- // zero = not used
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- // one = step length in assigning indices
- //"CERE-OFF";
- fprintf(pFlukaVmcInp,"OPT-PROD %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",zero,zero,zero,matMin,matMax,one);
- }
- else {
- fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('CKOV',?) call.\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- } // end of else if (strncmp(proc->GetName(),"CKOV",4) == 0)
-
- // Compton scattering
- // G3 default value: 1
- // G4 processes: G4ComptonScattering,
- // G4LowEnergyCompton,
- // G4PolarizedComptonScattering
- // Particles: gamma
- // Physics: EM
- // flag = 0 no Compton scattering
- // flag = 1 Compton scattering, electron processed
- // flag = 2 Compton scattering, no electron stored
- // gMC ->SetProcess("COMP",1); // EMFCUT -1. 0. 0. 3. lastmat 0. PHOT-THR
- else if (strncmp(proc->GetName(),"COMP",4) == 0) {
- if (proc->Flag() == 1 || proc->Flag() == 2) {
- fprintf(pFlukaVmcInp,"*\n*Energy threshold (GeV) for Compton scattering - resets to default=0.\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('COMP',1);\n");
- // - one = energy threshold (GeV) for Compton scattering - resets to default=0.
- // zero = not used
- // zero = not used
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- // one = step length in assigning indices
- //"PHOT-THR";
- fprintf(pFlukaVmcInp,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",-one,zero,zero,matMin,matMax,one);
- }
- else if (proc->Flag() == 0) {
- fprintf(pFlukaVmcInp,"*\n*No Compton scattering - no FLUKA card generated\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('COMP',0)\n");
- }
- else {
- fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('COMP',?) call.\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- } // end of else if (strncmp(proc->GetName(),"COMP",4) == 0)
-
- // decay
- // G3 default value: 1
- // G4 process: G4Decay
- //
- // Particles: all which decay is applicable for
- // Physics: General
- // flag = 0 no decays
- // flag = 1 decays, secondaries processed
- // flag = 2 decays, no secondaries stored
- //gMC ->SetProcess("DCAY",0); // not available
- else if ((strncmp(proc->GetName(),"DCAY",4) == 0) && proc->Flag() == 0)
- cout << "SetProcess for flag =" << proc->GetName() << " value=" << proc->Flag() << " not avaliable!" << endl;
- else if ((strncmp(proc->GetName(),"DCAY",4) == 0) && proc->Flag() == 1) {
- // Nothing to do decays are switched on by default
- }
-
-
- // delta-ray
- // G3 default value: 2
- // !! G4 treats delta rays in different way
- // G4 processes: G4eIonisation/G4IeIonization,
- // G4MuIonisation/G4IMuIonization,
- // G4hIonisation/G4IhIonisation
- // Particles: charged
- // Physics: EM
- // flag = 0 no energy loss
- // flag = 1 restricted energy loss fluctuations
- // flag = 2 complete energy loss fluctuations
- // flag = 3 same as 1
- // flag = 4 no energy loss fluctuations
- // gMC ->SetProcess("DRAY",0); // DELTARAY 1.E+6 0. 0. 3. lastmat 0.
- else if (strncmp(proc->GetName(),"DRAY",4) == 0) {
- if (proc->Flag() == 0 || proc->Flag() == 4) {
- fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('DRAY',0) or SetProcess('DRAY',4)\n");
- fprintf(pFlukaVmcInp,"*No delta ray production by muons - threshold set artificially high\n");
- Double_t emin = 1.0e+6; // kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
- // zero = ignored
- // zero = ignored
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- // one = step length in assigning indices
- fprintf(pFlukaVmcInp,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n",emin,zero,zero,matMin,matMax,one);
- }
- else if (proc->Flag() == 1 || proc->Flag() == 2 || proc->Flag() == 3) {
- fprintf(pFlukaVmcInp,"*\n*Kinetic energy threshold (GeV) for delta ray production\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('DRAY',flag), flag=1,2,3\n");
- fprintf(pFlukaVmcInp,"*Delta ray production by muons switched on\n");
- fprintf(pFlukaVmcInp,"*Energy threshold set by call SetCut('DCUTM',cut) or set to 1.0e+6.\n");
- fCut = 1.0e+6;
- nextc.Reset();
- while ((cut = (TFlukaConfigOption*)nextc())) {
- if (strncmp(cut->GetName(), "DCUTM", 5) == 0 &&
- cut->Medium() == proc->Medium()) fCut = cut->Cut();
- }
- // fCut = kinetic energy threshold (GeV) for delta ray production (discrete energy transfer)
- // zero = ignored
- // zero = ignored
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- // one = step length in assigning indices
- fprintf(pFlukaVmcInp,"DELTARAY %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n",fCut,zero,zero,matMin,matMax,one);
- }
- else {
- fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('DRAY',?) call.\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- } // end of else if (strncmp(proc->GetName(),"DRAY",4) == 0)
-
- // hadronic process
- // G3 default value: 1
- // G4 processes: all defined by TG4PhysicsConstructorHadron
- //
- // Particles: hadrons
- // Physics: Hadron
- // flag = 0 no multiple scattering
- // flag = 1 hadronic interactions, secondaries processed
- // flag = 2 hadronic interactions, no secondaries stored
- // gMC ->SetProcess("HADR",1); // ??? hadronic process
- //Select pure GEANH (HADR 1) or GEANH/NUCRIN (HADR 3) ?????
- else if (strncmp(proc->GetName(),"HADR",4) == 0) {
- if (proc->Flag() == 1 || proc->Flag() == 2) {
- fprintf(pFlukaVmcInp,"*\n*Hadronic interaction is ON by default in FLUKA\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- else if (proc->Flag() == 0) {
- fprintf(pFlukaVmcInp,"*\n*Hadronic interaction is set OFF\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('HADR',0);\n");
- fprintf(pFlukaVmcInp,"*Switching off hadronic interactions not foreseen in FLUKA\n");
- fprintf(pFlukaVmcInp,"THRESHOL %10.1f%10.1f%10.1f%10.1e%10.1f\n",zero, zero, zero, 1.e10, zero);
- }
- else {
- fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('HADR',?) call.\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- } // end of else if (strncmp(proc->GetName(),"HADR",4) == 0)
-
-
- // energy loss
- // G3 default value: 2
- // G4 processes: G4eIonisation/G4IeIonization,
- // G4MuIonisation/G4IMuIonization,
- // G4hIonisation/G4IhIonisation
- //
- // Particles: charged
- // Physics: EM
- // flag=0 no energy loss
- // flag=1 restricted energy loss fluctuations
- // flag=2 complete energy loss fluctuations
- // flag=3 same as 1
- // flag=4 no energy loss fluctuations
- // If the value ILOSS is changed, then (in G3) cross-sections and energy
- // loss tables must be recomputed via the command 'PHYSI'
- // gMC ->SetProcess("LOSS",2); // ??? IONFLUCT ? energy loss
- else if (strncmp(proc->GetName(),"LOSS",4) == 0) {
- if (proc->Flag() == 2) { // complete energy loss fluctuations
- fprintf(pFlukaVmcInp,"*\n*Complete energy loss fluctuations do not exist in FLUKA\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('LOSS',2);\n");
- fprintf(pFlukaVmcInp,"*flag=2=complete energy loss fluctuations\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- else if (proc->Flag() == 1 || proc->Flag() == 3) { // restricted energy loss fluctuations
- fprintf(pFlukaVmcInp,"*\n*Restricted energy loss fluctuations\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('LOSS',1) or SetProcess('LOSS',3)\n");
- // one = restricted energy loss fluctuations (for hadrons and muons) switched on
- // one = restricted energy loss fluctuations (for e+ and e-) switched on
- // one = minimal accuracy
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // upper bound of the material indices in which the respective thresholds apply
- fprintf(pFlukaVmcInp,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,one,one,matMin,matMax);
- }
- else if (proc->Flag() == 4) { // no energy loss fluctuations
- fprintf(pFlukaVmcInp,"*\n*No energy loss fluctuations\n");
- fprintf(pFlukaVmcInp,"*\n*Generated from call: SetProcess('LOSS',4)\n");
- // - one = restricted energy loss fluctuations (for hadrons and muons) switched off
- // - one = restricted energy loss fluctuations (for e+ and e-) switched off
- // one = minimal accuracy
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- fprintf(pFlukaVmcInp,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,-one,one,matMin,matMax);
- }
- else {
- fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('LOSS',?) call.\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- } // end of else if (strncmp(proc->GetName(),"LOSS",4) == 0)
-
-
- // multiple scattering
- // G3 default value: 1
- // G4 process: G4MultipleScattering/G4IMultipleScattering
- //
- // Particles: charged
- // Physics: EM
- // flag = 0 no multiple scattering
- // flag = 1 Moliere or Coulomb scattering
- // flag = 2 Moliere or Coulomb scattering
- // flag = 3 Gaussian scattering
- // gMC ->SetProcess("MULS",1); // MULSOPT multiple scattering
- else if (strncmp(proc->GetName(),"MULS",4) == 0) {
- if (proc->Flag() == 1 || proc->Flag() == 2 || proc->Flag() == 3) {
- fprintf(pFlukaVmcInp,"*\n*Multiple scattering is ON by default for e+e- and for hadrons/muons\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- else if (proc->Flag() == 0) {
- fprintf(pFlukaVmcInp,"*\n*Multiple scattering is set OFF\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('MULS',0);\n");
- // zero = ignored
- // three = multiple scattering for hadrons and muons is completely suppressed
- // three = multiple scattering for e+ and e- is completely suppressed
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- fprintf(pFlukaVmcInp,"MULSOPT %10.1f%10.1f%10.1f%10.1f%10.1f\n",zero,three,three,matMin,matMax);
- }
- else {
- fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('MULS',?) call.\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- } // end of else if (strncmp(proc->GetName(),"MULS",4) == 0)
-
-
- // muon nuclear interaction
- // G3 default value: 0
- // G4 processes: G4MuNuclearInteraction,
- // G4MuonMinusCaptureAtRest
- //
- // Particles: mu
- // Physics: Not set
- // flag = 0 no muon-nuclear interaction
- // flag = 1 nuclear interaction, secondaries processed
- // flag = 2 nuclear interaction, secondaries not processed
- // gMC ->SetProcess("MUNU",1); // MUPHOTON 1. 0. 0. 3. lastmat
- else if (strncmp(proc->GetName(),"MUNU",4) == 0) {
- if (proc->Flag() == 1) {
- fprintf(pFlukaVmcInp,"*\n*Muon nuclear interactions with production of secondary hadrons\n");
- fprintf(pFlukaVmcInp,"*\n*Generated from call: SetProcess('MUNU',1);\n");
- // one = full simulation of muon nuclear interactions and production of secondary hadrons
- // zero = ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
- // zero = fraction of rho-like interactions ( must be < 1) - Default = 0.75.
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- fprintf(pFlukaVmcInp,"MUPHOTON %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax);
- }
- else if (proc->Flag() == 2) {
- fprintf(pFlukaVmcInp,"*\n*Muon nuclear interactions without production of secondary hadrons\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('MUNU',2);\n");
- // two = full simulation of muon nuclear interactions and production of secondary hadrons
- // zero = ratio of longitudinal to transverse virtual photon cross-section - Default = 0.25.
- // zero = fraction of rho-like interactions ( must be < 1) - Default = 0.75.
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- fprintf(pFlukaVmcInp,"MUPHOTON %10.1f%10.1f%10.1f%10.1f%10.1f\n",two,zero,zero,matMin,matMax);
- }
- else if (proc->Flag() == 0) {
- fprintf(pFlukaVmcInp,"*\n*No muon nuclear interaction - no FLUKA card generated\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('MUNU',0)\n");
- }
- else {
- fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('MUNU',?) call.\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- } // end of else if (strncmp(proc->GetName(),"MUNU",4) == 0)
-
-
- // photofission
- // G3 default value: 0
- // G4 process: ??
- //
- // Particles: gamma
- // Physics: ??
- // gMC ->SetProcess("PFIS",0); // PHOTONUC -1. 0. 0. 3. lastmat 0.
- // flag = 0 no photon fission
- // flag = 1 photon fission, secondaries processed
- // flag = 2 photon fission, no secondaries stored
- else if (strncmp(proc->GetName(),"PFIS",4) == 0) {
- if (proc->Flag() == 0) {
- fprintf(pFlukaVmcInp,"*\n*No photonuclear interactions\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PFIS',0);\n");
- // - one = no photonuclear interactions
- // zero = not used
- // zero = not used
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- fprintf(pFlukaVmcInp,"PHOTONUC %10.1f%10.1f%10.1f%10.1f%10.1f\n",-one,zero,zero,matMin,matMax);
- }
- else if (proc->Flag() == 1) {
- fprintf(pFlukaVmcInp,"*\n*Photon nuclear interactions are activated at all energies\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PFIS',1);\n");
- // one = photonuclear interactions are activated at all energies
- // zero = not used
- // zero = not used
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- fprintf(pFlukaVmcInp,"PHOTONUC %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,zero,zero,matMin,matMax);
- }
- else if (proc->Flag() == 0) {
- fprintf(pFlukaVmcInp,"*\n*No photofission - no FLUKA card generated\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PFIS',0)\n");
- }
- else {
- fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('PFIS',?) call.\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- }
-
-
- // photo electric effect
- // G3 default value: 1
- // G4 processes: G4PhotoElectricEffect
- // G4LowEnergyPhotoElectric
- // Particles: gamma
- // Physics: EM
- // flag = 0 no photo electric effect
- // flag = 1 photo electric effect, electron processed
- // flag = 2 photo electric effect, no electron stored
- // gMC ->SetProcess("PHOT",1); // EMFCUT 0. -1. 0. 3. lastmat 0. PHOT-THR
- else if (strncmp(proc->GetName(),"PHOT",4) == 0) {
- if (proc->Flag() == 1 || proc->Flag() == 2) {
- fprintf(pFlukaVmcInp,"*\n*Photo electric effect is activated\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PHOT',1);\n");
- // zero = ignored
- // - one = resets to default=0.
- // zero = ignored
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- // one = step length in assigning indices
- //"PHOT-THR";
- fprintf(pFlukaVmcInp,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",zero,-one,zero,matMin,matMax,one);
- }
- else if (proc->Flag() == 0) {
- fprintf(pFlukaVmcInp,"*\n*No photo electric effect - no FLUKA card generated\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('PHOT',0)\n");
- }
- else {
- fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('PHOT',?) call.\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- } // else if (strncmp(proc->GetName(),"PHOT",4) == 0)
-
-
- // Rayleigh scattering
- // G3 default value: 0
- // G4 process: G4OpRayleigh
- //
- // Particles: optical photon
- // Physics: Optical
- // flag = 0 Rayleigh scattering off
- // flag = 1 Rayleigh scattering on
- //xx gMC ->SetProcess("RAYL",1);
- else if (strncmp(proc->GetName(),"RAYL",4) == 0) {
- if (proc->Flag() == 1) {
- fprintf(pFlukaVmcInp,"*\n*Rayleigh scattering is ON by default in FLUKA\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- else if (proc->Flag() == 0) {
- fprintf(pFlukaVmcInp,"*\n*Rayleigh scattering is set OFF\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('RAYL',0);\n");
- // - one = no Rayleigh scattering and no binding corrections for Compton
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- fprintf(pFlukaVmcInp,"EMFRAY %10.1f%10.1f%10.1f%10.1f\n",-one,three,matMin,matMax);
- }
- else {
- fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('RAYL',?) call.\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- } // end of else if (strncmp(proc->GetName(),"RAYL",4) == 0)
-
-
- // synchrotron radiation in magnetic field
- // G3 default value: 0
- // G4 process: G4SynchrotronRadiation
- //
- // Particles: ??
- // Physics: Not set
- // flag = 0 no synchrotron radiation
- // flag = 1 synchrotron radiation
- //xx gMC ->SetProcess("SYNC",1); // synchrotron radiation generation
- else if (strncmp(proc->GetName(),"SYNC",4) == 0) {
- fprintf(pFlukaVmcInp,"*\n*Synchrotron radiation generation is NOT implemented in FLUKA\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
-
-
- // Automatic calculation of tracking medium parameters
- // flag = 0 no automatic calculation
- // flag = 1 automatic calculation
- //xx gMC ->SetProcess("AUTO",1); // ??? automatic computation of the tracking medium parameters
- else if (strncmp(proc->GetName(),"AUTO",4) == 0) {
- fprintf(pFlukaVmcInp,"*\n*Automatic calculation of tracking medium parameters is always ON in FLUKA\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
-
-
- // To control energy loss fluctuation model
- // flag = 0 Urban model
- // flag = 1 PAI model
- // flag = 2 PAI+ASHO model (not active at the moment)
- //xx gMC ->SetProcess("STRA",1); // ??? energy fluctuation model
- else if (strncmp(proc->GetName(),"STRA",4) == 0) {
- if (proc->Flag() == 0 || proc->Flag() == 2 || proc->Flag() == 3) {
- fprintf(pFlukaVmcInp,"*\n*Ionization energy losses calculation is activated\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetProcess('STRA',n);, n=0,1,2\n");
- // one = restricted energy loss fluctuations (for hadrons and muons) switched on
- // one = restricted energy loss fluctuations (for e+ and e-) switched on
- // one = minimal accuracy
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- fprintf(pFlukaVmcInp,"IONFLUCT %10.1f%10.1f%10.1f%10.1f%10.1f\n",one,one,one,matMin,matMax);
- }
- else {
- fprintf(pFlukaVmcInp,"*\n*Illegal flag value in SetProcess('STRA',?) call.\n");
- fprintf(pFlukaVmcInp,"*No FLUKA card generated\n");
- }
- } // else if (strncmp(proc->GetName(),"STRA",4) == 0)
-
-
-
-
- else { // processes not yet treated
-
- // light photon absorption (Cerenkov photons)
- // it is turned on when Cerenkov process is turned on
- // G3 default value: 0
- // G4 process: G4OpAbsorption, G4OpBoundaryProcess
- //
- // Particles: optical photon
- // Physics: Optical
- // flag = 0 no absorption of Cerenkov photons
- // flag = 1 absorption of Cerenkov photons
- // gMC ->SetProcess("LABS",2); // ??? Cerenkov light absorption
-
-
-
- cout << "SetProcess for flag=" << proc->GetName() << " value=" << proc->Flag() << " not yet implemented!" << endl;
- }
- } //end of loop number of SetProcess calls
+// Open files
+ if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
+ Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
+ exit(1);
+ }
+ if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
+ Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcTmp.Data());
+ exit(1);
+ }
+ if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
+ Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcInp.Data());
+ exit(1);
+ }
-
-// Loop over number of SetCut calls
-
- nextc.Reset();
- while ((cut = (TFlukaConfigOption*)nextc())) {
- Float_t matMin = three;
- Float_t matMax = fLastMaterial;
- Bool_t global = kTRUE;
- if (cut->Medium() != -1) {
- Int_t mat;
- if ((mat = cut->Medium()) >= GetFlukaMaterials()->GetEntries()) continue;
- matMin = Float_t(mat);
- matMax = matMin;
- global = kFALSE;
- TGeoMaterial* material = (TGeoMaterial*) (GetFlukaMaterials())->At(GetMaterialIndex(mat));
- fprintf(pFlukaVmcInp,"*\n*Material specific cut setting for #%8d %s %s %13.3e\n", mat, material->GetName(), cut->GetName(), cut->Cut());
-
- }
-
- // cuts handled in SetProcess calls
- if (strncmp(cut->GetName(),"BCUTM",5) == 0) continue;
- else if (strncmp(cut->GetName(),"BCUTE",5) == 0) continue;
- else if (strncmp(cut->GetName(),"DCUTM",5) == 0) continue;
- else if (strncmp(cut->GetName(),"PPCUTM",6) == 0) continue;
-
- // delta-rays by electrons
- // G4 particles: "e-"
- // G3 default value: 10**4 GeV
- // gMC ->SetCut("DCUTE",cut); // cut for deltarays by electrons
- else if (strncmp(cut->GetName(),"DCUTE",5) == 0) {
- fprintf(pFlukaVmcInp,"*Cut for delta rays by electrons\n");
- // -cut->Cut();
- // zero = ignored
- // zero = ignored
- // matMin = lower bound of the material indices in which the respective thresholds apply
- // matMax = upper bound of the material indices in which the respective thresholds apply
- // loop over materials for EMFCUT FLUKA cards
- for (j=0; j < matMax-matMin+1; j++) {
- Int_t nreg, imat, *reglist;
- Float_t ireg;
- imat = (Int_t) matMin + j;
- reglist = fGeom->GetMaterialList(imat, nreg);
- // loop over regions of a given material
- for (k=0; k<nreg; k++) {
- ireg = reglist[k];
- fprintf(pFlukaVmcInp,"EMFCUT %10.4g%10.1f%10.1f%10.1f%10.1f\n",-cut->Cut(),zero,zero,ireg,ireg);
- }
+// Copy core input file
+ Char_t sLine[255];
+ Float_t fEventsPerRun;
+
+ while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
+ if (strncmp(sLine,"GEOEND",6) != 0)
+ fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
+ else {
+ fprintf(pFlukaVmcInp,"GEOEND\n"); // add GEOEND card
+ goto flukamat;
}
- fprintf(pFlukaVmcInp,"DELTARAY %10.4g%10.3f%10.3f%10.1f%10.1f%10.1f\n",cut->Cut(), 100., 1.03, matMin, matMax, 1.0);
- } // end of if for delta-rays by electrons
+ } // end of while until GEOEND card
- // gammas
- // G4 particles: "gamma"
- // G3 default value: 0.001 GeV
- // gMC ->SetCut("CUTGAM",cut); // cut for gammas
-
- else if (strncmp(cut->GetName(),"CUTGAM",6) == 0 && global) {
- fprintf(pFlukaVmcInp,"*\n*Cut for gamma\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTGAM',cut);\n");
- // -cut->Cut();
- // 7.0 = lower bound of the particle id-numbers to which the cut-off
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f\n",-cut->Cut(),7.0);
- }
- else if (strncmp(cut->GetName(),"CUTGAM",6) == 0 && !global) {
- // loop over materials for EMFCUT FLUKA cards
- for (j=0; j < matMax-matMin+1; j++) {
- Int_t nreg, imat, *reglist;
- Float_t ireg;
- imat = (Int_t) matMin + j;
- reglist = fGeom->GetMaterialList(imat, nreg);
- // loop over regions of a given material
- for (Int_t k = 0; k < nreg; k++) {
- ireg = reglist[k];
- fprintf(pFlukaVmcInp,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", zero, cut->Cut(), zero, ireg, ireg, one);
- }
- }
- } // end of else if for gamma
-
-
- // electrons
- // G4 particles: "e-"
- // ?? positrons
- // G3 default value: 0.001 GeV
- //gMC ->SetCut("CUTELE",cut); // cut for e+,e-
- else if (strncmp(cut->GetName(),"CUTELE",6) == 0 && global) {
- fprintf(pFlukaVmcInp,"*\n*Cut for electrons\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTELE',cut);\n");
- // -cut->Cut();
- // three = lower bound of the particle id-numbers to which the cut-off
- // 4.0 = upper bound of the particle id-numbers to which the cut-off
- // one = step length in assigning numbers
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f%10.1f\n",-cut->Cut(),three,4.0,one);
- }
- else if (strncmp(cut->GetName(),"CUTELE",6) == 0 && !global) {
- // loop over materials for EMFCUT FLUKA cards
- for (j=0; j < matMax-matMin+1; j++) {
- Int_t nreg, imat, *reglist;
- Float_t ireg;
- imat = (Int_t) matMin + j;
- reglist = fGeom->GetMaterialList(imat, nreg);
- // loop over regions of a given material
- for (k=0; k<nreg; k++) {
- ireg = reglist[k];
- fprintf(pFlukaVmcInp,"EMFCUT %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -cut->Cut(), zero, zero, ireg, ireg, one);
- }
+ flukamat:
+ while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
+ fprintf(pFlukaVmcInp,"%s\n",sLine);
+ }
+
+ while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
+ if (strncmp(sLine,"START",5) != 0)
+ fprintf(pFlukaVmcInp,"%s\n",sLine);
+ else {
+ sscanf(sLine+10,"%10f",&fEventsPerRun);
+ goto fin;
}
- } // end of else if for electrons
+ } //end of while until START card
+
+ fin:
- // neutral hadrons
- // G4 particles: of type "baryon", "meson", "nucleus" with zero charge
- // G3 default value: 0.01 GeV
- //gMC ->SetCut("CUTNEU",cut); // cut for neutral hadrons
- else if (strncmp(cut->GetName(),"CUTNEU",6) == 0 && global) {
- fprintf(pFlukaVmcInp,"*\n*Cut for neutral hadrons\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTNEU',cut);\n");
-
- // 8.0 = Neutron
- // 9.0 = Antineutron
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),8.0,9.0);
-
- // 12.0 = Kaon zero long
- // 12.0 = Kaon zero long
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),12.0,12.0);
-
- // 17.0 = Lambda, 18.0 = Antilambda
- // 19.0 = Kaon zero short
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),17.0,19.0);
-
- // 22.0 = Sigma zero, Pion zero, Kaon zero
- // 25.0 = Antikaon zero
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),22.0,25.0);
-
- // 32.0 = Antisigma zero
- // 32.0 = Antisigma zero
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),32.0,32.0);
-
- // 34.0 = Xi zero
- // 35.0 = AntiXi zero
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),34.0,35.0);
-
- // 47.0 = D zero
- // 48.0 = AntiD zero
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),47.0,48.0);
-
- // 53.0 = Xi_c zero
- // 53.0 = Xi_c zero
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),53.0,53.0);
-
- // 55.0 = Xi'_c zero
- // 56.0 = Omega_c zero
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),55.0,56.0);
-
- // 59.0 = AntiXi_c zero
- // 59.0 = AntiXi_c zero
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),59.0,59.0);
-
- // 61.0 = AntiXi'_c zero
- // 62.0 = AntiOmega_c zero
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),61.0,62.0);
- }
-
- // charged hadrons
- // G4 particles: of type "baryon", "meson", "nucleus" with non-zero charge
- // G3 default value: 0.01 GeV
- //gMC ->SetCut("CUTHAD",cut); // cut for charged hadrons
- else if (strncmp(cut->GetName(),"CUTHAD",6) == 0 && global) {
- fprintf(pFlukaVmcInp,"*\n*Cut for charged hadrons\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTHAD',cut);\n");
-
- // 1.0 = Proton
- // 2.0 = Antiproton
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),1.0,2.0);
-
- // 13.0 = Positive Pion, Negative Pion, Positive Kaon
- // 16.0 = Negative Kaon
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),13.0,16.0);
-
- // 20.0 = Negative Sigma
- // 21.0 = Positive Sigma
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),20.0,21.0);
-
- // 31.0 = Antisigma minus
- // 33.0 = Antisigma plus
- // 2.0 = step length
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f%10.1f\n",-cut->Cut(),31.0,33.0,2.0);
-
- // 36.0 = Negative Xi, Positive Xi, Omega minus
- // 39.0 = Antiomega
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),36.0,39.0);
-
- // 45.0 = D plus
- // 46.0 = D minus
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),45.0,46.0);
-
- // 49.0 = D_s plus, D_s minus, Lambda_c plus
- // 52.0 = Xi_c plus
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),49.0,52.0);
-
- // 54.0 = Xi'_c plus
- // 60.0 = AntiXi'_c minus
- // 6.0 = step length
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f%10.1f\n",-cut->Cut(),54.0,60.0,6.0);
-
- // 57.0 = Antilambda_c minus
- // 58.0 = AntiXi_c minus
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),57.0,58.0);
- }
+// Pass information to configuration objects
+
+ Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
+ TFlukaConfigOption::SetStaticInfo(pFlukaVmcInp, 3, fLastMaterial, fGeom);
+
+ TIter next(fUserConfig);
+ TFlukaConfigOption* proc;
+ while((proc = dynamic_cast<TFlukaConfigOption*> (next()))) proc->WriteFlukaInputCards();
+//
+// Process Fluka specific scoring options
+//
+ TFlukaScoringOption::SetStaticInfo(pFlukaVmcInp, fGeom);
+ Float_t loginp = -49.0;
+ Int_t inp = 0;
+ Int_t nscore = fUserScore->GetEntries();
+
+ TFlukaScoringOption *mopo = 0;
+ TFlukaScoringOption *mopi = 0;
- // muons
- // G4 particles: "mu+", "mu-"
- // G3 default value: 0.01 GeV
- //gMC ->SetCut("CUTMUO",cut); // cut for mu+, mu-
- else if (strncmp(cut->GetName(),"CUTMUO",6)== 0 && global) {
- fprintf(pFlukaVmcInp,"*\n*Cut for muons\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetCut('CUTMUO',cut);\n");
- // 10.0 = Muon+
- // 11.0 = Muon-
- fprintf(pFlukaVmcInp,"PART-THR %10.4g%10.1f%10.1f\n",-cut->Cut(),10.0,11.0);
- }
-
- //
- // time of flight cut in seconds
- // G4 particles: all
- // G3 default value: 0.01 GeV
- //gMC ->SetCut("TOFMAX",tofmax); // time of flight cuts in seconds
- else if (strncmp(cut->GetName(),"TOFMAX",6) == 0) {
- fprintf(pFlukaVmcInp,"*\n*Time of flight cuts in seconds\n");
- fprintf(pFlukaVmcInp,"*Generated from call: SetCut('TOFMAX',tofmax);\n");
- // zero = ignored
- // zero = ignored
- // -6.0 = lower bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
- // 64.0 = upper bound of the particle numbers for which the transport time cut-off and/or the start signal is to be applied
- fprintf(pFlukaVmcInp,"TIME-CUT %10.4g%10.1f%10.1f%10.1f%10.1f\n",cut->Cut()*1.e9,zero,zero,-6.0,64.0);
- }
-
- else if (global){
- cout << "SetCut for flag=" << cut->GetName() << " value=" << cut->Cut() << " not yet implemented!" << endl;
- }
- else {
- cout << "SetCut for flag=" << cut->GetName() << " value=" << cut->Cut() << " (material specific) not yet implemented!" << endl;
- }
-
- } //end of loop over SetCut calls
-
+ for (Int_t isc = 0; isc < nscore; isc++)
+ {
+ mopo = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isc));
+ char* fileName = mopo->GetFileName();
+ Int_t size = strlen(fileName);
+ Float_t lun = -1.;
+//
+// Check if new output file has to be opened
+ for (Int_t isci = 0; isci < isc; isci++) {
+
+
+ mopi = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isci));
+ if(strncmp(mopi->GetFileName(), fileName, size)==0) {
+ //
+ // No, the file already exists
+ lun = mopi->GetLun();
+ mopo->SetLun(lun);
+ break;
+ }
+ } // inner loop
+
+ if (lun == -1.) {
+ // Open new output file
+ inp++;
+ mopo->SetLun(loginp + inp);
+ mopo->WriteOpenFlukaFile();
+ }
+ mopo->WriteFlukaInputCards();
+ }
+
+// Add RANDOMIZ card
+ fprintf(pFlukaVmcInp,"RANDOMIZ %10.1f%10.0f\n", 1., Float_t(gRandom->GetSeed()));
// Add START and STOP card
- fprintf(pFlukaVmcInp,"START %10.1f\n",fEventsPerRun);
- fprintf(pFlukaVmcInp,"STOP \n");
+ fprintf(pFlukaVmcInp,"START %10.1f\n",fEventsPerRun);
+ fprintf(pFlukaVmcInp,"STOP \n");
// Close files
-
fclose(pFlukaVmcCoreInp);
fclose(pFlukaVmcFlukaMat);
fclose(pFlukaVmcInp);
-
+
+
+//
+// Initialisation needed for Cerenkov photon production and transport
+ TObjArray *matList = GetFlukaMaterials();
+ Int_t nmaterial = matList->GetEntriesFast();
+ fMaterials = new Int_t[nmaterial+25];
+
+ for (Int_t im = 0; im < nmaterial; im++)
+ {
+ TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
+ Int_t idmat = material->GetIndex();
+ fMaterials[idmat] = im;
+ }
} // end of InitPhysics
void TFluka::SetMaxStep(Double_t step)
{
// Set the maximum step size
- if (step > 1.e4) return;
+// if (step > 1.e4) return;
- Int_t mreg, latt;
- fGeom->GetCurrentRegion(mreg, latt);
+// Int_t mreg=0, latt=0;
+// fGeom->GetCurrentRegion(mreg, latt);
+ Int_t mreg = fGeom->GetCurrentRegion();
STEPSZ.stepmx[mreg - 1] = step;
}
// TRACKR.xtrack = x-position of the last point
// TRACKR.ytrack = y-position of the last point
// TRACKR.ztrack = z-position of the last point
- Int_t caller = GetCaller();
- if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
- position.SetX(GetXsco());
- position.SetY(GetYsco());
- position.SetZ(GetZsco());
- position.SetT(TRACKR.atrack);
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kENDRAW || caller == kUSDRAW ||
+ caller == kBXExiting || caller == kBXEntering ||
+ caller == kUSTCKV) {
+ position.SetX(GetXsco());
+ position.SetY(GetYsco());
+ position.SetZ(GetZsco());
+ position.SetT(TRACKR.atrack);
}
- else if (caller == 4) { // mgdraw
- position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
- position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
- position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
- position.SetT(TRACKR.atrack);
+ else if (caller == kMGDRAW) {
+ Int_t i = -1;
+ if ((i = fPrimaryElectronIndex) > -1) {
+ // Primary Electron Ionisation
+ Double_t x, y, z, t;
+ GetPrimaryElectronPosition(i, x, y, z, t);
+ position.SetX(x);
+ position.SetY(y);
+ position.SetZ(z);
+ position.SetT(t);
+ } else {
+ position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
+ position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
+ position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
+ position.SetT(TRACKR.atrack);
+ }
}
- else if (caller == 5) { // sodraw
- position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
- position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
- position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
- position.SetT(0);
+ else if (caller == kSODRAW) {
+ Int_t ist = FLKSTK.npflka;
+ position.SetX(FLKSTK.xflk[ist]);
+ position.SetY(FLKSTK.yflk[ist]);
+ position.SetZ(FLKSTK.zflk[ist]);
+ position.SetT(FLKSTK.agestk[ist]);
+ } else if (caller == kMGResumedTrack) {
+ position.SetX(TRACKR.spausr[0]);
+ position.SetY(TRACKR.spausr[1]);
+ position.SetZ(TRACKR.spausr[2]);
+ position.SetT(TRACKR.spausr[3]);
}
else
- Warning("TrackPosition","position not available");
+ Warning("TrackPosition","position not available");
}
//______________________________________________________________________________
// TRACKR.xtrack = x-position of the last point
// TRACKR.ytrack = y-position of the last point
// TRACKR.ztrack = z-position of the last point
- Int_t caller = GetCaller();
- if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
- x = GetXsco();
- y = GetYsco();
- z = GetZsco();
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kENDRAW || caller == kUSDRAW ||
+ caller == kBXExiting || caller == kBXEntering ||
+ caller == kUSTCKV) {
+ x = GetXsco();
+ y = GetYsco();
+ z = GetZsco();
+ }
+ else if (caller == kMGDRAW) {
+ Int_t i = -1;
+ if ((i = fPrimaryElectronIndex) > -1) {
+ Double_t t;
+ GetPrimaryElectronPosition(i, x, y, z, t);
+ } else {
+ x = TRACKR.xtrack[TRACKR.ntrack];
+ y = TRACKR.ytrack[TRACKR.ntrack];
+ z = TRACKR.ztrack[TRACKR.ntrack];
+ }
+ }
+ else if (caller == kSODRAW) {
+ Int_t ist = FLKSTK.npflka;
+ x = FLKSTK.xflk[ist];
+ y = FLKSTK.yflk[ist];
+ z = FLKSTK.zflk[ist];
}
- else if (caller == 4 || caller == 5) { // mgdraw, sodraw
- x = TRACKR.xtrack[TRACKR.ntrack];
- y = TRACKR.ytrack[TRACKR.ntrack];
- z = TRACKR.ztrack[TRACKR.ntrack];
+ else if (caller == kMGResumedTrack) {
+ x = TRACKR.spausr[0];
+ y = TRACKR.spausr[1];
+ z = TRACKR.spausr[2];
}
else
- Warning("TrackPosition","position not available");
+ Warning("TrackPosition","position not available");
}
//______________________________________________________________________________
// TRACKR.etrack = total energy of the particle
// TRACKR.jtrack = identity number of the particle
// PAPROP.am[TRACKR.jtrack] = particle mass in gev
- Int_t caller = GetCaller();
- if (caller != 2) { // not eedraw
- if (TRACKR.ptrack >= 0) {
- momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
- momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
- momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
- momentum.SetE(TRACKR.etrack);
- return;
- }
- else {
- Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
- momentum.SetPx(p*TRACKR.cxtrck);
- momentum.SetPy(p*TRACKR.cytrck);
- momentum.SetPz(p*TRACKR.cztrck);
- momentum.SetE(TRACKR.etrack);
+ FlukaCallerCode_t caller = GetCaller();
+ FlukaProcessCode_t icode = GetIcode();
+
+ if (caller != kEEDRAW &&
+ caller != kMGResumedTrack &&
+ caller != kSODRAW &&
+ caller != kUSDRAW &&
+ (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
+ if (TRACKR.ptrack >= 0) {
+ momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
+ momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
+ momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
+ momentum.SetE(TRACKR.etrack);
+ return;
+ }
+ else {
+ Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
+ momentum.SetPx(p*TRACKR.cxtrck);
+ momentum.SetPy(p*TRACKR.cytrck);
+ momentum.SetPz(p*TRACKR.cztrck);
+ momentum.SetE(TRACKR.etrack);
+ return;
+ }
+ } else if (caller == kMGResumedTrack) {
+ momentum.SetPx(TRACKR.spausr[4]);
+ momentum.SetPy(TRACKR.spausr[5]);
+ momentum.SetPz(TRACKR.spausr[6]);
+ momentum.SetE (TRACKR.spausr[7]);
return;
- }
+ } else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
+ momentum.SetPx(0.);
+ momentum.SetPy(0.);
+ momentum.SetPz(0.);
+ momentum.SetE(TrackMass());
+
+ } else if (caller == kSODRAW) {
+ Int_t ist = FLKSTK.npflka;
+ Double_t p = FLKSTK.pmoflk[ist];
+ Int_t ifl = FLKSTK.iloflk[ist];
+ Double_t m = PAPROP.am[ifl + 6];
+ Double_t e = TMath::Sqrt(p * p + m * m);
+ momentum.SetPx(p * FLKSTK.txflk[ist]);
+ momentum.SetPy(p * FLKSTK.tyflk[ist]);
+ momentum.SetPz(p * FLKSTK.tzflk[ist]);
+ momentum.SetE(e);
+ } else if (caller == kUSDRAW) {
+ if (icode == kEMFSCObrems ||
+ icode == kEMFSCOmoller ||
+ icode == kEMFSCObhabha ||
+ icode == kEMFSCOcompton )
+ {
+ momentum.SetPx(fPint[0]);
+ momentum.SetPy(fPint[1]);
+ momentum.SetPz(fPint[2]);
+ momentum.SetE(fPint[3]);
+ } else if (icode == kKASKADdray ||
+ icode == kKASKADbrems ||
+ icode == kKASKADpair) {
+ momentum.SetPx(GENSTK.plr[0] * GENSTK.cxr[0]);
+ momentum.SetPy(GENSTK.plr[0] * GENSTK.cyr[0]);
+ momentum.SetPz(GENSTK.plr[0] * GENSTK.czr[0]);
+ momentum.SetE (GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6]);
+ } else {
+ Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack
+ - ParticleMassFPC(TRACKR.jtrack)
+ * ParticleMassFPC(TRACKR.jtrack));
+ momentum.SetPx(p*TRACKR.cxtrck);
+ momentum.SetPy(p*TRACKR.cytrck);
+ momentum.SetPz(p*TRACKR.cztrck);
+ momentum.SetE(TRACKR.etrack);
+ }
}
else
Warning("TrackMomentum","momentum not available");
// TRACKR.etrack = total energy of the particle
// TRACKR.jtrack = identity number of the particle
// PAPROP.am[TRACKR.jtrack] = particle mass in gev
- Int_t caller = GetCaller();
- if (caller != 2) { // not eedraw
+ FlukaCallerCode_t caller = GetCaller();
+ FlukaProcessCode_t icode = GetIcode();
+ if (caller != kEEDRAW &&
+ caller != kMGResumedTrack &&
+ caller != kSODRAW &&
+ caller != kUSDRAW &&
+ (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
if (TRACKR.ptrack >= 0) {
px = TRACKR.ptrack*TRACKR.cxtrck;
py = TRACKR.ptrack*TRACKR.cytrck;
pz = TRACKR.ptrack*TRACKR.cztrck;
- e = TRACKR.etrack;
+ e = TRACKR.etrack;
return;
}
else {
- Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
+ Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
px = p*TRACKR.cxtrck;
py = p*TRACKR.cytrck;
pz = p*TRACKR.cztrck;
- e = TRACKR.etrack;
+ e = TRACKR.etrack;
return;
}
+ } else if (caller == kMGResumedTrack) {
+ px = TRACKR.spausr[4];
+ py = TRACKR.spausr[5];
+ pz = TRACKR.spausr[6];
+ e = TRACKR.spausr[7];
+ return;
+ } else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
+ px = 0.;
+ py = 0.;
+ pz = 0.;
+ e = TrackMass();
+ } else if (caller == kSODRAW) {
+ Int_t ist = FLKSTK.npflka;
+ Double_t p = FLKSTK.pmoflk[ist];
+ Int_t ifl = FLKSTK.iloflk[ist];
+ Double_t m = PAPROP.am[ifl + 6];
+ e = TMath::Sqrt(p * p + m * m);
+ px = p * FLKSTK.txflk[ist];
+ py = p * FLKSTK.tyflk[ist];
+ pz = p * FLKSTK.tzflk[ist];
+ } else if (caller == kUSDRAW) {
+ if (icode == kEMFSCObrems ||
+ icode == kEMFSCOmoller ||
+ icode == kEMFSCObhabha ||
+ icode == kEMFSCOcompton )
+ {
+ px = fPint[0];
+ py = fPint[1];
+ pz = fPint[2];
+ e = fPint[3];
+ } else if (icode == kKASKADdray ||
+ icode == kKASKADbrems ||
+ icode == kKASKADpair) {
+ px = GENSTK.plr[0] * GENSTK.cxr[0];
+ py = GENSTK.plr[0] * GENSTK.cyr[0];
+ pz = GENSTK.plr[0] * GENSTK.czr[0];
+ e = GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6];
+ } else {
+ Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
+ px = p*TRACKR.cxtrck;
+ py = p*TRACKR.cytrck;
+ pz = p*TRACKR.cztrck;
+ e = TRACKR.etrack;
+ }
}
else
- Warning("TrackMomentum","momentum not available");
+ Warning("TrackMomentum","momentum not available");
}
//______________________________________________________________________________
{
// Return the length in centimeters of the current step
// TRACKR.ctrack = total curved path
- Int_t caller = GetCaller();
- if (caller == 11 || caller==12 || caller == 3 || caller == 6) //bxdraw,endraw,usdraw
- return 0.0;
- else if (caller == 4) //mgdraw
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kBXEntering || caller == kBXExiting ||
+ caller == kENDRAW || caller == kUSDRAW ||
+ caller == kUSTCKV || caller == kMGResumedTrack ||
+ caller == kSODRAW)
+ return 0.0;
+ else if (caller == kMGDRAW)
return TRACKR.ctrack;
- else
- return -1.0;
+ else {
+ Warning("TrackStep", "track step not available");
+ return 0.0;
+ }
}
//______________________________________________________________________________
Double_t TFluka::TrackLength() const
{
// TRACKR.cmtrck = cumulative curved path since particle birth
- Int_t caller = GetCaller();
- if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kBXEntering || caller == kBXExiting ||
+ caller == kENDRAW || caller == kUSDRAW || caller == kMGDRAW ||
+ caller == kUSTCKV)
return TRACKR.cmtrck;
- else
- return -1.0;
+ else if (caller == kMGResumedTrack)
+ return TRACKR.spausr[8];
+ else if (caller == kSODRAW)
+ return 0.0;
+ else {
+ Warning("TrackLength", "track length not available for caller %5d \n", caller);
+ return 0.0;
+ }
}
//______________________________________________________________________________
{
// Return the current time of flight of the track being transported
// TRACKR.atrack = age of the particle
- Int_t caller = GetCaller();
- if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kMGDRAW) {
+ Int_t i;
+ if ((i = fPrimaryElectronIndex) > -1) {
+ Double_t x, y, z, t;
+ GetPrimaryElectronPosition(i, x, y, z, t);
+ return t;
+ } else {
+ return TRACKR.atrack;
+ }
+ } else if (caller == kBXEntering || caller == kBXExiting ||
+ caller == kENDRAW || caller == kUSDRAW ||
+ caller == kUSTCKV)
return TRACKR.atrack;
- else
- return -1;
+ else if (caller == kMGResumedTrack)
+ return TRACKR.spausr[3];
+ else if (caller == kSODRAW) {
+ return (FLKSTK.agestk[FLKSTK.npflka]);
+ }
+ else {
+ Warning("TrackTime", "track time not available");
+ return 0.0;
+ }
}
//______________________________________________________________________________
// TRACKR.dtrack = energy deposition of the jth deposition event
// If coming from bxdraw we have 2 steps of 0 length and 0 edep
- Int_t caller = GetCaller();
- if (caller == 11 || caller==12) return 0.0;
+ // If coming from usdraw we just signal particle production - no edep
+ // If just first time after resuming, no edep for the primary
+ FlukaCallerCode_t caller = GetCaller();
+
+ if (caller == kBXExiting || caller == kBXEntering ||
+ caller == kUSDRAW || caller == kMGResumedTrack ||
+ caller == kSODRAW)
+ return 0.0;
Double_t sum = 0;
- for ( Int_t j=0;j<TRACKR.mtrack;j++) {
- sum +=TRACKR.dtrack[j];
- }
- if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
- return fRull + sum;
- else {
- return sum;
+ Int_t i = -1;
+
+ // Material with primary ionisation activated but number of primary electrons nprim = 0
+ if (fPrimaryElectronIndex == -2) return 0.0;
+ // nprim > 0
+ if ((i = fPrimaryElectronIndex) > -1) {
+ // Primary ionisation
+ sum = GetPrimaryElectronKineticEnergy(i);
+ if (sum > 100.) {
+ printf("edep > 100. %d %d %f \n", i, ALLDLT.nalldl, sum);
+ }
+ return sum;
+ } else {
+ // Normal ionisation
+ if (TRACKR.mtrack > 1) printf("Edep: %6d\n", TRACKR.mtrack);
+
+ for ( Int_t j=0;j<TRACKR.mtrack;j++) {
+ sum +=TRACKR.dtrack[j];
+ }
+ if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
+ return fRull + sum;
+ else {
+ return sum;
+ }
}
}
+//______________________________________________________________________________
+Int_t TFluka::CorrectFlukaId() const
+{
+ // since we don't put photons and e- created bellow transport cut on the vmc stack
+ // and there is a call to endraw for energy deposition for each of them
+ // and they have the track number of their parent, but different identity (pdg)
+ // so we want to assign also their parent identity.
+
+ if( (IsTrackStop())
+ && TRACKR.ispusr[mkbmx2 - 4] == TRACKR.ispusr[mkbmx2 - 1]
+ && TRACKR.jtrack != TRACKR.ispusr[mkbmx2 - 3] ) {
+ if (fVerbosityLevel >=3)
+ cout << "CorrectFlukaId() for icode=" << GetIcode()
+ << " track=" << TRACKR.ispusr[mkbmx2 - 1]
+ << " current PDG=" << PDGFromId(TRACKR.jtrack)
+ << " assign parent PDG=" << PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
+ return TRACKR.ispusr[mkbmx2 - 3]; // assign parent identity
+ }
+ if (TRACKR.jtrack <= 64){
+ return TRACKR.jtrack;
+ } else {
+ return TRACKR.j0trck;
+ }
+}
+
+
//______________________________________________________________________________
Int_t TFluka::TrackPid() const
{
// Return the id of the particle transported
// TRACKR.jtrack = identity number of the particle
- Int_t caller = GetCaller();
- if (caller != 2) { // not eedraw
- return PDGFromId(TRACKR.jtrack);
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller != kEEDRAW && caller != kSODRAW) {
+ return PDGFromId( CorrectFlukaId() );
+ }
+ else if (caller == kSODRAW) {
+ return PDGFromId(FLKSTK.iloflk[FLKSTK.npflka]);
}
else
return -1000;
// Return charge of the track currently transported
// PAPROP.ichrge = electric charge of the particle
// TRACKR.jtrack = identity number of the particle
- Int_t caller = GetCaller();
- if (caller != 2) // not eedraw
- return PAPROP.ichrge[TRACKR.jtrack+6];
+
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller != kEEDRAW && caller != kSODRAW)
+ return PAPROP.ichrge[CorrectFlukaId() + 6];
+ else if (caller == kSODRAW) {
+ Int_t ifl = PDGFromId(FLKSTK.iloflk[FLKSTK.npflka]);
+ return PAPROP.ichrge[ifl + 6];
+ }
else
return -1000.0;
}
{
// PAPROP.am = particle mass in GeV
// TRACKR.jtrack = identity number of the particle
- Int_t caller = GetCaller();
- if (caller != 2) // not eedraw
- return PAPROP.am[TRACKR.jtrack+6];
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller != kEEDRAW && caller != kSODRAW)
+ return PAPROP.am[CorrectFlukaId()+6];
+ else if (caller == kSODRAW) {
+ Int_t ifl = FLKSTK.iloflk[FLKSTK.npflka];
+ return PAPROP.am[ifl + 6];
+ }
else
return -1000.0;
}
Double_t TFluka::Etot() const
{
// TRACKR.etrack = total energy of the particle
- Int_t caller = GetCaller();
- if (caller != 2) // not eedraw
- return TRACKR.etrack;
- else
- return -1000.0;
+ FlukaCallerCode_t caller = GetCaller();
+ FlukaProcessCode_t icode = GetIcode();
+ if (caller != kEEDRAW && caller != kSODRAW && caller != kUSDRAW)
+ {
+ return TRACKR.etrack;
+ } else if (caller == kUSDRAW) {
+ if (icode == kEMFSCObrems ||
+ icode == kEMFSCOmoller ||
+ icode == kEMFSCObhabha ||
+ icode == kEMFSCOcompton ) {
+ return fPint[3];
+ }
+ else if (icode == kKASKADdray ||
+ icode == kKASKADbrems ||
+ icode == kKASKADpair) {
+ return (GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6]);
+ }
+ }
+ else if (caller == kSODRAW) {
+ Int_t ist = FLKSTK.npflka;
+ Double_t p = FLKSTK.pmoflk[ist];
+ Int_t ifl = FLKSTK.iloflk[ist];
+ Double_t m = PAPROP.am[ifl + 6];
+ Double_t e = TMath::Sqrt(p * p + m * m);
+ return e;
+ }
+
+ return -1000.0;
}
//
// If the step would go behind the region of one material,
// it will be shortened to reach only the boundary.
// Therefore IsTrackInside() is always true.
- Int_t caller = GetCaller();
- if (caller == 11 || caller==12) // bxdraw
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kBXEntering || caller == kBXExiting)
return 0;
else
return 1;
{
// True if this is the first step of the track in the current volume
- Int_t caller = GetCaller();
- if (caller == 11) // bxdraw entering
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kBXEntering)
return 1;
else return 0;
}
{
// True if track is exiting volume
//
- Int_t caller = GetCaller();
- if (caller == 12) // bxdraw exiting
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kBXExiting)
return 1;
else return 0;
}
{
// True if the track is out of the setup
// means escape
-// Icode = 14: escape - call from Kaskad
-// Icode = 23: escape - call from Emfsco
-// Icode = 32: escape - call from Kasneu
-// Icode = 40: escape - call from Kashea
-// Icode = 51: escape - call from Kasoph
- if (fIcode == 14 ||
- fIcode == 23 ||
- fIcode == 32 ||
- fIcode == 40 ||
- fIcode == 51) return 1;
+ FlukaProcessCode_t icode = GetIcode();
+
+ if (icode == kKASKADescape ||
+ icode == kEMFSCOescape ||
+ icode == kKASNEUescape ||
+ icode == kKASHEAescape ||
+ icode == kKASOPHescape)
+ return 1;
else return 0;
}
//______________________________________________________________________________
Bool_t TFluka::IsTrackDisappeared() const
{
-// means all inelastic interactions and decays
+// All inelastic interactions and decays
// fIcode from usdraw
- if (fIcode == 101 || // inelastic interaction
- fIcode == 102 || // particle decay
- fIcode == 103 || // delta ray generation by hadron
- fIcode == 104 || // direct pair production
- fIcode == 105 || // bremsstrahlung (muon)
- fIcode == 208 || // bremsstrahlung (electron)
- fIcode == 214 || // in-flight annihilation
- fIcode == 215 || // annihilation at rest
- fIcode == 217 || // pair production
- fIcode == 219 || // Compton scattering
- fIcode == 221 || // Photoelectric effect
- fIcode == 300 || // hadronic interaction
- fIcode == 400 // delta-ray
+ FlukaProcessCode_t icode = GetIcode();
+ if (icode == kKASKADinelint || // inelastic interaction
+ icode == kKASKADdecay || // particle decay
+ icode == kKASKADdray || // delta ray generation by hadron
+ icode == kKASKADpair || // direct pair production
+ icode == kKASKADbrems || // bremsstrahlung (muon)
+ icode == kEMFSCObrems || // bremsstrahlung (electron)
+ icode == kEMFSCOmoller || // Moller scattering
+ icode == kEMFSCObhabha || // Bhaba scattering
+ icode == kEMFSCOanniflight || // in-flight annihilation
+ icode == kEMFSCOannirest || // annihilation at rest
+ icode == kEMFSCOpair || // pair production
+ icode == kEMFSCOcompton || // Compton scattering
+ icode == kEMFSCOphotoel || // Photoelectric effect
+ icode == kKASNEUhadronic || // hadronic interaction
+ icode == kKASHEAdray // delta-ray
) return 1;
else return 0;
}
{
// True if the track energy has fallen below the threshold
// means stopped by signal or below energy threshold
-// Icode = 12: stopping particle - call from Kaskad
-// Icode = 15: time kill - call from Kaskad
-// Icode = 21: below threshold, iarg=1 - call from Emfsco
-// Icode = 22: below threshold, iarg=2 - call from Emfsco
-// Icode = 24: time kill - call from Emfsco
-// Icode = 31: below threshold - call from Kasneu
-// Icode = 33: time kill - call from Kasneu
-// Icode = 41: time kill - call from Kashea
-// Icode = 52: time kill - call from Kasoph
- if (fIcode == 12 ||
- fIcode == 15 ||
- fIcode == 21 ||
- fIcode == 22 ||
- fIcode == 24 ||
- fIcode == 31 ||
- fIcode == 33 ||
- fIcode == 41 ||
- fIcode == 52) return 1;
+ FlukaProcessCode_t icode = GetIcode();
+ if (icode == kKASKADstopping || // stopping particle
+ icode == kKASKADtimekill || // time kill
+ icode == kEMFSCOstopping1 || // below user-defined cut-off
+ icode == kEMFSCOstopping2 || // below user cut-off
+ icode == kEMFSCOtimekill || // time kill
+ icode == kKASNEUstopping || // neutron below threshold
+ icode == kKASNEUtimekill || // time kill
+ icode == kKASHEAtimekill || // time kill
+ icode == kKASOPHtimekill) return 1; // time kill
else return 0;
}
{
// Number of secondary particles generated in the current step
-// FINUC.np = number of secondaries except light and heavy ions
+// GENSTK.np = number of secondaries except light and heavy ions
// FHEAVY.npheav = number of secondaries for light and heavy secondary ions
- Int_t caller = GetCaller();
- if (caller == 6) // valid only after usdraw
- return FINUC.np + FHEAVY.npheav;
- else if (caller == 50) {
- // Cerenkov Photon production
- return fNCerenkov;
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kUSDRAW) // valid only after usdraw
+ return GENSTK.np + FHEAVY.npheav;
+ else if (caller == kUSTCKV) {
+ // Cerenkov Photon production
+ return fNCerenkov;
}
return 0;
} // end of NSecondaries
//______________________________________________________________________________
void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
- TLorentzVector& position, TLorentzVector& momentum)
+ TLorentzVector& position, TLorentzVector& momentum)
{
// Copy particles from secondary stack to vmc stack
//
- Int_t caller = GetCaller();
- if (caller == 6) { // valid only after usdraw
- if (FINUC.np > 0) {
- // Hadronic interaction
- if (isec >= 0 && isec < FINUC.np) {
- particleId = PDGFromId(FINUC.kpart[isec]);
- position.SetX(fXsco);
- position.SetY(fYsco);
- position.SetZ(fZsco);
- position.SetT(TRACKR.atrack);
- momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
- momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
- momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
- momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
- }
- else if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
- Int_t jsec = isec - FINUC.np;
- particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
- position.SetX(fXsco);
- position.SetY(fYsco);
- position.SetZ(fZsco);
- position.SetT(TRACKR.atrack);
- momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
- momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
- momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
- if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
- momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
- else if (FHEAVY.tkheav[jsec] > 6)
- momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
- }
- else
- Warning("GetSecondary","isec out of range");
- }
- } else if (caller == 50) {
- Int_t index = OPPHST.lstopp - isec;
- position.SetX(OPPHST.xoptph[index]);
- position.SetY(OPPHST.yoptph[index]);
- position.SetZ(OPPHST.zoptph[index]);
- position.SetT(OPPHST.agopph[index]);
- Double_t p = OPPHST.poptph[index];
-
- momentum.SetPx(p * OPPHST.txopph[index]);
- momentum.SetPy(p * OPPHST.tyopph[index]);
- momentum.SetPz(p * OPPHST.tzopph[index]);
- momentum.SetE(p);
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kUSDRAW) { // valid only after usdraw
+ if (GENSTK.np > 0) {
+ // Hadronic interaction
+ if (isec >= 0 && isec < GENSTK.np) {
+ particleId = PDGFromId(GENSTK.kpart[isec]);
+ position.SetX(fXsco);
+ position.SetY(fYsco);
+ position.SetZ(fZsco);
+ position.SetT(TRACKR.atrack);
+ momentum.SetPx(GENSTK.plr[isec]*GENSTK.cxr[isec]);
+ momentum.SetPy(GENSTK.plr[isec]*GENSTK.cyr[isec]);
+ momentum.SetPz(GENSTK.plr[isec]*GENSTK.czr[isec]);
+ momentum.SetE(GENSTK.tki[isec] + PAPROP.am[GENSTK.kpart[isec]+6]);
+ }
+ else if (isec >= GENSTK.np && isec < GENSTK.np + FHEAVY.npheav) {
+ Int_t jsec = isec - GENSTK.np;
+ particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
+ position.SetX(fXsco);
+ position.SetY(fYsco);
+ position.SetZ(fZsco);
+ position.SetT(TRACKR.atrack);
+ momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
+ momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
+ momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
+ if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
+ momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
+ else if (FHEAVY.tkheav[jsec] > 6)
+ momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
+ }
+ else
+ Warning("GetSecondary","isec out of range");
+ }
+ } else if (caller == kUSTCKV) {
+ Int_t index = OPPHST.lstopp - isec;
+ position.SetX(OPPHST.xoptph[index]);
+ position.SetY(OPPHST.yoptph[index]);
+ position.SetZ(OPPHST.zoptph[index]);
+ position.SetT(OPPHST.agopph[index]);
+ Double_t p = OPPHST.poptph[index];
+
+ momentum.SetPx(p * OPPHST.txopph[index]);
+ momentum.SetPy(p * OPPHST.tyopph[index]);
+ momentum.SetPz(p * OPPHST.tzopph[index]);
+ momentum.SetE(p);
}
else
- Warning("GetSecondary","no secondaries available");
+ Warning("GetSecondary","no secondaries available");
} // end of GetSecondary
// Name of the process that has produced the secondary particles
// in the current step
- Int_t mugamma = (TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11);
-
- if (fIcode == 102) return kPDecay;
- else if (fIcode == 104 || fIcode == 217) return kPPair;
- else if (fIcode == 219) return kPCompton;
- else if (fIcode == 221) return kPPhotoelectric;
- else if (fIcode == 105 || fIcode == 208) return kPBrem;
- else if (fIcode == 103 || fIcode == 400) return kPDeltaRay;
- else if (fIcode == 210 || fIcode == 212) return kPDeltaRay;
- else if (fIcode == 214 || fIcode == 215) return kPAnnihilation;
- else if (fIcode == 101) return kPHadronic;
- else if (fIcode == 101) {
- if (!mugamma) return kPHadronic;
- else if (TRACKR.jtrack == 7) return kPPhotoFission;
- else return kPMuonNuclear;
+ Int_t mugamma = (TRACKR.jtrack == kFLUKAphoton ||
+ TRACKR.jtrack == kFLUKAmuplus ||
+ TRACKR.jtrack == kFLUKAmuminus);
+ FlukaProcessCode_t icode = GetIcode();
+
+ if (icode == kKASKADdecay) return kPDecay;
+ else if (icode == kKASKADpair || icode == kEMFSCOpair) return kPPair;
+ else if (icode == kEMFSCOcompton) return kPCompton;
+ else if (icode == kEMFSCOphotoel) return kPPhotoelectric;
+ else if (icode == kKASKADbrems || icode == kEMFSCObrems) return kPBrem;
+ else if (icode == kKASKADdray || icode == kKASHEAdray) return kPDeltaRay;
+ else if (icode == kEMFSCOmoller || icode == kEMFSCObhabha) return kPDeltaRay;
+ else if (icode == kEMFSCOanniflight || icode == kEMFSCOannirest) return kPAnnihilation;
+ else if (icode == kKASKADinelint) {
+ if (!mugamma) return kPHadronic;
+ else if (TRACKR.jtrack == kFLUKAphoton) return kPPhotoFission;
+ else return kPMuonNuclear;
}
- else if (fIcode == 225) return kPRayleigh;
+ else if (icode == kEMFSCOrayleigh) return kPRayleigh;
// Fluka codes 100, 300 and 400 still to be investigasted
- else return kPNoProcess;
+ else return kPNoProcess;
}
//
// Return processes active in the current step
//
+ FlukaProcessCode_t icode = GetIcode();
+ FlukaCallerCode_t caller = GetCaller();
+
proc.Set(1);
TMCProcess iproc;
- switch (fIcode) {
- case 15:
- case 24:
- case 33:
- case 41:
- case 52:
- iproc = kPTOFlimit;
- break;
- case 12:
- case 14:
- case 21:
- case 22:
- case 23:
- case 31:
- case 32:
- case 40:
- case 51:
- iproc = kPStop;
- break;
- case 50:
- iproc = kPLightAbsorption;
- break;
- case 20:
- iproc = kPPhotoelectric;
- break;
- default:
- iproc = ProdProcess(0);
+ if (caller == kBXEntering || caller == kBXExiting || caller == kEEDRAW || caller == kSODRAW) {
+ iproc = kPTransportation;
}
+ else if (caller == kUSTCKV) {
+ iproc = kPCerenkov;
+ } else {
+ switch (icode) {
+ case kEMFSCO:
+ if (Edep() > 0.) {
+ iproc = kPEnergyLoss;
+ } else {
+ iproc = kPTransportation;
+ }
+ break;
+ case kKASKAD:
+ if (Edep() > 0.) {
+ iproc = kPEnergyLoss;
+ } else {
+ iproc = kPTransportation;
+ }
+ break;
+ case kKASHEA:
+ case kKASNEU:
+ case kKASOPH:
+ case kKASKADescape:
+ case kEMFSCOescape:
+ case kKASNEUescape:
+ case kKASHEAescape:
+ case kKASOPHescape:
+ iproc = kPTransportation;
+ break;
+ case kKASKADtimekill:
+ case kEMFSCOtimekill:
+ case kKASNEUtimekill:
+ case kKASHEAtimekill:
+ case kKASOPHtimekill:
+ iproc = kPTOFlimit;
+ break;
+ case kKASKADstopping:
+ case kEMFSCOstopping1:
+ case kEMFSCOstopping2:
+ case kKASNEUstopping:
+ iproc = kPStop;
+ break;
+ case kKASOPHabsorption:
+ iproc = kPLightAbsorption;
+ break;
+ case kKASOPHrefraction:
+ iproc = kPLightRefraction;
+ break;
+ case kEMFSCOlocaldep :
+ iproc = kPPhotoelectric;
+ break;
+ default:
+ iproc = ProdProcess(0);
+ }
+ }
+
proc[0] = iproc;
return 1;
}
return fMCGeo->VolName(id);
}
+Int_t TFluka::MediumId(const Text_t* mediumName) const
+{
+ //
+ // Return the unique medium id for medium with name mediumName
+ TList *medlist = gGeoManager->GetListOfMedia();
+ TGeoMedium* med = (TGeoMedium*) medlist->FindObject(mediumName);
+ if (med) {
+ return (med->GetId());
+ } else {
+ return (-1);
+ }
+}
+
//______________________________________________________________________________
Int_t TFluka::VolId(const Text_t* volName) const
{
// Time consuming. (Only used during set-up)
// Could be replaced by hash-table
//
- return fMCGeo->VolId(volName);
+ char sname[20];
+ Int_t len;
+ strncpy(sname, volName, len = strlen(volName));
+ sname[len] = 0;
+ while (sname[len - 1] == ' ') sname[--len] = 0;
+ return fMCGeo->VolId(sname);
}
//______________________________________________________________________________
//
// Return the current volume name
//
- if (gGeoManager->IsOutside()) return 0;
+ if (gGeoManager->IsOutside()) return "OutOfWorld";
return gGeoManager->GetCurrentVolume()->GetName();
}
return node->GetVolume()->GetName();
}
+const char* TFluka::CurrentVolPath() {
+ // Return the current volume path
+ return gGeoManager->GetPath();
+}
//______________________________________________________________________________
-Int_t TFluka::CurrentMaterial(Float_t & /*a*/, Float_t & /*z*/,
- Float_t & /*dens*/, Float_t & /*radl*/, Float_t & /*absl*/) const
+Int_t TFluka::CurrentMaterial(Float_t & a, Float_t & z,
+ Float_t & dens, Float_t & radl, Float_t & absl) const
{
//
-// Return the current medium number ??? what about material properties
+// Return the current medium number and material properties
//
Int_t copy;
Int_t id = TFluka::CurrentVolID(copy);
Int_t med = TFluka::VolId2Mate(id);
+ TGeoVolume* vol = gGeoManager->GetCurrentVolume();
+ TGeoMaterial* mat = vol->GetMaterial();
+ a = mat->GetA();
+ z = mat->GetZ();
+ dens = mat->GetDensity();
+ radl = mat->GetRadLen();
+ absl = mat->GetIntLen();
+
return med;
}
//______________________________________________________________________________
void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
{
+//
+// See Gmtod(Float_t*, Float_t*, Int_t)
+//
if (iflag == 1) gGeoManager->MasterToLocal(xm,xd);
else gGeoManager->MasterToLocalVect(xm,xd);
}
//______________________________________________________________________________
void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
{
+//
+// See Gdtom(Float_t*, Float_t*, Int_t)
+//
if (iflag == 1) gGeoManager->LocalToMaster(xd,xm);
else gGeoManager->LocalToMasterVect(xd,xm);
}
//______________________________________________________________________________
TObjArray *TFluka::GetFlukaMaterials()
{
+//
+// Get array of Fluka materials
return fGeom->GetMatList();
}
//______________________________________________________________________________
-void TFluka::SetMreg(Int_t l)
+void TFluka::SetMreg(Int_t l, Int_t lttc)
{
// Set current fluka region
fCurrentFlukaRegion = l;
- fGeom->SetMreg(l);
+ fGeom->SetMreg(l,lttc);
}
+//______________________________________________________________________________
TString TFluka::ParticleName(Int_t pdg) const
{
// Return particle name for particle with pdg code pdg.
Int_t ifluka = IdFromPDG(pdg);
- return TString((CHPPRP.btype[ifluka+6]), 8);
+ return TString((CHPPRP.btype[ifluka - kFLUKAcodemin]), 8);
}
+//______________________________________________________________________________
Double_t TFluka::ParticleMass(Int_t pdg) const
{
// Return particle mass for particle with pdg code pdg.
Int_t ifluka = IdFromPDG(pdg);
- return (PAPROP.am[ifluka+6]);
+ return (PAPROP.am[ifluka - kFLUKAcodemin]);
}
+//______________________________________________________________________________
+Double_t TFluka::ParticleMassFPC(Int_t fpc) const
+{
+ // Return particle mass for particle with Fluka particle code fpc
+ return (PAPROP.am[fpc - kFLUKAcodemin]);
+}
+
+//______________________________________________________________________________
Double_t TFluka::ParticleCharge(Int_t pdg) const
{
// Return particle charge for particle with pdg code pdg.
Int_t ifluka = IdFromPDG(pdg);
- return Double_t(PAPROP.ichrge[ifluka+6]);
+ return Double_t(PAPROP.ichrge[ifluka - kFLUKAcodemin]);
}
+//______________________________________________________________________________
Double_t TFluka::ParticleLifeTime(Int_t pdg) const
{
// Return particle lifetime for particle with pdg code pdg.
Int_t ifluka = IdFromPDG(pdg);
- return (PAPROP.thalf[ifluka+6]);
+ return (PAPROP.tmnlf[ifluka - kFLUKAcodemin]);
}
+//______________________________________________________________________________
void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t& charge, Float_t& tlife)
{
// Retrieve particle properties for particle with pdg code pdg.
tlife = ParticleLifeTime(pdg);
}
+//______________________________________________________________________________
+void TFluka::PrintHeader()
+{
+ //
+ // Print a header
+ printf("\n");
+ printf("\n");
+ printf("------------------------------------------------------------------------------\n");
+ printf("- You are using the TFluka Virtual Monte Carlo Interface to FLUKA. -\n");
+ printf("- Please see the file fluka.out for FLUKA output and licensing information. -\n");
+ printf("------------------------------------------------------------------------------\n");
+ printf("\n");
+ printf("\n");
+}
-#define pushcerenkovphoton pushcerenkovphoton_
-#define usersteppingckv usersteppingckv_
+#define pshckp pshckp_
+#define ustckv ustckv_
extern "C" {
- void pushcerenkovphoton(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
- Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
- Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
+ void pshckp(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
+ Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
+ Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
+ {
+ //
+ // Pushes one cerenkov photon to the stack
+ //
+
+ TFluka* fluka = (TFluka*) gMC;
+ TVirtualMCStack* cppstack = fluka->GetStack();
+ Int_t parent = TRACKR.ispusr[mkbmx2-1];
+ cppstack->PushTrack(0, parent, 50000050,
+ px, py, pz, e,
+ vx, vy, vz, tof,
+ polx, poly, polz,
+ kPCerenkov, ntr, wgt, 0);
+ if (fluka->GetVerbosityLevel() >= 3)
+ printf("pshckp: track=%d parent=%d lattc=%d %s\n", ntr, parent, TRACKR.lt1trk, fluka->CurrentVolName());
+ }
+
+ void ustckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
{
- //
- // Pushes one cerenkov photon to the stack
- //
-
- TFluka* fluka = (TFluka*) gMC;
- TVirtualMCStack* cppstack = fluka->GetStack();
- Int_t parent = TRACKR.ispusr[mkbmx2-1];
- cppstack->PushTrack(0, parent, 50000050,
- px, py, pz, e,
- vx, vy, vz, tof,
- polx, poly, polz,
- kPCerenkov, ntr, wgt, 0);
+ //
+ // Calls stepping in order to signal cerenkov production
+ //
+ TFluka *fluka = (TFluka*)gMC;
+ fluka->SetMreg(mreg, TRACKR.lt1trk); //LTCLCM.mlatm1);
+ fluka->SetXsco(x);
+ fluka->SetYsco(y);
+ fluka->SetZsco(z);
+ fluka->SetNCerenkov(nphot);
+ fluka->SetCaller(kUSTCKV);
+ if (fluka->GetVerbosityLevel() >= 3)
+ printf("ustckv: %10d mreg=%d lattc=%d newlat=%d (%f, %f, %f) edep=%f vol=%s\n",
+ nphot, mreg, TRACKR.lt1trk, LTCLCM.newlat, x, y, z, fluka->Edep(), fluka->CurrentVolName());
+
+ // check region lattice consistency (debug Ernesto)
+ // *****************************************************
+ Int_t nodeId;
+ Int_t volId = fluka->CurrentVolID(nodeId);
+ Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
+
+ if( mreg != volId && !gGeoManager->IsOutside() ) {
+ cout << " ustckv: track=" << TRACKR.ispusr[mkbmx2-1] << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
+ << " icode=" << fluka->GetIcode() << " gNstep=" << fluka->GetNstep() << endl
+ << " fluka mreg=" << mreg << " mlttc=" << TRACKR.lt1trk << endl
+ << " TGeo volId=" << volId << " crtlttc=" << crtlttc << endl
+ << " common TRACKR lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
+ << " common LTCLCM newlat=" << LTCLCM.newlat << " mlatld=" << LTCLCM.mlatld << endl
+ << " mlatm1=" << LTCLCM.mlatm1 << " mltsen=" << LTCLCM.mltsen << endl
+ << " mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
+ if( TRACKR.lt1trk == crtlttc ) cout << " *************************************************************" << endl;
}
+ // *****************************************************
- void usersteppingckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
- {
- //
- // Calls stepping in order to signal cerenkov production
- //
- TFluka *fluka = (TFluka*)gMC;
- fluka->SetMreg(mreg);
- fluka->SetXsco(x);
- fluka->SetYsco(y);
- fluka->SetZsco(z);
- fluka->SetNCerenkov(nphot);
- fluka->SetCaller(50);
- printf("userstepping ckv: %10d %10d %13.3f %13.3f %13.2f\n", nphot, mreg, x, y, z);
+
+
+ (TVirtualMCApplication::Instance())->Stepping();
+ }
+}
+
+//______________________________________________________________________________
+void TFluka::AddParticlesToPdgDataBase() const
+{
+
+//
+// Add particles to the PDG data base
+
+ TDatabasePDG *pdgDB = TDatabasePDG::Instance();
+
+ const Double_t kAu2Gev = 0.9314943228;
+ const Double_t khSlash = 1.0545726663e-27;
+ const Double_t kErg2Gev = 1/1.6021773349e-3;
+ const Double_t khShGev = khSlash*kErg2Gev;
+ const Double_t kYear2Sec = 3600*24*365.25;
+//
+// Ions
+//
+ pdgDB->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE,
+ 0,3,"Ion",GetIonPdg(1,2));
+ pdgDB->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE,
+ khShGev/(12.33*kYear2Sec),3,"Ion",GetIonPdg(1,3));
+ pdgDB->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,
+ khShGev/(12.33*kYear2Sec),6,"Ion",GetIonPdg(2,4));
+ pdgDB->AddParticle("HE3","HE3",3*kAu2Gev+14.931e-3,kFALSE,
+ 0,6,"Ion",GetIonPdg(2,3));
+//
+//
+//
+// Special particles
+//
+ pdgDB->AddParticle("Cherenkov","Cherenkov",0,kFALSE,
+ 0,0,"Special",GetSpecialPdg(50));
+ pdgDB->AddParticle("FeedbackPhoton","FeedbackPhoton",0,kFALSE,
+ 0,0,"Special",GetSpecialPdg(51));
+}
+
+//
+// Info about primary ionization electrons
+//
+
+//______________________________________________________________________________
+Int_t TFluka::GetNPrimaryElectrons()
+{
+ // Get number of primary electrons
+ return ALLDLT.nalldl;
+}
+
+//______________________________________________________________________________
+Double_t TFluka::GetPrimaryElectronKineticEnergy(Int_t i) const
+{
+ // Returns kinetic energy of primary electron i
+
+ Double_t ekin = -1.;
+
+ if (i >= 0 && i < ALLDLT.nalldl) {
+ ekin = ALLDLT.talldl[i];
+ } else {
+ Warning("GetPrimaryElectronKineticEnergy",
+ "Primary electron index out of range %d %d \n",
+ i, ALLDLT.nalldl);
+ }
+ return ekin;
+}
+
+void TFluka::GetPrimaryElectronPosition(Int_t i, Double_t& x, Double_t& y, Double_t& z, Double_t& t) const
+{
+ // Returns position of primary electron i
+ if (i >= 0 && i < ALLDLT.nalldl) {
+ x = ALLDLT.xalldl[i];
+ y = ALLDLT.yalldl[i];
+ z = ALLDLT.zalldl[i];
+ t = ALLDLT.talldl[i];
+ return;
+ } else {
+ Warning("GetPrimaryElectronPosition",
+ "Primary electron index out of range %d %d \n",
+ i, ALLDLT.nalldl);
+ return;
+ }
+ return;
+}
+
+Int_t TFluka::GetIonPdg(Int_t z, Int_t a, Int_t i) const
+{
+// Acording to
+// http://cepa.fnal.gov/psm/stdhep/pdg/montecarlorpp-2006.pdf
+
+ return 1000000000 + 10*1000*z + 10*a + i;
+}
+
+//__________________________________________________________________
+Int_t TFluka::GetSpecialPdg(Int_t number) const
+{
+// Numbering for special particles
+
+ return 50000000 + number;
+}
+
+
+void TFluka::PrimaryIonisationStepping(Int_t nprim)
+{
+// Call Stepping for primary ionisation electrons
+ Int_t i;
+// Protection against nprim > mxalld
+
+// Multiple steps for nprim > 0
+ if (nprim > 0) {
+ for (i = 0; i < nprim; i++) {
+ SetCurrentPrimaryElectronIndex(i);
+ (TVirtualMCApplication::Instance())->Stepping();
+ if (i == 0) SetTrackIsNew(kFALSE);
+ }
+ } else {
+ // No primary electron ionisation
+ // Call Stepping anyway but flag nprim = 0 as index = -2
+ SetCurrentPrimaryElectronIndex(-2);
(TVirtualMCApplication::Instance())->Stepping();
}
+ // Reset the index
+ SetCurrentPrimaryElectronIndex(-1);
}
+//______________________________________________________________________
+Float_t* TFluka::CreateFloatArray(Double_t* array, Int_t size) const
+{
+// Converts Double_t* array to Float_t*,
+// !! The new array has to be deleted by user.
+// ---
+
+ Float_t* floatArray;
+ if (size>0) {
+ floatArray = new Float_t[size];
+ for (Int_t i=0; i<size; i++)
+ if (array[i] >= FLT_MAX )
+ floatArray[i] = FLT_MAX/100.;
+ else
+ floatArray[i] = array[i];
+ }
+ else {
+ //floatArray = 0;
+ floatArray = new Float_t[1];
+ }
+ return floatArray;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff