TargetCount.cxx
Go to the documentation of this file.
2 
4 
5 #include "TGeoManager.h"
6 #include "TGeoMaterial.h"
7 #include "TMath.h"
8 #include "TH1D.h"
9 #include "TCanvas.h"
10 #include "TRandom3.h"
11 
12 #include <cassert>
13 #include <iomanip>
14 #include <stdio.h>
15 #include <iostream>
16 
17 namespace ana
18 {
19  //----------------------------------------------------------------------
20 
21  //----------------------------------------------------------------------
22 
23 
24  /// Recurse through the geometry to count the mass
25  double TargetCount::CountMass(const TGeoVolume* vol, int Z, double &errorfrac)
26  {
27  double mass = 0;
28  double children_capacity = 0;
29 
30  // Recurse through all the children and count their mass. Also keep track
31  // of their volume ("capacity") so we don't double-count our own.
32  for(int d = 0; d < vol->GetNdaughters(); ++d){
33  const TGeoNode* node = vol->GetNode(d);
34  std::string name = node->GetName();
35  name.erase(name.find('_'));
36  TGeoVolume* child_vol = fGeoManager->FindVolumeFast(name.c_str());
37  children_capacity += child_vol->Capacity();
38  mass += CountMass(child_vol, Z, errorfrac);
39  }
40 
41  // Get detailed properties of the material
42  TGeoMaterial* mat = vol->GetMaterial();
43  // How much of this volume to count. If we're asking for all elements
44  // (Z=-1) it's all of it, if we're asking for a specific element it's
45  // zero unless we find some.
46  double w = (Z > 0) ? 0 : 1;
47  for(int i = 0; i < mat->GetNelements(); ++i){
48  double matA, matZ, matW;
49  mat->GetElementProp(matA, matZ, matW, i);
50  // We found the requested element, use its mass fraction
51  if(matZ == Z) w = matW;
52  }
53 
54  // g -> kg
55  double m = .001 * w * (vol->Capacity()-children_capacity) * vol->GetMaterial()->GetDensity();
56  mass += m;
57 
58  errorfrac += m*kMassError[mat->GetName()];
59 
60  return mass;
61  }
62 
63 
64  //----------------------------------------------------------------------
65 
66  double TargetCount::CountMass(int Z, TVector3 min, TVector3 max, int npoints, double& errorfrac)
67  {
68  TVector3 extent = max - min;
69  double scale = 1.0;
70 
71  const double cell_volume = (scale*extent.X())*(scale*extent.Y())*
72  (scale*extent.Z()) / double(npoints);
73 
74  std::map <std::string, double> density;
75  std::map <std::string, double> volume;
76  std::map <std::string, double> mass;
77  std::map <int,double> zerror;
78  fZMassMap.clear();
79  fZIsotopeMassMap.clear();
80 
81  TRandom3 random(0);
82 
83  for(int ipoint = 0; ipoint < npoints; ++ipoint){
84  double x = random.Uniform(min.X(), max.X());
85  double y = random.Uniform(min.Y(), max.Y());
86  double z = random.Uniform(min.Z(), max.Z());
87 
88  const TGeoNode* node = fGeoManager->FindNode(x,y,z);
89  TGeoMaterial* current_material = node->GetMedium()->GetMaterial();
90 
91  const std::string current_material_name = current_material->GetName();
92 
93  if(density.find(current_material_name) == density.end()){
94  density [current_material_name] = current_material->GetDensity(); // in g/cm^3
95  // Initialize these to zero
96  volume [current_material_name] = 0.0;
97  mass [current_material_name] = 0.0;
98 
99  }// end of searching current_material_name
100 
101  if( x > min.X() && x < max.X() &&
102  y > min.Y() && y < max.Y() &&
103  z > min.Z() && z < max.Z() ){
104 
105  volume[current_material_name] += cell_volume;
106  mass [current_material_name] += 0.001 * cell_volume * density[current_material_name];
107 
108  for(int i = 0; i < current_material->GetNelements(); ++i){
109  // atomic mass, number and mass fraction of ith element in
110  // current material
111  double matA, matZ, matW;
112  current_material->GetElementProp(matA, matZ, matW, i);
113  long double m = 0.001 * cell_volume * density[current_material_name]* matW;
114  fZMassMap[matZ] += m;
115  fZIsotopeMassMap[matZ][matA] += m;
116  zerror[matZ] += m * kMassError[current_material_name];
117  }
118  }
119  }// end loop over npoints
120 
121 
122  std::cout<<"\n"<<std::setw(20)<<"Material"
123  <<std::setw(20)<<"Density [g/cm3]"
124  <<std::setw(15)<<"Volume [cm3]"
125  <<std::setw(15)<<"Mass [kg]"
126  <<std::setw(20)<<"Fraction of Total"
127  <<std::setw(20)<<"Error in mass [kg] \n";
128 
129  double totmass = 0;
130  for( auto&& imap : mass)
131  totmass += imap.second;
132 
133  errorfrac = 0;
134  for(std::map<std::string, double>::iterator iter = density.begin(); iter != density.end(); ++iter) {
135 
136  const std::string key_str = iter->first;
137 
138  std::cout<<std::setw(20)<<key_str
139  <<std::setw(20)<<density[key_str]
140  <<std::setw(15)<<volume[key_str]
141  <<std::setw(15)<<mass[key_str]
142  <<std::setw(20)<<mass[key_str]/totmass
143  <<std::setw(20)<<mass[key_str]*kMassError[key_str]<<"\n";
144  errorfrac += (mass[key_str]*kMassError[key_str]*mass[key_str]*kMassError[key_str]);
145  }// end of printing the mass table for each detector material
146 
147  errorfrac = TMath::Sqrt(errorfrac)/totmass;
148  std::cout<<"total mass : "<<totmass<<" +/- "<<errorfrac*totmass<<"("<<errorfrac*100<<"%)\n";
149 
150 
151  std::cout<<"-------------------------------------\n";
152 
153  double masstot = -5;
154  if(Z>0)
155  return fZMassMap[Z];
156  else{
157  for (std::map<int, double>::iterator it=fZMassMap.begin(); it!=fZMassMap.end(); ++it)
158  masstot += it->second;
159  }
160 
161  std::cout<<"\n"<<std::setw(5)<<"Z"
162  <<std::setw(20)<<"Mass [kg]"
163  <<std::setw(20)<<"Fraction of Total"
164  <<std::setw(20)<<"Uncertainty [kg]\n";
165  for(std::map<int, double>::iterator iter = fZMassMap.begin(); iter != fZMassMap.end(); ++iter) {
166  std::cout<<std::setw(5)<<iter->first
167  <<std::setw(20)<<iter->second
168  <<std::setw(20)<<iter->second/masstot
169  <<std::setw(20)<<zerror[iter->first]<<"\n";
170 
171  }
172  double toterr = 0;
173  for(auto& err : zerror)
174  toterr += err.second*err.second;
175  toterr = TMath::Sqrt(toterr);
176  std::cout<<"total mass : "<<masstot<<" +/- "<<toterr<<"("<<toterr*100/masstot<<"%)\n";
177  return masstot;
178  }
179 
180  //----------------------------------------------------------------------
181 
183  {
184  const int old_verbosity = gGeoManager->GetVerboseLevel();
185  gGeoManager->SetVerboseLevel(0);
186 
187  const std::string gdmlName = "neardet-3x3-8block-xtru-vacuum-stagger.gdml";
188 
189  std::vector<std::string> path = Wildcard("$SRT_PRIVATE_CONTEXT/Geometry/gdml/"+gdmlName);
190 
191  if(path.empty()) path = Wildcard("$SRT_PUBLIC_CONTEXT/Geometry/gdml/"+gdmlName);
192  assert(path.size() == 1); // Must find in one or the other
193 
194  fGeoManager = TGeoManager::Import(path[0].c_str());
195 
196  // Error estimates are from Geometry tech note by Matt Strait
197  // NOvA DocDB 15421, Page 25.
198  // errors are in % by mass
199  // The tech note says the scintillator error is at least 0.3%, but we double that,
200  // to account for air bubbles.
201  kMassError["Scintillator"] = 0.006;
202  kMassError["Glue"] = 0.05; // 5%
203  kMassError["PVC"] = 0.016; // 1.6%
204  kMassError["Air"] = 0; // The error in air is folded into scintillator.
205 
206  return old_verbosity;
207 
208  }
209 
210  //----------------------------------------------------------------------
211 
213  {
214 
215  const int old_verbosity = SetupGeoManager();
216  TGeoVolume* vol = fGeoManager->FindVolumeFast("vSuperBlockND");
217  assert(vol);
218 
219  // fMass = vol->Weight(); // Count the whole volume
220  double errorfrac;
221  fMass = CountMass(vol, Z, errorfrac); // Only some subset
222  fErrorFrac = errorfrac;
223  fNNucleons = 1000 * fMass * TMath::Na();
224 
225  fGeoManager->SetVerboseLevel(old_verbosity);
226  }
227 
228  //----------------------------------------------------------------------
229 
230  TargetCount::TargetCount(const TVector3& min,
231  const TVector3& max,
232  const int npoints,
233  int Z)
234  {
235 
236 
237  const int old_verbosity = SetupGeoManager();
238 
239  double errorfrac;
240  fMass = CountMass(Z, min, max, npoints, errorfrac);
241  fErrorFrac = errorfrac;
242  fNNucleons = 1000 * fMass * TMath::Na();
243 
244  fGeoManager->SetVerboseLevel(old_verbosity);
245 
246  }
247 
248  //----------------------------------------------------------------------
249  double TargetCount::NNuclei(int Z, int A) const
250  {
251  double nnuclei = 0;
252  for (const auto & isotopePair : fZIsotopeMassMap)
253  {
254  int z = isotopePair.first;
255  if (Z > 0 && Z != z)
256  continue;
257 
258  for (const auto & massPair : isotopePair.second)
259  {
260  int a = massPair.first;
261  if (A > 0 && A != a)
262  continue;
263 
264  nnuclei += (1000 * massPair.second * TMath::Na()) / a;
265  }
266  }
267 
268  return nnuclei;
269  }
270 
271 }
const XML_Char * name
Definition: expat.h:151
TGeoManager * fGeoManager
Definition: TargetCount.h:41
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
set< int >::iterator it
std::map< int, double > fZMassMap
Definition: TargetCount.h:46
std::map< int, std::map< int, double > > fZIsotopeMassMap
Definition: TargetCount.h:47
const int Na
std::map< std::string, double > kMassError
Definition: TargetCount.h:48
Float_t Z
Definition: plot.C:38
Double_t scale
Definition: plot.C:25
std::vector< std::string > Wildcard(const std::string &wildcardString)
Find files matching a UNIX glob, plus expand environment variables.
Definition: UtilsExt.cxx:268
const double a
double NNuclei(int Z=-1, int A=-1) const
Number of nuclei, optionally restricted to a particular element or isotope.
double CountMass(int Z, TVector3 min, TVector3 max, int npoints, double &errofrac)
Definition: TargetCount.cxx:66
Float_t d
Definition: plot.C:236
Definition: structs.h:1
z
Definition: test.py:28
TargetCount(int Z=-1)
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
const std::string path
Definition: plot_BEN.C:43
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
static const double A
Definition: Units.h:82
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
assert(nhit_max >=nhit_nbins)
tuple density
Definition: barite.py:33
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
Float_t w
Definition: plot.C:20
Eigen::MatrixXd mat