SiStripRecHitConverter/src/StripCPE.cc

0001 #include "RecoLocalTracker/SiStripRecHitConverter/interface/StripCPE.h"
0002 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0003 #include "Geometry/CommonTopologies/interface/TkRadialStripTopology.h"
0004 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetType.h"
0005 #include <algorithm>
0006 #include <cmath>
0007 #include <cassert>
0008
0009 StripCPE::StripCPE(edm::ParameterSet& conf,
0010                    const MagneticField& mag,
0011                    const TrackerGeometry& geom,
0012                    const SiStripLorentzAngle& LorentzAngle,
0013                    const SiStripBackPlaneCorrection& BackPlaneCorrection,
0014                    const SiStripConfObject& confObj,
0015                    const SiStripLatency& latency)
0016     : peakMode_(latency.singleReadOutMode() == 1),
0017       geom_(geom),
0018       magfield_(mag),
0019       LorentzAngleMap_(LorentzAngle),
0020       BackPlaneCorrectionMap_(BackPlaneCorrection) {
0021   typedef std::map<std::string, SiStripModuleGeometry> map_t;
0022   map_t modules;
0023   modules["IB1"] = SiStripModuleGeometry::IB1;
0024   modules["IB2"] = SiStripModuleGeometry::IB2;
0025   modules["OB1"] = SiStripModuleGeometry::OB1;
0026   modules["OB2"] = SiStripModuleGeometry::OB2;
0027   modules["W1A"] = SiStripModuleGeometry::W1A;
0028   modules["W2A"] = SiStripModuleGeometry::W2A;
0029   modules["W3A"] = SiStripModuleGeometry::W3A;
0030   modules["W1B"] = SiStripModuleGeometry::W1B;
0031   modules["W2B"] = SiStripModuleGeometry::W2B;
0032   modules["W3B"] = SiStripModuleGeometry::W3B;
0033   modules["W4"] = SiStripModuleGeometry::W4;
0034   modules["W5"] = SiStripModuleGeometry::W5;
0035   modules["W6"] = SiStripModuleGeometry::W6;
0036   modules["W7"] = SiStripModuleGeometry::W7;
0037
0038   const unsigned size =
0039       static_cast<unsigned int>(
0040           max_element(modules.begin(), modules.end(), [&](auto& arg1, auto& arg2) { return arg1.second < arg2.second; })
0041               ->second) +
0042       1;
0043   xtalk1.resize(size);
0044   xtalk2.resize(size);
0045
0046   for (map_t::const_iterator it = modules.begin(); it != modules.end(); it++) {
0047     const std::string modeS(peakMode_ ? "Peak" : "Deco"), xtalk1S("xtalk1_" + it->first + modeS),
0048         xtalk2S("xtalk2_" + it->first + modeS);
0049
0050     if (!confObj.isParameter(xtalk1S))
0051       throw cms::Exception("SiStripConfObject does not contain: ") << xtalk1S;
0052     if (!confObj.isParameter(xtalk2S))
0053       throw cms::Exception("SiStripConfObject does not contain: ") << xtalk2S;
0054
0055     xtalk1[static_cast<int>(it->second)] = confObj.get<double>(xtalk1S);
0056     xtalk2[static_cast<int>(it->second)] = confObj.get<double>(xtalk2S);
0057   }
0058
0059   fillParams();
0060 }
0061
0062 StripClusterParameterEstimator::LocalValues StripCPE::localParameters(const SiStripCluster& cluster,
0063                                                                       const GeomDetUnit& det) const {
0064   StripCPE::Param const& p = param(det);
0065   const float barycenter = cluster.barycenter();
0066   const float fullProjection =
0067       p.coveredStrips(p.drift + LocalVector(0, 0, -p.thickness), p.topology->localPosition(barycenter));
0068   const float strip = barycenter - 0.5f * (1.f - p.backplanecorrection) * fullProjection;
0069
0070   return std::make_pair(p.topology->localPosition(strip), p.topology->localError(strip, 1.f / 12.f));
0071 }
0072
0073 float StripCPE::Param::coveredStrips(const LocalVector& lvec, const LocalPoint& lpos) const {
0074   return topology->coveredStrips(lpos + 0.5f * lvec, lpos - 0.5f * lvec);
0075 }
0076
0077 LocalVector StripCPE::driftDirection(const StripGeomDetUnit* det) const {
0078   LocalVector lbfield = (det->surface()).toLocal(magfield_.inTesla(det->surface().position()));
0079
0080   float tanLorentzAnglePerTesla = LorentzAngleMap_.getLorentzAngle(det->geographicalId().rawId());
0081
0082   float dir_x = -tanLorentzAnglePerTesla * lbfield.y();
0083   float dir_y = tanLorentzAnglePerTesla * lbfield.x();
0084   float dir_z = 1.f;  // E field always in z direction
0085
0086   return LocalVector(dir_x, dir_y, dir_z);
0087 }
0088
0089 void StripCPE::fillParams() {
0090   auto const& dus = geom_.detUnits();
0091   m_off = dus.size();
0092   for (unsigned int i = 1; i < 7; ++i) {
0093     LogDebug("LookingForFirstStrip")
0094         << "Subdetector " << i << " GeomDetEnumerator " << GeomDetEnumerators::tkDetEnum[i] << " offset "
0095         << geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) << " is it strip? "
0096         << (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) != dus.size()
0097                 ? dus[geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i])]->type().isTrackerStrip()
0098                 : false);
0099     if (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) != dus.size() &&
0100         dus[geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i])]->type().isTrackerStrip()) {
0101       if (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) < m_off)
0102         m_off = geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]);
0103     }
0104   }
0105   LogDebug("LookingForFirstStrip") << " Chosen offset: " << m_off;
0106   m_Params.resize(dus.size() - m_off);
0107   for (auto i = m_off; i != dus.size(); ++i) {
0108     auto& p = m_Params[i - m_off];
0109     const StripGeomDetUnit* stripdet = (const StripGeomDetUnit*)(dus[i]);
0110     assert(stripdet->index() == int(i));
0111     //    assert(stripdet->geographicalId().subdetId()>1); // not pixel..
0112     assert(stripdet->type().isTrackerStrip());  // not pixel
0113
0114     const Bounds& bounds = stripdet->specificSurface().bounds();
0115     p.maxLength = std::sqrt(std::pow(bounds.length(), 2.f) + std::pow(bounds.width(), 2.f));
0116     p.thickness = bounds.thickness();
0117     p.invThickness = 1.f / p.thickness;
0118     p.drift = driftDirection(stripdet) * p.thickness;
0119     p.topology = (StripTopology*)(&stripdet->topology());
0120     p.nstrips = p.topology->nstrips();
0121     p.moduleGeom = SiStripDetId(stripdet->geographicalId()).moduleGeometry();
0122     p.backplanecorrection = BackPlaneCorrectionMap_.getBackPlaneCorrection(stripdet->geographicalId().rawId());
0123
0124     const TkRadialStripTopology* rtop =
0125         dynamic_cast<const TkRadialStripTopology*>(&stripdet->specificType().specificTopology());
0126     p.pitch_rel_err2 =
0127         (rtop)
0128             ? pow(0.5f * rtop->angularWidth() * rtop->stripLength() / rtop->localPitch(LocalPoint(0, 0, 0)), 2.f) / 12.f
0129             : 0.f;
0130   }
0131 }