Belle II Software  release-05-01-25
BkgSensitiveDetector.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 and Marko Petric *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <simulation/background/BkgSensitiveDetector.h>
12 #include <simulation/background/BkgNeutronWeight.h>
13 #include <simulation/dataobjects/BeamBackHit.h>
14 #include <framework/gearbox/Unit.h>
15 
16 #include <G4Track.hh>
17 #include <G4Step.hh>
18 
19 #include <framework/datastore/StoreArray.h>
20 #include <framework/datastore/RelationArray.h>
21 #include <TVector3.h>
22 
23 using namespace std;
24 
25 namespace Belle2 {
32  BkgSensitiveDetector::BkgSensitiveDetector(const char* subDett, int iden):
33  Simulation::SensitiveDetectorBase("BKG", Const::invalidDetector), m_trackID(0), m_startPos(0., 0., 0.), m_startMom(0., 0., 0.),
34  m_startTime(0), m_startEnergy(0), m_energyDeposit(0), m_trackLength(0.)
35  {
36  // registration of store arrays and relations
37 
38  StoreArray<MCParticle> mcParticles;
39  StoreArray<BeamBackHit> beamBackHits;
40 
41  beamBackHits.registerInDataStore();
42  mcParticles.registerRelationTo(beamBackHits);
43 
44  // additional registration of MCParticle relation
45 
46  RelationArray relation(mcParticles, beamBackHits);
48 
49  m_eclrepscale = 0;
50  std::string subDet = subDett;
51  if (subDet == "IR") m_subDet = 0;
52  else if (subDet == "PXD") m_subDet = 1;
53  else if (subDet == "SVD") m_subDet = 2;
54  else if (subDet == "CDC") m_subDet = 3;
55  else if (subDet == "ARICH") m_subDet = 4;
56  else if (subDet == "TOP") m_subDet = 5;
57  else if (subDet == "ECL") m_subDet = 6;
58  else if (subDet == "EKLM") m_subDet = 7;
59  else if (subDet == "BKLM") m_subDet = 8;
60  else m_subDet = 99;
61 
62  if (m_subDet == 6) {
63  const int nc[] = {3, 4, 6, 9, 2, 9, 6, 4}; // the number of crystals in a sector per ring
64  const int indx[] = {96, 288, 864, 1151, 7776, 8064, 8544, 8736}; // regions with the same number of crystals in a ring
65  std::vector<int> indxv(indx, indx + sizeof(indx) / sizeof(int));
66  m_eclrepscale = nc[upper_bound(indxv.begin(), indxv.end(), iden) - indxv.begin()];
67  }
68 
69  m_identifier = iden;
70 
71  }
72 
73 
74  G4bool BkgSensitiveDetector::step(G4Step* aStep, G4TouchableHistory*)
75  {
76 
77  const G4StepPoint& preStep = *aStep->GetPreStepPoint();
78  G4Track& track = *aStep->GetTrack();
79 
80  if (m_trackID != track.GetTrackID()) {
81  //TrackID changed, store track informations
82  m_trackID = track.GetTrackID();
83  //Get world position
84  const G4ThreeVector& worldPosition = preStep.GetPosition();
85  m_startPos.SetXYZ(worldPosition.x() * Unit::mm / Unit::cm , worldPosition.y() * Unit::mm / Unit::cm,
86  worldPosition.z() * Unit::mm / Unit::cm);
87  //Get momentum
88  const G4ThreeVector& momentum = preStep.GetMomentum() ;
89  m_startMom.SetXYZ(momentum.x() * Unit::MeV, momentum.y() * Unit::MeV ,
90  momentum.z() * Unit::MeV);
91  //Get time
92  m_startTime = preStep.GetGlobalTime();
93  //Get energy
94  m_startEnergy = preStep.GetKineticEnergy() * Unit::MeV;
95  //Reset energy deposit;
96  m_energyDeposit = 0;
97  //Reset track lenght;
98  m_trackLength = 0;
99  }
100  //Update energy deposit
101  m_energyDeposit += aStep->GetTotalEnergyDeposit() * Unit::MeV;
102 
103  m_trackLength += aStep->GetStepLength() * Unit::mm;
104 
105  //Save Hit if track leaves volume or is killed
106  if (track.GetNextVolume() != track.GetVolume() || track.GetTrackStatus() >= fStopAndKill) {
107  int pdgCode = track.GetDefinition()->GetPDGEncoding();
108  double endEnergy = track.GetKineticEnergy() * Unit::MeV;
109  double neutWeight = 0;
110  if (pdgCode == 2112) {
112  neutWeight = wt.getWeight(m_startEnergy / Unit::MeV);
113  }
114  StoreArray<BeamBackHit> beamBackHits;
115  int nentr = beamBackHits.getEntries();
116  int id = m_identifier;
117  if (m_subDet == 6) { // ECL
118  int sector = preStep.GetTouchable()->GetCopyNumber(1);
119  if ((sector == 0) && (m_eclrepscale == 2) && ((id & 1) != 0)) id += 144; // barrel odd-numbered diodes in sector 0 -> sector 144
120  id += sector * m_eclrepscale;
121  }
122  beamBackHits.appendNew(m_subDet, id, pdgCode, m_trackID, m_startPos, m_startMom, m_startTime, m_startEnergy, endEnergy,
123  m_energyDeposit, m_trackLength, neutWeight);
124 
125 
126  // create relation to MCParticle
127  StoreArray<MCParticle> mcParticles;
128  RelationArray relBeamBackHitToMCParticle(mcParticles, beamBackHits);
129  relBeamBackHitToMCParticle.add(m_trackID, nentr);
130 
131  //Reset TrackID
132  m_trackID = 0;
133  }
134 
135  return true;
136  }
137 
138 
140 } // 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::Unit::cm
static const double cm
Standard units with the value = 1.
Definition: Unit.h:57
Belle2::BkgSensitiveDetector::m_startEnergy
double m_startEnergy
particle energy at the entrance in volume
Definition: BkgSensitiveDetector.h:61
Belle2::BkgSensitiveDetector::m_eclrepscale
int m_eclrepscale
replica (=sector) scale in ECL
Definition: BkgSensitiveDetector.h:64
Belle2::RelationArray
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:72
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::BkgSensitiveDetector::m_startPos
TVector3 m_startPos
particle position at the entrance in volume
Definition: BkgSensitiveDetector.h:58
Belle2::BkgSensitiveDetector::m_subDet
int m_subDet
subdetector id number
Definition: BkgSensitiveDetector.h:55
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::BkgSensitiveDetector::m_startMom
TVector3 m_startMom
particle momentum at the entrance in volume
Definition: BkgSensitiveDetector.h:59
Belle2::BkgSensitiveDetector::m_trackID
int m_trackID
track id
Definition: BkgSensitiveDetector.h:57
Belle2::Unit::MeV
static const double MeV
[megaelectronvolt]
Definition: Unit.h:124
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::BkgSensitiveDetector::m_energyDeposit
double m_energyDeposit
energy deposited in volume
Definition: BkgSensitiveDetector.h:62
Belle2::BkgNeutronWeight::getInstance
static BkgNeutronWeight & getInstance()
Return a reference to the singleton BkgNeutronWeight instance.
Definition: BkgNeutronWeight.cc:30
Belle2::BkgSensitiveDetector::step
bool step(G4Step *aStep, G4TouchableHistory *) override
Process each step and calculate variables defined in PXDSimHit.
Definition: BkgSensitiveDetector.cc:74
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::Unit::mm
static const double mm
[millimeters]
Definition: Unit.h:80
Belle2::BkgNeutronWeight::getWeight
double getWeight(double ke)
Get weighting factor to convert a neutron to its 1-MeV equivalent.
Definition: BkgNeutronWeight.cc:36
Belle2::StoreArray< MCParticle >
Belle2::Const
This class provides a set of constants for the framework.
Definition: Const.h:36
Belle2::BkgNeutronWeight
The class to get the weighting factor for a 1-MeV-equivalent neutron flux on silicon.
Definition: BkgNeutronWeight.h:33
Belle2::BkgSensitiveDetector::m_startTime
double m_startTime
global time
Definition: BkgSensitiveDetector.h:60
Belle2::BkgSensitiveDetector::m_trackLength
double m_trackLength
length of the track in the volume
Definition: BkgSensitiveDetector.h:63
Belle2::BkgSensitiveDetector::m_identifier
int m_identifier
identifier of subdetector component
Definition: BkgSensitiveDetector.h:56
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226