Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/hadronic/Hadr01/src/HistoManager.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/extended/hadronic/Hadr01/src/HistoManager.cc (Version 11.3.0) and /examples/extended/hadronic/Hadr01/src/HistoManager.cc (Version 9.6)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 /// \file hadronic/Hadr01/src/HistoManager.cc      26 /// \file hadronic/Hadr01/src/HistoManager.cc
 27 /// \brief Implementation of the HistoManager      27 /// \brief Implementation of the HistoManager class
 28 //                                                 28 //
                                                   >>  29 // $Id$
                                                   >>  30 //
 29 //--------------------------------------------     31 //---------------------------------------------------------------------------
 30 //                                                 32 //
 31 // ClassName:   HistoManager                       33 // ClassName:   HistoManager
 32 //                                                 34 //
 33 //                                                 35 //
 34 // Author:      V.Ivanchenko 30/01/01              36 // Author:      V.Ivanchenko 30/01/01
 35 //                                                 37 //
 36 // Modified:                                       38 // Modified:
 37 // 04.06.2006 Adoptation of Hadr01 (V.Ivanchen     39 // 04.06.2006 Adoptation of Hadr01 (V.Ivanchenko)
 38 // 16.11.2006 Add beamFlag (V.Ivanchenko)          40 // 16.11.2006 Add beamFlag (V.Ivanchenko)
 39 //                                                 41 //
 40 //--------------------------------------------     42 //----------------------------------------------------------------------------
 41 //                                                 43 //
 42                                                    44 
 43 //....oooOO0OOooo........oooOO0OOooo........oo     45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 44 //....oooOO0OOooo........oooOO0OOooo........oo     46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 45                                                    47 
 46 #include "HistoManager.hh"                         48 #include "HistoManager.hh"
 47                                                <<  49 #include "G4Track.hh"
 48 #include "Histo.hh"                            <<  50 #include "G4Step.hh"
 49                                                <<  51 #include "G4UnitsTable.hh"
 50 #include "G4Alpha.hh"                          <<  52 #include "G4Neutron.hh"
                                                   >>  53 #include "G4Proton.hh"
 51 #include "G4AntiProton.hh"                         54 #include "G4AntiProton.hh"
 52 #include "G4Deuteron.hh"                       << 
 53 #include "G4Electron.hh"                       << 
 54 #include "G4Gamma.hh"                              55 #include "G4Gamma.hh"
 55 #include "G4He3.hh"                            <<  56 #include "G4Electron.hh"
 56 #include "G4KaonMinus.hh"                      <<  57 #include "G4Positron.hh"
 57 #include "G4KaonPlus.hh"                       << 
 58 #include "G4KaonZeroLong.hh"                   << 
 59 #include "G4KaonZeroShort.hh"                  << 
 60 #include "G4MuonMinus.hh"                      << 
 61 #include "G4MuonPlus.hh"                           58 #include "G4MuonPlus.hh"
 62 #include "G4Neutron.hh"                        <<  59 #include "G4MuonMinus.hh"
 63 #include "G4PhysicalConstants.hh"              << 
 64 #include "G4PionMinus.hh"                      << 
 65 #include "G4PionPlus.hh"                           60 #include "G4PionPlus.hh"
                                                   >>  61 #include "G4PionMinus.hh"
 66 #include "G4PionZero.hh"                           62 #include "G4PionZero.hh"
 67 #include "G4Positron.hh"                       <<  63 #include "G4KaonPlus.hh"
 68 #include "G4Proton.hh"                         <<  64 #include "G4KaonMinus.hh"
 69 #include "G4Step.hh"                           <<  65 #include "G4KaonZeroShort.hh"
 70 #include "G4SystemOfUnits.hh"                  <<  66 #include "G4KaonZeroLong.hh"
 71 #include "G4Track.hh"                          <<  67 #include "G4Deuteron.hh"
 72 #include "G4Triton.hh"                             68 #include "G4Triton.hh"
 73 #include "G4UnitsTable.hh"                     <<  69 #include "G4He3.hh"
                                                   >>  70 #include "G4Alpha.hh"
                                                   >>  71 #include "Histo.hh"
 74 #include "globals.hh"                              72 #include "globals.hh"
                                                   >>  73 #include "G4SystemOfUnits.hh"
 75                                                    74 
 76 //....oooOO0OOooo........oooOO0OOooo........oo     75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 77                                                    76 
 78 HistoManager* HistoManager::fManager = nullptr <<  77 HistoManager* HistoManager::fManager = 0;
 79                                                    78 
 80 //....oooOO0OOooo........oooOO0OOooo........oo     79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 81                                                    80 
 82 HistoManager* HistoManager::GetPointer()           81 HistoManager* HistoManager::GetPointer()
 83 {                                                  82 {
 84   if (nullptr == fManager) {                   <<  83   if(!fManager) {
 85     static HistoManager manager;                   84     static HistoManager manager;
 86     fManager = &manager;                           85     fManager = &manager;
 87   }                                                86   }
 88   return fManager;                                 87   return fManager;
 89 }                                                  88 }
 90                                                    89 
 91 //....oooOO0OOooo........oooOO0OOooo........oo     90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 92                                                    91 
 93 HistoManager::HistoManager()                       92 HistoManager::HistoManager()
 94   : fPrimaryDef(nullptr),                      << 
 95     fRadius(10 * CLHEP::cm),                   << 
 96     fLength(300. * CLHEP::mm),                 << 
 97     fEdepMax(1.0 * CLHEP::GeV)                 << 
 98 {                                                  93 {
 99   fHisto = new Histo();                        <<  94   fVerbose=  0;
                                                   >>  95   fNSlices   = 300;
                                                   >>  96   fNBinsE    = 100;
                                                   >>  97   fNHisto    = 25;
                                                   >>  98   fLength    = 300.*mm;
                                                   >>  99   fEdepMax   = 1.0*GeV;
                                                   >> 100   fBeamFlag  = true;
                                                   >> 101   fHistoBooked = false;
                                                   >> 102   fHisto     = new Histo();
100   fHisto->SetVerbose(fVerbose);                   103   fHisto->SetVerbose(fVerbose);
101   fNeutron = G4Neutron::Neutron();             << 104   fNeutron   = G4Neutron::Neutron();
                                                   >> 105   fPrimaryKineticEnergy = 0.0;
                                                   >> 106   fPrimaryDef = 0;
102 }                                                 107 }
103                                                   108 
104 //....oooOO0OOooo........oooOO0OOooo........oo    109 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105                                                   110 
106 HistoManager::~HistoManager()                     111 HistoManager::~HistoManager()
107 {                                                 112 {
108   delete fHisto;                                  113   delete fHisto;
109 }                                                 114 }
110                                                   115 
111 //....oooOO0OOooo........oooOO0OOooo........oo    116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
112                                                   117 
113 void HistoManager::BookHisto()                 << 118 void HistoManager::bookHisto()
114 {                                                 119 {
115   fHistoBooked = true;                            120   fHistoBooked = true;
116   fHisto->Add1D("1", "Energy deposition (MeV/m << 121   fHisto->Add1D("1","Energy deposition (MeV/mm/event) in the target",
117                 MeV / mm);                     << 122                fNSlices,0.0,fLength/mm,MeV/mm);
118   fHisto->Add1D("2", "Log10 Energy (MeV) of ga << 123   fHisto->Add1D("2","Log10 Energy (MeV) of gammas",fNBinsE,-5.,5.,1.0);
119   fHisto->Add1D("3", "Log10 Energy (MeV) of el << 124   fHisto->Add1D("3","Log10 Energy (MeV) of electrons",fNBinsE,-5.,5.,1.0);
120   fHisto->Add1D("4", "Log10 Energy (MeV) of po << 125   fHisto->Add1D("4","Log10 Energy (MeV) of positrons",fNBinsE,-5.,5.,1.0);
121   fHisto->Add1D("5", "Log10 Energy (MeV) of pr << 126   fHisto->Add1D("5","Log10 Energy (MeV) of protons",fNBinsE,-5.,5.,1.0);
122   fHisto->Add1D("6", "Log10 Energy (MeV) of ne << 127   fHisto->Add1D("6","Log10 Energy (MeV) of neutrons",fNBinsE,-5.,5.,1.0);
123   fHisto->Add1D("7", "Log10 Energy (MeV) of ch << 128   fHisto->Add1D("7","Log10 Energy (MeV) of charged pions",fNBinsE,-4.,6.,1.0);
124   fHisto->Add1D("8", "Log10 Energy (MeV) of pi << 129   fHisto->Add1D("8","Log10 Energy (MeV) of pi0",fNBinsE,-4.,6.,1.0);
125   fHisto->Add1D("9", "Log10 Energy (MeV) of ch << 130   fHisto->Add1D("9","Log10 Energy (MeV) of charged kaons",fNBinsE,-4.,6.,1.0);
126   fHisto->Add1D("10", "Log10 Energy (MeV) of n << 131   fHisto->Add1D("10","Log10 Energy (MeV) of neutral kaons",fNBinsE,-4.,6.,1.0);
127   fHisto->Add1D("11", "Log10 Energy (MeV) of d << 132   fHisto->Add1D("11","Log10 Energy (MeV) of deuterons and tritons",fNBinsE,-5.,5.,1.0);
128   fHisto->Add1D("12", "Log10 Energy (MeV) of H << 133   fHisto->Add1D("12","Log10 Energy (MeV) of He3 and alpha",fNBinsE,-5.,5.,1.0);
129   fHisto->Add1D("13", "Log10 Energy (MeV) of G << 134   fHisto->Add1D("13","Log10 Energy (MeV) of Generic Ions",fNBinsE,-5.,5.,1.0);
130   fHisto->Add1D("14", "Log10 Energy (MeV) of m << 135   fHisto->Add1D("14","Log10 Energy (MeV) of muons",fNBinsE,-4.,6.,1.0);
131   fHisto->Add1D("15", "log10 Energy (MeV) of s << 136   fHisto->Add1D("15","log10 Energy (MeV) of side-leaked neutrons",fNBinsE,-5.,5.,1.0);
132   fHisto->Add1D("16", "log10 Energy (MeV) of f << 137   fHisto->Add1D("16","log10 Energy (MeV) of forward-leaked neutrons",fNBinsE,-5.,5.,1.0);
133   fHisto->Add1D("17", "log10 Energy (MeV) of b << 138   fHisto->Add1D("17","log10 Energy (MeV) of backward-leaked neutrons",fNBinsE,-5.,5.,1.0);
134   fHisto->Add1D("18", "log10 Energy (MeV) of l << 139   fHisto->Add1D("18","log10 Energy (MeV) of leaking protons",fNBinsE,-4.,6.,1.0);
135   fHisto->Add1D("19", "log10 Energy (MeV) of l << 140   fHisto->Add1D("19","log10 Energy (MeV) of leaking charged pions",fNBinsE,-4.,6.,1.0);
136   fHisto->Add1D("20", "Log10 Energy (MeV) of p << 141   fHisto->Add1D("20","Log10 Energy (MeV) of pi+",fNBinsE,-4.,6.,1.0);
137   fHisto->Add1D("21", "Log10 Energy (MeV) of p << 142   fHisto->Add1D("21","Log10 Energy (MeV) of pi-",fNBinsE,-4.,6.,1.0);
138   fHisto->Add1D("22", "Energy deposition in th << 143   fHisto->Add1D("22","Energy deposition in the target normalized to beam energy",
139                 1.0);                          << 144                110,0.0,1.1,1.0);
140   fHisto->Add1D("23", "EM energy deposition in << 145   fHisto->Add1D("23","EM energy deposition in the target normalized to beam energy",
141                 1.0);                          << 146                110,0.0,1.1,1.0);
142   fHisto->Add1D("24", "Pion energy deposition  << 147   fHisto->Add1D("24","Pion energy deposition in the target normalized to beam energy",
143                 1.1, 1.0);                     << 148                110,0.0,1.1,1.0);
144   fHisto->Add1D("25", "Proton energy depositio << 149   fHisto->Add1D("25","Proton energy deposition in the target normalized to beam energy",
145                 1.1, 1.0);                     << 150                110,0.0,1.1,1.0);
146   fHisto->Add1D("26", "Energy (MeV) of pi+", f << 
147   fHisto->Add1D("27", "Energy (MeV) of pi-", f << 
148   fHisto->Add1D("28", "Energy (MeV) of pi0", f << 
149 }                                                 151 }
150                                                   152 
151 //....oooOO0OOooo........oooOO0OOooo........oo    153 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
152                                                   154 
153 void HistoManager::BeginOfRun()                   155 void HistoManager::BeginOfRun()
154 {                                                 156 {
155   fR2 = fRadius * fRadius;                     << 157   fAbsZ0       = -0.5*fLength;
156   fAbsZ0 = 0.5 * fLength;                      << 158   fNevt       = 0;
157   fNevt = 0;                                   << 159   fNelec      = 0;
158   fNelec = 0;                                  << 160   fNposit     = 0;
159   fNposit = 0;                                 << 161   fNgam       = 0;
160   fNgam = 0;                                   << 162   fNstep      = 0;
161   fNstep = 0;                                  << 
162   fNprot_leak = 0;                                163   fNprot_leak = 0;
163   fNpiofNleak = 0;                                164   fNpiofNleak = 0;
164   fNions = 0;                                  << 165   fNions      = 0;
165   fNdeut = 0;                                  << 166   fNdeut      = 0;
166   fNalpha = 0;                                 << 167   fNalpha     = 0;
167   fNkaons = 0;                                 << 168   fNkaons     = 0;
168   fNmuons = 0;                                 << 169   fNmuons     = 0;
169   fNcpions = 0;                                << 170   fNcpions    = 0;
170   fNpi0 = 0;                                   << 171   fNpi0       = 0;
171   fNneutron = 0;                               << 172   fNneutron   = 0;
172   fNproton = 0;                                << 173   fNproton    = 0;
173   fNaproton = 0;                               << 174   fNaproton   = 0;
174   fNneu_forw = 0;                              << 175   fNneu_forw  = 0;
175   fNneu_leak = 0;                              << 176   fNneu_leak  = 0;
176   fNneu_back = 0;                              << 177   fNneu_back  = 0;
177                                                   178 
178   fEdepSum = 0.0;                              << 179   fEdepSum     = 0.0;
179   fEdepSum2 = 0.0;                             << 180   fEdepSum2    = 0.0;
180                                                   181 
181   if (!fHistoBooked) {                         << 182   if(!fHistoBooked) { bookHisto(); }
182     BookHisto();                               << 
183   }                                            << 
184   fHisto->Book();                                 183   fHisto->Book();
185                                                   184 
186   if (fVerbose > 0) {                          << 185   if(fVerbose > 0) {
187     G4cout << "HistoManager: Histograms are bo << 186     G4cout << "HistoManager: Histograms are booked and run has been started"
                                                   >> 187            <<G4endl;
188   }                                               188   }
189 }                                                 189 }
190                                                   190 
191 //....oooOO0OOooo........oooOO0OOooo........oo    191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
192                                                   192 
193 void HistoManager::EndOfRun()                     193 void HistoManager::EndOfRun()
194 {                                                 194 {
                                                   >> 195 
195   G4cout << "HistoManager: End of run actions     196   G4cout << "HistoManager: End of run actions are started" << G4endl;
196                                                   197 
197   // Average values                               198   // Average values
198   G4cout << "================================= << 199   G4cout<<"========================================================"<<G4endl;
199                                                   200 
200   G4double x = (G4double)fNevt;                   201   G4double x = (G4double)fNevt;
201   if (fNevt > 0) {                             << 202   if(fNevt > 0) { x = 1.0/x; }
202     x = 1.0 / x;                               << 
203   }                                            << 
204                                                   203 
205   G4double xe = x * (G4double)fNelec;          << 204   G4double xe = x*(G4double)fNelec;
206   G4double xg = x * (G4double)fNgam;           << 205   G4double xg = x*(G4double)fNgam;
207   G4double xp = x * (G4double)fNposit;         << 206   G4double xp = x*(G4double)fNposit;
208   G4double xs = x * (G4double)fNstep;          << 207   G4double xs = x*(G4double)fNstep;
209   G4double xn = x * (G4double)fNneutron;       << 208   G4double xn = x*(G4double)fNneutron;
210   G4double xpn = x * (G4double)fNproton;       << 209   G4double xpn = x*(G4double)fNproton;
211   G4double xap = x * (G4double)fNaproton;      << 210   G4double xap = x*(G4double)fNaproton;
212   G4double xnf = x * (G4double)fNneu_forw;     << 211   G4double xnf = x*(G4double)fNneu_forw;
213   G4double xnb = x * (G4double)fNneu_leak;     << 212   G4double xnb = x*(G4double)fNneu_leak;
214   G4double xnbw = x * (G4double)fNneu_back;    << 213   G4double xnbw= x*(G4double)fNneu_back;
215   G4double xpl = x * (G4double)fNprot_leak;    << 214   G4double xpl = x*(G4double)fNprot_leak;
216   G4double xal = x * (G4double)fNpiofNleak;    << 215   G4double xal = x*(G4double)fNpiofNleak;
217   G4double xpc = x * (G4double)fNcpions;       << 216   G4double xpc = x*(G4double)fNcpions;
218   G4double xp0 = x * (G4double)fNpi0;          << 217   G4double xp0 = x*(G4double)fNpi0;
219   G4double xpk = x * (G4double)fNkaons;        << 218   G4double xpk = x*(G4double)fNkaons;
220   G4double xpm = x * (G4double)fNmuons;        << 219   G4double xpm = x*(G4double)fNmuons;
221   G4double xid = x * (G4double)fNdeut;         << 220   G4double xid = x*(G4double)fNdeut;
222   G4double xia = x * (G4double)fNalpha;        << 221   G4double xia = x*(G4double)fNalpha;
223   G4double xio = x * (G4double)fNions;         << 222   G4double xio = x*(G4double)fNions;
224                                                   223 
225   fEdepSum *= x;                               << 224   fEdepSum  *= x;
226   fEdepSum2 *= x;                                 225   fEdepSum2 *= x;
227   fEdepSum2 -= fEdepSum * fEdepSum;            << 226   fEdepSum2 -= fEdepSum*fEdepSum;
228   if (fEdepSum2 > 0.0) {                       << 227   if(fEdepSum2 > 0.0) { fEdepSum2 = std::sqrt(fEdepSum2); }
229     fEdepSum2 = std::sqrt(fEdepSum2);          << 228   else                { fEdepSum2 = 0.0; }
230   }                                            << 229 
231   else {                                       << 230   G4cout                         << "Beam particle                        "
232     fEdepSum2 = 0.0;                           << 231                                  << fPrimaryDef->GetParticleName() <<G4endl;
233   }                                            << 232   G4cout                         << "Beam Energy(MeV)                     " 
234                                                << 233                                  << fPrimaryKineticEnergy/MeV <<G4endl;
235   G4cout << "Beam particle                     << 234   G4cout                         << "Number of events                     " << fNevt <<G4endl;
236   G4cout << "Beam Energy(MeV)                  << 235   G4cout << std::setprecision(4) << "Average energy deposit (MeV)         " << fEdepSum/MeV 
237   G4cout << "Number of events                  << 236          << "   RMS(MeV) " << fEdepSum2/MeV << G4endl;
238   G4cout << std::setprecision(4) << "Average e << 
239          << "   RMS(MeV) " << fEdepSum2 / MeV  << 
240   G4cout << std::setprecision(4) << "Average n    237   G4cout << std::setprecision(4) << "Average number of steps              " << xs << G4endl;
241   G4cout << std::setprecision(4) << "Average n    238   G4cout << std::setprecision(4) << "Average number of gamma              " << xg << G4endl;
242   G4cout << std::setprecision(4) << "Average n    239   G4cout << std::setprecision(4) << "Average number of e-                 " << xe << G4endl;
243   G4cout << std::setprecision(4) << "Average n    240   G4cout << std::setprecision(4) << "Average number of e+                 " << xp << G4endl;
244   G4cout << std::setprecision(4) << "Average n    241   G4cout << std::setprecision(4) << "Average number of neutrons           " << xn << G4endl;
245   G4cout << std::setprecision(4) << "Average n    242   G4cout << std::setprecision(4) << "Average number of protons            " << xpn << G4endl;
246   G4cout << std::setprecision(4) << "Average n    243   G4cout << std::setprecision(4) << "Average number of antiprotons        " << xap << G4endl;
247   G4cout << std::setprecision(4) << "Average n    244   G4cout << std::setprecision(4) << "Average number of pi+ & pi-          " << xpc << G4endl;
248   G4cout << std::setprecision(4) << "Average n    245   G4cout << std::setprecision(4) << "Average number of pi0                " << xp0 << G4endl;
249   G4cout << std::setprecision(4) << "Average n    246   G4cout << std::setprecision(4) << "Average number of kaons              " << xpk << G4endl;
250   G4cout << std::setprecision(4) << "Average n    247   G4cout << std::setprecision(4) << "Average number of muons              " << xpm << G4endl;
251   G4cout << std::setprecision(4) << "Average n    248   G4cout << std::setprecision(4) << "Average number of deuterons+tritons  " << xid << G4endl;
252   G4cout << std::setprecision(4) << "Average n    249   G4cout << std::setprecision(4) << "Average number of He3+alpha          " << xia << G4endl;
253   G4cout << std::setprecision(4) << "Average n    250   G4cout << std::setprecision(4) << "Average number of ions               " << xio << G4endl;
254   G4cout << std::setprecision(4) << "Average n    251   G4cout << std::setprecision(4) << "Average number of forward neutrons   " << xnf << G4endl;
255   G4cout << std::setprecision(4) << "Average n    252   G4cout << std::setprecision(4) << "Average number of reflected neutrons " << xnb << G4endl;
256   G4cout << std::setprecision(4) << "Average n    253   G4cout << std::setprecision(4) << "Average number of leaked neutrons    " << xnbw << G4endl;
257   G4cout << std::setprecision(4) << "Average n    254   G4cout << std::setprecision(4) << "Average number of proton leak        " << xpl << G4endl;
258   G4cout << std::setprecision(4) << "Average n    255   G4cout << std::setprecision(4) << "Average number of pion leak          " << xal << G4endl;
259   G4cout << "================================= << 256   G4cout<<"========================================================"<<G4endl;
260   G4cout << G4endl;                            << 257   G4cout<<G4endl;
261                                                   258 
262   // normalise histograms                         259   // normalise histograms
263   for (G4int i = 0; i < fNHisto; i++) {        << 260   for(G4int i=0; i<fNHisto; i++) { 
264     fHisto->ScaleH1(i, x);                     << 261     fHisto->ScaleH1(i,x);
265   }                                               262   }
266                                                   263 
267   fHisto->Save();                                 264   fHisto->Save();
268 }                                                 265 }
269                                                   266 
270 //....oooOO0OOooo........oooOO0OOooo........oo    267 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
271                                                   268 
272 void HistoManager::BeginOfEvent()                 269 void HistoManager::BeginOfEvent()
273 {                                                 270 {
274   fEdepEvt = 0.0;                                 271   fEdepEvt = 0.0;
275   fEdepEM = 0.0;                               << 272   fEdepEM  = 0.0;
276   fEdepPI = 0.0;                               << 273   fEdepPI  = 0.0;
277   fEdepP = 0.0;                                << 274   fEdepP   = 0.0;
278 }                                                 275 }
279                                                   276 
280 //....oooOO0OOooo........oooOO0OOooo........oo    277 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
281                                                   278 
282 void HistoManager::EndOfEvent()                   279 void HistoManager::EndOfEvent()
283 {                                                 280 {
284   fEdepSum += fEdepEvt;                        << 281   fEdepSum  += fEdepEvt;
285   fEdepSum2 += fEdepEvt * fEdepEvt;            << 282   fEdepSum2 += fEdepEvt*fEdepEvt;
286   fHisto->Fill(21, fEdepEvt / fPrimaryKineticE << 283   fHisto->Fill(21,fEdepEvt/fPrimaryKineticEnergy,1.0);
287   fHisto->Fill(22, fEdepEM / fPrimaryKineticEn << 284   fHisto->Fill(22,fEdepEM/fPrimaryKineticEnergy,1.0);
288   fHisto->Fill(23, fEdepPI / fPrimaryKineticEn << 285   fHisto->Fill(23,fEdepPI/fPrimaryKineticEnergy,1.0);
289   fHisto->Fill(24, fEdepP / fPrimaryKineticEne << 286   fHisto->Fill(24,fEdepP/fPrimaryKineticEnergy,1.0);
290 }                                                 287 }
291                                                   288 
292 //....oooOO0OOooo........oooOO0OOooo........oo    289 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
293                                                   290 
294 void HistoManager::ScoreNewTrack(const G4Track    291 void HistoManager::ScoreNewTrack(const G4Track* track)
295 {                                                 292 {
296   const G4ParticleDefinition* pd = track->GetD    293   const G4ParticleDefinition* pd = track->GetDefinition();
297   const G4String name = pd->GetParticleName(); << 294   G4String name = pd->GetParticleName();
298   const G4double ee = track->GetKineticEnergy( << 295   G4double e = track->GetKineticEnergy();
299   G4double e = ee;                             << 
300                                                   296 
301   // Primary track                                297   // Primary track
302   if (0 == track->GetParentID()) {             << 298   if(0 == track->GetParentID()) {
                                                   >> 299 
303     fNevt++;                                      300     fNevt++;
304     fPrimaryKineticEnergy = e;                    301     fPrimaryKineticEnergy = e;
305     fPrimaryDef = pd;                             302     fPrimaryDef = pd;
306     G4ThreeVector dir = track->GetMomentumDire    303     G4ThreeVector dir = track->GetMomentumDirection();
307     if (1 < fVerbose)                          << 304     if(1 < fVerbose) 
308       G4cout << "### Primary " << name << " ki << 305       G4cout << "### Primary " << name 
309              << "; m(MeV)= " << pd->GetPDGMass << 306              << " kinE(MeV)= " << e/MeV
310              << ";  dir= " << track->GetMoment << 307              << "; m(MeV)= " << pd->GetPDGMass()/MeV
                                                   >> 308              << "; pos(mm)= " << track->GetPosition()/mm 
                                                   >> 309              << ";  dir= " << track->GetMomentumDirection() 
                                                   >> 310              << G4endl;
311                                                   311 
312     // Secondary track                            312     // Secondary track
313   }                                            << 313   } else {
314   else {                                       << 314     if(1 < fVerbose) 
315     if (1 < fVerbose)                          << 315       G4cout << "=== Secondary " << name 
316       G4cout << "=== Secondary " << name << "  << 316              << " kinE(MeV)= " << e/MeV
317              << "; m(MeV)= " << pd->GetPDGMass << 317              << "; m(MeV)= " << pd->GetPDGMass()/MeV
318              << ";  dir= " << track->GetMoment << 318              << "; pos(mm)= " << track->GetPosition()/mm 
319     e = (e > 0.0) ? std::log10(e / MeV) : -100 << 319              << ";  dir= " << track->GetMomentumDirection() 
320     if (pd == G4Gamma::Gamma()) {              << 320              << G4endl;
                                                   >> 321     e = std::log10(e/MeV);
                                                   >> 322     if(pd == G4Gamma::Gamma()) {
321       fNgam++;                                    323       fNgam++;
322       fHisto->Fill(1, e, 1.0);                 << 324       fHisto->Fill(1,e,1.0);
323     }                                          << 325     } else if ( pd == G4Electron::Electron()) {
324     else if (pd == G4Electron::Electron()) {   << 
325       fNelec++;                                   326       fNelec++;
326       fHisto->Fill(2, e, 1.0);                 << 327       fHisto->Fill(2,e,1.0);
327     }                                          << 328     } else if ( pd == G4Positron::Positron()) {
328     else if (pd == G4Positron::Positron()) {   << 
329       fNposit++;                                  329       fNposit++;
330       fHisto->Fill(3, e, 1.0);                 << 330       fHisto->Fill(3,e,1.0);
331     }                                          << 331     } else if ( pd == G4Proton::Proton()) {
332     else if (pd == G4Proton::Proton()) {       << 
333       fNproton++;                                 332       fNproton++;
334       fHisto->Fill(4, e, 1.0);                 << 333       fHisto->Fill(4,e,1.0);
335     }                                          << 334     } else if ( pd == fNeutron) {
336     else if (pd == fNeutron) {                 << 
337       fNneutron++;                                335       fNneutron++;
338       fHisto->Fill(5, e, 1.0);                 << 336       fHisto->Fill(5,e,1.0);
339     }                                          << 337     } else if ( pd == G4AntiProton::AntiProton()) {
340     else if (pd == G4AntiProton::AntiProton()) << 
341       fNaproton++;                                338       fNaproton++;
342     }                                          << 339     } else if ( pd == G4PionPlus::PionPlus() ) {
343     else if (pd == G4PionPlus::PionPlus()) {   << 
344       fNcpions++;                                 340       fNcpions++;
345       fHisto->Fill(6, e, 1.0);                 << 341       fHisto->Fill(6,e,1.0);
346       fHisto->Fill(19, e, 1.0);                << 342       fHisto->Fill(19,e,1.0);
347       fHisto->Fill(25, ee, 1.0);               << 343 
348     }                                          << 344     } else if ( pd == G4PionMinus::PionMinus()) {
349     else if (pd == G4PionMinus::PionMinus()) { << 
350       fNcpions++;                                 345       fNcpions++;
351       fHisto->Fill(6, e, 1.0);                 << 346       fHisto->Fill(6,e,1.0);
352       fHisto->Fill(20, e, 1.0);                << 347       fHisto->Fill(20,e,1.0);
353       fHisto->Fill(26, ee, 1.0);               << 348 
354     }                                          << 349     } else if ( pd == G4PionZero::PionZero()) {
355     else if (pd == G4PionZero::PionZero()) {   << 
356       fNpi0++;                                    350       fNpi0++;
357       fHisto->Fill(7, e, 1.0);                 << 351       fHisto->Fill(7,e,1.0);
358       fHisto->Fill(27, ee, 1.0);               << 352     } else if ( pd == G4KaonPlus::KaonPlus() || pd == G4KaonMinus::KaonMinus()) {
359     }                                          << 
360     else if (pd == G4KaonPlus::KaonPlus() || p << 
361       fNkaons++;                                  353       fNkaons++;
362       fHisto->Fill(8, e, 1.0);                 << 354       fHisto->Fill(8,e,1.0);
363     }                                          << 355     } else if ( pd == G4KaonZeroShort::KaonZeroShort() || pd == G4KaonZeroLong::KaonZeroLong()) {
364     else if (pd == G4KaonZeroShort::KaonZeroSh << 
365       fNkaons++;                                  356       fNkaons++;
366       fHisto->Fill(9, e, 1.0);                 << 357       fHisto->Fill(9,e,1.0);
367     }                                          << 358     } else if ( pd == G4Deuteron::Deuteron() || pd == G4Triton::Triton()) {
368     else if (pd == G4Deuteron::Deuteron() || p << 
369       fNdeut++;                                   359       fNdeut++;
370       fHisto->Fill(10, e, 1.0);                << 360       fHisto->Fill(10,e,1.0);
371     }                                          << 361     } else if ( pd == G4He3::He3() || pd == G4Alpha::Alpha()) {
372     else if (pd == G4He3::He3() || pd == G4Alp << 
373       fNalpha++;                                  362       fNalpha++;
374       fHisto->Fill(11, e, 1.0);                << 363       fHisto->Fill(11,e,1.0);
375     }                                          << 364     } else if ( pd->GetParticleType() == "nucleus") {
376     else if (pd->GetParticleType() == "nucleus << 
377       fNions++;                                   365       fNions++;
378       fHisto->Fill(12, e, 1.0);                << 366       fHisto->Fill(12,e,1.0);
379     }                                          << 367     } else if ( pd == G4MuonPlus::MuonPlus() || pd == G4MuonMinus::MuonMinus()) {
380     else if (pd == G4MuonPlus::MuonPlus() || p << 
381       fNmuons++;                                  368       fNmuons++;
382       fHisto->Fill(13, e, 1.0);                << 369       fHisto->Fill(13,e,1.0);    
383     }                                             370     }
384   }                                               371   }
385 }                                                 372 }
386                                                   373 
387 //....oooOO0OOooo........oooOO0OOooo........oo    374 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
388                                                   375 
389 void HistoManager::AddTargetStep(const G4Step*    376 void HistoManager::AddTargetStep(const G4Step* step)
390 {                                                 377 {
391   ++fNstep;                                    << 378   fNstep++;
392   const G4double fEdep = step->GetTotalEnergyD << 379   G4double fEdep = step->GetTotalEnergyDeposit();
393   const G4Track* track = step->GetTrack();     << 380   if(fEdep >= DBL_MIN) { 
394   const G4ParticleDefinition* pd = track->GetD << 381     const G4Track* track = step->GetTrack();
395   if (1 < fVerbose) {                          << 382 
396     G4cout << "TargetSD::ProcessHits: " << pd- << 383     G4ThreeVector pos = 
397            << " E(MeV)=" << track->GetKineticE << 384       (step->GetPreStepPoint()->GetPosition() +
398            << " beta1= " << step->GetPreStepPo << 385        step->GetPostStepPoint()->GetPosition())*0.5;
399            << " beta2= " << step->GetPostStepP << 
400            << " weight= " << step->GetTrack()- << 
401            << " t(ns)=" << track->GetGlobalTim << 
402   }                                            << 
403   if (fEdep > 0.0) {                           << 
404     G4ThreeVector pos =                        << 
405       (step->GetPreStepPoint()->GetPosition()  << 
406                                                   386 
407     G4double z = pos.z() + fAbsZ0;             << 387     G4double z = pos.z() - fAbsZ0;
408                                                   388 
409     // scoring                                    389     // scoring
410     fEdepEvt += fEdep;                            390     fEdepEvt += fEdep;
411     fHisto->Fill(0, z, fEdep);                 << 391     fHisto->Fill(0,z,fEdep);
                                                   >> 392     const G4ParticleDefinition* pd = track->GetDefinition();
412                                                   393 
413     if (pd == G4Gamma::Gamma() || pd == G4Elec << 394     if(pd == G4Gamma::Gamma() || pd == G4Electron::Electron() 
                                                   >> 395        || pd == G4Positron::Positron()) {
414       fEdepEM += fEdep;                           396       fEdepEM += fEdep;
415     }                                          << 397     } else if ( pd == G4PionPlus::PionPlus() || pd == G4PionMinus::PionMinus()) {
416     else if (pd == G4PionPlus::PionPlus() || p << 
417       fEdepPI += fEdep;                           398       fEdepPI += fEdep;
418     }                                          << 399     } else if ( pd == G4Proton::Proton() || pd == G4AntiProton::AntiProton()) {
419     else if (pd == G4Proton::Proton() || pd == << 400       fEdepP  += fEdep;
420       fEdepP += fEdep;                         << 
421     }                                             401     }
422                                                   402 
423     if (1 < fVerbose) {                        << 403     if(1 < fVerbose) {
424       G4cout << "HistoManager::AddEnergy: e(ke << 404       G4cout << "HistoManager::AddEnergy: e(keV)= " << fEdep/keV
425              << "; step(mm)= " << step->GetSte << 405              << "; z(mm)= " << z/mm
426              << " E(MeV)= " << track->GetKinet << 406              << "; step(mm)= " << step->GetStepLength()/mm
                                                   >> 407              << " by " << pd->GetParticleName()
                                                   >> 408              << " E(MeV)= " << track->GetKineticEnergy()/MeV
                                                   >> 409              << G4endl;
427     }                                             410     }
428   }                                               411   }
429 }                                                 412 }
430                                                   413 
431 //....oooOO0OOooo........oooOO0OOooo........oo    414 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
432                                                   415 
433 void HistoManager::AddLeakingParticle(const G4    416 void HistoManager::AddLeakingParticle(const G4Track* track)
434 {                                                 417 {
435   const G4ParticleDefinition* pd = track->GetD << 418   const G4ParticleDefinition* pd = track->GetDefinition(); 
436   const G4StepPoint* sp = track->GetStep()->Ge << 419   G4double e = std::log10(track->GetKineticEnergy()/MeV);
437   const G4double ekin = sp->GetKineticEnergy() << 
438   G4double e = -100.;                          << 
439   if (ekin > 0.0) {                            << 
440     e = std::log10(ekin / CLHEP::MeV);         << 
441   }                                            << 
442   else {                                       << 
443     G4cout << "### HistoManager::AddLeakingPar << 
444            << "  Eprestep(MeV)=" << ekin << "  << 
445     return;                                    << 
446   }                                            << 
447                                                   420 
448   const G4ThreeVector& pos = sp->GetPosition() << 421   G4ThreeVector pos = track->GetPosition();
449   const G4ThreeVector& dir = sp->GetMomentumDi << 422   G4ThreeVector dir = track->GetMomentumDirection();
450   G4double x = pos.x();                           423   G4double x = pos.x();
451   G4double y = pos.y();                           424   G4double y = pos.y();
452   G4double z = pos.z();                           425   G4double z = pos.z();
453   G4double vx = dir.x();                       << 426  
454   G4double vy = dir.y();                       << 
455   G4double vz = dir.z();                       << 
456                                                << 
457   G4bool isLeaking = false;                       427   G4bool isLeaking = false;
458                                                   428 
459   // Forward                                   << 429   // Forward 
460   const G4double del = 0.001 * CLHEP::mm;      << 430   if(z > -fAbsZ0 && dir.z() > 0.0) {
461   if (std::abs(z - fAbsZ0) < del && vz > 0.0)  << 
462     isLeaking = true;                             431     isLeaking = true;
463     if (pd == fNeutron) {                      << 432     if(pd == fNeutron) {
464       ++fNneu_forw;                               433       ++fNneu_forw;
465       fHisto->Fill(15, e, 1.0);                << 434       fHisto->Fill(15,e,1.0);
466     }                                          << 435     } else isLeaking = true;
467     else                                       << 
468       isLeaking = true;                        << 
469                                                   436 
470     // Backward                                   437     // Backward
471   }                                            << 438   } else if (z < fAbsZ0 && dir.z() < 0.0) {
472   else if (std::abs(z + fAbsZ0) < del && vz <  << 
473     isLeaking = true;                             439     isLeaking = true;
474     if (pd == fNeutron) {                      << 440     if(pd == fNeutron) {
475       ++fNneu_back;                               441       ++fNneu_back;
476       fHisto->Fill(16, e, 1.0);                << 442       fHisto->Fill(16,e,1.0);
477     }                                          << 443     } else isLeaking = true;
478     else                                       << 
479       isLeaking = true;                        << 
480                                                   444 
481     // Side                                       445     // Side
482   }                                            << 446   } else if (std::abs(z) <= -fAbsZ0 && x*dir.x() + y*dir.y() > 0.0) {
483   else if (std::abs(z) <= fAbsZ0 + del && x *  << 
484            && std::abs(x * x + y * y - fR2) <  << 
485   {                                            << 
486     isLeaking = true;                             447     isLeaking = true;
487     if (pd == fNeutron) {                      << 448     if(pd == fNeutron) {
488       ++fNneu_leak;                               449       ++fNneu_leak;
489       fHisto->Fill(14, e, 1.0);                << 450       fHisto->Fill(14,e,1.0);
490     }                                          << 451     } else isLeaking = true;
491     else                                       << 
492       isLeaking = true;                        << 
493   }                                               452   }
494                                                   453 
495   // protons and pions                            454   // protons and pions
496   if (isLeaking) {                             << 455   if(isLeaking) {
497     if (pd == G4Proton::Proton()) {            << 456     if(pd == G4Proton::Proton()) {
498       fHisto->Fill(17, e, 1.0);                << 457       fHisto->Fill(17,e,1.0);
499       ++fNprot_leak;                              458       ++fNprot_leak;
500     }                                          << 459     } else if (pd == G4PionPlus::PionPlus() || pd == G4PionMinus::PionMinus()) {
501     else if (pd == G4PionPlus::PionPlus() || p << 460       fHisto->Fill(18,e,1.0);
502       fHisto->Fill(18, e, 1.0);                << 
503       ++fNpiofNleak;                              461       ++fNpiofNleak;
504     }                                             462     }
505   }                                               463   }
506 }                                                 464 }
507                                                   465 
508 //....oooOO0OOooo........oooOO0OOooo........oo    466 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
509                                                   467 
510 void HistoManager::SetVerbose(G4int val)       << 468 void HistoManager::SetVerbose(G4int val)        
511 {                                                 469 {
512   fVerbose = val;                              << 470   fVerbose = val; 
513   fHisto->SetVerbose(val);                        471   fHisto->SetVerbose(val);
514 }                                                 472 }
515                                                   473 
516 //....oooOO0OOooo........oooOO0OOooo........oo    474 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
517                                                   475 
518 void HistoManager::Fill(G4int id, G4double x,     476 void HistoManager::Fill(G4int id, G4double x, G4double w)
519 {                                                 477 {
520   fHisto->Fill(id, x, w);                         478   fHisto->Fill(id, x, w);
521 }                                                 479 }
522                                                   480 
523 //....oooOO0OOooo........oooOO0OOooo........oo    481 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 482 
524                                                   483