File indexing completed on 2025-02-23 09:19:40
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030 #include <iostream>
0031
0032 #include "CCalSteppingAction.hh"
0033 #include "CCalRunAction.hh"
0034
0035 #include "G4AnalysisManager.hh"
0036 #include "G4SystemOfUnits.hh"
0037 #include "G4SDManager.hh"
0038 #include "G4StepPoint.hh"
0039 #include "G4ThreeVector.hh"
0040 #include "G4RunManager.hh"
0041
0042 CCalSteppingAction::CCalSteppingAction()
0043 {
0044 timeHistoMaxBin = 200;
0045 for (G4int i=0; i<200; i++) {timeDeposit[i] = 0.f;}
0046 for (G4int i=0; i<70; i++) {LateralProfile[i] = 0.f;}
0047 }
0048
0049 CCalSteppingAction::~CCalSteppingAction(){
0050 }
0051
0052 void CCalSteppingAction::UserSteppingAction(const G4Step* aStep)
0053 {
0054 G4double de = aStep->GetTotalEnergyDeposit();
0055 if(de < CLHEP::eV) { return; }
0056
0057 const G4StepPoint* PostStepPoint= aStep->GetPostStepPoint();
0058 const G4StepPoint* PreStepPoint= aStep->GetPreStepPoint();
0059 G4double time =
0060 (PostStepPoint) ? PostStepPoint->GetGlobalTime()/nanosecond : 0.;
0061
0062 G4int it = (G4int)time;
0063 it = std::min(it, 10000);
0064
0065 G4int TSliceID = std::max(0,std::min(it,timeHistoMaxBin-1));
0066
0067 G4float fde = (G4float)(de/CLHEP::GeV);
0068
0069
0070 timeDeposit[TSliceID] += fde;
0071
0072 G4ThreeVector HitPoint = 0.5*(PostStepPoint->GetPosition()+
0073 PreStepPoint->GetPosition());
0074
0075
0076 G4double perp =
0077 std::sqrt(HitPoint.y()*HitPoint.y()+HitPoint.z()*HitPoint.z())/cm;
0078 G4int ir = (G4int)perp;
0079
0080 G4int radialPosition = std::max(0,std::min(69,ir));
0081 LateralProfile[radialPosition] += fde;
0082
0083 }
0084
0085 void CCalSteppingAction::endOfEvent(){
0086
0087 G4AnalysisManager* man = G4AnalysisManager::Instance();
0088 #ifdef debug
0089 G4double totalFilledProfileHcal = 0.0;
0090 #endif
0091 static G4int IDlateralProfile = -1;
0092 if (IDlateralProfile < 0) {
0093 IDlateralProfile = man->GetH1Id("h500");
0094 }
0095 for (G4int i=0; i<70; ++i) {
0096 man->FillH1(IDlateralProfile+i,LateralProfile[i]);
0097 #ifdef debug
0098 G4cout << "Fill Profile Hcal histo " << i << " with " << LateralProfile[i] << G4endl;
0099 totalFilledProfileHcal += LateralProfile[i];
0100 #endif
0101 }
0102
0103 #ifdef debug
0104 G4cout << "CCalAnalysis::InsertLateralProfile: Total filled Profile Hcal"
0105 << " histo " << totalFilledProfileHcal << G4endl;
0106 #endif
0107
0108 static G4int IDTimeHist = -1;
0109 if (IDTimeHist < 0) {
0110 IDTimeHist = man->GetH1Id("h300");
0111 }
0112 static G4int IDTimeProfile = -1;
0113 if (IDTimeProfile < 0) {
0114 IDTimeProfile = man->GetH1Id("h901");
0115 }
0116 #ifdef debug
0117 G4double totalFilledTimeProfile = 0.0;
0118 #endif
0119 for (G4int j=0; j<timeHistoMaxBin; ++j) {
0120 man->FillH1(IDTimeHist+j,timeDeposit[j]);
0121 #ifdef debug
0122 G4cout << "Fill Time slice histo " << j << " with " << timeDeposit[j] << G4endl;
0123 totalFilledTimeProfile += timeDeposit[j];
0124 #endif
0125
0126 G4double t = j + 0.5;
0127 man->FillH1(IDTimeProfile+1,t,timeDeposit[j]);
0128 #ifdef debug
0129 G4cout << "Fill Time profile histo 1 with " << t << " " << timeDeposit[j] << G4endl;
0130 #endif
0131 }
0132 #ifdef debug
0133 G4cout << "CCalAnalysis::InsertTime: Total filled Time profile histo "
0134 << totalFilledTimeProfile << G4endl;
0135 #endif
0136
0137 for (G4int i=0; i<70; i++){LateralProfile[i] = 0.f;}
0138 for (G4int i=0; i<200; i++){timeDeposit[i] = 0.f;}
0139
0140 G4cout << " --- End of event --- " << G4endl;
0141
0142 }