Belle II Software development
SensitiveDetector.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 <beast/dosi/simulation/SensitiveDetector.h>
10#include <beast/dosi/dataobjects/DosiSimHit.h>
11
12#include <framework/datastore/StoreArray.h>
13#include <framework/datastore/RelationArray.h>
14
15#include <G4Track.hh>
16#include <G4Step.hh>
17
18namespace Belle2 {
24 namespace dosi {
25
27 Simulation::SensitiveDetectorBase("DosiSensitiveDetector", Const::invalidDetector)
28 {
29 m_hitNum = 0;
30 m_EvnetNumber = 0;
32 m_trackID = 0;
33 m_startTime = 0;
34 m_endTime = 0;
35 m_WightedTime = 0;
36 m_startEnergy = 0;
38 m_trackLength = 0;
39 iECLCell = 0;
40 TimeIndex = 0;
41 local_pos = 0;
42 T_ave = 0;
43 firstcall = 0;
44 m_phiID = 0;
45 m_thetaID = 0;
46 m_cellID = 0;
47
48 //Make sure all collections are registered
49 StoreArray<MCParticle> mcParticles;
51 RelationArray relMCSimHit(mcParticles, simHits);
52
53 //Register all collections we want to modify and require those we want to use
54 mcParticles.registerInDataStore();
55 simHits.registerInDataStore();
56 relMCSimHit.registerInDataStore();
57
58 //Register the Relation so that the TrackIDs get replaced by the actual
59 //MCParticle indices after simulating the events. This is needed as
60 //secondary particles might not be stored so everything relating to those
61 //particles will be attributed to the last saved mother particle
62 registerMCParticleRelation(relMCSimHit);
63 }
64
66 {
67
68 }
69
70 bool SensitiveDetector::step(G4Step* step, G4TouchableHistory*)
71 {
72 const G4StepPoint& preStep = *step->GetPreStepPoint();
73 const G4StepPoint& postStep = *step->GetPostStepPoint();
74
75 G4Track& track = *step->GetTrack();
76 if (m_trackID != track.GetTrackID()) {
77 //TrackID changed, store track informations
78 m_trackID = track.GetTrackID();
79 //Get momentum
80 m_momentum = preStep.GetMomentum() ;
81 //Get energy
82 m_startEnergy = preStep.GetKineticEnergy() ;
83 //Reset energy deposit;
85 //Reset Wighted Time;
86 m_WightedTime = 0;
87 //Reset m_WightedPos;
88 m_WightedPos.SetXYZ(0, 0, 0);
89
90 }
91 //Update energy deposit
92 m_energyDeposit += step->GetTotalEnergyDeposit() ;
93
94 m_startTime = preStep.GetGlobalTime();
95 m_endTime = postStep.GetGlobalTime();
96 m_WightedTime += (m_startTime + m_endTime) / 2 * (step->GetTotalEnergyDeposit());
97
98 m_startPos = preStep.GetPosition();
99 m_endPos = postStep.GetPosition();
100 ROOT::Math::XYZVector position((m_startPos.getX() + m_endPos.getX()) / 2 / CLHEP::cm,
101 (m_startPos.getY() + m_endPos.getY()) / 2 / CLHEP::cm,
102 (m_startPos.getZ() + m_endPos.getZ()) / 2 / CLHEP::cm);
103 m_WightedPos += position * (step->GetTotalEnergyDeposit());
104
105 //Save Hit if track leaves volume or is killed
106 if (track.GetNextVolume() != track.GetVolume() || track.GetTrackStatus() >= fStopAndKill) {
107 int pdgCode = track.GetDefinition()->GetPDGEncoding();
108
109 //const G4VPhysicalVolume& v = * track.GetVolume();
110 //G4ThreeVector posCell = v.GetTranslation();
111 // Get layer ID
112
113 //if (v.GetName().find("Crystal") != std::string::npos) {
114 //CsiGeometryPar* eclp = CsiGeometryPar::Instance();
115 m_cellID = step->GetTrack()->GetVolume()->GetCopyNo();
116
117 double dTotalEnergy = 1 / m_energyDeposit; //avoid the error no match for 'operator/'
118 if (m_energyDeposit > 0.) {
119
121 m_energyDeposit, m_momentum, m_WightedPos * dTotalEnergy);
122 //}
123 }
124
125 //Reset TrackID
126 m_trackID = 0;
127 }
128
129 return true;
130 }
131 /*
132 //Get Track information
133 const G4Track& track = *step->GetTrack();
134 const int trackID = track.GetTrackID();
135 const double depEnergy = step->GetTotalEnergyDeposit() * CLHEP::MeV;
136 const double nielEnergy = step->GetNonIonizingEnergyDeposit() * CLHEP::MeV;
137 const G4ThreeVector G4tkPos = step->GetTrack()->GetPosition();
138 float tkPos[3];
139 tkPos[0] = G4tkPos.x() * CLHEP::cm;
140 tkPos[1] = G4tkPos.y() * CLHEP::cm;
141 tkPos[2] = G4tkPos.z() * CLHEP::cm;
142 const G4ThreeVector G4tkMom = step->GetTrack()->GetMomentum();
143 float tkMom[3];
144 tkMom[0] = G4tkMom.x() * CLHEP::MeV;
145 tkMom[1] = G4tkMom.y() * CLHEP::MeV;
146 tkMom[2] = G4tkMom.z() * CLHEP::MeV;
147 const G4ThreeVector G4tkMomDir = step->GetTrack()->GetMomentumDirection();
148 float tkMomDir[3];
149 tkMomDir[0] = G4tkMomDir.x() * CLHEP::MeV;
150 tkMomDir[1] = G4tkMomDir.y() * CLHEP::MeV;
151 tkMomDir[2] = G4tkMomDir.z() * CLHEP::MeV;
152 const int tkPDG = step->GetTrack()->GetDefinition()->GetPDGEncoding();
153 const double tkKEnergy = step->GetTrack()->GetKineticEnergy();
154 const int detNb = step->GetTrack()->GetVolume()->GetCopyNo();
155 const double GlTime = step->GetPreStepPoint()->GetGlobalTime();
156 //Ignore everything below 1eV
157 if (depEnergy < CLHEP::eV) return false;
158
159 //Get the datastore arrays
160 StoreArray<MCParticle> mcParticles;
161 StoreArray<DosiSimHit> simHits;
162 RelationArray relMCSimHit(mcParticles, simHits);
163
164 StoreArray<DosiSimHit> DosiHits;
165 if (!DosiHits.isValid()) DosiHits.create();
166 DosiSimHit* hit = DosiHits.appendNew(
167 depEnergy,
168 nielEnergy,
169 tkPDG,
170 tkKEnergy,
171 detNb,
172 GlTime,
173 tkPos,
174 tkMom,
175 tkMomDir
176 );
177
178 //Add Relation between SimHit and MCParticle with a weight of 1. Since
179 //the MCParticle index is not yet defined we use the trackID from Geant4
180 relMCSimHit.add(trackID, hit->getArrayIndex(), 1.0);
181
182 return true;
183 }
184 */
185
187 const G4int cellId,
188 const G4int trackID,
189 const G4int pid,
190 const G4double tof,
191 const G4double edep,
192 G4ThreeVector mom,
193 ROOT::Math::XYZVector posAve)
194 {
195
196 //Get the datastore arraus
197 StoreArray<MCParticle> mcParticles;
199 RelationArray relMCSimHit(mcParticles, simHits);
200
201 ROOT::Math::XYZVector momentum(mom.getX() / CLHEP::GeV, mom.getY() / CLHEP::GeV, mom.getZ() / CLHEP::GeV);
202 simHits.appendNew(cellId, trackID, pid, tof / CLHEP::ns, edep / CLHEP::GeV, momentum, posAve);
203 int simhitNumber = simHits.getEntries() - 1;
204 B2DEBUG(150, "HitNumber: " << simhitNumber);
205 relMCSimHit.add(trackID, simhitNumber, 1.0);
206 return (simhitNumber);
207 }//saveSimHit
208
209
210 } //dosi namespace
212} //Belle2 namespace
This class provides a set of constants for the framework.
Definition: Const.h:34
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:62
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.
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
int m_thetaID
The current phi ID in an event.
double m_energyDeposit
particle energy at the entrance in volume
int m_EvnetNumber
The current number of created hits in an event.
int saveSimHit(const G4int cellId, const G4int trackID, const G4int pid, const G4double tof, const G4double edep, G4ThreeVector mom, ROOT::Math::XYZVector WightedPos)
Save DosiSimHit into datastore.
int TimeIndex
Hit Energy of StoreArray.
int m_trackID
The current number of created hits in an event.
int m_hitNum
members of SensitiveDetector
int firstcall
flight time to diode sensor
double T_ave
position alongthe vector of crystal axis
double local_pos
vector of crystal axis
G4ThreeVector m_momentum
Wighted step Position.
int m_cellID
The current theta ID in an event.
double m_startTime
momentum of track
int iECLCell
length of the track in the volume
double m_trackLength
energy deposited in volume
int m_oldEvnetNumber
The current number of created hits in an event.
G4ThreeVector m_startPos
track id
ROOT::Math::XYZVector m_WightedPos
Position of poststep.
G4ThreeVector m_endPos
Position of prestep.
bool step(G4Step *step, G4TouchableHistory *) override
Step processing method.
Abstract base class for different kinds of events.