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/beamabort/simulation/SensitiveDetector.h>
12 #include <beast/beamabort/dataobjects/BeamabortSimHit.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 beamabort {
27 
29  Simulation::SensitiveDetectorBase("BeamabortSensitiveDetector", 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;
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 
67  SensitiveDetector::~SensitiveDetector()
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  saveSimHit(m_cellID, m_trackID, pdgCode, m_WightedTime / m_energyDeposit,
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<BeamabortSimHit> simHits;
162  RelationArray relMCSimHit(mcParticles, simHits);
163 
164  StoreArray<BeamabortSimHit> BeamabortHits;
165  if (!BeamabortHits.isValid()) BeamabortHits.create();
166  BeamabortSimHit* hit = BeamabortHits.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 
186  int SensitiveDetector::saveSimHit(
187  const G4int cellId,
188  const G4int trackID,
189  const G4int pid,
190  const G4double tof,
191  const G4double edep,
192  G4ThreeVector mom,
193  TVector3 posAve)
194  {
195 
196  //Get the datastore arraus
197  StoreArray<MCParticle> mcParticles;
199  RelationArray relMBeamabortmHit(mcParticles, simHits);
200 
201  StoreArray<BeamabortSimHit> BeamabortHits;
202  RelationArray beamabortSimHitRel(mcParticles, BeamabortHits);
203  TVector3 momentum(mom.getX() / CLHEP::GeV, mom.getY() / CLHEP::GeV, mom.getZ() / CLHEP::GeV);
204  BeamabortHits.appendNew(cellId, trackID, pid, tof / CLHEP::ns, edep / CLHEP::GeV, momentum, posAve);
205  int simhitNumber = BeamabortHits.getEntries() - 1;
206  B2DEBUG(150, "HitNumber: " << simhitNumber);
207  beamabortSimHitRel.add(trackID, simhitNumber);
208  return (simhitNumber);
209  }//saveSimHit
210 
211 
212  } //beamabort namespace
214 } //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::beamabort::SensitiveDetector::SensitiveDetector
SensitiveDetector()
Constructor.
Definition: SensitiveDetector.cc:36
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::arich::SensitiveDetector::step
bool step(G4Step *aStep, G4TouchableHistory *) override
Process each step and calculate variables defined in PXDSimHit.
Definition: SensitiveDetector.cc:51
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
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::StoreArray< MCParticle >
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226