137 MsgStream&
msg)
const {
139 double mantissa[4] = {0};
148 msg << MSG::VERBOSE <<
"called JetSignalStateCnv::compress() for :" <<
endmsg;
149 msg << MSG::VERBOSE <<
" raw momentum ( px | py | pz | m ) : ( "
150 << momRaw.
m_px <<
" | "
151 << momRaw.
m_py <<
" | "
152 << momRaw.
m_pz <<
" | "
154 msg << MSG::VERBOSE <<
" cal momentum ( px | py | pz | m ) : ( "
155 << momCal.
m_px <<
" | "
156 << momCal.
m_py <<
" | "
157 << momCal.
m_pz <<
" | "
160 bool forceNoCompression =
false;
163 if (momRaw.
m_m == 0 || momCal.
m_m == 0) {
164 forceNoCompression =
true;
166 if (momRaw.
m_px == 0 && momRaw.
m_py == 0) {
167 forceNoCompression =
true;
169 if (momCal.
m_px == 0 && momCal.
m_py == 0) {
170 forceNoCompression =
true;
173 if (forceNoCompression) {
174 msg << MSG::DEBUG <<
"M or PT of calibrated or raw signal state of jet zero !"
175 <<
" Switching off compression for this jet !!!" <<
endmsg;
181 double angleC = std::atan2(momCal.
m_py, momCal.
m_px);
182 double angleR = std::atan2(momRaw.
m_py, momRaw.
m_px);
183 double p_traC = std::sqrt(momCal.
m_px * momCal.
m_px + momCal.
m_py * momCal.
m_py);
184 double p_traR = std::sqrt(momRaw.
m_px * momRaw.
m_px + momRaw.
m_py * momRaw.
m_py);
185 mantissa[0] = std::frexp(momCal.
m_m / momRaw.
m_m - 1, &exponent[0]);
186 mantissa[1] = std::frexp(angleC - angleR, &exponent[1]);
187 mantissa[2] = std::frexp(momCal.
eta() - momRaw.
eta(), &exponent[2]);
188 mantissa[3] = std::frexp(p_traC / p_traR - 1, &exponent[3]);
190 delta0 = momCal.
m_m / momRaw.
m_m - 1;
191 delta1 = angleC - angleR;
192 delta2 = momCal.
eta() - momRaw.
eta();
193 delta3 = p_traC / p_traR - 1;
196 unsigned long tmp0(0);
197 unsigned long tmp1(0);
198 unsigned long tmp2(0);
199 unsigned long tmp3(0);
201 unsigned short vec0(0);
202 unsigned short vec1(0);
203 unsigned short vec2(0);
209 if (exponent[0] > 0x7) {
212 if (exponent[0] < 0x0) {
213 factor = -exponent[0];
216 tmp0 = int(std::fabs(mantissa[0]) * 0x10) >> factor;
217 tmp0 |= (exponent[0] & 0x7) << 5;
218 if (mantissa[0] < 0) {
224 if (exponent[1] > 0x7) {
227 if (exponent[1] < 0x0) {
228 factor = -exponent[1];
231 tmp1 = int(std::fabs(mantissa[1]) * 0x10) >> factor;
232 tmp1 |= (exponent[1] & 0x7) << 5;
233 if (mantissa[1] < 0) {
239 if (exponent[2] > 0x7) {
242 if (exponent[2] < 0x0) {
243 factor = -exponent[2];
246 tmp2 = int(std::fabs(mantissa[2]) * 0x10) >> factor;
247 tmp2 |= (exponent[2] & 0x7) << 5;
248 if (mantissa[2] < 0) {
254 if (exponent[3] > 0x7) {
257 if (exponent[3] < 0x0) {
258 factor = -exponent[3];
261 tmp3 = int(std::fabs(mantissa[3]) * 0x10) >> factor;
262 tmp3 |= (exponent[3] & 0x7) << 5;
263 if (mantissa[3] < 0) {
278 if (exponent[0] > 0xF) {
281 if (exponent[0] < 0x0) {
282 factor = -exponent[0];
285 tmp0 = int(std::fabs(mantissa[0]) * 0x80) >> factor;
286 if (mantissa[0] < 0) {
292 if (exponent[1] > 0xF) {
295 if (exponent[1] < 0x0) {
296 factor = -exponent[1];
299 tmp1 = int(std::fabs(mantissa[1]) * 0x80) >> factor;
300 if (mantissa[1] < 0) {
306 if (exponent[2] > 0xF) {
309 if (exponent[2] < 0x0) {
310 factor = -exponent[2];
313 tmp2 = int(std::fabs(mantissa[2]) * 0x80) >> factor;
314 if (mantissa[2] < 0) {
320 if (exponent[3] > 0xF) {
323 if (exponent[3] < 0x0) {
324 factor = -exponent[3];
327 tmp3 = int(std::fabs(mantissa[3]) * 0x80) >> factor;
328 if (mantissa[3] < 0) {
332 vec0 = (tmp0 & 0xFF);
333 vec0 += (tmp1 & 0xFF) << 8;
334 vec1 = (tmp2 & 0xFF);
335 vec1 += (tmp3 & 0xFF) << 8;
337 vec2 |= exponent[1] << 4;
338 vec2 |= exponent[2] << 8;
339 vec2 |= exponent[3] << 12;
348 if (exponent[0] > 31) {
351 if (exponent[0] < 0) {
352 factor = -exponent[0];
355 tmp0 = int(std::fabs(mantissa[0]) * 0x400) >> factor;
356 tmp0 |= exponent[0] << 11;
357 if (mantissa[0] < 0) {
363 if (exponent[1] > 31) {
366 if (exponent[1] < 0) {
367 factor = -exponent[1];
370 tmp1 = int(std::fabs(mantissa[1]) * 0x400) >> factor;
371 tmp1 |= exponent[1] << 11;
372 if (mantissa[1] < 0) {
378 if (exponent[2] > 31) {
381 if (exponent[2] < 0) {
382 factor = -exponent[2];
385 tmp2 = int(std::fabs(mantissa[2]) * 0x400) >> factor;
386 tmp2 |= exponent[2] << 11;
387 if (mantissa[2] < 0) {
393 if (exponent[3] > 31) {
396 if (exponent[3] < 0) {
397 factor = -exponent[3];
400 tmp3 = int(std::fabs(mantissa[3]) * 0x400) >> factor;
401 tmp3 |= exponent[3] << 11;
402 if (mantissa[3] < 0) {
428 ps.push_back(mass[0]);
429 ps.push_back(mass[1]);
431 msg << MSG::VERBOSE <<
" compress x : " << momRaw.
m_px
432 <<
" = " << px[0] <<
" = " << px[1] <<
endmsg;
433 msg << MSG::VERBOSE <<
" compress y : " << momRaw.
m_py
434 <<
" = " << py[0] <<
" = " << py[1] <<
endmsg;
435 msg << MSG::VERBOSE <<
" compress z : " << momRaw.
m_pz
436 <<
" = " << pz[0] <<
" = " << pz[1] <<
endmsg;
437 msg << MSG::VERBOSE <<
" compress m : " << momRaw.
m_m
438 <<
" = " << mass[0] <<
" = " << mass[1] <<
endmsg;
443 if (
msg.level() <= MSG::VERBOSE) {
444 msg << MSG::VERBOSE <<
" compress # ps : " << ps.size() <<
" : ";
445 for (JetConverterTypes::signalState_pers_t::const_iterator it = ps.begin(); it != ps.end(); ++it) {
446 msg << std::hex << *it <<
" ";
449 msg << MSG::VERBOSE <<
" compress x : " << momRaw.
m_px <<
endmsg;
450 msg << MSG::VERBOSE <<
" compress y : " << momRaw.
m_py <<
endmsg;
451 msg << MSG::VERBOSE <<
" compress z : " << momRaw.
m_pz <<
endmsg;
452 msg << MSG::VERBOSE <<
" compress m : " << momRaw.
m_m <<
endmsg;
475 msg << MSG::DEBUG <<
" compress M : " << delta0 <<
" = "
476 << mantissa[0] <<
" *2^ " << exponent[0] - d0 <<
endmsg;
477 msg << MSG::DEBUG <<
" compress phi : " << delta1 <<
" = "
478 << mantissa[1] <<
" *2^ " << exponent[1] - d1 <<
endmsg;
479 msg << MSG::DEBUG <<
" compress eta : " << delta2 <<
" = "
480 << mantissa[2] <<
" *2^ " << exponent[2] - d1 <<
endmsg;
481 msg << MSG::DEBUG <<
" compress pT : " << delta3 <<
" = "
482 << mantissa[3] <<
" *2^ " << exponent[3] - d1 <<
endmsg;
490 MsgStream&
msg)
const {
492 int exponent[4] = {0};
493 double mantissa[4] = {0};
497 exponent[0] = ((ps[0] >> 5) & 0x7) - 2;
498 exponent[1] = ((ps[0] >> 13) & 0x7) - 5;
499 exponent[2] = ((ps[1] >> 5) & 0x7) - 5;
500 exponent[3] = ((ps[1] >> 13) & 0x7) - 5;
502 mantissa[0] = double( ps[0] & 0x1F) / 0x10;
503 mantissa[1] = double((ps[0] >> 8) & 0x1F) / 0x10;
504 mantissa[2] = double( ps[1] & 0x1F) / 0x10;
505 mantissa[3] = double((ps[1] >> 8) & 0x1F) / 0x10;
507 mantissa[0] = -mantissa[0];
509 if (ps[0] & 0x1000) {
510 mantissa[1] = -mantissa[1];
513 mantissa[2] = -mantissa[2];
515 if (ps[1] & 0x1000) {
516 mantissa[3] = -mantissa[3];
521 exponent[0] = ( ps[2] & 0xF) - 3;
522 exponent[1] = ((ps[2] >> 4) & 0xF) - 6;
523 exponent[2] = ((ps[2] >> 8) & 0xF) - 6;
524 exponent[3] = ((ps[2] >> 12) & 0xF) - 6;
526 mantissa[0] = double( ps[0] & 0x7F) / 0x80;
527 mantissa[1] = double((ps[0] >> 8) & 0x7F) / 0x80;
528 mantissa[2] = double( ps[1] & 0x7F) / 0x80;
529 mantissa[3] = double((ps[1] >> 8) & 0x7F) / 0x80;
531 mantissa[0] = -mantissa[0];
533 if (ps[0] & 0x8000) {
534 mantissa[1] = -mantissa[1];
537 mantissa[2] = -mantissa[2];
539 if (ps[1] & 0x8000) {
540 mantissa[3] = -mantissa[3];
545 exponent[0] = ((ps[0] >> 11) & 0x1F) - 5;
546 exponent[1] = ((ps[1] >> 11) & 0x1F) - 7;
547 exponent[2] = ((ps[2] >> 11) & 0x1F) - 7;
548 exponent[3] = ((ps[3] >> 11) & 0x1F) - 7;
550 mantissa[0] = double(ps[0] & 0x3FF) / 0x400;
551 mantissa[1] = double(ps[1] & 0x3FF) / 0x400;
552 mantissa[2] = double(ps[2] & 0x3FF) / 0x400;
553 mantissa[3] = double(ps[3] & 0x3FF) / 0x400;
555 mantissa[0] = -mantissa[0];
558 mantissa[1] = -mantissa[1];
561 mantissa[2] = -mantissa[2];
564 mantissa[3] = -mantissa[3];
572 <<
"Bad signal state size " << ps.size()
573 <<
" in SignalStateCnv::decompress; returning default momentum."
583 msg << MSG::VERBOSE <<
" RS x : " << momRaw.
m_px
584 <<
" = " << ps[0] <<
" = " << ps[1] <<
endmsg;
585 msg << MSG::VERBOSE <<
" RS y : " << momRaw.
m_py
586 <<
" = " << ps[2] <<
" = " << ps[3] <<
endmsg;
587 msg << MSG::VERBOSE <<
" RS z : " << momRaw.
m_pz
588 <<
" = " << ps[4] <<
" = " << ps[5] <<
endmsg;
589 msg << MSG::VERBOSE <<
" RS m : " << momRaw.
m_m
590 <<
" = " << ps[6] <<
" = " << ps[7] <<
endmsg;
595 double delta0 = std::ldexp(mantissa[0], exponent[0]);
596 double delta1 = std::ldexp(mantissa[1], exponent[1]);
597 double delta2 = std::ldexp(mantissa[2], exponent[2]);
598 double delta3 = std::ldexp(mantissa[3], exponent[3]);
600 if ((delta0 == -1) || (delta3 == -1)) {
602 <<
"A jet was badly decompressed. Returning null to avoid fpe in SignalStateCnv::decompress "
607 msg << MSG::VERBOSE <<
" RS 0 : " << delta0 <<
" = " << mantissa[0] <<
" = " << exponent[0] <<
endmsg;
608 msg << MSG::VERBOSE <<
" RS 1 : " << delta1 <<
" = " << mantissa[1] <<
" = " << exponent[1] <<
endmsg;
609 msg << MSG::VERBOSE <<
" RS 2 : " << delta2 <<
" = " << mantissa[2] <<
" = " << exponent[2] <<
endmsg;
610 msg << MSG::VERBOSE <<
" RS 3 : " << delta3 <<
" = " << mantissa[3] <<
" = " << exponent[3] <<
endmsg;
612 double angleC = std::atan2(momCal.
m_py, momCal.
m_px);
613 double p_traC = std::sqrt(momCal.
m_px * momCal.
m_px + momCal.
m_py * momCal.
m_py);
615 double angleR = angleC - delta1;
616 double eta_R = momCal.
eta() - delta2;
617 double p_traR = p_traC / (delta3 + 1);
618 momRaw.
m_m = momCal.
m_m / (delta0 + 1);
619 momRaw.
m_px = p_traR * std::cos(angleR);
620 momRaw.
m_py = p_traR * std::sin(angleR);
621 momRaw.
m_pz = p_traR * std::sinh(eta_R);
624 msg << MSG::VERBOSE <<
" RS # ps : " << ps.size() <<
" : ";
625 for (JetConverterTypes::signalState_pers_t::const_iterator it = ps.begin(); it != ps.end(); ++it) {
626 msg << std::hex << *it <<
" ";
632 msg << MSG::VERBOSE <<
" RS m : " << momRaw.
m_m <<
endmsg;