Belle II Software  release-05-01-25
SensitiveBar.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2010 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Marko Petric, Marko Staric *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <top/simulation/SensitiveBar.h>
12 
13 #include <simulation/kernel/UserInfo.h>
14 #include <G4Step.hh>
15 #include <G4Track.hh>
16 #include <G4UnitsTable.hh>
17 #include <G4ParticleDefinition.hh>
18 #include <G4ParticleTypes.hh>
19 
20 #include <framework/logging/Logger.h>
21 #include <framework/gearbox/Unit.h>
22 
23 #include <TVector3.h>
24 #include <TRandom.h>
25 #include <cmath>
26 
27 using namespace std;
28 
29 namespace Belle2 {
34  namespace TOP {
35 
36  SensitiveBar::SensitiveBar():
37  Simulation::SensitiveDetectorBase("TOP", Const::TOP)
38  {
39 
40  // registration of store arrays and relations
41 
42  m_barHits.registerInDataStore();
44 
45  // additional registration of MCParticle relation (required for correct relations)
46 
48 
49  const auto* geo = m_topgp->getGeometry();
50  m_trackIDs.resize(geo->getNumModules(), 0);
51 
52  }
53 
54 
55  bool SensitiveBar::step(G4Step* aStep, G4TouchableHistory*)
56  {
57 
58  // get track and particle definition
59 
60  G4Track* aTrack = aStep->GetTrack();
61  G4ParticleDefinition* particle = aTrack->GetDefinition();
62 
63  // if optical photon, apply QE and return false
64 
65  if (particle == G4OpticalPhoton::OpticalPhotonDefinition()) {
66  auto* info = dynamic_cast<Simulation::TrackInfo*>(aTrack->GetUserInformation());
67  if (!info) return false;
68  if (info->getStatus() < 2) {
69  double energy = aTrack->GetKineticEnergy() * Unit::MeV / Unit::eV;
70  double qeffi = m_topgp->getPMTEfficiencyEnvelope(energy);
71  double fraction = info->getFraction();
72  if (qeffi == 0 or gRandom->Uniform() * fraction > qeffi) {
73  aTrack->SetTrackStatus(fStopAndKill);
74  return false;
75  }
76  info->setStatus(2);
77  info->setFraction(qeffi);
78  }
79  return false;
80  }
81 
82  // continue for other particles, if they enter the volume for the first time
83 
84  G4StepPoint* PrePosition = aStep->GetPreStepPoint();
85  if (PrePosition->GetStepStatus() != fGeomBoundary) return false; // not on boundary
86 
87  if (m_barHits.getEntries() == 0) {
88  for (auto& trackID : m_trackIDs) trackID = -1; // reset on new event
89  }
90 
91  int moduleID = PrePosition->GetTouchableHandle()->GetReplicaNumber(m_replicaDepth);
92  const auto* geo = m_topgp->getGeometry();
93  if (!geo->isModuleIDValid(moduleID)) {
94  B2ERROR("TOP::SensitiveBar: undefined module ID."
95  << LogVar("moduleID", moduleID));
96  return false;
97  }
98 
99  int trackID = aTrack->GetTrackID();
100  if (trackID == m_trackIDs[moduleID - 1]) return false; // not first time on boundary
101  m_trackIDs[moduleID - 1] = trackID;
102 
103  // particle other than optical photon entered the volume for the first time
104  // convert to basf2 units and write-out the hit
105 
106  G4ThreeVector worldPosition = PrePosition->GetPosition();
107  double tracklength = aTrack->GetTrackLength() - aStep->GetStepLength();
108  double globalTime = PrePosition->GetGlobalTime();
109  G4ThreeVector momentum = PrePosition->GetMomentum();
110 
111  TVector3 TPosition(worldPosition.x(), worldPosition.y(), worldPosition.z());
112  TVector3 TMomentum(momentum.x(), momentum.y(), momentum.z());
113  TVector3 TOrigin(aTrack->GetVertexPosition().x(),
114  aTrack->GetVertexPosition().y(),
115  aTrack->GetVertexPosition().z());
116 
117  // convert to Basf2 units
118  TPosition = TPosition * Unit::mm;
119  TMomentum = TMomentum * Unit::MeV;
120  TOrigin = TOrigin * Unit::mm;
121  tracklength = tracklength * Unit::mm;
122 
123  const auto& module = geo->getModule(moduleID);
124  TVector3 locPosition = module.pointToLocal(TPosition);
125  TVector3 locMomentum = module.momentumToLocal(TMomentum);
126  double theta = locMomentum.Theta();
127  double phi = locMomentum.Phi();
128 
129  int PDG = (int)(particle->GetPDGEncoding());
130 
131  // write hit to datastore
132  auto* hit = m_barHits.appendNew(moduleID, PDG, TOrigin, TPosition, TMomentum,
133  globalTime, tracklength, locPosition,
134  theta, phi);
135 
136  // set the relation
137  m_relParticleHit.add(trackID, hit->getArrayIndex());
138 
139  return true;
140  }
141 
142 
143  } // end of namespace top
145 } // end of namespace Belle2
Belle2::TOP::SensitiveBar::m_barHits
StoreArray< TOPBarHit > m_barHits
collection of entrance-to-bar hits
Definition: SensitiveBar.h:70
Belle2::TOP::SensitiveBar::m_replicaDepth
int m_replicaDepth
replica depth of module volume
Definition: SensitiveBar.h:65
Belle2::StoreArray::registerRelationTo
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:150
Belle2::TOP::SensitiveBar::m_mcParticles
StoreArray< MCParticle > m_mcParticles
collection of MC particles
Definition: SensitiveBar.h:69
Belle2::TOP::SensitiveBar::m_trackIDs
std::vector< int > m_trackIDs
track ID's
Definition: SensitiveBar.h:67
Belle2::Simulation::UserInfo
UserInfo class which is used to attach additional information to Geant4 particles and tracks.
Definition: UserInfo.h:46
Belle2::TOP::TOPGeometryPar::getPMTEfficiencyEnvelope
double getPMTEfficiencyEnvelope(double energy) const
Returns PMT efficiency envelope, e.g.
Definition: TOPGeometryPar.cc:179
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::TOP::SensitiveBar::m_relParticleHit
RelationArray m_relParticleHit
relations
Definition: SensitiveBar.h:71
Belle2::Unit::MeV
static const double MeV
[megaelectronvolt]
Definition: Unit.h:124
Belle2::Unit::eV
static const double eV
[electronvolt]
Definition: Unit.h:122
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::RelationArray::c_deleteElement
@ c_deleteElement
Delete the whole relation element if the original element got re-attributed.
Definition: RelationArray.h:91
LogVar
Class to store variables with their name which were sent to the logging service.
Definition: LogVariableStream.h:24
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::TOP::SensitiveBar::step
bool step(G4Step *aStep, G4TouchableHistory *) override
Process each step and fill TOPBarHits.
Definition: SensitiveBar.cc:55
Belle2::Unit::mm
static const double mm
[millimeters]
Definition: Unit.h:80
Belle2::TOP::SensitiveBar::m_topgp
TOPGeometryPar * m_topgp
geometry parameters
Definition: SensitiveBar.h:66
Belle2::TOP::TOPGeometryPar::getGeometry
const TOPGeometry * getGeometry() const
Returns pointer to geometry object using Basf2 units.
Definition: TOPGeometryPar.cc:167
Belle2::Const
This class provides a set of constants for the framework.
Definition: Const.h:36