Belle II Software  release-05-01-25
eclcovmat_bit.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2011 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Calculation some parameters for files that download *
7  * at eclectronics. These parameters are discribe bit capacity *
8  * of the weighting coefficients for time and amplitude calculation *
9  * *
10  * Contributors: Alexander Bobrov (a.v.bobrov@inp.nsk.su) , *
11  * Guglielmo De Nardo *
12  * *
13  * This software is provided "as is" without any warranty. *
14  **************************************************************************/
15 
16 #include <stdio.h>
17 #include <iostream>
18 #include <cassert>
19 #include <cmath>
20 #include <cstring>
21 #include <string>
22 #include <ecl/digitization/algorithms.h>
23 #include <ecl/digitization/WrapArray2D.h>
24 
25 using namespace std;
26 using namespace Belle2;
27 using namespace ECL;
28 
29 int myPow(int x, int p)
30 {
31  if (p == 0) return 1;
32  if (p == 1) return x;
33  return x * myPow(x, p - 1);
34 }
35 
36 short int bit_val(double f, int idd, short int bit, int a, int b, int lim)
37 {
38  short int i;
39  i = bit;
40 
41  int flag;
42  flag = 0;
43  if (i < 11)
44  i = 11;
45 
46  do {
47  const int R(myPow(2, i));
48  const int val((int)(a * R * f / idd / b + R + 0.5) - R);
49 
50  if ((val < lim && val > -lim)) {
51  flag = 1;
52  } else {
53  i = i - 1;
54 
55  }
56 
57  } while (i > 8 && flag == 0);
58 
59 
60  if (i > bit) {i = bit;}
61  return i;
62 
63 }
64 
65 int PhiC(int k)
66 {
67  int R;
68  if (k > -1 && k < 8736) {
69  if (k < 96)R = 48;
70  else if (k < 288)R = 64;
71  else if (k < 864)R = 96;
72  else if (k < 8064)R = 144;
73  else if (k < 8544)R = 96;
74  else R = 64;
75  } else {R = -1;}
76  return R;
77 }
78 
79 
80 int Rest(int k)
81 {
82  int R;
83  if (k > -1 && k < 8736) {
84  if (k < 96)R = 0;
85  else if (k < 288)R = 96;
86  else if (k < 864)R = 288;
87  else if (k < 8064)R = 864;
88  else if (k < 8544)R = 8064;
89  else R = 8544;
90  } else {R = -1;}
91  return R;
92 }
93 
94 
95 int Rest1(int k)
96 {
97  int R;
98  if (k > -1 && k < 8736) {
99  if (k < 96)R = 0;
100  else if (k < 288)R = 2;
101  else if (k < 864)R = 5;
102  else if (k < 8064)R = 11;
103  else if (k < 8544)R = 61;
104  else R = 66;
105  } else {R = -1;}
106  return R;
107 }
108 
109 
110 
111 
112 int ThetR(int k)
113 {
114  int R;
115 
116  if (k > -1 && k < 8736) {
117  const int l = Rest(k);
118  const int b = PhiC(k);
119  R = (k - l) / b + Rest1(k);
120 
121  } else
122  R = -1;
123 
124  return R;
125 }
126 
127 
128 
129 
130 
131 void bit(int num, char* inputPath)
132 {
133 
134 
135  FILE* BMcoIN;
136 
137  // varible for one channel
138 
139  double ss1[16][16];
140 
141  int ChN;
142 
143 
144  double g1g1[192], gg[192], gg1[192];//, dgg[192]; dgg is never used (AH)
145  WrapArray2D<double> sg1(16, 192), sg(16, 192), sg2(16, 192);
146  double gg2[192], g1g2[192], g2g2[192];
147  double dgg1[192], dgg2[192];
148 
149  WrapArray2D<double> f(192, 16);
150  WrapArray2D<double> f1(192, 16);
151  WrapArray2D<double> fg31(192, 16);
152  WrapArray2D<double> fg32(192, 16);
153  WrapArray2D<double> fg33(192, 16);
154 
155  WrapArray2D<double> fg41(24, 16);
156  WrapArray2D<double> fg43(24, 16);
157 
158 
159  WrapArray2D<int> m_f(192, 16);
160  WrapArray2D<int> m_f1(192, 16);
161  WrapArray2D<int> m_fg31(192, 16);
162  WrapArray2D<int> m_fg32(192, 16);
163  WrapArray2D<int> m_fg33(192, 16);
164 
165  WrapArray2D<int> m_fg41(24, 16);
166  WrapArray2D<int> m_fg43(24, 16);
167 
168 
169  double val_f;
170 
171  memset(g1g1, 0, sizeof(g1g1));
172  memset(gg, 0, sizeof(gg));
173  memset(gg1, 0, sizeof(gg1));
174 // memset(dgg, 0, sizeof(dgg)); dgg is never used (AH)
175  memset(sg1, 0, sizeof(sg1));
176  memset(sg, 0, sizeof(sg));
177  memset(sg2, 0, sizeof(sg2));
178  memset(gg2, 0, sizeof(gg2));
179  memset(g1g2, 0, sizeof(g1g2));
180  memset(g2g2, 0, sizeof(g2g2));
181  memset(dgg1, 0, sizeof(dgg1));
182  memset(dgg2, 0, sizeof(dgg2));
183 
184 
185  double dt;
186 
187  double fdd;
188 
189 
190  int j1, j, i;
191  double adt, ddt, tc1;
192  double tmd, tpd, ssssj, ssssj1, ssssi, ssssi1;
193  double svp, svm;
194 
195 
196  int idd;
197  int ia, i16;
198  int mxf, mxf1, mxfg31, mxfg32, mxfg33, mxfg41, mxfg43;
199  mxf = 0;
200  mxf1 = 0;
201  mxfg31 = 0;
202  mxfg32 = 0;
203  mxfg33 = 0;
204  mxfg41 = 0;
205  mxfg43 = 0;
206 
207  int k;
208 
209  short int bitf;
210  short int bitf1;
211  short int bitfg31;
212  short int bitfg32;
213  short int bitfg33;
214  short int bitfg41;
215  short int bitfg43;
216 
217 
218  bitf = 16;
219  bitf1 = 16;
220  bitfg31 = 16;
221  bitfg41 = 16;
222 
223 
224  bitfg32 = 16;
225 
226 
227  bitfg33 = 19;
228  bitfg43 = 19;
229 
230 
231  const int mapsize = 217;
232 
233  short int bfg31[mapsize];
234 
235  short int bfg32[mapsize];
236  short int bfg33[mapsize];
237  short int bfg41[mapsize];
238 
239  short int bfg43[mapsize];
240  short int bf[mapsize];
241  short int bf1[mapsize];
242  short int bad[mapsize];
243 
244  int badN;
245  badN = 0;
246  for (j = 0; j < mapsize; j++) {
247  bfg31[j] = 16;
248  bfg32[j] = 16;
249  bfg33[j] = 19;
250  bf[j] = 16;
251  bf1[j] = 16;
252 
253  bfg41[j] = 16;
254  bfg43[j] = 19;
255  bad[j] = 0;
256  }
257  double xf;
258  double xf1;
259  double xfg31;
260  double xfg32;
261  double xfg33;
262  double xfg41;
263  double xfg43;
264 
265 
266 
267  xf = 0.;
268  xf1 = 0.;
269  xfg31 = 0.;
270  xfg32 = 0.;
271  xfg33 = 0.;
272  xfg41 = 0.;
273  xfg43 = 0.;
274 
275 
276  int cf;
277  int cf1;
278  int cf31;
279  int cf32;
280  int cf33;
281  int cf41;
282  int cf43;
283 
284  cf = 12;
285  cf1 = 12;
286  cf31 = 12;
287  cf32 = 12;
288  cf33 = 12;
289  cf41 = 12;
290  cf43 = 12;
291 
292  dt = 0.5;
293 
294 
295  const double ndt = 96.;
296  adt = 1. / ndt;
297  const int endc = 2 * ndt;
298 
299 
300  const double dt0 = adt * dt;
301  double t0 = -dt0 * ndt;
302 
303 
304 
305 
306  //%%%%%%%%%%%%%%%%%%%%%%%%%55
307 
308 
309  const double del = 0.;
310  const double ts0 = 0.5;
311  dt = 0.5;
312 
313  for (j1 = 0; j1 < endc; j1++) {
314  t0 = t0 + dt0;
315  double t = t0 - dt - del;
316 
317 // const double tin(t + dt); //tin value unused (AH)
318 
319  for (j = 0; j < 16; j++) {
320  t = t + dt;
321 
322  if (t > 0) {
323 
324 
325  f[j1][j] = ShaperDSP(t / 0.881944444);
326 
327  ddt = 0.005 * dt;
328  tc1 = t - ts0;
329  tpd = t + ddt;
330  tmd = t - ddt;
331 
332 
333  if (tc1 > ddt) {
334  svp = ShaperDSP(tpd / 0.881944444);
335  svm = ShaperDSP(tmd / 0.881944444);
336 
337  f1[j1][j] = (svp - svm) / 2. / ddt;
338 
339  } else {
340 
341 
342  f1[j1][j] = ShaperDSP(tpd / 0.881944444) / (ddt + tc1);
343 
344 
345  }// else tc1>ddt
346  // if(fabs(f1[ChN][j1][j])>100.||(j==2&&j1==156)){
347 
348  if (fabs(f1[j1][j]) > 100. || (j == 2 && j1 == 156)) {
349 
350  }
351  } //if t>0
352  else {
353  f[j1][j] = 0.;
354  f1[j1][j] = 0.;
355 
356  }
357 
358  } //for j
359 
360 
361 
362 
363 
364 
365  }
366 
367 
368 
369 
371 
372 
373  badN = 0;
374  for (ChN = 0; ChN < mapsize; ChN++) {
375  if (ChN % (mapsize / 50) == 0) {printf("CnN=%d badN=%d \n", ChN, badN);}
376  //read matrix
377 
378  {
379  string filename(inputPath);
380  //filename += "/corr";
381  filename += to_string(num);
382  filename += "/Binmcor";
383  filename += to_string(ChN);
384  filename += "_L.dat";
385  // sprintf(BMin, "/hsm/belle/bdata2/users/avbobrov/belle2/corr%d/Binmcor%d_L.dat", num, ChN);
386 
387  if ((BMcoIN = fopen(filename.c_str(), "rb")) == NULL) {
388  printf(" file %s is absent \n", filename.c_str());
389  exit(1);
390  }
391 
392  // 2
393  if (fread(ss1, sizeof(double), 256, BMcoIN) != 256) {
394  printf("Error writing file %s \n", filename.c_str());
395  exit(1);
396  }
397 
398  fclose(BMcoIN);
399 
400  }
401 
402 
403 
404  // shape function parameters
405 
406 
407 
408 
409  for (j1 = 0; j1 < endc; j1++) {
410 
411  gg[j1] = 0.;
412  gg1[j1] = 0.;
413  g1g1[j1] = 0.;
414  gg2[j1] = 0.;
415  g1g2[j1] = 0.;
416  g2g2[j1] = 0.;
417  for (j = 0; j < 16; j++) {
418  sg[j][j1] = 0.;
419  sg1[j][j1] = 0.;
420  sg2[j][j1] = 0.;
421 
422  ssssj1 = f1[j1][j];
423  ssssj = f[j1][j];
424 
425  for (i = 0; i < 16; i++) {
426 
427  sg[j][j1] = sg[j][j1] + ss1[j][i] * f[j1][i];
428  sg1[j][j1] = sg1[j][j1] + ss1[j][i] * f1[j1][i];
429 
430  sg2[j][j1] = sg2[j][j1] + ss1[j][i];
431 
432  ssssi = f[j1][i];
433  ssssi1 = f1[j1][i];
434 
435  gg[j1] = gg[j1] + ss1[j][i] * ssssj * ssssi;
436 
437 
438 
439  gg1[j1] = gg1[j1] + ss1[j][i] * ssssj * ssssi1;
440  g1g1[j1] = g1g1[j1] + ss1[j][i] * ssssi1 * ssssj1;
441 
442  gg2[j1] = gg2[j1] + ss1[j][i] * ssssj;
443  g1g2[j1] = g1g2[j1] + ss1[j][i] * ssssj1;
444  g2g2[j1] = g2g2[j1] + ss1[j][i];
445 
446 
447  } // for i
448 
449 
450  } //for j
451 
452 // dgg[j1] = gg[j1] * g1g1[j1] - gg1[j1] * gg1[j1]; dgg is never used (AH)
453  dgg1[j1] = gg[j1] * g2g2[j1] - gg2[j1] * gg2[j1];
454  dgg2[j1] = gg[j1] * g1g1[j1] * g2g2[j1] - gg1[j1] * gg1[j1] * g2g2[j1] + 2 * gg1[j1] * g1g2[j1] * gg2[j1] - gg2[j1] * gg2[j1] *
455  g1g1[j1] - g1g2[j1] * g1g2[j1] * gg[j1];
456 
457 
458  for (i = 0; i < 16; i++) {
459  if (dgg2[j1] != 0) {
460 
461  fg31[j1][i] = ((g1g1[j1] * g2g2[j1] - g1g2[j1] * g1g2[j1]) * sg[i][j1] + (g1g2[j1] * gg2[j1] - gg1[j1] * g2g2[j1]) * sg1[i][j1] +
462  (gg1[j1] * g1g2[j1] - g1g1[j1] * gg2[j1]) * sg2[i][j1]) / dgg2[j1];
463 
464 
465  fg32[j1][i] = ((g1g2[j1] * gg2[j1] - gg1[j1] * g2g2[j1]) * sg[i][j1] + (gg[j1] * g2g2[j1] - gg2[j1] * gg2[j1]) * sg1[i][j1] +
466  (gg1[j1] * gg2[j1] - gg[j1] * g1g2[j1]) * sg2[i][j1]) / dgg2[j1];
467 
468  fg33[j1][i] = ((gg1[j1] * g1g2[j1] - g1g1[j1] * gg2[j1]) * sg[i][j1] + (gg1[j1] * gg2[j1] - gg[j1] * g1g2[j1]) * sg1[i][j1] +
469  (gg[j1] * g1g1[j1] - gg1[j1] * gg1[j1]) * sg2[i][j1]) / dgg2[j1];
470 
471 
472  }
473 
474 
475  //cold
476  if (dgg1[j1] != 0) {
477 
478  fg41[j1][i] = (g2g2[j1] * sg[i][j1] - gg2[j1] * sg2[i][j1]) / dgg1[j1];
479  fg43[j1][i] = (gg[j1] * sg2[i][j1] - gg2[j1] * sg[i][j1]) / dgg1[j1];
480 
481 
482  }
483  //cold
484 
485 
486  } // for i
487 
488 
489 
490  } // for j1 <endc
491 
492 
493 
494 
495  //%%%%%%%%%%%%%%%%%%%adduction to integer
496 
497  constexpr int n16 = 16;
498 
499 
500  const int ipardsp13 = 14 + 14 * 256;
501  const int ipardsp14 = 0 * 256 + 17;
502 
503  const int ibb = ipardsp13 / 256;
504  const int iaa = ipardsp13 - ibb * 256;
505  idd = ipardsp14 / 256;
506 // const int icc = ipardsp14 - idd * 256; //Value never used (AH)
507 
508 
509 
510 
511 
512  ia = myPow(2, iaa);
513 // const int ib = myPow(2, ibb); //ib value unused (AH)
514 
515 // const int ic = myPow(2, icc); //ic value unused (AH)
516  i16 = myPow(2, 15);
517 // const int ilim = myPow(2, 15); //ilim value unused (AH)
518 
519  for (i = 0; i < 16; i++) {
520  if (i == 0) {
521  idd = n16;
522  fdd = 16.;
523  } else {
524  idd = 1;
525  fdd = 1.;
526  }
527  for (k = 0; k < 192; k++) {
528 
529  val_f = f[k][i];
530  bf[ChN] = bit_val(val_f, idd, bf[ChN], 1, 1, i16);
531  if (bf[ChN] < bitf) {bitf = bf[ChN];}
532  val_f = f1[k][i];
533 
534  bf1[ChN] = bit_val(val_f, idd, bf1[ChN], 4, 3, i16);
535  if (bf1[ChN] < bitf1) {bitf1 = bf1[ChN];}
536  val_f = fg31[k][i];
537  bfg31[ChN] = bit_val(val_f, idd, bfg31[ChN], 1, 1, i16);
538  if (bfg31[ChN] < bitfg31) {bitfg31 = bfg31[ChN];}
539  val_f = fg32[k][i];
540  bfg32[ChN] = bit_val(val_f, idd, bfg32[ChN], 3, 4, i16);
541  if (bfg32[ChN] < bitfg32) {bitfg32 = bfg32[ChN];}
542  val_f = fg33[k][i];
543  bfg33[ChN] = bit_val(val_f, idd, bfg33[ChN], 1, 1, i16);
544  if (bfg33[ChN] < bitfg33) {bitfg33 = bfg33[ChN];}
545  if (k < 24) {
546  val_f = fg41[k][i];
547  bfg41[ChN] = bit_val(val_f, idd, bfg41[ChN], 1, 1, i16);
548  if (bfg41[ChN] < bitfg41) {bitfg41 = bfg41[ChN];}
549  val_f = fg43[k][i];
550  bfg43[ChN] = bit_val(val_f, idd, bfg43[ChN], 1, 1, i16);
551  if (bfg43[ChN] < bitfg43) {bitfg43 = bfg43[ChN];}
552  }
553 
554 
555 
556  if (fabs(f[k][i] / fdd) > xf) {xf = fabs(f[k][i] / fdd);}
557  if (fabs(f1[k][i] / fdd) > xf1) {xf1 = fabs(f1[k][i] / fdd);}
558  if (fabs(fg31[k][i] / fdd) > xfg31) {xfg31 = fabs(fg31[k][i] / fdd);}
559  if (fabs(fg32[k][i] / fdd) > xfg32) {xfg32 = fabs(fg32[k][i] / fdd);}
560  if (fabs(fg33[k][i] / fdd) > xfg33) {xfg33 = fabs(fg33[k][i] / fdd);}
561  if (k < 24) {
562  if (fabs(fg41[k][i] / fdd) > xfg41) {xfg41 = fabs(fg41[k][i] / fdd);}
563  if (fabs(fg43[k][i] / fdd) > xfg43) {xfg43 = fabs(fg43[k][i] / fdd);}
564  }
565 
566 
567  } // for k
568 
569 
570  } // for i
571 
572 
573  if (bf[ChN] < cf || bf1[ChN] < cf1 || bfg31[ChN] < cf31 || bfg32[ChN] < cf32 || bfg33[ChN] < cf33 || bfg41[ChN] < cf41
574  || bfg43[ChN] < cf43) {
575  badN = badN + 1;
576 
577  if (bf[ChN] < cf) {bad[ChN] = bad[ChN] + 1;}
578  if (bf1[ChN] < cf1) {bad[ChN] = bad[ChN] + 2;}
579  if (bfg31[ChN] < cf31) {bad[ChN] = bad[ChN] + 4;}
580  if (bfg32[ChN] < cf32) {bad[ChN] = bad[ChN] + 8;}
581  if (bfg33[ChN] < cf33) {bad[ChN] = bad[ChN] + 16;}
582  if (bfg41[ChN] < cf41) {bad[ChN] = bad[ChN] + 32;}
583  if (bfg43[ChN] < cf43) {bad[ChN] = bad[ChN] + 64;}
584 
585  }
586  }
587 
588  // printf("mxf=%d mxf1=%d mxfg31=%d mxfg32=%d mxfg33=%d mxfg41=%d mxfg43=%d \n",mxf,mxf1,mxfg31,mxfg32,mxfg33,mxfg41,mxfg43);
589 
590 
591  printf("bxf=%d bxf1=%d bxfg31=%d bxfg32=%d bxfg33=%d bxfg41=%d bxfg43=%d \n", bitf, bitf1, bitfg31, bitfg32, bitfg33, bitfg41,
592  bitfg43);
593  printf("xf=%lf xf1=%lf xfg31=%lf xfg32=%lf xfg33=%lf xfg41=%lf xfg43=%lf \n", xf, xf1, xfg31, xfg32, xfg33, xfg41, xfg43);
594 
595  ia = myPow(2, bitf);
596  mxf = (int)(xf * ia / idd + ia + 0.5) - ia;
597 
598  ia = myPow(2, bitf1);
599  mxf1 = (int)(4 * xf1 * ia / idd / 3 + ia + 0.5) - ia;
600  ia = myPow(2, bitfg31);
601  mxfg31 = (int)(xfg31 * ia / idd + ia + 0.5) - ia;
602  ia = myPow(2, bitfg32);
603  mxfg32 = (int)(3 * xfg32 * ia / idd / 4 + ia + 0.5) - ia;
604  ia = myPow(2, bitfg33);
605  mxfg33 = (int)(xfg33 * ia / idd + ia + 0.5) - ia;
606  ia = myPow(2, bitfg41);
607  mxfg41 = (int)(xfg41 * ia / idd + ia + 0.5) - ia;
608  ia = myPow(2, bitfg43);
609  mxfg43 = (int)(xfg43 * ia / idd + ia + 0.5) - ia;
610 
611  printf("mf=%d mf1=%d mfg31=%d mfg32=%d mfg33=%d mfg41=%d mfg43=%d \n", mxf, mxf1, mxfg31, mxfg32, mxfg33, mxfg41, mxfg43);
612 
613 
614  printf("Rf=%lf Rf1=%lf Rfg31=%lf Rfg32=%lf Rfg33=%lf Rfg41=%lf Rfg43=%lf \n", (double)mxf / (double)i16, (double)mxf1 / (double)i16,
615  (double)mxfg31 / (double)i16, (double)mxfg32 / (double)i16, (double)mxfg33 / (double)i16, (double)mxfg41 / (double)i16,
616  (double)mxfg43 / (double)i16);
617  printf("\n");
618 
619 
620  i16 = myPow(2, 15);
621 
622  if (bitf == 16) {
623  ia = myPow(2, bitf);
624  mxf = (int)(xf * ia / idd + ia + 0.5) - ia;
625 
626  printf("f at limit %lf\n", (double)mxf / (double)i16);
627  if ((double)mxf / (double)i16 > 1) {printf("wronf bit for f\n");}
628  } else {
629  ia = myPow(2, bitf);
630  mxf = (int)(xf * ia / idd + ia + 0.5) - ia;
631 
632  printf("f %lf \n", (double)mxf / (double)i16);
633  if ((double)mxf / (double)i16 > 1) {printf("wronf bit for f\n");}
634 
635 
636  ia = myPow(2, bitf + 1);
637  mxf = (int)(xf * ia / idd + ia + 0.5) - ia;
638 
639  printf("f+1 %lf\n", (double)mxf / (double)i16);
640 
641  ia = myPow(2, bitf - 1);
642  mxf = (int)(xf * ia / idd + ia + 0.5) - ia;
643 
644  printf("f-1 %lf\n", (double)mxf / (double)i16);
645 
646 
647  }
648 
649 
650  printf("\n");
651 
652  if (bitf1 == 16) {
653  ia = myPow(2, bitf1);
654  mxf1 = (int)(4 * xf1 * ia / idd / 3 + ia + 0.5) - ia;
655 
656  printf("f at limit %lf\n", (double)mxf1 / (double)i16);
657  if ((double)mxf / (double)i16 > 1) {printf("wronf bit for f1\n");}
658  } else {
659  ia = myPow(2, bitf1);
660  mxf1 = (int)(4 * xf1 * ia / idd / 3 + ia + 0.5) - ia;
661 
662  printf("f1 %lf\n", (double)mxf1 / (double)i16);
663  if ((double)mxf1 / (double)i16 > 1) {printf("wronf bit for f1\n");}
664 
665 
666  ia = myPow(2, bitf1 + 1);
667  mxf1 = (int)(4 * xf1 * ia / idd / 3 + ia + 0.5) - ia;
668 
669  printf("f1+1 %lf \n", (double)mxf1 / (double)i16);
670 
671  ia = myPow(2, bitf1 - 1);
672  printf("ia=%d i16=%d \n", ia, i16);
673 
674  mxf1 = (int)(4 * xf1 * ia / idd / 3 + ia + 0.5) - ia;
675 
676  printf("f1-1 %lf\n", (double)mxf1 / (double)i16);
677 
678 
679  }
680 
681 
682  printf("\n");
683 
684  if (bitfg31 == 16) {
685  ia = myPow(2, bitfg31);
686  mxfg31 = (int)(xfg31 * ia / idd + ia + 0.5) - ia;
687 
688  printf("fg31 at limit %lf\n", (double)mxfg31 / (double)i16);
689  if ((double)mxfg31 / (double)i16 > 1) {printf("wronf bit for fg31\n");}
690  } else {
691  ia = myPow(2, bitfg31);
692  mxfg31 = (int)(xfg31 * ia / idd + ia + 0.5) - ia;
693 
694  printf("fg31 %lf\n", (double)mxfg31 / (double)i16);
695  if ((double)mxfg31 / (double)i16 > 1) {printf("wronf bit for fg31\n");}
696 
697 
698  ia = myPow(2, bitfg31 + 1);
699  mxfg31 = (int)(xfg31 * ia / idd + ia + 0.5) - ia;
700 
701  printf("fg31+1 %lf \n", (double)mxfg31 / (double)i16);
702 
703  ia = myPow(2, bitfg31 - 1);
704  mxfg31 = (int)(xfg31 * ia / idd + ia + 0.5) - ia;
705 
706  printf("fg31-1 %lf\n", (double)mxfg31 / (double)i16);
707 
708 
709  }
710 
711 
712  printf("\n");
713 
714  if (bitfg32 == 16) {
715  ia = myPow(2, bitfg32);
716  mxfg31 = (int)(3 * xfg32 * ia / idd / 4 + ia + 0.5) - ia;
717 
718  printf("fg32 at limit %lf\n", (double)mxfg32 / (double)i16);
719  if ((double)mxfg32 / (double)i16 > 1) {printf("wronf bit for f32\n");}
720  } else {
721  ia = myPow(2, bitfg32);
722  mxfg32 = (int)(3 * xfg32 * ia / idd / 4 + ia + 0.5) - ia;
723 
724  printf("fg32 %lf\n", (double)mxfg32 / (double)i16);
725  if ((double)mxfg32 / (double)i16 > 1) {printf("wronf bit for f\n");}
726 
727 
728  ia = myPow(2, bitfg32 + 1);
729  mxfg31 = (int)(3 * xfg32 * ia / idd / 4 + ia + 0.5) - ia;
730 
731  printf("fg32+1 %lf \n", (double)mxfg31 / (double)i16);
732 
733  ia = myPow(2, bitfg31 - 1);
734  mxfg31 = (int)(3 * xfg32 * ia / idd / 4 + ia + 0.5) - ia;
735 
736  printf("fg32-1 %lf\n", (double)mxfg31 / (double)i16);
737 
738 
739  }
740 
741 
742 
743 
744  printf("\n");
745 
746  if (bitfg33 == 19) {
747  ia = myPow(2, bitfg33);
748  mxfg33 = (int)(xfg33 * ia / idd + ia + 0.5) - ia;
749 
750  printf("fg33 at limit %lf\n", (double)mxfg33 / (double)i16);
751  if ((double)mxfg33 / (double)i16 > 1) {printf("wronf bit for fg33\n");}
752  } else {
753  ia = myPow(2, bitfg33);
754  mxfg33 = (int)(xfg33 * ia / idd + ia + 0.5) - ia;
755 
756  printf("fg33 %lf\n", (double)mxfg33 / (double)i16);
757  if ((double)mxfg33 / (double)i16 > 1) {printf("wronf bit for fg33\n");}
758 
759 
760  ia = myPow(2, bitfg33 + 1);
761  mxfg33 = (int)(xfg33 * ia / idd + ia + 0.5) - ia;
762 
763  printf("fg33+1 %lf \n", (double)mxfg33 / (double)i16);
764 
765  ia = myPow(2, bitfg33 - 1);
766  mxfg33 = (int)(xfg33 * ia / idd + ia + 0.5) - ia;
767 
768  printf("fg33-1 %lf\n", (double)mxfg33 / (double)i16);
769 
770 
771  }
772 
773 
774  printf("\n");
775 
776 
777  if (bitfg41 == 16) {
778  ia = myPow(2, bitfg41);
779  mxfg41 = (int)(xfg41 * ia / idd + ia + 0.5) - ia;
780 
781  printf("fg41 at limit %lf\n", (double)mxfg41 / (double)i16);
782  if ((double)mxfg41 / (double)i16 > 1) {printf("wronf bit for fg41\n");}
783  } else {
784  ia = myPow(2, bitfg41);
785  mxfg41 = (int)(xfg41 * ia / idd + ia + 0.5) - ia;
786 
787  printf("fg41 %lf\n", (double)mxfg41 / (double)i16);
788  if ((double)mxfg41 / (double)i16 > 1) {printf("wronf bit for fg41\n");}
789 
790 
791  ia = myPow(2, bitfg41 + 1);
792  mxfg41 = (int)(xfg41 * ia / idd + ia + 0.5) - ia;
793 
794  printf("fg41+1 %lf \n", (double)mxfg41 / (double)i16);
795 
796  ia = myPow(2, bitfg41 - 1);
797  mxfg41 = (int)(xfg41 * ia / idd + ia + 0.5) - ia;
798 
799  printf("fg41-1 %lf\n", (double)mxfg41 / (double)i16);
800 
801 
802  }
803 
804 
805  printf("\n");
806 
807  if (bitfg43 == 19) {
808  ia = myPow(2, bitfg43);
809  mxfg43 = (int)(xfg43 * ia / idd + ia + 0.5) - ia;
810 
811  printf("fg43 at limit %lf\n", (double)mxfg43 / (double)i16);
812  if ((double)mxfg43 / (double)i16 > 1) {printf("wronf bit for fg33\n");}
813  } else {
814  ia = myPow(2, bitfg43);
815  mxfg43 = (int)(xfg43 * ia / idd + ia + 0.5) - ia;
816 
817  printf("fg43 %lf\n", (double)mxfg43 / (double)i16);
818  if ((double)mxfg33 / (double)i16 > 1) {printf("wronf bit for fg33\n");}
819 
820 
821  ia = myPow(2, bitfg43 + 1);
822  mxfg43 = (int)(xfg43 * ia / idd + ia + 0.5) - ia;
823 
824  printf("fg43+1 %lf \n", (double)mxfg43 / (double)i16);
825 
826  ia = myPow(2, bitfg43 - 1);
827  mxfg43 = (int)(xfg43 * ia / idd + ia + 0.5) - ia;
828 
829  printf("fg43-1 %lf\n", (double)mxfg43 / (double)i16);
830 
831 
832  }
833 
834 
835 
836 
837  printf("\n");
838 
839 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
840 
841  if (bitf < 13) {printf("problem bitf=%d \n", bitf);}
842  if (bitf1 < 13) {printf("problem bitf1=%d \n", bitf1);}
843  if (bitfg31 < 13) {printf("problem bitfg31=%d \n", bitfg31);}
844  if (bitfg32 < 13) {printf("problem bitfg32=%d \n", bitfg32);}
845  if (bitfg33 < 16) {printf("problem bitfg33=%d \n", bitfg33);}
846 
847  if (bitfg41 < 13) {printf("problem bitfg41=%d \n", bitfg41);}
848  if (bitfg43 < 16) {printf("problem bitfg43=%d \n", bitfg43);}
849 
850 
851  printf("\n");
852 
853  printf("bxf=%d bxf1=%d bxfg31=%d bxfg32=%d bxfg33=%d bxfg41=%d bxfg43=%d \n", bitf, bitf1, bitfg31, bitfg32, bitfg33, bitfg41,
854  bitfg43);
855 
856  printf("badN=%d \n", badN);
857 
858 // int RF[69]; RF is never used (AH)
859 // int RF1[69]; RF1 is never used (AH)
860  int RFG31[69];
861  int RFG32[69];
862  int RFG33[69];
863 // int RFG41[69]; RFG41 is never used (AH)
864  int RFG43[69];
865 
866 
867  int BF[20];
868  int BF1[20];
869  int BFG31[20];
870  int BFG32[20];
871  int BFG33[20];
872  int BFG41[20];
873  int BFG43[20];
874 
875  for (i = 0; i < 69; i++) {
876  if (i < 20) {
877  BF[i] = 0;
878  BF1[i] = 0;
879  BFG31[i] = 0;
880  BFG32[i] = 0;
881  BFG32[i] = 0;
882  BFG31[i] = 0;
883  BFG43[i] = 0;
884  }
885 
886 // RF[i] = 0; RF is never used (AH)
887 // RF1[i] = 0; RF is never used (AH)
888  RFG31[i] = 0;
889  RFG32[i] = 0;
890  RFG32[i] = 0;
891  RFG31[i] = 0;
892  RFG43[i] = 0;
893 
894  }
895 
896  for (j = 0; j < mapsize; j++) {
897 
898  for (i = 10; i < 20; i++) {
899  if (bf[j] == i) {BF[i] = BF[i] + 1;}
900  if (bf1[j] == i) {BF1[i] = BF1[i] + 1;}
901  if (bfg31[j] == i) {BFG31[i] = BFG31[i] + 1;}
902  if (bfg32[j] == i) {BFG32[i] = BFG32[i] + 1;}
903  if (bfg33[j] == i) {BFG33[i] = BFG33[i] + 1;}
904  if (bfg41[j] == i) {BFG41[i] = BFG41[i] + 1;}
905  if (bfg43[j] == i) {BFG43[i] = BFG43[i] + 1;}
906  }
907  if (bf[j] < cf || bf1[j] < cf1 || bfg31[j] < cf31 || bfg32[j] < cf32 || bfg33[j] < cf33 || bfg41[j] < cf41 || bfg43[j] < cf43) {
908  if (bfg31[j] < cf31) { //printf("fg31 ring=%d ",ThetR(j));
909  RFG31[ThetR(j)] = RFG31[ThetR(j)] + 1;
910  }
911  if (bfg32[j] < cf32) { //printf("fg32 ring=%d ",ThetR(j));
912  RFG32[ThetR(j)] = RFG32[ThetR(j)] + 1;
913  }
914  if (bfg33[j] < cf33) { //printf("fg33 ring=%d ",ThetR(j));
915  RFG33[ThetR(j)] = RFG33[ThetR(j)] + 1;
916  }
917  if (bfg43[j] < cf43) { //printf("fg41 ring=%d ",ThetR(j));
918  RFG43[ThetR(j)] = RFG43[ThetR(j)] + 1;
919  }
920  // printf("\n");
921  }
922  }
923 
924  for (i = 10; i < 20; i++) {
925  if (BF[i] > 0 || BF1[i] > 0 || BFG31[i] > 0 || BFG32[i] > 0 || BFG33[i] > 0 || BFG41[i] > 0 || BFG43[i] > 0) {
926  printf("bit=%d BF=%d BF1=%d BFG31=%d BFG32=%d BFG33=%d BFG41=%d BFG43=%d\n", i, BF[i], BF1[i], BFG31[i], BFG32[i], BFG33[i],
927  BFG41[i], BFG43[i]);
928  }
929 
930  }
931 
932  char BMin[256];
933  sprintf(BMin, "bitst%d.dat", num);
934 
935  if ((BMcoIN = fopen(BMin, "w")) == NULL) {
936  printf(" file %s is absent \n", BMin);
937  exit(1);
938  }
939 
940  for (i = 0; i < mapsize; i++) {
941  fprintf(BMcoIN, "%d %d %d %d %d %d %d %d \n", i, bf[i], bf1[i], bfg31[i], bfg32[i], bfg33[i], bfg41[i], bfg43[i]);
942  }
943  fclose(BMcoIN);
944 
945 }
946 
947 
948 int main(int argc, char** argv)
949 {
950  assert(argc == 1 || argc == 3);
951  if (argc == 1) {
952  cout << "Usage:" << endl;
953  cout << argv[0] << " <type> <covmat_path>" << endl;
954  cout << "type is a numerical integer value to define the source of the calibration" << endl;
955  cout << "covmat_path is the path to look for the directory corr<type> that contains the covariance matrix files Binmcorr<type>_L.dat"
956  << endl;
957  cout << " An output file with the translation to integers in written in the work directory" << endl;
958  } else bit(atoi(argv[1]), argv[2]);
959 }
Belle2::ECL::WrapArray2D
class to replace POD 2D array to save stack usage since it just allocates memory dynamically.
Definition: WrapArray2D.h:33
main
int main(int argc, char **argv)
Run all tests.
Definition: test_main.cc:77
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
alignment.constraints_generator.filename
filename
File name.
Definition: constraints_generator.py:224