001
002
003
004
005
006
007
008
009
010
011 #include "GaudiKernel/MsgStream.h"
012 #include "JetTagTools/SecVtxTag.h"
013 #include "TrkTrack/Track.h"
014 #include "TrkParametersBase/ParametersBase.h"
015 #include "TrkParameters/MeasuredPerigee.h"
016 #include "JetEvent/Jet.h"
017 #include "JetTagInfo/SecVtxInfo.h"
018 #include "CLHEP/GenericFunctions/CumulativeChiSquare.hh"
019 #include "CLHEP/Vector/LorentzVector.h"
020 #include "TrkVertexFitterInterfaces/IVertexFitter.h"
021 #include "Particle/TrackParticle.h"
022
023 #include "JetTagTools/LikelihoodTool.h"
024 #include "JetTagTools/HistoHelperRoot.h"
025 #include "VxVertex/VxCandidate.h"
026 #include "VxVertex/RecVertex.h"
027 #include "VxVertex/VxTrackAtVertex.h"
028 #include "ITrackToVertex/ITrackToVertex.h"
029
030 #include "GaudiKernel/ITHistSvc.h"
031 #include "JetTagInfo/TruthInfo.h"
032 #include "JetTagInfo/SvxSummary.h"
033
034 #include "JetTagEvent/ISvxAssociation.h"
035 #include "VxSecVertex/VxSecVertexInfo.h"
036 #include "VxSecVertex/VxSecVKalVertexInfo.h"
037 #include "TrkParticleBase/LinkToTrackParticleBase.h"
038 #include "TrkParticleBase/TrackParticleBase.h"
039
040 #include "JetTagEvent/TrackAssociation.h"
041
042 namespace Analysis
043 {
044 SecVtxTag::SecVtxTag(const std::string& t, const std::string& n, const IInterface* p)
045 : AlgTool(t,n,p),
046
047 m_runModus("analysis"),
048
049 m_histoHelperSig(0),
050 m_histoHelperBkg(0),
051
052 m_priVtx(0),
053
054 m_fitTool("Trk::FastVertexFitter"),
055 m_trackToVertexTool("Reco::TrackToVertex"),
056 m_likelihoodTool("Analysis::LikelihoodTool"),
057
058 m_useVKalTool(true),
059
060 m_useVariables(std::vector<std::string>()),
061
062 m_vertexFindingMode("BU"),
063
064 m_minpt(1.*GeV),
065 m_maxd0wrtPrimVtx(1.*mm),
066 m_maxz0wrtPrimVtx(1.5*mm),
067 m_min2DSigToPrimVtx(3.),
068
069 m_precisionHits(9),
070 m_PixelHits(2),
071 m_BLayerHits(1),
072
073 m_buMaxChi2OfTrack(4.),
074 m_minDistToPrimVtx(0.2),
075 m_secVxFinderName("InDetVKalVxInJetTool")
076 {
077 declareInterface<ITagTool>(this);
078 declareProperty("Runmodus", m_runModus);
079 declareProperty("useVariables", m_useVariables);
080
081 declareProperty("VertexFitterTool" , m_fitTool);
082 declareProperty("TrackToVertexTool" , m_trackToVertexTool);
083 declareProperty("LikelihoodTool" , m_likelihoodTool);
084
085 declareProperty("VertexFindingMode", m_vertexFindingMode);
086
087 declareProperty("minpt", m_minpt);
088 declareProperty("maxd0wrtPrimVtx", m_maxd0wrtPrimVtx);
089 declareProperty("maxz0wrtPrimVtx", m_maxz0wrtPrimVtx);
090 declareProperty("min2DSigToPrimVtx", m_min2DSigToPrimVtx);
091 declareProperty("PrecisionHits", m_precisionHits);
092 declareProperty("PixelHits", m_PixelHits);
093 declareProperty("BLayerHits", m_BLayerHits);
094
095 declareProperty("buMaxChi2OfTrack", m_buMaxChi2OfTrack);
096 declareProperty("minDistToPrimVtx", m_minDistToPrimVtx);
097
098 declareProperty("SecVxFinderName",m_secVxFinderName);
099
100 ITHistSvc* myHistoSvc(0);
101 MsgStream mLog(msgSvc(), name());
102
103 if( service( "THistSvc", myHistoSvc ).isSuccess() )
104 {
105 mLog << MSG::DEBUG << name() << ": HistoSvc loaded successfully." << endreq;
106 m_histoHelperSig = new HistoHelperRoot(myHistoSvc);
107 m_histoHelperBkg = new HistoHelperRoot(myHistoSvc);
108 m_histoHelperSig->setCheckOverflows(true);
109 m_histoHelperBkg->setCheckOverflows(true);
110 }
111 else mLog << MSG::ERROR << name() << ": HistoSvc could NOT be loaded." << endreq;
112 }
113
114 SecVtxTag::~SecVtxTag()
115 {
116 delete m_histoHelperSig;
117 delete m_histoHelperBkg;
118 }
119
120 StatusCode SecVtxTag::initialize()
121 {
122 MsgStream mlog(msgSvc(), name());
123 m_printParameterSettings();
124
125
126
127
128 if (!m_useVKalTool)
129 {
130 if ( m_fitTool.retrieve().isFailure() ) {
131 mlog << MSG::FATAL << "failed to retrieve tool " << m_fitTool << endreq;
132 return StatusCode::FAILURE;
133 } else {
134 mlog << MSG::INFO << "retrieved tool " << m_fitTool << endreq;
135 }
136
137
138 if ( m_trackToVertexTool.retrieve().isFailure() ) {
139 mlog << MSG::FATAL << "failed to retrieve tool " << m_trackToVertexTool << endreq;
140 return StatusCode::FAILURE;
141 } else {
142 mlog << MSG::INFO << "retrieved tool " << m_trackToVertexTool << endreq;
143 }
144 }
145
146
147
148
149 if (m_runModus == "analysis")
150 {
151
152 if ( m_likelihoodTool.retrieve().isFailure() ) {
153 mlog << MSG::FATAL << "failed to retrieve tool " << m_likelihoodTool << endreq;
154 return StatusCode::FAILURE;
155 } else {
156 mlog << MSG::INFO << "retrieved tool " << m_likelihoodTool << endreq;
157 }
158
159
160 m_likelihoodTool->readInHistosFromFile("SigCali/");
161 m_likelihoodTool->readInHistosFromFile("BkgCali/");
162 m_likelihoodTool->printStatus();
163 }
164
165
166
167
168 if (m_runModus=="reference")
169 {
170 m_histoHelperSig->bookHisto("/SigCali/vtxmassBU", "", 40, 0., 5.);
171 m_histoHelperBkg->bookHisto("/BkgCali/vtxmassBU", "", 40, 0., 5.);
172 m_histoHelperSig->bookHisto("/SigCali/vtxmultBU", "", 10, 0.1, 10.1);
173 m_histoHelperBkg->bookHisto("/BkgCali/vtxmultBU", "", 10, 0.1, 10.1);
174 m_histoHelperSig->bookHisto("/SigCali/vtxdistBU", "", 80, 0.1, 40.1);
175 m_histoHelperBkg->bookHisto("/BkgCali/vtxdistBU", "", 80, 0.1, 40.1);
176 m_histoHelperSig->bookHisto("/SigCali/vtxenergyFractionBU", "", 40, 0., 1.);
177 m_histoHelperBkg->bookHisto("/BkgCali/vtxenergyFractionBU", "", 40, 0., 1.);
178 m_histoHelperSig->bookHisto("/SigCali/vtxNGTVtxBU", "", 21, 0.1, 21.1);
179 m_histoHelperBkg->bookHisto("/BkgCali/vtxNGTVtxBU", "", 21, 0.1, 21.1);
180 }
181
182 mlog << MSG::DEBUG << name() << ": Initializing successfully." << endreq;
183 return StatusCode::SUCCESS;
184 }
185
186 StatusCode SecVtxTag::finalize()
187 {
188 MsgStream mlog(msgSvc(), name());
189
190 return StatusCode::SUCCESS;
191 }
192
193 void SecVtxTag::tagJet(Jet& jetToTag)
194 {
195 MsgStream mlog(msgSvc(), name());
196
197 mlog << MSG::DEBUG << "in tagJet()" << endreq;
198
199 std::string instanceName(name());
200 SecVtxInfo * secVtxInfo = new SecVtxInfo(instanceName.erase(0,8));
201 jetToTag.addInfo(secVtxInfo);
202
203 TrackVec tracksInSV;
204 const Trk::VxCandidate* myVxCandidate(0);
205 Trk::RecVertex VertInJet;
206 double Prob(-1.);
207 double RPhiDistance(-1.);
208 double VtxMass(-1000.);
209 double energyFraction(-1.);
210 int NGood2TrackVertices(-1);
211 int NTracksInJet(-1);
212 int NTracksInSV(-1);
213 double distance(-1.);
214
215 mlog << MSG::VERBOSE << "PVX x = " << m_priVtx->recVertex().position().x() << " y = " << m_priVtx->recVertex().position().y() << " z = " << m_priVtx->recVertex().position().z() << endreq;
216
217
218 if (!m_useVKalTool)
219 {
220
221
222 std::vector<const Rec::TrackParticle*> myJetToTag;
223
224 std::vector<const Rec::TrackParticle*>* trackVector = jetToTag.getAssociation<TrackAssociation>("Tracks")->tracks();
225 for (std::vector<const Rec::TrackParticle*>::iterator jetItr = trackVector->begin(); jetItr != trackVector->end() ; ++jetItr )
226 {
227 if (m_preselect(*jetItr)==true) myJetToTag.push_back(*jetItr);
228 }
229 delete trackVector; trackVector=0;
230
231 mlog << MSG::VERBOSE << myJetToTag.size() << " tracks after preselect for sec vtx fit." << endreq;
232 if (myJetToTag.size()<2)
233 {
234 mlog << MSG::VERBOSE << "Less than two tracks in the jet. Returning ..." << endreq;
235 secVtxInfo->setFitType(SecVtxInfo::NoFit);
236 return;
237 }
238
239
240
241
242 if (m_vertexFindingMode == "BU")
243 {
244 myVxCandidate=m_fitBuildUpVertex(myJetToTag, tracksInSV);
245 }
246 else if (m_vertexFindingMode == "TD")
247 {
248 myVxCandidate=m_fitTearDownVertex(myJetToTag, tracksInSV);
249 }
250 else
251 {
252 mlog << MSG::DEBUG << "Wrong vertex fitting method chosen ..." << m_vertexFindingMode << endreq;
253 return;
254 }
255
256 if (myVxCandidate!=0)
257 {
258 const Trk::RecVertex& vertexInJet = myVxCandidate->recVertex();
259 int ndof = vertexInJet.fitQuality().numberDoF();
260 double chi2 = vertexInJet.fitQuality().chiSquared();
261 if(ndof>0 && chi2>=0.)
262 {
263 Genfun::CumulativeChiSquare myCumulativeChiSquare(ndof);
264 Prob = 1.-myCumulativeChiSquare(chi2);
265 }
266
267
268 if (Prob >= 0.)
269 {
270 if (m_vertexFindingMode == "BU")
271 secVtxInfo->setFitType(SecVtxInfo::BuildUp);
272 else if (m_vertexFindingMode == "TD")
273 secVtxInfo->setFitType(SecVtxInfo::TearDown);
274
275 mlog << MSG::VERBOSE << m_vertexFindingMode << " Vtx fitted with " << tracksInSV.size() << " tracks. Probability = " << Prob << endreq;
276 VtxMass=m_massOfVertex(tracksInSV);
277 mlog << MSG::VERBOSE << m_vertexFindingMode << " Vtx mass = " << VtxMass/1000. << " GeV" << endreq;
278 if (VtxMass<0.)
279 VtxMass=0.;
280
281 double xMom(0.);
282 double yMom(0.);
283 double zMom(0.);
284 double Energy(0.);
285 for (TrackIter i=tracksInSV.begin(); i!=tracksInSV.end(); i++)
286 {
287 double x((*i)->measuredPerigee()->momentum()[Trk::px]);
288 double y((*i)->measuredPerigee()->momentum()[Trk::py]);
289 double z((*i)->measuredPerigee()->momentum()[Trk::pz]);
290 xMom += x;
291 yMom += y;
292 zMom += z;
293 Energy += sqrt(x*x + y*y + z*z);
294 }
295 energyFraction = Energy/jetToTag.e();
296
297 RPhiDistance = sqrt(
298 pow( m_priVtx->recVertex().position()[0] - vertexInJet.position()[0] , 2 ) +
299 pow( m_priVtx->recVertex().position()[1] - vertexInJet.position()[1] , 2 ) );
300
301 distance = sqrt(
302 pow( RPhiDistance , 2 ) +
303 pow( m_priVtx->recVertex().position()[2] - vertexInJet.position()[2] , 2 ) );
304 }
305
306 else if (Prob < 0.)
307 {
308 mlog << MSG::DEBUG << "Fitted secondary vertex probability < 0!" << endreq;
309 secVtxInfo->setFitType(SecVtxInfo::NoFit);
310 tracksInSV.clear();
311 return;
312 }
313 }
314 else
315 {
316 mlog << MSG::DEBUG << "Fitted secondary vertex pointer is zero!" << endreq;
317 secVtxInfo->setFitType(SecVtxInfo::NoFit);
318 tracksInSV.clear();
319 return;
320 }
321 }
322 else
323 {
324
325
326
327 std::vector<const Trk::TrackParticleBase*> TrkFromV0;
328
329 bool vertexIsThere(true);
330
331 const ISvxAssociation* newSvxAssociation=jetToTag.getAssociation<ISvxAssociation>(m_secVxFinderName);
332
333 if (newSvxAssociation==0) {
334 mlog << MSG::DEBUG << " No VKalVrt vertex found, attached to the Jet as association. Not going on with tagging..." << endreq;
335 return;
336 }
337
338 const Trk::VxSecVertexInfo* myVertexInfo=newSvxAssociation->vertexInfo();
339 if (myVertexInfo==0) {
340 vertexIsThere=false;
341 }
342
343 if (vertexIsThere)
344 {
345
346 const Trk::VxSecVKalVertexInfo* myVKalVertexInfo=dynamic_cast<const Trk::VxSecVKalVertexInfo*>(myVertexInfo);
347
348 if (myVKalVertexInfo==0) {
349 mlog << MSG::WARNING << " Could not cast to a VKalVrt Info object " << endreq;
350 return;
351 }
352
353 const std::vector<Trk::VxCandidate*> & myVertices=myVKalVertexInfo->vertices();
354
355 if (myVertices.size()>0) {
356
357 VtxMass=myVKalVertexInfo->mass();
358 energyFraction=myVKalVertexInfo->energyFraction();
359 NGood2TrackVertices=(long int)myVKalVertexInfo->n2trackvertices();
360 NTracksInJet=-1;
361
362 NTracksInSV=0;
363
364 const std::vector<Trk::VxCandidate*>::const_iterator verticesBegin=myVertices.begin();
365 const std::vector<Trk::VxCandidate*>::const_iterator verticesEnd=myVertices.end();
366
367 for (std::vector<Trk::VxCandidate*>::const_iterator verticesIter=verticesBegin;
368 verticesIter!=verticesEnd;verticesIter++) {
369 if ((*verticesIter)==0)
370 {
371 mlog << MSG::ERROR << " Secondary vertex from InDetVKalVxInJetFinder has zero pointer. Skipping... " << endreq;
372 continue;
373 }
374 const std::vector<Trk::VxTrackAtVertex*>* myTracks=(*verticesIter)->vxTrackAtVertex();
375 if (myTracks!=0) {
376 NTracksInSV+=myTracks->size();
377 }
378 }
379 }
380
381 if (myVertices[0]!=0) {
382 const Trk::RecVertex & myVertex=myVertices[0]->recVertex();
383 VertInJet = myVertex;
384
385 if (m_priVtx)
386 {
387 distance=(myVertex.position()-m_priVtx->recVertex().position()).mag();
388 RPhiDistance=(myVertex.position()-m_priVtx->recVertex().position()).perp();
389 }
390 else
391 {
392 mlog << MSG::WARNING << "Tagging requested, but no primary vertex supplied." << endreq;
393 }
394 }
395 else
396 {
397 mlog << MSG::VERBOSE << "No vertex. Cannot calculate distances." << endreq;
398 }
399
400 }
401 }
402
403 mlog << MSG::VERBOSE << "SecVtx Position x = " << VertInJet.position().x() << " y = " << VertInJet.position().y() << " z = " << VertInJet.position().z() << endreq;
404 mlog << MSG::VERBOSE << "(norm)Distance: " <<distance<< endreq;
405 mlog << MSG::VERBOSE << "RPhi distance to PV: " << RPhiDistance << endreq;
406 mlog << MSG::VERBOSE << "N tracks in jet from VKalVrt "<<NTracksInJet<< " in SVX "<<NTracksInSV<<endreq;
407 mlog << MSG::VERBOSE << "Svx Mass = "<< VtxMass/1000. << " Svx E frac = " << energyFraction << " NGood2TrackVertices = " << NGood2TrackVertices << endreq;
408
409
410 if (m_useVKalTool) secVtxInfo->setSecVtx(VertInJet, Prob, tracksInSV);
411 else secVtxInfo->setSecVtx(myVxCandidate->recVertex(), Prob, tracksInSV);
412 secVtxInfo->setNumSelTracksForFit(NTracksInJet);
413 secVtxInfo->setMult(NTracksInSV);
414 secVtxInfo->setMass(VtxMass/1000.);
415 secVtxInfo->setRPhiDist(RPhiDistance);
416 secVtxInfo->setDist(distance);
417 secVtxInfo->setEnergyFraction(energyFraction);
418 secVtxInfo->setNumberOfG2TrackVertices(NGood2TrackVertices);
419
420
421
422 if (m_runModus=="reference")
423 {
424 if (jetToTag.pt() > 15.*GeV && fabs(jetToTag.eta()) < 2.5)
425 {
426
427 const TruthInfo* truthInfo = jetToTag.tagInfo<TruthInfo>("TruthInfo");
428 if (truthInfo != 0)
429 {
430 HistoHelperRoot* histoHelper(0);
431 std::string hDir;
432 double rmin(0.);
433 if (truthInfo->jetTruthLabel() == "N/A")
434 {
435
436 rmin = truthInfo->deltaRMinTo("C");
437 if (rmin < truthInfo->deltaRMinTo("B")) rmin = truthInfo->deltaRMinTo("B");
438 if (rmin > 0.8) {
439 histoHelper = m_histoHelperBkg;
440 hDir = "/BkgCali/";
441 }
442 }
443 if (truthInfo->jetTruthLabel()=="B") {
444 histoHelper = m_histoHelperSig;
445 hDir = "/SigCali/";
446 }
447 if (histoHelper != 0)
448 {
449 histoHelper->fillHisto(hDir+"vtxmassBU", VtxMass/1000.);
450 histoHelper->fillHisto(hDir+"vtxdistBU", distance);
451 histoHelper->fillHisto(hDir+"vtxmultBU", (double)tracksInSV.size());
452 histoHelper->fillHisto(hDir+"vtxenergyFractionBU", energyFraction);
453 histoHelper->fillHisto(hDir+"vtxNGTVtxBU", (double)NGood2TrackVertices);
454 }
455 } else
456 {
457 mlog << MSG::WARNING << "SecVtxTag TruthLabel == 0!!!" << endreq;
458 }
459 }
460 }
461
462 if (m_likelihoodTool!=0)
463 {
464
465 m_likelihoodTool->setLhVariableValue("vtxmass"+m_vertexFindingMode, VtxMass/1000.);
466 m_likelihoodTool->setLhVariableValue("vtxdist"+m_vertexFindingMode, distance);
467 m_likelihoodTool->setLhVariableValue("vtxmult"+m_vertexFindingMode, (double)NTracksInSV);
468 m_likelihoodTool->setLhVariableValue("vtxenergyFraction"+m_vertexFindingMode, energyFraction);
469 m_likelihoodTool->setLhVariableValue("vtxNGTVtx"+m_vertexFindingMode, (double)NGood2TrackVertices);
470 secVtxInfo->setTagLikelihood(m_likelihoodTool->calculateLikelihood());
471 secVtxInfo->setWeight(m_likelihoodTool->calculateWeight());
472 }
473
474
475 secVtxInfo->makeValid();
476 mlog << MSG::DEBUG << "finished" << endreq;
477 return;
478 }
479
480 void SecVtxTag::finalizeHistos()
481 {
482 if (m_runModus=="reference")
483 {
484 m_histoHelperSig->normalizeHistos();
485 m_histoHelperBkg->normalizeHistos();
486
487
488
489
490
491 }
492 return;
493 }
494
495
496
497
498
499 Hep3Vector SecVtxTag::m_jetMomentum(const TrackVec& tracksInJet)
500 {
501 Hep3Vector jetDir;
502 double px=0., py=0., pz=0.;
503
504
505 for (TrackIter i=tracksInJet.begin(); i!=tracksInJet.end(); ++i)
506 {
507 px += (*i)->measuredPerigee()->momentum()[Trk::px];
508 py += (*i)->measuredPerigee()->momentum()[Trk::py];
509 pz += (*i)->measuredPerigee()->momentum()[Trk::pz];
510 }
511 jetDir.set(px,py,pz);
512 return jetDir;
513 }
514
515
516
517 bool SecVtxTag::m_preselect(const Rec::TrackParticle *trkPrt)
518 {
519 MsgStream mLog(msgSvc(), name());
520 const Trk::MeasuredPerigee* perigee(0);
521 const Trk::MeasuredPerigee* perigeeAtGlobal = trkPrt->measuredPerigee();
522 const Trk::MeasuredPerigee* perigeeAtVertex = m_trackToVertexTool->perigeeAtVertex(*trkPrt, m_priVtx->recVertex().position());
523
524 double d0wrtPriVtx( perigeeAtVertex->parameters()[Trk::d0] );
525
526
527
528
529
530
531
532
533
534
535
536
537 if (d0wrtPriVtx == 0.)
538 {
539 return false;
540 }
541
542 bool print(false);
543
544
545
546
547
548
549 if (trkPrt->pt()<m_minpt)
550 {
551 mLog << MSG::VERBOSE << "Skipping track because of pt " << trkPrt->pt() << " < " << m_minpt << endreq;
552 return false;
553 }
554
555 perigee=perigeeAtVertex;
556 if (fabs(perigee->parameters()[Trk::d0])>m_maxd0wrtPrimVtx)
557 {
558 mLog << MSG::VERBOSE << "Skipping track because of d0 " << fabs(perigee->parameters()[Trk::d0]) << " > " << m_maxd0wrtPrimVtx << endreq;
559 delete perigeeAtVertex;
560 return false;
561 }
562 if (fabs(perigee->parameters()[Trk::z0])>m_maxz0wrtPrimVtx)
563 {
564 mLog << MSG::VERBOSE << "Skipping track because of z0 " << fabs(perigee->parameters()[Trk::z0]) << " > " << m_maxz0wrtPrimVtx << endreq;
565 delete perigeeAtVertex;
566 return false;
567 }
568
569
570 perigee=perigeeAtGlobal;
571 double signif2D(perigee->parameters()[Trk::d0]/perigee->localErrorMatrix().error(Trk::d0));
572 if ( fabs(signif2D) < m_min2DSigToPrimVtx )
573 {
574 mLog << MSG::VERBOSE << "Skipping track because of d0/sigd0 " << signif2D << " < " << m_min2DSigToPrimVtx << endreq;
575 delete perigeeAtVertex;
576 return false;
577 }
578 delete perigeeAtVertex;
579
580 int SCTHits = trkPrt->trackSummary()->get
581 (Trk::numberOfSCTHits);
582 int PixelHits = trkPrt->trackSummary()->get
583 (Trk::numberOfPixelHits);
584 if (print)
585 {
586 std::cout <<"SCTHits: " << SCTHits << std::endl;
587 std::cout <<"PixelHits: " << PixelHits << std::endl;
588 }
589 if ((SCTHits+PixelHits)<m_precisionHits)
590 return false;
591 if (PixelHits<m_PixelHits)
592 return false;
593
594 int BLayerHits= trkPrt->trackSummary()->get
595 (Trk::numberOfBLayerHits);
596 if (print)
597 std::cout <<"BLayerHits: " << BLayerHits << std::endl;
598 if (BLayerHits<m_BLayerHits)
599 return false;
600 return true;
601 }
602
603 const Trk::VxCandidate* SecVtxTag::m_fitTearDownVertex(TrackVec tracksInJet, TrackVec& tracksInSV)
604 {
605 MsgStream mlog(msgSvc(), name());
606
607 mlog << MSG::DEBUG << "starting m_fitTearDown" << endreq;
608
609
610 std::vector<const Trk::ParametersBase*> vecOfMeasuredPerigee;
611 for (TrackVec::iterator trkPrtItr = tracksInJet.begin();
612 trkPrtItr != tracksInJet.end(); ++ trkPrtItr)
613 {
614 vecOfMeasuredPerigee.push_back((*trkPrtItr)->measuredPerigee());
615 }
616
617 const Trk::VxCandidate* myVxCandidate(0);
618 if (vecOfMeasuredPerigee.size()>1)
619 {
620 Hep3Vector startpoint(0.,0.,0.);
621
622 std::vector<int> indexOfSortedChi2;
623 myVxCandidate = m_fitTool->fit(vecOfMeasuredPerigee, Trk::Vertex(startpoint));
624 if (myVxCandidate==0) return 0;
625 m_sortTracksInChi2(indexOfSortedChi2, myVxCandidate);
626
627 startpoint=myVxCandidate->recVertex().position();
628
629
630
631 while ( (*(myVxCandidate->vxTrackAtVertex()))[*(indexOfSortedChi2.end()-1)]->trackQuality().chiSquared() > 5. &&
632 vecOfMeasuredPerigee.size() > 2 )
633 {
634 delete myVxCandidate;
635
636 tracksInJet.erase( tracksInJet.begin()+(*(indexOfSortedChi2.end()-1)) );
637 vecOfMeasuredPerigee.erase( vecOfMeasuredPerigee.begin()+(*(indexOfSortedChi2.end()-1)) );
638
639 myVxCandidate = m_fitTool->fit(vecOfMeasuredPerigee, Trk::Vertex(startpoint));
640 if (myVxCandidate==0) return 0;
641 m_sortTracksInChi2(indexOfSortedChi2, myVxCandidate);
642 startpoint=myVxCandidate->recVertex().position();
643 }
644 for (unsigned int i = 0; i < tracksInJet.size(); ++i)
645 {
646 tracksInSV.push_back(tracksInJet[indexOfSortedChi2[i]]);
647 }
648 }
649 return myVxCandidate;
650 }
651
652
653
654 const Trk::VxCandidate* SecVtxTag::m_fitBuildUpVertex(TrackVec tracksInJet, TrackVec& tracksInSV)
655 {
656 MsgStream mlog(msgSvc(), name());
657 mlog << MSG::DEBUG << "starting m_fitBuildUp" << endreq;
658
659 std::vector<const Trk::ParametersBase*> trktracksInBUSV;
660 double fitProb(-1.);
661
662 std::vector<const Trk::ParametersBase*> vecOfMeasuredPerigee;
663 for (TrackVec::iterator trkPrtItr = tracksInJet.begin();
664 trkPrtItr != tracksInJet.end(); ++trkPrtItr)
665 {
666 vecOfMeasuredPerigee.push_back((*trkPrtItr)->measuredPerigee());
667 }
668
669 Hep3Vector startpoint(0.,0.,0.);
670 const Trk::VxCandidate* myVxCandidate(0);
671
672 if (vecOfMeasuredPerigee.size()>1)
673 {
674
675 for (std::vector<const Trk::ParametersBase*>::iterator i=vecOfMeasuredPerigee.begin(); i!=vecOfMeasuredPerigee.end()-1; i++)
676 {
677 for (std::vector<const Trk::ParametersBase*>::iterator j=i+1; j!=vecOfMeasuredPerigee.end(); j++)
678 {
679 std::vector<const Trk::ParametersBase*> tracksToFit;
680 tracksToFit.push_back((*i));
681 tracksToFit.push_back((*j));
682 const Trk::VxCandidate* tmpVxCandidate(m_fitTool->fit(tracksToFit, Trk::Vertex(startpoint)));
683 if (tmpVxCandidate==0)
684 {
685 if (myVxCandidate!=0) delete myVxCandidate;
686 return 0;
687 }
688 const Trk::RecVertex& tempVtx(tmpVxCandidate->recVertex());
689 int ndof = tempVtx.fitQuality().numberDoF();
690 double chi2 = tempVtx.fitQuality().chiSquared();
691 double tempProb(-1.);
692 if(ndof>0 && chi2>=0.)
693 {
694 Genfun::CumulativeChiSquare myCumulativeChiSquare(ndof);
695 tempProb = 1.-myCumulativeChiSquare(chi2);
696 }
697
698 double tempDist(sqrt(pow(tempVtx.position()[0]-m_priVtx->recVertex().position()[0],2)+
699 pow(tempVtx.position()[1]-m_priVtx->recVertex().position()[1],2)+
700 pow(tempVtx.position()[2]-m_priVtx->recVertex().position()[2],2)
701 ));
702 if ((tempProb>fitProb) && (tempDist>m_minDistToPrimVtx))
703 {
704 if (myVxCandidate!=0)
705 delete myVxCandidate;
706 fitProb=tempProb;
707 trktracksInBUSV.clear();
708 trktracksInBUSV.push_back((*i));
709 trktracksInBUSV.push_back((*j));
710 myVxCandidate=tmpVxCandidate;
711 }
712 else
713 delete tmpVxCandidate;
714 }
715 }
716
717 if (trktracksInBUSV.size() == 0)
718 return 0;
719
720
721 for (std::vector<const Trk::ParametersBase*>::iterator i=vecOfMeasuredPerigee.begin(); i!=vecOfMeasuredPerigee.end(); i++)
722 {
723 if ((*i)==trktracksInBUSV[0] || (*i)==trktracksInBUSV[1])
724 continue;
725 trktracksInBUSV.push_back((*i));
726 startpoint=myVxCandidate->recVertex().position();
727 const Trk::VxCandidate* tmpVxCandidate(m_fitTool->fit(trktracksInBUSV, Trk::Vertex(startpoint)));
728 if (tmpVxCandidate==0)
729 {
730 if (myVxCandidate!=0) delete myVxCandidate;
731 return 0;
732 }
733 const Trk::RecVertex& tempVtx(tmpVxCandidate->recVertex());
734 int ndof = tempVtx.fitQuality().numberDoF();
735 double tchi2 = tempVtx.fitQuality().chiSquared();
736 double tempProb(-1.);
737 if(ndof>0 && tchi2>=0.)
738 {
739 Genfun::CumulativeChiSquare myCumulativeChiSquare(ndof);
740 tempProb = 1.-myCumulativeChiSquare(tchi2);
741 }
742 double chi2( (*(tmpVxCandidate->vxTrackAtVertex()))[trktracksInBUSV.size()-1]->trackQuality().chiSquared() );
743 if ( (chi2 > m_buMaxChi2OfTrack) && (trktracksInBUSV.size() > 1) )
744 {
745 trktracksInBUSV.pop_back();
746 delete tmpVxCandidate;
747 }
748 else
749 {
750 delete myVxCandidate;
751 myVxCandidate=tmpVxCandidate;
752 fitProb=tempProb;
753 }
754 }
755
756 for (TrackVec::iterator trkPrtItr = tracksInJet.begin();
757 trkPrtItr != tracksInJet.end(); ++ trkPrtItr)
758 {
759 for (std::vector<const Trk::ParametersBase*>::iterator i=trktracksInBUSV.begin(); i!=trktracksInBUSV.end(); i++)
760 {
761 if ( (*trkPrtItr)->measuredPerigee() == (*i) )
762 tracksInSV.push_back((*trkPrtItr));
763 }
764 }
765 }
766
767 return myVxCandidate;
768 }
769
770 double SecVtxTag::m_massOfVertex(const TrackVec& TPInVtx)
771 {
772 if (TPInVtx.size()<2)
773 return -1.;
774 double xMom(0.);
775 double yMom(0.);
776 double zMom(0.);
777 double Energy(0.);
778 for (TrackIter i=TPInVtx.begin(); i!=TPInVtx.end(); i++)
779 {
780 double x((*i)->measuredPerigee()->momentum()[Trk::px]);
781 double y((*i)->measuredPerigee()->momentum()[Trk::py]);
782 double z((*i)->measuredPerigee()->momentum()[Trk::pz]);
783 xMom += x;
784 yMom += y;
785 zMom += z;
786 Energy += sqrt(x*x + y*y + z*z);
787 }
788
789 HepLorentzVector final(xMom, yMom, zMom, Energy);
790 return final.m();
791 }
792
793 void SecVtxTag::m_sortTracksInChi2(std::vector<int>& indexOfSortedChi2, const Trk::VxCandidate* myVxCandidate)
794 {
795
796
797
798 std::map<double, int> mapOfChi2;
799 for (unsigned int i=0; i < myVxCandidate->vxTrackAtVertex()->size(); ++i)
800 {
801
802 double chi2PerTrack = (*(myVxCandidate->vxTrackAtVertex()))[i]->trackQuality().chiSquared();
803 mapOfChi2.insert( std::map<double, int>::value_type( chi2PerTrack, i ) );
804 }
805 indexOfSortedChi2.clear();
806 std::map<double, int>::const_iterator mItr = mapOfChi2.begin();
807 for ( ; mItr != mapOfChi2.end() ; ++mItr )
808 {
809 indexOfSortedChi2.push_back( (*mItr).second );
810 }
811 }
812
813 void SecVtxTag::m_printParameterSettings()
814 {
815 MsgStream log(msgSvc(), name());
816 log << MSG::INFO << name() << " initialize(): Parametersettings " << endreq;
817 log << MSG::INFO << "I am in " << m_runModus << " modus." << endreq;
818 log << MSG::INFO << "I am using the following variables:" << endreq;
819 for (std::vector<std::string>::iterator i = m_useVariables.begin();
820 i!=m_useVariables.end() ; ++i)
821 {
822 log << MSG::INFO << (*i) << endreq;
823 }
824 log << MSG::INFO << "Track quality cuts:" << endreq;
825 log << MSG::INFO << "pt > " << m_minpt << "\t MeV " << endreq;
826 log << MSG::INFO << "d0wrtPrimVtx < " << m_maxd0wrtPrimVtx << "\t mm " << endreq;
827 log << MSG::INFO << "z0wrtPrimVtx < " << m_maxz0wrtPrimVtx << "\t mm " << endreq;
828 log << MSG::INFO << "PrecisionHits >= " << m_precisionHits << endreq;
829 log << MSG::INFO << "PixelHits >= " << m_PixelHits << endreq;
830 log << MSG::INFO << "BLayerHits >= " << m_BLayerHits << endreq;
831 }
832
833 }
834
| 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.
|
|