001
002
003
004 #include "JetTagTools/SVTag.h"
005
006 #include "CLHEP/Vector/ThreeVector.h"
007 #include "JetEvent/Jet.h"
008 #include "GaudiKernel/IToolSvc.h"
009 #include "JetTagInfo/TruthInfo.h"
010 #include "JetTagInfo/SVInfoBase.h"
011 #include "JetTagInfo/SVInfoPlus.h"
012 #include "Navigation/NavigationToken.h"
013 #include "JetTagTools/NewLikelihoodTool.h"
014 #include "GaudiKernel/ITHistSvc.h"
015 #include "JetTagTools/HistoHelperRoot.h"
016 #include "JetTagTools/LikelihoodComponents.h"
017
018 #include "JetTagEvent/ISvxAssociation.h"
019 #include "VxSecVertex/VxSecVertexInfo.h"
020 #include "VxSecVertex/VxSecVKalVertexInfo.h"
021 #include "TrkParticleBase/LinkToTrackParticleBase.h"
022 #include "TrkParticleBase/TrackParticleBase.h"
023
024 #include "VxVertex/RecVertex.h"
025 #include "VxVertex/VxTrackAtVertex.h"
026
027 namespace Analysis
028 {
029
030 SVTag::SVTag(const std::string& t, const std::string& n, const IInterface* p)
031 : AthAlgTool(t,n,p),
032 m_likelihoodTool("Analysis::NewLikelihoodTool"),
033 m_histoHelper(0),
034 m_secVxFinderName("InDetVKalVxInJetTool")
035 {
036 declareInterface<ITagTool>(this);
037 declareProperty("Runmodus", m_runModus= "reference");
038 declareProperty("referenceType", m_refType = "ALL");
039 declareProperty("SVAlgType", m_SVmode);
040 declareProperty("jetPtMinRef", m_pTjetmin = 15.*GeV);
041 declareProperty("LikelihoodTool", m_likelihoodTool);
042 declareProperty("checkOverflows", m_checkOverflows = false);
043 declareProperty("purificationDeltaR", m_purificationDeltaR = 0.8);
044 declareProperty("UseBinInterpol", m_UseBinInterpol = false);
045
046 declareProperty("writeInfoBase", m_writeInfoBase = true);
047 declareProperty("infoPlusName", m_infoPlusName = "SVInfoPlus");
048 declareProperty("originalTPCollectionName", m_originalTPCollectionName = "TrackParticleCandidate");
049 declareProperty("jetCollectionList", m_jetCollectionList);
050 declareProperty("jetWithInfoPlus", m_jetWithInfoPlus);
051 m_jetWithInfoPlus.push_back("Cone4H1Tower");
052 declareProperty("useForcedCalibration", m_doForcedCalib = false);
053 declareProperty("ForcedCalibrationName", m_ForcedCalibName = "Cone4H1Tower");
054
055 declareProperty("SecVxFinderName",m_secVxFinderName);
056
057 declareProperty("UseDRJetPvSv", m_useDRJPVSV = true);
058
059 }
060
061 SVTag::~SVTag() {
062 delete m_histoHelper;
063 }
064
065 StatusCode SVTag::initialize() {
066 m_printParameterSettings();
067
068
069 m_hypotheses.push_back("B");
070 m_hypotheses.push_back("U");
071
072
073 m_SaveSVTrackInfo = false;
074
075
076 ITHistSvc* myHistoSvc;
077 if( service( "THistSvc", myHistoSvc ).isSuccess() ) {
078 ATH_MSG_DEBUG("#BTAG#" << name() << ": HistoSvc loaded successfully.");
079 m_histoHelper = new HistoHelperRoot(myHistoSvc);
080 m_histoHelper->setCheckOverflows(m_checkOverflows);
081 } else ATH_MSG_ERROR("#BTAG#" << name() << ": HistoSvc could NOT bo loaded.");
082
083
084 if( m_runModus == "analysis" && m_SVmode != "SV0" ) {
085 if ( m_likelihoodTool.retrieve().isFailure() ) {
086 ATH_MSG_FATAL("#BTAG# Failed to retrieve tool " << m_likelihoodTool);
087 return StatusCode::FAILURE;
088 } else {
089 ATH_MSG_INFO("#BTAG# Retrieved tool " << m_likelihoodTool);
090 }
091 m_likelihoodTool->defineHypotheses(m_hypotheses);
092
093 for(unsigned int ih=0;ih<m_hypotheses.size();ih++) {
094 std::string hName;
095 if(m_SVmode=="SV1") {
096 hName = m_hypotheses[ih]+"/N2T";
097 m_likelihoodTool->defineHistogram(hName);
098 hName = m_hypotheses[ih]+"/BidimME";
099 m_likelihoodTool->defineHistogram(hName);
100 if(m_useDRJPVSV) {
101 hName = m_hypotheses[ih]+"/DRJPVSV";
102 m_likelihoodTool->defineHistogram(hName);
103 }
104 }
105 if(m_SVmode=="SV2") {
106 hName = m_hypotheses[ih]+"/TridimMEN2T";
107 m_likelihoodTool->defineHistogram(hName);
108 }
109 }
110
111 for(unsigned int ih=0;ih<m_hypotheses.size();ih++) {
112 std::string hName = m_hypotheses[ih]+"/N2TEff"+m_SVmode;
113 m_likelihoodTool->defineHistogram(hName);
114 hName = m_hypotheses[ih]+"/N2TNorm"+m_SVmode;
115 m_likelihoodTool->defineHistogram(hName);
116 }
117 m_likelihoodTool->printStatus();
118 }
119
120 if (m_runModus == "analysis") m_SaveSVTrackInfo = true;
121
122
123
124
125 if (m_runModus=="reference" && m_SVmode!="SV0") {
126
127
128
129 ATH_MSG_VERBOSE("#BTAG# Booking histos...");
130 std::vector<double> xb;
131 double xbi[10] = {1., 2., 3., 4., 5., 6., 8., 10., 20., 50.};
132 for(uint ijc=0;ijc<m_jetCollectionList.size();ijc++) {
133
134 for(uint ih=0;ih<m_hypotheses.size();ih++) {
135
136 if (m_SVmode == "SV1") {
137 std::string hDir = "/RefFileSV1"+m_jetCollectionList[ijc]+"/"+m_hypotheses[ih]+"/";
138 m_histoHelper->bookHisto(hDir+"N2T", "Number of Good Two Track Vertices",9,xbi);
139 m_histoHelper->bookHisto(hDir+"N2TEffSV1", "Number of Good Two Track Vertices",9,xbi);
140 m_histoHelper->bookHisto(hDir+"N2TNormSV1", "Number of Good Two Track Vertices",30,0.,30.);
141 m_histoHelper->bookHisto(hDir+"BidimME", "(E fraction)**0.7 vs Mass/(Mass+1)" ,50,0.218261406,1.,50,0.,1.);
142 m_histoHelper->bookHisto(hDir+"DRJPVSV", "DeltaR between jet axis and (PV,SV) axis",100,0.,0.5);
143 } else if (m_SVmode == "SV2") {
144 std::string hDir = "/RefFileSV2"+m_jetCollectionList[ijc]+"/"+m_hypotheses[ih]+"/";
145
146 m_histoHelper->bookHisto(hDir+"N2TEffSV2", "Number of Good Two Track Vertices",9,xbi);
147 m_histoHelper->bookHisto(hDir+"N2TNormSV2", "Number of Good Two Track Vertices",30,0.,30.);
148 m_histoHelper->bookHisto(hDir+"TridimMEN2T", "ln(N) vs (E fraction)**0.7 vs Mass/(Mass+1)" ,20,0.,1.,20,0.,1.,7,0.,3.8);
149 if(ih==0) {
150
151 hDir = "/RefFileSV2"+m_jetCollectionList[ijc]+"/controlSV/";
152 m_histoHelper->bookHisto(hDir+"eta","eta",60,-3.,3.);
153 m_histoHelper->bookHisto(hDir+"phi","phi",64,-3.2,3.2);
154 m_histoHelper->bookHisto(hDir+"pt","pt",50,0.,300.);
155 }
156 }
157 }
158 }
159 m_histoHelper->print();
160 }
161
162
163 c_mom = 1000.;
164
165 m_expos = 0.7;
166
167 m_nbjet = 0;
168 m_nljet = 0;
169 return StatusCode::SUCCESS;
170 }
171
172
173 StatusCode SVTag::finalize()
174 {
175 if( m_runModus == "reference" ) {
176 ATH_MSG_INFO("#BTAG# Number of jets used for calibration for reference "
177 << m_refType << " : #b= " << m_nbjet << " #light= " << m_nljet
178 );
179 }
180 return StatusCode::SUCCESS;
181 }
182
183 void SVTag::tagJet(Jet& jetToTag)
184 {
185
186 std::string author = jetToTag.jetAuthor();
187 if (m_doForcedCalib) author = m_ForcedCalibName;
188 ATH_MSG_VERBOSE("#BTAG# Using jet type " << author << " for calibrations.");
189
190
191
192 bool keepInfoPlus = false;
193 if (std::find( m_jetWithInfoPlus.begin(),
194 m_jetWithInfoPlus.end(),
195 author ) != m_jetWithInfoPlus.end()) keepInfoPlus = true;
196
197
198 if(m_runModus == "analysis" && keepInfoPlus) {
199 bool ppb = true;
200 if (m_SaveSVTrackInfo) {
201 StatusCode sc = evtStore()->retrieve(m_originalTPCollection, m_originalTPCollectionName);
202 if (sc.isFailure()) {
203 ATH_MSG_ERROR("#BTAG# " << m_originalTPCollectionName << " not found in StoreGate.");
204 } else {
205 ATH_MSG_VERBOSE("#BTAG# TrackParticleContainer " << m_originalTPCollectionName
206 << " found.");
207 ppb = false;
208 }
209 }
210 if(ppb) {
211 ATH_MSG_VERBOSE("#BTAG# Not able to persistify infos ! Exiting...");
212 return;
213 }
214 }
215
216
217 double jeteta = jetToTag.eta(), jetphi = jetToTag.phi(), jetpt = jetToTag.pt();
218 ATH_MSG_VERBOSE("#BTAG# Jet properties : eta = "<<jeteta
219 <<" phi = "<<jetphi
220 <<" pT = "<<jetpt/1.e3);
221
222
223 if (m_runModus=="reference" && m_SVmode == "SV2") {
224 if (fabs(jeteta) <= 2.5) {
225 m_histoHelper->fillHisto("/RefFileSV2"+author+"/controlSV/eta",(double)jeteta);
226 m_histoHelper->fillHisto("/RefFileSV2"+author+"/controlSV/phi",(double)jetphi);
227 m_histoHelper->fillHisto("/RefFileSV2"+author+"/controlSV/pt",(double)jetpt/1.e3);
228 }
229 }
230
231
232
233 double ambtot = -1., xratio = -1., distnrm = 0., drJPVSV = 0.;
234 long int NSVPair = -1, npsec = -1, npprm = -1;
235
236 Trk::RecVertex myVert;
237 std::vector<const Trk::TrackParticleBase*> TrkList;
238
239
240
241
242 bool vertexIsThere(true);
243
244 const ISvxAssociation* newSvxAssociation=jetToTag.getAssociation<ISvxAssociation>(m_secVxFinderName);
245
246 if (newSvxAssociation==0) {
247 ATH_MSG_DEBUG( "#BTAG# No VKalVrt vertex found, attached to the Jet as association. Not going on with tagging...");
248 return;
249 } else {
250 ATH_MSG_DEBUG( "#BTAG# VKalVrt vertex info " << m_secVxFinderName << " found.");
251 }
252
253 const Trk::VxSecVertexInfo* myVertexInfo=newSvxAssociation->vertexInfo();
254 if (myVertexInfo==0) {
255 vertexIsThere=false;
256
257
258 }
259
260 if (vertexIsThere)
261 {
262
263 const Trk::VxSecVKalVertexInfo* myVKalVertexInfo=dynamic_cast<const Trk::VxSecVKalVertexInfo*>(myVertexInfo);
264
265 if (myVKalVertexInfo==0) {
266 ATH_MSG_WARNING( "#BTAG# Could not cast to a VKalVrt Info object ");
267 return;
268 }
269
270 const std::vector<Trk::VxCandidate*> & myVertices=myVKalVertexInfo->vertices();
271
272 std::vector<const Trk::RecVertex*> secondaryRecVertices;
273
274 if (myVertices.size()>0) {
275
276 ambtot = myVKalVertexInfo->mass()/c_mom;
277 xratio = myVKalVertexInfo->energyFraction();
278 NSVPair = (long int)(myVKalVertexInfo->n2trackvertices());
279
280
281
282
283
284 if (myVertices[0]!=0) {
285 Hep3Vector jetDir = jetToTag.hlv().vect();
286 Hep3Vector PvSvDir( myVertices[0]->recVertex().position().x() - m_priVtx->recVertex().position().x(),
287 myVertices[0]->recVertex().position().y() - m_priVtx->recVertex().position().y(),
288 myVertices[0]->recVertex().position().z() - m_priVtx->recVertex().position().z() );
289 drJPVSV = jetDir.deltaR(PvSvDir);
290 }else{
291 ATH_MSG_VERBOSE("#BTAG# No secondary vertex.");
292 }
293
294 const std::vector<Trk::VxCandidate*>::const_iterator verticesBegin=myVertices.begin();
295 const std::vector<Trk::VxCandidate*>::const_iterator verticesEnd=myVertices.end();
296
297 npsec=0;
298
299 for (std::vector<Trk::VxCandidate*>::const_iterator verticesIter=verticesBegin;
300 verticesIter!=verticesEnd;verticesIter++) {
301 if ((*verticesIter)==0)
302 {
303 ATH_MSG_ERROR("#BTAG# Secondary vertex from InDetVKalVxInJetFinder has zero pointer. Skipping... ");
304 continue;
305 }
306 const std::vector<Trk::VxTrackAtVertex*>* myTracks=(*verticesIter)->vxTrackAtVertex();
307 secondaryRecVertices.push_back(&((*verticesIter)->recVertex()));
308 if (myTracks!=0) {
309 npsec+=myTracks->size();
310 const std::vector<Trk::VxTrackAtVertex*>::const_iterator tracksBegin=myTracks->begin();
311 const std::vector<Trk::VxTrackAtVertex*>::const_iterator tracksEnd=myTracks->end();
312 for (std::vector<Trk::VxTrackAtVertex*>::const_iterator tracksIter=tracksBegin;
313 tracksIter!=tracksEnd;++tracksIter) {
314 const Trk::ITrackLink* myTrackLink=((*tracksIter)->trackOrParticleLink());
315
316 if (myTrackLink==0) {
317 ATH_MSG_WARNING("#BTAG# No ITrack Link attached to the track at the sec vertex... ");
318 }
319 const Trk::LinkToTrackParticleBase* myTPLink=
320 dynamic_cast<const Trk::LinkToTrackParticleBase*>(myTrackLink);
321 if (myTPLink!=0) {
322 TrkList.push_back(**myTPLink);
323 } else {
324 ATH_MSG_WARNING("#BTAG# Cannot cast to TrackParticlBase. Running on tracks not yet foreseen... ");
325 }
326 }
327 }
328 }
329
330 if (myVertices[0]!=0) {
331 const Trk::RecVertex & myVertex=myVertices[0]->recVertex();
332 myVert = myVertex;
333
334 if (m_priVtx)
335 {
336 distnrm=get3DSignificance(m_priVtx->recVertex(),
337 secondaryRecVertices,
338 jetToTag.hlv().vect());
339 }
340 else
341 {
342 ATH_MSG_WARNING("#BTAG# Tagging requested, but no primary vertex supplied.");
343 distnrm=0.;
344 }
345 }
346 else
347 {
348 ATH_MSG_VERBOSE("#BTAG# No vertex. Cannot calculate normalized distance.");
349 distnrm=0.;
350 }
351
352 ATH_MSG_VERBOSE("#BTAG# SVX x = " << myVert.position().x() << " y = " << myVert.position().y() << " z = " << myVert.position().z());
353 ATH_MSG_VERBOSE("#BTAG# The SVX prop. sign3d: " << distnrm);
354 ATH_MSG_VERBOSE("#BTAG# N tracks in SVX "<<npsec);
355 ATH_MSG_VERBOSE("#BTAG# Svx Mass = "<< ambtot << " Svx E frac = " << xratio << " NGood2TrackVertices = " << NSVPair);
356 } else {
357 keepInfoPlus = false;
358 ATH_MSG_VERBOSE("#BTAG# No SVx !");
359 }
360 }
361 else
362 {
363 keepInfoPlus = false;
364 ATH_MSG_VERBOSE("#BTAG# No SVx !");
365 }
366
367
368
369 SVInfoBase *mySVInfoBase(0);
370 SVInfoPlus *mySVInfoPlus(0);
371 bool isIPInfoPlusAlreadyThere = false;
372 if (m_runModus=="analysis") {
373 if (m_writeInfoBase) {
374 mySVInfoBase = new SVInfoBase(m_SVmode);
375 }
376
377 if(keepInfoPlus) {
378 ATH_MSG_VERBOSE("#BTAG# Filling SVInfoPlus...");
379 mySVInfoPlus = const_cast<SVInfoPlus*>(jetToTag.tagInfo<SVInfoPlus>(m_infoPlusName));
380 if(mySVInfoPlus) {
381 ATH_MSG_VERBOSE("#BTAG# Previous SVInfoPlus " << m_infoPlusName << " found...");
382 isIPInfoPlusAlreadyThere = true;
383 } else {
384 mySVInfoPlus = new SVInfoPlus(m_infoPlusName);
385 if(mySVInfoPlus) {
386 ATH_MSG_VERBOSE("#BTAG# " << m_infoPlusName << " not found. New SVInfoPlus created...");
387 mySVInfoPlus->setRecSvx(myVert);
388 mySVInfoPlus->setNGTrackInJet(npprm);
389 mySVInfoPlus->setNGTrackInSvx(npsec);
390 mySVInfoPlus->setN2T(NSVPair);
391 mySVInfoPlus->setMass(ambtot);
392 mySVInfoPlus->setEnergyFraction(xratio);
393 mySVInfoPlus->setNormDist(distnrm);
394 for (uint itk = 0; itk < TrkList.size(); itk++) {
395 const Rec::TrackParticle* myTP=dynamic_cast<const Rec::TrackParticle*>(TrkList[itk]);
396 if (myTP==0) {
397 ATH_MSG_WARNING("#BTAG# Could not cast back to TrackParticle ");
398 } else {
399 SVTrackInfo tinfo(m_originalTPCollection, myTP);
400 mySVInfoPlus->addTrackInfo(tinfo);
401 }
402 }
403 }
404 }
405 }
406 }
407
408
409 ATH_MSG_VERBOSE("#BTAG# SV mode = " << m_SVmode);
410
411 if(m_SVmode != "SV0" ) {
412 float ambtotp = ambtot != -1 ? ambtot/(1.+ambtot): 0.;
413 float xratiop = xratio > 0. ? (float)pow(xratio,m_expos) : 0.;
414 std::string pref = "";
415 if (m_runModus=="reference") {
416 if (jetpt >= m_pTjetmin && fabs(jeteta) <= 2.5) {
417 std::string label = "N/A";
418 double deltaRmin(0.);
419 const TruthInfo* mcTrueInfo = jetToTag.tagInfo<TruthInfo>("TruthInfo");
420 if (mcTrueInfo) {
421 label = mcTrueInfo->jetTruthLabel();
422 ATH_MSG_VERBOSE("#BTAG# label found : ");
423
424 double deltaRtoClosestB = mcTrueInfo->deltaRMinTo("B");
425 double deltaRtoClosestC = mcTrueInfo->deltaRMinTo("C");
426 deltaRmin = deltaRtoClosestB < deltaRtoClosestC ? deltaRtoClosestB : deltaRtoClosestC;
427 } else {
428 ATH_MSG_ERROR("#BTAG# No Label ! Cannot run on reference mode !");
429 return;
430 }
431 if ( ( "B"==m_refType && "B"==label ) ||
432 ( "UDSG"==m_refType && "N/A"==label ) ||
433 ( "ALL"==m_refType &&
434 ( "B"==label || ( "N/A"==label && deltaRmin > m_purificationDeltaR ) ) )
435 ) {
436 if ("B"==label) {
437 pref = m_hypotheses[0];
438 m_nbjet++;
439 } else if ("N/A"==label) {
440 pref = m_hypotheses[1];
441 m_nljet++;
442 }
443 }
444 if (pref == "B" || pref == "U") {
445 std::string hDir = "/RefFile"+m_SVmode+author+"/"+pref+"/";
446
447
448 if (m_SVmode == "SV1") m_histoHelper->fillHisto(hDir+"N2TNormSV1",(float)NSVPair);
449 if (m_SVmode == "SV2") m_histoHelper->fillHisto(hDir+"N2TNormSV2",(float)NSVPair);
450 if (NSVPair > 0 && ambtot > 0.) {
451 if (xratiop == 1.) xratiop = 0.999999;
452 if (m_SVmode == "SV1") {
453 m_histoHelper->fillHisto(hDir+"N2T",(float)NSVPair);
454 m_histoHelper->fillHisto(hDir+"N2TEffSV1",(float)NSVPair);
455 m_histoHelper->fillHisto(hDir+"BidimME",ambtotp,xratiop);
456 m_histoHelper->fillHisto(hDir+"DRJPVSV",(float)drJPVSV);
457 }
458 if (m_SVmode == "SV2") {
459 m_histoHelper->fillHisto(hDir+"N2TEffSV2",(float)NSVPair);
460 m_histoHelper->fillHisto(hDir+"TridimMEN2T",ambtotp,xratiop,log((float)NSVPair));
461 }
462 }
463 }
464 }
465 } else if (m_runModus=="analysis") {
466 std::vector<double> probi;
467
468 double effb = m_likelihoodTool->getEff(m_hypotheses[0], author+"#N2T", m_SVmode);
469 double effu = m_likelihoodTool->getEff(m_hypotheses[1], author+"#N2T", m_SVmode);
470 ATH_MSG_DEBUG( "#BTAG# EFF b,u= " << effb << " " << effu);
471 if (NSVPair>0 && ambtot > 0.) {
472 std::vector<Slice> nslices;
473 AtomicProperty atom2(ambtotp,"SecVtx Transformed Mass");
474 AtomicProperty atom3(xratiop,"SecVtx Transformed Energy Fraction");
475 if (m_SVmode == "SV1") {
476 AtomicProperty atom1(NSVPair,"Number of Two Track Vertices");
477 std::string compoName(author+"#");
478 Composite compo1(compoName+"N2T");
479 Composite compo2(compoName+"BidimME");
480 compo1.atoms.push_back(atom1);
481 compo2.atoms.push_back(atom2);
482 compo2.atoms.push_back(atom3);
483 Slice slice1("SV1");
484 slice1.composites.push_back(compo1);
485 slice1.composites.push_back(compo2);
486 if(m_useDRJPVSV) {
487 AtomicProperty atom4(drJPVSV,"DeltaR between jet axis and (PV,SV) axis");
488 Composite compo3(compoName+"DRJPVSV");
489 compo3.atoms.push_back(atom4);
490 slice1.composites.push_back(compo3);
491 }
492 nslices.push_back(slice1);
493 } else if (m_SVmode == "SV2") {
494 AtomicProperty atom1(log((float)NSVPair),"log(Number of Two Track Vertices)");
495 std::string compoName(author+"#");
496 Composite compo(compoName+"TridimMEN2T");
497 compo.atoms.push_back(atom2);
498 compo.atoms.push_back(atom3);
499 compo.atoms.push_back(atom1);
500 Slice slice1("SV2");
501 slice1.composites.push_back(compo);
502 nslices.push_back(slice1);
503 }
504 m_likelihoodTool->setLhVariableValue(nslices);
505 probi = m_likelihoodTool->calculateLikelihood();
506 ATH_MSG_DEBUG( "#BTAG# WEIGHT: pb, pu= "
507 << probi[0] << " " << probi[1]);
508 if (probi.size() >= 2) {
509 probi[0] *= effb;
510 probi[1] *= effu;
511 } else {
512 ATH_MSG_ERROR("#BTAG# Missing number in jet probabilities ! "<<probi.size());
513 }
514 } else {
515
516
517 probi.push_back((1.-effb));
518 probi.push_back((1.-effu));
519 }
520 if (mySVInfoBase) mySVInfoBase->setTagLikelihood(probi);
521 m_likelihoodTool->clear();
522 }
523
524
525 } else {
526 ATH_MSG_VERBOSE("#BTAG# SV0 Lxy3D significance = " << distnrm);
527 std::vector<double> probi;
528 probi.push_back(distnrm);
529 if (mySVInfoBase) mySVInfoBase->setTagLikelihood(probi);
530 }
531
532
533
534 if(m_runModus=="analysis") {
535 if (mySVInfoBase) {
536 jetToTag.addInfo(mySVInfoBase);
537 mySVInfoBase->makeValid();
538 }
539 if (mySVInfoPlus && !isIPInfoPlusAlreadyThere) jetToTag.addInfo(mySVInfoPlus);
540 }
541 ATH_MSG_VERBOSE("#BTAG# Finalyzing... ");
542
543 return;
544 }
545
546 void SVTag::finalizeHistos() {
547 }
548
549 void SVTag::m_printParameterSettings() {
550 ATH_MSG_INFO("#BTAG# " << name() << "Parameter settings ");
551 ATH_MSG_INFO("#BTAG# I am in " << m_runModus << " modus.");
552 ATH_MSG_INFO("#BTAG# The method is " << m_SVmode);
553 if (m_runModus == "reference") ATH_MSG_INFO("#BTAG# Preparing "<< m_refType<< "-jet probability density functions...");
554 }
555
556 double SVTag::get3DSignificance(const Trk::RecVertex & priVertex,
557 std::vector<const Trk::RecVertex* > & secVertex,
558 const Hep3Vector jetDirection) {
559
560 std::vector<HepVector> positions;
561 std::vector<const HepSymMatrix*> weightMatrices;
562
563 std::vector<const Trk::RecVertex* >::const_iterator secBegin=secVertex.begin();
564 std::vector<const Trk::RecVertex* >::const_iterator secEnd=secVertex.end();
565
566 for (std::vector<const Trk::RecVertex* >::const_iterator secIter=secBegin;
567 secIter!=secEnd;++secIter)
568 {
569 HepVector position(3);
570 position[0]=(*secIter)->position().x();
571 position[1]=(*secIter)->position().y();
572 position[2]=(*secIter)->position().z();
573 positions.push_back(position);
574
575 weightMatrices.push_back(&((*secIter)->errorPosition().weight()));
576
577 }
578
579 HepVector weightTimesPosition(3,0);
580 HepSymMatrix sumWeights(3,0);
581
582 int count=0;
583 for (std::vector<const Trk::RecVertex* >::const_iterator secIter=secBegin;
584 secIter!=secEnd;++secIter)
585 {
586
587 weightTimesPosition+=(*weightMatrices[count])*positions[count];
588 sumWeights+=(*weightMatrices[count]);
589
590 count+=1;
591 }
592
593 int failed(0);
594 HepSymMatrix meanCovariance=sumWeights.inverse(failed);
595 if (failed!=0) {
596 ATH_MSG_ERROR("#BTAG# Could not invert sum of sec vtx matrices");
597 return 0.;
598 }
599
600 HepVector meanPosition=meanCovariance*weightTimesPosition;
601
602 HepVector Diff(3);
603 Diff[0]=meanPosition[0]-priVertex.position().x();
604 Diff[1]=meanPosition[1]-priVertex.position().y();
605 Diff[2]=meanPosition[2]-priVertex.position().z();
606
607 HepSymMatrix CovLaterWeightSum=meanCovariance
608 +priVertex.errorPosition().covariance();
609
610
611 failed=0;
612 CovLaterWeightSum.invert(failed);
613 if (failed!=0) {
614 ATH_MSG_ERROR("#BTAG# Could not invert prima+seco vertex matrix");
615 return 0;
616 }
617
618 Hep3Vector difference(Diff[0],Diff[1],Diff[2]);
619
620 return std::sqrt(CovLaterWeightSum.similarity(Diff))*difference.unit().dot(jetDirection.unit());
621 }
622 }
623
| 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.
|
|