35 #include "DetectorConstruction.hh"
36 #include "PrimaryGeneratorAction.hh"
37 #include "HistoManager.hh"
50 fParticle(0), fEkin(0.)
89 const Run* localRun =
static_cast<const Run*
>(run);
139 if (TotNbofEvents == 0)
return;
142 EnergyBalance /= TotNbofEvents;
146 if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep/TotNbofEvents);
151 if (rmsTLCh>0.) rmsTLCh = std::sqrt(rmsTLCh/TotNbofEvents);
156 if (rmsTLNe>0.) rmsTLNe = std::sqrt(rmsTLNe/TotNbofEvents);
161 if (rmsStCh>0.) rmsStCh = std::sqrt(rmsStCh/TotNbofEvents);
166 if (rmsStNe>0.) rmsStNe = std::sqrt(rmsStNe/TotNbofEvents);
174 transmit[0] = 100.*
fTransmit[0]/TotNbofEvents;
175 transmit[1] = 100.*
fTransmit[1]/TotNbofEvents;
178 reflect[0] = 100.*
fReflect[0]/TotNbofEvents;
179 reflect[1] = 100.*
fReflect[1]/TotNbofEvents;
185 if (rmsMsc > 0.) { rmsMsc = std::sqrt(rmsMsc); }
191 fEnergyLeak[0] /= TotNbofEvents;
fEnergyLeak2[0] /= TotNbofEvents;
193 if (rmsEl0>0.) rmsEl0 = std::sqrt(rmsEl0/TotNbofEvents);
196 fEnergyLeak[1] /= TotNbofEvents;
fEnergyLeak2[1] /= TotNbofEvents;
198 if (rmsEl1>0.) rmsEl1 = std::sqrt(rmsEl1/TotNbofEvents);
210 G4double dEdxTable = 0., dEdxFull = 0.;
215 G4double stopTable = dEdxTable/density;
216 G4double stopFull = dEdxFull /density;
220 G4double meandEdx = fEnergyDeposit/length;
221 G4double stopPower = meandEdx/density;
223 G4cout <<
"\n ======================== run summary ======================\n";
227 G4cout <<
"\n The run was " << TotNbofEvents <<
" " << partName <<
" of "
230 << material->
GetName() <<
" (density: "
235 G4cout <<
"\n Total energy deposit in absorber per event = "
236 <<
G4BestUnit(fEnergyDeposit,
"Energy") <<
" +- "
240 G4cout <<
"\n -----> Mean dE/dx = " << meandEdx/(
MeV/
cm) <<
" MeV/cm"
241 <<
"\t(" << stopPower/(
MeV*
cm2/
g) <<
" MeV*cm2/g)"
245 G4cout <<
" restricted dEdx = " << dEdxTable/(
MeV/
cm) <<
" MeV/cm"
246 <<
"\t(" << stopTable/(
MeV*
cm2/
g) <<
" MeV*cm2/g)"
249 G4cout <<
" full dEdx = " << dEdxFull/(
MeV/
cm) <<
" MeV/cm"
250 <<
"\t(" << stopFull/(
MeV*
cm2/
g) <<
" MeV*cm2/g)"
253 G4cout <<
"\n Leakage : primary = "
254 <<
G4BestUnit(fEnergyLeak[0],
"Energy") <<
" +- "
257 <<
G4BestUnit(fEnergyLeak[1],
"Energy") <<
" +- "
261 G4cout <<
" Energy balance : edep + eleak = "
265 G4cout <<
"\n Total track length (charged) in absorber per event = "
266 <<
G4BestUnit(fTrakLenCharged,
"Length") <<
" +- "
269 G4cout <<
" Total track length (neutral) in absorber per event = "
270 <<
G4BestUnit(fTrakLenNeutral,
"Length") <<
" +- "
273 G4cout <<
"\n Number of steps (charged) in absorber per event = "
274 << fNbStepsCharged <<
" +- " << rmsStCh <<
G4endl;
276 G4cout <<
" Number of steps (neutral) in absorber per event = "
277 << fNbStepsNeutral <<
" +- " << rmsStNe <<
G4endl;
279 G4cout <<
"\n Number of secondaries per event : Gammas = " << Gamma
280 <<
"; electrons = " << Elect
281 <<
"; positrons = " << Posit <<
G4endl;
283 G4cout <<
"\n Number of events with the primary particle transmitted = "
284 << transmit[1] <<
" %" <<
G4endl;
286 G4cout <<
" Number of events with at least 1 particle transmitted "
287 <<
"(same charge as primary) = " << transmit[0] <<
" %" <<
G4endl;
289 G4cout <<
"\n Number of events with the primary particle reflected = "
290 << reflect[1] <<
" %" <<
G4endl;
292 G4cout <<
" Number of events with at least 1 particle reflected "
293 <<
"(same charge as primary) = " << reflect[0] <<
" %" <<
G4endl;
297 G4cout <<
"\n MultipleScattering:"
298 <<
"\n rms proj angle of transmit primary particle = "
299 << rmsMsc/
mrad <<
" mrad (central part only)" <<
G4endl;
301 G4cout <<
" computed theta0 (Highland formula) = "
304 G4cout <<
" central part defined as +- "
306 <<
" Tail ratio = " << tailMsc <<
" %" <<
G4endl;
313 G4double binWidth = analysisManager->GetH1Width(ih);
315 analysisManager->ScaleH1(ih,fac);
318 binWidth = analysisManager->GetH1Width(ih);
319 fac = 1./(TotNbofEvents*binWidth);
320 analysisManager->ScaleH1(ih,fac);
323 analysisManager->ScaleH1(ih,1./TotNbofEvents);
346 G4double teta0 = 13.6*
MeV*z*std::sqrt(t)*(1.+0.038*std::log(t))/bpc;
void SetPrimary(G4ParticleDefinition *particle, G4double energy)
DetectorConstruction * fDetector
static constexpr double MeV
std::vector< G4double > fEnergyDeposit[kMaxAbsor]
G4double fMscProjecTheta2
G4double GetAbsorberThickness()
G4double GetRadlen() const
G4double GetDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
const G4String & GetParticleName() const
G4double fTrakLenCharged2
G4double GetPDGCharge() const
G4double GetPDGMass() const
G4double fTrakLenNeutral2
G4double fMscThetaCentral
G4Material * GetAbsorberMaterial()
const G4String & GetName() const
static constexpr double g
static constexpr double cm2
G4double ComputeTotalDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, G4double cut=DBL_MAX)
virtual void Merge(const G4Run *)
G4CsvAnalysisManager G4AnalysisManager
G4ParticleDefinition * fParticle
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
static constexpr double eplus
static const G4double fac
G4double ComputeMscHighland()
static constexpr double cm
G4GLOB_DLL std::ostream G4cout
static constexpr double mrad
G4double fNbStepsCharged2
virtual void Merge(const G4Run *)
Simple detector construction with a box volume placed in a world.
G4double fNbStepsNeutral2
G4double GetDensity() const