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: Luka Santelj *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <arich/simulation/SensitiveDetector.h>
12 #include <arich/dataobjects/ARICHSimHit.h>
13 
14 // geant4
15 #include <G4OpBoundaryProcess.hh>
16 #include <G4ProcessManager.hh>
17 
18 
19 #include <simulation/kernel/UserInfo.h>
20 #include <G4Track.hh>
21 #include <G4Step.hh>
22 
23 #include <framework/datastore/StoreArray.h>
24 #include <framework/datastore/RelationArray.h>
25 #include <TVector2.h>
26 #include <TRandom3.h>
27 
28 using namespace std;
29 
30 namespace Belle2 {
35  namespace arich {
36 
37  SensitiveDetector::SensitiveDetector():
38  Simulation::SensitiveDetectorBase("ARICH", Const::ARICH)
39  {
40 
41  StoreArray<MCParticle> particles;
43  RelationArray relation(particles, hits);
45  hits.registerInDataStore();
46  particles.registerRelationTo(hits);
47 
48  }
49 
50 
51  G4bool SensitiveDetector::step(G4Step* aStep, G4TouchableHistory*)
52  {
53  //Get particle ID
54  G4Track& track = *aStep->GetTrack();
55  if (track.GetDefinition()->GetParticleName() != "opticalphoton") return false;
56 
57  //Get time (check for proper global time)
58  const G4double globalTime = track.GetGlobalTime();
59 
60  const G4StepPoint& postStep = *aStep->GetPostStepPoint();
61  const G4ThreeVector& postworldPosition = postStep.GetPosition();
62 
63  const G4StepPoint& preStep = *aStep->GetPreStepPoint();
64  const G4ThreeVector& preworldPosition = preStep.GetPosition();
65 
66  // direction of photon (+1 is into hapd)
67  int dir = -1;
68  if (postworldPosition.z() - preworldPosition.z() > 0) dir = +1;
69 
70  // trigger only on window geometrical boundary
71  if (dir == 1 && postStep.GetStepStatus() != fGeomBoundary) { return false;}
72  if (dir == -1 && preStep.GetStepStatus() != fGeomBoundary) { return false;}
73 
74  if ((track.GetNextVolume())->GetName() == "ARICH.HAPDWall") {
75  track.SetTrackStatus(fStopAndKill);
76  return false;
77  }
78 
79  // Check if photon is internally reflected in HAPD window
80  G4OpBoundaryProcessStatus theStatus = Undefined;
81 
82  G4ProcessManager* OpManager = G4OpticalPhoton::OpticalPhoton()->GetProcessManager();
83 
84  if (OpManager) {
85  G4int MAXofPostStepLoops =
86  OpManager->GetPostStepProcessVector()->entries();
87  G4ProcessVector* fPostStepDoItVector =
88  OpManager->GetPostStepProcessVector(typeDoIt);
89  for (G4int i = 0; i < MAXofPostStepLoops; i++) {
90  G4VProcess* fCurrentProcess = (*fPostStepDoItVector)[i];
91  G4OpBoundaryProcess* opProcess = dynamic_cast<G4OpBoundaryProcess*>(fCurrentProcess);
92  if (opProcess) { theStatus = opProcess->GetStatus(); break;}
93  }
94  }
95 
96  // if photon is internally reflected and going backward, do nothing
97  if (theStatus == 3 && dir < 0) return 0;
98  if (theStatus == 4 && dir < 0) { if (gRandom->Uniform() < 0.5) track.SetTrackStatus(fStopAndKill); return 0; }
99 
100  // apply quantum efficiency if not yet done
101  bool applyQE = true;
102  Simulation::TrackInfo* info =
103  dynamic_cast<Simulation::TrackInfo*>(track.GetUserInformation());
104  if (info) applyQE = info->getStatus() < 2;
105 
106  if (applyQE) {
107 
108  double energy = track.GetKineticEnergy() / CLHEP::eV;
109  double qeffi = m_simPar->getQE(energy) * m_simPar->getColEff();
110 
111  double fraction = info->getFraction();
112  //correct Q.E. for internally reflected photons
113  if (theStatus == 3) qeffi *= m_simPar->getQEScaling();
114  if (gRandom->Uniform() * fraction > qeffi) {
115  // apply possible absorbtion in HAPD window (for internally reflected photons only)
116  if (theStatus == 3 && gRandom->Uniform() < m_simPar->getWindowAbsorbtion()) track.SetTrackStatus(fStopAndKill);
117  return false;
118  }
119  }
120 
121  //Get module ID number
122  const G4int moduleID = dir > 0 ? postStep.GetTouchableHandle()->GetReplicaNumber(1) :
123  preStep.GetTouchableHandle()->GetReplicaNumber(2);
124 
125  if (moduleID <= 0) return false;
126  //std::cout << "moduleID " << moduleID << std::endl;
127  // std::cout << m_arichgp->getColEff() << " " << info->getFraction() << " " << m_arichgp->getQEScaling() <<" " << m_arichgp->QE(2.0) << std::endl;
128  //Transform to local position
129  G4ThreeVector localPosition = dir > 0 ? postStep.GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(
130  postworldPosition) :
131  preStep.GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(preworldPosition);
132 
133  //Get photon energy
134  const G4double energy = track.GetKineticEnergy() / CLHEP::eV;
135 
136  //------------------------------------------------------------
137  // Create ARICHSimHit and save it to datastore
138  //------------------------------------------------------------
139 
140  TVector2 locpos(localPosition.x() / CLHEP::cm, localPosition.y() / CLHEP::cm);
141  StoreArray<ARICHSimHit> arichSimHits;
142  ARICHSimHit* simHit = arichSimHits.appendNew(moduleID, locpos, globalTime, energy);
143 
144  // add relation to MCParticle
145  StoreArray<MCParticle> mcParticles;
146  RelationArray arichSimHitRel(mcParticles, arichSimHits);
147  arichSimHitRel.add(track.GetParentID(), simHit->getArrayIndex());
148 
149  // after detection photon track is killed
150  track.SetTrackStatus(fStopAndKill);
151  return true;
152  }
153 
154 
155  } // end of namespace arich
157 } // end of namespace Belle2
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::Simulation::UserInfo
UserInfo class which is used to attach additional information to Geant4 particles and tracks.
Definition: UserInfo.h:46
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::arich::SensitiveDetector::m_simPar
DBObjPtr< ARICHSimulationPar > m_simPar
simulation parameters from the DB
Definition: SensitiveDetector.h:56
Belle2::ARICHSimHit
Class ARICHSimHit - Geant4 simulated hit for ARICH.
Definition: ARICHSimHit.h:39
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::RelationsInterface::getArrayIndex
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
Definition: RelationsObject.h:387
Belle2::StoreArray< MCParticle >
Belle2::Const
This class provides a set of constants for the framework.
Definition: Const.h:36