Belle II Software light-2406-ragdoll
B2BIIFixMdstModule_muid.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8// -*- C++ -*-
9//
10// Package: B2BIIFixMdst
11// Module: B2BIIFixMdst_muid
12//
13// Description: Probability density provider.
14//
15// Usage: This is the class definition of the object that provides
16// the probability densities used in muon / hadron discrimination.
17//
18// Author: Leo Piilonen
19// Created: 31-Jan-2001
20//
21// $Id: B2BIIFixMdst_muid.cc 9944 2006-11-29 07:36:07Z katayama $
22//
23// Revision history
24//
25// $Log$
26//
27// Revision 2.0 2015/03/11 tkeck
28// Conversion to Belle II
29//
30// Revision 1.1 2002/08/13 05:18:45 hitoshi
31// muid_user is absorped in B2BIIFixMdst (by Piilonen).
32//
33// Revision 1.00 2002/08/08 14:00:00 piilonen
34// Copied from muid_user
35//
36// my include file -- should always be first!
37
38#include <b2bii/modules/B2BIIMdstInput/B2BIIFixMdstModule.h>
39#include "belle_legacy/panther/panther.h"
40#include "belle_legacy/tables/mdst.h"
41#include "belle_legacy/tables/belletdf.h"
42#include "belle_legacy/belleutil/belutil.h"
43#include <framework/utilities/Spline.h>
44
45// system include files
46#include <cmath>
47#include <algorithm>
48#include <vector>
49#include <string>
50#include <iostream>
51#include <fstream>
52#include <cstdio>
53
54using namespace Belle2;
55
56const double MuidProb::kRchisqMax = 10.0 ;
57const double MuidProb::kEEclMax = 1.0 ;
58const double MuidProb::kPTrkMax = 5.0 ;
59
60// Start of muid routines in B2BIIFixMdst __________________________________
61
63{
64
65 m_old_expno = (m_expno != 0) ? m_expno : 5; // expt #5, by default
67 if (m_mapped_expno > 0) {
68 m_muonprob = new MuidProb("muon", m_mapped_expno);
69 m_pionprob = new MuidProb("pion", m_mapped_expno);
70 m_kaonprob = new MuidProb("kaon", m_mapped_expno);
71
72 B2INFO("%B2BIIFixMdst muid is initialized by basf with"
73 << std::endl << " Endcap_MX_layer = " << m_eklm_max_layer
74 << ((m_eklm_max_layer == 0) ?
75 " ==> use value in data file" : "")
76 << ((m_eklm_max_layer == 13) ? " (if possible)" : "")
77 << std::endl << " UseECL = "
78 << ((m_use_ecl == 0) ? "false" : "true")
79 << std::endl << " ExpNo = " << m_expno
80 << ((m_expno == 0) ?
81 " ==> use value in data file" : ""));
82 } else {
83 B2INFO("%B2BIIFixMdst muid is disabled by user (ExpNo = "
84 << m_expno << " is negative)");
85 }
86
87 return;
88
89}
90
92{
93
94 if (m_mapped_expno > 0) {
95 delete m_muonprob;
96 delete m_pionprob;
97 delete m_kaonprob;
98 }
99
100 return;
101
102}
103
104void B2BIIFixMdstModule::Muid_begin_run(int expmc, int exp, int run)
105{
106
107 if (m_mapped_expno > 0) {
108 int expno = (m_expno != 0) ? m_expno : exp;
109 if (expno != m_old_expno) {
110 B2INFO("Reading DB constants for "
111 << ((expmc == 1) ? "data" : "MC") << ", Exp #" << expno
112 << ", Run #" << run);
113 m_old_expno = expno;
118 }
119 }
120
121 return;
122
123}
124
126{
127
128 return;
129
130}
131
133{
134
135 Belle::Mdst_klm_mu_ex_Manager& muexMgr = Belle::Mdst_klm_mu_ex_Manager::get_manager();
136 // In well formed C++ references must always be legal but panther sometimes
137 // returns references to NULL. To catch this reliably we check the address
138 // of the reference and to make sure the compiler does not optimize this
139 // away we cast to a void pointer first so he doesn't know it's the address
140 // of a reference anymore.
141 if ((void*)&muexMgr == NULL) {
142 B2ERROR("%B2BIIFixMdst muid did not find Belle::Mdst_klm_mu_ex table "
143 << "in this event, which implies" << std::endl
144 << " that muid_set and muid_dec were not run yet. "
145 << "Required order of modules is " << std::endl
146 << " ... muid_set muid_dec ... rec2mdst ... B2BIIFixMdst ..."
147 << std::endl << " This program is terminated.");
148 exit(-1);
149 }
150 Belle::Mdst_ecl_trk_Manager& eclTrkMgr = Belle::Mdst_ecl_trk_Manager::get_manager();
151 Belle::Mdst_ecl_trk_Index eclTrkIdx = eclTrkMgr.index("trk");
152 eclTrkIdx.update();
153
154// For each Belle::Mdst_klm_mu_ex, refigure the likelihoods by using the
155// most up-to-date probability density functions of muons, pions, and
156// kaons for the KLM range and reduced chi-squared. Optionally, fold
157// in the probability density function of muons, pions, and kaons for
158// the ECL energy associated with the corresponding charged track.
159// The chain of pointers among panther tables is a little byzantine.
160
161 for (std::vector<Belle::Mdst_klm_mu_ex>::iterator iMuex = muexMgr.begin();
162 iMuex < muexMgr.end(); ++iMuex) {
163 Belle::Mdst_charged& chg = iMuex->pMDST_Charged();
164 Belle::Mdst_muid& muid = chg.muid();
165 int ECMaxLyr = (muid.quality() & 0x400000) ? 1 : 0;
166 if (m_eklm_max_layer == 11) {
167 ECMaxLyr = 1;
168 } else {
169 if ((m_eklm_max_layer == 13) && (ECMaxLyr == 0)) {
170 ECMaxLyr = 0;
171 }
172 }
173 double px = chg.px();
174 double py = chg.py();
175 double pz = chg.pz();
176 double p = sqrt(px * px + py * py + pz * pz);
177 double muon = iMuex->Muon_likelihood();
178 double pion = iMuex->Pion_likelihood();
179 double kaon = iMuex->Kaon_likelihood();
180 double miss = iMuex->Miss_likelihood();
181 double junk = iMuex->Junk_likelihood();
182 int outcome = iMuex->Outcome();
183 int lyrPasMax = iMuex->Layer_trk_brl() +
184 iMuex->Layer_trk_end() + 1;
185 int lyrPicMax = iMuex->Layer_hit_brl();
186 if (iMuex->Layer_hit_end() != -1) {
187 lyrPicMax = iMuex->Layer_trk_brl() +
188 iMuex->Layer_hit_end() + 1;
189 }
190 int lyrDiff = lyrPasMax - lyrPicMax;
191 int lyrExt = iMuex->Layer_trk_brl();
192 if ((outcome == 2) || (outcome == 4)) {
193 lyrExt = iMuex->Layer_trk_end();
194 if ((ECMaxLyr == 1) && (lyrExt > 11)) {
195 outcome = 4; // unconditionally an escape from endcap
196 lyrExt = 11;
197 int lyrEnd = std::min(iMuex->Layer_hit_end(), 11);
198 lyrPasMax = iMuex->Layer_trk_brl() + lyrExt + 1;
199 lyrPicMax = iMuex->Layer_trk_brl() + lyrEnd + 1;
200 lyrDiff = lyrPasMax - lyrPicMax;
201 }
202 }
203 int nPoints = iMuex->N_hits() * 2;
204 double rChi2 = 1.0E10;
205 if (nPoints > 0) {
206 rChi2 = iMuex->Chi_2() / nPoints;
207 }
208
209 if (outcome > 0) {
210 muon = m_muonprob->prob(ECMaxLyr, outcome, lyrExt, lyrDiff, rChi2);
211 pion = m_pionprob->prob(ECMaxLyr, outcome, lyrExt, lyrDiff, rChi2);
212 kaon = m_kaonprob->prob(ECMaxLyr, outcome, lyrExt, lyrDiff, rChi2);
213 double denom = muon + pion + kaon;
214 if (denom < 1.0E-10) {
215 denom = 1.0;
216 muon = pion = kaon = 0.0;
217 junk = 1.0; // should happen very rarely
218 }
219 muon /= denom;
220 pion /= denom;
221 kaon /= denom;
222 }
223 if (m_use_ecl) {
224 Belle::Mdst_trk& trk = chg.trk();
225 std::vector<Belle::Mdst_ecl_trk> eclTrk = point_from(trk.get_ID(), eclTrkIdx);
226 std::vector<Belle::Mdst_ecl_trk>::iterator iEclTrk = eclTrk.begin();
227 if (iEclTrk != eclTrk.end()) {
228 Belle::Mdst_ecl& ecl = iEclTrk->ecl();
229 if (ecl) {
230 if (ecl.match() >= 1) {
231 double eEcl = ecl.energy();
232 if ((eEcl > 0.0) && (p > 0.25)) {
233 double muonE = m_muonprob->probECL(eEcl, p);
234 double pionE = m_pionprob->probECL(eEcl, p);
235 double kaonE = m_kaonprob->probECL(eEcl, p);
236 double denom = muonE + pionE + kaonE;
237 if (denom > 1.0E-10) {
238 if (miss + junk > 0.0) {
239 muon = muonE;
240 pion = pionE;
241 kaon = kaonE;
242 miss = 0.0;
243 junk = 0.0;
244 outcome = 5;
245 } else {
246 if (nPoints == 0) {
247 muon *= muonE;
248 pion *= pionE;
249 kaon *= kaonE;
250 }
251 }
252 }
253 denom = muon + pion + kaon;
254 if (denom < 1.0E-10) {
255 denom = 1.0;
256 muon = pion = kaon = 0.0;
257 junk = 1.0; // should happen very rarely
258 }
259 muon /= denom;
260 pion /= denom;
261 kaon /= denom;
262 } // if ( (eEcl > 0.0) && (p > 0.25) )
263 } // if ( ecl.match() >= 1 )
264 } // if ( ecl )
265 } // if ( iEclTrk != eclTrk.end() )
266 } // if ( m_use_ecl )
267
268 iMuex->Muon_likelihood(muon);
269 iMuex->Pion_likelihood(pion);
270 iMuex->Kaon_likelihood(kaon);
271 iMuex->Miss_likelihood(miss);
272 iMuex->Junk_likelihood(junk);
273 iMuex->Outcome(outcome);
274
275 } // for iMuex loop
276
277 return;
278
279}
280
281// End of muid routines in B2BIIFixMdst ____________________________________
282
283// Start of class MuidProb _____________________________________________
284
285// Default constructor: open database file and read in the tables
286
287MuidProb::MuidProb(const char* dbtemplate, int& expNo)
288{
289
290 readDB(dbtemplate, expNo);
291
292}
293
294void MuidProb::readDB(const char* dbtemplate, int& expNo)
295{
296
297 double dx; // bin size
298
299 if (expNo == 0) return; // do nothing if pdf file was not found previously
300
301 if ((expNo < 5) || ((expNo & 1) == 0)) {
302 expNo = 5; // Runs 1-4 and all even-numbered runs map to Exp #5
303 }
304
305 char dbname[] = "muid_xxxx_e000000.dat";
306 std::string pathname = "";
307 bool tmp = Belle::set_belfnm_verbose(false);
308 while ((pathname == "") && (expNo >= 5)) {
309 std::sprintf(dbname, "%s%s%s%06d%s", "muid_", dbtemplate,
310 "_e", expNo, ".dat");
311 pathname = Belle::belfnm(dbname, 0, "share/belle_legacy/data-files/muid");
312 expNo -= 2;
313 }
314 (void)Belle::set_belfnm_verbose(tmp);
315 expNo += 2;
316 if (pathname == "") {
317 B2ERROR("%MuidProb: Failed to open database file."
318 << "Module B2BIIFixMdst muid will do nothing.");
319 expNo = 0;
320 return;
321 }
322
323 std::ifstream probDB;
324 char ident[90];
325 probDB.open(pathname.c_str());
326 probDB.get(ident, 80, '\n');
327 B2INFO("%MuidProb: Reading file " << pathname
328 << std::endl << "%MuidProb: " << ident
329 << std::endl << "%MuidProb: ");
330
331 for (int l = 0; l < 2; l++) {
332 B2INFO("range ");
333 for (int k = 0; k < 4; k++) {
334 for (int j = 0; j < 15; j++) {
335 for (int i = 0; i < kRange; i++) {
336 probDB >> fRange[l][k][j][i];
337 }
338 }
339 }
340 B2INFO("chi**2 ");
341 dx = kRchisqMax / kRchisq;
342 for (int k = 0; k < 4; k++) {
343 for (int i = 0; i < kRchisq; i++) {
344 probDB >> fRchisq[l][k][i];
345 }
346 Spline::muidSpline(kRchisq, dx, &fRchisq[l][k][0], &fRchisqD1[l][k][0],
347 &fRchisqD2[l][k][0], &fRchisqD3[l][k][0]);
348 }
349 for (int k = 0; k < 4; k++) {
350 for (int j = 0; j < 15; j++) {
351 probDB >> fRchisqN[l][k][j];
352 }
353 }
354 }
355 B2INFO("eEcl ");
356 dx = kEEclMax / kEEcl;
357 for (int k = 0; k < kPTrk; k++) {
358 double probSum = 0.0; // normalization denominator
359 for (int i = 0; i < kEEcl; i++) {
360 probDB >> fEEcl[k][i];
361 probSum += fEEcl[k][i];
362 }
363 fEEcl[k][kEEcl] = 1.0 - probSum; // overflow bin
364 if (fEEcl[k][kEEcl] < 0.0) { fEEcl[k][kEEcl] = 0.0; }
365 if (fEEcl[k][kEEcl] > 1.0) { fEEcl[k][kEEcl] = 1.0; }
366 Spline::muidSpline(kEEcl, dx, &fEEcl[k][0], &fEEclD1[k][0],
367 &fEEclD2[k][0], &fEEclD3[k][0]);
368 }
369
370 probDB.close();
371
372}
373
374//____________________________________________________________________________
375
376double MuidProb::prob(int ECMaxLyr, int outcome,
377 int lyr_ext, int lyr_dif, double chi2_red) const
378{
379
380 // ECMaxLyr: 0 if Endcap_MX_layer == 13; 1 if Endcap_MX_layer == 11
381 // outcome: 0=??, 1=Barrel Stop, 2=Endcap Stop, 3=Barrel Exit, 4=Endcap Exit
382 // lyr_ext: last layer that Ext track touched
383 // lyr_dif: difference between last Ext layer and last hit layer
384 // chi2_red: reduced chi**2 of the transverse deviations of all associated
385 // hits from the corresponding Ext track crossings
386
387 return probRange(ECMaxLyr, outcome, lyr_ext, lyr_dif) *
388 probRchisq(ECMaxLyr, outcome, lyr_ext, chi2_red);
389
390}
391
392//____________________________________________________________________________
393
394double MuidProb::probRange(int ECMaxLyr, int outcome, int lyr_ext,
395 int lyr_dif) const
396{
397
398
399 // ECMaxLyr: 0 if max endcap layer is 13; 1 if it's 11.
400 // outcome: 0=??, 1=Barrel Stop, 2=Endcap Stop, 3=Barrel Exit, 4=Endcap Exit
401 // lyr_ext: last layer that Ext track touched
402 // lyr_dif: difference between last Ext layer and last hit layer
403
404 if (outcome <= 0) { return 0.0; }
405 if (outcome > 4) { return 0.0; }
406 if (lyr_ext < 0) { return 0.0; }
407 if (lyr_ext > 14) { return 0.0; }
408 if (lyr_dif < 0) { return 0.0; }
409 if (lyr_dif >= kRange) { lyr_dif = kRange - 1; }
410
411 return fRange[ECMaxLyr][outcome - 1][lyr_ext][lyr_dif];
412
413}
414
415//____________________________________________________________________________
416
417double MuidProb::probRchisq(int ECMaxLyr, int outcome, int lyr_ext,
418 double chi2_red) const
419{
420
421 // ECMaxLyr: 0 if max endcap layer is 13; 1 if it's 11.
422 // outcome: 0=??, 1=Barrel Stop, 2=Endcap Stop, 3=Barrel Exit, 4=Endcap Exit
423 // lyr_ext: last layer that Ext track touched
424 // chi2_red: reduced chi**2 of the transverse deviations of all associated
425 // hits from the corresponding Ext track crossings
426
427 // Extract the probability density for a particular value of reduced
428 // chi-squared by spline interpolation of the binned histogram values.
429 // This avoids binning artifacts that were seen when the histogram
430 // was sampled directly.
431
432 if (outcome <= 0) { return 0.0; }
433 if (outcome > 4) { return 0.0; }
434 if (lyr_ext < 0) { return 0.0; }
435 if (lyr_ext > 14) { return 0.0; }
436 if (chi2_red < 0.0) { return 0.0; }
437
438 double prob;
439 double area = fRchisqN[ECMaxLyr][outcome - 1][lyr_ext];
440 if (chi2_red >= kRchisqMax) {
441 prob = area * fRchisq[ECMaxLyr][outcome - 1][kRchisq - 1] + (1.0 - area);
442 } else {
443 int i = (int)(chi2_red / (kRchisqMax / kRchisq));
444 double dx = chi2_red - i * (kRchisqMax / kRchisq);
445 prob = fRchisq[ECMaxLyr][outcome - 1][i] +
446 dx * (fRchisqD1[ECMaxLyr][outcome - 1][i] +
447 dx * (fRchisqD2[ECMaxLyr][outcome - 1][i] +
448 dx * fRchisqD3[ECMaxLyr][outcome - 1][i]));
449 prob = (prob < 0.0) ? 0.0 : area * prob;
450 }
451 return prob;
452}
453
454//____________________________________________________________________________
455
456double MuidProb::probECL(double eEcl, double p) const
457{
458
459 // eEcl: ECL energy (GeV)
460 // p: CDC momentum (GeV/c)
461
462 double prob;
463 int jTrk = (p < kPTrkMax) ? (int)(p / (kPTrkMax / kPTrk)) : kPTrk - 1;
464 if (eEcl >= kEEclMax) {
465 prob = fEEcl[jTrk][kEEcl]; // overflow bin
466 } else {
467 int i = (int)(eEcl / (kEEclMax / kEEcl));
468 double dx = eEcl - i * (kEEclMax / kEEcl);
469 prob = fEEcl[jTrk][i] +
470 dx * (fEEclD1[jTrk][i] +
471 dx * (fEEclD2[jTrk][i] +
472 dx * fEEclD3[jTrk][i]));
473 if (prob < 0.0) { prob = 0.0; } // don't let spline fit go negative
474 }
475 return prob;
476
477}
void Muid_end_run(void)
Called when the current run ends.
int m_expno
parameter ExpNo: Experiment number for muid (default=0) 0 ==> use experiment number stored in data fi...
void Muid_init(void)
Initialize the Muid module.
void Muid_begin_run(const int, const int, const int)
Called for each new run.
MuidProb * m_kaonprob
Pointer to kaons' prob-density object.
int m_use_ecl
parameter UseECL: Use (1) or don't use (0) ECL in muid (default=0)
void Muid_event(void)
Called for each event.
MuidProb * m_pionprob
Pointer to pions' prob-density object.
void Muid_term(void)
Terminate the Muid module.
int m_mapped_expno
Mapped value of m_old_exp_no.
MuidProb * m_muonprob
Pointer to muons' prob-density object.
int m_old_expno
Most recently used experiment # in muid.
int m_eklm_max_layer
parameter Endcap_MX_layer: Max layer # of KLM endcap (default=11) 0 ==> use same value as was used by...
Class computes probability density for Muid calculation.
double fRchisq[2][4][kRchisq+1]
Reduced chi-squared pdf.
double fRchisqD2[2][4][kRchisq+1]
Second derivatives of Reduced chi-squared pdf.
static const int kRange
Array size of range.
double fRange[2][4][15][kRange]
Range pdf.
double probRange(int, int, int, int) const
Compute probability density for range.
double probECL(double, double) const
Compute probability density for ECL energy deposit.
double fRchisqD1[2][4][kRchisq+1]
First derivatives of Reduced chi-squared pdf.
double prob(int, int, int, int, double) const
Compute probability density.
double probRchisq(int, int, int, double) const
Compute probability density for reduced chi-squared.
MuidProb(const char *, int &)
Constructor.
static const int kRchisq
Array size of reduced chi-squared.
static const double kEEclMax
Overflow value of ECL energy.
double fEEclD1[kPTrk][kEEcl+1]
First derivatives of ECL energy pdf.
double fEEclD2[kPTrk][kEEcl+1]
Second derivatives of ECL energy pdf.
static const double kPTrkMax
Overflow value of CDC momentum.
double fEEclD3[kPTrk][kEEcl+1]
Third derivatives of ECL energy pdf.
double fRchisqN[2][4][15]
Non-overflow normalization.
static const int kEEcl
Array size of ECL energy.
double fRchisqD3[2][4][kRchisq+1]
Third derivatives of Reduced chi-squared pdf.
static const int kPTrk
Array size of CDC momentum.
static const double kRchisqMax
Overflow value of reduced chi-squared.
void readDB(const char *, int &)
Read in probability density functions from database.
double fEEcl[kPTrk][kEEcl+1]
ECL energy pdf.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24