Belle II Software  release-05-01-25
SensitivePMT.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/SensitivePMT.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 
26 using namespace std;
27 
28 namespace Belle2 {
33  namespace TOP {
34 
35  SensitivePMT::SensitivePMT():
36  Simulation::SensitiveDetectorBase("TOP", Const::TOP)
37  {
38 
39  m_simHits.registerInDataStore();
41 
42  m_simPhotons.registerInDataStore(DataStore::c_DontWriteOut);
43  m_simHits.registerRelationTo(m_simPhotons, DataStore::c_Event,
45 
47 
48  }
49 
50 
51  G4bool SensitivePMT::step(G4Step* aStep, G4TouchableHistory*)
52  {
53  // photon track
54  G4Track& photon = *aStep->GetTrack();
55 
56  // check if the track is an optical photon
57  if (photon.GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) return false;
58 
59  // photon energy in [eV]
60  double energy = photon.GetKineticEnergy() * Unit::MeV / Unit::eV;
61 
62  // pmt and module ID
63  int pmtID = photon.GetTouchableHandle()->GetReplicaNumber(m_pmtReplicaDepth);
64  int moduleID = photon.GetTouchableHandle()->GetReplicaNumber(m_moduleReplicaDepth);
65 
66  // hit position in local frame, converted to Basf units
67  G4ThreeVector localPosition = photon.GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(photon.GetPosition());
68  double xLocal = localPosition.x() * Unit::mm;
69  double yLocal = localPosition.y() * Unit::mm;
70 
71  // apply quantum efficiency if not yet done
72  bool applyQE = true;
73  double fraction = 1;
74  auto* info = dynamic_cast<Simulation::TrackInfo*>(photon.GetUserInformation());
75  if (info) {
76  applyQE = info->getStatus() < 3;
77  fraction = info->getFraction();
78  }
79  if (applyQE) {
80  double qeffi = m_topgp->getPMTEfficiency(energy, moduleID, pmtID, xLocal, yLocal);
81  if (qeffi == 0 or gRandom->Uniform() * fraction > qeffi) {
82  photon.SetTrackStatus(fStopAndKill);
83  return false;
84  }
85  info->setStatus(3);
86  info->setFraction(qeffi);
87  }
88 
89  // photon at detection
90  const G4ThreeVector& g_detPoint = photon.GetPosition();
91  const G4ThreeVector& g_detMomDir = photon.GetMomentumDirection();
92  TVector3 detPoint(g_detPoint.x(), g_detPoint.y(), g_detPoint.z());
93  TVector3 detMomDir(g_detMomDir.x(), g_detMomDir.y(), g_detMomDir.z());
94  double detTime = photon.GetGlobalTime();
95  double length = photon.GetTrackLength();
96 
97  // photon at emission
98  const G4ThreeVector& g_emiPoint = photon.GetVertexPosition();
99  const G4ThreeVector& g_emiMomDir = photon.GetVertexMomentumDirection();
100  TVector3 emiPoint(g_emiPoint.x(), g_emiPoint.y(), g_emiPoint.z());
101  TVector3 emiMomDir(g_emiMomDir.x(), g_emiMomDir.y(), g_emiMomDir.z());
102  double emiTime = photon.GetGlobalTime() - photon.GetLocalTime();
103 
104  // convert to Basf2 units
105  emiPoint = emiPoint * Unit::mm;
106  detPoint = detPoint * Unit::mm;
107  length = length * Unit::mm;
108 
109  // write to store arrays; add relations
110 
111  auto* simHit = m_simHits.appendNew(moduleID, pmtID, xLocal, yLocal,
112  detTime, energy);
113 
114  int parentID = photon.GetParentID();
115  if (parentID == 0) parentID = photon.GetTrackID();
116  m_relParticleHit.add(parentID, simHit->getArrayIndex());
117 
118  const auto* geo = m_topgp->getGeometry();
119  if (geo->isModuleIDValid(moduleID)) {
120  // transform to local frame
121  const auto& module = geo->getModule(moduleID);
122  emiPoint = module.pointToLocal(emiPoint);
123  detPoint = module.pointToLocal(detPoint);
124  emiMomDir = module.momentumToLocal(emiMomDir);
125  detMomDir = module.momentumToLocal(detMomDir);
126  } else {
127  B2ERROR("TOP::SensitivePMT: undefined module ID."
128  << LogVar("moduleID", moduleID));
129  }
130  auto* simPhoton = m_simPhotons.appendNew(moduleID,
131  emiPoint, emiMomDir, emiTime,
132  detPoint, detMomDir, detTime,
133  length, energy);
134 
135  simHit->addRelationTo(simPhoton);
136 
137  // kill photon after detection
138  photon.SetTrackStatus(fStopAndKill);
139 
140  return true;
141  }
142 
143 
144  } // end of namespace top
146 } // end of namespace Belle2
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::SensitivePMT::m_topgp
TOPGeometryPar * m_topgp
geometry parameters
Definition: SensitivePMT.h:72
Belle2::TOP::SensitivePMT::m_simHits
StoreArray< TOPSimHit > m_simHits
collection of simulated hits
Definition: SensitivePMT.h:75
Belle2::TOP::SensitivePMT::m_mcParticles
StoreArray< MCParticle > m_mcParticles
collection of MC particles
Definition: SensitivePMT.h:74
Belle2::Simulation::UserInfo
UserInfo class which is used to attach additional information to Geant4 particles and tracks.
Definition: UserInfo.h:46
Belle2::DataStore::c_DontWriteOut
@ c_DontWriteOut
Object/array should be NOT saved by output modules.
Definition: DataStore.h:73
Belle2::TOP::SensitivePMT::step
bool step(G4Step *aStep, G4TouchableHistory *) override
Process each step, fill TOPSimHits and TOPSimPhotons.
Definition: SensitivePMT.cc:51
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::SensitivePMT::m_simPhotons
StoreArray< TOPSimPhoton > m_simPhotons
collection of simulated photons
Definition: SensitivePMT.h:76
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::TOP::SensitivePMT::m_relParticleHit
RelationArray m_relParticleHit
relations
Definition: SensitivePMT.h:77
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::SensitivePMT::m_moduleReplicaDepth
int m_moduleReplicaDepth
replica depth of module volume
Definition: SensitivePMT.h:71
Belle2::Unit::mm
static const double mm
[millimeters]
Definition: Unit.h:80
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
Belle2::TOP::TOPGeometryPar::getPMTEfficiency
double getPMTEfficiency(double energy, int moduleID, int pmtID, double x, double y) const
Returns PMT pixel efficiency, a product of quantum and collection efficiency.
Definition: TOPGeometryPar.cc:191
Belle2::DataStore::c_Event
@ c_Event
Different object in each event, all objects/arrays are invalidated after event() function has been ca...
Definition: DataStore.h:61
Belle2::TOP::SensitivePMT::m_pmtReplicaDepth
int m_pmtReplicaDepth
replica depth of PMT volume
Definition: SensitivePMT.h:70