35 #include "DetectorConstruction.hh"
36 #include "PrimaryGeneratorAction.hh"
37 #include "HistoManager.hh"
50 fParticle(0), fEkin(0.)
52 fEnergyDeposit = fEnergyDeposit2 = 0.;
53 fTrakLenCharged = fTrakLenCharged2 = 0.;
54 fTrakLenNeutral = fTrakLenNeutral2 = 0.;
55 fNbStepsCharged = fNbStepsCharged2 = 0.;
56 fNbStepsNeutral = fNbStepsNeutral2 = 0.;
57 fMscProjecTheta = fMscProjecTheta2 = 0.;
58 fMscThetaCentral = 0.;
60 fNbGamma = fNbElect = fNbPosit = 0;
62 fTransmit[0] = fTransmit[1] = fReflect[0] = fReflect[1] = 0;
66 fEnergyLeak[0] = fEnergyLeak[1] = fEnergyLeak2[0] = fEnergyLeak2[1] = 0.;
89 const Run* localRun =
static_cast<const Run*
>(run);
92 fParticle = localRun->fParticle;
93 fEkin = localRun->fEkin;
95 fMscThetaCentral = localRun->fMscThetaCentral;
99 fEnergyDeposit += localRun->fEnergyDeposit;
100 fEnergyDeposit2 += localRun->fEnergyDeposit2;
101 fTrakLenCharged += localRun->fTrakLenCharged;
102 fTrakLenCharged2 += localRun->fTrakLenCharged2;
103 fTrakLenNeutral += localRun->fTrakLenNeutral;
104 fTrakLenNeutral2 += localRun->fTrakLenNeutral2;
105 fNbStepsCharged += localRun->fNbStepsCharged;
106 fNbStepsCharged2 += localRun->fNbStepsCharged2;
107 fNbStepsNeutral += localRun->fNbStepsNeutral;
108 fNbStepsNeutral2 += localRun->fNbStepsNeutral2;
109 fMscProjecTheta += localRun->fMscProjecTheta;
110 fMscProjecTheta2 += localRun->fMscProjecTheta2;
113 fNbGamma += localRun->fNbGamma;
114 fNbElect += localRun->fNbElect;
115 fNbPosit += localRun->fNbPosit;
117 fTransmit[0] += localRun->fTransmit[0];
118 fTransmit[1] += localRun->fTransmit[1];
119 fReflect[0] += localRun->fReflect[0];
120 fReflect[1] += localRun->fReflect[1];
122 fMscEntryCentral += localRun->fMscEntryCentral;
124 fEnergyLeak[0] += localRun->fEnergyLeak[0];
125 fEnergyLeak[1] += localRun->fEnergyLeak[1];
126 fEnergyLeak2[0] += localRun->fEnergyLeak2[0];
127 fEnergyLeak2[1] += localRun->fEnergyLeak2[1];
139 if (TotNbofEvents == 0)
return;
141 G4double EnergyBalance = fEnergyDeposit + fEnergyLeak[0] + fEnergyLeak[1];
142 EnergyBalance /= TotNbofEvents;
144 fEnergyDeposit /= TotNbofEvents; fEnergyDeposit2 /= TotNbofEvents;
145 G4double rmsEdep = fEnergyDeposit2 - fEnergyDeposit*fEnergyDeposit;
146 if (rmsEdep>0.) rmsEdep = std::sqrt(rmsEdep/TotNbofEvents);
149 fTrakLenCharged /= TotNbofEvents; fTrakLenCharged2 /= TotNbofEvents;
150 G4double rmsTLCh = fTrakLenCharged2 - fTrakLenCharged*fTrakLenCharged;
151 if (rmsTLCh>0.) rmsTLCh = std::sqrt(rmsTLCh/TotNbofEvents);
154 fTrakLenNeutral /= TotNbofEvents; fTrakLenNeutral2 /= TotNbofEvents;
155 G4double rmsTLNe = fTrakLenNeutral2 - fTrakLenNeutral*fTrakLenNeutral;
156 if (rmsTLNe>0.) rmsTLNe = std::sqrt(rmsTLNe/TotNbofEvents);
159 fNbStepsCharged /= TotNbofEvents; fNbStepsCharged2 /= TotNbofEvents;
160 G4double rmsStCh = fNbStepsCharged2 - fNbStepsCharged*fNbStepsCharged;
161 if (rmsStCh>0.) rmsStCh = std::sqrt(rmsStCh/TotNbofEvents);
164 fNbStepsNeutral /= TotNbofEvents; fNbStepsNeutral2 /= TotNbofEvents;
165 G4double rmsStNe = fNbStepsNeutral2 - fNbStepsNeutral*fNbStepsNeutral;
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;
182 if (fMscEntryCentral > 0) {
183 fMscProjecTheta /= fMscEntryCentral; fMscProjecTheta2 /= fMscEntryCentral;
184 rmsMsc = fMscProjecTheta2 - fMscProjecTheta*fMscProjecTheta;
185 if (rmsMsc > 0.) { rmsMsc = std::sqrt(rmsMsc); }
186 if(fTransmit[1] > 0.0) {
187 tailMsc = 100.- (100.*fMscEntryCentral)/(2*fTransmit[1]);
191 fEnergyLeak[0] /= TotNbofEvents; fEnergyLeak2[0] /= TotNbofEvents;
192 G4double rmsEl0 = fEnergyLeak2[0] - fEnergyLeak[0]*fEnergyLeak[0];
193 if (rmsEl0>0.) rmsEl0 = std::sqrt(rmsEl0/TotNbofEvents);
196 fEnergyLeak[1] /= TotNbofEvents; fEnergyLeak2[1] /= TotNbofEvents;
197 G4double rmsEl1 = fEnergyLeak2[1] - fEnergyLeak[1]*fEnergyLeak[1];
198 if (rmsEl1>0.) rmsEl1 = std::sqrt(rmsEl1/TotNbofEvents);
210 G4double dEdxTable = 0., dEdxFull = 0.;
212 dEdxTable = emCalculator.
GetDEDX(fEkin,fParticle,material);
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 +- "
305 << fMscThetaCentral/
mrad <<
" mrad; "
306 <<
" Tail ratio = " << tailMsc <<
" %" <<
G4endl;
313 G4double binWidth = analysisManager->GetH1Width(ih);
314 G4double unit = analysisManager->GetH1Unit(ih);
316 analysisManager->ScaleH1(ih,fac);
319 binWidth = analysisManager->GetH1Width(ih);
320 unit = analysisManager->GetH1Unit(ih);
321 fac = unit/(TotNbofEvents*binWidth);
322 analysisManager->ScaleH1(ih,fac);
325 analysisManager->ScaleH1(ih,1./TotNbofEvents);
348 G4double teta0 = 13.6*
MeV*z*std::sqrt(t)*(1.+0.038*std::log(t))/bpc;
virtual void Merge(const G4Run *)
static constexpr double cm2
const G4String & GetName() const
G4double ComputeMscHighland()
G4double GetDensity() const
G4double ComputeTotalDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, G4double cut=DBL_MAX)
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
const G4String & GetParticleName() const
G4GLOB_DLL std::ostream G4cout
G4Material * GetAbsorberMaterial()
G4double GetDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
static constexpr double cm
static constexpr double eplus
static constexpr double mrad
G4double GetAbsorberThickness()
G4double GetRadlen() const
G4double GetPDGMass() const
G4double energy(const ThreeVector &p, const G4double m)
static const G4double fac
Detector construction class to define materials and geometry.
static constexpr double MeV
void SetPrimary(G4ParticleDefinition *particle, G4double energy)
G4double GetPDGCharge() const
virtual void Merge(const G4Run *)
G4CsvAnalysisManager G4AnalysisManager