X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TFluka%2FTFlukaMCGeometry.cxx;h=1199d3cf89de1b5615922d88916ec9d4d7990652;hb=9968e86cb4f48f01d9997d8922c3820d8d46f268;hp=166d6dbdb84100db29822c5a1f6fe18b72e30127;hpb=784c3bdd1a97531f3e6d5257ee01cda602b3b700;p=u%2Fmrichter%2FAliRoot.git diff --git a/TFluka/TFlukaMCGeometry.cxx b/TFluka/TFlukaMCGeometry.cxx index 166d6dbdb84..1199d3cf89d 100644 --- a/TFluka/TFlukaMCGeometry.cxx +++ b/TFluka/TFlukaMCGeometry.cxx @@ -24,14 +24,18 @@ // Author: Andrei Gheata 10/07/2003 #include "Riostream.h" +#include "TList.h" #include "TCallf77.h" #include "TFluka.h" +#include "TSystem.h" #include "TFlukaMCGeometry.h" #include "TFlukaConfigOption.h" #include "TGeoManager.h" #include "TGeoVolume.h" #include "TObjString.h" -#include "Fepisor.h" +#include "Fsourcm.h" +#include "Ftrackr.h" +#include "Fstepsz.h" //(STEPSZ) fluka common #ifndef WIN32 # define idnrwr idnrwr_ @@ -78,37 +82,37 @@ extern "C" 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; @@ -117,43 +121,53 @@ Int_t gNstep = 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; } //_____________________________________________________________________________ @@ -172,14 +186,6 @@ TFlukaMCGeometry::~TFlukaMCGeometry() // // private methods // -//_____________________________________________________________________________ -TFlukaMCGeometry::TFlukaMCGeometry(const TFlukaMCGeometry &) - : TNamed() -{ - // - // Copy constructor - // -} //_____________________________________________________________________________ Double_t* TFlukaMCGeometry::CreateDoubleArray(Float_t* array, Int_t size) const @@ -219,11 +225,13 @@ Int_t TFlukaMCGeometry::GetMedium() 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; } @@ -238,14 +246,16 @@ Int_t *TFlukaMCGeometry::GetRegionList(Int_t imed, Int_t &nreg) 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) @@ -257,14 +267,16 @@ 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 @@ -296,10 +308,10 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) } 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 @@ -310,7 +322,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) 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); @@ -325,14 +337,14 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) 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); @@ -345,7 +357,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) 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); @@ -367,7 +379,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) 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); @@ -378,7 +390,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) 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); @@ -410,7 +422,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) 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); @@ -423,7 +435,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) 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); @@ -440,14 +452,14 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) 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); @@ -460,7 +472,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) 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); @@ -473,7 +485,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) 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); @@ -486,7 +498,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) 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); @@ -499,7 +511,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) 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); @@ -511,7 +523,7 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) 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); @@ -520,31 +532,31 @@ TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) 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(); @@ -565,8 +577,7 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname) // program load the right cross sections, and equal to the names included in // the .pemf. Otherwise the user must define the LOW-MAT CARDS, and make his // own .pemf, in order to get the right cross sections loaded in memory. - - + // Materials defined by FLUKA TString sname; gGeoManager->Export("flgeom.root"); if (fname) sname = fname; @@ -582,21 +593,19 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname) 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 ... element = table->GetElement(15); element->SetTitle("PHOSPHO"); // same story ... -// element = table->GetElement(10); -// element->SetTitle("ARGON"); // NEON not in neutron xsec table Int_t nelements = table->GetNelements(); TList *matlist = gGeoManager->GetListOfMaterials(); -// TList *medlist = gGeoManager->GetListOfMedia(); -// Int_t nmed = medlist->GetSize(); TIter next(matlist); Int_t nmater = matlist->GetSize(); - Int_t nfmater = 0; + Int_t nfmater = 0; + Int_t newind = 26; // here non predefined materials start + TGeoMaterial *mat; TGeoMixture *mix = 0; TString matname; @@ -611,7 +620,14 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname) rho = 0.999; mat = new TGeoMaterial(matname, element->A(), element->Z(), rho); - mat->SetIndex(nfmater+3); + // Check if the element has been predefined by FLUKA + Int_t pmid = GetPredefinedMaterialId(Int_t(element->Z())); + if (pmid > 0) { + mat->SetIndex(pmid); + } else { + mat->SetIndex(newind++); + } + mat->SetUsed(kTRUE); fMatList->Add(mat); objstr = new TObjString(matname.Data()); @@ -620,11 +636,10 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname) } fIndmat = nfmater; -// TGeoMedium *med; // Adjust material names and add them to FLUKA list for (i=0; iAt(i); - if (!mat->IsUsed()) continue; + if (!mat->IsUsed()) continue; z = mat->GetZ(); a = mat->GetA(); rho = mat->GetDensity(); @@ -634,28 +649,8 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname) } 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; jAt(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); + + mat->SetIndex(newind++); objstr = new TObjString(matname.Data()); fMatList->Add(mat); fMatNames->Add(objstr); @@ -775,9 +770,10 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname) 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; jAt(j); @@ -798,13 +794,20 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname) 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(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; - } } void TFlukaMCGeometry::CreatePemfFile() @@ -824,32 +827,32 @@ void TFlukaMCGeometry::CreatePemfFile() TString sname; for (i = fIndmat; i < fLastMaterial - 2; i++) { - printf("Write Peg Files %d\n", i); - - mat = (TGeoMaterial*)fMatList->At(i); - if (!mat->IsUsed()) continue; - sname = "mat"; - sprintf(number, "%d", i); - sname.Append(number); - cout << endl; - cout << endl; - cout << "******************************************************************************" << endl; - cout << "******************************************************************************" << endl; - cout << endl; - WritePegFile(i, &countNoStern, &countElemError, &countMixError, &countGas); - sname.Prepend("$FLUPRO/pemf/rpemf peg/"); - gSystem->Exec(sname.Data()); - - // check if the pemf file was created - TString sname = Form("peg/mat%d.pemf", i); - ifstream in( sname.Data() ); - if ( in ) { - countMatOK++; - in.close(); - } else { - cout << "ERROR Fail to create the pemf file " << sname << endl; - countPemfError++; - } + printf("Write Peg Files %d\n", i); + + mat = (TGeoMaterial*)fMatList->At(i); + if (!mat->IsUsed()) continue; + sname = "mat"; + sprintf(number, "%d", i); + sname.Append(number); + cout << endl; + cout << endl; + cout << "******************************************************************************" << endl; + cout << "******************************************************************************" << endl; + cout << endl; + WritePegFile(i, &countNoStern, &countElemError, &countMixError, &countGas); + sname.Prepend("$FLUPRO/pemf/rpemf peg/"); + gSystem->Exec(sname.Data()); + + // check if the pemf file was created + TString sname = Form("peg/mat%d.pemf", i); + ifstream in( sname.Data() ); + if ( in ) { + countMatOK++; + in.close(); + } else { + cout << "ERROR Fail to create the pemf file " << sname << endl; + countPemfError++; + } } cout << "Materials (pemf created) " << countMatOK << endl; cout << "Not Sternheimer par. found " << countNoStern << endl; @@ -1022,11 +1025,11 @@ void TFlukaMCGeometry::WritePegFile(Int_t imat, Int_t *NoStern, Int_t *ElemError while((proc = (TFlukaConfigOption*)next())) { if (proc->Medium() == mat->GetIndex()) { - ap = proc->Cut(kCUTGAM); - ae = proc->Cut(kCUTELE); - if (ap == -1.) ap = TFlukaConfigOption::DefaultCut(kCUTGAM); - if (ae == -1.) ae = TFlukaConfigOption::DefaultCut(kCUTELE); - break; + ap = proc->Cut(kCUTGAM); + ae = proc->Cut(kCUTELE); + if (ap == -1.) ap = TFlukaConfigOption::DefaultCut(kCUTGAM); + if (ae == -1.) ae = TFlukaConfigOption::DefaultCut(kCUTELE); + break; } } @@ -1070,6 +1073,10 @@ Double_t * TFlukaMCGeometry::GetISSB(Double_t rho, Int_t nElem, Double_t *zelem, 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; @@ -1219,25 +1226,31 @@ Int_t TFlukaMCGeometry::GetElementIndex(Int_t z) const } //_____________________________________________________________________________ -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); + while (gGeoManager->GetCurrentVolume()->IsAssembly()) gGeoManager->CdUp(); } //_____________________________________________________________________________ @@ -1278,6 +1291,18 @@ void TFlukaMCGeometry::ToFlukaString(TString &str) const //_____________________________________________________________________________ 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; @@ -1306,7 +1331,7 @@ void TFlukaMCGeometry::FlukaMatName(TString &str) const } } } - + //______________________________________________________________________________ void TFlukaMCGeometry::Vname(const char *name, char *vname) const { @@ -1322,8 +1347,40 @@ void TFlukaMCGeometry::Vname(const char *name, char *vname) const } -// FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers +//______________________________________________________________________________ +Int_t TFlukaMCGeometry::GetNstep() +{ + // return gNstep for debug propose + return gNstep; +} + + +Int_t TFlukaMCGeometry::GetPredefinedMaterialId(Int_t z) const +{ +// Get predifined material id from Z if present in list + const Int_t kMax = 25; + + static Int_t idFluka[kMax] = + {-1, -1, 1, 2, 4, 6, 7, 8, 12, 13, + 26, 29, 47, 14, 79, 80, 82, 73, 11, 18, + 20, 50, 74, 22, 28}; + + Int_t id = -1; + for (Int_t i = 0; i < kMax; i++) + { + if (z == idFluka[i]) { + id = i + 1; + break; + } + + } + + return id; +} + +// FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers +// //_____________________________________________________________________________ Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/) { @@ -1342,29 +1399,38 @@ Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/) void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz, Double_t *pV, Int_t &oldReg , const Int_t &oldLttc, Double_t &propStep, Int_t &/*nascFlag*/, Double_t &retStep, Int_t &newReg, - Double_t &saf, Int_t &newLttc, Int_t <tcFlag, + Double_t &saf, Int_t &newLttc, Int_t <tcFlag, Double_t *sLt, Int_t *jrLt) { // Initialize FLUKa point and direction; + static Int_t ierr = 0; gNstep++; +// gMCGeom->SetDebugMode(kTRUE); + NORLAT.xn[0] = pSx; NORLAT.xn[1] = pSy; NORLAT.xn[2] = pSz; - 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) { + 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; + cout << " point: (" << pSx << ", " << pSy << ", " << pSz << ") dir: (" + << pV[0] << ", " << pV[1] << ", " << pV[2] << ")" << endl; + } + + + Int_t ccreg=0,cclat=0; + gMCGeom->GetCurrentRegion(ccreg,cclat); + Bool_t crossed = (ccreg==oldReg && cclat==olttc)?kFALSE:kTRUE; + gMCGeom->SetCurrentRegion(oldReg, olttc); // Initialize default return values lttcFlag = 0; @@ -1374,35 +1440,27 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz, sLt[lttcFlag+1] = 0.; newReg = oldReg; newLttc = olttc; - // check if dummy boundary flag is set + 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 = (olttc == 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 gGeoManager->SetOutside(kFALSE); - - // Reset dummy boundary flag - gFluka->SetDummyBoundary(0); - + curLttc = gGeoManager->GetCurrentNodeId()+1; curReg = gGeoManager->GetCurrentVolume()->GetNumber(); if (olttc != curLttc) { @@ -1411,106 +1469,143 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz, 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 - 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()); - gGeoManager->SetCurrentPoint(pSx+extra*pV[0], pSy+extra*pV[1], pSz+extra*pV[2]); + gGeoManager->SetCurrentPoint(pSx, pSy, pSz); gGeoManager->SetCurrentDirection(pV); - gGeoManager->FindNextBoundary(-propStep); - Double_t snext = gGeoManager->GetStep(); - - if (snext<=0.0) { - saf = 0.0; - newReg = -3; - sLt[lttcFlag] = 0.0; - if (gMCGeom->IsDebugging()) printf("BACK SCATTERING\n"); + Double_t pt[3], local[3], ldir[3]; + Int_t pid = TRACKR.jtrack; + pt[0] = pSx; + pt[1] = pSy; + pt[2] = pSz; + gGeoManager->MasterToLocal(pt,local); + gGeoManager->MasterToLocalVect(pV,ldir); +/* + Bool_t valid = gGeoManager->GetCurrentVolume()->Contains(local); + if (!valid) { + printf("location not valid in %s pid=%i\n", gGeoManager->GetPath(),pid); + printf("local:(%f, %f, %f) ldir:(%f, %f, %f)\n", local[0],local[1],local[2],ldir[0],ldir[1],ldir[2]); +// gGeoManager->FindNode(); +// printf(" -> actual location: %s\n", gGeoManager->GetPath()); + } +*/ + Double_t pstep = propStep; + Double_t snext = propStep; + const Double_t epsil = 0.9999999 * TGeoShape::Tolerance(); + // This should never happen !!! + if (pstepFindNextBoundaryAndStep(pstep); + snext = gGeoManager->GetStep(); + saf = 0.; + if (snext <= TGeoShape::Tolerance()) { +// printf("FLUKA crossed boundary but snext=%g\n", snext); + ierr++; + snext = epsil; + } else { + snext += TGeoShape::Tolerance(); + ierr = 0; + } + } else { + gGeoManager->FindNextBoundaryAndStep(pstep, kTRUE); + snext = gGeoManager->GetStep(); + saf = gGeoManager->GetSafeDistance(); + if (snext <= TGeoShape::Tolerance()) { +// printf("FLUKA put particle on bondary without crossing\n"); + ierr++; + snext = epsil; + saf = 0.; + } else { + snext += TGeoShape::Tolerance(); + ierr = 0; + } + if (saf<0) saf=0.; + saf -= saf*3.0e-09; + } +// if (ierr>1) { +// printf("%d snext=%g\n", ierr, snext); +// } + if (ierr == 10) { +// printf("Too many null steps - sending error code -33...\n"); + newReg = -33; // Error code + ierr = 0; return; - } - - snext += extra; - saf = gGeoManager->GetSafeDistance(); - saf -= extra; - if (saf<0) saf=0.0; - else saf -= saf*3.0e-09; + } + + PAREM.dist = snext; 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) { + 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. + if (saf>snext) saf = snext; // Safety should be less than the proposed step if a boundary will be crossed 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+1E-6)*dir[i]; - // locate next region - gGeoManager->FindNode(); newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1; - gGeoManager->SetCurrentPoint(pt); - 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 && newLttc!=olttc) { + + 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!=olttc) { - if (gGeoManager->IsOutside()) gGeoManager->SetOutside(kFALSE); - gGeoManager->CdNode(olttc-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; iIsDebugging()) 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; } @@ -1519,7 +1614,7 @@ void inihwr(Int_t &intHist) if (gMCGeom->IsDebugging()) { printf(" --- current path: %s\n", gGeoManager->GetPath()); printf("<= INIHWR\n"); - } + } } //_____________________________________________________________________________ @@ -1530,100 +1625,120 @@ void jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/ // 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) + 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: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); } - lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,newReg,flagErr,newLttc); + 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) + 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: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); } - lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,newReg,flagErr,newLttc); + 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) + 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: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); } - lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,newReg,flagErr,newLttc); + 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) + 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: 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); + flagErr = 0; + Double_t epsil = 1000.*TGeoShape::Tolerance(); + TGeoNode *node = gGeoManager->FindNode(pSx+epsil*pV[0], pSy+epsil*pV[1], pSz+epsil*pV[2]); if (gGeoManager->IsOutside()) { - newReg = gMCGeom->NofVolumes()+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 && oldLttcCdNode(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); } gGeoManager->SetCurrentPoint(NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]); - gGeoManager->SetCurrentDirection(NORLAT.wn[0], NORLAT.wn[1], NORLAT.wn[2]); + gGeoManager->SetCurrentDirection(pVx, pVy, pVz); Double_t *dnorm = gGeoManager->FindNormalFast(); flagErr = 0; if (!dnorm) { @@ -1632,14 +1747,17 @@ void nrmlwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, norml[0] = -pVx; norml[1] = -pVy; norml[2] = -pVz; - } - norml[0] = -dnorm[0]; - norml[1] = -dnorm[1]; - norml[2] = -dnorm[2]; + } else { + norml[0] = -dnorm[0]; + norml[1] = -dnorm[1]; + norml[2] = -dnorm[2]; + } + if (gMCGeom->IsDebugging()) { printf(" normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]); printf("<= NRMLWR\n"); - } + } + } //_____________________________________________________________________________ @@ -1665,7 +1783,7 @@ Int_t isvhwr(const Int_t &check, const Int_t & intHist) // For TGeo, just return the current node ID. No copy need to be made. - if (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());