TrackTrigger/src/StubPtConsistency.cc

0001 #include <iostream>
0002 #include <memory>
0003
0004 #include "L1Trigger/TrackTrigger/interface/StubPtConsistency.h"
0005 #include "CLHEP/Units/PhysicalConstants.h"
0006
0007 namespace StubPtConsistency {
0008
0009   float getConsistency(TTTrack<Ref_Phase2TrackerDigi_> aTrack,
0010                        const TrackerGeometry* theTrackerGeom,
0011                        const TrackerTopology* tTopo,
0012                        double mMagneticFieldStrength,
0013                        int nPar) {
0014     if (!(nPar == 4 || nPar == 5))
0015       throw cms::Exception("IncorrectInput") << "Not a valid nPar option!";
0016
0017     double trk_bendchi2 = 0.0;
0018     double bend_resolution = 0.483;
0019     // B*c/2E11 - converts q/pt to track angle at some radius from beamline
0020     float speedOfLightConverted = CLHEP::c_light / 1.0E5;
0021
0022     // Need the pT signed in order to determine if bend is positive or negative
0023     // P(MeV/c) = (c/10^9)·Q·B(kG)·R(cm)
0024     float trk_signedPt = speedOfLightConverted * mMagneticFieldStrength / aTrack.rInv();
0025
0026     // loop over stubs
0027     const auto& stubRefs = aTrack.getStubRefs();
0028     int nStubs = stubRefs.size();
0029
0030     for (const auto& stubRef : stubRefs) {
0031       DetId detIdStub = theTrackerGeom->idToDet((stubRef->clusterRef(0))->getDetId())->geographicalId();
0032       MeasurementPoint coords = stubRef->clusterRef(0)->findAverageLocalCoordinatesCentered();
0033       const GeomDet* theGeomDet = theTrackerGeom->idToDet(detIdStub);
0034       Global3DPoint posStub = theGeomDet->surface().toGlobal(theGeomDet->topology().localPosition(coords));
0035
0036       float stub_r = posStub.perp();
0037       float stub_z = posStub.z();
0038
0039       bool isBarrel = (detIdStub.subdetId() == StripSubdetector::TOB);
0040
0041       const GeomDetUnit* det0 = theTrackerGeom->idToDetUnit(detIdStub);
0042       const GeomDetUnit* det1 = theTrackerGeom->idToDetUnit(tTopo->partnerDetId(detIdStub));
0043       const PixelGeomDetUnit* unit = reinterpret_cast<const PixelGeomDetUnit*>(det0);
0044       const PixelTopology& topo = unit->specificTopology();
0045
0046       // Calculation of snesor spacing obtained from TMTT: https://github.com/CMS-TMTT/cmssw/blob/TMTT_938/L1Trigger/TrackFindingTMTT/src/Stub.cc#L138-L146
0047       float stripPitch = topo.pitch().first;
0048
0049       float modMinR = std::min(det0->position().perp(), det1->position().perp());
0050       float modMaxR = std::max(det0->position().perp(), det1->position().perp());
0051       float modMinZ = std::min(det0->position().z(), det1->position().z());
0052       float modMaxZ = std::max(det0->position().z(), det1->position().z());
0053       float sensorSpacing = sqrt((modMaxR - modMinR) * (modMaxR - modMinR) + (modMaxZ - modMinZ) * (modMaxZ - modMinZ));
0054
0055       // Approximation of phiOverBendCorrection, from TMTT: https://github.com/CMS-TMTT/cmssw/blob/TMTT_938/L1Trigger/TrackFindingTMTT/src/Stub.cc#L440-L448
0056       bool tiltedBarrel = (isBarrel && tTopo->tobSide(detIdStub) != 3);
0057       float gradient = 0.886454;
0058       float intercept = 0.504148;
0059       float correction;
0060       if (tiltedBarrel)
0061         correction = gradient * std::abs(stub_z) / stub_r + intercept;
0062       else if (isBarrel)
0063         correction = 1;
0064       else
0065         correction = std::abs(stub_z) / stub_r;
0066
0067       float stubBend = stubRef->bendFE();
0068       if (!isBarrel && stub_z < 0.0)
0069         stubBend = -stubBend;  // flip sign of bend if in negative end cap
0070
0071       float trackBend = -(sensorSpacing * stub_r * mMagneticFieldStrength * (speedOfLightConverted / 2)) /
0072                         (stripPitch * trk_signedPt * correction);
0073       float bendDiff = trackBend - stubBend;
0074
0075       trk_bendchi2 += (bendDiff * bendDiff) / (bend_resolution * bend_resolution);
0076     }  // end loop over stubs
0077
0078     float bendchi2 = trk_bendchi2 / nStubs;
0079     return bendchi2;
0080   }  //end getConsistency()
0081 }  // namespace StubPtConsistency