34 #include "RunAction.hh"
36 #include "PrimaryGeneratorAction.hh"
37 #include "RunActionMessenger.hh"
38 #include "HistoManager.hh"
39 #include "EmAcceptance.hh"
61 :fDetector(det), fPrimary(prim), fRunMessenger(0), fHistoManager(0)
67 fChargedStep = fNeutralStep = 0.0;
69 for (
G4int k=0; k<
MaxAbsor; k++) { fEdeptrue[k] = fRmstrue[k] = 1.;
95 fSumEAbs[k] = fSum2EAbs[k] = fSumLAbs[k] = fSum2LAbs[k] = 0.;
96 fEnergyDeposit[k].clear();
99 fChargedStep = fNeutralStep = 0.0;
108 fEnergyFlow.resize(nbPlanes);
109 fLateralEleak.resize(nbPlanes);
110 for (
G4int k=0; k<nbPlanes; k++) {fEnergyFlow[k] = fLateralEleak[k] = 0.; }
115 if (analysis->IsActive()) analysis->OpenFile();
128 if(fApplyLimit) fEnergyDeposit[
kAbs].push_back(EAbs);
129 fSumEAbs[
kAbs] += EAbs; fSum2EAbs[
kAbs] += EAbs*EAbs;
130 fSumLAbs[
kAbs] += LAbs; fSum2LAbs[
kAbs] += LAbs*LAbs;
140 if(norm > 0) norm = 1./
norm;
143 fChargedStep *=
norm;
144 fNeutralStep *=
norm;
151 G4double MeanEAbs,MeanEAbs2,rmsEAbs,resolution,rmsres;
152 G4double MeanLAbs,MeanLAbs2,rmsLAbs;
154 std::ios::fmtflags mode =
G4cout.flags();
156 G4cout <<
"\n------------------------------------------------------------\n";
157 G4cout << std::setw(14) <<
"material"
158 << std::setw(17) <<
"Edep RMS"
159 << std::setw(33) <<
"sqrt(E0(GeV))*rmsE/Emean"
160 << std::setw(23) <<
"total tracklen \n \n";
164 MeanEAbs = fSumEAbs[k]*
norm;
165 MeanEAbs2 = fSum2EAbs[k]*
norm;
166 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
175 for(
G4int i=0; i<nEvt; i++) {
177 if(std::abs(e - MeanEAbs) < lim) {
184 if(norm1 > 0.0) norm1 = 1.0/norm1;
185 MeanEAbs = sume*norm1;
186 MeanEAbs2 = sume2*norm1;
187 rmsEAbs = std::sqrt(std::abs(MeanEAbs2 - MeanEAbs*MeanEAbs));
190 resolution= 100.*sqbeam*rmsEAbs/MeanEAbs;
191 rmsres = resolution*qnorm;
194 fSumEAbs[k] = MeanEAbs;
195 fSum2EAbs[k] = rmsEAbs;
197 MeanLAbs = fSumLAbs[k]*
norm;
198 MeanLAbs2 = fSum2LAbs[k]*
norm;
199 rmsLAbs = std::sqrt(std::abs(MeanLAbs2 - MeanLAbs*MeanLAbs));
205 << std::setprecision(5)
206 << std::setw(6) <<
G4BestUnit(MeanEAbs,
"Energy") <<
" : "
207 << std::setprecision(4)
208 << std::setw(5) <<
G4BestUnit( rmsEAbs,
"Energy")
209 << std::setw(10) << resolution <<
" +- "
210 << std::setw(5) << rmsres <<
" %"
211 << std::setprecision(3)
212 << std::setw(10) <<
G4BestUnit(MeanLAbs,
"Length") <<
" +- "
213 << std::setw(4) <<
G4BestUnit( rmsLAbs,
"Length")
216 G4cout <<
"\n------------------------------------------------------------\n";
218 G4cout <<
" Beam particle "
220 GetParticleDefinition()->GetParticleName()
225 G4cout << std::setprecision(6)
226 <<
" Mean number of charged steps " << fChargedStep <<
G4endl;
227 G4cout <<
" Mean number of neutral steps " << fNeutralStep <<
G4endl;
228 G4cout <<
"------------------------------------------------------------\n";
234 for (
G4int Id=1; Id<=Idmax+1; Id++) {
235 analysis->FillH1(2*MaxAbsor+1, (
G4double)Id, fEnergyFlow[Id]);
236 analysis->FillH1(2*MaxAbsor+2, (
G4double)Id, fLateralEleak[Id]);
245 for (
G4int Id=1; Id<=Idmax; Id++) {
246 G4int iAbsor = Id%nbOfAbsor;
if (iAbsor==0) iAbsor = nbOfAbsor;
247 EdepTot [iAbsor] += (fEnergyFlow[Id] - fEnergyFlow[Id+1] - fLateralEleak[Id]);
250 G4cout << std::setprecision(3)
251 <<
"\n Energy deposition from Energy flow balance : \n"
252 << std::setw(10) <<
" material \t Total Edep \n \n";
255 for (
G4int k=1; k<=nbOfAbsor; k++) {
258 <<
"\t " <<
G4BestUnit(EdepTot [k],
"Energy") <<
"\n";
261 G4cout <<
"\n------------------------------------------------------------\n"
264 G4cout.setf(mode,std::ios::floatfield);
276 MeanEAbs = fSumEAbs[j];
277 rmsEAbs = fSum2EAbs[j];
280 fEdeptrue[j], fRmstrue[j], fLimittrue[j]);
282 fRmstrue[j], fRmstrue[j], 2.0*fLimittrue[j]);
290 analysis->ScaleH1(ih,norm/
MeV);
294 if (analysis->IsActive()) {
296 analysis->CloseFile();
314 const G4int ncolumn = 5;
318 const G4double dp = std::log10(tkmax/tkmin)/nbin;
319 const G4double dt = std::pow(10.,dp);
321 for (
G4int i=1; i<nbin; ++i) tk[i] = tk[i-1]*dt;
325 std::ios::fmtflags mode =
G4cout.flags();
326 G4cout.setf(std::ios::fixed,std::ios::floatfield);
329 G4cout <<
"\n kinetic energies \n ";
330 for (
G4int j=0; j<nbin; ++j) {
332 if ((j+1)%ncolumn == 0)
G4cout <<
"\n ";
338 G4cout.setf(std::ios::scientific,std::ios::floatfield);
352 for (
size_t i=0; i<numOfCouples; i++) {
358 <<
" in " << mat ->
GetName() <<
"\nC";
362 G4cout <<
"\nERAN " << tkmin/
GeV <<
" (ekmin)\t"
363 << tkmax/
GeV <<
" (ekmax)\t"
364 << nbin <<
" (nekbin)";
367 if (cutgam < tkmin) cutgam = tkmin;
368 if (cutgam > tkmax) cutgam = tkmax;
371 if (cutele < tkmin) cutele = tkmin;
372 if (cutele > tkmax) cutele = tkmax;
373 G4cout <<
"\nCUTS " << cutgam/
GeV <<
" (cutgam)\t"
374 << cutele/
GeV <<
" (cutele)";
378 for (
G4int l=0;l<nbin; ++l)
383 if ((l+1)%ncolumn == 0)
G4cout <<
"\n ";
389 G4cout.setf(mode,std::ios::floatfield);
406 if (i>=0 && i<MaxAbsor) {
407 fEdeptrue [i] =
edep;