#include "TList.h"
#include "TCallf77.h"
#include "TFluka.h"
+#include "TSystem.h"
#include "TFlukaMCGeometry.h"
#include "TFlukaConfigOption.h"
#include "TGeoManager.h"
// program load the right cross sections, and equal to the names included in
// the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his
// own .pemf, in order to get the right cross sections loaded in memory.
-
-
+ // Materials defined by FLUKA
TString sname;
gGeoManager->Export("flgeom.root");
if (fname) sname = fname;
element->SetTitle("ALUMINUM"); // this is how FLUKA likes it ...
element = table->GetElement(15);
element->SetTitle("PHOSPHO"); // same story ...
-// element = table->GetElement(10);
-// element->SetTitle("ARGON"); // NEON not in neutron xsec table
Int_t nelements = table->GetNelements();
TList *matlist = gGeoManager->GetListOfMaterials();
-// TList *medlist = gGeoManager->GetListOfMedia();
-// Int_t nmed = medlist->GetSize();
TIter next(matlist);
Int_t nmater = matlist->GetSize();
- Int_t nfmater = 0;
+ Int_t nfmater = 0;
+ Int_t newind = 26; // here non predefined materials start
+
TGeoMaterial *mat;
TGeoMixture *mix = 0;
TString matname;
rho = 0.999;
mat = new TGeoMaterial(matname, element->A(), element->Z(), rho);
- mat->SetIndex(nfmater+3);
+ // Check if the element has been predefined by FLUKA
+ Int_t pmid = GetPredefinedMaterialId(Int_t(element->Z()));
+ if (pmid > 0) {
+ mat->SetIndex(pmid);
+ } else {
+ mat->SetIndex(newind++);
+ }
+
mat->SetUsed(kTRUE);
fMatList->Add(mat);
objstr = new TObjString(matname.Data());
}
fIndmat = nfmater;
-// TGeoMedium *med;
// Adjust material names and add them to FLUKA list
for (i=0; i<nmater; i++) {
mat = (TGeoMaterial*)matlist->At(i);
matname = mat->GetName();
FlukaMatName(matname);
- mat->SetIndex(nfmater+3);
+ mat->SetIndex(newind++);
objstr = new TObjString(matname.Data());
fMatList->Add(mat);
fMatNames->Add(objstr);
Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
if (crtlttc == lttc) return;
gGeoManager->CdNode(lttc-1);
+ while (gGeoManager->GetCurrentVolume()->IsAssembly()) gGeoManager->CdUp();
}
//_____________________________________________________________________________
return gNstep;
}
-// FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers
+Int_t TFlukaMCGeometry::GetPredefinedMaterialId(Int_t z) const
+{
+// Get predifined material id from Z if present in list
+ const Int_t kMax = 25;
+
+ static Int_t idFluka[kMax] =
+ {-1, -1, 1, 2, 4, 6, 7, 8, 12, 13,
+ 26, 29, 47, 14, 79, 80, 82, 73, 11, 18,
+ 20, 50, 74, 22, 28};
+
+ Int_t id = -1;
+
+ for (Int_t i = 0; i < kMax; i++)
+ {
+ if (z == idFluka[i]) {
+ id = i + 1;
+ break;
+ }
+
+ }
+
+ return id;
+}
+
+// FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers
+//
//_____________________________________________________________________________
Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/)
{
{
// Initialize FLUKa point and direction;
+ static Int_t ierr = 0;
gNstep++;
-
-// if( (gNstep > 43912170 && gNstep < 43912196 ) ||
-// (gNstep > 47424560 && gNstep < 47424581 ) ||
-// (gNstep > 54388266 && gNstep < 54388319 )
-// ) gMCGeom->SetDebugMode();
-// else gMCGeom->SetDebugMode(kFALSE);
+// gMCGeom->SetDebugMode(kTRUE);
NORLAT.xn[0] = pSx;
NORLAT.xn[1] = pSy;
oldReg = gGeoManager->GetCurrentVolume()->GetNumber();
}
- if (gMCGeom->IsDebugging()) cout << "g1wr gNstep=" << gNstep
- << " oldReg="<< oldReg <<" olttc="<< olttc
- << " track=" << TRACKR.ispusr[mkbmx2-1] << endl;
-
- Int_t ccreg,cclat;
+ if (gMCGeom->IsDebugging()) {
+ cout << "g1wr gNstep=" << gNstep << " oldReg="<< oldReg <<" olttc="<< olttc
+ << " track=" << TRACKR.ispusr[mkbmx2-1] << endl;
+ cout << " point: (" << pSx << ", " << pSy << ", " << pSz << ") dir: ("
+ << pV[0] << ", " << pV[1] << ", " << pV[2] << ")" << endl;
+ }
+
+
+ Int_t ccreg=0,cclat=0;
gMCGeom->GetCurrentRegion(ccreg,cclat);
Bool_t crossed = (ccreg==oldReg && cclat==olttc)?kFALSE:kTRUE;
gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
gGeoManager->SetCurrentDirection(pV);
-
+ Double_t pt[3], local[3], ldir[3];
+ Int_t pid = TRACKR.jtrack;
+ pt[0] = pSx;
+ pt[1] = pSy;
+ pt[2] = pSz;
+ gGeoManager->MasterToLocal(pt,local);
+ gGeoManager->MasterToLocalVect(pV,ldir);
+/*
+ Bool_t valid = gGeoManager->GetCurrentVolume()->Contains(local);
+ if (!valid) {
+ printf("location not valid in %s pid=%i\n", gGeoManager->GetPath(),pid);
+ printf("local:(%f, %f, %f) ldir:(%f, %f, %f)\n", local[0],local[1],local[2],ldir[0],ldir[1],ldir[2]);
+// gGeoManager->FindNode();
+// printf(" -> actual location: %s\n", gGeoManager->GetPath());
+ }
+*/
+ Double_t pstep = propStep;
+ Double_t snext = propStep;
+ const Double_t epsil = 0.9999999 * TGeoShape::Tolerance();
+ // This should never happen !!!
+ if (pstep<TGeoShape::Tolerance()) {
+ printf("Proposed step is 0 !!!\n");
+ pstep = 2.*TGeoShape::Tolerance();
+ }
if (crossed) {
- gGeoManager->FindNextBoundaryAndStep(propStep);
- saf = 0.0;
+ gGeoManager->FindNextBoundaryAndStep(pstep);
+ snext = gGeoManager->GetStep();
+ saf = 0.;
+ if (snext <= TGeoShape::Tolerance()) {
+// printf("FLUKA crossed boundary but snext=%g\n", snext);
+ ierr++;
+ snext = epsil;
+ } else {
+ snext += TGeoShape::Tolerance();
+ ierr = 0;
+ }
} else {
- gGeoManager->FindNextBoundaryAndStep(propStep, kTRUE);
+ gGeoManager->FindNextBoundaryAndStep(pstep, kTRUE);
+ snext = gGeoManager->GetStep();
saf = gGeoManager->GetSafeDistance();
- if (saf<0) saf=0.0;
+ if (snext <= TGeoShape::Tolerance()) {
+// printf("FLUKA put particle on bondary without crossing\n");
+ ierr++;
+ snext = epsil;
+ saf = 0.;
+ } else {
+ snext += TGeoShape::Tolerance();
+ ierr = 0;
+ }
+ if (saf<0) saf=0.;
saf -= saf*3.0e-09;
}
-
- Double_t snext = gGeoManager->GetStep();
-
- if (snext<=0.0) snext = TGeoShape::Tolerance();
-
+// if (ierr>1) {
+// printf("%d snext=%g\n", ierr, snext);
+// }
+ if (ierr == 10) {
+// printf("Too many null steps - sending error code -33...\n");
+ newReg = -33; // Error code
+ ierr = 0;
+ return;
+ }
+
PAREM.dist = snext;
NORLAT.distn = snext;
NORLAT.xn[0] += snext*pV[0];
// We really crossed the boundary, but is it the same region ?
gMCGeom->SetNextRegion(newReg, newLttc);
- Int_t pid = TRACKR.jtrack;
if ( ((newReg==oldReg && newLttc!=olttc) || (oldReg!=newReg && olttc==newLttc) ) && pid!=-1) {
// Virtual boundary between replicants
newReg = gFluka->GetDummyRegion();
printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
}
flagErr = 0;
- TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
+ Double_t epsil = 1000.*TGeoShape::Tolerance();
+ TGeoNode *node = gGeoManager->FindNode(pSx+epsil*pV[0], pSy+epsil*pV[1], pSz+epsil*pV[2]);
if (gGeoManager->IsOutside()) {
newReg = gMCGeom->NofVolumes()+2;
newLttc = TFlukaMCGeometry::kLttcOutside;
if( oldReg!=newReg && oldLttc==newLttc ) {
// this should not happen!! ??? Ernesto
- cout << " lkwr oldReg!=newReg ("<< oldReg <<"!="<< newReg
- << ") && oldLttc==newLttc ("<< newLttc <<") !!!!!!!!!!!!!!!!!" << endl;
+// cout << " lkwr oldReg!=newReg ("<< oldReg <<"!="<< newReg
+// << ") && oldLttc==newLttc ("<< newLttc <<") !!!!!!!!!!!!!!!!!" << endl;
newReg = dummy;
newLttc = TFlukaMCGeometry::kLttcVirtual;
flagErr = newReg;
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff