181 double currentComponentTime = storeBgMetaData->getRealTime();
183 B2FATAL(
"Mismatch in component times:\n"
185 <<
"Background file: " << currentComponentTime);
187 VxdID currentSensorID(0);
188 double currentSensorThickness(0);
189 double currentSensorArea(0);
193 B2DEBUG(100,
"Expo and dose");
194 currentSensorID.
setID(0);
195 double currentSensorMass(0);
197 for (
const SVDSimHit& hit : storeSimHits) {
199 VxdID sensorID = hit.getSensorID();
200 if (sensorID != currentSensorID) {
201 currentSensorID = sensorID;
209 (hitEnergy /
Unit::J) / (currentSensorMass / 1000) * (
c_smy / currentComponentTime);
211 m_sensorData[currentSensorID].m_expo += hitEnergy / currentSensorArea / (currentComponentTime /
Unit::s);
213 const ROOT::Math::XYZVector localPos = hit.getPosIn();
214 const ROOT::Math::XYZVector globalPos =
pointToGlobal(currentSensorID, localPos);
215 float globalPosXYZ[3];
216 globalPos.GetCoordinates(globalPosXYZ);
219 hit.getPDGcode(), hit.getGlobalTime(),
220 localPos.X(), localPos.Y(), globalPosXYZ, hitEnergy,
221 (hitEnergy /
Unit::J) / (currentSensorMass / 1000) / (currentComponentTime /
Unit::s),
222 (hitEnergy /
Unit::J) / currentSensorArea / (currentComponentTime /
Unit::s)
230 B2DEBUG(100,
"Neutron flux");
231 currentSensorID.
setID(0);
233 VxdID sensorID = hit.getSensorID();
235 if (sensorID != currentSensorID) {
236 currentSensorID = sensorID;
242 ROOT::Math::XYZVector entryPos(hit.getEntryU(), hit.getEntryV(), hit.getEntryW());
243 ROOT::Math::XYZVector exitPos(hit.getExitU(), hit.getExitV(), hit.getExitW());
244 double stepLength = (exitPos - entryPos).
R();
250 double minDistance = 1.0e10;
251 for (
const SVDSimHit& related : storeSimHits) {
252 double distance = (entryPos - related.getPosIn()).R();
253 if (distance < minDistance) {
254 minDistance = distance;
260 B2WARNING(
"No related SVDSimHit found");
267 ROOT::Math::XYZVector hitMomentum(hit.getMomentum());
268 hitMomentum.SetX(std::isfinite(hitMomentum.X()) ? hitMomentum.X() : 0.0);
269 hitMomentum.SetY(std::isfinite(hitMomentum.Y()) ? hitMomentum.Y() : 0.0);
270 hitMomentum.SetZ(std::isfinite(hitMomentum.Z()) ? hitMomentum.Z() : 0.0);
272 double kineticEnergy(0.0);
273 double nielWeight(0.0);
276 kineticEnergy =
sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
281 kineticEnergy =
sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
286 kineticEnergy =
sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
291 kineticEnergy =
sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
296 kineticEnergy =
sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
300 nielWeight = std::isfinite(nielWeight) ? nielWeight : 0.0;
301 m_sensorData[currentSensorID].m_neutronFlux += nielWeight * stepLength / currentSensorThickness / currentSensorArea /
302 currentComponentTime *
c_smy;
306 ROOT::Math::XYZVector localPos(hit.getU(), hit.getV(), hit.getW());
307 const ROOT::Math::XYZVector globalPos =
pointToGlobal(currentSensorID, localPos);
308 float globalPosXYZ[3];
309 globalPos.GetCoordinates(globalPosXYZ);
310 ROOT::Math::XYZVector localMom = hit.getMomentum();
311 const ROOT::Math::XYZVector globalMom =
vectorToGlobal(currentSensorID, localMom);
312 float globalMomXYZ[3];
313 globalMom.GetCoordinates(globalMomXYZ);
317 hit.getU(), hit.getV(), globalPosXYZ, globalMomXYZ, kineticEnergy,
318 stepLength, nielWeight,
319 stepLength / currentSensorThickness / currentSensorArea / (currentComponentTime /
Unit::s),
320 nielWeight * stepLength / currentSensorThickness / currentSensorArea / (currentComponentTime /
Unit::s)
328 B2DEBUG(100,
"Fired strips");
329 currentSensorID.
setID(0);
330 double currentSensorUCut = 0;
331 double currentSensorVCut = 0;
333 std::map<VxdID, std::multiset<unsigned short> > firedStrips;
337 VxdID sensorID = digit.getSensorID();
338 if (sensorID != currentSensorID) {
339 currentSensorID = sensorID;
341 currentSensorUCut = eToADU(3.0 * info.getElectronicNoiseU());
342 currentSensorVCut = eToADU(3.0 * info.getElectronicNoiseV());
344 B2DEBUG(30,
"MaxCharge: " << digit.getMaxADCCounts() <<
" threshold: " << (digit.isUStrip() ? currentSensorUCut :
346 if (digit.getMaxADCCounts() < (digit.isUStrip() ? currentSensorUCut : currentSensorVCut))
continue;
347 B2DEBUG(30,
"Passed.");
349 VxdID writeID(sensorID);
350 if (digit.isUStrip())
354 firedStrips[writeID].insert(digit.getCellID());
357 for (
auto idAndSet : firedStrips) {
358 bool isUStrip = (idAndSet.first.getSegmentNumber() == 0);
359 VxdID sensorID = idAndSet.first;
362 int nFired_APV = idAndSet.second.size();
364 for (
auto it = idAndSet.second.begin();
365 it != idAndSet.second.end();
366 it = idAndSet.second.upper_bound(*it)) nFired++;
367 double fired = nFired / (currentComponentTime /
Unit::s) / sensorArea;
393 B2DEBUG(100,
"Occupancy");
394 currentSensorID.
setID(0);
395 double currentNoiseU = 0;
396 double currentNoiseV = 0;
399 for (
auto cluster : storeClsuters) {
400 VxdID sensorID = cluster.getSensorID();
401 if (currentSensorID != sensorID) {
402 currentSensorID = sensorID;
404 currentNoiseU = eToADU(info.getElectronicNoiseU());
405 currentNoiseV = eToADU(info.getElectronicNoiseV());
406 nStripsU = info.getUCells();
407 nStripsV = info.getVCells();
409 bool isU = cluster.isUCluster();
410 double snr = (isU) ? cluster.getCharge() / currentNoiseU : cluster.getCharge() / currentNoiseV;
411 int nStrips = (isU) ? nStripsU : nStripsV;
412 double tau_error = 45 / snr *
Unit::ns;
415 double w_acceptance = tau_acceptance / currentComponentTime;
417 double occupancy = 1.0 / nStrips * cluster.getSize();
419 m_sensorData[sensorID].m_occupancyU += w_acceptance * occupancy;
420 m_sensorData[sensorID].m_occupancyU_APV += w_acceptance_APV * occupancy;
422 m_sensorData[sensorID].m_occupancyV += w_acceptance * occupancy;
423 m_sensorData[sensorID].m_occupancyV_APV += w_acceptance_APV * occupancy;
429 cluster.isUCluster(), cluster.getPosition(), cluster.getSize(),
430 cluster.getCharge(), snr, w_acceptance, w_acceptance * occupancy,
431 w_acceptance_APV * occupancy