Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
ana::TargetCount Class Reference

Count number of targets within the main part of the ND (vSuperBlockND) More...

#include "/cvmfs/nova-development.opensciencegrid.org/novasoft/releases/N20-11-30/CAFAna/XSec/TargetCount.h"

Public Member Functions

 TargetCount (int Z=-1)
 
 TargetCount (const TVector3 &min, const TVector3 &max, const int npoints, int Z=-1)
 
double Mass () const
 Mass in kg. More...
 
double NNucleons () const
 Number of nucleons (mass * avogadro's number) More...
 
double NNuclei (int Z=-1, int A=-1) const
 Number of nuclei, optionally restricted to a particular element or isotope. More...
 
double ErrorFraction () const
 Fractional error on mass or number of nucleons. More...
 
std::map< int, double > ZMassMap () const
 Map of mass for each Z. More...
 

Protected Member Functions

int SetupGeoManager ()
 
double CountMass (int Z, TVector3 min, TVector3 max, int npoints, double &errofrac)
 
double CountMass (const TGeoVolume *vol, int Z, double &errorfrac)
 Recurse through the geometry to count the mass. More...
 

Protected Attributes

TGeoManager * fGeoManager
 
double fMass
 
double fErrorFrac
 
double fNNucleons
 
std::map< int, double > fZMassMap
 
std::map< int, std::map< int, double > > fZIsotopeMassMap
 
std::map< std::string, double > kMassError
 

Detailed Description

Count number of targets within the main part of the ND (vSuperBlockND)

Definition at line 10 of file TargetCount.h.

Constructor & Destructor Documentation

ana::TargetCount::TargetCount ( int  Z = -1)
Parameters
ZCount instances of this element. -1 = all elements

Definition at line 212 of file TargetCount.cxx.

References ana::assert(), CountMass(), fErrorFrac, fGeoManager, fMass, fNNucleons, Na, and SetupGeoManager().

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  }
TGeoManager * fGeoManager
Definition: TargetCount.h:41
const int Na
Float_t Z
Definition: plot.C:38
double CountMass(int Z, TVector3 min, TVector3 max, int npoints, double &errofrac)
Definition: TargetCount.cxx:66
assert(nhit_max >=nhit_nbins)
ana::TargetCount::TargetCount ( const TVector3 &  min,
const TVector3 &  max,
const int  npoints,
int  Z = -1 
)

Count the number of nuclear targets of atomic number Z, and inside the fiducial bounds specified by min and max The method works by random-polling and is based on calculateMassesLong function in Geometry service

Parameters
ZCount instances of this element. -1 = all elements
npointsis the number of samples in the random poll This value should be 1e6 or more.

Definition at line 230 of file TargetCount.cxx.

References CountMass(), fErrorFrac, fGeoManager, fMass, fNNucleons, Na, and SetupGeoManager().

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  }
TGeoManager * fGeoManager
Definition: TargetCount.h:41
const int Na
Float_t Z
Definition: plot.C:38
double CountMass(int Z, TVector3 min, TVector3 max, int npoints, double &errofrac)
Definition: TargetCount.cxx:66
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68

Member Function Documentation

double ana::TargetCount::CountMass ( int  Z,
TVector3  min,
TVector3  max,
int  npoints,
double &  errofrac 
)
protected

Definition at line 66 of file TargetCount.cxx.

References om::cout, barite::density, visualisationForPaperMasterPlot::extent, fGeoManager, fZIsotopeMassMap, fZMassMap, MECModelEnuComparisons::i, it, kMassError, m, min(), file_size_ana::node, scale, ana::Sqrt(), string, submit_syst::x, submit_syst::y, test::z, and Z.

Referenced by CountMass(), and TargetCount().

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  }
TGeoManager * fGeoManager
Definition: TargetCount.h:41
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
std::map< std::string, double > kMassError
Definition: TargetCount.h:48
Float_t Z
Definition: plot.C:38
Double_t scale
Definition: plot.C:25
Definition: structs.h:1
z
Definition: test.py:28
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
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
enum BeamMode string
double ana::TargetCount::CountMass ( const TGeoVolume *  vol,
int  Z,
double &  errorfrac 
)
protected

Recurse through the geometry to count the mass.

Definition at line 25 of file TargetCount.cxx.

References CountMass(), d, fGeoManager, MECModelEnuComparisons::i, kMassError, m, file_size_ana::node, string, and w.

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  }
const XML_Char * name
Definition: expat.h:151
TGeoManager * fGeoManager
Definition: TargetCount.h:41
std::map< std::string, double > kMassError
Definition: TargetCount.h:48
Float_t Z
Definition: plot.C:38
double CountMass(int Z, TVector3 min, TVector3 max, int npoints, double &errofrac)
Definition: TargetCount.cxx:66
Float_t d
Definition: plot.C:236
Float_t w
Definition: plot.C:20
Eigen::MatrixXd mat
enum BeamMode string
double ana::TargetCount::ErrorFraction ( ) const
inline

Fractional error on mass or number of nucleons.

Definition at line 35 of file TargetCount.h.

References fErrorFrac.

35 {return fErrorFrac;}
double ana::TargetCount::Mass ( ) const
inline

Mass in kg.

Definition at line 29 of file TargetCount.h.

References fMass.

Referenced by test_flux().

29 {return fMass;}
double ana::TargetCount::NNuclei ( int  Z = -1,
int  A = -1 
) const

Number of nuclei, optionally restricted to a particular element or isotope.

Definition at line 249 of file TargetCount.cxx.

References a, fZIsotopeMassMap, Na, and test::z.

Referenced by NNucleons(), and ana::CrossSectionAnalysis::NucleusCount().

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  }
std::map< int, std::map< int, double > > fZIsotopeMassMap
Definition: TargetCount.h:47
const int Na
Float_t Z
Definition: plot.C:38
const double a
z
Definition: test.py:28
static const double A
Definition: Units.h:82
double ana::TargetCount::NNucleons ( ) const
inline

Number of nucleons (mass * avogadro's number)

Definition at line 31 of file TargetCount.h.

References genie::units::A, fNNucleons, NNuclei(), and Z.

Referenced by CalculateXSec(), ana::DeriveFlux(), GeniePredictionToRoot(), ana::ICrossSectionAnalysis::NTargets(), ana::CrossSectionAnalysis::NucleonCount(), and test_flux().

31 {return fNNucleons;}
int ana::TargetCount::SetupGeoManager ( )
protected

Definition at line 182 of file TargetCount.cxx.

References ana::assert(), fGeoManager, kMassError, path, string, and ana::Wildcard().

Referenced by TargetCount(), and ZMassMap().

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  }
TGeoManager * fGeoManager
Definition: TargetCount.h:41
std::map< std::string, double > kMassError
Definition: TargetCount.h:48
std::vector< std::string > Wildcard(const std::string &wildcardString)
Find files matching a UNIX glob, plus expand environment variables.
Definition: UtilsExt.cxx:268
const std::string path
Definition: plot_BEN.C:43
assert(nhit_max >=nhit_nbins)
enum BeamMode string
std::map<int, double> ana::TargetCount::ZMassMap ( ) const
inline

Map of mass for each Z.

Definition at line 37 of file TargetCount.h.

References fZMassMap, and SetupGeoManager().

37 {return fZMassMap;}
std::map< int, double > fZMassMap
Definition: TargetCount.h:46

Member Data Documentation

double ana::TargetCount::fErrorFrac
protected

Definition at line 43 of file TargetCount.h.

Referenced by ErrorFraction(), and TargetCount().

TGeoManager* ana::TargetCount::fGeoManager
protected

Definition at line 41 of file TargetCount.h.

Referenced by CountMass(), SetupGeoManager(), and TargetCount().

double ana::TargetCount::fMass
protected

Definition at line 42 of file TargetCount.h.

Referenced by Mass(), and TargetCount().

double ana::TargetCount::fNNucleons
protected

Definition at line 44 of file TargetCount.h.

Referenced by NNucleons(), and TargetCount().

std::map<int, std::map<int, double> > ana::TargetCount::fZIsotopeMassMap
protected

Definition at line 47 of file TargetCount.h.

Referenced by CountMass(), and NNuclei().

std::map<int, double> ana::TargetCount::fZMassMap
protected

Definition at line 46 of file TargetCount.h.

Referenced by CountMass(), and ZMassMap().

std::map<std::string, double> ana::TargetCount::kMassError
protected

Definition at line 48 of file TargetCount.h.

Referenced by CountMass(), and SetupGeoManager().


The documentation for this class was generated from the following files: