Belle II Software development
BkgSensitiveDetector.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <simulation/background/BkgSensitiveDetector.h>
10#include <simulation/background/BkgNeutronWeight.h>
11#include <simulation/dataobjects/BeamBackHit.h>
12#include <framework/gearbox/Unit.h>
13#include <framework/gearbox/Const.h>
14
15#include <G4Track.hh>
16#include <G4Step.hh>
17
18#include <framework/datastore/StoreArray.h>
19#include <framework/datastore/RelationArray.h>
20
21using namespace std;
22
23namespace Belle2 {
28
29
30 BkgSensitiveDetector::BkgSensitiveDetector(const char* subDett, int iden):
31 Simulation::SensitiveDetectorBase("BKG", Const::invalidDetector), m_trackID(0), m_startPos(0., 0., 0.), m_startMom(0., 0., 0.),
33 {
34 // registration of store arrays and relations
35
36 StoreArray<MCParticle> mcParticles;
37 StoreArray<BeamBackHit> beamBackHits;
38
39 beamBackHits.registerInDataStore();
40 mcParticles.registerRelationTo(beamBackHits);
41
42 // additional registration of MCParticle relation
43
44 RelationArray relation(mcParticles, beamBackHits);
46
47 m_eclrepscale = 0;
48 std::string subDet = subDett;
49 if (subDet == "IR") m_subDet = 0;
50 else if (subDet == "PXD") m_subDet = 1;
51 else if (subDet == "SVD") m_subDet = 2;
52 else if (subDet == "CDC") m_subDet = 3;
53 else if (subDet == "ARICH") m_subDet = 4;
54 else if (subDet == "TOP") m_subDet = 5;
55 else if (subDet == "ECL") m_subDet = 6;
56 else if (subDet == "EKLM") m_subDet = 7;
57 else if (subDet == "BKLM") m_subDet = 8;
58 else m_subDet = 99;
59
60 if (m_subDet == 6) {
61 const int nc[] = {3, 4, 6, 9, 2, 9, 6, 4}; // the number of crystals in a sector per ring
62 const int indx[] = {96, 288, 864, 1151, 7776, 8064, 8544, 8736}; // regions with the same number of crystals in a ring
63 std::vector<int> indxv(indx, indx + sizeof(indx) / sizeof(int));
64 m_eclrepscale = nc[upper_bound(indxv.begin(), indxv.end(), iden) - indxv.begin()];
65 }
66
67 m_identifier = iden;
68
69 }
70
71
72 G4bool BkgSensitiveDetector::step(G4Step* aStep, G4TouchableHistory*)
73 {
74
75 const G4StepPoint& preStep = *aStep->GetPreStepPoint();
76 G4Track& track = *aStep->GetTrack();
77
78 if (m_trackID != track.GetTrackID()) {
79 //TrackID changed, store track informations
80 m_trackID = track.GetTrackID();
81 //Get world position
82 const G4ThreeVector& worldPosition = preStep.GetPosition();
83 m_startPos.SetXYZ(worldPosition.x() * Unit::mm / Unit::cm, worldPosition.y() * Unit::mm / Unit::cm,
84 worldPosition.z() * Unit::mm / Unit::cm);
85 //Get momentum
86 const G4ThreeVector& momentum = preStep.GetMomentum() ;
87 m_startMom.SetXYZ(momentum.x() * Unit::MeV, momentum.y() * Unit::MeV,
88 momentum.z() * Unit::MeV);
89 //Get time
90 m_startTime = preStep.GetGlobalTime();
91 //Get energy
92 m_startEnergy = preStep.GetKineticEnergy() * Unit::MeV;
93 //Reset energy deposit;
95 //Reset track lenght;
96 m_trackLength = 0;
97 }
98 //Update energy deposit
99 m_energyDeposit += aStep->GetTotalEnergyDeposit() * Unit::MeV;
100
101 m_trackLength += aStep->GetStepLength() * Unit::mm;
102
103 //Save Hit if track leaves volume or is killed
104 if (track.GetNextVolume() != track.GetVolume() || track.GetTrackStatus() >= fStopAndKill) {
105 int pdgCode = track.GetDefinition()->GetPDGEncoding();
106 double endEnergy = track.GetKineticEnergy() * Unit::MeV;
107 double neutWeight = 0;
108 if (pdgCode == Const::neutron.getPDGCode()) {
110 neutWeight = wt.getWeight(m_startEnergy / Unit::MeV);
111 }
112 StoreArray<BeamBackHit> beamBackHits;
113 int nentr = beamBackHits.getEntries();
114 int id = m_identifier;
115 if (m_subDet == 6) { // ECL
116 int sector = preStep.GetTouchable()->GetCopyNumber(1);
117 if ((sector == 0) && (m_eclrepscale == 2) && ((id & 1) != 0)) id += 144; // barrel odd-numbered diodes in sector 0 -> sector 144
118 id += sector * m_eclrepscale;
119 }
120 beamBackHits.appendNew(m_subDet, id, pdgCode, m_trackID, m_startPos, m_startMom, m_startTime, m_startEnergy, endEnergy,
121 m_energyDeposit, m_trackLength, neutWeight);
122
123
124 // create relation to MCParticle
125 StoreArray<MCParticle> mcParticles;
126 RelationArray relBeamBackHitToMCParticle(mcParticles, beamBackHits);
127 relBeamBackHitToMCParticle.add(m_trackID, nentr);
128
129 //Reset TrackID
130 m_trackID = 0;
131 }
132
133 return true;
134 }
135
136
138} // end of namespace Belle2
The class to get the weighting factor for a 1-MeV-equivalent neutron flux on silicon.
double m_energyDeposit
Energy deposited in volume.
ROOT::Math::XYZVector m_startPos
Particle position at the entrance in volume.
double m_startEnergy
Particle energy at the entrance in volume.
double m_trackLength
Length of the track in the volume.
int m_eclrepscale
Replica (=sector) scale in ECL.
int m_subDet
Subdetector id number.
ROOT::Math::XYZVector m_startMom
Particle momentum at the entrance in volume.
int m_identifier
Identifier of subdetector component.
This class provides a set of constants for the framework.
Definition Const.h:34
static const ParticleType neutron
neutron particle
Definition Const.h:675
Low-level class to create/modify relations between StoreArrays.
void add(index_type from, index_type to, weight_type weight=1.0)
Add a new element to the relation.
static void registerMCParticleRelation(const std::string &name, RelationArray::EConsolidationAction ignoreAction=RelationArray::c_negativeWeight)
Register an relation involving MCParticles.
SensitiveDetectorBase(const std::string &name, Const::EDetector subdetector)
Create a new Sensitive detecor with a given name and belonging to a given subdetector.
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
Accessor to arrays stored in the data store.
Definition StoreArray.h:113
T * appendNew()
Construct a new T object at the end of the array.
Definition StoreArray.h:246
int getEntries() const
Get the number of objects in the array.
Definition StoreArray.h:216
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:140
static const double mm
[millimeters]
Definition Unit.h:70
static const double MeV
[megaelectronvolt]
Definition Unit.h:114
static const double cm
Standard units with the value = 1.
Definition Unit.h:47
double getWeight(double ke)
Get weighting factor to convert a neutron to its 1-MeV equivalent.
static BkgNeutronWeight & getInstance()
Return a reference to the singleton BkgNeutronWeight instance.
BkgSensitiveDetector(const char *subDet, int iden=0)
Constructor.
bool step(G4Step *aStep, G4TouchableHistory *) override
Process each step and calculate variables defined in PXDSimHit.
Abstract base class for different kinds of events.
STL namespace.