001 #include "JetTagTools/TrackSelector.h"
002
003 #include "Particle/TrackParticle.h"
004 #include "Particle/TrackParticleContainer.h"
005 #include "TrkTrackSummary/TrackSummary.h"
006 #include "TrkTrack/Track.h"
007 #include "TrkEventPrimitives/FitQuality.h"
008 #include "CLHEP/GenericFunctions/CumulativeChiSquare.hh"
009 #include "CLHEP/Vector/ThreeVector.h"
010 #include "GaudiKernel/IToolSvc.h"
011 #include "ITrackToVertex/ITrackToVertex.h"
012 #include <TMath.h>
013 #include <string>
014 #include <bitset>
015 #include <algorithm>
016
017 #include "InDetTestBLayer/InDetTestBLayerTool.h"
018
019 namespace Analysis {
020
021 static const InterfaceID IID_ITrackSelector("Analysis::TrackSelector", 1, 0);
022
023 TrackSelector::TrackSelector(const std::string& type,
024 const std::string& name, const IInterface* parent) :
025 AthAlgTool(type, name, parent),
026 m_primaryVertex(Hep3Vector()),
027 m_ntri(0),
028 m_ntrf(0),
029 m_trackToVertexTool("Reco::TrackToVertex"),
030 m_inDetTestBLayerTool("InDet::InDetTestBLayerTool") {
031
032 declareInterface<TrackSelector>(this);
033 declareProperty("trackToVertexTool", m_trackToVertexTool);
034 declareProperty("inDetTestBLayerTool", m_inDetTestBLayerTool);
035 declareProperty("useBLayerHitPrediction", m_useBLayerHitPrediction = false);
036 declareProperty("usePerigeeParameters", m_usePerigeeParameters = false);
037 declareProperty("pTMin", m_pTMin = 1.*GeV);
038 declareProperty("d0Max", m_d0Max = 1.*mm);
039 declareProperty("z0Max", m_z0Max = 1.5*mm);
040 declareProperty("sigd0Max",m_sigd0Max = 999.*mm);
041 declareProperty("sigz0Max",m_sigz0Max = 999.*mm);
042 declareProperty("etaMax", m_etaMax = 9999.);
043 declareProperty("useTrackSummaryInfo", m_useTrackSummaryInfo = true);
044 declareProperty("nHitBLayer", m_nHitBLayer = 1);
045 declareProperty("nHitPix", m_nHitPix = 2);
046 declareProperty("nHitSct", m_nHitSct = 0);
047 declareProperty("nHitSi", m_nHitSi = 7);
048 declareProperty("nHitTrt", m_nHitTrt = 0);
049 declareProperty("nHitTrtHighE", m_nHitTrtHighE = 0);
050 declareProperty("useDeadPixInfo", m_useDeadPixInfo = true);
051 declareProperty("useDeadSctInfo", m_useDeadSctInfo = true);
052 declareProperty("useTrackQualityInfo", m_useTrackQualityInfo = true);
053 declareProperty("fitChi2", m_fitChi2 = 99999.);
054 declareProperty("fitProb", m_fitProb = -1.);
055 declareProperty("fitChi2OnNdfMax",m_fitChi2OnNdfMax=999.);
056 declareProperty("inputTrackCollection", m_inputTrackCollection);
057 declareProperty("outputTrackCollection", m_outputTrackCollection);
058 }
059
060 TrackSelector::~TrackSelector() {
061 }
062
063 const InterfaceID& TrackSelector::interfaceID() {
064 return IID_ITrackSelector;
065 }
066
067 StatusCode TrackSelector::initialize() {
068 for(int i=0;i<16;i++) m_ntrc[i]=0;
069
070 IToolSvc* toolSvc;
071 StatusCode sc = service("ToolSvc", toolSvc);
072 if (StatusCode::SUCCESS != sc) {
073 ATH_MSG_ERROR("#BTAG# Can't get ToolSvc");
074 return StatusCode::FAILURE;
075 }
076
077 if ( m_trackToVertexTool.retrieve().isFailure() ) {
078 ATH_MSG_FATAL("#BTAG# Failed to retrieve tool " << m_trackToVertexTool);
079 return StatusCode::FAILURE;
080 } else {
081 ATH_MSG_DEBUG("#BTAG# Retrieved tool " << m_trackToVertexTool);
082 }
083
084
085 if ( m_useBLayerHitPrediction ) {
086 if ( m_inDetTestBLayerTool.retrieve().isFailure() ) {
087 ATH_MSG_FATAL("#BTAG# Failed to retrieve tool " << m_inDetTestBLayerTool );
088 return StatusCode::FAILURE;
089 } else {
090 ATH_MSG_DEBUG("#BTAG# Retrieved tool " << m_inDetTestBLayerTool );
091 }
092 }
093
094
095
096 if (msgLvl(MSG::DEBUG)) {
097 msg(MSG::DEBUG) << "#BTAG# TrackSelector " << name() << " cuts: " << endreq;
098 msg(MSG::DEBUG) << "#BTAG# - pT >= " << m_pTMin << endreq;
099 msg(MSG::DEBUG) << "#BTAG# - |eta| <= " << m_etaMax << endreq;
100 msg(MSG::DEBUG) << "#BTAG# - |d0| <= " << m_d0Max << endreq;
101 msg(MSG::DEBUG) << "#BTAG# - |z0| <= " << m_z0Max << endreq;
102 msg(MSG::DEBUG) << "#BTAG# - |sigd0| <= " << m_sigd0Max << endreq;
103 msg(MSG::DEBUG) << "#BTAG# - |sigz0| <= " << m_sigz0Max << endreq;
104 if(m_useTrackSummaryInfo) {
105 msg(MSG::DEBUG) << "#BTAG# - nbHitsBLayer >= " << m_nHitBLayer << endreq;
106 if(m_useDeadPixInfo) {
107 msg(MSG::DEBUG) << "#BTAG# - nbHitsPix+nbDeadPix >= " << m_nHitPix << endreq;
108 } else {
109 msg(MSG::DEBUG) << "#BTAG# - nbHitsPix >= " << m_nHitPix << endreq;
110 }
111 if(m_useBLayerHitPrediction)
112 msg(MSG::DEBUG) << "#BTAG# using conddb for b-layer hit requirements " << endreq;
113
114 if(m_useDeadSctInfo) {
115 msg(MSG::DEBUG) << "#BTAG# - nbHitsSct+nbDeadSct >= " << m_nHitSct << endreq;
116 } else {
117 msg(MSG::DEBUG) << "#BTAG# - nbHitsSct >= " << m_nHitSct << endreq;
118 }
119 if(m_useDeadPixInfo) {
120 if(m_useDeadSctInfo) {
121 msg(MSG::DEBUG) << "#BTAG# - nbHitsSi+nbDeadPix+nbDeadSct >= " << m_nHitSi << endreq;
122 } else {
123 msg(MSG::DEBUG) << "#BTAG# - nbHitsSi+nbDeadPix >= " << m_nHitSi << endreq;
124 }
125 } else {
126 if(m_useDeadSctInfo) {
127 msg(MSG::DEBUG) << "#BTAG# - nbHitsSi+nbDeadSct >= " << m_nHitSi << endreq;
128 } else {
129 msg(MSG::DEBUG) << "#BTAG# - nbHitsSi >= " << m_nHitSi << endreq;
130 }
131 }
132 msg(MSG::DEBUG) << "#BTAG# - nbHitsTrt >= " << m_nHitTrt << endreq;
133 msg(MSG::DEBUG) << "#BTAG# - nbHitsTrtHighE >= " << m_nHitTrtHighE << endreq;
134 }
135 msg(MSG::DEBUG) << "#BTAG# - fit chi2 <= " << m_fitChi2 << endreq;
136 msg(MSG::DEBUG) << "#BTAG# - fit proba >= " << m_fitProb << endreq;
137 msg(MSG::DEBUG) << "#BTAG# - fit chi2 / ndf <= " << m_fitChi2OnNdfMax << endreq;
138 }
139
140 return StatusCode::SUCCESS;
141 }
142
143
144 StatusCode TrackSelector::finalize() {
145 ATH_MSG_VERBOSE("#BTAG# tracks selected: In= " << m_ntri);
146 for(int i=0;i<16;i++) ATH_MSG_VERBOSE("#BTAG# cut" << i << "= " << m_ntrc[i]);
147 ATH_MSG_VERBOSE("#BTAG# Out= " << m_ntrf);
148 return StatusCode::SUCCESS;
149 }
150
151
152 bool TrackSelector::selectTrack(const Rec::TrackParticle* track) {
153
154
155 enum Cuts { pTMin, d0Max, z0Max, sigd0Max, sigz0Max, etaMax,
156 nHitBLayer, deadBLayer, nHitPix, nHitSct, nHitSi, nHitTrt, nHitTrtHighE,
157 fitChi2, fitProb,fitChi2OnNdfMax,
158 numCuts };
159 std::bitset<numCuts> failedCuts;
160
161 if(m_primaryVertex==Hep3Vector()) {
162 ATH_MSG_WARNING("#BTAG# primary vertex is not set !");
163 }
164
165 double trackD0;
166 double trackZ0;
167 double tracksigD0;
168 double tracksigZ0;
169 if(m_usePerigeeParameters) {
170 trackD0 = track->measuredPerigee()->parameters()[Trk::d0];
171 trackZ0 = track->measuredPerigee()->parameters()[Trk::z0] - m_primaryVertex.z();
172 tracksigD0 = TMath::Sqrt(track->measuredPerigee()->localErrorMatrix().covariance()[Trk::d0][Trk::d0]);
173 tracksigZ0 = TMath::Sqrt(track->measuredPerigee()->localErrorMatrix().covariance()[Trk::z0][Trk::z0]);
174 } else {
175
176 const Trk::MeasuredPerigee* perigee = m_trackToVertexTool->perigeeAtVertex(*track, m_primaryVertex);
177 if (perigee==0) {
178 ATH_MSG_WARNING("#BTAG# Extrapolation failed. Rejecting track... ");
179 return false;
180 }
181 trackD0 = perigee->parameters()[Trk::d0];
182 trackZ0 = perigee->parameters()[Trk::z0];
183 tracksigD0 = TMath::Sqrt(perigee->localErrorMatrix().covariance()[Trk::d0][Trk::d0]);
184 tracksigZ0 = TMath::Sqrt(perigee->localErrorMatrix().covariance()[Trk::z0][Trk::z0]);
185 delete perigee;
186 }
187 msg(MSG::VERBOSE) << "#BTAG# Track "
188 << " Eta= " << track->eta() << " Phi= " << track->phi() << " pT= " <<track->pt()
189 << " d0= " << trackD0
190 << " z0= " << trackZ0 << " sigd0= " << tracksigD0 << " sigz0: " << tracksigZ0 << endreq;
191
192
193 bool pass = true;
194 if(track->pt()<m_pTMin) {
195 pass = false;
196 failedCuts.set(pTMin);
197 }
198 if(fabs(trackD0)>m_d0Max) {
199 pass = false;
200 failedCuts.set(d0Max);
201 }
202 if(fabs(trackZ0*track->sinTh())>m_z0Max) {
203 pass = false;
204 failedCuts.set(z0Max);
205 }
206 if (tracksigD0>m_sigd0Max) {
207 pass = false;
208 failedCuts.set(sigd0Max);
209 }
210 if (tracksigZ0>m_sigz0Max) {
211 pass = false;
212 failedCuts.set(sigz0Max);
213 }
214 if(fabs(track->eta())>m_etaMax) {
215 pass = false;
216 failedCuts.set(etaMax);
217 }
218 if(m_useTrackSummaryInfo) {
219 const Trk::TrackSummary* trks = track->trackSummary();
220 if(trks) {
221 int nb = trks->get(Trk::numberOfBLayerHits);
222 if(nb<0) nb=0;
223 if(nb < m_nHitBLayer) {
224 failedCuts.set(nHitBLayer);
225 if(!m_useBLayerHitPrediction) {
226 pass = false;
227 failedCuts.set(deadBLayer);
228 } else {
229 if(m_inDetTestBLayerTool->expectHitInBLayer(track)) {
230 pass = false;
231 failedCuts.set(deadBLayer);
232 }
233 }
234 }
235 int np = trks->get(Trk::numberOfPixelHits);
236 if(np<0) np=0;
237 if(m_useDeadPixInfo) np += std::max(trks->get(Trk::numberOfPixelDeadSensors),0);
238 if(np < m_nHitPix) {
239 pass = false;
240 failedCuts.set(nHitPix);
241 }
242 int ns = trks->get(Trk::numberOfSCTHits);
243 if(ns<0) ns=0;
244 if(m_useDeadSctInfo) ns += std::max(trks->get(Trk::numberOfSCTDeadSensors),0);
245 if(ns < m_nHitSct) {
246 pass = false;
247 failedCuts.set(nHitSct);
248 }
249 if((np+ns) < m_nHitSi) {
250 pass = false;
251 failedCuts.set(nHitSi);
252 }
253 int nh = trks->get(Trk::numberOfTRTHits);
254 if(nh<0) nh=0;
255 if(nh < m_nHitTrt) {
256 pass = false;
257 failedCuts.set(nHitTrt);
258 }
259 int nhe = trks->get(Trk::numberOfTRTHighThresholdHits);
260 if(nhe<0) nhe=0;
261 if(nhe < m_nHitTrtHighE) {
262 pass = false;
263 failedCuts.set(nHitTrtHighE);
264 }
265 }
266 }
267
268 double chi2 = track->fitQuality()->chiSquared();
269 int ndf = track->fitQuality()->numberDoF();
270 double proba = 1.;
271 if(ndf>0 && chi2>=0.) {
272 Genfun::CumulativeChiSquare myCumulativeChiSquare(ndf);
273 proba = 1.-myCumulativeChiSquare(chi2);
274 }
275 if(chi2>m_fitChi2) {
276 pass = false;
277 failedCuts.set(fitChi2);
278 }
279 if(proba<m_fitProb) {
280 pass = false;
281 failedCuts.set(fitProb);
282 }
283 if(ndf>0) {
284 if(chi2/double(ndf)>m_fitChi2OnNdfMax) {
285 pass = false;
286 failedCuts.set(fitChi2OnNdfMax);
287 }
288 } else {
289 pass = false;
290 failedCuts.set(fitChi2OnNdfMax);
291 }
292
293
294 ATH_MSG_VERBOSE("#BTAG# passedCuts for track ");
295 for(int i=0;i<numCuts;i++) {
296 int passl = ~failedCuts[(Cuts)i];
297 if(passl) m_ntrc[i]++;
298 msg(MSG::VERBOSE) << passl;
299 }
300 msg(MSG::VERBOSE) << endreq;
301
302 m_passedCuts = ~failedCuts;
303 return pass;
304 }
305
306
307 StatusCode TrackSelector::selectAllTracks() {
308
309
310 const Rec::TrackParticleContainer* inputTracks(0);
311 StatusCode sc = evtStore()->retrieve(inputTracks, m_inputTrackCollection);
312 if (sc.isFailure()) {
313 ATH_MSG_ERROR("#BTAG# TrackParticleContainer " << m_inputTrackCollection << " not found.");
314 return sc;
315 }
316 ATH_MSG_VERBOSE("#BTAG# TrackParticleContainer " << m_inputTrackCollection << " found.");
317
318
319 Rec::TrackParticleContainer* outputTracks = new Rec::TrackParticleContainer();
320
321
322 Rec::TrackParticleContainer::const_iterator nextTrk(inputTracks->begin());
323 Rec::TrackParticleContainer::const_iterator lastTrk(inputTracks->end());
324 for (; nextTrk!=lastTrk; nextTrk++) {
325 if( this->selectTrack( (*nextTrk) ) ) {
326 ATH_MSG_VERBOSE("#BTAG# track selected");
327 outputTracks->push_back( new Rec::TrackParticle(**nextTrk) );
328 }
329 }
330
331
332 ATH_MSG_VERBOSE("#BTAG# recording to StoreGate: " << m_outputTrackCollection);
333 sc = evtStore()->record(outputTracks, m_outputTrackCollection, false);
334 if (sc.isFailure()) {
335 ATH_MSG_ERROR("#BTAG# unable to store TrackParticleContainer " << m_outputTrackCollection);
336 return sc;
337 }
338 ATH_MSG_VERBOSE("#BTAG# TrackParticleContainer " << m_outputTrackCollection << " stored.");
339
340 return StatusCode::SUCCESS;
341 }
342
343 }
| Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems |
|
This page was automatically generated by the
LXR engine.
|
|