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/dosi/simulation/SensitiveDetector.h>
12 #include <beast/dosi/dataobjects/DosiSimHit.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 dosi {
27 
29  Simulation::SensitiveDetectorBase("DosiSensitiveDetector", 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<DosiSimHit> 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<DosiSimHit> simHits;
163  RelationArray relMCSimHit(mcParticles, simHits);
164 
165  StoreArray<DosiSimHit> DosiHits;
166  if (!DosiHits.isValid()) DosiHits.create();
167  DosiSimHit* hit = DosiHits.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<DosiSimHit> simHits;
200  RelationArray relMCSimHit(mcParticles, simHits);
201 
202  TVector3 momentum(mom.getX() / CLHEP::GeV, mom.getY() / CLHEP::GeV, mom.getZ() / CLHEP::GeV);
203  simHits.appendNew(cellId, trackID, pid, tof / CLHEP::ns, edep / CLHEP::GeV, momentum, posAve);
204  int simhitNumber = simHits.getEntries() - 1;
205  B2DEBUG(150, "HitNumber: " << simhitNumber);
206  relMCSimHit.add(trackID, simhitNumber, 1.0);
207  return (simhitNumber);
208  }//saveSimHit
209 
210 
211  } //dosi namespace
213 } //Belle2 namespace
Belle2::StoreArray::appendNew
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:256
Belle2::RelationArray
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:72
Belle2::csi::SensitiveDetector::saveSimHit
int saveSimHit(const G4int cellId, const G4int trackID, const G4int pid, const G4double tof, const G4double edep, G4ThreeVector mom, G4ThreeVector pos, G4ThreeVector PosCell, G4ThreeVector VecCell)
Save ECLSimHits of the event into datastore.
Definition: SensitiveDetector.cc:165
Belle2::csi::SensitiveDetector::m_EvnetNumber
int m_EvnetNumber
The current number of created hits in an event.
Definition: SensitiveDetector.h:77
Belle2::csi::SensitiveDetector::m_trackLength
double m_trackLength
energy deposited in volume
Definition: SensitiveDetector.h:89
Belle2::csi::SensitiveDetector::step
bool step(G4Step *aStep, G4TouchableHistory *) override
Step processing method.
Definition: SensitiveDetector.cc:97
Belle2::dosi::SensitiveDetector::SensitiveDetector
SensitiveDetector()
Constructor.
Definition: SensitiveDetector.cc:36
Belle2::csi::SensitiveDetector::m_WightedTime
double m_WightedTime
global time
Definition: SensitiveDetector.h:86
Belle2::csi::SensitiveDetector::firstcall
int firstcall
flight time to diode sensor
Definition: SensitiveDetector.h:98
Belle2::csi::SensitiveDetector::m_trackID
int m_trackID
The current number of created hits in an event.
Definition: SensitiveDetector.h:79
Belle2::csi::SensitiveDetector::m_startEnergy
double m_startEnergy
global time
Definition: SensitiveDetector.h:87
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::csi::SensitiveDetector::m_startTime
double m_startTime
momentum of track
Definition: SensitiveDetector.h:84
Belle2::csi::SensitiveDetector::m_cellID
int m_cellID
The current theta ID in an event.
Definition: SensitiveDetector.h:101
Belle2::csi::SensitiveDetector::iECLCell
int iECLCell
Hit index of StoreArray.
Definition: SensitiveDetector.h:92
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::csi::SensitiveDetector::m_momentum
G4ThreeVector m_momentum
Wighted step Position.
Definition: SensitiveDetector.h:83
Belle2::csi::SensitiveDetector::m_startPos
G4ThreeVector m_startPos
track id
Definition: SensitiveDetector.h:80
Belle2::csi::SensitiveDetector::m_endPos
G4ThreeVector m_endPos
Position of prestep.
Definition: SensitiveDetector.h:81
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::csi::SensitiveDetector::m_thetaID
int m_thetaID
The current phi ID in an event.
Definition: SensitiveDetector.h:100
Belle2::csi::SensitiveDetector::m_oldEvnetNumber
int m_oldEvnetNumber
The current number of created hits in an event.
Definition: SensitiveDetector.h:78
Belle2::csi::SensitiveDetector::T_ave
double T_ave
position alongthe vector of crystal axis
Definition: SensitiveDetector.h:97
Belle2::csi::SensitiveDetector::~SensitiveDetector
~SensitiveDetector()
Destructor.
Definition: SensitiveDetector.cc:90
Belle2::csi::SensitiveDetector::local_pos
double local_pos
Hit Time of StoreArray.
Definition: SensitiveDetector.h:96
Belle2::StoreArray< MCParticle >
Belle2::csi::SensitiveDetector::TimeIndex
int TimeIndex
Hit Energy of StoreArray.
Definition: SensitiveDetector.h:93
Belle2::csi::SensitiveDetector::m_phiID
int m_phiID
flag of first call
Definition: SensitiveDetector.h:99
Belle2::csi::SensitiveDetector::m_hitNum
int m_hitNum
members of SensitiveDetector
Definition: SensitiveDetector.h:76
Belle2::csi::SensitiveDetector::m_energyDeposit
double m_energyDeposit
particle energy at the entrance in volume
Definition: SensitiveDetector.h:88
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226
Belle2::csi::SensitiveDetector::m_WightedPos
G4ThreeVector m_WightedPos
Position of poststep.
Definition: SensitiveDetector.h:82
Belle2::csi::SensitiveDetector::m_endTime
double m_endTime
global time
Definition: SensitiveDetector.h:85