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