10 #include <ecl/simulation/SensitiveDetector.h>
11 #include <ecl/geometry/ECLGeometryPar.h>
12 #include <framework/gearbox/Unit.h>
13 #include <simulation/background/BkgNeutronWeight.h>
16 using namespace Belle2::ECL;
20 SensitiveDetector::SensitiveDetector(G4String name, G4double UNUSED(thresholdEnergyDeposit),
21 G4double UNUSED(thresholdKineticEnergy)):
22 Simulation::SensitiveDetectorBase(name, Const::ECL),
23 m_eclSimHitRel(m_mcParticles, m_eclSimHits),
24 m_eclHitRel(m_mcParticles, m_eclHits),
25 m_ECLHadronComponentEmissionFunction()
32 m_hadronenergyDeposit = 0;
37 m_eclSimHits.registerInDataStore();
38 m_eclHits.registerInDataStore();
40 m_mcParticles.registerRelationTo(m_eclSimHits);
41 m_mcParticles.registerRelationTo(m_eclHits);
43 m_HadronEmissionFunction = m_ECLHadronComponentEmissionFunction->getHadronComponentEmissionFunction();
46 SensitiveDetector::~SensitiveDetector()
50 void SensitiveDetector::Initialize(G4HCofThisEvent*)
56 double SensitiveDetector::GetHadronIntensityFromDEDX(
double x)
59 if (x > 232)
return m_HadronEmissionFunction->Eval(232);
60 return m_HadronEmissionFunction->Eval(x);
65 double GetCsITlScintillationEfficiency(
double x)
67 const double p0 = 1.52;
68 const double p1 = 0.00344839;
69 const double p2 = 2.0;
72 val = (x * p0) / (x + (p1 * x * x) + p2);
82 double preKineticEnergy = aStep->GetPreStepPoint()->GetKineticEnergy();
83 double postKineticEnergy = aStep->GetPostStepPoint()->GetKineticEnergy();
84 double avgKineticEnergy = 0.5 * (preKineticEnergy + postKineticEnergy);
85 const G4ParticleDefinition* StepParticleDefinition = aStep->GetTrack()->GetParticleDefinition();
86 G4Material* StepMaterial = aStep->GetTrack()->GetMaterial();
87 const double CsIDensity = 4.51;
88 double ELE_DEDX = m_emCal.ComputeDEDX(avgKineticEnergy, StepParticleDefinition,
"eIoni" ,
89 StepMaterial) / CLHEP::MeV * CLHEP::cm / (CsIDensity);
90 double MU_DEDX = m_emCal.ComputeDEDX(avgKineticEnergy, StepParticleDefinition,
"muIoni" ,
91 StepMaterial) / CLHEP::MeV * CLHEP::cm / (CsIDensity);
92 double HAD_DEDX = m_emCal.ComputeDEDX(avgKineticEnergy, StepParticleDefinition,
"hIoni" ,
93 StepMaterial) / CLHEP::MeV * CLHEP::cm / (CsIDensity);
94 double ION_DEDX = m_emCal.ComputeDEDX(avgKineticEnergy, StepParticleDefinition,
"ionIoni",
95 StepMaterial) / CLHEP::MeV * CLHEP::cm / (CsIDensity);
96 G4double DEDX_val = ELE_DEDX + MU_DEDX + HAD_DEDX + ION_DEDX;
98 G4double edep = aStep->GetTotalEnergyDeposit();
99 double LightOutputCorrection = GetCsITlScintillationEfficiency(DEDX_val);
100 edep *= LightOutputCorrection;
101 const G4StepPoint& s0 = *aStep->GetPreStepPoint();
102 const G4StepPoint& s1 = *aStep->GetPostStepPoint();
103 const G4Track& track = *aStep->GetTrack();
105 if (m_trackID != track.GetTrackID()) {
106 m_trackID = track.GetTrackID();
107 m_momentum = s0.GetMomentum();
111 m_WeightedPos.set(0, 0, 0);
112 m_hadronenergyDeposit = 0;
115 G4double hedep = 0.5 * edep;
116 m_energyDeposit += edep;
117 m_WeightedTime += (s0.GetGlobalTime() + s1.GetGlobalTime()) * hedep;
118 m_WeightedPos += (s0.GetPosition() + s1.GetPosition()) * hedep;
121 m_hadronenergyDeposit += (edep / CLHEP::GeV) * GetHadronIntensityFromDEDX(DEDX_val);
124 if (track.GetNextVolume() != track.GetVolume() || track.GetTrackStatus() >= fStopAndKill) {
125 if (m_energyDeposit > 0.0) {
128 G4double dTotalEnergy = 1 / m_energyDeposit;
129 m_WeightedTime *= dTotalEnergy;
130 m_WeightedPos *= dTotalEnergy;
131 int pdgCode = track.GetDefinition()->GetPDGEncoding();
132 saveSimHit(cellID, m_trackID, pdgCode, m_WeightedTime, m_energyDeposit, m_momentum, m_WeightedPos, m_hadronenergyDeposit);
140 void SensitiveDetector::EndOfEvent(G4HCofThisEvent*)
142 struct hit_t {
double e, t; };
145 struct thit_t {
int tid, hid; };
150 for (
const ECLSimHit& t : m_eclSimHits) {
151 double tof = t.getFlightTime();
152 if (tof >= 8000.0 || tof < -8000.0)
continue;
153 int TimeIndex = (tof + 8000.0) * (1. / 100);
154 int cellId = t.getCellId() - 1;
155 int key = cellId * 160 + TimeIndex;
156 double edep = t.getEnergyDep();
157 double tsen = tof + eclp->time2sensor(cellId, t.getPosition());
161 int trkid = t.getTrackId();
162 tr.push_back({trkid, key});
164 double old_edep = h.e, old_tsen = h.t;
165 double new_edep = old_edep + edep;
168 h.t = (old_edep * old_tsen + edep * tsen) / new_edep;
171 int hitNum = m_eclHits.getEntries();
173 for (
const pair<int, hit_t>& t : a) {
174 int key = t.first, cellId = key / 160;
175 for (
const thit_t& s : tr)
176 if (s.hid == key) m_eclHitRel.add(s.tid, hitNum);
177 const hit_t& h = t.second;
178 m_eclHits.appendNew(cellId + 1, h.e, h.t); hitNum++;
182 int SensitiveDetector::saveSimHit(G4int cellId, G4int trackID, G4int pid, G4double tof, G4double edep,
183 const G4ThreeVector& mom,
const G4ThreeVector& pos,
double Hadronedep)
185 int simhitNumber = m_eclSimHits.getEntries();
186 m_eclSimHitRel.add(trackID, simhitNumber);
187 tof *= 1 / CLHEP::ns;
188 edep *= 1 / CLHEP::GeV;
189 m_eclSimHits.appendNew(cellId + 1, trackID, pid, tof, edep, mom * (1 / CLHEP::GeV), pos * (1 / CLHEP::cm), Hadronedep);
194 Simulation::SensitiveDetectorBase(name,
Const::ECL)
196 m_eclHits.registerInDataStore(
"ECLDiodeHits");
218 const G4StepPoint& s0 = *aStep->GetPreStepPoint();
219 const G4StepPoint& s1 = *aStep->GetPostStepPoint();
220 const G4Track& track = *aStep->GetTrack();
229 G4double edep = aStep->GetTotalEnergyDeposit();
231 m_tsum += (s0.GetGlobalTime() + s1.GetGlobalTime()) * edep;
233 if (track.GetNextVolume() != track.GetVolume() ||
234 track.GetTrackStatus() >= fStopAndKill) {
242 const G4ThreeVector& p = s0.GetPosition();
245 double dz = (p - cp) * cv;
246 double rho = ((p - cp) - dz * cv).perp();
247 cout <<
"diff " << cellID <<
" " << dz <<
" " << rho << endl;
248 if (abs(dz - 150) > 1) {
249 const G4VTouchable* touch = s0.GetTouchable();
250 const G4NavigationHistory* h = touch->GetHistory();
251 int hd = h->GetDepth();
252 int i1 = h->GetReplicaNo(hd - 1);
253 int i2 = h->GetReplicaNo(hd - 2);
254 const G4String& vname = touch->GetVolume()->GetName();
255 cout << vname <<
" " << hd <<
" " << i1 <<
" " << i2 << endl;
256 for (
int i = 0; i <= hd; i++) {
257 G4VPhysicalVolume* v = h->GetVolume(i);
258 cout << v->GetName() << endl;
271 struct dep_t {
double e, t;};
276 if (cellId == h.cellId) v.push_back({h.e, h.t});
279 double esum = 0, tsum = 0;
280 for (
const dep_t& t : v) {
287 for (
const dep_t& t : v) trms += pow(t.t - tsum, 2) * t.e ;
290 tsum *= 1 / CLHEP::ns;
291 esum *= 1 / CLHEP::GeV;
292 m_eclHits.appendNew(cellId + 1, esum, tsum);
306 Simulation::SensitiveDetectorBase(name,
Const::invalidDetector),
307 m_eclBeamBkgHitRel(m_mcParticles, m_eclBeamBkgHits)
326 const G4StepPoint& s0 = *aStep->GetPreStepPoint();
327 const G4Track& track = *aStep->GetTrack();
332 const G4ThreeVector& worldPosition = s0.GetPosition();
334 m_startPos.SetXYZ(worldPosition.x() * mm2cm , worldPosition.y() * mm2cm, worldPosition.z() * mm2cm);
336 const G4ThreeVector& momentum = s0.GetMomentum() ;
351 if (track.GetNextVolume() != track.GetVolume() ||
352 track.GetTrackStatus() >= fStopAndKill) {
353 int pdgCode = track.GetDefinition()->GetPDGEncoding();
354 double endEnergy = track.GetKineticEnergy() *
Unit::MeV;
355 double neutWeight = 0;
356 if (pdgCode == 2112) {