File indexing completed on 2025-04-04 01:27:02
0001 #include "SimMuon/Neutron/interface/SubsystemNeutronWriter.h"
0002 #include "FWCore/AbstractServices/interface/RandomNumberGenerator.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 #include "FWCore/Utilities/interface/EDMException.h"
0007 #include "SimMuon/Neutron/src/AsciiNeutronWriter.h"
0008 #include "SimMuon/Neutron/src/RootNeutronWriter.h"
0009 #include "SimMuon/Neutron/src/NeutronWriter.h"
0010 #include "SimMuon/Neutron/src/EDMNeutronWriter.h"
0011 #include "FWCore/ServiceRegistry/interface/Service.h"
0012
0013 #include "CLHEP/Random/RandFlat.h"
0014
0015 #include <algorithm>
0016
0017 using namespace std;
0018
0019 class SortByTime {
0020 public:
0021 bool operator()(const PSimHit& h1, const PSimHit& h2) { return (h1.tof() < h2.tof()); }
0022 };
0023
0024 SubsystemNeutronWriter::SubsystemNeutronWriter(edm::ParameterSet const& pset)
0025 : theHitWriter(nullptr),
0026 useRandFlat(false),
0027 theInputTag(pset.getParameter<edm::InputTag>("input")),
0028 theNeutronTimeCut(pset.getParameter<double>("neutronTimeCut")),
0029 theTimeWindow(pset.getParameter<double>("timeWindow")),
0030 theT0(pset.getParameter<double>("t0")),
0031 hitToken_(consumes<edm::PSimHitContainer>(theInputTag)),
0032 theNEvents(0),
0033 initialized(false),
0034 useLocalDetId_(true) {
0035 string writer = pset.getParameter<string>("writer");
0036 string output = pset.getParameter<string>("output");
0037 if (writer == "ASCII") {
0038 theHitWriter = new AsciiNeutronWriter(output);
0039 } else if (writer == "ROOT") {
0040 theHitWriter = new RootNeutronWriter(output);
0041 } else if (writer == "EDM") {
0042 produces<edm::PSimHitContainer>();
0043 theHitWriter = new EDMNeutronWriter();
0044
0045 useLocalDetId_ = false;
0046
0047 edm::Service<edm::RandomNumberGenerator> rng;
0048 if (!rng.isAvailable()) {
0049 throw cms::Exception("Configuration")
0050 << "SubsystemNeutronWriter requires the RandomNumberGeneratorService\n"
0051 "which is not present in the configuration file. You must add the service\n"
0052 "in the configuration file or remove the modules that require it.";
0053 }
0054 useRandFlat = true;
0055 } else {
0056 throw cms::Exception("NeutronWriter") << "Bad writer: " << writer;
0057 }
0058 }
0059
0060 SubsystemNeutronWriter::~SubsystemNeutronWriter() {
0061 printStats();
0062 delete theHitWriter;
0063 }
0064
0065 void SubsystemNeutronWriter::printStats() {
0066 edm::LogInfo("SubsystemNeutronWriter") << "SubsystemNeutronWriter Statistics:\n";
0067 for (map<int, int>::iterator mapItr = theCountPerChamberType.begin(); mapItr != theCountPerChamberType.end();
0068 ++mapItr) {
0069 edm::LogInfo("SubsystemNeutronWriter")
0070 << " theEventOccupancy[" << mapItr->first << "] = " << mapItr->second << " / " << theNEvents << " / NCT \n";
0071 }
0072 }
0073
0074 void SubsystemNeutronWriter::produce(edm::Event& e, edm::EventSetup const& c) {
0075 CLHEP::HepRandomEngine* engine = nullptr;
0076 if (useRandFlat) {
0077 edm::Service<edm::RandomNumberGenerator> rng;
0078 engine = &rng->getEngine(e.streamID());
0079 }
0080 theHitWriter->beginEvent(e, c);
0081 ++theNEvents;
0082 const edm::Handle<edm::PSimHitContainer>& hits = e.getHandle(hitToken_);
0083
0084
0085 std::map<int, edm::PSimHitContainer> hitsByChamber;
0086 for (edm::PSimHitContainer::const_iterator hitItr = hits->begin(); hitItr != hits->end(); ++hitItr) {
0087 int chamberIndex = chamberId(hitItr->detUnitId());
0088 hitsByChamber[chamberIndex].push_back(*hitItr);
0089 }
0090
0091
0092 for (std::map<int, edm::PSimHitContainer>::iterator hitsByChamberItr = hitsByChamber.begin();
0093 hitsByChamberItr != hitsByChamber.end();
0094 ++hitsByChamberItr) {
0095 int chambertype = chamberType(hitsByChamberItr->first);
0096 writeHits(chambertype, hitsByChamberItr->second, engine);
0097 }
0098 theHitWriter->endEvent();
0099 }
0100
0101 void SubsystemNeutronWriter::initialize(int chamberType) {
0102
0103
0104 theHitWriter->initialize(chamberType);
0105 }
0106
0107 void SubsystemNeutronWriter::writeHits(int chamberType,
0108 edm::PSimHitContainer& chamberHits,
0109 CLHEP::HepRandomEngine* engine) {
0110 sort(chamberHits.begin(), chamberHits.end(), SortByTime());
0111 edm::PSimHitContainer cluster;
0112 float startTime = -1000.;
0113 float smearing = 0.;
0114 for (size_t i = 0; i < chamberHits.size(); ++i) {
0115 PSimHit hit = chamberHits[i];
0116 float tof = hit.tof();
0117 LogDebug("SubsystemNeutronWriter") << "found hit from part type " << hit.particleType() << " at tof " << tof
0118 << " p " << hit.pabs() << " on det " << hit.detUnitId() << " chamber type "
0119 << chamberType;
0120 if (tof > theNeutronTimeCut) {
0121 if (tof > (startTime + theTimeWindow)) {
0122 startTime = tof;
0123
0124 smearing = theT0;
0125 if (useRandFlat) {
0126 smearing += CLHEP::RandFlat::shoot(engine, 25.);
0127 }
0128 if (!cluster.empty()) {
0129 LogDebug("SubsystemNeutronWriter") << "filling old cluster";
0130 writeCluster(chamberType, cluster);
0131 cluster.clear();
0132 }
0133 LogDebug("SubsystemNeutronWriter")
0134 << "starting neutron cluster at time " << startTime << " on detType " << chamberType;
0135 }
0136 adjust(hit, -1. * startTime, smearing);
0137 cluster.push_back(hit);
0138 }
0139 }
0140
0141 if (!cluster.empty()) {
0142 writeCluster(chamberType, cluster);
0143 }
0144 }
0145
0146 void SubsystemNeutronWriter::writeCluster(int chamberType, const edm::PSimHitContainer& cluster) {
0147 if (accept(cluster)) {
0148 theHitWriter->writeCluster(chamberType, cluster);
0149 updateCount(chamberType);
0150 }
0151 }
0152
0153 void SubsystemNeutronWriter::adjust(PSimHit& h, float timeOffset, float smearing) {
0154 unsigned int detId = useLocalDetId_ ? localDetId(h.detUnitId()) : h.detUnitId();
0155 float htime = h.timeOfFlight() + timeOffset + smearing;
0156
0157 if (h.timeOfFlight() > 1.E+6) {
0158 htime = smearing;
0159 }
0160 h = PSimHit(h.entryPoint(),
0161 h.exitPoint(),
0162 h.pabs(),
0163 htime,
0164 h.energyLoss(),
0165 h.particleType(),
0166 detId,
0167 h.trackId(),
0168 h.momentumAtEntry().theta(),
0169 h.momentumAtEntry().phi());
0170 }
0171
0172 void SubsystemNeutronWriter::updateCount(int chamberType) {
0173 map<int, int>::iterator entry = theCountPerChamberType.find(chamberType);
0174 if (entry == theCountPerChamberType.end()) {
0175 theCountPerChamberType.insert(pair<int, int>(chamberType, 1));
0176 } else {
0177 ++(entry->second);
0178 }
0179 }