SiStripHitResolution/src/HitResol.cc

0001 ////////////////////////////////////////////////////////////////////////////////
0002 // Package:          CalibTracker/SiStripHitResolution
0003 // Class:            HitResol
0004 // Original Authors: Denis Gele and Kathryn Coldham (adapted from HitEff)
0005 //                   modified by Khawla Jaffel for CPE studies
0006 //                   ported to cmssw by M. Musich
0007 //
0008 ///////////////////////////////////////////////////////////////////////////////
0009
0010 // system include files
0011 #include <memory>
0012 #include <string>
0013 #include <iostream>
0014
0015 // user includes
0016 #include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h"
0017 #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h"
0018 #include "CalibTracker/Records/interface/SiStripDetCablingRcd.h"
0019 #include "CalibTracker/Records/interface/SiStripQualityRcd.h"
0020 #include "CalibTracker/SiStripHitEfficiency/interface/SiStripHitEfficiencyHelpers.h"
0021 #include "CalibTracker/SiStripHitResolution/interface/HitResol.h"
0022 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0023 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0024 #include "DataFormats/Common/interface/Handle.h"
0025 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0026 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0027 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0028 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0029 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0030 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0031 #include "DataFormats/TrackReco/interface/Track.h"
0032 #include "DataFormats/TrackReco/interface/TrackBase.h"
0033 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0034 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0035 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0036 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0037 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0038 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0039 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0040 #include "FWCore/Framework/interface/ESHandle.h"
0041 #include "FWCore/Framework/interface/Event.h"
0042 #include "FWCore/Framework/interface/EventSetup.h"
0043 #include "FWCore/Framework/interface/Frameworkfwd.h"
0044 #include "FWCore/Framework/interface/MakerMacros.h"
0045 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0046 #include "FWCore/ServiceRegistry/interface/Service.h"
0047 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0048 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0049 #include "Geometry/CommonDetUnit/interface/GluedGeomDet.h"
0050 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0051 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0052 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0053 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0054 #include "RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h"
0055 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0056 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0057 #include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
0058 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCurvilinearToLocal.h"
0059 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h"
0060 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0061 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0062
0063 // ROOT includes
0064 #include "TMath.h"
0065 #include "TH1F.h"
0066
0067 //
0068 // constructors and destructor
0069 //
0070 using namespace std;
0071 HitResol::HitResol(const edm::ParameterSet& conf)
0072     : scalerToken_(consumes<LumiScalersCollection>(conf.getParameter<edm::InputTag>("lumiScalers"))),
0073       combinatorialTracks_token_(
0074           consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("combinatorialTracks"))),
0075       tjToken_(consumes<std::vector<Trajectory> >(conf.getParameter<edm::InputTag>("trajectories"))),
0076       topoToken_(esConsumes()),
0077       geomToken_(esConsumes()),
0078       cpeToken_(esConsumes(edm::ESInputTag("", "StripCPEfromTrackAngle"))),
0079       siStripQualityToken_(esConsumes()),
0080       magFieldToken_(esConsumes()),
0081       addLumi_(conf.getUntrackedParameter<bool>("addLumi", false)),
0082       DEBUG_(conf.getParameter<bool>("Debug")),
0083       cutOnTracks_(conf.getUntrackedParameter<bool>("cutOnTracks", false)),
0084       momentumCut_(conf.getUntrackedParameter<double>("MomentumCut", 3.)),
0085       compSettings_(conf.getUntrackedParameter<int>("CompressionSettings", -1)),
0086       usePairsOnly_(conf.getUntrackedParameter<unsigned int>("UsePairsOnly", 1)),
0087       layers_(conf.getParameter<int>("Layer")),
0088       trackMultiplicityCut_(conf.getUntrackedParameter<unsigned int>("trackMultiplicity", 100)) {
0089   usesResource(TFileService::kSharedResource);
0090 }
0091
0092 void HitResol::beginJob() {
0093   edm::Service<TFileService> fs;
0094   if (compSettings_ > 0) {
0095     edm::LogInfo("SiStripHitResolution:HitResol") << "the compressions settings are:" << compSettings_ << std::endl;
0096     fs->file().SetCompressionSettings(compSettings_);
0097   }
0098
0099   reso = fs->make<TTree>("reso", "tree hit pairs for resolution studies");
0100   reso->Branch("momentum", &mymom, "momentum/F");
0101   reso->Branch("numHits", &numHits, "numHits/I");
0102   reso->Branch("trackChi2", &ProbTrackChi2, "trackChi2/F");
0103   reso->Branch("detID1", &iidd1, "detID1/I");
0104   reso->Branch("pitch1", &mypitch1, "pitch1/F");
0105   reso->Branch("clusterW1", &clusterWidth, "clusterW1/I");
0106   reso->Branch("expectedW1", &expWidth, "expectedW1/F");
0107   reso->Branch("atEdge1", &atEdge, "atEdge1/F");
0108   reso->Branch("simpleRes", &simpleRes, "simpleRes/F");
0109   reso->Branch("detID2", &iidd2, "detID2/I");
0110   reso->Branch("clusterW2", &clusterWidth_2, "clusterW2/I");
0111   reso->Branch("expectedW2", &expWidth_2, "expectedW2/F");
0112   reso->Branch("atEdge2", &atEdge_2, "atEdge2/F");
0113   reso->Branch("pairPath", &pairPath, "pairPath/F");
0114   reso->Branch("hitDX", &hitDX, "hitDX/F");
0115   reso->Branch("trackDX", &trackDX, "trackDX/F");
0116   reso->Branch("trackDXE", &trackDXE, "trackDXE/F");
0117   reso->Branch("trackParamX", &trackParamX, "trackParamX/F");
0118   reso->Branch("trackParamY", &trackParamY, "trackParamY/F");
0119   reso->Branch("trackParamDXDZ", &trackParamDXDZ, "trackParamDXDZ/F");
0120   reso->Branch("trackParamDYDZ", &trackParamDYDZ, "trackParamDYDZ/F");
0121   reso->Branch("trackParamXE", &trackParamXE, "trackParamXE/F");
0122   reso->Branch("trackParamYE", &trackParamYE, "trackParamYE/F");
0123   reso->Branch("trackParamDXDZE", &trackParamDXDZE, "trackParamDXDZE/F");
0124   reso->Branch("trackParamDYDZE", &trackParamDYDZE, "trackParamDYDZE/F");
0125   reso->Branch("pairsOnly", &pairsOnly, "pairsOnly/I");
0126   treso = fs->make<TTree>("treso", "tree tracks  for resolution studies");
0127   treso->Branch("track_momentum", &track_momentum, "track_momentum/F");
0128   treso->Branch("track_pt", &track_pt, "track_pt/F");
0129   treso->Branch("track_eta", &track_eta, "track_eta/F");
0130   treso->Branch("track_phi", &track_phi, "track_phi/F");
0131   treso->Branch("track_trackChi2", &track_trackChi2, "track_trackChi2/F");
0132   treso->Branch("track_width", &expWidth, "track_width/F");  // from 1D HIT
0133   treso->Branch("NumberOf_tracks", &NumberOf_tracks, "NumberOf_tracks/I");
0134
0135   events = 0;
0136   EventTrackCKF = 0;
0137
0138   histos2d_["track_phi_vs_eta"] = new TH2F("track_phi_vs_eta", ";track phi;track eta", 60, -3.5, 3.5, 60, -3., 3.);
0139   histos2d_["residual_vs_trackMomentum"] = new TH2F("residual_vs_trackMomentum",
0140                                                     ";track momentum [GeV]; x_{pred_track} - x_{reco_hit} [#mum]",
0141                                                     60,
0142                                                     0.,
0143                                                     10.,
0144                                                     60,
0145                                                     0.,
0146                                                     200.);
0147   histos2d_["residual_vs_trackPt"] = new TH2F(
0148       "residual_vs_trackPt", ";track p_{T}[GeV];x_{pred_track} - x_{reco_hit} [#mum]", 60, 0., 10., 60, 0., 200.);
0149   histos2d_["residual_vs_trackEta"] =
0150       new TH2F("residual_vs_trackEta", ";track #eta;x_{pred_track} - x_{reco_hit} [#mum]", 60, 0., 3., 60, 0., 200.);
0151   histos2d_["residual_vs_trackPhi"] =
0152       new TH2F("residual_vs_trackPhi", ";track #phi;x_{pred_track} - x_{reco_hit} [#mum]", 60, 0., 3.5, 60, 0., 200.);
0153   histos2d_["residual_vs_expectedWidth"] = new TH2F(
0154       "residual_vs_expectedWidth", ";track Width;x_{pred_track} - x_{reco_hit} [#mum]", 3, 0., 3., 60, 0., 200.);
0155   histos2d_["numHits_vs_residual"] =
0156       new TH2F("numHits_vs_residual", ";x_{pred_track} - x_{reco_hit} [#mum];N Hits", 60, 0., 200., 15, 0., 15.);
0157 }
0158
0159 void HitResol::analyze(const edm::Event& e, const edm::EventSetup& es) {
0160   //Retrieve tracker topology from geometry
0161   const TrackerTopology* tTopo = &es.getData(topoToken_);
0162
0163   LogDebug("SiStripHitResolution:HitResol") << "beginning analyze from HitResol" << endl;
0164
0165   using namespace edm;
0166   using namespace reco;
0167
0168   // Step A: Get Inputs
0169
0170   int run_nr = e.id().run();
0171   int ev_nr = e.id().event();
0172
0173   // get the tracks
0174   edm::Handle<reco::TrackCollection> trackCollectionCKF;
0175   e.getByToken(combinatorialTracks_token_, trackCollectionCKF);
0176   const reco::TrackCollection* tracksCKF = trackCollectionCKF.product();
0177
0178   // get the trajectory collection
0179   edm::Handle<std::vector<Trajectory> > trajectoryCollectionHandle;
0180   e.getByToken(tjToken_, trajectoryCollectionHandle);
0181   const TrajectoryCollection* trajectoryCollection = trajectoryCollectionHandle.product();
0182
0183   //get tracker geometry
0184   edm::ESHandle<TrackerGeometry> tracker = es.getHandle(geomToken_);
0185   const TrackerGeometry* tkgeom = &(*tracker);
0186
0187   //get Cluster Parameter Estimator
0188   edm::ESHandle<StripClusterParameterEstimator> parameterestimator = es.getHandle(cpeToken_);
0189   const StripClusterParameterEstimator& stripcpe(*parameterestimator);
0190
0191   // get the SiStripQuality records
0192   edm::ESHandle<SiStripQuality> SiStripQuality_ = es.getHandle(siStripQualityToken_);
0193
0194   // get the magnetic field
0195   const MagneticField* magField_ = &es.getData(magFieldToken_);
0196
0197   events++;
0198
0199   // List of variables for SiStripHitResolution ntuple
0200   mymom = 0;
0201   numHits = 0;
0202   ProbTrackChi2 = 0;
0203   iidd1 = 0;
0204   mypitch1 = 0;
0205   clusterWidth = 0;
0206   expWidth = 0;
0207   atEdge = 0;
0208   simpleRes = 0;
0209   iidd2 = 0;
0210   clusterWidth_2 = 0;
0211   expWidth_2 = 0;
0212   atEdge_2 = 0;
0213   pairPath = 0;
0214   hitDX = 0;
0215   trackDX = 0;
0216   trackDXE = 0;
0217   trackParamX = 0;
0218   trackParamY = 0;
0219   trackParamDXDZ = 0;
0220   trackParamDYDZ = 0;
0221   trackParamXE = 0;
0222   trackParamYE = 0;
0223   trackParamDXDZE = 0;
0224   trackParamDYDZE = 0;
0225   pairsOnly = 0;
0226
0227   LogDebug("HitResol") << "Starting analysis, nrun nevent, tracksCKF->size(): " << run_nr << " " << ev_nr << " "
0228                        << tracksCKF->size() << std::endl;
0229
0230   for (unsigned int iT = 0; iT < tracksCKF->size(); ++iT) {
0231     track_momentum = tracksCKF->at(iT).p();
0232     track_pt = tracksCKF->at(iT).p();
0233     track_eta = tracksCKF->at(iT).eta();
0234     track_phi = tracksCKF->at(iT).phi();
0235     track_trackChi2 = ChiSquaredProbability((double)(tracksCKF->at(iT).chi2()), (double)(tracksCKF->at(iT).ndof()));
0236     treso->Fill();
0237   }
0238
0239   histos2d_["track_phi_vs_eta"]->Fill(track_phi, track_eta);
0240
0241   // loop over trajectories from refit
0242   for (const auto& traj : *trajectoryCollection) {
0243     const auto& TMeas = traj.measurements();
0244     // Loop on each measurement and take it into consideration
0245     //--------------------------------------------------------
0246     for (auto itm = TMeas.cbegin(); itm != TMeas.cend(); ++itm) {
0247       if (!itm->updatedState().isValid()) {
0248         LogDebug("HitResol") << "trajectory measurement not valid" << std::endl;
0249         continue;
0250       }
0251
0252       const TransientTrackingRecHit::ConstRecHitPointer mypointhit = itm->recHit();
0253       const TrackingRecHit* myhit = (*itm->recHit()).hit();
0254
0255       ProbTrackChi2 = 0;
0256       numHits = 0;
0257
0258       LogDebug("HitResol") << "TrackChi2 =  "
0259                            << ChiSquaredProbability((double)(traj.chiSquared()), (double)(traj.ndof(false))) << "\n"
0260                            << "itm->updatedState().globalMomentum().perp(): "
0261                            << itm->updatedState().globalMomentum().perp() << "\n"
0262                            << "numhits " << traj.foundHits() << std::endl;
0263
0264       numHits = traj.foundHits();
0265       ProbTrackChi2 = ChiSquaredProbability((double)(traj.chiSquared()), (double)(traj.ndof(false)));
0266
0267       mymom = itm->updatedState().globalMomentum().perp();
0268
0269       //Now for the first hit
0270       TrajectoryStateOnSurface mytsos = itm->updatedState();
0271       const auto hit1 = itm->recHit();
0272       DetId id1 = hit1->geographicalId();
0273       if (id1.subdetId() < StripSubdetector::TIB || id1.subdetId() > StripSubdetector::TEC)
0274         continue;
0275
0276       if (hit1->isValid() && mymom > momentumCut_ &&
0277           (id1.subdetId() >= StripSubdetector::TIB && id1.subdetId() <= StripSubdetector::TEC)) {
0278         const auto stripdet = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(hit1->geographicalId()));
0279         const StripTopology& Topo = stripdet->specificTopology();
0280         int Nstrips = Topo.nstrips();
0281         mypitch1 = stripdet->surface().bounds().width() / Topo.nstrips();
0282
0283         const auto det = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(mypointhit->geographicalId()));
0284
0285         TrajectoryStateOnSurface mytsos = itm->updatedState();
0286         LocalVector trackDirection = mytsos.localDirection();
0287         LocalVector drift = stripcpe.driftDirection(stripdet);
0288
0289         const auto hit1d = dynamic_cast<const SiStripRecHit1D*>(myhit);
0290
0291         if (hit1d) {
0292           getSimHitRes(det, trackDirection, *hit1d, expWidth, &mypitch1, drift);
0293           clusterWidth = hit1d->cluster()->amplitudes().size();
0294           uint16_t firstStrip = hit1d->cluster()->firstStrip();
0295           uint16_t lastStrip = firstStrip + (hit1d->cluster()->amplitudes()).size() - 1;
0296           atEdge = (firstStrip == 0 || lastStrip == (Nstrips - 1));
0297         }
0298
0299         const auto hit2d = dynamic_cast<const SiStripRecHit2D*>(myhit);
0300
0301         if (hit2d) {
0302           getSimHitRes(det, trackDirection, *hit2d, expWidth, &mypitch1, drift);
0303           clusterWidth = hit2d->cluster()->amplitudes().size();
0304           uint16_t firstStrip = hit2d->cluster()->firstStrip();
0305           uint16_t lastStrip = firstStrip + (hit2d->cluster()->amplitudes()).size() - 1;
0306           atEdge = (firstStrip == 0 || lastStrip == (Nstrips - 1));
0307         }
0308
0309         simpleRes =
0310             getSimpleRes(&(*itm));  // simple resolution by using the track re-fit forward and backward predicted state
0311
0312         histos2d_["residual_vs_trackMomentum"]->Fill(itm->updatedState().globalMomentum().mag(),
0313                                                      simpleRes * 10000);   // reso in cm *10000 == micro-meter
0314         histos2d_["residual_vs_trackPt"]->Fill(mymom, simpleRes * 10000);  // reso in cm *10000 == micro-meter
0315         histos2d_["residual_vs_trackEta"]->Fill(itm->updatedState().globalMomentum().eta(), simpleRes * 10000);
0316         histos2d_["residual_vs_trackPhi"]->Fill(itm->updatedState().globalMomentum().phi(), simpleRes * 10000);
0317         histos2d_["residual_vs_expectedWidth"]->Fill(expWidth, simpleRes * 10000);
0318         histos2d_["numHits_vs_residual"]->Fill(simpleRes * 10000, numHits);
0319
0320         // Now to see if there is a match - pair method - hit in overlapping sensors
0321         vector<TrajectoryMeasurement>::const_iterator itTraj2 = TMeas.end();  // last hit along the fitted track
0322
0323         for (auto itmCompare = itm - 1;
0324              // start to compare from the 5th hit
0325              itmCompare >= TMeas.cbegin() && itmCompare > itm - 4;
0326              --itmCompare) {
0327           const auto hit2 = itmCompare->recHit();
0328           if (!hit2->isValid())
0329             continue;
0330           DetId id2 = hit2->geographicalId();
0331
0332           //must be from the same detector and layer
0333           iidd1 = hit1->geographicalId().rawId();
0334           iidd2 = hit2->geographicalId().rawId();
0335           if (id1.subdetId() != id2.subdetId() || ::checkLayer(iidd1, tTopo) != ::checkLayer(iidd2, tTopo))
0336             break;
0337           //must both be stereo if one is
0338           if (tTopo->isStereo(id1) != tTopo->isStereo(id2))
0339             continue;
0340           //A check i dont completely understand but might as well keep there
0341           if (tTopo->glued(id1) == id1.rawId())
0342             LogDebug("HitResol") << "BAD GLUED: Have glued layer with id = " << id1.rawId()
0343                                  << " and glued id = " << tTopo->glued(id1) << "  and stereo = " << tTopo->isStereo(id1)
0344                                  << endl;
0345           if (tTopo->glued(id2) == id2.rawId())
0346             LogDebug("HitResol") << "BAD GLUED: Have glued layer with id = " << id2.rawId()
0347                                  << " and glued id = " << tTopo->glued(id2) << "  and stereo = " << tTopo->isStereo(id2)
0348                                  << endl;
0349
0350           itTraj2 = itmCompare;
0351           break;
0352         }
0353
0354         if (itTraj2 == TMeas.cend()) {
0355         } else {
0356           LogDebug("HitResol") << "Found overlapping sensors " << std::endl;
0357           pairsOnly = usePairsOnly_;
0358
0359           //We found one....let's fill in the truth info!
0360           TrajectoryStateOnSurface tsos_2 = itTraj2->updatedState();
0361           LocalVector trackDirection_2 = tsos_2.localDirection();
0362           const auto myhit2 = itTraj2->recHit();
0363           const auto myhit_2 = (*itTraj2->recHit()).hit();
0364           const auto stripdet_2 = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(myhit2->geographicalId()));
0365           const StripTopology& Topo_2 = stripdet_2->specificTopology();
0366           int Nstrips_2 = Topo_2.nstrips();
0367           float mypitch_2 = stripdet_2->surface().bounds().width() / Topo_2.nstrips();
0368
0369           if (mypitch1 != mypitch_2)
0370             return;  // for PairsOnly
0371
0372           const auto det_2 = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(myhit2->geographicalId()));
0373
0374           LocalVector drift_2 = stripcpe.driftDirection(stripdet_2);
0375
0376           const auto hit1d_2 = dynamic_cast<const SiStripRecHit1D*>(myhit_2);
0377           if (hit1d_2) {
0378             getSimHitRes(det_2, trackDirection_2, *hit1d_2, expWidth_2, &mypitch_2, drift_2);
0379             clusterWidth_2 = hit1d_2->cluster()->amplitudes().size();
0380             uint16_t firstStrip_2 = hit1d_2->cluster()->firstStrip();
0381             uint16_t lastStrip_2 = firstStrip_2 + (hit1d_2->cluster()->amplitudes()).size() - 1;
0382             atEdge_2 = (firstStrip_2 == 0 || lastStrip_2 == (Nstrips_2 - 1));
0383           }
0384
0385           const auto hit2d_2 = dynamic_cast<const SiStripRecHit2D*>(myhit_2);
0386           if (hit2d_2) {
0387             getSimHitRes(det_2, trackDirection_2, *hit2d_2, expWidth_2, &mypitch_2, drift_2);
0388             clusterWidth_2 = hit2d_2->cluster()->amplitudes().size();
0389             uint16_t firstStrip_2 = hit2d_2->cluster()->firstStrip();
0390             uint16_t lastStrip_2 = firstStrip_2 + (hit2d_2->cluster()->amplitudes()).size() - 1;
0391             atEdge_2 = (firstStrip_2 == 0 || lastStrip_2 == (Nstrips_2 - 1));
0392           }
0393
0394           // if(pairsOnly && (pitch != pitch2) ) fill = false;
0395
0396           // Make AnalyticalPropagator to use in getPairParameters
0397           AnalyticalPropagator mypropagator(magField_, anyDirection);
0398
0399           if (!getPairParameters(&(*magField_),
0400                                  mypropagator,
0401                                  &(*itTraj2),
0402                                  &(*itm),
0403                                  pairPath,
0404                                  hitDX,
0405                                  trackDX,
0406                                  trackDXE,
0407                                  trackParamX,
0408                                  trackParamY,
0409                                  trackParamDXDZ,
0410                                  trackParamDYDZ,
0411                                  trackParamXE,
0412                                  trackParamYE,
0413                                  trackParamDXDZE,
0414                                  trackParamDYDZE)) {
0415           } else {
0416             LogDebug("HitResol") << " \n\n\n"
0417                                  << " momentum       " << mymom << "\n"
0418                                  << " numHits        " << numHits << "\n"
0419                                  << " trackChi2      " << ProbTrackChi2 << "\n"
0420                                  << " detID1         " << iidd1 << "\n"
0421                                  << " pitch1         " << mypitch1 << "\n"
0422                                  << " clusterW1      " << clusterWidth << "\n"
0423                                  << " expectedW1     " << expWidth << "\n"
0424                                  << " atEdge1        " << atEdge << "\n"
0425                                  << " simpleRes      " << simpleRes << "\n"
0426                                  << " detID2         " << iidd2 << "\n"
0427                                  << " clusterW2      " << clusterWidth_2 << "\n"
0428                                  << " expectedW2     " << expWidth_2 << "\n"
0429                                  << " atEdge2        " << atEdge_2 << "\n"
0430                                  << " pairPath       " << pairPath << "\n"
0431                                  << " hitDX          " << hitDX << "\n"
0432                                  << " trackDX        " << trackDX << "\n"
0433                                  << " trackDXE       " << trackDXE << "\n"
0434                                  << " trackParamX      " << trackParamX << "\n"
0435                                  << " trackParamY      " << trackParamY << "\n"
0436                                  << " trackParamDXDZ " << trackParamDXDZ << "\n"
0437                                  << " trackParamDYDZ " << trackParamDYDZ << "\n"
0438                                  << " trackParamXE   " << trackParamXE << "\n"
0439                                  << " trackParamYE   " << trackParamYE << "\n"
0440                                  << " trackParamDXDZE" << trackParamDXDZE << "\n"
0441                                  << " trackParamDYDZE" << trackParamDYDZE << std::endl;
0442             reso->Fill();
0443           }
0444         }  //itTraj2 != TMeas.end()
0445       }  //hit1->isValid()....
0446     }  // itm
0447   }  // it
0448 }
0449
0450 void HitResol::endJob() {
0451   LogDebug("SiStripHitResolution:HitResol") << " Events Analysed             " << events << endl;
0452   LogDebug("SiStripHitResolution:HitResol") << " Number Of Tracked events    " << EventTrackCKF << endl;
0453
0454   reso->GetDirectory()->cd();
0455   reso->Write();
0456   treso->Write();
0457 }
0458
0459 double HitResol::checkConsistency(const StripClusterParameterEstimator::LocalValues& parameters,
0460                                   double xx,
0461                                   double xerr) {
0462   double error = sqrt(parameters.second.xx() + xerr * xerr);
0463   double separation = abs(parameters.first.x() - xx);
0464   double consistency = separation / error;
0465   return consistency;
0466 }
0467
0468 void HitResol::getSimHitRes(const GeomDetUnit* det,
0469                             const LocalVector& trackdirection,
0470                             const TrackingRecHit& recHit,
0471                             float& trackWidth,
0472                             float* pitch,
0473                             LocalVector& drift) {
0474   const auto stripdet = dynamic_cast<const StripGeomDetUnit*>(det);
0475   const auto& topol = dynamic_cast<const StripTopology&>(stripdet->topology());
0476
0477   LocalPoint position = recHit.localPosition();
0478   (*pitch) = topol.localPitch(position);
0479
0480   float anglealpha = 0;
0481   if (trackdirection.z() != 0) {
0482     anglealpha = atan(trackdirection.x() / trackdirection.z()) * TMath::RadToDeg();
0483   }
0484
0485   //  LocalVector drift = stripcpe.driftDirection(stripdet);
0486   float thickness = stripdet->surface().bounds().thickness();
0487   float tanalpha = tan(anglealpha * TMath::DegToRad());
0488   float tanalphaL = drift.x() / drift.z();
0489   (trackWidth) = fabs((thickness / (*pitch)) * tanalpha - (thickness / (*pitch)) * tanalphaL);
0490 }
0491
0492 double HitResol::getSimpleRes(const TrajectoryMeasurement* traj1) {
0493   TrajectoryStateOnSurface theCombinedPredictedState;
0494
0495   if (traj1->backwardPredictedState().isValid())
0496     theCombinedPredictedState =
0497         TrajectoryStateCombiner().combine(traj1->forwardPredictedState(), traj1->backwardPredictedState());
0498   else
0499     theCombinedPredictedState = traj1->forwardPredictedState();
0500
0501   if (!theCombinedPredictedState.isValid()) {
0502     return -100;
0503   }
0504
0505   const TransientTrackingRecHit::ConstRecHitPointer& firstRecHit = traj1->recHit();
0506   double recHitX_1 = firstRecHit->localPosition().x();
0507   return (theCombinedPredictedState.localPosition().x() - recHitX_1);
0508 }
0509
0510 //traj1 is the matched trajectory...traj2 is the original
0511 bool HitResol::getPairParameters(const MagneticField* magField_,
0512                                  AnalyticalPropagator& propagator,
0513                                  const TrajectoryMeasurement* traj1,
0514                                  const TrajectoryMeasurement* traj2,
0515                                  float& pairPath,
0516                                  float& hitDX,
0517                                  float& trackDX,
0518                                  float& trackDXE,
0519                                  float& trackParamX,
0520                                  float& trackParamY,
0521                                  float& trackParamDXDZ,
0522                                  float& trackParamDYDZ,
0523                                  float& trackParamXE,
0524                                  float& trackParamYE,
0525                                  float& trackParamDXDZE,
0526                                  float& trackParamDYDZE) {
0527   pairPath = 0;
0528   hitDX = 0;
0529   trackDX = 0;
0530   trackDXE = 0;
0531
0532   trackParamX = 0;
0533   trackParamY = 0;
0534   trackParamDXDZ = 0;
0535   trackParamDYDZ = 0;
0536   trackParamXE = 0;
0537   trackParamYE = 0;
0538   trackParamDXDZE = 0;
0539   trackParamDYDZE = 0;
0540
0541   TrajectoryStateCombiner combiner_;
0542
0543   // backward predicted state at module 1
0544   const TrajectoryStateOnSurface& bwdPred1 = traj1->backwardPredictedState();
0545   if (!bwdPred1.isValid())
0546     return false;
0547   LogDebug("HitResol") << "momentum from backward predicted state = " << bwdPred1.globalMomentum().mag() << endl;
0548   // forward predicted state at module 2
0549   const TrajectoryStateOnSurface& fwdPred2 = traj2->forwardPredictedState();
0550   LogDebug("HitResol") << "momentum from forward predicted state = " << fwdPred2.globalMomentum().mag() << endl;
0551   if (!fwdPred2.isValid())
0552     return false;
0553   // extrapolate fwdPred2 to module 1
0554   TrajectoryStateOnSurface fwdPred2At1 = propagator.propagate(fwdPred2, bwdPred1.surface());
0555   if (!fwdPred2At1.isValid())
0556     return false;
0557   // combine fwdPred2At1 with bwdPred1 (ref. state, best estimate without hits 1 and 2)
0558   TrajectoryStateOnSurface comb1 = combiner_.combine(bwdPred1, fwdPred2At1);
0559   if (!comb1.isValid())
0560     return false;
0561
0562   //
0563   // propagation of reference parameters to module 2
0564   //
0565   std::pair<TrajectoryStateOnSurface, double> tsosWithS = propagator.propagateWithPath(comb1, fwdPred2.surface());
0566   TrajectoryStateOnSurface comb1At2 = tsosWithS.first;
0567   if (!comb1At2.isValid())
0568     return false;
0569   //distance of propagation from one surface to the next==could cut here
0570   pairPath = tsosWithS.second;
0571   if (TMath::Abs(pairPath) > 15)
0572     return false;  //cut to remove hit pairs > 15 cm apart
0573
0574   // local parameters and errors on module 1
0575   AlgebraicVector5 pars = comb1.localParameters().vector();
0576   AlgebraicSymMatrix55 errs = comb1.localError().matrix();
0577   //number 3 is predX
0578   double predX1 = pars[3];
0579   //track fitted parameters in local coordinates for position 0
0580   (trackParamX) = pars[3];
0581   (trackParamY) = pars[4];
0582   (trackParamDXDZ) = pars[1];
0583   (trackParamDYDZ) = pars[2];
0584   (trackParamXE) = TMath::Sqrt(errs(3, 3));
0585   (trackParamYE) = TMath::Sqrt(errs(4, 4));
0586   (trackParamDXDZE) = TMath::Sqrt(errs(1, 1));
0587   (trackParamDYDZE) = TMath::Sqrt(errs(2, 2));
0588
0589   // local parameters and errors on module 2
0590   pars = comb1At2.localParameters().vector();
0591   errs = comb1At2.localError().matrix();
0592   double predX2 = pars[3];
0593
0594   ////
0595   //// jacobians (local-to-global@1,global 1-2,global-to-local@2)
0596   ////
0597   JacobianLocalToCurvilinear jacLocToCurv(comb1.surface(), comb1.localParameters(), *magField_);
0598   AnalyticalCurvilinearJacobian jacCurvToCurv(
0599       comb1.globalParameters(), comb1At2.globalPosition(), comb1At2.globalMomentum(), tsosWithS.second);
0600   JacobianCurvilinearToLocal jacCurvToLoc(comb1At2.surface(), comb1At2.localParameters(), *magField_);
0601   // combined jacobian local-1-to-local-2
0602   AlgebraicMatrix55 jacobian = jacLocToCurv.jacobian() * jacCurvToCurv.jacobian() * jacCurvToLoc.jacobian();
0603   // covariance on module 1
0604   AlgebraicSymMatrix55 covComb1 = comb1.localError().matrix();
0605   // variance and correlations for predicted local_x on modules 1 and 2
0606   double c00 = covComb1(3, 3);
0607   double c10(0.);
0608   double c11(0.);
0609   for (int i = 1; i < 5; ++i) {
0610     c10 += jacobian(3, i) * covComb1(i, 3);
0611     for (int j = 1; j < 5; ++j)
0612       c11 += jacobian(3, i) * covComb1(i, j) * jacobian(3, j);
0613   }
0614   // choose relative sign in order to minimize error on difference
0615   double diff = c00 - 2 * fabs(c10) + c11;
0616   diff = diff > 0 ? sqrt(diff) : -sqrt(-diff);
0617   (trackDXE) = diff;
0618   double relativeXSign_ = c10 > 0 ? -1 : 1;
0619
0620   (trackDX) = predX1 + relativeXSign_ * predX2;
0621
0622   double recHitX_1 = traj1->recHit()->localPosition().x();
0623   double recHitX_2 = traj2->recHit()->localPosition().x();
0624
0625   (hitDX) = recHitX_1 + relativeXSign_ * recHitX_2;
0626
0627   return true;
0628 }
0629
0630 void HitResol::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0631   edm::ParameterSetDescription desc;
0632   desc.add<edm::InputTag>("lumiScalers", edm::InputTag("scalersRawToDigi"));
0633   desc.add<edm::InputTag>("combinatorialTracks", edm::InputTag("generalTracks"));
0634   desc.add<edm::InputTag>("trajectories", edm::InputTag("generalTracks"));
0635   desc.addUntracked<int>("CompressionSettings", -1);
0636   desc.add<int>("Layer", 0);
0637   desc.add<bool>("Debug", false);
0638   desc.addUntracked<bool>("addLumi", false);
0639   desc.addUntracked<bool>("cutOnTracks", false);
0640   desc.addUntracked<unsigned int>("trackMultiplicity", 100);
0641   desc.addUntracked<double>("MomentumCut", 3.);
0642   desc.addUntracked<unsigned int>("UsePairsOnly", 1);
0643   descriptions.addWithDefaultLabel(desc);
0644 }
0645
0646 //define this as a plug-in
0647 DEFINE_FWK_MODULE(HitResol);