Neutron/src/SubsystemNeutronWriter.cc

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     // write out the real DetId, not the local one
0045     useLocalDetId_ = false;
0046     // smear the times
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   // sort hits by chamber
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   // now write out each chamber's contents
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   // should instantiate one of every chamber type, just so
0103   // ROOT knows what file to associate them with
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)) {  // 1st in cluster
0122         startTime = tof;
0123         // set the time to be [t0, t0+25] at start of event
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   // get any leftover cluster
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   // prevent float precision loss
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 }