Neutron/src/SubsystemNeutronReader.cc

0001 #include "SimMuon/Neutron/interface/SubsystemNeutronReader.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "FWCore/ParameterSet/interface/FileInPath.h"
0004 #include "CLHEP/Random/RandFlat.h"
0005 #include "CLHEP/Random/RandPoissonQ.h"
0006 #include "SimMuon/Neutron/src/NeutronReader.h"
0007 #include "SimMuon/Neutron/src/RootNeutronReader.h"
0008 #include "SimMuon/Neutron/src/AsciiNeutronReader.h"
0009 #include <algorithm>
0010
0011 using namespace std;
0012
0013 SubsystemNeutronReader::SubsystemNeutronReader(const edm::ParameterSet& pset)
0014     : theHitReader(nullptr),
0015       theLuminosity(pset.getParameter<double>("luminosity")),  // in units of 10^34
0016       theStartTime(pset.getParameter<double>("startTime")),
0017       theEndTime(pset.getParameter<double>("endTime")),
0018       theEventOccupancy(pset.getParameter<vector<double> >("eventOccupancy"))  // TODO make map
0019 {
0020   // 17.3 collisions per live bx, 79.5% of bx live
0021   float collisionsPerCrossing = 13.75 * theLuminosity;
0022   int windowSize = (int)((theEndTime - theStartTime) / 25.);
0023   theEventsInWindow = collisionsPerCrossing * windowSize;
0024   string reader = pset.getParameter<string>("reader");
0025   edm::FileInPath input = pset.getParameter<edm::FileInPath>("input");
0026   if (reader == "ASCII") {
0027     theHitReader = new AsciiNeutronReader(input.fullPath());
0028   } else if (reader == "ROOT") {
0029     theHitReader = new RootNeutronReader(input.fullPath());
0030   }
0031 }
0032
0033 SubsystemNeutronReader::~SubsystemNeutronReader() { delete theHitReader; }
0034
0035 void SubsystemNeutronReader::generateChamberNoise(int chamberType,
0036                                                   int chamberIndex,
0037                                                   edm::PSimHitContainer& result,
0038                                                   CLHEP::HepRandomEngine* engine) {
0039   // make sure this chamber hasn't been done before
0040   if (find(theChambersDone.begin(), theChambersDone.end(), chamberIndex) == theChambersDone.end()) {
0041     float meanNumberOfEvents = theEventOccupancy[chamberType - 1] * theEventsInWindow;
0042     CLHEP::RandPoissonQ randPoissonQ(*engine, meanNumberOfEvents);
0043     int nEventsToAdd = randPoissonQ.fire();
0044     //    LogDebug("NeutronReader") << "Number of neutron events to add: "
0045     //std::cout << "Number of neutron events to add for chamber type " << chamberType << " : "
0046     // << nEventsToAdd <<  " mean " << meanNumberOfEvents << std::endl;
0047     //                   << nEventsToAdd <<  " mean " << meanNumberOfEvents;
0048
0049     for (int i = 0; i < nEventsToAdd; ++i) {
0050       // find the time for this event
0051       float timeOffset = CLHEP::RandFlat::shoot(engine, theStartTime, theEndTime);
0052       vector<PSimHit> neutronHits;
0053       theHitReader->readNextEvent(chamberType, neutronHits);
0054
0055       for (vector<PSimHit>::const_iterator neutronHitItr = neutronHits.begin(); neutronHitItr != neutronHits.end();
0056            ++neutronHitItr) {
0057         const PSimHit& rawHit = *neutronHitItr;
0058         // do the time offset and local det id
0059         int det = detId(chamberIndex, rawHit.detUnitId());
0060         PSimHit hit(rawHit.entryPoint(),
0061                     rawHit.exitPoint(),
0062                     rawHit.pabs(),
0063                     rawHit.tof() + timeOffset,
0064                     rawHit.energyLoss(),
0065                     rawHit.particleType(),
0066                     det,
0067                     rawHit.trackId(),
0068                     rawHit.thetaAtEntry(),
0069                     rawHit.phiAtEntry(),
0070                     rawHit.processType());
0071         //std::cout << "NEWHIT " << hit << std::endl;
0072         result.push_back(hit);
0073       }
0074     }
0075     theChambersDone.push_back(chamberIndex);
0076   }
0077 }