#include "TCallf77.h"
#include "TFluka.h"
#include "TFlukaMCGeometry.h"
+#include "TFlukaConfigOption.h"
#include "TGeoManager.h"
#include "TGeoVolume.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*/);
void type_of_call 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*/);
+ Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
void type_of_call 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*/);
+ 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 & /*newReg*/, Int_t & /*flagErr*/, 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*/);
-};
+}
// TFluka global pointer
TFluka *gFluka = 0;
ClassImp(TFlukaMCGeometry)
-TFlukaMCGeometry* TFlukaMCGeometry::fgInstance=0;
+TFlukaMCGeometry* TFlukaMCGeometry::fgInstance= NULL;
//_____________________________________________________________________________
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;
- fNextRegion = 0;
- fNextLattice = 0;
- fRegionList = 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;
- fNextRegion = 0;
- fNextLattice = 0;
- fRegionList = 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
Int_t TFlukaMCGeometry::GetFlukaMaterial(Int_t imed) const
{
// Returns FLUKA material index for medium IMED
- TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1);
+ TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1);
if (!med) {
Error("GetFlukaMaterial", "MEDIUM %i nor found", imed);
return -1;
}
+ TGeoMaterial* mat = med->GetMaterial();
+ if (!mat->IsUsed()) return -1;
Int_t imatfl = med->GetMaterial()->GetIndex();
return imatfl;
}
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
}
TGeoMixture *mix = 0;
TGeoElement *element;
- TGeoElementTable *table = TGeoElementTable::Instance();
+ TGeoElementTable *table = gGeoManager->GetElementTable();
switch (ind) {
case 0: // AIR
- mix = new TGeoMixture("AIR", 4, 0.001205);
+ mix = new TGeoMixture("TPC_AIR", 4, 0.001205);
element = table->GetElement(6); // C
mix->DefineElement(0, element, 0.000124);
element = table->GetElement(7); // N
mix->DefineElement(3, element, 0.012827);
break;
case 1: //SDD SI CHIP
- mix = new TGeoMixture("SDD_SI", 6, 2.4485);
+ mix = new TGeoMixture("ITS_SDD_SI", 6, 2.4485);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.004367771);
element = table->GetElement(6);
mix->DefineElement(5, element, 0.09814344903);
break;
case 2: // WATER
- mix = new TGeoMixture("WATER", 2, 1.0);
+ mix = new TGeoMixture("ITS_WATER", 2, 1.0);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.111898344);
element = table->GetElement(8);
mix->DefineElement(1, element, 0.888101656);
break;
case 3: // CERAMICS
- mix = new TGeoMixture("CERAMICS", 5, 3.6);
+ mix = new TGeoMixture("ITS_CERAMICS", 5, 3.6);
element = table->GetElement(8);
mix->DefineElement(0, element, 0.59956);
element = table->GetElement(13);
mix->DefineElement(4, element, 0.0115);
break;
case 4: // EPOXY
- mix = new TGeoMixture("G10FR4", 4, 1.8);
+ mix = new TGeoMixture("MUON_G10FR4", 4, 1.8);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.19);
element = table->GetElement(6);
mix->DefineElement(3, element, 0.28);
break;
case 6: // KAPTON
- mix = new TGeoMixture("KAPTON", 4, 1.3);
+ mix = new TGeoMixture("ITS_KAPTON", 4, 1.3);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.026363415);
element = table->GetElement(6);
mix->DefineElement(3, element, 0.209238060);
break;
case 7: // INOX
- mix = new TGeoMixture("INOX", 9, 7.9);
+ mix = new TGeoMixture("ITS_INOX", 9, 7.9);
element = table->GetElement(6);
mix->DefineElement(0, element, 0.0003);
element = table->GetElement(14);
mix->DefineElement(3, element, 0.19137459);
break;
case 9: // SDD-C-AL
- mix = new TGeoMixture("SDD-C-AL", 5, 1.9837);
+ mix = new TGeoMixture("ITS_SDD-C-AL", 5, 1.9837);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.022632);
element = table->GetElement(6);
mix->DefineElement(4, element, 0.1);
break;
case 10: // X7R-CAP
- mix = new TGeoMixture("X7R-CAP", 7, 6.72);
+ mix = new TGeoMixture("ITS_X7R-CAP", 7, 6.72);
element = table->GetElement(8);
mix->DefineElement(0, element, 0.085975822);
element = table->GetElement(22);
mix->DefineElement(6, element, 0.2081768);
break;
case 11: // SDD ruby sph. Al2O3
- mix = new TGeoMixture("AL2O3", 2, 3.97);
+ mix = new TGeoMixture("ITS_AL2O3", 2, 3.97);
element = table->GetElement(8);
mix->DefineElement(0, element, 0.5293);
element = table->GetElement(13);
mix->DefineElement(1, element, 0.4707);
break;
case 12: // SDD HV microcable
- mix = new TGeoMixture("HV-CABLE", 5, 1.6087);
+ mix = new TGeoMixture("ITS_HV-CABLE", 5, 1.6087);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.01983871336);
element = table->GetElement(6);
mix->DefineElement(4, element, 0.247536);
break;
case 13: //SDD LV+signal cable
- mix = new TGeoMixture("LV-CABLE", 5, 2.1035);
+ mix = new TGeoMixture("ITS_LV-CABLE", 5, 2.1035);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.0082859922);
element = table->GetElement(6);
mix->DefineElement(4, element, 0.68572);
break;
case 14: //SDD hybrid microcab
- mix = new TGeoMixture("HYB-CAB", 5, 2.0502);
+ mix = new TGeoMixture("ITS_HYB-CAB", 5, 2.0502);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.00926228815);
element = table->GetElement(6);
mix->DefineElement(4, element, 0.64869);
break;
case 15: //SDD anode microcab
- mix = new TGeoMixture("ANOD-CAB", 5, 1.7854);
+ mix = new TGeoMixture("ITS_ANOD-CAB", 5, 1.7854);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.0128595919215);
element = table->GetElement(6);
mix->DefineElement(4, element, 0.431909);
break;
case 16: // inox/alum
- mix = new TGeoMixture("INOX-AL", 5, 3.0705);
+ mix = new TGeoMixture("ITS_INOX-AL", 5, 3.0705);
element = table->GetElement(13);
mix->DefineElement(0, element, 0.816164);
element = table->GetElement(14);
element = table->GetElement(28);
mix->DefineElement(4, element, 0.0183836);
case 17: // MYLAR
- mix = new TGeoMixture("MYLAR", 3, 1.39);
+ mix = new TGeoMixture("TPC_MYLAR", 3, 1.39);
element = table->GetElement(1);
mix->DefineElement(0, element, 0.0416667);
element = table->GetElement(6);
mix->DefineElement(2, element, 0.333333);
break;
case 18: // SPDBUS(AL+KPT+EPOX) - unknown composition
- mix = new TGeoMixture("SPDBUS", 1, 1.906);
+ mix = new TGeoMixture("ITS_SPDBUS", 1, 1.906);
element = table->GetElement(9);
mix->DefineElement(0, element, 1.);
z = element->Z();
break;
case 19: // SDD/SSD rings - unknown composition
- mix = new TGeoMixture("SDDRINGS", 1, 1.8097);
+ mix = new TGeoMixture("ITS_SDDRINGS", 1, 1.8097);
element = table->GetElement(6);
mix->DefineElement(0, element, 1.);
z = element->Z();
break;
case 20: // SPD end ladder - unknown composition
- mix = new TGeoMixture("SPDEL", 1, 3.6374);
+ mix = new TGeoMixture("ITS_SPDEL", 1, 3.6374);
element = table->GetElement(22);
mix->DefineElement(0, element, 1.);
z = element->Z();
break;
case 21: // SDD end ladder - unknown composition
- mix = new TGeoMixture("SDDEL", 1, 0.3824);
+ mix = new TGeoMixture("ITS_SDDEL", 1, 0.3824);
element = table->GetElement(30);
mix->DefineElement(0, element, 1.);
z = element->Z();
break;
case 22: // SSD end ladder - unknown composition
- mix = new TGeoMixture("SSDEL", 1, 0.68);
+ mix = new TGeoMixture("ITS_SSDEL", 1, 0.68);
element = table->GetElement(16);
mix->DefineElement(0, element, 1.);
z = element->Z();
Int_t i,j,idmat;
Int_t counttothree, nelem;
Double_t a,z,rho, w;
- TGeoElementTable *table = TGeoElementTable::Instance();
+ TGeoElementTable *table = gGeoManager->GetElementTable();
TGeoElement *element;
element = table->GetElement(13);
element->SetTitle("ALUMINUM"); // this is how FLUKA likes it ...
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++) {
mat = (TGeoMaterial*)matlist->At(i);
- if (!mat->IsUsed()) continue;
+ if (!mat->IsUsed()) continue;
z = mat->GetZ();
a = mat->GetA();
rho = mat->GetDensity();
}
matname = mat->GetName();
FlukaMatName(matname);
-/*
- // material with one element: create it as mixture since it can be duplicated
- if (!mat->IsMixture()) {
- // normal materials
- mix = new TGeoMixture(matname.Data(), 1, rho);
- mix->DefineElement(0, mat->GetElement(), 1.);
- mat->SetIndex(nfmater+3);
- for (j=0; j<nmed; j++) {
- med = (TGeoMedium*)medlist->At(j);
- if (med->GetMaterial() == mat) {
- med->SetMaterial(mix);
- if (mat->GetCerenkovProperties()) {
- mix->SetCerenkovProperties(mat->GetCerenkovProperties());
- mat->SetCerenkovProperties(0);
- }
- break;
- }
- }
- mat = (TGeoMaterial*)mix;
- }
-*/
+
mat->SetIndex(nfmater+3);
objstr = new TObjString(matname.Data());
fMatList->Add(mat);
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);
}
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);
out << setw(10) << "0.0";
out << endl;
}
+ // dummy region
+ idmat = 2; // vacuum
+ fDummyRegion = nvols+1;
+ out << "* Dummy region: " << endl;
+ out << setw(10) << "ASSIGNMAT ";
+ out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
+ out << setw(10) << setiosflags(ios::fixed) << idmat;
+ out << setw(10) << setiosflags(ios::fixed) << fDummyRegion;
+ out << setw(10) << "0.0";
+ out << setw(10) << "0.0";
+ out << setw(10) << "0.0";
+ out << setw(10) << "0.0" << endl;
out.close();
fLastMaterial = nfmater+2;
-
- if (!gFluka->IsGeneratePemf()) {
- 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
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;
}
//_____________________________________________________________________________
-void TFlukaMCGeometry::SetMreg(Int_t mreg)
+void TFlukaMCGeometry::SetMreg(Int_t mreg, Int_t lttc)
{
// Update if needed next history;
- if (gFluka->GetDummyBoundary()==2) {
- gGeoManager->CdNode(fNextLattice-1);
- return;
- }
- Int_t curreg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
- if (mreg==curreg) return;
- if (mreg==fNextRegion) {
- if (fNextLattice!=999999999) gGeoManager->CdNode(fNextLattice-1);
+// if (gFluka->GetDummyBoundary()==2) {
+// gGeoManager->CdNode(fNextLattice-1);
+// return;
+// }
+ if (lttc == TFlukaMCGeometry::kLttcOutside) {
+ fCurrentRegion = NofVolumes()+2;
+ fCurrentLattice = lttc;
+ gGeoManager->CdTop();
+ gGeoManager->SetOutside(kTRUE);
+ }
+ if (lttc == TFlukaMCGeometry::kLttcVirtual) return;
+ if (lttc <=0) {
+ Error("TFlukaMCGeometry::SetMreg","Invalide lattice %i",lttc);
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);
+ }
+ fCurrentRegion = mreg;
+ fCurrentLattice = lttc;
+
+ Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
+ if (crtlttc == lttc) return;
+ gGeoManager->CdNode(lttc-1);
}
//_____________________________________________________________________________
//_____________________________________________________________________________
void TFlukaMCGeometry::FlukaMatName(TString &str) const
{
+// Strip the detector name
+
+ 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);
+ }
+ tokens->Clear();
+ delete tokens;
+
// Convert a name to upper case 8 chars.
ToFlukaString(str);
Int_t ilast;
}
}
}
-
+
//______________________________________________________________________________
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
//_____________________________________________________________________________
//_____________________________________________________________________________
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,
+ 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);
-*/
- 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);
- }
- gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
- gGeoManager->SetCurrentDirection(pV);
- gMCGeom->SetCurrentRegion(oldReg, oldLttc);
+// if( (gNstep > 43912170 && gNstep < 43912196 ) ||
+// (gNstep > 47424560 && gNstep < 47424581 ) ||
+// (gNstep > 54388266 && gNstep < 54388319 )
+// ) gMCGeom->SetDebugMode();
+// else gMCGeom->SetDebugMode(kFALSE);
+
+ NORLAT.xn[0] = pSx;
+ NORLAT.xn[1] = pSy;
+ NORLAT.xn[2] = pSz;
+
+ Int_t olttc = oldLttc;
+ if (oldLttc<=0) {
+ gGeoManager->FindNode(pSx,pSy,pSz);
+ olttc = gGeoManager->GetCurrentNodeId()+1;
+ oldReg = gGeoManager->GetCurrentVolume()->GetNumber();
+ }
+
+ if (gMCGeom->IsDebugging()) cout << "g1wr gNstep=" << gNstep
+ << " oldReg="<< oldReg <<" olttc="<< olttc
+ << " track=" << TRACKR.ispusr[mkbmx2-1] << endl;
+
+ Int_t ccreg,cclat;
+ gMCGeom->GetCurrentRegion(ccreg,cclat);
+ Bool_t crossed = (ccreg==oldReg && cclat==olttc)?kFALSE:kTRUE;
+
+ 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;
- // check if dummy boundary flag is set
+ newLttc = olttc;
+ Bool_t crossedDummy = (oldReg == gFluka->GetDummyRegion())?kTRUE:kFALSE;
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;
- if (crossedDummy) {
- // FLUKA crossed the dummy boundary - update new region/history
- retStep = 0.;
- saf = 0.;
- gMCGeom->GetNextRegion(newReg, newLttc);
- gMCGeom->SetMreg(newReg);
- if (gMCGeom->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.;
- gFluka->SetDummyBoundary(0);
- return;
- }
- }
-
+ if (crossedDummy) {
+ // FLUKA crossed the dummy boundary - update new region/history
+ retStep = TGeoShape::Tolerance();
+ saf = 0.;
+ gMCGeom->GetNextRegion(newReg, newLttc);
+ gMCGeom->SetMreg(newReg, newLttc);
+ sLt[lttcFlag] = TGeoShape::Tolerance(); // null step in current region
+ lttcFlag++;
+ jrLt[lttcFlag] = newLttc;
+ 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
- if (gGeoManager->IsOutside()) {
- gGeoManager->SetOutside(kFALSE);
- gGeoManager->CdTop();
- }
-
- // Reset dummy boundary flag
- gFluka->SetDummyBoundary(0);
-
+ gGeoManager->SetOutside(kFALSE);
+
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
- if (gMCGeom->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 (gMCGeom->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 (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) {
+ // Now the current TGeo state reflects the FLUKA state
+
+ gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
+ gGeoManager->SetCurrentDirection(pV);
+
+ if (crossed) {
+ gGeoManager->FindNextBoundaryAndStep(propStep);
+ saf = 0.0;
+ } else {
+ gGeoManager->FindNextBoundaryAndStep(propStep, kTRUE);
+ saf = gGeoManager->GetSafeDistance();
+ if (saf<0) saf=0.0;
+ saf -= saf*3.0e-09;
+ }
+
+ 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[1] += snext*pV[1];
+ NORLAT.xn[2] += snext*pV[2];
+ if (!gGeoManager->IsOnBoundary()) {
// Next boundary further than proposed step, which is approved
+ if (saf>propStep) saf = propStep;
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();
+ if (saf>snext) saf = snext; // Safety should be less than the proposed step if a boundary will be crossed
+ gGeoManager->SetCurrentPoint(pSx,pSy,pSz);
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);
+ newReg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+2):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) {
+
+ Int_t pid = TRACKR.jtrack;
+ if ( ((newReg==oldReg && newLttc!=olttc) || (oldReg!=newReg && olttc==newLttc) ) && pid!=-1) {
// Virtual boundary between replicants
- if (gMCGeom->IsDebugging()) printf(" DUMMY boundary\n");
- newReg = 1; // cheat FLUKA telling it it crossed the TOP region
+ newReg = gFluka->GetDummyRegion();
newLttc = TFlukaMCGeometry::kLttcVirtual;
- // mark that next boundary is virtual
- gFluka->SetDummyBoundary(1);
- }
+ 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.;
- 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);
- }
+ 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]);
- }
- if (gMCGeom->IsDebugging()) printf("<= G1WR (in: %s)\n", gGeoManager->GetPath());
+ }
}
//_____________________________________________________________________________
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");
- }
+ }
}
//_____________________________________________________________________________
// number of regions (volumes in TGeo)
// build application geometry
if (gMCGeom->IsDebugging()) printf("========== Inside JOMIWR\n");
- flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast();
+ flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast()+1;
if (gMCGeom->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)
+ 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 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,flagErr,newReg,newLttc);
}
//_____________________________________________________________________________
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)
+ 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 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,flagErr,newReg,newLttc);
}
//_____________________________________________________________________________
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)
+ 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,flagErr,newReg,newLttc);
}
//_____________________________________________________________________________
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)
+ 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 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);
}
+ flagErr = 0;
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;
+ newReg = gMCGeom->NofVolumes()+2;
+ newLttc = TFlukaMCGeometry::kLttcOutside;
+ gGeoManager->SetOutside(kFALSE);
+ if (oldLttc>0 && oldLttc<newLttc) gGeoManager->CdNode(oldLttc-1);
return;
}
+ gGeoManager->SetOutside(kFALSE);
newReg = node->GetVolume()->GetNumber();
newLttc = gGeoManager->GetCurrentNodeId()+1;
- gMCGeom->SetNextRegion(newReg, newLttc);
- flagErr = newReg;
+ if (oldLttc==TFlukaMCGeometry::kLttcOutside || oldLttc==0) return;
+
+ Int_t dummy = gFluka->GetDummyRegion();
+ if (oldReg==dummy) {
+ Int_t newreg1, newlttc1;
+ gMCGeom->GetNextRegion(newreg1, newlttc1);
+ 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 = 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(" out: newReg=%i newLttc=%i in %s\n", newReg, newLttc, gGeoManager->GetPath());
- printf("<= LKWR\n");
+ 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);
}
-// 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(pVx, pVy, pVz);
Double_t *dnorm = gGeoManager->FindNormalFast();
flagErr = 0;
if (!dnorm) {
norml[0] = -pVx;
norml[1] = -pVy;
norml[2] = -pVz;
- }
- 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;
+ } else {
+ norml[0] = -dnorm[0];
+ norml[1] = -dnorm[1];
+ norml[2] = -dnorm[2];
+ }
+
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");
- }
+ }
+
}
//_____________________________________________________________________________
// 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