composite_calorimeter/src/CCalMaterialFactory.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 ///////////////////////////////////////////////////////////////////////////////
0027 // File: CCalMaterialFactory.cc
0028 // Description: CCalMaterialFactory is a factory class to vuild G4Material
0029 //              from CCalMaterial and CCalAmaterial
0030 ///////////////////////////////////////////////////////////////////////////////
0031 #include <fstream>
0032 #include <stdlib.h>
0033
0034 #include "CCalMaterialFactory.hh"
0035 #include "CCalutils.hh"
0036
0037 #include "G4PhysicalConstants.hh"
0038 #include "G4SystemOfUnits.hh"
0039 #include "G4Material.hh"
0040
0041 //#define ddebug
0042 //#define debug
0043
0044 typedef CCalMaterial*  ptrCCalMaterial;
0045 typedef CCalAMaterial* ptrCCalAMaterial;
0046
0047
0048 CCalMaterialFactory * CCalMaterialFactory::instance = 0;
0049 G4String CCalMaterialFactory::elementfile = "";
0050 G4String CCalMaterialFactory::mixturefile = "";
0051
0052
0053 CCalMaterialFactory* CCalMaterialFactory::getInstance(const G4String& matfile,
0054                                                       const G4String& mixfile){
0055   if ((matfile=="" || matfile==elementfile) &&
0056       (mixfile=="" || mixfile==mixturefile))
0057     return getInstance();
0058   else if ((matfile != "" && elementfile != "" && matfile != elementfile) ||
0059            (mixfile != "" && mixturefile != "" && mixfile != mixturefile)) {
0060     G4cerr << "ERROR: Trying to get materials from " << matfile << " and "
0061          << mixfile << " while previously were retrieved from "
0062          << elementfile << " and " << mixturefile << "." << G4endl;
0063     return 0;
0064   } else {
0065     if (elementfile == "")
0066       elementfile=matfile;
0067     if (mixturefile == "")
0068       mixturefile=mixfile;
0069     return getInstance();
0070   }
0071 }
0072
0073
0074 CCalMaterialFactory* CCalMaterialFactory::getInstance(const G4String& matfile){
0075   return getInstance(matfile,matfile);
0076 }
0077
0078
0079 CCalMaterialFactory* CCalMaterialFactory::getInstance(){
0080   if (elementfile=="" || mixturefile=="") {
0081     G4cerr << "ERROR: You haven't defined files to be used for materials in "
0082          << "CCalMaterialFactory::getInstance(const G4String&,const G4String&)"
0083          << G4endl;
0084     return 0;
0085   }
0086
0087   if (instance==0) {
0088     instance = new CCalMaterialFactory;
0089     return instance;
0090   }
0091   else
0092     return instance;
0093 }
0094
0095
0096 CCalMaterialFactory::~CCalMaterialFactory(){
0097   CCalMaterialTable::iterator ite;
0098   for(ite = theCCalMaterials.begin(); ite != theCCalMaterials.end(); ite++ ){
0099     delete *ite;
0100   }
0101   theCCalMaterials.clear();
0102   CCalAMaterialTable::iterator itea;
0103   for(itea = theCCalAMaterials.begin(); itea != theCCalAMaterials.end();
0104       itea++ ){
0105     delete *itea;
0106   }
0107   theCCalAMaterials.clear();
0108 }
0109
0110
0111 G4Material* CCalMaterialFactory::findMaterial(const G4String & mat) const {
0112   G4Material* theMat=findG4Material(mat);
0113
0114   if (theMat) {
0115 #ifdef ddebug
0116     G4cout << "Material " << mat << " already defined. Returning previous "
0117          << "instance." << G4endl;
0118 #endif
0119     return theMat;
0120   } else {
0121     CCalMaterial* CCalmat=findCCalMaterial(mat);
0122     if (CCalmat){
0123       G4Material* G4Mat = new G4Material(CCalmat->Name(),
0124                                          CCalmat->Density()*g/cm3,
0125                                          CCalmat->NElements());
0126       for(G4int i=0; i<CCalmat->NElements(); i++) {
0127         G4Element* elem = findElement(CCalmat->Element(i));
0128         if (!elem) {
0129           G4ExceptionDescription ed;
0130           ed << "   Could not build material " << mat << "." << G4endl;
0131           G4Exception("CCalMaterialFactory::findMaterial()","ccal001",
0132                       FatalException,ed);
0133         }
0134         G4Mat->AddElement(elem, CCalmat->Weight(i));
0135       }
0136 #ifdef ddebug
0137     G4cout << "Material " << mat << " has been built successfully." << G4endl;
0138 #endif
0139       return G4Mat;
0140     } else {
0141       G4cerr << "ERROR: Material " << mat << " not found in CCal database!!!"
0142            << G4endl;
0143       return 0;
0144     }
0145   }
0146 }
0147
0148
0149 G4Element* CCalMaterialFactory::findElement(const G4String & mat) const {
0150   const G4ElementTable  theElements = *(G4Element::GetElementTable());
0151   for (unsigned int i=0; i<theElements.size(); i++)
0152     if (theElements[i]->GetName()==mat){
0153 #ifdef ddebug
0154       G4cout << "Element " << mat << " found!" << G4endl;
0155 #endif
0156       return theElements[i];
0157     }
0158   return 0;
0159 }
0160
0161
0162 G4Element* CCalMaterialFactory::addElement(const G4String & name,
0163                                            const G4String & symbol,
0164                                            G4double Z, G4double A,
0165                                            G4double density) {
0166
0167   G4Element* theEl = new G4Element(name, symbol, Z, A*g/mole);
0168   //Make it also as a material.
0169   CCalAMaterial* theMat = new CCalAMaterial(name,A,density);
0170   theCCalAMaterials.push_back(theMat);
0171
0172 #ifdef ddebug
0173   G4cout << "Element " << name << " created!" << G4endl;
0174 #endif
0175   return theEl;
0176 }
0177
0178
0179 G4Material* CCalMaterialFactory::addMaterial(const G4String& name,
0180                                              G4double density,
0181                                              G4int nconst,
0182                                              G4String mats[],
0183                                              G4double prop[],
0184                                              MatDescription md){
0185   addCCalMaterial(name, density, nconst, mats, prop, md);
0186   return findMaterial(name);
0187 }
0188
0189
0190 void CCalMaterialFactory::readElements(const G4String& matfile) {
0191
0192   G4String path = "NULL";
0193   if (std::getenv("CCAL_GLOBALPATH"))
0194     path = std::getenv("CCAL_GLOBALPATH");
0195
0196   G4cout << " ==> Opening file " << matfile << " to read elements..." << G4endl;
0197   std::ifstream is;
0198   G4bool ok = openGeomFile(is, path, matfile);
0199   if (!ok) {
0200     G4cerr << "ERROR: Could not open file " << matfile << G4endl;
0201     return;
0202   }
0203
0204   // Find *DO GMAT
0205   findDO(is, G4String("GMAT"));
0206
0207   readElements(is);
0208
0209   is.close();
0210 }
0211
0212
0213 void CCalMaterialFactory::readMaterials(const G4String& matfile) {
0214
0215   G4String path = "NULL";
0216   if (std::getenv("CCAL_GLOBALPATH"))
0217     path = std::getenv("CCAL_GLOBALPATH");
0218
0219   G4cout << " ==> Opening file " << matfile << " to read materials..." << G4endl;
0220   std::ifstream is;
0221   bool ok = openGeomFile(is, path, matfile);
0222   if (!ok) {
0223     G4cerr << "ERROR: Could not open file " << matfile << G4endl;
0224     return;
0225   }
0226
0227   // Find *DO GMIX
0228   findDO(is, G4String("GMIX"));
0229
0230   readMaterials(is);
0231
0232   is.close();
0233 }
0234
0235
0236 //===========================================================================
0237 // Protected & private methods ==============================================
0238
0239
0240 G4Material* CCalMaterialFactory::findG4Material(const G4String & mat) const {
0241   const G4MaterialTable theG4Materials = *(G4Material::GetMaterialTable());
0242   for (unsigned int i=0; i<theG4Materials.size(); i++) {
0243     if (theG4Materials[i]->GetName()==mat){
0244       return theG4Materials[i];
0245     }
0246   }
0247   return 0;
0248 }
0249
0250
0251 CCalMaterial* CCalMaterialFactory::findCCalMaterial(const G4String & mat)
0252   const {
0253   for (unsigned int i=0; i<theCCalMaterials.size(); i++)
0254     if (theCCalMaterials[i]->Name()==mat){
0255 #ifdef ddebug
0256       G4cout << "CCalMaterial " << mat << " found!" << G4endl;
0257 #endif
0258       return theCCalMaterials[i];
0259     }
0260   return (CCalMaterial*) findCCalAMaterial(mat);
0261 }
0262
0263
0264 CCalAMaterial* CCalMaterialFactory::findCCalAMaterial(const G4String & mat)
0265   const {
0266   for (unsigned int i=0; i<theCCalAMaterials.size(); i++)
0267     if (theCCalAMaterials[i]->Name()==mat){
0268 #ifdef ddebug
0269       G4cout << "CCalMaterial " << mat << " found!" << G4endl;
0270 #endif
0271       return theCCalAMaterials[i];
0272     }
0273   return 0;
0274 }
0275
0276
0277 CCalMaterial* CCalMaterialFactory::addCCalMaterial(const G4String& name,
0278                                                    G4double density,
0279                                                    G4int nconst,
0280                                                    G4String mats[],
0281                                                    G4double prop[],
0282                                                    MatDescription md){
0283   ptrCCalMaterial* matcol=0;
0284   ptrCCalAMaterial* amatcol=0;
0285
0286   if (md==byAtomic)
0287     amatcol = new ptrCCalAMaterial[nconst];
0288   else
0289     matcol = new ptrCCalMaterial[nconst];
0290
0291   for (G4int i=0; i<nconst; i++){
0292     if (md==byAtomic) {
0293       CCalAMaterial* amat = findCCalAMaterial(mats[i]);
0294       if (amat)
0295         amatcol[i]=amat;
0296       else {
0297         G4cerr << "ERROR: Trying to build" << name << " out of unknown "
0298              << mats[i] << "." << G4endl
0299              << "Skiping this material!" << G4endl;
0300         delete[] amatcol;
0301         return 0;
0302       }
0303     } //by Atomic fractions
0304     else {
0305       CCalMaterial* mat = findCCalMaterial(mats[i]);
0306       if (mat)
0307         matcol[i]=mat;
0308       else {
0309         G4cerr << "ERROR: Trying to build" <<name << " out of unknown "
0310              << mats[i] << "." << G4endl
0311              << "Skiping this material!" << G4endl;
0312         delete[] matcol;
0313         return 0;
0314       }
0315     }
0316   } //for
0317
0318   //Let's do the CCalMaterial!
0319   if (md==byAtomic) {
0320     CCalAMaterial* amaterial = new CCalAMaterial(name, density, nconst,
0321                                                  amatcol, prop);
0322     delete[] amatcol;
0323     theCCalAMaterials.push_back(amaterial);
0324 #ifdef ddebug
0325     G4cout << *amaterial << G4endl;
0326 #endif
0327     return amaterial;
0328   } else {
0329     CCalMaterial::FractionType ft;
0330     if (md == byWeight)
0331       ft=CCalMaterial::FTWeight;
0332     else
0333       ft=CCalMaterial::FTVolume;
0334     CCalMaterial* material = new CCalMaterial(name, density, nconst,
0335                                             matcol, prop, ft);
0336     delete[] matcol;
0337     theCCalMaterials.push_back(material);
0338 #ifdef ddebug
0339     G4cout << *material << G4endl;
0340 #endif
0341     return material;
0342   }
0343 }
0344
0345
0346 void CCalMaterialFactory::readElements(std::ifstream& is){
0347   G4String name, symbol;
0348
0349   G4cout << "     ==> Reading elements... " << G4endl;
0350 #ifdef debug
0351   G4cout << "       Element    \tsymbol\tA\tZ\tdensity\tX_0          abs_l"<< G4endl;
0352 #endif
0353   //There should come the list of materials. #. Defines a comment
0354   //*DO defines the beguining of the Mixes block.
0355
0356   readName(is,name);
0357   while (name != "*ENDDO") {
0358     //It should be an element definition
0359     G4double A, Z, density;
0360     is >> symbol >> A >> Z >> density >> jump;
0361 #ifdef debug
0362     G4cout << "       " << name << "    \t" << symbol << "\t"
0363          << A << "\t" << Z << "\t" << density << G4endl;
0364 #endif
0365     addElement(name, symbol, Z, A, density);
0366     readName(is,name);
0367   };
0368   G4cout << "     " << G4Element::GetElementTable()->size()
0369        << " elements read from file" << G4endl << G4endl;
0370 }
0371
0372
0373 void CCalMaterialFactory::readMaterials(std::ifstream& is){
0374   G4String name, matname;
0375
0376   G4cout << "     ==> Reading materials... " << G4endl;
0377
0378   //Take into account the special case of vacuum...
0379 #ifdef debug
0380   G4cout <<"       \"Vacuum\"" << G4endl;
0381 #endif
0382   G4double density  = universe_mean_density;   //from PhysicalConstants.h
0383   G4double pressure = 1.E-19*pascal;
0384   G4double temperature = 0.1*kelvin;
0385   new G4Material("Vacuum", /*Z=*/ 1., /*A=*/ 1.01*g/mole,
0386                  density, kStateGas, temperature, pressure);
0387
0388   //There should come the list of materials. #. Defines a comment
0389   //*ENDDO defines the block.
0390   readName(is,name);
0391   while (name != "*ENDDO") {
0392     //It should be a material definition
0393     matname=name;
0394     G4int nElem;
0395     G4double dens;
0396     is >> nElem >> dens >> jump;
0397
0398 #ifdef debug
0399     G4cout <<"       " << matname
0400            << " made of " << nElem
0401            << " elements. Density=" << dens
0402            << G4endl;
0403 #endif
0404
0405     G4int absnelem = std::abs(nElem);
0406
0407     G4String* mats    = new G4String[absnelem];
0408     G4double* weights = new G4double[absnelem];
0409
0410     G4double prop;
0411     for(int i=0; i<absnelem; i++) {
0412       readName(is, name);
0413       is >> prop >> jump;
0414       mats[i]=name;
0415       weights[i]=std::abs(prop);
0416     } //for...
0417     MatDescription md;
0418     if (nElem>0 && prop<0)
0419       md = byAtomic;
0420     else if (nElem>0)
0421       md = byWeight;
0422     else
0423       md = byVolume;
0424
0425     addCCalMaterial(matname, dens, absnelem, mats, weights, md);
0426     delete[] mats;
0427     delete[] weights;
0428
0429     readName(is,name);
0430   };  //while
0431
0432   G4cout << "     " << theCCalMaterials.size() << " materials read from "
0433        << mixturefile << G4endl << G4endl;
0434 }
0435
0436
0437 CCalMaterialFactory::CCalMaterialFactory() {
0438   readElements (elementfile);
0439   readMaterials(mixturefile);
0440 }