File indexing completed on 2022-12-03 01:10:34
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
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.;
0089 static constexpr double etacutREC_ = 3.;
0090 static constexpr double pTcut_ = 0.7;
0091 static constexpr double deltaZcut_ = 0.1;
0092 static constexpr double deltaPTcut_ = 0.05;
0093 static constexpr double deltaDRcut_ = 0.03;
0094 static constexpr double depositBTLthreshold_ = 1;
0095 static constexpr double depositETLthreshold_ = 0.001;
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
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
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
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
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
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
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
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
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
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 }
0429 if (!tpIsAssoc && goodCell) {
0430 meUnassDeposit_->Fill(0.5);
0431 }
0432
0433 } else {
0434
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 }
0469 if (!tpIsAssoc && goodCell) {
0470 meUnassDeposit_->Fill(1.5);
0471 }
0472 }
0473 }
0474 }
0475 }
0476
0477 unsigned int index = 0;
0478
0479
0480 const bool isGoodVtx = std::abs(primRecoVtx.z() - zsim) < deltaZcut_ || primRecoVtx.isFake();
0481
0482
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
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
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 }
0546
0547 else {
0548
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
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
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
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
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 }
0729 }
0730 }
0731
0732
0733 void MtdTracksValidation::bookHistograms(DQMStore::IBooker& ibook, edm::Run const& run, edm::EventSetup const& iSetup) {
0734 ibook.setCurrentFolder(folder_);
0735
0736
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
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);
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
0973 if (found == r2s_->end())
0974 return nullptr;
0975
0976
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
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);
