CTPPSSimHitProducer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    FastSimulation/CTPPSSimHitProducer
// Class:      CTPPSSimHitProducer
// 
//
// Original Author:  Dilson De Jesus Damiao
//         Created:  Mon, 05 Sep 2016 18:49:10 GMT
//
//

// system include files
#include <memory>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/StreamID.h"

// SimHitContainer
#include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"

// STL headers 
#include <vector>
#include <iostream>

// HepMC headers
#include "HepMC/GenEvent.h"
#include "HepMC/GenVertex.h"
#include "HepMC/GenParticle.h"
#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/Forward/interface/LHCTransportLink.h"

#include "DataFormats/Math/interface/LorentzVector.h"

#include "DataFormats/GeometryVector/interface/LocalPoint.h"
#include "DataFormats/GeometryVector/interface/LocalVector.h"

#include "Utilities/PPS/interface/PPSUnitConversion.h"

//
// class declaration
//

class CTPPSSimHitProducer : public edm::stream::EDProducer<> {
    public:
        explicit CTPPSSimHitProducer(const edm::ParameterSet&);
        ~CTPPSSimHitProducer() override;

    private:
        void beginStream(edm::StreamID) override;
        void produce(edm::Event&, const edm::EventSetup&) override;
        void endStream() override;

        edm::EDGetTokenT< edm::HepMCProduct > mcEventToken; // label of MC event
        edm::Handle< edm::HepMCProduct > EvtHandle ;
        // ----------member data ---------------------------
    double fz_tracker1, fz_tracker2, fz_timing;
 


};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
CTPPSSimHitProducer::CTPPSSimHitProducer(const edm::ParameterSet& iConfig)
{

    produces<edm::PSimHitContainer>("CTPPSHits");
    // consumes
    mcEventToken = mayConsume<edm::HepMCProduct>(iConfig.getUntrackedParameter<edm::InputTag>("MCEvent",std::string("")));

    // Read the position of the trackers and timing detectors
    fz_tracker1   = iConfig.getParameter<double>("Z_Tracker1");
    fz_tracker2   = iConfig.getParameter<double>("Z_Tracker2"); 
    fz_timing =  iConfig.getParameter<double>("Z_Timing");  
}

CTPPSSimHitProducer::~CTPPSSimHitProducer()
{

}


//
// member functions
//

// ------------ method called to produce the data  ------------
    void
CTPPSSimHitProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
{
    using namespace edm;

    std::vector<PSimHit> theCTPPSHits;
    iEvent.getByToken( mcEventToken , EvtHandle );

    const HepMC::GenEvent* Evt = EvtHandle->GetEvent() ;
    std::vector<math::XYZTLorentzVector> protonCTPPS;
    protonCTPPS.clear();
    for(HepMC::GenEvent::vertex_const_iterator ivtx = Evt->vertices_begin();ivtx!=Evt->vertices_end();ivtx++) {
        if ((*ivtx)->id()!=0) continue;
        double prim_vtxZ=(*ivtx)->position().z()*mm_to_m; //in meters
        // Get the vertices at the entrance of CTPPS and get the protons coming out of them (propagated by Hector)
        for(HepMC::GenVertex::particles_out_const_iterator i=(*ivtx)->particles_out_const_begin(); i != (*ivtx)->particles_out_const_end();++i) {
            int pid = (*i)->pdg_id();
            if(pid!=2212) continue;

            HepMC::GenVertex* pv = (*i)->production_vertex();
            const HepMC::FourVector& vertex = pv->position();
            const HepMC::FourVector p((*i)->momentum());
            protonCTPPS.push_back(math::XYZTLorentzVector(p.x(),p.y(),p.z(),p.t()));

            LocalPoint initialPosition_tr1, initialPosition_tr2, initialPosition_tof;

            double t0_tr1 = 0., t0_tr2 = 0., t0_tof = 0.;
            //  t0_* has dimensions of mm
            //  Convert to ns for internal calculations.
            const double c_light_s = 2.99792458e+11;// mm/s;
            const double s_to_ns = 1.e9;
            const double m_to_mm = 1.e3;
            double x_tr1 = 0. , x_tr2 = 0. , x_tof = 0., y_tr1 = 0. , y_tr2 = 0. , y_tof = 0., z_tr1 = 0.; 
            double z_tr2 = fz_tracker2; //m
            double z_tof = fz_timing; //m
            int Direction =0;

            // Read the vertex made by SimTransport at the entrance of det1
            if(std::abs(vertex.eta())>8. && (*i)->status() == 1){
                if (vertex.z()>0) Direction = 1; 
                else if (vertex.z() < 0 ) Direction = -1;

                //Get the global coordinates at Tracker1, equal to those of the vertex at CTPPS
                x_tr1 = vertex.x(); 
                y_tr1 = vertex.y();
                z_tr1 =  vertex.z()*mm_to_m;//move z from mm to meters 
                Local3DPoint xyzzy_tr1(x_tr1, y_tr1, z_tr1);
                initialPosition_tr1 = xyzzy_tr1;
                t0_tr1 = vertex.t()/c_light_s*s_to_ns; 
                //Get the global coordinates at Tracker2 by propagating as a straight line from Tracker1
                t0_tr2 = z_tr2*m_to_mm/c_light_s*s_to_ns;//discuss latter if needs to be corrected with vertex position
                t0_tr2 = t0_tr1+(z_tr2-z_tr1)*m_to_mm/c_light_s*s_to_ns;//corrected with vertex position
                z_tr2 *= Direction; 
                x_tr2 = x_tr1 + (p.x()/p.z())*(z_tr2-z_tr1)*m_to_mm;
                y_tr2 = y_tr1 + (p.y()/p.z())*(z_tr2-z_tr1)*m_to_mm;
                Local3DPoint xyzzy_tr2(x_tr2, y_tr2, z_tr2);
                initialPosition_tr2 = xyzzy_tr2;
                //Propagate as a straight line from Tracker1
                t0_tof = z_tof*m_to_mm/c_light_s*s_to_ns;//discuss latter if needs to be corrected with vertex position
                t0_tof = (z_tof-prim_vtxZ)*m_to_mm/c_light_s*s_to_ns;//corrected with vertex position
                z_tof *= Direction; 
                x_tof = x_tr1 + (p.x()/p.z())*(z_tof-z_tr1)*m_to_mm;
                y_tof = y_tr1 + (p.y()/p.z())*(z_tof-z_tr1)*m_to_mm;
                Local3DPoint xyzzy_tof(x_tof, y_tof, z_tof);
                initialPosition_tof = xyzzy_tof;

                // DetId codification for PSimHit from CTPPSPixel- It will be replaced by CTPPSDetId
                // 2014314496 -> Tracker1 zPositive
                // 2014838784 -> Tracker2 zPositive
                // 2046820352 -> Timing   zPositive
                // 2031091712 -> Tracker1 zNegative
                // 2031616000 -> Tracker2 zNegative
                // 2063597568 -> Timing   zNegative

                if(Direction > 0.)  {
                    PSimHit hit_tr1(xyzzy_tr1,xyzzy_tr1,0.,t0_tr1,0.,pid,2014314496,0,0.,0.,2);
                    PSimHit hit_tr2(xyzzy_tr2,xyzzy_tr2,0.,t0_tr2,0.,pid,2014838784,0,0.,0.,2);
                    PSimHit hit_tof(xyzzy_tof,xyzzy_tof,0.,t0_tof,0.,pid,2046820352,0,0.,0.,2);
                    theCTPPSHits.push_back(hit_tr1);
                    theCTPPSHits.push_back(hit_tr2);
                    theCTPPSHits.push_back(hit_tof);
                }
                if(Direction < 0.)  {
                    PSimHit hit_tr1(xyzzy_tr1,xyzzy_tr1,0.,t0_tr1,0.,pid,2031091712,0,0.,0.,2);
                    PSimHit hit_tr2(xyzzy_tr2,xyzzy_tr2,0.,t0_tr2,0.,pid,2031616000,0,0.,0.,2);
                    PSimHit hit_tof(xyzzy_tof,xyzzy_tof,0.,t0_tof,0.,pid,2063597568,0,0.,0.,2);
                    theCTPPSHits.push_back(hit_tr1);
                    theCTPPSHits.push_back(hit_tr2);
                    theCTPPSHits.push_back(hit_tof);
                }
            }
        }
    }
    std::unique_ptr<edm::PSimHitContainer> pctpps(new edm::PSimHitContainer);
    int  n = 0;
    for ( std::vector<PSimHit>::const_iterator i = theCTPPSHits.begin();
            i != theCTPPSHits.end(); i++ ) {
        pctpps->push_back(*i);
        n += 1;
    }
    iEvent.put(std::move(pctpps),"CTPPSHits");
}

// ------------ method called once each stream before processing any runs, lumis or events  ------------
    void
CTPPSSimHitProducer::beginStream(edm::StreamID)
{
}

// ------------ method called once each stream after processing all runs, lumis and events  ------------
void
CTPPSSimHitProducer::endStream() {
}


//define this as a plug-in
DEFINE_FWK_MODULE(CTPPSSimHitProducer);