#include "TFluka.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
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;
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);
}
}
{
// 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;
+ return GENSTK.np + FHEAVY.npheav;
else if (caller == 50) {
// Cerenkov Photon production
return fNCerenkov;
Int_t caller = GetCaller();
if (caller == 6) { // valid only after usdraw
- if (FINUC.np > 0) {
+ 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);
{
// Return particle lifetime for particle with pdg code pdg.
Int_t ifluka = IdFromPDG(pdg);
- return (PAPROP.thalf[ifluka+6]);
+ return (PAPROP.tmnlf[ifluka+6]);
}
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
fluka->SetNCerenkov(nphot);
fluka->SetCaller(50);
if (fluka->GetVerbosityLevel() >= 3)
- printf("userstepping ckv: %10d %10d %13.3f %13.3f %13.2f %s\n", nphot, mreg, x, y, z, fluka->CurrentVolName());
(TVirtualMCApplication::Instance())->Stepping();
+
}
}
-
void TFlukaConfigOption::SetProcess(const char* flagname, Int_t flag)
{
- printf("SetProcess %s %5d %5d \n", flagname, flag, fMedium);
-
// Set a process flag
const TString process[15] =
{"DCAY", "PAIR", "COMP", "PHOT", "PFIS", "DRAY", "ANNI", "BREM", "MUNU", "CKOV",
Int_t i;
for (i = 0; i < 15; i++) {
if (process[i].CompareTo(flagname) == 0) {
- printf("flag %5d\n", i);
-
fProcessFlag[i] = flag;
if (fMedium == -1) fgDProcessFlag[i] = flag;
break;
// Write the FLUKA input cards for the set of process flags and cuts
//
//
- if (fMedium > -1) {
+ if (fMedium != -1) {
TFluka* fluka = (TFluka*) gMC;
+ fMedium = fgGeom->GetFlukaMaterial(fMedium);
+ //
+ // Check if material is actually used
+ Int_t* reglist;
+ Int_t nreg;
+ reglist = fgGeom->GetMaterialList(fMedium, nreg);
+ if (nreg == 0) {
+ printf("Material not used !\n");
+ return;
+ }
+
TObjArray *matList = fluka->GetFlukaMaterials();
Int_t nmaterial = matList->GetEntriesFast();
TGeoMaterial* material = 0;
void TFlukaConfigOption::ProcessPAIR()
{
// Process PAIR option
- fprintf(fgFile,"*\n* --- PAIR --- Pair production by gammas, muons and hadrons. Flag = %5d, PPCUTM = %13.4g \n",
- fProcessFlag[kPAIR], fCutValue[kPPCUTM]);
+ fprintf(fgFile,"*\n* --- PAIR --- Pair production by gammas, muons and hadrons. Flag = %5d, PPCUTM = %13.4g, PPCUTE = %13.4g\n",
+ fProcessFlag[kPAIR], fCutValue[kCUTELE], fCutValue[kPPCUTM]);
//
// gamma -> e+ e-
//
if (fProcessFlag[kPAIR] > 0) {
- fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0., fCMatMin, fCMatMax, 1.);
+ fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0.0, fCMatMin, fCMatMax, 1.);
} else {
fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 1e10, fCMatMin, fCMatMax, 1.);
}
//
if (fProcessFlag[kPHOT] > 0) {
- fprintf(fgFile,"EMFCUT %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0. , 0., 0., fCMatMin, fCMatMax, 1.);
+ fprintf(fgFile,"EMFCUT %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0., fCMatMin, fCMatMax, 1.);
} else {
fprintf(fgFile,"EMFCUT %10.1f%10.4g%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0., 1.e10, 0., fCMatMin, fCMatMax, 1.);
}
#include "TGeoManager.h"
#include "TGeoVolume.h"
#include "TObjString.h"
-#include "Fepisor.h"
+#include "Fsourcm.h"
#ifndef WIN32
# define idnrwr idnrwr_
ClassImp(TFlukaMCGeometry)
-TFlukaMCGeometry* TFlukaMCGeometry::fgInstance=0;
+TFlukaMCGeometry* TFlukaMCGeometry::fgInstance= NULL;
//_____________________________________________________________________________
TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title)
//
fDebug = kFALSE;
fLastMaterial = 0;
+ fCurrentRegion = 0;
+ fCurrentLattice = 0;
fNextRegion = 0;
fNextLattice = 0;
fRegionList = 0;
+ fIndmat = 0;
gFluka = (TFluka*)gMC;
gMCGeom = this;
gNstep = 0;
//
fDebug = kFALSE;
fLastMaterial = 0;
+ fCurrentRegion = 0;
+ fCurrentLattice = 0;
fNextRegion = 0;
fNextLattice = 0;
fRegionList = 0;
+ fIndmat = 0;
gFluka = (TFluka*)gMC;
gMCGeom = this;
gNstep = 0;
Int_t TFlukaMCGeometry::GetFlukaMaterial(Int_t imed) const
{
// Returns FLUKA material index for medium IMED
- TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1);
+ TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1);
if (!med) {
Error("GetFlukaMaterial", "MEDIUM %i nor found", imed);
return -1;
}
+ TGeoMaterial* mat = med->GetMaterial();
+ if (!mat->IsUsed()) return -1;
Int_t imatfl = med->GetMaterial()->GetIndex();
return imatfl;
}
}
TGeoMixture *mix = 0;
TGeoElement *element;
- TGeoElementTable *table = TGeoElementTable::Instance();
+ TGeoElementTable *table = gGeoManager->GetElementTable();
switch (ind) {
case 0: // AIR
- mix = new TGeoMixture("AIR", 4, 0.001205);
+ mix = new TGeoMixture("TPC_AIR", 4, 0.001205);
element = table->GetElement(6); // C
mix->DefineElement(0, element, 0.000124);
element = table->GetElement(7); // N
mix->DefineElement(3, element, 0.012827);
break;
case 1: //SDD SI CHIP
- mix = new TGeoMixture("SDD_SI", 6, 2.4485);
+ mix = new TGeoMixture("ITS_SDD_SI", 6, 2.4485);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.004367771);
element = table->GetElement(6);
mix->DefineElement(5, element, 0.09814344903);
break;
case 2: // WATER
- mix = new TGeoMixture("WATER", 2, 1.0);
+ mix = new TGeoMixture("ITS_WATER", 2, 1.0);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.111898344);
element = table->GetElement(8);
mix->DefineElement(1, element, 0.888101656);
break;
case 3: // CERAMICS
- mix = new TGeoMixture("CERAMICS", 5, 3.6);
+ mix = new TGeoMixture("ITS_CERAMICS", 5, 3.6);
element = table->GetElement(8);
mix->DefineElement(0, element, 0.59956);
element = table->GetElement(13);
mix->DefineElement(4, element, 0.0115);
break;
case 4: // EPOXY
- mix = new TGeoMixture("G10FR4", 4, 1.8);
+ mix = new TGeoMixture("MUON_G10FR4", 4, 1.8);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.19);
element = table->GetElement(6);
mix->DefineElement(3, element, 0.28);
break;
case 6: // KAPTON
- mix = new TGeoMixture("KAPTON", 4, 1.3);
+ mix = new TGeoMixture("ITS_KAPTON", 4, 1.3);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.026363415);
element = table->GetElement(6);
mix->DefineElement(3, element, 0.209238060);
break;
case 7: // INOX
- mix = new TGeoMixture("INOX", 9, 7.9);
+ mix = new TGeoMixture("ITS_INOX", 9, 7.9);
element = table->GetElement(6);
mix->DefineElement(0, element, 0.0003);
element = table->GetElement(14);
mix->DefineElement(3, element, 0.19137459);
break;
case 9: // SDD-C-AL
- mix = new TGeoMixture("SDD-C-AL", 5, 1.9837);
+ mix = new TGeoMixture("ITS_SDD-C-AL", 5, 1.9837);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.022632);
element = table->GetElement(6);
mix->DefineElement(4, element, 0.1);
break;
case 10: // X7R-CAP
- mix = new TGeoMixture("X7R-CAP", 7, 6.72);
+ mix = new TGeoMixture("ITS_X7R-CAP", 7, 6.72);
element = table->GetElement(8);
mix->DefineElement(0, element, 0.085975822);
element = table->GetElement(22);
mix->DefineElement(6, element, 0.2081768);
break;
case 11: // SDD ruby sph. Al2O3
- mix = new TGeoMixture("AL2O3", 2, 3.97);
+ mix = new TGeoMixture("ITS_AL2O3", 2, 3.97);
element = table->GetElement(8);
mix->DefineElement(0, element, 0.5293);
element = table->GetElement(13);
mix->DefineElement(1, element, 0.4707);
break;
case 12: // SDD HV microcable
- mix = new TGeoMixture("HV-CABLE", 5, 1.6087);
+ mix = new TGeoMixture("ITS_HV-CABLE", 5, 1.6087);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.01983871336);
element = table->GetElement(6);
mix->DefineElement(4, element, 0.247536);
break;
case 13: //SDD LV+signal cable
- mix = new TGeoMixture("LV-CABLE", 5, 2.1035);
+ mix = new TGeoMixture("ITS_LV-CABLE", 5, 2.1035);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.0082859922);
element = table->GetElement(6);
mix->DefineElement(4, element, 0.68572);
break;
case 14: //SDD hybrid microcab
- mix = new TGeoMixture("HYB-CAB", 5, 2.0502);
+ mix = new TGeoMixture("ITS_HYB-CAB", 5, 2.0502);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.00926228815);
element = table->GetElement(6);
mix->DefineElement(4, element, 0.64869);
break;
case 15: //SDD anode microcab
- mix = new TGeoMixture("ANOD-CAB", 5, 1.7854);
+ mix = new TGeoMixture("ITS_ANOD-CAB", 5, 1.7854);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.0128595919215);
element = table->GetElement(6);
mix->DefineElement(4, element, 0.431909);
break;
case 16: // inox/alum
- mix = new TGeoMixture("INOX-AL", 5, 3.0705);
+ mix = new TGeoMixture("ITS_INOX-AL", 5, 3.0705);
element = table->GetElement(13);
mix->DefineElement(0, element, 0.816164);
element = table->GetElement(14);
element = table->GetElement(28);
mix->DefineElement(4, element, 0.0183836);
case 17: // MYLAR
- mix = new TGeoMixture("MYLAR", 3, 1.39);
+ mix = new TGeoMixture("TPC_MYLAR", 3, 1.39);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.0416667);
element = table->GetElement(6);
mix->DefineElement(2, element, 0.333333);
break;
case 18: // SPDBUS(AL+KPT+EPOX) - unknown composition
- mix = new TGeoMixture("SPDBUS", 1, 1.906);
+ mix = new TGeoMixture("ITS_SPDBUS", 1, 1.906);
element = table->GetElement(9);
mix->DefineElement(0, element, 1.);
z = element->Z();
break;
case 19: // SDD/SSD rings - unknown composition
- mix = new TGeoMixture("SDDRINGS", 1, 1.8097);
+ mix = new TGeoMixture("ITS_SDDRINGS", 1, 1.8097);
element = table->GetElement(6);
mix->DefineElement(0, element, 1.);
z = element->Z();
break;
case 20: // SPD end ladder - unknown composition
- mix = new TGeoMixture("SPDEL", 1, 3.6374);
+ mix = new TGeoMixture("ITS_SPDEL", 1, 3.6374);
element = table->GetElement(22);
mix->DefineElement(0, element, 1.);
z = element->Z();
break;
case 21: // SDD end ladder - unknown composition
- mix = new TGeoMixture("SDDEL", 1, 0.3824);
+ mix = new TGeoMixture("ITS_SDDEL", 1, 0.3824);
element = table->GetElement(30);
mix->DefineElement(0, element, 1.);
z = element->Z();
break;
case 22: // SSD end ladder - unknown composition
- mix = new TGeoMixture("SSDEL", 1, 0.68);
+ mix = new TGeoMixture("ITS_SSDEL", 1, 0.68);
element = table->GetElement(16);
mix->DefineElement(0, element, 1.);
z = element->Z();
Int_t i,j,idmat;
Int_t counttothree, nelem;
Double_t a,z,rho, w;
- TGeoElementTable *table = TGeoElementTable::Instance();
+ TGeoElementTable *table = gGeoManager->GetElementTable();
TGeoElement *element;
element = table->GetElement(13);
element->SetTitle("ALUMINUM"); // this is how FLUKA likes it ...
// Adjust material names and add them to FLUKA list
for (i=0; i<nmater; i++) {
mat = (TGeoMaterial*)matlist->At(i);
- if (!mat->IsUsed()) continue;
+ if (!mat->IsUsed()) continue;
z = mat->GetZ();
a = mat->GetA();
rho = mat->GetDensity();
}
matname = mat->GetName();
FlukaMatName(matname);
-/*
- // material with one element: create it as mixture since it can be duplicated
- if (!mat->IsMixture()) {
- // normal materials
- mix = new TGeoMixture(matname.Data(), 1, rho);
- mix->DefineElement(0, mat->GetElement(), 1.);
- mat->SetIndex(nfmater+3);
- for (j=0; j<nmed; j++) {
- med = (TGeoMedium*)medlist->At(j);
- if (med->GetMaterial() == mat) {
- med->SetMaterial(mix);
- if (mat->GetCerenkovProperties()) {
- mix->SetCerenkovProperties(mat->GetCerenkovProperties());
- mat->SetCerenkovProperties(0);
- }
- break;
- }
- }
- mat = (TGeoMaterial*)mix;
- }
-*/
+
mat->SetIndex(nfmater+3);
objstr = new TObjString(matname.Data());
fMatList->Add(mat);
}
out.close();
fLastMaterial = nfmater+2;
-
- if (!gFluka->IsGeneratePemf()) {
- if (gSystem->AccessPathName("FlukaVmc.pemf")) Fatal("CreateFlukaMatFile", "No pemf file in working directory");
- return;
- }
}
void TFlukaMCGeometry::CreatePemfFile()
//_____________________________________________________________________________
void TFlukaMCGeometry::FlukaMatName(TString &str) const
{
+// Strip the detector name
+
+ TObjArray * tokens = str.Tokenize("_");
+ Int_t ntok = tokens->GetEntries();
+ if (ntok > 0) {
+ TString head = ((TObjString*) tokens->At(0))->GetString();
+ Int_t nhead = head.Length();
+ str = str.Remove(0, nhead + 1);
+ }
+ tokens->Clear();
+ delete tokens;
+
// Convert a name to upper case 8 chars.
ToFlukaString(str);
Int_t ilast;
}
+
// FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers
//_____________________________________________________________________________
printf(" oldReg=%i oldLttc=%i pstep=%f\n",oldReg, oldLttc, propStep);
}
Int_t olttc = oldLttc;
- if (oldLttc<0) {
+ if (oldLttc<=0) {
gGeoManager->FindNode(pSx,pSy,pSz);
olttc = gGeoManager->GetCurrentNodeId()+1;
+ if (gMCGeom->IsDebugging()) {
+ printf("WOOPS: old reg/latt = %i/%i\n",oldReg,oldLttc);
+ printf("point: (%16.12f, %16.12f, %16.12f) in lttc=%i\n", pSx,pSy,pSz,olttc);
+ }
}
+ Int_t ccreg,cclat;
+ gMCGeom->GetCurrentRegion(ccreg,cclat);
gMCGeom->SetCurrentRegion(oldReg, olttc);
// Initialize default return values
lttcFlag = 0;
Double_t extra = 1.E-10;
-// printf("ERROR: (%f, %f, %f)\n",pSx,pSy,pSz);
if (gMCGeom->IsDebugging()) printf(" current path: %s\n", gGeoManager->GetPath());
gGeoManager->SetCurrentPoint(pSx+extra*pV[0], pSy+extra*pV[1], pSz+extra*pV[2]);
gGeoManager->SetCurrentDirection(pV);
- gGeoManager->FindNextBoundary(-propStep);
+// gGeoManager->FindNextBoundary(-propStep);
+ gGeoManager->FindNextBoundary(-100000.);
Double_t snext = gGeoManager->GetStep();
if (snext<=0.0) {
- saf = 0.0;
- newReg = -3;
- sLt[lttcFlag] = 0.0;
- if (gMCGeom->IsDebugging()) printf("BACK SCATTERING\n");
- return;
+ snext = 0.0;
+ /*
+ if (!crossed) {
+ // artefact due to MS
+ saf = 0.0;
+ newReg = -3;
+ sLt[lttcFlag] = 0.0;
+ if (gMCGeom->IsDebugging()) {
+ printf("BACK SCATTERING (reg/lttc = %i/%i)\n", oldReg,oldLttc);
+ printf("point: (%16.11f, %16.11f, %16.11f)\n", pSx,pSy,pSz);
+ printf("dir : (%16.11f, %16.11f, %16.11f)\n", pV[0],pV[1],pV[2]);
+ }
+ return;
+ }
+ */
+ // else ...
+ // FLUKA native detects some extra errors even when a boundary was crossed, returning a -3
+ // If the particle turns back to the original region before crossing at the first step,
+ // we just return the distance to the boundary, not issuing an error (not due to MS)
}
-
snext += extra;
saf = gGeoManager->GetSafeDistance();
saf -= extra;
+ if (saf>snext) printf("ERROR: saf=%f .GT. snext=%f\n", saf,snext);
if (saf<0) saf=0.0;
else saf -= saf*3.0e-09;
+// saf *= 0.3;
+ PAREM.dist = snext;
NORLAT.distn = snext;
NORLAT.xn[0] += snext*pV[0];
NORLAT.xn[1] += snext*pV[1];
NORLAT.xn[2] += snext*pV[2];
-// saf = 0.0; // !!! TEMPORARY FOR TESTING MAGSPHF - TO BE REMOVED
if (snext>propStep) {
// Next boundary further than proposed step, which is approved
retStep = propStep;
return;
}
// The next boundary is closer. We try to cross it.
+ if (saf>propStep) saf = propStep; // Safety should be less than the proposed step if a boundary will be crossed
+// saf = 0.0; // !!! SAFETY SHOULD BE 0 IF THE BOUNDARY WILL BE CROSSED ???
gGeoManager->SetCurrentPoint(pSx,pSy,pSz);
Double_t *point = gGeoManager->GetCurrentPoint();
Double_t *dir = gGeoManager->GetCurrentDirection();
memcpy(pt, point, 3*sizeof(Double_t));
Int_t i;
- for (i=0;i<3;i++) point[i] += (snext+1E-6)*dir[i];
+ extra = 1.E-13;
+ for (i=0;i<3;i++) point[i] += (snext+extra)*dir[i];
// locate next region
gGeoManager->FindNode();
newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
+ while (newLttc==olttc) {
+ extra *= 10.;
+ if (extra>1.E-5) break;
+ for (i=0;i<3;i++) point[i] += extra*dir[i];
+ gGeoManager->FindNode();
+ newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
+ }
gGeoManager->SetCurrentPoint(pt);
newReg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
if (gMCGeom->IsDebugging()) printf(" newReg=%i newLttc=%i\n", newReg, newLttc);
+
// We really crossed the boundary, but is it the same region ?
gMCGeom->SetNextRegion(newReg, newLttc);
if (newReg==oldReg && newLttc!=olttc) {
printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
}
+/*
NORLAT.xn[0] = pSx;
NORLAT.xn[1] = pSy;
NORLAT.xn[2] = pSz;
NORLAT.wn[0] = pV[0];
NORLAT.wn[1] = pV[1];
NORLAT.wn[2] = pV[2];
+*/
TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
if (gGeoManager->IsOutside()) {
newReg = gMCGeom->NofVolumes()+1;
extern "C" {
Double_t abscff(Double_t& wvlngt, Double_t& /*omgpho*/, Int_t& mmat)
{
+// printf("abscff%f %d\n", wvlngt, mmat);
+
//
// Return absorption length for given photon energy and material
//
#include "Fdblprc.h" //(DBLPRC) fluka common
#include "Ftrackr.h" //(TRACKR) fluka common
#include "Fopphst.h" //(OPPHST) fluka common
-#include "Fstack.h" //(STACK) fluka common
+#include "Fflkstk.h" //(FLKSTK) fluka common
#ifndef WIN32
# define mgdraw mgdraw_
if (TRACKR.jtrack == -1) {
trackId = OPPHST.louopp[OPPHST.lstopp];
if (trackId == 0) {
- trackId = STACK.ispark[STACK.lstack][mkbmx2-1];
+ trackId = FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1];
}
} else {
trackId = TRACKR.ispusr[mkbmx2-1];
if (!TRACKR.ispusr[mkbmx2 - 2]) {
//
// Single step
- if (verbosityLevel >= 3) {
+ if (TRACKR.jtrack == -1 && trackId == 109340) {
cout << endl << " !!! I am in mgdraw - calling Stepping(): " << icode << endl;
cout << endl << " Track Id = " << trackId << " region = " << mreg << endl;
+ printf("Stepsize %13.5e \n", fluka->TrackStep());
}
+
+
+
(TVirtualMCApplication::Instance())->Stepping();
fluka->SetTrackIsNew(kFALSE);
// Fluka commons
#include "Fdblprc.h" //(DBLPRC) fluka common
#include "Fdimpar.h" //(DIMPAR) fluka parameters
-#include "Fepisor.h" //(EPISOR) fluka common
-#include "Fstack.h" //(STACK) fluka common
-#include "Fstars.h" //(STARS) fluka common
-#include "Fbeam.h" //(BEAM) fluka common
+#include "Fsourcm.h" //(EPISOR) fluka common
+#include "Fflkstk.h" //(FLKSTK) fluka common
+#include "Fsumcou.h" //(SUMCOU) fluka common
#include "Fpaprop.h" //(PAPROP) fluka common
#include "Fltclcm.h" //(LTCLCM) fluka common
#include "Fopphst.h" //(OPPHST) fluka common
-//#include "Fcaslim.h" //(CASLIM) fluka common
//Virutal MC
#include "TFluka.h"
Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
if (debug) {
cout << "==> source(" << nomore << ")" << endl;
- cout << "\t* EPISOR.lsouit = " << (EPISOR.lsouit?'T':'F') << endl;
+ cout << "\t* SOURCM.lsouit = " << (SOURCM.lsouit?'T':'F') << endl;
}
static Bool_t lfirst = true;
}
if (lfirst) {
- EPISOR.tkesum = zerzer;
+ SOURCM.tkesum = zerzer;
lfirst = false;
- EPISOR.lussrc = true;
+ SOURCM.lussrc = true;
} else {
//
// Post-track actions for primary track
if (itrack<0) {
nomore = 1;
- EPISOR.lsouit = false;
+ SOURCM.lsouit = false;
if (debug) {
- cout << "\t* EPISOR.lsouit = " << (EPISOR.lsouit?'T':'F') << endl;
+ cout << "\t* SOURCM.lsouit = " << (SOURCM.lsouit?'T':'F') << endl;
cout << "\t* No more particles. Exiting..." << endl;
cout << "<== source(" << nomore << ")" << endl;
}
printf("Event has been stopped by user !");
fluka->SetStopEvent(kFALSE);
nomore = 1;
- EPISOR.lsouit = false;
+ SOURCM.lsouit = false;
return;
}
<< particle->Energy() << " GeV "
<< particle->GetMass() << " GeV " << endl;
}
- /* Lstack is the stack counter: of course any time source is called it
+ /* Npflka is the stack counter: of course any time source is called it
* must be =0
*/
/* Cosines (tx,ty,tz)*/
Double_t cosz = TMath::Sqrt(oneone - cosx*cosx - cosy*cosy);
if (particle->Pz() < 0.) cosz = -cosz;
- //STACK.ilo[STACK.lstack] = BEAM.ijbeam;
if (pdg != 50000050 && pdg != 50000051) {
- STACK.lstack++;
+ FLKSTK.npflka++;
Int_t ifl = fluka-> IdFromPDG(pdg);
- STACK.ilo[STACK.lstack] = ifl;
- /* Wt is the weight of the particle*/
- STACK.wt[STACK.lstack] = oneone;
- STARS.weipri += STACK.wt[STACK.lstack];
+ FLKSTK.iloflk[FLKSTK.npflka] = ifl;
+ /* Wtflk is the weight of the particle*/
+ FLKSTK.wtflk[FLKSTK.npflka] = oneone;
+ SUMCOU.weipri += FLKSTK.wtflk[FLKSTK.npflka];
- STACK.lo[STACK.lstack] = 1;
+ FLKSTK.loflk[FLKSTK.npflka] = 1;
/* User dependent flag:*/
- STACK.louse[STACK.lstack] = 0;
+ FLKSTK.louse[FLKSTK.npflka] = 0;
/* User dependent spare variables:*/
Int_t ispr = 0;
for (ispr = 0; ispr < mkbmx1; ispr++)
- STACK.sparek[STACK.lstack][ispr] = zerzer;
+ FLKSTK.sparek[FLKSTK.npflka][ispr] = zerzer;
/* User dependent spare flags:*/
for (ispr = 0; ispr < mkbmx2; ispr++)
- STACK.ispark[STACK.lstack][ispr] = 0;
+ FLKSTK.ispark[FLKSTK.npflka][ispr] = 0;
/* Save the track number of the stack particle:*/
- STACK.ispark[STACK.lstack][mkbmx2-1] = itrack;
- STACK.nparma++;
- STACK.numpar[STACK.lstack] = STACK.nparma;
- STACK.nevent[STACK.lstack] = 0;
- STACK.dfnear[STACK.lstack] = +zerzer;
+ FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1] = itrack;
+ FLKSTK.nparma++;
+ FLKSTK.numpar[FLKSTK.npflka] = FLKSTK.nparma;
+ FLKSTK.nevent[FLKSTK.npflka] = 0;
+ FLKSTK.dfnear[FLKSTK.npflka] = +zerzer;
/* Particle age (s)*/
- STACK.agestk[STACK.lstack] = +zerzer;
- STACK.cmpath[STACK.lstack] = +zerzer;
- STACK.aknshr[STACK.lstack] = -twotwo;
+ FLKSTK.agestk[FLKSTK.npflka] = +zerzer;
+ FLKSTK.cmpath[FLKSTK.npflka] = +zerzer;
+ FLKSTK.aknshr[FLKSTK.npflka] = -twotwo;
/* Group number for "low" energy neutrons, set to 0 anyway*/
- STACK.igroup[STACK.lstack] = 0;
+ FLKSTK.igroup[FLKSTK.npflka] = 0;
/* Kinetic energy */
Double_t p = particle->P();
Double_t mass = PAPROP.am[ifl + 6];
- STACK.tke[STACK.lstack] = TMath::Sqrt( p * p + mass * mass) - mass;
+ FLKSTK.tkeflk[FLKSTK.npflka] = TMath::Sqrt( p * p + mass * mass) - mass;
/* Particle momentum*/
- STACK.pmom [STACK.lstack] = p;
+ FLKSTK.pmoflk [FLKSTK.npflka] = p;
- STACK.tx [STACK.lstack] = cosx;
- STACK.ty [STACK.lstack] = cosy;
- STACK.tz [STACK.lstack] = cosz;
+ FLKSTK.txflk [FLKSTK.npflka] = cosx;
+ FLKSTK.tyflk [FLKSTK.npflka] = cosy;
+ FLKSTK.tzflk [FLKSTK.npflka] = cosz;
/* Polarization cosines:*/
if (polarisation.Mag()) {
Double_t cospolx = polarisation.Px() / polarisation.Mag();
Double_t cospoly = polarisation.Py() / polarisation.Mag();
Double_t cospolz = sqrt(oneone - cospolx * cospolx - cospoly * cospoly);
- STACK.txpol [STACK.lstack] = cospolx;
- STACK.typol [STACK.lstack] = cospoly;
- STACK.tzpol [STACK.lstack] = cospolz;
+ FLKSTK.txpol [FLKSTK.npflka] = cospolx;
+ FLKSTK.typol [FLKSTK.npflka] = cospoly;
+ FLKSTK.tzpol [FLKSTK.npflka] = cospolz;
}
else {
- STACK.txpol [STACK.lstack] = -twotwo;
- STACK.typol [STACK.lstack] = +zerzer;
- STACK.tzpol [STACK.lstack] = +zerzer;
+ FLKSTK.txpol [FLKSTK.npflka] = -twotwo;
+ FLKSTK.typol [FLKSTK.npflka] = +zerzer;
+ FLKSTK.tzpol [FLKSTK.npflka] = +zerzer;
}
/* Particle coordinates*/
// Vertext coordinates;
- STACK.xa [STACK.lstack] = particle->Vx();
- STACK.ya [STACK.lstack] = particle->Vy();
- STACK.za [STACK.lstack] = particle->Vz();
+ FLKSTK.xflk [FLKSTK.npflka] = particle->Vx();
+ FLKSTK.yflk [FLKSTK.npflka] = particle->Vy();
+ FLKSTK.zflk [FLKSTK.npflka] = particle->Vz();
/* Calculate the total kinetic energy of the primaries: don't change*/
- Int_t st_ilo = STACK.ilo[STACK.lstack];
+ Int_t st_ilo = FLKSTK.iloflk[FLKSTK.npflka];
if ( st_ilo != 0 )
- EPISOR.tkesum +=
- ((STACK.tke[STACK.lstack] + PAPROP.amdisc[st_ilo+6])
- * STACK.wt[STACK.lstack]);
+ SOURCM.tkesum +=
+ ((FLKSTK.tkeflk[FLKSTK.npflka] + PAPROP.amdisc[st_ilo+6])
+ * FLKSTK.wtflk[FLKSTK.npflka]);
else
- EPISOR.tkesum += (STACK.tke[STACK.lstack] * STACK.wt[STACK.lstack]);
+ SOURCM.tkesum += (FLKSTK.tkeflk[FLKSTK.npflka] * FLKSTK.wtflk[FLKSTK.npflka]);
/* Here we ask for the region number of the hitting point.
- * NREG (LSTACK) = ...
+ * NRGFLK (LFLKSTK) = ...
* The following line makes the starting region search much more
* robust if particles are starting very close to a boundary:
*/
- geocrs( STACK.tx[STACK.lstack],
- STACK.ty[STACK.lstack],
- STACK.tz[STACK.lstack] );
+ geocrs( FLKSTK.txflk[FLKSTK.npflka],
+ FLKSTK.tyflk[FLKSTK.npflka],
+ FLKSTK.tzflk[FLKSTK.npflka] );
Int_t idisc;
- georeg ( STACK.xa[STACK.lstack],
- STACK.ya[STACK.lstack],
- STACK.za[STACK.lstack],
- STACK.nreg[STACK.lstack],
+ georeg ( FLKSTK.xflk[FLKSTK.npflka],
+ FLKSTK.yflk[FLKSTK.npflka],
+ FLKSTK.zflk[FLKSTK.npflka],
+ FLKSTK.nrgflk[FLKSTK.npflka],
idisc);//<-- dummy return variable not used
/* Do not change these cards:*/
Int_t igeohsm1 = 1;
Int_t igeohsm2 = -11;
- geohsm ( STACK.nhspnt[STACK.lstack], igeohsm1, igeohsm2, LTCLCM.mlattc );
- STACK.nlattc[STACK.lstack] = LTCLCM.mlattc;
+ geohsm ( FLKSTK.nhspnt[FLKSTK.npflka], igeohsm1, igeohsm2, LTCLCM.mlattc );
+ FLKSTK.nlattc[FLKSTK.npflka] = LTCLCM.mlattc;
soevsv();
} else {
//
// Ckeck transport cut first
Int_t ireg = EMFSTK.iremf[kp];
- Double_t cut = (TMath::Abs(EMFSTK.ichemf[kp]) == 1) ? EFMRGN.ecut[ireg-1] : EFMRGN.pcut[ireg-1];
+ Double_t cut = (TMath::Abs(EMFSTK.ichemf[kp]) == 1) ? EFMRGN.elethr[ireg-1] : EFMRGN.phothr[ireg-1];
Double_t e = EMFSTK.etemf[kp];
if ((e < cut)
&& (
(EMFSTK.ichemf[kp] == 0) ||
(EMFSTK.ichemf[kp] == -1) ||
- (EMFSTK.ichemf[kp] == 1 && EFMRGN.pcut[ireg-1] > emassmev)
+ (EMFSTK.ichemf[kp] == 1 && EFMRGN.phothr[ireg-1] > emassmev)
)
)
{
e *= emvgev;
Int_t pdg = fluka->PDGFromId(flukaid);
Double_t p = sqrt(e * e - PAPROP.am[flukaid+6] * PAPROP.am[flukaid+6]);
- Double_t px = p * EMFSTK.u[kp];
- Double_t pz = p * EMFSTK.v[kp];
- Double_t py = p * EMFSTK.w[kp];
+ Double_t px = p * EMFSTK.uemf[kp];
+ Double_t pz = p * EMFSTK.vemf[kp];
+ Double_t py = p * EMFSTK.wemf[kp];
Double_t tof = EMFSTK.agemf[kp];
Double_t polx = EMFSTK.upol[kp];
Double_t poly = EMFSTK.vpol[kp];
Double_t polz = EMFSTK.wpol[kp];
- Double_t vx = EMFSTK.x[kp];
- Double_t vy = EMFSTK.y[kp];
- Double_t vz = EMFSTK.z[kp];
+ Double_t vx = EMFSTK.xemf[kp];
+ Double_t vy = EMFSTK.yemf[kp];
+ Double_t vz = EMFSTK.zemf[kp];
Double_t weight = EMFSTK.wtemf[kp];
Int_t ntr;
#include "Fdblprc.h" //(DBLPRC) fluka common
#include "Fevtflg.h" //(EVTFLG) fluka common
#include "Fpaprop.h" //(PAPROP) fluka common
-#include "Fstack.h" //(STACK) fluka common
+#include "Fflkstk.h" //(FLKSTK) fluka common
#include "Ftrackr.h" //(TRACKR) fluka common
-#include "Ffinuc.h" //(FINUC) fluka common
+#include "Fgenstk.h" //(GENSTK) fluka common
//Virtual MC
#include "TFluka.h"
//* *
//*----------------------------------------------------------------------*
-// STACK.lstack = stack pointer
-// STACK.louse = user flag
+// FLKSTK.npflka = stack pointer
+// FLKSTK.louse = user flag
// TRACKR.llouse = user defined flag for the current particle
- STACK.louse[STACK.lstack] = TRACKR.llouse;
+ FLKSTK.louse[FLKSTK.npflka] = TRACKR.llouse;
// mkbmx1 = dimension for kwb real spare array in fluka stack in DIMPAR
// mkbmx2 = dimension for kwb int. spare array in fluka stack in DIMPAR
-// STACK.sparek = spare real variables available for k.w.burn
-// STACK.ispark = spare integer variables available for k.w.burn
+// FLKSTK.sparek = spare real variables available for k.w.burn
+// FLKSTK.ispark = spare integer variables available for k.w.burn
// TRACKR.spausr = user defined spare variables for the current particle
// TRACKR.ispusr = user defined spare flags for the current particle
Int_t ispr;
for (ispr = 0; ispr <= mkbmx1 - 1; ispr++) {
- STACK.sparek[STACK.lstack][ispr] = TRACKR.spausr[ispr];
+ FLKSTK.sparek[FLKSTK.npflka][ispr] = TRACKR.spausr[ispr];
}
for (ispr = 0; ispr <= mkbmx2 - 1; ispr++) {
- STACK.ispark[STACK.lstack][ispr] = TRACKR.ispusr[ispr];
+ FLKSTK.ispark[FLKSTK.npflka][ispr] = TRACKR.ispusr[ispr];
}
// Get the pointer to the VMC
Int_t done = 0;
Int_t parent = TRACKR.ispusr[mkbmx2-1];
- Int_t kpart = FINUC.kpart[numsec-1];
+ Int_t kpart = GENSTK.kpart[numsec-1];
if (kpart < -6) return;
Int_t pdg = fluka->PDGFromId(kpart);
- Double_t px = FINUC.plr[numsec-1] * FINUC.cxr[numsec-1];
- Double_t py = FINUC.plr[numsec-1] * FINUC.cyr[numsec-1];
- Double_t pz = FINUC.plr[numsec-1] * FINUC.czr[numsec-1];
- Double_t e = FINUC.tki[numsec-1] + PAPROP.am[FINUC.kpart[numsec-1]+6];
+ Double_t px = GENSTK.plr[numsec-1] * GENSTK.cxr[numsec-1];
+ Double_t py = GENSTK.plr[numsec-1] * GENSTK.cyr[numsec-1];
+ Double_t pz = GENSTK.plr[numsec-1] * GENSTK.czr[numsec-1];
+ Double_t e = GENSTK.tki[numsec-1] + PAPROP.am[GENSTK.kpart[numsec-1]+6];
Double_t vx = xx;
Double_t vy = yy;
Double_t vz = zz;
Double_t tof = TRACKR.atrack;
- Double_t polx = FINUC.cxrpol[numsec-1];
- Double_t poly = FINUC.cyrpol[numsec-1];
- Double_t polz = FINUC.czrpol[numsec-1];
+ Double_t polx = GENSTK.cxrpol[numsec-1];
+ Double_t poly = GENSTK.cyrpol[numsec-1];
+ Double_t polz = GENSTK.czrpol[numsec-1];
TMCProcess mech = kPHadronic;
}
- Double_t weight = FINUC.wei[numsec-1];
+ Double_t weight = GENSTK.wei[numsec-1];
Int_t is = 0;
Int_t ntr;
//
<< numsec << "npprmr " << npprmr << endl;
//
// Save current track number
- STACK.ispark[STACK.lstack][mkbmx2-1] = ntr;
- STACK.ispark[STACK.lstack][mkbmx2-2] = 0;
+ FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1] = ntr;
+ FLKSTK.ispark[FLKSTK.npflka][mkbmx2-2] = 0;
} // end of if (numsec-1 > npprmr)
} // end of stuprf
} // end of extern "C"
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff