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