#include <Riostream.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 "TVirtualMC.h"
#include "TMCProcess.h"
#include "TFlukaScoringOption.h"
#include "TLorentzVector.h"
#include "TArrayI.h"
+#include "TArrayD.h"
// Fluka methods that may be needed.
#ifndef WIN32
fGeneratePemf = kFALSE;
fNVolumes = 0;
fCurrentFlukaRegion = -1;
+ fNewReg = -1;
fGeom = 0;
fMCGeo = 0;
fMaterials = 0;
SetGeneratePemf(kFALSE);
fNVolumes = 0;
fCurrentFlukaRegion = -1;
+ fNewReg = -1;
fDummyBoundary = 0;
fFieldFlag = 1;
fGeneratePemf = kFALSE;
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
-
- // Now we have TGeo geometry created and we have to patch FlukaVmc.inp
- // with the material mapping file FlukaMat.inp
-
+ 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) {
// Prepare input file with the current physics settings
InitPhysics();
-
- cout << "\t* InitPhysics() - Prepare input file was called" << endl;
-
- if (fVerbosityLevel >=2)
- cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
- << ") in fluka..." << endl;
- GLOBAL.lfdrtr = true;
-
- if (fVerbosityLevel >=2)
- cout << "\t* Opening file " << fInputFileName << endl;
+// Open fortran files
const char* fname = fInputFileName;
-
fluka_openinp(lunin, PASSCHARA(fname));
fluka_openout(11, PASSCHARA("fluka.out"));
-
- if (fVerbosityLevel >=2)
- cout << "\t* Calling flukam..." << endl;
+// Read input cards
+ GLOBAL.lfdrtr = true;
flukam(1);
-
- if (fVerbosityLevel >=2)
- cout << "\t* Closing file " << fInputFileName << 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");
+ Warning("ProcessEvent", "User Run Abortion: No more events handled !\n");
fNEvent += 1;
return;
}
if (fVerbosityLevel >=3)
cout << "==> TFluka::ProcessEvent() called." << endl;
fApplication->GeneratePrimaries();
- EPISOR.lsouit = true;
+ SOURCM.lsouit = true;
flukam(1);
if (fVerbosityLevel >=3)
cout << "<== TFluka::ProcessEvent() called." << endl;
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;
}
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(fGeom->GetFlukaMaterial(itmed), nreg);
- if (nreg == 0) {
- return;
- }
-
//
Bool_t process = kFALSE;
if (strncmp(param, "DCAY", 4) == 0 ||
{
process = kTRUE;
}
+
if (process) {
- SetProcess(param, Int_t (parval), fGeom->GetFlukaMaterial(itmed));
+ SetProcess(param, Int_t (parval), itmed);
} else {
- SetCut(param, parval, fGeom->GetFlukaMaterial(itmed));
+ SetCut(param, parval, itmed);
}
}
// 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) {
Int_t /*number*/, Int_t /*nlevel*/) {
//
// Not with TGeo
- Warning("WriteEuclid", "Not implemented with TGeo");
+ Warning("WriteEuclid", "Not implemented !");
}
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
Int_t idSpecial[6] = {10020040, 10020030, 10010030, 10010020, 10000000, 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;
}
// Error id
- if (id == 0 || id < -6 || id > 250) {
+ if (id == 0 || id < kFLUKAcodemin || id > kFLUKAcodemax) {
if (fVerbosityLevel >= 3)
printf("PDGFromId: Error id = 0\n");
return -1;
}
if (fVerbosityLevel >= 3)
printf("mpdgha called with %d %d \n", id, intfluka);
- // MPDGHA() goes from fluka internal to pdg.
return mpdgha(intfluka);
} else {
// ions and optical photons
- return idSpecial[id + 6];
+ return idSpecial[id - kFLUKAcodemin];
}
}
//______________________________________________________________________________
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");
-
// Construct file names
FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
TString sFlukaVmcCoreInp = getenv("ALICE_ROOT");
// Open files
if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
- printf("\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
+ Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
exit(1);
}
if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
- printf("\nCannot open file %s\n",sFlukaVmcTmp.Data());
+ Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcTmp.Data());
exit(1);
}
if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
- printf("\nCannot open file %s\n",sFlukaVmcInp.Data());
+ Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcInp.Data());
exit(1);
}
Int_t inp = 0;
Int_t nscore = fUserScore->GetEntries();
- TFlukaScoringOption *mopo = 0x0;
- TFlukaScoringOption *mopi = 0x0;
+ TFlukaScoringOption *mopo = 0;
+ TFlukaScoringOption *mopi = 0;
for (Int_t isc = 0; isc < nscore; isc++)
{
}
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");
// 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 || caller == 50) { //bxdraw,endraw,usdraw,ckov
+ 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,mgdraw resuming
+ else if (caller == kMGDRAW) {
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
+ else if (caller == kSODRAW) {
position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
position.SetT(0);
- } else if (caller == 40) { // mgdraw resuming transport
+ } else if (caller == kMGResumedTrack) {
position.SetX(TRACKR.spausr[0]);
position.SetY(TRACKR.spausr[1]);
position.SetZ(TRACKR.spausr[2]);
// 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 || caller == 50) { //bxdraw,endraw,usdraw,ckov
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kENDRAW || caller == kUSDRAW ||
+ caller == kBXExiting || caller == kBXEntering ||
+ caller == kUSTCKV) {
x = GetXsco();
y = GetYsco();
z = GetZsco();
}
- else if (caller == 4 || caller == 5) { // mgdraw, sodraw, mgdraw resuming
+ else if (caller == kMGDRAW || caller == kSODRAW) {
x = TRACKR.xtrack[TRACKR.ntrack];
y = TRACKR.ytrack[TRACKR.ntrack];
z = TRACKR.ztrack[TRACKR.ntrack];
}
- else if (caller == 40) { // mgdraw resuming transport
+ else if (caller == kMGResumedTrack) {
x = TRACKR.spausr[0];
y = TRACKR.spausr[1];
z = TRACKR.spausr[2];
// 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 && caller != 40) { // not eedraw or mgdraw resuming
+ FlukaCallerCode_t caller = GetCaller();
+ FlukaProcessCode_t icode = GetIcode();
+
+ if (caller != kEEDRAW && caller != kMGResumedTrack &&
+ (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
if (TRACKR.ptrack >= 0) {
momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
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));
momentum.SetPx(p*TRACKR.cxtrck);
momentum.SetPy(p*TRACKR.cytrck);
momentum.SetPz(p*TRACKR.cztrck);
momentum.SetE(TRACKR.etrack);
return;
}
- } else if (caller == 40) { // mgdraw resuming transport
+ } 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
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 && caller != 40) { // not eedraw and not mgdraw resuming
+ FlukaCallerCode_t caller = GetCaller();
+ FlukaProcessCode_t icode = GetIcode();
+ if (caller != kEEDRAW && caller != kMGResumedTrack &&
+ (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 == 40) { // mgdraw resuming transport
+ } 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
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 || caller == 50 || caller == 40) //bxdraw,endraw,usdraw, ckov
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kBXEntering || caller == kBXExiting ||
+ caller == kENDRAW || caller == kUSDRAW ||
+ caller == kUSTCKV || caller == kMGResumedTrack)
return 0.0;
- else if (caller == 4) //mgdraw
+ else if (caller == kMGDRAW)
return TRACKR.ctrack;
else {
Warning("TrackStep", "track step not available");
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 || caller == 50) //bxdraw,endraw,mgdraw,usdraw,ckov
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kBXEntering || caller == kBXExiting ||
+ caller == kENDRAW || caller == kUSDRAW || caller == kMGDRAW ||
+ caller == kUSTCKV)
return TRACKR.cmtrck;
- else if (caller == 40) // mgdraw resuming transport
+ else if (caller == kMGResumedTrack)
return TRACKR.spausr[8];
else {
Warning("TrackLength", "track length not available");
{
// 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 || caller == 50) //bxdraw,endraw,mgdraw,usdraw,ckov
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kBXEntering || caller == kBXExiting ||
+ caller == kENDRAW || caller == kUSDRAW || caller == kMGDRAW ||
+ caller == kUSTCKV)
return TRACKR.atrack;
- else if (caller == 40)
+ else if (caller == kMGResumedTrack)
return TRACKR.spausr[3];
else {
Warning("TrackTime", "track time not available");
// If coming from bxdraw we have 2 steps of 0 length and 0 edep
// If coming from usdraw we just signal particle production - no edep
// If just first time after resuming, no edep for the primary
- Int_t caller = GetCaller();
- if (caller == 11 || caller==12 || caller==6 || caller == 40) return 0.0;
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kBXExiting || caller == kBXEntering ||
+ caller == kUSDRAW || caller == kMGResumedTrack) return 0.0;
Double_t sum = 0;
for ( Int_t j=0;j<TRACKR.mtrack;j++) {
sum +=TRACKR.dtrack[j];
{
// Return the id of the particle transported
// TRACKR.jtrack = identity number of the particle
- Int_t caller = GetCaller();
- if (caller != 2) { // not eedraw
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller != kEEDRAW) {
return PDGFromId(TRACKR.jtrack);
}
else
// 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
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller != kEEDRAW)
return PAPROP.ichrge[TRACKR.jtrack+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
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller != kEEDRAW)
return PAPROP.am[TRACKR.jtrack+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
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller != kEEDRAW)
return TRACKR.etrack;
else
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 ||
+ icode == kKASKADtimekill ||
+ icode == kEMFSCOstopping1 ||
+ icode == kEMFSCOstopping2 ||
+ icode == kEMFSCOtimekill ||
+ icode == kKASNEUstopping ||
+ icode == kKASNEUtimekill ||
+ icode == kKASHEAtimekill ||
+ icode == kKASOPHtimekill) return 1;
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) {
+ 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;
}
// Copy particles from secondary stack to vmc stack
//
- Int_t caller = GetCaller();
- if (caller == 6) { // valid only after usdraw
- if (FINUC.np > 0) {
+ FlukaCallerCode_t caller = GetCaller();
+ if (caller == kUSDRAW) { // valid only after usdraw
+ if (GENSTK.np > 0) {
// Hadronic interaction
- if (isec >= 0 && isec < FINUC.np) {
- particleId = PDGFromId(FINUC.kpart[isec]);
+ 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(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]);
+ 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 >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
- Int_t jsec = isec - FINUC.np;
+ 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);
else
Warning("GetSecondary","isec out of range");
}
- } else if (caller == 50) {
+ } else if (caller == kUSTCKV) {
Int_t index = OPPHST.lstopp - isec;
position.SetX(OPPHST.xoptph[index]);
position.SetY(OPPHST.yoptph[index]);
// 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();
proc.Set(1);
TMCProcess iproc;
- switch (fIcode) {
- case 15:
- case 24:
- case 33:
- case 41:
- case 52:
+ switch (icode) {
+ case kKASKADtimekill:
+ case kEMFSCOtimekill:
+ case kKASNEUtimekill:
+ case kKASHEAtimekill:
+ case kKASOPHtimekill:
iproc = kPTOFlimit;
break;
- case 12:
- case 14:
- case 21:
- case 22:
- case 23:
- case 31:
- case 32:
- case 40:
- case 51:
+ case kKASKADstopping:
+ case kKASKADescape:
+ case kEMFSCOstopping1:
+ case kEMFSCOstopping2:
+ case kEMFSCOescape:
+ case kKASNEUstopping:
+ case kKASNEUescape:
+ case kKASHEAescape:
+ case kKASOPHescape:
iproc = kPStop;
break;
- case 50:
+ case kKASOPHabsorption:
iproc = kPLightAbsorption;
break;
- case 59:
+ case kKASOPHrefraction:
iproc = kPLightRefraction;
- case 20:
+ case kEMSCOlocaledep :
iproc = kPPhotoelectric;
break;
default:
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
//______________________________________________________________________________
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);
}
{
// 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);
}
{
// 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)
}
-
-#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)
- {
- //
- // 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);
- }
-
- void usersteppingckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
+ 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);
+ }
+
+ void ustckv(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->SetMreg(mreg,LTCLCM.mlatm1);
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 %s\n", nphot, mreg, x, y, z, fluka->CurrentVolName());
+ fluka->SetCaller(kUSTCKV);
+ if (fluka->GetVerbosityLevel() >= 3)
(TVirtualMCApplication::Instance())->Stepping();
+
}
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff