Functions | Variables
get_target_mass.C File Reference
#include <iostream>
#include <string>
#include <iomanip>
#include "Conventions/Units.h"
#include "PDG/PDGCodes.h"
#include "PDG/PDGUtils.h"
#include "PDG/PDGLibrary.h"

Go to the source code of this file.

Functions

Int_t set_top_volume (string name)
 
void get_mass (double length_unit, double density_unit)
 
void get_target_mass (string geometry_file, string topvolname="", double length_unit=units::cm, double density_unit=units::g_cm3, Bool_t checkoverlaps=kFALSE)
 
void get_mass (Double_t length_unit, Double_t density_unit)
 

Variables

string gFileName
 

Function Documentation

void get_mass ( double  length_unit,
double  density_unit 
)
void get_mass ( Double_t  length_unit,
Double_t  density_unit 
)

Definition at line 91 of file get_target_mass.C.

92 {
93  //tables of Z and A
94  const Int_t lcin_Z = 150;
95  const Int_t lcin_A = 300;
96 
97  // calc unit conversion factors
98  Double_t density_unit_to_SI = density_unit / units::kg_m3;
99  Double_t length_unit_to_SI = length_unit / units::m;
100  Double_t volume_unit_to_SI = TMath::Power(length_unit_to_SI, 3.);
101 #ifdef _debug_
102  cout << "Input density unit --> kg/m^3 : x" << density_unit_to_SI << endl;
103  cout << "Input length unit --> m : x" << length_unit_to_SI << endl;
104 #endif
105 
106  // get materials in geometry
107  TList *matlist = gGeoManager->GetListOfMaterials();
108  if (!matlist ) {
109  cout << "Null list of materials!" << endl;
110  return;
111  } else {
112 #ifdef _debug_
113  matlist->Print();
114 #endif
115  }
116 
117  int max_idx = 0; // number of mixtures in geometry
118  Int_t nmat = matlist->GetEntries();
119  for( Int_t imat = 0; imat < nmat; imat++ )
120  {
121  Int_t idx = gGeoManager->GetMaterial(imat)->GetIndex();
122  max_idx = TMath::Max(max_idx, idx);
123  }
124 
125  //check if material index is unique
126  Int_t * checkindex = new Int_t[max_idx+1];
127  for( Int_t i = 0; i<max_idx+1; i++ ) checkindex[i] = 0;
128  for( Int_t imat = 0; imat < nmat; imat++ )
129  {
130  if( !checkindex[imat] ) checkindex[imat] = 1;
131  else
132  {
133  cout << "material index is not unique" << endl;
134  return;
135  }
136  }
137 
138 #ifdef _debug_
139  cout << "max_idx = " << max_idx << endl;
140  cout << "nmat = " << nmat << endl;
141 #endif
142 
143  TGeoVolume * topvol = gGeoManager->GetTopVolume(); //get top volume
144  if (!topvol) {
145  cout << "volume does not exist" << endl;
146  return;
147  }
148 
149  TGeoIterator NodeIter(topvol);
150  TGeoNode *node;
151  NodeIter.SetType(0); // include all daughters
152 
153  Double_t * volume = new Double_t[max_idx+1];
154  Double_t * mass = new Double_t[max_idx+1];
155 
156  for( Int_t i = 0; i<max_idx+1; i++ ){ volume[i]=0.; mass[i]=0.; } // IMPORTANT! force empty arrays, allows repated calls without ending ROOT session
157 
158  volume[ topvol->GetMaterial()->GetIndex() ] = topvol->Capacity() * volume_unit_to_SI; //iterator does not include topvolume
159 
160  while ( (node=NodeIter()) )
161  {
162  Int_t momidx = node->GetMotherVolume()->GetMaterial()->GetIndex() ;
163  Int_t idx = node->GetVolume() ->GetMaterial()->GetIndex() ;
164 
165  Double_t node_vol = node->GetVolume()->Capacity() * volume_unit_to_SI;
166 
167  volume[ momidx ] -= node_vol; //substract subvolume from mother
168  volume[ idx ] += node_vol;
169  }
170 
171 
172  Double_t larr_MassIsotopes[lcin_Z][lcin_A] = {0.}; //[Z][A], no map in pure ROOT
173  Double_t larr_VolumeIsotopes[lcin_Z][lcin_A] = {0.}; //[Z][A], no map in pure ROOT
174 
175  for( Int_t i=0; i<gGeoManager->GetListOfMaterials()->GetEntries(); i++ )
176  {
177  TGeoMaterial *lgeo_Mat = gGeoManager->GetMaterial(i);
178  Int_t idx = gGeoManager->GetMaterial(i)->GetIndex();
179 
180  if( lgeo_Mat->IsMixture() )
181  {
182  TGeoMixture * lgeo_Mix = dynamic_cast <TGeoMixture*> ( lgeo_Mat );
183  Int_t lint_Nelements = lgeo_Mix->GetNelements();
184 
185  for ( Int_t j=0; j<lint_Nelements; j++)
186  {
187  Int_t lint_Z = TMath::Nint( (Double_t) lgeo_Mix->GetZmixt()[j] );
188  Int_t lint_A = TMath::Nint( (Double_t) lgeo_Mix->GetAmixt()[j] );
189  Double_t ldou_Fraction = lgeo_Mix->GetWmixt()[j];
190  Double_t ldou_Density = lgeo_Mix->GetDensity() * density_unit_to_SI;
191 
192  larr_MassIsotopes[ lint_Z ][ lint_A ] += volume[idx] * ldou_Fraction * ldou_Density;
193  larr_VolumeIsotopes[ lint_Z ][ lint_A ] += volume[idx] * ldou_Fraction;
194  }
195  }
196  }
197 
198  //
199  // print out volume/mass for each `material'
200  //
201 
202  Double_t ldou_MinimumVolume = 1e-20;
203 
204  cout << endl
205  << " Geometry: \"" << gFileName << "\"" << endl
206  << " TopVolume: \"" << topvol->GetName() << "\""
207  << endl;
208 
209  cout <<endl << "materials:" << endl;
210  cout << setw(5) << "index"
211  << setw(15) << "name"
212  << setprecision(6)
213  << setw(14) << "volume (m^3)"
214  << setw(14) << "mass (kg)"
215  << setw(14) << "mass (%)"
216  << endl;
217 
218  double total_mass_materials = 0;
219  for( Int_t i=0; i<gGeoManager->GetListOfMaterials()->GetEntries(); i++ )
220  {
221  Int_t idx = gGeoManager->GetMaterial(i)->GetIndex();
222  Double_t density = gGeoManager->GetMaterial(i)->GetDensity() * density_unit_to_SI;
223  Double_t mass_material = density * volume[idx];
224  if ( volume[idx] > ldou_MinimumVolume ) {
225  total_mass_materials += mass_material;
226  }
227  }
228 
229 
230  for( Int_t i=0; i<gGeoManager->GetListOfMaterials()->GetEntries(); i++ )
231  {
232  Int_t idx = gGeoManager->GetMaterial(i)->GetIndex();
233  Double_t density = gGeoManager->GetMaterial(i)->GetDensity() * density_unit_to_SI;
234 
235  mass[idx] = density * volume[idx];
236 
237  if( volume[idx] > ldou_MinimumVolume ) {
238  cout << setw(5) << i
239  << setw(15) << gGeoManager->GetMaterial(i)->GetName()
240  << setprecision(6)
241  << setw(14) << volume[idx]
242  << setw(14) << mass[idx]
243  << setw(14) << mass[idx]*100./total_mass_materials
244  << endl;
245  }
246  }
247 
248 
249  //
250  // print out mass contribution for each nuclear target
251  //
252  PDGLibrary* pdglib = PDGLibrary::Instance();
253 
254  cout <<endl << "isotopes:" << endl;
255  cout << setw(4) << "Z"
256  << setw(4) << "A"
257  << setw(14) << "PDG isotope"
258  << setw(6) << " "
259  << setprecision(6)
260  << setw(14) << "volume (m^3)"
261  << setw(14) << "mass (kg)"
262  << setw(14) << "mass (%)"
263  << endl;
264 
265  double total_mass_isotopes = 0;
266  for( Int_t i=0; i<lcin_Z; i++ ) {
267  for( Int_t j=0; j<lcin_A; j++ ) {
268  if( larr_VolumeIsotopes[ i ][ j ] > ldou_MinimumVolume ) {
269  total_mass_isotopes += larr_MassIsotopes[ i ][ j ];
270  }
271  }
272  }
273 
274  for( Int_t i=0; i<lcin_Z; i++ )
275  {
276  for( Int_t j=0; j<lcin_A; j++ )
277  {
278  if( larr_VolumeIsotopes[ i ][ j ] > ldou_MinimumVolume ) {
279  int pdgcode = 1000000000 + i*10000 + j*10;
280  cout << setw(4) << i
281  << setw(4)<< j
282  << setw(14) << pdgcode
283  << setw(6) << pdglib->Find(pdgcode)->GetName()
284  << setprecision(6)
285  << setw(14) << larr_VolumeIsotopes[ i ][ j ]
286  << setw(14) << larr_MassIsotopes[ i ][ j ]
287  << setw(14) << larr_MassIsotopes[ i ][ j ]*100.0/total_mass_isotopes
288  << endl;
289  }
290  else if ( larr_VolumeIsotopes[ i ][ j ] < -ldou_MinimumVolume ) {
291  cout << "negative volume, check geometry " << larr_VolumeIsotopes[ i ][ j ] << endl;
292  }
293  }
294  }
295 
296  cout << endl << " mass totals: " << total_mass_materials << " " << total_mass_isotopes
297  << endl << endl;
298 
299  delete [] volume;
300  delete [] mass;
301 
302 }
genie::units::m
static const double m
Definition: Units.h:79
i
size_t i(0)
cout
the ParameterSet object passed in for the configuration of a destination should be the only source that can affect the behavior of that destination This is to eliminate the dependencies of configuring a destination from multiple mostly from the defaults It suppresses possible glitches about changing the configuration file somewhere outside of a destination segament might still affect the behavior of that destination In the previous configuration for a specific the value of a certain e may come from following and have been suppressed It the configuring ParameterSet object for each destination will be required to carry a parameter list as complete as possible If a parameter still cannot be found in the ParameterSet the configuration code will go look for a hardwired default directly The model is a great simplicity comparing with the previous especially when looking for default values Another great advantage is most of the parameters now have very limited places that allows to appear Usually they can only appear at one certain level in a configuration file For in the old configuring model or in a default ParameterSet object inside of a or in a category or in a severity object This layout of multiple sources for a single parameter brings great confusion in both writing a configuration and in processing the configuration file Under the new the only allowed place for the parameter limit to appear is inside of a category which is inside of a destination object Other improvements simplify the meaning of a destination name In the old a destination name has multiple folds of meanings the e cout and cerr have the special meaning of logging messages to standard output or standard error the name also serves as the output filename if the destination is a file these names are also references to look up for detailed configurations in configuring the MessageFacility The multi purpose of the destination name might cause some unwanted behavior in either writing or parsing the configuration file To amend in the new model the destination name is now merely a name for a which might represent the literal purpose of this or just an id All other meanings of the destinations names now go into the destination ParameterSet as individual such as the type parameter and filename parameter Following is the deatiled rule for the new configuring Everything that is related with MessageFacility configuration must be wrapped in a single ParameterSet object with the name MessageFacility The MessageFacility ParameterSet object contains a series of top level parameters These parameters can be chosen a vector of string listing the name of debug enabled models Or use *to enable debug messages in all modules a vector of string a vector of string a vector of string a ParameterSet object containing the list of all destinations The destinations ParameterSet object is a combination of ParameterSet objects for individual destinations There are two types of destinations that you can insert in the destinations ParameterSet ordinary including cout
Definition: MessageFacilityAsStandAlonePackage.txt:82
density
double density[4]
Definition: LAr.C:78
e
const double e
Definition: gUpMuFluxGen.cxx:165
IMap::volume
static const char * volume[nvol]
Definition: CommonIMapIncludes.h:78
genie::units::kg_m3
static const double kg_m3
Definition: Units.h:153
gFileName
string gFileName
Definition: get_target_mass.C:45
genie::PDGLibrary
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:27
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:63
IMap::matlist
static const char * matlist[nspecialmat]
Definition: CommonIMapIncludes.h:74
void get_target_mass ( string  geometry_file,
string  topvolname = "",
double  length_unit = units::cm,
double  density_unit = units::g_cm3,
Bool_t  checkoverlaps = kFALSE 
)

Definition at line 48 of file get_target_mass.C.

52 {
53  PDGLibrary* pdglib = PDGLibrary::Instance(); // get message out of the way
54 
55  gFileName = geometry_file;
56  // import geometry
57  TGeoManager *tgeo = new TGeoManager("TGeo","TGeo");
58  tgeo->Import( geometry_file.c_str() );
59 
60  // set top volume
61  if( ! set_top_volume(topvolname) ) return;
62 
63  // draw
64  gGeoManager->GetTopVolume()->Draw();
65 
66  // method fails if there are volume overlaps, check it if the user wants to... (takes some time)
67  if(checkoverlaps) {
68  Double_t overlap_acc = 0.05;
69  gGeoManager->CheckOverlaps( overlap_acc );
70  }
71 
72  get_mass(length_unit, density_unit);
73 }
get_mass
void get_mass(double length_unit, double density_unit)
set_top_volume
Int_t set_top_volume(string name)
Definition: get_target_mass.C:76
gFileName
string gFileName
Definition: get_target_mass.C:45
genie::PDGLibrary
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:27
Int_t set_top_volume ( string  name)

Definition at line 76 of file get_target_mass.C.

77 {
78  // no user input, set to overal top volume
79  if( topvolname=="" ) {
80  topvolname = gGeoManager->GetTopVolume()->GetName();
81  }
82  TGeoVolume * topvol = gGeoManager->GetVolume( topvolname.c_str() );
83  if (!topvol) {
84  cout << "top volume does not exist" << endl;
85  return 0;
86  }
87  gGeoManager->SetTopVolume(topvol);
88  return 1;
89 }
cout
the ParameterSet object passed in for the configuration of a destination should be the only source that can affect the behavior of that destination This is to eliminate the dependencies of configuring a destination from multiple mostly from the defaults It suppresses possible glitches about changing the configuration file somewhere outside of a destination segament might still affect the behavior of that destination In the previous configuration for a specific the value of a certain e may come from following and have been suppressed It the configuring ParameterSet object for each destination will be required to carry a parameter list as complete as possible If a parameter still cannot be found in the ParameterSet the configuration code will go look for a hardwired default directly The model is a great simplicity comparing with the previous especially when looking for default values Another great advantage is most of the parameters now have very limited places that allows to appear Usually they can only appear at one certain level in a configuration file For in the old configuring model or in a default ParameterSet object inside of a or in a category or in a severity object This layout of multiple sources for a single parameter brings great confusion in both writing a configuration and in processing the configuration file Under the new the only allowed place for the parameter limit to appear is inside of a category which is inside of a destination object Other improvements simplify the meaning of a destination name In the old a destination name has multiple folds of meanings the e cout and cerr have the special meaning of logging messages to standard output or standard error the name also serves as the output filename if the destination is a file these names are also references to look up for detailed configurations in configuring the MessageFacility The multi purpose of the destination name might cause some unwanted behavior in either writing or parsing the configuration file To amend in the new model the destination name is now merely a name for a which might represent the literal purpose of this or just an id All other meanings of the destinations names now go into the destination ParameterSet as individual such as the type parameter and filename parameter Following is the deatiled rule for the new configuring Everything that is related with MessageFacility configuration must be wrapped in a single ParameterSet object with the name MessageFacility The MessageFacility ParameterSet object contains a series of top level parameters These parameters can be chosen a vector of string listing the name of debug enabled models Or use *to enable debug messages in all modules a vector of string a vector of string a vector of string a ParameterSet object containing the list of all destinations The destinations ParameterSet object is a combination of ParameterSet objects for individual destinations There are two types of destinations that you can insert in the destinations ParameterSet ordinary including cout
Definition: MessageFacilityAsStandAlonePackage.txt:82

Variable Documentation

string gFileName

Definition at line 45 of file get_target_mass.C.

Generated by   doxygen 1.8.11