MtdValidation/plugins/BtlLocalRecoValidation.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    Validation/MtdValidation
0004 // Class:      BtlLocalRecoValidation
0005 //
0006 /**\class BtlLocalRecoValidation BtlLocalRecoValidation.cc Validation/MtdValidation/plugins/BtlLocalRecoValidation.cc
0007
0008  Description: BTL RECO hits and clusters validation
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013
0014 #include <string>
0015
0016 #include "FWCore/Framework/interface/Frameworkfwd.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/MakerMacros.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020
0021 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0022 #include "DQMServices/Core/interface/DQMStore.h"
0023
0024 #include "DataFormats/Common/interface/ValidHandle.h"
0025 #include "DataFormats/Math/interface/GeantUnits.h"
0026 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0027 #include "DataFormats/FTLRecHit/interface/FTLRecHitCollections.h"
0028 #include "DataFormats/FTLRecHit/interface/FTLClusterCollections.h"
0029 #include "DataFormats/TrackerRecHit2D/interface/MTDTrackingRecHit.h"
0030
0031 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
0032 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0033 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0034
0035 #include "Geometry/Records/interface/MTDDigiGeometryRecord.h"
0036 #include "Geometry/Records/interface/MTDTopologyRcd.h"
0037 #include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h"
0038 #include "Geometry/MTDNumberingBuilder/interface/MTDTopology.h"
0039
0040 #include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h"
0041 #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
0042
0043 #include "Geometry/MTDCommonData/interface/MTDTopologyMode.h"
0044
0045 #include "MTDHit.h"
0046
0047 class BtlLocalRecoValidation : public DQMEDAnalyzer {
0048 public:
0049   explicit BtlLocalRecoValidation(const edm::ParameterSet&);
0050   ~BtlLocalRecoValidation() override;
0051
0052   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0053
0054 private:
0055   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0056
0057   void analyze(const edm::Event&, const edm::EventSetup&) override;
0058
0059   bool isSameCluster(const FTLCluster&, const FTLCluster&);
0060
0061   // ------------ member data ------------
0062
0063   const std::string folder_;
0064   const double hitMinEnergy_;
0065   const bool optionalPlots_;
0066   const bool uncalibRecHitsPlots_;
0067   const double hitMinAmplitude_;
0068
0069   edm::EDGetTokenT<FTLRecHitCollection> btlRecHitsToken_;
0070   edm::EDGetTokenT<FTLUncalibratedRecHitCollection> btlUncalibRecHitsToken_;
0071   edm::EDGetTokenT<CrossingFrame<PSimHit> > btlSimHitsToken_;
0072   edm::EDGetTokenT<FTLClusterCollection> btlRecCluToken_;
0073   edm::EDGetTokenT<MTDTrackingDetSetVector> mtdTrackingHitToken_;
0074
0075   edm::ESGetToken<MTDGeometry, MTDDigiGeometryRecord> mtdgeoToken_;
0076   edm::ESGetToken<MTDTopology, MTDTopologyRcd> mtdtopoToken_;
0077
0078   // --- histograms declaration
0079
0080   MonitorElement* meNevents_;
0081
0082   MonitorElement* meNhits_;
0083
0084   MonitorElement* meHitEnergy_;
0085   MonitorElement* meHitLogEnergy_;
0086   MonitorElement* meHitTime_;
0087   MonitorElement* meHitTimeError_;
0088
0089   MonitorElement* meOccupancy_;
0090
0091   //local position monitoring
0092   MonitorElement* meLocalOccupancy_;
0093   MonitorElement* meHitXlocal_;
0094   MonitorElement* meHitYlocal_;
0095   MonitorElement* meHitZlocal_;
0096
0097   MonitorElement* meHitX_;
0098   MonitorElement* meHitY_;
0099   MonitorElement* meHitZ_;
0100   MonitorElement* meHitPhi_;
0101   MonitorElement* meHitEta_;
0102
0103   MonitorElement* meHitTvsE_;
0104   MonitorElement* meHitEvsPhi_;
0105   MonitorElement* meHitEvsEta_;
0106   MonitorElement* meHitEvsZ_;
0107   MonitorElement* meHitTvsPhi_;
0108   MonitorElement* meHitTvsEta_;
0109   MonitorElement* meHitTvsZ_;
0110   MonitorElement* meHitLongPos_;
0111   MonitorElement* meHitLongPosErr_;
0112
0113   MonitorElement* meTimeRes_;
0114   MonitorElement* meEnergyRes_;
0115
0116   MonitorElement* meLongPosPull_;
0117   MonitorElement* meLongPosPullvsE_;
0118   MonitorElement* meLongPosPullvsEta_;
0119
0120   MonitorElement* meTPullvsE_;
0121   MonitorElement* meTPullvsEta_;
0122
0123   MonitorElement* meNclusters_;
0124
0125   MonitorElement* meCluTime_;
0126   MonitorElement* meCluTimeError_;
0127   MonitorElement* meCluEnergy_;
0128   MonitorElement* meCluPhi_;
0129   MonitorElement* meCluEta_;
0130   MonitorElement* meCluHits_;
0131   MonitorElement* meCluZvsPhi_;
0132   MonitorElement* meCluEnergyvsEta_;
0133   MonitorElement* meCluHitsvsEta_;
0134
0135   MonitorElement* meCluTimeRes_;
0136   MonitorElement* meCluEnergyRes_;
0137   MonitorElement* meCluTResvsE_;
0138   MonitorElement* meCluTResvsEta_;
0139   MonitorElement* meCluTPullvsE_;
0140   MonitorElement* meCluTPullvsEta_;
0141   MonitorElement* meCluRhoRes_;
0142   MonitorElement* meCluPhiRes_;
0143   MonitorElement* meCluLocalXRes_;
0144   MonitorElement* meCluLocalYRes_;
0145   MonitorElement* meCluZRes_;
0146   MonitorElement* meCluLocalXPull_;
0147   MonitorElement* meCluLocalYPull_;
0148   MonitorElement* meCluZPull_;
0149   MonitorElement* meCluYXLocal_;
0150   MonitorElement* meCluYXLocalSim_;
0151   MonitorElement* meCluXLocalErr_;
0152   MonitorElement* meCluYLocalErr_;
0153
0154   MonitorElement* meUnmatchedCluEnergy_;
0155
0156   // --- UncalibratedRecHits histograms
0157
0158   static constexpr int nBinsQ_ = 20;
0159   static constexpr float binWidthQ_ = 30.;
0160   static constexpr int nBinsQEta_ = 3;
0161   static constexpr float binsQEta_[nBinsQEta_ + 1] = {0., 0.65, 1.15, 1.55};
0162
0163   MonitorElement* meTimeResQ_[nBinsQ_];
0164   MonitorElement* meTimeResQvsEta_[nBinsQ_][nBinsQEta_];
0165
0166   static constexpr int nBinsEta_ = 31;
0167   static constexpr float binWidthEta_ = 0.05;
0168   static constexpr int nBinsEtaQ_ = 7;
0169   static constexpr float binsEtaQ_[nBinsEtaQ_ + 1] = {0., 30., 60., 90., 120., 150., 360., 600.};
0170
0171   MonitorElement* meTimeResEta_[nBinsEta_];
0172   MonitorElement* meTimeResEtavsQ_[nBinsEta_][nBinsEtaQ_];
0173 };
0174
0175 bool BtlLocalRecoValidation::isSameCluster(const FTLCluster& clu1, const FTLCluster& clu2) {
0176   return clu1.id() == clu2.id() && clu1.size() == clu2.size() && clu1.x() == clu2.x() && clu1.y() == clu2.y() &&
0177          clu1.time() == clu2.time();
0178 }
0179
0180 // ------------ constructor and destructor --------------
0181 BtlLocalRecoValidation::BtlLocalRecoValidation(const edm::ParameterSet& iConfig)
0182     : folder_(iConfig.getParameter<std::string>("folder")),
0183       hitMinEnergy_(iConfig.getParameter<double>("HitMinimumEnergy")),
0184       optionalPlots_(iConfig.getParameter<bool>("optionalPlots")),
0185       uncalibRecHitsPlots_(iConfig.getParameter<bool>("UncalibRecHitsPlots")),
0186       hitMinAmplitude_(iConfig.getParameter<double>("HitMinimumAmplitude")) {
0187   btlRecHitsToken_ = consumes<FTLRecHitCollection>(iConfig.getParameter<edm::InputTag>("recHitsTag"));
0188   if (uncalibRecHitsPlots_)
0189     btlUncalibRecHitsToken_ =
0190         consumes<FTLUncalibratedRecHitCollection>(iConfig.getParameter<edm::InputTag>("uncalibRecHitsTag"));
0191   btlSimHitsToken_ = consumes<CrossingFrame<PSimHit> >(iConfig.getParameter<edm::InputTag>("simHitsTag"));
0192   btlRecCluToken_ = consumes<FTLClusterCollection>(iConfig.getParameter<edm::InputTag>("recCluTag"));
0193   mtdTrackingHitToken_ = consumes<MTDTrackingDetSetVector>(iConfig.getParameter<edm::InputTag>("trkHitTag"));
0194
0195   mtdgeoToken_ = esConsumes<MTDGeometry, MTDDigiGeometryRecord>();
0196   mtdtopoToken_ = esConsumes<MTDTopology, MTDTopologyRcd>();
0197 }
0198
0199 BtlLocalRecoValidation::~BtlLocalRecoValidation() {}
0200
0201 // ------------ method called for each event  ------------
0202 void BtlLocalRecoValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0203   using namespace edm;
0204   using namespace std;
0205   using namespace geant_units::operators;
0206
0207   auto geometryHandle = iSetup.getTransientHandle(mtdgeoToken_);
0208   const MTDGeometry* geom = geometryHandle.product();
0209
0210   auto topologyHandle = iSetup.getTransientHandle(mtdtopoToken_);
0211   const MTDTopology* topology = topologyHandle.product();
0212
0213   auto btlRecHitsHandle = makeValid(iEvent.getHandle(btlRecHitsToken_));
0214   auto btlSimHitsHandle = makeValid(iEvent.getHandle(btlSimHitsToken_));
0215   auto btlRecCluHandle = makeValid(iEvent.getHandle(btlRecCluToken_));
0216   auto mtdTrkHitHandle = makeValid(iEvent.getHandle(mtdTrackingHitToken_));
0217   MixCollection<PSimHit> btlSimHits(btlSimHitsHandle.product());
0218
0219 #ifdef EDM_ML_DEBUG
0220   for (const auto& hits : *mtdTrkHitHandle) {
0221     if (MTDDetId(hits.id()).mtdSubDetector() == MTDDetId::MTDType::BTL) {
0222       LogDebug("BtlLocalRecoValidation") << "MTD cluster DetId " << hits.id() << " # cluster " << hits.size();
0223       for (const auto& hit : hits) {
0224         LogDebug("BtlLocalRecoValidation")
0225             << "MTD_TRH: " << hit.localPosition().x() << "," << hit.localPosition().y() << " : "
0226             << hit.localPositionError().xx() << "," << hit.localPositionError().yy() << " : " << hit.time() << " : "
0227             << hit.timeError();
0228       }
0229     }
0230   }
0231 #endif
0232
0233   // --- Loop over the BTL SIM hits
0234   std::unordered_map<uint32_t, MTDHit> m_btlSimHits;
0235   for (auto const& simHit : btlSimHits) {
0236     // --- Use only hits compatible with the in-time bunch-crossing
0237     if (simHit.tof() < 0 || simHit.tof() > 25.)
0238       continue;
0239
0240     DetId id = simHit.detUnitId();
0241
0242     auto simHitIt = m_btlSimHits.emplace(id.rawId(), MTDHit()).first;
0243
0244     // --- Accumulate the energy (in MeV) of SIM hits in the same detector cell
0245     (simHitIt->second).energy += convertUnitsTo(0.001_MeV, simHit.energyLoss());
0246
0247     // --- Get the time of the first SIM hit in the cell
0248     if ((simHitIt->second).time == 0 || simHit.tof() < (simHitIt->second).time) {
0249       (simHitIt->second).time = simHit.tof();
0250
0251       auto hit_pos = simHit.entryPoint();
0252       (simHitIt->second).x = hit_pos.x();
0253       (simHitIt->second).y = hit_pos.y();
0254       (simHitIt->second).z = hit_pos.z();
0255     }
0256
0257   }  // simHit loop
0258
0259   // --- Loop over the BTL RECO hits
0260   unsigned int n_reco_btl = 0;
0261   for (const auto& recHit : *btlRecHitsHandle) {
0262     BTLDetId detId = recHit.id();
0263     DetId geoId = detId.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topology->getMTDTopologyMode()));
0264     const MTDGeomDet* thedet = geom->idToDet(geoId);
0265     if (thedet == nullptr)
0266       throw cms::Exception("BtlLocalRecoValidation") << "GeographicalID: " << std::hex << geoId.rawId() << " ("
0267                                                      << detId.rawId() << ") is invalid!" << std::dec << std::endl;
0268     const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
0269     const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0270
0271     Local3DPoint local_point(0., 0., 0.);
0272     local_point = topo.pixelToModuleLocalPoint(local_point, detId.row(topo.nrows()), detId.column(topo.nrows()));
0273     const auto& global_point = thedet->toGlobal(local_point);
0274
0275     meHitEnergy_->Fill(recHit.energy());
0276     meHitLogEnergy_->Fill(log10(recHit.energy()));
0277     meHitTime_->Fill(recHit.time());
0278     meHitTimeError_->Fill(recHit.timeError());
0279     meHitLongPos_->Fill(recHit.position());
0280     meHitLongPosErr_->Fill(recHit.positionError());
0281
0282     meOccupancy_->Fill(global_point.z(), global_point.phi());
0283
0284     if (optionalPlots_) {
0285       meLocalOccupancy_->Fill(local_point.x() + recHit.position(), local_point.y());
0286       meHitXlocal_->Fill(local_point.x());
0287       meHitYlocal_->Fill(local_point.y());
0288       meHitZlocal_->Fill(local_point.z());
0289     }
0290     meHitX_->Fill(global_point.x());
0291     meHitY_->Fill(global_point.y());
0292     meHitZ_->Fill(global_point.z());
0293     meHitPhi_->Fill(global_point.phi());
0294     meHitEta_->Fill(global_point.eta());
0295
0296     meHitTvsE_->Fill(recHit.energy(), recHit.time());
0297     meHitEvsPhi_->Fill(global_point.phi(), recHit.energy());
0298     meHitEvsEta_->Fill(global_point.eta(), recHit.energy());
0299     meHitEvsZ_->Fill(global_point.z(), recHit.energy());
0300     meHitTvsPhi_->Fill(global_point.phi(), recHit.time());
0301     meHitTvsEta_->Fill(global_point.eta(), recHit.time());
0302     meHitTvsZ_->Fill(global_point.z(), recHit.time());
0303
0304     // Resolution histograms
0305     if (m_btlSimHits.count(detId.rawId()) == 1 && m_btlSimHits[detId.rawId()].energy > hitMinEnergy_) {
0306       float longpos_res = recHit.position() - convertMmToCm(m_btlSimHits[detId.rawId()].x);
0307       float time_res = recHit.time() - m_btlSimHits[detId.rawId()].time;
0308       float energy_res = recHit.energy() - m_btlSimHits[detId.rawId()].energy;
0309
0310       Local3DPoint local_point_sim(convertMmToCm(m_btlSimHits[detId.rawId()].x),
0311                                    convertMmToCm(m_btlSimHits[detId.rawId()].y),
0312                                    convertMmToCm(m_btlSimHits[detId.rawId()].z));
0313       local_point_sim =
0314           topo.pixelToModuleLocalPoint(local_point_sim, detId.row(topo.nrows()), detId.column(topo.nrows()));
0315       const auto& global_point_sim = thedet->toGlobal(local_point_sim);
0316
0317       meTimeRes_->Fill(time_res);
0318       meEnergyRes_->Fill(energy_res);
0319
0320       meLongPosPull_->Fill(longpos_res / recHit.positionError());
0321       meLongPosPullvsEta_->Fill(std::abs(global_point_sim.eta()), longpos_res / recHit.positionError());
0322       meLongPosPullvsE_->Fill(m_btlSimHits[detId.rawId()].energy, longpos_res / recHit.positionError());
0323
0324       meTPullvsEta_->Fill(std::abs(global_point_sim.eta()), time_res / recHit.timeError());
0325       meTPullvsE_->Fill(m_btlSimHits[detId.rawId()].energy, time_res / recHit.timeError());
0326     }
0327
0328     n_reco_btl++;
0329
0330   }  // recHit loop
0331
0332   if (n_reco_btl > 0)
0333     meNhits_->Fill(log10(n_reco_btl));
0334
0335   // --- Loop over the BTL RECO clusters ---
0336   unsigned int n_clus_btl(0);
0337   for (const auto& DetSetClu : *btlRecCluHandle) {
0338     for (const auto& cluster : DetSetClu) {
0339       if (cluster.energy() < hitMinEnergy_)
0340         continue;
0341       BTLDetId cluId = cluster.id();
0342       DetId detIdObject(cluId);
0343       const auto& genericDet = geom->idToDetUnit(detIdObject);
0344       if (genericDet == nullptr) {
0345         throw cms::Exception("BtlLocalRecoValidation")
0346             << "GeographicalID: " << std::hex << cluId << " is invalid!" << std::dec << std::endl;
0347       }
0348       n_clus_btl++;
0349
0350       const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(genericDet->topology());
0351       const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0352
0353       // --- Cluster position in the module reference frame
0354       Local3DPoint local_point(topo.localX(cluster.x()), topo.localY(cluster.y()), 0.);
0355       const auto& global_point = genericDet->toGlobal(local_point);
0356
0357       meCluEnergy_->Fill(cluster.energy());
0358       meCluTime_->Fill(cluster.time());
0359       meCluTimeError_->Fill(cluster.timeError());
0360       meCluPhi_->Fill(global_point.phi());
0361       meCluEta_->Fill(global_point.eta());
0362       meCluZvsPhi_->Fill(global_point.z(), global_point.phi());
0363       meCluHits_->Fill(cluster.size());
0364
0365       // --- Get the SIM hits associated to the cluster and calculate
0366       //     the cluster SIM energy, time and position
0367
0368       double cluEneSIM = 0.;
0369       double cluTimeSIM = 0.;
0370       double cluLocXSIM = 0.;
0371       double cluLocYSIM = 0.;
0372       double cluLocZSIM = 0.;
0373
0374       for (int ihit = 0; ihit < cluster.size(); ++ihit) {
0375         int hit_row = cluster.minHitRow() + cluster.hitOffset()[ihit * 2];
0376         int hit_col = cluster.minHitCol() + cluster.hitOffset()[ihit * 2 + 1];
0377
0378         // Match the RECO hit to the corresponding SIM hit
0379         for (const auto& recHit : *btlRecHitsHandle) {
0380           BTLDetId hitId(recHit.id().rawId());
0381
0382           if (m_btlSimHits.count(hitId.rawId()) == 0)
0383             continue;
0384
0385           // Check the hit position
0386           if (hitId.mtdSide() != cluId.mtdSide() || hitId.mtdRR() != cluId.mtdRR() || recHit.row() != hit_row ||
0387               recHit.column() != hit_col)
0388             continue;
0389
0390           // Check the hit energy and time
0391           if (recHit.energy() != cluster.hitENERGY()[ihit] || recHit.time() != cluster.hitTIME()[ihit])
0392             continue;
0393
0394           // SIM hit's position in the module reference frame
0395           Local3DPoint local_point_sim(convertMmToCm(m_btlSimHits[recHit.id().rawId()].x),
0396                                        convertMmToCm(m_btlSimHits[recHit.id().rawId()].y),
0397                                        convertMmToCm(m_btlSimHits[recHit.id().rawId()].z));
0398           local_point_sim =
0399               topo.pixelToModuleLocalPoint(local_point_sim, hitId.row(topo.nrows()), hitId.column(topo.nrows()));
0400
0401           // Calculate the SIM cluster's position in the module reference frame
0402           cluLocXSIM += local_point_sim.x() * m_btlSimHits[recHit.id().rawId()].energy;
0403           cluLocYSIM += local_point_sim.y() * m_btlSimHits[recHit.id().rawId()].energy;
0404           cluLocZSIM += local_point_sim.z() * m_btlSimHits[recHit.id().rawId()].energy;
0405
0406           // Calculate the SIM cluster energy and time
0407           cluEneSIM += m_btlSimHits[recHit.id().rawId()].energy;
0408           cluTimeSIM += m_btlSimHits[recHit.id().rawId()].time * m_btlSimHits[recHit.id().rawId()].energy;
0409
0410           break;
0411
0412         }  // recHit loop
0413
0414       }  // ihit loop
0415
0416       // Find the MTDTrackingRecHit corresponding to the cluster
0417       const MTDTrackingRecHit* comp(nullptr);
0418       bool matchClu = false;
0419       const auto& trkHits = (*mtdTrkHitHandle)[detIdObject];
0420       for (const auto& trkHit : trkHits) {
0421         if (isSameCluster(trkHit.mtdCluster(), cluster)) {
0422           comp = trkHit.clone();
0423           matchClu = true;
0424           break;
0425         }
0426       }
0427       if (!matchClu) {
0428         edm::LogWarning("BtlLocalRecoValidation")
0429             << "No valid TrackingRecHit corresponding to cluster, detId = " << detIdObject.rawId();
0430       }
0431
0432       // --- Fill the cluster resolution histograms
0433       if (cluTimeSIM > 0. && cluEneSIM > 0.) {
0434         cluTimeSIM /= cluEneSIM;
0435
0436         Local3DPoint cluLocalPosSIM(cluLocXSIM / cluEneSIM, cluLocYSIM / cluEneSIM, cluLocZSIM / cluEneSIM);
0437         const auto& cluGlobalPosSIM = genericDet->toGlobal(cluLocalPosSIM);
0438
0439         float time_res = cluster.time() - cluTimeSIM;
0440         float energy_res = cluster.energy() - cluEneSIM;
0441         meCluTimeRes_->Fill(time_res);
0442         meCluEnergyRes_->Fill(energy_res);
0443
0444         float rho_res = global_point.perp() - cluGlobalPosSIM.perp();
0445         float phi_res = global_point.phi() - cluGlobalPosSIM.phi();
0446
0447         meCluRhoRes_->Fill(rho_res);
0448         meCluPhiRes_->Fill(phi_res);
0449
0450         float xlocal_res = local_point.x() - cluLocalPosSIM.x();
0451         float ylocal_res = local_point.y() - cluLocalPosSIM.y();
0452         float z_res = global_point.z() - cluGlobalPosSIM.z();
0453
0454         meCluZRes_->Fill(z_res);
0455
0456         if (optionalPlots_) {
0457           if (matchClu && comp != nullptr) {
0458             meCluLocalXRes_->Fill(xlocal_res);
0459             meCluLocalYRes_->Fill(ylocal_res);
0460             meCluLocalXPull_->Fill(xlocal_res / std::sqrt(comp->localPositionError().xx()));
0461             meCluLocalYPull_->Fill(ylocal_res / std::sqrt(comp->localPositionError().yy()));
0462             meCluZPull_->Fill(z_res / std::sqrt(comp->globalPositionError().czz()));
0463             meCluXLocalErr_->Fill(std::sqrt(comp->localPositionError().xx()));
0464             meCluYLocalErr_->Fill(std::sqrt(comp->localPositionError().yy()));
0465           }
0466
0467           meCluYXLocal_->Fill(local_point.x(), local_point.y());
0468           meCluYXLocalSim_->Fill(cluLocalPosSIM.x(), cluLocalPosSIM.y());
0469         }
0470
0471         meCluEnergyvsEta_->Fill(std::abs(cluGlobalPosSIM.eta()), cluster.energy());
0472         meCluHitsvsEta_->Fill(std::abs(cluGlobalPosSIM.eta()), cluster.size());
0473
0474         meCluTResvsEta_->Fill(std::abs(cluGlobalPosSIM.eta()), time_res);
0475         meCluTResvsE_->Fill(cluEneSIM, time_res);
0476
0477         meCluTPullvsEta_->Fill(std::abs(cluGlobalPosSIM.eta()), time_res / cluster.timeError());
0478         meCluTPullvsE_->Fill(cluEneSIM, time_res / cluster.timeError());
0479
0480       }  // if ( cluTimeSIM > 0. &&  cluEneSIM > 0. )
0481       else {
0482         meUnmatchedCluEnergy_->Fill(std::log10(cluster.energy()));
0483       }
0484
0485     }  // cluster loop
0486
0487   }  // DetSetClu loop
0488
0489   if (n_clus_btl > 0)
0490     meNclusters_->Fill(log10(n_clus_btl));
0491
0492   // --- This is to count the number of processed events, needed in the harvesting step
0493   meNevents_->Fill(0.5);
0494
0495   // --- Loop over the BTL Uncalibrated RECO hits
0496   if (uncalibRecHitsPlots_) {
0497     auto btlUncalibRecHitsHandle = makeValid(iEvent.getHandle(btlUncalibRecHitsToken_));
0498
0499     for (const auto& uRecHit : *btlUncalibRecHitsHandle) {
0500       BTLDetId detId = uRecHit.id();
0501
0502       // --- Skip UncalibratedRecHits not matched to SimHits
0503       if (m_btlSimHits.count(detId.rawId()) != 1)
0504         continue;
0505
0506       DetId geoId = detId.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topology->getMTDTopologyMode()));
0507       const MTDGeomDet* thedet = geom->idToDet(geoId);
0508       if (thedet == nullptr)
0509         throw cms::Exception("BtlLocalRecoValidation") << "GeographicalID: " << std::hex << geoId.rawId() << " ("
0510                                                        << detId.rawId() << ") is invalid!" << std::dec << std::endl;
0511       const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
0512       const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0513
0514       Local3DPoint local_point(0., 0., 0.);
0515       local_point = topo.pixelToModuleLocalPoint(local_point, detId.row(topo.nrows()), detId.column(topo.nrows()));
0516       const auto& global_point = thedet->toGlobal(local_point);
0517
0518       // --- Combine the information from the left and right BTL cell sides
0519
0520       float nHits = 0.;
0521       float hit_amplitude = 0.;
0522       float hit_time = 0.;
0523
0524       // left side:
0525       if (uRecHit.amplitude().first > 0.) {
0526         hit_amplitude += uRecHit.amplitude().first;
0527         hit_time += uRecHit.time().first;
0528         nHits += 1.;
0529       }
0530       // right side:
0531       if (uRecHit.amplitude().second > 0.) {
0532         hit_amplitude += uRecHit.amplitude().second;
0533         hit_time += uRecHit.time().second;
0534         nHits += 1.;
0535       }
0536
0537       hit_amplitude /= nHits;
0538       hit_time /= nHits;
0539
0540       // --- Fill the histograms
0541
0542       if (hit_amplitude < hitMinAmplitude_)
0543         continue;
0544
0545       float time_res = hit_time - m_btlSimHits[detId.rawId()].time;
0546
0547       // amplitude histograms
0548
0549       int qBin = (int)(hit_amplitude / binWidthQ_);
0550       if (qBin > nBinsQ_ - 1)
0551         qBin = nBinsQ_ - 1;
0552
0553       meTimeResQ_[qBin]->Fill(time_res);
0554
0555       int etaBin = 0;
0556       for (int ibin = 1; ibin < nBinsQEta_; ++ibin)
0557         if (fabs(global_point.eta()) >= binsQEta_[ibin] && fabs(global_point.eta()) < binsQEta_[ibin + 1])
0558           etaBin = ibin;
0559
0560       meTimeResQvsEta_[qBin][etaBin]->Fill(time_res);
0561
0562       // eta histograms
0563
0564       etaBin = (int)(fabs(global_point.eta()) / binWidthEta_);
0565       if (etaBin > nBinsEta_ - 1)
0566         etaBin = nBinsEta_ - 1;
0567
0568       meTimeResEta_[etaBin]->Fill(time_res);
0569
0570       qBin = 0;
0571       for (int ibin = 1; ibin < nBinsEtaQ_; ++ibin)
0572         if (hit_amplitude >= binsEtaQ_[ibin] && hit_amplitude < binsEtaQ_[ibin + 1])
0573           qBin = ibin;
0574
0575       meTimeResEtavsQ_[etaBin][qBin]->Fill(time_res);
0576
0577     }  // uRecHit loop
0578   }
0579 }
0580
0581 // ------------ method for histogram booking ------------
0582 void BtlLocalRecoValidation::bookHistograms(DQMStore::IBooker& ibook,
0583                                             edm::Run const& run,
0584                                             edm::EventSetup const& iSetup) {
0585   ibook.setCurrentFolder(folder_);
0586
0587   // --- histograms booking
0588
0589   meNevents_ = ibook.book1D("BtlNevents", "Number of events", 1, 0., 1.);
0590
0591   meNhits_ = ibook.book1D("BtlNhits", "Number of BTL RECO hits;log_{10}(N_{RECO})", 100, 0., 5.25);
0592
0593   meHitEnergy_ = ibook.book1D("BtlHitEnergy", "BTL RECO hits energy;E_{RECO} [MeV]", 100, 0., 20.);
0594   meHitLogEnergy_ = ibook.book1D("BtlHitLogEnergy", "BTL RECO hits energy;log_{10}(E_{RECO} [MeV])", 25, -1., 1.5);
0595   meHitTime_ = ibook.book1D("BtlHitTime", "BTL RECO hits ToA;ToA_{RECO} [ns]", 100, 0., 25.);
0596   meHitTimeError_ = ibook.book1D("BtlHitTimeError", "BTL RECO hits ToA error;#sigma^{ToA}_{RECO} [ns]", 50, 0., 0.1);
0597   meOccupancy_ = ibook.book2D(
0598       "BtlOccupancy", "BTL RECO hits occupancy;Z_{RECO} [cm]; #phi_{RECO} [rad]", 65, -260., 260., 126, -3.2, 3.2);
0599   if (optionalPlots_) {
0600     meLocalOccupancy_ = ibook.book2D(
0601         "BtlLocalOccupancy", "BTL RECO hits local occupancy;X_{RECO} [cm]; Y_{RECO} [cm]", 100, 10., 10., 60, -3., 3.);
0602     meHitXlocal_ = ibook.book1D("BtlHitXlocal", "BTL RECO local X;X_{RECO}^{LOC} [cm]", 100, -10., 10.);
0603     meHitYlocal_ = ibook.book1D("BtlHitYlocal", "BTL RECO local Y;Y_{RECO}^{LOC} [cm]", 60, -3, 3);
0604     meHitZlocal_ = ibook.book1D("BtlHitZlocal", "BTL RECO local z;z_{RECO}^{LOC} [cm]", 10, -1, 1);
0605   }
0606   meHitX_ = ibook.book1D("BtlHitX", "BTL RECO hits X;X_{RECO} [cm]", 60, -120., 120.);
0607   meHitY_ = ibook.book1D("BtlHitY", "BTL RECO hits Y;Y_{RECO} [cm]", 60, -120., 120.);
0608   meHitZ_ = ibook.book1D("BtlHitZ", "BTL RECO hits Z;Z_{RECO} [cm]", 100, -260., 260.);
0609   meHitPhi_ = ibook.book1D("BtlHitPhi", "BTL RECO hits #phi;#phi_{RECO} [rad]", 126, -3.2, 3.2);
0610   meHitEta_ = ibook.book1D("BtlHitEta", "BTL RECO hits #eta;#eta_{RECO}", 100, -1.55, 1.55);
0611   meHitTvsE_ =
0612       ibook.bookProfile("BtlHitTvsE", "BTL RECO ToA vs energy;E_{RECO} [MeV];ToA_{RECO} [ns]", 50, 0., 20., 0., 100.);
0613   meHitEvsPhi_ = ibook.bookProfile(
0614       "BtlHitEvsPhi", "BTL RECO energy vs #phi;#phi_{RECO} [rad];E_{RECO} [MeV]", 50, -3.2, 3.2, 0., 100.);
0615   meHitEvsEta_ = ibook.bookProfile(
0616       "BtlHitEvsEta", "BTL RECO energy vs #eta;#eta_{RECO};E_{RECO} [MeV]", 50, -1.55, 1.55, 0., 100.);
0617   meHitEvsZ_ =
0618       ibook.bookProfile("BtlHitEvsZ", "BTL RECO energy vs Z;Z_{RECO} [cm];E_{RECO} [MeV]", 50, -260., 260., 0., 100.);
0619   meHitTvsPhi_ = ibook.bookProfile(
0620       "BtlHitTvsPhi", "BTL RECO ToA vs #phi;#phi_{RECO} [rad];ToA_{RECO} [ns]", 50, -3.2, 3.2, 0., 100.);
0621   meHitTvsEta_ =
0622       ibook.bookProfile("BtlHitTvsEta", "BTL RECO ToA vs #eta;#eta_{RECO};ToA_{RECO} [ns]", 50, -1.6, 1.6, 0., 100.);
0623   meHitTvsZ_ =
0624       ibook.bookProfile("BtlHitTvsZ", "BTL RECO ToA vs Z;Z_{RECO} [cm];ToA_{RECO} [ns]", 50, -260., 260., 0., 100.);
0625   meHitLongPos_ = ibook.book1D("BtlLongPos", "BTL RECO hits longitudinal position;long. pos._{RECO}", 100, -10, 10);
0626   meHitLongPosErr_ =
0627       ibook.book1D("BtlLongPosErr", "BTL RECO hits longitudinal position error; long. pos. error_{RECO}", 100, -1, 1);
0628   meTimeRes_ = ibook.book1D("BtlTimeRes", "BTL time resolution;T_{RECO}-T_{SIM}", 100, -0.5, 0.5);
0629   meEnergyRes_ = ibook.book1D("BtlEnergyRes", "BTL energy resolution;E_{RECO}-E_{SIM}", 100, -0.5, 0.5);
0630   meLongPosPull_ = ibook.book1D("BtlLongPosPull",
0631                                 "BTL longitudinal position pull;X^{loc}_{RECO}-X^{loc}_{SIM}/#sigma_{xloc_{RECO}}",
0632                                 100,
0633                                 -5.,
0634                                 5.);
0635   meLongPosPullvsE_ = ibook.bookProfile(
0636       "BtlLongposPullvsE",
0637       "BTL longitudinal position pull vs E;E_{SIM} [MeV];X^{loc}_{RECO}-X^{loc}_{SIM}/#sigma_{xloc_{RECO}}",
0638       20,
0639       0.,
0640       20.,
0641       -5.,
0642       5.,
0643       "S");
0644   meLongPosPullvsEta_ = ibook.bookProfile(
0645       "BtlLongposPullvsEta",
0646       "BTL longitudinal position pull vs #eta;|#eta_{RECO}|;X^{loc}_{RECO}-X^{loc}_{SIM}/#sigma_{xloc_{RECO}}",
0647       32,
0648       0,
0649       1.55,
0650       -5.,
0651       5.,
0652       "S");
0653   meTPullvsE_ = ibook.bookProfile(
0654       "BtlTPullvsE", "BTL time pull vs E;E_{SIM} [MeV];(T_{RECO}-T_{SIM})/#sigma_{T_{RECO}}", 20, 0., 20., -5., 5., "S");
0655   meTPullvsEta_ = ibook.bookProfile("BtlTPullvsEta",
0656                                     "BTL time pull vs #eta;|#eta_{RECO}|;(T_{RECO}-T_{SIM})/#sigma_{T_{RECO}}",
0657                                     30,
0658                                     0,
0659                                     1.55,
0660                                     -5.,
0661                                     5.,
0662                                     "S");
0663
0664   meNclusters_ = ibook.book1D("BtlNclusters", "Number of BTL RECO clusters;log_{10}(N_{RECO})", 100, 0., 5.25);
0665   meCluTime_ = ibook.book1D("BtlCluTime", "BTL cluster time ToA;ToA [ns]", 250, 0, 25);
0666   meCluTimeError_ = ibook.book1D("BtlCluTimeError", "BTL cluster time error;#sigma_{t} [ns]", 100, 0, 0.1);
0667   meCluEnergy_ = ibook.book1D("BtlCluEnergy", "BTL cluster energy;E_{RECO} [MeV]", 100, 0, 20);
0668   meCluPhi_ = ibook.book1D("BtlCluPhi", "BTL cluster #phi;#phi_{RECO} [rad]", 144, -3.2, 3.2);
0669   meCluEta_ = ibook.book1D("BtlCluEta", "BTL cluster #eta;#eta_{RECO}", 100, -1.55, 1.55);
0670   meCluHits_ = ibook.book1D("BtlCluHitNumber", "BTL hits per cluster; Cluster size", 10, 0, 10);
0671   meCluZvsPhi_ = ibook.book2D(
0672       "BtlOccupancy", "BTL cluster Z vs #phi;Z_{RECO} [cm]; #phi_{RECO} [rad]", 144, -260., 260., 50, -3.2, 3.2);
0673   meCluEnergyvsEta_ = ibook.bookProfile(
0674       "BtlCluEnergyVsEta", "BTL cluster energy vs #eta; |#eta_{RECO}|; E_{RECO} [cm]", 30, 0., 1.55, 0., 20., "S");
0675   meCluHitsvsEta_ = ibook.bookProfile(
0676       "BtlCluHitsVsEta", "BTL hits per cluster vs #eta; |#eta_{RECO}|;Cluster size", 30, 0., 1.55, 0., 10., "S");
0677
0678   meCluTimeRes_ = ibook.book1D("BtlCluTimeRes", "BTL cluster time resolution;T_{RECO}-T_{SIM} [ns]", 100, -0.5, 0.5);
0679   meCluEnergyRes_ =
0680       ibook.book1D("BtlCluEnergyRes", "BTL cluster energy resolution;E_{RECO}-E_{SIM} [MeV]", 100, -0.5, 0.5);
0681   meCluTResvsE_ = ibook.bookProfile("BtlCluTResvsE",
0682                                     "BTL cluster time resolution vs E;E_{SIM} [MeV];(T_{RECO}-T_{SIM}) [ns]",
0683                                     20,
0684                                     0.,
0685                                     20.,
0686                                     -0.5,
0687                                     0.5,
0688                                     "S");
0689   meCluTResvsEta_ = ibook.bookProfile("BtlCluTResvsEta",
0690                                       "BTL cluster time resolution vs #eta;|#eta_{RECO}|;(T_{RECO}-T_{SIM}) [ns]",
0691                                       30,
0692                                       0,
0693                                       1.55,
0694                                       -0.5,
0695                                       0.5,
0696                                       "S");
0697   meCluTPullvsE_ = ibook.bookProfile("BtlCluTPullvsE",
0698                                      "BTL cluster time pull vs E;E_{SIM} [MeV];(T_{RECO}-T_{SIM})/#sigma_{T_{RECO}}",
0699                                      20,
0700                                      0.,
0701                                      20.,
0702                                      -5.,
0703                                      5.,
0704                                      "S");
0705   meCluTPullvsEta_ =
0706       ibook.bookProfile("BtlCluTPullvsEta",
0707                         "BTL cluster time pull vs #eta;|#eta_{RECO}|;(T_{RECO}-T_{SIM})/#sigma_{T_{RECO}}",
0708                         30,
0709                         0,
0710                         1.55,
0711                         -5.,
0712                         5.,
0713                         "S");
0714   meCluRhoRes_ =
0715       ibook.book1D("BtlCluRhoRes", "BTL cluster #rho resolution;#rho_{RECO}-#rho_{SIM} [cm]", 100, -0.5, 0.5);
0716   meCluPhiRes_ =
0717       ibook.book1D("BtlCluPhiRes", "BTL cluster #phi resolution;#phi_{RECO}-#phi_{SIM} [rad]", 100, -0.03, 0.03);
0718   meCluZRes_ = ibook.book1D("BtlCluZRes", "BTL cluster Z resolution;Z_{RECO}-Z_{SIM} [cm]", 100, -0.2, 0.2);
0719   if (optionalPlots_) {
0720     meCluLocalXRes_ =
0721         ibook.book1D("BtlCluLocalXRes", "BTL cluster local X resolution;X_{RECO}-X_{SIM} [cm]", 100, -3.1, 3.1);
0722     meCluLocalYRes_ =
0723         ibook.book1D("BtlCluLocalYRes", "BTL cluster local Y resolution;Y_{RECO}-Y_{SIM} [cm]", 100, -3.1, 3.1);
0724     meCluLocalXPull_ =
0725         ibook.book1D("BtlCluLocalXPull", "BTL cluster local X pull;X_{RECO}-X_{SIM}/sigmaX_[RECO]", 100, -5., 5.);
0726     meCluLocalYPull_ =
0727         ibook.book1D("BtlCluLocalYPull", "BTL cluster local Y pull;Y_{RECO}-Y_{SIM}/sigmaY_[RECO]", 100, -5., 5.);
0728     meCluZPull_ = ibook.book1D("BtlCluZPull", "BTL cluster Z pull;Z_{RECO}-Z_{SIM}/sigmaZ_[RECO]", 100, -5., 5.);
0729     meCluYXLocal_ = ibook.book2D("BtlCluYXLocal",
0730                                  "BTL cluster local Y vs X;X^{local}_{RECO} [cm];Y^{local}_{RECO} [cm]",
0731                                  200,
0732                                  -9.5,
0733                                  9.5,
0734                                  200,
0735                                  -2.8,
0736                                  2.8);
0737     meCluYXLocalSim_ = ibook.book2D("BtlCluYXLocalSim",
0738                                     "BTL cluster local Y vs X;X^{local}_{SIM} [cm];Y^{local}_{SIM} [cm]",
0739                                     200,
0740                                     -9.5,
0741                                     9.5,
0742                                     200,
0743                                     -2.8,
0744                                     2.8);
0745     meCluXLocalErr_ = ibook.book1D("BtlCluXLocalErr", "BTL cluster X local error;sigmaX_{RECO,loc} [cm]", 30, 0., 3.);
0746     meCluYLocalErr_ = ibook.book1D("BtlCluYLocalErr", "BTL cluster Y local error;sigmaY_{RECO,loc} [cm]", 30, 0., 0.9);
0747   }
0748   meUnmatchedCluEnergy_ =
0749       ibook.book1D("BtlUnmatchedCluEnergy", "BTL unmatched cluster log10(energy);log10(E_{RECO} [MeV])", 5, -3, 2);
0750
0751   // --- UncalibratedRecHits histograms
0752
0753   if (uncalibRecHitsPlots_) {
0754     for (unsigned int ihistoQ = 0; ihistoQ < nBinsQ_; ++ihistoQ) {
0755       std::string hname = Form("TimeResQ_%d", ihistoQ);
0756       std::string htitle = Form("BTL time resolution (Q bin = %d);T_{RECO} - T_{SIM} [ns]", ihistoQ);
0757       meTimeResQ_[ihistoQ] = ibook.book1D(hname, htitle, 200, -0.3, 0.7);
0758
0759       for (unsigned int ihistoEta = 0; ihistoEta < nBinsQEta_; ++ihistoEta) {
0760         hname = Form("TimeResQvsEta_%d_%d", ihistoQ, ihistoEta);
0761         htitle = Form("BTL time resolution (Q bin = %d, |#eta| bin = %d);T_{RECO} - T_{SIM} [ns]", ihistoQ, ihistoEta);
0762         meTimeResQvsEta_[ihistoQ][ihistoEta] = ibook.book1D(hname, htitle, 200, -0.3, 0.7);
0763
0764       }  // ihistoEta loop
0765
0766     }  // ihistoQ loop
0767
0768     for (unsigned int ihistoEta = 0; ihistoEta < nBinsEta_; ++ihistoEta) {
0769       std::string hname = Form("TimeResEta_%d", ihistoEta);
0770       std::string htitle = Form("BTL time resolution (|#eta| bin = %d);T_{RECO} - T_{SIM} [ns]", ihistoEta);
0771       meTimeResEta_[ihistoEta] = ibook.book1D(hname, htitle, 200, -0.3, 0.7);
0772
0773       for (unsigned int ihistoQ = 0; ihistoQ < nBinsEtaQ_; ++ihistoQ) {
0774         hname = Form("TimeResEtavsQ_%d_%d", ihistoEta, ihistoQ);
0775         htitle = Form("BTL time resolution (|#eta| bin = %d, Q bin = %d);T_{RECO} - T_{SIM} [ns]", ihistoEta, ihistoQ);
0776         meTimeResEtavsQ_[ihistoEta][ihistoQ] = ibook.book1D(hname, htitle, 200, -0.3, 0.7);
0777
0778       }  // ihistoQ loop
0779
0780     }  // ihistoEta loop
0781   }
0782 }
0783
0784 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0785 void BtlLocalRecoValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0786   edm::ParameterSetDescription desc;
0787
0788   desc.add<std::string>("folder", "MTD/BTL/LocalReco");
0789   desc.add<edm::InputTag>("recHitsTag", edm::InputTag("mtdRecHits", "FTLBarrel"));
0790   desc.add<edm::InputTag>("uncalibRecHitsTag", edm::InputTag("mtdUncalibratedRecHits", "FTLBarrel"));
0791   desc.add<edm::InputTag>("simHitsTag", edm::InputTag("mix", "g4SimHitsFastTimerHitsBarrel"));
0792   desc.add<edm::InputTag>("recCluTag", edm::InputTag("mtdClusters", "FTLBarrel"));
0793   desc.add<edm::InputTag>("trkHitTag", edm::InputTag("mtdTrackingRecHits"));
0794   desc.add<double>("HitMinimumEnergy", 1.);  // [MeV]
0795   desc.add<bool>("optionalPlots", false);
0796   desc.add<bool>("UncalibRecHitsPlots", false);
0797   desc.add<double>("HitMinimumAmplitude", 30.);  // [pC]
0798
0799   descriptions.add("btlLocalRecoValid", desc);
0800 }
0801
0802 DEFINE_FWK_MODULE(BtlLocalRecoValidation);