Belle II Software  release-05-01-25
SensitiveDetector.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Poyuan Chen *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 #include <ecl/simulation/SensitiveDetector.h>
11 #include <ecl/geometry/ECLGeometryPar.h>
12 #include <framework/gearbox/Unit.h>
13 #include <simulation/background/BkgNeutronWeight.h>
14 
15 using namespace std;
16 using namespace Belle2::ECL;
17 
18 #define UNUSED(x)
19 
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()
26  // m_thresholdEnergyDeposit(thresholdEnergyDeposit),
27  // m_thresholdKineticEnergy(thresholdKineticEnergy)
28 {
29  m_trackID = 0;
30  m_WeightedTime = 0;
31  m_energyDeposit = 0;
32  m_hadronenergyDeposit = 0;
33 
34  registerMCParticleRelation(m_eclSimHitRel);
35  registerMCParticleRelation(m_eclHitRel);
36 
37  m_eclSimHits.registerInDataStore();
38  m_eclHits.registerInDataStore();
39 
40  m_mcParticles.registerRelationTo(m_eclSimHits);
41  m_mcParticles.registerRelationTo(m_eclHits);
42 
43  m_HadronEmissionFunction = m_ECLHadronComponentEmissionFunction->getHadronComponentEmissionFunction();
44 }
45 
46 SensitiveDetector::~SensitiveDetector()
47 {
48 }
49 
50 void SensitiveDetector::Initialize(G4HCofThisEvent*)
51 {
52 }
53 
54 //Returns percent of scintillation emission from hadron component for given ionization dEdx value.
55 //See section 5 of S. Longo and J. M. Roney 2018 JINST 13 P03018 for additional details.
56 double SensitiveDetector::GetHadronIntensityFromDEDX(double x)
57 {
58  if (x < 2) return 0;
59  if (x > 232) return m_HadronEmissionFunction->Eval(232);
60  return m_HadronEmissionFunction->Eval(x);
61 }
62 
63 //Returns total scintillation efficiency, normalized to 1 for photons, for CsI(Tl) given ionization dEdx value.
64 //See section 5 of S. Longo and J. M. Roney 2018 JINST 13 P03018 for additional details.
65 double GetCsITlScintillationEfficiency(double x)
66 {
67  const double p0 = 1.52;
68  const double p1 = 0.00344839;
69  const double p2 = 2.0;
70  double val = 1;
71  if (x > 3.9406) {
72  val = (x * p0) / (x + (p1 * x * x) + p2);
73  }
74  return val;
75 }
76 //-----------------------------------------------------
77 // Method invoked for every step in sensitive detector
78 //-----------------------------------------------------
79 bool SensitiveDetector::step(G4Step* aStep, G4TouchableHistory*)
80 {
81  //
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; //gcm^-3
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; //Ionization dE/dx for any particle type
97 
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();
104 
105  if (m_trackID != track.GetTrackID()) { //TrackID changed, store track informations
106  m_trackID = track.GetTrackID();
107  m_momentum = s0.GetMomentum(); // Get momentum
108  //Reset track parameters
109  m_energyDeposit = 0;
110  m_WeightedTime = 0;
111  m_WeightedPos.set(0, 0, 0);
112  m_hadronenergyDeposit = 0;
113  }
114  //Update energy deposit
115  G4double hedep = 0.5 * edep; // half of energy deposition -- for averaging purpose
116  m_energyDeposit += edep;
117  m_WeightedTime += (s0.GetGlobalTime() + s1.GetGlobalTime()) * hedep;
118  m_WeightedPos += (s0.GetPosition() + s1.GetPosition()) * hedep;
119  //
120  //Calculate hadronic scintillation component intensity from dEdx
121  m_hadronenergyDeposit += (edep / CLHEP::GeV) * GetHadronIntensityFromDEDX(DEDX_val);
122  //
123  //Save Hit if track leaves volume or is killed
124  if (track.GetNextVolume() != track.GetVolume() || track.GetTrackStatus() >= fStopAndKill) {
125  if (m_energyDeposit > 0.0) {
127  G4int cellID = eclp->TouchableToCellID(s0.GetTouchable());
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);
133  }
134  //Reset TrackID
135  m_trackID = 0;
136  }
137  return true;
138 }
139 
140 void SensitiveDetector::EndOfEvent(G4HCofThisEvent*)
141 {
142  struct hit_t { double e, t; };
143  map<int, hit_t> a;
144 
145  struct thit_t { int tid, hid; };
146  vector<thit_t> tr;
147 
149 
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()); // flight time to diode sensor
158 
159  hit_t& h = a[key];
160 
161  int trkid = t.getTrackId();
162  tr.push_back({trkid, key});
163 
164  double old_edep = h.e, old_tsen = h.t;
165  double new_edep = old_edep + edep;
166 
167  h.e = new_edep;
168  h.t = (old_edep * old_tsen + edep * tsen) / new_edep;
169  }
170 
171  int hitNum = m_eclHits.getEntries();
172  // assert(hitNum == 0);
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++;
179  }
180 }
181 
182 int SensitiveDetector::saveSimHit(G4int cellId, G4int trackID, G4int pid, G4double tof, G4double edep,
183  const G4ThreeVector& mom, const G4ThreeVector& pos, double Hadronedep)
184 {
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);
190  return simhitNumber;
191 }
192 
193 SensitiveDiode::SensitiveDiode(const G4String& name):
194  Simulation::SensitiveDetectorBase(name, Const::ECL)
195 {
196  m_eclHits.registerInDataStore("ECLDiodeHits");
198  m_trackID = -1;
199  m_tsum = 0;
200  m_esum = 0;
201 }
202 
204 {
205 }
206 
207 void SensitiveDiode::Initialize(G4HCofThisEvent*)
208 {
209  m_hits.clear();
210  m_cells.clear();
211 }
212 
213 //-----------------------------------------------------
214 // Method invoked for every step in sensitive detector
215 //-----------------------------------------------------
216 bool SensitiveDiode::step(G4Step* aStep, G4TouchableHistory*)
217 {
218  const G4StepPoint& s0 = *aStep->GetPreStepPoint();
219  const G4StepPoint& s1 = *aStep->GetPostStepPoint();
220  const G4Track& track = *aStep->GetTrack();
221 
222  if (m_trackID != track.GetTrackID()) { //TrackID changed, store track informations
223  m_trackID = track.GetTrackID();
224  //Reset track parameters
225  m_esum = 0;
226  m_tsum = 0;
227  }
228  //Update energy deposit
229  G4double edep = aStep->GetTotalEnergyDeposit();
230  m_esum += edep;
231  m_tsum += (s0.GetGlobalTime() + s1.GetGlobalTime()) * edep;
232  //Save Hit if track leaves volume or is killed
233  if (track.GetNextVolume() != track.GetVolume() ||
234  track.GetTrackStatus() >= fStopAndKill) {
235  if (m_esum > 0.0) {
236  int cellID = m_eclp->TouchableDiodeToCellID(s0.GetTouchable());
237  m_tsum /= 2 * m_esum; // average
238  m_hits.push_back({cellID, m_esum, m_tsum});
239  auto it = find(m_cells.begin(), m_cells.end(), cellID);
240  if (it == m_cells.end())m_cells.push_back(cellID);
241 #if 0
242  const G4ThreeVector& p = s0.GetPosition();
243  G4ThreeVector cp = 10 * m_eclp->getCrystalPos(cellID);
244  G4ThreeVector cv = m_eclp->getCrystalVec(cellID);
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); // index of each volume is set at physical volume creation
253  int i2 = h->GetReplicaNo(hd - 2); // go up in volume hierarchy
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;
259  }
260  }
261 #endif
262  }
263  //Reset TrackID
264  m_trackID = 0;
265  }
266  return true;
267 }
268 
269 void SensitiveDiode::EndOfEvent(G4HCofThisEvent*)
270 {
271  struct dep_t {double e, t;};
272  vector<dep_t> v;
273  for (int cellId : m_cells) {
274  v.clear();
275  for (const hit_t& h : m_hits) {
276  if (cellId == h.cellId) v.push_back({h.e, h.t});
277  }
278 
279  double esum = 0, tsum = 0;
280  for (const dep_t& t : v) {
281  esum += t.e;
282  tsum += t.t * t.e;
283  }
284  tsum /= esum;
285 
286  double trms = 0;
287  for (const dep_t& t : v) trms += pow(t.t - tsum, 2) * t.e ;
288  trms /= esum;
289 
290  tsum *= 1 / CLHEP::ns;
291  esum *= 1 / CLHEP::GeV;
292  m_eclHits.appendNew(cellId + 1, esum, tsum);
293 
294 
295  // cout<<"DD: "<<cellId<<" "<<esum<<" "<<tsum<<" +- "<<sqrt(tsum2)<<endl;
296 
297  // sort(v.begin(),v.end(),[](const dep_t &a, const dep_t &b){return a.t<b.t;});
298  // cout<<"DD: "<<cellId<<" ";
299  // for(const dep_t &t: v)cout<<t.e<<" MeV "<<t.t<<" nsec, ";
300  // cout<<"\n";
301 
302  }
303 }
304 
306  Simulation::SensitiveDetectorBase(name, Const::invalidDetector),
307  m_eclBeamBkgHitRel(m_mcParticles, m_eclBeamBkgHits)
308 {
310  m_eclBeamBkgHits.registerInDataStore();
312 
314  m_energyDeposit = 0;
315  m_startEnergy = 0;
316  m_startTime = 0;
317  m_trackLength = 0;
318  m_trackID = -1;
319 }
320 
321 //-----------------------------------------------------
322 // Method invoked for every step in sensitive detector
323 //-----------------------------------------------------
324 bool BkgSensitiveDiode::step(G4Step* aStep, G4TouchableHistory*)
325 {
326  const G4StepPoint& s0 = *aStep->GetPreStepPoint();
327  const G4Track& track = *aStep->GetTrack();
328 
329  if (m_trackID != track.GetTrackID()) { //TrackID changed, store track informations
330  m_trackID = track.GetTrackID();
331  //Get world position
332  const G4ThreeVector& worldPosition = s0.GetPosition();
333  const double mm2cm = Unit::mm / Unit::cm;
334  m_startPos.SetXYZ(worldPosition.x() * mm2cm , worldPosition.y() * mm2cm, worldPosition.z() * mm2cm);
335  //Get momentum
336  const G4ThreeVector& momentum = s0.GetMomentum() ;
337  m_startMom.SetXYZ(momentum.x() * Unit::MeV, momentum.y() * Unit::MeV, momentum.z() * Unit::MeV);
338  //Get time
339  m_startTime = s0.GetGlobalTime();
340  //Get energy
341  m_startEnergy = s0.GetKineticEnergy() * Unit::MeV;
342  //Reset energy deposit;
343  m_energyDeposit = 0;
344  //Reset track lenght;
345  m_trackLength = 0;
346  }
347  //Update energy deposit
348  m_energyDeposit += aStep->GetTotalEnergyDeposit() * Unit::MeV;
349  m_trackLength += aStep->GetStepLength() * Unit::mm;
350  //Save Hit if track leaves volume or is killed
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) {
358  neutWeight = wt.getWeight(m_startEnergy / Unit::MeV);
359  }
360 
361  int bkgHitNumber = m_eclBeamBkgHits.getEntries();
362  m_eclBeamBkgHitRel.add(m_trackID, bkgHitNumber);
363 
364  int cellID = m_eclp->TouchableDiodeToCellID(s0.GetTouchable());
365  m_eclBeamBkgHits.appendNew(6, cellID, pdgCode, m_trackID, m_startPos, m_startMom, m_startTime, endEnergy, m_startEnergy,
366  m_energyDeposit, m_trackLength, neutWeight);
367 
368  //Reset TrackID
369  m_trackID = 0;
370  }
371  return true;
372 }
Belle2::Unit::cm
static const double cm
Standard units with the value = 1.
Definition: Unit.h:57
Belle2::ECL::SensitiveDiode::m_esum
double m_esum
total energy deposited in a volume by a track
Definition: SensitiveDetector.h:113
Belle2::ECL::BkgSensitiveDiode::m_eclBeamBkgHitRel
RelationArray m_eclBeamBkgHitRel
MCParticle to BeamBackHit relation array.
Definition: SensitiveDetector.h:139
Belle2::ECL::SensitiveDiode::EndOfEvent
void EndOfEvent(G4HCofThisEvent *eventHC) override
Do what you want to do at the end of each event.
Definition: SensitiveDetector.cc:269
Belle2::StoreArray::registerRelationTo
bool registerRelationTo(const StoreArray< TO > &toArray, DataStore::EDurability durability=DataStore::c_Event, DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut, const std::string &namedRelation="") const
Register a relation to the given StoreArray.
Definition: StoreArray.h:150
Belle2::ECL::BkgSensitiveDiode::m_energyDeposit
double m_energyDeposit
energy deposited in volume
Definition: SensitiveDetector.h:134
Belle2::ECL::ECLGeometryPar::getCrystalVec
G4ThreeVector getCrystalVec(int cid)
The direction of crystal.
Definition: ECLGeometryPar.h:95
Belle2::ECLSimHit
ClassECLSimHit - Geant4 simulated hit for the ECL.
Definition: ECLSimHit.h:42
Belle2::ECL::BkgSensitiveDiode::m_eclp
ECLGeometryPar * m_eclp
pointer to ECLGeometryPar
Definition: SensitiveDetector.h:136
Belle2::ECL::BkgSensitiveDiode::step
bool step(G4Step *aStep, G4TouchableHistory *history) override
Process each step and calculate variables defined in ECLHit.
Definition: SensitiveDetector.cc:324
Belle2::ECL::SensitiveDiode::step
bool step(G4Step *aStep, G4TouchableHistory *history) override
Process each step and calculate variables defined in ECLHit.
Definition: SensitiveDetector.cc:216
Belle2::ECL::SensitiveDiode::m_eclp
ECLGeometryPar * m_eclp
pointer to ECLGeometryPar
Definition: SensitiveDetector.h:104
Belle2::Simulation::SensitiveDetectorBase::registerMCParticleRelation
static void registerMCParticleRelation(const std::string &name, RelationArray::EConsolidationAction ignoreAction=RelationArray::c_negativeWeight)
Register an relation involving MCParticles.
Definition: SensitiveDetectorBase.cc:24
Belle2::ECL::SensitiveDiode::hit_t
simple hit structure
Definition: SensitiveDetector.h:106
Belle2::ECL::SensitiveDiode::m_eclHits
StoreArray< ECLHit > m_eclHits
ECLHit array.
Definition: SensitiveDetector.h:117
Belle2::ECL::BkgSensitiveDiode::m_trackLength
double m_trackLength
length of the track in the volume
Definition: SensitiveDetector.h:135
Belle2::ECL::ECLGeometryPar::TouchableToCellID
int TouchableToCellID(const G4VTouchable *)
The same as above but without sanity checks.
Definition: ECLGeometryPar.cc:398
Belle2::ECL::SensitiveDiode::~SensitiveDiode
~SensitiveDiode()
Destructor.
Definition: SensitiveDetector.cc:203
Belle2::ECL::ECLGeometryPar::Instance
static ECLGeometryPar * Instance()
Static method to get a reference to the ECLGeometryPar instance.
Definition: ECLGeometryPar.cc:151
Belle2::ECL::ECLGeometryPar::TouchableDiodeToCellID
int TouchableDiodeToCellID(const G4VTouchable *)
Mapping from G4VTouchable copyNumbers to Crystal CellID.
Definition: ECLGeometryPar.cc:393
Belle2::Unit::MeV
static const double MeV
[megaelectronvolt]
Definition: Unit.h:124
Belle2::ECL::BkgSensitiveDiode::m_startEnergy
double m_startEnergy
particle energy at the entrance in volume
Definition: SensitiveDetector.h:133
Belle2::ECL::BkgSensitiveDiode::m_eclBeamBkgHits
StoreArray< BeamBackHit > m_eclBeamBkgHits
BeamBackHit array.
Definition: SensitiveDetector.h:138
Belle2::srsensor::SensitiveDetector::step
bool step(G4Step *step, G4TouchableHistory *) override
Step processing method.
Definition: SensitiveDetector.cc:56
Belle2::ECL::SensitiveDiode::m_trackID
int m_trackID
current track id
Definition: SensitiveDetector.h:111
Belle2::BkgNeutronWeight::getInstance
static BkgNeutronWeight & getInstance()
Return a reference to the singleton BkgNeutronWeight instance.
Definition: BkgNeutronWeight.cc:30
Belle2::ECL::BkgSensitiveDiode::m_startPos
TVector3 m_startPos
particle position at the entrance in volume
Definition: SensitiveDetector.h:130
Belle2::ECL::SensitiveDiode::m_hits
std::vector< hit_t > m_hits
array of hits
Definition: SensitiveDetector.h:115
Belle2::RelationArray::add
void add(index_type from, index_type to, weight_type weight=1.0)
Add a new element to the relation.
Definition: RelationArray.h:291
Belle2::ECL::SensitiveDiode::m_tsum
double m_tsum
average track time weighted by energy deposition
Definition: SensitiveDetector.h:112
Belle2::ECL::SensitiveDiode::SensitiveDiode
SensitiveDiode(const G4String &)
Constructor.
Definition: SensitiveDetector.cc:193
Belle2::ECL::ECLGeometryPar::getCrystalPos
G4ThreeVector getCrystalPos(int cid)
The Position of crystal.
Definition: ECLGeometryPar.h:88
Belle2::ECL::BkgSensitiveDiode::BkgSensitiveDiode
BkgSensitiveDiode(const G4String &)
Constructor.
Definition: SensitiveDetector.cc:305
Belle2::Unit::mm
static const double mm
[millimeters]
Definition: Unit.h:80
Belle2::ECL::BkgSensitiveDiode::m_startMom
TVector3 m_startMom
particle momentum at the entrance in volume
Definition: SensitiveDetector.h:131
Belle2::ECL::BkgSensitiveDiode::m_mcParticles
StoreArray< MCParticle > m_mcParticles
MCParticle array.
Definition: SensitiveDetector.h:137
Belle2::ECL::SensitiveDiode::Initialize
void Initialize(G4HCofThisEvent *HCTE) override
Register ECL hits collection into G4HCofThisEvent.
Definition: SensitiveDetector.cc:207
Belle2::BkgNeutronWeight::getWeight
double getWeight(double ke)
Get weighting factor to convert a neutron to its 1-MeV equivalent.
Definition: BkgNeutronWeight.cc:36
Belle2::ECL::ECLGeometryPar
The Class for ECL Geometry Parameters.
Definition: ECLGeometryPar.h:45
Belle2::ECL::BkgSensitiveDiode::m_trackID
int m_trackID
track id
Definition: SensitiveDetector.h:129
Belle2::Const
This class provides a set of constants for the framework.
Definition: Const.h:36
Belle2::BkgNeutronWeight
The class to get the weighting factor for a 1-MeV-equivalent neutron flux on silicon.
Definition: BkgNeutronWeight.h:33
Belle2::ECL::BkgSensitiveDiode::m_startTime
double m_startTime
global time
Definition: SensitiveDetector.h:132
Belle2::ECL::SensitiveDiode::m_cells
std::vector< int > m_cells
array of hitted crystals
Definition: SensitiveDetector.h:116