get_target_mass.C
Go to the documentation of this file.
1 //****************************************************************************
2 //*
3 //* get_target_mass()
4 //*
5 //* Reads-in a ROOT geometry and prints-out its mass composition.
6 //*
7 //* Input arguments:
8 //* - geometry_file
9 //* - topvolname (default: "")
10 //* - length_unit (default: units::cm)
11 //* - density_unit (default: units::g_cm3)
12 //* - checkoverlaps (default: kFALSE)
13 //*
14 //*
15 //* T. Ferber (University Hamburg)
16 //* Luruper Chaussee 149
17 //* 20259 Hamburg
18 //* torben.ferber@desy.de
19 //*
20 //* @ Aug 30, 2010 - TF
21 //* added list of isotope masses, added index check, added arrary initialization, corrected iterator initialization
22 //*
23 //****************************************************************************
24 
25 
26 #include <iostream>
27 #include <string>
28 #include <iostream>
29 #include <iomanip>
30 
31 #include "Conventions/Units.h"
32 
33 #include "PDG/PDGCodes.h"
34 #include "PDG/PDGUtils.h"
35 #include "PDG/PDGLibrary.h"
36 
37 // #define _debug_
38 
39 using namespace std;
40 using namespace genie;
41 
42 Int_t set_top_volume(string name);
43 void get_mass (double length_unit, double density_unit);
44 
45 string gFileName;
46 
47 //____________________________________________________________________________
49  string geometry_file, string topvolname = "",
50  double length_unit=units::cm, double density_unit=units::g_cm3,
51  Bool_t checkoverlaps = kFALSE)
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 }
74 
75 //____________________________________________________________________________
76 Int_t set_top_volume(string topvolname)
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 }
90 //____________________________________________________________________________
91 void get_mass(Double_t length_unit, Double_t density_unit)
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 }
303 //____________________________________________________________________________
304 
const XML_Char * name
Definition: expat.h:151
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
static const double g_cm3
Definition: Units.h:154
Definition: structs.h:1
static const double kg_m3
Definition: Units.h:153
const double j
Definition: BetheBloch.cxx:29
void get_mass(double length_unit, double density_unit)
OStream cout
Definition: OStream.cxx:6
Int_t set_top_volume(string name)
string gFileName
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:30
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
void get_target_mass(string geometry_file, string topvolname="", double length_unit=units::cm, double density_unit=units::g_cm3, Bool_t checkoverlaps=kFALSE)
tuple density
Definition: barite.py:33
Float_t e
Definition: plot.C:35
static constexpr Double_t cm
Definition: Munits.h:140