1 #include <svd/modules/svdBackground/SVDBackgroundModule.h>
3 #include <framework/datastore/DataStore.h>
4 #include <framework/datastore/StoreObjPtr.h>
5 #include <framework/datastore/StoreArray.h>
6 #include <framework/datastore/RelationArray.h>
8 #include <framework/dataobjects/FileMetaData.h>
9 #include <framework/dataobjects/BackgroundMetaData.h>
10 #include <mdst/dataobjects/MCParticle.h>
11 #include <svd/dataobjects/SVDSimHit.h>
12 #include <svd/dataobjects/SVDTrueHit.h>
13 #include <svd/dataobjects/SVDShaperDigit.h>
14 #include <svd/dataobjects/SVDCluster.h>
15 #include <svd/dataobjects/SVDEnergyDepositionEvent.h>
16 #include <svd/dataobjects/SVDNeutronFluxEvent.h>
17 #include <svd/dataobjects/SVDOccupancyEvent.h>
22 #include <boost/format.hpp>
34 double eToADU(
double charge)
36 double minADC = -96000;
37 double maxADC = 288000;
38 double unitADC = (maxADC - minADC) / 1024.0;
39 return round(std::min(maxADC, std::max(minADC, charge)) / unitADC);
53 Module(), m_outputDirectoryName(""),
54 m_doseReportingLevel(c_reportNTuple),
55 m_nfluxReportingLevel(c_reportNTuple),
56 m_occupancyReportingLevel(c_reportNTuple),
57 m_componentName(""), m_componentTime(0),
58 m_triggerWidth(5), m_acceptanceWidth(2.5),
59 m_nielNeutrons(new
TNiel(c_niel_neutronFile)),
60 m_nielProtons(new
TNiel(c_niel_protonFile)),
61 m_nielPions(new
TNiel(c_niel_pionFile)),
62 m_nielElectrons(new
TNiel(c_niel_electronFile))
65 setDescription(
"SVD background module");
66 setPropertyFlags(c_ParallelProcessingCertified);
68 addParam(
"componentName", m_componentName,
"Background component name to process", m_componentName);
69 addParam(
"componentTime", m_componentTime,
"Background component time", m_componentTime);
70 addParam(
"triggerWidth", m_triggerWidth,
"RMS of trigger time estimate in ns", m_triggerWidth);
71 addParam(
"acceptanceWidth", m_acceptanceWidth,
72 "A hit is accepted if arrived within +/- accpetanceWidth * RMS(hit time - trigger time) of trigger; in ns", m_acceptanceWidth);
73 addParam(
"doseReportingLevel", m_doseReportingLevel,
"0 - no data, 1 - summary only, 2 - summary + ntuple", m_doseReportingLevel);
74 addParam(
"nfluxReportingLevel", m_nfluxReportingLevel,
"0 - no data, 1 - summary only, 2 - summary + ntuple",
75 m_nfluxReportingLevel);
76 addParam(
"occupancyReportingLevel", m_occupancyReportingLevel,
"0 - no data, 1 - summary only, 2 - summary + ntuple",
77 m_occupancyReportingLevel);
78 addParam(
"outputDirectory", m_outputDirectoryName,
"Name of output directory", m_outputDirectoryName);
81 const TVector3& SVDBackgroundModule::pointToGlobal(
VxdID sensorID,
const TVector3& local)
83 static TVector3 result(0, 0, 0);
86 result = info.pointToGlobal(local);
90 const TVector3& SVDBackgroundModule::vectorToGlobal(
VxdID sensorID,
const TVector3& local)
92 static TVector3 result(0, 0, 0);
95 result = info.vectorToGlobal(local);
99 SVDBackgroundModule::~SVDBackgroundModule()
103 void SVDBackgroundModule::initialize()
114 RelationArray relDigitsMCParticles(storeDigits, storeMCParticles);
116 RelationArray relMCParticlesTrueHits(storeMCParticles, storeTrueHits);
117 RelationArray relTrueHitsSimHits(storeTrueHits, storeSimHits);
121 storeEnergyDeposits.registerInDataStore();
123 storeNeutronFluxes.registerInDataStore();
125 storeOccupancyEvents.registerInDataStore();
128 m_storeFileMetaDataName = storeFileMetaData.getName();
129 m_storeBgMetaDataName = storeBgMetaData.getName();
130 m_storeMCParticlesName = storeMCParticles.getName();
131 m_storeSimHitsName = storeSimHits.getName();
132 m_storeTrueHitsName = storeTrueHits.getName();
133 m_storeDigitsName = storeDigits.getName();
134 m_relDigitsMCParticlesName = relDigitsMCParticles.
getName();
135 m_relDigitsTrueHitsName = relDigitsTrueHits.
getName();
136 m_relParticlesTrueHitsName = relMCParticlesTrueHits.
getName();
137 m_relTrueHitsSimHitsName = relTrueHitsSimHits.
getName();
138 m_storeEnergyDepositsName = storeEnergyDeposits.getName();
139 m_storeNeutronFluxesName = storeNeutronFluxes.getName();
141 m_componentTime *= Unit::us;
142 m_acceptanceWidth *= Unit::ns;
143 m_triggerWidth *= Unit::ns;
146 void SVDBackgroundModule::beginRun()
150 void SVDBackgroundModule::event()
167 RelationArray relDigitsMCParticles(storeDigits, storeMCParticles, m_relDigitsMCParticlesName);
168 RelationArray relDigitsTrueHits(storeDigits, storeTrueHits, m_relDigitsTrueHitsName);
169 RelationArray relTrueHitsSimHits(storeTrueHits, storeSimHits, m_relTrueHitsSimHitsName);
170 RelationArray relTrueHitsMCParticles(storeMCParticles, storeTrueHits, m_relParticlesTrueHitsName);
173 double currentComponentTime = storeBgMetaData->getRealTime();
174 if (currentComponentTime != m_componentTime)
175 B2FATAL(
"Mismatch in component times:\n"
176 <<
"Steering file: " << m_componentTime <<
"\n"
177 <<
"Background file: " << currentComponentTime);
179 VxdID currentSensorID(0);
180 double currentSensorThickness(0);
181 double currentSensorMass(0);
182 double currentSensorArea(0);
185 if (m_doseReportingLevel > c_reportNone) {
186 B2DEBUG(100,
"Expo and dose");
187 currentSensorID.
setID(0);
188 for (
const SVDSimHit& hit : storeSimHits) {
190 VxdID sensorID = hit.getSensorID();
191 if (sensorID != currentSensorID) {
192 currentSensorID = sensorID;
193 currentSensorThickness = getSensorThickness(currentSensorID);
194 currentSensorMass = getSensorMass(currentSensorID);
195 currentSensorArea = getSensorArea(currentSensorID);
197 double hitEnergy = hit.getElectrons() * Const::ehEnergy;
199 m_sensorData[currentSensorID].m_dose +=
200 (hitEnergy / Unit::J) / (currentSensorMass / 1000) * (c_smy / currentComponentTime);
202 m_sensorData[currentSensorID].m_expo += hitEnergy / currentSensorArea / (currentComponentTime / Unit::s);
203 if (m_doseReportingLevel == c_reportNTuple) {
204 const TVector3 localPos = hit.getPosIn();
205 const TVector3 globalPos = pointToGlobal(currentSensorID, localPos);
206 float globalPosXYZ[3];
207 globalPos.GetXYZ(globalPosXYZ);
210 hit.getPDGcode(), hit.getGlobalTime(),
211 localPos.X(), localPos.Y(), globalPosXYZ, hitEnergy,
212 (hitEnergy / Unit::J) / (currentSensorMass / 1000) / (currentComponentTime / Unit::s),
213 (hitEnergy / Unit::J) / currentSensorArea / (currentComponentTime / Unit::s)
220 if (m_nfluxReportingLevel > c_reportNone) {
221 B2DEBUG(100,
"Neutron flux");
222 currentSensorID.
setID(0);
224 VxdID sensorID = hit.getSensorID();
226 if (sensorID != currentSensorID) {
227 currentSensorID = sensorID;
228 currentSensorThickness = getSensorThickness(currentSensorID);
229 currentSensorMass = getSensorMass(currentSensorID);
230 currentSensorArea = getSensorArea(currentSensorID);
233 TVector3 entryPos(hit.getEntryU(), hit.getEntryV(), hit.getEntryW());
234 TVector3 exitPos(hit.getExitU(), hit.getExitV(), hit.getExitW());
235 double stepLength = (exitPos - entryPos).Mag();
241 double minDistance = 1.0e10;
242 for (
const SVDSimHit& related : storeSimHits) {
243 double distance = (entryPos - related.getPosIn()).Mag();
244 if (distance < minDistance) {
245 minDistance = distance;
254 TVector3 hitMomentum(hit.getMomentum());
255 hitMomentum.SetX(std::isfinite(hitMomentum.X()) ? hitMomentum.X() : 0.0);
256 hitMomentum.SetY(std::isfinite(hitMomentum.Y()) ? hitMomentum.Y() : 0.0);
257 hitMomentum.SetZ(std::isfinite(hitMomentum.Z()) ? hitMomentum.Z() : 0.0);
259 double kineticEnergy(0.0);
260 double nielWeight(0.0);
263 kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
264 nielWeight = m_nielNeutrons->getNielFactor(kineticEnergy / Unit::MeV);
268 kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
269 nielWeight = m_nielProtons->getNielFactor(kineticEnergy / Unit::MeV);
271 if (pdg == 111 || pdg == 211) {
273 kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
274 nielWeight = m_nielPions->getNielFactor(kineticEnergy / Unit::MeV);
277 double m0 = 0.511e-3;
278 kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
279 nielWeight = m_nielElectrons->getNielFactor(kineticEnergy / Unit::MeV);
283 kineticEnergy = sqrt(hitMomentum.Mag2() + m0 * m0) - m0;
287 nielWeight = std::isfinite(nielWeight) ? nielWeight : 0.0;
288 m_sensorData[currentSensorID].m_neutronFlux += nielWeight * stepLength / currentSensorThickness / currentSensorArea /
289 currentComponentTime * c_smy;
292 if (m_nfluxReportingLevel == c_reportNTuple) {
293 TVector3 localPos(hit.getU(), hit.getV(), hit.getW());
294 const TVector3 globalPos = pointToGlobal(currentSensorID, localPos);
295 float globalPosXYZ[3];
296 globalPos.GetXYZ(globalPosXYZ);
297 TVector3 localMom = hit.getMomentum();
298 const TVector3 globalMom = vectorToGlobal(currentSensorID, localMom);
299 float globalMomXYZ[3];
300 globalMom.GetXYZ(globalMomXYZ);
304 hit.getU(), hit.getV(), globalPosXYZ, globalMomXYZ, kineticEnergy,
305 stepLength, nielWeight,
306 stepLength / currentSensorThickness / currentSensorArea / (currentComponentTime / Unit::s),
307 nielWeight * stepLength / currentSensorThickness / currentSensorArea / (currentComponentTime / Unit::s)
314 if (m_occupancyReportingLevel > c_reportNone) {
315 B2DEBUG(100,
"Fired strips");
316 currentSensorID.
setID(0);
317 double currentSensorUCut = 0;
318 double currentSensorVCut = 0;
320 std::map<VxdID, std::multiset<unsigned short> > firedStrips;
324 VxdID sensorID = digit.getSensorID();
325 if (sensorID != currentSensorID) {
326 currentSensorID = sensorID;
327 auto info = getInfo(sensorID);
328 currentSensorUCut = eToADU(3.0 * info.getElectronicNoiseU());
329 currentSensorVCut = eToADU(3.0 * info.getElectronicNoiseV());
331 B2DEBUG(30,
"MaxCharge: " << digit.getMaxADCCounts() <<
" threshold: " << (digit.isUStrip() ? currentSensorUCut :
333 if (digit.getMaxADCCounts() < (digit.isUStrip() ? currentSensorUCut : currentSensorVCut))
continue;
334 B2DEBUG(30,
"Passed.");
336 VxdID writeID(sensorID);
337 if (digit.isUStrip())
341 firedStrips[writeID].insert(digit.getCellID());
344 for (
auto idAndSet : firedStrips) {
345 bool isUStrip = (idAndSet.first.getSegmentNumber() == 0);
346 VxdID sensorID = idAndSet.first;
348 double sensorArea = getSensorArea(sensorID);
349 int nFired_APV = idAndSet.second.size();
351 for (
auto it = idAndSet.second.begin();
352 it != idAndSet.second.end();
353 it = idAndSet.second.upper_bound(*it)) nFired++;
354 double fired = nFired / (currentComponentTime / Unit::s) / sensorArea;
355 double fired_t = nFired_APV * c_APVCycleTime / (currentComponentTime / Unit::s) / sensorArea;
357 m_sensorData[sensorID].m_firedU += fired;
358 m_sensorData[sensorID].m_firedU_t += fired_t;
360 m_sensorData[sensorID].m_firedV += fired;
361 m_sensorData[sensorID].m_firedV_t += fired_t;
380 B2DEBUG(100,
"Occupancy");
381 currentSensorID.
setID(0);
382 double currentNoiseU = 0;
383 double currentNoiseV = 0;
386 for (
auto cluster : storeClsuters) {
387 VxdID sensorID = cluster.getSensorID();
388 if (currentSensorID != sensorID) {
389 currentSensorID = sensorID;
390 auto info = getInfo(sensorID);
391 currentNoiseU = eToADU(info.getElectronicNoiseU());
392 currentNoiseV = eToADU(info.getElectronicNoiseV());
393 nStripsU = info.getUCells();
394 nStripsV = info.getVCells();
396 bool isU = cluster.isUCluster();
397 double snr = (isU) ? cluster.getCharge() / currentNoiseU : cluster.getCharge() / currentNoiseV;
398 int nStrips = (isU) ? nStripsU : nStripsV;
399 double tau_error = 45 / snr * Unit::ns;
400 tau_error = sqrt(m_triggerWidth * m_triggerWidth + tau_error * tau_error);
401 double tau_acceptance = 2 * m_acceptanceWidth * tau_error;
402 double w_acceptance = tau_acceptance / currentComponentTime;
403 double w_acceptance_APV = c_APVCycleTime / currentComponentTime;
404 double occupancy = 1.0 / nStrips * cluster.getSize();
406 m_sensorData[sensorID].m_occupancyU += w_acceptance * occupancy;
407 m_sensorData[sensorID].m_occupancyU_APV += w_acceptance_APV * occupancy;
409 m_sensorData[sensorID].m_occupancyV += w_acceptance * occupancy;
410 m_sensorData[sensorID].m_occupancyV_APV += w_acceptance_APV * occupancy;
412 if (m_occupancyReportingLevel == c_reportNTuple) {
416 cluster.isUCluster(), cluster.getPosition(), cluster.getSize(),
417 cluster.getCharge(), snr, w_acceptance, w_acceptance * occupancy,
418 w_acceptance_APV * occupancy
425 void SVDBackgroundModule::endRun()
430 void SVDBackgroundModule::terminate()
434 string outfileName(m_outputDirectoryName + m_componentName +
"_summary.txt");
435 outfile.open(outfileName.c_str(), ios::out | ios::trunc);
436 outfile <<
"component_name\t"
437 <<
"component_time\t"
450 <<
"occupancy_u_APV\t"
453 double componentTime = m_componentTime / Unit::us;
454 for (
auto vxdSensor : m_sensorData) {
455 outfile << m_componentName.c_str() <<
"\t"
456 << componentTime <<
"\t"
457 << vxdSensor.first.getLayerNumber() <<
"\t"
458 << vxdSensor.first.getLadderNumber() <<
"\t"
459 << vxdSensor.first.getSensorNumber() <<
"\t"
460 << vxdSensor.second.m_dose <<
"\t"
461 << vxdSensor.second.m_expo <<
"\t"
462 << vxdSensor.second.m_neutronFlux <<
"\t"
463 << vxdSensor.second.m_firedU <<
"\t"
464 << vxdSensor.second.m_firedV <<
"\t"
465 << vxdSensor.second.m_firedU_t <<
"\t"
466 << vxdSensor.second.m_firedV_t <<
"\t"
467 << vxdSensor.second.m_occupancyU <<
"\t"
468 << vxdSensor.second.m_occupancyV <<
"\t"
469 << vxdSensor.second.m_occupancyU_APV <<
"\t"
470 << vxdSensor.second.m_occupancyV_APV