Geant4 Cross Reference |
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 electromagnetic/TestEm2/src/Run.cc 26 /// \file electromagnetic/TestEm2/src/Run.cc 27 /// \brief Implementation of the Run class 27 /// \brief Implementation of the Run class 28 // 28 // 29 // << 29 // $Id: Run.cc 75577 2013-11-04 12:03:26Z vnivanch $ >> 30 // 30 //....oooOO0OOooo........oooOO0OOooo........oo 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 31 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 32 33 33 #include "Run.hh" 34 #include "Run.hh" 34 35 35 #include "EmAcceptance.hh" << 36 #include "DetectorConstruction.hh" 36 #include "PrimaryGeneratorAction.hh" 37 #include "PrimaryGeneratorAction.hh" >> 38 #include "EmAcceptance.hh" 37 39 38 #include "G4Run.hh" 40 #include "G4Run.hh" 39 #include "G4SystemOfUnits.hh" << 40 #include "G4UnitsTable.hh" 41 #include "G4UnitsTable.hh" 41 << 42 #include "G4SystemOfUnits.hh" 42 #include <iomanip> 43 #include <iomanip> 43 44 44 //....oooOO0OOooo........oooOO0OOooo........oo 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 45 46 46 Run::Run(DetectorConstruction* det, PrimaryGen << 47 Run::Run(DetectorConstruction* det, PrimaryGeneratorAction* kin) >> 48 :G4Run(),fDet(det),fKin(kin), >> 49 f_nLbin(MaxBin),f_nRbin(MaxBin) 47 { 50 { 48 Reset(); 51 Reset(); 49 } 52 } 50 53 51 //....oooOO0OOooo........oooOO0OOooo........oo 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 52 55 53 void Run::Reset() 56 void Run::Reset() 54 { << 57 { 55 f_nLbin = fDet->GetnLtot(); 58 f_nLbin = fDet->GetnLtot(); 56 f_dEdL.resize(f_nLbin); 59 f_dEdL.resize(f_nLbin); 57 fSumELongit.resize(f_nLbin); 60 fSumELongit.resize(f_nLbin); 58 fSumELongitCumul.resize(f_nLbin); 61 fSumELongitCumul.resize(f_nLbin); 59 fSumE2Longit.resize(f_nLbin); 62 fSumE2Longit.resize(f_nLbin); 60 fSumE2LongitCumul.resize(f_nLbin); 63 fSumE2LongitCumul.resize(f_nLbin); 61 64 62 f_nRbin = fDet->GetnRtot(); 65 f_nRbin = fDet->GetnRtot(); 63 f_dEdR.resize(f_nRbin); 66 f_dEdR.resize(f_nRbin); 64 fSumERadial.resize(f_nRbin); 67 fSumERadial.resize(f_nRbin); 65 fSumERadialCumul.resize(f_nRbin); 68 fSumERadialCumul.resize(f_nRbin); 66 fSumE2Radial.resize(f_nRbin); 69 fSumE2Radial.resize(f_nRbin); 67 fSumE2RadialCumul.resize(f_nRbin); 70 fSumE2RadialCumul.resize(f_nRbin); 68 << 71 69 fChargedStep = 0.0; 72 fChargedStep = 0.0; 70 fNeutralStep = 0.0; << 73 fNeutralStep = 0.0; 71 << 74 72 fVerbose = 0; 75 fVerbose = 0; 73 << 76 74 // initialize arrays of cumulative energy de << 77 //initialize arrays of cumulative energy deposition 75 // 78 // 76 for (G4int i = 0; i < f_nLbin; ++i) { << 79 for (G4int i=0; i<f_nLbin; i++) 77 fSumELongit[i] = fSumE2Longit[i] = fSumELo << 80 fSumELongit[i]=fSumE2Longit[i]=fSumELongitCumul[i]=fSumE2LongitCumul[i]=0.; 78 } << 81 79 for (G4int j = 0; j < f_nRbin; ++j) { << 82 for (G4int j=0; j<f_nRbin; j++) 80 fSumERadial[j] = fSumE2Radial[j] = fSumERa << 83 fSumERadial[j]=fSumE2Radial[j]=fSumERadialCumul[j]=fSumE2RadialCumul[j]=0.; 81 } << 84 82 // initialize track length << 85 //initialize track length 83 fSumChargTrLength = fSum2ChargTrLength = fSu << 86 fSumChargTrLength=fSum2ChargTrLength=fSumNeutrTrLength=fSum2NeutrTrLength=0.; >> 87 84 } 88 } 85 89 86 //....oooOO0OOooo........oooOO0OOooo........oo 90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 87 91 88 void Run::InitializePerEvent() << 92 Run::~Run() 89 { << 93 {} 90 // initialize arrays of energy deposit per b << 91 for (G4int i = 0; i < f_nLbin; ++i) { << 92 f_dEdL[i] = 0.; << 93 } << 94 94 95 for (G4int j = 0; j < f_nRbin; ++j) { << 95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 96 f_dEdR[j] = 0.; << 97 } << 98 96 99 // initialize tracklength << 97 void Run::InitializePerEvent() >> 98 { >> 99 //initialize arrays of energy deposit per bin >> 100 for (G4int i=0; i<f_nLbin; i++) >> 101 { f_dEdL[i] = 0.; } >> 102 >> 103 for (G4int j=0; j<f_nRbin; j++) >> 104 { f_dEdR[j] = 0.; } >> 105 >> 106 //initialize tracklength 100 fChargTrLength = fNeutrTrLength = 0.; 107 fChargTrLength = fNeutrTrLength = 0.; 101 } 108 } 102 109 103 //....oooOO0OOooo........oooOO0OOooo........oo 110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 104 111 105 void Run::FillPerEvent() 112 void Run::FillPerEvent() 106 { 113 { 107 // accumulate statistic << 114 //accumulate statistic 108 // 115 // 109 G4double dLCumul = 0.; 116 G4double dLCumul = 0.; 110 for (G4int i = 0; i < f_nLbin; ++i) { << 117 for (G4int i=0; i<f_nLbin; i++) 111 fSumELongit[i] += f_dEdL[i]; << 118 { 112 fSumE2Longit[i] += f_dEdL[i] * f_dEdL[i]; << 119 fSumELongit[i] += f_dEdL[i]; 113 dLCumul += f_dEdL[i]; << 120 fSumE2Longit[i] += f_dEdL[i]*f_dEdL[i]; 114 fSumELongitCumul[i] += dLCumul; << 121 dLCumul += f_dEdL[i]; 115 fSumE2LongitCumul[i] += dLCumul * dLCumul; << 122 fSumELongitCumul[i] += dLCumul; 116 } << 123 fSumE2LongitCumul[i] += dLCumul*dLCumul; >> 124 } 117 125 118 G4double dRCumul = 0.; 126 G4double dRCumul = 0.; 119 for (G4int j = 0; j < f_nRbin; j++) { << 127 for (G4int j=0; j<f_nRbin; j++) 120 fSumERadial[j] += f_dEdR[j]; << 128 { 121 fSumE2Radial[j] += f_dEdR[j] * f_dEdR[j]; << 129 fSumERadial[j] += f_dEdR[j]; 122 dRCumul += f_dEdR[j]; << 130 fSumE2Radial[j] += f_dEdR[j]*f_dEdR[j]; 123 fSumERadialCumul[j] += dRCumul; << 131 dRCumul += f_dEdR[j]; 124 fSumE2RadialCumul[j] += dRCumul * dRCumul; << 132 fSumERadialCumul[j] += dRCumul; 125 } << 133 fSumE2RadialCumul[j] += dRCumul*dRCumul; >> 134 } 126 135 127 fSumChargTrLength += fChargTrLength; << 136 fSumChargTrLength += fChargTrLength; 128 fSum2ChargTrLength += fChargTrLength * fChar << 137 fSum2ChargTrLength += fChargTrLength*fChargTrLength; 129 fSumNeutrTrLength += fNeutrTrLength; << 138 fSumNeutrTrLength += fNeutrTrLength; 130 fSum2NeutrTrLength += fNeutrTrLength * fNeut << 139 fSum2NeutrTrLength += fNeutrTrLength*fNeutrTrLength; 131 140 132 // fill histograms << 141 //fill histograms 133 // 142 // 134 143 135 G4double Ekin = fKin->GetParticleGun()->GetP << 144 G4double Ekin=fKin->GetParticleGun()->GetParticleEnergy(); 136 G4double mass = fKin->GetParticleGun()->GetP << 145 G4double mass=fKin->GetParticleGun()->GetParticleDefinition()->GetPDGMass(); 137 G4double radl = fDet->GetMaterial()->GetRadl << 146 G4double radl=fDet->GetMaterial()->GetRadlen(); 138 147 139 G4AnalysisManager* analysisManager = G4Analy 148 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); 140 analysisManager->FillH1(1, 100. * dLCumul / << 149 analysisManager->FillH1(1, 100.*dLCumul/(Ekin+mass)); 141 analysisManager->FillH1(2, fChargTrLength / << 150 analysisManager->FillH1(2, fChargTrLength/radl); 142 analysisManager->FillH1(3, fNeutrTrLength / << 151 analysisManager->FillH1(3, fNeutrTrLength/radl); 143 << 152 144 // profiles << 153 //profiles 145 G4double norm = 100. / (Ekin + mass); << 154 G4double norm = 100./(Ekin+mass); 146 G4double dLradl = fDet->GetdLradl(); << 155 G4double dLradl = fDet->GetdLradl(); 147 for (G4int i = 0; i < f_nLbin; ++i) { << 156 for (G4int i=0; i<f_nLbin; i++) { 148 G4double bin = (i + 0.5) * dLradl; << 157 G4double bin = (i+0.5)*dLradl; 149 analysisManager->FillP1(0, bin, norm * f_d << 158 analysisManager->FillP1(0, bin, norm*f_dEdL[i]/dLradl); 150 } << 159 } 151 G4double dRradl = fDet->GetdRradl(); << 160 G4double dRradl = fDet->GetdRradl(); 152 for (G4int j = 0; j < f_nRbin; ++j) { << 161 for (G4int j=0; j<f_nRbin; j++) { 153 G4double bin = (j + 0.5) * dRradl; << 162 G4double bin = (j+0.5)*dRradl; 154 analysisManager->FillP1(1, bin, norm * f_d << 163 analysisManager->FillP1(1, bin, norm*f_dEdR[j]/dRradl); 155 } << 164 } 156 } 165 } 157 166 158 //....oooOO0OOooo........oooOO0OOooo........oo 167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 159 168 160 void Run::Merge(const G4Run* run) 169 void Run::Merge(const G4Run* run) 161 { 170 { 162 const Run* localRun = static_cast<const Run* 171 const Run* localRun = static_cast<const Run*>(run); 163 172 164 fChargedStep += localRun->fChargedStep; 173 fChargedStep += localRun->fChargedStep; 165 fNeutralStep += localRun->fNeutralStep; << 174 fNeutralStep += localRun->fNeutralStep; 166 175 167 for (G4int i = 0; i < f_nLbin; ++i) { << 176 for (G4int i=0; i<f_nLbin; i++) { 168 fSumELongit[i] += localRun->fSumELongit[i] 177 fSumELongit[i] += localRun->fSumELongit[i]; 169 fSumE2Longit[i] += localRun->fSumE2Longit[ 178 fSumE2Longit[i] += localRun->fSumE2Longit[i]; 170 fSumELongitCumul[i] += localRun->fSumELong 179 fSumELongitCumul[i] += localRun->fSumELongitCumul[i]; 171 fSumE2LongitCumul[i] += localRun->fSumE2Lo 180 fSumE2LongitCumul[i] += localRun->fSumE2LongitCumul[i]; 172 } 181 } 173 182 174 for (G4int j = 0; j < f_nRbin; ++j) { << 183 for (G4int j=0; j<f_nRbin; j++) { 175 fSumERadial[j] += localRun->fSumERadial[j] 184 fSumERadial[j] += localRun->fSumERadial[j]; 176 fSumE2Radial[j] += localRun->fSumE2Radial[ 185 fSumE2Radial[j] += localRun->fSumE2Radial[j]; 177 fSumERadialCumul[j] += localRun->fSumERadi 186 fSumERadialCumul[j] += localRun->fSumERadialCumul[j]; 178 fSumE2RadialCumul[j] += localRun->fSumE2Ra 187 fSumE2RadialCumul[j] += localRun->fSumE2RadialCumul[j]; 179 } 188 } 180 189 181 fSumChargTrLength += localRun->fSumChargTrLe << 190 fSumChargTrLength += localRun->fSumChargTrLength; 182 fSum2ChargTrLength += localRun->fSum2ChargTr 191 fSum2ChargTrLength += localRun->fSum2ChargTrLength; 183 fSumNeutrTrLength += localRun->fSumNeutrTrLe << 192 fSumNeutrTrLength += localRun->fSumNeutrTrLength; 184 fSum2NeutrTrLength += localRun->fSum2NeutrTr 193 fSum2NeutrTrLength += localRun->fSum2NeutrTrLength; 185 194 186 G4Run::Merge(run); << 195 G4Run::Merge(run); 187 } 196 } 188 197 189 //....oooOO0OOooo........oooOO0OOooo........oo 198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 190 199 191 void Run::EndOfRun(G4double edep, G4double rms << 200 void Run::EndOfRun(G4double edep, G4double rms, G4double& limit) 192 { 201 { 193 G4int NbOfEvents = GetNumberOfEvent(); 202 G4int NbOfEvents = GetNumberOfEvent(); 194 203 195 G4double kinEnergy = fKin->GetParticleGun()- 204 G4double kinEnergy = fKin->GetParticleGun()->GetParticleEnergy(); 196 assert(NbOfEvents * kinEnergy > 0); << 205 assert(NbOfEvents*kinEnergy > 0); 197 206 198 fChargedStep /= G4double(NbOfEvents); 207 fChargedStep /= G4double(NbOfEvents); 199 fNeutralStep /= G4double(NbOfEvents); << 208 fNeutralStep /= G4double(NbOfEvents); 200 << 201 G4double mass = fKin->GetParticleGun()->GetP << 202 G4double norme = 100. / (NbOfEvents * (kinEn << 203 209 204 // longitudinal << 210 G4double mass=fKin->GetParticleGun()->GetParticleDefinition()->GetPDGMass(); >> 211 G4double norme = 100./(NbOfEvents*(kinEnergy+mass)); >> 212 >> 213 //longitudinal 205 // 214 // 206 G4double dLradl = fDet->GetdLradl(); 215 G4double dLradl = fDet->GetdLradl(); 207 216 208 MyVector MeanELongit(f_nLbin), rmsELongit(f_ << 217 MyVector MeanELongit(f_nLbin), rmsELongit(f_nLbin); 209 MyVector MeanELongitCumul(f_nLbin), rmsELong 218 MyVector MeanELongitCumul(f_nLbin), rmsELongitCumul(f_nLbin); 210 219 211 G4AnalysisManager* analysisManager = G4Analy 220 G4AnalysisManager* analysisManager = G4AnalysisManager::Instance(); 212 221 213 G4int i; 222 G4int i; 214 for (i = 0; i < f_nLbin; ++i) { << 223 for (i=0; i<f_nLbin; i++) { 215 MeanELongit[i] = norme * fSumELongit[i]; << 224 MeanELongit[i] = norme*fSumELongit[i]; 216 rmsELongit[i] = << 225 rmsELongit[i] = 217 norme * std::sqrt(std::abs(NbOfEvents * << 226 norme*std::sqrt(std::abs(NbOfEvents*fSumE2Longit[i] 218 << 227 - fSumELongit[i]*fSumELongit[i])); 219 MeanELongitCumul[i] = norme * fSumELongitC << 228 220 rmsELongitCumul[i] = norme << 229 MeanELongitCumul[i] = norme*fSumELongitCumul[i]; 221 * std::sqrt(std::abs( << 230 rmsELongitCumul[i] = norme*std::sqrt(std::abs(NbOfEvents* 222 << 231 fSumE2LongitCumul[i] - fSumELongitCumul[i]*fSumELongitCumul[i])); 223 G4double bin = (i + 0.5) * dLradl; << 232 G4double bin = (i+0.5)*dLradl; 224 analysisManager->FillH1(4, bin, MeanELongi << 233 analysisManager->FillH1(4, bin,MeanELongit[i]/dLradl); 225 analysisManager->FillH1(5, bin, rmsELongit << 234 analysisManager->FillH1(5, bin, rmsELongit[i]/dLradl); 226 bin = (i + 1) * dLradl; << 235 bin = (i+1)*dLradl; 227 analysisManager->FillH1(6, bin, MeanELongi << 236 analysisManager->FillH1(6, bin,MeanELongitCumul[i]); 228 analysisManager->FillH1(7, bin, rmsELongit 237 analysisManager->FillH1(7, bin, rmsELongitCumul[i]); 229 } 238 } 230 239 231 // radial << 240 //radial 232 // 241 // 233 G4double dRradl = fDet->GetdRradl(); 242 G4double dRradl = fDet->GetdRradl(); 234 243 235 MyVector MeanERadial(f_nRbin), rmsERadial(f_ << 244 MyVector MeanERadial(f_nRbin), rmsERadial(f_nRbin); 236 MyVector MeanERadialCumul(f_nRbin), rmsERadi 245 MyVector MeanERadialCumul(f_nRbin), rmsERadialCumul(f_nRbin); 237 246 238 for (i = 0; i < f_nRbin; ++i) { << 247 for (i=0; i<f_nRbin; i++) { 239 MeanERadial[i] = norme * fSumERadial[i]; << 248 MeanERadial[i] = norme*fSumERadial[i]; 240 rmsERadial[i] = << 249 rmsERadial[i] = norme*std::sqrt(std::abs(NbOfEvents*fSumE2Radial[i] 241 norme * std::sqrt(std::abs(NbOfEvents * << 250 - fSumERadial[i]*fSumERadial[i])); 242 << 251 243 MeanERadialCumul[i] = norme * fSumERadialC << 252 MeanERadialCumul[i] = norme*fSumERadialCumul[i]; 244 rmsERadialCumul[i] = norme << 253 rmsERadialCumul[i] = 245 * std::sqrt(std::abs( << 254 norme*std::sqrt(std::abs(NbOfEvents*fSumE2RadialCumul[i] 246 << 255 - fSumERadialCumul[i]*fSumERadialCumul[i])); 247 << 256 248 G4double bin = (i + 0.5) * dRradl; << 257 G4double bin = (i+0.5)*dRradl; 249 analysisManager->FillH1(8, bin, MeanERadia << 258 analysisManager->FillH1(8, bin,MeanERadial[i]/dRradl); 250 analysisManager->FillH1(9, bin, rmsERadial << 259 analysisManager->FillH1(9, bin, rmsERadial[i]/dRradl); 251 bin = (i + 1) * dRradl; << 260 bin = (i+1)*dRradl; 252 analysisManager->FillH1(10, bin, MeanERadi << 261 analysisManager->FillH1(10, bin,MeanERadialCumul[i]); 253 analysisManager->FillH1(11, bin, rmsERadia << 262 analysisManager->FillH1(11, bin, rmsERadialCumul[i]); 254 } 263 } 255 264 256 // find Moliere confinement << 265 //find Moliere confinement 257 // 266 // 258 const G4double EMoliere = 90.; 267 const G4double EMoliere = 90.; 259 G4double iMoliere = 0.; 268 G4double iMoliere = 0.; 260 if ((MeanERadialCumul[0] <= EMoliere) && (Me << 269 if ((MeanERadialCumul[0] <= EMoliere) && >> 270 (MeanERadialCumul[f_nRbin-1] >= EMoliere)) { 261 G4int imin = 0; 271 G4int imin = 0; 262 while ((imin < f_nRbin - 1) && (MeanERadia << 272 while( (imin < f_nRbin-1) && (MeanERadialCumul[imin] < EMoliere) ) 263 ++imin; << 273 { ++imin; } 264 } << 274 G4double ratio = (EMoliere - MeanERadialCumul[imin]) / 265 << 275 (MeanERadialCumul[imin+1] - MeanERadialCumul[imin]); 266 G4double del = MeanERadialCumul[imin + 1] << 267 G4double ratio = (del > 0.0) ? (EMoliere - << 268 iMoliere = 1. + imin + ratio; 276 iMoliere = 1. + imin + ratio; 269 } << 277 } 270 << 278 271 // track length << 279 //track length 272 // 280 // 273 norme = 1. / (NbOfEvents * (fDet->GetMateria << 281 norme = 1./(NbOfEvents*(fDet->GetMaterial()->GetRadlen())); 274 G4double MeanChargTrLength = norme * fSumCha << 282 G4double MeanChargTrLength = norme*fSumChargTrLength; 275 G4double rmsChargTrLength = << 283 G4double rmsChargTrLength = 276 norme << 284 norme*std::sqrt(std::fabs(NbOfEvents*fSum2ChargTrLength 277 * std::sqrt(std::abs(NbOfEvents * fSum2Cha << 285 - fSumChargTrLength*fSumChargTrLength)); 278 << 286 279 G4double MeanNeutrTrLength = norme * fSumNeu << 287 G4double MeanNeutrTrLength = norme*fSumNeutrTrLength; 280 G4double rmsNeutrTrLength = << 288 G4double rmsNeutrTrLength = 281 norme << 289 norme*std::sqrt(std::fabs(NbOfEvents*fSum2NeutrTrLength 282 * std::sqrt(std::abs(NbOfEvents * fSum2Neu << 290 - fSumNeutrTrLength*fSumNeutrTrLength)); 283 291 284 // print << 292 //print 285 std::ios::fmtflags mode = G4cout.flags(); 293 std::ios::fmtflags mode = G4cout.flags(); 286 G4cout.setf(std::ios::fixed, std::ios::float << 294 G4cout.setf(std::ios::fixed,std::ios::floatfield); 287 G4int prec = G4cout.precision(2); << 295 G4int prec = G4cout.precision(2); 288 << 296 289 if (fVerbose) { 297 if (fVerbose) { >> 298 290 G4cout << " LOGITUDINAL PR 299 G4cout << " LOGITUDINAL PROFILE " 291 << " CUMULATIVE LOGITUDINAL PR 300 << " CUMULATIVE LOGITUDINAL PROFILE" << G4endl << G4endl; 292 301 293 G4cout << " bin " << 302 G4cout << " bin " << " Mean rms " 294 << " Mean rms << 303 << " bin " << " Mean rms \n" << G4endl; 295 << " bin " << 304 296 << " Mean rms \n" << 305 for (i=0; i<f_nLbin; i++) { 297 << G4endl; << 306 G4double inf=i*dLradl, sup=inf+dLradl; 298 << 307 299 for (i = 0; i < f_nLbin; ++i) { << 308 G4cout << std::setw(8) << inf << "->" 300 G4double inf = i * dLradl, sup = inf + d << 309 << std::setw(5) << sup << " radl: " 301 << 310 << std::setw(7) << MeanELongit[i] << "% " 302 G4cout << std::setw(8) << inf << "->" << << 311 << std::setw(9) << rmsELongit[i] << "% " 303 << MeanELongit[i] << "% " << std << 312 << " 0->" << std::setw(5) << sup << " radl: " 304 << " 0->" << std::setw(5) << << 313 << std::setw(7) << MeanELongitCumul[i] << "% " 305 << MeanELongitCumul[i] << "% " < << 314 << std::setw(7) << rmsELongitCumul[i] << "% " 306 << G4endl; << 315 <<G4endl; 307 } 316 } 308 317 309 G4cout << G4endl << G4endl << G4endl; 318 G4cout << G4endl << G4endl << G4endl; 310 319 311 G4cout << " RADIAL PROFIL 320 G4cout << " RADIAL PROFILE " 312 << " CUMULATIVE RADIAL PROFIL 321 << " CUMULATIVE RADIAL PROFILE" << G4endl << G4endl; 313 322 314 G4cout << " bin " << 323 G4cout << " bin " << " Mean rms " 315 << " Mean rms << 324 << " bin " << " Mean rms \n" << G4endl; 316 << " bin " << 325 317 << " Mean rms \n" << 326 for (i=0; i<f_nRbin; i++) { 318 << G4endl; << 327 G4double inf=i*dRradl, sup=inf+dRradl; 319 << 320 for (i = 0; i < f_nRbin; ++i) { << 321 G4double inf = i * dRradl, sup = inf + d << 322 << 323 G4cout << std::setw(8) << inf << "->" << << 324 << MeanERadial[i] << "% " << std << 325 << " 0->" << std::setw(5) << << 326 << MeanERadialCumul[i] << "% " < << 327 << G4endl; << 328 } << 329 } << 330 328 >> 329 G4cout << std::setw(8) << inf << "->" >> 330 << std::setw(5) << sup << " radl: " >> 331 << std::setw(7) << MeanERadial[i] << "% " >> 332 << std::setw(9) << rmsERadial[i] << "% " >> 333 << " 0->" << std::setw(5) << sup << " radl: " >> 334 << std::setw(7) << MeanERadialCumul[i] << "% " >> 335 << std::setw(7) << rmsERadialCumul[i] << "% " >> 336 <<G4endl; >> 337 } >> 338 } >> 339 331 G4cout << "\n ===== SUMMARY ===== \n" << G4e 340 G4cout << "\n ===== SUMMARY ===== \n" << G4endl; 332 341 333 G4cout << " Total number of events: " 342 G4cout << " Total number of events: " << NbOfEvents << "\n" 334 << " Mean number of charged steps: " 343 << " Mean number of charged steps: " << fChargedStep << G4endl; 335 G4cout << " Mean number of neutral steps: " << 344 G4cout << " Mean number of neutral steps: " << fNeutralStep >> 345 << "\n" << G4endl; 336 346 337 G4cout << " energy deposit : " << std::setw( << 347 G4cout << " energy deposit : " 338 << std::setw(7) << rmsELongitCumul[f_ << 348 << std::setw(7) << MeanELongitCumul[f_nLbin-1] << " % E0 +- " 339 G4cout << " charged traklen: " << std::setw( << 349 << std::setw(7) << rmsELongitCumul[f_nLbin-1] << " % E0" << G4endl; 340 << rmsChargTrLength << " radl" << G4e << 350 G4cout << " charged traklen: " 341 G4cout << " neutral traklen: " << std::setw( << 351 << std::setw(7) << MeanChargTrLength << " radl +- " 342 << rmsNeutrTrLength << " radl" << G4e << 352 << std::setw(7) << rmsChargTrLength << " radl" << G4endl; 343 << 353 G4cout << " neutral traklen: " 344 if (iMoliere > 0.) { << 354 << std::setw(7) << MeanNeutrTrLength << " radl +- " 345 G4double RMoliere1 = iMoliere * fDet->Getd << 355 << std::setw(7) << rmsNeutrTrLength << " radl" << G4endl; 346 G4double RMoliere2 = iMoliere * fDet->Getd << 356 347 G4cout << "\n " << EMoliere << " % confine << 357 if (iMoliere > 0. ) { 348 << G4BestUnit(RMoliere2, "Length") << 358 G4double RMoliere1 = iMoliere*fDet->GetdRradl(); 349 << "\n" << 359 G4double RMoliere2 = iMoliere*fDet->GetdRlength(); 350 << G4endl; << 360 G4cout << "\n " << EMoliere << " % confinement: radius = " 351 } << 361 << RMoliere1 << " radl (" >> 362 << G4BestUnit( RMoliere2, "Length") << ")" << "\n" << G4endl; >> 363 } 352 364 353 G4cout.setf(mode, std::ios::floatfield); << 365 G4cout.setf(mode,std::ios::floatfield); 354 G4cout.precision(prec); 366 G4cout.precision(prec); 355 367 356 // Acceptance 368 // Acceptance 357 << 369 358 G4int nLbin = fDet->GetnLtot(); << 370 G4int nLbin = fDet->GetnLtot(); 359 if (limit < DBL_MAX) { << 371 if (limit < DBL_MAX) { 360 EmAcceptance acc; 372 EmAcceptance acc; 361 acc.BeginOfAcceptance("Total Energy in Abs << 373 acc.BeginOfAcceptance("Total Energy in Absorber",NbOfEvents); 362 G4double e = MeanELongitCumul[nLbin - 1] / << 374 G4double e = MeanELongitCumul[nLbin-1]/100.; 363 G4double r = rmsELongitCumul[nLbin - 1] / << 375 G4double r = rmsELongitCumul[nLbin-1]/100.; 364 acc.EmAcceptanceGauss("Edep", NbOfEvents, << 376 acc.EmAcceptanceGauss("Edep",NbOfEvents,e,edep,rms,limit); 365 acc.EmAcceptanceGauss("Erms", NbOfEvents, << 377 acc.EmAcceptanceGauss("Erms",NbOfEvents,r,rms,rms,2.0*limit); 366 acc.EndOfAcceptance(); 378 acc.EndOfAcceptance(); 367 } 379 } 368 limit = DBL_MAX; 380 limit = DBL_MAX; 369 } 381 } 370 382 371 //....oooOO0OOooo........oooOO0OOooo........oo 383 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 372 384