EgammaPhotonAlgos/src/ConversionSeedFinder.cc

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "FWCore/Framework/interface/ConsumesCollector.h"
0003 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0004 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionSeedFinder.h"
0005 // Field
0006 // Geometry
0007 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0008 //
0009 #include "RecoTracker/MeasurementDet/interface/LayerCollector.h"
0010 //
0011
0012 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0013 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0014 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0015
0016 ConversionSeedFinder::ConversionSeedFinder(const edm::ParameterSet& config, edm::ConsumesCollector& iC)
0017     :  //  conf_(config),
0018       theUpdator_() {
0019   measurementTrkToken_ = iC.consumes<MeasurementTrackerEvent>(
0020       edm::InputTag("MeasurementTrackerEvent"));  //hardcoded because the original was and no time to fix (sigh)
0021   beamSpotToken_ = iC.consumes<reco::BeamSpot>(
0022       edm::InputTag("offlineBeamSpot"));  //hardcoded because the original was and no time to fix (sigh)
0023
0024   theMFToken_ = iC.esConsumes();
0025   theGeomSearchTrackerToken_ = iC.esConsumes();
0026   thePropagatorAlongMomentumToken_ = iC.esConsumes(edm::ESInputTag("", "alongMomElePropagator"));
0027   thePropagatorOppositeToMomentumToken_ = iC.esConsumes(edm::ESInputTag("", "oppositeToMomElePropagator"));
0028   LogDebug("ConversionSeedFinder") << " CTOR "
0029                                    << "\n";
0030
0031   auto measurementTrackerName = config.getParameter<std::string>("MeasurementTrackerName");
0032   theMeasurementTrackerToken_ = iC.esConsumes(edm::ESInputTag("", measurementTrackerName));
0033 }
0034
0035 void ConversionSeedFinder::setEvent(const edm::Event& evt) {
0036   theTrackerGeom_ = this->getMeasurementTracker()->geomTracker();
0037
0038   //get the BeamSpot
0039   edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
0040   evt.getByToken(beamSpotToken_, recoBeamSpotHandle);
0041   theBeamSpot_ = *recoBeamSpotHandle;
0042
0043   evt.getByToken(measurementTrkToken_, theTrackerData_);
0044 }
0045
0046 void ConversionSeedFinder::setEventSetup(const edm::EventSetup& es) {
0047   theGeomSearchTracker_ = es.getHandle(theGeomSearchTrackerToken_);
0048   theMF_ = es.getHandle(theMFToken_);
0049
0050   theMeasurementTracker_ = &es.getData(theMeasurementTrackerToken_);
0051
0052   thePropagatorAlongMomentum_ = &es.getData(thePropagatorAlongMomentumToken_);
0053
0054   thePropagatorOppositeToMomentum_ = &es.getData(thePropagatorOppositeToMomentumToken_);
0055 }
0056
0057 void ConversionSeedFinder::findLayers() {
0058   int charge;
0059   //List the DetLayers crossed by a straight line from the centre of the
0060   //detector to the supercluster position
0061   //  GlobalPoint  vertex(0.,0.,0.);
0062   GlobalPoint vertex(theBeamSpot_.position().x(), theBeamSpot_.position().y(), theBeamSpot_.position().z());
0063   charge = -1;
0064   FreeTrajectoryState theStraightLineFTS = trackStateFromClusters(charge, vertex, alongMomentum, 1.);
0065
0066   findLayers(theStraightLineFTS);
0067 }
0068
0069 FreeTrajectoryState ConversionSeedFinder::trackStateFromClusters(int charge,
0070                                                                  const GlobalPoint& theOrigin,
0071                                                                  PropagationDirection dir,
0072                                                                  float scaleFactor) const {
0073   double caloEnergy = theSCenergy_ * scaleFactor;
0074
0075   GlobalVector radiusCalo = theSCPosition_ - theOrigin;
0076
0077   GlobalVector momentumWithoutCurvature = radiusCalo.unit() * caloEnergy;
0078
0079   GlobalTrajectoryParameters gtp;
0080   if (dir == alongMomentum) {
0081     gtp = GlobalTrajectoryParameters(theOrigin, momentumWithoutCurvature, charge, &(*theMF_));
0082   } else {
0083     gtp = GlobalTrajectoryParameters(theSCPosition_, momentumWithoutCurvature, charge, &(*theMF_));
0084   }
0085
0086   // now create error matrix
0087   // dpos = 4mm/sqrt(E), dtheta = move vertex by 1sigma
0088   float dpos = 0.4 / sqrt(theSCenergy_);
0089   dpos *= 2.;
0090   float dphi = dpos / theSCPosition_.perp();
0091   //  float dp = 0.03 * sqrt(theCaloEnergy);
0092   //  float dp = theCaloEnergy / sqrt(12.); // for fun
0093   float theta1 = theSCPosition_.theta();
0094   float theta2 = atan2(double(theSCPosition_.perp()), theSCPosition_.z() - 5.5);
0095   float dtheta = theta1 - theta2;
0096   AlgebraicSymMatrix55 m;
0097   m[0][0] = 1.;
0098   m[1][1] = dpos * dpos;
0099   m[2][2] = dpos * dpos;
0100   m[3][3] = dphi * dphi;
0101   m[4][4] = dtheta * dtheta;
0102
0103   FreeTrajectoryState fts(gtp, CurvilinearTrajectoryError(m));
0104
0105   return fts;
0106 }
0107
0108 void ConversionSeedFinder::findLayers(const FreeTrajectoryState& traj) {
0109   theLayerList_.clear();
0110
0111   StraightLinePropagator prop(&(*theMF_), alongMomentum);
0112
0113   LayerCollector collector(theNavigationSchool_, &prop, this->getMeasurementTracker(), 5., 5.);
0114
0115   theLayerList_ = collector.allLayers(traj);
0116
0117   for (unsigned int i = 0; i < theLayerList_.size(); ++i) {
0118     printLayer(i);
0119   }
0120 }
0121
0122 void ConversionSeedFinder::printLayer(int i) const {
0123   const DetLayer* layer = theLayerList_[i];
0124   if (layer->location() == GeomDetEnumerators::barrel) {
0125     //    const BarrelDetLayer * barrelLayer = dynamic_cast<const BarrelDetLayer*>(layer);
0126     //float r = barrelLayer->specificSurface().radius();
0127     //    std::cout   <<  " barrel layer radius " << r << " " << barrelLayer->specificSurface().bounds().length()/2. << "\n";
0128
0129   } else {
0130     //    const ForwardDetLayer * forwardLayer = dynamic_cast<const ForwardDetLayer*>(layer);
0131     // float z =  fabs(forwardLayer->surface().position().z());
0132     //    std::cout   << " forward layer position " << z << " " << forwardLayer->specificSurface().innerRadius() << " " << forwardLayer->specificSurface().outerRadius() << "\n";
0133   }
0134 }