| 1 |
#include "edwenger/HiTrkEffAnalyzer/interface/HiHackedAnalyticalTrackSelector.h"
|
| 2 |
|
| 3 |
#include <Math/DistFunc.h>
|
| 4 |
#include "TMath.h"
|
| 5 |
|
| 6 |
using reco::modules::HiHackedAnalyticalTrackSelector;
|
| 7 |
|
| 8 |
HiHackedAnalyticalTrackSelector::HiHackedAnalyticalTrackSelector( const edm::ParameterSet & cfg ) :
|
| 9 |
src_( cfg.getParameter<edm::InputTag>( "src" ) ),
|
| 10 |
beamspot_( cfg.getParameter<edm::InputTag>( "beamspot" ) ),
|
| 11 |
vertices_( cfg.getParameter<edm::InputTag>( "vertices" ) ),
|
| 12 |
copyExtras_(cfg.getUntrackedParameter<bool>("copyExtras", false)),
|
| 13 |
copyTrajectories_(cfg.getUntrackedParameter<bool>("copyTrajectories", false)),
|
| 14 |
keepAllTracks_( cfg.exists("keepAllTracks") ?
|
| 15 |
cfg.getParameter<bool>("keepAllTracks") :
|
| 16 |
false ), // as this is what you expect from a well behaved selector
|
| 17 |
setQualityBit_( false ),
|
| 18 |
qualityToSet_( TrackBase::undefQuality ),
|
| 19 |
min_relpterr_( cfg.getParameter<double>("min_relpterr") ),
|
| 20 |
min_nhits_( cfg.getParameter<uint32_t>("min_nhits") ),
|
| 21 |
vtxNumber_( cfg.getParameter<int32_t>("vtxNumber") ),
|
| 22 |
vtxTracks_( cfg.getParameter<uint32_t>("vtxTracks") ),
|
| 23 |
vtxChi2Prob_( cfg.getParameter<double>("vtxChi2Prob") ),
|
| 24 |
// parameters for adapted optimal cuts on chi2 and primary vertex compatibility
|
| 25 |
res_par_(cfg.getParameter< std::vector<double> >("res_par") ),
|
| 26 |
chi2n_par_( cfg.getParameter<double>("chi2n_par") ),
|
| 27 |
d0_par1_(cfg.getParameter< std::vector<double> >("d0_par1")),
|
| 28 |
dz_par1_(cfg.getParameter< std::vector<double> >("dz_par1")),
|
| 29 |
d0_par2_(cfg.getParameter< std::vector<double> >("d0_par2")),
|
| 30 |
dz_par2_(cfg.getParameter< std::vector<double> >("dz_par2")),
|
| 31 |
// Boolean indicating if adapted primary vertex compatibility cuts are to be applied.
|
| 32 |
applyAdaptedPVCuts_(cfg.getParameter<bool>("applyAdaptedPVCuts")),
|
| 33 |
// Impact parameter absolute cuts.
|
| 34 |
max_d0_(cfg.getParameter<double>("max_d0")),
|
| 35 |
max_z0_(cfg.getParameter<double>("max_z0")),
|
| 36 |
// Cuts on numbers of layers with hits/3D hits/lost hits.
|
| 37 |
min_layers_(cfg.getParameter<uint32_t>("minNumberLayers") ),
|
| 38 |
min_3Dlayers_(cfg.getParameter<uint32_t>("minNumber3DLayers") ),
|
| 39 |
max_lostLayers_(cfg.getParameter<uint32_t>("maxNumberLostLayers") )
|
| 40 |
{
|
| 41 |
if (cfg.exists("qualityBit")) {
|
| 42 |
std::string qualityStr = cfg.getParameter<std::string>("qualityBit");
|
| 43 |
if (qualityStr != "") {
|
| 44 |
setQualityBit_ = true;
|
| 45 |
qualityToSet_ = TrackBase::qualityByName(cfg.getParameter<std::string>("qualityBit"));
|
| 46 |
}
|
| 47 |
}
|
| 48 |
if (keepAllTracks_ && !setQualityBit_) throw cms::Exception("Configuration") <<
|
| 49 |
"If you set 'keepAllTracks' to true, you must specify which qualityBit to set.\n";
|
| 50 |
if (setQualityBit_ && (qualityToSet_ == TrackBase::undefQuality)) throw cms::Exception("Configuration") <<
|
| 51 |
"You can't set the quality bit " << cfg.getParameter<std::string>("qualityBit") << " as it is 'undefQuality' or unknown.\n";
|
| 52 |
|
| 53 |
std::string alias( cfg.getParameter<std::string>( "@module_label" ) );
|
| 54 |
produces<reco::TrackCollection>().setBranchAlias( alias + "Tracks");
|
| 55 |
if (copyExtras_) {
|
| 56 |
produces<reco::TrackExtraCollection>().setBranchAlias( alias + "TrackExtras");
|
| 57 |
produces<TrackingRecHitCollection>().setBranchAlias( alias + "RecHits");
|
| 58 |
}
|
| 59 |
if (copyTrajectories_) {
|
| 60 |
produces< std::vector<Trajectory> >().setBranchAlias( alias + "Trajectories");
|
| 61 |
produces< TrajTrackAssociationCollection >().setBranchAlias( alias + "TrajectoryTrackAssociations");
|
| 62 |
}
|
| 63 |
|
| 64 |
}
|
| 65 |
|
| 66 |
HiHackedAnalyticalTrackSelector::~HiHackedAnalyticalTrackSelector() {
|
| 67 |
}
|
| 68 |
|
| 69 |
void HiHackedAnalyticalTrackSelector::produce( edm::Event& evt, const edm::EventSetup& es )
|
| 70 |
{
|
| 71 |
using namespace std;
|
| 72 |
using namespace edm;
|
| 73 |
using namespace reco;
|
| 74 |
|
| 75 |
Handle<TrackCollection> hSrcTrack;
|
| 76 |
Handle< vector<Trajectory> > hTraj;
|
| 77 |
Handle< vector<Trajectory> > hTrajP;
|
| 78 |
Handle< TrajTrackAssociationCollection > hTTAss;
|
| 79 |
|
| 80 |
bool isTrajThere = evt.getByLabel(src_, hTraj);
|
| 81 |
|
| 82 |
// looking for the beam spot
|
| 83 |
edm::Handle<reco::BeamSpot> hBsp;
|
| 84 |
evt.getByLabel(beamspot_, hBsp);
|
| 85 |
reco::BeamSpot vertexBeamSpot;
|
| 86 |
vertexBeamSpot = *hBsp;
|
| 87 |
|
| 88 |
// Select good primary vertices for use in subsequent track selection
|
| 89 |
edm::Handle<reco::VertexCollection> hVtx;
|
| 90 |
evt.getByLabel(vertices_, hVtx);
|
| 91 |
std::vector<Point> points;
|
| 92 |
std::vector<double> vterr, vzerr;
|
| 93 |
selectVertices(*hVtx, points, vterr, vzerr);
|
| 94 |
|
| 95 |
// Get tracks
|
| 96 |
evt.getByLabel( src_, hSrcTrack );
|
| 97 |
|
| 98 |
selTracks_ = auto_ptr<TrackCollection>(new TrackCollection());
|
| 99 |
rTracks_ = evt.getRefBeforePut<TrackCollection>();
|
| 100 |
if (copyExtras_) {
|
| 101 |
selTrackExtras_ = auto_ptr<TrackExtraCollection>(new TrackExtraCollection());
|
| 102 |
selHits_ = auto_ptr<TrackingRecHitCollection>(new TrackingRecHitCollection());
|
| 103 |
rHits_ = evt.getRefBeforePut<TrackingRecHitCollection>();
|
| 104 |
rTrackExtras_ = evt.getRefBeforePut<TrackExtraCollection>();
|
| 105 |
}
|
| 106 |
|
| 107 |
if (isTrajThere && copyTrajectories_) trackRefs_.resize(hSrcTrack->size());
|
| 108 |
|
| 109 |
// Loop over tracks
|
| 110 |
size_t current = 0;
|
| 111 |
for (TrackCollection::const_iterator it = hSrcTrack->begin(), ed = hSrcTrack->end(); it != ed; ++it, ++current) {
|
| 112 |
const Track & trk = * it;
|
| 113 |
// Check if this track passes cuts
|
| 114 |
bool ok = select(vertexBeamSpot, trk, points, vterr, vzerr);
|
| 115 |
if (!ok) {
|
| 116 |
if (isTrajThere && copyTrajectories_) trackRefs_[current] = reco::TrackRef();
|
| 117 |
if (!keepAllTracks_) continue;
|
| 118 |
}
|
| 119 |
selTracks_->push_back( Track( trk ) ); // clone and store
|
| 120 |
if (ok && setQualityBit_) selTracks_->back().setQuality(qualityToSet_);
|
| 121 |
if (copyExtras_) {
|
| 122 |
// TrackExtras
|
| 123 |
selTrackExtras_->push_back( TrackExtra( trk.outerPosition(), trk.outerMomentum(), trk.outerOk(),
|
| 124 |
trk.innerPosition(), trk.innerMomentum(), trk.innerOk(),
|
| 125 |
trk.outerStateCovariance(), trk.outerDetId(),
|
| 126 |
trk.innerStateCovariance(), trk.innerDetId(),
|
| 127 |
trk.seedDirection(), trk.seedRef() ) );
|
| 128 |
selTracks_->back().setExtra( TrackExtraRef( rTrackExtras_, selTrackExtras_->size() - 1) );
|
| 129 |
TrackExtra & tx = selTrackExtras_->back();
|
| 130 |
tx.setResiduals(trk.residuals());
|
| 131 |
// TrackingRecHits
|
| 132 |
for( trackingRecHit_iterator hit = trk.recHitsBegin(); hit != trk.recHitsEnd(); ++ hit ) {
|
| 133 |
selHits_->push_back( (*hit)->clone() );
|
| 134 |
tx.add( TrackingRecHitRef( rHits_, selHits_->size() - 1) );
|
| 135 |
}
|
| 136 |
}
|
| 137 |
if (isTrajThere && copyTrajectories_) {
|
| 138 |
trackRefs_[current] = TrackRef(rTracks_, selTracks_->size() - 1);
|
| 139 |
}
|
| 140 |
}
|
| 141 |
if (isTrajThere && copyTrajectories_ ) {
|
| 142 |
Handle< vector<Trajectory> > hTraj;
|
| 143 |
Handle< TrajTrackAssociationCollection > hTTAss;
|
| 144 |
evt.getByLabel(src_, hTTAss);
|
| 145 |
evt.getByLabel(src_, hTraj); // it's gotten up there
|
| 146 |
selTrajs_ = auto_ptr< vector<Trajectory> >(new vector<Trajectory>());
|
| 147 |
rTrajectories_ = evt.getRefBeforePut< vector<Trajectory> >();
|
| 148 |
selTTAss_ = auto_ptr< TrajTrackAssociationCollection >(new TrajTrackAssociationCollection());
|
| 149 |
for (size_t i = 0, n = hTraj->size(); i < n; ++i) {
|
| 150 |
Ref< vector<Trajectory> > trajRef(hTraj, i);
|
| 151 |
TrajTrackAssociationCollection::const_iterator match = hTTAss->find(trajRef);
|
| 152 |
if (match != hTTAss->end()) {
|
| 153 |
const Ref<TrackCollection> &trkRef = match->val;
|
| 154 |
short oldKey = static_cast<short>(trkRef.key());
|
| 155 |
if (trackRefs_[oldKey].isNonnull()) {
|
| 156 |
selTrajs_->push_back( Trajectory(*trajRef) );
|
| 157 |
selTTAss_->insert ( Ref< vector<Trajectory> >(rTrajectories_, selTrajs_->size() - 1), trackRefs_[oldKey] );
|
| 158 |
}
|
| 159 |
}
|
| 160 |
}
|
| 161 |
}
|
| 162 |
|
| 163 |
static const std::string emptyString;
|
| 164 |
evt.put(selTracks_);
|
| 165 |
if (copyExtras_ ) {
|
| 166 |
evt.put(selTrackExtras_);
|
| 167 |
evt.put(selHits_);
|
| 168 |
}
|
| 169 |
if (isTrajThere && copyTrajectories_ ) {
|
| 170 |
evt.put(selTrajs_);
|
| 171 |
evt.put(selTTAss_);
|
| 172 |
}
|
| 173 |
}
|
| 174 |
|
| 175 |
|
| 176 |
bool HiHackedAnalyticalTrackSelector::select(const reco::BeamSpot &vertexBeamSpot, const reco::Track &tk, const std::vector<Point> &points, std::vector<double> &vterr, std::vector<double> &vzerr) {
|
| 177 |
// Decide if the given track passes selection cuts.
|
| 178 |
|
| 179 |
using namespace std;
|
| 180 |
|
| 181 |
// Cuts on numbers of layers with hits/3D hits/lost hits.
|
| 182 |
uint32_t nlayers = tk.hitPattern().trackerLayersWithMeasurement();
|
| 183 |
uint32_t nlayers3D = tk.hitPattern().pixelLayersWithMeasurement() +
|
| 184 |
tk.hitPattern().numberOfValidStripLayersWithMonoAndStereo();
|
| 185 |
uint32_t nlayersLost = tk.hitPattern().trackerLayersWithoutMeasurement();
|
| 186 |
if (nlayers < min_layers_) return false;
|
| 187 |
if (nlayers3D < min_3Dlayers_) return false;
|
| 188 |
if (nlayersLost > max_lostLayers_) return false;
|
| 189 |
|
| 190 |
// Get track parameters
|
| 191 |
double pt = tk.pt(),eta = tk.eta(), chi2n = tk.normalizedChi2();
|
| 192 |
double d0 = -tk.dxy(vertexBeamSpot.position()), d0E = tk.d0Error(),
|
| 193 |
dz = tk.dz(vertexBeamSpot.position()), dzE = tk.dzError();
|
| 194 |
|
| 195 |
double relpterr = tk.ptError()/pt;
|
| 196 |
uint32_t nhits = tk.numberOfValidHits();
|
| 197 |
|
| 198 |
if(relpterr > min_relpterr_) return false;
|
| 199 |
if(nhits < min_nhits_) return false;
|
| 200 |
|
| 201 |
// optimized cuts adapted to the track nlayers, pt, eta:
|
| 202 |
// cut on chiquare/ndof
|
| 203 |
if (chi2n > chi2n_par_*nlayers) return false;
|
| 204 |
|
| 205 |
|
| 206 |
// parametrized d0 resolution for the track pt
|
| 207 |
double nomd0E = sqrt(res_par_[0]*res_par_[0]+(res_par_[1]/max(pt,1e-9))*(res_par_[1]/max(pt,1e-9)));
|
| 208 |
// parametrized z0 resolution for the track pt and eta
|
| 209 |
double nomdzE = nomd0E*(std::cosh(eta));
|
| 210 |
|
| 211 |
|
| 212 |
// ---- PrimaryVertex compatibility cut
|
| 213 |
bool primaryVertexZCompatibility(false);
|
| 214 |
bool primaryVertexD0Compatibility(false);
|
| 215 |
|
| 216 |
if (points.empty()) { //If not primaryVertices are reconstructed, check just the compatibility with the BS
|
| 217 |
//z0 within three sigma of the beam spot z, if no good vertex is found
|
| 218 |
if ( abs(dz) < (vertexBeamSpot.sigmaZ()*3) ) primaryVertexZCompatibility = true;
|
| 219 |
|
| 220 |
// d0 compatibility with beam line
|
| 221 |
if (abs(d0) < pow(d0_par1_[0]*nlayers,d0_par1_[1])*nomd0E &&
|
| 222 |
abs(d0) < pow(d0_par2_[0]*nlayers,d0_par2_[1])*d0E) primaryVertexD0Compatibility = true;
|
| 223 |
}
|
| 224 |
|
| 225 |
|
| 226 |
int iv=0;
|
| 227 |
for (std::vector<Point>::const_iterator point = points.begin(), end = points.end(); point != end; ++point) {
|
| 228 |
if(primaryVertexZCompatibility && primaryVertexD0Compatibility) break;
|
| 229 |
double dzPV = tk.dz(*point); //re-evaluate the dz with respect to the vertex position
|
| 230 |
double d0PV = tk.dxy(*point); //re-evaluate the dxy with respect to the vertex position
|
| 231 |
|
| 232 |
/*
|
| 233 |
if (abs(dzPV) < pow(dz_par1_[0]*nlayers,dz_par1_[1])*nomdzE &&
|
| 234 |
abs(dzPV) < pow(dz_par2_[0]*nlayers,dz_par2_[1])*dzE ) primaryVertexZCompatibility = true;
|
| 235 |
|
| 236 |
if (abs(d0PV) < pow(d0_par1_[0]*nlayers,d0_par1_[1])*nomd0E &&
|
| 237 |
abs(d0PV) < pow(d0_par2_[0]*nlayers,d0_par2_[1])*d0E) primaryVertexD0Compatibility = true;
|
| 238 |
*/
|
| 239 |
|
| 240 |
// Edward's temporary changes
|
| 241 |
double dzErrPV = sqrt(dzE*dzE+vzerr[iv]*vzerr[iv]);
|
| 242 |
double d0ErrPV = sqrt(d0E*d0E+vterr[iv]*vterr[iv]);
|
| 243 |
iv++;
|
| 244 |
|
| 245 |
// Wei's temporary changes (plus Edward's max cuts)
|
| 246 |
if (abs(dzPV) < dz_par1_[0]*pow(nlayers,dz_par1_[1])*nomdzE &&
|
| 247 |
abs(dzPV) < dz_par2_[0]*pow(nlayers,dz_par2_[1])*dzErrPV &&
|
| 248 |
abs(dzPV) < max_z0_) primaryVertexZCompatibility = true;
|
| 249 |
|
| 250 |
if (abs(d0PV) < d0_par1_[0]*pow(nlayers,d0_par1_[1])*nomd0E &&
|
| 251 |
abs(d0PV) < d0_par2_[0]*pow(nlayers,d0_par2_[1])*d0ErrPV &&
|
| 252 |
abs(d0PV) < max_d0_) primaryVertexD0Compatibility = true;
|
| 253 |
|
| 254 |
|
| 255 |
}
|
| 256 |
|
| 257 |
|
| 258 |
// Absolute cuts on all tracks impact parameters with respect to beam-spot.
|
| 259 |
// If BS is not compatible, verify if at least the reco-vertex is compatible (useful for incorrect BS settings)
|
| 260 |
if (abs(d0) > max_d0_ && !primaryVertexD0Compatibility) return false;
|
| 261 |
if (abs(dz) > max_z0_ && !primaryVertexZCompatibility) return false;
|
| 262 |
|
| 263 |
reco::TrackBase::TrackQuality trackQualityTight = TrackBase::qualityByName("tight");
|
| 264 |
if(!tk.quality(trackQualityTight)) return false;
|
| 265 |
|
| 266 |
if (applyAdaptedPVCuts_) {
|
| 267 |
return (primaryVertexD0Compatibility && primaryVertexZCompatibility);
|
| 268 |
} else {
|
| 269 |
return true;
|
| 270 |
}
|
| 271 |
|
| 272 |
|
| 273 |
}
|
| 274 |
void HiHackedAnalyticalTrackSelector::selectVertices(const reco::VertexCollection &vtxs, std::vector<Point> &points, std::vector<double> &vterr, std::vector<double> &vzerr) {
|
| 275 |
// Select good primary vertices
|
| 276 |
using namespace reco;
|
| 277 |
int32_t toTake = vtxNumber_;
|
| 278 |
for (VertexCollection::const_iterator it = vtxs.begin(), ed = vtxs.end(); it != ed; ++it) {
|
| 279 |
if ((it->tracksSize() >= vtxTracks_) &&
|
| 280 |
( (it->chi2() == 0.0) || (TMath::Prob(it->chi2(), static_cast<int32_t>(it->ndof()) ) >= vtxChi2Prob_) ) ) {
|
| 281 |
points.push_back(it->position());
|
| 282 |
vterr.push_back(sqrt(it->yError()*it->xError()));
|
| 283 |
vzerr.push_back(it->zError());
|
| 284 |
toTake--; if (toTake == 0) break;
|
| 285 |
}
|
| 286 |
}
|
| 287 |
}
|
| 288 |
|
