CTPPSRecHitProducer/plugins/CTPPSRecHitProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    FastSimulation/CTPPSRecHitProducer
0004 // Class:      CTPPSRecHitProducer
0005 //
0006 /**\class CTPPSRecHitProducer CTPPSRecHitProducer.cc FastSimulation/CTPPSRecHitProducer/plugins/CTPPSRecHitProducer.cc
0007
0008 Description: [one line class summary]
0009
0010 Implementation:
0011 [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Dilson De Jesus Damiao
0015 //         Created:  Mon, 19 Sep 2016 18:10:30 GMT
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/stream/EDProducer.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/Utilities/interface/StreamID.h"
0031 /////
0032 #include <cmath>
0033
0034 #include <CLHEP/Random/RandGauss.h>
0035 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0036 #include "FWCore/Utilities/interface/Exception.h"
0037 #include "FWCore/ServiceRegistry/interface/Service.h"
0038
0039 #include "FWCore/Utilities/interface/InputTag.h"
0040 #include "FWCore/Utilities/interface/EDGetToken.h"
0041 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0042 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0043 #include <vector>
0044 #include "FWCore/Framework/interface/ConsumesCollector.h"
0045
0046 #include "FastSimDataFormats/CTPPSFastSim/interface/CTPPSFastRecHitContainer.h"
0047 #include "FastSimulation/CTPPSFastGeometry/interface/CTPPSToFDetector.h"
0048 #include "Utilities/PPS/interface/PPSUnitConversion.h"
0049
0050 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0051 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0052 #include <iostream>
0053
0054 #include <string>
0055
0056 //Mixing Collection and CF info
0057 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0058 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
0059 #include "HepMC/GenEvent.h"
0060
0061 //
0062 // class declaration
0063 //
0064 class TRandom3;
0065 class CTPPSToFDetector;
0066
0067 class CTPPSRecHitProducer : public edm::stream::EDProducer<> {
0068 public:
0069   explicit CTPPSRecHitProducer(const edm::ParameterSet&);
0070   ~CTPPSRecHitProducer() override;
0071
0072 private:
0073   void beginStream(edm::StreamID) override;
0074   void produce(edm::Event&, const edm::EventSetup&) override;
0075   void endStream() override;
0076
0077   // ----------member data ---------------------------
0078   typedef std::vector<PSimHit> PSimHitContainer;
0079
0080   //    std::string mix_;
0081   //    std::string collection_for_XF;
0082
0083   edm::EDGetTokenT<CrossingFrame<PSimHit> > cf_token;
0084
0085   double fTrackerWidth, fTrackerHeight, fTrackerInsertion, fBeamXRMS_Trk1, fBeamXRMS_Trk2, fTrk1XOffset, fTrk2XOffset;
0086   double fHitSigmaX, fHitSigmaY, fHitSigmaZ;
0087   std::vector<double> fToFCellWidth;
0088   double fToFCellHeight, fToFPitchX, fToFPitchY;
0089   int fToFNCellX, fToFNCellY;
0090   double fToFInsertion, fBeamXRMS_ToF, fToFXOffset, fTimeSigma;
0091 };
0092
0093 CTPPSRecHitProducer::CTPPSRecHitProducer(const edm::ParameterSet& iConfig) {
0094   produces<edm::CTPPSFastRecHitContainer>("CTPPSFastRecHits");
0095
0096   //Name of Collection use for create the XF
0097   std::string mix_ = iConfig.getParameter<std::string>("mixLabel");
0098   std::string collection_for_XF = iConfig.getParameter<std::string>("InputCollection");
0099   cf_token = consumes<CrossingFrame<PSimHit> >(edm::InputTag(mix_, collection_for_XF));
0100
0101   // Read the detector parameters
0102   fTrackerWidth = iConfig.getParameter<double>("TrackerWidth");
0103   fTrackerHeight = iConfig.getParameter<double>("TrackerHeight");
0104   fTrackerInsertion = iConfig.getParameter<double>("TrackerInsertion");
0105   fBeamXRMS_Trk1 = iConfig.getParameter<double>("BeamXRMS_Trk1");
0106   fBeamXRMS_Trk2 = iConfig.getParameter<double>("BeamXRMS_Trk2");
0107   fTrk1XOffset = iConfig.getParameter<double>("Trk1XOffset");
0108   fTrk2XOffset = iConfig.getParameter<double>("Trk2XOffset");
0109   fHitSigmaX = iConfig.getParameter<double>("HitSigmaX");
0110   fHitSigmaY = iConfig.getParameter<double>("HitSigmaY");
0111   fHitSigmaZ = iConfig.getParameter<double>("HitSigmaZ");
0112   fToFCellWidth = iConfig.getUntrackedParameter<std::vector<double> >("ToFCellWidth");
0113   fToFCellHeight = iConfig.getParameter<double>("ToFCellHeight");
0114   fToFPitchX = iConfig.getParameter<double>("ToFPitchX");
0115   fToFPitchY = iConfig.getParameter<double>("ToFPitchY");
0116   fToFNCellX = iConfig.getParameter<int>("ToFNCellX");
0117   fToFNCellY = iConfig.getParameter<int>("ToFNCellY");
0118   fToFInsertion = iConfig.getParameter<double>("ToFInsertion");
0119   fBeamXRMS_ToF = iConfig.getParameter<double>("BeamXRMS_ToF");
0120   fToFXOffset = iConfig.getParameter<double>("ToFXOffset");
0121   fTimeSigma = iConfig.getParameter<double>("TimeSigma");
0122
0123   edm::Service<edm::RandomNumberGenerator> rng;
0124   if (!rng.isAvailable()) {
0125     throw cms::Exception("Configuration")
0126         << "CTPPSRecHitProducer requires the RandomNumberGeneratorService\n"
0127            "which is not present in the configuration file.  You must add the service\n"
0128            "in the configuration file or remove the modules that require it.";
0129   }
0130 }
0131
0132 CTPPSRecHitProducer::~CTPPSRecHitProducer() {}
0133
0134 void CTPPSRecHitProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0135   using namespace edm;
0136
0137   edm::Service<edm::RandomNumberGenerator> rng;
0138   CLHEP::HepRandomEngine* engine = &rng->getEngine(iEvent.streamID());
0139   if (engine->name() != "TRandom3") {
0140     throw cms::Exception("Configuration") << "The TRandom3 engine type must be used with CTPPSRecHitProducer, Random "
0141                                              "Number Generator Service not correctly configured!";
0142   }
0143
0144   //Mix the PSimHits for pileup
0145   Handle<CrossingFrame<PSimHit> > xFrame;
0146   iEvent.getByToken(cf_token, xFrame);
0147
0148   std::unique_ptr<MixCollection<PSimHit> > simHits(new MixCollection<PSimHit>(xFrame.product()));
0149
0150   std::vector<CTPPSFastRecHit> theCTPPSFastRecHit;
0151
0152   for (MixCollection<PSimHit>::MixItr ihit = simHits->begin(); ihit != simHits->end(); ihit++)
0153
0154   {
0155     const PSimHit* simHit = &(*ihit);
0156     unsigned int detlayerId = simHit->detUnitId();
0157     // DetId layer codification for PSimHit
0158     // 2014314496 -> Tracker1 zPositive
0159     // 2014838784 -> Tracker2 zPositive
0160     // 2046820352 -> Timing   zPositive
0161     // 2031091712 -> Tracker1 zNegative
0162     // 2031616000 -> Tracker2 zNegative
0163     // 2063597568 -> Timing   zNegative
0164     bool outside = false;
0165
0166     //Make Tracker RecHits by smearing the SimHits positions and applying fiducial cuts
0167     CTPPSFastRecHit rechit;  // (const Local3DPoint& entry, unsigned int detId, float tof, unsigned int cellId)
0168     if (detlayerId == 2014314496 || detlayerId == 2031091712) {
0169       // Apply position smearing
0170       float x_trk1 = simHit->entryPoint().x() + CLHEP::RandGauss::shoot(engine, 0, fHitSigmaX * um_to_mm);
0171       float y_trk1 = simHit->entryPoint().y() + CLHEP::RandGauss::shoot(engine, 0, fHitSigmaY * um_to_mm);
0172       float z_trk1 = simHit->entryPoint().z();
0173
0174       // Apply fiducial cuts
0175       double pos_trk1 = fTrackerInsertion * fBeamXRMS_Trk1 + fTrk1XOffset;
0176       if (x_trk1 > 0 || fabs(x_trk1) < pos_trk1 || fabs(x_trk1) > (fTrackerWidth + pos_trk1) ||
0177           fabs(y_trk1) > fTrackerHeight / 2.)
0178         outside = true;
0179
0180       float tof = 0.0;
0181       unsigned int cellId = 1;
0182       Local3DPoint xyzzy = Local3DPoint(x_trk1, y_trk1, z_trk1);
0183       if (!outside) {
0184         rechit.setLocal3DPoint(xyzzy);
0185         rechit.setTof(tof);
0186         rechit.setDetUnitId(detlayerId);
0187         rechit.setCellId(cellId);
0188         theCTPPSFastRecHit.push_back(rechit);
0189       }
0190     }
0191     if (detlayerId == 2014838784 || detlayerId == 2031616000) {
0192       // Apply position smearing
0193       double x_trk2 = simHit->entryPoint().x() + CLHEP::RandGauss::shoot(engine, 0, fHitSigmaX * um_to_mm);
0194       double y_trk2 = simHit->entryPoint().y() + CLHEP::RandGauss::shoot(engine, 0, fHitSigmaY * um_to_mm);
0195       double z_trk2 = simHit->entryPoint().z();
0196       // Apply fiducial cuts
0197       double pos_trk2 = fTrackerInsertion * fBeamXRMS_Trk2 + fTrk2XOffset;
0198       if (x_trk2 > 0 || fabs(x_trk2) < pos_trk2 || fabs(x_trk2) > (fTrackerWidth + pos_trk2) ||
0199           fabs(y_trk2) > fTrackerHeight / 2.)
0200         outside = true;
0201
0202       float tof = 0.0;
0203       unsigned int cellId = 2;
0204       Local3DPoint xyzzy = Local3DPoint(x_trk2, y_trk2, z_trk2);
0205       if (!outside) {
0206         rechit.setLocal3DPoint(xyzzy);
0207         rechit.setTof(tof);
0208         rechit.setDetUnitId(detlayerId);
0209         rechit.setCellId(cellId);
0210         theCTPPSFastRecHit.push_back(rechit);
0211       }
0212     }
0213     //Make Timing RecHits by smearing the SimHits time of flight and checking the cell of the hit
0214     //The RecHit position is the centre of the cell
0215     if (detlayerId == 2046820352 || detlayerId == 2063597568) {
0216       float t = simHit->tof();
0217       unsigned int cellId = 0;
0218       float tof = CLHEP::RandGauss::shoot(engine, t, fTimeSigma);
0219       double x_tof = simHit->entryPoint().x();
0220       double y_tof = simHit->entryPoint().y();
0221
0222       double pos_tof = fToFInsertion * fBeamXRMS_ToF + fToFXOffset;
0223
0224       std::vector<double> vToFCellWidth;
0225       vToFCellWidth.reserve(8);
0226       for (int i = 0; i < 8; i++) {
0227         vToFCellWidth.push_back(fToFCellWidth[i]);
0228       }
0229       CTPPSToFDetector* ToFDet = new CTPPSToFDetector(
0230           fToFNCellX, fToFNCellY, vToFCellWidth, fToFCellHeight, fToFPitchX, fToFPitchY, pos_tof, fTimeSigma);
0231       cellId = ToFDet->findCellId(x_tof, y_tof);
0232       if (cellId > 0) {
0233         double xc_tof = 0., yc_tof = 0.;
0234         ToFDet->get_CellCenter(cellId, xc_tof, yc_tof);
0235         Local3DPoint xyzzy = Local3DPoint(xc_tof, yc_tof, simHit->entryPoint().z());
0236         rechit.setLocal3DPoint(xyzzy);
0237         rechit.setTof(tof);
0238         rechit.setDetUnitId(detlayerId);
0239         rechit.setCellId(cellId);
0240         theCTPPSFastRecHit.push_back(rechit);
0241       }
0242     }
0243   }
0244
0245   std::unique_ptr<CTPPSFastRecHitContainer> output_recHits(new edm::CTPPSFastRecHitContainer);
0246   output_recHits->reserve(simHits->size());
0247
0248   for (std::vector<CTPPSFastRecHit>::const_iterator i = theCTPPSFastRecHit.begin(); i != theCTPPSFastRecHit.end();
0249        i++) {
0250     output_recHits->push_back(*i);
0251   }
0252
0253   iEvent.put(std::move(output_recHits), "CTPPSFastRecHits");
0254 }
0255
0256 // ------------ method called once each stream before processing any runs, lumis or events  ------------
0257 void CTPPSRecHitProducer::beginStream(edm::StreamID) {}
0258
0259 // ------------ method called once each stream after processing all runs, lumis and events  ------------
0260 void CTPPSRecHitProducer::endStream() {}
0261
0262 //define this as a plug-in
0263 DEFINE_FWK_MODULE(CTPPSRecHitProducer);