/*****************************************************************************
* Copyright (C) 2009-2014 this file is part of the NPTool Project *
* *
* For the licensing terms see $NPTOOL/Licence/NPTool_Licence *
* For the list of contributors see $NPTOOL/Licence/Contributors *
*****************************************************************************/
/*****************************************************************************
* Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr *
* *
* Creation Date : October 2014 *
* Last update : *
*---------------------------------------------------------------------------*
* Decription: *
* This singleton class contain a librairy of material available for *
* the geometry. Instantiate the needed material on demand, and generate the *
* associate DEDX table. *
* *
*---------------------------------------------------------------------------*
* Comment: *
* *
*****************************************************************************/
// NPS
#include "MaterialManager.hh"
// Geant4
#include "G4Element.hh"
#include "G4EmCalculator.hh"
#include "G4Box.hh"
#include "G4PVPlacement.hh"
#include "G4VisAttributes.hh"
#include "G4NistManager.hh"
// STL
#include <iostream>
#include <string>
using namespace std;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
MaterialManager* MaterialManager::instance = 0 ;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
MaterialManager* MaterialManager::getInstance(){
if(instance == 0) instance = new MaterialManager();
return instance;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
MaterialManager::~MaterialManager(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
MaterialManager::MaterialManager(){
m_D = NULL;
m_T = NULL;
m_He3 = NULL;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MaterialManager::Destroy(){
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MaterialManager::ClearMaterialLibrary(){
map<string,G4Material*>::iterator it;
for(it = m_Material.begin() ; it != m_Material.end() ; it++){
delete it->second;
}
m_Material.clear();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4Material* MaterialManager::GetMaterialFromLibrary(string Name){
// Search if the material is already instantiate
map<string,G4Material*>::iterator it;
it = m_Material.find(Name);
// The element is not found
if(it==m_Material.end()){
//Usual compound
if(Name == "Vacuum"){
G4Material* material = new G4Material(Name, 0.000000001 * mg / cm3,2);
material->AddElement(GetElementFromLibrary("N"),7);
material->AddElement(GetElementFromLibrary("O"),3);
m_Material[Name]=material;
return material;
}
else if(Name == "PCB"){
// Actually taken value fron Epoxy
G4Material* material = new G4Material(Name, 1.85 * g / cm3,3);
material->AddElement(GetElementFromLibrary("H"),.475);
material->AddElement(GetElementFromLibrary("C"),.45);
material->AddElement(GetElementFromLibrary("O"),.075);
m_Material[Name]=material;
return material;
}
else if(Name == "Mylar"){
G4Material* material = new G4Material(Name, 1.397*g/cm3,3);
material->AddElement(GetElementFromLibrary("H"),8);
material->AddElement(GetElementFromLibrary("C"),10);
material->AddElement(GetElementFromLibrary("O"),4);
m_Material[Name]=material;
return material;
}
else if(Name == "Kapton"){
G4Material* material = new G4Material(Name, 1.39*g/cm3,3);
material->AddElement(GetElementFromLibrary("H"),4);
material->AddElement(GetElementFromLibrary("C"),5);
material->AddElement(GetElementFromLibrary("O"),2);
m_Material[Name]=material;
return material;
}
else if(Name == "Havar"){
G4Material* material = new G4Material(Name, 8.3*g / cm3,5);
material->AddElement(GetElementFromLibrary("Co"),42);
material->AddElement(GetElementFromLibrary("Cr"),20);
material->AddElement(GetElementFromLibrary("Ni"),13);
material->AddElement(GetElementFromLibrary("Fe"),19);
material->AddElement(GetElementFromLibrary("W"),1);
m_Material[Name]=material;
return material;
}
// Usual Target
else if(Name == "CD2"){
G4Material* material = new G4Material(Name, 1.06*g/cm3,2);
material->AddElement(GetElementFromLibrary("C"),1);
material->AddElement(GetElementFromLibrary("D"),2);
m_Material[Name]=material;
return material;
}
else if(Name == "CH2"){
G4Material* material = new G4Material(Name, 0.93*g/cm3,2);
material->AddElement(GetElementFromLibrary("C"),1);
material->AddElement(GetElementFromLibrary("D"),2);
m_Material[Name]=material;
return material;
}
else if(Name == "Cu"){
G4Material* material = new G4Material(Name, 8.96*g/cm3,1);
material->AddElement(GetElementFromLibrary("Cu"),1);
m_Material[Name]=material;
return material;
}
else if(Name == "C"){ // Graphite
G4Material* material = new G4Material(Name, 2.267*g/cm3,1);
material->AddElement(GetElementFromLibrary("C"),1);
m_Material[Name]=material;
return material;
}
else if(Name == "Pb"){
G4Material* material = new G4Material(Name, 11.342*g / cm3,1);
material->AddElement(GetElementFromLibrary("Pb"),1);
m_Material[Name]=material;
return material;
}
else if(Name == "D2_solid"){ // density should be checked
G4Material* material = new G4Material(Name, 0.0715*g/cm3,1);
material->AddElement(GetElementFromLibrary("D"),2);
m_Material[Name]=material;
return material;
}
else if(Name == "H2_solid"){
G4Material* material = new G4Material(Name, 0.0715*g/cm3,1);
material->AddElement(GetElementFromLibrary("H"),2);
m_Material[Name]=material;
return material;
}
// Usual detector material
else if(Name == "Si"){
G4Material* material = new G4Material(Name,2.321*g/cm3 ,1);
material->AddElement(GetElementFromLibrary("Si"),1);
m_Material[Name]=material;
return material;
}
else if(Name == "Ge"){
G4Material* material = new G4Material(Name,5.323*g/cm3 ,1);
material->AddElement(GetElementFromLibrary("Ge"),1);
m_Material[Name]=material;
return material;
}
else if(Name == "BC400"){
G4Material* material = new G4Material(Name,1.032*g/cm3 ,2);
material->AddElement(GetElementFromLibrary("H"),10);
material->AddElement(GetElementFromLibrary("C"),9);
m_Material[Name]=material;
return material;
}
else if(Name == "NaI"){
G4Material* material = new G4Material(Name, 3.67*g/cm3,2);
material->AddElement(GetElementFromLibrary("Na"),1);
material->AddElement(GetElementFromLibrary("I"),1);
m_Material[Name]=material;
return material;
}
else if(Name == "CsI"){
G4Material* material = new G4Material(Name, 4.51*g/cm3,2);
material->AddElement(GetElementFromLibrary("Cs"),1);
material->AddElement(GetElementFromLibrary("I"),1);
m_Material[Name]=material;
return material;
}
else if(Name == "LaBr3"){
G4Material* material = new G4Material(Name, 5.29*g/cm3 , 2);
material->AddElement(GetElementFromLibrary("La"),1);
material->AddElement(GetElementFromLibrary("Br"),3);
m_Material[Name]=material;
return material;
}
else if(Name == "BaF2"){
G4Material* material = new G4Material(Name, 4.89*g/cm3 , 2);
material->AddElement(GetElementFromLibrary("Ba"),1);
material->AddElement(GetElementFromLibrary("F"),2);
m_Material[Name]=material;
return material;
}
// Misc
else if(Name == "Al"){
G4Material* material = new G4Material(Name, 2.702*g/cm3,1);
material->AddElement(GetElementFromLibrary("Al"),1);
m_Material[Name]=material;
return material;
}
else if(Name == "Fe"){
G4Material* material = new G4Material(Name, 7.874*g/cm3,1);
material->AddElement(GetElementFromLibrary("Fe"),1);
m_Material[Name]=material;
return material;
}
else{
G4cout << "ERROR: Material requested \""<< Name <<"\" is not available in the Material Librairy" << G4endl;
exit(1);
}
}
else return it->second;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MaterialManager::AddMaterialToLibrary(G4Material* material){
m_Material[material->GetName()]=material;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4Element* MaterialManager::GetElementFromLibrary(string Name){
if(Name=="D" || Name=="d"){
if(!m_D){
m_D = new G4Element(Name.c_str(), Name.c_str(),1);
G4Isotope* isotope = new G4Isotope(Name.c_str(), 1,2,2.01410178*g/mole);
m_D->AddIsotope(isotope,1);
}
return m_D;
}
else if(Name=="T" || Name=="t"){
if(!m_T){
m_T = new G4Element(Name.c_str(), Name.c_str(),1);
G4Isotope* isotope = new G4Isotope(Name.c_str(), 1,3,3.0160492*g/mole);
m_T->AddIsotope(isotope,1);
}
return m_T;
}
else if(Name=="He3" || Name=="3He"){
if(!m_He3){
m_He3 = new G4Element(Name.c_str(), Name.c_str(),1);
G4Isotope* isotope = new G4Isotope(Name.c_str(), 2,1,3.0160293*g/mole);
m_He3->AddIsotope(isotope,1);
}
return m_He3;
}
G4NistManager* man = G4NistManager::Instance();
return man->FindOrBuildElement(Name.c_str());
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Generate a DEDX file table using the material used in the geometry
void MaterialManager::WriteDEDXTable(G4ParticleDefinition* Particle ,G4double Emin,G4double Emax){
map<string,G4Material*>::iterator it;
for(it = m_Material.begin() ; it != m_Material.end() ; it++){
// Opening hte output file
G4String GlobalPath = getenv("NPTOOL");
G4String Path = GlobalPath + "/Inputs/EnergyLoss/" + Particle->GetParticleName() + "_" + it->second->GetName() + ".G4table";
ofstream File ;
File.open(Path) ;
if(!File) return ;
File << "Table from Geant4 generate using NPSimulation \t"
<< "Particle: " << Particle->GetParticleName() << "\tMaterial: " << it->second->GetName() << G4endl ;
// G4cout << Particle->GetParticleName() << "\tMaterial: " << it->second->GetName() <<endl;
G4EmCalculator emCalculator;
G4double dedx ;
// Tipical Range needed, if Emax is larger, then adapted
if(Emax < 1*TeV) Emax = 1*TeV;
double step = 1*keV;
double before = 0 ;
for (G4double E=Emin; E < Emax; E+=step){
dedx = emCalculator.ComputeTotalDEDX(E, Particle, it->second)/(MeV/micrometer);
if(before){
if(abs(before-dedx)/abs(before)<0.01) step*=10;
}
before = dedx;
File << E/MeV << "\t" << dedx << G4endl ;
}
File.close();
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void MaterialManager::CreateSampleVolumes(G4LogicalVolume* world_log){
// Crate a micrometer big cube for each material
G4double SampleSize = 10 * cm;
G4double WorldSize = 10.0 * m ;
G4Box* sample_box = new G4Box("sample_box",SampleSize ,SampleSize ,SampleSize);
G4int i = 1;
G4double Coord1 = WorldSize-SampleSize;
G4double Coord2 = 0 ;
map<string,G4Material*>::iterator it;
for(it = m_Material.begin() ; it != m_Material.end() ; it++){
G4LogicalVolume* sample_log = new G4LogicalVolume(sample_box, it->second, "sample_log", 0, 0, 0);
sample_log->SetVisAttributes(G4VisAttributes::Invisible);
Coord2 = WorldSize-i*SampleSize;
i++;
new G4PVPlacement(0, G4ThreeVector(Coord1,Coord2,-Coord1), sample_log, "sample", world_log, false, 0);
}
}