// Author: Andrei Gheata 10/07/2003
#include "TObjString.h"
-#include "TFluka.h"
+#include "TFlukaGeo.h"
//#include "TVirtualMCApplication.h"
#include "TFlukaMCGeometry.h"
#include "TGeoManager.h"
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 & /*newReg*/, Int_t & /*flagErr*/, 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*/);
+// 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*/);
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*/,
// TFluka global pointer
TFluka *fluka = 0;
TFlukaMCGeometry *mcgeom = 0;
+Int_t kNstep = 0;
ClassImp(TFlukaMCGeometry)
//
// Standard constructor
//
+ fDebug = kFALSE;
+ fLastMaterial = 0;
+ fNextRegion = 0;
+ fNextLattice = 0;
+ fRegionList = 0;
fluka = (TFluka*)gMC;
mcgeom = this;
+ kNstep = 0;
}
//_____________________________________________________________________________
//
// Default constructor
//
+ fDebug = kFALSE;
+ fLastMaterial = 0;
+ fNextRegion = 0;
+ fNextLattice = 0;
+ fRegionList = 0;
fluka = (TFluka*)gMC;
mcgeom = this;
+ kNstep = 0;
}
//_____________________________________________________________________________
// Destructor
//
fgInstance=0;
+ if (fRegionList) delete [] fRegionList;
if (gGeoManager) delete gGeoManager;
}
Float_t &dens, Float_t &radl, Float_t &absl,
Float_t* /*ubuf*/, Int_t& /*nbuf*/)
{
- printf("Gfmate %i\n", imat);
+ if (fDebug) printf("Gfmate %i\n", imat);
TGeoMaterial *mat;
TIter next (gGeoManager->GetListOfMaterials());
while ((mat = (TGeoMaterial*)next())) {
dens = mat->GetDensity();
radl = mat->GetRadLen();
absl = mat->GetIntLen();
- printf(" ->material found : %s a=%g, z=%g, dens=%g, radl=%g, absl=%g\n", name, a,z,dens,radl,absl);
+ if (fDebug) printf(" ->material found : %s a=%g, z=%g, dens=%g, radl=%g, absl=%g\n", name, a,z,dens,radl,absl);
}
//_____________________________________________________________________________
Double_t &dens, Double_t &radl, Double_t &absl,
Double_t* /*ubuf*/, Int_t& /*nbuf*/)
{
- printf("Gfmate %i\n", imat);
- TGeoMaterial *mat;
+ if (fDebug) printf("Gfmate %i\n", imat);
+ TGeoMaterial *mat;
TIter next (gGeoManager->GetListOfMaterials());
while ((mat = (TGeoMaterial*)next())) {
if (mat->GetUniqueID() == (UInt_t)imat) break;
dens = mat->GetDensity();
radl = mat->GetRadLen();
absl = mat->GetIntLen();
- printf(" ->material found : %s a=%g, z=%g, dens=%g, radl=%g, absl=%g\n", name, a,z,dens,radl,absl);
+ if (fDebug) printf(" ->material found : %s a=%g, z=%g, dens=%g, radl=%g, absl=%g\n", name, a,z,dens,radl,absl);
}
//_____________________________________________________________________________
kmat = gGeoManager->GetListOfMaterials()->GetSize();
gGeoManager->Material(name, a, z, dens, kmat, radl, absl);
- printf("Material %s: kmat=%i, a=%g, z=%g, dens=%g\n", name, kmat, a, z, dens);
+ if (fDebug) printf("Material %s: kmat=%i, a=%g, z=%g, dens=%g\n", name, kmat, a, z, dens);
}
//_____________________________________________________________________________
}
}
kmat = gGeoManager->GetListOfMaterials()->GetSize();
- printf("Mixture %s with %i elem: kmat=%i, dens=%g\n", name, nlmat, kmat, dens);
- for (Int_t j=0; j<nlmat; j++) printf(" Elem %i: z=%g a=%g w=%g\n",j,z[j],a[j],wmat[j]);
+ if (fDebug) {
+ printf("Mixture %s with %i elem: kmat=%i, dens=%g\n", name, nlmat, kmat, dens);
+ for (Int_t j=0; j<nlmat; j++) printf(" Elem %i: z=%g a=%g w=%g\n",j,z[j],a[j],wmat[j]);
+ }
gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
}
//_____________________________________________________________________________
TGeoNode *node = gGeoManager->GetCurrentNode();
if (!node) imed = gGeoManager->GetTopNode()->GetVolume()->GetMedium()->GetId();
else imed = node->GetVolume()->GetMedium()->GetId();
- printf("GetMedium=%i\n", imed);
+ if (fDebug) printf("GetMedium=%i\n", imed);
return imed;
}
+//_____________________________________________________________________________
+Int_t TFlukaMCGeometry::GetFlukaMaterial(Int_t imed) const
+{
+// Returns FLUKA material index for medium IMED
+ TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1);
+ if (!med) {
+ Error("GetFlukaMaterial", "MEDIUM %i nor found", imed);
+ return -1;
+ }
+ Int_t imatfl = med->GetMaterial()->GetIndex();
+ return imatfl;
+}
+
+//_____________________________________________________________________________
+Int_t *TFlukaMCGeometry::GetRegionList(Int_t imed, Int_t &nreg)
+{
+// Get an ordered list of regions matching a given medium number
+ nreg = 0;
+ if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
+ TIter next(gGeoManager->GetListOfUVolumes());
+ TGeoVolume *vol;
+ Int_t imedium, ireg;
+ while ((vol = (TGeoVolume*)next())) {
+ imedium = vol->GetMedium()->GetId();
+ if (imedium == imed) {
+ ireg = vol->GetNumber();
+ fRegionList[nreg++] = ireg;
+ }
+ }
+ return fRegionList;
+}
+
+//_____________________________________________________________________________
+Int_t *TFlukaMCGeometry::GetMaterialList(Int_t imat, Int_t &nreg)
+{
+// Get an ordered list of regions matching a given medium number
+ nreg = 0;
+ if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1];
+ TIter next(gGeoManager->GetListOfUVolumes());
+ TGeoVolume *vol;
+ Int_t imaterial, ireg;
+ while ((vol = (TGeoVolume*)next())) {
+ imaterial = vol->GetMedium()->GetMaterial()->GetIndex();
+ if (imaterial == imat) {
+ ireg = vol->GetNumber();
+ fRegionList[nreg++] = ireg;
+ }
+ }
+ return fRegionList;
+}
//_____________________________________________________________________________
void TFlukaMCGeometry::Medium(Int_t& kmed, const char* name, Int_t nmat, Int_t isvol,
Int_t ifield, Double_t fieldm, Double_t tmaxfd,
// performed with g3helix; ifield = 3 if tracking performed with g3helx3.
//
- kmed = gGeoManager->GetListOfMedia()->GetSize()+3; // !!! in FLUKA they start with 3
+ kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
gGeoManager->Medium(name,kmed,nmat, isvol, ifield, fieldm, tmaxfd, stemax,deemax, epsil, stmin);
- printf("Medium %s: kmed=%i, nmat=%i, isvol=%i\n", name, kmed, nmat,isvol);
+ if (fDebug) printf("Medium %s: kmed=%i, nmat=%i, isvol=%i\n", name, kmed, nmat,isvol);
}
//_____________________________________________________________________________
krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
gGeoManager->Matrix(krot, thex, phix, they, phiy, thez, phiz);
- printf("Rotation %i defined\n", krot);
+ if (fDebug) printf("Rotation %i defined\n", krot);
}
//_____________________________________________________________________________
Vname(shape,vshape);
TGeoVolume* vol = gGeoManager->Volume(vname, shape, nmed, upar, npar);
- printf("Volume %s: id=%i shape=%s, nmed=%i\n", vname, vol->GetNumber(), shape, nmed);
+ if (fDebug) printf("Volume %s: id=%i shape=%s, nmed=%i\n", vname, vol->GetNumber(), shape, nmed);
return vol->GetNumber();
}
Vname(mother,vmother);
gGeoManager->Division(vname, vmother, iaxis, ndiv, 0, 0, 0, "n");
- printf("Division %s: mother=%s iaxis=%i ndiv=%i\n", vname, vmother, iaxis, ndiv);
+ if (fDebug) printf("Division %s: mother=%s iaxis=%i ndiv=%i\n", vname, vmother, iaxis, ndiv);
}
//_____________________________________________________________________________
Double_t *upar=0;
gGeoManager->Node(vname, nr, vmother, x, y, z, irot, isOnly, upar);
- printf("Adding daughter %s to %s: cpy=%i irot=%i only=%s\n", vname,vmother,nr,irot,only.Data());
+ if (fDebug) printf("Adding daughter %s to %s: cpy=%i irot=%i only=%s\n", vname,vmother,nr,irot,only.Data());
}
//_____________________________________________________________________________
Vname(mother,vmother);
gGeoManager->Node(vname,nr,vmother, x,y,z,irot,isOnly,upar,np);
- printf("Adding daughter(s) %s to %s: cpy=%i irot=%i only=%s\n", vname,vmother,nr,irot,only.Data());
+ if (fDebug) printf("Adding daughter(s) %s to %s: cpy=%i irot=%i only=%s\n", vname,vmother,nr,irot,only.Data());
}
//_____________________________________________________________________________
printf("VolId: Volume %s not found\n",name);
return 0;
}
- printf("VolId for %s: %i\n", name, uid);
+ if (fDebug) printf("VolId for %s: %i\n", name, uid);
return uid;
}
Error("VolName","volume with id=%d does not exist",id);
return "NULL";
}
- printf("VolName for id=%i: %s\n", id, volume->GetName());
+ if (fDebug) printf("VolName for id=%i: %s\n", id, volume->GetName());
return volume->GetName();
}
}
TGeoMedium *med = volume->GetMedium();
if (!med) return 0;
- printf("VolId2Mate id=%i: idmed=%i\n", id, med->GetId());
+ if (fDebug) printf("VolId2Mate id=%i: idmed=%i\n", id, med->GetId());
return med->GetId();
}
TGeoNode *node = gGeoManager->GetCurrentNode();
copyNo = node->GetNumber();
Int_t id = node->GetVolume()->GetNumber();
- printf("CurrentVolId(cpy=%i) = %i\n", copyNo, id);
+ if (fDebug) printf("CurrentVolId(cpy=%i) = %i\n", copyNo, id);
return id;
}
TGeoNode *node = gGeoManager->GetMother(off);
if (!node) return 0;
copyNo = node->GetNumber();
- printf("CurrentVolOffId(off=%i,cpy=%i) = %i\n", off,copyNo,node->GetVolume()->GetNumber() );
+ if (fDebug) printf("CurrentVolOffId(off=%i,cpy=%i) = %i\n", off,copyNo,node->GetVolume()->GetNumber() );
return node->GetVolume()->GetNumber();
}
// FLUKA specific
// Returns the current volume name
//
if (gGeoManager->IsOutside()) return 0;
- printf("CurrentVolName : %s\n", gGeoManager->GetCurrentVolume()->GetName());
+ if (fDebug) printf("CurrentVolName : %s\n", gGeoManager->GetCurrentVolume()->GetName());
return gGeoManager->GetCurrentVolume()->GetName();
}
//_____________________________________________________________________________
if (off==0) return CurrentVolName();
TGeoNode *node = gGeoManager->GetMother(off);
if (!node) return 0;
- printf("CurrentVolOffName(off=%i) : %s\n", off,node->GetVolume()->GetName());
+ if (fDebug) printf("CurrentVolOffName(off=%i) : %s\n", off,node->GetVolume()->GetName());
return node->GetVolume()->GetName();
}
}
//_____________________________________________________________________________
-void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname) const
+void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname)
{
// ==== from FLUGG ====
// NAMES OF ELEMENTS AND COMPOUNDS: the names must be written in upper case,
Double_t zel, ael, wel, rho;
char elname[8] = {' ',' ','_', 'E','L','E','M','\0'};
char digit[3];
+ Bool_t found = kFALSE;
- printf("Creating materials and compounds\n");
+ if (fDebug) printf("Creating materials and compounds\n");
for (i=0; i<nmater; i++) {
mat = (TGeoMaterial*)matlist->At(i);
- printf("material: %s index=%i\n", mat->GetName(), mat->GetIndex());
+ if (mat->GetZ()<1E-1) {
+ mat->SetIndex(2); // vacuum, built-in inside FLUKA
+ continue;
+ }
+// printf("material: %s index=%i: Z=%f A=%f rho=%f\n", mat->GetName(), mat->GetIndex(),mat->GetZ(),mat->GetA(),mat->GetDensity());
matorig = gGeoManager->FindDuplicateMaterial(mat);
if (matorig) {
idmat = matorig->GetIndex();
mat->SetIndex(idmat);
- printf(" -> found a duplicate: %s with index %i\n", matorig->GetName(), idmat);
+// printf(" -> found a duplicate: %s with index %i\n", matorig->GetName(), idmat);
matorig = 0;
} else {
- printf(" Adding to temp list with index %i\n", nfmater+3);
+// printf(" Adding to temp list with index %i\n", nfmater+3);
listfluka->Add(mat);
mat->SetIndex(nfmater+3);
matorig = mat;
}
}
}
- printf(" newmat name: %s\n", matname.Data());
+// printf(" newmat name: %s\n", matname.Data());
}
// now we have unique materials with unique names in the lists
-
if (matorig && matorig->IsMixture()) {
// create dummy materials for elements
rho = 0.999;
mix = (TGeoMixture*)matorig;
nelem = mix->GetNelements();
- printf(" material is a MIXTURE with %i elements:\n", nelem);
+// printf(" material is a MIXTURE with %i elements:\n", nelem);
for (j=0; j<nelem; j++) {
+ found = kFALSE;
zel = (mix->GetZmixt())[j];
- printf(" Zelem[%i] = %g\n",j,zel);
+ ael = (mix->GetAmixt())[j];
+// printf(" Zelem[%i] = %g\n",j,zel);
if ((zel-Int_t(zel))>0.01) {
- Warning("CreateFlukaMatFile", "element with Z=%f\n", zel);
+ TGeoMaterial *mat1;
+ for (Int_t imat=0; imat<nfmater; imat++) {
+ mat1 = (TGeoMaterial*)listfluka->At(imat);
+ if (TMath::Abs(mat1->GetZ()-zel)>1E-4) continue;
+ if (TMath::Abs(mat1->GetA()-ael)>1E-4) continue;
+ found = kTRUE;
+ break;
+ }
+ if (!found) Warning("CreateFlukaMatFile", "element with Z=%f\n", zel);
}
- if (!zelem[Int_t(zel)]) {
+ if (!zelem[Int_t(zel)] && !found) {
// write fluka element
memcpy(elname, &elNames[2*Int_t(zel-1)], 2);
zelem[Int_t(zel)] = 1;
- ael = (mix->GetAmixt())[j];
mat = new TGeoMaterial(elname, ael, zel, rho);
mat->SetIndex(nfmater+3);
- printf(" element not in list: new material %s at index=%i, Z=%g, A=%g, dummyrho=%g\n",
- elname,nfmater+3,zel,ael,rho);
+// printf(" element not in list: new material %s at index=%i, Z=%g, A=%g, dummyrho=%g\n",
+// elname,nfmater+3,zel,ael,rho);
listfluka->Add(mat);
objstr = new TObjString(elname);
listflukanames->Add(objstr);
}
}
// now dump materials in the file
- printf("DUMPING %i materials\n", nfmater);
+// printf("DUMPING %i materials\n", nfmater);
for (i=0; i<nfmater; i++) {
mat = (TGeoMaterial*)listfluka->At(i);
out << setw(10) << "MATERIAL ";
out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
- matname = mat->GetName();
+// matname = mat->GetName();
+ objstr = (TObjString*)listflukanames->At(i);
+ matname = objstr->GetString();
ToFlukaString(matname);
zmat = mat->GetZ();
+ if (zmat-Int_t(zmat)>0.01) {
+ if (zmat-Int_t(zmat)>0.5) zmat = Int_t(zmat)+1.;
+ else zmat = Int_t(zmat);
+ }
amat = mat->GetA();
rhomat = mat->GetDensity();
// write material card
out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(i+3);
out << setw(10) << " ";
- if (mat->IsMixture())
- out << setw(10) << mix->GetNelements();
- else
+ out << setw(10) << " ";
+ out << setw(8) << matname.Data() << endl;
+ // add LOW-MAT crd
+ if (!mat->IsMixture()) {
+ out << setw(10) << "LOW-MAT ";
+ out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
+ out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(i+3);
out << setw(10) << " ";
- out << setw(10) << matname.Data() << endl;
+ out << setw(10) << " ";
+ out << setw(10) << " ";
+ out << setw(10) << " ";
+ out << setw(10) << " ";
+ out << setw(8) << matname.Data() << endl;
+ }
}
// write mixture header
PrintHeader(out, "COMPOUNDS");
nelem = mix->GetNelements();
objstr = (TObjString*)listflukanames->At(i);
matname = objstr->GetString();
- printf("MIXTURE %s with index %i having %i elements\n", matname.Data(), mat->GetIndex(),nelem);
+// printf("MIXTURE %s with index %i having %i elements\n", matname.Data(), mat->GetIndex(),nelem);
for (j=0; j<nelem; j++) {
// dump mixture cards
- printf(" #elem %i: Z=%g, A=%g, W=%g\n", j, (mix->GetZmixt())[j],
- (mix->GetAmixt())[j],(mix->GetWmixt())[j]);
+// printf(" #elem %i: Z=%g, A=%g, W=%g\n", j, (mix->GetZmixt())[j],
+// (mix->GetAmixt())[j],(mix->GetWmixt())[j]);
wel = (mix->GetWmixt())[j];
zel = (mix->GetZmixt())[j];
- memcpy(elname, &elNames[2*Int_t(zel-1)], 2);
- element = (TGeoMaterial*)listfluka->FindObject(elname);
+ ael = (mix->GetAmixt())[j];
+ if (zel-Int_t(zel)>0.01) {
+ // loop the temporary list
+ element = 0;
+ TGeoMaterial *mat1;
+ for (Int_t imat=0; imat<i; imat++) {
+ mat1 = (TGeoMaterial*)listfluka->At(imat);
+ if (TMath::Abs(mat1->GetZ()-zel)>1E-4) continue;
+ if (TMath::Abs(mat1->GetA()-ael)>1E-4) continue;
+ element = mat1;
+ break;
+ }
+ } else {
+ memcpy(elname, &elNames[2*Int_t(zel-1)], 2);
+ element = (TGeoMaterial*)listfluka->FindObject(elname);
+ }
if (!element) {
Error("CreateFlukaMatFile", "Element Z=%g %s not found", zel, elname);
return;
}
idmat = element->GetIndex();
- printf("element %s , index=%i\n", element->GetName(), idmat);
+// printf("element %s , index=%i\n", element->GetName(), idmat);
out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
out << setw(10) << setiosflags(ios::fixed) << setprecision(6) << -wel;
out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
}
*/
// Now print the material assignments
- Double_t flagfield;
+ Double_t flagfield = 0.;
+ printf("#############################################################\n");
+ if (fluka->IsFieldEnabled()) {
+ flagfield = 1.;
+ printf("Magnetic field enabled\n");
+ } else printf("Magnetic field disabled\n");
+ printf("#############################################################\n");
+
PrintHeader(out, "TGEO MATERIAL ASSIGNMENTS");
for (i=1; i<=nvols; i++) {
vol = gGeoManager->GetVolume(i);
mat = vol->GetMedium()->GetMaterial();
idmat = mat->GetIndex();
- flagfield = (vol->GetField())?1.:0.;
out << setw(10) << "ASSIGNMAT ";
out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
out << setw(10) << setiosflags(ios::fixed) << Double_t(idmat);
out << setw(10) << setiosflags(ios::fixed) << Double_t(i);
out << setw(10) << "0.0";
+ out << setw(10) << "0.0";
out << setw(10) << setiosflags(ios::fixed) << flagfield;
+ out << setw(10) << "0.0";
out << endl;
}
delete listfluka;
listflukanames->Delete();
delete listflukanames;
out.close();
+ fLastMaterial = nfmater+2;
}
//_____________________________________________________________________________
if (gGeoManager->IsOutside()) return 0;
return gGeoManager->GetCurrentNode()->GetUniqueID();
}
+//_____________________________________________________________________________
+void TFlukaMCGeometry::SetMreg(Int_t mreg)
+{
+// Update if needed next history;
+ if (fluka->GetDummyBoundary()==2) {
+ gGeoManager->CdNode(fNextLattice-1);
+ return;
+ }
+ Int_t curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
+ if (mreg==curreg) return;
+ if (mreg==fNextRegion) {
+ if (fNextLattice!=999999999) gGeoManager->CdNode(fNextLattice-1);
+ return;
+ } else {
+ if (mreg == fCurrentRegion) {
+ if (fCurrentLattice!=999999999) gGeoManager->CdNode(fCurrentLattice-1);
+ return;
+ }
+ }
+ if (fDebug) printf("ERROR: mreg=%i neither current nor next region\n", mreg);
+}
+
+//_____________________________________________________________________________
+void TFlukaMCGeometry::SetCurrentRegion(Int_t mreg, Int_t latt)
+{
+// Set index/history for next entered region
+ fCurrentRegion = mreg;
+ fCurrentLattice = latt;
+}
+
+//_____________________________________________________________________________
+void TFlukaMCGeometry::SetNextRegion(Int_t mreg, Int_t latt)
+{
+// Set index/history for next entered region
+ fNextRegion = mreg;
+ fNextLattice = latt;
+}
//_____________________________________________________________________________
void TFlukaMCGeometry::ToFlukaString(TString &str) const
// FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers
-//_____________________________________________________________________________
-void jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
- Int_t &flukaReg)
-{
-// Geometry initialization wrapper called by FLUKAM. Provides to FLUKA the
-// number of regions (volumes in TGeo)
- // build application geometry
- printf("=> Inside JOMIWR\n");
- flukaReg = gGeoManager->GetListOfUVolumes()->GetEntries()+1;
-}
-
//_____________________________________________________________________________
Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/)
{
// card in fluka input), returns 1 if user wants Fluka always to
// use DNEAR (in this case, be sure that GEANT4 DNEAR is unique,
// coming from all directions!!!)
+ if (mcgeom->IsDebugging()) printf("========== Dummy IDNRWR\n");
return 0;
}
+//_____________________________________________________________________________
+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 *sLt, Int_t *jrLt)
+
+{
+ // Initialize FLUKa point and direction;
+ kNstep++;
+/*
+ if (kNstep>0) {
+ mcgeom->SetDebugMode(kTRUE);
+ fluka->SetVerbosityLevel(3);
+ }
+ if (kNstep>6520) {
+ mcgeom->SetDebugMode(kFALSE);
+ fluka->SetVerbosityLevel(0);
+ }
+ if ((kNstep%10)==0) printf("step %i\n", kNstep);
+*/
+
+ if (mcgeom->IsDebugging()) {
+ printf("========== Inside G1WR\n");
+ printf(" point/dir:(%14.9f, %14.9f, %14.9f, %g, %g, %g)\n", pSx,pSy,pSz,pV[0],pV[1],pV[2]);
+ printf(" oldReg=%i oldLttc=%i pstep=%f\n",oldReg, oldLttc, propStep);
+ }
+ gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
+ gGeoManager->SetCurrentDirection(pV);
+ mcgeom->SetCurrentRegion(oldReg, oldLttc);
+ // Initialize default return values
+ lttcFlag = 0;
+ jrLt[lttcFlag] = oldLttc;
+ sLt[lttcFlag] = propStep;
+ jrLt[lttcFlag+1] = -1;
+ sLt[lttcFlag+1] = 0.;
+ newReg = oldReg;
+ newLttc = oldLttc;
+ // check if dummy boundary flag is set
+ Int_t curLttc, curReg;
+ if (fluka->IsDummyBoundary()) {
+ // printf("Dummy boundary intercepted. Point is: %f, %f, %f\n", pSx, pSy, pSz);
+ Bool_t crossedDummy = (oldLttc == TFlukaMCGeometry::kLttcVirtual)?kTRUE:kFALSE;
+ if (crossedDummy) {
+ // FLUKA crossed the dummy boundary - update new region/history
+ retStep = 0.;
+ saf = 0.;
+ mcgeom->GetNextRegion(newReg, newLttc);
+ mcgeom->SetMreg(newReg);
+ if (mcgeom->IsDebugging()) printf(" virtual newReg=%i newLttc=%i\n", newReg, newLttc);
+ sLt[lttcFlag] = 0.; // null step in current region
+ lttcFlag++;
+ jrLt[lttcFlag] = newLttc;
+ sLt[lttcFlag] = 0.; // null step in next region
+ jrLt[lttcFlag+1] = -1;
+ sLt[lttcFlag+1] = 0.;
+ fluka->SetDummyBoundary(0);
+ return;
+ }
+ }
+
+ // Reset outside flag
+ if (gGeoManager->IsOutside()) {
+ gGeoManager->SetOutside(kFALSE);
+ gGeoManager->CdTop();
+ }
+
+ // Reset dummy boundary flag
+ fluka->SetDummyBoundary(0);
+
+ curLttc = gGeoManager->GetCurrentNodeId()+1;
+ curReg = gGeoManager->GetCurrentVolume()->GetNumber();
+ if (oldLttc != curLttc) {
+ // FLUKA crossed the boundary : we trust that the given point is really there,
+ // so we just update TGeo state
+ gGeoManager->CdNode(oldLttc-1);
+ curLttc = gGeoManager->GetCurrentNodeId()+1;
+ curReg = gGeoManager->GetCurrentVolume()->GetNumber();
+ if (mcgeom->IsDebugging()) printf(" re-initialized point: curReg=%i curLttc=%i\n", curReg, curLttc);
+ }
+ // Now the current TGeo state reflects the FLUKA state
+ if (mcgeom->IsDebugging()) printf(" current path: %s\n", gGeoManager->GetPath());
+ Double_t extra = 1E-6;
+ Double_t tmpStep = propStep + extra;
+ gGeoManager->FindNextBoundary(-tmpStep);
+ Double_t snext = gGeoManager->GetStep();
+ // !!!!!
+ if (snext<=0) {
+ // FLUKA is in the wrong region, notify it
+ if (mcgeom->IsDebugging()) printf("ERROR: snext=%f\n", snext);
+// newReg = -3;
+// return;
+ snext = extra;
+ }
+ saf = gGeoManager->GetSafeDistance();
+ Bool_t cross = kFALSE;
+ Bool_t onBound = kFALSE;
+ if (snext<tmpStep) {
+ // We have some boundary in the way
+ Double_t dd = snext-propStep;
+ if (dd < 0) {
+ cross = kTRUE;
+ dd = -dd;
+ }
+ if (dd < 1E-8) onBound = kTRUE;
+ }
+ snext += 1.E-8;
+ if (mcgeom->IsDebugging()) {
+ if (!cross) printf(" physical step approved: %f\n", propStep);
+ else printf(" boundary crossing at: %f\n", snext);
+ if (onBound) printf(" step on boundary limit ! NASC=%i\n", nascFlag);
+ }
+ if (!cross) {
+ // Next boundary further than proposed step, which is approved
+ retStep = propStep;
+ sLt[lttcFlag] = propStep;
+ return;
+ }
+ // The next boundary is closer. We try to cross it.
+ Double_t *point = gGeoManager->GetCurrentPoint();
+ Double_t *dir = gGeoManager->GetCurrentDirection();
+ Double_t pt[3];
+ memcpy(pt, point, 3*sizeof(Double_t));
+
+ Int_t i;
+ for (i=0;i<3;i++) point[i] += snext*dir[i];
+ gGeoManager->FindNode();
+ newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
+ if (newLttc == oldLttc) {
+ // brute force ...
+ // Just try a fast extra small step
+ snext += 1E-6;
+ for (i=0;i<3;i++) point[i] = pt[i]+snext*dir[i];
+ gGeoManager->FindNode();
+ newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
+ if (newLttc == oldLttc) {
+ // check if state changes at the end of the proposed step
+ for (i=0;i<3;i++) point[i] = pt[i]+propStep*dir[i];
+ gGeoManager->FindNode();
+ newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
+ if (newLttc==oldLttc) {
+ // approve step
+ retStep = propStep;
+ sLt[lttcFlag] = propStep;
+ return;
+ }
+ // snext is underestimated - we will create a virtual one to overcome the error
+// printf("some boundary in the way...\n");
+ }
+ }
+ gGeoManager->SetCurrentPoint(pt);
+// newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
+ newReg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
+ if (mcgeom->IsDebugging()) printf(" newReg=%i newLttc=%i\n", newReg, newLttc);
+
+ // We really crossed the boundary, but is it the same region ?
+ mcgeom->SetNextRegion(newReg, newLttc);
+ if (newReg == oldReg) {
+ // Virtual boundary between replicants
+ if (mcgeom->IsDebugging()) printf(" DUMMY boundary\n");
+ newReg = 1; // cheat FLUKA telling it it crossed the TOP region
+ newLttc = TFlukaMCGeometry::kLttcVirtual;
+ // mark that next boundary is virtual
+ fluka->SetDummyBoundary(1);
+ }
+ retStep = snext;
+ sLt[lttcFlag] = snext;
+ lttcFlag++;
+ jrLt[lttcFlag] = newLttc;
+ sLt[lttcFlag] = snext;
+ jrLt[lttcFlag+1] = -1;
+ sLt[lttcFlag+1] = 0.;
+
+ if (newLttc!=oldLttc) {
+ if (gGeoManager->IsOutside()) {
+ gGeoManager->SetOutside(kFALSE);
+ gGeoManager->CdTop();
+ }
+ gGeoManager->CdTop();
+ if (!gGeoManager->GetCurrentMatrix()->IsIdentity()) printf("ERROR at step %i\n", kNstep);
+ gGeoManager->CdNode(oldLttc-1);
+ }
+ if (mcgeom->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]);
+ }
+ if (mcgeom->IsDebugging()) printf("<= G1WR (in: %s)\n", gGeoManager->GetPath());
+}
+
+//_____________________________________________________________________________
+void g1rtwr()
+{
+ if (mcgeom->IsDebugging()) printf("========== Dummy G1RTWR\n");
+}
+
+//_____________________________________________________________________________
+void conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/)
+{
+ if (mcgeom->IsDebugging()) printf("========== Dummy CONHWR\n");
+}
+
+//_____________________________________________________________________________
+void inihwr(Int_t &intHist)
+{
+ if (mcgeom->IsDebugging()) printf("========== Inside INIHWR -> reinitializing history: %i\n", intHist);
+ if (gGeoManager->IsOutside()) gGeoManager->CdTop();
+ if (intHist<=0) {
+// printf("=== wrong history number\n");
+ return;
+ }
+ if (intHist==0) gGeoManager->CdTop();
+ else gGeoManager->CdNode(intHist-1);
+ if (mcgeom->IsDebugging()) {
+ printf(" --- current path: %s\n", gGeoManager->GetPath());
+ printf("<= INIHWR\n");
+ }
+}
+
+//_____________________________________________________________________________
+void jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
+ Int_t &flukaReg)
+{
+// Geometry initialization wrapper called by FLUKAM. Provides to FLUKA the
+// number of regions (volumes in TGeo)
+ // build application geometry
+ if (mcgeom->IsDebugging()) printf("========== Inside JOMIWR\n");
+ flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast();
+ if (mcgeom->IsDebugging()) printf("<= JOMIWR: last region=%i\n", flukaReg);
+}
+
+//_____________________________________________________________________________
+void lkdbwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
+ Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
+ Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
+{
+ if (mcgeom->IsDebugging()) {
+ printf("========== Inside LKDBWR (%f, %f, %f)\n",pSx, pSy, pSz);
+// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
+ printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
+ }
+ TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
+ if (gGeoManager->IsOutside()) {
+ newReg = mcgeom->NofVolumes()+1;
+// newLttc = gGeoManager->GetCurrentNodeId();
+ newLttc = 999999999;
+ if (mcgeom->IsDebugging()) {
+ printf("OUTSIDE\n");
+ printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
+ printf("<= LKMGWR\n");
+ }
+ flagErr = newReg;
+ return;
+ }
+ newReg = node->GetVolume()->GetNumber();
+ newLttc = gGeoManager->GetCurrentNodeId()+1;
+ mcgeom->SetNextRegion(newReg, newLttc);
+ flagErr = newReg;
+ if (mcgeom->IsDebugging()) {
+ printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
+ printf("<= LKDBWR\n");
+ }
+}
+
+//_____________________________________________________________________________
+void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
+ Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
+ Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
+{
+ if (mcgeom->IsDebugging()) {
+ printf("========== Inside LKFXWR (%f, %f, %f)\n",pSx, pSy, pSz);
+// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
+ printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
+ }
+ TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
+ if (gGeoManager->IsOutside()) {
+ newReg = mcgeom->NofVolumes()+1;
+// newLttc = gGeoManager->GetCurrentNodeId();
+ newLttc = 999999999;
+ if (mcgeom->IsDebugging()) {
+ printf("OUTSIDE\n");
+ printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
+ printf("<= LKMGWR\n");
+ }
+ flagErr = newReg;
+ return;
+ }
+ newReg = node->GetVolume()->GetNumber();
+ newLttc = gGeoManager->GetCurrentNodeId()+1;
+ mcgeom->SetNextRegion(newReg, newLttc);
+ flagErr = newReg;
+ if (mcgeom->IsDebugging()) {
+ printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
+ printf("<= LKFXWR\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)
+{
+ if (mcgeom->IsDebugging()) {
+ printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz);
+// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
+ printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
+ }
+ TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
+ if (gGeoManager->IsOutside()) {
+ newReg = mcgeom->NofVolumes()+1;
+// newLttc = gGeoManager->GetCurrentNodeId();
+ newLttc = 999999999;
+ if (mcgeom->IsDebugging()) {
+ printf("OUTSIDE\n");
+ printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
+ printf("<= LKMGWR\n");
+ }
+ flagErr = newReg;
+ return;
+ }
+ newReg = node->GetVolume()->GetNumber();
+ newLttc = gGeoManager->GetCurrentNodeId()+1;
+ mcgeom->SetNextRegion(newReg, newLttc);
+ flagErr = newReg;
+ if (mcgeom->IsDebugging()) {
+ printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
+ printf("<= LKMGWR\n");
+ }
+}
+
+//_____________________________________________________________________________
+void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
+ Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
+ Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
+{
+ if (mcgeom->IsDebugging()) {
+ printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz);
+// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
+ printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
+ }
+ TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
+ if (gGeoManager->IsOutside()) {
+ newReg = mcgeom->NofVolumes()+1;
+// newLttc = gGeoManager->GetCurrentNodeId();
+ newLttc = 999999999;
+ if (mcgeom->IsDebugging()) {
+ printf("OUTSIDE\n");
+ printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
+ printf("<= LKMGWR\n");
+ }
+ flagErr = newReg;
+ return;
+ }
+ newReg = node->GetVolume()->GetNumber();
+ newLttc = gGeoManager->GetCurrentNodeId()+1;
+ mcgeom->SetNextRegion(newReg, newLttc);
+ flagErr = newReg;
+ if (mcgeom->IsDebugging()) {
+ printf(" out: newReg=%i newLttc=%i in %s\n", newReg, newLttc, gGeoManager->GetPath());
+ printf("<= LKWR\n");
+ }
+}
+
+//_____________________________________________________________________________
+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)
+{
+ if (mcgeom->IsDebugging()) {
+ printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz);
+ printf(" oldReg=%i, newReg=%i\n", oldReg,newReg);
+ }
+// Int_t curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
+// Int_t curLttc = gGeoManager->GetCurrentNodeId()+1;
+// if (mcgeom->IsDebugging()) printf(" curReg=%i, curLttc=%i in: %s\n", curreg, curLttc, gGeoManager->GetPath());
+// Bool_t regsame = (curreg==oldReg)?kTRUE:kFALSE;
+ gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
+ gGeoManager->SetCurrentDirection(pVx,pVy,pVz);
+/*
+ if (!regsame) {
+ if (mcgeom->IsDebugging()) printf(" REGIONS DOEN NOT MATCH\n");
+ gGeoManager->FindNode();
+ curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
+ curLttc = gGeoManager->GetCurrentNodeId()+1;
+ if (mcgeom->IsDebugging()) printf(" re-initialized point: curReg=%i curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath());
+ }
+*/
+ Double_t *dnorm = gGeoManager->FindNormalFast();
+ flagErr = 0;
+ if (!dnorm) {
+ printf(" ERROR: Cannot compute fast normal\n");
+ flagErr = 1;
+ norml[0] = -pVx;
+ norml[1] = -pVy;
+ norml[2] = -pVz;
+ }
+ norml[0] = -dnorm[0];
+ norml[1] = -dnorm[1];
+ norml[2] = -dnorm[2];
+ if (mcgeom->IsDebugging()) printf(" normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]);
+// curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
+// curLttc = gGeoManager->GetCurrentNodeId()+1;
+ if (mcgeom->IsDebugging()) {
+// printf(" final location: curReg=%i, curLttc=%i in %s\n", curreg,curLttc,gGeoManager->GetPath());
+ printf("<= NRMLWR\n");
+ }
+}
+
+//_____________________________________________________________________________
+void rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
+ Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/)
+{
+ if (mcgeom->IsDebugging()) printf("=> Dummy RGRPWR\n");
+}
+
//_____________________________________________________________________________
Int_t isvhwr(const Int_t &check, const Int_t & intHist)
{
// For TGeo, just return the current node ID. No copy need to be made.
+ if (mcgeom->IsDebugging()) printf("=> Inside ISVHWR\n");
if (check<0) return intHist;
- Int_t histInt = gGeoManager->GetCurrentNodeId();
+ Int_t histInt = gGeoManager->GetCurrentNodeId()+1;
+ if (mcgeom->IsDebugging()) printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath());
return histInt;
}
-//_____________________________________________________________________________
-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 *sLt, Int_t *jrLt)
-{
-// from FLUGG:
-// Wrapper for geometry tracking: returns approved step of
-// particle and all variables that fluka G1 computes.
-
- // Initialize current point/direction
- gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
- gGeoManager->SetCurrentDirection(pV);
-
- // Initialize old path (FLUKA lattice history)
- if (oldLttc != jrLt[lttcFlag])
- printf("Woops: old history not matching jrLt[%i]. Checking other histories.\n",lttcFlag);
-
- gGeoManager->CdNode(oldLttc);
- TGeoVolume *oldvol = (gGeoManager->IsOutside())?0:gGeoManager->GetCurrentVolume();
- Int_t oldregion = (oldvol)?(mcgeom->NofVolumes()+1):oldvol->GetNumber(); // should it be 0?
- if (oldregion != oldReg) {
- while (lttcFlag>=0) {
- gGeoManager->CdNode(jrLt[lttcFlag]);
- oldvol = (gGeoManager->IsOutside())?0:gGeoManager->GetCurrentVolume();
- oldregion = (oldvol)?(mcgeom->NofVolumes()+1):oldvol->GetNumber();
- if (oldregion == oldReg) break;
- // bad history -> clean up jrLt[lttcFlag], sLt[lttcFlag]
- sLt[lttcFlag] = 0.;
- jrLt[lttcFlag] = -1;
- lttcFlag--;
- }
-
- if (oldregion != oldReg) {
- printf("Error: g1wr: history not found\n");
- printf(" relocating current point (%f, %f, %f)\n", pSx, pSy, pSz);
- gGeoManager->FindNode();
- oldvol = (gGeoManager->IsOutside())?0:gGeoManager->GetCurrentVolume();
- oldregion = (oldvol)?(mcgeom->NofVolumes()+1):oldvol->GetNumber();
- lttcFlag = 0;
- jrLt[lttcFlag]=isvhwr(0,0);
- }
- }
- sLt[lttcFlag] = 0.;
- // now 'oldregion' contains the real region, matching or not the old history
-
- // Compute geometry step/safety within physical step limit
- newReg = oldregion;
- Double_t steptot = 0.;
- Double_t snext = 0.;
- Int_t istep = 0;
- Bool_t done = kFALSE;
- while (!done) {
- gGeoManager->FindNextBoundary(propStep-steptot);
- snext = gGeoManager->GetStep();
- if (steptot == 0) saf = gGeoManager->GetSafeDistance();
- if (snext<propStep) {
- // There is a boundary on the way.
- // Make a step=snext+1E-6 to force boundary crossing
- steptot += snext;
- sLt[lttcFlag] = snext;
- retStep = snext;
- gGeoManager->Step();
- // Hopefully we end-up in a new region, else we do few small steps
- if (!gGeoManager->IsEntering()) {
-// sameregion = kTRUE;
- gGeoManager->SetStep(0.);
- istep = 0;
- }
- while (!gGeoManager->IsEntering() && steptot<propStep) {
- gGeoManager->Step();
- sLt[lttcFlag] += 1E-6;
- retStep = sLt[lttcFlag];
- steptot += 1E-6;
- istep++;
- if (istep>1000) {
- // we already made 10 extra microns and nothing
- printf("Woops: g1wr: extra 10 microns and no boundary...\n");
- gGeoManager->SetStep(propStep-steptot-1E-6);
- gGeoManager->Step();
- if (gGeoManager->IsEntering()) {
- retStep = sLt[lttcFlag];
- lttcFlag++;
- sLt[lttcFlag] = propStep-steptot;
- newReg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
- } else {
- sLt[lttcFlag] += propStep-steptot;
- }
- return;
- }
- }
- if (steptot>propStep) return;
- // we managed to cross the boundary -> in which region
- newReg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
- lttcFlag++;
- if (gGeoManager->IsOutside()) return;
-
- }
- }
-}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff