001 #include "JetTagTools/IPTag.h"
002 #include "CLHEP/Vector/ThreeVector.h"
003 #include "CLHEP/Vector/LorentzVector.h"
004 #include "GaudiKernel/ITHistSvc.h"
005 #include "GaudiKernel/IToolSvc.h"
006 #include "JetEvent/Jet.h"
007 #include "JetTagEvent/TrackAssociation.h"
008 #include "JetTagInfo/IPInfoBase.h"
009 #include "JetTagInfo/IPInfoPlus.h"
010 #include "JetTagInfo/TruthInfo.h"
011 #include "JetTagInfo/SvxSummary.h"
012 #include "JetTagTools/NewLikelihoodTool.h"
013 #include "JetTagTools/LikelihoodComponents.h"
014 #include "JetTagTools/HistoHelperRoot.h"
015 #include "JetTagTools/TrackSelector.h"
016 #include "JetTagTools/GradedTrack.h"
017 #include "JetTagTools/SVForIPTool.h"
018 #include "JetTagInfo/TrackGrade.h"
019 #include "JetTagInfo/TrackGradesDefinition.h"
020 #include "JetTagTools/ITrackGradeFactory.h"
021 #include "Particle/TrackParticle.h"
022 #include "Navigation/NavigationToken.h"
023 #include "ITrackToVertex/ITrackToVertex.h"
024 #include "TrkParticleBase/TrackParticleBase.h"
025 #include "TrkVertexFitterInterfaces/ITrackToVertexIPEstimator.h"
026 #include "TH1.h"
027 #include <cmath>
028 #include <sstream>
029 #include <algorithm>
030 #include <vector>
031
032 namespace Analysis {
033
034
035
036
037
038
039
040 typedef std::vector<double> FloatVec;
041 typedef std::vector<double>::iterator FloatVecIter;
042
043 IPTag::IPTag(const std::string& t, const std::string& n, const IInterface* p)
044 : AthAlgTool(t,n,p),
045 m_runModus("analysis"),
046 m_histoHelper(0),
047 m_trackToVertexTool("Reco::TrackToVertex"),
048 m_trackSelectorTool("Analysis::TrackSelector"),
049 m_likelihoodTool("Analysis::NewLikelihoodTool"),
050 m_secVxFinderNameForV0Removal("InDetVKalVxInJetTool"),
051 m_secVxFinderNameForIPSign("InDetVKalVxInJetTool"),
052 m_SVForIPTool("Analysis::SVForIPTool"),
053 m_trackGradeFactory("Analysis::BasicTrackGradeFactory"),
054 m_unbiasIPEstimation(false) {
055
056 declareInterface<ITagTool>(this);
057
058 declareProperty("Runmodus", m_runModus);
059 declareProperty("hypotheses", m_hypotheses);
060 m_hypotheses.push_back("B");
061 m_hypotheses.push_back("U");
062 declareProperty("useVariables", m_useVariables);
063 declareProperty("impactParameterView", m_impactParameterView = "2D");
064 declareProperty("SignWithSvx", m_SignWithSvx = false);
065 declareProperty("SVForIPTool",m_SVForIPTool);
066
067
068 declareProperty("trackGradePartitions", m_trackGradePartitionsDefinition);
069 m_trackGradePartitionsDefinition.push_back("Good");
070 declareProperty("RejectBadTracks", m_RejectBadTracks = false);
071
072
073 declareProperty("LikelihoodTool", m_likelihoodTool);
074 declareProperty("trackSelectorTool", m_trackSelectorTool);
075 declareProperty("trackToVertexTool", m_trackToVertexTool);
076 declareProperty("trackGradeFactory",m_trackGradeFactory);
077
078
079 declareProperty("originalTPCollectionName", m_originalTPCollectionName = "TrackParticleCandidate");
080 declareProperty("writeInfoBase", m_writeInfoBase = true);
081 declareProperty("infoPlusName", m_infoPlusName = "IPInfoPlus");
082
083
084 declareProperty("referenceType", m_referenceType = "ALL");
085 declareProperty("truthMatchingName", m_truthMatchingName = "TruthInfo");
086 declareProperty("checkOverflows", m_checkOverflows = false);
087 declareProperty("purificationDeltaR", m_purificationDeltaR = 0.8);
088 declareProperty("jetPtMinRef", m_jetPtMinRef = 15.*GeV);
089
090 declareProperty("jetCollectionList", m_jetCollectionList);
091 declareProperty("jetWithInfoPlus", m_jetWithInfoPlus);
092 m_jetWithInfoPlus.push_back("Cone4H1Tower");
093 declareProperty("useForcedCalibration", m_doForcedCalib = false);
094 declareProperty("ForcedCalibrationName", m_ForcedCalibName = "Cone4H1Tower");
095
096 declareProperty("SecVxFinderNameForV0Removal",m_secVxFinderNameForV0Removal);
097 declareProperty("SecVxFinderNameForIPSign",m_secVxFinderNameForIPSign);
098
099 declareProperty("TrackToVertexIPEstimator",m_trackToVertexIPEstimator);
100
101 declareProperty("unbiasIPEstimation",m_unbiasIPEstimation);
102
103 }
104
105 IPTag::~IPTag() {
106 delete m_histoHelper;
107 for (size_t i = 0; i < m_trackGradePartitions.size(); i++)
108 delete m_trackGradePartitions[i];
109 }
110
111
112 StatusCode IPTag::initialize() {
113
114
115 if ( m_trackToVertexTool.retrieve().isFailure() ) {
116 ATH_MSG_FATAL("#BTAG# Failed to retrieve tool " << m_trackToVertexTool);
117 return StatusCode::FAILURE;
118 } else {
119 ATH_MSG_DEBUG("#BTAG# Retrieved tool " << m_trackToVertexTool);
120 }
121
122 if (m_SVForIPTool.retrieve().isFailure() ) {
123 ATH_MSG_FATAL("#BTAG# Failed to retrieve tool " << m_SVForIPTool);
124 return StatusCode::FAILURE;
125 } else {
126 ATH_MSG_DEBUG("#BTAG# Retrieved tool " << m_SVForIPTool);
127 }
128
129 if (m_trackToVertexIPEstimator.retrieve().isFailure() ) {
130 ATH_MSG_FATAL("#BTAG# Failed to retrieve tool " << m_trackToVertexIPEstimator);
131 return StatusCode::FAILURE;
132 } else {
133 ATH_MSG_DEBUG("#BTAG# Retrieved tool " << m_trackToVertexIPEstimator);
134 }
135
136
137
138 if ( m_trackSelectorTool.retrieve().isFailure() ) {
139 ATH_MSG_FATAL("#BTAG# Failed to retrieve tool " << m_trackSelectorTool);
140 return StatusCode::FAILURE;
141 } else {
142 ATH_MSG_DEBUG("#BTAG# Retrieved tool " << m_trackSelectorTool);
143 }
144
145
146 ITHistSvc* myHistoSvc;
147 if( service( "THistSvc", myHistoSvc ).isSuccess() ) {
148 ATH_MSG_DEBUG("#BTAG# HistoSvc loaded successfully.");
149 m_histoHelper = new HistoHelperRoot(myHistoSvc);
150 m_histoHelper->setCheckOverflows(m_checkOverflows);
151 } else {
152 ATH_MSG_ERROR("#BTAG# HistoSvc could NOT bo loaded.");
153 }
154
155
156 if ( m_trackGradeFactory.retrieve().isFailure() ) {
157 ATH_MSG_FATAL("#BTAG# Failed to retrieve tool " << m_trackGradeFactory);
158 return StatusCode::FAILURE;
159 } else {
160 ATH_MSG_DEBUG("#BTAG# Retrieved tool " << m_trackGradeFactory);
161 }
162
163
164
165 int nbPart = m_trackGradePartitionsDefinition.size();
166 ATH_MSG_INFO("#BTAG# Defining " << nbPart <<" track partitions: ");
167 for(int i=0;i<nbPart;i++) {
168 TrackGradePartition* part(0);
169 try {
170 part = new TrackGradePartition(m_trackGradePartitionsDefinition[i],
171 *m_trackGradeFactory);
172 }
173 catch (std::string error) {
174 ATH_MSG_ERROR("#BTAG# Reported error " << error);
175 ATH_MSG_ERROR("#BTAG# List of categories provided to IPTag by jO : ");
176 for (int l=0;l<nbPart;l++) {
177 ATH_MSG_ERROR("#BTAG# string " << m_trackGradePartitionsDefinition[l]);
178 }
179 ATH_MSG_ERROR("#BTAG# List of categories provided by the TrackGradeFactory " << m_trackGradeFactory);
180
181 const TrackGradesDefinition & trackFactoryGradesDefinition = m_trackGradeFactory->getTrackGradesDefinition();
182 const std::vector<TrackGrade> & gradeList = trackFactoryGradesDefinition.getList();
183 std::vector<TrackGrade>::const_iterator listIter=gradeList.begin();
184 std::vector<TrackGrade>::const_iterator listEnd=gradeList.end();
185 if (msgLvl(MSG::ERROR)) {
186 for ( ; listIter !=listEnd ; ++listIter ) {
187 ATH_MSG_ERROR("#BTAG# n. " << (*listIter).gradeNumber() << " string " << (*listIter).gradeString());
188 }
189 }
190 ATH_MSG_ERROR("#BTAG# Terminating now... ");
191 return StatusCode::FAILURE;
192 }
193 ATH_MSG_INFO((*part));
194 m_trackGradePartitions.push_back(part);
195 }
196
197
198 if( m_runModus == "analysis" ) {
199 if ( m_likelihoodTool.retrieve().isFailure() ) {
200 ATH_MSG_FATAL("#BTAG# Failed to retrieve tool " << m_likelihoodTool);
201 return StatusCode::FAILURE;
202 } else {
203 ATH_MSG_DEBUG("#BTAG# Retrieved tool " << m_likelihoodTool);
204 }
205 m_likelihoodTool->defineHypotheses(m_hypotheses);
206
207 for(unsigned int i=0;i<m_trackGradePartitions.size();i++) {
208 for(unsigned int ih=0;ih<m_hypotheses.size();ih++) {
209 std::string hName;
210 if(m_impactParameterView=="1D")
211 hName = m_hypotheses[ih]+"/"+m_trackGradePartitions[i]->suffix()+"/SipZ0";
212 if(m_impactParameterView=="2D")
213 hName = m_hypotheses[ih]+"/"+m_trackGradePartitions[i]->suffix()+"/SipA0";
214 if(m_impactParameterView=="3D")
215 hName = m_hypotheses[ih]+"/"+m_trackGradePartitions[i]->suffix()+"/Sip3D";
216 m_likelihoodTool->defineHistogram(hName);
217 }
218 }
219 m_likelihoodTool->printStatus();
220 }
221
222 const TrackGradesDefinition & trackFactoryGradesDefinition = m_trackGradeFactory->getTrackGradesDefinition();
223 const std::vector<TrackGrade> & gradeList = trackFactoryGradesDefinition.getList();
224 std::vector<TrackGrade>::const_iterator listBegin=gradeList.begin();
225 std::vector<TrackGrade>::const_iterator listEnd=gradeList.end();
226
227
228 if( m_runModus == "reference" ) {
229 for(uint j=0;j<m_jetCollectionList.size();j++) {
230
231
232
233 for (std::vector<TrackGrade>::const_iterator listIter=listBegin ; listIter !=listEnd ; ++listIter ) {
234 const TrackGrade & grd = (*listIter);
235 if(m_impactParameterView=="1D") {
236 for(uint ih=0;ih<m_hypotheses.size();ih++) {
237 std::string hName = "/RefFileIP1D" + m_jetCollectionList[j] + "/"
238 + m_hypotheses[ih] + "/" + (std::string)grd + "/SipZ0";
239 ATH_MSG_VERBOSE("#BTAG# booking for IP1D: " << hName);
240 m_histoHelper->bookHisto(hName,"Signed Impact Parameter (z)",240,-20.,40.);
241 }
242 }
243 if(m_impactParameterView=="2D") {
244 for(uint ih=0;ih<m_hypotheses.size();ih++) {
245 std::string hName = "/RefFileIP2D" + m_jetCollectionList[j] + "/"
246 + m_hypotheses[ih] + "/" + (std::string)grd + "/SipA0";
247 ATH_MSG_VERBOSE("#BTAG# booking for IP2D: " << hName);
248 m_histoHelper->bookHisto(hName,"Signed Impact Parameter (rphi)",240,-20.,40.);
249 }
250 }
251 if(m_impactParameterView=="3D") {
252 const int nbx = 32;
253 const int nby = 14;
254 double xbi[32] = {-20.0,-10.0,-8.0,-6.0,-5.5,-5.0,-4.5,-4.0,-3.5,-3.0,-2.5,-2.0,-1.5,-1.0,-0.5,
255 0.0,0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0,8.0,10.0,20.0,40.0};
256 double ybi[14] = {-20.,-8.,-4.,-3.,-2.,-1.,0.,1.,2.,3.,4.,8.,12.,40.};
257 for(uint ih=0;ih<m_hypotheses.size();ih++) {
258 std::string hName = "/RefFileIP3D" + m_jetCollectionList[j] + "/"
259 + m_hypotheses[ih] + "/" + (std::string)grd + "/Sip3D";
260 ATH_MSG_VERBOSE("#BTAG# booking for IP3D: " << hName);
261 m_histoHelper->bookHisto(hName,"Signed IP Z0 vs (rphi)",nbx-1,xbi,nby-1,ybi);
262 }
263 }
264 }
265 }
266 m_histoHelper->print();
267 }
268 m_nbjet = 0;
269 m_nljet = 0;
270 return StatusCode::SUCCESS;
271 }
272
273
274 StatusCode IPTag::finalize() {
275 if( m_runModus == "reference" ) {
276 ATH_MSG_INFO("#BTAG# Number of jets used for calibration for reference "
277 << m_referenceType << " : #b= " << m_nbjet << " #light= " << m_nljet );
278 }
279 return StatusCode::SUCCESS;
280 }
281
282 void IPTag::tagJet(Jet& jetToTag) {
283
284
285 std::string author = jetToTag.jetAuthor();
286 if (m_doForcedCalib) {
287 author = m_ForcedCalibName;
288 }
289 ATH_MSG_VERBOSE("#BTAG# Using jet type " << author << " for calibrations.");
290
291
292 bool keepInfoPlus = false;
293 if(std::find( m_jetWithInfoPlus.begin(),
294 m_jetWithInfoPlus.end(),
295 author ) != m_jetWithInfoPlus.end()) keepInfoPlus = true;
296 if(keepInfoPlus) ATH_MSG_VERBOSE("#BTAG# Writing infoPlus for this jet type.");
297
298
299 if(m_runModus == "analysis" && keepInfoPlus) {
300 StatusCode sc;
301 sc = evtStore()->retrieve(m_originalTPCollection, m_originalTPCollectionName);
302 if (sc.isFailure()) {
303 ATH_MSG_ERROR("#BTAG# TrackParticleContainer " << m_originalTPCollectionName << " not found.");
304 return;
305 } else {
306 ATH_MSG_VERBOSE("#BTAG# TrackParticleContainer " << m_originalTPCollectionName << " found.");
307 }
308 }
309
310
311 std::string label = "N/A";
312 std::string pref = "";
313 if( m_runModus == "reference" ) {
314
315 if( jetToTag.pt()>m_jetPtMinRef && fabs(jetToTag.eta())<2.5 ) {
316
317 const TruthInfo* mcinfo = jetToTag.tagInfo<TruthInfo>("TruthInfo");
318 double deltaRmin(0.);
319 if( mcinfo ) {
320 label = mcinfo->jetTruthLabel();
321
322 double deltaRtoClosestB = mcinfo->deltaRMinTo("B");
323 double deltaRtoClosestC = mcinfo->deltaRMinTo("C");
324 deltaRmin = deltaRtoClosestB < deltaRtoClosestC ? deltaRtoClosestB : deltaRtoClosestC;
325 double deltaRtoClosestT = mcinfo->deltaRMinTo("T");
326 deltaRmin = deltaRtoClosestT < deltaRmin ? deltaRtoClosestT : deltaRmin;
327 } else {
328 ATH_MSG_ERROR("#BTAG# No TruthInfo ! Cannot run on reference mode !");
329 return;
330 }
331 if ( ( "B"==m_referenceType && "B"==label ) ||
332 ( "UDSG"==m_referenceType && "N/A"==label ) ||
333 ( "ALL"==m_referenceType &&
334 ( "B"==label || ( "N/A"==label && deltaRmin > m_purificationDeltaR ) ) )
335 ) {
336 if ("B"==label) {
337 pref = m_hypotheses[0];
338 m_nbjet++;
339 } else if ("N/A"==label) {
340 pref = m_hypotheses[1];
341 m_nljet++;
342 }
343 }
344 }
345 }
346
347 m_tracksInJet.clear();
348 int nbPart = m_trackGradePartitionsDefinition.size();
349
350 std::vector<const Trk::TrackParticleBase*> TrkFromV0;
351 Hep3Vector SvxDirection;
352 bool canUseSvxDirection=false;
353
354 if (m_SignWithSvx) {
355 m_SVForIPTool->getDirectionFromSecondaryVertexInfo(SvxDirection,canUseSvxDirection,
356 jetToTag,m_secVxFinderNameForIPSign,(m_priVtx)->recVertex());
357 }
358
359
360 m_SVForIPTool->getTrkFromV0FromSecondaryVertexInfo(TrkFromV0,
361 jetToTag,m_secVxFinderNameForV0Removal);
362 ATH_MSG_VERBOSE("#BTAG# TrkFromV0 : number of reconstructed bad tracks: " << TrkFromV0.size());
363
364
365 int nbTrak = 0;
366 m_trackSelectorTool->primaryVertex(m_priVtx->recVertex().position());
367 m_trackSelectorTool->prepare();
368 std::vector<const Rec::TrackParticle*>* trackVector = jetToTag.getAssociation<TrackAssociation>("Tracks")->tracks();
369 std::vector<const Rec::TrackParticle*>::iterator jetItr;
370 for( jetItr = trackVector->begin(); jetItr != trackVector->end() ; ++jetItr ) {
371 const Rec::TrackParticle * aTemp = *jetItr;
372 nbTrak++;
373 if( m_trackSelectorTool->selectTrack(aTemp) ) {
374 TrackGrade * theGrade = m_trackGradeFactory->getGrade(*aTemp,
375 jetToTag.hlv() );
376 ATH_MSG_VERBOSE("#BTAG# result of selectTrack is OK, grade= " << (std::string)(*theGrade) );
377
378 bool tobeUsed = false;
379 for(int i=0;i<nbPart;i++) {
380 if (std::find( (m_trackGradePartitions[i]->grades()).begin(),
381 (m_trackGradePartitions[i]->grades()).end(),
382 *theGrade )
383 != (m_trackGradePartitions[i]->grades()).end()) tobeUsed = true;
384 }
385
386 if (std::find(TrkFromV0.begin(),TrkFromV0.end(),aTemp) != TrkFromV0.end()) {
387 ATH_MSG_VERBOSE("#BTAG# Bad track in jet, pt = " << aTemp->pt() << " eta = "
388 << aTemp->eta() << " phi = " << aTemp->phi() );
389 if (m_RejectBadTracks) tobeUsed = false;
390 }
391 if (tobeUsed) m_tracksInJet.push_back(GradedTrack(aTemp, *theGrade));
392 delete theGrade;
393 theGrade=0;
394 }
395 }
396 delete trackVector;
397 ATH_MSG_VERBOSE("#BTAG# #tracks = " << nbTrak);
398 ATH_MSG_VERBOSE("#BTAG# the z of the primary = " << m_priVtx->recVertex().position().z());
399
400
401 Hep3Vector jetDirection(jetToTag.px(),jetToTag.py(),jetToTag.pz());
402 Hep3Vector unit = jetDirection.unit();
403 if (m_SignWithSvx && canUseSvxDirection) {
404 unit = SvxDirection.unit();
405 ATH_MSG_DEBUG("#BTAG# Using direction from sec vertex finder: " <<
406 " phi: " << unit.phi() << " theta: " << unit.theta() <<
407 " instead of jet direction phi: " << jetDirection.phi() <<
408 " theta: " << jetDirection.theta() );
409 }
410
411 FloatVec vectD0;
412 FloatVec vectD0Signi;
413 FloatVec vectZ0;
414 FloatVec vectZ0Signi;
415 std::vector<TrackGrade> vectGrades;
416 std::vector<const Rec::TrackParticle*> vectTP;
417 std::vector<bool> vectFromV0;
418
419 const int nbTrackMean = 5;
420 vectD0.reserve(nbTrackMean);
421 vectD0Signi.reserve(nbTrackMean);
422 vectZ0.reserve(nbTrackMean);
423 vectZ0Signi.reserve(nbTrackMean);
424 vectGrades.reserve(nbTrackMean);
425 vectTP.reserve(nbTrackMean);
426 vectFromV0.reserve(nbTrackMean);
427
428 for (std::vector<GradedTrack>::iterator trkItr = m_tracksInJet.begin();
429 trkItr != m_tracksInJet.end(); ++trkItr) {
430 const Rec::TrackParticle* trk = (*trkItr).track;
431 bool isFromV0 = (std::find(TrkFromV0.begin(),TrkFromV0.end(),trk) != TrkFromV0.end());
432
433
434 double d0wrtPriVtx(0.);
435 double z0wrtPriVtx(0.);
436 double d0ErrwrtPriVtx(1.);
437 double z0ErrwrtPriVtx(1.);
438
439
440 const Trk::ImpactParametersAndSigma* myIPandSigma(0);
441 if (m_trackToVertexIPEstimator) {
442 myIPandSigma = m_trackToVertexIPEstimator->estimate(trk,m_priVtx,m_unbiasIPEstimation);
443 }
444 if(0==myIPandSigma) {
445 ATH_MSG_WARNING("#BTAG# IPTAG: trackToVertexIPEstimator failed !");
446 } else {
447 d0wrtPriVtx=myIPandSigma->IPd0;
448 d0ErrwrtPriVtx=myIPandSigma->sigmad0;
449 z0wrtPriVtx=myIPandSigma->IPz0SinTheta;
450 z0ErrwrtPriVtx=myIPandSigma->sigmaz0SinTheta;
451 delete myIPandSigma;
452 myIPandSigma=0;
453 }
454
455
456
457 double signOfIP(1.);
458 if(m_impactParameterView=="2D" || m_impactParameterView=="1D") {
459 signOfIP=m_trackToVertexIPEstimator->get2DLifetimeSignOfTrack(trk->definingParameters(),
460 unit,m_priVtx->recVertex());
461 }
462
463 if(m_impactParameterView=="3D") {
464 signOfIP=m_trackToVertexIPEstimator->get3DLifetimeSignOfTrack(trk->definingParameters(),
465 unit,m_priVtx->recVertex());
466 }
467 double signOfZIP = m_trackToVertexIPEstimator->getZLifetimeSignOfTrack(trk->definingParameters(),
468 unit,m_priVtx->recVertex());
469
470
471 double sIP = signOfIP*fabs(d0wrtPriVtx);
472 double significance= signOfIP*fabs(d0wrtPriVtx/d0ErrwrtPriVtx);
473 double szIP = signOfZIP*fabs(z0wrtPriVtx);
474 double z0Sig = signOfZIP*fabs(z0wrtPriVtx/z0ErrwrtPriVtx);
475
476 vectD0.push_back(sIP);
477 vectD0Signi.push_back(significance);
478 vectZ0.push_back(szIP);
479 vectZ0Signi.push_back(z0Sig);
480 vectGrades.push_back((*trkItr).grade);
481 vectTP.push_back(trk);
482 vectFromV0.push_back(isFromV0);
483
484 msg(MSG::VERBOSE) << "#BTAG# IPTAG: Trk: grade= " << (std::string)(*trkItr).grade
485 << " Eta= " << trk->eta() << " Phi= " << trk->phi() << " pT= " << trk->pt()
486 << " d0= " << sIP
487 << "+-" << d0ErrwrtPriVtx
488 << " sigd0= " << significance
489 << " z0= " << szIP
490 << "+-" << z0ErrwrtPriVtx
491 << " sigz0= " << z0Sig << endreq;
492
493
494
495 if( m_runModus == "reference" ) {
496 if( pref != "" ) {
497 ATH_MSG_DEBUG("#BTAG# filling ref histo for " << pref);
498 const TrackGradesDefinition & trackFactoryGradesDefinition = m_trackGradeFactory->getTrackGradesDefinition();
499 const std::vector<TrackGrade> & gradeList = trackFactoryGradesDefinition.getList();
500 std::vector<TrackGrade>::const_iterator listIter=gradeList.begin();
501 std::vector<TrackGrade>::const_iterator listEnd=gradeList.end();
502 for ( ; listIter !=listEnd ; ++listIter ) {
503 const TrackGrade & grd = (*listIter);
504 ATH_MSG_DEBUG("#BTAG# filling ref histo for grade " << (std::string)grd );
505 if( grd==(*trkItr).grade ) {
506 ATH_MSG_DEBUG("#BTAG# track is of required grade ");
507 const std::string suffix = "_" + (std::string)grd;
508 if( m_impactParameterView=="1D" ) {
509 std::string hName = "/RefFileIP1D" + author + "/"
510 + pref + "/" + (std::string)grd + "/SipZ0";
511 ATH_MSG_DEBUG("#BTAG# histo 1D: " << hName);
512 m_histoHelper->fillHisto(hName,z0Sig);
513 }
514 if( m_impactParameterView=="2D" ) {
515 std::string hName = "/RefFileIP2D" + author + "/"
516 + pref + "/" + (std::string)grd + "/SipA0";
517 ATH_MSG_DEBUG("#BTAG# histo 2D: " << hName);
518 m_histoHelper->fillHisto(hName,significance);
519 }
520 if( m_impactParameterView=="3D" ) {
521 std::string hName = "/RefFileIP3D" + author + "/"
522 + pref + "/" + (std::string)grd + "/Sip3D";
523 ATH_MSG_DEBUG("#BTAG# histo 3D: " << hName);
524 m_histoHelper->fillHisto(hName,significance,z0Sig);
525 }
526 }
527 }
528 }
529 }
530
531 }
532
533
534 if(m_runModus=="analysis") {
535 std::string instanceName("IP"+m_impactParameterView);
536 IPInfoBase* infoBase(0);
537 IPInfoPlus* infoPlus(0);
538
539
540 if(m_writeInfoBase) {
541 infoBase = new IPInfoBase(instanceName);
542 infoBase->nbTracks(vectD0.size());
543 }
544
545 bool isIPInfoPlusAlreadyThere = false;
546 uint nbt = vectTP.size();
547 if(keepInfoPlus && nbt > 0) {
548 ATH_MSG_VERBOSE("#BTAG# Filling IPInfoPlus...");
549 infoPlus = const_cast<IPInfoPlus*>(jetToTag.tagInfo<IPInfoPlus>(m_infoPlusName));
550 if(infoPlus) {
551 ATH_MSG_VERBOSE("#BTAG# Previous IPInfoPlus " << m_infoPlusName << " found...");
552 isIPInfoPlusAlreadyThere = true;
553
554 for(uint i=0;i<nbt;i++) {
555 infoPlus->updateTrackWeight(vectTP[i], m_impactParameterView,
556 this->trackWeight(author,vectGrades[i],vectD0Signi[i],vectZ0Signi[i]));
557 }
558 } else {
559 infoPlus = new IPInfoPlus(m_infoPlusName);
560 if(infoPlus) {
561 ATH_MSG_VERBOSE("#BTAG#" << m_infoPlusName << " not found. New IPInfoPlus created...");
562 for(uint i=0;i<nbt;i++) {
563 IPTrackInfo tinfo(m_originalTPCollection, vectTP[i],
564 vectGrades[i], vectFromV0[i],
565 vectD0[i], vectD0Signi[i],
566 vectZ0[i], vectZ0Signi[i]);
567 tinfo.resetW();
568 infoPlus->addTrackInfo(tinfo);
569 infoPlus->updateTrackWeight(vectTP[i], m_impactParameterView,
570 this->trackWeight(author,vectGrades[i],vectD0Signi[i],vectZ0Signi[i]));
571 }
572 }
573 }
574 }
575
576
577 std::vector<Slice> slices;
578 for(unsigned int i=0; i<vectD0Signi.size(); i++) {
579 if(m_impactParameterView=="3D") {
580 AtomicProperty atom1(vectD0Signi[i],"Significance of IP (rphi)");
581 AtomicProperty atom2(vectZ0Signi[i],"Significance of Z0");
582 std::string compoName(author+"#");
583 compoName += (std::string)vectGrades[i];
584 compoName += "/Sip3D";
585 Composite compo1(compoName);
586 compo1.atoms.push_back(atom1);
587 compo1.atoms.push_back(atom2);
588 Slice slice1("IP3D");
589 slice1.composites.push_back(compo1);
590 slices.push_back(slice1);
591 }
592 if(m_impactParameterView=="2D") {
593 AtomicProperty atom1(vectD0Signi[i],"Significance of IP (rphi)");
594 std::string compoName(author+"#");
595 compoName += (std::string)vectGrades[i];
596 compoName += "/SipA0";
597 Composite compo1(compoName);
598 compo1.atoms.push_back(atom1);
599 Slice slice1("IP2D");
600 slice1.composites.push_back(compo1);
601 slices.push_back(slice1);
602 }
603 if(m_impactParameterView=="1D") {
604 AtomicProperty atom1(vectZ0Signi[i],"Significance of IP (z)");
605 std::string compoName(author+"#");
606 compoName += (std::string)vectGrades[i];
607 compoName += "/SipZ0";
608 Composite compo1(compoName);
609 compo1.atoms.push_back(atom1);
610 Slice slice1("IP1D");
611 slice1.composites.push_back(compo1);
612 slices.push_back(slice1);
613 }
614 }
615 m_likelihoodTool->setLhVariableValue(slices);
616 if(vectD0Signi.size()>0) {
617 std::vector<double> tmp = m_likelihoodTool->calculateLikelihood();
618 if(infoBase) infoBase->setTagLikelihood(tmp);
619 } else {
620 std::vector<double> null;
621 null.push_back(1.);
622 null.push_back(1.e9);
623 if(infoBase) infoBase->setTagLikelihood(null);
624 }
625 std::vector<double> newvecp = m_likelihoodTool->tagLikelihood();
626 if (newvecp.size() == 2) {
627 double pb = newvecp[0], pu = newvecp[1];
628 if (pu != 0) {
629 ATH_MSG_DEBUG("#BTAG# (new) WEIGHT " << pb << " " << pu << " " << log(pb/pu) );
630 } else {
631 ATH_MSG_VERBOSE("#BTAG# light jet proba = 0 !");
632 }
633 } else {
634 ATH_MSG_VERBOSE("#BTAG# No likelihood !");
635 }
636
637
638 if(infoBase) {
639 jetToTag.addInfo(infoBase);
640 if (vectD0Signi.size()>0) infoBase->makeValid();
641 }
642 if(infoPlus && !isIPInfoPlusAlreadyThere) jetToTag.addInfo(infoPlus);
643 m_likelihoodTool->clear();
644 }
645 m_tracksInJet.clear();
646 }
647
648 double IPTag::trackWeight(std::string author, TrackGrade grade, double sa0, double sz0) {
649
650 std::vector<Slice> slices;
651 if(m_impactParameterView=="3D") {
652 AtomicProperty atom1(sa0,"Significance of IP (rphi)");
653 AtomicProperty atom2(sz0,"Significance of Z0");
654 std::string compoName(author+"#");
655 compoName += (std::string)grade;
656 compoName += "/Sip3D";
657 Composite compo1(compoName);
658 compo1.atoms.push_back(atom1);
659 compo1.atoms.push_back(atom2);
660 Slice slice1("IP3D");
661 slice1.composites.push_back(compo1);
662 slices.push_back(slice1);
663 }
664 if(m_impactParameterView=="2D") {
665 AtomicProperty atom1(sa0,"Significance of IP (rphi)");
666 std::string compoName(author+"#");
667 compoName += (std::string)grade;
668 compoName += "/SipA0";
669 Composite compo1(compoName);
670 compo1.atoms.push_back(atom1);
671 Slice slice1("IP2D");
672 slice1.composites.push_back(compo1);
673 slices.push_back(slice1);
674 }
675 if(m_impactParameterView=="1D") {
676 AtomicProperty atom1(sz0,"Significance of IP (z)");
677 std::string compoName(author+"#");
678 compoName += (std::string)grade;
679 compoName += "/SipZ0";
680 Composite compo1(compoName);
681 compo1.atoms.push_back(atom1);
682 Slice slice1("IP1D");
683 slice1.composites.push_back(compo1);
684 slices.push_back(slice1);
685 }
686 m_likelihoodTool->setLhVariableValue(slices);
687 std::vector<double> tmp = m_likelihoodTool->calculateLikelihood();
688 m_likelihoodTool->clear();
689 double pb = tmp[0];
690 double pu = tmp[1];
691 double w = -100.;
692 if(pb<=0.) {
693 w = -30.;
694 } else if (pu<=0.) {
695 w = +100.;
696 } else {
697 w = log(pb/pu);
698 }
699 return w;
700 }
701
702
703
704 }
705
706
707
| 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.
|
|