Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:27:24

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