**************************************************************************/
// $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 "TCallf77.h"
#include "TFluka.h"
#include "TFlukaMCGeometry.h"
+#include "TFlukaConfigOption.h"
#include "TGeoManager.h"
#include "TGeoVolume.h"
-
-#include "TCallf77.h"
+#include "TObjString.h"
+#include "Fepisor.h"
#ifndef WIN32
# define idnrwr idnrwr_
// TFluka global pointer
TFluka *gFluka = 0;
TFlukaMCGeometry *gMCGeom = 0;
-Int_t kNstep = 0;
+Int_t gNstep = 0;
ClassImp(TFlukaMCGeometry)
fRegionList = 0;
gFluka = (TFluka*)gMC;
gMCGeom = this;
- kNstep = 0;
+ gNstep = 0;
fMatList = new TObjArray(256);
fMatNames = new TObjArray(256);
}
fRegionList = 0;
gFluka = (TFluka*)gMC;
gMCGeom = this;
- kNstep = 0;
+ gNstep = 0;
fMatList = 0;
fMatNames = 0;
}
mix = new TGeoMixture("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);
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);
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);
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);
element = table->GetElement(16);
mix->DefineElement(0, element, 1.);
+ z = element->Z();
break;
}
mix->SetZ(z);
// 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();
+// TList *medlist = gGeoManager->GetListOfMedia();
+// Int_t nmed = medlist->GetSize();
TIter next(matlist);
Int_t nmater = matlist->GetSize();
Int_t nfmater = 0;
matname = element->GetTitle();
ToFlukaString(matname);
rho = 0.999;
+
mat = new TGeoMaterial(matname, element->A(), element->Z(), rho);
mat->SetIndex(nfmater+3);
mat->SetUsed(kTRUE);
fMatNames->Add(objstr);
nfmater++;
}
- Int_t indmat = nfmater;
- TGeoMedium *med;
+
+ fIndmat = nfmater;
+// TGeoMedium *med;
// Adjust material names and add them to FLUKA list
for (i=0; i<nmater; i++) {
mat = (TGeoMaterial*)matlist->At(i);
- if (!mat->IsUsed()) continue;
+ if (!mat->IsUsed()) continue;
z = mat->GetZ();
a = mat->GetA();
rho = mat->GetDensity();
}
matname = mat->GetName();
FlukaMatName(matname);
+/*
// material with one element: create it as mixture since it can be duplicated
if (!mat->IsMixture()) {
// normal materials
}
mat = (TGeoMaterial*)mix;
}
+*/
mat->SetIndex(nfmater+3);
objstr = new TObjString(matname.Data());
fMatList->Add(mat);
mat = (TGeoMaterial*)fMatList->At(i);
// mat->SetUsed(kFALSE);
mix = 0;
-// out << "* " << mat->GetName() << endl;
out << setw(10) << "MATERIAL ";
out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
objstr = (TObjString*)fMatNames->At(i);
out << setw(8) << matname.Data() << endl;
if (!mix) {
// add LOW-MAT card for NEON to associate with ARGON neutron xsec
- if (z==10 && matname.Contains("NEON")) {
+ if (z==10) {
out << setw(10) << "LOW-MAT ";
out.setf(static_cast<std::ios::fmtflags>(0),std::ios::floatfield);
out << setw(10) << setiosflags(ios::fixed) << setprecision(1) << Double_t(mat->GetIndex());
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<std::ios::fmtflags>(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;
PrintHeader(out, "TGEO MATERIAL ASSIGNMENTS");
for (i=1; i<=nvols; i++) {
+
vol = gGeoManager->GetVolume(i);
mat = vol->GetMedium()->GetMaterial();
-// mat->SetUsed(kTRUE);
idmat = mat->GetIndex();
for (Int_t j=0; j<nfmater; j++) {
mat = (TGeoMaterial*)fMatList->At(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<std::ios::fmtflags>(0),std::ios::floatfield);
out << setw(10) << setiosflags(ios::fixed) << Double_t(idmat);
out << setw(10) << setiosflags(ios::fixed) << Double_t(i);
out << setw(10) << "0.0";
out << setw(10) << "0.0";
- out << setw(10) << setiosflags(ios::fixed) << flagfield;
+ out << setw(10) << setiosflags(ios::fixed) << hasfield;
out << setw(10) << "0.0";
out << endl;
}
out.close();
fLastMaterial = nfmater+2;
-/*
- TGeoMaterial *mat1 = 0;
- for (i=1; i<=nvols; i++) {
- vol = gGeoManager->GetVolume(i);
- med = vol->GetMedium();
- mat = med->GetMaterial();
- printf("Region %d: %s\n", i, vol->GetName());
- printf(" medium %d: %s\n", med->GetId(), med->GetName());
- for (j=0; j<nfmater; j++) {
- mat1 = (TGeoMaterial*)fMatList->At(j);
- if (mat1 != mat) continue;
- objstr = (TObjString*)fMatNames->At(j);
- matname = objstr->GetString();
- break;
- }
- if (mat1 != mat) printf(" (*) material not found in Fluka list\n");
- printf(" material %s (at ind=%d): FlukaID=%d FlukaName=%s\n",
- mat->GetName(), j, mat->GetIndex(), matname.Data());
- if (mat->GetCerenkovProperties()) printf(" Cerenkov properties found\n");
- }
-*/
- if (!gFluka->IsGeneratePemf()) return;
- // Write peg files
- char number[20];
- for (i=indmat; i<nfmater; i++) {
- mat = (TGeoMaterial*)fMatList->At(i);
- if (!mat->IsUsed()) continue;
- sname = "mat";
- sprintf(number, "%d", i);
- sname.Append(number);
- WritePegFile(i);
- sname.Prepend("$FLUPRO/pemf/rpemf peg/");
- gSystem->Exec(sname.Data());
+ if (!gFluka->IsGeneratePemf()) {
+ if (gSystem->AccessPathName("FlukaVmc.pemf")) Fatal("CreateFlukaMatFile", "No pemf file in working directory");
+ return;
}
- sname = "cat peg/*.pemf > peg/alice.pemf";
- gSystem->Exec(sname.Data());
- sname = "mv peg/alice.pemf alice.pemf";
- gSystem->Exec(sname.Data());
+}
+
+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) const
+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
+ // Write the .peg file for one material
+
TGeoMaterial *mat = (TGeoMaterial*)fMatList->At(imat);
TString name = ((TObjString*)fMatNames->At(imat))->GetString();
TString line;
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;
- if (dens<0.01) out << " GASP=1." << 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=" <<issb_parm[3] <<", AFACT="<<issb_parm[4] <<", SK="
+ << issb_parm[5] << ", DELTA0=" << issb_parm[6];
+ }
+ else {
+ cout << "WARNING: Strange element, Sternheimer parameters not found" << endl;
+ (*NoStern)++;
+ }
+
+ if (dens<0.01) {
+ (*countGas)++;
+ out << " GASP=1." << endl;
+ }
+
out << " &END" << endl;
out << name.Data() << endl;
out << elem->GetName() << endl;
- } else {
+
+ }
+ 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; j<nel; j++) {
+ zt[j] = (mix->GetZmixt())[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=";
- line = "";
- for (i=0; i<nel; i++) {
- sprintf(number, "%f", mix->GetWmixt()[i]);
- line += number;
- line += ", ";
- if (line.Length() > 30) {
- out << line.Data() << endl;
- line = "";
- }
+ 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<nel; j++) {
+ out << " " << wt[j] << ",";
+ line += Form(" %g,", wt[j] );
+ if( line.Length() > 60 ) { out << endl; line = ""; }
}
- if (line.Length()) out << " " << line.Data() << endl;
- out << " RHO=" << dens;
- if (dens<0.01) out << ", GASP=1." << endl;
- out << " &END" << endl;
+ 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; i<nel; i++) {
elem = mix->GetElement(i);
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=0.56099906, UE=3000000., AP=.03, UP=3000000. $END" << endl;
+ out << " $INP AE=" << ae << ", UE=" << ue <<", AP=" << ap << ", UP=" << up << " $END" << endl;
out << "PWLF" << endl;
out << " $INP NALE=300, NALG=400, NALR=100 $END" << endl;
out << "DECK" << endl;
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
+ };
+
+ 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].nelems; i++) {
+ in >> 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<sternDataArray[is].nelems; j++) {
+ if( sternDataArray[is].Z[j] == Int_t(zelem[j]) &&
+ TMath::Abs( (sternDataArray[is].wt[j] - welem[j])/sternDataArray[is].wt[j] ) < 0.1 )
+ nmatch++;
+ }
+ }
+
+ if( sternDataArray[is].nelems > 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
{
//_____________________________________________________________________________
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;
//_____________________________________________________________________________
void TFlukaMCGeometry::FlukaMatName(TString &str) const
{
+// Convert a name to upper case 8 chars.
ToFlukaString(str);
Int_t ilast;
for (ilast=7; ilast>0; ilast--) if (str(ilast)!=' ') break;
//_____________________________________________________________________________
void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
Double_t *pV, Int_t &oldReg , const Int_t &oldLttc, Double_t &propStep,
- Int_t &nascFlag, Double_t &retStep, Int_t &newReg,
+ Int_t &/*nascFlag*/, Double_t &retStep, Int_t &newReg,
Double_t &saf, Int_t &newLttc, Int_t <tcFlag,
Double_t *sLt, Int_t *jrLt)
{
// Initialize FLUKa point and direction;
- kNstep++;
-/*
- if (kNstep>0) {
- gMCGeom->SetDebugMode(kTRUE);
- gFluka->SetVerbosityLevel(3);
- }
- if (kNstep>6520) {
- gMCGeom->SetDebugMode(kFALSE);
- gFluka->SetVerbosityLevel(0);
- }
- if ((kNstep%10)==0) printf("step %i\n", kNstep);
-*/
+ gNstep++;
+ NORLAT.xn[0] = pSx;
+ NORLAT.xn[1] = pSy;
+ NORLAT.xn[2] = pSz;
+ NORLAT.wn[0] = pV[0];
+ NORLAT.wn[1] = pV[1];
+ NORLAT.wn[2] = pV[2];
if (gMCGeom->IsDebugging()) {
printf("========== Inside G1WR\n");
printf(" point/dir:(%14.9f, %14.9f, %14.9f, %g, %g, %g)\n", pSx,pSy,pSz,pV[0],pV[1],pV[2]);
printf(" oldReg=%i oldLttc=%i pstep=%f\n",oldReg, oldLttc, propStep);
+ }
+ Int_t olttc = oldLttc;
+ if (oldLttc<0) {
+ gGeoManager->FindNode(pSx,pSy,pSz);
+ olttc = gGeoManager->GetCurrentNodeId()+1;
}
- gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
- gGeoManager->SetCurrentDirection(pV);
- gMCGeom->SetCurrentRegion(oldReg, oldLttc);
+ gMCGeom->SetCurrentRegion(oldReg, olttc);
// Initialize default return values
lttcFlag = 0;
- jrLt[lttcFlag] = oldLttc;
+ jrLt[lttcFlag] = olttc;
sLt[lttcFlag] = propStep;
jrLt[lttcFlag+1] = -1;
sLt[lttcFlag+1] = 0.;
newReg = oldReg;
- newLttc = oldLttc;
+ newLttc = olttc;
// check if dummy boundary flag is set
Int_t curLttc, curReg;
if (gFluka->IsDummyBoundary()) {
- // printf("Dummy boundary intercepted. Point is: %f, %f, %f\n", pSx, pSy, pSz);
- Bool_t crossedDummy = (oldLttc == TFlukaMCGeometry::kLttcVirtual)?kTRUE:kFALSE;
+// printf("Dummy boundary intercepted. Point is: %f, %f, %f\n", pSx, pSy, pSz);
+ Bool_t crossedDummy = (olttc == TFlukaMCGeometry::kLttcVirtual)?kTRUE:kFALSE;
if (crossedDummy) {
// FLUKA crossed the dummy boundary - update new region/history
retStep = 0.;
}
// Reset outside flag
- if (gGeoManager->IsOutside()) {
- gGeoManager->SetOutside(kFALSE);
- gGeoManager->CdTop();
- }
+ gGeoManager->SetOutside(kFALSE);
// Reset dummy boundary flag
gFluka->SetDummyBoundary(0);
curLttc = gGeoManager->GetCurrentNodeId()+1;
curReg = gGeoManager->GetCurrentVolume()->GetNumber();
- if (oldLttc != curLttc) {
+ if (olttc != curLttc) {
// FLUKA crossed the boundary : we trust that the given point is really there,
// so we just update TGeo state
- gGeoManager->CdNode(oldLttc-1);
+ gGeoManager->CdNode(olttc-1);
curLttc = gGeoManager->GetCurrentNodeId()+1;
curReg = gGeoManager->GetCurrentVolume()->GetNumber();
if (gMCGeom->IsDebugging()) printf(" re-initialized point: curReg=%i curLttc=%i\n", curReg, curLttc);
}
- // Now the current TGeo state reflects the FLUKA state
+ // Now the current TGeo state reflects the FLUKA state
+ Double_t extra = 1.E-10;
+
+
+// printf("ERROR: (%f, %f, %f)\n",pSx,pSy,pSz);
if (gMCGeom->IsDebugging()) printf(" current path: %s\n", gGeoManager->GetPath());
- Double_t extra = 1E-6;
- Double_t tmpStep = propStep + extra;
- gGeoManager->FindNextBoundary(-tmpStep);
+ gGeoManager->SetCurrentPoint(pSx+extra*pV[0], pSy+extra*pV[1], pSz+extra*pV[2]);
+ gGeoManager->SetCurrentDirection(pV);
+ gGeoManager->FindNextBoundary(-propStep);
Double_t snext = gGeoManager->GetStep();
- // !!!!!
- if (snext<=0) {
- // FLUKA is in the wrong region, notify it
- if (gMCGeom->IsDebugging()) printf("ERROR: snext=%f\n", snext);
-// newReg = -3;
-// return;
- snext = extra;
- }
+
+ if (snext<=0.0) {
+ saf = 0.0;
+ newReg = -3;
+ sLt[lttcFlag] = 0.0;
+ if (gMCGeom->IsDebugging()) printf("BACK SCATTERING\n");
+ return;
+ }
+
+ snext += extra;
saf = gGeoManager->GetSafeDistance();
- Bool_t cross = kFALSE;
- Bool_t onBound = kFALSE;
- if (snext<tmpStep) {
- // We have some boundary in the way
- Double_t dd = snext-propStep;
- if (dd < 0) {
- cross = kTRUE;
- dd = -dd;
- }
- if (dd < 1E-8) onBound = kTRUE;
- }
- snext += 1.E-8;
- if (gMCGeom->IsDebugging()) {
- if (!cross) printf(" physical step approved: %f\n", propStep);
- else printf(" boundary crossing at: %f\n", snext);
- if (onBound) printf(" step on boundary limit ! NASC=%i\n", nascFlag);
- }
- if (!cross) {
+ saf -= extra;
+ if (saf<0) saf=0.0;
+ else saf -= saf*3.0e-09;
+ NORLAT.distn = snext;
+ NORLAT.xn[0] += snext*pV[0];
+ NORLAT.xn[1] += snext*pV[1];
+ NORLAT.xn[2] += snext*pV[2];
+// saf = 0.0; // !!! TEMPORARY FOR TESTING MAGSPHF - TO BE REMOVED
+ if (snext>propStep) {
// Next boundary further than proposed step, which is approved
retStep = propStep;
sLt[lttcFlag] = propStep;
return;
}
// The next boundary is closer. We try to cross it.
+ gGeoManager->SetCurrentPoint(pSx,pSy,pSz);
Double_t *point = gGeoManager->GetCurrentPoint();
Double_t *dir = gGeoManager->GetCurrentDirection();
Double_t pt[3];
memcpy(pt, point, 3*sizeof(Double_t));
Int_t i;
- for (i=0;i<3;i++) point[i] += snext*dir[i];
+ for (i=0;i<3;i++) point[i] += (snext+1E-6)*dir[i];
+ // locate next region
gGeoManager->FindNode();
newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
- if (newLttc == oldLttc) {
- // brute force ...
- // Just try a fast extra small step
- snext += 1E-6;
- for (i=0;i<3;i++) point[i] = pt[i]+snext*dir[i];
- gGeoManager->FindNode();
- newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
- if (newLttc == oldLttc) {
- // check if state changes at the end of the proposed step
- for (i=0;i<3;i++) point[i] = pt[i]+propStep*dir[i];
- gGeoManager->FindNode();
- newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
- if (newLttc==oldLttc) {
- // approve step
- retStep = propStep;
- sLt[lttcFlag] = propStep;
- return;
- }
- // snext is underestimated - we will create a virtual one to overcome the error
-// printf("some boundary in the way...\n");
- }
- }
gGeoManager->SetCurrentPoint(pt);
-// newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
newReg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
if (gMCGeom->IsDebugging()) printf(" newReg=%i newLttc=%i\n", newReg, newLttc);
// We really crossed the boundary, but is it the same region ?
gMCGeom->SetNextRegion(newReg, newLttc);
- if (newReg == oldReg) {
+ if (newReg==oldReg && newLttc!=olttc) {
// Virtual boundary between replicants
if (gMCGeom->IsDebugging()) printf(" DUMMY boundary\n");
newReg = 1; // cheat FLUKA telling it it crossed the TOP region
jrLt[lttcFlag+1] = -1;
sLt[lttcFlag+1] = 0.;
- if (newLttc!=oldLttc) {
- if (gGeoManager->IsOutside()) {
- gGeoManager->SetOutside(kFALSE);
- gGeoManager->CdTop();
- }
- gGeoManager->CdTop();
- if (!gGeoManager->GetCurrentMatrix()->IsIdentity()) printf("ERROR at step %i\n", kNstep);
- gGeoManager->CdNode(oldLttc-1);
+ if (newLttc!=olttc) {
+ if (gGeoManager->IsOutside()) gGeoManager->SetOutside(kFALSE);
+ gGeoManager->CdNode(olttc-1);
}
if (gMCGeom->IsDebugging()) {
printf("=> snext=%g safe=%g\n", snext, saf);
//_____________________________________________________________________________
void lkdbwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
- Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
+ Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
{
if (gMCGeom->IsDebugging()) {
printf("========== Inside LKDBWR (%f, %f, %f)\n",pSx, pSy, pSz);
-// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
+ printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
}
- TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
- if (gGeoManager->IsOutside()) {
- newReg = gMCGeom->NofVolumes()+1;
-// newLttc = gGeoManager->GetCurrentNodeId();
- newLttc = 999999999;
- if (gMCGeom->IsDebugging()) {
- printf("OUTSIDE\n");
- printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
- printf("<= LKMGWR\n");
- }
- flagErr = newReg;
- return;
- }
- newReg = node->GetVolume()->GetNumber();
- newLttc = gGeoManager->GetCurrentNodeId()+1;
- gMCGeom->SetNextRegion(newReg, newLttc);
- flagErr = newReg;
- if (gMCGeom->IsDebugging()) {
- printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
- printf("<= LKDBWR\n");
- }
+ lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,newReg,flagErr,newLttc);
}
//_____________________________________________________________________________
void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
- Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
+ Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
{
if (gMCGeom->IsDebugging()) {
printf("========== Inside LKFXWR (%f, %f, %f)\n",pSx, pSy, pSz);
-// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
+ printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
}
- TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
- if (gGeoManager->IsOutside()) {
- newReg = gMCGeom->NofVolumes()+1;
-// newLttc = gGeoManager->GetCurrentNodeId();
- newLttc = 999999999;
- if (gMCGeom->IsDebugging()) {
- printf("OUTSIDE\n");
- printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
- printf("<= LKMGWR\n");
- }
- flagErr = newReg;
- return;
- }
- newReg = node->GetVolume()->GetNumber();
- newLttc = gGeoManager->GetCurrentNodeId()+1;
- gMCGeom->SetNextRegion(newReg, newLttc);
- flagErr = newReg;
- if (gMCGeom->IsDebugging()) {
- printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
- printf("<= LKFXWR\n");
- }
+ lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,newReg,flagErr,newLttc);
}
//_____________________________________________________________________________
void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
- Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
+ Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
{
if (gMCGeom->IsDebugging()) {
printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz);
-// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
+ printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
}
- TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
- if (gGeoManager->IsOutside()) {
- newReg = gMCGeom->NofVolumes()+1;
-// newLttc = gGeoManager->GetCurrentNodeId();
- newLttc = 999999999;
- if (gMCGeom->IsDebugging()) {
- printf("OUTSIDE\n");
- printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
- printf("<= LKMGWR\n");
- }
- flagErr = newReg;
- return;
- }
- newReg = node->GetVolume()->GetNumber();
- newLttc = gGeoManager->GetCurrentNodeId()+1;
- gMCGeom->SetNextRegion(newReg, newLttc);
- flagErr = newReg;
- if (gMCGeom->IsDebugging()) {
- printf(" out: newReg=%i newLttc=%i\n", newReg, newLttc);
- printf("<= LKMGWR\n");
- }
+ lkwr(pSx,pSy,pSz,pV,oldReg,oldLttc,newReg,flagErr,newLttc);
}
//_____________________________________________________________________________
void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
- Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
+ Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
{
if (gMCGeom->IsDebugging()) {
printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz);
-// printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
+ printf(" in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
printf(" in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
}
+ NORLAT.xn[0] = pSx;
+ NORLAT.xn[1] = pSy;
+ NORLAT.xn[2] = pSz;
+ NORLAT.wn[0] = pV[0];
+ NORLAT.wn[1] = pV[1];
+ NORLAT.wn[2] = pV[2];
TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
if (gGeoManager->IsOutside()) {
newReg = gMCGeom->NofVolumes()+1;
printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz);
printf(" oldReg=%i, newReg=%i\n", oldReg,newReg);
}
-// Int_t curreg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
-// Int_t curLttc = gGeoManager->GetCurrentNodeId()+1;
-// if (gMCGeom->IsDebugging()) printf(" curReg=%i, curLttc=%i in: %s\n", curreg, curLttc, gGeoManager->GetPath());
-// Bool_t regsame = (curreg==oldReg)?kTRUE:kFALSE;
- gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
- gGeoManager->SetCurrentDirection(pVx,pVy,pVz);
-/*
- if (!regsame) {
- if (gMCGeom->IsDebugging()) printf(" REGIONS DOEN NOT MATCH\n");
- gGeoManager->FindNode();
- curreg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
- curLttc = gGeoManager->GetCurrentNodeId()+1;
- if (gMCGeom->IsDebugging()) printf(" re-initialized point: curReg=%i curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath());
- }
-*/
+ gGeoManager->SetCurrentPoint(NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]);
+ gGeoManager->SetCurrentDirection(NORLAT.wn[0], NORLAT.wn[1], NORLAT.wn[2]);
Double_t *dnorm = gGeoManager->FindNormalFast();
flagErr = 0;
if (!dnorm) {
norml[0] = -dnorm[0];
norml[1] = -dnorm[1];
norml[2] = -dnorm[2];
- if (gMCGeom->IsDebugging()) printf(" normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]);
-// curreg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
-// curLttc = gGeoManager->GetCurrentNodeId()+1;
if (gMCGeom->IsDebugging()) {
-// printf(" final location: curReg=%i, curLttc=%i in %s\n", curreg,curLttc,gGeoManager->GetPath());
+ printf(" normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]);
printf("<= NRMLWR\n");
}
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff