Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Analyser.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: G4Analyser.cc 66241 2012-12-13 18:34:42Z gunter $
27 //
28 // 20100726 M. Kelsey -- Use references for fetched lists
29 // 20101010 M. Kelsey -- Migrate to integer A and Z
30 // 20101019 M. Kelsey -- CoVerity report, unitialized constructor
31 
32 #include "G4Analyser.hh"
33 #include <cmath>
34 #include <iomanip>
35 
37  : verboseLevel(0), eventNumber(0.0), averageMultiplicity(0.0),
38  averageProtonNumber(0.0), averageNeutronNumber(0.0),
39  averagePionNumber(0.0), averageNucleonKinEnergy(0.0),
40  averageProtonKinEnergy(0.0), averageNeutronKinEnergy(0.0),
41  averagePionKinEnergy(0.0), averageExitationEnergy(0.0),
42  averageOutgoingNuclei(0.0), fissy_prob(0.0), averagePionPl(0.0),
43  averagePionMin(0.0), averagePion0(0.0), averageA(0.0), averageZ(0.0),
44  inel_csec(0.0), withNuclei(false) {
45  if (verboseLevel > 3) {
46  G4cout << " >>> G4Analyser::G4Analyser" << G4endl;
47  }
48 }
49 
51 
52  if (verboseLevel > 3) {
53  G4cout << " >>> G4Analyser::setInelCsec" << G4endl;
54  }
55 
56  inel_csec = csec; // mb
57  withNuclei = withn;
58 
59  if (verboseLevel > 3) {
60  G4cout << " total inelastic " << inel_csec << G4endl;
61  }
62 }
63 
64 void G4Analyser::setWatchers(const std::vector<G4NuclWatcher>& watchers) {
65 
66  if (verboseLevel > 3) {
67  G4cout << " >>> G4Analyser::setWatchers" << G4endl;
68  }
69 
70  ana_watchers = watchers;
71  if (verboseLevel > 3) {
72  G4cout << " watchers set " << watchers.size() << G4endl;
73  }
74 }
75 
77 
78  if (verboseLevel > 3) {
79  G4cout << " >>> G4Analyser::try_watchers" << G4endl;
80  }
81 
82  for (G4int iw = 0; iw < G4int(ana_watchers.size()); iw++) {
83 
84  if (if_nucl) {
85 
86  if (ana_watchers[iw].look_forNuclei()) ana_watchers[iw].watch(a, z);
87 
88  } else {
89 
90  if (!ana_watchers[iw].look_forNuclei()) ana_watchers[iw].watch(a, z);
91  }
92  }
93 }
94 
95 
97 
98  if (verboseLevel > 3) {
99  G4cout << " >>> G4Analyser::analyse" << G4endl;
100  }
101 
102  if (withNuclei) {
103  const std::vector<G4InuclNuclei>& nucleus = output.getOutgoingNuclei();
104 
105  // if (nucleus.size() >= 0) {
106  if (nucleus.size() > 0) {
107  G4int nbig = 0;
108  averageOutgoingNuclei += nucleus.size();
109 
110  for (G4int in = 0; in < G4int(nucleus.size()); in++) {
111  averageExitationEnergy += nucleus[in].getExitationEnergy();
112 
113  G4int a = nucleus[in].getA();
114  G4int z = nucleus[in].getZ();
115 
116  if (in == 0) {
117  averageA += a;
118  averageZ += z;
119  };
120 
121  if (a > 10) nbig++;
122  try_watchers(a, z, true);
123  };
124 
125  if (nbig > 1) fissy_prob += 1.0;
126  eventNumber += 1.0;
127  const std::vector<G4InuclElementaryParticle>& particles =
128  output.getOutgoingParticles();
129  averageMultiplicity += particles.size();
130 
131  for (G4int i = 0; i < G4int(particles.size()); i++) {
132  G4int ap = 0;
133  G4int zp = 0;
134 
135  if (particles[i].nucleon()) {
136  averageNucleonKinEnergy += particles[i].getKineticEnergy();
137 
138  if (particles[i].type() == 1) {
139  zp = 1;
140  ap = 1;
141  averageProtonNumber += 1.0;
142  averageProtonKinEnergy += particles[i].getKineticEnergy();
143 
144  } else {
145  ap = 1;
146  zp = 0;
147  averageNeutronNumber += 1.0;
148  averageNeutronKinEnergy += particles[i].getKineticEnergy();
149  };
150 
151  } else if (particles[i].pion()) {
152  averagePionKinEnergy += particles[i].getKineticEnergy();
153  averagePionNumber += 1.0;
154  ap = 0;
155 
156  if (particles[i].type() == 3) {
157  zp = 1;
158  averagePionPl += 1.0;
159 
160  } else if (particles[i].type() == 5) {
161  zp = -1;
162  averagePionMin += 1.0;
163 
164  } else if (particles[i].type() == 7) {
165  zp = 0;
166  averagePion0 += 1.0;
167  };
168  };
169  try_watchers(ap, zp, false);
170  };
171  };
172 
173  } else {
174  eventNumber += 1.0;
175  const std::vector<G4InuclElementaryParticle>& particles =
176  output.getOutgoingParticles();
177  averageMultiplicity += particles.size();
178 
179  for (G4int i = 0; i < G4int(particles.size()); i++) {
180 
181  if (particles[i].nucleon()) {
182  averageNucleonKinEnergy += particles[i].getKineticEnergy();
183 
184  if (particles[i].type() == 1) {
185  averageProtonNumber += 1.0;
186  averageProtonKinEnergy += particles[i].getKineticEnergy();
187 
188  } else {
189  averageNeutronNumber += 1.0;
190  averageNeutronKinEnergy += particles[i].getKineticEnergy();
191  }
192 
193  } else if (particles[i].pion()) {
194  averagePionKinEnergy += particles[i].getKineticEnergy();
195  averagePionNumber += 1.0;
196  }
197  }
198  }
199 }
200 
202 
203  if (verboseLevel > 3) {
204  G4cout << " >>> G4Analyser::printResultsSimple" << G4endl;
205  }
206 
207  G4cout << " Number of events " << G4int(eventNumber + 0.1) << G4endl
208  << " average multiplicity " << averageMultiplicity / eventNumber << G4endl
209  << " average proton number " << averageProtonNumber / eventNumber << G4endl
210  << " average neutron number " << averageNeutronNumber / eventNumber << G4endl
211  << " average nucleon Ekin " << averageNucleonKinEnergy /
212  (averageProtonNumber + averageNeutronNumber) << G4endl
213  << " average proton Ekin " << averageProtonKinEnergy / (averageProtonNumber +
214  1.0e-10) << G4endl
215  << " average neutron Ekin " << averageNeutronKinEnergy / (averageNeutronNumber +
216  1.0e-10) << G4endl
217  << " average pion number " << averagePionNumber / eventNumber << G4endl
218  << " average pion Ekin " << averagePionKinEnergy / (averagePionNumber +
219  1.0e-10) << G4endl;
220  if (withNuclei) {
221  G4cout
222  << " average Exitation Energy " <<
223  averageExitationEnergy / averageOutgoingNuclei << G4endl
224  << " average num of fragments " << averageOutgoingNuclei / eventNumber << G4endl;
225  G4cout << " fission prob. " << fissy_prob / eventNumber << " c.sec " <<
226  inel_csec * fissy_prob / eventNumber << G4endl;
227  }
228 }
229 
231 
232  if (verboseLevel > 3) {
233  G4cout << " >>> G4Analyser::printResults" << G4endl;
234  }
235 
236  G4cout << " Number of events " << G4int(eventNumber + 0.1) << G4endl
237  << " average multiplicity " << averageMultiplicity / eventNumber << G4endl
238  << " average proton number " << averageProtonNumber / eventNumber << G4endl
239  << " average neutron number " << averageNeutronNumber / eventNumber << G4endl
240  << " average nucleon Ekin " << averageNucleonKinEnergy /
241  (averageProtonNumber + averageNeutronNumber) << G4endl
242  << " average proton Ekin " << averageProtonKinEnergy / (averageProtonNumber +
243  1.0e-10) << G4endl
244  << " average neutron Ekin " << averageNeutronKinEnergy / (averageNeutronNumber +
245  1.0e-10) << G4endl
246  << " average pion number " << averagePionNumber / eventNumber << G4endl
247  << " average pion Ekin " << averagePionKinEnergy / (averagePionNumber +
248  1.0e-10) << G4endl
249  << " average pi+ " << averagePionPl / eventNumber << G4endl
250  << " average pi- " << averagePionMin / eventNumber << G4endl
251  << " average pi0 " << averagePion0 / eventNumber << G4endl;
252 
253  if (withNuclei) {
254  G4cout
255  << " average A " << averageA / eventNumber << G4endl
256  << " average Z " << averageZ / eventNumber << G4endl
257  << " average Exitation Energy " <<
258  averageExitationEnergy / averageOutgoingNuclei << G4endl
259  << " average num of fragments " << averageOutgoingNuclei / eventNumber << G4endl;
260  G4cout << " fission prob. " << fissy_prob / eventNumber << " c.sec " <<
261  inel_csec * fissy_prob / eventNumber << G4endl;
263  }
264 }
265 
267 
268  if (verboseLevel > 3) {
269  G4cout << " >>> G4Analyser::handleWatcherStatistics" << G4endl;
270  }
271 
272  // const G4double small = 1.0e-10;
273 
274  if (verboseLevel > 3) {
275  G4cout << " >>>Izotop analysis:" << G4endl;
276  };
277 
278  G4double fgr = 0.0;
279  G4double averat = 0.0;
280  G4double ave_err = 0.0;
281  G4double gl_chsq = 0.0;
282  G4double tot_exper = 0.0;
283  G4double tot_exper_err = 0.0;
284  G4double tot_inucl = 0.0;
285  G4double tot_inucl_err = 0.0;
286  G4double checked = 0.0;
287 
288  for (G4int iw = 0; iw < G4int(ana_watchers.size()); iw++) {
289  ana_watchers[iw].setInuclCs(inel_csec, G4int(eventNumber));
290  ana_watchers[iw].print();
291 
292  if (ana_watchers[iw].to_check()) {
293  std::pair<G4double, G4double> rat_err = ana_watchers[iw].getAverageRatio();
294  averat += rat_err.first;
295  ave_err += rat_err.second;
296  gl_chsq += ana_watchers[iw].getChsq();
297  std::pair<G4double, G4double> cs_err = ana_watchers[iw].getExpCs();
298  tot_exper += cs_err.first;
299  tot_exper_err += cs_err.second;
300  std::pair<G4double, G4double> inucl_cs_err = ana_watchers[iw].getInuclCs();
301  tot_inucl += inucl_cs_err.first;
302  tot_inucl_err += inucl_cs_err.second;
303  G4double iz_checked = ana_watchers[iw].getNmatched();
304 
305  if (iz_checked > 0.0) {
306  fgr += ana_watchers[iw].getLhood();
307  checked += iz_checked;
308  };
309  };
310  };
311 
312  if (checked > 0.0) {
313  gl_chsq = std::sqrt(gl_chsq) / checked;
314  averat /= checked;
315  ave_err /= checked;
316  fgr = std::pow(10.0, std::sqrt(fgr / checked));
317  };
318 
319  if (verboseLevel > 3) {
320  G4cout << " total exper c.s. " << tot_exper << " err " << tot_exper_err <<
321  " tot inucl c.s. " << tot_inucl << " err " << tot_inucl_err << G4endl;
322  G4cout << " checked total " << checked << " lhood " << fgr << G4endl
323  << " average ratio " << averat << " err " << ave_err << G4endl
324  << " global chsq " << gl_chsq << G4endl;
325  }
326 }
327 
329 
330  if (verboseLevel > 3) {
331  G4cout << " >>> G4Analyser::printResultsNtuple" << G4endl;
332  }
333 
334  // Create one line of ACII data.
335  // Several runs should create ntuple for data-analysis
336  G4cout <<
337  std::setw(15) << int(eventNumber + 0.1) <<
338  std::setw(15) << averageMultiplicity / eventNumber <<
339  std::setw(15) << averageProtonNumber / eventNumber <<
340  std::setw(15) << averageNeutronNumber / eventNumber << " " <<
341  std::setw(15) << averageNucleonKinEnergy / (averageProtonNumber + averageNeutronNumber) << " " <<
342  std::setw(15) << averageProtonKinEnergy / (averageProtonNumber + 1.0e-10) << " " <<
343  std::setw(15) << averageNeutronKinEnergy / (averageNeutronNumber + 1.0e-10) << " " <<
344  std::setw(15) << averagePionNumber / eventNumber << " " <<
345  std::setw(15) << averagePionKinEnergy / (averagePionNumber + 1.0e-10) << G4endl;
346 }
G4bool pion(G4int ityp)
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
void analyse(const G4CollisionOutput &output)
Definition: G4Analyser.cc:96
void handleWatcherStatistics()
Definition: G4Analyser.cc:266
void printResultsNtuple()
Definition: G4Analyser.cc:328
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
void try_watchers(G4int a, G4int z, G4bool if_nucl)
Definition: G4Analyser.cc:76
G4bool nucleon(G4int ityp)
void printResultsSimple()
Definition: G4Analyser.cc:201
bool G4bool
Definition: G4Types.hh:79
typedef int(XMLCALL *XML_NotStandaloneHandler)(void *userData)
void printResults()
Definition: G4Analyser.cc:230
const std::vector< G4InuclNuclei > & getOutgoingNuclei() const
const std::vector< G4InuclElementaryParticle > & getOutgoingParticles() const
tuple z
Definition: test.py:28
void setInelCsec(G4double csec, G4bool withn)
Definition: G4Analyser.cc:50
#define G4endl
Definition: G4ios.hh:61
void setWatchers(const std::vector< G4NuclWatcher > &watchers)
Definition: G4Analyser.cc:64
double G4double
Definition: G4Types.hh:76