CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/SimMuon/Neutron/src/SubsystemNeutronReader.cc

Go to the documentation of this file.
00001 #include "SimMuon/Neutron/interface/SubsystemNeutronReader.h"
00002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00003 #include "FWCore/ParameterSet/interface/FileInPath.h"
00004 #include "CLHEP/Random/RandFlat.h"
00005 #include "CLHEP/Random/RandPoissonQ.h"
00006 #include "SimMuon/Neutron/src/NeutronReader.h"
00007 #include "SimMuon/Neutron/src/RootNeutronReader.h"
00008 #include "SimMuon/Neutron/src/AsciiNeutronReader.h"
00009 #include <algorithm>
00010 
00011 using namespace std;
00012 
00013 SubsystemNeutronReader::SubsystemNeutronReader(const edm::ParameterSet & pset)
00014 : theHitReader(0),
00015   theRandFlat(0), 
00016   theRandPoisson(0),
00017   theLuminosity(pset.getParameter<double>("luminosity")), // in units of 10^34
00018   theStartTime(pset.getParameter<double>("startTime")), 
00019   theEndTime(pset.getParameter<double>("endTime")),
00020   theEventOccupancy(pset.getParameter<vector<double> >("eventOccupancy")) // TODO make map
00021 {
00022   // 17.3 collisions per live bx, 79.5% of bx live
00023   float collisionsPerCrossing = 13.75 * theLuminosity;
00024   int windowSize = (int)((theEndTime-theStartTime)/25.);
00025   theEventsInWindow = collisionsPerCrossing * windowSize;
00026   string reader = pset.getParameter<string>("reader");
00027   edm::FileInPath input = pset.getParameter<edm::FileInPath>("input");
00028   if(reader == "ASCII")
00029   {
00030     theHitReader = new AsciiNeutronReader(input.fullPath());
00031   }
00032   else if (reader == "ROOT")
00033   {
00034     theHitReader = new RootNeutronReader(input.fullPath());
00035   }
00036 }
00037 
00038 
00039 SubsystemNeutronReader::~SubsystemNeutronReader() {
00040   delete theHitReader;
00041   delete theRandFlat;
00042   delete theRandPoisson;
00043 }
00044 
00045 
00046 void SubsystemNeutronReader::setRandomEngine(CLHEP::HepRandomEngine & engine)
00047 {
00048   theRandFlat = new CLHEP::RandFlat(engine);
00049   theRandPoisson = new CLHEP::RandPoissonQ(engine);
00050 }
00051 
00052 
00053 void
00054 SubsystemNeutronReader::generateChamberNoise(int chamberType, int chamberIndex, 
00055                                              edm::PSimHitContainer & result) 
00056 {
00057   // make sure this chamber hasn't been done before
00058   if(find(theChambersDone.begin(), theChambersDone.end(), chamberIndex) 
00059      == theChambersDone.end()) 
00060   {
00061     float meanNumberOfEvents = theEventOccupancy[chamberType-1] 
00062                              * theEventsInWindow;
00063     int nEventsToAdd = theRandPoisson->fire(meanNumberOfEvents);
00064 //    LogDebug("NeutronReader") << "Number of neutron events to add: " 
00065 //std::cout << "Number of neutron events to add for chamber type " << chamberType << " : " 
00066 // << nEventsToAdd <<  " mean " << meanNumberOfEvents << std::endl;
00067 //                   << nEventsToAdd <<  " mean " << meanNumberOfEvents;
00068 
00069     for(int i = 0; i < nEventsToAdd; ++i) {
00070       // find the time for this event
00071       float timeOffset = theRandFlat->fire(theStartTime, theEndTime);
00072       vector<PSimHit> neutronHits;
00073       theHitReader->readNextEvent(chamberType, neutronHits);
00074 
00075       for( vector<PSimHit>::const_iterator neutronHitItr = neutronHits.begin();
00076            neutronHitItr != neutronHits.end(); ++neutronHitItr)
00077       {
00078          const PSimHit & rawHit = *neutronHitItr;
00079          // do the time offset and local det id
00080          int det = detId(chamberIndex, rawHit.detUnitId());
00081          PSimHit hit(rawHit.entryPoint(), rawHit.exitPoint(), rawHit.pabs(),
00082                      rawHit.tof()+timeOffset,
00083                      rawHit.energyLoss(), rawHit.particleType(),
00084                      det, rawHit.trackId(),
00085                      rawHit.thetaAtEntry(),  rawHit.phiAtEntry(), rawHit.processType());
00086 //std::cout << "NEWHIT " << hit << std::endl;
00087          result.push_back(hit);
00088       }
00089 
00090     }
00091     theChambersDone.push_back(chamberIndex);
00092   }
00093 }
00094