CTPPSSimHitProducer/plugins/CTPPSSimHitProducer.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    FastSimulation/CTPPSSimHitProducer
0004 // Class:      CTPPSSimHitProducer
0005 //
0006 /**\class CTPPSSimHitProducer CTPPSSimHitProducer.cc FastSimulation/CTPPSSimHitProducer/plugins/CTPPSSimHitProducer.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, 05 Sep 2016 18:49:10 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 // SimHitContainer
0033 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0034
0035 // STL headers
0036 #include <vector>
0037 #include <iostream>
0038
0039 // HepMC headers
0040 #include "HepMC/GenEvent.h"
0041 #include "HepMC/GenVertex.h"
0042 #include "HepMC/GenParticle.h"
0043 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0044 #include "SimDataFormats/Forward/interface/LHCTransportLink.h"
0045
0046 #include "DataFormats/Math/interface/LorentzVector.h"
0047
0048 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0049 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0050
0051 #include "Utilities/PPS/interface/PPSUnitConversion.h"
0052
0053 //
0054 // class declaration
0055 //
0056
0057 class CTPPSSimHitProducer : public edm::stream::EDProducer<> {
0058 public:
0059   explicit CTPPSSimHitProducer(const edm::ParameterSet&);
0060   ~CTPPSSimHitProducer() override;
0061
0062 private:
0063   void beginStream(edm::StreamID) override;
0064   void produce(edm::Event&, const edm::EventSetup&) override;
0065   void endStream() override;
0066
0067   edm::EDGetTokenT<edm::HepMCProduct> mcEventToken;  // label of MC event
0068   edm::Handle<edm::HepMCProduct> EvtHandle;
0069   // ----------member data ---------------------------
0070   double fz_tracker1, fz_tracker2, fz_timing;
0071 };
0072
0073 //
0074 // constants, enums and typedefs
0075 //
0076
0077 //
0078 // static data member definitions
0079 //
0080
0081 //
0082 // constructors and destructor
0083 //
0084 CTPPSSimHitProducer::CTPPSSimHitProducer(const edm::ParameterSet& iConfig) {
0085   produces<edm::PSimHitContainer>("CTPPSHits");
0086   // consumes
0087   mcEventToken =
0088       mayConsume<edm::HepMCProduct>(iConfig.getUntrackedParameter<edm::InputTag>("MCEvent", std::string("")));
0089
0090   // Read the position of the trackers and timing detectors
0091   fz_tracker1 = iConfig.getParameter<double>("Z_Tracker1");
0092   fz_tracker2 = iConfig.getParameter<double>("Z_Tracker2");
0093   fz_timing = iConfig.getParameter<double>("Z_Timing");
0094 }
0095
0096 CTPPSSimHitProducer::~CTPPSSimHitProducer() {}
0097
0098 //
0099 // member functions
0100 //
0101
0102 // ------------ method called to produce the data  ------------
0103 void CTPPSSimHitProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0104   using namespace edm;
0105
0106   std::vector<PSimHit> theCTPPSHits;
0107   iEvent.getByToken(mcEventToken, EvtHandle);
0108
0109   const HepMC::GenEvent* Evt = EvtHandle->GetEvent();
0110   std::vector<math::XYZTLorentzVector> protonCTPPS;
0111   protonCTPPS.clear();
0112   for (HepMC::GenEvent::vertex_const_iterator ivtx = Evt->vertices_begin(); ivtx != Evt->vertices_end(); ivtx++) {
0113     if ((*ivtx)->id() != 0)
0114       continue;
0115     double prim_vtxZ = (*ivtx)->position().z() * mm_to_m;  //in meters
0116     // Get the vertices at the entrance of CTPPS and get the protons coming out of them (propagated by Hector)
0117     for (HepMC::GenVertex::particles_out_const_iterator i = (*ivtx)->particles_out_const_begin();
0118          i != (*ivtx)->particles_out_const_end();
0119          ++i) {
0120       int pid = (*i)->pdg_id();
0121       if (pid != 2212)
0122         continue;
0123
0124       HepMC::GenVertex* pv = (*i)->production_vertex();
0125       const HepMC::FourVector& vertex = pv->position();
0126       const HepMC::FourVector p((*i)->momentum());
0127       protonCTPPS.push_back(math::XYZTLorentzVector(p.x(), p.y(), p.z(), p.t()));
0128
0129       LocalPoint initialPosition_tr1, initialPosition_tr2, initialPosition_tof;
0130
0131       double t0_tr1 = 0., t0_tr2 = 0., t0_tof = 0.;
0132       //  t0_* has dimensions of mm
0133       //  Convert to ns for internal calculations.
0134       const double c_light_s = 2.99792458e+11;  // mm/s;
0135       const double s_to_ns = 1.e9;
0136       const double m_to_mm = 1.e3;
0137       double x_tr1 = 0., x_tr2 = 0., x_tof = 0., y_tr1 = 0., y_tr2 = 0., y_tof = 0., z_tr1 = 0.;
0138       double z_tr2 = fz_tracker2;  //m
0139       double z_tof = fz_timing;    //m
0140       int Direction = 0;
0141
0142       // Read the vertex made by SimTransport at the entrance of det1
0143       if (std::abs(vertex.eta()) > 8. && (*i)->status() == 1) {
0144         if (vertex.z() > 0)
0145           Direction = 1;
0146         else if (vertex.z() < 0)
0147           Direction = -1;
0148
0149         //Get the global coordinates at Tracker1, equal to those of the vertex at CTPPS
0150         x_tr1 = vertex.x();
0151         y_tr1 = vertex.y();
0152         z_tr1 = vertex.z() * mm_to_m;  //move z from mm to meters
0153         Local3DPoint xyzzy_tr1(x_tr1, y_tr1, z_tr1);
0154         initialPosition_tr1 = xyzzy_tr1;
0155         t0_tr1 = vertex.t() / c_light_s * s_to_ns;
0156         //Get the global coordinates at Tracker2 by propagating as a straight line from Tracker1
0157         t0_tr2 = z_tr2 * m_to_mm / c_light_s * s_to_ns;  //discuss latter if needs to be corrected with vertex position
0158         t0_tr2 = t0_tr1 + (z_tr2 - z_tr1) * m_to_mm / c_light_s * s_to_ns;  //corrected with vertex position
0159         z_tr2 *= Direction;
0160         x_tr2 = x_tr1 + (p.x() / p.z()) * (z_tr2 - z_tr1) * m_to_mm;
0161         y_tr2 = y_tr1 + (p.y() / p.z()) * (z_tr2 - z_tr1) * m_to_mm;
0162         Local3DPoint xyzzy_tr2(x_tr2, y_tr2, z_tr2);
0163         initialPosition_tr2 = xyzzy_tr2;
0164         //Propagate as a straight line from Tracker1
0165         t0_tof = z_tof * m_to_mm / c_light_s * s_to_ns;  //discuss latter if needs to be corrected with vertex position
0166         t0_tof = (z_tof - prim_vtxZ) * m_to_mm / c_light_s * s_to_ns;  //corrected with vertex position
0167         z_tof *= Direction;
0168         x_tof = x_tr1 + (p.x() / p.z()) * (z_tof - z_tr1) * m_to_mm;
0169         y_tof = y_tr1 + (p.y() / p.z()) * (z_tof - z_tr1) * m_to_mm;
0170         Local3DPoint xyzzy_tof(x_tof, y_tof, z_tof);
0171         initialPosition_tof = xyzzy_tof;
0172
0173         // DetId codification for PSimHit from CTPPSPixel- It will be replaced by CTPPSDetId
0174         // 2014314496 -> Tracker1 zPositive
0175         // 2014838784 -> Tracker2 zPositive
0176         // 2046820352 -> Timing   zPositive
0177         // 2031091712 -> Tracker1 zNegative
0178         // 2031616000 -> Tracker2 zNegative
0179         // 2063597568 -> Timing   zNegative
0180
0181         if (Direction > 0.) {
0182           PSimHit hit_tr1(xyzzy_tr1, xyzzy_tr1, 0., t0_tr1, 0., pid, 2014314496, 0, 0., 0., 2);
0183           PSimHit hit_tr2(xyzzy_tr2, xyzzy_tr2, 0., t0_tr2, 0., pid, 2014838784, 0, 0., 0., 2);
0184           PSimHit hit_tof(xyzzy_tof, xyzzy_tof, 0., t0_tof, 0., pid, 2046820352, 0, 0., 0., 2);
0185           theCTPPSHits.push_back(hit_tr1);
0186           theCTPPSHits.push_back(hit_tr2);
0187           theCTPPSHits.push_back(hit_tof);
0188         }
0189         if (Direction < 0.) {
0190           PSimHit hit_tr1(xyzzy_tr1, xyzzy_tr1, 0., t0_tr1, 0., pid, 2031091712, 0, 0., 0., 2);
0191           PSimHit hit_tr2(xyzzy_tr2, xyzzy_tr2, 0., t0_tr2, 0., pid, 2031616000, 0, 0., 0., 2);
0192           PSimHit hit_tof(xyzzy_tof, xyzzy_tof, 0., t0_tof, 0., pid, 2063597568, 0, 0., 0., 2);
0193           theCTPPSHits.push_back(hit_tr1);
0194           theCTPPSHits.push_back(hit_tr2);
0195           theCTPPSHits.push_back(hit_tof);
0196         }
0197       }
0198     }
0199   }
0200   std::unique_ptr<edm::PSimHitContainer> pctpps(new edm::PSimHitContainer);
0201   for (std::vector<PSimHit>::const_iterator i = theCTPPSHits.begin(); i != theCTPPSHits.end(); i++) {
0202     pctpps->push_back(*i);
0203   }
0204   iEvent.put(std::move(pctpps), "CTPPSHits");
0205 }
0206
0207 // ------------ method called once each stream before processing any runs, lumis or events  ------------
0208 void CTPPSSimHitProducer::beginStream(edm::StreamID) {}
0209
0210 // ------------ method called once each stream after processing all runs, lumis and events  ------------
0211 void CTPPSSimHitProducer::endStream() {}
0212
0213 //define this as a plug-in
0214 DEFINE_FWK_MODULE(CTPPSSimHitProducer);