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
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.;
0117 static constexpr double etacutREC_ = 3.;
0118 static constexpr double pTcut_ = 0.7;
0119 static constexpr double deltaZcut_ = 0.1;
0120 static constexpr double deltaPTcut_ = 0.05;
0121 static constexpr double deltaDRcut_ = 0.03;
0122 static constexpr double depositBTLthreshold_ = 1;
0123 static constexpr double depositETLthreshold_ = 0.001;
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;
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
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
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
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
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
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
0343 (simHitIt->second).energy += convertUnitsTo(0.001_MeV, simHit.energyLoss());
0344 }
0345
0346
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
0387 const bool isGoodVtx = std::abs(primRecoVtx.z() - zsim) < deltaZcut_ || primRecoVtx.isFake();
0388
0389
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
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
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 }
0458
0459 else {
0460
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
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
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 }
0624 }
0625
0626 if (isGoodVtx) {
0627 const bool vtxFake = primRecoVtx.isFake();
0628
0629 if (mvaRecSel(trackGen, primRecoVtx, t0Safe[trackref], Sigmat0Safe[trackref])) {
0630
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
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 }
0671 }
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
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
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
0765
0766
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
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;
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
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
0836
0837
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
0862 void MtdTracksValidation::bookHistograms(DQMStore::IBooker& ibook, edm::Run const& run, edm::EventSetup const& iSetup) {
0863 ibook.setCurrentFolder(folder_);
0864
0865
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
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);
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
1102 if (found == r2s_->end())
1103 return nullptr;
1104
1105
1106 for (const auto& tp : found->val) {
1107 if (tp.first->eventId().bunchCrossing() == 0)
1108 return &tp.first;
1109 }
1110
1111
1112 return nullptr;
1113 }
1114
1115 DEFINE_FWK_MODULE(MtdTracksValidation);
