31 #include "RunAction.hh"
33 #include "PrimaryGeneratorAction.hh"
34 #include "HistoManager.hh"
47 :detector(det), primary(prim), histoManager(hist)
71 G4int size = totalEnergy.size();
72 if (size < nbPixels) {
73 visibleEnergy.resize(nbPixels);
74 totalEnergy.resize(nbPixels);
76 visibleEnergy2.resize(nbPixels);
77 totalEnergy2.resize(nbPixels);
80 for (
G4int k=0; k<nbPixels; k++) {
81 visibleEnergy[k] = visibleEnergy2[k] = totalEnergy[k]= totalEnergy2[k] = 0.0;
85 size = layerEtot.size();
87 layerEvis.resize(n1pxl);
88 layerEtot.resize(n1pxl);
89 layerEvis2.resize(n1pxl);
90 layerEtot2.resize(n1pxl);
93 for (
G4int k=0; k<n1pxl; k++) {
94 layerEvis[k] = layerEvis2[k] = layerEtot[k]= layerEtot2[k] = 0.0;
98 calorEvis = calorEvis2 = calorEtot = calorEtot2 = Eleak = Eleak2 = 0.;
99 EdLeak[0] = EdLeak[1] = EdLeak[2] = 0.;
100 nbRadLen = nbRadLen2 = 0.;
104 histoManager->
book();
117 visibleEnergy[pixel] += Evis; visibleEnergy2[pixel] += Evis*Evis;
118 totalEnergy[pixel] += Etot; totalEnergy2[pixel] += Etot*Etot;
127 layerEvis[layer] += Evis; layerEvis2[layer] += Evis*Evis;
128 layerEtot[layer] += Etot; layerEtot2[layer] += Etot*Etot;
139 calorEvis += calEvis; calorEvis2 += calEvis*calEvis;
140 calorEtot += calEtot; calorEtot2 += calEtot*calEtot;
141 Eleak += eleak; Eleak2 += eleak*eleak;
159 nbRadLen += dn; nbRadLen2 += dn*dn;
180 G4cout <<
" The run was " << nbEvents <<
" " << partName <<
" of "
183 G4cout <<
"------------------------------------------------------------"
188 if (nbEvents == 0)
return;
192 std::ios::fmtflags mode =
G4cout.flags();
197 G4double meanNbRadL = nbRadLen/ nbEvents;
198 G4double meanNbRadL2 = nbRadLen2/nbEvents;
199 G4double varNbRadL = meanNbRadL2 - meanNbRadL*meanNbRadL;
201 if (varNbRadL > 0.) rmsNbRadL = std::sqrt(varNbRadL);
205 <<
"\n Calor : mean number of Rad Length = "
206 << meanNbRadL <<
" +- "<< rmsNbRadL
207 <<
" --> Effective Rad Length = "
210 G4cout <<
"------------------------------------------------------------"
216 <<
"visible Energy (rms/mean) "
217 <<
"total Energy (rms/mean)" <<
G4endl;
219 G4double meanEvis,meanEvis2,varianceEvis,rmsEvis,resEvis;
220 G4double meanEtot,meanEtot2,varianceEtot,rmsEtot,resEtot;
224 for (
G4int i1=0; i1<n1pxl; i1++) {
226 meanEvis = layerEvis[i1] /nbEvents;
227 meanEvis2 = layerEvis2[i1]/nbEvents;
228 varianceEvis = meanEvis2 - meanEvis*meanEvis;
230 if (varianceEvis > 0.) rmsEvis = std::sqrt(varianceEvis);
231 resEvis = 100*rmsEvis/meanEvis;
232 histoManager->
FillHisto(3, i1+0.5, meanEvis);
235 meanEtot = layerEtot[i1] /nbEvents;
236 meanEtot2 = layerEtot2[i1]/nbEvents;
237 varianceEtot = meanEtot2 - meanEtot*meanEtot;
239 if (varianceEtot > 0.) rmsEtot = std::sqrt(varianceEtot);
240 resEtot = 100*rmsEtot/meanEtot;
241 histoManager->
FillHisto(4, i1+0.5, meanEtot);
246 <<
"\n layer " << i1 <<
": "
247 << std::setprecision(5)
248 << std::setw(6) <<
G4BestUnit(meanEvis,
"Energy") <<
" +- "
249 << std::setprecision(4)
250 << std::setw(5) <<
G4BestUnit( rmsEvis,
"Energy") <<
" ("
251 << std::setprecision(2)
252 << std::setw(3) << resEvis <<
" %)"
254 << std::setprecision(5)
255 << std::setw(6) <<
G4BestUnit(meanEtot,
"Energy") <<
" +- "
256 << std::setprecision(4)
257 << std::setw(5) <<
G4BestUnit( rmsEtot,
"Energy") <<
" ("
258 << std::setprecision(2)
259 << std::setw(3) << resEtot <<
" %)";
264 meanEvis = calorEvis /nbEvents;
265 meanEvis2 = calorEvis2/nbEvents;
266 varianceEvis = meanEvis2 - meanEvis*meanEvis;
268 if (varianceEvis > 0.) rmsEvis = std::sqrt(varianceEvis);
269 resEvis = 100*rmsEvis/meanEvis;
272 meanEtot = calorEtot /nbEvents;
273 meanEtot2 = calorEtot2/nbEvents;
274 varianceEtot = meanEtot2 - meanEtot*meanEtot;
276 if (varianceEtot > 0.) rmsEtot = std::sqrt(varianceEtot);
277 resEtot = 100*rmsEtot/meanEtot;
282 <<
"\n total calor : "
283 << std::setprecision(5)
284 << std::setw(6) <<
G4BestUnit(meanEvis,
"Energy") <<
" +- "
285 << std::setprecision(4)
286 << std::setw(5) <<
G4BestUnit( rmsEvis,
"Energy") <<
" ("
287 << std::setprecision(2)
288 << std::setw(3) << resEvis <<
" %)"
290 << std::setprecision(5)
291 << std::setw(6) <<
G4BestUnit(meanEtot,
"Energy") <<
" +- "
292 << std::setprecision(4)
293 << std::setw(5) <<
G4BestUnit( rmsEtot,
"Energy") <<
" ("
294 << std::setprecision(2)
295 << std::setw(3) << resEtot <<
" %)";
297 G4cout <<
"\n------------------------------------------------------------"
301 G4double meanEleak,meanEleak2,varianceEleak,rmsEleak,ratio;
302 meanEleak = Eleak /nbEvents;
303 meanEleak2 = Eleak2/nbEvents;
304 varianceEleak = meanEleak2 - meanEleak*meanEleak;
306 if (varianceEleak > 0.) rmsEleak = std::sqrt(varianceEleak);
307 ratio = 100*meanEleak/
energy;
316 << std::setprecision(5)
317 << std::setw(6) <<
G4BestUnit(meanEleak,
"Energy") <<
" +- "
318 << std::setprecision(4)
319 << std::setw(5) <<
G4BestUnit( rmsEleak,
"Energy")
320 <<
"\n Eleak/Ebeam ="
321 << std::setprecision(3)
322 << std::setw(4) << ratio <<
" % ( forward ="
323 << std::setw(4) << forward <<
" %; backward ="
324 << std::setw(4) << bakward <<
" %; lateral ="
325 << std::setw(4) << lateral <<
" %)"
328 G4cout.setf(mode,std::ios::floatfield);
332 histoManager->
save();
347 G4String fileName = name +
".pixels.ascii";
349 std::ofstream File(fileName, std::ios::out);
356 File << noEvents <<
" " << n1pxl <<
" " << n2pxl <<
" " << n1shift
void fillNbRadLen(G4double)
void BeginOfRunAction(const G4Run *)
void FillHisto(G4int id, G4double bin, G4double weight=1.0)
G4double GetCalorThickness()
#define G4BestUnit(a, b)
#define G4_USE_G4BESTUNIT_FOR_VERBOSE 1
void SetRandomNumberStore(G4bool flag)
void fillPerEvent_1(G4int, G4double, G4double)
void fillPerEvent_2(G4int, G4double, G4double)
const G4Run * GetCurrentRun() const
const G4String & GetParticleName() const
G4GLOB_DLL std::ostream G4cout
void EndOfRunAction(const G4Run *)
static void showEngineStatus()
static G4RunManager * GetRunManager()
G4int GetSizeVectorPixels()
void fillDetailedLeakage(G4int, G4double)
G4ParticleGun * GetParticleGun()
G4ParticleDefinition * GetParticleDefinition() const
static G4Geantino * Geantino()
G4int GetNumberOfEventToBeProcessed() const
void fillPerEvent_3(G4double, G4double, G4double)
void PrintCalorParameters()
G4double GetParticleEnergy() const