Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NuclWatcher.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4NuclWatcher.cc 66241 2012-12-13 18:34:42Z gunter $
27 //
28 // 20100202 M. Kelsey -- Move most code here from .hh file, clean up
29 // 20100405 M. Kelsey -- Pass const-ref std::vector<>
30 // 20101010 M. Kelsey -- Migrate to integer A and Z
31 // 20101019 M. Kelsey -- CoVerity report: "!true" should be "!here"
32 
33 #include "G4NuclWatcher.hh"
34 #include "globals.hh"
35 
36 #include <algorithm>
37 #include <vector>
38 #include <cmath>
39 
41  const std::vector<G4double>& expa,
42  const std::vector<G4double>& expcs,
43  const std::vector<G4double>& experr,
44  G4bool check,
45  G4bool nucl)
46  : nuclz(z), izotop_chsq(0.), average_ratio(0.), aver_rat_err(0.),
47  aver_lhood(0.), aver_matched(0.), exper_as(expa), exper_cs(expcs),
48  exper_err(experr), checkable(check), nucleable(nucl) {}
49 
51  const G4double small = 0.001;
52 
53  if (std::abs(z-nuclz) >= small) return;
54 
55  G4bool here = false; // Increment specified nucleus count
56  G4int simulatedAsSize = simulated_as.size();
57  for (G4int i = 0; i<simulatedAsSize && !here; i++) {
58  if (std::abs(simulated_as[i] - a) < small) {
59  simulated_cs[i] += 1.0;
60  here = true; // Terminates loop
61  }
62  }
63 
64  if (!here) { // Nothing found, so add new entry
65  simulated_as.push_back(a);
66  simulated_cs.push_back(1.0);
67  }
68 }
69 
71  G4int simulatedAsSize = simulated_as.size();
72  for(G4int i = 0; i < simulatedAsSize ; i++) {
73  double err = std::sqrt(simulated_cs[i]) / simulated_cs[i];
74 
75  simulated_prob.push_back(simulated_cs[i] / nev);
76  simulated_cs[i] *= csec / nev;
77  simulated_errors.push_back(simulated_cs[i] * err);
78  }
79 }
80 
81 std::pair<G4double, G4double> G4NuclWatcher::getExpCs() const {
82  G4double cs = 0.0;
83  G4double err = 0.0;
84 
85  G4int experAsSize = exper_as.size();
86  for(G4int iz = 0; iz < experAsSize; iz++) {
87  cs += exper_cs[iz];
88  err += exper_err[iz];
89  }
90 
91  return std::pair<G4double, G4double>(cs, err);
92 }
93 
94 std::pair<G4double, G4double> G4NuclWatcher::getInuclCs() const {
95  G4double cs = 0.0;
96  G4double err = 0.0;
97  G4int simulatedAsSize = simulated_as.size();
98  for(G4int iz = 0; iz < simulatedAsSize; iz++) {
99  cs += simulated_cs[iz];
100  err += simulated_errors[iz];
101  }
102 
103  return std::pair<G4double, G4double>(cs, err);
104 }
105 
107  const G4double small = 0.001;
108 
109  G4cout << "\n ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ "
110  << "\n **** izotop Z **** " << nuclz << G4endl;
111 
112  izotop_chsq = 0.0;
113  average_ratio = 0.0;
114  aver_rat_err = 0.0;
115  G4double exp_cs = 0.0;
116  G4double exp_cs_err = 0.0;
117  G4double inucl_cs = 0.0;
118  G4double inucl_cs_err = 0.0;
119  std::vector<G4bool> not_used(simulated_cs.size(), true);
120  G4int nmatched = exper_as.size();
121  G4int nused = simulated_cs.size();
122  G4double lhood = 0.0;
123  G4int experAsSize = exper_as.size();
124 
125  for (G4int iz = 0; iz < experAsSize; iz++) {
126  G4double a = exper_as[iz];
127 
128  exp_cs += exper_cs[iz];
129  exp_cs_err += exper_err[iz];
130 
131  G4bool found = false;
132  G4int simulatedAsSize = simulated_as.size();
133  for (G4int i = 0; i<simulatedAsSize && !found; i++) {
134  if (std::fabs(simulated_as[i] - a) < small) {
135  G4double rat = simulated_cs[i] / exper_cs[iz];
136 
137  lhood += std::log10(rat) * std::log10(rat);
138 
139  G4double rat_err
140  = std::sqrt(simulated_errors[i]*simulated_errors[i] +
141  exper_err[iz]*exper_err[iz] * rat*rat) / exper_cs[iz];
142  average_ratio += rat;
143  aver_rat_err += rat_err;
144 
145  G4cout << " A " << a << " exp.cs " << exper_cs[iz] << " err "
146  << exper_err[iz] << G4endl << " sim. cs " << simulated_cs[i]
147  << " err " << simulated_errors[i] << G4endl
148  << " ratio " << rat << " err " << rat_err << G4endl
149  << " simulated production rate " << simulated_prob[i] << G4endl;
150 
151  not_used[i] = false;
152  izotop_chsq += (rat - 1.0) * (rat - 1.0) / rat_err / rat_err;
153  found = true;
154  nused--;
155  }
156  }
157 
158  if (found) nmatched--;
159  else
160  G4cout << " not found exper.: A " << a << " exp.cs " << exper_cs[iz]
161  << " err " << exper_err[iz] << G4endl;
162  }
163 
164  G4cout << " not found in simulations " << nmatched << G4endl
165  << " not found in exper: " << nused << G4endl;
166 
167  G4int simulatedAsSize = simulated_as.size();
168  for(G4int i = 0; i < simulatedAsSize; i++) {
169  inucl_cs += simulated_cs[i];
170  inucl_cs_err += simulated_errors[i];
171 
172  if(not_used[i])
173  G4cout << " extra simul.: A " << simulated_as[i] << " sim. cs "
174  << simulated_cs[i] << " err " << simulated_errors[i] << G4endl;
175 
176  G4cout << " simulated production rate " << simulated_prob[i] << G4endl;
177  }
178 
179  G4int matched = exper_as.size() - nmatched;
180 
181  if (matched > 0) {
182  aver_lhood = lhood;
183  aver_matched = matched;
184  lhood = std::pow(10.0, std::sqrt(lhood/matched));
185 
186  G4cout << " matched " << matched << " CHSQ " << std::sqrt(izotop_chsq) / matched
187  << G4endl
188  << " raw chsq " << izotop_chsq << G4endl
189  << " average ratio " << average_ratio / matched
190  << " err " << aver_rat_err / matched << G4endl
191  << " lhood " << lhood << G4endl;
192  }
193  else {
194  izotop_chsq = 0.0;
195  aver_lhood = 0.0;
196  }
197 
198  G4cout << " exper. cs " << exp_cs << " err " << exp_cs_err << G4endl
199  << " inucl. cs " << inucl_cs << " err " << inucl_cs_err << G4endl
200  << " ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ "
201  << G4endl;
202 }
G4NuclWatcher(G4int z, const std::vector< G4double > &expa, const std::vector< G4double > &expcs, const std::vector< G4double > &experr, G4bool check, G4bool nucl)
std::pair< G4double, G4double > getInuclCs() const
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
int G4int
Definition: G4Types.hh:78
void setInuclCs(G4double csec, G4int nev)
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
std::pair< G4double, G4double > getExpCs() const
tuple z
Definition: test.py:28
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void watch(G4int a, G4int z)