GEMDigitizer/plugins/ME0ReDigiProducer.cc

0001 #ifndef SimMuon_GEMDigitizer_ME0ReDigiProducer_h
0002 #define SimMuon_GEMDigitizer_ME0ReDigiProducer_h
0003
0004 /*
0005  * This module smears and discretizes the timing and position of the
0006  * ME0 pseudo digis.
0007  */
0008
0009 #include "FWCore/Framework/interface/stream/EDProducer.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 #include "FWCore/Utilities/interface/ESGetToken.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/Framework/interface/ConsumesCollector.h"
0016 #include "FWCore/ServiceRegistry/interface/Service.h"
0017 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0018 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/Utilities/interface/Exception.h"
0021
0022 #include "DataFormats/GEMDigi/interface/ME0DigiPreRecoCollection.h"
0023 #include "DataFormats/Common/interface/Handle.h"
0024
0025 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0026 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0027 #include "Geometry/GEMGeometry/interface/ME0Geometry.h"
0028 #include "Geometry/GEMGeometry/interface/ME0EtaPartitionSpecs.h"
0029
0030 #include "CLHEP/Random/RandGaussQ.h"
0031 #include "CLHEP/Random/RandFlat.h"
0032 #include "CLHEP/Units/PhysicalConstants.h"
0033
0034 #include <sstream>
0035 #include <string>
0036 #include <map>
0037 #include <vector>
0038
0039 typedef class Point3DBase<float, LocalTag> LocalPoint;
0040
0041 class ME0ReDigiProducer : public edm::stream::EDProducer<> {
0042 private:
0043   //Class used to define custom geometry, if required
0044   //assume that all ME0 chambers have the same dimension
0045   //and for the same layer have the same radial and Z position
0046   //Good for now, can build in support for more varied geos later
0047   //if necessary
0048   class TemporaryGeometry {
0049   public:
0050     TemporaryGeometry(const ME0Geometry* geometry,
0051                       const unsigned int numberOfStrips,
0052                       const unsigned int numberOfPartitions);
0053     ~TemporaryGeometry();
0054     unsigned int findEtaPartition(float locY) const;
0055     const TrapezoidalStripTopology* getTopo(const unsigned int partIdx) const { return stripTopos[partIdx]; }
0056     float getPartCenter(const unsigned int partIdx) const;  //position of part. in chamber
0057     float getCentralTOF(const ME0DetId& me0Id, unsigned int partIdx) const {
0058       return tofs[me0Id.layer() - 1][partIdx];
0059     }  //in detId layer numbers stat at 1
0060     unsigned int numLayers() const { return tofs.size(); }
0061
0062   private:
0063     TrapezoidalStripTopology* buildTopo(const std::vector<float>& _p) const;
0064
0065   private:
0066     float middleDistanceFromBeam;                       // radiusOfMainPartitionInCenter;
0067     std::vector<TrapezoidalStripTopology*> stripTopos;  // vector of Topos, one for each part
0068     std::vector<std::vector<double>> tofs;              //TOF to center of the partition:  [layer][part]
0069     std::vector<float> partitionTops;                   //Top of each partition in the chamber's local coords
0070   };
0071
0072 public:
0073   explicit ME0ReDigiProducer(const edm::ParameterSet& ps);
0074
0075   ~ME0ReDigiProducer() override;
0076
0077   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0078
0079   void produce(edm::Event&, const edm::EventSetup&) override;
0080
0081   void buildDigis(const ME0DigiPreRecoCollection&,
0082                   ME0DigiPreRecoCollection&,
0083                   ME0DigiPreRecoMap&,
0084                   CLHEP::HepRandomEngine* engine);
0085
0086 private:
0087   void fillCentralTOFs();
0088   void getStripProperties(const ME0EtaPartition* etaPart,
0089                           const ME0DigiPreReco* inDigi,
0090                           float& tof,
0091                           int& strip,
0092                           LocalPoint& digiLocalPoint,
0093                           LocalError& digiLocalError) const;
0094   int getCustomStripProperties(const ME0DetId& detId,
0095                                const ME0DigiPreReco* inDigi,
0096                                float& tof,
0097                                int& strip,
0098                                LocalPoint& digiLocalPoint,
0099                                LocalError& digiLocalError) const;
0100
0101   typedef std::tuple<unsigned int, unsigned int, unsigned int> DigiIndicies;
0102   typedef std::map<DigiIndicies, unsigned int> ChamberDigiMap;
0103   //fills map...returns -1 if digi is not already in the map
0104   unsigned int fillDigiMap(
0105       ChamberDigiMap& chDigiMap, unsigned int bx, unsigned int part, unsigned int strip, unsigned int currentIDX) const;
0106
0107   //paramters
0108   const float bxWidth;              // ns
0109   bool useCusGeoFor1PartGeo;        //Use custom strips and partitions for digitization for single partition geometry
0110   bool usePads;                     //sets strip granularity to x2 coarser
0111   unsigned int numberOfStrips;      // Custom number of strips per partition
0112   unsigned int numberOfPartitions;  // Custom number of partitions per chamber
0113   double neutronAcceptance;         // fraction of neutron events to keep in event (>= 1 means no filtering)
0114   double timeResolution;            // smear time by gaussian with this sigma (in ns)....negative for no smearing
0115   int minBXReadout;                 // Minimum BX to readout
0116   int maxBXReadout;                 // Maximum BX to readout
0117   std::vector<int> layerReadout;  // Don't readout layer if entry is 0 (Layer number 1 in the numbering scheme is idx 0)
0118   bool mergeDigis;                // Keep only one digi at the same chamber, strip, partition, and BX
0119   edm::EDGetTokenT<ME0DigiPreRecoCollection> token;
0120   edm::ESGetToken<ME0Geometry, MuonGeometryRecord> geom_token_;
0121
0122   bool useBuiltinGeo;
0123   const ME0Geometry* geometry;
0124   TemporaryGeometry* tempGeo;
0125   std::vector<std::vector<double>> tofs;  //used for built in geo
0126 };
0127
0128 ME0ReDigiProducer::TemporaryGeometry::TemporaryGeometry(const ME0Geometry* geometry,
0129                                                         const unsigned int numberOfStrips,
0130                                                         const unsigned int numberOfPartitions) {
0131   //First test geometry to make sure that it is compatible with our assumptions
0132   const auto& chambers = geometry->chambers();
0133   if (chambers.empty())
0134     throw cms::Exception("Setup")
0135         << "ME0ReDigiProducer::TemporaryGeometry::TemporaryGeometry() - No ME0Chambers in geometry.";
0136   const auto* mainChamber = chambers.front();
0137   const unsigned int nLayers = chambers.front()->nLayers();
0138   if (!nLayers)
0139     throw cms::Exception("Setup")
0140         << "ME0ReDigiProducer::TemporaryGeometry::TemporaryGeometry() - ME0Chamber has no layers.";
0141   const auto* mainLayer = mainChamber->layers()[0];
0142   if (!mainLayer->nEtaPartitions())
0143     throw cms::Exception("Setup")
0144         << "ME0ReDigiProducer::TemporaryGeometry::TemporaryGeometry() - ME0Layer has no partitions.";
0145   if (mainLayer->nEtaPartitions() != 1)
0146     throw cms::Exception("Setup") << "ME0ReDigiProducer::TemporaryGeometry::TemporaryGeometry() - This module is only "
0147                                      "compatitble with geometries that contain only one partition per ME0Layer.";
0148
0149   const auto* mainPartition = mainLayer->etaPartitions()[0];
0150   const TrapezoidalStripTopology* mainTopo = dynamic_cast<const TrapezoidalStripTopology*>(&mainPartition->topology());
0151   if (!mainTopo)
0152     throw cms::Exception("Setup") << "ME0ReDigiProducer::TemporaryGeometry::TemporaryGeometry() - ME0 strip topology "
0153                                      "must be of type TrapezoidalStripTopology. This module cannot be used";
0154
0155   for (const auto& chamber : geometry->chambers()) {
0156     if (chamber->nLayers() != int(nLayers))
0157       throw cms::Exception("Setup") << "ME0ReDigiProducer::TemporaryGeometry::TemporaryGeometry() - Not all "
0158                                        "ME0Chambers have the same number of layers. This module cannot be used.";
0159     for (unsigned int iL = 0; iL < nLayers; ++iL) {
0160       if (chamber->layers()[iL]->nEtaPartitions() != mainLayer->nEtaPartitions())
0161         throw cms::Exception("Setup") << "ME0ReDigiProducer::TemporaryGeometry::TemporaryGeometry() - Not all "
0162                                          "ME0Layers have the same number of partitions. This module cannot be used.";
0163       if (chamber->layers()[iL]->etaPartitions()[0]->specs()->parameters() != mainPartition->specs()->parameters())
0164         throw cms::Exception("Setup") << "ME0ReDigiProducer::TemporaryGeometry::TemporaryGeometry() - Not all ME0 ETA "
0165                                          "partitions have the same properties. This module cannot be used.";
0166       if (std::fabs(chamber->layers()[iL]->etaPartitions()[0]->position().z()) !=
0167           std::fabs(mainChamber->layers()[iL]->etaPartitions()[0]->position().z()))
0168         throw cms::Exception("Setup")
0169             << "ME0ReDigiProducer::TemporaryGeometry::TemporaryGeometry() - Not all ME0 ETA partitions in a single "
0170                "layer have the same Z position. This module cannot be used.";
0171     }
0172   }
0173
0174   //Calculate radius to center of partition
0175   middleDistanceFromBeam = mainTopo->radius();
0176
0177   //calculate the top of each eta partition, assuming equal distance in eta between partitions
0178   const auto localTop = LocalPoint(0, mainTopo->stripLength() / 2);
0179   const auto localBottom = LocalPoint(0, -1 * mainTopo->stripLength() / 2);
0180   const auto globalTop = mainPartition->toGlobal(localTop);
0181   const auto globalBottom = mainPartition->toGlobal(localBottom);
0182   const double etaTop = globalTop.eta();
0183   const double etaBottom = globalBottom.eta();
0184   const double zBottom = globalBottom.z();
0185
0186   //Build topologies
0187   partitionTops.reserve(numberOfPartitions);
0188   stripTopos.reserve(numberOfPartitions);
0189   const auto& mainPars = mainPartition->specs()->parameters();
0190   for (unsigned int iP = 0; iP < numberOfPartitions; ++iP) {
0191     const double eta = (etaTop - etaBottom) * double(iP + 1) / double(numberOfPartitions) + etaBottom;
0192     const double distFromBeam = std::fabs(zBottom / std::sinh(eta));
0193     partitionTops.push_back(distFromBeam - middleDistanceFromBeam);
0194     LogDebug("ME0ReDigiProducer::TemporaryGeometry") << "Top of new partition: " << partitionTops.back() << std::endl;
0195
0196     std::vector<float> params(4, 0);
0197
0198     //half width of trapezoid at local coordinate Y
0199     auto getWidth = [&](float locY) -> float {
0200       return (mainPars[2] * (mainPars[1] + mainPars[0]) + locY * (mainPars[1] - mainPars[0])) / (2 * mainPars[2]);
0201     };
0202
0203     params[0] = iP == 0 ? mainPars[0] : getWidth(partitionTops[iP - 1]);  // Half width of bottom of chamber
0204     params[1] =
0205         iP + 1 == numberOfPartitions ? mainPars[1] : getWidth(partitionTops[iP]);  // Half width of top of chamber
0206     params[2] = ((iP + 1 == numberOfPartitions ? localTop.y() : partitionTops[iP]) -
0207                  (iP == 0 ? localBottom.y() : partitionTops[iP - 1])) /
0208                 2;  // Half width of height of chamber
0209     params[3] = numberOfStrips;
0210
0211     stripTopos.push_back(buildTopo(params));
0212   }
0213
0214   //Get TOF at center of each partition
0215   tofs.resize(nLayers);
0216   LogDebug("ME0ReDigiProducer::TemporaryGeometry") << "TOF numbers [layer][partition]: ";
0217   for (unsigned int iL = 0; iL < nLayers; ++iL) {
0218     tofs[iL].resize(numberOfPartitions);
0219     for (unsigned int iP = 0; iP < numberOfPartitions; ++iP) {
0220       const LocalPoint partCenter(0., getPartCenter(iP), 0.);
0221       const GlobalPoint centralGP(mainChamber->layers()[iL]->etaPartitions()[0]->toGlobal(partCenter));
0222       tofs[iL][iP] = (centralGP.mag() / (CLHEP::c_light / CLHEP::cm));  //speed of light [cm/ns]
0223       LogDebug("ME0ReDigiProducer::TemporaryGeometry")
0224           << "[" << iL << "][" << iP << "]=" << tofs[iL][iP] << " " << std::endl;
0225     }
0226   }
0227 }
0228
0229 unsigned int ME0ReDigiProducer::TemporaryGeometry::findEtaPartition(float locY) const {
0230   unsigned int etaPart = stripTopos.size() - 1;
0231   for (unsigned int iP = 0; iP < stripTopos.size(); ++iP) {
0232     if (locY < partitionTops[iP]) {
0233       etaPart = iP;
0234       break;
0235     }
0236   }
0237   return etaPart;
0238 }
0239
0240 float ME0ReDigiProducer::TemporaryGeometry::getPartCenter(const unsigned int partIdx) const {
0241   return stripTopos[partIdx]->radius() - middleDistanceFromBeam;
0242 }
0243
0244 ME0ReDigiProducer::TemporaryGeometry::~TemporaryGeometry() {
0245   for (auto* p : stripTopos) {
0246     delete p;
0247   }
0248 }
0249
0250 TrapezoidalStripTopology* ME0ReDigiProducer::TemporaryGeometry::buildTopo(const std::vector<float>& _p) const {
0251   float b = _p[0];
0252   float B = _p[1];
0253   float h = _p[2];
0254   float r0 = h * (B + b) / (B - b);
0255   float striplength = h * 2;
0256   float strips = _p[3];
0257   float pitch = (b + B) / strips;
0258   int nstrip = static_cast<int>(strips);
0259
0260   LogDebug("ME0ReDigiProducer::TemporaryGeometry")
0261       << "New partition parameters: "
0262       << "bottom width(" << 2 * b << ") top width(" << 2 * B << ") height(" << 2 * h << ") radius to center(" << r0
0263       << ") nStrips(" << strips << ") pitch(" << pitch << ")" << std::endl;
0264
0265   return new TrapezoidalStripTopology(nstrip, pitch, striplength, r0);
0266 }
0267
0268 ME0ReDigiProducer::ME0ReDigiProducer(const edm::ParameterSet& ps)
0269     : bxWidth(25.0),
0270       useCusGeoFor1PartGeo(ps.getParameter<bool>("useCusGeoFor1PartGeo")),
0271       usePads(ps.getParameter<bool>("usePads")),
0272       numberOfStrips(ps.getParameter<unsigned int>("numberOfStrips")),
0273       numberOfPartitions(ps.getParameter<unsigned int>("numberOfPartitions")),
0274       neutronAcceptance(ps.getParameter<double>("neutronAcceptance")),
0275       timeResolution(ps.getParameter<double>("timeResolution")),
0276       minBXReadout(ps.getParameter<int>("minBXReadout")),
0277       maxBXReadout(ps.getParameter<int>("maxBXReadout")),
0278       layerReadout(ps.getParameter<std::vector<int>>("layerReadout")),
0279       mergeDigis(ps.getParameter<bool>("mergeDigis")),
0280       token(consumes<ME0DigiPreRecoCollection>(edm::InputTag(ps.getParameter<std::string>("inputCollection")))),
0281       geom_token_(esConsumes<ME0Geometry, MuonGeometryRecord, edm::Transition::BeginRun>()) {
0282   produces<ME0DigiPreRecoCollection>();
0283   produces<ME0DigiPreRecoMap>();
0284
0285   edm::Service<edm::RandomNumberGenerator> rng;
0286   if (!rng.isAvailable()) {
0287     throw cms::Exception("Configuration")
0288         << "ME0ReDigiProducer::ME0PreRecoDigiProducer() - RandomNumberGeneratorService is not present in configuration "
0289            "file.\n"
0290         << "Add the service in the configuration file or remove the modules that require it.";
0291   }
0292   geometry = nullptr;
0293   tempGeo = nullptr;
0294   useBuiltinGeo = true;
0295
0296   if (useCusGeoFor1PartGeo) {
0297     if (usePads)
0298       numberOfStrips = numberOfStrips / 2;
0299     if (numberOfStrips == 0)
0300       throw cms::Exception("Setup")
0301           << "ME0ReDigiProducer::ME0PreRecoDigiProducer() - Must have at least one strip if using custom geometry.";
0302     if (numberOfPartitions == 0)
0303       throw cms::Exception("Setup")
0304           << "ME0ReDigiProducer::ME0PreRecoDigiProducer() - Must have at least one partition if using custom geometry.";
0305   }
0306
0307   if (neutronAcceptance < 0)
0308     throw cms::Exception("Setup") << "ME0ReDigiProducer::ME0PreRecoDigiProducer() - neutronAcceptance must be >= 0.";
0309 }
0310
0311 ME0ReDigiProducer::~ME0ReDigiProducer() {
0312   if (tempGeo)
0313     delete tempGeo;
0314 }
0315
0316 void ME0ReDigiProducer::beginRun(const edm::Run&, const edm::EventSetup& eventSetup) {
0317   // set geometry
0318   edm::ESHandle<ME0Geometry> hGeom = eventSetup.getHandle(geom_token_);
0319   geometry = &*hGeom;
0320
0321   const auto& chambers = geometry->chambers();
0322   if (chambers.empty())
0323     throw cms::Exception("Setup") << "ME0ReDigiProducer::beginRun() - No ME0Chambers in geometry.";
0324
0325   const unsigned int nLayers = chambers.front()->nLayers();
0326   if (!nLayers)
0327     throw cms::Exception("Setup") << "ME0ReDigiProducer::beginRun() - No layers in ME0 geometry.";
0328
0329   const unsigned int nPartitions = chambers.front()->layers()[0]->nEtaPartitions();
0330
0331   if (useCusGeoFor1PartGeo && nPartitions == 1) {
0332     useBuiltinGeo = false;
0333   }
0334
0335   if (useBuiltinGeo) {
0336     if (nLayers != layerReadout.size())
0337       throw cms::Exception("Configuration")
0338           << "ME0ReDigiProducer::beginRun() - The geometry has " << nLayers << " layers, but the readout of "
0339           << layerReadout.size() << " were specified with the layerReadout parameter.";
0340     fillCentralTOFs();
0341   } else {
0342     LogDebug("ME0ReDigiProducer") << "Building temporary geometry:" << std::endl;
0343     tempGeo = new TemporaryGeometry(geometry, numberOfStrips, numberOfPartitions);
0344     LogDebug("ME0ReDigiProducer") << "Done building temporary geometry!" << std::endl;
0345
0346     if (tempGeo->numLayers() != layerReadout.size())
0347       throw cms::Exception("Configuration")
0348           << "ME0ReDigiProducer::beginRun() - The geometry has " << tempGeo->numLayers()
0349           << " layers, but the readout of " << layerReadout.size()
0350           << " were specified with the layerReadout parameter.";
0351   }
0352 }
0353
0354 void ME0ReDigiProducer::produce(edm::Event& e, const edm::EventSetup& eventSetup) {
0355   edm::Service<edm::RandomNumberGenerator> rng;
0356   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0357
0358   edm::Handle<ME0DigiPreRecoCollection> input_digis;
0359   e.getByToken(token, input_digis);
0360
0361   std::unique_ptr<ME0DigiPreRecoCollection> output_digis(new ME0DigiPreRecoCollection());
0362   std::unique_ptr<ME0DigiPreRecoMap> output_digimap(new ME0DigiPreRecoMap());
0363
0364   // build the digis
0365   buildDigis(*(input_digis.product()), *output_digis, *output_digimap, engine);
0366
0367   // store them in the event
0368   e.put(std::move(output_digis));
0369   e.put(std::move(output_digimap));
0370 }
0371
0372 void ME0ReDigiProducer::buildDigis(const ME0DigiPreRecoCollection& input_digis,
0373                                    ME0DigiPreRecoCollection& output_digis,
0374                                    ME0DigiPreRecoMap& output_digimap,
0375                                    CLHEP::HepRandomEngine* engine) {
0376   /*
0377       Starting form the incoming pseudo-digi, which has perfect time and position resolution:
0378       1A. Smear time using sigma_t by some value
0379       1B. Correct the smeared time with the central arrival time for partition
0380       1C. Apply discretization: if the smeared time is outside the BX window (-12.5ns;+12.5ns),
0381       the hit should be assigned to the next (or previous) BX
0382
0383       2A. Find strip that the digi belongs to
0384       2B. Get the center of this strip and the error on the position assuming the geometry
0385
0386       3A. Filter event if a digi at this partition/strip/BX already exists
0387       3B. Add to collection
0388     */
0389
0390   LogDebug("ME0ReDigiProducer::buildDigis") << "Begin building digis." << std::endl;
0391   ME0DigiPreRecoCollection::DigiRangeIterator me0dgIt;
0392   for (me0dgIt = input_digis.begin(); me0dgIt != input_digis.end(); ++me0dgIt) {
0393     const auto& me0Id = (*me0dgIt).first;
0394     LogTrace("ME0ReDigiProducer::buildDigis") << "Starting with roll: " << me0Id << std::endl;
0395
0396     //setup map for this chamber/eta partition
0397     ChamberDigiMap chDigiMap;
0398
0399     int newDigiIdx = 0;
0400     const ME0DigiPreRecoCollection::Range& range = (*me0dgIt).second;
0401     for (ME0DigiPreRecoCollection::const_iterator digi = range.first; digi != range.second; digi++) {
0402       LogTrace("ME0ReDigiProducer::buildDigis") << std::endl
0403                                                 << "(" << digi->x() << "," << digi->y() << "," << digi->tof() << ","
0404                                                 << digi->pdgid() << "," << digi->prompt() << ")-> ";
0405
0406       //If we don't readout this layer skip
0407       if (!layerReadout[me0Id.layer() - 1]) {
0408         output_digimap.insertDigi(me0Id, -1);
0409         continue;
0410       }
0411
0412       //if neutron and we are filtering skip
0413       if (!digi->prompt() && neutronAcceptance < 1.0)
0414         if (CLHEP::RandFlat::shoot(engine) > neutronAcceptance) {
0415           output_digimap.insertDigi(me0Id, -1);
0416           continue;
0417         }
0418
0419       //smear TOF if necessary
0420       float tof = digi->tof() + (timeResolution < 0 ? 0.0 : CLHEP::RandGaussQ::shoot(engine, 0, timeResolution));
0421
0422       //Values used to fill objet
0423       int mapPartIDX = me0Id.roll() - 1;
0424       int strip = 0;
0425       LocalPoint digiLocalPoint;
0426       LocalError digiLocalError;
0427       if (useBuiltinGeo) {
0428         getStripProperties(geometry->etaPartition(me0Id), &*digi, tof, strip, digiLocalPoint, digiLocalError);
0429       } else {
0430         mapPartIDX = getCustomStripProperties(me0Id, &*digi, tof, strip, digiLocalPoint, digiLocalError);
0431       }
0432
0433       //filter if outside of readout window
0434       const int bxIdx = std::round(tof / bxWidth);
0435       LogTrace("ME0ReDigiProducer::buildDigis") << tof << "(" << bxIdx << ") ";
0436       if (bxIdx < minBXReadout) {
0437         output_digimap.insertDigi(me0Id, -1);
0438         continue;
0439       }
0440       if (bxIdx > maxBXReadout) {
0441         output_digimap.insertDigi(me0Id, -1);
0442         continue;
0443       }
0444       tof = bxIdx * bxWidth;
0445
0446       //If we are merging check to see if it already exists
0447       LogTrace("ME0ReDigiProducer::buildDigis") << "(" << bxIdx << "," << mapPartIDX << "," << strip << ") ";
0448       if (mergeDigis) {
0449         int matchIDX = fillDigiMap(chDigiMap, bxIdx, mapPartIDX, strip, newDigiIdx);
0450         if (matchIDX >= 0) {
0451           output_digimap.insertDigi(me0Id, matchIDX);
0452           continue;
0453         }
0454       }
0455
0456       //Digis store sigmaX,sigmaY, correlationCoef
0457       const float sigmaX = std::sqrt(digiLocalError.xx());
0458       const float sigmaY = std::sqrt(digiLocalError.yy());
0459       const float corrCoef = digiLocalError.xy() / (sigmaX * sigmaY);
0460
0461       //Fill in the new collection
0462       ME0DigiPreReco out_digi(
0463           digiLocalPoint.x(), digiLocalPoint.y(), sigmaX, sigmaY, corrCoef, tof, digi->pdgid(), digi->prompt());
0464       output_digis.insertDigi(me0Id, out_digi);
0465
0466       // store index of previous detid and digi
0467       output_digimap.insertDigi(me0Id, newDigiIdx);
0468       newDigiIdx++;
0469
0470       LogTrace("ME0ReDigiProducer::buildDigis") << "(" << digiLocalPoint.x() << "," << digiLocalPoint.y() << ","
0471                                                 << sigmaX << "," << sigmaY << "," << tof << ") ";
0472     }
0473
0474     chDigiMap.clear();
0475   }
0476 }
0477
0478 void ME0ReDigiProducer::fillCentralTOFs() {
0479   const auto* mainChamber = geometry->chambers().front();
0480   const unsigned int nLayers = mainChamber->nLayers();
0481   //Get TOF at center of each partition
0482   tofs.clear();
0483   tofs.resize(nLayers);
0484   LogDebug("ME0ReDigiProducer::fillCentralTOFs()") << "TOF numbers [layer][partition]: ";
0485   for (unsigned int iL = 0; iL < nLayers; ++iL) {
0486     const auto* layer = mainChamber->layers()[iL];
0487     const unsigned int mapLayIDX = layer->id().layer() - 1;
0488     const unsigned int nPartitions = layer->nEtaPartitions();
0489     if (!nPartitions)
0490       throw cms::Exception("Setup") << "ME0ReDigiProducer::fillCentralTOFs() - ME0Layer has no partitions.";
0491     tofs[mapLayIDX].resize(nPartitions);
0492     for (unsigned int iP = 0; iP < nPartitions; ++iP) {
0493       const unsigned int mapPartIDX = layer->etaPartitions()[iP]->id().roll() - 1;
0494       const GlobalPoint centralGP(layer->etaPartitions()[iP]->position());
0495       tofs[mapLayIDX][mapPartIDX] = (centralGP.mag() / (CLHEP::c_light / CLHEP::cm));  //speed of light [cm/ns]
0496       LogDebug("ME0ReDigiProducer::fillCentralTOFs()")
0497           << "[" << mapLayIDX << "][" << mapPartIDX << "]=" << tofs[mapLayIDX][mapPartIDX] << " " << std::endl;
0498     }
0499   }
0500 }
0501 int ME0ReDigiProducer::getCustomStripProperties(const ME0DetId& detId,
0502                                                 const ME0DigiPreReco* inDigi,
0503                                                 float& tof,
0504                                                 int& strip,
0505                                                 LocalPoint& digiLocalPoint,
0506                                                 LocalError& digiLocalError) const {
0507   const unsigned int partIdx = tempGeo->findEtaPartition(inDigi->y());
0508   LogTrace("ME0ReDigiProducer::buildDigis") << partIdx << " ";
0509   const float partMeanTof = tempGeo->getCentralTOF(detId, partIdx);
0510
0511   //convert to relative to partition
0512   tof -= partMeanTof;
0513
0514   //get coordinates and errors
0515   const float partCenter = tempGeo->getPartCenter(partIdx);
0516   const auto* topo = tempGeo->getTopo(partIdx);
0517
0518   //find channel
0519   const LocalPoint partLocalPoint(inDigi->x(), inDigi->y() - partCenter, 0.);
0520   strip = topo->channel(partLocalPoint);
0521   const float stripF = float(strip) + 0.5;
0522
0523   //get digitized location
0524   LocalPoint digiPartLocalPoint = topo->localPosition(stripF);
0525   digiLocalError = topo->localError(
0526       stripF,
0527       1. /
0528           sqrt(
0529               12.));  //std dev. flat distribution with length L is L/sqrt(12). The strip topology expects the error in units of strips.
0530   digiLocalPoint = LocalPoint(digiPartLocalPoint.x(), digiPartLocalPoint.y() + partCenter, 0.0);
0531   return partIdx;
0532 }
0533 void ME0ReDigiProducer::getStripProperties(const ME0EtaPartition* etaPart,
0534                                            const ME0DigiPreReco* inDigi,
0535                                            float& tof,
0536                                            int& strip,
0537                                            LocalPoint& digiLocalPoint,
0538                                            LocalError& digiLocalError) const {
0539   //convert to relative to partition
0540   tof -= tofs[etaPart->id().layer() - 1][etaPart->id().roll() - 1];
0541
0542   const TrapezoidalStripTopology* origTopo = (const TrapezoidalStripTopology*)(&etaPart->specificTopology());
0543   TrapezoidalStripTopology padTopo(
0544       origTopo->nstrips() / 2, origTopo->pitch() * 2, origTopo->stripLength(), origTopo->radius());
0545   const auto& topo = usePads ? padTopo : etaPart->specificTopology();
0546
0547   //find channel
0548   const LocalPoint partLocalPoint(inDigi->x(), inDigi->y(), 0.);
0549   strip = topo.channel(partLocalPoint);
0550   const float stripF = float(strip) + 0.5;
0551
0552   //get digitized location
0553   digiLocalPoint = topo.localPosition(stripF);
0554   digiLocalError = topo.localError(stripF, 1. / sqrt(12.));
0555 }
0556
0557 unsigned int ME0ReDigiProducer::fillDigiMap(
0558     ChamberDigiMap& chDigiMap, unsigned int bx, unsigned int part, unsigned int strip, unsigned int currentIDX) const {
0559   DigiIndicies newIDX(bx, part, strip);
0560   auto it1 = chDigiMap.find(newIDX);
0561   if (it1 == chDigiMap.end()) {
0562     chDigiMap[newIDX] = currentIDX;
0563     return -1;
0564   }
0565   return it1->second;
0566 }
0567
0568 DEFINE_FWK_MODULE(ME0ReDigiProducer);
0569 #endif