File indexing completed on 2023-03-17 11:22:10
0001 #ifndef CkfDebugger_H
0002 #define CkfDebugger_H
0003
0004 #include "DataFormats/Common/interface/Handle.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Utilities/interface/ESGetToken.h"
0009
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011
0012 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0013 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0014 #include "MagneticField/Engine/interface/MagneticField.h"
0015
0016 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
0017 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0018 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0019 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0020 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0021 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0022 #include "TrackingTools/PatternTools/interface/MeasurementExtractor.h"
0023 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0024 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0025
0026 #include <vector>
0027 #include <iostream>
0028 #include <TFile.h>
0029 #include <TH1F.h>
0030 #include <TH2F.h>
0031
0032 class Trajectory;
0033 class TrajectoryMeasurement;
0034 class PSimHit;
0035 class MeasurementTrackerEvent;
0036 class TrajectoryStateOnSurface;
0037 class MagneticField;
0038 class Chi2MeasurementEstimator;
0039 class Propagator;
0040 class NavigationSchool;
0041 class TrackerDigiGeometryRecord;
0042 class IdealMagneticFieldRecord;
0043 class IdealGeometryRecord;
0044 class NavigationSchoolRecord;
0045
0046 typedef TransientTrackingRecHit::ConstRecHitPointer CTTRHp;
0047
0048 class CkfDebugger {
0049 public:
0050 CkfDebugger(edm::ConsumesCollector iC);
0051
0052 ~CkfDebugger();
0053
0054 void setConditions(edm::EventSetup const& es);
0055
0056 void printSimHits(const edm::Event& iEvent);
0057
0058 void countSeed() { totSeeds++; }
0059
0060 void fillSeedHist(CTTRHp h1, CTTRHp h2, TrajectoryStateOnSurface t) {
0061
0062 hchi2seedAll->Fill(testSeed(h1, h2, t));
0063 }
0064
0065 bool analyseCompatibleMeasurements(const Trajectory&,
0066 const std::vector<TrajectoryMeasurement>&,
0067 const MeasurementTrackerEvent*,
0068 const Propagator*,
0069 const Chi2MeasurementEstimatorBase*,
0070 const TransientTrackingRecHitBuilder*);
0071
0072 void deleteHitAssociator() { delete hitAssociator; }
0073
0074 private:
0075 typedef TrajectoryMeasurement TM;
0076 typedef TrajectoryStateOnSurface TSOS;
0077
0078 class SimHit {
0079 public:
0080 SimHit(const PSimHit* phit, const GeomDetUnit* gdu) : thePHit(phit), theDet(gdu) {}
0081 LocalPoint localPosition() const { return thePHit->localPosition(); }
0082 GlobalPoint globalPosition() const { return theDet->toGlobal(thePHit->localPosition()); }
0083 const GeomDetUnit* det() const { return theDet; }
0084 unsigned int trackId() const { return thePHit->trackId(); }
0085 LocalVector localDirection() const { return thePHit->localDirection(); }
0086 Geom::Theta<float> thetaAtEntry() const { return thePHit->thetaAtEntry(); }
0087 Geom::Phi<float> phiAtEntry() const { return thePHit->phiAtEntry(); }
0088 float pabs() const { return thePHit->pabs(); }
0089 float timeOfFlight() const { return thePHit->timeOfFlight(); }
0090 float energyLoss() const { return thePHit->energyLoss(); }
0091 int particleType() const { return thePHit->particleType(); }
0092 unsigned int detUnitId() const { return thePHit->detUnitId(); }
0093 unsigned short processType() const { return thePHit->processType(); }
0094 const PSimHit& psimHit() const { return *thePHit; }
0095
0096 private:
0097 const PSimHit* thePHit;
0098 const GeomDetUnit* theDet;
0099 };
0100
0101 const TrackerGeometry* theTrackerGeom;
0102 const MagneticField* theMagField;
0103 const GeometricSearchTracker* theGeomSearchTracker;
0104 const MeasurementEstimator* theChi2;
0105 const Propagator* theForwardPropagator;
0106 edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> theTrackerToken;
0107 edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> theFieldToken;
0108 edm::ESGetToken<TrackerTopology, IdealGeometryRecord> theTopoHandToken;
0109 edm::ESGetToken<NavigationSchool, NavigationSchoolRecord> theNavToken;
0110
0111 TrackerHitAssociator::Config trackerHitAssociatorConfig_;
0112 TrackerHitAssociator* hitAssociator;
0113 const MeasurementTrackerEvent* theMeasurementTracker;
0114 const TransientTrackingRecHitBuilder* theTTRHBuilder;
0115 const TrackerTopology* theTopo;
0116 NavigationSchool const* theNavSchool;
0117
0118 std::map<unsigned int, std::vector<PSimHit*> > idHitsMap;
0119
0120 void dumpSimHit(const SimHit& hit) const;
0121
0122 bool correctTrajectory(const Trajectory&, unsigned int&) const;
0123
0124 int assocTrackId(CTTRHp rechit) const;
0125
0126
0127 std::vector<const PSimHit*> nextCorrectHits(const Trajectory&, unsigned int&);
0128
0129 bool associated(CTTRHp rechit, const PSimHit& sh) const;
0130
0131 bool goodSimHit(const PSimHit& sh) const;
0132
0133 bool correctMeas(const TM& tm, const PSimHit* correctHit) const;
0134
0135 std::pair<CTTRHp, double> analyseRecHitExistance(const PSimHit& sh, const TSOS& startingState);
0136
0137 int analyseRecHitNotFound(const Trajectory&, CTTRHp);
0138
0139 double testSeed(CTTRHp, CTTRHp, TrajectoryStateOnSurface);
0140
0141 const PSimHit* pSimHit(unsigned int tkId, DetId detId);
0142
0143 bool hasDelta(const PSimHit* correctHit) {
0144 bool delta = false;
0145 for (std::vector<PSimHit>::iterator isim = hitAssociator->SimHitMap[correctHit->detUnitId()].begin();
0146 isim != hitAssociator->SimHitMap[correctHit->detUnitId()].end();
0147 ++isim) {
0148
0149
0150 if (isim->processType() == 9)
0151 delta = true;
0152 }
0153 return delta;
0154 }
0155
0156 Global3DPoint position(const PSimHit* sh) const {
0157 return theTrackerGeom->idToDetUnit(DetId(sh->detUnitId()))->toGlobal(sh->localPosition());
0158 };
0159 const GeomDetUnit* det(const PSimHit* sh) const { return theTrackerGeom->idToDetUnit(DetId(sh->detUnitId())); };
0160
0161 int layer(const GeomDet* det) {
0162
0163 return theTopo->layer(det->geographicalId());
0164 }
0165
0166 template <unsigned int D>
0167 std::pair<double, double> computePulls(CTTRHp recHit, TSOS startingState) {
0168 typedef typename AlgebraicROOTObject<D>::Vector VecD;
0169 typedef typename AlgebraicROOTObject<D, D>::SymMatrix SMatDD;
0170 TSOS detState = theForwardPropagator->propagate(startingState, recHit->det()->surface());
0171 LogTrace("CkfDebugger") << "parameters=" << recHit->parameters();
0172 LogTrace("CkfDebugger") << "parametersError=" << recHit->parametersError();
0173 MeasurementExtractor me(detState);
0174 VecD r = asSVector<D>(recHit->parameters()) - me.measuredParameters<D>(*recHit);
0175 LogTrace("CkfDebugger") << "me.measuredParameters=" << me.measuredParameters<D>(*recHit);
0176 LogTrace("CkfDebugger") << "me.measuredError=" << me.measuredError<D>(*recHit);
0177 SMatDD R = asSMatrix<D>(recHit->parametersError()) + me.measuredError<D>(*recHit);
0178 LogTrace("CkfDebugger") << "r=" << r;
0179 LogTrace("CkfDebugger") << "R=" << R;
0180 R.Invert();
0181 LogTrace("CkfDebugger") << "R(-1)=" << R;
0182 LogTrace("CkfDebugger") << "chi2=" << ROOT::Math::Similarity(r, R);
0183 double pullX = (-r[0]) * sqrt(R(0, 0));
0184 double r_1 = 0;
0185 if (VecD::Dim() >= 2) {
0186 r_1 = r[1];
0187 }
0188 double pullY = (-r_1) * sqrt(R(1, 1));
0189 LogTrace("CkfDebugger") << "pullX=" << pullX;
0190 LogTrace("CkfDebugger") << "pullY=" << pullY;
0191 return std::pair<double, double>(pullX, pullY);
0192 }
0193 std::pair<double, double> computePulls(CTTRHp recHit, TSOS startingState) {
0194 switch (recHit->dimension()) {
0195 case 1:
0196 return computePulls<1>(recHit, startingState);
0197 case 2:
0198 return computePulls<2>(recHit, startingState);
0199 case 3:
0200 return computePulls<3>(recHit, startingState);
0201 case 4:
0202 return computePulls<4>(recHit, startingState);
0203 case 5:
0204 return computePulls<5>(recHit, startingState);
0205 }
0206 throw cms::Exception("CkfDebugger error: rechit of dimension not 1,2,3,4,5");
0207 }
0208
0209 std::vector<int> dump;
0210 std::map<std::pair<int, int>, int> dump2;
0211 std::map<std::pair<int, int>, int> dump3;
0212 std::map<std::pair<int, int>, int> dump4;
0213 std::map<std::pair<int, int>, int> dump5;
0214 std::map<std::pair<int, int>, int> dump6;
0215
0216 TFile* file;
0217 TH1F *hchi2seedAll, *hchi2seedProb;
0218
0219 std::map<std::string, TH1F*> hPullX_shrh;
0220 std::map<std::string, TH1F*> hPullY_shrh;
0221 std::map<std::string, TH1F*> hPullX_shst;
0222 std::map<std::string, TH1F*> hPullY_shst;
0223 std::map<std::string, TH1F*> hPullX_strh;
0224 std::map<std::string, TH1F*> hPullY_strh;
0225
0226 std::map<std::string, TH1F*> hPullM_shrh;
0227 std::map<std::string, TH1F*> hPullS_shrh;
0228 std::map<std::string, TH1F*> hPullM_shst;
0229 std::map<std::string, TH1F*> hPullS_shst;
0230 std::map<std::string, TH1F*> hPullM_strh;
0231 std::map<std::string, TH1F*> hPullS_strh;
0232
0233 std::map<std::string, TH1F*> hPullGP_X_shst;
0234 std::map<std::string, TH1F*> hPullGP_Y_shst;
0235 std::map<std::string, TH1F*> hPullGP_Z_shst;
0236
0237 TH2F* hPullGPXvsGPX_shst;
0238 TH2F* hPullGPXvsGPY_shst;
0239 TH2F* hPullGPXvsGPZ_shst;
0240 TH2F* hPullGPXvsGPr_shst;
0241 TH2F* hPullGPXvsGPeta_shst;
0242 TH2F* hPullGPXvsGPphi_shst;
0243
0244 int seedWithDelta;
0245 int problems;
0246 int no_sim_hit;
0247 int no_layer;
0248 int layer_not_found;
0249 int det_not_found;
0250 int chi2gt30;
0251 int chi2gt30delta;
0252 int chi2gt30deltaSeed;
0253 int chi2ls30;
0254 int simple_hit_not_found;
0255 int no_component;
0256 int only_one_component;
0257 int matched_not_found;
0258 int matched_not_associated;
0259 int partner_det_not_fuond;
0260 int glued_det_not_fuond;
0261 int propagation;
0262 int other;
0263 int totchi2gt30;
0264
0265 int totSeeds;
0266 };
0267
0268 #endif
