#include "TCallf77.h"
#include "TFluka.h"
#include "TFlukaMCGeometry.h"
+#include "TFlukaConfigOption.h"
#include "TGeoManager.h"
#include "TGeoVolume.h"
#include "TObjString.h"
-
+#include "Fepisor.h"
#ifndef WIN32
# define idnrwr idnrwr_
fMatNames->Add(objstr);
nfmater++;
}
- Int_t indmat = nfmater;
+
+ fIndmat = nfmater;
// TGeoMedium *med;
// Adjust material names and add them to FLUKA list
for (i=0; i<nmater; i++) {
}
Float_t hasfield = (vol->GetMedium()->GetParam(1) > 0) ? flagfield : 0.;
-
+ out << "* Assigning material: " << vol->GetMedium()->GetMaterial()->GetName() << " to Volume: " << vol->GetName();
+ out << endl;
+
out << setw(10) << "ASSIGNMAT ";
out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
out << setw(10) << setiosflags(ios::fixed) << Double_t(idmat);
if (gSystem->AccessPathName("FlukaVmc.pemf")) Fatal("CreateFlukaMatFile", "No pemf file in working directory");
return;
}
-
- // Write peg files
- char number[20];
- Int_t countMatOK = 0;
- Int_t countElemError = 0;
- Int_t countNoStern = 0;
- Int_t countMixError = 0;
- Int_t countGas = 0;
- // Int_t countGasError = 0;
- Int_t countPemfError = 0;
- for (i=indmat; i<nfmater; 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;
- cout << "Elements with error definitions (Z not integer) " << countElemError << endl;
- cout << "Mixtures with error definitions (Z not integer) " << countMixError << endl;
- cout << "Posible Gas (rho < 0.01) " << countGas << endl;
- // cout << "Posible Gas (without pressure information) " << countGasError << endl;
+}
+
+void TFlukaMCGeometry::CreatePemfFile()
+{
+ //
+ // Steering routine to write and process peg files producing the pemf input
+ //
+ char number[20];
+ Int_t countMatOK = 0;
+ Int_t countElemError = 0;
+ Int_t countNoStern = 0;
+ Int_t countMixError = 0;
+ Int_t countGas = 0;
+ Int_t countPemfError = 0;
+ Int_t i;
+ TGeoMaterial* mat = 0x0;
+ 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++;
+ }
+ }
+ cout << "Materials (pemf created) " << countMatOK << endl;
+ cout << "Not Sternheimer par. found " << countNoStern << endl;
+ cout << "Elements with error definitions (Z not integer) " << countElemError << endl;
+ cout << "Mixtures with error definitions (Z not integer) " << countMixError << endl;
+ cout << "Posible Gas (rho < 0.01) " << countGas << endl;
+ // cout << "Posible Gas (without pressure information) " << countGasError << endl;
cout << "Pemf files Error " << countPemfError << endl;
- cout << endl << endl;
-
- sname = "cat peg/*.pemf > peg/FlukaVmc.pemf";
- gSystem->Exec(sname.Data());
- sname = "mv peg/FlukaVmc.pemf FlukaVmc.pemf";
- gSystem->Exec(sname.Data());
+ cout << endl << endl;
+
+ sname = "cat peg/*.pemf > peg/FlukaVmc.pemf";
+ gSystem->Exec(sname.Data());
+ sname = "mv peg/FlukaVmc.pemf FlukaVmc.pemf";
+ gSystem->Exec(sname.Data());
}
//_____________________________________________________________________________
delete [] wt;
}
+ Double_t ue = 3000000.; // [MeV]
+ Double_t up = 3000000.; // [MeV]
+ Double_t ae = -1.;
+ Double_t ap = -1.;
+
+
+ TObjArray* cutList = ((TFluka*) gMC)->GetListOfUserConfigs();
+ TIter next(cutList);
+ TFlukaConfigOption* proc;
+
+ 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;
+ }
+ }
+
+ if (ap == -1.) ap = TFlukaConfigOption::DefaultCut(kCUTGAM);
+ if (ae == -1.) ae = TFlukaConfigOption::DefaultCut(kCUTELE);
+
+ ap *= 1000.; // [MeV]
+ ae = (ae + 0.00051099906) * 1000.; // [MeV]
+
out << "ENER" << endl;
- out << " $INP AE=0.56099906, UE=3000000., AP=.03, UP=3000000. $END" << endl;
+ out << " $INP AE=" << ae << ", UE=" << ue <<", AP=" << ap << ", UP=" << up << " $END" << endl;
out << "PWLF" << endl;
out << " $INP NALE=300, NALG=400, NALR=100 $END" << endl;
out << "DECK" << endl;
Double_t wt[20]; //[nelems] weight fraction
Double_t density; // g/cm3
Double_t iev; // Average Ion potential (eV)
- // **** Sternheimer parameters ****
+ // **** Sternheimer parameters ****
Double_t cbar; // CBAR
Double_t x0; // X0
Double_t x1; // X1
//_____________________________________________________________________________
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,
+ 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;
gNstep++;
-/*
- if (kNstep>0) {
- gMCGeom->SetDebugMode(kTRUE);
- gFluka->SetVerbosityLevel(3);
- }
- if (kNstep>6520) {
- gMCGeom->SetDebugMode(kFALSE);
- gFluka->SetVerbosityLevel(0);
- }
- if ((kNstep%10)==0) printf("step %i\n", kNstep);
-*/
+ NORLAT.xn[0] = pSx;
+ NORLAT.xn[1] = pSy;
+ NORLAT.xn[2] = pSz;
+ NORLAT.wn[0] = pV[0];
+ NORLAT.wn[1] = pV[1];
+ NORLAT.wn[2] = pV[2];
if (gMCGeom->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);
+ }
+ Int_t olttc = oldLttc;
+ if (oldLttc<0) {
+ gGeoManager->FindNode(pSx,pSy,pSz);
+ olttc = gGeoManager->GetCurrentNodeId()+1;
}
- gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
- gGeoManager->SetCurrentDirection(pV);
- gMCGeom->SetCurrentRegion(oldReg, oldLttc);
+ gMCGeom->SetCurrentRegion(oldReg, olttc);
// Initialize default return values
lttcFlag = 0;
- jrLt[lttcFlag] = oldLttc;
+ jrLt[lttcFlag] = olttc;
sLt[lttcFlag] = propStep;
jrLt[lttcFlag+1] = -1;
sLt[lttcFlag+1] = 0.;
newReg = oldReg;
- newLttc = oldLttc;
+ newLttc = olttc;
// check if dummy boundary flag is set
Int_t curLttc, curReg;
if (gFluka->IsDummyBoundary()) {
- // printf("Dummy boundary intercepted. Point is: %f, %f, %f\n", pSx, pSy, pSz);
- Bool_t crossedDummy = (oldLttc == TFlukaMCGeometry::kLttcVirtual)?kTRUE:kFALSE;
+// printf("Dummy boundary intercepted. Point is: %f, %f, %f\n", pSx, pSy, pSz);
+ Bool_t crossedDummy = (olttc == TFlukaMCGeometry::kLttcVirtual)?kTRUE:kFALSE;
if (crossedDummy) {
// FLUKA crossed the dummy boundary - update new region/history
retStep = 0.;
}
// Reset outside flag
- if (gGeoManager->IsOutside()) {
- gGeoManager->SetOutside(kFALSE);
- gGeoManager->CdTop();
- }
+ gGeoManager->SetOutside(kFALSE);
// Reset dummy boundary flag
gFluka->SetDummyBoundary(0);
curLttc = gGeoManager->GetCurrentNodeId()+1;
curReg = gGeoManager->GetCurrentVolume()->GetNumber();
- if (oldLttc != curLttc) {
+ if (olttc != curLttc) {
// FLUKA crossed the boundary : we trust that the given point is really there,
// so we just update TGeo state
- gGeoManager->CdNode(oldLttc-1);
+ gGeoManager->CdNode(olttc-1);
curLttc = gGeoManager->GetCurrentNodeId()+1;
curReg = gGeoManager->GetCurrentVolume()->GetNumber();
if (gMCGeom->IsDebugging()) printf(" re-initialized point: curReg=%i curLttc=%i\n", curReg, curLttc);
}
- // Now the current TGeo state reflects the FLUKA state
+ // Now the current TGeo state reflects the FLUKA state
+ Double_t extra = 1.E-10;
+
+
+// printf("ERROR: (%f, %f, %f)\n",pSx,pSy,pSz);
if (gMCGeom->IsDebugging()) printf(" current path: %s\n", gGeoManager->GetPath());
- Double_t extra = 1E-6;
- Double_t tmpStep = propStep + extra;
- gGeoManager->FindNextBoundary(-tmpStep);
+ gGeoManager->SetCurrentPoint(pSx+extra*pV[0], pSy+extra*pV[1], pSz+extra*pV[2]);
+ gGeoManager->SetCurrentDirection(pV);
+ gGeoManager->FindNextBoundary(-propStep);
Double_t snext = gGeoManager->GetStep();
- // !!!!!
- if (snext<=0) {
- // FLUKA is in the wrong region, notify it
- if (gMCGeom->IsDebugging()) printf("ERROR: snext=%f\n", snext);
-// newReg = -3;
-// return;
- snext = extra;
- }
+
+ if (snext<=0.0) {
+ saf = 0.0;
+ newReg = -3;
+ sLt[lttcFlag] = 0.0;
+ if (gMCGeom->IsDebugging()) printf("BACK SCATTERING\n");
+ 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 (gMCGeom->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) {
+ saf -= extra;
+ if (saf<0) saf=0.0;
+ else saf -= saf*3.0e-09;
+ NORLAT.distn = snext;
+ NORLAT.xn[0] += snext*pV[0];
+ NORLAT.xn[1] += snext*pV[1];
+ NORLAT.xn[2] += snext*pV[2];
+// saf = 0.0; // !!! TEMPORARY FOR TESTING MAGSPHF - TO BE REMOVED
+ if (snext>propStep) {
// Next boundary further than proposed step, which is approved
retStep = propStep;
sLt[lttcFlag] = propStep;
return;
}
// The next boundary is closer. We try to cross it.
+ gGeoManager->SetCurrentPoint(pSx,pSy,pSz);
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];
+ for (i=0;i<3;i++) point[i] += (snext+1E-6)*dir[i];
+ // locate next region
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())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
if (gMCGeom->IsDebugging()) printf(" newReg=%i newLttc=%i\n", newReg, newLttc);
// We really crossed the boundary, but is it the same region ?
gMCGeom->SetNextRegion(newReg, newLttc);
- if (newReg == oldReg) {
+ if (newReg==oldReg && newLttc!=olttc) {
// Virtual boundary between replicants
if (gMCGeom->IsDebugging()) printf(" DUMMY boundary\n");
newReg = 1; // cheat FLUKA telling it it crossed the TOP region
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", gNstep);
- gGeoManager->CdNode(oldLttc-1);
+ if (newLttc!=olttc) {
+ if (gGeoManager->IsOutside()) gGeoManager->SetOutside(kFALSE);
+ gGeoManager->CdNode(olttc-1);
}
if (gMCGeom->IsDebugging()) {
printf("=> snext=%g safe=%g\n", snext, saf);
//_____________________________________________________________________________
void lkdbwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
- Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
+ Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
{
if (gMCGeom->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: 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 = gMCGeom->NofVolumes()+1;
-// newLttc = gGeoManager->GetCurrentNodeId();
- newLttc = 999999999;
- if (gMCGeom->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;
- gMCGeom->SetNextRegion(newReg, newLttc);
- flagErr = newReg;
- if (gMCGeom->IsDebugging()) {
- printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
- printf("<= LKDBWR\n");
- }
+ lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,newReg,flagErr,newLttc);
}
//_____________________________________________________________________________
void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
- Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
+ Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
{
if (gMCGeom->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: 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 = gMCGeom->NofVolumes()+1;
-// newLttc = gGeoManager->GetCurrentNodeId();
- newLttc = 999999999;
- if (gMCGeom->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;
- gMCGeom->SetNextRegion(newReg, newLttc);
- flagErr = newReg;
- if (gMCGeom->IsDebugging()) {
- printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
- printf("<= LKFXWR\n");
- }
+ lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,newReg,flagErr,newLttc);
}
//_____________________________________________________________________________
void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
- Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
+ Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
{
if (gMCGeom->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: 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 = gMCGeom->NofVolumes()+1;
-// newLttc = gGeoManager->GetCurrentNodeId();
- newLttc = 999999999;
- if (gMCGeom->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;
- gMCGeom->SetNextRegion(newReg, newLttc);
- flagErr = newReg;
- if (gMCGeom->IsDebugging()) {
- printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
- printf("<= LKMGWR\n");
- }
+ lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,newReg,flagErr,newLttc);
}
//_____________________________________________________________________________
void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
- Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
+ Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
{
if (gMCGeom->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: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
}
+ NORLAT.xn[0] = pSx;
+ NORLAT.xn[1] = pSy;
+ NORLAT.xn[2] = pSz;
+ NORLAT.wn[0] = pV[0];
+ NORLAT.wn[1] = pV[1];
+ NORLAT.wn[2] = pV[2];
TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
if (gGeoManager->IsOutside()) {
newReg = gMCGeom->NofVolumes()+1;
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())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
-// Int_t curLttc = gGeoManager->GetCurrentNodeId()+1;
-// if (gMCGeom->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 (gMCGeom->IsDebugging()) printf(" REGIONS DOEN NOT MATCH\n");
- gGeoManager->FindNode();
- curreg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
- curLttc = gGeoManager->GetCurrentNodeId()+1;
- if (gMCGeom->IsDebugging()) printf(" re-initialized point: curReg=%i curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath());
- }
-*/
+ gGeoManager->SetCurrentPoint(NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]);
+ gGeoManager->SetCurrentDirection(NORLAT.wn[0], NORLAT.wn[1], NORLAT.wn[2]);
Double_t *dnorm = gGeoManager->FindNormalFast();
flagErr = 0;
if (!dnorm) {
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]);
-// curreg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
-// curLttc = gGeoManager->GetCurrentNodeId()+1;
if (gMCGeom->IsDebugging()) {
-// printf(" final location: curReg=%i, curLttc=%i in %s\n", curreg,curLttc,gGeoManager->GetPath());
+ printf(" normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]);
printf("<= NRMLWR\n");
}
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff