Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:30:52 

0001 /****************************************************************************
0002  * Authors:
0003  *   Jan Kašpar
0004  *   Laurent Forthomme
0005  ****************************************************************************/
0006 
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/stream/EDProducer.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 #include "FWCore/Framework/interface/ESWatcher.h"
0013 #include "FWCore/Framework/interface/MakerMacros.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 
0016 #include "FWCore/Utilities/interface/ESInputTag.h"
0017 #include "FWCore/Utilities/interface/Exception.h"
0018 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0019 
0020 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0021 #include "DataFormats/CTPPSDetId/interface/CTPPSPixelDetId.h"
0022 #include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h"
0023 
0024 #include "DataFormats/Common/interface/DetSetVector.h"
0025 #include "DataFormats/CTPPSReco/interface/TotemRPRecHit.h"
0026 #include "DataFormats/CTPPSReco/interface/CTPPSDiamondRecHit.h"
0027 #include "DataFormats/CTPPSReco/interface/CTPPSPixelRecHit.h"
0028 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
0029 
0030 #include "CondFormats/DataRecord/interface/CTPPSInterpolatedOpticsRcd.h"
0031 #include "CondFormats/PPSObjects/interface/LHCInterpolatedOpticalFunctionsSetCollection.h"
0032 
0033 #include "CondFormats/PPSObjects/interface/CTPPSBeamParameters.h"
0034 #include "CondFormats/DataRecord/interface/CTPPSBeamParametersRcd.h"
0035 
0036 #include "CondFormats/RunInfo/interface/LHCInfo.h"
0037 #include "CondFormats/DataRecord/interface/LHCInfoRcd.h"
0038 
0039 #include "CondFormats/DataRecord/interface/PPSDirectSimulationDataRcd.h"
0040 #include "CondFormats/PPSObjects/interface/PPSDirectSimulationData.h"
0041 
0042 #include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h"
0043 #include "Geometry/Records/interface/VeryForwardMisalignedGeometryRecord.h"
0044 #include "Geometry/VeryForwardRPTopology/interface/RPTopology.h"
0045 #include "CondFormats/PPSObjects/interface/PPSPixelTopology.h"
0046 #include "CondFormats/DataRecord/interface/PPSPixelTopologyRcd.h"
0047 
0048 #include "FWCore/ServiceRegistry/interface/Service.h"
0049 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0050 
0051 #include "CLHEP/Random/RandGauss.h"
0052 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0053 
0054 #include <unordered_map>
0055 
0056 #include "TMath.h"
0057 #include "TMatrixD.h"
0058 #include "TVectorD.h"
0059 #include "TF1.h"
0060 #include "TF2.h"
0061 #include "TFile.h"
0062 #include "CLHEP/Random/RandFlat.h"
0063 
0064 //----------------------------------------------------------------------------------------------------
0065 
0066 class PPSDirectProtonSimulation : public edm::stream::EDProducer<> {
0067 public:
0068   explicit PPSDirectProtonSimulation(const edm::ParameterSet &);
0069   ~PPSDirectProtonSimulation() override {}
0070 
0071   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0072 
0073 private:
0074   void produce(edm::Event &, const edm::EventSetup &) override;
0075 
0076   void processProton(const HepMC::GenVertex *in_vtx,
0077                      const HepMC::GenParticle *in_trk,
0078                      const CTPPSGeometry &geometry,
0079                      const LHCInfo &lhcInfo,
0080                      const CTPPSBeamParameters &beamParameters,
0081                      const PPSPixelTopology &ppt,
0082                      const LHCInterpolatedOpticalFunctionsSetCollection &opticalFunctions,
0083                      CLHEP::HepRandomEngine *rndEngine,
0084 
0085                      std::vector<CTPPSLocalTrackLite> &out_tracks,
0086 
0087                      edm::DetSetVector<TotemRPRecHit> &out_strip_hits,
0088                      edm::DetSetVector<CTPPSPixelRecHit> &out_pixel_hits,
0089                      edm::DetSetVector<CTPPSDiamondRecHit> &out_diamond_hits,
0090 
0091                      std::map<int, edm::DetSetVector<TotemRPRecHit>> &out_strip_hits_per_particle,
0092                      std::map<int, edm::DetSetVector<CTPPSPixelRecHit>> &out_pixel_hits_per_particle,
0093                      std::map<int, edm::DetSetVector<CTPPSDiamondRecHit>> &out_diamond_hits_per_particle) const;
0094 
0095   // ------------ config file parameters ------------
0096 
0097   // conditions
0098   edm::ESGetToken<LHCInfo, LHCInfoRcd> tokenLHCInfo_;
0099   edm::ESGetToken<CTPPSBeamParameters, CTPPSBeamParametersRcd> tokenBeamParameters_;
0100   edm::ESGetToken<PPSPixelTopology, PPSPixelTopologyRcd> pixelTopologyToken_;
0101   edm::ESGetToken<LHCInterpolatedOpticalFunctionsSetCollection, CTPPSInterpolatedOpticsRcd> tokenOpticalFunctions_;
0102   edm::ESGetToken<CTPPSGeometry, VeryForwardMisalignedGeometryRecord> tokenGeometry_;
0103   edm::ESGetToken<PPSDirectSimulationData, PPSDirectSimulationDataRcd> tokenDirectSimuData_;
0104 
0105   edm::ESWatcher<PPSDirectSimulationDataRcd> directSimuDataRcdWatcher_;
0106 
0107   // input
0108   edm::EDGetTokenT<edm::HepMCProduct> hepMCToken_;
0109 
0110   // flags what output to be produced
0111   bool produceScoringPlaneHits_;
0112   bool produceRecHits_;
0113 
0114   // settings of LHC aperture limitations (high xi)
0115   bool useEmpiricalApertures_;
0116   std::unique_ptr<TF2> empiricalAperture45_;
0117   std::unique_ptr<TF2> empiricalAperture56_;
0118 
0119   // efficiency flags
0120   bool useTrackingEfficiencyPerRP_;
0121   bool useTimingEfficiencyPerRP_;
0122   bool useTrackingEfficiencyPerPlane_;
0123   bool useTimingEfficiencyPerPlane_;
0124 
0125   // efficiency maps
0126   std::map<unsigned int, std::unique_ptr<TH2F>> efficiencyMapsPerRP_;
0127   std::map<unsigned int, std::unique_ptr<TH2F>> efficiencyMapsPerPlane_;
0128 
0129   // other parameters
0130   bool produceHitsRelativeToBeam_;
0131   bool roundToPitch_;
0132   bool checkIsHit_;
0133 
0134   double pitchStrips_;              ///< strip pitch in mm
0135   double insensitiveMarginStrips_;  ///< size of insensitive margin at sensor's edge facing the beam, in mm
0136 
0137   double pitchPixelsHor_;
0138   double pitchPixelsVer_;
0139 
0140   unsigned int verbosity_;
0141 
0142   std::unique_ptr<TF1> timeResolutionDiamonds45_, timeResolutionDiamonds56_;
0143 
0144   // ------------ internal parameters ------------
0145 
0146   /// internal variable: v position of strip 0, in mm
0147   double stripZeroPosition_;
0148 };
0149 
0150 //----------------------------------------------------------------------------------------------------
0151 
0152 PPSDirectProtonSimulation::PPSDirectProtonSimulation(const edm::ParameterSet &iConfig)
0153     : tokenLHCInfo_(esConsumes(edm::ESInputTag{"", iConfig.getParameter<std::string>("lhcInfoLabel")})),
0154       tokenBeamParameters_(esConsumes()),
0155       pixelTopologyToken_(esConsumes()),
0156       tokenOpticalFunctions_(esConsumes(edm::ESInputTag{"", iConfig.getParameter<std::string>("opticsLabel")})),
0157       tokenGeometry_(esConsumes()),
0158       tokenDirectSimuData_(esConsumes()),
0159 
0160       hepMCToken_(consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("hepMCTag"))),
0161 
0162       produceScoringPlaneHits_(iConfig.getParameter<bool>("produceScoringPlaneHits")),
0163       produceRecHits_(iConfig.getParameter<bool>("produceRecHits")),
0164 
0165       useEmpiricalApertures_(iConfig.getParameter<bool>("useEmpiricalApertures")),
0166 
0167       useTrackingEfficiencyPerRP_(iConfig.getParameter<bool>("useTrackingEfficiencyPerRP")),
0168       useTimingEfficiencyPerRP_(iConfig.getParameter<bool>("useTimingEfficiencyPerRP")),
0169       useTrackingEfficiencyPerPlane_(iConfig.getParameter<bool>("useTrackingEfficiencyPerPlane")),
0170       useTimingEfficiencyPerPlane_(iConfig.getParameter<bool>("useTimingEfficiencyPerPlane")),
0171 
0172       produceHitsRelativeToBeam_(iConfig.getParameter<bool>("produceHitsRelativeToBeam")),
0173       roundToPitch_(iConfig.getParameter<bool>("roundToPitch")),
0174       checkIsHit_(iConfig.getParameter<bool>("checkIsHit")),
0175 
0176       pitchStrips_(iConfig.getParameter<double>("pitchStrips")),
0177       insensitiveMarginStrips_(iConfig.getParameter<double>("insensitiveMarginStrips")),
0178 
0179       pitchPixelsHor_(iConfig.getParameter<double>("pitchPixelsHor")),
0180       pitchPixelsVer_(iConfig.getParameter<double>("pitchPixelsVer")),
0181 
0182       verbosity_(iConfig.getUntrackedParameter<unsigned int>("verbosity", 0)) {
0183   if (produceScoringPlaneHits_)
0184     produces<std::vector<CTPPSLocalTrackLite>>();
0185 
0186   if (produceRecHits_) {
0187     produces<edm::DetSetVector<TotemRPRecHit>>();
0188     produces<edm::DetSetVector<CTPPSDiamondRecHit>>();
0189     produces<edm::DetSetVector<CTPPSPixelRecHit>>();
0190 
0191     produces<std::map<int, edm::DetSetVector<TotemRPRecHit>>>();
0192     produces<std::map<int, edm::DetSetVector<CTPPSDiamondRecHit>>>();
0193     produces<std::map<int, edm::DetSetVector<CTPPSPixelRecHit>>>();
0194   }
0195 
0196   // check user input
0197   if (useTrackingEfficiencyPerRP_ && useTrackingEfficiencyPerPlane_)
0198     throw cms::Exception("PPS")
0199         << "useTrackingEfficiencyPerRP and useTrackingEfficiencyPerPlane should not be simultaneously set true.";
0200 
0201   if (useTimingEfficiencyPerRP_ && useTimingEfficiencyPerPlane_)
0202     throw cms::Exception("PPS")
0203         << "useTimingEfficiencyPerRP and useTimingEfficiencyPerPlane should not be simultaneously set true.";
0204 
0205   // v position of strip 0
0206   stripZeroPosition_ = RPTopology::last_strip_to_border_dist_ + (RPTopology::no_of_strips_ - 1) * RPTopology::pitch_ -
0207                        RPTopology::y_width_ / 2.;
0208 }
0209 
0210 //----------------------------------------------------------------------------------------------------
0211 
0212 void PPSDirectProtonSimulation::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0213   edm::ParameterSetDescription desc;
0214   desc.addUntracked<unsigned int>("verbosity", 0);
0215 
0216   desc.add<std::string>("lhcInfoLabel", "")->setComment("label of the LHCInfo record");
0217   desc.add<std::string>("opticsLabel", "")->setComment("label of the optics records");
0218 
0219   desc.add<edm::InputTag>("hepMCTag", edm::InputTag("generator", "unsmeared"));
0220 
0221   desc.add<bool>("produceScoringPlaneHits", true);
0222   desc.add<bool>("produceRecHits", true);
0223 
0224   desc.add<bool>("useEmpiricalApertures", true);
0225 
0226   desc.add<bool>("useTrackingEfficiencyPerRP", false);
0227   desc.add<bool>("useTimingEfficiencyPerRP", false);
0228   desc.add<bool>("useTrackingEfficiencyPerPlane", false);
0229   desc.add<bool>("useTimingEfficiencyPerPlane", false);
0230 
0231   desc.add<bool>("produceHitsRelativeToBeam", true);
0232   desc.add<bool>("roundToPitch", true);
0233   desc.add<bool>("checkIsHit", true);
0234   desc.add<double>("pitchStrips", 66.e-3);              // in mm
0235   desc.add<double>("insensitiveMarginStrips", 34.e-3);  // in mm
0236 
0237   desc.add<double>("pitchPixelsHor", 100.e-3);
0238   desc.add<double>("pitchPixelsVer", 150.e-3);
0239 
0240   descriptions.add("ppsDirectProtonSimulation", desc);
0241 }
0242 
0243 //----------------------------------------------------------------------------------------------------
0244 
0245 void PPSDirectProtonSimulation::produce(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0246   // get input
0247   edm::Handle<edm::HepMCProduct> hepmc_prod;
0248   iEvent.getByToken(hepMCToken_, hepmc_prod);
0249 
0250   // get conditions
0251   auto const &lhcInfo = iSetup.getData(tokenLHCInfo_);
0252   auto const &beamParameters = iSetup.getData(tokenBeamParameters_);
0253   auto const &ppt = iSetup.getData(pixelTopologyToken_);
0254   auto const &opticalFunctions = iSetup.getData(tokenOpticalFunctions_);
0255   auto const &geometry = iSetup.getData(tokenGeometry_);
0256   auto const &directSimuData = iSetup.getData(tokenDirectSimuData_);
0257 
0258   if (directSimuDataRcdWatcher_.check(iSetup)) {
0259     timeResolutionDiamonds45_ =
0260         std::make_unique<TF1>(TF1("timeResolutionDiamonds45", directSimuData.getTimeResolutionDiamonds45().c_str()));
0261     timeResolutionDiamonds56_ =
0262         std::make_unique<TF1>(TF1("timeResolutionDiamonds56", directSimuData.getTimeResolutionDiamonds56().c_str()));
0263 
0264     empiricalAperture45_ =
0265         std::make_unique<TF2>(TF2("empiricalAperture45", directSimuData.getEmpiricalAperture45().c_str()));
0266     empiricalAperture56_ =
0267         std::make_unique<TF2>(TF2("empiricalAperture56", directSimuData.getEmpiricalAperture56().c_str()));
0268 
0269     // load the efficiency maps
0270     if (useTrackingEfficiencyPerRP_ || useTimingEfficiencyPerRP_)
0271       efficiencyMapsPerRP_ = directSimuData.loadEffeciencyHistogramsPerRP();
0272     if (useTrackingEfficiencyPerPlane_ || useTimingEfficiencyPerPlane_)
0273       efficiencyMapsPerPlane_ = directSimuData.loadEffeciencyHistogramsPerPlane();
0274   }
0275 
0276   // prepare outputs
0277   std::unique_ptr<edm::DetSetVector<TotemRPRecHit>> pStripRecHits(new edm::DetSetVector<TotemRPRecHit>());
0278   std::unique_ptr<edm::DetSetVector<CTPPSDiamondRecHit>> pDiamondRecHits(new edm::DetSetVector<CTPPSDiamondRecHit>());
0279   std::unique_ptr<edm::DetSetVector<CTPPSPixelRecHit>> pPixelRecHits(new edm::DetSetVector<CTPPSPixelRecHit>());
0280 
0281   auto pStripRecHitsPerParticle = std::make_unique<std::map<int, edm::DetSetVector<TotemRPRecHit>>>();
0282   auto pDiamondRecHitsPerParticle = std::make_unique<std::map<int, edm::DetSetVector<CTPPSDiamondRecHit>>>();
0283   auto pPixelRecHitsPerParticle = std::make_unique<std::map<int, edm::DetSetVector<CTPPSPixelRecHit>>>();
0284 
0285   std::unique_ptr<std::vector<CTPPSLocalTrackLite>> pTracks(new std::vector<CTPPSLocalTrackLite>());
0286 
0287   // get random engine
0288   edm::Service<edm::RandomNumberGenerator> rng;
0289   CLHEP::HepRandomEngine *engine = &rng->getEngine(iEvent.streamID());
0290 
0291   // loop over event vertices
0292   auto evt = hepmc_prod->GetEvent();
0293   for (auto it_vtx = evt->vertices_begin(); it_vtx != evt->vertices_end(); ++it_vtx) {
0294     auto vtx = *(it_vtx);
0295 
0296     // loop over outgoing particles
0297     for (auto it_part = vtx->particles_out_const_begin(); it_part != vtx->particles_out_const_end(); ++it_part) {
0298       auto part = *(it_part);
0299 
0300       // accept only stable protons
0301       if (part->pdg_id() != 2212)
0302         continue;
0303 
0304       if (part->status() != 1 && part->status() < 83)
0305         continue;
0306 
0307       processProton(vtx,
0308                     part,
0309                     geometry,
0310                     lhcInfo,
0311                     beamParameters,
0312                     ppt,
0313                     opticalFunctions,
0314                     engine,
0315                     *pTracks,
0316                     *pStripRecHits,
0317                     *pPixelRecHits,
0318                     *pDiamondRecHits,
0319                     *pStripRecHitsPerParticle,
0320                     *pPixelRecHitsPerParticle,
0321                     *pDiamondRecHitsPerParticle);
0322     }
0323   }
0324 
0325   if (produceScoringPlaneHits_)
0326     iEvent.put(std::move(pTracks));
0327 
0328   if (produceRecHits_) {
0329     iEvent.put(std::move(pStripRecHits));
0330     iEvent.put(std::move(pPixelRecHits));
0331     iEvent.put(std::move(pDiamondRecHits));
0332 
0333     iEvent.put(std::move(pStripRecHitsPerParticle));
0334     iEvent.put(std::move(pPixelRecHitsPerParticle));
0335     iEvent.put(std::move(pDiamondRecHitsPerParticle));
0336   }
0337 }
0338 
0339 //----------------------------------------------------------------------------------------------------
0340 
0341 void PPSDirectProtonSimulation::processProton(
0342     const HepMC::GenVertex *in_vtx,
0343     const HepMC::GenParticle *in_trk,
0344     const CTPPSGeometry &geometry,
0345     const LHCInfo &lhcInfo,
0346     const CTPPSBeamParameters &beamParameters,
0347     const PPSPixelTopology &ppt,
0348     const LHCInterpolatedOpticalFunctionsSetCollection &opticalFunctions,
0349     CLHEP::HepRandomEngine *rndEngine,
0350     std::vector<CTPPSLocalTrackLite> &out_tracks,
0351     edm::DetSetVector<TotemRPRecHit> &out_strip_hits,
0352     edm::DetSetVector<CTPPSPixelRecHit> &out_pixel_hits,
0353     edm::DetSetVector<CTPPSDiamondRecHit> &out_diamond_hits,
0354     std::map<int, edm::DetSetVector<TotemRPRecHit>> &out_strip_hits_per_particle,
0355     std::map<int, edm::DetSetVector<CTPPSPixelRecHit>> &out_pixel_hits_per_particle,
0356     std::map<int, edm::DetSetVector<CTPPSDiamondRecHit>> &out_diamond_hits_per_particle) const {
0357   /// xi is positive for diffractive protons, thus proton momentum p = (1-xi) * p_nom
0358   /// horizontal component of proton momentum: p_x = th_x * (1-xi) * p_nom
0359 
0360   std::stringstream ssLog;
0361 
0362   // vectors in CMS convention
0363   const HepMC::FourVector &vtx_cms = in_vtx->position();  // in mm
0364   const HepMC::FourVector &mom_cms = in_trk->momentum();
0365 
0366   // transformation to LHC/TOTEM convention
0367   HepMC::FourVector vtx_lhc(-vtx_cms.x(), vtx_cms.y(), -vtx_cms.z(), vtx_cms.t());
0368   HepMC::FourVector mom_lhc(-mom_cms.x(), mom_cms.y(), -mom_cms.z(), mom_cms.t());
0369 
0370   // determine the LHC arm and related parameters
0371   unsigned int arm = 3;
0372   double z_sign;
0373   double beamMomentum = 0.;
0374   double xangle = 0.;
0375   const std::unique_ptr<TF2> *empiricalAperture;
0376   if (mom_lhc.z() < 0)  // sector 45
0377   {
0378     arm = 0;
0379     z_sign = -1;
0380     beamMomentum = beamParameters.getBeamMom45();
0381     xangle = beamParameters.getHalfXangleX45();
0382     empiricalAperture = &empiricalAperture45_;
0383   } else {  // sector 56
0384     arm = 1;
0385     z_sign = +1;
0386     beamMomentum = beamParameters.getBeamMom56();
0387     xangle = beamParameters.getHalfXangleX56();
0388     empiricalAperture = &empiricalAperture56_;
0389   }
0390 
0391   // calculate effective RP arrival time
0392   // effective time mimics the timing calibration -> effective times are distributed about 0
0393   // units:
0394   //    vertex: all components in mm
0395   //    c_light: in mm/ns
0396   //    time_eff: in ns
0397   const double time_eff = (vtx_lhc.t() - z_sign * vtx_lhc.z()) / CLHEP::c_light;
0398 
0399   // calculate kinematics for optics parametrisation
0400   const double p = mom_lhc.rho();
0401   const double xi = 1. - p / beamMomentum;
0402   const double th_x_phys = mom_lhc.x() / p;
0403   const double th_y_phys = mom_lhc.y() / p;
0404   const double vtx_lhc_eff_x = vtx_lhc.x() - vtx_lhc.z() * (mom_lhc.x() / mom_lhc.z() + xangle);
0405   const double vtx_lhc_eff_y = vtx_lhc.y() - vtx_lhc.z() * (mom_lhc.y() / mom_lhc.z());
0406 
0407   if (verbosity_) {
0408     ssLog << "simu: xi = " << xi << ", th_x_phys = " << th_x_phys << ", th_y_phys = " << th_y_phys
0409           << ", vtx_lhc_eff_x = " << vtx_lhc_eff_x << ", vtx_lhc_eff_y = " << vtx_lhc_eff_y << std::endl;
0410   }
0411 
0412   // check empirical aperture
0413   if (useEmpiricalApertures_) {
0414     const auto &xangle = lhcInfo.crossingAngle();
0415     (*empiricalAperture)->SetParameter("xi", xi);
0416     (*empiricalAperture)->SetParameter("xangle", xangle);
0417     const double th_x_th = (*empiricalAperture)->EvalPar(nullptr);
0418 
0419     if (th_x_th > th_x_phys) {
0420       if (verbosity_) {
0421         ssLog << "stop because of empirical appertures";
0422         edm::LogInfo("PPS") << ssLog.str();
0423       }
0424 
0425       return;
0426     }
0427   }
0428 
0429   // transport the proton into each pot/scoring plane
0430   for (const auto &ofp : opticalFunctions) {
0431     CTPPSDetId rpId(ofp.first);
0432     const unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0433 
0434     // first check the arm
0435     if (rpId.arm() != arm)
0436       continue;
0437 
0438     if (verbosity_)
0439       ssLog << "  RP " << rpDecId << std::endl;
0440 
0441     // transport proton
0442     LHCInterpolatedOpticalFunctionsSet::Kinematics k_in = {
0443         vtx_lhc_eff_x * 1E-1, th_x_phys, vtx_lhc_eff_y * 1E-1, th_y_phys, xi};  // conversions: mm -> cm
0444     LHCInterpolatedOpticalFunctionsSet::Kinematics k_out;
0445     ofp.second.transport(k_in, k_out, true);
0446 
0447     double b_x = k_out.x * 1E1, b_y = k_out.y * 1E1;  // conversions: cm -> mm
0448     double a_x = k_out.th_x, a_y = k_out.th_y;
0449 
0450     // if needed, subtract beam position and angle
0451     if (produceHitsRelativeToBeam_) {
0452       // determine beam position
0453       LHCInterpolatedOpticalFunctionsSet::Kinematics k_be_in = {0., 0., 0., 0., 0.};
0454       LHCInterpolatedOpticalFunctionsSet::Kinematics k_be_out;
0455       ofp.second.transport(k_be_in, k_be_out, true);
0456 
0457       a_x -= k_be_out.th_x;
0458       a_y -= k_be_out.th_y;
0459       b_x -= k_be_out.x * 1E1;
0460       b_y -= k_be_out.y * 1E1;  // conversions: cm -> mm
0461     }
0462 
0463     const double z_scoringPlane = ofp.second.getScoringPlaneZ() * 1E1;  // conversion: cm --> mm
0464 
0465     if (verbosity_) {
0466       ssLog << "    proton transported: a_x = " << a_x << " rad, a_y = " << a_y << " rad, b_x = " << b_x
0467             << " mm, b_y = " << b_y << " mm, z = " << z_scoringPlane << " mm" << std::endl;
0468     }
0469 
0470     // RP type
0471     const bool isTrackingRP =
0472         (rpId.subdetId() == CTPPSDetId::sdTrackingStrip || rpId.subdetId() == CTPPSDetId::sdTrackingPixel);
0473     const bool isTimingRP = (rpId.subdetId() == CTPPSDetId::sdTimingDiamond);
0474 
0475     // apply per-RP efficiency
0476     if ((useTimingEfficiencyPerRP_ && isTimingRP) || (useTrackingEfficiencyPerRP_ && isTrackingRP)) {
0477       const auto it = efficiencyMapsPerRP_.find(rpId);
0478 
0479       if (it != efficiencyMapsPerRP_.end()) {
0480         const double r = CLHEP::RandFlat::shoot(rndEngine, 0., 1.);
0481         auto *effMap = it->second.get();
0482         const double eff = effMap->GetBinContent(effMap->FindBin(b_x, b_y));
0483         if (r > eff) {
0484           if (verbosity_)
0485             ssLog << "    stop due to per-RP efficiency" << std::endl;
0486           continue;
0487         }
0488       }
0489     }
0490 
0491     // save scoring plane hit
0492     if (produceScoringPlaneHits_)
0493       out_tracks.emplace_back(
0494           rpId.rawId(), b_x, 0., b_y, 0., 0., 0., 0., 0., 0., CTPPSpixelLocalTrackReconstructionInfo::invalid, 0, 0., 0.);
0495 
0496     // stop if rec hits are not to be produced
0497     if (!produceRecHits_)
0498       continue;
0499 
0500     // loop over all sensors in the RP
0501     for (const auto &detIdInt : geometry.sensorsInRP(rpId)) {
0502       CTPPSDetId detId(detIdInt);
0503 
0504       // determine the track impact point (in global coordinates)
0505       // !! this assumes that local axes (1, 0, 0) and (0, 1, 0) describe the sensor surface
0506       const auto &gl_o = geometry.localToGlobal(detId, CTPPSGeometry::Vector(0, 0, 0));
0507       const auto &gl_a1 = geometry.localToGlobal(detId, CTPPSGeometry::Vector(1, 0, 0)) - gl_o;
0508       const auto &gl_a2 = geometry.localToGlobal(detId, CTPPSGeometry::Vector(0, 1, 0)) - gl_o;
0509 
0510       const double gl_o_z = gl_o.z();
0511 
0512       TMatrixD A(3, 3);
0513       TVectorD B(3);
0514       A(0, 0) = a_x;
0515       A(0, 1) = -gl_a1.x();
0516       A(0, 2) = -gl_a2.x();
0517       B(0) = gl_o.x() - b_x;
0518       A(1, 0) = a_y;
0519       A(1, 1) = -gl_a1.y();
0520       A(1, 2) = -gl_a2.y();
0521       B(1) = gl_o.y() - b_y;
0522       A(2, 0) = z_sign;
0523       A(2, 1) = -gl_a1.z();
0524       A(2, 2) = -gl_a2.z();
0525       B(2) = gl_o_z - z_scoringPlane;
0526       TMatrixD Ai(3, 3);
0527       Ai = A.Invert();
0528       TVectorD P(3);
0529       P = Ai * B;
0530 
0531       double ze = P(0);
0532       const CTPPSGeometry::Vector h_glo(a_x * ze + b_x, a_y * ze + b_y, z_sign * ze + z_scoringPlane);
0533 
0534       // hit in local coordinates
0535       auto h_loc = geometry.globalToLocal(detId, h_glo);
0536 
0537       if (verbosity_) {
0538         ssLog << std::endl
0539               << "    de z = " << P(0) << " mm, p1 = " << P(1) << " mm, p2 = " << P(2) << " mm" << std::endl
0540               << "    h_glo: x = " << h_glo.x() << " mm, y = " << h_glo.y() << " mm, z = " << h_glo.z() << " mm"
0541               << std::endl
0542               << "    h_loc: c1 = " << h_loc.x() << " mm, c2 = " << h_loc.y() << " mm, c3 = " << h_loc.z() << " mm"
0543               << std::endl;
0544       }
0545 
0546       // apply per-plane efficiency
0547       if ((useTimingEfficiencyPerPlane_ && isTimingRP) || (useTrackingEfficiencyPerPlane_ && isTrackingRP)) {
0548         const auto it = efficiencyMapsPerPlane_.find(detId);
0549 
0550         if (it != efficiencyMapsPerPlane_.end()) {
0551           const double r = CLHEP::RandFlat::shoot(rndEngine, 0., 1.);
0552           auto *effMap = it->second.get();
0553           const double eff = effMap->GetBinContent(effMap->FindBin(h_glo.x(), h_glo.y()));
0554           if (r > eff) {
0555             if (verbosity_)
0556               ssLog << "    stop due to per-plane efficiency" << std::endl;
0557             continue;
0558           }
0559         }
0560       }
0561 
0562       // strips
0563       if (detId.subdetId() == CTPPSDetId::sdTrackingStrip) {
0564         double u = h_loc.x();
0565         double v = h_loc.y();
0566 
0567         if (verbosity_ > 5)
0568           ssLog << "            u=" << u << ", v=" << v;
0569 
0570         // is it within detector?
0571         if (checkIsHit_ && !RPTopology::IsHit(u, v, insensitiveMarginStrips_)) {
0572           if (verbosity_ > 5)
0573             ssLog << " | no hit" << std::endl;
0574           continue;
0575         }
0576 
0577         // round the measurement
0578         if (roundToPitch_) {
0579           double m = stripZeroPosition_ - v;
0580           signed int strip = (int)floor(m / pitchStrips_ + 0.5);
0581 
0582           v = stripZeroPosition_ - pitchStrips_ * strip;
0583 
0584           if (verbosity_ > 5)
0585             ssLog << " | strip=" << strip;
0586         }
0587 
0588         double sigma = pitchStrips_ / sqrt(12.);
0589 
0590         if (verbosity_ > 5)
0591           ssLog << " | m=" << v << ", sigma=" << sigma << std::endl;
0592 
0593         edm::DetSet<TotemRPRecHit> &hits = out_strip_hits.find_or_insert(detId);
0594         hits.emplace_back(v, sigma);
0595 
0596         edm::DetSet<TotemRPRecHit> &hits_per_particle =
0597             out_strip_hits_per_particle[in_trk->barcode()].find_or_insert(detId);
0598         hits_per_particle.emplace_back(v, sigma);
0599       }
0600 
0601       // diamonds
0602       if (detId.subdetId() == CTPPSDetId::sdTimingDiamond) {
0603         CTPPSDiamondDetId diamondDetId(detIdInt);
0604 
0605         // check acceptance
0606         const auto *dg = geometry.sensor(detIdInt);
0607         const auto &diamondDimensions = dg->getDiamondDimensions();
0608         const auto x_half_width = diamondDimensions.xHalfWidth;
0609         const auto y_half_width = diamondDimensions.yHalfWidth;
0610         const auto z_half_width = diamondDimensions.zHalfWidth;
0611 
0612         if (h_loc.x() < -x_half_width || h_loc.x() > +x_half_width || h_loc.y() < -y_half_width ||
0613             h_loc.y() > +y_half_width)
0614           continue;
0615 
0616         // timing information
0617         const double time_resolution = (diamondDetId.arm() == 0) ? timeResolutionDiamonds45_->Eval(h_glo.x())
0618                                                                  : timeResolutionDiamonds56_->Eval(h_glo.x());
0619 
0620         const double t0 = time_eff + CLHEP::RandGauss::shoot(rndEngine, 0., time_resolution);
0621         const double tot = 1.23456;
0622         const double ch_t_precis = time_resolution;
0623         const int time_slice = 0;
0624 
0625         // build rec hit
0626         const bool multiHit = false;
0627 
0628         CTPPSDiamondRecHit rc(gl_o.x(),
0629                               2. * x_half_width,
0630                               gl_o.y(),
0631                               2. * y_half_width,
0632                               gl_o_z,
0633                               2. * z_half_width,
0634                               t0,
0635                               tot,
0636                               ch_t_precis,
0637                               time_slice,
0638                               HPTDCErrorFlags(),
0639                               multiHit);
0640 
0641         edm::DetSet<CTPPSDiamondRecHit> &hits = out_diamond_hits.find_or_insert(detId);
0642         hits.push_back(rc);
0643 
0644         edm::DetSet<CTPPSDiamondRecHit> &hits_per_particle =
0645             out_diamond_hits_per_particle[in_trk->barcode()].find_or_insert(detId);
0646         hits_per_particle.push_back(rc);
0647       }
0648 
0649       // pixels
0650       if (detId.subdetId() == CTPPSDetId::sdTrackingPixel) {
0651         if (verbosity_) {
0652           CTPPSPixelDetId pixelDetId(detIdInt);
0653           ssLog << "    pixel plane " << pixelDetId.plane() << ": local hit x = " << h_loc.x()
0654                 << " mm, y = " << h_loc.y() << " mm, z = " << h_loc.z() << " mm" << std::endl;
0655         }
0656 
0657         bool module3By2 = (geometry.sensor(detIdInt)->sensorType() != DDD_CTPPS_PIXELS_SENSOR_TYPE_2x2);
0658         if (checkIsHit_ && !ppt.isPixelHit(h_loc.x(), h_loc.y(), module3By2))
0659           continue;
0660 
0661         if (roundToPitch_) {
0662           h_loc.SetX(pitchPixelsHor_ * floor(h_loc.x() / pitchPixelsHor_ + 0.5));
0663           h_loc.SetY(pitchPixelsVer_ * floor(h_loc.y() / pitchPixelsVer_ + 0.5));
0664         }
0665 
0666         if (verbosity_ > 5)
0667           ssLog << "            hit accepted: m1 = " << h_loc.x() << " mm, m2 = " << h_loc.y() << " mm" << std::endl;
0668 
0669         const double sigmaHor = pitchPixelsHor_ / sqrt(12.);
0670         const double sigmaVer = pitchPixelsVer_ / sqrt(12.);
0671 
0672         const LocalPoint lp(h_loc.x(), h_loc.y(), h_loc.z());
0673         const LocalError le(sigmaHor, 0., sigmaVer);
0674 
0675         edm::DetSet<CTPPSPixelRecHit> &hits = out_pixel_hits.find_or_insert(detId);
0676         hits.emplace_back(lp, le);
0677 
0678         edm::DetSet<CTPPSPixelRecHit> &hits_per_particle =
0679             out_pixel_hits_per_particle[in_trk->barcode()].find_or_insert(detId);
0680         hits_per_particle.emplace_back(lp, le);
0681       }
0682     }
0683   }
0684 
0685   if (verbosity_)
0686     edm::LogInfo("PPS") << ssLog.str();
0687 }
0688 
0689 //----------------------------------------------------------------------------------------------------
0690 
0691 DEFINE_FWK_MODULE(PPSDirectProtonSimulation);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!