//
// Standard constructor
//
+ fDebug = kFALSE;
fLastMaterial = 0;
+ fNextRegion = 0;
+ fNextLattice = 0;
+ fRegionList = 0;
fluka = (TFluka*)gMC;
mcgeom = this;
}
//
// Default constructor
//
+ fDebug = kFALSE;
fLastMaterial = 0;
+ fNextRegion = 0;
+ fNextLattice = 0;
+ fRegionList = 0;
fluka = (TFluka*)gMC;
mcgeom = this;
}
// 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();
}
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);
if (mat->GetZ()<1E-1) {
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) {
vol = gGeoManager->GetVolume(i);
mat = vol->GetMedium()->GetMaterial();
idmat = mat->GetIndex();
- flagfield = (vol->GetField())?1.:0.;
+// flagfield = (vol->GetField())?1.:0.;
+ flagfield = 1.;
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;
if (gGeoManager->IsOutside()) return 0;
return gGeoManager->GetCurrentNode()->GetUniqueID();
}
+//_____________________________________________________________________________
+void TFlukaMCGeometry::SetMreg(Int_t mreg)
+{
+// Update if needed next history;
+ 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;
+ }
+ if (fDebug) printf("ERROR: mreg=%i neither current nor next region\n", mreg);
+}
+
+//_____________________________________________________________________________
+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
// 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!!!)
- printf("=> Inside IDNRWR\n");
+ 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 &saf, Int_t &newLttc, Int_t <tcFlag,
Double_t *sLt, Int_t *jrLt)
{
// from FLUGG:
// particle and all variables that fluka G1 computes.
// Initialize current point/direction
- printf("=> Inside G1WR\n");
+ 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]);
+ }
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]-1);
- 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);
- }
- }
+ if (mcgeom->IsDebugging()) printf(" oldReg=%i oldLttc=%i pstep=%f\n",oldReg, oldLttc, propStep);
+ if (oldLttc==999999999) printf("WOOPS - wrong old lattice\n");
+ if (gGeoManager->IsOutside()) {
+ gGeoManager->SetOutside(kFALSE);
+ gGeoManager->CdTop();
+ }
+ Int_t curLttc = gGeoManager->GetCurrentNodeId()+1;
+ Int_t curreg = gGeoManager->GetCurrentVolume()->GetNumber();
+ if (mcgeom->IsDebugging()) printf(" curReg=%i curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath());
+ Bool_t regsame = (curreg==oldReg)?kTRUE:kFALSE;
+ if (!regsame && mcgeom->IsDebugging()) printf(" REGIONS DOES NOT MATCH\n");
+ if (oldLttc != curLttc) {
+ if (mcgeom->IsDebugging()) printf(" HISTORIES DOES NOT MATCH\n");
+ gGeoManager->CdNode(oldLttc-1);
+ curLttc = gGeoManager->GetCurrentNodeId()+1;
+ curreg = gGeoManager->GetCurrentVolume()->GetNumber();
+ if (mcgeom->IsDebugging()) printf(" re-initialized point: curReg=%i curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath());
+ }
+ lttcFlag = 0;
sLt[lttcFlag] = 0.;
+ jrLt[lttcFlag] = curLttc;
// now 'oldregion' contains the real region, matching or not the old history
// Compute geometry step/safety within physical step limit
- newReg = oldregion;
+// newReg = oldregion;
+ Double_t *point = gGeoManager->GetCurrentPoint();
+ Double_t *dir = gGeoManager->GetCurrentDirection();
Double_t steptot = 0.;
Double_t snext = 0.;
Int_t istep = 0;
Bool_t done = kFALSE;
+ Double_t pst;
+ Int_t i;
while (!done) {
- gGeoManager->FindNextBoundary(propStep-steptot);
+ gGeoManager->FindNextBoundary(-propStep);
snext = gGeoManager->GetStep();
- if (steptot == 0) saf = gGeoManager->GetSafeDistance();
+ if (mcgeom->IsDebugging()) printf(" FindNextBoundary(%g) snext=%g\n", propStep, snext);
+ if (steptot == 0) {
+ saf = gGeoManager->GetSafeDistance();
+ if (mcgeom->IsDebugging()) printf(" Safety: %g\n", saf);
+ }
+ sLt[lttcFlag] = propStep;
+ jrLt[lttcFlag] = gGeoManager->GetCurrentNodeId()+1;
+ lttcFlag++; //1
+ sLt[lttcFlag] = 0.;
+ jrLt[lttcFlag] = -1;
+ newReg = curreg;
+ newLttc = oldLttc;
if (snext<propStep) {
// There is a boundary on the way.
// Make a step=snext+1E-6 to force boundary crossing
+ lttcFlag--; // 0
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();
+// lttcFlag++;
+ // make the step to get into the next region
+ for (i=0;i<3;i++) point[i]+=(snext+1E-6)*dir[i];
+ gGeoManager->FindNode();
+ istep = 0;
+ if (mcgeom->IsDebugging()) printf(" boundary: step made %g\n", snext);
+ while (gGeoManager->IsSameLocation() && steptot<propStep) {
+ if (istep>1E3) {
+ printf("Geometry error: could not cross boundary after extra 10 microns\n");
+ return;
+ }
+ for (i=0;i<3;i++) point[i]+=1E-6*dir[i];
+ gGeoManager->FindNode();
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;
+ }
+ if (steptot>propStep) {printf("Error\n");return;}
// we managed to cross the boundary -> in which region
newReg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
- lttcFlag++;
- if (gGeoManager->IsOutside()) return;
-
- }
+ lttcFlag++; //1
+ newLttc = (gGeoManager->IsOutside())?999999999:gGeoManager->GetCurrentNodeId()+1;
+ sLt[lttcFlag] = snext; // at 1
+ jrLt[lttcFlag] = newLttc;
+ sLt[lttcFlag+1] = 0.;
+ jrLt[lttcFlag+1] = -1;
+ // !!!!!!!!!!
+
+ while (newReg==oldReg && steptot<propStep) {
+ if (mcgeom->IsDebugging()) printf(" Entered SAME region... continue\n");
+ pst = propStep-steptot;
+ gGeoManager->FindNextBoundary(-pst);
+ snext = gGeoManager->GetStep();
+ steptot += snext;
+ if (snext<pst) {
+ printf("Found new boundary\n");
+ sLt[lttcFlag] = snext;
+ retStep = steptot; // ???
+ for (i=0;i<3;i++) point[i]+=(snext+1E-6)*dir[i];
+ steptot += 1E-6;
+ gGeoManager->FindNode();
+ if (gGeoManager->IsSameLocation()) {
+ printf("Cannot cross boundary\n");
+ break;
+ }
+ newReg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
+ newLttc = (gGeoManager->IsOutside())?999999999:gGeoManager->GetCurrentNodeId()+1;
+ if (mcgeom->IsDebugging()) printf("Found newreg=%i, newLttc=%i, lttFlag is: %i\n", newReg, newLttc, lttcFlag);
+ sLt[lttcFlag-1] += snext; // correct step in old region
+ sLt[lttcFlag] = propStep-snext;
+ jrLt[lttcFlag] = newLttc;
+ sLt[lttcFlag+1] = 0.;
+ jrLt[lttcFlag+1] = -1;
+ if (newReg != oldReg) break; // lttcFlag=1
+ lttcFlag++;
+ } else {
+ if (mcgeom->IsDebugging()) printf("Not crossing next\n");
+ lttcFlag--; //0
+ retStep=steptot;
+ sLt[lttcFlag] = retStep;
+ sLt[lttcFlag+1] = 0.;
+ jrLt[lttcFlag+1] = -1;
+ done = kTRUE;
+ }
+ }
+
+ lttcFlag++; //2
+ if (mcgeom->IsDebugging()) {
+ if (!gGeoManager->IsOutside()) {
+ printf(" ENTERED region %i, newLttc=%i in: %s\n", newReg,newLttc,gGeoManager->GetPath());
+ } else printf(" EXIT GEOMETRY: BLKHOLE reg=%i\n", newReg);
+ }
+ }
+ // no boundary within proposed step
+ lttcFlag--;
+ done = kTRUE;
+ }
+ if (mcgeom->IsDebugging()) printf("=> newReg=%i newLttc=%i lttcFlag=%i\n", newReg, newLttc, lttcFlag);
+ mcgeom->SetNextRegion(newReg, newLttc);
+ 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 (newLttc!=oldLttc) {
+ if (gGeoManager->IsOutside()) {
+ gGeoManager->SetOutside(kFALSE);
+ gGeoManager->CdTop();
+ }
+ gGeoManager->CdNode(oldLttc-1);
}
+ if (mcgeom->IsDebugging()) printf("<= G1WR (in: %s)\n", gGeoManager->GetPath());
}
//_____________________________________________________________________________
void g1rtwr()
{
- printf("=> Inside G1RTWR\n");
+ if (mcgeom->IsDebugging()) printf("========== Dummy G1RTWR\n");
}
//_____________________________________________________________________________
void conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/)
{
- printf("=> Inside CONHWR\n");
+ if (mcgeom->IsDebugging()) printf("========== Dummy CONHWR\n");
}
//_____________________________________________________________________________
-void inihwr(Int_t & /*intHist*/)
+void inihwr(Int_t &intHist)
{
- printf("=> Inside INIHWR\n");
+ 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");
+ }
}
//_____________________________________________________________________________
// 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()->GetEntriesFast()+1;
+ if (mcgeom->IsDebugging()) printf("========== Inside JOMIWR\n");
+ flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast();
+ if (mcgeom->IsDebugging()) printf("<= JOMIWR: last region=%i\n", flukaReg);
}
//_____________________________________________________________________________
Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
{
- printf("=> Inside LKDBWR\n");
- printf("Locate :(%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);
+ 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 = 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;
flagErr = newReg;
- printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
+ if (mcgeom->IsDebugging()) {
+ printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
+ printf("<= LKDBWR\n");
+ }
}
//_____________________________________________________________________________
Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
{
- printf("=> Inside LKFXWR\n");
- printf("Locate :(%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);
+ 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 = 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;
flagErr = newReg;
- printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
+ if (mcgeom->IsDebugging()) {
+ printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
+ printf("<= LKFXWR\n");
+ }
}
//_____________________________________________________________________________
Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
{
- printf("=> Inside LKMGWR\n");
- printf("Locate :(%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);
+ 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 = 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;
flagErr = newReg;
- printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
+ if (mcgeom->IsDebugging()) {
+ printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
+ printf("<= LKMGWR\n");
+ }
}
//_____________________________________________________________________________
Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
{
- printf("=> Inside LKWR\n");
- printf("Locate :(%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);
+ 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 = 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;
flagErr = newReg;
- printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
+ 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*/)
+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)
{
- printf("=> Inside NRMLWR\n");
+ 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*/)
{
- printf("=> Inside RGRPWR\n");
+ if (mcgeom->IsDebugging()) printf("=> Dummy RGRPWR\n");
}
//_____________________________________________________________________________
// For TGeo, just return the current node ID. No copy need to be made.
- printf("=> Inside ISVHWR\n");
+ 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;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff