AlignmentRelative/plugins/PPSFastLocalSimulation.cc

0001 /****************************************************************************
0002 * Authors:
0003 *  Jan Kašpar (jan.kaspar@gmail.com)
0004 ****************************************************************************/
0005
0006 #include "FWCore/Framework/interface/stream/EDProducer.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 #include "FWCore/ServiceRegistry/interface/Service.h"
0011 #include "FWCore/Utilities/interface/ESGetToken.h"
0012 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/Framework/interface/ESHandle.h"
0015
0016 #include "CLHEP/Random/RandomEngine.h"
0017 #include "CLHEP/Random/RandGauss.h"
0018 #include "CLHEP/Random/RandExponential.h"
0019
0020 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0021
0022 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0023 #include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h"
0024
0025 #include "DataFormats/Common/interface/DetSetVector.h"
0026 #include "DataFormats/CTPPSReco/interface/TotemRPRecHit.h"
0027 #include "DataFormats/CTPPSReco/interface/CTPPSDiamondRecHit.h"
0028 #include "DataFormats/CTPPSReco/interface/CTPPSPixelRecHit.h"
0029
0030 #include "Geometry/VeryForwardRPTopology/interface/RPTopology.h"
0031 #include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h"
0032 #include "Geometry/Records/interface/VeryForwardMisalignedGeometryRecord.h"
0033
0034 #include "CalibPPS/AlignmentRelative/interface/Utilities.h"
0035
0036 #include "TMath.h"
0037 #include "TMatrixD.h"
0038 #include "TVectorD.h"
0039
0040 /**
0041  *\brief Fast (no G4) proton simulation in within one station.
0042  * Uses misaligned geometry.
0043  */
0044 class PPSFastLocalSimulation : public edm::stream::EDProducer<> {
0045 public:
0046   PPSFastLocalSimulation(const edm::ParameterSet &);
0047   ~PPSFastLocalSimulation() override;
0048
0049 protected:
0050   /// verbosity level
0051   unsigned int verbosity_;
0052
0053   /// whether a HepMC description of the proton shall be saved in the event
0054   bool makeHepMC_;
0055
0056   /// whether the hits of the proton shall be calculated and saved
0057   bool makeHits_;
0058
0059   /// the list of RPs to simulate
0060   std::vector<unsigned int> RPs_;
0061
0062   /// number of particles to generate per event
0063   unsigned int particlesPerEvent_;
0064
0065   /// particle energy and momentum
0066   double particle_E_, particle_p_;
0067
0068   /// the "origin" of tracks, in mm
0069   double z0_;
0070
0071   /// whether measurement values shall be rounded to the nearest strip
0072   bool roundToPitch_;
0073
0074   /// in mm
0075   double pitchStrips_, pitchDiamonds_, pitchPixels_;
0076
0077   /// size of insensitive margin at sensor's edge facing the beam, in mm
0078   double insensitiveMarginStrips_;
0079
0080   struct Distribution {
0081     enum Type { dtBox, dtGauss, dtGaussLimit } type_;
0082     double x_mean_, x_width_, x_min_, x_max_;
0083     double y_mean_, y_width_, y_min_, y_max_;
0084
0085     Distribution(const edm::ParameterSet &);
0086
0087     void Generate(CLHEP::HepRandomEngine &rndEng, double &x, double &y);
0088   };
0089
0090   /// position parameters in mm
0091   Distribution position_dist_;
0092
0093   /// angular parameters in rad
0094   Distribution angular_dist_;
0095
0096   //---------- internal parameters ----------
0097
0098   /// v position of strip 0, in mm
0099   double stripZeroPosition_;
0100
0101   edm::ESGetToken<CTPPSGeometry, VeryForwardMisalignedGeometryRecord> esTokenGeometry_;
0102
0103   void GenerateTrack(unsigned int pi,
0104                      CLHEP::HepRandomEngine &rndEng,
0105                      HepMC::GenEvent *gEv,
0106                      std::unique_ptr<edm::DetSetVector<TotemRPRecHit>> &stripHitColl,
0107                      std::unique_ptr<edm::DetSetVector<CTPPSDiamondRecHit>> &diamondHitColl,
0108                      std::unique_ptr<edm::DetSetVector<CTPPSPixelRecHit>> &pixelHitColl,
0109                      const CTPPSGeometry &geometry);
0110
0111   //---------- framework methods ----------
0112
0113   void produce(edm::Event &, const edm::EventSetup &) override;
0114 };
0115
0116 //----------------------------------------------------------------------------------------------------
0117
0118 using namespace edm;
0119 using namespace std;
0120 using namespace CLHEP;
0121 using namespace HepMC;
0122
0123 //----------------------------------------------------------------------------------------------------
0124
0125 PPSFastLocalSimulation::Distribution::Distribution(const edm::ParameterSet &ps) {
0126   // get type
0127   string typeName = ps.getParameter<string>("type");
0128   if (!typeName.compare("box"))
0129     type_ = dtBox;
0130   else if (!typeName.compare("gauss"))
0131     type_ = dtGauss;
0132   else if (!typeName.compare("gauss-limit"))
0133     type_ = dtGaussLimit;
0134   else
0135     throw cms::Exception("PPS") << "Unknown distribution type `" << typeName << "'.";
0136
0137   x_mean_ = ps.getParameter<double>("x_mean");
0138   x_width_ = ps.getParameter<double>("x_width");
0139   x_min_ = ps.getParameter<double>("x_min");
0140   x_max_ = ps.getParameter<double>("x_max");
0141
0142   y_mean_ = ps.getParameter<double>("y_mean");
0143   y_width_ = ps.getParameter<double>("y_width");
0144   y_min_ = ps.getParameter<double>("y_min");
0145   y_max_ = ps.getParameter<double>("y_max");
0146 }
0147
0148 //----------------------------------------------------------------------------------------------------
0149
0150 void PPSFastLocalSimulation::Distribution::Generate(CLHEP::HepRandomEngine &rndEng, double &x, double &y) {
0151   switch (type_) {
0152     case dtBox:
0153       x = x_mean_ + x_width_ * (rndEng.flat() - 0.5);
0154       y = y_mean_ + y_width_ * (rndEng.flat() - 0.5);
0155       break;
0156
0157     case dtGauss:
0158       x = x_mean_ + RandGauss::shoot(&rndEng) * x_width_;
0159       y = y_mean_ + RandGauss::shoot(&rndEng) * y_width_;
0160       break;
0161
0162     case dtGaussLimit: {
0163       const double u_x = rndEng.flat(), u_y = rndEng.flat();
0164
0165       const double cdf_x_min = (1. + TMath::Erf((x_min_ - x_mean_) / x_width_ / sqrt(2.))) / 2.;
0166       const double cdf_x_max = (1. + TMath::Erf((x_max_ - x_mean_) / x_width_ / sqrt(2.))) / 2.;
0167       const double a_x = cdf_x_max - cdf_x_min, b_x = cdf_x_min;
0168
0169       const double cdf_y_min = (1. + TMath::Erf((y_min_ - y_mean_) / y_width_ / sqrt(2.))) / 2.;
0170       const double cdf_y_max = (1. + TMath::Erf((y_max_ - y_mean_) / y_width_ / sqrt(2.))) / 2.;
0171       const double a_y = cdf_y_max - cdf_y_min, b_y = cdf_y_min;
0172
0173       x = x_mean_ + x_width_ * sqrt(2.) * TMath::ErfInverse(2. * (a_x * u_x + b_x) - 1.);
0174       y = y_mean_ + y_width_ * sqrt(2.) * TMath::ErfInverse(2. * (a_y * u_y + b_y) - 1.);
0175     }
0176
0177     break;
0178
0179     default:
0180       x = y = 0.;
0181   }
0182 }
0183
0184 //----------------------------------------------------------------------------------------------------
0185
0186 PPSFastLocalSimulation::PPSFastLocalSimulation(const edm::ParameterSet &ps)
0187     : verbosity_(ps.getUntrackedParameter<unsigned int>("verbosity", 0)),
0188
0189       makeHepMC_(ps.getParameter<bool>("makeHepMC")),
0190       makeHits_(ps.getParameter<bool>("makeHits")),
0191
0192       RPs_(ps.getParameter<vector<unsigned int>>("RPs")),
0193
0194       particlesPerEvent_(ps.getParameter<unsigned int>("particlesPerEvent")),
0195       particle_E_(ps.getParameter<double>("particle_E")),
0196       particle_p_(ps.getParameter<double>("particle_p")),
0197       z0_(ps.getParameter<double>("z0")),
0198
0199       roundToPitch_(ps.getParameter<bool>("roundToPitch")),
0200       pitchStrips_(ps.getParameter<double>("pitchStrips")),
0201       pitchDiamonds_(ps.getParameter<double>("pitchDiamonds")),
0202       pitchPixels_(ps.getParameter<double>("pitchPixels")),
0203
0204       insensitiveMarginStrips_(ps.getParameter<double>("insensitiveMarginStrips")),
0205
0206       position_dist_(ps.getParameterSet("position_distribution")),
0207       angular_dist_(ps.getParameterSet("angular_distribution")),
0208
0209       esTokenGeometry_(esConsumes()) {
0210   // v position of strip 0
0211   stripZeroPosition_ = RPTopology::last_strip_to_border_dist_ + (RPTopology::no_of_strips_ - 1) * RPTopology::pitch_ -
0212                        RPTopology::y_width_ / 2.;
0213
0214   // register the output
0215   if (makeHepMC_)
0216     produces<HepMCProduct>();
0217
0218   if (makeHits_) {
0219     produces<DetSetVector<TotemRPRecHit>>();
0220     produces<DetSetVector<CTPPSDiamondRecHit>>();
0221     produces<DetSetVector<CTPPSPixelRecHit>>();
0222   }
0223 }
0224
0225 //----------------------------------------------------------------------------------------------------
0226
0227 PPSFastLocalSimulation::~PPSFastLocalSimulation() {}
0228
0229 //----------------------------------------------------------------------------------------------------
0230
0231 void PPSFastLocalSimulation::GenerateTrack(unsigned int idx,
0232                                            CLHEP::HepRandomEngine &rndEng,
0233                                            HepMC::GenEvent *gEv,
0234                                            unique_ptr<edm::DetSetVector<TotemRPRecHit>> &stripHitColl,
0235                                            unique_ptr<edm::DetSetVector<CTPPSDiamondRecHit>> &diamondHitColl,
0236                                            unique_ptr<edm::DetSetVector<CTPPSPixelRecHit>> &pixelHitColl,
0237                                            const CTPPSGeometry &geometry) {
0238   // generate track
0239   double bx = 0., by = 0., ax = 0., ay = 0.;
0240   position_dist_.Generate(rndEng, bx, by);
0241   angular_dist_.Generate(rndEng, ax, ay);
0242
0243   if (verbosity_ > 5)
0244     printf("\tax = %.3f mrad, bx = %.3f mm, ay = %.3f mrad, by = %.3f mm, z0 = %.3f m\n",
0245            ax * 1E3,
0246            bx,
0247            ay * 1E3,
0248            by,
0249            z0_ * 1E-3);
0250
0251   // add HepMC track description
0252   if (makeHepMC_) {
0253     GenVertex *gVx = new GenVertex(HepMC::FourVector(bx, by, z0_, 0.));
0254     gEv->add_vertex(gVx);
0255
0256     GenParticle *gPe;
0257     double az = sqrt(1. - ax * ax - ay * ay);
0258     gPe = new GenParticle(HepMC::FourVector(particle_p_ * ax, particle_p_ * ay, particle_p_ * az, particle_E_),
0259                           2212,
0260                           1);  // add a proton in final state
0261     gPe->suggest_barcode(idx + 1);
0262     gVx->add_particle_out(gPe);
0263   }
0264
0265   if (makeHits_) {
0266     // check all sensors known to geometry
0267     for (CTPPSGeometry::mapType::const_iterator it = geometry.beginSensor(); it != geometry.endSensor(); ++it) {
0268       // get RP decimal id
0269       CTPPSDetId detId(it->first);
0270       unsigned int decRPId = detId.arm() * 100 + detId.station() * 10 + detId.rp();
0271
0272       // stop if the RP is not selected
0273       if (find(RPs_.begin(), RPs_.end(), decRPId) == RPs_.end())
0274         continue;
0275
0276       // keep only 1 diamond channel to represent 1 plane
0277       if (detId.subdetId() == CTPPSDetId::sdTimingDiamond) {
0278         CTPPSDiamondDetId channelId(it->first);
0279         if (channelId.channel() != 0)
0280           continue;
0281       }
0282
0283       if (verbosity_ > 5) {
0284         printf("        ");
0285         printId(it->first);
0286         printf(": ");
0287       }
0288
0289       // determine the track impact point (in global coordinates)
0290       // !! this assumes that local axes (1, 0, 0) and (0, 1, 0) describe the sensor surface
0291       const auto gl_o = geometry.localToGlobal(detId, CTPPSGeometry::Vector(0, 0, 0));
0292       const auto gl_a1 = geometry.localToGlobal(detId, CTPPSGeometry::Vector(1, 0, 0)) - gl_o;
0293       const auto gl_a2 = geometry.localToGlobal(detId, CTPPSGeometry::Vector(0, 1, 0)) - gl_o;
0294
0295       TMatrixD A(3, 3);
0296       TVectorD B(3);
0297       A(0, 0) = ax;
0298       A(0, 1) = -gl_a1.x();
0299       A(0, 2) = -gl_a2.x();
0300       B(0) = gl_o.x() - bx;
0301       A(1, 0) = ay;
0302       A(1, 1) = -gl_a1.y();
0303       A(1, 2) = -gl_a2.y();
0304       B(1) = gl_o.y() - by;
0305       A(2, 0) = 1.;
0306       A(2, 1) = -gl_a1.z();
0307       A(2, 2) = -gl_a2.z();
0308       B(2) = gl_o.z() - z0_;
0309       TMatrixD Ai(3, 3);
0310       Ai = A.Invert();
0311       TVectorD P(3);
0312       P = Ai * B;
0313
0314       double de_z = P(0);
0315       CTPPSGeometry::Vector h_glo(ax * de_z + bx, ay * de_z + by, de_z + z0_);
0316
0317       // hit in local coordinates
0318       CTPPSGeometry::Vector h_loc = geometry.globalToLocal(detId, h_glo);
0319
0320       // strips
0321       if (detId.subdetId() == CTPPSDetId::sdTrackingStrip) {
0322         double u = h_loc.x();
0323         double v = h_loc.y();
0324
0325         if (verbosity_ > 5)
0326           printf("            u=%+8.4f, v=%+8.4f", u, v);
0327
0328         // is it within detector?
0329         if (!RPTopology::IsHit(u, v, insensitiveMarginStrips_)) {
0330           if (verbosity_ > 5)
0331             printf(" | no hit\n");
0332           continue;
0333         }
0334
0335         // round the measurement
0336         if (roundToPitch_) {
0337           double m = stripZeroPosition_ - v;
0338           signed int strip = (int)floor(m / pitchStrips_ + 0.5);
0339
0340           v = stripZeroPosition_ - pitchStrips_ * strip;
0341
0342           if (verbosity_ > 5)
0343             printf(" | strip=%+4i", strip);
0344         }
0345
0346         double sigma = pitchStrips_ / sqrt(12.);
0347
0348         if (verbosity_ > 5)
0349           printf(" | m=%+8.4f, sigma=%+8.4f\n", v, sigma);
0350
0351         DetSet<TotemRPRecHit> &hits = stripHitColl->find_or_insert(detId);
0352         hits.emplace_back(v, sigma);
0353       }
0354
0355       // diamonds
0356       if (detId.subdetId() == CTPPSDetId::sdTimingDiamond) {
0357         if (roundToPitch_) {
0358           h_loc.SetX(pitchDiamonds_ * floor(h_loc.x() / pitchDiamonds_ + 0.5));
0359         }
0360
0361         if (verbosity_ > 5)
0362           printf("            m = %.3f\n", h_loc.x());
0363
0364         const double width = pitchDiamonds_;
0365
0366         DetSet<CTPPSDiamondRecHit> &hits = diamondHitColl->find_or_insert(detId);
0367         hits.emplace_back(h_loc.x(), width, 0., 0., 0., 0., 0., 0., 0., 0, HPTDCErrorFlags(), false);
0368       }
0369
0370       // pixels
0371       if (detId.subdetId() == CTPPSDetId::sdTrackingPixel) {
0372         if (roundToPitch_) {
0373           h_loc.SetX(pitchPixels_ * floor(h_loc.x() / pitchPixels_ + 0.5));
0374           h_loc.SetY(pitchPixels_ * floor(h_loc.y() / pitchPixels_ + 0.5));
0375         }
0376
0377         if (verbosity_ > 5)
0378           printf("            m1 = %.3f, m2 = %.3f\n", h_loc.x(), h_loc.y());
0379
0380         const double sigma = pitchPixels_ / sqrt(12.);
0381
0382         const LocalPoint lp(h_loc.x(), h_loc.y(), h_loc.z());
0383         const LocalError le(sigma, 0., sigma);
0384
0385         DetSet<CTPPSPixelRecHit> &hits = pixelHitColl->find_or_insert(detId);
0386         hits.emplace_back(lp, le);
0387       }
0388     }
0389   }
0390 }
0391
0392 //----------------------------------------------------------------------------------------------------
0393
0394 void PPSFastLocalSimulation::produce(edm::Event &event, const edm::EventSetup &es) {
0395   if (verbosity_ > 2)
0396     printf(">> PPSFastLocalSimulation::produce > event %llu\n", event.id().event());
0397
0398   Service<edm::RandomNumberGenerator> rng;
0399   HepRandomEngine &rndEng = rng->getEngine(event.streamID());
0400
0401   if (verbosity_ > 5)
0402     printf("\tseed = %li\n", rndEng.getSeed());
0403
0404   // get geometry
0405   auto const &geometry = es.getData(esTokenGeometry_);
0406
0407   // initialize products
0408   GenEvent *gEv = new GenEvent();
0409   gEv->set_event_number(event.id().event());
0410
0411   unique_ptr<DetSetVector<TotemRPRecHit>> stripHitColl(new DetSetVector<TotemRPRecHit>());
0412   unique_ptr<DetSetVector<CTPPSDiamondRecHit>> diamondHitColl(new DetSetVector<CTPPSDiamondRecHit>());
0413   unique_ptr<DetSetVector<CTPPSPixelRecHit>> pixelHitColl(new DetSetVector<CTPPSPixelRecHit>());
0414
0415   // run particle loop
0416   for (unsigned int pi = 0; pi < particlesPerEvent_; pi++) {
0417     if (verbosity_ > 5)
0418       printf("    generating track %u\n", pi);
0419
0420     GenerateTrack(pi, rndEng, gEv, stripHitColl, diamondHitColl, pixelHitColl, geometry);
0421   }
0422
0423   // save products
0424   if (makeHepMC_) {
0425     unique_ptr<HepMCProduct> hepMCoutput(new HepMCProduct());
0426     hepMCoutput->addHepMCData(gEv);
0427     event.put(std::move(hepMCoutput));
0428   }
0429
0430   if (makeHits_) {
0431     event.put(std::move(stripHitColl));
0432     event.put(std::move(diamondHitColl));
0433     event.put(std::move(pixelHitColl));
0434   }
0435 }
0436
0437 //----------------------------------------------------------------------------------------------------
0438
0439 DEFINE_FWK_MODULE(PPSFastLocalSimulation);