g4sim: src/MaterialSvc.cc Source File

00001 //---------------------------------------------------------------------------//
00002 //Description: Setup Material
00003 //Author: Wu Chen(wuchen@mail.ihep.ac.cn)
00004 //Created: 17 Oct, 2012
00005 //Comment:
00006 //---------------------------------------------------------------------------//
00007 
00008 #include "MaterialSvc.hh"
00009 
00010 #include "G4Isotope.hh"
00011 #include "G4Material.hh"
00012 #include "G4UnitsTable.hh"
00013 #include "myglobals.hh"
00014 
00015 #include <sstream>
00016 #include <fstream>
00017 #include <string.h>
00018 #include <iostream>
00019 #include <iomanip>
00020 
00021 #include "DEBUG.hh"
00022 
00023 MaterialSvc* MaterialSvc::fMaterialSvc = 0;
00024 
00025 MaterialSvc::MaterialSvc()
00026 {
00027         if (fMaterialSvc){
00028                 G4Exception("MaterialSvc::MaterialSvc()","Run0031",
00029                                 FatalException, "MaterialSvc constructed twice.");
00030         }
00031         fMaterialSvc = this;
00032 }
00033 
00034 MaterialSvc::~MaterialSvc()
00035 {
00036         printf("~MaterialSvc\n");
00037 }
00038 
00039 MaterialSvc* MaterialSvc::GetMaterialSvc(){
00040         if ( !fMaterialSvc ){
00041                 fMaterialSvc = new MaterialSvc;
00042         }
00043         return fMaterialSvc;
00044 }
00045 
00046 void MaterialSvc::SetMaterial( G4String file_name ){
00047         if(file_name[0] != '/'){ // Relative Dir
00048                 G4String dir_name = getenv("CONFIGUREROOT");
00049                 if (dir_name[dir_name.size()-1] != '/') dir_name.append("/");
00050                 file_name = dir_name + file_name;
00051         }
00052         MAT_LINEVAR(file_name) // DEBUG.hh
00053         std::ifstream fin_card(file_name);
00054         if ( !fin_card ){
00055                 std::cout<<"material list "<<file_name<<" is not available!!!"<<std::endl;
00056                 G4Exception("MaterialSvc::SetMaterial()","Run0031",
00057                                 FatalException, "material list is not available.");
00058         }
00059         std::stringstream buf_card;
00060         std::string s_card;
00061         std::cout<<"########IN MaterialSvc::SetMaterial##########"<<std::endl;
00062         std::cout<<std::setiosflags(std::ios::left)<<std::setw(20) << "Material Name"
00063                                                  <<std::setiosflags(std::ios::left)<<std::setw(15)<< "Z"
00064                                                  <<std::setiosflags(std::ios::left)<<std::setw(15)<< "A"
00065                                                  <<std::setiosflags(std::ios::left)<<std::setw(15)<< "Ionization"
00066                                                  <<std::setiosflags(std::ios::left)<<std::setw(15)<< "Density"
00067                                                  <<std::setiosflags(std::ios::left)<<std::setw(15)<< "Radlen"
00068                                                  <<std::endl;
00069         std::cout<<std::setiosflags(std::ios::left)<<std::setw(20) << ""
00070                                                  <<std::setiosflags(std::ios::left)<<std::setw(15)<< ""
00071                                                  <<std::setiosflags(std::ios::left)<<std::setw(15)<< "g/mole"
00072                                                  <<std::setiosflags(std::ios::left)<<std::setw(15)<< "eV"
00073                                                  <<std::setiosflags(std::ios::left)<<std::setw(15)<< "g/cm3"
00074                                                  <<std::setiosflags(std::ios::left)<<std::setw(15)<< "mm"
00075                                                  <<std::endl;
00076         while(getline(fin_card,s_card)){
00077                 buf_card.str("");
00078                 buf_card.clear();
00079                 buf_card<<s_card;
00080                 const char* c = s_card.c_str();
00081                 int length = strlen(c);
00082                 int offset = 0;
00083                 for ( ; offset < length; offset++ ){
00084                         if ( c[offset] != ' ' ) break;
00085                 }
00086                 if ( c[offset] == '#' ) continue;
00087                 else if ( c[offset] == '/' && c[offset+1] == '/' ) continue;
00088                 else if ( length - offset == 0 ) continue;
00089                 std::string keyword;
00090                 buf_card>>keyword;
00091                 if ( keyword == "TYPE" ){
00092                         buf_card>>fMode;
00093                         continue;
00094                 }
00095                 AddMaterial(s_card);
00096         }
00097         MAT_LINECONT("Finished adding") // DEBUG.hh
00098         fin_card.close();
00099         buf_card.str("");
00100         buf_card.clear();
00101         // print table
00102         //
00103         G4cout << *(G4Element::GetElementTable()) << G4endl;
00104         G4cout << *(G4Material::GetMaterialTable()) << G4endl;
00105 }
00106 
00107 void MaterialSvc::AddMaterial( G4String content ){
00108         G4Material *aMaterial = 0;
00109         G4Element *aElement = 0;
00110         G4String symbol;
00111         G4String name;
00112         G4String s_temp;
00113         double z = 0;
00114         double n = 0;
00115         double a = 0;
00116         double density = 0;
00117         int ncomponents = 0;
00118         int natoms[50];
00119         double comFrac[50];
00120         double rel_dens = 0;
00121         G4String material[50];
00122         G4String element[50];
00123         std::stringstream buf_card;
00124         buf_card.str("");
00125         buf_card.clear();
00126         buf_card<<content;
00127         MAT_LINEVAR(fMode) // DEBUG.hh
00128         if ( fMode == "Isotopes"){
00129                 buf_card>>symbol;
00130                 buf_card>>z;
00131                 buf_card>>n;
00132                 buf_card>>a;
00133                 //    std::cout<<"symbol = "<<symbol<<", a = "<<a<<", z = "<<z<<", n = "<<n<<std::endl;
00134                 new G4Isotope(symbol.c_str(), z, n, a*g/mole );
00135         }
00136         else if ( fMode == "Elements"){
00137                 buf_card>>symbol;
00138                 buf_card>>name;
00139                 buf_card>>z;
00140                 buf_card>>a;
00141                 //    std::cout<<"symbol = "<<symbol<<", name = "<<name<<", a = "<<a<<", z = "<<z<<std::endl;
00142                 new G4Element(name.c_str(), symbol.c_str(), z, a*g/mole );
00143         }
00144         else if ( fMode == "Comp_Elements"){
00145                 buf_card>>symbol;
00146                 buf_card>>name;
00147                 buf_card>>ncomponents;
00148                 double sum_frac = 0;
00149                 for ( int i = 0; i < ncomponents; i++ ){
00150                         buf_card>>element[i];
00151                         buf_card>>comFrac[i];
00152                         sum_frac = sum_frac + comFrac[i];
00153                 }
00154                 //    std::cout<<"symbol = "<<symbol<<", name= "<<name<<", ncomponents = "<<ncomponents<<std::endl;
00155                 //    for ( int i = 0; i < ncomponents; i++ ){
00156                 //      std::cout<<"  "<<i<<": "<<"isotope = "<<element[i]<<", comFrac = "<<comFrac[i]<<std::endl;
00157                 //    }
00158                 if ( sum_frac == 0 ){
00159                         std::cout<<"Please check mass fractions for "<<name<<std::endl;
00160                         G4Exception("MaterialSvc::AddMaterial()","Run0031",
00161                                         FatalException, "total fraction is zero.");
00162                 }
00163                 if ( sum_frac != 1 ){
00164                         std::cout<<"the total fraction for "<<name<<" is not 1!!!"<<std::endl;
00165                         std::cout<<"MaterialSvc will normalize it to 1!"<<std::endl;
00166                         for ( int i = 0; i < ncomponents; i++ ){
00167                                 comFrac[i] = comFrac[i]/sum_frac;
00168                         }
00169                 }
00170                 aElement = new G4Element(name.c_str(), symbol.c_str(), ncomponents);
00171                 for ( int i = 0; i < ncomponents; i++ ){
00172                         G4Isotope* new_iso = G4Isotope::GetIsotope(element[i]);
00173                         aElement->AddIsotope(new_iso, comFrac[i]);
00174                 }
00175         }
00176         else if ( fMode == "Simple_Materials" ){
00177                 buf_card>>name;
00178                 buf_card>>z;
00179                 buf_card>>a;
00180                 buf_card>>density;
00181                 buf_card>>rel_dens;
00182                 density = rel_dens*density*g/cm3;
00183                 //    std::cout<<"name = "<<name<<", a = "<<a<<", z = "<<z<<", density = "<<density<<", rel_dens = "<<rel_dens<<std::endl;
00184                 aMaterial = new G4Material(name.c_str(), z, a*g/mole, density);
00185         }
00186         else if ( fMode == "Molecule_Materials" ){
00187                 buf_card>>name;
00188                 buf_card>>density;
00189                 buf_card>>rel_dens;
00190                 density = rel_dens*density*g/cm3;
00191                 buf_card>>ncomponents;
00192                 for ( int i = 0; i < ncomponents; i++ ){
00193                         buf_card>>element[i];
00194                         buf_card>>natoms[i];
00195                 }
00196                 //    std::cout<<"name = "<<name<<", density = "<<density<<", rel_dens = "<<rel_dens<<", ncomponents = "<<ncomponents<<std::endl;
00197                 //    for ( int i = 0; i < ncomponents; i++ ){
00198                 //      std::cout<<"  "<<i<<": "<<"element = "<<element[i]<<", natoms = "<<natoms[i]<<std::endl;
00199                 //    }
00200                 aMaterial = new G4Material(name.c_str(), density, ncomponents);
00201                 for ( int i = 0; i < ncomponents; i++ ){
00202                         G4Element* new_ele = G4Element::GetElement(element[i]);
00203                         aMaterial->AddElement(new_ele, natoms[i]);
00204                 }
00205         }
00206         else if ( fMode == "MixEle_Materials" ){
00207                 buf_card>>name;
00208                 buf_card>>density;
00209                 buf_card>>rel_dens;
00210                 density = rel_dens*density*g/cm3;
00211                 buf_card>>ncomponents;
00212                 double sum_frac = 0;
00213                 for ( int i = 0; i < ncomponents; i++ ){
00214                         buf_card>>element[i];
00215                         buf_card>>comFrac[i];
00216                         sum_frac = sum_frac + comFrac[i];
00217                 }
00218                 //    std::cout<<"name = "<<name<<", density = "<<density<<", rel_dens = "<<rel_dens<<", ncomponents = "<<ncomponents<<std::endl;
00219                 //    for ( int i = 0; i < ncomponents; i++ ){
00220                 //      std::cout<<"  "<<i<<": "<<"element = "<<element[i]<<", comFrac = "<<comFrac[i]<<std::endl;
00221                 //    }
00222                 if ( sum_frac == 0 ){
00223                         std::cout<<"Please check mass fractions for "<<name<<std::endl;
00224                         G4Exception("MaterialSvc::AddMaterial()","Run0031",
00225                                         FatalException, "total fraction is zero.");
00226                 }
00227                 if ( sum_frac != 1 ){
00228                         std::cout<<"the total fraction for "<<name<<" is not 1!!!"<<std::endl;
00229                         std::cout<<"MaterialSvc will normalize it to 1!"<<std::endl;
00230                         for ( int i = 0; i < ncomponents; i++ ){
00231                                 comFrac[i] = comFrac[i]/sum_frac;
00232                         }
00233                 }
00234                 aMaterial = new G4Material(name.c_str(), density, ncomponents);
00235                 for ( int i = 0; i < ncomponents; i++ ){
00236                         G4Element* new_ele = G4Element::GetElement(element[i]);
00237                         aMaterial->AddElement(new_ele, comFrac[i]);
00238                 }
00239         }
00240         else if ( fMode == "Mixture_Materials" ){
00241                 buf_card>>name;
00242                 buf_card>>rel_dens;
00243                 buf_card>>ncomponents;
00244                 double sum_frac = 0;
00245                 for ( int i = 0; i < ncomponents; i++ ){
00246                         buf_card>>material[i];
00247                         buf_card>>comFrac[i];
00248                         sum_frac += comFrac[i];
00249                 }
00250                 //    std::cout<<"name = "<<name<<", rel_dens = "<<rel_dens<<", ncomponents = "<<ncomponents<<std::endl;
00251                 if ( sum_frac == 0 ){
00252                         std::cout<<"Please check mass fractions for "<<name<<std::endl;
00253                         G4Exception("MaterialSvc::AddMaterial()","Run0031",
00254                                         FatalException, "total fraction is zero.");
00255                 }
00256                 if ( sum_frac != 1 ){
00257                         std::cout<<"the total fraction for "<<name<<" is not 1!!!"<<std::endl;
00258                         std::cout<<"MaterialSvc will normalize it to 1!"<<std::endl;
00259                         for ( int i = 0; i < ncomponents; i++ ){
00260                                 comFrac[i] = comFrac[i]/sum_frac;
00261                         }
00262                 }
00263                 for ( int i = 0; i < ncomponents; i++ ){
00264                         G4Material* new_mat_com = G4Material::GetMaterial(material[i]);
00265                         G4double i_density = new_mat_com->GetDensity();
00266                         //      std::cout<<"  "<<i<<": "<<"material = "<<material[i]<<", comFrac = "<<comFrac[i]<<", density = "<<(i_density/(g/cm3))<<std::endl;
00267                         density += comFrac[i]/(i_density/(g/cm3));
00268                 }
00269                 if ( density!=0 ){
00270                         density = rel_dens/density*g/cm3;
00271                 }
00272                 aMaterial = new G4Material(name.c_str(), density, ncomponents);
00273                 for ( int i = 0; i < ncomponents; i++ ){
00274                         G4Material* new_mat_com = G4Material::GetMaterial(material[i]);
00275                         aMaterial->AddMaterial(new_mat_com, comFrac[i]);
00276                 }
00277         }
00278         else if ( fMode == "VolMix_Materials" ){
00279                 buf_card>>name;
00280                 buf_card>>rel_dens;
00281                 buf_card>>ncomponents;
00282                 double sum_frac = 0;
00283                 for ( int i = 0; i < ncomponents; i++ ){
00284                         buf_card>>material[i];
00285                         buf_card>>comFrac[i];
00286                         sum_frac += comFrac[i];
00287                 }
00288                 //    std::cout<<"name = "<<name<<", rel_dens = "<<rel_dens<<", ncomponents = "<<ncomponents<<std::endl;
00289                 if ( sum_frac == 0 ){
00290                         std::cout<<"Please check mass fractions for "<<name<<std::endl;
00291                         G4Exception("MaterialSvc::AddMaterial()","Run0031",
00292                                         FatalException, "total fraction is zero.");
00293                 }
00294                 if ( sum_frac != 1 ){
00295                         std::cout<<"the total fraction for "<<name<<" is not 1!!!"<<std::endl;
00296                         std::cout<<"MaterialSvc will normalize it to 1!"<<std::endl;
00297                         for ( int i = 0; i < ncomponents; i++ ){
00298                                 comFrac[i] = comFrac[i]/sum_frac;
00299                         }
00300                 }
00301                 for ( int i = 0; i < ncomponents; i++ ){
00302                         G4Material* new_mat_com = G4Material::GetMaterial(material[i]);
00303                         G4double i_density = new_mat_com->GetDensity();
00304                         //      std::cout<<"  "<<i<<": "<<"material = "<<material[i]<<", comFrac = "<<comFrac[i]<<", density = "<<(i_density/(g/cm3))<<std::endl;
00305                         density += comFrac[i]*(i_density);
00306                 }
00307                 density = rel_dens*density;
00308                 aMaterial = new G4Material(name.c_str(), density, ncomponents);
00309                 for ( int i = 0; i < ncomponents; i++ ){
00310                         G4Material* new_mat_com = G4Material::GetMaterial(material[i]);
00311                         G4double i_density = new_mat_com->GetDensity();
00312                         aMaterial->AddMaterial(new_mat_com, comFrac[i]*i_density/density);
00313                 }
00314         }
00315         else{
00316                 std::cout<<"mode "<<fMode<<" is not defined yet!!!"<<std::endl;
00317                 G4Exception("MaterialSvc::AddMaterial()","Run0031",
00318                                 FatalException, "unknown mode type.");
00319         }
00320         // Display useful information
00321         if (aMaterial){
00322                 double Z(0.),A(0.),Ionization(0.),Density(0.),Radlen(0.);
00323                 for(int i=0; i<aMaterial->GetElementVector()->size(); i++){
00324                         Z += (aMaterial->GetElement(i)->GetZ())*
00325                                 (aMaterial->GetFractionVector()[i]);
00326                         A += (aMaterial->GetElement(i)->GetA())*
00327                                 (aMaterial->GetFractionVector()[i]);
00328                 }
00329                 Ionization = aMaterial->GetIonisation()->GetMeanExcitationEnergy();
00330                 Density = aMaterial->GetDensity();
00331                 Radlen = aMaterial->GetRadlen();
00332                 std::cout<<std::setiosflags(std::ios::left)<<std::setw(20) << aMaterial->GetName()
00333                                                  <<std::setiosflags(std::ios::left)<<std::setw(15)<< Z
00334                                                  <<std::setiosflags(std::ios::left)<<std::setw(15)<< A/(g/mole)
00335                                                  <<std::setiosflags(std::ios::left)<<std::setw(15)<< Ionization/eV
00336                                                  <<std::setiosflags(std::ios::left)<<std::setw(15)<< Density/(g/cm3)
00337                                                  <<std::setiosflags(std::ios::left)<<std::setw(15)<< Radlen/mm
00338                                                  <<std::endl;
00339         }
00340 }