Belle II Software  release-05-01-25
ECLDspUtilities.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2019 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Sergei Gribanov, Mikhail Remnev *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 // ECL
12 #include <ecl/utility/ECLDspUtilities.h>
13 #include <ecl/utility/ECLChannelMapper.h>
14 #include <ecl/dbobjects/ECLDspData.h>
15 #include <ecl/dbobjects/ECLCrystalCalib.h>
16 #include <ecl/dataobjects/ECLDigit.h>
17 // Framework
18 #include <framework/logging/Logger.h>
19 #include <framework/database/DBObjPtr.h>
20 #include <framework/database/DBArray.h>
21 #include <framework/utilities/FileSystem.h>
22 // ROOT
23 #include <TFile.h>
24 #include <TTree.h>
25 
26 using namespace Belle2;
27 using namespace ECL;
28 
30 const short ECLDSP_FORMAT_VERSION = 1;
31 
32 int ECLDspUtilities::pedestal_fit_initialized = 0;
33 float ECLDspUtilities::pedfit_fg31[768] = {};
34 float ECLDspUtilities::pedfit_fg32[768] = {};
35 
40 void readEclDspCoefs(FILE* fl, void* ptr, int word_size, int word_count)
41 {
42  int read_items = fread(ptr, word_size, word_count, fl);
43  if (read_items != word_count) {
44  B2ERROR("Error reading DSP coefficients");
45  }
46 }
47 
48 ECLDspData* ECLDspUtilities::readEclDsp(const char* filename, int boardNumber)
49 {
50  FILE* fl;
51  fl = fopen(filename, "rb");
52  if (fl == NULL) {
53  B2ERROR("Can't open file " << filename);
54  }
55 
56  ECLDspData* data = new ECLDspData(boardNumber);
57 
58  //== Do not fill data for non-existent shapers
59  int crateNum = (boardNumber - 1) / 12 + 1;
60  int shaperNum = (boardNumber - 1) % 12;
61 
62  if (shaperNum >= 8) {
63  if (crateNum >= 45) {
64  fclose(fl);
65  return data;
66  }
67  if (shaperNum >= 10) {
68  if (crateNum >= 37 && crateNum <= 44) {
69  fclose(fl);
70  return data;
71  }
72  }
73  }
74 
75  // Word size
76  const int nsiz = 2;
77  // Size of fragment to read (in words)
78  int nsiz1 = 256;
79  short int id[256];
80  int size = fread(id, nsiz, nsiz1, fl);
81  if (size != nsiz1) {
82  B2ERROR("Error reading header of DSP file");
83  }
84 
85  std::vector<short int> extraData;
86  extraData.push_back(ECLDSP_FORMAT_VERSION);
87  extraData.push_back(data->getPackerVersion());
88  data->setExtraData(extraData);
89 
90  data->setverMaj(id[8] & 0xFF);
91  data->setverMin(id[8] >> 8);
92  data->setkb(id[13] >> 8);
93  data->setka(id[13] - 256 * data->getkb());
94  data->sety0Startr(id[14] >> 8);
95  data->setkc(id[14] - 256 * data->gety0Startr());
96  data->setchiThresh(id[15]);
97  data->setk2(id[16] >> 8);
98  data->setk1(id[16] - 256 * data->getk2());
99  data->sethT(id[64]);
100  data->setlAT(id[128]);
101  data->setsT(id[192]);
102  data->setaAT(id[208]);
103 
104  std::vector<short int> f(49152), f1(49152), f31(49152),
105  f32(49152), f33(49152), f41(6144), f43(6144);
106 
107  for (int i = 0; i < 16; ++i) {
108  nsiz1 = 384;
109  readEclDspCoefs(fl, &(*(f41.begin() + i * nsiz1)), nsiz, nsiz1);
110  nsiz1 = 3072;
111  readEclDspCoefs(fl, &(*(f31.begin() + i * nsiz1)), nsiz, nsiz1);
112  readEclDspCoefs(fl, &(*(f32.begin() + i * nsiz1)), nsiz, nsiz1);
113  readEclDspCoefs(fl, &(*(f33.begin() + i * nsiz1)), nsiz, nsiz1);
114  nsiz1 = 384;
115  readEclDspCoefs(fl, &(*(f43.begin() + i * nsiz1)), nsiz, nsiz1);
116  nsiz1 = 3072;
117  readEclDspCoefs(fl, &(*(f.begin() + i * nsiz1)), nsiz, nsiz1);
118  readEclDspCoefs(fl, &(*(f1.begin() + i * nsiz1)), nsiz, nsiz1);
119  }
120  fclose(fl);
121 
122  data->setF41(f41);
123  data->setF31(f31);
124  data->setF32(f32);
125  data->setF33(f33);
126  data->setF43(f43);
127  data->setF(f);
128  data->setF1(f1);
129 
130  return data;
131 }
132 
133 void ECLDspUtilities::writeEclDsp(const char* filename, ECLDspData* data)
134 {
135  // Default header for DSP file.
136  // Words '0xABCD' are overwitten with current parameters.
137  const unsigned short DEFAULT_HEADER[256] {
138  // ECLDSP FILE.....
139  0x4543, 0x4c44, 0x5350, 0x2046, 0x494c, 0x4500, 0x0000, 0x0000,
140  0xABCD, 0xffff, 0x0000, 0x0000, 0x0000, 0xABCD, 0xABCD, 0xABCD,
141  0xABCD, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
142  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
143  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
144  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
145  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
146  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
147  0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD,
148  0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD,
149  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
150  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
151  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
152  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
153  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
154  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
155  0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD,
156  0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD,
157  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
158  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
159  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
160  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
161  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
162  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
163  0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD,
164  0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD,
165  0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD,
166  0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD, 0xABCD,
167  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
168  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
169  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
170  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
171  };
172 
173  FILE* fl;
174  fl = fopen(filename, "wb");
175  // Word size
176  const int nsiz = 2;
177  // Size of fragment to write (in words)
178  int nsiz1 = 256;
179  // modification of DEFAULT_HEADER
180  unsigned short header[256];
181 
182  int format_version = data->getExtraData()[0];
183  if (format_version > ECLDSP_FORMAT_VERSION) {
184  B2WARNING("Format version " << format_version << " is not fully supported, some data might be discarded.");
185  }
186 
187  for (int i = 0; i < 256; i++) {
188  if (i == 8) {
189  header[i] = (data->getverMin() << 8) | data->getverMaj();
190  } else if (i == 13) {
191  header[i] = (data->getkb() << 8) | data->getka();
192  } else if (i == 14) {
193  header[i] = (data->gety0Startr() << 8) | data->getkc();
194  } else if (i == 15) {
195  header[i] = data->getchiThresh();
196  } else if (i == 16) {
197  header[i] = (data->getk2() << 8) | data->getk1();
198  } else if (i >= 64 && i < 80) {
199  header[i] = data->gethT();
200  } else if (i >= 128 && i < 144) {
201  header[i] = data->getlAT();
202  } else if (i >= 192 && i < 208) {
203  header[i] = data->getsT();
204  } else if (i >= 208 && i < 224) {
205  header[i] = data->getaAT();
206  } else {
207  // Reverse bytes for header.
208  int high = (DEFAULT_HEADER[i] & 0xFF00) >> 8;
209  int low = (DEFAULT_HEADER[i] & 0x00FF);
210  header[i] = (low << 8) + high;
211  }
212  }
213 
214  int size = fwrite(header, nsiz, nsiz1, fl);
215  if (size != nsiz1) {
216  B2FATAL("Error writing header of DSP file " << filename);
217  }
218 
219  std::vector<short int> f(49152), f1(49152), f31(49152),
220  f32(49152), f33(49152), f41(6144), f43(6144);
221  data->getF41(f41);
222  data->getF31(f31);
223  data->getF32(f32);
224  data->getF33(f33);
225  data->getF43(f43);
226  data->getF(f);
227  data->getF1(f1);
228 
229  for (int i = 0; i < 16; ++i) {
230  nsiz1 = 384;
231  fwrite(&(*(f41.begin() + i * nsiz1)), nsiz, nsiz1, fl);
232  nsiz1 = 3072;
233  fwrite(&(*(f31.begin() + i * nsiz1)), nsiz, nsiz1, fl);
234  fwrite(&(*(f32.begin() + i * nsiz1)), nsiz, nsiz1, fl);
235  fwrite(&(*(f33.begin() + i * nsiz1)), nsiz, nsiz1, fl);
236  nsiz1 = 384;
237  fwrite(&(*(f43.begin() + i * nsiz1)), nsiz, nsiz1, fl);
238  nsiz1 = 3072;
239  fwrite(&(*(f.begin() + i * nsiz1)), nsiz, nsiz1, fl);
240  fwrite(&(*(f1.begin() + i * nsiz1)), nsiz, nsiz1, fl);
241  }
242  fclose(fl);
243 }
244 
245 /**********************************************/
246 /* SHAPE FITTER */
247 
248 short int* vectorsplit(std::vector<short int>& vectorFrom, int channel)
249 {
250  size_t size = vectorFrom.size();
251  if (size % 16) B2ERROR("Split is impossible!" << "Vector size" << size);
252  return (vectorFrom.data() + (size / 16) * (channel - 1));
253 }
254 
255 ECLShapeFit ECLDspUtilities::shapeFitter(int cid, std::vector<int> adc, int ttrig)
256 {
257  ECLChannelMapper mapper;
258  mapper.initFromDB();
259 
260  //== Get mapping data
261 
262  int crate = mapper.getCrateID(cid);
263  int shaper = mapper.getShaperPosition(cid);
264  int channel = mapper.getShaperChannel(cid);
265 
266  //== Get DSP data
267 
268  int dsp_group = (crate - 1) / 18;
269  int dsp_id = (crate - 1) % 18;
270  std::string payload_name = "ECLDSPPars" + std::to_string(dsp_group);
271  DBArray<ECLDspData> dsp_data(payload_name);
272  const ECLDspData* data = dsp_data[dsp_id * 12 + shaper - 1];
273 
274  std::vector<short> vec_f; data->getF(vec_f);
275  std::vector<short> vec_f1; data->getF1(vec_f1);
276  std::vector<short> vec_fg31; data->getF31(vec_fg31);
277  std::vector<short> vec_fg32; data->getF32(vec_fg32);
278  std::vector<short> vec_fg33; data->getF33(vec_fg33);
279  std::vector<short> vec_fg41; data->getF41(vec_fg41);
280  std::vector<short> vec_fg43; data->getF43(vec_fg43);
281 
282  short* f = vectorsplit(vec_f, channel);
283  short* f1 = vectorsplit(vec_f1, channel);
284  short* fg31 = vectorsplit(vec_fg31, channel);
285  short* fg32 = vectorsplit(vec_fg32, channel);
286  short* fg33 = vectorsplit(vec_fg33, channel);
287  short* fg41 = vectorsplit(vec_fg41, channel);
288  short* fg43 = vectorsplit(vec_fg43, channel);
289 
290  int k_a = data->getka();
291  int k_b = data->getkb();
292  int k_c = data->getkc();
293  int k_1 = data->getk1();
294  int k_2 = data->getk2();
295  int k_16 = data->gety0Startr();
296  int chi_thres = data->getchiThresh();
297 
298  //== Get ADC data
299  int* y = adc.data();
300 
301  //== Get trigger time
302  int ttrig2 = ttrig - 2 * (ttrig / 8);
303 
304  //== Get thresholds
305  DBObjPtr<ECLCrystalCalib> thr_LowAmp("ECL_FPGA_LowAmp");
306  DBObjPtr<ECLCrystalCalib> thr_HitThresh("ECL_FPGA_HitThresh");
307  DBObjPtr<ECLCrystalCalib> thr_StoreDigit("ECL_FPGA_StoreDigit");
308 
309  int A0 = thr_LowAmp->getCalibVector()[cid - 1];
310  int Ahard = thr_HitThresh->getCalibVector()[cid - 1];
311  int Askip = thr_StoreDigit->getCalibVector()[cid - 1];
312 
313  //== Perform fit
314  auto result = lftda_(f, f1, fg41, fg43, fg31, fg32, fg33, y, ttrig2, A0,
315  Ahard, Askip, k_a, k_b, k_c, k_16, k_1, k_2,
316  chi_thres);
317 
318  return result;
319 }
320 
321 void ECLDspUtilities::initPedestalFit()
322 {
323  std::string path = FileSystem::findFile("ecl/data/ecl_pedestal_peak_fit.root");
324  TFile* file = new TFile(path.c_str(), "read");
325  if (!file->IsOpen()) {
326  B2FATAL("Unable to load coefficients for ECL pedestal fit");
327  }
328  TTree* tree = (TTree*)file->Get("dsp_coefs");
329  int nentries = tree->GetEntries();
330  float fg31_i, fg32_i;
331  tree->SetBranchAddress("fg31", &fg31_i);
332  tree->SetBranchAddress("fg32", &fg32_i);
333  //== Load DSP coefficients used in pedestal fitting
334  for (int i = 0; i < nentries; i++) {
335  tree->GetEntry(i);
336  pedfit_fg31[i] = fg31_i;
337  pedfit_fg32[i] = fg32_i;
338  }
339  file->Close();
340 
341  pedestal_fit_initialized = 1;
342 }
343 
344 ECLPedestalFit ECLDspUtilities::pedestalFit(std::vector<int> adc)
345 {
346  if (!pedestal_fit_initialized) {
347  initPedestalFit();
348  }
349 
350  ECLPedestalFit result;
351  float amp;
352 
353  //== Find maximum in the pedestal
354 
355  int ped_max = 0;
356  int ped_max_pos = 0;
357  for (int i = 0; i < 16; i++) {
358  if (adc[i] > ped_max) {
359  ped_max = adc[i];
360  ped_max_pos = i;
361  }
362  }
363 
364  // Skip fit if the peak position is on the edge of the fit region
365  if (ped_max_pos < 2 || ped_max_pos > 13) {
366  result.amp = -128;
367  return result;
368  }
369 
370  //== Get first position estimate from maximum location
371 
372  int time_index = ped_max_pos * 4 - 8;
373  if (time_index > 47) time_index = 47;
374  if (time_index < 0) time_index = 0;
375 
376  //== Run two iterations of chi2 minimization
377  for (int iter = 0; iter < 2; iter++) {
378  amp = 0;
379  float tim = 0;
380  for (int j = 0; j < 16; j++) {
381  amp += pedfit_fg31[j + time_index * 16] * adc[j];
382  tim += pedfit_fg32[j + time_index * 16] * adc[j];
383  }
384  tim = tim / amp;
385  time_index -= tim * 4;
386  if (time_index > 47) time_index = 47;
387  if (time_index < 0) time_index = 0;
388  }
389 
390  result.amp = amp;
391 
392  return result;
393 }
394 
Belle2::ECL::ECLDspUtilities::writeEclDsp
static void writeEclDsp(const char *raw_file, ECLDspData *obj)
Convert ECLDspData from Root object to *.dat file.
Definition: ECLDspUtilities.cc:133
Belle2::ECL::ECLDspUtilities::pedfit_fg32
static float pedfit_fg32[768]
DSP coefficients used to determine time in pedestalFit.
Definition: ECLDspUtilities.h:108
Belle2::ECL::ECLChannelMapper::getCrateID
int getCrateID(int iCOPPERNode, int iFINESSE)
get crate number by given COPPER node number and FINESSE number
Definition: ECLChannelMapper.cc:211
Belle2::ECL::ECLChannelMapper::getShaperPosition
int getShaperPosition(int cellID)
get position of the shaper in the crate by given CellId
Definition: ECLChannelMapper.cc:289
Belle2::DBArray
Class for accessing arrays of objects in the database.
Definition: DBArray.h:36
Belle2::ECL::ECLDspUtilities::shapeFitter
static ECLShapeFit shapeFitter(int cid, std::vector< int > adc, int ttrig)
Emulate shape fitting algorithm from ShaperDSP using algorithm from ecl/utility/src/ECLDspEmulator....
Definition: ECLDspUtilities.cc:255
Belle2::ECL::ECLChannelMapper::getShaperChannel
int getShaperChannel(int cellID)
get number of DSP channel in the shaper by given number of CellId
Definition: ECLChannelMapper.cc:296
Belle2::ECL::ECLDspUtilities::initPedestalFit
static void initPedestalFit()
Load DSP coefficients used in the pedestal fit function.
Definition: ECLDspUtilities.cc:321
Belle2::ECL::ECLShapeFit
ShaperDSP fit results from _lftda function.
Definition: ECLDspEmulator.h:30
Belle2::ECL::ECLDspUtilities::pedestalFit
static ECLPedestalFit pedestalFit(std::vector< int > adc)
Fit pedestal part of the signal waveform (first 16 samples) This method will fit the first 16 samples...
Definition: ECLDspUtilities.cc:344
Belle2::DBObjPtr
Class for accessing objects in the database.
Definition: DBObjPtr.h:31
Belle2::ECLDspData
This object contains ECL DSP coefs – electromagnetic calorimeter digital signal processing coefficien...
Definition: ECLDspData.h:44
Belle2::ECL::ECLChannelMapper
This class provides access to ECL channel map that is either a) Loaded from the database (see ecl/dbo...
Definition: ECLChannelMapper.h:36
Belle2::ECL::ECLPedestalFit
This struct is returned by the pedestalFit method that fits the first 16 samples of the waveform (ped...
Definition: ECLDspUtilities.h:42
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::ECL::ECLDspUtilities::pedfit_fg31
static float pedfit_fg31[768]
DSP coefficients used to determine amplitude in pedestalFit.
Definition: ECLDspUtilities.h:106
Belle2::ECL::ECLChannelMapper::initFromDB
bool initFromDB()
Initialize channel mapper from the conditions database.
Definition: ECLChannelMapper.cc:105
Belle2::ECL::ECLDspUtilities::pedestal_fit_initialized
static int pedestal_fit_initialized
Flag indicating whether arrays fg31,fg32 are filled.
Definition: ECLDspUtilities.h:104
Belle2::ECL::ECLDspUtilities::readEclDsp
static ECLDspData * readEclDsp(const char *raw_file, int boardNumber)
Convert ECLDspData from *.dat file to Root object.
Definition: ECLDspUtilities.cc:48
Belle2::FileSystem::findFile
static std::string findFile(const std::string &path, bool silent=false)
Search for given file or directory in local or central release directory, and return absolute path if...
Definition: FileSystem.cc:147