X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TFluka%2FTFlukaMCGeometry.cxx;h=1199d3cf89de1b5615922d88916ec9d4d7990652;hb=9968e86cb4f48f01d9997d8922c3820d8d46f268;hp=6c8fcb96a22cd86b372029a407e112a5fc069c40;hpb=2573ac89dc45ded8d40502ea8a1feff8fdd43b01;p=u%2Fmrichter%2FAliRoot.git diff --git a/TFluka/TFlukaMCGeometry.cxx b/TFluka/TFlukaMCGeometry.cxx index 6c8fcb96a22..1199d3cf89d 100644 --- a/TFluka/TFlukaMCGeometry.cxx +++ b/TFluka/TFlukaMCGeometry.cxx @@ -1,14 +1,41 @@ -// @(#):$Name$:$Id$ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * Permission to use, copy, modify and distribute this software and its * + * documentation strictly for non-commercial purposes is hereby granted * + * without fee, provided that the above copyright notice appears in all * + * copies and that both the copyright notice and this permission notice * + * appear in the supporting documentation. The authors make no claims * + * about the suitability of this software for any purpose. It is * + * provided "as is" without express or implied warranty. * + **************************************************************************/ + +// $Id$ + +// Class TFlukaMCGeometry +// -------------------- +// Implementation of the TVirtualMCGeometry interface +// for defining and using TGeo geometry with FLUKA. +// This allows the FLUKA MonteCarlo to run with the TGeo +// geometrical modeller // Author: Andrei Gheata 10/07/2003 -#include "TObjString.h" +#include "Riostream.h" +#include "TList.h" +#include "TCallf77.h" #include "TFluka.h" -//#include "TVirtualMCApplication.h" +#include "TSystem.h" #include "TFlukaMCGeometry.h" +#include "TFlukaConfigOption.h" #include "TGeoManager.h" #include "TGeoVolume.h" - -#include "TCallf77.h" +#include "TObjString.h" +#include "Fsourcm.h" +#include "Ftrackr.h" +#include "Fstepsz.h" //(STEPSZ) fluka common #ifndef WIN32 # define idnrwr idnrwr_ @@ -55,68 +82,92 @@ 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 *fluka = 0; -TFlukaMCGeometry *mcgeom = 0; +TFluka *gFluka = 0; +TFlukaMCGeometry *gMCGeom = 0; +Int_t gNstep = 0; ClassImp(TFlukaMCGeometry) -TFlukaMCGeometry* TFlukaMCGeometry::fgInstance=0; +TFlukaMCGeometry* TFlukaMCGeometry::fgInstance= NULL; //_____________________________________________________________________________ TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title) - : TVirtualMCGeometry(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 // - fLastMaterial = 0; - fluka = (TFluka*)gMC; - mcgeom = this; + gFluka = (TFluka*)gMC; + gMCGeom = this; + gNstep = 0; } //_____________________________________________________________________________ TFlukaMCGeometry::TFlukaMCGeometry() - : TVirtualMCGeometry() -{ + :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 // - fLastMaterial = 0; - fluka = (TFluka*)gMC; - mcgeom = this; + gFluka = (TFluka*)gMC; + gMCGeom = this; + gNstep = 0; } //_____________________________________________________________________________ @@ -126,20 +177,15 @@ TFlukaMCGeometry::~TFlukaMCGeometry() // Destructor // fgInstance=0; + if (fRegionList) delete [] fRegionList; + if (fMatList) delete fMatList; + if (fMatNames) {fMatNames->Delete(); delete fMatNames;} if (gGeoManager) delete gGeoManager; } // // private methods // -//_____________________________________________________________________________ -TFlukaMCGeometry::TFlukaMCGeometry(const TFlukaMCGeometry &) - : TVirtualMCGeometry() -{ - // - // Copy constructor - // -} //_____________________________________________________________________________ Double_t* TFlukaMCGeometry::CreateDoubleArray(Float_t* array, Int_t size) const @@ -161,174 +207,8 @@ Double_t* TFlukaMCGeometry::CreateDoubleArray(Float_t* array, Int_t size) const } // // public methods -//_____________________________________________________________________________ -void TFlukaMCGeometry::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z, - Float_t &dens, Float_t &radl, Float_t &absl, - Float_t* /*ubuf*/, Int_t& /*nbuf*/) -{ - printf("Gfmate %i\n", imat); - TGeoMaterial *mat; - TIter next (gGeoManager->GetListOfMaterials()); - while ((mat = (TGeoMaterial*)next())) { - if (mat->GetUniqueID() == (UInt_t)imat) break; - } - if (!mat) { - Error("Gfmate", "no material with index %i found", imat); - return; - } - sprintf(name, "%s", mat->GetName()); - a = mat->GetA(); - z = mat->GetZ(); - dens = mat->GetDensity(); - radl = mat->GetRadLen(); - absl = mat->GetIntLen(); - printf(" ->material found : %s a=%g, z=%g, dens=%g, radl=%g, absl=%g\n", name, a,z,dens,radl,absl); -} - -//_____________________________________________________________________________ -void TFlukaMCGeometry::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z, - Double_t &dens, Double_t &radl, Double_t &absl, - Double_t* /*ubuf*/, Int_t& /*nbuf*/) -{ - printf("Gfmate %i\n", imat); - TGeoMaterial *mat; - TIter next (gGeoManager->GetListOfMaterials()); - while ((mat = (TGeoMaterial*)next())) { - if (mat->GetUniqueID() == (UInt_t)imat) break; - } - if (!mat) { - Error("Gfmate", "no material with index %i found", imat); - return; - } - sprintf(name, "%s", mat->GetName()); - a = mat->GetA(); - z = mat->GetZ(); - dens = mat->GetDensity(); - radl = mat->GetRadLen(); - absl = mat->GetIntLen(); - printf(" ->material found : %s a=%g, z=%g, dens=%g, radl=%g, absl=%g\n", name, a,z,dens,radl,absl); -} - -//_____________________________________________________________________________ -void TFlukaMCGeometry::Material(Int_t& kmat, const char* name, Double_t a, Double_t z, - Double_t dens, Double_t radl, Double_t absl, Float_t* buf, - Int_t nwbuf) -{ - // - // Defines a Material - // - // kmat number assigned to the material - // name material name - // a atomic mass in au - // z atomic number - // dens density in g/cm3 - // absl absorbtion length in cm - // if >=0 it is ignored and the program - // calculates it, if <0. -absl is taken - // radl radiation length in cm - // if >=0 it is ignored and the program - // calculates it, if <0. -radl is taken - // buf pointer to an array of user words - // nbuf number of user words - // - - Double_t* dbuf = CreateDoubleArray(buf, nwbuf); - Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf); - delete [] dbuf; -} - -//_____________________________________________________________________________ -void TFlukaMCGeometry::Material(Int_t& kmat, const char* name, Double_t a, Double_t z, - Double_t dens, Double_t radl, Double_t absl, Double_t* /*buf*/, - Int_t /*nwbuf*/) -{ - // - // Defines a Material - // - // kmat number assigned to the material - // name material name - // a atomic mass in au - // z atomic number - // dens density in g/cm3 - // absl absorbtion length in cm - // if >=0 it is ignored and the program - // calculates it, if <0. -absl is taken - // radl radiation length in cm - // if >=0 it is ignored and the program - // calculates it, if <0. -radl is taken - // buf pointer to an array of user words - // nbuf number of user words - // - kmat = gGeoManager->GetListOfMaterials()->GetSize(); - gGeoManager->Material(name, a, z, dens, kmat, radl, absl); - printf("Material %s: kmat=%i, a=%g, z=%g, dens=%g\n", name, kmat, a, z, dens); -} -//_____________________________________________________________________________ -void TFlukaMCGeometry::Mixture(Int_t& kmat, const char* name, Float_t* a, Float_t* z, - Double_t dens, Int_t nlmat, Float_t* wmat) -{ - // - // Defines mixture OR COMPOUND IMAT as composed by - // THE BASIC NLMAT materials defined by arrays A,Z and WMAT - // - // If NLMAT > 0 then wmat contains the proportion by - // weights of each basic material in the mixture. - // - // If nlmat < 0 then WMAT contains the number of atoms - // of a given kind into the molecule of the COMPOUND - // In this case, WMAT in output is changed to relative - // weigths. - // - - Double_t* da = CreateDoubleArray(a, TMath::Abs(nlmat)); - Double_t* dz = CreateDoubleArray(z, TMath::Abs(nlmat)); - Double_t* dwmat = CreateDoubleArray(wmat, TMath::Abs(nlmat)); - - Mixture(kmat, name, da, dz, dens, nlmat, dwmat); - for (Int_t i=0; i 0 then wmat contains the proportion by - // weights of each basic material in the mixture. - // - // If nlmat < 0 then WMAT contains the number of atoms - // of a given kind into the molecule of the COMPOUND - // In this case, WMAT in output is changed to relative - // weigths. - // - - if (nlmat < 0) { - nlmat = - nlmat; - Double_t amol = 0; - Int_t i; - for (i=0;iGetListOfMaterials()->GetSize(); - printf("Mixture %s with %i elem: kmat=%i, dens=%g\n", name, nlmat, kmat, dens); - for (Int_t j=0; jMixture(name, a, z, dens, nlmat, wmat, kmat); -} //_____________________________________________________________________________ Int_t TFlukaMCGeometry::GetMedium() const { @@ -337,345 +217,67 @@ Int_t TFlukaMCGeometry::GetMedium() const TGeoNode *node = gGeoManager->GetCurrentNode(); if (!node) imed = gGeoManager->GetTopNode()->GetVolume()->GetMedium()->GetId(); else imed = node->GetVolume()->GetMedium()->GetId(); - printf("GetMedium=%i\n", imed); + if (fDebug) printf("GetMedium=%i\n", imed); return imed; } //_____________________________________________________________________________ -void TFlukaMCGeometry::Medium(Int_t& kmed, const char* name, Int_t nmat, Int_t isvol, - Int_t ifield, Double_t fieldm, Double_t tmaxfd, - Double_t stemax, Double_t deemax, Double_t epsil, - Double_t stmin, Float_t* ubuf, Int_t nbuf) +Int_t TFlukaMCGeometry::GetFlukaMaterial(Int_t imed) const { - // - // kmed tracking medium number assigned - // name tracking medium name - // nmat material number - // isvol sensitive volume flag - // ifield magnetic field - // fieldm max. field value (kilogauss) - // tmaxfd max. angle due to field (deg/step) - // stemax max. step allowed - // deemax max. fraction of energy lost in a step - // epsil tracking precision (cm) - // stmin min. step due to continuous processes (cm) - // - // ifield = 0 if no magnetic field; ifield = -1 if user decision in guswim; - // ifield = 1 if tracking performed with g3rkuta; ifield = 2 if tracking - // performed with g3helix; ifield = 3 if tracking performed with g3helx3. - // - - //printf("Creating mediuma: %s, numed=%d, nmat=%d\n",name,kmed,nmat); - Double_t* dubuf = CreateDoubleArray(ubuf, nbuf); - Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, epsil, - stmin, dubuf, nbuf); - delete [] dubuf; +// Returns FLUKA material index for medium IMED + TGeoMedium *med = (TGeoMedium*)gGeoManager->GetListOfMedia()->At(imed-1); + if (!med) { + Error("GetFlukaMaterial", "MEDIUM %i nor found", imed); + return -1; + } + TGeoMaterial* mat = med->GetMaterial(); + if (!mat->IsUsed()) return -1; + Int_t imatfl = med->GetMaterial()->GetIndex(); + return imatfl; } //_____________________________________________________________________________ -void TFlukaMCGeometry::Medium(Int_t& kmed, const char* name, Int_t nmat, Int_t isvol, - Int_t ifield, Double_t fieldm, Double_t tmaxfd, - Double_t stemax, Double_t deemax, Double_t epsil, - Double_t stmin, Double_t* /*ubuf*/, Int_t /*nbuf*/) +Int_t *TFlukaMCGeometry::GetRegionList(Int_t imed, Int_t &nreg) { - // - // kmed tracking medium number assigned - // name tracking medium name - // nmat material number - // isvol sensitive volume flag - // ifield magnetic field - // fieldm max. field value (kilogauss) - // tmaxfd max. angle due to field (deg/step) - // stemax max. step allowed - // deemax max. fraction of energy lost in a step - // epsil tracking precision (cm) - // stmin min. step due to continuos processes (cm) - // - // ifield = 0 if no magnetic field; ifield = -1 if user decision in guswim; - // ifield = 1 if tracking performed with g3rkuta; ifield = 2 if tracking - // performed with g3helix; ifield = 3 if tracking performed with g3helx3. - // - - kmed = gGeoManager->GetListOfMedia()->GetSize()+3; // !!! in FLUKA they start with 3 - gGeoManager->Medium(name,kmed,nmat, isvol, ifield, fieldm, tmaxfd, stemax,deemax, epsil, stmin); - printf("Medium %s: kmed=%i, nmat=%i, isvol=%i\n", name, kmed, nmat,isvol); +// Get an ordered list of regions matching a given medium number + nreg = 0; + if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1]; + TIter next(gGeoManager->GetListOfUVolumes()); + TGeoVolume *vol; + Int_t imedium, ireg; + while ((vol = (TGeoVolume*)next())) { + TGeoMedium* med; + if ((med = vol->GetMedium()) == 0) continue; + imedium = med->GetId(); + if (imedium == imed) { + ireg = vol->GetNumber(); + fRegionList[nreg++] = ireg; + } + } + return fRegionList; } //_____________________________________________________________________________ -void TFlukaMCGeometry::Matrix(Int_t& krot, Double_t thex, Double_t phix, Double_t they, - Double_t phiy, Double_t thez, Double_t phiz) +Int_t *TFlukaMCGeometry::GetMaterialList(Int_t imat, Int_t &nreg) { - // - // krot rotation matrix number assigned - // theta1 polar angle for axis i - // phi1 azimuthal angle for axis i - // theta2 polar angle for axis ii - // phi2 azimuthal angle for axis ii - // theta3 polar angle for axis iii - // phi3 azimuthal angle for axis iii - // - // it defines the rotation matrix number irot. - // - - krot = gGeoManager->GetListOfMatrices()->GetEntriesFast(); - gGeoManager->Matrix(krot, thex, phix, they, phiy, thez, phiz); - printf("Rotation %i defined\n", krot); +// Get an ordered list of regions matching a given medium number + nreg = 0; + if (!fRegionList) fRegionList = new Int_t[NofVolumes()+1]; + TIter next(gGeoManager->GetListOfUVolumes()); + TGeoVolume *vol; + Int_t imaterial, ireg; + while ((vol = (TGeoVolume*)next())) { + 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::Gsvolu(const char *name, const char *shape, Int_t nmed, - Float_t *upar, Int_t npar) -{ - // - // NAME Volume name - // SHAPE Volume type - // NUMED Tracking medium number - // NPAR Number of shape parameters - // UPAR Vector containing shape parameters - // - // It creates a new volume in the JVOLUM data structure. - // - - Double_t* dupar = CreateDoubleArray(upar, npar); - Int_t id = Gsvolu(name, shape, nmed, dupar, npar); - delete [] dupar; - return id; -} - -//_____________________________________________________________________________ -Int_t TFlukaMCGeometry::Gsvolu(const char *name, const char *shape, Int_t nmed, - Double_t *upar, Int_t npar) -{ - // - // NAME Volume name - // SHAPE Volume type - // NUMED Tracking medium number - // NPAR Number of shape parameters - // UPAR Vector containing shape parameters - // - // It creates a new volume in the JVOLUM data structure. - // - char vname[5]; - Vname(name,vname); - char vshape[5]; - Vname(shape,vshape); - - TGeoVolume* vol = gGeoManager->Volume(vname, shape, nmed, upar, npar); - printf("Volume %s: id=%i shape=%s, nmed=%i\n", vname, vol->GetNumber(), shape, nmed); - return vol->GetNumber(); -} - -//_____________________________________________________________________________ -void TFlukaMCGeometry::Gsdvn(const char *name, const char *mother, Int_t ndiv, - Int_t iaxis) -{ - // - // Create a new volume by dividing an existing one - // - // NAME Volume name - // MOTHER Mother volume name - // NDIV Number of divisions - // IAXIS Axis value - // - // X,Y,Z of CAXIS will be translated to 1,2,3 for IAXIS. - // It divides a previously defined volume. - // - char vname[5]; - Vname(name,vname); - char vmother[5]; - Vname(mother,vmother); - - gGeoManager->Division(vname, vmother, iaxis, ndiv, 0, 0, 0, "n"); - printf("Division %s: mother=%s iaxis=%i ndiv=%i\n", vname, vmother, iaxis, ndiv); -} - -//_____________________________________________________________________________ -void TFlukaMCGeometry::Gsdvn2(const char *name, const char *mother, Int_t ndiv, - Int_t iaxis, Double_t c0i, Int_t numed) -{ - // - // Create a new volume by dividing an existing one - // - // Divides mother into ndiv divisions called name - // along axis iaxis starting at coordinate value c0. - // the new volume created will be medium number numed. - // - char vname[5]; - Vname(name,vname); - char vmother[5]; - Vname(mother,vmother); - - gGeoManager->Division(vname, vmother, iaxis, ndiv, c0i, 0, numed, "nx"); -} -//_____________________________________________________________________________ -void TFlukaMCGeometry::Gsdvt(const char *name, const char *mother, Double_t step, - Int_t iaxis, Int_t numed, Int_t /*ndvmx*/) -{ - // - // Create a new volume by dividing an existing one - // - // Divides MOTHER into divisions called NAME along - // axis IAXIS in steps of STEP. If not exactly divisible - // will make as many as possible and will centre them - // with respect to the mother. Divisions will have medium - // number NUMED. If NUMED is 0, NUMED of MOTHER is taken. - // NDVMX is the expected maximum number of divisions - // (If 0, no protection tests are performed) - // - char vname[5]; - Vname(name,vname); - char vmother[5]; - Vname(mother,vmother); - - gGeoManager->Division(vname, vmother, iaxis, 0, 0, step, numed, "s"); -} - -//_____________________________________________________________________________ -void TFlukaMCGeometry::Gsdvt2(const char *name, const char *mother, Double_t step, - Int_t iaxis, Double_t c0, Int_t numed, Int_t /*ndvmx*/) -{ - // - // Create a new volume by dividing an existing one - // - // Divides MOTHER into divisions called NAME along - // axis IAXIS starting at coordinate value C0 with step - // size STEP. - // The new volume created will have medium number NUMED. - // If NUMED is 0, NUMED of mother is taken. - // NDVMX is the expected maximum number of divisions - // (If 0, no protection tests are performed) - // - char vname[5]; - Vname(name,vname); - char vmother[5]; - Vname(mother,vmother); - - gGeoManager->Division(vname, vmother, iaxis, 0, c0, step, numed, "sx"); -} - -//_____________________________________________________________________________ -void TFlukaMCGeometry::Gsord(const char * /*name*/, Int_t /*iax*/) -{ - // - // Flags volume CHNAME whose contents will have to be ordered - // along axis IAX, by setting the search flag to -IAX - // IAX = 1 X axis - // IAX = 2 Y axis - // IAX = 3 Z axis - // IAX = 4 Rxy (static ordering only -> GTMEDI) - // IAX = 14 Rxy (also dynamic ordering -> GTNEXT) - // IAX = 5 Rxyz (static ordering only -> GTMEDI) - // IAX = 15 Rxyz (also dynamic ordering -> GTNEXT) - // IAX = 6 PHI (PHI=0 => X axis) - // IAX = 7 THETA (THETA=0 => Z axis) - // - - // TBC - keep this function - // nothing to be done for TGeo //xx -} - -//_____________________________________________________________________________ -void TFlukaMCGeometry::Gspos(const char *name, Int_t nr, const char *mother, Double_t x, - Double_t y, Double_t z, Int_t irot, const char *konly) -{ - // - // Position a volume into an existing one - // - // NAME Volume name - // NUMBER Copy number of the volume - // MOTHER Mother volume name - // X X coord. of the volume in mother ref. sys. - // Y Y coord. of the volume in mother ref. sys. - // Z Z coord. of the volume in mother ref. sys. - // IROT Rotation matrix number w.r.t. mother ref. sys. - // ONLY ONLY/MANY flag - // - // It positions a previously defined volume in the mother. - // - - TString only = konly; - only.ToLower(); - Bool_t isOnly = kFALSE; - if (only.Contains("only")) isOnly = kTRUE; - char vname[5]; - Vname(name,vname); - char vmother[5]; - Vname(mother,vmother); - - Double_t *upar=0; - gGeoManager->Node(vname, nr, vmother, x, y, z, irot, isOnly, upar); - printf("Adding daughter %s to %s: cpy=%i irot=%i only=%s\n", vname,vmother,nr,irot,only.Data()); -} -//_____________________________________________________________________________ -void TFlukaMCGeometry::Gsposp(const char *name, Int_t nr, const char *mother, - Double_t x, Double_t y, Double_t z, Int_t irot, - const char *konly, Float_t *upar, Int_t np ) -{ - // - // Place a copy of generic volume NAME with user number - // NR inside MOTHER, with its parameters UPAR(1..NP) - // - - Double_t* dupar = CreateDoubleArray(upar, np); - Gsposp(name, nr, mother, x, y, z, irot, konly, dupar, np); - delete [] dupar; -} - -//_____________________________________________________________________________ -void TFlukaMCGeometry::Gsposp(const char *name, Int_t nr, const char *mother, - Double_t x, Double_t y, Double_t z, Int_t irot, - const char *konly, Double_t *upar, Int_t np ) -{ - // - // Place a copy of generic volume NAME with user number - // NR inside MOTHER, with its parameters UPAR(1..NP) - // - - TString only = konly; - only.ToLower(); - Bool_t isOnly = kFALSE; - if (only.Contains("only")) isOnly = kTRUE; - char vname[5]; - Vname(name,vname); - char vmother[5]; - Vname(mother,vmother); - - gGeoManager->Node(vname,nr,vmother, x,y,z,irot,isOnly,upar,np); - printf("Adding daughter(s) %s to %s: cpy=%i irot=%i only=%s\n", vname,vmother,nr,irot,only.Data()); -} - -//_____________________________________________________________________________ -Int_t TFlukaMCGeometry::VolId(const Text_t *name) const -{ - // - // Return the unique numeric identifier for volume name - // - - Int_t uid = gGeoManager->GetUID(name); - if (uid<0) { - printf("VolId: Volume %s not found\n",name); - return 0; - } - printf("VolId for %s: %i\n", name, uid); - return uid; -} - -//_____________________________________________________________________________ -const char* TFlukaMCGeometry::VolName(Int_t id) const -{ - // - // Return the volume name given the volume identifier - // - TGeoVolume *volume = gGeoManager->GetVolume(id); - if (!volume) { - Error("VolName","volume with id=%d does not exist",id); - return "NULL"; - } - printf("VolName for id=%i: %s\n", id, volume->GetName()); - return volume->GetName(); -} - //_____________________________________________________________________________ Int_t TFlukaMCGeometry::NofVolumes() const { @@ -685,151 +287,286 @@ Int_t TFlukaMCGeometry::NofVolumes() const return gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1; } - -//_____________________________________________________________________________ -Int_t TFlukaMCGeometry::VolId2Mate(Int_t id) const -{ - // - // Return material number for a given volume id - // - TGeoVolume *volume = gGeoManager->GetVolume(id); - if (!volume) { - Error("VolId2Mate","volume with id=%d does not exist",id); - return 0; - } - TGeoMedium *med = volume->GetMedium(); - if (!med) return 0; - printf("VolId2Mate id=%i: idmed=%i\n", id, med->GetId()); - return med->GetId(); -} - -//_____________________________________________________________________________ -Int_t TFlukaMCGeometry::CurrentVolID(Int_t& copyNo) const -{ - // Returns the current volume ID and copy number - if (gGeoManager->IsOutside()) return 0; - TGeoNode *node = gGeoManager->GetCurrentNode(); - copyNo = node->GetNumber(); - Int_t id = node->GetVolume()->GetNumber(); - printf("CurrentVolId(cpy=%i) = %i\n", copyNo, id); - return id; -} - -//_____________________________________________________________________________ -Int_t TFlukaMCGeometry::CurrentVolOffID(Int_t off, Int_t& copyNo) const -{ - // Return the current volume "off" upward in the geometrical tree - // ID and copy number - if (off<0 || off>gGeoManager->GetLevel()) return 0; - if (off==0) return CurrentVolID(copyNo); - TGeoNode *node = gGeoManager->GetMother(off); - if (!node) return 0; - copyNo = node->GetNumber(); - printf("CurrentVolOffId(off=%i,cpy=%i) = %i\n", off,copyNo,node->GetVolume()->GetNumber() ); - return node->GetVolume()->GetNumber(); -} -// FLUKA specific - -//_____________________________________________________________________________ -const char* TFlukaMCGeometry::CurrentVolName() const -{ - // - // Returns the current volume name - // - if (gGeoManager->IsOutside()) return 0; - printf("CurrentVolName : %s\n", gGeoManager->GetCurrentVolume()->GetName()); - return gGeoManager->GetCurrentVolume()->GetName(); -} -//_____________________________________________________________________________ -const char* TFlukaMCGeometry::CurrentVolOffName(Int_t off) const -{ - // - // Return the current volume "off" upward in the geometrical tree - // ID, name and copy number - // if name=0 no name is returned - // - if (off<0 || off>gGeoManager->GetLevel()) return 0; - if (off==0) return CurrentVolName(); - TGeoNode *node = gGeoManager->GetMother(off); - if (!node) return 0; - printf("CurrentVolOffName(off=%i) : %s\n", off,node->GetVolume()->GetName()); - return node->GetVolume()->GetName(); -} //_____________________________________________________________________________ -void TFlukaMCGeometry::Gsatt(const char *name, const char *att, Int_t val) -{ - // - // NAME Volume name - // IOPT Name of the attribute to be set - // IVAL Value to which the attribute is to be set - // see: TFluka::Gsatt - char vname[5]; - Vname(name,vname); - char vatt[5]; - Vname(att,vatt); - gGeoManager->SetVolumeAttribute(vname, vatt, val); -} - -//_____________________________________________________________________________ -void TFlukaMCGeometry::Gdtom(Float_t *xd, Float_t *xm, Int_t iflag) -{ - // - // Computes coordinates XM (Master Reference System - // knowing the coordinates XD (Detector Ref System) - // The local reference system can be initialized by - // - the tracking routines and GDTOM used in GUSTEP - // - a call to GSCMED(NLEVEL,NAMES,NUMBER) - // (inverse routine is GMTOD) - // - // If IFLAG=1 convert coordinates - // IFLAG=2 convert direction cosinus - // - Double_t XM[3], XD[3]; - Int_t i; - for (i=0;i<3;i++) XD[i] = xd[i]; - if (iflag == 1) gGeoManager->LocalToMaster(XD,XM); - else gGeoManager->LocalToMasterVect(XD,XM); - for (i=0;i<3;i++) xm[i]=XM[i]; +TGeoMaterial * TFlukaMCGeometry::GetMakeWrongMaterial(Double_t z) +{ +// Try to replace a wrongly-defined material + static Double_t kz[23] = {7.3, 17.8184, 7.2167, 10.856, 8.875, 8.9, 7.177, + 25.72, 6.2363, 7.1315, 47.7056, 10.6467, 7.8598, 2.10853, 10.6001, 9.1193, + 15.3383, 4.55, 9.6502, 6.4561, 21.7963, 29.8246, 15.4021}; + + Int_t ind; + Double_t dz; + for (ind=0; ind<23; ind++) { + dz = TMath::Abs(z-kz[ind]); + if (dz<1E-4) break; + } + if (ind>22) { + printf("Cannot patch material with Z=%g\n", z); + return 0; + } + TGeoMixture *mix = 0; + TGeoElement *element; + TGeoElementTable *table = gGeoManager->GetElementTable(); + switch (ind) { + case 0: // AIR + 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(1, element, 0.755267); + element = table->GetElement(8); // O + mix->DefineElement(2, element, 0.231781); + element = table->GetElement(18); // AR + mix->DefineElement(3, element, 0.012827); + break; + case 1: //SDD SI CHIP + 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(1, element, 0.039730642); + element = table->GetElement(7); + mix->DefineElement(2, element, 0.001396798); + element = table->GetElement(8); + mix->DefineElement(3, element, 0.01169634); + element = table->GetElement(14); + mix->DefineElement(4, element, 0.844665); + element = table->GetElement(47); + mix->DefineElement(5, element, 0.09814344903); + break; + case 2: // WATER + 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("ITS_CERAMICS", 5, 3.6); + element = table->GetElement(8); + mix->DefineElement(0, element, 0.59956); + element = table->GetElement(13); + mix->DefineElement(1, element, 0.3776); + element = table->GetElement(14); + mix->DefineElement(2, element, 0.00933); + element = table->GetElement(24); + mix->DefineElement(3, element, 0.002); + element = table->GetElement(25); + mix->DefineElement(4, element, 0.0115); + break; + case 4: // EPOXY + mix = new TGeoMixture("MUON_G10FR4", 4, 1.8); + element = table->GetElement(1); + mix->DefineElement(0, element, 0.19); + element = table->GetElement(6); + mix->DefineElement(1, element, 0.18); + element = table->GetElement(8); + mix->DefineElement(2, element, 0.35); + element = table->GetElement(14); + mix->DefineElement(3, element, 0.28); + break; + case 5: // EPOXY + mix = new TGeoMixture("G10FR4", 4, 1.8); + element = table->GetElement(1); + mix->DefineElement(0, element, 0.19); + element = table->GetElement(6); + mix->DefineElement(1, element, 0.18); + element = table->GetElement(8); + mix->DefineElement(2, element, 0.35); + element = table->GetElement(14); + mix->DefineElement(3, element, 0.28); + break; + case 6: // KAPTON + mix = new TGeoMixture("ITS_KAPTON", 4, 1.3); + element = table->GetElement(1); + mix->DefineElement(0, element, 0.026363415); + element = table->GetElement(6); + mix->DefineElement(1, element, 0.6911272); + element = table->GetElement(7); + mix->DefineElement(2, element, 0.073271325); + element = table->GetElement(8); + mix->DefineElement(3, element, 0.209238060); + break; + case 7: // INOX + mix = new TGeoMixture("ITS_INOX", 9, 7.9); + element = table->GetElement(6); + mix->DefineElement(0, element, 0.0003); + element = table->GetElement(14); + mix->DefineElement(1, element, 0.01); + element = table->GetElement(15); + mix->DefineElement(2, element, 0.00045); + element = table->GetElement(16); + mix->DefineElement(3, element, 0.0003); + element = table->GetElement(24); + mix->DefineElement(4, element, 0.17); + element = table->GetElement(25); + mix->DefineElement(5, element, 0.02); + element = table->GetElement(26); + mix->DefineElement(6, element, 0.654); + element = table->GetElement(28); + mix->DefineElement(7, element, 0.12); + element = table->GetElement(42); + mix->DefineElement(8, element, 0.025); + break; + case 8: // ROHACELL + mix = new TGeoMixture("ROHACELL", 4, 0.05); + element = table->GetElement(1); + mix->DefineElement(0, element, 0.07836617); + element = table->GetElement(6); + mix->DefineElement(1, element, 0.64648941); + element = table->GetElement(7); + mix->DefineElement(2, element, 0.08376983); + element = table->GetElement(8); + mix->DefineElement(3, element, 0.19137459); + break; + case 9: // SDD-C-AL + 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(1, element, 0.8176579); + element = table->GetElement(7); + mix->DefineElement(2, element, 0.0093488); + element = table->GetElement(8); + mix->DefineElement(3, element, 0.0503618); + element = table->GetElement(13); + mix->DefineElement(4, element, 0.1); + break; + case 10: // X7R-CAP + 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(1, element, 0.084755042); + element = table->GetElement(28); + mix->DefineElement(2, element, 0.038244751); + element = table->GetElement(29); + mix->DefineElement(3, element, 0.009471271); + element = table->GetElement(50); + mix->DefineElement(4, element, 0.321736471); + element = table->GetElement(56); + mix->DefineElement(5, element, 0.251639432); + element = table->GetElement(82); + mix->DefineElement(6, element, 0.2081768); + break; + case 11: // SDD ruby sph. Al2O3 + 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("ITS_HV-CABLE", 5, 1.6087); + element = table->GetElement(1); + mix->DefineElement(0, element, 0.01983871336); + element = table->GetElement(6); + mix->DefineElement(1, element, 0.520088819984); + element = table->GetElement(7); + mix->DefineElement(2, element, 0.0551367996); + element = table->GetElement(8); + mix->DefineElement(3, element, 0.157399667056); + element = table->GetElement(13); + mix->DefineElement(4, element, 0.247536); + break; + case 13: //SDD LV+signal cable + 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(1, element, 0.21722436468); + element = table->GetElement(7); + mix->DefineElement(2, element, 0.023028867); + element = table->GetElement(8); + mix->DefineElement(3, element, 0.06574077612); + element = table->GetElement(13); + mix->DefineElement(4, element, 0.68572); + break; + case 14: //SDD hybrid microcab + 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(1, element, 0.24281879711); + element = table->GetElement(7); + mix->DefineElement(2, element, 0.02574224025); + element = table->GetElement(8); + mix->DefineElement(3, element, 0.07348667449); + element = table->GetElement(13); + mix->DefineElement(4, element, 0.64869); + break; + case 15: //SDD anode microcab + 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(1, element, 0.392653705471); + element = table->GetElement(7); + mix->DefineElement(2, element, 0.041626868025); + element = table->GetElement(8); + mix->DefineElement(3, element, 0.118832707289); + element = table->GetElement(13); + mix->DefineElement(4, element, 0.431909); + break; + case 16: // inox/alum + mix = new TGeoMixture("ITS_INOX-AL", 5, 3.0705); + element = table->GetElement(13); + mix->DefineElement(0, element, 0.816164); + element = table->GetElement(14); + mix->DefineElement(1, element, 0.000919182); + element = table->GetElement(24); + mix->DefineElement(2, element, 0.0330906); + element = table->GetElement(26); + mix->DefineElement(3, element, 0.131443); + element = table->GetElement(28); + mix->DefineElement(4, element, 0.0183836); + case 17: // MYLAR + mix = new TGeoMixture("TPC_MYLAR", 3, 1.39); + element = table->GetElement(1); + mix->DefineElement(0, element, 0.0416667); + element = table->GetElement(6); + mix->DefineElement(1, element, 0.625); + element = table->GetElement(8); + mix->DefineElement(2, element, 0.333333); + break; + case 18: // SPDBUS(AL+KPT+EPOX) - unknown composition + 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("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("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("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("ITS_SSDEL", 1, 0.68); + element = table->GetElement(16); + mix->DefineElement(0, element, 1.); + z = element->Z(); + break; + } + mix->SetZ(z); + printf("Patched with mixture %s\n", mix->GetName()); + return mix; } -//_____________________________________________________________________________ -void TFlukaMCGeometry::Gdtom(Double_t *xd, Double_t *xm, Int_t iflag) -{ - if (iflag == 1) gGeoManager->LocalToMaster(xd,xm); - else gGeoManager->LocalToMasterVect(xd,xm); -} - -//_____________________________________________________________________________ -void TFlukaMCGeometry::Gmtod(Float_t *xm, Float_t *xd, Int_t iflag) -{ - // - // Computes coordinates XD (in DRS) - // from known coordinates XM in MRS - // The local reference system can be initialized by - // - the tracking routines and GMTOD used in GUSTEP - // - a call to GMEDIA(XM,NUMED,CHECK) - // - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER) - // (inverse routine is GDTOM) - // - // If IFLAG=1 convert coordinates - // IFLAG=2 convert direction cosinus - // - Double_t XM[3], XD[3]; - Int_t i; - for (i=0;i<3;i++) XM[i]=xm[i]; - if (iflag == 1) gGeoManager->MasterToLocal(XM,XD); - else gGeoManager->MasterToLocalVect(XM,XD); - for (i=0;i<3;i++) xd[i] = XD[i]; -} - -//_____________________________________________________________________________ -void TFlukaMCGeometry::Gmtod(Double_t *xm, Double_t *xd, Int_t iflag) -{ - if (iflag == 1) gGeoManager->MasterToLocal(xm,xd); - else gGeoManager->MasterToLocalVect(xm,xd); -} - //_____________________________________________________________________________ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname) { @@ -840,19 +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. - - Int_t zelem[128]; - static char elNames[220] = { - // 1 ============================= 5 ==================================== 10 ===================================== 15 === - 'H','_','H','E','L','I','B','E','B','_','C','_','N','_','O','_','F','_','N','E','N','A','M','G','A','L','S','I','P','_', - 'S','_','C','L','A','R','K','_','C','A','S','C','T','I','V','_','C','R','M','N','F','E','C','O','N','I','C','U','Z','N', - 'G','A','G','E','A','S','S','E','B','R','K','R','R','B','S','R','Y','_','Z','R','N','B','M','O','T','C','R','U','R','H', - 'P','D','A','G','C','D','I','N','S','N','S','B','T','E','I','_','X','E','C','S','B','A','L','A','C','E','P','R','N','D', - 'P','M','S','M','E','U','G','D','T','B','D','Y','H','O','E','R','T','M','Y','B','L','U','H','F','T','A','W','_','R','E', - 'O','S','I','R','P','T','A','U','H','G','T','L','P','B','B','I','P','O','A','T','R','N','F','R','R','A','A','C','T','H', - 'P','A','U','_','N','P','P','U','A','M','C','M','B','K','C','F','E','S','F','M','M','D','N','O','L','R','R','F','D','B', - 'S','G','B','H','H','S','M','T','D','S'}; - memset(zelem, 0, 128*sizeof(Int_t)); + // Materials defined by FLUKA TString sname; gGeoManager->Export("flgeom.root"); if (fname) sname = fname; @@ -865,201 +590,156 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname) } PrintHeader(out, "MATERIALS AND COMPOUNDS"); PrintHeader(out, "MATERIALS"); + Int_t i,j,idmat; + Int_t counttothree, nelem; + Double_t a,z,rho, w; + 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 ... + Int_t nelements = table->GetNelements(); TList *matlist = gGeoManager->GetListOfMaterials(); TIter next(matlist); Int_t nmater = matlist->GetSize(); - Int_t nfmater = 0; - TObjArray *listfluka = new TObjArray(nmater+50); - TObjArray *listflukanames = new TObjArray(nmater+50); - TGeoMaterial *mat, *matorig; + Int_t nfmater = 0; + Int_t newind = 26; // here non predefined materials start + + TGeoMaterial *mat; TGeoMixture *mix = 0; TString matname; - TObjString *objstr, *objstrother; - Int_t i,j,k,idmat; - Bool_t done; - Int_t nelem, nidmat; - Double_t amat,zmat,rhomat; - Double_t zel, ael, wel, rho; - char elname[8] = {' ',' ','_', 'E','L','E','M','\0'}; - char digit[3]; - Bool_t found = kFALSE; + TObjString *objstr; + // Create all needed elements + for (Int_t i=1; iGetElement(i); + // skip elements which are not defined + if (!element->IsUsed() && !element->IsDefined()) continue; + matname = element->GetTitle(); + ToFlukaString(matname); + rho = 0.999; + + mat = new TGeoMaterial(matname, element->A(), element->Z(), rho); + // 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()); + fMatNames->Add(objstr); + nfmater++; + } - printf("Creating materials and compounds\n"); + fIndmat = nfmater; + // Adjust material names and add them to FLUKA list for (i=0; iAt(i); - if (mat->GetZ()<1E-1) { + if (!mat->IsUsed()) continue; + z = mat->GetZ(); + a = mat->GetA(); + rho = mat->GetDensity(); + if (mat->GetZ()<0.001) { mat->SetIndex(2); // vacuum, built-in inside FLUKA continue; - } -// printf("material: %s index=%i: Z=%f A=%f rho=%f\n", mat->GetName(), mat->GetIndex(),mat->GetZ(),mat->GetA(),mat->GetDensity()); - matorig = gGeoManager->FindDuplicateMaterial(mat); - if (matorig) { - idmat = matorig->GetIndex(); - mat->SetIndex(idmat); -// printf(" -> found a duplicate: %s with index %i\n", matorig->GetName(), idmat); - matorig = 0; - } else { -// printf(" Adding to temp list with index %i\n", nfmater+3); - listfluka->Add(mat); - mat->SetIndex(nfmater+3); - matorig = mat; - objstr = new TObjString(mat->GetName()); - listflukanames->Add(objstr); - nfmater++; - // look if name is existing - nidmat = 0; - matname = objstr->GetString(); - ToFlukaString(matname); - objstr->SetString(matname.Data()); - done = kFALSE; - while (!done) { - if (nfmater == 1) break; - for (j=0; jAt(j); - if (objstr->IsEqual(objstrother)) { - // we have to change the name - if (nidmat>98) { - Error("CreateFlukaMatFile", "too many materials having same name"); - return; - } - nidmat++; - k = matname.Index(" "); - if (k<0 || k>6) k=6; - if (nidmat>9) { - sprintf(digit, "%d", nidmat); - } else { - digit[0] = '0'; - sprintf(&digit[1], "%d", nidmat); - } - matname.Insert(k,digit); - matname.Remove(8); - objstr->SetString(matname.Data()); - break; - } - if (j == nfmater-2) { - done = kTRUE; - break; - } - } - } -// printf(" newmat name: %s\n", matname.Data()); - } - // now we have unique materials with unique names in the lists - - if (matorig && matorig->IsMixture()) { - // create dummy materials for elements - rho = 0.999; - mix = (TGeoMixture*)matorig; - nelem = mix->GetNelements(); -// printf(" material is a MIXTURE with %i elements:\n", nelem); - for (j=0; jGetZmixt())[j]; - ael = (mix->GetAmixt())[j]; -// printf(" Zelem[%i] = %g\n",j,zel); - if ((zel-Int_t(zel))>0.01) { - TGeoMaterial *mat1; - for (Int_t imat=0; imatAt(imat); - if (TMath::Abs(mat1->GetZ()-zel)>1E-4) continue; - if (TMath::Abs(mat1->GetA()-ael)>1E-4) continue; - found = kTRUE; - break; - } - if (!found) Warning("CreateFlukaMatFile", "element with Z=%f\n", zel); - } - if (!zelem[Int_t(zel)] && !found) { - // write fluka element - memcpy(elname, &elNames[2*Int_t(zel-1)], 2); - zelem[Int_t(zel)] = 1; - mat = new TGeoMaterial(elname, ael, zel, rho); - mat->SetIndex(nfmater+3); -// printf(" element not in list: new material %s at index=%i, Z=%g, A=%g, dummyrho=%g\n", -// elname,nfmater+3,zel,ael,rho); - listfluka->Add(mat); - objstr = new TObjString(elname); - listflukanames->Add(objstr); - nfmater++; - } - } - } - } - // now dump materials in the file -// printf("DUMPING %i materials\n", nfmater); + } + matname = mat->GetName(); + FlukaMatName(matname); + + mat->SetIndex(newind++); + objstr = new TObjString(matname.Data()); + fMatList->Add(mat); + fMatNames->Add(objstr); + nfmater++; + } + + // Dump all elements with MATERIAL cards for (i=0; iAt(i); + mat = (TGeoMaterial*)fMatList->At(i); +// mat->SetUsed(kFALSE); + mix = 0; out << setw(10) << "MATERIAL "; out.setf(static_cast(0),std::ios::floatfield); - matname = mat->GetName(); - ToFlukaString(matname); - zmat = mat->GetZ(); - if (zmat-Int_t(zmat)>0.01) { - if (zmat-Int_t(zmat)>0.5) zmat = Int_t(zmat)+1.; - else zmat = Int_t(zmat); - } - amat = mat->GetA(); - rhomat = mat->GetDensity(); - // write material card + objstr = (TObjString*)fMatNames->At(i); + matname = objstr->GetString(); + z = mat->GetZ(); + a = mat->GetA(); + rho = mat->GetDensity(); if (mat->IsMixture()) { out << setw(10) << " "; out << setw(10) << " "; mix = (TGeoMixture*)mat; } else { - out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << zmat; - out << setw(10) << setprecision(3) << amat; + out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << z; + out << setw(10) << setprecision(3) << a; } out.setf(static_cast(0),std::ios::floatfield); - out << setw(10) << setiosflags(ios::scientific) << setprecision(3) << rhomat; + out << setw(10) << setiosflags(ios::scientific) << setprecision(3) << rho; out.setf(static_cast(0),std::ios::floatfield); - out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(i+3); + out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex()); out << setw(10) << " "; out << setw(10) << " "; out << setw(8) << matname.Data() << endl; - } - // write mixture header - PrintHeader(out, "COMPOUNDS"); - Int_t counttothree; - TGeoMaterial *element; - for (i=0; iAt(i); - if (!mat->IsMixture()) continue; - mix = (TGeoMixture*)mat; + if (!mix) { + // add LOW-MAT card for NEON to associate with ARGON neutron xsec + if (z==10) { + out << setw(10) << "LOW-MAT "; + out.setf(static_cast(0),std::ios::floatfield); + out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex()); + out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << 18.; + out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << -2.; + out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << 293.; + out << setw(10) << " "; + out << setw(10) << " "; +// out << setw(8) << matname.Data() << endl; + out << setw(8) << " " << endl; + } + else { + element = table->GetElement((int)z); + TString elename = element->GetTitle(); + ToFlukaString(elename); + if( matname.CompareTo( elename ) != 0 ) { + out << setw(10) << "LOW-MAT "; + out.setf(static_cast(0),std::ios::floatfield); + out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex()); + out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << z; + out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << " "; + out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << " "; + out << setw(10) << " "; + out << setw(10) << " "; + // missing material at Low Energy Cross Section Table + if( (int)z==10 || (int)z==21 || (int)z==34 || (int)z==37 || (int)z==39 || (int)z==44 || + (int)z==45 || (int)z==46 || (int)z==52 || (int)z==57 || (int)z==59 || (int)z==60 || + (int)z==61 || (int)z==65 || (int)z==66 || (int)z==67 || (int)z==68 || (int)z==69 || + (int)z==70 || (int)z==71 || (int)z==72 || (int)z==76 || (int)z==77 || (int)z==78 || + (int)z==81 || (int)z==84 || (int)z==85 || (int)z==86 || (int)z==87 || (int)z==88 || + (int)z==89 || (int)z==91 ) + out << setw(8) << "UNKNOWN " << endl; + else + out << setw(8) << elename.Data() << endl; + // out << setw(8) << " " << endl; + } + } + continue; + } counttothree = 0; out << setw(10) << "COMPOUND "; nelem = mix->GetNelements(); - objstr = (TObjString*)listflukanames->At(i); + objstr = (TObjString*)fMatNames->At(i); matname = objstr->GetString(); -// printf("MIXTURE %s with index %i having %i elements\n", matname.Data(), mat->GetIndex(),nelem); for (j=0; jGetZmixt())[j], -// (mix->GetAmixt())[j],(mix->GetWmixt())[j]); - wel = (mix->GetWmixt())[j]; - zel = (mix->GetZmixt())[j]; - ael = (mix->GetAmixt())[j]; - if (zel-Int_t(zel)>0.01) { - // loop the temporary list - element = 0; - TGeoMaterial *mat1; - for (Int_t imat=0; imatAt(imat); - if (TMath::Abs(mat1->GetZ()-zel)>1E-4) continue; - if (TMath::Abs(mat1->GetA()-ael)>1E-4) continue; - element = mat1; - break; - } - } else { - memcpy(elname, &elNames[2*Int_t(zel-1)], 2); - element = (TGeoMaterial*)listfluka->FindObject(elname); - } - if (!element) { - Error("CreateFlukaMatFile", "Element Z=%g %s not found", zel, elname); - return; - } - idmat = element->GetIndex(); -// printf("element %s , index=%i\n", element->GetName(), idmat); + w = (mix->GetWmixt())[j]; + if (w<0.00001) w=0.00001; + z = (mix->GetZmixt())[j]; + a = (mix->GetAmixt())[j]; + idmat = GetElementIndex(Int_t(z)); + if (!idmat) Error("CreateFlukaMatFile", "element with Z=%f not found", z); out.setf(static_cast(0),std::ios::floatfield); - out << setw(10) << setiosflags(ios::fixed) << setprecision(6) << -wel; + out << setw(10) << setiosflags(ios::fixed) << setprecision(6) << -w; out.setf(static_cast(0),std::ios::floatfield); out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(idmat); counttothree++; @@ -1068,53 +748,449 @@ void TFlukaMCGeometry::CreateFlukaMatFile(const char *fname) out << endl; if ( (j+1) != nelem) out << setw(10) << "COMPOUND "; counttothree = 0; - } - } - //Unless we have 3, 6, 9... submaterials we need to put some empty - //space and the compound name + } + } if (nelem%3) { for (j=0; j<(3-(nelem%3)); j++) out << setw(10) << " " << setw(10) << " "; out << matname.Data(); out << endl; - } - } - - // Now print the list of regions (volumes in TGeo) + } + } Int_t nvols = gGeoManager->GetListOfUVolumes()->GetEntriesFast()-1; TGeoVolume *vol; -/* - PrintHeader(out, "TGEO VOLUMES"); - for (i=1; i<=nvols; i++) { - vol = gGeoManager->GetVolume(i); - out.setf(std::ios::left, std::ios::adjustfield); - out << setw(10) << i; - out << setw(20) << vol->GetName() << endl; - } -*/ // Now print the material assignments - Double_t flagfield; + Double_t flagfield = 0.; + printf("#############################################################\n"); + if (gFluka->IsFieldEnabled()) { + flagfield = 1.; + printf("Magnetic field enabled\n"); + } else printf("Magnetic field disabled\n"); + printf("#############################################################\n"); + PrintHeader(out, "TGEO MATERIAL ASSIGNMENTS"); for (i=1; i<=nvols; i++) { + TGeoMedium* med; vol = gGeoManager->GetVolume(i); - mat = vol->GetMedium()->GetMaterial(); + if ((med = vol->GetMedium()) == 0) continue; + mat = med->GetMaterial(); idmat = mat->GetIndex(); - flagfield = (vol->GetField())?1.:0.; + for (Int_t j=0; jAt(j); + if (mat->GetIndex() == idmat) mat->SetUsed(kTRUE); + } + + 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(0),std::ios::floatfield); out << setw(10) << setiosflags(ios::fixed) << Double_t(idmat); out << setw(10) << setiosflags(ios::fixed) << Double_t(i); out << setw(10) << "0.0"; - out << setw(10) << setiosflags(ios::fixed) << flagfield; + out << setw(10) << "0.0"; + out << setw(10) << setiosflags(ios::fixed) << hasfield; + out << setw(10) << "0.0"; out << endl; } - delete listfluka; - listflukanames->Delete(); - delete listflukanames; + // 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; } +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()); +} + +//_____________________________________________________________________________ +void TFlukaMCGeometry::WritePegFile(Int_t imat, Int_t *NoStern, Int_t *ElemError, + Int_t *MixError, Int_t *countGas) const +{ + // Write the .peg file for one material + + TGeoMaterial *mat = (TGeoMaterial*)fMatList->At(imat); + TString name = ((TObjString*)fMatNames->At(imat))->GetString(); + TString line; + char number[20]; + TGeoElement *elem = mat->GetElement(); + name = name.Strip(); + TString sname = "mat"; + sprintf(number, "%d", imat); + sname.Append(number); + sname.Append(".peg"); + sname.Prepend("peg/"); + ofstream out; + out.open(sname.Data(), ios::out); + if (!out.good()) return; + Double_t dens = mat->GetDensity(); + TGeoMixture *mix = 0; + Int_t nel = 1; + Int_t i; + if (mat->IsMixture()) { + mix = (TGeoMixture*)mat; + nel = mix->GetNelements(); + } + + if (nel==1) { + cout << "( Element ) " << name << " Z=" << mat->GetZ() << " Rho " << mat->GetDensity() << endl; + + Double_t zel = mat->GetZ(); + if( (zel-Int_t(zel))>0.001 || zel < 1 ) { + cout << " ERROR: A Element with not integer Z=" << zel << endl; + cout << endl; + (*ElemError)++; + return; + } + + out << "ELEM" << endl; + out << " &INP IRAYL=1, RHO=" << dens << ", " << endl; + + // check for the Sternheimer parameters + Double_t *issb_parm = GetISSB( mat->GetDensity(), 1, &zel, 0 ); + if( issb_parm[0] > 0 && issb_parm[1] > 0 ) { + cout << "Sternheimer parameters found" << endl; + out << ", ISSB=1, IEV=" << issb_parm[0] << ", CBAR=" << issb_parm[1] + << ", X0=" << issb_parm[2] << "," << endl; + out << "X1=" <GetName() << endl; + + } + else { + + cout << "( Mixture ) " << name << " Rho " << dens << " nElem " << nel << endl; + + Double_t *zt = new Double_t[nel]; + Double_t *wt = new Double_t[nel]; + for (int j=0; jGetZmixt())[j]; + wt[j] = (mix->GetWmixt())[j]; + if( (zt[j]-Int_t(zt[j])) > 0.001 || zt[j] < 1 ) { + cout << "ERROR Mixture " << name << " with an element with not integer Z=" << zt[j] << endl; + cout << endl; + (*MixError)++; + // just continue since the mixtures are not patch, + // but the final release should include the return + // return; + } + } + Double_t *issb_parm = GetISSB( mat->GetDensity(), nel, zt, wt ); + out << "MIXT" << endl; + out << " &INP IRAYL=1, NE=" << nel << ", RHOZ=" << wt[0] << ","; + line = Form(" &INP IRAYL=1, NE=%d RHOZ=%g", nel, wt[0]); + for(int j=1; j 60 ) { out << endl; line = ""; } + } + out << " RHO=" << mat->GetDensity() << ", "; + line += Form(" RHO=%g, ", mat->GetDensity()); + if( line.Length() > 60 ) { out << endl; line = ""; } + + if( issb_parm[0] > 0 && issb_parm[1] > 0 ) { + cout << "Sternheimer parameters found" << endl; + out << " ISSB=1, IEV=" << issb_parm[0] << ","; + line += Form(" ISSB=1, IEV=%g,", issb_parm[0]); + if( line.Length() > 60 ) { out << endl; line = ""; } + out << " CBAR=" << issb_parm[1] << ","; + line += Form(" CBAR=%g,",issb_parm[1]); + if( line.Length() > 60 ) { out << endl; line = ""; } + out << " X0=" << issb_parm[2] << ","; + line += Form(" X0=%g,", issb_parm[2]); + if( line.Length() > 60 ) { out << endl; line = ""; } + out << " X1=" << issb_parm[3] << ","; + line += Form(" X1=%g,", issb_parm[3]); + if( line.Length() > 60 ) { out << endl; line = ""; } + out << " AFACT="<< issb_parm[4] << ","; + line += Form(" AFACT=%g,", issb_parm[4]); + if( line.Length() > 60 ) { out << endl; line = ""; } + out << " SK=" << issb_parm[5] << ","; + line += Form(" SK=%g,", issb_parm[5]); + if( line.Length() > 60 ) { out << endl; line = ""; } + } + else { + cout << "Sternheimer parameters not found" << endl; + (*NoStern)++; + } + + if (dens<0.01){ + (*countGas)++; + out << " GASP=1." << endl; + } + + out << " &END" << endl; + out << name.Data() << endl; + for (i=0; iGetElement(i); + line = elem->GetName(); + if (line.Length()==1) line.Append(" "); + out << line.Data() << " "; + } + out << endl; + + delete [] zt; + 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=" << 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; + out << " $INP $END" << endl; + out << "TEST" << endl; + out << " $INP $END" << endl; + out.close(); +} + +Double_t * TFlukaMCGeometry::GetISSB(Double_t rho, Int_t nElem, Double_t *zelem, Double_t *welem ) const +{ + // Read the density effect parameters + // from R.M. Sternheimer et al. Atomic Data + // and Nuclear Data Tables, Vol. 30 No. 2 + // + // return the parameters if the element/mixture match with one of the list + // otherwise returns the parameters set to 0 + + struct sternheimerData { + TString longname; // element/mixture name + Int_t nelems; // number of constituents N + Int_t Z[20]; //[nelems] Z + Double_t wt[20]; //[nelems] weight fraction + Double_t density; // g/cm3 + Double_t iev; // Average Ion potential (eV) + // **** 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; + TString formula; + Int_t num; + char state; + + static Double_t parameters[7]; + memset( parameters, 0, sizeof(Double_t) ); + + static sternheimerData sternDataArray[300]; + static Bool_t isFileRead = kFALSE; + + // Read the data file if is needed + if( isFileRead == kFALSE ) { + TString sSternheimerInp = getenv("ALICE_ROOT"); + sSternheimerInp +="/TFluka/input/Sternheimer.data"; + + ifstream in(sSternheimerInp); + char line[100]; + in.getline(line, 100); + in.getline(line, 100); + in.getline(line, 100); + in.getline(line, 100); + in.getline(line, 100); + in.getline(line, 100); + + + Int_t is = 0; + while( !in.eof() ) { + in >> shortname >> num >> sternDataArray[is].nelems + >> sternDataArray[is].longname >> formula >> state; + if( in.eof() ) break; + for(int i=0; i> sternDataArray[is].Z[i] >> sternDataArray[is].wt[i]; + } + in >> sternDataArray[is].density; + in >> sternDataArray[is].iev; + in >> sternDataArray[is].cbar; + in >> sternDataArray[is].x0; + in >> sternDataArray[is].x1; + in >> sternDataArray[is].afact; + in >> sternDataArray[is].sk; + if( sternDataArray[is].nelems == 1 ) in >> sternDataArray[is].delta0; + is++; + } + isFileRead = kTRUE; + in.close(); + } + + Int_t is = 0; + while( is < 280 ) { + + // check for elements + if( sternDataArray[is].nelems == 1 && nElem == 1 + && sternDataArray[is].Z[0] == Int_t(*zelem) + && TMath::Abs( (sternDataArray[is].density - rho)/sternDataArray[is].density ) < 0.1 ) { + cout << sternDataArray[is].longname << " #elems:" << sternDataArray[is].nelems << " Rho:" + << sternDataArray[is].density << endl; + cout << sternDataArray[is].iev << " " + << sternDataArray[is].cbar << " " + << sternDataArray[is].x0 << " " + << sternDataArray[is].x1 << " " + << sternDataArray[is].afact << " " + << sternDataArray[is].sk << " " + << sternDataArray[is].delta0 << endl; + + parameters[0] = sternDataArray[is].iev; + parameters[1] = sternDataArray[is].cbar; + parameters[2] = sternDataArray[is].x0; + parameters[3] = sternDataArray[is].x1; + parameters[4] = sternDataArray[is].afact; + parameters[5] = sternDataArray[is].sk; + parameters[6] = sternDataArray[is].delta0; + return parameters; + } + + // check for mixture + int nmatch = 0; + if( sternDataArray[is].nelems > 1 && sternDataArray[is].nelems == nElem ) { + for(int j=0; j 1 && + TMath::Abs( (sternDataArray[is].density - rho)/sternDataArray[is].density ) < 0.1 + && nmatch == sternDataArray[is].nelems ) { + cout << sternDataArray[is].longname << " #elem:" << sternDataArray[is].nelems << " Rho:" + << sternDataArray[is].density << endl; + cout << sternDataArray[is].iev << " " + << sternDataArray[is].cbar << " " + << sternDataArray[is].x0 << " " + << sternDataArray[is].x1 << " " + << sternDataArray[is].afact << " " + << sternDataArray[is].sk << " " + << sternDataArray[is].delta0 << endl; + + parameters[0] = sternDataArray[is].iev; + parameters[1] = sternDataArray[is].cbar; + parameters[2] = sternDataArray[is].x0; + parameters[3] = sternDataArray[is].x1; + parameters[4] = sternDataArray[is].afact; + parameters[5] = sternDataArray[is].sk; + parameters[6] = 0; + return parameters; + } + is++; + } + return parameters; +} + //_____________________________________________________________________________ void TFlukaMCGeometry::PrintHeader(ofstream &out, const char *text) const { @@ -1135,6 +1211,64 @@ Int_t TFlukaMCGeometry::RegionId() const return gGeoManager->GetCurrentNode()->GetUniqueID(); } +//_____________________________________________________________________________ +Int_t TFlukaMCGeometry::GetElementIndex(Int_t z) const +{ +// Get index of a material having a given Z element. + TIter next(fMatList); + TGeoMaterial *mat; + Int_t index = 0; + while ((mat=(TGeoMaterial*)next())) { + if (mat->IsMixture()) continue; + if (mat->GetElement()->Z() == z) return mat->GetIndex(); + } + return index; +} + +//_____________________________________________________________________________ +void TFlukaMCGeometry::SetMreg(Int_t mreg, Int_t lttc) +{ +// Update if needed next history; +// 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; + } + fCurrentRegion = mreg; + fCurrentLattice = lttc; + + Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1; + if (crtlttc == lttc) return; + gGeoManager->CdNode(lttc-1); + while (gGeoManager->GetCurrentVolume()->IsAssembly()) gGeoManager->CdUp(); +} + +//_____________________________________________________________________________ +void TFlukaMCGeometry::SetCurrentRegion(Int_t mreg, Int_t latt) +{ +// Set index/history for next entered region + fCurrentRegion = mreg; + fCurrentLattice = latt; +} + +//_____________________________________________________________________________ +void TFlukaMCGeometry::SetNextRegion(Int_t mreg, Int_t latt) +{ +// Set index/history for next entered region + fNextRegion = mreg; + fNextLattice = latt; +} + //_____________________________________________________________________________ void TFlukaMCGeometry::ToFlukaString(TString &str) const { @@ -1153,6 +1287,51 @@ void TFlukaMCGeometry::ToFlukaString(TString &str) const if (str(pos)==' ') str.Replace(pos,1,"_",1); return; } + +//_____________________________________________________________________________ +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; + for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break; + if (ilast>5) ilast = 5; + char number[3]; + TIter next(fMatNames); + TObjString *objstr; + TString matname; + Int_t index = 0; + while ((objstr=(TObjString*)next())) { + matname = objstr->GetString(); + if (matname == str) { + index++; + if (index<10) { + number[0] = '0'; + sprintf(&number[1], "%d", index); + } else if (index<100) { + sprintf(number, "%d", index); + } else { + Error("FlukaMatName", "Too many materials %s", str.Data()); + return; + } + str.Replace(ilast+1, 2, number); + str.Remove(8); + } + } +} + //______________________________________________________________________________ void TFlukaMCGeometry::Vname(const char *name, char *vname) const { @@ -1168,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*/) { @@ -1180,145 +1391,230 @@ Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/) // card in fluka input), returns 1 if user wants Fluka always to // use DNEAR (in this case, be sure that GEANT4 DNEAR is unique, // coming from all directions!!!) - printf("========== Dummy IDNRWR\n"); + if (gMCGeom->IsDebugging()) printf("========== Dummy IDNRWR\n"); return 0; } //_____________________________________________________________________________ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz, - Double_t *pV, Int_t &oldReg , const Int_t &oldLttc, Double_t & propStep, - Int_t & /*nascFlag*/, Double_t &retStep, Int_t &newReg, - Double_t &saf, Int_t &newLttc, Int_t <tcFlag, + Double_t *pV, Int_t &oldReg , const Int_t &oldLttc, Double_t &propStep, + Int_t &/*nascFlag*/, Double_t &retStep, Int_t &newReg, + Double_t &saf, Int_t &newLttc, Int_t <tcFlag, Double_t *sLt, Int_t *jrLt) + { -// from FLUGG: -// Wrapper for geometry tracking: returns approved step of -// particle and all variables that fluka G1 computes. + // Initialize 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; + + 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; + 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; + jrLt[lttcFlag] = olttc; + sLt[lttcFlag] = propStep; + jrLt[lttcFlag+1] = -1; + sLt[lttcFlag+1] = 0.; + newReg = oldReg; + newLttc = olttc; + Bool_t crossedDummy = (oldReg == gFluka->GetDummyRegion())?kTRUE:kFALSE; + Int_t curLttc, curReg; + 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); + + curLttc = gGeoManager->GetCurrentNodeId()+1; + curReg = gGeoManager->GetCurrentVolume()->GetNumber(); + if (olttc != curLttc) { + // FLUKA crossed the boundary : we trust that the given point is really there, + // so we just update TGeo state + gGeoManager->CdNode(olttc-1); + curLttc = gGeoManager->GetCurrentNodeId()+1; + curReg = gGeoManager->GetCurrentVolume()->GetNumber(); + } + // Now the current TGeo state reflects the FLUKA state - // Initialize current point/direction - printf("========== Inside G1WR\n"); - printf(" point/dir:(%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pV[0],pV[1],pV[2]); gGeoManager->SetCurrentPoint(pSx, pSy, pSz); gGeoManager->SetCurrentDirection(pV); - printf(" oldReg=%i oldLttc=%i pstep=%f\n",oldReg, oldLttc, propStep); - Int_t curLttc = gGeoManager->GetCurrentNodeId()+1; - Int_t curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); - printf(" curReg=%i curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath()); - Bool_t regsame = (curreg==oldReg)?kTRUE:kFALSE; - if (!regsame) printf(" REGIONS DOES NOT MATCH\n"); - if (oldLttc != curLttc) { - printf(" HISTORIES DOES NOT MATCH\n"); - gGeoManager->CdNode(oldLttc-1); - curLttc = gGeoManager->GetCurrentNodeId()+1; - curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); - printf(" re-initialized point: curReg=%i curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath()); - } - lttcFlag = 0; - sLt[lttcFlag] = 0.; - jrLt[lttcFlag] = curLttc; - // now 'oldregion' contains the real region, matching or not the old history - - // Compute geometry step/safety within physical step limit -// newReg = oldregion; - Double_t *point = gGeoManager->GetCurrentPoint(); - Double_t *dir = gGeoManager->GetCurrentDirection(); - Double_t steptot = 0.; - Double_t snext = 0.; - Int_t istep = 0; - Bool_t done = kFALSE; - Int_t i; - while (!done) { - gGeoManager->FindNextBoundary(-propStep); - snext = gGeoManager->GetStep(); - printf(" FindNextBoundary(%g) snext=%g\n", propStep, snext); - if (steptot == 0) { - saf = gGeoManager->GetSafeDistance(); - printf(" Safety: %g\n", saf); - } - sLt[lttcFlag] = propStep; - jrLt[lttcFlag] = gGeoManager->GetCurrentNodeId()+1; - lttcFlag++; //1 - newReg = curreg; - newLttc = oldLttc; - if (snextFindNode(); - istep = 0; - printf(" boundary: step made %g\n", snext); - while (gGeoManager->IsSameLocation() && steptot1E3) { - printf("Geometry error: could not cross boundary after extra 10 microns\n"); - return; - } - for (i=0;i<3;i++) point[i]+=1E-6*dir[i]; - gGeoManager->FindNode(); - sLt[lttcFlag] += 1E-6; - retStep = sLt[lttcFlag]; - steptot += 1E-6; - istep++; - } - lttcFlag++; //1 - if (steptot>propStep) {printf("Error\n");return;} - // we managed to cross the boundary -> in which region - newReg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); - sLt[lttcFlag] = propStep-steptot; - newLttc = (gGeoManager->IsOutside())?999999999:gGeoManager->GetCurrentNodeId()+1; - jrLt[lttcFlag] = newLttc; - if (!gGeoManager->IsOutside()) { - lttcFlag++; //2 - printf(" ENTERED region %i, newLttc=%i in: %s\n", newReg,newLttc,gGeoManager->GetPath()); - } else printf(" EXIT GEOMETRY: BLKHOLE reg=%i\n", newReg); - } - // no boundary within proposed step - lttcFlag--; - done = kTRUE; + 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()); } - printf("=> newReg=%i newLttc=%i lttcFlag=%i\n", newReg, newLttc, lttcFlag); - printf("=> snext=%g safe=%g\n", snext, saf); - for (Int_t i=0; iIsOutside()) { - gGeoManager->SetOutside(kFALSE); - gGeoManager->CdTop(); - } - gGeoManager->CdNode(oldLttc-1); +*/ + 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; } - printf("<= G1WR (in: %s)\n", gGeoManager->GetPath()); + + 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; + } + 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; + 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) || (oldReg!=newReg && olttc==newLttc) ) && pid!=-1) { + // Virtual boundary between replicants + newReg = gFluka->GetDummyRegion(); + newLttc = TFlukaMCGeometry::kLttcVirtual; + 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.; + 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("========== Dummy G1RTWR\n"); } //_____________________________________________________________________________ -void conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/) +void conhwr(Int_t & intHist, Int_t & incrCount) { - 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) { - 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 (intHist==0) gGeoManager->CdTop(); else gGeoManager->CdNode(intHist-1); - printf(" --- current path: %s\n", gGeoManager->GetPath()); - printf("<= INIHWR\n"); + if (gMCGeom->IsDebugging()) { + printf(" --- current path: %s\n", gGeoManager->GetPath()); + printf("<= INIHWR\n"); + } } //_____________________________________________________________________________ @@ -1328,112 +1624,121 @@ void jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/ // Geometry initialization wrapper called by FLUKAM. Provides to FLUKA the // number of regions (volumes in TGeo) // build application geometry - printf("========== Inside JOMIWR\n"); - flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast(); - printf("<= JOMIWR: last region=%i\n", flukaReg); + if (gMCGeom->IsDebugging()) printf("========== Inside JOMIWR\n"); + 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) { - printf("========== Inside LKDBWR (%f, %f, %f)\n",pSx, pSy, pSz); -// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); - printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); - TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz); - if (gGeoManager->IsOutside()) { - newReg = mcgeom->NofVolumes()+1; - newLttc = gGeoManager->GetCurrentNodeId(); - } - newReg = node->GetVolume()->GetNumber(); - newLttc = gGeoManager->GetCurrentNodeId()+1; - flagErr = newReg; - printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); - printf("<= LKDBWR\n"); + 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,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) { - printf("========== Inside LKFXWR (%f, %f, %f)\n",pSx, pSy, pSz); -// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); - printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); - TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz); - if (gGeoManager->IsOutside()) { - newReg = mcgeom->NofVolumes()+1; - newLttc = gGeoManager->GetCurrentNodeId(); - } - newReg = node->GetVolume()->GetNumber(); - newLttc = gGeoManager->GetCurrentNodeId()+1; - flagErr = newReg; - printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); - printf("<= LKFXWR\n"); + 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,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) { - printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz); -// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); - printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); - TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz); - if (gGeoManager->IsOutside()) { - newReg = mcgeom->NofVolumes()+1; - newLttc = gGeoManager->GetCurrentNodeId(); - } - newReg = node->GetVolume()->GetNumber(); - newLttc = gGeoManager->GetCurrentNodeId()+1; - flagErr = newReg; - printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc); - printf("<= LKMGWR\n"); + 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,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) { - printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz); -// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]); - printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc); - TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz); + if (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); + } + 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 = mcgeom->NofVolumes()+1; - newLttc = gGeoManager->GetCurrentNodeId(); + 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; - flagErr = newReg; - printf(" out: newReg=%i newLttc=%i in %s\n", newReg, newLttc, gGeoManager->GetPath()); - printf("<= LKWR\n"); + newLttc = gGeoManager->GetCurrentNodeId()+1; + 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(" 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) { - printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz); - printf(" oldReg=%i, newReg=%i\n", oldReg,newReg); - Int_t curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); - Int_t curLttc = gGeoManager->GetCurrentNodeId()+1; - 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) { - printf(" REGIONS DOEN NOT MATCH\n"); - gGeoManager->FindNode(); - curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); - curLttc = gGeoManager->GetCurrentNodeId()+1; - printf(" re-initialized point: curReg=%i curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath()); - } + 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(pVx, pVy, pVz); Double_t *dnorm = gGeoManager->FindNormalFast(); flagErr = 0; if (!dnorm) { @@ -1442,22 +1747,24 @@ void nrmlwr(Double_t &pSx, Double_t &pSy, Double_t &pSz, norml[0] = -pVx; norml[1] = -pVy; norml[2] = -pVz; + } 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"); } - norml[0] = -dnorm[0]; - norml[1] = -dnorm[1]; - norml[2] = -dnorm[2]; - printf(" normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]); - curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber(); - curLttc = gGeoManager->GetCurrentNodeId()+1; - printf(" final location: curReg=%i, curLttc=%i in %s\n", curreg,curLttc,gGeoManager->GetPath()); - printf("<= NRMLWR\n"); + } //_____________________________________________________________________________ void rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/, Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/) { - printf("=> Dummy RGRPWR\n"); + if (gMCGeom->IsDebugging()) printf("=> Dummy RGRPWR\n"); } //_____________________________________________________________________________ @@ -1476,10 +1783,10 @@ 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. - 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; - printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath()); + if (gMCGeom->IsDebugging()) printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath()); return histInt; }