/
githubmirror
/
cmssw
Обзор
Документация
Войти
/
githubmirror
/
cmssw
Код
Запросы
0
Пакеты
0
Релизы
0
Аналитика
Безопасность
master
CalibTracker/SiStripHitResolution/src/HitResol.cc
647 строк
30 KB
Shahzad Malik Muzaffar
[ALCA-DB] Drop Geometry/CommonDetUnit package
28 май 2026, 12:57
28 май 2026, 12:57
993d9d5
Код
Авторство
О чём код?
//////////////////////////////////////////////////////////////////////////////// // Package: CalibTracker/SiStripHitResolution // Class: HitResol // Original Authors: Denis Gele and Kathryn Coldham (adapted from HitEff) // modified by Khawla Jaffel for CPE studies // ported to cmssw by M. Musich // /////////////////////////////////////////////////////////////////////////////// // system include files #include <memory> #include <string> #include <iostream> // user includes #include "CalibFormats/SiStripObjects/interface/SiStripDetCabling.h" #include "CalibFormats/SiStripObjects/interface/SiStripQuality.h" #include "CalibTracker/Records/interface/SiStripDetCablingRcd.h" #include "CalibTracker/Records/interface/SiStripQualityRcd.h" #include "CalibTracker/SiStripHitEfficiency/interface/SiStripHitEfficiencyHelpers.h" #include "CalibTracker/SiStripHitResolution/interface/HitResol.h" #include "CommonTools/Statistics/interface/ChiSquaredProbability.h" #include "CommonTools/UtilAlgos/interface/TFileService.h" #include "DataFormats/Common/interface/Handle.h" #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h" #include "DataFormats/GeometryVector/interface/GlobalPoint.h" #include "DataFormats/GeometryVector/interface/GlobalVector.h" #include "DataFormats/GeometryVector/interface/LocalVector.h" #include "DataFormats/SiStripCluster/interface/SiStripCluster.h" #include "DataFormats/SiStripDetId/interface/StripSubdetector.h" #include "DataFormats/TrackReco/interface/Track.h" #include "DataFormats/TrackReco/interface/TrackBase.h" #include "DataFormats/TrackReco/interface/TrackExtra.h" #include "DataFormats/TrackReco/interface/TrackFwd.h" #include "DataFormats/TrackerCommon/interface/TrackerTopology.h" #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h" #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h" #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h" #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h" #include "FWCore/Framework/interface/ESHandle.h" #include "FWCore/Framework/interface/Event.h" #include "FWCore/Framework/interface/EventSetup.h" #include "FWCore/Framework/interface/Frameworkfwd.h" #include "FWCore/Framework/interface/MakerMacros.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" #include "FWCore/ServiceRegistry/interface/Service.h" #include "Geometry/CommonTopologies/interface/GeomDet.h" #include "Geometry/CommonTopologies/interface/GeomDetType.h" #include "Geometry/CommonTopologies/interface/GluedGeomDet.h" #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" #include "Geometry/Records/interface/TrackerTopologyRcd.h" #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" #include "RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h" #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h" #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h" #include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h" #include "TrackingTools/AnalyticalJacobians/interface/JacobianCurvilinearToLocal.h" #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h" #include "TrackingTools/Records/interface/TransientRecHitRecord.h" #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h" // ROOT includes #include "TMath.h" #include "TH1F.h" // // constructors and destructor // using namespace std; HitResol::HitResol(const edm::ParameterSet& conf) : scalerToken_(consumes<LumiScalersCollection>(conf.getParameter<edm::InputTag>("lumiScalers"))), combinatorialTracks_token_( consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("combinatorialTracks"))), tjToken_(consumes<std::vector<Trajectory> >(conf.getParameter<edm::InputTag>("trajectories"))), topoToken_(esConsumes()), geomToken_(esConsumes()), cpeToken_(esConsumes(edm::ESInputTag("", "StripCPEfromTrackAngle"))), siStripQualityToken_(esConsumes()), magFieldToken_(esConsumes()), addLumi_(conf.getUntrackedParameter<bool>("addLumi", false)), DEBUG_(conf.getParameter<bool>("Debug")), cutOnTracks_(conf.getUntrackedParameter<bool>("cutOnTracks", false)), momentumCut_(conf.getUntrackedParameter<double>("MomentumCut", 3.)), compSettings_(conf.getUntrackedParameter<int>("CompressionSettings", -1)), usePairsOnly_(conf.getUntrackedParameter<unsigned int>("UsePairsOnly", 1)), layers_(conf.getParameter<int>("Layer")), trackMultiplicityCut_(conf.getUntrackedParameter<unsigned int>("trackMultiplicity", 100)) { usesResource(TFileService::kSharedResource); } void HitResol::beginJob() { edm::Service<TFileService> fs; if (compSettings_ > 0) { edm::LogInfo("SiStripHitResolution:HitResol") << "the compressions settings are:" << compSettings_ << std::endl; fs->file().SetCompressionSettings(compSettings_); } reso = fs->make<TTree>("reso", "tree hit pairs for resolution studies"); reso->Branch("momentum", &mymom, "momentum/F"); reso->Branch("numHits", &numHits, "numHits/I"); reso->Branch("trackChi2", &ProbTrackChi2, "trackChi2/F"); reso->Branch("detID1", &iidd1, "detID1/I"); reso->Branch("pitch1", &mypitch1, "pitch1/F"); reso->Branch("clusterW1", &clusterWidth, "clusterW1/I"); reso->Branch("expectedW1", &expWidth, "expectedW1/F"); reso->Branch("atEdge1", &atEdge, "atEdge1/F"); reso->Branch("simpleRes", &simpleRes, "simpleRes/F"); reso->Branch("detID2", &iidd2, "detID2/I"); reso->Branch("clusterW2", &clusterWidth_2, "clusterW2/I"); reso->Branch("expectedW2", &expWidth_2, "expectedW2/F"); reso->Branch("atEdge2", &atEdge_2, "atEdge2/F"); reso->Branch("pairPath", &pairPath, "pairPath/F"); reso->Branch("hitDX", &hitDX, "hitDX/F"); reso->Branch("trackDX", &trackDX, "trackDX/F"); reso->Branch("trackDXE", &trackDXE, "trackDXE/F"); reso->Branch("trackParamX", &trackParamX, "trackParamX/F"); reso->Branch("trackParamY", &trackParamY, "trackParamY/F"); reso->Branch("trackParamDXDZ", &trackParamDXDZ, "trackParamDXDZ/F"); reso->Branch("trackParamDYDZ", &trackParamDYDZ, "trackParamDYDZ/F"); reso->Branch("trackParamXE", &trackParamXE, "trackParamXE/F"); reso->Branch("trackParamYE", &trackParamYE, "trackParamYE/F"); reso->Branch("trackParamDXDZE", &trackParamDXDZE, "trackParamDXDZE/F"); reso->Branch("trackParamDYDZE", &trackParamDYDZE, "trackParamDYDZE/F"); reso->Branch("pairsOnly", &pairsOnly, "pairsOnly/I"); treso = fs->make<TTree>("treso", "tree tracks for resolution studies"); treso->Branch("track_momentum", &track_momentum, "track_momentum/F"); treso->Branch("track_pt", &track_pt, "track_pt/F"); treso->Branch("track_eta", &track_eta, "track_eta/F"); treso->Branch("track_phi", &track_phi, "track_phi/F"); treso->Branch("track_trackChi2", &track_trackChi2, "track_trackChi2/F"); treso->Branch("track_width", &expWidth, "track_width/F"); // from 1D HIT treso->Branch("NumberOf_tracks", &NumberOf_tracks, "NumberOf_tracks/I"); events = 0; EventTrackCKF = 0; histos2d_["track_phi_vs_eta"] = new TH2F("track_phi_vs_eta", ";track phi;track eta", 60, -3.5, 3.5, 60, -3., 3.); histos2d_["residual_vs_trackMomentum"] = new TH2F("residual_vs_trackMomentum", ";track momentum [GeV]; x_{pred_track} - x_{reco_hit} [#mum]", 60, 0., 10., 60, 0., 200.); histos2d_["residual_vs_trackPt"] = new TH2F( "residual_vs_trackPt", ";track p_{T}[GeV];x_{pred_track} - x_{reco_hit} [#mum]", 60, 0., 10., 60, 0., 200.); histos2d_["residual_vs_trackEta"] = new TH2F("residual_vs_trackEta", ";track #eta;x_{pred_track} - x_{reco_hit} [#mum]", 60, 0., 3., 60, 0., 200.); histos2d_["residual_vs_trackPhi"] = new TH2F("residual_vs_trackPhi", ";track #phi;x_{pred_track} - x_{reco_hit} [#mum]", 60, 0., 3.5, 60, 0., 200.); histos2d_["residual_vs_expectedWidth"] = new TH2F( "residual_vs_expectedWidth", ";track Width;x_{pred_track} - x_{reco_hit} [#mum]", 3, 0., 3., 60, 0., 200.); histos2d_["numHits_vs_residual"] = new TH2F("numHits_vs_residual", ";x_{pred_track} - x_{reco_hit} [#mum];N Hits", 60, 0., 200., 15, 0., 15.); } void HitResol::analyze(const edm::Event& e, const edm::EventSetup& es) { //Retrieve tracker topology from geometry const TrackerTopology* tTopo = &es.getData(topoToken_); LogDebug("SiStripHitResolution:HitResol") << "beginning analyze from HitResol" << endl; using namespace edm; using namespace reco; // Step A: Get Inputs int run_nr = e.id().run(); int ev_nr = e.id().event(); // get the tracks edm::Handle<reco::TrackCollection> trackCollectionCKF; e.getByToken(combinatorialTracks_token_, trackCollectionCKF); const reco::TrackCollection* tracksCKF = trackCollectionCKF.product(); // get the trajectory collection edm::Handle<std::vector<Trajectory> > trajectoryCollectionHandle; e.getByToken(tjToken_, trajectoryCollectionHandle); const TrajectoryCollection* trajectoryCollection = trajectoryCollectionHandle.product(); //get tracker geometry edm::ESHandle<TrackerGeometry> tracker = es.getHandle(geomToken_); const TrackerGeometry* tkgeom = &(*tracker); //get Cluster Parameter Estimator edm::ESHandle<StripClusterParameterEstimator> parameterestimator = es.getHandle(cpeToken_); const StripClusterParameterEstimator& stripcpe(*parameterestimator); // get the SiStripQuality records edm::ESHandle<SiStripQuality> SiStripQuality_ = es.getHandle(siStripQualityToken_); // get the magnetic field const MagneticField* magField_ = &es.getData(magFieldToken_); events++; // List of variables for SiStripHitResolution ntuple mymom = 0; numHits = 0; ProbTrackChi2 = 0; iidd1 = 0; mypitch1 = 0; clusterWidth = 0; expWidth = 0; atEdge = 0; simpleRes = 0; iidd2 = 0; clusterWidth_2 = 0; expWidth_2 = 0; atEdge_2 = 0; pairPath = 0; hitDX = 0; trackDX = 0; trackDXE = 0; trackParamX = 0; trackParamY = 0; trackParamDXDZ = 0; trackParamDYDZ = 0; trackParamXE = 0; trackParamYE = 0; trackParamDXDZE = 0; trackParamDYDZE = 0; pairsOnly = 0; LogDebug("HitResol") << "Starting analysis, nrun nevent, tracksCKF->size(): " << run_nr << " " << ev_nr << " " << tracksCKF->size() << std::endl; for (unsigned int iT = 0; iT < tracksCKF->size(); ++iT) { track_momentum = tracksCKF->at(iT).p(); track_pt = tracksCKF->at(iT).p(); track_eta = tracksCKF->at(iT).eta(); track_phi = tracksCKF->at(iT).phi(); track_trackChi2 = ChiSquaredProbability((double)(tracksCKF->at(iT).chi2()), (double)(tracksCKF->at(iT).ndof())); treso->Fill(); } histos2d_["track_phi_vs_eta"]->Fill(track_phi, track_eta); // loop over trajectories from refit for (const auto& traj : *trajectoryCollection) { const auto& TMeas = traj.measurements(); // Loop on each measurement and take it into consideration //-------------------------------------------------------- for (auto itm = TMeas.cbegin(); itm != TMeas.cend(); ++itm) { if (!itm->updatedState().isValid()) { LogDebug("HitResol") << "trajectory measurement not valid" << std::endl; continue; } const TransientTrackingRecHit::ConstRecHitPointer mypointhit = itm->recHit(); const TrackingRecHit* myhit = (*itm->recHit()).hit(); ProbTrackChi2 = 0; numHits = 0; LogDebug("HitResol") << "TrackChi2 = " << ChiSquaredProbability((double)(traj.chiSquared()), (double)(traj.ndof(false))) << "\n" << "itm->updatedState().globalMomentum().perp(): " << itm->updatedState().globalMomentum().perp() << "\n" << "numhits " << traj.foundHits() << std::endl; numHits = traj.foundHits(); ProbTrackChi2 = ChiSquaredProbability((double)(traj.chiSquared()), (double)(traj.ndof(false))); mymom = itm->updatedState().globalMomentum().perp(); //Now for the first hit TrajectoryStateOnSurface mytsos = itm->updatedState(); const auto hit1 = itm->recHit(); DetId id1 = hit1->geographicalId(); if (id1.subdetId() < StripSubdetector::TIB || id1.subdetId() > StripSubdetector::TEC) continue; if (hit1->isValid() && mymom > momentumCut_ && (id1.subdetId() >= StripSubdetector::TIB && id1.subdetId() <= StripSubdetector::TEC)) { const auto stripdet = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(hit1->geographicalId())); const StripTopology& Topo = stripdet->specificTopology(); int Nstrips = Topo.nstrips(); mypitch1 = stripdet->surface().bounds().width() / Topo.nstrips(); const auto det = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(mypointhit->geographicalId())); TrajectoryStateOnSurface mytsos = itm->updatedState(); LocalVector trackDirection = mytsos.localDirection(); LocalVector drift = stripcpe.driftDirection(stripdet); const auto hit1d = dynamic_cast<const SiStripRecHit1D*>(myhit); if (hit1d) { getSimHitRes(det, trackDirection, *hit1d, expWidth, &mypitch1, drift); clusterWidth = hit1d->cluster()->amplitudes().size(); uint16_t firstStrip = hit1d->cluster()->firstStrip(); uint16_t lastStrip = firstStrip + (hit1d->cluster()->amplitudes()).size() - 1; atEdge = (firstStrip == 0 || lastStrip == (Nstrips - 1)); } const auto hit2d = dynamic_cast<const SiStripRecHit2D*>(myhit); if (hit2d) { getSimHitRes(det, trackDirection, *hit2d, expWidth, &mypitch1, drift); clusterWidth = hit2d->cluster()->amplitudes().size(); uint16_t firstStrip = hit2d->cluster()->firstStrip(); uint16_t lastStrip = firstStrip + (hit2d->cluster()->amplitudes()).size() - 1; atEdge = (firstStrip == 0 || lastStrip == (Nstrips - 1)); } simpleRes = getSimpleRes(&(*itm)); // simple resolution by using the track re-fit forward and backward predicted state histos2d_["residual_vs_trackMomentum"]->Fill(itm->updatedState().globalMomentum().mag(), simpleRes * 10000); // reso in cm *10000 == micro-meter histos2d_["residual_vs_trackPt"]->Fill(mymom, simpleRes * 10000); // reso in cm *10000 == micro-meter histos2d_["residual_vs_trackEta"]->Fill(itm->updatedState().globalMomentum().eta(), simpleRes * 10000); histos2d_["residual_vs_trackPhi"]->Fill(itm->updatedState().globalMomentum().phi(), simpleRes * 10000); histos2d_["residual_vs_expectedWidth"]->Fill(expWidth, simpleRes * 10000); histos2d_["numHits_vs_residual"]->Fill(simpleRes * 10000, numHits); // Now to see if there is a match - pair method - hit in overlapping sensors vector<TrajectoryMeasurement>::const_iterator itTraj2 = TMeas.end(); // last hit along the fitted track for (auto itmCompare = itm - 1; // start to compare from the 5th hit itmCompare >= TMeas.cbegin() && itmCompare > itm - 4; --itmCompare) { const auto hit2 = itmCompare->recHit(); if (!hit2->isValid()) continue; DetId id2 = hit2->geographicalId(); //must be from the same detector and layer iidd1 = hit1->geographicalId().rawId(); iidd2 = hit2->geographicalId().rawId(); if (id1.subdetId() != id2.subdetId() || ::checkLayer(iidd1, tTopo) != ::checkLayer(iidd2, tTopo)) break; //must both be stereo if one is if (tTopo->isStereo(id1) != tTopo->isStereo(id2)) continue; //A check i dont completely understand but might as well keep there if (tTopo->glued(id1) == id1.rawId()) LogDebug("HitResol") << "BAD GLUED: Have glued layer with id = " << id1.rawId() << " and glued id = " << tTopo->glued(id1) << " and stereo = " << tTopo->isStereo(id1) << endl; if (tTopo->glued(id2) == id2.rawId()) LogDebug("HitResol") << "BAD GLUED: Have glued layer with id = " << id2.rawId() << " and glued id = " << tTopo->glued(id2) << " and stereo = " << tTopo->isStereo(id2) << endl; itTraj2 = itmCompare; break; } if (itTraj2 == TMeas.cend()) { } else { LogDebug("HitResol") << "Found overlapping sensors " << std::endl; pairsOnly = usePairsOnly_; //We found one....let's fill in the truth info! TrajectoryStateOnSurface tsos_2 = itTraj2->updatedState(); LocalVector trackDirection_2 = tsos_2.localDirection(); const auto myhit2 = itTraj2->recHit(); const auto myhit_2 = (*itTraj2->recHit()).hit(); const auto stripdet_2 = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(myhit2->geographicalId())); const StripTopology& Topo_2 = stripdet_2->specificTopology(); int Nstrips_2 = Topo_2.nstrips(); float mypitch_2 = stripdet_2->surface().bounds().width() / Topo_2.nstrips(); if (mypitch1 != mypitch_2) return; // for PairsOnly const auto det_2 = dynamic_cast<const StripGeomDetUnit*>(tkgeom->idToDetUnit(myhit2->geographicalId())); LocalVector drift_2 = stripcpe.driftDirection(stripdet_2); const auto hit1d_2 = dynamic_cast<const SiStripRecHit1D*>(myhit_2); if (hit1d_2) { getSimHitRes(det_2, trackDirection_2, *hit1d_2, expWidth_2, &mypitch_2, drift_2); clusterWidth_2 = hit1d_2->cluster()->amplitudes().size(); uint16_t firstStrip_2 = hit1d_2->cluster()->firstStrip(); uint16_t lastStrip_2 = firstStrip_2 + (hit1d_2->cluster()->amplitudes()).size() - 1; atEdge_2 = (firstStrip_2 == 0 || lastStrip_2 == (Nstrips_2 - 1)); } const auto hit2d_2 = dynamic_cast<const SiStripRecHit2D*>(myhit_2); if (hit2d_2) { getSimHitRes(det_2, trackDirection_2, *hit2d_2, expWidth_2, &mypitch_2, drift_2); clusterWidth_2 = hit2d_2->cluster()->amplitudes().size(); uint16_t firstStrip_2 = hit2d_2->cluster()->firstStrip(); uint16_t lastStrip_2 = firstStrip_2 + (hit2d_2->cluster()->amplitudes()).size() - 1; atEdge_2 = (firstStrip_2 == 0 || lastStrip_2 == (Nstrips_2 - 1)); } // if(pairsOnly && (pitch != pitch2) ) fill = false; // Make AnalyticalPropagator to use in getPairParameters AnalyticalPropagator mypropagator(magField_, anyDirection); if (!getPairParameters(&(*magField_), mypropagator, &(*itTraj2), &(*itm), pairPath, hitDX, trackDX, trackDXE, trackParamX, trackParamY, trackParamDXDZ, trackParamDYDZ, trackParamXE, trackParamYE, trackParamDXDZE, trackParamDYDZE)) { } else { LogDebug("HitResol") << " \n\n\n" << " momentum " << mymom << "\n" << " numHits " << numHits << "\n" << " trackChi2 " << ProbTrackChi2 << "\n" << " detID1 " << iidd1 << "\n" << " pitch1 " << mypitch1 << "\n" << " clusterW1 " << clusterWidth << "\n" << " expectedW1 " << expWidth << "\n" << " atEdge1 " << atEdge << "\n" << " simpleRes " << simpleRes << "\n" << " detID2 " << iidd2 << "\n" << " clusterW2 " << clusterWidth_2 << "\n" << " expectedW2 " << expWidth_2 << "\n" << " atEdge2 " << atEdge_2 << "\n" << " pairPath " << pairPath << "\n" << " hitDX " << hitDX << "\n" << " trackDX " << trackDX << "\n" << " trackDXE " << trackDXE << "\n" << " trackParamX " << trackParamX << "\n" << " trackParamY " << trackParamY << "\n" << " trackParamDXDZ " << trackParamDXDZ << "\n" << " trackParamDYDZ " << trackParamDYDZ << "\n" << " trackParamXE " << trackParamXE << "\n" << " trackParamYE " << trackParamYE << "\n" << " trackParamDXDZE" << trackParamDXDZE << "\n" << " trackParamDYDZE" << trackParamDYDZE << std::endl; reso->Fill(); } } //itTraj2 != TMeas.end() } //hit1->isValid().... } // itm } // it } void HitResol::endJob() { LogDebug("SiStripHitResolution:HitResol") << " Events Analysed " << events << endl; LogDebug("SiStripHitResolution:HitResol") << " Number Of Tracked events " << EventTrackCKF << endl; reso->GetDirectory()->cd(); reso->Write(); treso->Write(); } double HitResol::checkConsistency(const StripClusterParameterEstimator::LocalValues& parameters, double xx, double xerr) { double error = sqrt(parameters.second.xx() + xerr * xerr); double separation = abs(parameters.first.x() - xx); double consistency = separation / error; return consistency; } void HitResol::getSimHitRes(const GeomDetUnit* det, const LocalVector& trackdirection, const TrackingRecHit& recHit, float& trackWidth, float* pitch, LocalVector& drift) { const auto stripdet = dynamic_cast<const StripGeomDetUnit*>(det); const auto& topol = dynamic_cast<const StripTopology&>(stripdet->topology()); LocalPoint position = recHit.localPosition(); (*pitch) = topol.localPitch(position); float anglealpha = 0; if (trackdirection.z() != 0) { anglealpha = atan(trackdirection.x() / trackdirection.z()) * TMath::RadToDeg(); } // LocalVector drift = stripcpe.driftDirection(stripdet); float thickness = stripdet->surface().bounds().thickness(); float tanalpha = tan(anglealpha * TMath::DegToRad()); float tanalphaL = drift.x() / drift.z(); (trackWidth) = fabs((thickness / (*pitch)) * tanalpha - (thickness / (*pitch)) * tanalphaL); } double HitResol::getSimpleRes(const TrajectoryMeasurement* traj1) { TrajectoryStateOnSurface theCombinedPredictedState; if (traj1->backwardPredictedState().isValid()) theCombinedPredictedState = TrajectoryStateCombiner().combine(traj1->forwardPredictedState(), traj1->backwardPredictedState()); else theCombinedPredictedState = traj1->forwardPredictedState(); if (!theCombinedPredictedState.isValid()) { return -100; } const TransientTrackingRecHit::ConstRecHitPointer& firstRecHit = traj1->recHit(); double recHitX_1 = firstRecHit->localPosition().x(); return (theCombinedPredictedState.localPosition().x() - recHitX_1); } //traj1 is the matched trajectory...traj2 is the original bool HitResol::getPairParameters(const MagneticField* magField_, AnalyticalPropagator& propagator, const TrajectoryMeasurement* traj1, const TrajectoryMeasurement* traj2, float& pairPath, float& hitDX, float& trackDX, float& trackDXE, float& trackParamX, float& trackParamY, float& trackParamDXDZ, float& trackParamDYDZ, float& trackParamXE, float& trackParamYE, float& trackParamDXDZE, float& trackParamDYDZE) { pairPath = 0; hitDX = 0; trackDX = 0; trackDXE = 0; trackParamX = 0; trackParamY = 0; trackParamDXDZ = 0; trackParamDYDZ = 0; trackParamXE = 0; trackParamYE = 0; trackParamDXDZE = 0; trackParamDYDZE = 0; TrajectoryStateCombiner combiner_; // backward predicted state at module 1 const TrajectoryStateOnSurface& bwdPred1 = traj1->backwardPredictedState(); if (!bwdPred1.isValid()) return false; LogDebug("HitResol") << "momentum from backward predicted state = " << bwdPred1.globalMomentum().mag() << endl; // forward predicted state at module 2 const TrajectoryStateOnSurface& fwdPred2 = traj2->forwardPredictedState(); LogDebug("HitResol") << "momentum from forward predicted state = " << fwdPred2.globalMomentum().mag() << endl; if (!fwdPred2.isValid()) return false; // extrapolate fwdPred2 to module 1 TrajectoryStateOnSurface fwdPred2At1 = propagator.propagate(fwdPred2, bwdPred1.surface()); if (!fwdPred2At1.isValid()) return false; // combine fwdPred2At1 with bwdPred1 (ref. state, best estimate without hits 1 and 2) TrajectoryStateOnSurface comb1 = combiner_.combine(bwdPred1, fwdPred2At1); if (!comb1.isValid()) return false; // // propagation of reference parameters to module 2 // std::pair<TrajectoryStateOnSurface, double> tsosWithS = propagator.propagateWithPath(comb1, fwdPred2.surface()); TrajectoryStateOnSurface comb1At2 = tsosWithS.first; if (!comb1At2.isValid()) return false; //distance of propagation from one surface to the next==could cut here pairPath = tsosWithS.second; if (TMath::Abs(pairPath) > 15) return false; //cut to remove hit pairs > 15 cm apart // local parameters and errors on module 1 AlgebraicVector5 pars = comb1.localParameters().vector(); AlgebraicSymMatrix55 errs = comb1.localError().matrix(); //number 3 is predX double predX1 = pars[3]; //track fitted parameters in local coordinates for position 0 (trackParamX) = pars[3]; (trackParamY) = pars[4]; (trackParamDXDZ) = pars[1]; (trackParamDYDZ) = pars[2]; (trackParamXE) = TMath::Sqrt(errs(3, 3)); (trackParamYE) = TMath::Sqrt(errs(4, 4)); (trackParamDXDZE) = TMath::Sqrt(errs(1, 1)); (trackParamDYDZE) = TMath::Sqrt(errs(2, 2)); // local parameters and errors on module 2 pars = comb1At2.localParameters().vector(); errs = comb1At2.localError().matrix(); double predX2 = pars[3]; //// //// jacobians (local-to-global@1,global 1-2,global-to-local@2) //// JacobianLocalToCurvilinear jacLocToCurv(comb1.surface(), comb1.localParameters(), *magField_); AnalyticalCurvilinearJacobian jacCurvToCurv( comb1.globalParameters(), comb1At2.globalPosition(), comb1At2.globalMomentum(), tsosWithS.second); JacobianCurvilinearToLocal jacCurvToLoc(comb1At2.surface(), comb1At2.localParameters(), *magField_); // combined jacobian local-1-to-local-2 AlgebraicMatrix55 jacobian = jacLocToCurv.jacobian() * jacCurvToCurv.jacobian() * jacCurvToLoc.jacobian(); // covariance on module 1 AlgebraicSymMatrix55 covComb1 = comb1.localError().matrix(); // variance and correlations for predicted local_x on modules 1 and 2 double c00 = covComb1(3, 3); double c10(0.); double c11(0.); for (int i = 1; i < 5; ++i) { c10 += jacobian(3, i) * covComb1(i, 3); for (int j = 1; j < 5; ++j) c11 += jacobian(3, i) * covComb1(i, j) * jacobian(3, j); } // choose relative sign in order to minimize error on difference double diff = c00 - 2 * fabs(c10) + c11; diff = diff > 0 ? sqrt(diff) : -sqrt(-diff); (trackDXE) = diff; double relativeXSign_ = c10 > 0 ? -1 : 1; (trackDX) = predX1 + relativeXSign_ * predX2; double recHitX_1 = traj1->recHit()->localPosition().x(); double recHitX_2 = traj2->recHit()->localPosition().x(); (hitDX) = recHitX_1 + relativeXSign_ * recHitX_2; return true; } void HitResol::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; desc.add<edm::InputTag>("lumiScalers", edm::InputTag("scalersRawToDigi")); desc.add<edm::InputTag>("combinatorialTracks", edm::InputTag("generalTracks")); desc.add<edm::InputTag>("trajectories", edm::InputTag("generalTracks")); desc.addUntracked<int>("CompressionSettings", -1); desc.add<int>("Layer", 0); desc.add<bool>("Debug", false); desc.addUntracked<bool>("addLumi", false); desc.addUntracked<bool>("cutOnTracks", false); desc.addUntracked<unsigned int>("trackMultiplicity", 100); desc.addUntracked<double>("MomentumCut", 3.); desc.addUntracked<unsigned int>("UsePairsOnly", 1); descriptions.addWithDefaultLabel(desc); } //define this as a plug-in DEFINE_FWK_MODULE(HitResol);