g4simulation/g4gdml/PHG4GDMLWriteMaterials.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 // $Id: PHG4GDMLWriteMaterials.cc 70764 2013-06-05 12:54:37Z gcosmo $
0028 //
0029 // class PHG4GDMLWriteMaterials Implementation
0030 //
0031 // Original author: Zoltan Torzsok, November 2007
0032 //
0033 // --------------------------------------------------------------------
0034
0035 #include "PHG4GDMLWriteMaterials.hh"
0036
0037 #include <Geant4/G4Element.hh>
0038 #include <Geant4/G4Isotope.hh>
0039 #include <Geant4/G4Material.hh>
0040 #include <Geant4/G4MaterialPropertiesTable.hh>
0041 #include <Geant4/G4PhysicalConstants.hh>
0042 #include <Geant4/G4PhysicsFreeVector.hh>
0043 #include <Geant4/G4SystemOfUnits.hh>
0044
0045 #include <sstream>
0046
0047 PHG4GDMLWriteMaterials::PHG4GDMLWriteMaterials()
0048   :
0049    materialsElement(nullptr)
0050 {
0051 }
0052
0053 void PHG4GDMLWriteMaterials::
0054     AtomWrite(xercesc::DOMElement* element, const G4double& a)
0055 {
0056   xercesc::DOMElement* atomElement = NewElement("atom");
0057   atomElement->setAttributeNode(NewAttribute("unit", "g/mole"));
0058   atomElement->setAttributeNode(NewAttribute("value", a * mole / g));
0059   element->appendChild(atomElement);
0060 }
0061
0062 void PHG4GDMLWriteMaterials::
0063     DWrite(xercesc::DOMElement* element, const G4double& d)
0064 {
0065   xercesc::DOMElement* DElement = NewElement("D");
0066   DElement->setAttributeNode(NewAttribute("unit", "g/cm3"));
0067   DElement->setAttributeNode(NewAttribute("value", d * cm3 / g));
0068   element->appendChild(DElement);
0069 }
0070
0071 void PHG4GDMLWriteMaterials::
0072     PWrite(xercesc::DOMElement* element, const G4double& P)
0073 {
0074   xercesc::DOMElement* PElement = NewElement("P");
0075   PElement->setAttributeNode(NewAttribute("unit", "pascal"));
0076   PElement->setAttributeNode(NewAttribute("value", P / hep_pascal));
0077   element->appendChild(PElement);
0078 }
0079
0080 void PHG4GDMLWriteMaterials::
0081     TWrite(xercesc::DOMElement* element, const G4double& T)
0082 {
0083   xercesc::DOMElement* TElement = NewElement("T");
0084   TElement->setAttributeNode(NewAttribute("unit", "K"));
0085   TElement->setAttributeNode(NewAttribute("value", T / kelvin));
0086   element->appendChild(TElement);
0087 }
0088
0089 void PHG4GDMLWriteMaterials::
0090     MEEWrite(xercesc::DOMElement* element, const G4double& MEE)
0091 {
0092   xercesc::DOMElement* PElement = NewElement("MEE");
0093   PElement->setAttributeNode(NewAttribute("unit", "eV"));
0094   PElement->setAttributeNode(NewAttribute("value", MEE / electronvolt));
0095   element->appendChild(PElement);
0096 }
0097
0098 void PHG4GDMLWriteMaterials::
0099     IsotopeWrite(const G4Isotope* const isotopePtr)
0100 {
0101   const G4String name = GenerateName(isotopePtr->GetName(), isotopePtr);
0102
0103   xercesc::DOMElement* isotopeElement = NewElement("isotope");
0104   isotopeElement->setAttributeNode(NewAttribute("name", name));
0105   isotopeElement->setAttributeNode(NewAttribute("N", isotopePtr->GetN()));
0106   isotopeElement->setAttributeNode(NewAttribute("Z", isotopePtr->GetZ()));
0107   materialsElement->appendChild(isotopeElement);
0108   AtomWrite(isotopeElement, isotopePtr->GetA());
0109 }
0110
0111 void PHG4GDMLWriteMaterials::ElementWrite(const G4Element* const elementPtr)
0112 {
0113   const G4String name = GenerateName(elementPtr->GetName(), elementPtr);
0114
0115   xercesc::DOMElement* elementElement = NewElement("element");
0116   elementElement->setAttributeNode(NewAttribute("name", name));
0117
0118   const size_t NumberOfIsotopes = elementPtr->GetNumberOfIsotopes();
0119
0120   if (NumberOfIsotopes > 0)
0121   {
0122     const G4double* RelativeAbundanceVector =
0123         elementPtr->GetRelativeAbundanceVector();
0124     for (size_t i = 0; i < NumberOfIsotopes; i++)
0125     {
0126       G4String fractionref = GenerateName(elementPtr->GetIsotope(i)->GetName(),
0127                                           elementPtr->GetIsotope(i));
0128       xercesc::DOMElement* fractionElement = NewElement("fraction");
0129       fractionElement->setAttributeNode(NewAttribute("n",
0130                                                      RelativeAbundanceVector[i]));
0131       fractionElement->setAttributeNode(NewAttribute("ref", fractionref));
0132       elementElement->appendChild(fractionElement);
0133       AddIsotope(elementPtr->GetIsotope(i));
0134     }
0135   }
0136   else
0137   {
0138     elementElement->setAttributeNode(NewAttribute("Z", elementPtr->GetZ()));
0139     AtomWrite(elementElement, elementPtr->GetA());
0140   }
0141
0142   materialsElement->appendChild(elementElement);
0143   // Append the element AFTER all the possible components are appended!
0144 }
0145
0146 void PHG4GDMLWriteMaterials::MaterialWrite(const G4Material* const materialPtr)
0147 {
0148   G4String state_str("undefined");
0149   const G4State state = materialPtr->GetState();
0150   if (state == kStateSolid)
0151   {
0152     state_str = "solid";
0153   }
0154   else if (state == kStateLiquid)
0155   {
0156     state_str = "liquid";
0157   }
0158   else if (state == kStateGas)
0159   {
0160     state_str = "gas";
0161   }
0162
0163   const G4String name = GenerateName(materialPtr->GetName(), materialPtr);
0164
0165   xercesc::DOMElement* materialElement = NewElement("material");
0166   materialElement->setAttributeNode(NewAttribute("name", name));
0167   materialElement->setAttributeNode(NewAttribute("state", state_str));
0168
0169   // Write any property attached to the material...
0170   //
0171   if (materialPtr->GetMaterialPropertiesTable())
0172   {
0173     PropertyWrite(materialElement, materialPtr);
0174   }
0175
0176   if (materialPtr->GetTemperature() != STP_Temperature)
0177   {
0178     TWrite(materialElement, materialPtr->GetTemperature());
0179   }
0180   if (materialPtr->GetPressure() != STP_Pressure)
0181   {
0182     PWrite(materialElement, materialPtr->GetPressure());
0183   }
0184
0185   // Write Ionisation potential (mean excitation energy)
0186   MEEWrite(materialElement, materialPtr->GetIonisation()->GetMeanExcitationEnergy());
0187
0188   DWrite(materialElement, materialPtr->GetDensity());
0189
0190   const size_t NumberOfElements = materialPtr->GetNumberOfElements();
0191
0192   if ((NumberOfElements > 1) || (materialPtr->GetElement(0) && materialPtr->GetElement(0)->GetNumberOfIsotopes() > 1))
0193   {
0194     const G4double* MassFractionVector = materialPtr->GetFractionVector();
0195
0196     for (size_t i = 0; i < NumberOfElements; i++)
0197     {
0198       const G4String fractionref =
0199           GenerateName(materialPtr->GetElement(i)->GetName(),
0200                        materialPtr->GetElement(i));
0201       xercesc::DOMElement* fractionElement = NewElement("fraction");
0202       fractionElement->setAttributeNode(NewAttribute("n",
0203                                                      MassFractionVector[i]));
0204       fractionElement->setAttributeNode(NewAttribute("ref", fractionref));
0205       materialElement->appendChild(fractionElement);
0206       AddElement(materialPtr->GetElement(i));
0207     }
0208   }
0209   else
0210   {
0211     materialElement->setAttributeNode(NewAttribute("Z", materialPtr->GetZ()));
0212     AtomWrite(materialElement, materialPtr->GetA());
0213   }
0214
0215   // Append the material AFTER all the possible components are appended!
0216   //
0217   materialsElement->appendChild(materialElement);
0218 }
0219
0220 #if G4VERSION_NUMBER >= 1100
0221
0222 // --------------------------------------------------------------------
0223 void PHG4GDMLWriteMaterials::PropertyConstWrite(
0224     const G4String& key, const G4double pval,
0225     const G4MaterialPropertiesTable* ptable)
0226 {
0227   const G4String matrixref = GenerateName(key, ptable);
0228   xercesc::DOMElement* matrixElement = NewElement("matrix");
0229   matrixElement->setAttributeNode(NewAttribute("name", matrixref));
0230   matrixElement->setAttributeNode(NewAttribute("coldim", "1"));
0231   std::ostringstream pvalues;
0232
0233   pvalues << pval;
0234   matrixElement->setAttributeNode(NewAttribute("values", pvalues.str()));
0235
0236   defineElement->appendChild(matrixElement);
0237 }
0238
0239 void PHG4GDMLWriteMaterials::PropertyVectorWrite(
0240     const G4String& key, const G4PhysicsFreeVector* const pvec)
0241 {
0242   for (std::size_t i = 0; i < propertyList.size(); ++i)  // Check if property is
0243   {                                                      // already in the list!
0244     if (propertyList[i] == pvec)
0245     {
0246       return;
0247     }
0248   }
0249   propertyList.push_back(pvec);
0250
0251   const G4String matrixref = GenerateName(key, pvec);
0252   xercesc::DOMElement* matrixElement = NewElement("matrix");
0253   matrixElement->setAttributeNode(NewAttribute("name", matrixref));
0254   matrixElement->setAttributeNode(NewAttribute("coldim", "2"));
0255   std::ostringstream pvalues;
0256   for (std::size_t i = 0; i < pvec->GetVectorLength(); ++i)
0257   {
0258     if (i != 0)
0259     {
0260       pvalues << " ";
0261     }
0262     pvalues << pvec->Energy(i) << " " << (*pvec)[i];
0263   }
0264   matrixElement->setAttributeNode(NewAttribute("values", pvalues.str()));
0265
0266   defineElement->appendChild(matrixElement);
0267 }
0268
0269 void PHG4GDMLWriteMaterials::PropertyWrite(xercesc::DOMElement* matElement,
0270                                            const G4Material* const mat)
0271 {
0272   xercesc::DOMElement* propElement;
0273   G4MaterialPropertiesTable* ptable = mat->GetMaterialPropertiesTable();
0274
0275   auto pvec = ptable->GetProperties();
0276   auto cvec = ptable->GetConstProperties();
0277
0278   for (size_t i = 0; i < pvec.size(); ++i)
0279   {
0280     if (pvec[i] != nullptr)
0281     {
0282       propElement = NewElement("property");
0283       propElement->setAttributeNode(
0284           NewAttribute("name", ptable->GetMaterialPropertyNames()[i]));
0285       propElement->setAttributeNode(NewAttribute(
0286           "ref", GenerateName(ptable->GetMaterialPropertyNames()[i],
0287                               pvec[i])));
0288       PropertyVectorWrite(ptable->GetMaterialPropertyNames()[i],
0289                           pvec[i]);
0290       matElement->appendChild(propElement);
0291     }
0292   }
0293
0294   for (size_t i = 0; i < cvec.size(); ++i)
0295   {
0296     if (cvec[i].second == true)
0297     {
0298       propElement = NewElement("property");
0299       propElement->setAttributeNode(NewAttribute(
0300           "name", ptable->GetMaterialConstPropertyNames()[i]));
0301       propElement->setAttributeNode(NewAttribute(
0302           "ref", GenerateName(ptable->GetMaterialConstPropertyNames()[i],
0303                               ptable)));
0304       PropertyConstWrite(ptable->GetMaterialConstPropertyNames()[i],
0305                          cvec[i].first, ptable);
0306       matElement->appendChild(propElement);
0307     }
0308   }
0309 }
0310
0311 #else   // #if G4VERSION_NUMBER >= 1100
0312 void PHG4GDMLWriteMaterials::PropertyVectorWrite(const G4String& key,
0313                                                  const G4PhysicsOrderedFreeVector* const pvec)
0314 {
0315   const G4String matrixref = GenerateName(key, pvec);
0316   xercesc::DOMElement* matrixElement = NewElement("matrix");
0317   matrixElement->setAttributeNode(NewAttribute("name", matrixref));
0318   matrixElement->setAttributeNode(NewAttribute("coldim", "2"));
0319   std::ostringstream pvalues;
0320   for (size_t i = 0; i < pvec->GetVectorLength(); i++)
0321   {
0322     if (i != 0)
0323     {
0324       pvalues << " ";
0325     }
0326     pvalues << pvec->Energy(i) << " " << (*pvec)[i];
0327   }
0328   matrixElement->setAttributeNode(NewAttribute("values", pvalues.str()));
0329
0330   defineElement->appendChild(matrixElement);
0331 }
0332
0333 void PHG4GDMLWriteMaterials::PropertyWrite(xercesc::DOMElement* matElement,
0334                                            const G4Material* const mat)
0335 {
0336   xercesc::DOMElement* propElement;
0337   G4MaterialPropertiesTable* ptable = mat->GetMaterialPropertiesTable();
0338
0339   const std::map<G4int, G4PhysicsOrderedFreeVector*,
0340                  std::less<G4int> >* pmap = ptable->GetPropertyMap();
0341   const std::map<G4int, G4double,
0342                  std::less<G4int> >* cmap = ptable->GetConstPropertyMap();
0343   std::map<G4int, G4PhysicsOrderedFreeVector*,
0344            std::less<G4int> >::const_iterator mpos;
0345   std::map<G4int, G4double,
0346            std::less<G4int> >::const_iterator cpos;
0347
0348   for (mpos = pmap->begin(); mpos != pmap->end(); ++mpos)
0349   {
0350     propElement = NewElement("property");
0351     propElement->setAttributeNode(NewAttribute("name",
0352                                                ptable->GetMaterialPropertyNames()[mpos->first]));
0353     propElement->setAttributeNode(NewAttribute("ref",
0354                                                GenerateName(ptable->GetMaterialPropertyNames()[mpos->first],
0355                                                             mpos->second)));
0356     if (mpos->second)
0357     {
0358       PropertyVectorWrite(ptable->GetMaterialPropertyNames()[mpos->first],
0359                           mpos->second);
0360       matElement->appendChild(propElement);
0361     }
0362     else
0363     {
0364       G4String warn_message = "Null pointer for material property -" + ptable->GetMaterialPropertyNames()[mpos->first] + "- of material -" + mat->GetName() + "- !";
0365       G4Exception("G4GDMLWriteMaterials::PropertyWrite()", "NullPointer",
0366                   JustWarning, warn_message);
0367       continue;
0368     }
0369   }
0370
0371   for (cpos = cmap->begin(); cpos != cmap->end(); ++cpos)
0372   {
0373     propElement = NewElement("property");
0374     propElement->setAttributeNode(NewAttribute("name",
0375                                                ptable->GetMaterialConstPropertyNames()[cpos->first]));
0376     propElement->setAttributeNode(NewAttribute("ref",
0377                                                ptable->GetMaterialConstPropertyNames()[cpos->first]));
0378     xercesc::DOMElement* constElement = NewElement("constant");
0379     constElement->setAttributeNode(NewAttribute("name",
0380                                                 ptable->GetMaterialConstPropertyNames()[cpos->first]));
0381     constElement->setAttributeNode(NewAttribute("value", cpos->second));
0382     defineElement->appendChild(constElement);
0383     matElement->appendChild(propElement);
0384   }
0385 }
0386 #endif  // #if G4VERSION_NUMBER >= 1100
0387
0388 void PHG4GDMLWriteMaterials::MaterialsWrite(xercesc::DOMElement* element)
0389 {
0390   std::cout << "G4GDML: Writing materials..." << std::endl;
0391
0392   materialsElement = NewElement("materials");
0393   element->appendChild(materialsElement);
0394
0395   isotopeList.clear();
0396   elementList.clear();
0397   materialList.clear();
0398   propertyList.clear();
0399 }
0400
0401 void PHG4GDMLWriteMaterials::AddIsotope(const G4Isotope* const isotopePtr)
0402 {
0403   for (auto & i : isotopeList)  // Check if isotope is
0404   {                                                // already in the list!
0405     if (i == isotopePtr)
0406     {
0407       return;
0408     }
0409   }
0410   isotopeList.push_back(isotopePtr);
0411   IsotopeWrite(isotopePtr);
0412 }
0413
0414 void PHG4GDMLWriteMaterials::AddElement(const G4Element* const elementPtr)
0415 {
0416   for (auto & i : elementList)  // Check if element is
0417   {                                                // already in the list!
0418     if (i == elementPtr)
0419     {
0420       return;
0421     }
0422   }
0423   elementList.push_back(elementPtr);
0424   ElementWrite(elementPtr);
0425 }
0426
0427 void PHG4GDMLWriteMaterials::AddMaterial(const G4Material* const materialPtr)
0428 {
0429   for (auto & i : materialList)  // Check if material is
0430   {                                                 // already in the list!
0431     if (i == materialPtr)
0432     {
0433       return;
0434     }
0435   }
0436   materialList.push_back(materialPtr);
0437   MaterialWrite(materialPtr);
0438 }