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=" <<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 {
-
- 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=" << 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 = ""; }
- }
- 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);
- 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].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
{
git://git.uio.no / u/mrichter/AliRoot.git / commitdiff