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