ATLAS Offline Software
Loading...
Searching...
No Matches
JetSignalStateCnv.h
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3*/
4
5#ifndef JETEVENTTPCNV_SIGNALSTATESTORE_H
6#define JETEVENTTPCNV_SIGNALSTATESTORE_H
7
45
46#include "JetEventTPCnv/JetConverterBase.h"
47
48#include <cmath>
49
50class SignalStateCnv {
51 public:
52
53 SignalStateCnv(){
54 m_delta_p =1.0;
55 m_step_p = 1.0 / 127;
56 m_delta_m = 0.4;
57 m_step_m = 0.4 / 127;
58 m_center = 0.8;
59 }
60
61 // converting double ratio to char to store uncalibrated quantities :
62 // We will assume that ratio uncal/cal is in [1-delta_m, 1+delta_p]
63 double m_delta_p ;
64 double m_step_p ;
65 double m_delta_m ;
66 double m_step_m ;
67 double m_center ;
68
69
70 // round to closest integer :
71 char round(double d) const {return (char)floor(d+0.5);}
72
73 char char_from_ratio(double r) const {
74 double d= r -m_center ;
75 if( d>0){
76 if(d>m_delta_p) d = m_delta_p;
77 return round(d/m_step_p);
78 }else{
79 if(d<-m_delta_m) d = m_delta_m;
80 return (round(fabs(d)/m_step_m) | 128);
81 }
82
83 }
84
85 double ratio_from_char(char c) const {
86 bool isneg = c & 128; // are storing a negativ diff ?
87 int val = c & 127; // the actual diff value
88 if(isneg){
89 return (m_center - val*m_step_m);
90 }else{
91 return (m_center + val*m_step_p);
92 }
93
94 }
95
96 enum CompressionLevel { HIGH, MEDIUM, LOW, NO_COMPRESSION };
97
98 /* the new configurable converters */
99 JetConverterTypes::signalState_pers_t
100 compress( const JetConverterTypes::momentum & momCal,
101 const JetConverterTypes::momentum & momRaw,
102 CompressionLevel level, MsgStream& msg ) const
103 {
104 int exponent[4];
105 double mantissa[4] = {0};
106 JetConverterTypes::signalState_pers_t ps;
107
108 double delta0(0);
109 double delta1(0);
110 double delta2(0);
111 double delta3(0);
112 int factor;
113
114 msg << MSG::VERBOSE << "called JetSignalStateCnv::compress() for :" << endmsg;
115 msg << MSG::VERBOSE << " raw momentum ( px | py | pz | m ) : ( "
116 << momRaw.m_px << " | "
117 << momRaw.m_py << " | "
118 << momRaw.m_pz << " | "
119 << momRaw.m_m << " )" << endmsg;
120 msg << MSG::VERBOSE << " cal momentum ( px | py | pz | m ) : ( "
121 << momCal.m_px << " | "
122 << momCal.m_py << " | "
123 << momCal.m_pz << " | "
124 << momCal.m_m << " )" << endmsg;
125
126 bool forceNoCompression=false;
127
128 // switch off compression, if M or pT of raw or calibrated jet are zero...
129 if ( momRaw.m_m == 0 || momCal.m_m == 0 )
130 forceNoCompression=true;
131 if ( momRaw.m_px == 0 && momRaw.m_py == 0 )
132 forceNoCompression=true;
133 if ( momCal.m_px == 0 && momCal.m_py == 0 )
134 forceNoCompression=true;
135
136 if ( forceNoCompression )
137 {
138 msg << MSG::DEBUG << "M or PT of calibrated or raw signal state of jet zero !"
139 << " Switching off compression for this jet !!!" << endmsg;
140 level = NO_COMPRESSION;
141 }
142
143 // prepare for compression the ratios
144 if ( level != NO_COMPRESSION )
145 {
146 double angleC = atan2( momCal.m_py, momCal.m_px );
147 double angleR = atan2( momRaw.m_py, momRaw.m_px );
148 double p_traC = sqrt( momCal.m_px * momCal.m_px + momCal.m_py * momCal.m_py );
149 double p_traR = sqrt( momRaw.m_px * momRaw.m_px + momRaw.m_py * momRaw.m_py );
150 mantissa[0] = frexp( momCal.m_m / momRaw.m_m - 1, &exponent[0] );
151 mantissa[1] = frexp( angleC - angleR, &exponent[1] );
152 mantissa[2] = frexp( momCal.eta() - momRaw.eta(), &exponent[2] );
153 mantissa[3] = frexp( p_traC / p_traR - 1, &exponent[3] );
154
155 delta0 = momCal.m_m / momRaw.m_m - 1;
156 delta1 = angleC - angleR;
157 delta2 = momCal.eta() - momRaw.eta();
158 delta3 = p_traC / p_traR - 1;
159 }
160
161 unsigned long tmp0(0);
162 unsigned long tmp1(0);
163 unsigned long tmp2(0);
164 unsigned long tmp3(0);
165
166 // here, we'll store the
167 unsigned short vec0(0);
168 unsigned short vec1(0);
169 unsigned short vec2(0);
170
171 switch(level)
172 {
173 case HIGH: // HIGH compression means, store differnces in 2 shorts
174 factor = 0; // compress the mass ratio
175 exponent[0] += 2;
176 if ( exponent[0] > 0x7 ) exponent[0] = 0x7;
177 if ( exponent[0] < 0x0 )
178 {
179 factor = -exponent[0];
180 exponent[0] = 0x0;
181 }
182 tmp0 = int( fabs(mantissa[0]) * 0x10 ) >> factor;
183 tmp0 |= ( exponent[0] & 0x7 ) << 5;
184 if ( mantissa[0] < 0 ) tmp0 |= 0x10;
185
186 factor = 0; // compress delta phi
187 exponent[1] += 5;
188 if ( exponent[1] > 0x7 ) exponent[1] = 0x7;
189 if ( exponent[1] < 0x0 )
190 {
191 factor = -exponent[1];
192 exponent[1] = 0x0;
193 }
194 tmp1 = int( fabs(mantissa[1]) * 0x10 ) >> factor;
195 tmp1 |= ( exponent[1] & 0x7 ) << 5;
196 if ( mantissa[1] < 0 ) tmp1 |= 0x10;
197
198 factor = 0; // compresss delta eta
199 exponent[2] += 5;
200 if ( exponent[2] > 0x7 ) exponent[2] = 0x7;
201 if ( exponent[2] < 0x0 )
202 {
203 factor = -exponent[2];
204 exponent[2] = 0x0;
205 }
206 tmp2 = int( fabs(mantissa[2]) * 0x10 ) >> factor;
207 tmp2 |= ( exponent[2] & 0x7 ) << 5;
208 if ( mantissa[2] < 0 ) tmp2 |= 0x10;
209
210 factor = 0; // compresss the pT ratio
211 exponent[3] += 5;
212 if ( exponent[3] > 0x7 ) exponent[3] = 0x7;
213 if ( exponent[3] < 0x0 )
214 {
215 factor = -exponent[3];
216 exponent[3] = 0x0;
217 }
218 tmp3 = int( fabs(mantissa[3]) * 0x10 ) >> factor;
219 tmp3 |= ( exponent[3] & 0x7 ) << 5;
220 if ( mantissa[3] < 0 ) tmp3 |= 0x10;
221
222 vec0 = tmp0;
223 vec0 += tmp1 << 8;
224 vec1 = tmp2;
225 vec1 += tmp3 << 8;
226 ps.push_back(vec0);
227 ps.push_back(vec1);
228 break;
229
230 case MEDIUM: // MEDIUM compression means, store differences in 3 shorts
231 factor = 0; // compress the mass ratio
232 exponent[0] += 3;
233 if ( exponent[0] > 0xF ) exponent[0] = 0xF;
234 if ( exponent[0] < 0x0 )
235 {
236 factor = -exponent[0];
237 exponent[0] = 0;
238 }
239 tmp0 = int( fabs(mantissa[0]) * 0x80 ) >> factor;
240 if ( mantissa[0] < 0 ) tmp0 |= 0x80;
241
242 factor = 0; // compress delta phi
243 exponent[1] += 6;
244 if ( exponent[1] > 0xF ) exponent[1] = 0xF;
245 if ( exponent[1] < 0x0 )
246 {
247 factor = -exponent[1];
248 exponent[1] = 0;
249 }
250 tmp1 = int( fabs(mantissa[1]) * 0x80 ) >> factor;
251 if ( mantissa[1] < 0 ) tmp1 |= 0x80;
252
253 factor = 0; // compresss delta eta
254 exponent[2] += 6;
255 if ( exponent[2] > 0xF ) exponent[2] = 0xF;
256 if ( exponent[2] < 0x0 )
257 {
258 factor = -exponent[2];
259 exponent[2] = 0;
260 }
261 tmp2 = int( fabs(mantissa[2]) * 0x80 ) >> factor;
262 if ( mantissa[2] < 0 ) tmp2 |= 0x80;
263
264 factor = 0; // compresss the pT ratio
265 exponent[3] += 6;
266 if ( exponent[3] > 0xF ) exponent[3] = 0xF;
267 if ( exponent[3] < 0x0 )
268 {
269 factor = -exponent[3];
270 exponent[3] = 0;
271 }
272 tmp3 = int( fabs(mantissa[3]) * 0x80 ) >> factor;
273 if ( mantissa[3] < 0 ) tmp3 |= 0x80;
274
275 vec0 = ( tmp0 & 0xFF );
276 vec0 += ( tmp1 & 0xFF ) << 8;
277 vec1 = ( tmp2 & 0xFF );
278 vec1 += ( tmp3 & 0xFF ) << 8;
279 vec2 = exponent[0];
280 vec2 |= exponent[1] << 4;
281 vec2 |= exponent[2] << 8;
282 vec2 |= exponent[3] << 12;
283 ps.push_back(vec0);
284 ps.push_back(vec1);
285 ps.push_back(vec2);
286 break;
287
288 case LOW: // LOW compression means, store differences in 4 shorts
289 factor = 0; // compress the mass ratio
290 exponent[0] += 5;
291 if ( exponent[0] > 31 ) exponent[0] = 31;
292 if ( exponent[0] < 0 )
293 {
294 factor = -exponent[0];
295 exponent[0] = 0;
296 }
297 tmp0 = int( fabs(mantissa[0]) * 0x400 ) >> factor;
298 tmp0 |= exponent[0] << 11;
299 if ( mantissa[0] < 0 ) tmp0 |= 0x400;
300
301 factor = 0; // compress delta phi
302 exponent[1] += 7;
303 if ( exponent[1] > 31 ) exponent[1] = 31;
304 if ( exponent[1] < 0 )
305 {
306 factor = -exponent[1];
307 exponent[1] = 0;
308 }
309 tmp1 = int( fabs(mantissa[1]) * 0x400 ) >> factor;
310 tmp1 |= exponent[1] << 11;
311 if ( mantissa[1] < 0 ) tmp1 |= 0x400;
312
313 factor = 0; // compresss delta eta
314 exponent[2] += 7;
315 if ( exponent[2] > 31 ) exponent[2] = 31;
316 if ( exponent[2] < 0 )
317 {
318 factor = -exponent[2];
319 exponent[2] = 0;
320 }
321 tmp2 = int( fabs(mantissa[2]) * 0x400 ) >> factor;
322 tmp2 |= exponent[2] << 11;
323 if ( mantissa[2] < 0 ) tmp2 |= 0x400;
324
325 factor = 0; // compresss the pT ratio
326 exponent[3] += 7;
327 if ( exponent[3] > 31 ) exponent[3] = 31;
328 if ( exponent[3] < 0 )
329 {
330 factor = -exponent[3];
331 exponent[3] = 0;
332 }
333 tmp3 = int( fabs(mantissa[3]) * 0x400 ) >> factor;
334 tmp3 |= exponent[3] << 11;
335 if ( mantissa[3] < 0 ) tmp3 |= 0x400;
336
337 ps.push_back(tmp0);
338 ps.push_back(tmp1);
339 ps.push_back(tmp2);
340 ps.push_back(tmp3);
341 break;
342 case NO_COMPRESSION:
343 default:
344 union {
345 unsigned short s[2];
346 float f;
347 } m;
348
349 m.f = momRaw.m_px; ps.push_back(m.s[0]); ps.push_back(m.s[1]);
350 m.f = momRaw.m_py; ps.push_back(m.s[0]); ps.push_back(m.s[1]);
351 m.f = momRaw.m_pz; ps.push_back(m.s[0]); ps.push_back(m.s[1]);
352 m.f = momRaw.m_m; ps.push_back(m.s[0]); ps.push_back(m.s[1]);
353 msg << MSG::VERBOSE << " compress x : " << m.f << " = " << m.s[0] << " = " << m.s[1] << endmsg;
354 msg << MSG::VERBOSE << " compress y : " << m.f << " = " << m.s[0] << " = " << m.s[1] << endmsg;
355 msg << MSG::VERBOSE << " compress z : " << m.f << " = " << m.s[0] << " = " << m.s[1] << endmsg;
356 msg << MSG::VERBOSE << " compress m : " << m.f << " = " << m.s[0] << " = " << m.s[1] << endmsg;
357 break;
358 }
359 if ( msg.level() <= MSG::VERBOSE )
360 {
361 msg << MSG::VERBOSE << " compress # ps : " << ps.size() << " : ";
362 for ( JetConverterTypes::signalState_pers_t::const_iterator it=ps.begin(); it != ps.end(); ++it )
363 msg << std::hex << *it << " ";
364 msg << std::dec << endmsg;
365 msg << MSG::VERBOSE << " compress x : " << momRaw.m_px << endmsg;
366 msg << MSG::VERBOSE << " compress y : " << momRaw.m_py << endmsg;
367 msg << MSG::VERBOSE << " compress z : " << momRaw.m_pz << endmsg;
368 msg << MSG::VERBOSE << " compress m : " << momRaw.m_m << endmsg;
369
370 int d0(0);
371 int d1(0);
372
373 switch(level)
374 {
375 case HIGH: d0=2; d1=5; break;
376 case MEDIUM: d0=3; d1=6; break;
377 case LOW: d0=5; d1=7; break;
378 case NO_COMPRESSION:
379 default: break;
380 };
381
382 msg << MSG::DEBUG << " compress M : " << delta0 << " = "
383 << mantissa[0] << " *2^ " << exponent[0]-d0 << endmsg;
384 msg << MSG::DEBUG << " compress phi : " << delta1 << " = "
385 << mantissa[1] << " *2^ " << exponent[1]-d1 << endmsg;
386 msg << MSG::DEBUG << " compress eta : " << delta2 << " = "
387 << mantissa[2] << " *2^ " << exponent[2]-d1 << endmsg;
388 msg << MSG::DEBUG << " compress pT : " << delta3 << " = "
389 << mantissa[3] << " *2^ " << exponent[3]-d1 << endmsg;
390 }
391 return ps;
392 };
393
394 JetConverterTypes::momentum decompress( const JetConverterTypes::signalState_pers_t& ps,
395 JetConverterTypes::momentum momCal,
396 MsgStream& msg ) const
397 {
398 JetConverterTypes::momentum momRaw;
399 int exponent[4] = {0};
400 double mantissa[4] = {0};
401 //std::cout<< " decompressing momCal px="<< momCal.m_px << " "<< momCal.m_py << " "<< momCal.m_pz << " "<< momCal.m_m << std::endl;
402
403 switch( ps.size() )
404 {
405 case 2:
406 exponent[0] = ( ( ps[0] >> 5 ) & 0x7 ) - 2;
407 exponent[1] = ( ( ps[0] >> 13 ) & 0x7 ) - 5;
408 exponent[2] = ( ( ps[1] >> 5 ) & 0x7 ) - 5;
409 exponent[3] = ( ( ps[1] >> 13 ) & 0x7 ) - 5;
410
411 mantissa[0] = double( ps[0] & 0x1F ) / 0x10;
412 mantissa[1] = double( ( ps[0] >> 8 ) & 0x1F ) / 0x10;
413 mantissa[2] = double( ps[1] & 0x1F ) / 0x10;
414 mantissa[3] = double( ( ps[1] >> 8 ) & 0x1F ) / 0x10;
415 if ( ps[0] & 0x10 ) mantissa[0] = -mantissa[0];
416 if ( ps[0] & 0x1000 ) mantissa[1] = -mantissa[1];
417 if ( ps[1] & 0x10 ) mantissa[2] = -mantissa[2];
418 if ( ps[1] & 0x1000 ) mantissa[3] = -mantissa[3];
419 break;
420
421 case 3:
422 exponent[0] = ( ps[2] & 0xF ) - 3;
423 exponent[1] = ( ( ps[2] >> 4 ) & 0xF ) - 6;
424 exponent[2] = ( ( ps[2] >> 8 ) & 0xF ) - 6;
425 exponent[3] = ( ( ps[2] >> 12 ) & 0xF ) - 6;
426
427 mantissa[0] = double( ps[0] & 0x7F ) / 0x80;
428 mantissa[1] = double( ( ps[0] >> 8 ) & 0x7F ) / 0x80;
429 mantissa[2] = double( ps[1] & 0x7F ) / 0x80;
430 mantissa[3] = double( ( ps[1] >> 8 ) & 0x7F ) / 0x80;
431 if ( ps[0] & 0x80 ) mantissa[0] = -mantissa[0];
432 if ( ps[0] & 0x8000 ) mantissa[1] = -mantissa[1];
433 if ( ps[1] & 0x80 ) mantissa[2] = -mantissa[2];
434 if ( ps[1] & 0x8000 ) mantissa[3] = -mantissa[3];
435 break;
436
437 case 4:
438 exponent[0] = ( ( ps[0] >> 11 ) & 0x1F ) - 5;
439 exponent[1] = ( ( ps[1] >> 11 ) & 0x1F ) - 7;
440 exponent[2] = ( ( ps[2] >> 11 ) & 0x1F ) - 7;
441 exponent[3] = ( ( ps[3] >> 11 ) & 0x1F ) - 7;
442
443 mantissa[0] = double( ps[0] & 0x3FF ) / 0x400;
444 mantissa[1] = double( ps[1] & 0x3FF ) / 0x400;
445 mantissa[2] = double( ps[2] & 0x3FF ) / 0x400;
446 mantissa[3] = double( ps[3] & 0x3FF ) / 0x400;
447 if ( ps[0] & 0x400 ) mantissa[0] = -mantissa[0];
448 if ( ps[1] & 0x400 ) mantissa[1] = -mantissa[1];
449 if ( ps[2] & 0x400 ) mantissa[2] = -mantissa[2];
450 if ( ps[3] & 0x400 ) mantissa[3] = -mantissa[3];
451 break;
452
453 case 8:
454 default:
455 union {
456 unsigned short s[2];
457 float f;
458 } m;
459
460 m.s[0] = ps[0]; m.s[1] = ps[1]; momRaw.m_px = m.f;
461 m.s[0] = ps[2]; m.s[1] = ps[3]; momRaw.m_py = m.f;
462 m.s[0] = ps[4]; m.s[1] = ps[5]; momRaw.m_pz = m.f;
463 m.s[0] = ps[6]; m.s[1] = ps[7]; momRaw.m_m = m.f;
464 msg << MSG::VERBOSE << " RS x : " << m.f << " = " << m.s[0] << " = " << m.s[1] << endmsg;
465 msg << MSG::VERBOSE << " RS y : " << m.f << " = " << m.s[0] << " = " << m.s[1] << endmsg;
466 msg << MSG::VERBOSE << " RS z : " << m.f << " = " << m.s[0] << " = " << m.s[1] << endmsg;
467 msg << MSG::VERBOSE << " RS m : " << m.f << " = " << m.s[0] << " = " << m.s[1] << endmsg;
468 break;
469 };
470
471 if ( ps.size() < 8 )
472 {
473 double delta0 = ldexp( mantissa[0], exponent[0] );
474 double delta1 = ldexp( mantissa[1], exponent[1] );
475 double delta2 = ldexp( mantissa[2], exponent[2] );
476 double delta3 = ldexp( mantissa[3], exponent[3] );
477
478 //std::cout<< " decompressing "<< delta0 <<" "<< delta1 << " "<< delta2 <<" "<< delta3 << " "<< std::endl;
479
480 if( ( delta0==-1) || (delta3==-1)){
481 msg << MSG::WARNING << "A jet was badly decompressed. Returning null to avoid fpe in SignalStateCnv::decompress "<<endmsg;
482 return momRaw;
483 }
484
485 msg << MSG::VERBOSE << " RS 0 : " << delta0 << " = " << mantissa[0] << " = " << exponent[0] << endmsg;
486 msg << MSG::VERBOSE << " RS 1 : " << delta1 << " = " << mantissa[1] << " = " << exponent[1] << endmsg;
487 msg << MSG::VERBOSE << " RS 2 : " << delta2 << " = " << mantissa[2] << " = " << exponent[2] << endmsg;
488 msg << MSG::VERBOSE << " RS 3 : " << delta3 << " = " << mantissa[3] << " = " << exponent[3] << endmsg;
489
490 double angleC = atan2( momCal.m_py, momCal.m_px );
491 double p_traC = sqrt( momCal.m_px * momCal.m_px + momCal.m_py * momCal.m_py );
492
493 double angleR = angleC - delta1;
494 double eta_R = momCal.eta() - delta2;
495 double p_traR = p_traC / ( delta3 + 1 );
496 momRaw.m_m = momCal.m_m / ( delta0 + 1 );
497 momRaw.m_px = p_traR * cos(angleR);
498 momRaw.m_py = p_traR * sin(angleR);
499 momRaw.m_pz = p_traR * sinh(eta_R);
500 }
501 msg << MSG::VERBOSE << " RS # ps : " << ps.size() << " : ";
502 for ( JetConverterTypes::signalState_pers_t::const_iterator it=ps.begin(); it != ps.end(); ++it )
503 msg << std::hex << *it << " ";
504 msg << std::dec << endmsg;
505 msg << MSG::VERBOSE << " RS x : " << momRaw.m_px << endmsg;
506 msg << MSG::VERBOSE << " RS y : " << momRaw.m_py << endmsg;
507 msg << MSG::VERBOSE << " RS z : " << momRaw.m_pz << endmsg;
508 msg << MSG::VERBOSE << " RS m : " << momRaw.m_m << endmsg;
509
510 return momRaw;
511 };
512
513};
514#endif
endmsg
#define endmsg
Definition AnalysisConfig_Ntuple.cxx:63
SignalStateCnv::decompress
JetConverterTypes::momentum decompress(const JetConverterTypes::signalState_pers_t &ps, JetConverterTypes::momentum momCal, MsgStream &msg) const
Definition JetSignalStateCnv.h:394
SignalStateCnv::compress
JetConverterTypes::signalState_pers_t compress(const JetConverterTypes::momentum &momCal, const JetConverterTypes::momentum &momRaw, CompressionLevel level, MsgStream &msg) const
Definition JetSignalStateCnv.h:100
SignalStateCnv::ratio_from_char
double ratio_from_char(char c) const
Definition JetSignalStateCnv.h:85
SignalStateCnv::char_from_ratio
char char_from_ratio(double r) const
Definition JetSignalStateCnv.h:73
SignalStateCnv::round
char round(double d) const
Definition JetSignalStateCnv.h:71
r
int r
Definition globals.cxx:22
JetConverterTypes::signalState_pers_t
std::vector< unsigned short > signalState_pers_t
Definition JetConverterBase.h:12
msg
MsgStream & msg
Definition testRead.cxx:32
Generated on for ATLAS Offline Software by doxygen 1.14.0