// Author: Andrei Gheata 10/07/2003
#include "Riostream.h"
+#include "TList.h"
#include "TCallf77.h"
#include "TFluka.h"
+#include "TSystem.h"
#include "TFlukaMCGeometry.h"
#include "TFlukaConfigOption.h"
#include "TGeoManager.h"
#include "TObjString.h"
#include "Fsourcm.h"
#include "Ftrackr.h"
+#include "Fstepsz.h" //(STEPSZ) fluka common
#ifndef WIN32
# define idnrwr idnrwr_
void type_of_call g1wr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
Double_t * /*pV*/, Int_t & /*oldReg*/ , const Int_t & /*oldLttc*/, Double_t & /*propStep*/,
Int_t & /*nascFlag*/, Double_t & /*retStep*/, Int_t & /*newReg*/,
- Double_t & /*saf*/, Int_t & /*newLttc*/, Int_t & /*LttcFlag*/,
+ Double_t & /*saf*/, Int_t & /*newLttc*/, Int_t & /*LttcFlag*/,
Double_t *s /*Lt*/, Int_t * /*jrLt*/);
void type_of_call g1rtwr();
- void type_of_call conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/);
+ void type_of_call conhwr(Int_t & /*intHist*/, Int_t & /*incrCount*/);
void type_of_call inihwr(Int_t & /*intHist*/);
void type_of_call jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
Int_t & /*flukaReg*/);
Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
void type_of_call lkmgwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
- Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
+ Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
void type_of_call lkwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
- Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
+ Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
// void type_of_call magfld(const Double_t & /*pX*/, const Double_t & /*pY*/, const Double_t & /*pZ*/,
// Double_t & /*cosBx*/, Double_t & /*cosBy*/, Double_t & /*cosBz*/,
-// Double_t & /*Bmag*/, Int_t & /*reg*/, Int_t & /*idiscflag*/);
+// Double_t & /*Bmag*/, Int_t & /*reg*/, Int_t & /*idiscflag*/);
void type_of_call nrmlwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
Double_t & /*pVx*/, Double_t & /*pVy*/, Double_t & /*pVz*/,
- Double_t * /*norml*/, const Int_t & /*oldReg*/,
- const Int_t & /*newReg*/, Int_t & /*flagErr*/);
+ Double_t * /*norml*/, const Int_t & /*oldReg*/,
+ const Int_t & /*newReg*/, Int_t & /*flagErr*/);
void type_of_call rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/);
Int_t type_of_call isvhwr(const Int_t & /*fCheck*/, const Int_t & /*intHist*/);
//_____________________________________________________________________________
TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title)
- : TNamed(name, title)
+ :TNamed(name, title),
+ fDebug(kFALSE),
+ fLastMaterial(0),
+ fDummyRegion(0),
+ fCurrentRegion(0),
+ fCurrentLattice(0),
+ fNextRegion(0),
+ fNextLattice(0),
+ fRegionList(0),
+ fIndmat(0),
+ fMatList(new TObjArray(256)),
+ fMatNames(new TObjArray(256))
{
//
// Standard constructor
//
- fDebug = kFALSE;
- fLastMaterial = 0;
- fCurrentRegion = 0;
- fCurrentLattice = 0;
- fDummyRegion = 0;
- fNextRegion = 0;
- fNextLattice = 0;
- fRegionList = 0;
- fIndmat = 0;
gFluka = (TFluka*)gMC;
gMCGeom = this;
gNstep = 0;
- fMatList = new TObjArray(256);
- fMatNames = new TObjArray(256);
}
//_____________________________________________________________________________
TFlukaMCGeometry::TFlukaMCGeometry()
-{
+ :TNamed(),
+ fDebug(kFALSE),
+ fLastMaterial(0),
+ fDummyRegion(0),
+ fCurrentRegion(0),
+ fCurrentLattice(0),
+ fNextRegion(0),
+ fNextLattice(0),
+ fRegionList(0),
+ fIndmat(0),
+ fMatList(0),
+ fMatNames(0)
+
+{
//
// Default constructor
//
- fDebug = kFALSE;
- fLastMaterial = 0;
- fCurrentRegion = 0;
- fCurrentLattice = 0;
- fDummyRegion = 0;
- fNextRegion = 0;
- fNextLattice = 0;
- fRegionList = 0;
- fIndmat = 0;
gFluka = (TFluka*)gMC;
gMCGeom = this;
gNstep = 0;
- fMatList = 0;
- fMatNames = 0;
}
//_____________________________________________________________________________
//
// private methods
//
-//_____________________________________________________________________________
-TFlukaMCGeometry::TFlukaMCGeometry(const TFlukaMCGeometry &)
- : TNamed()
-{
- //
- // Copy constructor
- //
-}
//_____________________________________________________________________________
Double_t* TFlukaMCGeometry::CreateDoubleArray(Float_t* array, Int_t size) const
TGeoVolume *vol;
Int_t imedium, ireg;
while ((vol = (TGeoVolume*)next())) {
- imedium = vol->GetMedium()->GetId();
- if (imedium == imed) {
- ireg = vol->GetNumber();
- fRegionList[nreg++] = ireg;
- }
+ TGeoMedium* med;
+ if ((med = vol->GetMedium()) == 0) continue;
+ imedium = med->GetId();
+ if (imedium == imed) {
+ ireg = vol->GetNumber();
+ fRegionList[nreg++] = ireg;
+ }
}
return fRegionList;
-}
+}
//_____________________________________________________________________________
Int_t *TFlukaMCGeometry::GetMaterialList(Int_t imat, Int_t &nreg)
TGeoVolume *vol;
Int_t imaterial, ireg;
while ((vol = (TGeoVolume*)next())) {
- imaterial = vol->GetMedium()->GetMaterial()->GetIndex();
- if (imaterial == imat) {
- ireg = vol->GetNumber();
- fRegionList[nreg++] = ireg;
- }
+ TGeoMedium* med;
+ if ((med = vol->GetMedium()) == 0) continue;
+ imaterial = med->GetMaterial()->GetIndex();
+ if (imaterial == imat) {
+ ireg = vol->GetNumber();
+ fRegionList[nreg++] = ireg;
+ }
}
return fRegionList;
-}
+}
//_____________________________________________________________________________
Int_t TFlukaMCGeometry::NofVolumes() const
// 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);
PrintHeader(out, "TGEO MATERIAL ASSIGNMENTS");
for (i=1; i<=nvols; i++) {
-
+ TGeoMedium* med;
vol = gGeoManager->GetVolume(i);
- mat = vol->GetMedium()->GetMaterial();
+ if ((med = vol->GetMedium()) == 0) continue;
+ mat = med->GetMaterial();
idmat = mat->GetIndex();
for (Int_t j=0; j<nfmater; j++) {
mat = (TGeoMaterial*)fMatList->At(j);
TString sname;
for (i = fIndmat; i < fLastMaterial - 2; i++) {
- printf("Write Peg Files %d\n", i);
-
- mat = (TGeoMaterial*)fMatList->At(i);
- if (!mat->IsUsed()) continue;
- sname = "mat";
- sprintf(number, "%d", i);
- sname.Append(number);
- cout << endl;
- cout << endl;
- cout << "******************************************************************************" << endl;
- cout << "******************************************************************************" << endl;
- cout << endl;
- WritePegFile(i, &countNoStern, &countElemError, &countMixError, &countGas);
- sname.Prepend("$FLUPRO/pemf/rpemf peg/");
- gSystem->Exec(sname.Data());
-
- // check if the pemf file was created
- TString sname = Form("peg/mat%d.pemf", i);
- ifstream in( sname.Data() );
- if ( in ) {
- countMatOK++;
- in.close();
- } else {
- cout << "ERROR Fail to create the pemf file " << sname << endl;
- countPemfError++;
- }
+ printf("Write Peg Files %d\n", i);
+
+ mat = (TGeoMaterial*)fMatList->At(i);
+ if (!mat->IsUsed()) continue;
+ sname = "mat";
+ sprintf(number, "%d", i);
+ sname.Append(number);
+ cout << endl;
+ cout << endl;
+ cout << "******************************************************************************" << endl;
+ cout << "******************************************************************************" << endl;
+ cout << endl;
+ WritePegFile(i, &countNoStern, &countElemError, &countMixError, &countGas);
+ sname.Prepend("$FLUPRO/pemf/rpemf peg/");
+ gSystem->Exec(sname.Data());
+
+ // check if the pemf file was created
+ TString sname = Form("peg/mat%d.pemf", i);
+ ifstream in( sname.Data() );
+ if ( in ) {
+ countMatOK++;
+ in.close();
+ } else {
+ cout << "ERROR Fail to create the pemf file " << sname << endl;
+ countPemfError++;
+ }
}
cout << "Materials (pemf created) " << countMatOK << endl;
cout << "Not Sternheimer par. found " << countNoStern << endl;
while((proc = (TFlukaConfigOption*)next()))
{
if (proc->Medium() == mat->GetIndex()) {
- ap = proc->Cut(kCUTGAM);
- ae = proc->Cut(kCUTELE);
- if (ap == -1.) ap = TFlukaConfigOption::DefaultCut(kCUTGAM);
- if (ae == -1.) ae = TFlukaConfigOption::DefaultCut(kCUTELE);
- break;
+ ap = proc->Cut(kCUTGAM);
+ ae = proc->Cut(kCUTELE);
+ if (ap == -1.) ap = TFlukaConfigOption::DefaultCut(kCUTGAM);
+ if (ae == -1.) ae = TFlukaConfigOption::DefaultCut(kCUTELE);
+ break;
}
}
Double_t afact; // AFACT
Double_t sk; // SK
Double_t delta0; // DELTA0
+
+ sternheimerData():
+ longname(""), nelems(0), density(0), iev(0), cbar(0),
+ x0(0), x1(0), afact(0), sk(0), delta0(0) {}
};
TString shortname;
Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
if (crtlttc == lttc) return;
gGeoManager->CdNode(lttc-1);
+ while (gGeoManager->GetCurrentVolume()->IsAssembly()) gGeoManager->CdUp();
}
//_____________________________________________________________________________
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);
+ TString head = ((TObjString*) tokens->At(0))->GetString();
+ Int_t nhead = head.Length();
+ str = str.Remove(0, nhead + 1);
}
tokens->Clear();
delete tokens;
}
}
}
-
+
//______________________________________________________________________________
void TFlukaMCGeometry::Vname(const char *name, char *vname) const
{
}
+//______________________________________________________________________________
+Int_t TFlukaMCGeometry::GetNstep()
+{
+ // return gNstep for debug propose
+ 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*/)
{
void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
Double_t *pV, Int_t &oldReg , const Int_t &oldLttc, Double_t &propStep,
Int_t &/*nascFlag*/, Double_t &retStep, Int_t &newReg,
- Double_t &saf, Int_t &newLttc, Int_t <tcFlag,
+ Double_t &saf, Int_t &newLttc, Int_t <tcFlag,
Double_t *sLt, Int_t *jrLt)
{
// Initialize FLUKa point and direction;
+ static Int_t ierr = 0;
gNstep++;
+// gMCGeom->SetDebugMode(kTRUE);
+
NORLAT.xn[0] = pSx;
NORLAT.xn[1] = pSy;
NORLAT.xn[2] = pSz;
if (oldLttc<=0) {
gGeoManager->FindNode(pSx,pSy,pSz);
olttc = gGeoManager->GetCurrentNodeId()+1;
- }
- Int_t ccreg,cclat;
+ oldReg = gGeoManager->GetCurrentVolume()->GetNumber();
+ }
+
+ 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==oldLttc)?kFALSE:kTRUE;
+ Bool_t crossed = (ccreg==oldReg && cclat==olttc)?kFALSE:kTRUE;
+
gMCGeom->SetCurrentRegion(oldReg, olttc);
// Initialize default return values
lttcFlag = 0;
sLt[lttcFlag] = TGeoShape::Tolerance(); // null step in next region
jrLt[lttcFlag+1] = -1;
sLt[lttcFlag+1] = 0.; // null step in next region;
+ if (gMCGeom->IsDebugging()) printf("=> crossed dummy!! newReg=%i newLttc=%i\n", newReg, newLttc);
return;
- }
-
+ }
+
// Reset outside flag
gGeoManager->SetOutside(kFALSE);
-
+
curLttc = gGeoManager->GetCurrentNodeId()+1;
curReg = gGeoManager->GetCurrentVolume()->GetNumber();
if (olttc != curLttc) {
gGeoManager->CdNode(olttc-1);
curLttc = gGeoManager->GetCurrentNodeId()+1;
curReg = gGeoManager->GetCurrentVolume()->GetNumber();
- }
+ }
// Now the current TGeo state reflects the FLUKA state
-
+
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;
+ }
+// 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;
}
-
- Double_t snext = gGeoManager->GetStep();
-
- if (snext<=0.0) snext = TGeoShape::Tolerance();
-
+
PAREM.dist = snext;
NORLAT.distn = snext;
NORLAT.xn[0] += snext*pV[0];
NORLAT.xn[2] += snext*pV[2];
if (!gGeoManager->IsOnBoundary()) {
// Next boundary further than proposed step, which is approved
- if (saf>propStep) saf = propStep;
+ if (saf>propStep) saf = propStep;
retStep = propStep;
sLt[lttcFlag] = propStep;
return;
gGeoManager->SetCurrentPoint(pSx,pSy,pSz);
newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
newReg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+2):gGeoManager->GetCurrentVolume()->GetNumber();
- if (gMCGeom->IsDebugging()) printf(" newReg=%i newLttc=%i\n", newReg, newLttc);
+ 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);
- Int_t pid = TRACKR.jtrack;
- if (newReg==oldReg && newLttc!=olttc && pid!=-1) {
+ if ( ((newReg==oldReg && newLttc!=olttc) || (oldReg!=newReg && olttc==newLttc) ) && pid!=-1) {
// Virtual boundary between replicants
newReg = gFluka->GetDummyRegion();
newLttc = TFlukaMCGeometry::kLttcVirtual;
- }
-
+ if (gMCGeom->IsDebugging()) printf("=> virtual boundary!! newReg=%i newLttc=%i\n", newReg, newLttc);
+ }
+
retStep = snext;
sLt[lttcFlag] = snext;
lttcFlag++;
jrLt[lttcFlag] = newLttc;
sLt[lttcFlag] = snext;
jrLt[lttcFlag+1] = -1;
- sLt[lttcFlag+1] = 0.;
+
+ sLt[lttcFlag+1] = 0.;
gGeoManager->SetOutside(kFALSE);
gGeoManager->CdNode(olttc-1);
if (gMCGeom->IsDebugging()) {
printf("=> snext=%g safe=%g\n", snext, saf);
for (Int_t i=0; i<lttcFlag+1; i++) printf(" jrLt[%i]=%i sLt[%i]=%g\n", i,jrLt[i],i,sLt[i]);
- }
+ }
}
//_____________________________________________________________________________
void g1rtwr()
{
+
if (gMCGeom->IsDebugging()) printf("========== Dummy G1RTWR\n");
}
//_____________________________________________________________________________
-void conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/)
+void conhwr(Int_t & intHist, Int_t & incrCount)
{
- if (gMCGeom->IsDebugging()) printf("========== Dummy CONHWR\n");
+ if (gMCGeom->IsDebugging()) printf("========== Dummy CONHWR intHist=%d incrCount=%d currentNodeId=%d\n",
+ intHist, incrCount, gGeoManager->GetCurrentNodeId()+1 );
+// if( incrCount != -1 ) {
+// if (intHist==0) gGeoManager->CdTop();
+// else gGeoManager->CdNode(intHist-1);
+// }
+// intHist = gGeoManager->GetCurrentNodeId()+1;
}
//_____________________________________________________________________________
void inihwr(Int_t &intHist)
{
- if (gMCGeom->IsDebugging()) printf("========== Inside INIHWR -> reinitializing history: %i\n", intHist);
+ if (gMCGeom->IsDebugging())
+ printf("========== Inside INIHWR -> reinitializing history: %i \n", intHist);
if (gGeoManager->IsOutside()) gGeoManager->CdTop();
- if (intHist<=0) {
+ if (intHist<0) {
// printf("=== wrong history number\n");
return;
}
if (gMCGeom->IsDebugging()) {
printf(" --- current path: %s\n", gGeoManager->GetPath());
printf("<= INIHWR\n");
- }
+ }
}
//_____________________________________________________________________________
//_____________________________________________________________________________
void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
- Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
+ Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
{
if (gMCGeom->IsDebugging()) {
printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz);
//_____________________________________________________________________________
void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
- Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
+ Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
{
if (gMCGeom->IsDebugging()) {
printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz);
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 (newreg1==newReg && newlttc1==newLttc) {
newReg = dummy;
newLttc = TFlukaMCGeometry::kLttcVirtual;
- }
+ flagErr = newReg;
+ if (gMCGeom->IsDebugging()) printf(" virtual boundary (oldReg==dummy) !! newReg=%i newLttc=%i\n", newReg, newLttc);
+ }
return;
- }
+ }
if (oldReg==newReg && oldLttc!=newLttc) {
- newReg = gFluka->GetDummyRegion();
+ newReg = dummy;
newLttc = TFlukaMCGeometry::kLttcVirtual;
- }
-
+ if (gMCGeom->IsDebugging()) printf(" virtual boundary!! newReg=%i newLttc=%i\n", newReg, newLttc);
+ }
+
+ if( oldReg!=newReg && oldLttc==newLttc ) {
+ // this should not happen!! ??? Ernesto
+// cout << " lkwr oldReg!=newReg ("<< oldReg <<"!="<< newReg
+// << ") && oldLttc==newLttc ("<< newLttc <<") !!!!!!!!!!!!!!!!!" << endl;
+ newReg = dummy;
+ newLttc = TFlukaMCGeometry::kLttcVirtual;
+ flagErr = newReg;
+ }
+
if (gMCGeom->IsDebugging()) {
printf(" LKWR: newReg=%i newLttc=%i\n", newReg, newLttc);
}
//_____________________________________________________________________________
void nrmlwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
Double_t &pVx, Double_t &pVy, Double_t &pVz,
- Double_t *norml, const Int_t &oldReg,
- const Int_t &newReg, Int_t &flagErr)
+ Double_t *norml, const Int_t &oldReg,
+ const Int_t &newReg, Int_t &flagErr)
{
if (gMCGeom->IsDebugging()) {
printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz);
+ printf(" (%g, %g, %g)\n", NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]);
printf(" oldReg=%i, newReg=%i\n", oldReg,newReg);
}
gGeoManager->SetCurrentPoint(NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]);
norml[1] = -pVy;
norml[2] = -pVz;
} else {
- norml[0] = -dnorm[0];
- norml[1] = -dnorm[1];
- norml[2] = -dnorm[2];
+ norml[0] = -dnorm[0];
+ norml[1] = -dnorm[1];
+ norml[2] = -dnorm[2];
}
-
+
if (gMCGeom->IsDebugging()) {
printf(" normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]);
printf("<= NRMLWR\n");
- }
+ }
+
}
//_____________________________________________________________________________
// For TGeo, just return the current node ID. No copy need to be made.
- if (gMCGeom->IsDebugging()) printf("=> Inside ISVHWR\n");
+ if (gMCGeom->IsDebugging()) printf("=> Inside ISVHWR check=%d intHist=%d\n", check, intHist);
if (check<0) return intHist;
Int_t histInt = gGeoManager->GetCurrentNodeId()+1;
if (gMCGeom->IsDebugging()) printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath());
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff