ConversionSeedGenerators/plugins/SeedForPhotonConversion1Leg.cc

0001 #include "SeedForPhotonConversion1Leg.h"
0002
0003 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0004 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
0007 #include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
0008 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0009 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0010
0011 //#define mydebug_seed
0012 namespace {
0013   template <class T>
0014   T sqr(T t) {
0015     return t * t;
0016   }
0017 }  // namespace
0018
0019 const TrajectorySeed* SeedForPhotonConversion1Leg::trajectorySeed(TrajectorySeedCollection& seedCollection,
0020                                                                   const SeedingHitSet& hits,
0021                                                                   const GlobalPoint& vertex,
0022                                                                   const GlobalVector& vertexBounds,
0023                                                                   float ptmin,
0024                                                                   const edm::EventSetup& es,
0025                                                                   float cotTheta,
0026                                                                   std::stringstream& ss) {
0027   pss = &ss;
0028   if (hits.size() < 2)
0029     return nullptr;
0030
0031   GlobalTrajectoryParameters kine = initialKinematic(hits, vertex, es, cotTheta);
0032   float sinTheta = sin(kine.momentum().theta());
0033
0034   CurvilinearTrajectoryError error = initialError(vertexBounds, ptmin, sinTheta);
0035   FreeTrajectoryState fts(kine, error);
0036
0037   return buildSeed(seedCollection, hits, fts, es);
0038 }
0039
0040 GlobalTrajectoryParameters SeedForPhotonConversion1Leg::initialKinematic(const SeedingHitSet& hits,
0041                                                                          const GlobalPoint& vertexPos,
0042                                                                          const edm::EventSetup& es,
0043                                                                          const float cotTheta) const {
0044   GlobalTrajectoryParameters kine;
0045
0046   SeedingHitSet::ConstRecHitPointer tth1 = hits[0];
0047   SeedingHitSet::ConstRecHitPointer tth2 = hits[1];
0048
0049   // FIXME optimize: move outside loop
0050   const auto& bfield = es.getData(theBfieldToken);
0051   float nomField = bfield.nominalValue();
0052
0053   FastHelix helix(tth2->globalPosition(), tth1->globalPosition(), vertexPos, nomField, &bfield, vertexPos);
0054   kine = helix.stateAtVertex();
0055
0056   //force the pz/pt equal to the measured one
0057   if (fabs(cotTheta) < cotTheta_Max)
0058     kine = GlobalTrajectoryParameters(
0059         kine.position(),
0060         GlobalVector(kine.momentum().x(), kine.momentum().y(), kine.momentum().perp() * cotTheta),
0061         kine.charge(),
0062         &kine.magneticField());
0063   else
0064     kine = GlobalTrajectoryParameters(
0065         kine.position(),
0066         GlobalVector(kine.momentum().x(), kine.momentum().y(), kine.momentum().perp() * cotTheta_Max),
0067         kine.charge(),
0068         &kine.magneticField());
0069
0070 #ifdef mydebug_seed
0071   uint32_t detid;
0072   (*pss) << "[SeedForPhotonConversion1Leg] initialKinematic tth1 ";
0073   detid = tth1->geographicalId().rawId();
0074   po.print(*pss, detid);
0075   (*pss) << " \t " << detid << " " << tth1->localPosition() << " " << tth1->globalPosition();
0076   detid = tth2->geographicalId().rawId();
0077   (*pss) << " \n\t tth2 ";
0078   po.print(*pss, detid);
0079   (*pss) << " \t " << detid << " " << tth2->localPosition() << " " << tth2->globalPosition() << "\nhelix momentum "
0080          << kine.momentum() << " pt " << kine.momentum().perp() << " radius " << 1 / kine.transverseCurvature();
0081 #endif
0082
0083   bool isBOFF = (0 == nomField);
0084   ;
0085   if (isBOFF && (theBOFFMomentum > 0)) {
0086     kine =
0087         GlobalTrajectoryParameters(kine.position(), kine.momentum().unit() * theBOFFMomentum, kine.charge(), &bfield);
0088   }
0089   return kine;
0090 }
0091
0092 CurvilinearTrajectoryError SeedForPhotonConversion1Leg::initialError(const GlobalVector& vertexBounds,
0093                                                                      float ptMin,
0094                                                                      float sinTheta) const {
0095   // Set initial uncertainty on track parameters, using only P.V. constraint and no hit
0096   // information.
0097   GlobalError vertexErr(sqr(vertexBounds.x()), 0, sqr(vertexBounds.y()), 0, 0, sqr(vertexBounds.z()));
0098
0099   AlgebraicSymMatrix55 C = ROOT::Math::SMatrixIdentity();
0100
0101   // FIXME: minC00. Prevent apriori uncertainty in 1/P from being too small,
0102   // to avoid instabilities.
0103   // N.B. This parameter needs optimising ...
0104   float sin2th = sqr(sinTheta);
0105   float minC00 = 1.0;
0106   C[0][0] = std::max(sin2th / sqr(ptMin), minC00);
0107   float zErr = vertexErr.czz();
0108   float transverseErr = vertexErr.cxx();  // assume equal cxx cyy
0109   C[3][3] = transverseErr;
0110   C[4][4] = zErr * sin2th + transverseErr * (1 - sin2th);
0111
0112   return CurvilinearTrajectoryError(C);
0113 }
0114
0115 const TrajectorySeed* SeedForPhotonConversion1Leg::buildSeed(TrajectorySeedCollection& seedCollection,
0116                                                              const SeedingHitSet& hits,
0117                                                              const FreeTrajectoryState& fts,
0118                                                              const edm::EventSetup& es) const {
0119   // FIXME all this stuff shoould go in an initialized...
0120
0121   // get tracker
0122   const auto* tracker = &es.getData(theTrackerToken);
0123
0124   // get propagator
0125   const Propagator* propagator = &es.getData(thePropagatorToken);
0126
0127   auto builder = static_cast<TkTransientTrackingRecHitBuilder const*>(&es.getData(theTTRHBuilderToken));
0128   auto cloner = (*builder).cloner();
0129
0130   // get updator
0131   KFUpdator updator;
0132
0133   // Now update initial state track using information from seed hits.
0134
0135   TrajectoryStateOnSurface updatedState;
0136   edm::OwnVector<TrackingRecHit> seedHits;
0137
0138   const TrackingRecHit* hit = nullptr;
0139   for (unsigned int iHit = 0; iHit < hits.size() && iHit < 1; iHit++) {
0140     hit = hits[iHit];
0141     TrajectoryStateOnSurface state =
0142         (iHit == 0) ? propagator->propagate(fts, tracker->idToDet(hit->geographicalId())->surface())
0143                     : propagator->propagate(updatedState, tracker->idToDet(hit->geographicalId())->surface());
0144     if (!state.isValid())
0145       return nullptr;
0146
0147     SeedingHitSet::ConstRecHitPointer tth = hits[iHit];
0148
0149     std::unique_ptr<BaseTrackerRecHit> newtth(refitHit(tth, state, cloner));
0150
0151     if (!checkHit(state, &*newtth, es))
0152       return nullptr;
0153
0154     updatedState = updator.update(state, *newtth);
0155     if (!updatedState.isValid())
0156       return nullptr;
0157
0158     seedHits.push_back(newtth.release());
0159 #ifdef mydebug_seed
0160     uint32_t detid = hit->geographicalId().rawId();
0161     (*pss) << "\n[SeedForPhotonConversion1Leg] hit " << iHit;
0162     po.print(*pss, detid);
0163     (*pss) << " " << detid << "\t lp " << hit->localPosition() << " tth " << tth->localPosition() << " newtth "
0164            << newtth->localPosition() << " state " << state.globalMomentum().perp();
0165 #endif
0166   }
0167
0168   if (!hit)
0169     return nullptr;
0170
0171   PTrajectoryStateOnDet const& PTraj =
0172       trajectoryStateTransform::persistentState(updatedState, hit->geographicalId().rawId());
0173
0174   seedCollection.push_back(TrajectorySeed(PTraj, seedHits, alongMomentum));
0175   return &seedCollection.back();
0176 }
0177
0178 SeedingHitSet::RecHitPointer SeedForPhotonConversion1Leg::refitHit(SeedingHitSet::ConstRecHitPointer hit,
0179                                                                    const TrajectoryStateOnSurface& state,
0180                                                                    const TkClonerImpl& cloner) const {
0181   //const TransientTrackingRecHit* a= hit.get();
0182   //return const_cast<TransientTrackingRecHit*> (a);
0183   //This was modified otherwise the rechit will have just the local x component and local y=0
0184   // To understand how to modify for pixels
0185
0186   //const TSiStripRecHit2DLocalPos* b = dynamic_cast<const TSiStripRecHit2DLocalPos*>(a);
0187   //return const_cast<TSiStripRecHit2DLocalPos*>(b);
0188   return (SeedingHitSet::RecHitPointer)(cloner(*hit, state));
0189 }