34 #include "RunAction.hh"
36 #include "PrimaryGeneratorAction.hh"
37 #include "RunActionMessenger.hh"
38 #include "HistoManager.hh"
39 #include "EmAcceptance.hh"
68 fChargedStep = fNeutralStep = 0.0;
70 for (
G4int k=0; k<
MaxAbsor; k++) { fEdeptrue[k] = fRmstrue[k] = 1.;
96 fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = fSum2LAbs[k] = 0.;
97 fEnergyDeposit[k].clear();
100 fChargedStep = fNeutralStep = 0.0;
109 fEnergyFlow.resize(nbPlanes);
110 fLateralEleak.resize(nbPlanes);
111 for (
G4int k=0; k<nbPlanes; k++) {fEnergyFlow[k] = fLateralEleak[k] = 0.; }
116 if (analysis->IsActive()) analysis->OpenFile();
129 if(fApplyLimit) fEnergyDeposit[
kAbs].push_back(EAbs);
130 fSumEAbs[
kAbs] += EAbs; fSum2EAbs[
kAbs] += EAbs*EAbs;
131 fSumLAbs[
kAbs] += LAbs; fSum2LAbs[
kAbs] += LAbs*LAbs;
141 if(norm > 0) norm = 1./norm;
144 fChargedStep *= norm;
145 fNeutralStep *= norm;
152 G4double MeanEAbs,MeanEAbs2,rmsEAbs,resolution,rmsres;
153 G4double MeanLAbs,MeanLAbs2,rmsLAbs;
155 std::ios::fmtflags mode =
G4cout.flags();
157 G4cout <<
"\n------------------------------------------------------------\n";
158 G4cout << std::setw(14) <<
"material"
159 << std::setw(17) <<
"Edep RMS"
160 << std::setw(33) <<
"sqrt(E0(GeV))*rmsE/Emean"
161 << std::setw(23) <<
"total tracklen \n \n";
165 MeanEAbs = fSumEAbs[k]*norm;
166 MeanEAbs2 = fSum2EAbs[k]*norm;
167 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
176 for(
G4int i=0; i<nEvt; i++) {
177 G4double e = (fEnergyDeposit[k])[i];
178 if(std::abs(e - MeanEAbs) < lim) {
185 if(norm1 > 0.0) norm1 = 1.0/norm1;
186 MeanEAbs = sume*norm1;
187 MeanEAbs2 = sume2*norm1;
188 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
191 resolution= 100.*sqbeam*rmsEAbs/MeanEAbs;
192 rmsres = resolution*qnorm;
195 fSumEAbs[k] = MeanEAbs;
196 fSum2EAbs[k] = rmsEAbs;
198 MeanLAbs = fSumLAbs[k]*norm;
199 MeanLAbs2 = fSum2LAbs[k]*norm;
200 rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - MeanLAbs*MeanLAbs));
206 << std::setprecision(5)
207 << std::setw(6) <<
G4BestUnit(MeanEAbs,
"Energy") <<
" : "
208 << std::setprecision(4)
209 << std::setw(5) <<
G4BestUnit( rmsEAbs,
"Energy")
210 << std::setw(10) << resolution <<
" +- "
211 << std::setw(5) << rmsres <<
" %"
212 << std::setprecision(3)
213 << std::setw(10) <<
G4BestUnit(MeanLAbs,
"Length") <<
" +- "
214 << std::setw(4) <<
G4BestUnit( rmsLAbs,
"Length")
217 G4cout <<
"\n------------------------------------------------------------\n";
219 G4cout <<
" Beam particle "
221 GetParticleDefinition()->GetParticleName()
226 G4cout << std::setprecision(6)
227 <<
" Mean number of charged steps " << fChargedStep <<
G4endl;
228 G4cout <<
" Mean number of neutral steps " << fNeutralStep <<
G4endl;
229 G4cout <<
"------------------------------------------------------------\n";
235 for (
G4int Id=1; Id<=Idmax+1; Id++) {
236 analysis->FillH1(2*MaxAbsor+1, (
G4double)Id, fEnergyFlow[Id]);
237 analysis->FillH1(2*MaxAbsor+2, (
G4double)Id, fLateralEleak[Id]);
246 for (
G4int Id=1; Id<=Idmax; Id++) {
247 G4int iAbsor = Id%nbOfAbsor;
if (iAbsor==0) iAbsor = nbOfAbsor;
248 EdepTot [iAbsor] += (fEnergyFlow[Id] - fEnergyFlow[Id+1] - fLateralEleak[Id]);
251 G4cout << std::setprecision(3)
252 <<
"\n Energy deposition from Energy flow balance : \n"
253 << std::setw(10) <<
" material \t Total Edep \n \n";
256 for (
G4int k=1; k<=nbOfAbsor; k++) {
259 <<
"\t " <<
G4BestUnit(EdepTot [k],
"Energy") <<
"\n";
262 G4cout <<
"\n------------------------------------------------------------\n"
265 G4cout.setf(mode,std::ios::floatfield);
277 MeanEAbs = fSumEAbs[j];
278 rmsEAbs = fSum2EAbs[j];
281 fEdeptrue[j], fRmstrue[j], fLimittrue[j]);
283 fRmstrue[j], fRmstrue[j], 2.0*fLimittrue[j]);
291 analysis->ScaleH1(ih,norm/
MeV);
295 if (analysis->IsActive()) {
297 analysis->CloseFile();
315 const G4int ncolumn = 5;
319 const G4double dp = std::log10(tkmax/tkmin)/nbin;
320 const G4double dt = std::pow(10.,dp);
322 for (
G4int i=1; i<nbin; ++i) tk[i] = tk[i-1]*dt;
326 std::ios::fmtflags mode =
G4cout.flags();
327 G4cout.setf(std::ios::fixed,std::ios::floatfield);
330 G4cout <<
"\n kinetic energies \n ";
331 for (
G4int j=0; j<nbin; ++j) {
333 if ((j+1)%ncolumn == 0)
G4cout <<
"\n ";
339 G4cout.setf(std::ios::scientific,std::ios::floatfield);
353 for (
size_t i=0; i<numOfCouples; i++) {
355 if (couple->
GetMaterial() == mat) {index = i;
break;}
359 <<
" in " << mat ->
GetName() <<
"\nC";
363 G4cout <<
"\nERAN " << tkmin/
GeV <<
" (ekmin)\t"
364 << tkmax/
GeV <<
" (ekmax)\t"
365 << nbin <<
" (nekbin)";
368 if (cutgam < tkmin) cutgam = tkmin;
369 if (cutgam > tkmax) cutgam = tkmax;
372 if (cutele < tkmin) cutele = tkmin;
373 if (cutele > tkmax) cutele = tkmax;
374 G4cout <<
"\nCUTS " << cutgam/
GeV <<
" (cutgam)\t"
375 << cutele/
GeV <<
" (cutele)";
379 for (
G4int l=0;l<nbin; ++l)
384 if ((l+1)%ncolumn == 0)
G4cout <<
"\n ";
390 G4cout.setf(mode,std::ios::floatfield);
407 if (i>=0 && i<MaxAbsor) {
408 fEdeptrue [i] = edep;
G4ParticleDefinition * GetDefinition() const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
static G4LossTableManager * Instance()
void BeginOfRunAction(const G4Run *)
G4double GetDEDX(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
void FillPerEvent(G4int, G4double, G4double)
void EmAcceptanceGauss(const G4String &title, G4int stat, G4double avr, G4double avr0, G4double rms, G4double limit)
const G4String & GetName() const
void AddSecondaryTrack(const G4Track *)
void BeginOfAcceptance(const G4String &title, G4int stat)
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
const G4String & GetParticleName() const
G4Material * GetAbsorMaterial(G4int i)
G4GLOB_DLL std::ostream G4cout
size_t GetTableSize() const
G4int GetNumberOfEvent() const
void SetEdepAndRMS(G4ThreeVector)
void EndOfRunAction(const G4Run *)
ExG4HbookAnalysisManager G4AnalysisManager
static void showEngineStatus()
static G4ProductionCutsTable * GetProductionCutsTable()
static G4Positron * Positron()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
G4ParticleGun * GetParticleGun()
G4ParticleDefinition * GetParticleDefinition() const
static G4Electron * Electron()
G4double GetParticleEnergy() const
const G4Material * GetMaterial() const