source: trunk/source/processes/hadronic/models/cascade/cascade/src/G4Analyser.cc @ 819

Last change on this file since 819 was 819, checked in by garnier, 16 years ago

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