Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-30 04:06:16

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