ATLAS Offline Software
Loading...
Searching...
No Matches
TileCorrelation.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3*/
4
6//
7// Filename : TileCorrelation.cxx
8//
9// Author : Valencia TileCal group, cristobal.cuenca@cern.ch,
10// Mantained by Ximo Poveda, jpoveda@cern.ch
11//
12// Created : May, 2004
13// Moved to TileRecUtils on Jan'05
14//
16
19
24
25
26#include "CLHEP/Matrix/Matrix.h"
27//#include "TileConditions/TilePulseShapes.h"
28
29#include <cmath>
30#include <iostream>
31#include <iomanip>
32
33//using CLHEP::HepMatrix;
34
35
38 : AthMessaging ("TileCorrelation")
39 , m_SS()
40 , m_S()
41 , m_R()
42 , m_corr()
43 , m_corrSum()
44 , m_corrSum2()
45 , m_nCorr(0.0)
46 , m_N()
47 , m_jEntry(0)
48 , m_lag(0)
49 , m_nPairs()
50 , m_nD(0.0)
51 , m_S1()
52 , m_S2()
53 , m_S11()
54 , m_S12()
55 , m_S22()
56{
57}
58
59
63
66
67 ATH_MSG_DEBUG("SetCorrelationZero");
68
69 for (int ros = 0; ros < 4; ++ros)
70 for (unsigned int drawer = 0; drawer < TileCalibUtils::MAX_DRAWER ; ++drawer)
71 for (unsigned int channel = 0; channel < TileCalibUtils::MAX_CHAN; ++channel)
72 for (unsigned int gain = 0; gain < TileCalibUtils::MAX_GAIN; ++gain) {
73 m_N[ros][drawer][channel][gain] = 0;
74 //N[ros][drawer][channel][gain]=0;
75 for (int i = 0; i < dignum; ++i) {
76 m_S[ros][drawer][channel][gain][i] = 0.;
77 for (int j = 0; j < dignum; ++j) {
78 m_SS[ros][drawer][channel][gain][i][j] = 0.;
79 m_R[ros][drawer][channel][gain][i][j] = 0.;
80 }
81 }
82
83 for (m_lag = 0; m_lag < 9; ++m_lag) {
84 m_S1[ros][drawer][channel][gain][m_lag] = 0.;
85 m_S2[ros][drawer][channel][gain][m_lag] = 0.;
86 m_S11[ros][drawer][channel][gain][m_lag] = 0.;
87 m_S12[ros][drawer][channel][gain][m_lag] = 0.;
88 m_S22[ros][drawer][channel][gain][m_lag] = 0.;
89 m_nPairs[ros][drawer][channel][gain][m_lag] = 0;
90 m_corrSum[ros][drawer][channel][gain][m_lag] = 0.;
91 m_corrSum2[ros][drawer][channel][gain][m_lag] = 0.;
92 }
93 }
94
95 for (m_lag = 0; m_lag < 9; ++m_lag) {
96 m_corr[m_lag] = 0.;
97 }
98
99}
100
103
104 ATH_MSG_DEBUG("SetCorrelationDelta");
105
106 for (int ros = 0; ros < 4; ++ros)
107 for (unsigned int drawer = 0; drawer < TileCalibUtils::MAX_DRAWER ; ++drawer)
108 for (unsigned int channel = 0; channel < TileCalibUtils::MAX_CHAN; ++channel)
109 for (unsigned int gain = 0; gain < TileCalibUtils::MAX_GAIN; ++gain) {
110 m_N[ros][drawer][channel][gain] = 1;
111 for (int i = 0; i < dignum; ++i)
112 for (int j = 0; j < dignum; ++j)
113 if (i == j)
114 m_R[ros][drawer][channel][gain][i][j] = 1.;
115 else
116 m_R[ros][drawer][channel][gain][i][j] = 0.;
117 }
118}
119
121void TileCorrelation::sum(std::vector<double> &digits, int ros, int drawer, int channel, int gain, int &dignum) {
122
123 ATH_MSG_VERBOSE("Sum");
124
125 double N_d = 0.;
126 dignum = digits.size();
127
128 m_N[ros][drawer][channel][gain]++;
129 N_d = double(m_N[ros][drawer][channel][gain]);
130 for (int i = 0; i < dignum; ++i) {
131 m_S[ros][drawer][channel][gain][i] += digits[i];
132 for (int j = 0; j < dignum; ++j)
133 m_SS[ros][drawer][channel][gain][i][j] += digits[i] * digits[j];
134 //for (int j=0;j<i+1;j++) SS[ros][drawer][channel][gain][i][j]+=digits[i]*digits[j];
135 }
136
137 ATH_MSG_DEBUG(" Sum, ros=" << ros
138 << " drawer=" << drawer
139 << " channel=" << channel
140 << " gain=" << gain
141 << " N=" << m_N[ros][drawer][channel][gain]
142 << " Sum[1]=" << m_S[ros][drawer][channel][gain][1]
143 << " Sum[2]=" << m_S[ros][drawer][channel][gain][2]
144 << " Sum[1][2]=" << m_SS[ros][drawer][channel][gain][1][2]
145 << " Sum[1][1]=" << m_SS[ros][drawer][channel][gain][1][1]
146 << " Sum[2][2]=" << m_SS[ros][drawer][channel][gain][2][2]
147 << " B[1][2]=" << m_SS[ros][drawer][channel][gain][1][2] / N_d
148 - m_S[ros][drawer][channel][gain][1] / N_d * m_S[ros][drawer][channel][gain][2] / N_d
149 << " Correlation[1][2]=" << (N_d * m_SS[ros][drawer][channel][gain][1][2]
150 - m_S[ros][drawer][channel][gain][1] * m_S[ros][drawer][channel][gain][2])
151 / sqrt((N_d * m_SS[ros][drawer][channel][gain][1][1]
152 - m_S[ros][drawer][channel][gain][1]
153 * m_S[ros][drawer][channel][gain][1])
154 * (N_d * m_SS[ros][drawer][channel][gain][2][2]
155 - m_S[ros][drawer][channel][gain][2]
156 * m_S[ros][drawer][channel][gain][2]))
157 );
158
159}
160
163
164 ATH_MSG_DEBUG("CalcCorrelation");
165
166 for (int ros = 0; ros < 4; ++ros)
167 for (unsigned int drawer = 0; drawer < TileCalibUtils::MAX_DRAWER ; ++drawer)
168 for (unsigned int channel = 0; channel < TileCalibUtils::MAX_CHAN; ++channel)
169 for (unsigned int gain = 0; gain < TileCalibUtils::MAX_GAIN; ++gain) {
170 double N_d = double(m_N[ros][drawer][channel][gain]);
171 if (N_d > 0.)
172 ATH_MSG_VERBOSE(" CalcCorrelation, ros=" << ros
173 << " drawer=" << drawer
174 << " channel=" << channel
175 << " gain=" << gain
176 << " N_d=" << N_d);
177
178 for (int i = 0; i < dignum; ++i)
179 // for (int j=0;j<i+1;j++)
180 for (int j = 0; j < dignum; ++j) {
181 if (N_d > 0.) {
182 // std::cout<<"b i="<<i<<" j="<<j<<std::endl
183 // <<" R[ros][drawer][channel][gain][i][j]="<<R[ros][drawer][channel][gain][i][j]<<std::endl
184 // <<" N_d="<<N_d<<std::endl
185 // <<" S[ros][drawer][channel][gain][i]="<<S[ros][drawer][channel][gain][i]<<std::endl
186 // <<" S[ros][drawer][channel][gain][j]="<<S[ros][drawer][channel][gain][j]<<std::endl
187 // <<" SS[ros][drawer][channel][gain][i][j]="<<SS[ros][drawer][channel][gain][i][j]<<std::endl
188 // <<" SS[ros][drawer][channel][gain][i][i]="<<SS[ros][drawer][channel][gain][i][i]<<std::endl
189 // <<" SS[ros][drawer][channel][gain][j][j]="<<SS[ros][drawer][channel][gain][j][j]<<std::endl
190 //<<" S[ros][drawer][channel][gain][i]*S[ros][drawer][channel][gain][j]="<<S[ros][drawer][channel][gain][i]*S[ros][drawer][channel][gain][j]<<std::endl
191 //<<" N_d*SS[ros][drawer][channel][gain][i][j]="<<N_d*SS[ros][drawer][channel][gain][i][j]<<std::endl
192 // <<" N_d*SS[ros][drawer][channel][gain][i][j]-S[ros][drawer][channel][gain][i]*S[ros][drawer][channel][gain][j]="<<N_d*SS[ros][drawer][channel][gain][i][j]-S[ros][drawer][channel][gain][i]*S[ros][drawer][channel][gain][j]<<std::endl
193 // <<" N_d*SS[ros][drawer][channel][gain][j][j]-S[ros][drawer][channel][gain][j]*S[ros][drawer][channel][gain][j])="<<N_d*SS[ros][drawer][channel][gain][j][j]-S[ros][drawer][channel][gain][j]*S[ros][drawer][channel][gain][j]<<std::endl
194 // <<" N_d*SS[ros][drawer][channel][gain][i][i]-S[ros][drawer][channel][gain][i]*S[ros][drawer][channel][gain][i]="<<N_d*SS[ros][drawer][channel][gain][i][i]-S[ros][drawer][channel][gain][i]*S[ros][drawer][channel][gain][i]<<std::endl
195 // <<" sqrt((N_d*SS[ros][drawer][channel][gain][i][i]-S[ros][drawer][channel][gain][i]*S[ros][drawer][channel][gain][i])*(N_d*SS[ros][drawer][channel][gain][j][j]-S[ros][drawer][channel][gain][j]*S[ros][drawer][channel][gain][j]))="<< sqrt((N_d*SS[ros][drawer][channel][gain][i][i]-S[ros][drawer][channel][gain][i]*S[ros][drawer][channel][gain][i])*(N_d*SS[ros][drawer][channel][gain][j][j]-S[ros][drawer][channel][gain][j]*S[ros][drawer][channel][gain][j]))
196 // <<std::endl;
197 m_R[ros][drawer][channel][gain][i][j] = (N_d * m_SS[ros][drawer][channel][gain][i][j]
198 - m_S[ros][drawer][channel][gain][i] * m_S[ros][drawer][channel][gain][j])
199 / sqrt(
200 (N_d * m_SS[ros][drawer][channel][gain][i][i]
201 - m_S[ros][drawer][channel][gain][i] * m_S[ros][drawer][channel][gain][i])
202 * (N_d * m_SS[ros][drawer][channel][gain][j][j]
203 - m_S[ros][drawer][channel][gain][j] * m_S[ros][drawer][channel][gain][j]));
204 // std::cout<<"R[ros][drawer][channel][gain][i][j]="<<R[ros][drawer][channel][gain][i][j]<<std::endl;
205 } else
206 m_R[ros][drawer][channel][gain][i][j] = -1234.;
207 }
208 }
209}
210
212void TileCorrelation::runningCorrelation(std::vector<double> &digits, int ros, int drawer, int channel, int gain,
213 int &dignum, int chthres) {
214
215 ATH_MSG_VERBOSE("RunningCorrelation");
216
217 dignum = digits.size();
218
219 //chthres=10;
220 //update sums
221 m_N[ros][drawer][channel][gain]++;
222 m_jEntry = m_N[ros][drawer][channel][gain];
223 m_nD = double(m_jEntry);
224
225 if (ros == 1 && drawer == 1 && channel == 0 && gain == 1)
226 ATH_MSG_INFO("Computing RunningCorrelation for jentry=" << m_jEntry);
227
228 for (m_lag = 1; m_lag < dignum; ++m_lag) {
229 for (int i = 0; i < dignum - m_lag; ++i) {
230 m_S1[ros][drawer][channel][gain][m_lag - 1] += digits[i];
231 m_S2[ros][drawer][channel][gain][m_lag - 1] += digits[i + m_lag];
232 m_S12[ros][drawer][channel][gain][m_lag - 1] += digits[i] * digits[i + m_lag];
233 m_S11[ros][drawer][channel][gain][m_lag - 1] += digits[i] * digits[i];
234 m_S22[ros][drawer][channel][gain][m_lag - 1] += digits[i + m_lag] * digits[i + m_lag];
235 m_nPairs[ros][drawer][channel][gain][m_lag - 1]++;
236 }
237 if (m_lag == 1 && ros == 1 && drawer == 1 && channel == 0 && gain == 1)
238 ATH_MSG_VERBOSE(" jentry=" << m_jEntry
239 << " N=" << m_nPairs[ros][drawer][channel][gain][m_lag - 1]
240 << " S1=" << m_S1[ros][drawer][channel][gain][m_lag - 1]);
241
242
243 if (m_jEntry > chthres) {
244 m_nCorr = double(m_nPairs[ros][drawer][channel][gain][m_lag - 1]);
245 m_corr[m_lag - 1] = (m_nCorr * m_S12[ros][drawer][channel][gain][m_lag - 1]
246 - m_S1[ros][drawer][channel][gain][m_lag - 1] * m_S2[ros][drawer][channel][gain][m_lag - 1])
247 / sqrt(
248 (m_nCorr * m_S11[ros][drawer][channel][gain][m_lag - 1]
249 - m_S1[ros][drawer][channel][gain][m_lag - 1] * m_S1[ros][drawer][channel][gain][m_lag - 1])
250 * (m_nCorr * m_S22[ros][drawer][channel][gain][m_lag - 1]
251 - m_S2[ros][drawer][channel][gain][m_lag - 1] * m_S2[ros][drawer][channel][gain][m_lag - 1]));
252
253 if (m_lag == 1 && ros == 1 && drawer == 1 && channel == 0 && gain == 1)
254 ATH_MSG_DEBUG(" corr=" << m_corr[m_lag - 1]
255 << " corr_sum=" << m_corrSum[ros][drawer][channel][gain][m_lag - 1]
256 << " corr_sum_sq=" << m_corrSum2[ros][drawer][channel][gain][m_lag - 1]);
257
258
259 m_corrSum[ros][drawer][channel][gain][m_lag - 1] += m_corr[m_lag - 1];
260 m_corrSum2[ros][drawer][channel][gain][m_lag - 1] += m_corr[m_lag - 1] * m_corr[m_lag - 1];
261 // corr_mean=corr_sum[lag-1]/(chthres-jentry);
262 // corr_RMS=sqrt(corr_sum_sq[lag-1]*(chthres-jentry)-corr_sum[lag-1]*corr_sum[lag-1])/(chthres-jentry)
263
264 if (m_lag == 1 && ros == 1 && drawer == 1 && channel == 0 && gain == 1)
265 ATH_MSG_DEBUG(" jentry=" << m_jEntry
266 << " jentry-chthres=" << m_jEntry - chthres
267 << " lag=1, ros=1, drawer=1, channel=0, gain=1"
268 << " corr=" << m_corr[m_lag - 1]
269// <<" corr_mean="<<corr_sum[lag-1]/(jentry-chthres)
270 << " sum corr_mean=" << m_corrSum[ros][drawer][channel][gain][m_lag - 1]
271 << " corr_mean=" << m_corrSum[ros][drawer][channel][gain][m_lag - 1] / (m_jEntry - chthres)
272// <<" RMS="<<sqrt(corr_sum_sq[lag-1]*(jentry-chthres)-corr_sum[lag-1]*corr_sum[lag-1])/(jentry-chthres)
273 << " sum RMS=" << m_corrSum2[ros][drawer][channel][gain][m_lag - 1]
274 << " RMS=" << sqrt( m_corrSum2[ros][drawer][channel][gain][m_lag - 1] * (m_jEntry - chthres)
275 - m_corrSum[ros][drawer][channel][gain][m_lag - 1]
276 * m_corrSum[ros][drawer][channel][gain][m_lag - 1]) / (m_jEntry - chthres));
277
278 }
279 }
280}
281
283void TileCorrelation::calculateRunningCorrelation(int dignum, int chthres, bool is7to9) {
284
285 ATH_MSG_VERBOSE("CalcRunningCorrelation");
286
287 for (int ros = 0; ros < 4; ++ros)
288 for (unsigned int drawer = 0; drawer < TileCalibUtils::MAX_DRAWER ; ++drawer)
289 for (unsigned int channel = 0; channel < TileCalibUtils::MAX_CHAN; ++channel)
290 for (unsigned int gain = 0; TileCalibUtils::MAX_GAIN < 2; ++gain) {
291 m_jEntry = m_N[ros][drawer][channel][gain];
292 m_nCorr = double(m_jEntry - chthres);
293
294 if (m_jEntry > 0) {
295 if (is7to9 && dignum == 7) {
296 for (int i = 0; i < 9; ++i)
297 m_R[ros][drawer][channel][gain][i][i] = 1.;
298
299 for (m_lag = 1; m_lag < 9; ++m_lag)
300 for (int i = 0; i < 9 - m_lag; ++i) {
301 if (m_lag < 7) {
302 m_R[ros][drawer][channel][gain][i][i + m_lag] = m_corrSum[ros][drawer][channel][gain][m_lag - 1]
303 / m_nCorr;
304 m_R[ros][drawer][channel][gain][i + m_lag][i] = m_corrSum[ros][drawer][channel][gain][m_lag - 1]
305 / m_nCorr;
306 } else {
307 m_R[ros][drawer][channel][gain][i][i + m_lag] = 0.;
308 m_R[ros][drawer][channel][gain][i + m_lag][i] = 0.;
309 }
310
311 if (-1. > m_R[ros][drawer][channel][gain][i][i + m_lag]
312 || m_R[ros][drawer][channel][gain][i][i + m_lag] > 1.)
313 m_R[ros][drawer][channel][gain][i][i + m_lag] = 0.;
314 if (-1. > m_R[ros][drawer][channel][gain][i + m_lag][i]
315 || m_R[ros][drawer][channel][gain][i + m_lag][i] > 1.)
316 m_R[ros][drawer][channel][gain][i + m_lag][i] = 0.;
317
318 }
319 } else {
320 for (int i = 0; i < dignum; i++)
321 m_R[ros][drawer][channel][gain][i][i] = 1.;
322
323 for (m_lag = 1; m_lag < dignum; ++m_lag)
324 for (int i = 0; i < dignum - m_lag; ++i) {
325 m_R[ros][drawer][channel][gain][i][i + m_lag] = m_corrSum[ros][drawer][channel][gain][m_lag - 1]
326 / m_nCorr;
327 m_R[ros][drawer][channel][gain][i + m_lag][i] = m_corrSum[ros][drawer][channel][gain][m_lag - 1]
328 / m_nCorr;
329 if (-1. > m_R[ros][drawer][channel][gain][i][i + m_lag]
330 || m_R[ros][drawer][channel][gain][i][i + m_lag] > 1.)
331 m_R[ros][drawer][channel][gain][i][i + m_lag] = 0.;
332 if (-1. > m_R[ros][drawer][channel][gain][i + m_lag][i]
333 || m_R[ros][drawer][channel][gain][i + m_lag][i] > 1.)
334 m_R[ros][drawer][channel][gain][i + m_lag][i] = 0.;
335 }
336 }
337 }
338 }
339}
340
343
344 std::cout << " TileCorrelation::PrintCorrelation()..." << std::endl;
345 for (int ros = 0; ros < 1; ++ros) {
346 std::cout << " ros=" << ros << std::endl;
347 for (int drawer = 31; drawer < 32; ++drawer) {
348 std::cout << " drawer=" << drawer << std::endl;
349 for (int channel = 17; channel < 24; ++channel) {
350 std::cout << " channel=" << channel << std::endl;
351 for (unsigned int gain = 0; gain < TileCalibUtils::MAX_GAIN; ++gain) {
352 std::cout << " gain=" << gain << std::endl;
353 for (int i = 0; i < dignum; ++i) {
354 for (int j = 0; j < dignum; ++j) {
355 std::cout << " " << m_R[ros][drawer][channel][gain][i][j];
356 }
357 std::cout << std::endl;
358 }
359 std::cout << std::endl;
360 }
361 }
362 }
363 }
364
365}
366
368void TileCorrelation::saveCorrelationSumm(bool deltaCorrelation, const std::string& correlationSummOptFilterFile,
369 const TileHWID *tileHWID, int dignum) {
370
371 ATH_MSG_DEBUG("SaveCorrelationSumm");
372
373 CLHEP::HepMatrix correlation(dignum, 1, 0);
374
375 std::fstream *correlationFile = new std::fstream(correlationSummOptFilterFile.c_str(), std::fstream::out);
376 if (correlationFile->is_open()) ATH_MSG_INFO(correlationSummOptFilterFile << " file open");
377
378 if (deltaCorrelation) {
379 // for (int i=0;i<dignum;i++)
380 for (int j = 0; j < dignum; ++j) {
381 int i = 0;
382 if (m_R[0][0][0][0][i][j] > -100000. && m_R[0][0][0][0][i][j] < 100000.)
383 correlation[i][j] = m_R[0][0][0][0][i][j];
384 else
385 correlation[i][j] = 0.0;
386 }
387
388 *correlationFile << correlation.T() << std::endl;
389 } else {
390 for (int ros = 0; ros < 4; ++ros)
391 for (unsigned int drawer = 0; drawer < TileCalibUtils::MAX_DRAWER ; ++drawer) {
392 int frag = tileHWID->frag(ros + 1, drawer);
393 for (unsigned int channel = 0; channel < TileCalibUtils::MAX_CHAN; ++channel)
394 for (unsigned int gain = 0; gain < TileCalibUtils::MAX_GAIN; ++gain) {
395 ATH_MSG_VERBOSE("ros " << ros
396 << " drawer " << drawer << MSG::hex
397 << " frag0x " << frag << MSG::dec
398 << " channel " << channel
399 << " gain " << gain
400 << " N " << m_N[ros][drawer][channel][gain]);
401
402 if (m_N[ros][drawer][channel][gain] > 0) {
403 //for (int i=0;i<dignum;i++)
404 for (int j = 0; j < dignum; ++j) {
405 int i = 0;
406 if (m_R[ros][drawer][channel][gain][i][j] > -100000. && m_R[ros][drawer][channel][gain][i][j] < 100000.)
407 correlation[i][j] = m_R[ros][drawer][channel][gain][i][j];
408 else
409 correlation[i][j] = 0.0;
410 }
411
412 *correlationFile << "ros " << ros
413 << " drawer " << drawer << std::hex
414 << " frag0x " << frag << std::dec
415 << " channel " << channel
416 << " gain " << gain
417 << " N " << m_N[ros][drawer][channel][gain]
418 << correlation.T();
419 // <<M_correlation.T()<<std::endl;
420 // for (int i=0;i<dignum;i++)
421 // for (int j=0;j<dignum;j++)
422 // *f_correlation<<R[ros][drawer][channel][gain][i][j]<< M_correlation[i][j]<<std::endl;
423
424 }
425 }
426 }
427 }
428 correlationFile->close();
429}
430
432void TileCorrelation::saveCorrelationMatrix(bool deltaCorrelation, const std::string& correlationMatrixOptFilterFile,
433 const TileHWID *tileHWID, int dignum) {
434
435 ATH_MSG_DEBUG("SaveCorrelationMatrix");
436
437 CLHEP::HepMatrix correlation(dignum, dignum, 0);
438
439 std::fstream *correlationFile = new std::fstream(correlationMatrixOptFilterFile.c_str(), std::fstream::out);
440 if (correlationFile->is_open()) ATH_MSG_INFO(correlationMatrixOptFilterFile << " file open");
441
442 if (deltaCorrelation) {
443 for (int i = 0; i < dignum; ++i)
444 for (int j = 0; j < dignum; ++j) {
445 if (m_R[0][0][0][0][i][j] > -100000. && m_R[0][0][0][0][i][j] < 100000.)
446 correlation[i][j] = m_R[0][0][0][0][i][j];
447 else
448 correlation[i][j] = 0.0;
449 }
450
451 *correlationFile << correlation << std::endl;
452 } else {
453 for (int ros = 0; ros < 4; ++ros)
454 for (unsigned int drawer = 0; drawer < TileCalibUtils::MAX_DRAWER ; ++drawer) {
455 int frag = tileHWID->frag(ros + 1, drawer);
456 for (unsigned int channel = 0; channel < TileCalibUtils::MAX_CHAN; ++channel)
457 for (unsigned int gain = 0; gain < TileCalibUtils::MAX_GAIN; ++gain) {
458 ATH_MSG_VERBOSE("ros " << ros
459 << " drawer " << drawer << MSG::hex
460 << " frag0x " << frag << MSG::dec
461 << " channel " << channel
462 << " gain " << gain
463 << " N " << m_N[ros][drawer][channel][gain]);
464
465
466 if (m_N[ros][drawer][channel][gain] > 0) {
467 for (int i = 0; i < dignum; ++i)
468 for (int j = 0; j < dignum; ++j) {
469 if (m_R[ros][drawer][channel][gain][i][j] > -100000.
470 && m_R[ros][drawer][channel][gain][i][j] < 100000.)
471 correlation[i][j] = m_R[ros][drawer][channel][gain][i][j];
472 else
473 correlation[i][j] = 0.0;
474 }
475
476 *correlationFile << "ros " << ros
477 << " drawer " << drawer << std::hex
478 << " frag0x " << frag << std::dec
479 << " channel " << channel
480 << " gain " << gain
481 << " N " << m_N[ros][drawer][channel][gain]
482 << correlation << std::endl;
483 // for (int i=0;i<dignum;i++)
484 // for (int j=0;j<dignum;j++)
485 // *f_correlation<<R[ros][drawer][channel][gain][i][j]<< M_correlation[i][j]<<std::endl;
486
487 }
488 }
489 }
490 }
491 correlationFile->close();
492}
493
496 (bool deltaCorrelation,
497 const std::vector<double>& shapeFormLG,
498 const std::vector<double>& shapeFormHG,
499 const std::vector<double>& shapeFormDerivativeLG,
500 const std::vector<double>& shapeFormDerivativeHG,
501 const std::string& aiLoOptFilterFile,
502 const std::string& biLoOptFilterFile,
503 const std::string& aiHiOptFilterFile,
504 const std::string& biHiOptFilterFile,
505 const TileHWID *tileHWID, int dignum) {
506
507 ATH_MSG_DEBUG("CalcWeights");
508
509 CLHEP::HepMatrix correlation(dignum, dignum, 0);
510 CLHEP::HepMatrix inverse(dignum, dignum, 0);
511 CLHEP::HepMatrix zero(dignum, dignum, 0);
512 CLHEP::HepMatrix pulseShape(dignum, 1, 0);
513 CLHEP::HepMatrix pulseShapeDerivative(dignum, 1, 0);
514 CLHEP::HepMatrix a(dignum, 1, 0);
515 CLHEP::HepMatrix b(dignum, 1, 0);
516
517 std::fstream *aiLoFile = new std::fstream(aiLoOptFilterFile.c_str(), std::fstream::out);
518 std::fstream *biLoFile = new std::fstream(biLoOptFilterFile.c_str(), std::fstream::out);
519 std::fstream *aiHiFile = new std::fstream(aiHiOptFilterFile.c_str(), std::fstream::out);
520 std::fstream *biHiFile = new std::fstream(biHiOptFilterFile.c_str(), std::fstream::out);
521
522 //Open Weights files
523 if (aiLoFile->is_open() && aiLoFile->is_open() && aiLoFile->is_open() && aiLoFile->is_open()) {
524 ATH_MSG_INFO(" Weights files open");
525 } else {
526 ATH_MSG_INFO("Weights files didn't open successfully");
527 }
528
529 //pulse shape
530 // std::vector<double> new_shapeForm;
531 // double max=0.;
532 // int nmax=0;
533
534 // for (int i=0; i<int(m_shapeForm.size());i++)
535 // if (m_shapeForm[i]>max)
536 // {
537 // max=m_shapeForm[i];
538 // nmax=i;
539 // }
540 // new_shapeForm.resize(dignum*25,0.);
541
542 // if (lDebug)
543 // log<<MSG::DEBUG<<"m_shapeForm.size()="<<m_shapeForm.size()<<"m_shapeForm nmax="<<nmax<<" new_shapeForm.size()="<<new_shapeForm.size()<<" new_shapeForm nmax="<<dignum*25/2<<endmsg;
544
545 // for (int i=0;i<dignum*25;i++)
546 // {
547 // if (i<(dignum*25/2-nmax)) new_shapeForm[i]=0.;
548 // if (i>=(dignum*25/2-nmax) && i<(dignum*25/2-nmax+int(m_shapeForm.size()))) new_shapeForm[i]=m_shapeForm[i-dignum*25/2 +nmax];
549 // if (i>=(dignum*25/2-nmax+int(m_shapeForm.size()))) new_shapeForm[i]=0.;
550 // }
551
552 if (msgLvl(MSG::VERBOSE)) {
553 msg(MSG::VERBOSE) << "shapeFormLG, shapeFormDerivativeLG, shapeFormHG, shapeFormDerivativeHG" << endmsg;
554 for (int i = 0; i < int(shapeFormLG.size()); ++i)
555 msg(MSG::VERBOSE) << i << " " << std::setw(18) << std::setprecision(10)
556 << shapeFormLG[i] << " " << std::setw(18) << std::setprecision(10)
557 << shapeFormDerivativeLG[i] << " " << std::setw(18) << std::setprecision(10)
558 << shapeFormHG[i] << " " << std::setw(18) << std::setprecision(10)
559 << shapeFormDerivativeHG[i] << " " << endmsg;
560 }
561 //if (lDebug) {
562 // log<<MSG::DEBUG<<"m_HshapeForm"<<endmsg;
563 // for(int i=0;i<int(m_HshapeForm.size());i++) log<<MSG::DEBUG<<" "<<i<<" "<<m_HshapeForm[i]<<endmsg;
564 // log<<MSG::DEBUG<<"m_LdshapeForm"<<endmsg;
565 // for(int i=0;i<int(m_LdshapeForm.size());i++) log<<MSG::DEBUG<<" "<<i<<" "<<m_LdshapeForm[i]<<endmsg;
566 // log<<MSG::DEBUG<<"shapeFormDerivativeHG"<<endmsg;
567 // for(int i=0;i<int(m_HdshapeForm.size());i++) log<<MSG::DEBUG<<" "<<i<<" "<<m_HdshapeForm[i]<<endmsg;
568 //}
569
570 double Q1, Q2, Q3, Delta;
571 int ierr = 0;
572
573 if (deltaCorrelation) {
574 for (unsigned int gain = 0; gain < TileCalibUtils::MAX_GAIN; ++gain)
575 for (int pha = -12; pha < 13; ++pha) {
576 correlation = zero;
577 inverse = zero;
578
579 for (int i = 0; i < dignum; ++i)
580 for (int j = 0; j < dignum; ++j)
581 correlation[i][j] = m_R[0][0][0][0][i][j];
582
583 inverse = correlation.inverse(ierr);
584 if (ierr == 0) {
585 for (int i = 0; i < dignum; ++i) {
586 if (gain == 0) {
587 pulseShape[i][0] = shapeFormLG[i * 25 + 12 + pha];
588 pulseShapeDerivative[i][0] = shapeFormDerivativeLG[i * 25 + 12 + pha];
589 } else {
590 pulseShape[i][0] = shapeFormHG[i * 25 + 12 + pha];
591 pulseShapeDerivative[i][0] = shapeFormDerivativeHG[i * 25 + 12 + pha];
592 }
593 }
594
595 Q1 = ((pulseShape.T()) * inverse * pulseShape).determinant();
596 Q2 = ((pulseShapeDerivative.T()) * inverse * pulseShapeDerivative).determinant();
597 Q3 = ((pulseShapeDerivative.T()) * inverse * pulseShape).determinant();
598 Delta = Q1 * Q2 - Q3 * Q3;
599
600 a = Q2 / Delta * inverse * pulseShape - Q3 / Delta * inverse * pulseShapeDerivative;
601 b = Q3 / Delta * inverse * pulseShape - Q1 / Delta * inverse * pulseShapeDerivative;
602
603 if (gain == 0) {
604 *aiLoFile << std::setw(6) << pha;
605 for (int i = 0; i < dignum; ++i)
606 *aiLoFile << std::setw(18) << std::setprecision(10) << a(i + 1, 1);
607 *aiLoFile << std::endl;
608
609 *biLoFile << std::setw(6) << pha;
610 for (int i = 0; i < dignum; ++i)
611 *biLoFile << std::setw(18) << std::setprecision(10) << b(i + 1, 1);
612 *biLoFile << std::endl;
613 } else {
614 *aiHiFile << std::setw(6) << pha;
615 for (int i = 0; i < dignum; ++i)
616 *aiHiFile << std::setw(18) << std::setprecision(10) << a(i + 1, 1);
617 *aiHiFile << std::endl;
618
619 *biHiFile << std::setw(6) << pha;
620 for (int i = 0; i < dignum; ++i)
621 *biHiFile << std::setw(18) << std::setprecision(10) << b(i + 1, 1);
622 *biHiFile << std::endl;
623 }
624 }
625 }
626 } else {
627 for (int ros = 0; ros < 4; ++ros)
628 for (unsigned int drawer = 0; drawer < TileCalibUtils::MAX_DRAWER ; ++drawer) {
629 int frag = tileHWID->frag(ros + 1, drawer);
630 for (unsigned int channel = 0; channel < TileCalibUtils::MAX_CHAN; ++channel)
631 for (unsigned int gain = 0; gain < TileCalibUtils::MAX_GAIN; ++gain)
632 if (m_N[ros][drawer][channel][gain] > 0) {
633 if (gain == 0) {
634 *aiLoFile << "ros " << ros
635 << " drawer " << drawer << std::hex
636 << " frag0x " << frag << std::dec
637 << " channel " << channel
638 << " N " << m_N[ros][drawer][channel][0] << std::endl;
639
640 *biLoFile << "ros " << ros
641 << " drawer " << drawer << std::hex
642 << " frag0x " << frag << std::dec
643 << " channel " << channel
644 << " N " << m_N[ros][drawer][channel][0] << std::endl;
645 }
646 if (gain == 1) {
647 *aiHiFile << "ros " << ros
648 << " drawer " << drawer << std::hex
649 << " frag0x " << frag << std::dec
650 << " channel " << channel
651 << " N " << m_N[ros][drawer][channel][1] << std::endl;
652
653 *biHiFile << "ros " << ros
654 << " drawer " << drawer << std::hex
655 << " frag0x " << frag << std::dec
656 << " channel " << channel
657 << " N " << m_N[ros][drawer][channel][1] << std::endl;
658 }
659
660 for (int pha = -12; pha < 13; ++pha) {
661 correlation = zero;
662 inverse = zero;
663
664 for (int i = 0; i < dignum; ++i)
665 for (int j = 0; j < dignum; ++j)
666 correlation[i][j] = m_R[ros][drawer][channel][gain][i][j];
667
668 inverse = correlation.inverse(ierr);
669 if (ierr == 0) {
670 for (int i = 0; i < dignum; ++i) {
671 if (gain == 0) {
672 pulseShape[i][0] = shapeFormLG[i * 25 + 12 + pha];
673 pulseShapeDerivative[i][0] = shapeFormDerivativeLG[i * 25 + 12 + pha];
674 } else {
675 pulseShape[i][0] = shapeFormHG[i * 25 + 12 + pha];
676 pulseShapeDerivative[i][0] = shapeFormDerivativeHG[i * 25 + 12 + pha];
677 }
678
679 // PulseShape[i][0]=new_shapeForm[i*25+12+pha];
680 // if ((i*25+12+pha)>0 || (i*25+12+pha)<224)
681 // DPulseShape[i][0]=.5*(new_shapeForm[i*25+13+pha]-new_shapeForm[i*25+11+pha]);
682 // else DPulseShape[i][0]=0.;
683 }
684
685 //HepStd::std::cout<<" correlation "<<Correlation<<Hepstd::endl;
686 //f_weights<<" correlation "<<Correlation<<std::endl;
687 // HepStd::std::cout<<" inverse Matrix "<<Correlation.inverse(ierr)<<Hepstd::endl;
688 //HepStd::std::cout<<" inverse Matrix "<<Inverse<<Hepstd::endl;
689 //f_weights<<" inverse Matrix "<<Inverse<<std::endl;
690 // HepStd::std::cout<<" Product "<<Inverse*Correlation<<Hepstd::endl;
691
692 //std::cout<<" Q1 number of columns="<<((PulseShape.T())*Inverse*PulseShape).num_col()
693 //<<" number of rows="<<((PulseShape.T())*Inverse*PulseShape).num_row()<<std::endl;
694
695 Q1 = ((pulseShape.T()) * inverse * pulseShape).determinant();
696 Q2 = ((pulseShapeDerivative.T()) * inverse * pulseShapeDerivative).determinant();
697 Q3 = ((pulseShapeDerivative.T()) * inverse * pulseShape).determinant();
698 Delta = Q1 * Q2 - Q3 * Q3;
699
700 //std::cout<<" Coeffs: Q1="<<Q1<<" Q2="<<Q2<<" Q3="<<Q3<<" Delta="<<Delta<<std::endl;
701 a = Q2 / Delta * inverse * pulseShape - Q3 / Delta * inverse * pulseShapeDerivative;
702 b = Q3 / Delta * inverse * pulseShape - Q1 / Delta * inverse * pulseShapeDerivative;
703
704 //HepStd::std::cout<<" a Weights= "<<a<<Hepstd::endl;
705 //HepStd::std::cout<<" b Weights= "<<b<<Hepstd::endl;
706
707 if (gain == 0) {
708 *aiLoFile << std::setw(6) << pha;
709 for (int i = 0; i < dignum; ++i)
710 *aiLoFile << std::setw(18) << std::setprecision(10) << a(i + 1, 1);
711 *aiLoFile << std::endl;
712
713 *biLoFile << std::setw(6) << pha;
714 for (int i = 0; i < dignum; ++i)
715 *biLoFile << std::setw(18) << std::setprecision(10) << b(i + 1, 1);
716 *biLoFile << std::endl;
717 } else {
718 *aiHiFile << std::setw(6) << pha;
719 for (int i = 0; i < dignum; ++i)
720 *aiHiFile << std::setw(18) << std::setprecision(10) << a(i + 1, 1);
721 *aiHiFile << std::endl;
722
723 *biHiFile << std::setw(6) << pha;
724 for (int i = 0; i < dignum; ++i)
725 *biHiFile << std::setw(18) << std::setprecision(10) << b(i + 1, 1);
726 *biHiFile << std::endl;
727 }
728 }
729 }
730 }
731 }
732 }
733
734 aiLoFile->close();
735 biLoFile->close();
736 aiHiFile->close();
737 biHiFile->close();
738}
739
741void TileCorrelation::buildPulseShape(std::vector<double> &pulseShape, std::vector<double> &pulseShapeY,
742 std::vector<double> &pulseShapeT, int dignum) {
743
744 ATH_MSG_DEBUG("BuildPulseShape");
745
746 //1: set m_pulseShape
747 pulseShape.resize(dignum * 25);
748 ATH_MSG_DEBUG("Set dimension of m_pulseShape to dignum*25=" << dignum * 25);
749
750 //2: scan m_pulseShapeT for: tmin, tmax, nt0 and size: m_pulseShapeX[nt0]=1.0;
751 int nt0 = 0, size;
752 double tmin = 10000., tmax = -10000.;
753 size = pulseShapeT.size();
754 for (int i = 0; i < size; ++i) {
755 if (pulseShapeT[i] < tmin) tmin = pulseShapeT[i];
756 if (pulseShapeT[i] > tmax) tmax = pulseShapeT[i];
757 if (pulseShapeT[i] == 0) nt0 = i;
758 }
759
760 ATH_MSG_DEBUG("pulseShapeY & pulseShapeT size =" << size
761 << ", tmin=" << tmin
762 << ", tmax=" << tmax
763 << " central point=" << nt0
764 << " pulseShapeT[nt0]=" << pulseShapeT[nt0]
765 << " pulseShapeY[nt0]=" << pulseShapeY[nt0]);
766
767
768 //3: fill m_pulseShape
769 bool exact;
770 int nminn, nminp;
771 double minn, minp, tdist;
772 pulseShape[dignum * 25 / 2] = pulseShapeY[nt0];
773 for (int i = 1; i < dignum * 25 / 2 + 1; ++i) {
774 // negative times: 0->dignum*25/2
775 if (-i < tmin)
776 pulseShape[dignum * 25 / 2 - i] = 0.;
777 else {
778 exact = false;
779 minn = -10000.;
780 minp = 10000.;
781 nminn = 0;
782 nminp = size - 1;
783 for (int j = 0; j < nt0 + 1 && !exact; ++j) {
784 if (pulseShapeT[j] == double(-i)) {
785 pulseShape[dignum * 25 / 2 - i] = pulseShapeY[j];
786 exact = true;
787 } else {
788 tdist = pulseShapeT[j] - double(-i);
789 if (tdist < 0. && tdist > minn) {
790 minn = tdist;
791 nminn = j;
792 }
793 if (tdist > 0. && tdist < minp) {
794 minp = tdist;
795 nminp = j;
796 }
797 }
798 }
799
800 if (exact) {
801 ATH_MSG_VERBOSE("exact value found for time=" << -i
802 << " pulseShape=" << pulseShape[dignum * 25 / 2 - i]);
803
804 } else {
805
806 ATH_MSG_VERBOSE("exact value NOT found for time=" << -i
807 << " nminn=" << nminn
808 << " pulseShapeT=" << pulseShapeT[nminn]
809 << " pulseShapeY=" << pulseShapeY[nminn] << std::endl
810 << " nminp=" << nminp
811 << " pulseShapeT=" << pulseShapeT[nminp]
812 << " pulseShapeY=" << pulseShapeY[nminp]);
813
814
815 pulseShape[dignum * 25 / 2 - i] = pulseShapeY[nminn]
816 + (pulseShapeY[nminp] - pulseShapeY[nminn]) / (pulseShapeT[nminp] - pulseShapeT[nminn])
817 * (-i - pulseShapeT[nminn]);
818 }
819
820 }
821
822 // positive times: dignum*25/2->dignum*25
823 if (i > tmax)
824 pulseShape[dignum * 25 / 2 + i] = 0.;
825 else {
826 exact = false;
827 minn = -10000.;
828 minp = 10000.;
829 nminn = 0;
830 nminp = size;
831 for (int j = nt0; j < size && !exact; ++j) {
832 if (pulseShapeT[j] == double(i)) {
833 pulseShape[dignum * 25 / 2 + i] = pulseShapeY[j];
834 exact = true;
835 } else {
836 tdist = pulseShapeT[j] - double(i);
837 if (tdist < 0) if (tdist > minn) {
838 minn = tdist;
839 nminn = j;
840 }
841 if (tdist > 0) if (tdist < minp) {
842 minp = tdist;
843 nminp = j;
844 }
845 }
846 }
847 if (exact) {
848 ATH_MSG_VERBOSE("exact value found for time=" << i
849 << " pulseShape=" << pulseShape[dignum * 25 / 2 + i]);
850
851 } else {
852 ATH_MSG_VERBOSE("exact value NOT found for time=" << i
853 << " nminn=" << nminn
854 << " pulseShapeT=" << pulseShapeT[nminn]
855 << " pulseShapeY=" << pulseShapeY[nminn] << std::endl
856 << " nminp=" << nminp
857 << " pulseShapeT=" << pulseShapeT[nminp]
858 << " pulseShapeY=" << pulseShapeY[nminp]);
859
860
861 pulseShape[dignum * 25 / 2 + i] = pulseShapeY[nminn]
862 + (pulseShapeY[nminp] - pulseShapeY[nminn]) / (pulseShapeT[nminp] - pulseShapeT[nminn])
863 * (i - pulseShapeT[nminn]);
864 }
865 }
866 }
867}
#define endmsg
#define ATH_MSG_INFO(x)
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_DEBUG(x)
static Double_t a
MsgStream & msg() const
The standard message stream.
bool msgLvl(const MSG::Level lvl) const
Test the output level.
AthMessaging(IMessageSvc *msgSvc, const std::string &name)
Constructor.
static const unsigned int MAX_GAIN
Number of gains per channel.
static const unsigned int MAX_DRAWER
Number of drawers in ROS 1-4.
static const unsigned int MAX_CHAN
Number of channels in drawer.
void calculateRunningCorrelation(int dignum, int chthres, bool is7to9)
void calculateCorrelation(int dignum)
void sum(std::vector< double > &digits, int ros, int drawer, int channel, int gain, int &dignum)
int m_nPairs[4][64][48][2][9]
double m_R[4][64][48][2][9][9]
double m_S22[4][64][48][2][9]
double m_corrSum2[4][64][48][2][9]
double m_S[4][64][48][2][9]
double m_S2[4][64][48][2][9]
double m_S1[4][64][48][2][9]
void saveCorrelationMatrix(bool deltaCorrelation, const std::string &correlationMatrixOptFilterFile, const TileHWID *tileHWID, int dignum)
double m_corrSum[4][64][48][2][9]
int m_N[4][64][48][2]
void calculateWeights(bool deltaCorrelation, const std::vector< double > &shapeFormLG, const std::vector< double > &shapeFormHG, const std::vector< double > &shapeFormDerivativeLG, const std::vector< double > &shapeFormDerivativeHG, const std::string &aiLoOptFilterFile, const std::string &biLoOptFilterFile, const std::string &aiHiOptFilterFile, const std::string &biHiOptFilterFile, const TileHWID *tileHWID, int dignum)
void buildPulseShape(std::vector< double > &pulseShape, std::vector< double > &pulseShapeY, std::vector< double > &pulseShapeT, int dignum)
void runningCorrelation(std::vector< double > &digits, int ros, int drawer, int channel, int gain, int &dignum, int chthres)
void saveCorrelationSumm(bool deltaCorrelation, const std::string &correlationSummOptFilterFile, const TileHWID *tileHWID, int dignum)
double m_S11[4][64][48][2][9]
void setCorrelationZero(int dignum)
void printCorrelation(int dignum)
void setCorrelationDelta(int dignum)
double m_S12[4][64][48][2][9]
double m_SS[4][64][48][2][9][9]
Helper class for TileCal online (hardware) identifiers.
Definition TileHWID.h:49
int frag(const HWIdentifier &id) const
extract frag field from HW identifier
Definition TileHWID.h:181
void zero(TH2 *h)
zero the contents of a 2d histogram