ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NuclWatcher.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4NuclWatcher.cc
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 //
27 // 20100202 M. Kelsey -- Move most code here from .hh file, clean up
28 // 20100405 M. Kelsey -- Pass const-ref std::vector<>
29 // 20101010 M. Kelsey -- Migrate to integer A and Z
30 // 20101019 M. Kelsey -- CoVerity report: "!true" should be "!here"
31 
32 #include "G4NuclWatcher.hh"
33 #include "globals.hh"
34 
35 #include <algorithm>
36 #include <vector>
37 #include <cmath>
38 
40  const std::vector<G4double>& expa,
41  const std::vector<G4double>& expcs,
42  const std::vector<G4double>& experr,
43  G4bool check,
44  G4bool nucl)
45  : nuclz(z), izotop_chsq(0.), average_ratio(0.), aver_rat_err(0.),
46  aver_lhood(0.), aver_matched(0.), exper_as(expa), exper_cs(expcs),
47  exper_err(experr), checkable(check), nucleable(nucl) {}
48 
50  const G4double small = 0.001;
51 
52  if (std::abs(z-nuclz) >= small) return;
53 
54  G4bool here = false; // Increment specified nucleus count
55  G4int simulatedAsSize = simulated_as.size();
56  for (G4int i = 0; i<simulatedAsSize && !here; i++) {
57  if (std::abs(simulated_as[i] - a) < small) {
58  simulated_cs[i] += 1.0;
59  here = true; // Terminates loop
60  }
61  }
62 
63  if (!here) { // Nothing found, so add new entry
64  simulated_as.push_back(a);
65  simulated_cs.push_back(1.0);
66  }
67 }
68 
70  G4int simulatedAsSize = simulated_as.size();
71  for(G4int i = 0; i < simulatedAsSize ; i++) {
72  double err = std::sqrt(simulated_cs[i]) / simulated_cs[i];
73 
74  simulated_prob.push_back(simulated_cs[i] / nev);
75  simulated_cs[i] *= csec / nev;
76  simulated_errors.push_back(simulated_cs[i] * err);
77  }
78 }
79 
80 std::pair<G4double, G4double> G4NuclWatcher::getExpCs() const {
81  G4double cs = 0.0;
82  G4double err = 0.0;
83 
84  G4int experAsSize = exper_as.size();
85  for(G4int iz = 0; iz < experAsSize; iz++) {
86  cs += exper_cs[iz];
87  err += exper_err[iz];
88  }
89 
90  return std::pair<G4double, G4double>(cs, err);
91 }
92 
93 std::pair<G4double, G4double> G4NuclWatcher::getInuclCs() const {
94  G4double cs = 0.0;
95  G4double err = 0.0;
96  G4int simulatedAsSize = simulated_as.size();
97  for(G4int iz = 0; iz < simulatedAsSize; iz++) {
98  cs += simulated_cs[iz];
99  err += simulated_errors[iz];
100  }
101 
102  return std::pair<G4double, G4double>(cs, err);
103 }
104 
106  const G4double small = 0.001;
107 
108  G4cout << "\n ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ "
109  << "\n **** izotop Z **** " << nuclz << G4endl;
110 
111  izotop_chsq = 0.0;
112  average_ratio = 0.0;
113  aver_rat_err = 0.0;
114  G4double exp_cs = 0.0;
115  G4double exp_cs_err = 0.0;
116  G4double inucl_cs = 0.0;
117  G4double inucl_cs_err = 0.0;
118  std::vector<G4bool> not_used(simulated_cs.size(), true);
119  G4int nmatched = exper_as.size();
120  G4int nused = simulated_cs.size();
121  G4double lhood = 0.0;
122  G4int experAsSize = exper_as.size();
123 
124  for (G4int iz = 0; iz < experAsSize; iz++) {
125  G4double a = exper_as[iz];
126 
127  exp_cs += exper_cs[iz];
128  exp_cs_err += exper_err[iz];
129 
130  G4bool found = false;
131  G4int simulatedAsSize = simulated_as.size();
132  for (G4int i = 0; i<simulatedAsSize && !found; i++) {
133  if (std::fabs(simulated_as[i] - a) < small) {
134  G4double rat = simulated_cs[i] / exper_cs[iz];
135 
136  lhood += std::log10(rat) * std::log10(rat);
137 
138  G4double rat_err
139  = std::sqrt(simulated_errors[i]*simulated_errors[i] +
140  exper_err[iz]*exper_err[iz] * rat*rat) / exper_cs[iz];
141  average_ratio += rat;
142  aver_rat_err += rat_err;
143 
144  G4cout << " A " << a << " exp.cs " << exper_cs[iz] << " err "
145  << exper_err[iz] << G4endl << " sim. cs " << simulated_cs[i]
146  << " err " << simulated_errors[i] << G4endl
147  << " ratio " << rat << " err " << rat_err << G4endl
148  << " simulated production rate " << simulated_prob[i] << G4endl;
149 
150  not_used[i] = false;
151  izotop_chsq += (rat - 1.0) * (rat - 1.0) / rat_err / rat_err;
152  found = true;
153  nused--;
154  }
155  }
156 
157  if (found) nmatched--;
158  else
159  G4cout << " not found exper.: A " << a << " exp.cs " << exper_cs[iz]
160  << " err " << exper_err[iz] << G4endl;
161  }
162 
163  G4cout << " not found in simulations " << nmatched << G4endl
164  << " not found in exper: " << nused << G4endl;
165 
166  G4int simulatedAsSize = simulated_as.size();
167  for(G4int i = 0; i < simulatedAsSize; i++) {
168  inucl_cs += simulated_cs[i];
169  inucl_cs_err += simulated_errors[i];
170 
171  if(not_used[i])
172  G4cout << " extra simul.: A " << simulated_as[i] << " sim. cs "
173  << simulated_cs[i] << " err " << simulated_errors[i] << G4endl;
174 
175  G4cout << " simulated production rate " << simulated_prob[i] << G4endl;
176  }
177 
178  G4int matched = exper_as.size() - nmatched;
179 
180  if (matched > 0) {
181  aver_lhood = lhood;
182  aver_matched = matched;
183  lhood = std::pow(10.0, std::sqrt(lhood/matched));
184 
185  G4cout << " matched " << matched << " CHSQ " << std::sqrt(izotop_chsq) / matched
186  << G4endl
187  << " raw chsq " << izotop_chsq << G4endl
188  << " average ratio " << average_ratio / matched
189  << " err " << aver_rat_err / matched << G4endl
190  << " lhood " << lhood << G4endl;
191  }
192  else {
193  izotop_chsq = 0.0;
194  aver_lhood = 0.0;
195  }
196 
197  G4cout << " exper. cs " << exp_cs << " err " << exp_cs_err << G4endl
198  << " inucl. cs " << inucl_cs << " err " << inucl_cs_err << G4endl
199  << " ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ "
200  << G4endl;
201 }