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  * Alexandre Beaulieu
8  * *
9  * This software is provided "as is" without any warranty. *
10  **************************************************************************/
11 
12 #include <beast/csi/simulation/SensitiveDetector.h>
13 #include <beast/csi/dataobjects/CsiSimHit.h>
14 #include <beast/csi/dataobjects/CsiHit.h>
15 #include <framework/dataobjects/EventMetaData.h>
16 #include <beast/csi/geometry/CsiGeometryPar.h>
17 #include <framework/datastore/StoreArray.h>
18 #include <framework/datastore/RelationArray.h>
19 #include <framework/datastore/StoreObjPtr.h>
20 
21 #include <G4Track.hh>
22 #include <G4Step.hh>
23 
24 namespace Belle2 {
30  namespace csi {
31 
32  SensitiveDetector::SensitiveDetector():
33  Simulation::SensitiveDetectorBase("CsiSensitiveDetector", Const::invalidDetector)
34  {
35  m_hitNum = 0;
36  m_EvnetNumber = 0;
37  m_oldEvnetNumber = 0;
38  m_trackID = 0;
39  m_startTime = 0;
40  m_endTime = 0;
41  m_WightedTime = 0;
42  m_startEnergy = 0;
43  m_energyDeposit = 0;
44  m_trackLength = 0;
45  iECLCell = 0;
46  TimeIndex = 0;
47  local_pos = 0;
48  T_ave = 0;
49  firstcall = 0;
50  m_phiID = 0;
51  m_thetaID = 0;
52  m_cellID = 0;
53 
54 
55  // CsiHitIndex: array of 18 (Nb. crystals) x 80 (max time index (tof [ns]/100))
56  for (int i = 0; i < 18; i++) {for (int j = 0; j < 80; j++) { CsiHitIndex[i][j] = 0;}}
57 
58  //Make sure all collections are registered
59  StoreArray<MCParticle> mcParticles;
60  StoreArray<CsiSimHit> simHits;
61  StoreArray<CsiHit> csiHits;
62 
63  RelationArray simHitRel(mcParticles, simHits);
64  RelationArray csiHitRel(mcParticles, csiHits);
65 
66  //Register all collections we want to modify and require those we want to use
67  mcParticles.registerInDataStore();
68  simHits.registerInDataStore();
69  csiHits.registerInDataStore();
70  simHitRel.registerInDataStore();
71  csiHitRel.registerInDataStore();
72 
73  //Register the Relation so that the TrackIDs get replaced by the actual
74  //MCParticle indices after simulating the events. This is needed as
75  //secondary particles might not be stored so everything relating to those
76  //particles will be attributed to the last saved mother particle
77  registerMCParticleRelation(simHitRel);
78  registerMCParticleRelation(csiHitRel);
79  }
80 
81  SensitiveDetector::~SensitiveDetector()
82  {
83 
84  }
85 
86 
87 
88  bool SensitiveDetector::step(G4Step* aStep, G4TouchableHistory*)
89  {
90 
91  const G4StepPoint& preStep = *aStep->GetPreStepPoint();
92  const G4StepPoint& postStep = * aStep->GetPostStepPoint();
93 
94  G4Track& track = * aStep->GetTrack();
95  if (m_trackID != track.GetTrackID()) {
96  //TrackID changed, store track informations
97  m_trackID = track.GetTrackID();
98  //Get momentum
99  m_momentum = preStep.GetMomentum() ;
100  //Get energy
101  m_startEnergy = preStep.GetKineticEnergy() ;
102  //Reset energy deposit;
103  m_energyDeposit = 0;
104  //Reset Wighted Time;
105  m_WightedTime = 0;
106  //Reset m_WightedPos;
107  m_WightedPos.set(0, 0, 0);
108 
109  }
110  //Update energy deposit
111  m_energyDeposit += aStep->GetTotalEnergyDeposit() ;
112 
113  m_startTime = preStep.GetGlobalTime();
114  m_endTime = postStep.GetGlobalTime();
115  m_WightedTime += (m_startTime + m_endTime) / 2 * (aStep->GetTotalEnergyDeposit());
116 
117  m_startPos = preStep.GetPosition();
118  m_endPos = postStep.GetPosition();
119 
120  G4ThreeVector position(0.5 * (m_startPos + m_endPos));
121  m_WightedPos += position * (aStep->GetTotalEnergyDeposit());
122 
123  //Save Hit if track leaves volume or is killed
124  if (track.GetNextVolume() != track.GetVolume() || track.GetTrackStatus() >= fStopAndKill) {
125  int pdgCode = track.GetDefinition()->GetPDGEncoding();
126 
127  const G4VPhysicalVolume& v = * track.GetVolume();
128 
129  //Column z of the rotation matrix R is the result of R * zhat.
130  // This is what we want since G4Trap are created aligned along z
131  G4RotationMatrix rotCell = * v.GetObjectRotation();
132  // G4RotationMatrix rotCell = * v.GetRotation();
133  G4ThreeVector vecCell = rotCell.colZ();
134  G4ThreeVector posCell = v.GetTranslation();
135 
136  // Get layer ID
137  if (v.GetName().find("Crystal") != std::string::npos) {
138  CsiGeometryPar* csip = CsiGeometryPar::Instance();
139  m_cellID = csip->CsiVolNameToCellID(v.GetName());
140 
141  double dTotalEnergy = 1.0 / m_energyDeposit; //avoid the error no match for 'operator/'
142  if (m_energyDeposit > 0.) {
143  saveSimHit(m_cellID, m_trackID, pdgCode, m_WightedTime * dTotalEnergy,
144  m_energyDeposit, m_momentum, m_WightedPos * dTotalEnergy,
145  posCell, vecCell);
146  }
147  }
148 
149  //Reset TrackID
150  m_trackID = 0;
151  }
152 
153  return true;
154  }
155 
156  int SensitiveDetector::saveSimHit(
157  const G4int cellId,
158  const G4int trackID,
159  const G4int pid,
160  const G4double tof,
161  const G4double edep,
162  G4ThreeVector mom,
163  G4ThreeVector pos,
164  G4ThreeVector PosCell,
165  G4ThreeVector VecCell)
166  {
167 
168  //Get the datastore arrays
169  StoreArray<MCParticle> mcParticles;
170  StoreArray<CsiSimHit> simHits;
171 
172  RelationArray simHitRel(mcParticles, simHits);
173 
174  // convert Hep3Vectors to TVector3
175  TVector3 momentum(mom.getX() , mom.getY() , mom.getZ());
176  TVector3 position(pos.getX() , pos.getY() , pos.getZ());
177 
178  //Append
179  simHits.appendNew(cellId, trackID, pid, tof / CLHEP::ns, edep / CLHEP::GeV,
180  momentum * (1.0 / CLHEP::GeV), position * (1.0 / CLHEP::cm));
181 
182  int simhitNumber = simHits.getEntries() - 1;
183  simHitRel.add(trackID, simhitNumber);
184  B2DEBUG(150, "Saving CsiSimHit Number: " << simhitNumber);
185 
186  // Store CsiHit: total energy deposited per hit, with average time.
187 
188  // Old names:
189  /*
190  StoreArray<CsiHit> eclHitArray;
191  RelationArray eclHitRel(mcParticles, eclHitArray);
192  */
193 
194  // New names:
195 
196  StoreArray<CsiHit> csiHits;
197 
198  RelationArray csiHitRel(mcParticles, csiHits);
199 
200 
201  StoreObjPtr<EventMetaData> eventMetaDataPtr;
202  int m_currentEvnetNumber = eventMetaDataPtr->getEvent();
203 
204 
205  if (firstcall == 0 || m_currentEvnetNumber != m_oldEvnetNumber) {
206  m_oldEvnetNumber = m_currentEvnetNumber;
207  for (int iCSICell = 0; iCSICell < 18; iCSICell++) {
208  for (int TimeInd = 0; TimeInd < 80; TimeInd++) {
209  CsiHitIndex[iCSICell][TimeInd] = -1;
210  }
211  }
212  firstcall++;
213  }
214 
215 
216  if (m_currentEvnetNumber == m_oldEvnetNumber) {
217  if ((tof / CLHEP::ns) < 8000) {
218  TimeIndex = (int)((tof / CLHEP::ns) / 100);
219 
220  // local_pos is the distance between the hit and the PIN-diode end of the crystal.
221  local_pos = (15. - (pos - PosCell) * VecCell);
222 
223  // From the ECL package T_ave = 6.05 + 0.0749 * local_pos - 0.00112 * local_pos * local_pos + (tof / CLHEP::ns) ;
224 
225  //For CsI: we remove the offset (I hypothetize it if from photo-sensor delay)
226  // we remove the acceleration of the light
227  // we correct the velocity term using n_CsI=1.80
228  T_ave = 0.0600 * local_pos + (tof / CLHEP::ns) ;
229 
230 
231  double E_cell = (edep / CLHEP::GeV);
232  if (CsiHitIndex[cellId][TimeIndex] == -1) {
233  //m_hitNum = csiHits->GetLast() + 1;
234  //new(csiHits->AddrAt(m_hitNum)) CSIHit();
235  csiHits.appendNew();
236  m_hitNum = csiHits.getEntries() - 1;
237  CsiHitIndex[cellId][TimeIndex] = m_hitNum;
238  csiHits[m_hitNum]->setCellId(cellId);
239  csiHits[m_hitNum]->setEnergyDep(E_cell);
240  csiHits[m_hitNum]->setTimeAve(T_ave);
241  csiHits[m_hitNum]->setTimeVar(0);
242  csiHitRel.add(trackID, m_hitNum);
243 
244  } else {
245  m_hitNum = CsiHitIndex[cellId][TimeIndex];
246  double old_edep = csiHits[m_hitNum]->getEnergyDep();
247  double old_TimeAve = csiHits[m_hitNum]->getTimeAve();
248  double old_TimeVar = csiHits[m_hitNum]->getTimeVar();
249 
250  //double new_edep = E_cell + old_edep;
251  //double new_TimeAve = (old_edep * old_TimeAve + E_cell * T_ave) / new_edep;
252 
253  double temp = E_cell + old_edep;
254  double delta = T_ave - old_TimeAve;
255  double R = delta * E_cell / temp;
256  double new_TimeAve = old_TimeAve + R;
257  double new_TimeVar = old_TimeVar + old_edep * delta * R;
258  double new_edep = temp;
259 
260  //Algo from Wiki:
261  /*
262  for x, weight in dataWeightPairs: # Alternatively "for x, weight in zip(data, weights):"
263  temp = weight + sumweight
264  delta = x - mean
265  R = delta * weight / temp
266  mean = mean + R
267  M2 = M2 + sumweight * delta * R # Alternatively, "M2 = M2 + weight * delta * (x−mean)"
268  sumweight = temp
269 
270  variance_n = M2/sumweight
271  */
272 
273  csiHits[m_hitNum]->setEnergyDep(new_edep);
274  csiHits[m_hitNum]->setTimeAve(new_TimeAve);
275  csiHits[m_hitNum]->setTimeVar(new_TimeVar);
276  csiHitRel.add(trackID, m_hitNum);
277 
278  }
279  }
280  } else { B2ERROR("m_currentEvnetNumber ERROR: m_oldEvnetNumber==m_oldEvnetNumber"); }
281 
282  return (simhitNumber);
283  }//saveSimHit
284  } //csi namespace
286 } //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::csi::CsiGeometryPar
The Class for CSI Geometry Parameters.
Definition: CsiGeometryPar.h:54
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
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
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::csi::CsiGeometryPar::CsiVolNameToCellID
int CsiVolNameToCellID(const G4String VolumeName)
Get Cell Id.
Definition: CsiGeometryPar.cc:151
Belle2::StoreArray< MCParticle >
Belle2::claws::SensitiveDetector::step
bool step(G4Step *step, G4TouchableHistory *) override
Step processing method.
Definition: SensitiveDetector.cc:62
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition: StoreArray.h:226
Belle2::csi::SensitiveDetector::SensitiveDetector
SensitiveDetector()
Constructor.
Definition: SensitiveDetector.cc:41