Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-09-20 02:50:44

0001 #include <string>
0002 
0003 #include "FWCore/Framework/interface/Frameworkfwd.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "DataFormats/Common/interface/ValueMap.h"
0008 
0009 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
0010 #include "DQMServices/Core/interface/DQMStore.h"
0011 
0012 #include "DataFormats/Common/interface/ValidHandle.h"
0013 #include "DataFormats/Math/interface/deltaR.h"
0014 #include "DataFormats/Math/interface/GeantUnits.h"
0015 #include "DataFormats/Math/interface/angle_units.h"
0016 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0017 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0018 
0019 #include "DataFormats/Common/interface/Ptr.h"
0020 #include "DataFormats/Common/interface/PtrVector.h"
0021 #include "DataFormats/Common/interface/RefProd.h"
0022 #include "DataFormats/Common/interface/Ref.h"
0023 #include "DataFormats/Common/interface/RefVector.h"
0024 
0025 #include "DataFormats/TrackReco/interface/Track.h"
0026 #include "DataFormats/VertexReco/interface/Vertex.h"
0027 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0028 
0029 #include "Geometry/Records/interface/MTDDigiGeometryRecord.h"
0030 #include "Geometry/Records/interface/MTDTopologyRcd.h"
0031 #include "Geometry/MTDNumberingBuilder/interface/MTDTopology.h"
0032 #include "Geometry/MTDCommonData/interface/MTDTopologyMode.h"
0033 #include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h"
0034 #include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h"
0035 #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
0036 
0037 #include "RecoMTD/DetLayers/interface/MTDDetLayerGeometry.h"
0038 #include "RecoMTD/Records/interface/MTDRecoGeometryRecord.h"
0039 #include "MagneticField/Engine/interface/MagneticField.h"
0040 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0041 #include "RecoMTD/DetLayers/interface/MTDDetLayerGeometry.h"
0042 #include "RecoMTD/DetLayers/interface/MTDTrayBarrelLayer.h"
0043 #include "RecoMTD/DetLayers/interface/MTDDetTray.h"
0044 #include "RecoMTD/DetLayers/interface/MTDSectorForwardDoubleLayer.h"
0045 #include "RecoMTD/DetLayers/interface/MTDDetSector.h"
0046 #include "RecoMTD/Records/interface/MTDRecoGeometryRecord.h"
0047 
0048 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0049 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0050 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0051 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0052 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0053 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0054 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0055 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h"
0056 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h"
0057 #include "DataFormats/FTLRecHit/interface/FTLRecHitCollections.h"
0058 
0059 #include "DataFormats/Common/interface/OneToMany.h"
0060 #include "DataFormats/Common/interface/AssociationMap.h"
0061 
0062 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0063 #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h"
0064 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
0065 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
0066 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
0067 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0068 
0069 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0070 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0071 #include "HepMC/GenRanges.h"
0072 #include "CLHEP/Units/PhysicalConstants.h"
0073 #include "MTDHit.h"
0074 
0075 class MtdTracksValidation : public DQMEDAnalyzer {
0076 public:
0077   explicit MtdTracksValidation(const edm::ParameterSet&);
0078   ~MtdTracksValidation() override;
0079 
0080   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0081 
0082 private:
0083   void bookHistograms(DQMStore::IBooker&, edm::Run const&, edm::EventSetup const&) override;
0084 
0085   void analyze(const edm::Event&, const edm::EventSetup&) override;
0086 
0087   const std::pair<bool, bool> checkAcceptance(
0088       const reco::Track&, const edm::Event&, const edm::EventSetup&, size_t&, float&, float&, float&, float&);
0089 
0090   const bool mvaGenSel(const HepMC::GenParticle&, const float&);
0091   const bool mvaTPSel(const TrackingParticle&);
0092   const bool mvaRecSel(const reco::TrackBase&, const reco::Vertex&, const double&, const double&);
0093   const bool mvaGenRecMatch(const HepMC::GenParticle&, const double&, const reco::TrackBase&, const bool&);
0094   const edm::Ref<std::vector<TrackingParticle>>* getMatchedTP(const reco::TrackBaseRef&);
0095 
0096   const unsigned long int uniqueId(const uint32_t x, const EncodedEventId& y) {
0097     const uint64_t a = static_cast<uint64_t>(x);
0098     const uint64_t b = static_cast<uint64_t>(y.rawId());
0099 
0100     if (x < y.rawId())
0101       return (b << 32) | a;
0102     else
0103       return (a << 32) | b;
0104   }
0105 
0106   bool isETL(const double eta) const { return (std::abs(eta) > trackMinEtlEta_) && (std::abs(eta) < trackMaxEtlEta_); }
0107 
0108   // ------------ member data ------------
0109 
0110   const std::string folder_;
0111   const float trackMinPt_;
0112   const float trackMaxBtlEta_;
0113   const float trackMinEtlEta_;
0114   const float trackMaxEtlEta_;
0115 
0116   static constexpr double etacutGEN_ = 4.;               // |eta| < 4;
0117   static constexpr double etacutREC_ = 3.;               // |eta| < 3;
0118   static constexpr double pTcut_ = 0.7;                  // PT > 0.7 GeV
0119   static constexpr double deltaZcut_ = 0.1;              // dz separation 1 mm
0120   static constexpr double deltaPTcut_ = 0.05;            // dPT < 5%
0121   static constexpr double deltaDRcut_ = 0.03;            // DeltaR separation
0122   static constexpr double depositBTLthreshold_ = 1;      // threshold for energy deposit in BTL cell [MeV]
0123   static constexpr double depositETLthreshold_ = 0.001;  // threshold for energy deposit in ETL cell [MeV]
0124   static constexpr double rBTL_ = 110.0;
0125   static constexpr double zETL_ = 290.0;
0126   static constexpr double etaMatchCut_ = 0.05;
0127   static constexpr double cluDRradius_ = 0.05;  // to cluster rechits around extrapolated track
0128 
0129   const reco::RecoToSimCollection* r2s_;
0130   const reco::SimToRecoCollection* s2r_;
0131 
0132   edm::EDGetTokenT<reco::TrackCollection> GenRecTrackToken_;
0133   edm::EDGetTokenT<reco::TrackCollection> RecTrackToken_;
0134   edm::EDGetTokenT<std::vector<reco::Vertex>> RecVertexToken_;
0135 
0136   edm::EDGetTokenT<edm::HepMCProduct> HepMCProductToken_;
0137   edm::EDGetTokenT<TrackingParticleCollection> trackingParticleCollectionToken_;
0138   edm::EDGetTokenT<reco::SimToRecoCollection> simToRecoAssociationToken_;
0139   edm::EDGetTokenT<reco::RecoToSimCollection> recoToSimAssociationToken_;
0140   edm::EDGetTokenT<CrossingFrame<PSimHit>> btlSimHitsToken_;
0141   edm::EDGetTokenT<CrossingFrame<PSimHit>> etlSimHitsToken_;
0142   edm::EDGetTokenT<FTLRecHitCollection> btlRecHitsToken_;
0143   edm::EDGetTokenT<FTLRecHitCollection> etlRecHitsToken_;
0144 
0145   edm::EDGetTokenT<edm::ValueMap<int>> trackAssocToken_;
0146   edm::EDGetTokenT<edm::ValueMap<float>> pathLengthToken_;
0147 
0148   edm::EDGetTokenT<edm::ValueMap<float>> tmtdToken_;
0149   edm::EDGetTokenT<edm::ValueMap<float>> SigmatmtdToken_;
0150   edm::EDGetTokenT<edm::ValueMap<float>> t0SrcToken_;
0151   edm::EDGetTokenT<edm::ValueMap<float>> Sigmat0SrcToken_;
0152   edm::EDGetTokenT<edm::ValueMap<float>> t0PidToken_;
0153   edm::EDGetTokenT<edm::ValueMap<float>> Sigmat0PidToken_;
0154   edm::EDGetTokenT<edm::ValueMap<float>> t0SafePidToken_;
0155   edm::EDGetTokenT<edm::ValueMap<float>> Sigmat0SafePidToken_;
0156   edm::EDGetTokenT<edm::ValueMap<float>> trackMVAQualToken_;
0157 
0158   edm::ESGetToken<MTDGeometry, MTDDigiGeometryRecord> mtdgeoToken_;
0159   edm::ESGetToken<MTDTopology, MTDTopologyRcd> mtdtopoToken_;
0160   edm::ESGetToken<MTDDetLayerGeometry, MTDRecoGeometryRecord> mtdlayerToken_;
0161   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magfieldToken_;
0162   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> builderToken_;
0163   edm::ESGetToken<HepPDT::ParticleDataTable, edm::DefaultRecord> particleTableToken_;
0164 
0165   MonitorElement* meBTLTrackRPTime_;
0166   MonitorElement* meBTLTrackEffEtaTot_;
0167   MonitorElement* meBTLTrackEffPhiTot_;
0168   MonitorElement* meBTLTrackEffPtTot_;
0169   MonitorElement* meBTLTrackEffEtaMtd_;
0170   MonitorElement* meBTLTrackEffPhiMtd_;
0171   MonitorElement* meBTLTrackEffPtMtd_;
0172   MonitorElement* meBTLTrackPtRes_;
0173 
0174   MonitorElement* meETLTrackRPTime_;
0175   MonitorElement* meETLTrackEffEtaTot_[2];
0176   MonitorElement* meETLTrackEffPhiTot_[2];
0177   MonitorElement* meETLTrackEffPtTot_[2];
0178   MonitorElement* meETLTrackEffEtaMtd_[2];
0179   MonitorElement* meETLTrackEffPhiMtd_[2];
0180   MonitorElement* meETLTrackEffPtMtd_[2];
0181   MonitorElement* meETLTrackEffEta2Mtd_[2];
0182   MonitorElement* meETLTrackEffPhi2Mtd_[2];
0183   MonitorElement* meETLTrackEffPt2Mtd_[2];
0184   MonitorElement* meETLTrackPtRes_;
0185 
0186   MonitorElement* meTracktmtd_;
0187   MonitorElement* meTrackt0Src_;
0188   MonitorElement* meTrackSigmat0Src_;
0189   MonitorElement* meTrackt0Pid_;
0190   MonitorElement* meTrackSigmat0Pid_;
0191   MonitorElement* meTrackt0SafePid_;
0192   MonitorElement* meTrackSigmat0SafePid_;
0193   MonitorElement* meTrackNumHits_;
0194   MonitorElement* meTrackNumHitsNT_;
0195   MonitorElement* meTrackMVAQual_;
0196   MonitorElement* meTrackPathLenghtvsEta_;
0197 
0198   MonitorElement* meTrackPtTot_;
0199   MonitorElement* meMVATrackEffPtTot_;
0200   MonitorElement* meMVATrackMatchedEffPtTot_;
0201   MonitorElement* meMVATrackMatchedEffPtMtd_;
0202   MonitorElement* meExtraPtMtd_;
0203   MonitorElement* meExtraPtEtl2Mtd_;
0204   MonitorElement* meTrackMatchedTPEffPtTot_;
0205   MonitorElement* meTrackMatchedTPEffPtMtd_;
0206   MonitorElement* meTrackMatchedTPEffPtEtl2Mtd_;
0207   MonitorElement* meTrackMatchedTPmtdEffPtTot_;
0208   MonitorElement* meTrackMatchedTPmtdEffPtMtd_;
0209   MonitorElement* meTrackEtaTot_;
0210   MonitorElement* meMVATrackEffEtaTot_;
0211   MonitorElement* meMVATrackMatchedEffEtaTot_;
0212   MonitorElement* meMVATrackMatchedEffEtaMtd_;
0213   MonitorElement* meExtraEtaMtd_;
0214   MonitorElement* meExtraEtaEtl2Mtd_;
0215   MonitorElement* meTrackMatchedTPEffEtaTot_;
0216   MonitorElement* meTrackMatchedTPEffEtaMtd_;
0217   MonitorElement* meTrackMatchedTPEffEtaEtl2Mtd_;
0218   MonitorElement* meTrackMatchedTPmtdEffEtaTot_;
0219   MonitorElement* meTrackMatchedTPmtdEffEtaMtd_;
0220   MonitorElement* meMVATrackResTot_;
0221   MonitorElement* meMVATrackPullTot_;
0222   MonitorElement* meMVATrackZposResTot_;
0223 
0224   MonitorElement* meExtraPhiAtBTL_;
0225   MonitorElement* meExtraPhiAtBTLmatched_;
0226   MonitorElement* meExtraBTLeneInCone_;
0227   MonitorElement* meExtraMTDfailExtenderEta_;
0228   MonitorElement* meExtraMTDfailExtenderPt_;
0229 };
0230 
0231 // ------------ constructor and destructor --------------
0232 MtdTracksValidation::MtdTracksValidation(const edm::ParameterSet& iConfig)
0233     : folder_(iConfig.getParameter<std::string>("folder")),
0234       trackMinPt_(iConfig.getParameter<double>("trackMinimumPt")),
0235       trackMaxBtlEta_(iConfig.getParameter<double>("trackMaximumBtlEta")),
0236       trackMinEtlEta_(iConfig.getParameter<double>("trackMinimumEtlEta")),
0237       trackMaxEtlEta_(iConfig.getParameter<double>("trackMaximumEtlEta")) {
0238   GenRecTrackToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("inputTagG"));
0239   RecTrackToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("inputTagT"));
0240   RecVertexToken_ = consumes<std::vector<reco::Vertex>>(iConfig.getParameter<edm::InputTag>("inputTagV"));
0241   HepMCProductToken_ = consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("inputTagH"));
0242   trackingParticleCollectionToken_ =
0243       consumes<TrackingParticleCollection>(iConfig.getParameter<edm::InputTag>("SimTag"));
0244   simToRecoAssociationToken_ =
0245       consumes<reco::SimToRecoCollection>(iConfig.getParameter<edm::InputTag>("TPtoRecoTrackAssoc"));
0246   recoToSimAssociationToken_ =
0247       consumes<reco::RecoToSimCollection>(iConfig.getParameter<edm::InputTag>("TPtoRecoTrackAssoc"));
0248   btlSimHitsToken_ = consumes<CrossingFrame<PSimHit>>(iConfig.getParameter<edm::InputTag>("btlSimHits"));
0249   etlSimHitsToken_ = consumes<CrossingFrame<PSimHit>>(iConfig.getParameter<edm::InputTag>("etlSimHits"));
0250   btlRecHitsToken_ = consumes<FTLRecHitCollection>(iConfig.getParameter<edm::InputTag>("btlRecHits"));
0251   etlRecHitsToken_ = consumes<FTLRecHitCollection>(iConfig.getParameter<edm::InputTag>("etlRecHits"));
0252   trackAssocToken_ = consumes<edm::ValueMap<int>>(iConfig.getParameter<edm::InputTag>("trackAssocSrc"));
0253   pathLengthToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("pathLengthSrc"));
0254   tmtdToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tmtd"));
0255   SigmatmtdToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmatmtd"));
0256   t0SrcToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0Src"));
0257   Sigmat0SrcToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmat0Src"));
0258   t0PidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0PID"));
0259   Sigmat0PidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmat0PID"));
0260   t0SafePidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("t0SafePID"));
0261   Sigmat0SafePidToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("sigmat0SafePID"));
0262   trackMVAQualToken_ = consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("trackMVAQual"));
0263   mtdgeoToken_ = esConsumes<MTDGeometry, MTDDigiGeometryRecord>();
0264   mtdtopoToken_ = esConsumes<MTDTopology, MTDTopologyRcd>();
0265   mtdlayerToken_ = esConsumes<MTDDetLayerGeometry, MTDRecoGeometryRecord>();
0266   magfieldToken_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0267   builderToken_ = esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"));
0268   particleTableToken_ = esConsumes<HepPDT::ParticleDataTable, edm::DefaultRecord>();
0269 }
0270 
0271 MtdTracksValidation::~MtdTracksValidation() {}
0272 
0273 // ------------ method called for each event  ------------
0274 void MtdTracksValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0275   using namespace edm;
0276   using namespace geant_units::operators;
0277   using namespace std;
0278 
0279   auto GenRecTrackHandle = makeValid(iEvent.getHandle(GenRecTrackToken_));
0280   auto RecVertexHandle = makeValid(iEvent.getHandle(RecVertexToken_));
0281 
0282   std::unordered_map<uint32_t, MTDHit> m_btlHits;
0283   std::unordered_map<uint32_t, MTDHit> m_etlHits;
0284   std::unordered_map<uint32_t, std::set<unsigned long int>> m_btlTrkPerCell;
0285   std::unordered_map<uint32_t, std::set<unsigned long int>> m_etlTrkPerCell;
0286   std::map<TrackingParticleRef, std::vector<uint32_t>> m_tp2detid;
0287 
0288   const auto& tMtd = iEvent.get(tmtdToken_);
0289   const auto& SigmatMtd = iEvent.get(SigmatmtdToken_);
0290   const auto& t0Src = iEvent.get(t0SrcToken_);
0291   const auto& Sigmat0Src = iEvent.get(Sigmat0SrcToken_);
0292   const auto& t0Pid = iEvent.get(t0PidToken_);
0293   const auto& Sigmat0Pid = iEvent.get(Sigmat0PidToken_);
0294   const auto& t0Safe = iEvent.get(t0SafePidToken_);
0295   const auto& Sigmat0Safe = iEvent.get(Sigmat0SafePidToken_);
0296   const auto& mtdQualMVA = iEvent.get(trackMVAQualToken_);
0297   const auto& trackAssoc = iEvent.get(trackAssocToken_);
0298   const auto& pathLength = iEvent.get(pathLengthToken_);
0299 
0300   const auto& primRecoVtx = *(RecVertexHandle.product()->begin());
0301 
0302   // generator level information (HepMC format)
0303   auto GenEventHandle = makeValid(iEvent.getHandle(HepMCProductToken_));
0304   const HepMC::GenEvent* mc = GenEventHandle->GetEvent();
0305   double zsim = convertMmToCm((*(mc->vertices_begin()))->position().z());
0306   double tsim = (*(mc->vertices_begin()))->position().t() * CLHEP::mm / CLHEP::c_light;
0307 
0308   auto pdt = iSetup.getHandle(particleTableToken_);
0309   const HepPDT::ParticleDataTable* pdTable = pdt.product();
0310 
0311   auto simToRecoH = makeValid(iEvent.getHandle(simToRecoAssociationToken_));
0312   s2r_ = simToRecoH.product();
0313 
0314   auto recoToSimH = makeValid(iEvent.getHandle(recoToSimAssociationToken_));
0315   r2s_ = recoToSimH.product();
0316 
0317   //Fill maps with simhits accumulated per DetId
0318 
0319   auto btlSimHitsHandle = makeValid(iEvent.getHandle(btlSimHitsToken_));
0320   MixCollection<PSimHit> btlSimHits(btlSimHitsHandle.product());
0321   for (auto const& simHit : btlSimHits) {
0322     if (simHit.tof() < 0 || simHit.tof() > 25.)
0323       continue;
0324     DetId id = simHit.detUnitId();
0325     auto const thisHId = uniqueId(simHit.trackId(), simHit.eventId());
0326     m_btlTrkPerCell[id.rawId()].insert(thisHId);
0327     auto simHitIt = m_btlHits.emplace(id.rawId(), MTDHit()).first;
0328     // --- Accumulate the energy (in MeV) of SIM hits in the same detector cell
0329     (simHitIt->second).energy += convertUnitsTo(0.001_MeV, simHit.energyLoss());
0330   }
0331 
0332   auto etlSimHitsHandle = makeValid(iEvent.getHandle(etlSimHitsToken_));
0333   MixCollection<PSimHit> etlSimHits(etlSimHitsHandle.product());
0334   for (auto const& simHit : etlSimHits) {
0335     if (simHit.tof() < 0 || simHit.tof() > 25.) {
0336       continue;
0337     }
0338     DetId id = simHit.detUnitId();
0339     auto const thisHId = uniqueId(simHit.trackId(), simHit.eventId());
0340     m_etlTrkPerCell[id.rawId()].insert(thisHId);
0341     auto simHitIt = m_etlHits.emplace(id.rawId(), MTDHit()).first;
0342     // --- Accumulate the energy (in MeV) of SIM hits in the same detector cell
0343     (simHitIt->second).energy += convertUnitsTo(0.001_MeV, simHit.energyLoss());
0344   }
0345 
0346   //Fill map of DetId per ref to TP
0347 
0348   auto tpHandle = makeValid(iEvent.getHandle(trackingParticleCollectionToken_));
0349   TrackingParticleCollection tpColl = *(tpHandle.product());
0350   size_t tpindex(0);
0351   for (auto tp = tpColl.begin(); tp != tpColl.end(); tp++, ++tpindex) {
0352     TrackingParticleRef tpref(iEvent.getHandle(trackingParticleCollectionToken_), tpindex);
0353     if (tp->eventId().bunchCrossing() == 0 && tp->eventId().event() == 0) {
0354       if (!mvaTPSel(*tp))
0355         continue;
0356       for (const auto& simTrk : tp->g4Tracks()) {
0357         auto const thisTId = uniqueId(simTrk.trackId(), simTrk.eventId());
0358         for (auto const& cell : m_btlTrkPerCell) {
0359           if (m_btlHits[cell.first].energy < depositBTLthreshold_) {
0360             continue;
0361           }
0362           for (auto const& simtrack : cell.second) {
0363             if (thisTId == simtrack) {
0364               m_tp2detid[tpref].emplace_back(cell.first);
0365               break;
0366             }
0367           }
0368         }
0369         for (auto const& cell : m_etlTrkPerCell) {
0370           if (m_etlHits[cell.first].energy < depositETLthreshold_) {
0371             continue;
0372           }
0373           for (auto const& simtrack : cell.second) {
0374             if (thisTId == simtrack) {
0375               m_tp2detid[tpref].emplace_back(cell.first);
0376               break;
0377             }
0378           }
0379         }
0380       }
0381     }
0382   }
0383 
0384   unsigned int index = 0;
0385 
0386   // flag to select events with reco vertex close to true simulated primary vertex, or PV fake (particle guns)
0387   const bool isGoodVtx = std::abs(primRecoVtx.z() - zsim) < deltaZcut_ || primRecoVtx.isFake();
0388 
0389   // --- Loop over all RECO tracks ---
0390   for (const auto& trackGen : *GenRecTrackHandle) {
0391     const reco::TrackRef trackref(iEvent.getHandle(GenRecTrackToken_), index);
0392     index++;
0393 
0394     if (trackAssoc[trackref] == -1) {
0395       LogInfo("mtdTracks") << "Extended track not associated";
0396       continue;
0397     }
0398 
0399     const reco::TrackRef mtdTrackref = reco::TrackRef(iEvent.getHandle(RecTrackToken_), trackAssoc[trackref]);
0400     const reco::Track& track = *mtdTrackref;
0401 
0402     bool isBTL = false;
0403     bool isETL = false;
0404     bool twoETLdiscs = false;
0405     bool noCrack = std::abs(trackGen.eta()) < trackMaxBtlEta_ || std::abs(trackGen.eta()) > trackMinEtlEta_;
0406 
0407     if (track.pt() >= trackMinPt_ && std::abs(track.eta()) <= trackMaxEtlEta_) {
0408       meTracktmtd_->Fill(tMtd[trackref]);
0409       if (std::round(SigmatMtd[trackref] - Sigmat0Pid[trackref]) != 0) {
0410         LogWarning("mtdTracks")
0411             << "TimeError associated to refitted track is different from TimeError stored in tofPID "
0412                "sigmat0 ValueMap: this should not happen";
0413       }
0414 
0415       meTrackt0Src_->Fill(t0Src[trackref]);
0416       meTrackSigmat0Src_->Fill(Sigmat0Src[trackref]);
0417 
0418       meTrackt0Pid_->Fill(t0Pid[trackref]);
0419       meTrackSigmat0Pid_->Fill(Sigmat0Pid[trackref]);
0420       meTrackt0SafePid_->Fill(t0Safe[trackref]);
0421       meTrackSigmat0SafePid_->Fill(Sigmat0Safe[trackref]);
0422       meTrackMVAQual_->Fill(mtdQualMVA[trackref]);
0423 
0424       meTrackPathLenghtvsEta_->Fill(std::abs(track.eta()), pathLength[trackref]);
0425 
0426       if (std::abs(track.eta()) < trackMaxBtlEta_) {
0427         // --- all BTL tracks (with and without hit in MTD) ---
0428         meBTLTrackEffEtaTot_->Fill(track.eta());
0429         meBTLTrackEffPhiTot_->Fill(track.phi());
0430         meBTLTrackEffPtTot_->Fill(track.pt());
0431 
0432         bool MTDBtl = false;
0433         int numMTDBtlvalidhits = 0;
0434         for (const auto hit : track.recHits()) {
0435           if (hit->isValid() == false)
0436             continue;
0437           MTDDetId Hit = hit->geographicalId();
0438           if ((Hit.det() == 6) && (Hit.subdetId() == 1) && (Hit.mtdSubDetector() == 1)) {
0439             MTDBtl = true;
0440             numMTDBtlvalidhits++;
0441           }
0442         }
0443         meTrackNumHits_->Fill(numMTDBtlvalidhits);
0444 
0445         // --- keeping only tracks with last hit in MTD ---
0446         if (MTDBtl == true) {
0447           isBTL = true;
0448           meBTLTrackEffEtaMtd_->Fill(track.eta());
0449           meBTLTrackEffPhiMtd_->Fill(track.phi());
0450           meBTLTrackEffPtMtd_->Fill(track.pt());
0451           meBTLTrackRPTime_->Fill(track.t0());
0452           meBTLTrackPtRes_->Fill((trackGen.pt() - track.pt()) / trackGen.pt());
0453         }
0454         if (isBTL && Sigmat0Safe[trackref] < 0.) {
0455           meTrackNumHitsNT_->Fill(numMTDBtlvalidhits);
0456         }
0457       }  //loop over (geometrical) BTL tracks
0458 
0459       else {
0460         // --- all ETL tracks (with and without hit in MTD) ---
0461         if ((track.eta() < -trackMinEtlEta_) && (track.eta() > -trackMaxEtlEta_)) {
0462           meETLTrackEffEtaTot_[0]->Fill(track.eta());
0463           meETLTrackEffPhiTot_[0]->Fill(track.phi());
0464           meETLTrackEffPtTot_[0]->Fill(track.pt());
0465         }
0466 
0467         if ((track.eta() > trackMinEtlEta_) && (track.eta() < trackMaxEtlEta_)) {
0468           meETLTrackEffEtaTot_[1]->Fill(track.eta());
0469           meETLTrackEffPhiTot_[1]->Fill(track.phi());
0470           meETLTrackEffPtTot_[1]->Fill(track.pt());
0471         }
0472 
0473         bool MTDEtlZnegD1 = false;
0474         bool MTDEtlZnegD2 = false;
0475         bool MTDEtlZposD1 = false;
0476         bool MTDEtlZposD2 = false;
0477         int numMTDEtlvalidhits = 0;
0478         for (const auto hit : track.recHits()) {
0479           if (hit->isValid() == false)
0480             continue;
0481           MTDDetId Hit = hit->geographicalId();
0482           if ((Hit.det() == 6) && (Hit.subdetId() == 1) && (Hit.mtdSubDetector() == 2)) {
0483             isETL = true;
0484             ETLDetId ETLHit = hit->geographicalId();
0485 
0486             if ((ETLHit.zside() == -1) && (ETLHit.nDisc() == 1)) {
0487               MTDEtlZnegD1 = true;
0488               meETLTrackRPTime_->Fill(track.t0());
0489               meETLTrackPtRes_->Fill((trackGen.pt() - track.pt()) / trackGen.pt());
0490               numMTDEtlvalidhits++;
0491             }
0492             if ((ETLHit.zside() == -1) && (ETLHit.nDisc() == 2)) {
0493               MTDEtlZnegD2 = true;
0494               meETLTrackRPTime_->Fill(track.t0());
0495               meETLTrackPtRes_->Fill((trackGen.pt() - track.pt()) / trackGen.pt());
0496               numMTDEtlvalidhits++;
0497             }
0498             if ((ETLHit.zside() == 1) && (ETLHit.nDisc() == 1)) {
0499               MTDEtlZposD1 = true;
0500               meETLTrackRPTime_->Fill(track.t0());
0501               meETLTrackPtRes_->Fill((trackGen.pt() - track.pt()) / trackGen.pt());
0502               numMTDEtlvalidhits++;
0503             }
0504             if ((ETLHit.zside() == 1) && (ETLHit.nDisc() == 2)) {
0505               MTDEtlZposD2 = true;
0506               meETLTrackRPTime_->Fill(track.t0());
0507               meETLTrackPtRes_->Fill((trackGen.pt() - track.pt()) / trackGen.pt());
0508               numMTDEtlvalidhits++;
0509             }
0510           }
0511         }
0512         meTrackNumHits_->Fill(-numMTDEtlvalidhits);
0513         if (isETL && Sigmat0Safe[trackref] < 0.) {
0514           meTrackNumHitsNT_->Fill(-numMTDEtlvalidhits);
0515         }
0516 
0517         // --- keeping only tracks with last hit in MTD ---
0518         if ((track.eta() < -trackMinEtlEta_) && (track.eta() > -trackMaxEtlEta_)) {
0519           twoETLdiscs = (MTDEtlZnegD1 == true) && (MTDEtlZnegD2 == true);
0520           if ((MTDEtlZnegD1 == true) || (MTDEtlZnegD2 == true)) {
0521             meETLTrackEffEtaMtd_[0]->Fill(track.eta());
0522             meETLTrackEffPhiMtd_[0]->Fill(track.phi());
0523             meETLTrackEffPtMtd_[0]->Fill(track.pt());
0524             if (twoETLdiscs) {
0525               meETLTrackEffEta2Mtd_[0]->Fill(track.eta());
0526               meETLTrackEffPhi2Mtd_[0]->Fill(track.phi());
0527               meETLTrackEffPt2Mtd_[0]->Fill(track.pt());
0528             }
0529           }
0530         }
0531         if ((track.eta() > trackMinEtlEta_) && (track.eta() < trackMaxEtlEta_)) {
0532           twoETLdiscs = (MTDEtlZposD1 == true) && (MTDEtlZposD2 == true);
0533           if ((MTDEtlZposD1 == true) || (MTDEtlZposD2 == true)) {
0534             meETLTrackEffEtaMtd_[1]->Fill(track.eta());
0535             meETLTrackEffPhiMtd_[1]->Fill(track.phi());
0536             meETLTrackEffPtMtd_[1]->Fill(track.pt());
0537             if (twoETLdiscs) {
0538               meETLTrackEffEta2Mtd_[1]->Fill(track.eta());
0539               meETLTrackEffPhi2Mtd_[1]->Fill(track.phi());
0540               meETLTrackEffPt2Mtd_[1]->Fill(track.pt());
0541             }
0542           }
0543         }
0544       }
0545 
0546       LogDebug("MtdTracksValidation") << "Track p/pt = " << track.p() << " " << track.pt() << " eta " << track.eta()
0547                                       << " BTL " << isBTL << " ETL " << isETL << " 2disks " << twoETLdiscs;
0548 
0549       // TrackingParticle based matching
0550 
0551       const reco::TrackBaseRef trkrefb(trackref);
0552       auto tp_info = getMatchedTP(trkrefb);
0553 
0554       meTrackPtTot_->Fill(trackGen.pt());
0555       meTrackEtaTot_->Fill(std::abs(trackGen.eta()));
0556       if (tp_info != nullptr && mvaTPSel(**tp_info)) {
0557         const bool withMTD = (m_tp2detid.find(*tp_info) != m_tp2detid.end());
0558         LogDebug("MtdTracksValidation") << "Matched with selected TP, MTD sim hits association: " << withMTD;
0559         if (noCrack) {
0560           meTrackMatchedTPEffPtTot_->Fill(trackGen.pt());
0561           if (withMTD) {
0562             meTrackMatchedTPmtdEffPtTot_->Fill(trackGen.pt());
0563           }
0564         }
0565         meTrackMatchedTPEffEtaTot_->Fill(std::abs(trackGen.eta()));
0566         if (withMTD) {
0567           meTrackMatchedTPmtdEffEtaTot_->Fill(std::abs(trackGen.eta()));
0568         }
0569         if (isBTL || isETL) {
0570           if (noCrack) {
0571             meTrackMatchedTPEffPtMtd_->Fill(trackGen.pt());
0572             if (isBTL || twoETLdiscs) {
0573               meTrackMatchedTPEffPtEtl2Mtd_->Fill(trackGen.pt());
0574             }
0575             if (withMTD) {
0576               meTrackMatchedTPmtdEffPtMtd_->Fill(trackGen.pt());
0577             }
0578           }
0579           meTrackMatchedTPEffEtaMtd_->Fill(std::abs(trackGen.eta()));
0580           if (isBTL || twoETLdiscs) {
0581             meTrackMatchedTPEffEtaEtl2Mtd_->Fill(std::abs(trackGen.eta()));
0582           }
0583           if (withMTD) {
0584             meTrackMatchedTPmtdEffEtaMtd_->Fill(std::abs(trackGen.eta()));
0585           }
0586         }
0587 
0588         size_t nlayers(0);
0589         float extrho(0.);
0590         float exteta(0.);
0591         float extphi(0.);
0592         float selvar(0.);
0593         auto accept = checkAcceptance(trackGen, iEvent, iSetup, nlayers, extrho, exteta, extphi, selvar);
0594         if (accept.first && std::abs(exteta) < trackMaxBtlEta_) {
0595           meExtraPhiAtBTL_->Fill(angle_units::operators::convertRadToDeg(extphi));
0596           meExtraBTLeneInCone_->Fill(selvar);
0597         }
0598         if (accept.second) {
0599           if (std::abs(exteta) < trackMaxBtlEta_) {
0600             meExtraPhiAtBTLmatched_->Fill(angle_units::operators::convertRadToDeg(extphi));
0601           }
0602           if (noCrack) {
0603             meExtraPtMtd_->Fill(trackGen.pt());
0604             if (nlayers == 2) {
0605               meExtraPtEtl2Mtd_->Fill(trackGen.pt());
0606             }
0607           }
0608           meExtraEtaMtd_->Fill(std::abs(trackGen.eta()));
0609           if (nlayers == 2) {
0610             meExtraEtaEtl2Mtd_->Fill(std::abs(trackGen.eta()));
0611           }
0612           if (accept.first && accept.second && !(isBTL || isETL)) {
0613             edm::LogInfo("MtdTracksValidation")
0614                 << "MtdTracksValidation: extender fail in " << iEvent.id().run() << " " << iEvent.id().event()
0615                 << " pt= " << trackGen.pt() << " eta= " << trackGen.eta();
0616             meExtraMTDfailExtenderEta_->Fill(std::abs(trackGen.eta()));
0617             if (noCrack) {
0618               meExtraMTDfailExtenderPt_->Fill(trackGen.pt());
0619             }
0620           }
0621         }
0622 
0623       }  // TP matching
0624     }
0625 
0626     if (isGoodVtx) {
0627       const bool vtxFake = primRecoVtx.isFake();
0628 
0629       if (mvaRecSel(trackGen, primRecoVtx, t0Safe[trackref], Sigmat0Safe[trackref])) {
0630         // reco-gen matching used for MVA quality flag
0631 
0632         if (noCrack) {
0633           meMVATrackEffPtTot_->Fill(trackGen.pt());
0634         }
0635         meMVATrackEffEtaTot_->Fill(std::abs(trackGen.eta()));
0636 
0637         double dZ = trackGen.vz() - zsim;
0638         double dT(-9999.);
0639         double pullT(-9999.);
0640         if (Sigmat0Safe[trackref] != -1.) {
0641           dT = t0Safe[trackref] - tsim;
0642           pullT = dT / Sigmat0Safe[trackref];
0643         }
0644         for (const auto& genP : mc->particle_range()) {
0645           // select status 1 genParticles and match them to the reconstructed track
0646 
0647           float charge = pdTable->particle(HepPDT::ParticleID(genP->pdg_id())) != nullptr
0648                              ? pdTable->particle(HepPDT::ParticleID(genP->pdg_id()))->charge()
0649                              : 0.f;
0650           if (mvaGenSel(*genP, charge)) {
0651             if (mvaGenRecMatch(*genP, zsim, trackGen, vtxFake)) {
0652               meMVATrackZposResTot_->Fill(dZ);
0653               if (noCrack) {
0654                 meMVATrackMatchedEffPtTot_->Fill(trackGen.pt());
0655               }
0656               meMVATrackMatchedEffEtaTot_->Fill(std::abs(trackGen.eta()));
0657               if (isBTL || isETL) {
0658                 meMVATrackResTot_->Fill(dT);
0659                 meMVATrackPullTot_->Fill(pullT);
0660                 if (noCrack) {
0661                   meMVATrackMatchedEffPtMtd_->Fill(trackGen.pt());
0662                 }
0663                 meMVATrackMatchedEffEtaMtd_->Fill(std::abs(trackGen.eta()));
0664               }
0665               break;
0666             }
0667           }
0668         }
0669       }
0670     }  // MC truth matich analysis for good PV
0671   }    //RECO tracks loop
0672 }
0673 
0674 const std::pair<bool, bool> MtdTracksValidation::checkAcceptance(const reco::Track& track,
0675                                                                  const edm::Event& iEvent,
0676                                                                  edm::EventSetup const& iSetup,
0677                                                                  size_t& nlayers,
0678                                                                  float& extrho,
0679                                                                  float& exteta,
0680                                                                  float& extphi,
0681                                                                  float& selvar) {
0682   bool isMatched(false);
0683   nlayers = 0;
0684   extrho = 0.;
0685   exteta = -999.;
0686   extphi = -999.;
0687   selvar = 0.;
0688 
0689   auto geometryHandle = iSetup.getTransientHandle(mtdgeoToken_);
0690   const MTDGeometry* geom = geometryHandle.product();
0691   auto topologyHandle = iSetup.getTransientHandle(mtdtopoToken_);
0692   const MTDTopology* topology = topologyHandle.product();
0693 
0694   auto layerHandle = iSetup.getTransientHandle(mtdlayerToken_);
0695   const MTDDetLayerGeometry* layerGeo = layerHandle.product();
0696 
0697   auto magfieldHandle = iSetup.getTransientHandle(magfieldToken_);
0698   const MagneticField* mfield = magfieldHandle.product();
0699 
0700   auto ttrackBuilder = iSetup.getTransientHandle(builderToken_);
0701 
0702   auto tTrack = ttrackBuilder->build(track);
0703   TrajectoryStateOnSurface tsos = tTrack.outermostMeasurementState();
0704   float theMaxChi2 = 500.;
0705   float theNSigma = 10.;
0706   std::unique_ptr<MeasurementEstimator> theEstimator =
0707       std::make_unique<Chi2MeasurementEstimator>(theMaxChi2, theNSigma);
0708   SteppingHelixPropagator prop(mfield, anyDirection);
0709 
0710   auto btlRecHitsHandle = makeValid(iEvent.getHandle(btlRecHitsToken_));
0711   auto etlRecHitsHandle = makeValid(iEvent.getHandle(etlRecHitsToken_));
0712 
0713   edm::LogVerbatim("MtdTracksValidation")
0714       << "MtdTracksValidation: extrapolating track, pt= " << track.pt() << " eta= " << track.eta();
0715 
0716   //try BTL
0717   bool inBTL = false;
0718   float eneSum(0.);
0719   const std::vector<const DetLayer*>& layersBTL = layerGeo->allBTLLayers();
0720   for (const DetLayer* ilay : layersBTL) {
0721     std::pair<bool, TrajectoryStateOnSurface> comp = ilay->compatible(tsos, prop, *theEstimator);
0722     if (!comp.first)
0723       continue;
0724     if (!inBTL) {
0725       inBTL = true;
0726       extrho = comp.second.globalPosition().perp();
0727       exteta = comp.second.globalPosition().eta();
0728       extphi = comp.second.globalPosition().phi();
0729       edm::LogVerbatim("MtdTracksValidation") << "MtdTracksValidation: extrapolation at BTL surface, rho= " << extrho
0730                                               << " eta= " << exteta << " phi= " << extphi;
0731     }
0732     std::vector<DetLayer::DetWithState> compDets = ilay->compatibleDets(tsos, prop, *theEstimator);
0733     for (const auto& detWithState : compDets) {
0734       const auto& det = detWithState.first;
0735 
0736       // loop on compatible rechits and check energy in a fixed size cone around the extrapolation point
0737 
0738       edm::LogVerbatim("MtdTracksValidation")
0739           << "MtdTracksValidation: DetId= " << det->geographicalId().rawId()
0740           << " gp= " << detWithState.second.globalPosition().x() << " " << detWithState.second.globalPosition().y()
0741           << " " << detWithState.second.globalPosition().z() << " rho= " << detWithState.second.globalPosition().perp()
0742           << " eta= " << detWithState.second.globalPosition().eta()
0743           << " phi= " << detWithState.second.globalPosition().phi();
0744 
0745       for (const auto& recHit : *btlRecHitsHandle) {
0746         BTLDetId detId = recHit.id();
0747         DetId geoId = detId.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topology->getMTDTopologyMode()));
0748         const MTDGeomDet* thedet = geom->idToDet(geoId);
0749         if (thedet == nullptr)
0750           throw cms::Exception("MtdTracksValidation") << "GeographicalID: " << std::hex << geoId.rawId() << " ("
0751                                                       << detId.rawId() << ") is invalid!" << std::dec << std::endl;
0752         if (geoId == det->geographicalId()) {
0753           const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
0754           const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0755 
0756           Local3DPoint local_point(0., 0., 0.);
0757           local_point = topo.pixelToModuleLocalPoint(local_point, detId.row(topo.nrows()), detId.column(topo.nrows()));
0758           const auto& global_point = thedet->toGlobal(local_point);
0759           edm::LogVerbatim("MtdTracksValidation")
0760               << "MtdTracksValidation: Hit id= " << detId.rawId() << " ene= " << recHit.energy()
0761               << " dr= " << reco::deltaR(global_point, detWithState.second.globalPosition());
0762           if (reco::deltaR(global_point, detWithState.second.globalPosition()) < cluDRradius_) {
0763             eneSum += recHit.energy();
0764             //extrho = detWithState.second.globalPosition().perp();
0765             //exteta = detWithState.second.globalPosition().eta();
0766             //extphi = detWithState.second.globalPosition().phi();
0767           }
0768         }
0769       }
0770     }
0771     if (eneSum > depositBTLthreshold_) {
0772       nlayers++;
0773       selvar = eneSum;
0774       isMatched = true;
0775       edm::LogVerbatim("MtdTracksValidation")
0776           << "MtdTracksValidation: BTL matched, energy= " << eneSum << " #layers= " << nlayers;
0777     }
0778   }
0779   if (inBTL) {
0780     return std::make_pair(inBTL, isMatched);
0781   }
0782 
0783   //try ETL
0784   bool inETL = false;
0785   const std::vector<const DetLayer*>& layersETL = layerGeo->allETLLayers();
0786   for (const DetLayer* ilay : layersETL) {
0787     size_t hcount(0);
0788     const BoundDisk& disk = static_cast<const MTDSectorForwardDoubleLayer*>(ilay)->specificSurface();
0789     const double diskZ = disk.position().z();
0790     if (tsos.globalPosition().z() * diskZ < 0)
0791       continue;  // only propagate to the disk that's on the same side
0792     std::pair<bool, TrajectoryStateOnSurface> comp = ilay->compatible(tsos, prop, *theEstimator);
0793     if (!comp.first)
0794       continue;
0795     if (!inETL) {
0796       inETL = true;
0797       extrho = comp.second.globalPosition().perp();
0798       exteta = comp.second.globalPosition().eta();
0799       extphi = comp.second.globalPosition().phi();
0800     }
0801     edm::LogVerbatim("MtdTracksValidation") << "MtdTracksValidation: extrapolation at ETL surface, rho= " << extrho
0802                                             << " eta= " << exteta << " phi= " << extphi;
0803     std::vector<DetLayer::DetWithState> compDets = ilay->compatibleDets(tsos, prop, *theEstimator);
0804     for (const auto& detWithState : compDets) {
0805       const auto& det = detWithState.first;
0806 
0807       // loop on compatible rechits and check hits in a fixed size cone around the extrapolation point
0808 
0809       edm::LogVerbatim("MtdTracksValidation")
0810           << "MtdTracksValidation: DetId= " << det->geographicalId().rawId()
0811           << " gp= " << detWithState.second.globalPosition().x() << " " << detWithState.second.globalPosition().y()
0812           << " " << detWithState.second.globalPosition().z() << " rho= " << detWithState.second.globalPosition().perp()
0813           << " eta= " << detWithState.second.globalPosition().eta()
0814           << " phi= " << detWithState.second.globalPosition().phi();
0815 
0816       for (const auto& recHit : *etlRecHitsHandle) {
0817         ETLDetId detId = recHit.id();
0818         DetId geoId = detId.geographicalId();
0819         const MTDGeomDet* thedet = geom->idToDet(geoId);
0820         if (thedet == nullptr)
0821           throw cms::Exception("MtdTracksValidation") << "GeographicalID: " << std::hex << geoId.rawId() << " ("
0822                                                       << detId.rawId() << ") is invalid!" << std::dec << std::endl;
0823         if (geoId == det->geographicalId()) {
0824           const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
0825           const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0826 
0827           Local3DPoint local_point(topo.localX(recHit.row()), topo.localY(recHit.column()), 0.);
0828           const auto& global_point = thedet->toGlobal(local_point);
0829           edm::LogVerbatim("MtdTracksValidation")
0830               << "MtdTracksValidation: Hit id= " << detId.rawId() << " time= " << recHit.time()
0831               << " dr= " << reco::deltaR(global_point, detWithState.second.globalPosition());
0832           if (reco::deltaR(global_point, detWithState.second.globalPosition()) < cluDRradius_) {
0833             hcount++;
0834             if (hcount == 1) {
0835               //extrho = detWithState.second.globalPosition().perp();
0836               //exteta = detWithState.second.globalPosition().eta();
0837               //extphi = detWithState.second.globalPosition().phi();
0838             }
0839           }
0840         }
0841       }
0842     }
0843     if (hcount > 0) {
0844       nlayers++;
0845       selvar = (float)hcount;
0846       isMatched = true;
0847       edm::LogVerbatim("MtdTracksValidation")
0848           << "MtdTracksValidation: ETL matched, counts= " << hcount << " #layers= " << nlayers;
0849     }
0850   }
0851 
0852   if (!inBTL && !inETL) {
0853     edm::LogVerbatim("MtdTracksValidation")
0854         << "MtdTracksValidation: track not extrapolating to MTD: pt= " << track.pt() << " eta= " << track.eta()
0855         << " phi= " << track.phi() << " vz= " << track.vz()
0856         << " vxy= " << std::sqrt(track.vx() * track.vx() + track.vy() * track.vy());
0857   }
0858   return std::make_pair(inETL, isMatched);
0859 }
0860 
0861 // ------------ method for histogram booking ------------
0862 void MtdTracksValidation::bookHistograms(DQMStore::IBooker& ibook, edm::Run const& run, edm::EventSetup const& iSetup) {
0863   ibook.setCurrentFolder(folder_);
0864 
0865   // histogram booking
0866   meBTLTrackRPTime_ = ibook.book1D("TrackBTLRPTime", "Track t0 with respect to R.P.;t0 [ns]", 100, -1, 3);
0867   meBTLTrackEffEtaTot_ = ibook.book1D("TrackBTLEffEtaTot", "Track efficiency vs eta (Tot);#eta_{RECO}", 100, -1.6, 1.6);
0868   meBTLTrackEffPhiTot_ =
0869       ibook.book1D("TrackBTLEffPhiTot", "Track efficiency vs phi (Tot);#phi_{RECO} [rad]", 100, -3.2, 3.2);
0870   meBTLTrackEffPtTot_ = ibook.book1D("TrackBTLEffPtTot", "Track efficiency vs pt (Tot);pt_{RECO} [GeV]", 50, 0, 10);
0871   meBTLTrackEffEtaMtd_ = ibook.book1D("TrackBTLEffEtaMtd", "Track efficiency vs eta (Mtd);#eta_{RECO}", 100, -1.6, 1.6);
0872   meBTLTrackEffPhiMtd_ =
0873       ibook.book1D("TrackBTLEffPhiMtd", "Track efficiency vs phi (Mtd);#phi_{RECO} [rad]", 100, -3.2, 3.2);
0874   meBTLTrackEffPtMtd_ = ibook.book1D("TrackBTLEffPtMtd", "Track efficiency vs pt (Mtd);pt_{RECO} [GeV]", 50, 0, 10);
0875   meBTLTrackPtRes_ =
0876       ibook.book1D("TrackBTLPtRes", "Track pT resolution  ;pT_{Gentrack}-pT_{MTDtrack}/pT_{Gentrack} ", 100, -0.1, 0.1);
0877   meETLTrackRPTime_ = ibook.book1D("TrackETLRPTime", "Track t0 with respect to R.P.;t0 [ns]", 100, -1, 3);
0878   meETLTrackEffEtaTot_[0] =
0879       ibook.book1D("TrackETLEffEtaTotZneg", "Track efficiency vs eta (Tot) (-Z);#eta_{RECO}", 100, -3.2, -1.4);
0880   meETLTrackEffEtaTot_[1] =
0881       ibook.book1D("TrackETLEffEtaTotZpos", "Track efficiency vs eta (Tot) (+Z);#eta_{RECO}", 100, 1.4, 3.2);
0882   meETLTrackEffPhiTot_[0] =
0883       ibook.book1D("TrackETLEffPhiTotZneg", "Track efficiency vs phi (Tot) (-Z);#phi_{RECO} [rad]", 100, -3.2, 3.2);
0884   meETLTrackEffPhiTot_[1] =
0885       ibook.book1D("TrackETLEffPhiTotZpos", "Track efficiency vs phi (Tot) (+Z);#phi_{RECO} [rad]", 100, -3.2, 3.2);
0886   meETLTrackEffPtTot_[0] =
0887       ibook.book1D("TrackETLEffPtTotZneg", "Track efficiency vs pt (Tot) (-Z);pt_{RECO} [GeV]", 50, 0, 10);
0888   meETLTrackEffPtTot_[1] =
0889       ibook.book1D("TrackETLEffPtTotZpos", "Track efficiency vs pt (Tot) (+Z);pt_{RECO} [GeV]", 50, 0, 10);
0890   meETLTrackEffEtaMtd_[0] =
0891       ibook.book1D("TrackETLEffEtaMtdZneg", "Track efficiency vs eta (Mtd) (-Z);#eta_{RECO}", 100, -3.2, -1.4);
0892   meETLTrackEffEtaMtd_[1] =
0893       ibook.book1D("TrackETLEffEtaMtdZpos", "Track efficiency vs eta (Mtd) (+Z);#eta_{RECO}", 100, 1.4, 3.2);
0894   meETLTrackEffPhiMtd_[0] =
0895       ibook.book1D("TrackETLEffPhiMtdZneg", "Track efficiency vs phi (Mtd) (-Z);#phi_{RECO} [rad]", 100, -3.2, 3.2);
0896   meETLTrackEffPhiMtd_[1] =
0897       ibook.book1D("TrackETLEffPhiMtdZpos", "Track efficiency vs phi (Mtd) (+Z);#phi_{RECO} [rad]", 100, -3.2, 3.2);
0898   meETLTrackEffPtMtd_[0] =
0899       ibook.book1D("TrackETLEffPtMtdZneg", "Track efficiency vs pt (Mtd) (-Z);pt_{RECO} [GeV]", 50, 0, 10);
0900   meETLTrackEffPtMtd_[1] =
0901       ibook.book1D("TrackETLEffPtMtdZpos", "Track efficiency vs pt (Mtd) (+Z);pt_{RECO} [GeV]", 50, 0, 10);
0902   meETLTrackEffEta2Mtd_[0] =
0903       ibook.book1D("TrackETLEffEta2MtdZneg", "Track efficiency vs eta (Mtd 2 hit) (-Z);#eta_{RECO}", 100, -3.2, -1.4);
0904   meETLTrackEffEta2Mtd_[1] =
0905       ibook.book1D("TrackETLEffEta2MtdZpos", "Track efficiency vs eta (Mtd 2 hit) (+Z);#eta_{RECO}", 100, 1.4, 3.2);
0906   meETLTrackEffPhi2Mtd_[0] = ibook.book1D(
0907       "TrackETLEffPhi2MtdZneg", "Track efficiency vs phi (Mtd 2 hit) (-Z);#phi_{RECO} [rad]", 100, -3.2, 3.2);
0908   meETLTrackEffPhi2Mtd_[1] = ibook.book1D(
0909       "TrackETLEffPhi2MtdZpos", "Track efficiency vs phi (Mtd 2 hit) (+Z);#phi_{RECO} [rad]", 100, -3.2, 3.2);
0910   meETLTrackEffPt2Mtd_[0] =
0911       ibook.book1D("TrackETLEffPt2MtdZneg", "Track efficiency vs pt (Mtd 2 hit) (-Z);pt_{RECO} [GeV]", 50, 0, 10);
0912   meETLTrackEffPt2Mtd_[1] =
0913       ibook.book1D("TrackETLEffPt2MtdZpos", "Track efficiency vs pt (Mtd 2 hit) (+Z);pt_{RECO} [GeV]", 50, 0, 10);
0914   meETLTrackPtRes_ =
0915       ibook.book1D("TrackETLPtRes", "Track pT resolution;pT_{Gentrack}-pT_{MTDtrack}/pT_{Gentrack} ", 100, -0.1, 0.1);
0916 
0917   meTracktmtd_ = ibook.book1D("Tracktmtd", "Track time from TrackExtenderWithMTD;tmtd [ns]", 150, 1, 16);
0918   meTrackt0Src_ = ibook.book1D("Trackt0Src", "Track time from TrackExtenderWithMTD;t0Src [ns]", 100, -1.5, 1.5);
0919   meTrackSigmat0Src_ =
0920       ibook.book1D("TrackSigmat0Src", "Time Error from TrackExtenderWithMTD; #sigma_{t0Src} [ns]", 100, 0, 0.1);
0921 
0922   meTrackt0Pid_ = ibook.book1D("Trackt0Pid", "Track t0 as stored in TofPid;t0 [ns]", 100, -1, 1);
0923   meTrackSigmat0Pid_ = ibook.book1D("TrackSigmat0Pid", "Sigmat0 as stored in TofPid; #sigma_{t0} [ns]", 100, 0, 0.1);
0924   meTrackt0SafePid_ = ibook.book1D("Trackt0SafePID", "Track t0 Safe as stored in TofPid;t0 [ns]", 100, -1, 1);
0925   meTrackSigmat0SafePid_ =
0926       ibook.book1D("TrackSigmat0SafePID", "Sigmat0 Safe as stored in TofPid; #sigma_{t0} [ns]", 100, 0, 0.1);
0927   meTrackNumHits_ = ibook.book1D("TrackNumHits", "Number of valid MTD hits per track ; Number of hits", 10, -5, 5);
0928   meTrackNumHitsNT_ = ibook.book1D(
0929       "TrackNumHitsNT", "Number of valid MTD hits per track no time associated; Number of hits", 10, -5, 5);
0930   meTrackMVAQual_ = ibook.book1D("TrackMVAQual", "Track MVA Quality as stored in Value Map ; MVAQual", 100, 0, 1);
0931   meTrackPathLenghtvsEta_ = ibook.bookProfile(
0932       "TrackPathLenghtvsEta", "MTD Track pathlength vs MTD track Eta;|#eta|;Pathlength", 100, 0, 3.2, 100.0, 400.0, "S");
0933 
0934   meMVATrackEffPtTot_ = ibook.book1D("MVAEffPtTot", "Pt of tracks associated to LV; track pt [GeV] ", 110, 0., 11.);
0935   meMVATrackMatchedEffPtTot_ =
0936       ibook.book1D("MVAMatchedEffPtTot", "Pt of tracks associated to LV matched to GEN; track pt [GeV] ", 110, 0., 11.);
0937   meMVATrackMatchedEffPtMtd_ = ibook.book1D(
0938       "MVAMatchedEffPtMtd", "Pt of tracks associated to LV matched to GEN with time; track pt [GeV] ", 110, 0., 11.);
0939 
0940   meExtraPtMtd_ = ibook.book1D("ExtraPtMtd", "Pt of tracks extrapolated to hits; track pt [GeV] ", 110, 0., 11.);
0941   meExtraPtEtl2Mtd_ =
0942       ibook.book1D("ExtraPtEtl2Mtd", "Pt of tracks extrapolated to hits, 2 ETL layers; track pt [GeV] ", 110, 0., 11.);
0943 
0944   meTrackPtTot_ = ibook.book1D("TrackPtTot", "Pt of tracks ; track pt [GeV] ", 110, 0., 11.);
0945   meTrackMatchedTPEffPtTot_ =
0946       ibook.book1D("MatchedTPEffPtTot", "Pt of tracks  matched to TP; track pt [GeV] ", 110, 0., 11.);
0947   meTrackMatchedTPEffPtMtd_ =
0948       ibook.book1D("MatchedTPEffPtMtd", "Pt of tracks  matched to TP with time; track pt [GeV] ", 110, 0., 11.);
0949   meTrackMatchedTPEffPtEtl2Mtd_ = ibook.book1D(
0950       "MatchedTPEffPtEtl2Mtd", "Pt of tracks  matched to TP with time, 2 ETL hits; track pt [GeV] ", 110, 0., 11.);
0951 
0952   meTrackMatchedTPmtdEffPtTot_ =
0953       ibook.book1D("MatchedTPmtdEffPtTot", "Pt of tracks  matched to TP-mtd hit; track pt [GeV] ", 110, 0., 11.);
0954   meTrackMatchedTPmtdEffPtMtd_ = ibook.book1D(
0955       "MatchedTPmtdEffPtMtd", "Pt of tracks  matched to TP-mtd hit with time; track pt [GeV] ", 110, 0., 11.);
0956 
0957   meMVATrackEffEtaTot_ = ibook.book1D("MVAEffEtaTot", "Eta of tracks associated to LV; track eta ", 66, 0., 3.3);
0958   meMVATrackMatchedEffEtaTot_ =
0959       ibook.book1D("MVAMatchedEffEtaTot", "Eta of tracks associated to LV matched to GEN; track eta ", 66, 0., 3.3);
0960   meMVATrackMatchedEffEtaMtd_ = ibook.book1D(
0961       "MVAMatchedEffEtaMtd", "Eta of tracks associated to LV matched to GEN with time; track eta ", 66, 0., 3.3);
0962 
0963   meExtraEtaMtd_ = ibook.book1D("ExtraEtaMtd", "Eta of tracks extrapolated to hits; track eta ", 66, 0., 3.3);
0964   meExtraEtaEtl2Mtd_ =
0965       ibook.book1D("ExtraEtaEtl2Mtd", "Eta of tracks extrapolated to hits, 2 ETL layers; track eta ", 66, 0., 3.3);
0966 
0967   meTrackEtaTot_ = ibook.book1D("TrackEtaTot", "Eta of tracks ; track eta ", 66, 0., 3.3);
0968   meTrackMatchedTPEffEtaTot_ =
0969       ibook.book1D("MatchedTPEffEtaTot", "Eta of tracks  matched to TP; track eta ", 66, 0., 3.3);
0970   meMVATrackEffEtaTot_ = ibook.book1D("MVAEffEtaTot", "Eta of tracks ; track eta ", 66, 0., 3.3);
0971   meTrackMatchedTPEffEtaMtd_ =
0972       ibook.book1D("MatchedTPEffEtaMtd", "Eta of tracks  matched to TP with time; track eta ", 66, 0., 3.3);
0973   meTrackMatchedTPEffEtaEtl2Mtd_ = ibook.book1D(
0974       "MatchedTPEffEtaEtl2Mtd", "Eta of tracks  matched to TP with time, 2 ETL hits; track eta ", 66, 0., 3.3);
0975 
0976   meTrackMatchedTPmtdEffEtaTot_ =
0977       ibook.book1D("MatchedTPmtdEffEtaTot", "Eta of tracks  matched to TP-mtd hit; track eta ", 66, 0., 3.3);
0978   meTrackMatchedTPmtdEffEtaMtd_ =
0979       ibook.book1D("MatchedTPmtdEffEtaMtd", "Eta of tracks  matched to TP-mtd hit with time; track eta ", 66, 0., 3.3);
0980 
0981   meMVATrackResTot_ = ibook.book1D(
0982       "MVATrackRes", "t_{rec} - t_{sim} for LV associated tracks; t_{rec} - t_{sim} [ns] ", 120, -0.15, 0.15);
0983   meMVATrackPullTot_ =
0984       ibook.book1D("MVATrackPull", "Pull for associated tracks; (t_{rec}-t_{sim})/#sigma_{t}", 50, -5., 5.);
0985   meMVATrackZposResTot_ = ibook.book1D(
0986       "MVATrackZposResTot", "Z_{PCA} - Z_{sim} for associated tracks;Z_{PCA} - Z_{sim} [cm] ", 100, -0.1, 0.1);
0987 
0988   meExtraPhiAtBTL_ =
0989       ibook.book1D("ExtraPhiAtBTL", "Phi at BTL surface of extrapolated tracks; phi [deg]", 720, -180., 180.);
0990   meExtraPhiAtBTLmatched_ = ibook.book1D("ExtraPhiAtBTLmatched",
0991                                          "Phi at BTL surface of extrapolated tracksi matched with BTL hits; phi [deg]",
0992                                          720,
0993                                          -180.,
0994                                          180.);
0995   meExtraBTLeneInCone_ = ibook.book1D(
0996       "ExtraBTLeneInCone", "BTL reconstructed energy in cone arounnd extrapolated track; E [MeV]", 100, 0., 50.);
0997   meExtraMTDfailExtenderEta_ =
0998       ibook.book1D("ExtraMTDfailExtenderEta",
0999                    "Eta of tracks extrapolated to MTD with no track extender match to hits; track eta",
1000                    66,
1001                    0.,
1002                    3.3);
1003   ;
1004   meExtraMTDfailExtenderPt_ =
1005       ibook.book1D("ExtraMTDfailExtenderPt",
1006                    "Pt of tracks extrapolated to MTD with no track extender match to hits; track pt [GeV] ",
1007                    110,
1008                    0.,
1009                    11.);
1010 }
1011 
1012 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
1013 
1014 void MtdTracksValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
1015   edm::ParameterSetDescription desc;
1016 
1017   desc.add<std::string>("folder", "MTD/Tracks");
1018   desc.add<edm::InputTag>("inputTagG", edm::InputTag("generalTracks"));
1019   desc.add<edm::InputTag>("inputTagT", edm::InputTag("trackExtenderWithMTD"));
1020   desc.add<edm::InputTag>("inputTagV", edm::InputTag("offlinePrimaryVertices4D"));
1021   desc.add<edm::InputTag>("inputTagH", edm::InputTag("generatorSmeared"));
1022   desc.add<edm::InputTag>("SimTag", edm::InputTag("mix", "MergedTrackTruth"));
1023   desc.add<edm::InputTag>("TPtoRecoTrackAssoc", edm::InputTag("trackingParticleRecoTrackAsssociation"));
1024   desc.add<edm::InputTag>("btlSimHits", edm::InputTag("mix", "g4SimHitsFastTimerHitsBarrel"));
1025   desc.add<edm::InputTag>("etlSimHits", edm::InputTag("mix", "g4SimHitsFastTimerHitsEndcap"));
1026   desc.add<edm::InputTag>("btlRecHits", edm::InputTag("mtdRecHits", "FTLBarrel"));
1027   desc.add<edm::InputTag>("etlRecHits", edm::InputTag("mtdRecHits", "FTLEndcap"));
1028   desc.add<edm::InputTag>("tmtd", edm::InputTag("trackExtenderWithMTD:generalTracktmtd"));
1029   desc.add<edm::InputTag>("sigmatmtd", edm::InputTag("trackExtenderWithMTD:generalTracksigmatmtd"));
1030   desc.add<edm::InputTag>("t0Src", edm::InputTag("trackExtenderWithMTD:generalTrackt0"));
1031   desc.add<edm::InputTag>("sigmat0Src", edm::InputTag("trackExtenderWithMTD:generalTracksigmat0"));
1032   desc.add<edm::InputTag>("trackAssocSrc", edm::InputTag("trackExtenderWithMTD:generalTrackassoc"))
1033       ->setComment("Association between General and MTD Extended tracks");
1034   desc.add<edm::InputTag>("pathLengthSrc", edm::InputTag("trackExtenderWithMTD:generalTrackPathLength"));
1035   desc.add<edm::InputTag>("t0SafePID", edm::InputTag("tofPID:t0safe"));
1036   desc.add<edm::InputTag>("sigmat0SafePID", edm::InputTag("tofPID:sigmat0safe"));
1037   desc.add<edm::InputTag>("sigmat0PID", edm::InputTag("tofPID:sigmat0"));
1038   desc.add<edm::InputTag>("t0PID", edm::InputTag("tofPID:t0"));
1039   desc.add<edm::InputTag>("trackMVAQual", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA"));
1040   desc.add<double>("trackMinimumPt", 0.7);  // [GeV]
1041   desc.add<double>("trackMaximumBtlEta", 1.5);
1042   desc.add<double>("trackMinimumEtlEta", 1.6);
1043   desc.add<double>("trackMaximumEtlEta", 3.);
1044   desc.addUntracked<bool>("optionalPlots", true);
1045 
1046   descriptions.add("mtdTracksValid", desc);
1047 }
1048 
1049 const bool MtdTracksValidation::mvaGenSel(const HepMC::GenParticle& gp, const float& charge) {
1050   bool match = false;
1051   if (gp.status() != 1) {
1052     return match;
1053   }
1054   match = charge != 0.f && gp.momentum().perp() > pTcut_ && std::abs(gp.momentum().eta()) < etacutGEN_;
1055   return match;
1056 }
1057 
1058 const bool MtdTracksValidation::mvaTPSel(const TrackingParticle& tp) {
1059   bool match = false;
1060   if (tp.status() != 1) {
1061     return match;
1062   }
1063   auto x_pv = tp.parentVertex()->position().x();
1064   auto y_pv = tp.parentVertex()->position().y();
1065   auto z_pv = tp.parentVertex()->position().z();
1066 
1067   auto r_pv = std::sqrt(x_pv * x_pv + y_pv * y_pv);
1068 
1069   match = tp.charge() != 0 && tp.pt() > pTcut_ && std::abs(tp.eta()) < etacutGEN_ && r_pv < rBTL_ && z_pv < zETL_;
1070   return match;
1071 }
1072 
1073 const bool MtdTracksValidation::mvaRecSel(const reco::TrackBase& trk,
1074                                           const reco::Vertex& vtx,
1075                                           const double& t0,
1076                                           const double& st0) {
1077   bool match = false;
1078   match = trk.pt() > pTcut_ && std::abs(trk.eta()) < etacutREC_ &&
1079           (std::abs(trk.vz() - vtx.z()) <= deltaZcut_ || vtx.isFake());
1080   if (st0 > 0.) {
1081     match = match && std::abs(t0 - vtx.t()) < 3. * st0;
1082   }
1083   return match;
1084 }
1085 
1086 const bool MtdTracksValidation::mvaGenRecMatch(const HepMC::GenParticle& genP,
1087                                                const double& zsim,
1088                                                const reco::TrackBase& trk,
1089                                                const bool& vtxFake) {
1090   bool match = false;
1091   double dR = reco::deltaR(genP.momentum(), trk.momentum());
1092   double genPT = genP.momentum().perp();
1093   match = std::abs(genPT - trk.pt()) < trk.pt() * deltaPTcut_ && dR < deltaDRcut_ &&
1094           (std::abs(trk.vz() - zsim) < deltaZcut_ || vtxFake);
1095   return match;
1096 }
1097 
1098 const edm::Ref<std::vector<TrackingParticle>>* MtdTracksValidation::getMatchedTP(const reco::TrackBaseRef& recoTrack) {
1099   auto found = r2s_->find(recoTrack);
1100 
1101   // reco track not matched to any TP
1102   if (found == r2s_->end())
1103     return nullptr;
1104 
1105   //matched TP equal to any TP associated to in time events
1106   for (const auto& tp : found->val) {
1107     if (tp.first->eventId().bunchCrossing() == 0)
1108       return &tp.first;
1109   }
1110 
1111   // reco track not matched to any TP from vertex
1112   return nullptr;
1113 }
1114 
1115 DEFINE_FWK_MODULE(MtdTracksValidation);