001
002
003
004
005
006
007 #include "JetTagTools/SoftMuonTag.h"
008
009 #include "CLHEP/Vector/ThreeVector.h"
010 #include "FourMom/P4PxPyPzE.h"
011 #include "JetEvent/Jet.h"
012 #include "GaudiKernel/MsgStream.h"
013 #include "Navigation/NavigationToken.h"
014 #include "GaudiKernel/IToolSvc.h"
015 #include "ITrackToVertex/ITrackToVertex.h"
016
017 #include "JetTagInfo/TruthInfo.h"
018 #include "JetTagInfo/SoftMuonInfo.h"
019 #include "JetTagInfo/SoftLeptonTruthInfo.h"
020 #include "JetTagInfo/SLTrueInfo.h"
021
022 #include "JetTagTools/NewLikelihoodTool.h"
023
024 #include "JetTagTools/HistoHelperRoot.h"
025
026 #include "GaudiKernel/ITHistSvc.h"
027 #include "JetTagCalibration/CalibrationBroker.h"
028
029 #include "muonEvent/Muon.h"
030 #include "MuonIDEvent/MuonAssociation.h"
031
032 #include "StoreGate/StoreGateSvc.h"
033
034 namespace Analysis
035 {
036
037 SoftMuonTag::SoftMuonTag(const std::string& t, const std::string& n, const IInterface* p)
038 : AlgTool(t,n,p),
039 m_trackToVertexTool("Reco::TrackToVertex"),
040 m_likelihoodTool("Analysis::NewLikelihoodTool"),
041 m_histoHelper(0)
042 {
043 declareInterface<ITagTool>(this);
044 declareProperty("Runmodus", m_runModus = "analysis");
045 declareProperty("RefType", m_refType = "ALL");
046 declareProperty("TaggingAlgType", m_algMode = "L1D");
047 declareProperty("BTagJetPTmin", m_pTjetmin = 15.*GeV);
048 declareProperty("BTagJetEtamin", m_etajetmin = 2.7);
049 declareProperty("LikelihoodTool", m_likelihoodTool);
050 declareProperty("TrackToVertexTool", m_trackToVertexTool);
051 declareProperty("checkOverflows", m_checkOverflows = false);
052 declareProperty("purificationDeltaR", m_purificationDeltaR = 0.8);
053 declareProperty("muonIsolDeltaR", m_muonIsolDeltaR = 0.7);
054 declareProperty("UseBinInterpol", m_UseBinInterpol = true);
055
056 declareProperty("jetCollectionList", m_jetCollectionList);
057 m_jetCollectionList.push_back("Cone4CalTower");
058 m_jetCollectionList.push_back("Cone7CalTower");
059 m_jetCollectionList.push_back("Kt4CalTower");
060 m_jetCollectionList.push_back("Kt6CalTower");
061 m_jetCollectionList.push_back("Cone4Topo");
062 m_jetCollectionList.push_back("ConeTopo");
063 m_jetCollectionList.push_back("Kt4Topo");
064 m_jetCollectionList.push_back("Kt6Topo");
065 declareProperty("useForcedCalibration", m_doForcedCalib = false);
066 declareProperty("ForcedCalibrationName", m_ForcedCalibName = "Cone4CalTower");
067
068 declareProperty("RecAlgorithm", m_alg = 1);
069 declareProperty("CutD0", m_d0cut = 4.*mm);
070 declareProperty("CutPT", m_pTcut = 4.*GeV);
071 declareProperty("CutDR", m_DRcut = 0.5);
072 declareProperty("CutMatchChi2", m_MatchChi2cut = 10);
073
074 declareProperty("writeInfoPlus", m_writeInfoPlus = true);
075
076 declareProperty("originalMuCollectionName", m_originalMuCollectionName = "StacoMuonCollection");
077
078
079 m_hypothese.push_back("SoftMb");
080 m_hypothese.push_back("SoftMl");
081 m_hypothese.push_back("SoftMc");
082
083
084 m_histoname.push_back("pT");
085 m_histoname.push_back("pTrel");
086 m_histoname.push_back("pTLowPt");
087 m_histoname.push_back("pTrelLowPt");
088 m_histoname.push_back("pTpTrel");
089 m_histoname.push_back("pTpTrelLowPt");
090 m_histoname.push_back("JetETEffL1D");
091 m_histoname.push_back("JetETLowPtEffL1D");
092 m_histoname.push_back("JetETNormL1D");
093 m_histoname.push_back("JetETLowPtNormL1D");
094 m_histoname.push_back("JetETEffL1DL1D");
095 m_histoname.push_back("JetETLowPtEffL1DL1D");
096 m_histoname.push_back("JetETNormL1DL1D");
097 m_histoname.push_back("JetETLowPtNormL1DL1D");
098 m_histoname.push_back("JetETEffL2D");
099 m_histoname.push_back("JetETLowPtEffL2D");
100 m_histoname.push_back("JetETNormL2D");
101 m_histoname.push_back("JetETLowPtNormL2D");
102
103 }
104
105 SoftMuonTag::~SoftMuonTag()
106 {
107 }
108
109 StatusCode SoftMuonTag::initialize()
110 {
111 MsgStream mLog(msgSvc(), name());
112
113 mLog << MSG::INFO << "#BTAG# Initializing..." << endreq;
114 m_printParameterSettings();
115
116
117 IToolSvc* toolSvc;
118 StatusCode sc = service("ToolSvc", toolSvc);
119 if (StatusCode::SUCCESS != sc) {
120 mLog << MSG::ERROR << "#BTAG# Can't get ToolSvc" << endreq;
121 return StatusCode::FAILURE;
122 }
123
124
125 if ( m_trackToVertexTool.retrieve().isFailure() ) {
126 mLog << MSG::FATAL << "#BTAG# Failed to retrieve tool " << m_trackToVertexTool << endreq;
127 return StatusCode::FAILURE;
128 } else {
129 mLog << MSG::INFO << "#BTAG# Retrieved tool " << m_trackToVertexTool << endreq;
130 }
131
132
133
134
135 if (m_doForcedCalib) {
136 if (std::find( m_jetCollectionList.begin(),
137 m_jetCollectionList.end(),
138 m_ForcedCalibName ) == m_jetCollectionList.end()) {
139 mLog << MSG::ERROR << "#BTAG# Error, forced calibration to an unloaded one" << endreq;
140 return StatusCode::FAILURE;
141 }
142 }
143
144
145
146
147 if (m_runModus == "analysis") {
148 mLog << MSG::INFO <<"#BTAG# Reading histos..."<<endreq;
149
150 if ( m_likelihoodTool.retrieve().isFailure() ) {
151 mLog << MSG::FATAL << "#BTAG# Failed to retrieve tool " << m_likelihoodTool << endreq;
152 return StatusCode::FAILURE;
153 } else {
154 mLog << MSG::INFO << "#BTAG# Retrieved tool " << m_likelihoodTool << endreq;
155 }
156
157 if (3!=m_hypothese.size()) {
158 mLog << MSG::ERROR << "#BTAG# There should be 3 hypotheses instead of " << m_hypothese.size() << endreq;
159 return StatusCode::FAILURE;
160 }
161 m_likelihoodTool->defineHypotheses(m_hypothese);
162 for(uint ih=0;ih<m_hypothese.size();ih++) {
163
164 std::string hName = m_hypothese[ih]+"/JetETEff"+m_algMode;
165 m_likelihoodTool->defineHistogram(hName);
166 hName = m_hypothese[ih]+"/JetETNorm"+m_algMode;
167 m_likelihoodTool->defineHistogram(hName);
168 hName = m_hypothese[ih]+"/JetETLowPtEff"+m_algMode;
169 m_likelihoodTool->defineHistogram(hName);
170 hName = m_hypothese[ih]+"/JetETLowPtNorm"+m_algMode;
171 m_likelihoodTool->defineHistogram(hName);
172
173 std::string hBase = m_hypothese[ih]+"/";
174 if(m_algMode == "L1D"){
175 m_likelihoodTool->defineHistogram(hBase+"pTrel");
176 m_likelihoodTool->defineHistogram(hBase+"pTrelLowPt");
177 }
178 else if(m_algMode == "L1DL1D"){
179 m_likelihoodTool->defineHistogram(hBase+"pT");
180 m_likelihoodTool->defineHistogram(hBase+"pTLowPt");
181 m_likelihoodTool->defineHistogram(hBase+"pTrel");
182 m_likelihoodTool->defineHistogram(hBase+"pTrelLowPt");
183 }
184 else if(m_algMode == "L2D"){
185 m_likelihoodTool->defineHistogram(hBase+"pTpTrel");
186 m_likelihoodTool->defineHistogram(hBase+"pTpTrelLowPt");
187 }
188 }
189
190
191
192
193 m_likelihoodTool->printStatus();
194 }
195
196
197
198
199 if (m_runModus=="reference") {
200
201
202
203 ITHistSvc* myHistoSvc;
204 if( service( "THistSvc", myHistoSvc ).isSuccess() ) {
205 mLog << MSG::DEBUG << "#BTAG# "<< name() << ": HistoSvc loaded successfully." << endreq;
206 m_histoHelper = new HistoHelperRoot(myHistoSvc);
207 m_histoHelper->setCheckOverflows(m_checkOverflows);
208 } else mLog << MSG::ERROR << "#BTAG# " << name() << ": HistoSvc could NOT bo loaded." << endreq;
209
210 mLog << MSG::INFO <<"#BTAG# Booking histos..."<<endreq;
211 for(uint ijc=0;ijc<m_jetCollectionList.size();ijc++) {
212 for(uint ih=0;ih<m_hypothese.size();ih++) {
213 if(ih==0) {
214
215 std::string hDir = "/RefFileSoftMu"+m_jetCollectionList[ijc]+"/controlSoftMu/";
216 m_histoHelper->bookHisto(hDir+"eta","eta",60,-3.,3.);
217 m_histoHelper->bookHisto(hDir+"phi","phi",64,-3.2,3.2);
218 m_histoHelper->bookHisto(hDir+"pt","pt",50,0.,300.);
219 m_histoHelper->bookHisto(hDir+"smpt","Soft Muon pT",100,0.,100.);
220 }
221 std::string hDir = "/RefFileSoftMu"+m_jetCollectionList[ijc]+"/"+m_hypothese[ih]+"/";
222
223 m_histoHelper->bookHisto(hDir+"pT", "pT/(pT+5)",100,0.,1.);
224 m_histoHelper->bookHisto(hDir+"pTrel", "pTrel/(pTrel+0.5)",100,0.,1.);
225 m_histoHelper->bookHisto(hDir+"pTLowPt", "pT/(pT+5)",100,0.,1.);
226 m_histoHelper->bookHisto(hDir+"pTrelLowPt", "pTrel/(pTrel+0.5)",100,0.,1.);
227 m_histoHelper->bookHisto(hDir+"pTpTrel", "pT/(pT+5) vs pTrel/(pTrel+0.5)" ,100,0.,1.,100,0.,1.);
228 m_histoHelper->bookHisto(hDir+"pTpTrelLowPt", "pT/(pT+5) vs pTrel/(pTrel+0.5)" ,100,0.,1.,100,0.,1.);
229
230
231 m_histoHelper->bookHisto(hDir+"JetETEffL1D", "Jet Et",100,0,500);
232 m_histoHelper->bookHisto(hDir+"JetETLowPtEffL1D", "Jet Et",100,0,500);
233 m_histoHelper->bookHisto(hDir+"JetETNormL1D", "Jet Et",100,0,500);
234 m_histoHelper->bookHisto(hDir+"JetETLowPtNormL1D", "Jet Et",100,0,500);
235
236 m_histoHelper->bookHisto(hDir+"JetETEffL1DL1D", "Jet Et",100,0,500);
237 m_histoHelper->bookHisto(hDir+"JetETLowPtEffL1DL1D", "Jet Et",100,0,500);
238 m_histoHelper->bookHisto(hDir+"JetETNormL1DL1D", "Jet Et",100,0,500);
239 m_histoHelper->bookHisto(hDir+"JetETLowPtNormL1DL1D", "Jet Et",100,0,500);
240
241 m_histoHelper->bookHisto(hDir+"JetETEffL2D", "Jet Et",100,0,500);
242 m_histoHelper->bookHisto(hDir+"JetETLowPtEffL2D", "Jet Et",100,0,500);
243 m_histoHelper->bookHisto(hDir+"JetETNormL2D", "Jet Et",100,0,500);
244 m_histoHelper->bookHisto(hDir+"JetETLowPtNormL2D", "Jet Et",100,0,500);
245 }
246 }
247 m_histoHelper->print();
248
249 }
250
251 return StatusCode::SUCCESS;
252 }
253
254
255 StatusCode SoftMuonTag::finalize()
256 {
257 MsgStream mLog(msgSvc(), name());
258 return StatusCode::SUCCESS;
259 }
260
261 void SoftMuonTag::tagJet(Jet& jetToTag)
262 {
263 MsgStream mLog(msgSvc(), name());
264
265 mLog << MSG::DEBUG << "#BTAG# Starting tagJet" << endreq;
266
267
268 std::string author = jetToTag.jetAuthor();
269 if (m_doForcedCalib) {
270 author = m_ForcedCalibName;
271 } else {
272
273 if (std::find( m_jetCollectionList.begin(),
274 m_jetCollectionList.end(),
275 author ) == m_jetCollectionList.end()) {
276 mLog << MSG::DEBUG << "#BTAG# Jet Algorithm " << author << " not found in the standard list" << endreq;
277 mLog << MSG::DEBUG << "#BTAG# Trying to find a similar one..." << endreq;
278 if (author.find("Cone4",0) != std::string::npos) author = "Cone4CalTower";
279 else if (author.find("Cone7",0) != std::string::npos) author = "Cone7CalTower";
280 else if (author.find("Kt4",0) != std::string::npos) author = "Kt4CalTower";
281 else if (author.find("Kt6",0) != std::string::npos) author = "Kt6CalTower";
282 else {
283 mLog << MSG::DEBUG << "#BTAG# None found, taking " << m_ForcedCalibName << " calibration" << endreq;
284 author = m_ForcedCalibName;
285 }
286 }
287 }
288
289
290 double jeteta = jetToTag.eta(), jetphi = jetToTag.phi(), jetpt = jetToTag.pt();
291 mLog << MSG::DEBUG << "#BTAG# Jet properties : eta = "<<jeteta
292 <<" phi = "<<jetphi
293 <<" pT = "<<jetpt/1.e3
294 <<" author = " <<author
295 << endreq;
296
297
298 std::string pref = "";
299 std::string label = "N/A";
300 std::vector<Hep3Vector> hardMus;
301 if (m_runModus=="reference") {
302
303 bool hasHardMu(false);
304 const SoftLeptonTruthInfo* sltinfo = jetToTag.tagInfo<SoftLeptonTruthInfo>("SoftLeptonTruthInfo");
305 if (sltinfo) {
306 int nslt = sltinfo->numSLTrueInfo();
307 mLog << MSG::DEBUG << "#BTAG# SL truth info exist. Found " << nslt << " true leptons in jet" << endreq;
308 for (int islt = 0; islt < nslt; islt++) {
309 const SLTrueInfo slt = sltinfo->getSLTrueInfo(islt);
310 mLog << MSG::DEBUG << "#BTAG# SLT info " << slt.pdgId()
311 << " " << slt.momentum().perp()
312 << " " << slt.FromB() << " " << slt.FromD() << " " << slt.FromGH()
313 << endreq;
314 if ( abs(slt.pdgId()) == 13 &&
315 !(slt.FromB()) &&
316 !(slt.FromD()) &&
317 slt.FromGH()
318 ) {
319 Hep3Vector v(slt.momentum());
320 hardMus.push_back(v);
321 double dR = v.deltaR(jetToTag.hlv());
322 mLog << MSG::DEBUG << "#BTAG# DR info "
323 << v.eta() << " " << jetToTag.eta() << " "
324 << v.phi() << " " << jetToTag.phi() << " "
325 << dR
326 << endreq;
327 if(dR<m_muonIsolDeltaR) {
328 hasHardMu = true;
329 }
330 }
331 }
332 }
333
334 if(hasHardMu)return;
335
336 m_histoHelper->fillHisto("/RefFileSoftMu"+author+"/controlSoftMu/eta",(double)jeteta);
337 if (fabs(jeteta) <= m_etajetmin) {
338 m_histoHelper->fillHisto("/RefFileSoftMu"+author+"/controlSoftMu/phi",(double)jetphi);
339 m_histoHelper->fillHisto("/RefFileSoftMu"+author+"/controlSoftMu/pt",(double)jetpt/1.e3);
340
341 if(jetpt>m_pTjetmin)
342 {
343 const TruthInfo* mcTrueInfo = jetToTag.tagInfo<TruthInfo>("TruthInfo");
344 if (mcTrueInfo) {
345 label = mcTrueInfo->jetTruthLabel();
346 mLog << MSG::DEBUG << "#BTAG# Jet label="<<label<< endreq;
347
348 double deltaRtoClosestB = mcTrueInfo->deltaRMinTo("B");
349 double deltaRtoClosestC = mcTrueInfo->deltaRMinTo("C");
350 double deltaRmin = deltaRtoClosestB < deltaRtoClosestC ? deltaRtoClosestB : deltaRtoClosestC;
351
352 if ( ( "B"==m_refType && "B"==label ) ||
353 ( "C"==m_refType && "C"==label ) ||
354 ( "L"==m_refType && "N/A"==label ) ||
355 ( "ALL"==m_refType &&
356 (
357 ( "B"==label && deltaRtoClosestC > m_purificationDeltaR ) ||
358 ( "C"==label && deltaRtoClosestB > m_purificationDeltaR ) ||
359 ( "N/A"==label && deltaRmin > m_purificationDeltaR ) ) )
360 ) {
361 if ("B"==label) {
362 pref = m_hypothese[0];
363 } else if ("N/A"==label) {
364 pref = m_hypothese[1];
365 } else if ("C"==label) {
366 pref = m_hypothese[2];
367 }
368 std::string hDir = "/RefFileSoftMu"+author+"/"+pref+"/";
369 m_histoHelper->fillHisto(hDir+"JetETNormL1D",(double)jetpt/1.e3);
370 m_histoHelper->fillHisto(hDir+"JetETNormL1DL1D",(double)jetpt/1.e3);
371 m_histoHelper->fillHisto(hDir+"JetETNormL2D",(double)jetpt/1.e3);
372 m_histoHelper->fillHisto(hDir+"JetETLowPtNormL1D",(double)jetpt/1.e3);
373 m_histoHelper->fillHisto(hDir+"JetETLowPtNormL1DL1D",(double)jetpt/1.e3);
374 m_histoHelper->fillHisto(hDir+"JetETLowPtNormL2D",(double)jetpt/1.e3);
375 }
376 else return;
377 }
378 else {
379 mLog << MSG::ERROR << "#BTAG# No Label ! Cannot run on reference mode !"<<endreq;
380 return;
381 }
382 }
383 else return;
384 }
385 else return;
386 }
387
388 if(m_runModus=="analysis" && m_writeInfoPlus) {
389 bool ppb = true;
390 StoreGateSvc* m_StoreGate;
391 StatusCode sc = service("StoreGateSvc", m_StoreGate);
392 if (sc.isFailure()) {
393 mLog << MSG::ERROR << "#BTAG# StoreGate service not found !" << endreq;
394 } else {
395 sc = m_StoreGate->retrieve(m_originalMuCollection, m_originalMuCollectionName);
396 if (sc.isFailure()) {
397 mLog << MSG::WARNING << "#BTAG# " << m_originalMuCollectionName << " not found in StoreGate." << endreq;
398 } else {
399 mLog << MSG::DEBUG << "#BTAG# MuonContainer " << m_originalMuCollectionName << " found." << endreq;
400 ppb = false;
401 }
402 }
403 if(ppb) {
404 mLog << MSG::WARNING << "#BTAG# Not able to persistify infos ! Exiting..." << endreq;
405 return;
406 }
407 }
408
409 const MuonAssociation *mc = jetToTag.getAssociation<MuonAssociation>("Muons");
410 if (mc == 0) {
411 mLog << MSG::INFO << "#BTAG# No muon constituent" << endreq;
412 return;
413 }
414
415 SoftMuonInfo* softmInfo(0);
416
417 int muCounter = 0;
418 int muCounterL1D = 0;
419 int muCounterLowPt = 0;
420 int muCounterL1DLowPt = 0;
421 std::vector<double> bestMuProbi;
422 double bestMuWeight = 0;
423 double highestPT = 0;
424 double highestPTrel = 0;
425 bool highestIsLowP = 0;
426
427 for(Navigable<MuonContainer,double>::object_iter it = mc->begin(); it!= mc->end(); ++it) {
428 const Muon *m = (*it);
429 if (m != 0) {
430
431 bool closeToHardMu(false);
432 for(uint i=0;i<hardMus.size();i++)
433 {
434 Hep3Vector v = hardMus[i];
435 double dR = v.deltaR(m->hlv());
436 mLog << MSG::DEBUG << "#BTAG# DR(mu-mu) info "
437 << v.eta() << " " << m->eta() << " "
438 << v.phi() << " " << m->phi() << " "
439 << dR
440 << endreq;
441 if(dR<0.1) {
442 closeToHardMu = true;
443 }
444 }
445 if(closeToHardMu)
446 {
447 mLog << MSG::DEBUG << "#BTAG# skipping this muon" << endreq;
448 continue;
449 }
450
451 const MuonParameters::Author mAuthor = m->author();
452 mLog << MSG::DEBUG << "#BTAG# Muon Author=" << mAuthor << " "
453 << MuonParameters:: mediumPt << " "
454 << endreq;
455
456 if(mAuthor == MuonParameters::mediumPt ) continue;
457 bool isComb(m->isCombinedMuon());
458 bool isLowP(m->isLowPtReconstructedMuon());
459 bool acceptAlg(true);
460 if( 0==m_alg && ( isLowP || !isComb ) ) acceptAlg = false;
461 else if( 1==m_alg && (!isLowP && !isComb) ) acceptAlg = false;
462 if(!acceptAlg)continue;
463
464 mLog << MSG::DEBUG << "#BTAG# Muon Match Chi2=" << m->matchChi2OverDoF() << " " << m_MatchChi2cut << endreq;
465 if( m->matchChi2OverDoF()>m_MatchChi2cut && (isComb||isLowP) )continue;
466
467 bool passD0cut(true);
468 double d0wrtPriVtx(999.);
469 if( isLowP || isComb ) {
470 const Rec::TrackParticle* trk = m->inDetTrackParticle();
471 d0wrtPriVtx = trk->measuredPerigee()->parameters()[Trk::d0];
472 const Trk::MeasuredPerigee* perigee =
473 m_trackToVertexTool->perigeeAtVertex(*trk, m_priVtx->recVertex().position());
474 if (perigee) {
475 d0wrtPriVtx = perigee->parameters()[Trk::d0];
476 delete perigee;
477 }
478 if(fabs(d0wrtPriVtx)>m_d0cut)passD0cut = false;
479 }
480 if(!passD0cut)continue;
481
482 double dR = jetToTag.hlv().deltaR(m->hlv());
483 double pt(1./m->iPt());
484 double ptrel = m->hlv().perp((jetToTag.hlv()+m->hlv()).vect());
485
486 double ptN ( pt / (pt +5.*GeV) );
487 double ptrelN ( ptrel / (ptrel+ 0.5*GeV) );
488
489 mLog << MSG::DEBUG << "#BTAG# Found muon isLowP=" << isLowP << " isComb=" << isComb
490 << " acceptAlg=" << acceptAlg
491 << " pt=" << pt << "(cut=" << m_pTcut << ") "
492 << " d0=" << d0wrtPriVtx << "(cut=" << m_d0cut << ") "
493 << " dR=" << dR << "(cut=" << m_DRcut << ") "
494 << " ptrel=" << ptrel
495 << " cuts=" << m_d0cut << " " << m_DRcut << " " << m_pTcut
496 << endreq;
497
498 if(dR<m_DRcut)
499 {
500 if(isLowP)muCounterLowPt++;
501 else muCounter++;
502 if(m_runModus == "reference")
503 {
504 if (jetpt >= m_pTjetmin && fabs(jeteta) <= m_etajetmin) {
505 std::string hDir = "/RefFileSoftMu"+author+"/"+pref+"/";
506 if(1==muCounter){
507 m_histoHelper->fillHisto(hDir+"JetETEffL1DL1D",(double)jetpt/1.e3);
508 m_histoHelper->fillHisto(hDir+"JetETEffL2D",(double)jetpt/1.e3);
509 }
510 if(1==muCounterLowPt){
511 m_histoHelper->fillHisto(hDir+"JetETLowPtEffL1DL1D",(double)jetpt/1.e3);
512 m_histoHelper->fillHisto(hDir+"JetETLowPtEffL2D",(double)jetpt/1.e3);
513 }
514 if(pt>highestPT)
515 {
516 highestPT = pt;
517 highestPTrel = ptrel;
518 highestIsLowP = isLowP;
519 }
520 }
521 }
522 else if(m_runModus == "analysis")
523 {
524 std::vector<Slice> slices;
525 AtomicProperty atom2(ptrelN,"pTrel");
526 if(m_algMode == "L2D")
527 {
528 AtomicProperty atom1(ptN,"pT");
529 Slice slice("L2D");
530 if(isLowP){
531 Composite compo(author+"#pTpTrelLowPt");
532 compo.atoms.push_back(atom1);
533 compo.atoms.push_back(atom2);
534 slice.composites.push_back(compo);
535 }
536 else{
537 Composite compo(author+"#pTpTrel");
538 compo.atoms.push_back(atom1);
539 compo.atoms.push_back(atom2);
540 slice.composites.push_back(compo);
541 }
542 slices.push_back(slice);
543 }
544 else if(m_algMode == "L1DL1D")
545 {
546 AtomicProperty atom1(ptN,"pT");
547 Slice slice("L1DL1D");
548 if(isLowP){
549 Composite compo1(author+"#pTLowPt");
550 compo1.atoms.push_back(atom1);
551 slice.composites.push_back(compo1);
552 Composite compo2(author+"#pTrelLowPt");
553 compo2.atoms.push_back(atom2);
554 slice.composites.push_back(compo2);
555 }
556 else{
557 Composite compo1(author+"#pT");
558 compo1.atoms.push_back(atom1);
559 slice.composites.push_back(compo1);
560 Composite compo2(author+"#pTrel");
561 compo2.atoms.push_back(atom2);
562 slice.composites.push_back(compo2);
563 }
564 slices.push_back(slice);
565 }
566 else if(m_algMode == "L1D" && fabs(pt)>m_pTcut)
567 {
568 Slice slice("L1D");
569 if(isLowP){
570 muCounterL1DLowPt++;
571 Composite compo(author+"#pTrelLowPt");
572 compo.atoms.push_back(atom2);
573 slice.composites.push_back(compo);
574 }
575 else{
576 muCounterL1D++;
577 Composite compo(author+"#pTrel");
578 compo.atoms.push_back(atom2);
579 slice.composites.push_back(compo);
580 }
581 slices.push_back(slice);
582 }
583 std::vector<double> probi;
584 if(slices.size())
585 {
586 m_likelihoodTool->setLhVariableValue(slices);
587 probi = m_likelihoodTool->calculateLikelihood();
588 double w(1);
589 if (probi.size() >= 3)
590 {
591 mLog << MSG::VERBOSE << "#BTAG# SoftMu probabilities "
592 <<" p_b = "<<probi[0]
593 <<" p_c = "<<probi[2]
594 <<" p_l = "<<probi[1]
595 << endreq;
596 double effb = 0.5, effl = 0.5, effc = 0.5;
597 if(isLowP){
598 effb = m_likelihoodTool->getEff(m_hypothese[0],author+"#JetETLowPt",m_algMode);
599 effl = m_likelihoodTool->getEff(m_hypothese[1],author+"#JetETLowPt",m_algMode);
600 effc = m_likelihoodTool->getEff(m_hypothese[2],author+"#JetETLowPt",m_algMode);
601 mLog << MSG::VERBOSE << "#BTAG# SoftMu Low PT efficiencies for jetColl "<<author
602 <<" eps_b = "<<effb<<" eps_c = "<<effc<<" eps_l = "<<effl
603 << endreq;
604
605 }
606 else{
607 effb = m_likelihoodTool->getEff(m_hypothese[0],author+"#JetET",m_algMode);
608 effl = m_likelihoodTool->getEff(m_hypothese[1],author+"#JetET",m_algMode);
609 effc = m_likelihoodTool->getEff(m_hypothese[2],author+"#JetET",m_algMode);
610 mLog << MSG::VERBOSE << "#BTAG# SoftMu efficiencies for jetColl "<<author
611 <<" eps_b = "<<effb<<" eps_c = "<<effc<<" eps_l = "<<effl
612 << endreq;
613 }
614 probi[0] *= effb;
615 probi[1] *= effl;
616 probi[2] *= effc;
617 w = probi[0];
618 }
619 else
620 {
621 mLog << MSG::ERROR << "#BTAG# Missing number in jet probabilities ! "<<probi.size()<<endreq;
622 }
623 if(w>bestMuWeight)
624 {
625 bestMuWeight = w;
626 bestMuProbi = probi;
627 }
628 if(0==softmInfo)
629 {
630
631 std::string instanceName(name());
632 softmInfo = new SoftMuonInfo(instanceName.erase(0,8));
633 jetToTag.addInfo(softmInfo);
634 }
635
636 if(m_writeInfoPlus)
637 {
638 SMTrackInfo tinfo(m_originalMuCollection,m,d0wrtPriVtx,ptrel,probi);
639 softmInfo->addTrackInfo(tinfo);
640 }
641 m_likelihoodTool->clear();
642 }
643 }
644 }
645 }
646 mLog << MSG::DEBUG << "#BTAG# Done with muon " << m <<endreq;
647 }
648
649
650 if (muCounter+muCounterLowPt<1) {
651 mLog << MSG::DEBUG << "#BTAG# Jet does not contain any good muons" << endreq;
652 return;
653 }
654 else mLog << MSG::DEBUG << "#BTAG# Jet contains " << muCounter+muCounterLowPt << " muons" << endreq;
655
656 if (m_runModus == "reference" && highestPT>0)
657 {
658 std::string hDir = "/RefFileSoftMu"+author+"/"+pref+"/";
659 double ptN ( highestPT / (highestPT +5.*GeV) );
660 double ptrelN ( highestPTrel / (highestPTrel+ 0.5*GeV) );
661 if(highestIsLowP){
662 m_histoHelper->fillHisto(hDir+"pTLowPt",ptN);
663 m_histoHelper->fillHisto(hDir+"pTpTrelLowPt",ptN,ptrelN);
664 }
665 else{
666 m_histoHelper->fillHisto(hDir+"pT",ptN);
667 m_histoHelper->fillHisto(hDir+"pTpTrel",ptN,ptrelN);
668 }
669 if(fabs(highestPT)>m_pTcut){
670 if(highestIsLowP){
671 m_histoHelper->fillHisto(hDir+"JetETLowPtEffL1D",(double)jetpt/1.e3);
672 m_histoHelper->fillHisto(hDir+"pTrelLowPt",ptrelN);
673 }
674 else{
675 m_histoHelper->fillHisto(hDir+"JetETEffL1D",(double)jetpt/1.e3);
676 m_histoHelper->fillHisto(hDir+"pTrel",ptrelN);
677 }
678 hDir = "/RefFileSoftMu"+author+"/controlSoftMu/";
679 m_histoHelper->fillHisto(hDir+"smpt",highestPT*1e-3);
680 }
681 }
682 else if (softmInfo)
683 {
684 softmInfo->setTagLikelihood(bestMuProbi);
685 double w(-100);
686 if(bestMuProbi.size()>1)
687 {
688 double pb = bestMuProbi[0];
689 double pu = bestMuProbi[1];
690 if(pb<=0.) {
691 w = -30.;
692 } else if (pu<=0.) {
693 w = +100.;
694 } else {
695 w = log(pb/pu);
696 }
697 }
698 softmInfo->setWeight(w);
699
700 softmInfo->makeValid();
701 }
702
703 mLog << MSG::DEBUG << "#BTAG# tagJet is done" << endreq;
704
705 return;
706 }
707
708 void SoftMuonTag::finalizeHistos()
709 {
710 if (m_runModus=="reference") {
711 for(uint ijc=0;ijc<m_jetCollectionList.size();ijc++) {
712 for(uint ih=0;ih<m_hypothese.size();ih++) {
713 std::string hDir = "/BTAG/CALIB/"+m_hypothese[ih]+"/";
714 m_likelihoodTool
715 ->smoothAndNormalizeHistogram( m_histoHelper->getHisto1D(hDir+"pT") , "" );
716 m_likelihoodTool
717 ->smoothAndNormalizeHistogram( m_histoHelper->getHisto1D(hDir+"pTrel") , "" );
718 m_likelihoodTool
719 ->smoothAndNormalizeHistogram( m_histoHelper->getHisto1D(hDir+"pTLowPt") , "" );
720 m_likelihoodTool
721 ->smoothAndNormalizeHistogram( m_histoHelper->getHisto1D(hDir+"pTrelLowPt") , "" );
722 m_likelihoodTool
723 ->smoothAndNormalizeHistogram( m_histoHelper->getHisto1D(hDir+"pTpTrel") , "" );
724 m_likelihoodTool
725 ->smoothAndNormalizeHistogram( m_histoHelper->getHisto1D(hDir+"pTpTrelLowPt") , "" );
726 }
727 }
728 }
729 return;
730 }
731
732 void SoftMuonTag::m_printParameterSettings()
733 {
734 MsgStream log(msgSvc(), name());
735 log << MSG::INFO << "#BTAG# " << name() << "Parameter settings " << endreq;
736 log << MSG::INFO << "#BTAG# I am in " << m_runModus << " modus." << endreq;
737 log << MSG::INFO << "#BTAG# The method is "<<m_algMode<< endreq;
738 if (m_runModus == "reference") {
739 log << MSG::INFO << "#BTAG# Preparing "<< m_refType<< "-jet probability density functions..."<<endreq;
740 }
741 }
742
743 }
744
| 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.
|
|