Geant4 Cross Reference |
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 /// \file exoticphysics/dmparticle/src/Run.cc 27 /// \brief Implementation of the Run class 28 // 29 30 #include "Run.hh" 31 32 #include "TestParameters.hh" 33 34 #include "G4ElectronIonPair.hh" 35 #include "G4LossTableManager.hh" 36 #include "G4PhysicalConstants.hh" 37 #include "G4Run.hh" 38 #include "G4Step.hh" 39 #include "G4SystemOfUnits.hh" 40 #include "Randomize.hh" 41 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 43 44 Run::Run() : G4Run(), fParam(TestParameters::GetPointer()) 45 { 46 fTotStepGas = fOverflow = fTotEdep = fStepGas = fMaxEnergy = 0.0; 47 fEvt = fNbins = 0; 48 } 49 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 51 52 Run::~Run() {} 53 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 55 56 void Run::BeginOfRun() 57 { 58 // initilise scoring 59 fTotStepGas = fOverflow = fTotEdep = fStepGas = fMaxEnergy = 0.0; 60 fEvt = 0; 61 62 SetVerbose(1); 63 64 fNbins = fParam->GetNumberBins(); 65 fMaxEnergy = fParam->GetMaxEnergy(); 66 67 fEgas.resize(fNbins, 0.0); 68 fEdep.reset(); 69 70 if (fVerbose > 0) { 71 G4cout << " BinsE= " << fNbins << " Emax(keV)= " << fMaxEnergy / keV << G4endl; 72 } 73 } 74 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 76 77 void Run::EndOfRun() 78 { 79 G4int nEvt = GetNumberOfEvent(); 80 81 G4double norm = (nEvt > 0) ? 1.0 / (G4double)nEvt : 0.0; 82 83 fTotStepGas *= norm; 84 fOverflow *= norm; 85 86 G4double y1 = fEdep.mean(); 87 G4double y2 = fEdep.rms(); 88 89 // G4double de = fMaxEnergy/G4double(fNbins); 90 // G4double x1 = -de*0.5; 91 92 G4cout << " ====================================================" << G4endl; 93 G4cout << " Beam Particle: " << fParam->GetBeamParticle()->GetParticleName() << G4endl 94 << " Ekin(GeV) = " << fParam->GetBeamEnergy() / GeV << G4endl 95 << " Z(mm) = " << fParam->GetPositionZ() / mm << G4endl; 96 G4cout << " ================== run summary =====================" << G4endl; 97 G4int prec = G4cout.precision(5); 98 G4cout << " End of Run TotNbofEvents = " << nEvt << G4endl; 99 100 G4cout << G4endl; 101 G4cout << " Mean energy deposit in absorber = " << y1 / GeV << " +- " 102 << y2 * std::sqrt(norm) / GeV << " GeV; "; 103 if (y1 > 0.0) { 104 G4cout << " RMS/Emean = " << y2 / y1; 105 } 106 G4cout << G4endl; 107 G4cout << " Mean number of steps in absorber= " << fTotStepGas << G4endl; 108 G4cout << G4endl; 109 /* 110 G4cout << " ====== Energy deposit distribution Noverflows= " << fOverflow 111 << " ====== " << G4endl ; 112 G4cout << " bin nb Elow entries normalized " << G4endl; 113 114 x1 = 0.0; 115 116 for(G4int j=0; j<fNbins; ++j) 117 { 118 G4cout << std::setw(5) << j << std::setw(10) << x1/keV 119 << std::setw(12) << fEgas[j] << std::setw(12) << fEgas[j]*norm 120 << G4endl ; 121 x1 += de; 122 } 123 */ 124 G4cout.precision(prec); 125 126 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); 127 128 // normalize histograms 129 130 analysisManager->ScaleH1(1, norm); 131 132 G4cout << " ================== run end ==========================" << G4endl; 133 } 134 135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 136 137 void Run::BeginOfEvent() 138 { 139 fTotEdep = 0.0; 140 fStepGas = 0; 141 ++fEvt; 142 } 143 144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 145 146 void Run::EndOfEvent() 147 { 148 /* 149 G4cout << "fStepGas = " << fStepGas << " fTotStepGas= " << fTotStepGas 150 << " fTotEdep= " << fTotEdep << " fNbins= " << fNbins 151 << " fMaxEnergy= " << fMaxEnergy << G4endl; 152 */ 153 fTotStepGas += fStepGas; 154 155 G4int idx = G4lrint(fTotEdep * fNbins / fMaxEnergy); 156 157 if (idx < 0) { 158 fEgas[0] += 1.0; 159 } 160 if (idx >= fNbins) { 161 fOverflow += 1.0; 162 } 163 else { 164 fEgas[idx] += 1.0; 165 } 166 167 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); 168 169 // fill histo 170 171 analysisManager->FillH1(1, fTotEdep / GeV, 1.0); 172 fEdep.fill(fTotEdep, 1.0); 173 } 174 175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 176 // 177 // Allows one to accumulate data from different threads 178 179 void Run::Merge(const G4Run* run) 180 { 181 const Run* localRun = static_cast<const Run*>(run); 182 183 fTotStepGas += localRun->fTotStepGas; 184 fOverflow += localRun->fOverflow; 185 186 G4StatDouble* stat = const_cast<G4StatDouble*>(localRun->GetStat()); 187 188 fEdep.add(stat); 189 190 for (G4int j = 0; j < fNbins; ++j) { 191 fEgas[j] += localRun->fEgas[j]; 192 } 193 194 G4Run::Merge(run); 195 } 196 197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 198 // 199 // Is called from TargetSD 200 201 void Run::AddEnergy(G4double edep, const G4Step* step) 202 { 203 if (1 < fVerbose) { 204 G4cout << "Run::AddEnergy: e(keV)= " << edep / keV << G4endl; 205 } 206 fTotEdep += edep; 207 208 if (step) { 209 if (1 == step->GetTrack()->GetTrackID()) { 210 fStepGas += 1.0; 211 } 212 } 213 } 214 215 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 216