Back to home page

Project CMSSW displayed by LXR

 
 

    


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     //edm::LogVerbatim("CkfDebugger") << "CkfDebugger::fillSeedHist";
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   //const PSimHit* nextCorrectHit( const Trajectory&, unsigned int&) ;
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       /*       edm::LogVerbatim("CkfDebugger") << "SimHit on this det at pos="<< position(&*isim)  */
0149       /*         << " det=" << isim->detUnitId() << " process=" << isim->processType() ; */
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     //return ((int)(((det->geographicalId().rawId() >>16) & 0xF)));
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