BmnRoot
Loading...
Searching...
No Matches
BmnAdcProcessor.cxx
Go to the documentation of this file.
1#include "BmnAdcProcessor.h"
2
3#include "P4_F32vec4.h"
4#include "smmintrin.h"
5
6#include <BmnMath.h>
7
8BmnAdcProcessor::BmnAdcProcessor(uint32_t period, uint32_t run, TString det, Int_t nCh, const Int_t nSmpl)
9 : fVerbose(0)
10 , drawCnt(0)
11 , drawCnt2d(0)
12 , thrDif(0)
13 , FinalThr(0)
14 , fNIters(0)
15 , cmodcut(0)
16 , thrMax(0)
17 , fNSerials(0)
18 , fNSamples(nSmpl)
19 , fNChannels(nCh)
20 , fEvForPedestals(N_EV_FOR_PEDESTALS)
21 , fPedDat(nullptr)
22 , fSetup(kBMNSETUP)
23 , fApplyThreshold(true)
24 , fSigProf(nullptr)
25 , fNoisyChannels(nullptr)
26 , fAdc(nullptr)
27 , fPedRms(nullptr)
28 , fPedVal(nullptr)
29 , fPedValTemp(nullptr)
30 , fCMode(nullptr)
31 , fCMode0(nullptr)
32 , fSMode(nullptr)
33 , fSMode0(nullptr)
34 , fNoisyChipChannels(nullptr)
35 , fNvals(nullptr)
36 , fNvalsCMod(nullptr)
37 , fNvalsADC(nullptr)
38 , fPedCMod(nullptr)
39 , fPedCMod2(nullptr)
40 , fSumRmsV(nullptr)
41 , fDetName(det)
42 , fPeriod(period)
43 , fRun(run)
44{
46 ?
47 // #if CMAKE_BUILD_TYPE == Debug
48 // &BmnAdcProcessor::PrecalcEventMods : &BmnAdcProcessor::PrecalcEventModsOld;
49 // printf("\n\nDebug!!!\n\n");
50 // #else
51 // &BmnAdcProcessor::PrecalcEventMods_simd : &BmnAdcProcessor::PrecalcEventModsOld;
52 // printf("\n\nRelease!!!\n\n");
53 // #endif
54
55#ifdef BUILD_DEBUG
58 if (fVerbose)
59 printf("\n\nDebug!!!\n");
60#else
63 if (fVerbose)
64 printf("\n\nRelease!!!\n");
65#endif
66}
67
69{
70 auto serIter = fSerMap.find(serial);
71 if (serIter == fSerMap.end()) {
72 if (fVerbose > 1)
73 printf("adding %08X to the map, iSer = %d\n", serial, fNSerials);
74 fSerMap.insert(pair<UInt_t, Int_t>(serial, fNSerials++));
75 fAdcSerials.push_back(serial);
76 }
77}
78
80{
81 fPedVal = new Float_t**[fNSerials];
82 // fPedVal = new Double_t**[fNSerials];
83 fPedValTemp = new Double_t**[fNSerials];
84 // fNvals = new UInt_t*[fNSerials];
85 fNvals = new Float_t*[fNSerials];
86 fNvalsCMod = new UInt_t**[fNSerials];
87 fNvalsADC = new UInt_t**[fNSerials];
88 fPedRms = new Double_t**[fNSerials];
89 fPedCMod = new Double_t**[fNSerials];
90 fPedCMod2 = new Double_t**[fNSerials];
91 fSumRmsV = new Double_t*[fNSerials];
92 fNoisyChipChannels = new Bool_t**[fNSerials];
93 fCMode = new Float_t*[fNSerials];
94 fCMode0 = new Float_t*[fNSerials];
95 fSMode = new Float_t*[fNSerials];
96 fSMode0 = new Float_t*[fNSerials];
97 fAdcProfiles = new UInt_t**[fNSerials];
98 fAdc = new Float_t**[fNSerials];
99 // fAdc = new Double_t**[fNSerials];
100 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
101 // fPedVal[iCr] = new Double_t*[fNChannels];
102 fPedVal[iCr] = new Float_t*[fNChannels];
103 fPedValTemp[iCr] = new Double_t*[fNChannels];
104 // fNvals[iCr] = new UInt_t[fNChannels];
105 fNvals[iCr] = new Float_t[fNChannels];
106 memset(fNvals[iCr], 0, sizeof(Float_t) * fNChannels);
107 fNvalsCMod[iCr] = new UInt_t*[fNChannels];
108 fNvalsADC[iCr] = new UInt_t*[fNChannels];
109 fPedRms[iCr] = new Double_t*[fNChannels];
110 fPedCMod[iCr] = new Double_t*[fNChannels];
111 fPedCMod2[iCr] = new Double_t*[fNChannels];
112 fSumRmsV[iCr] = new Double_t[fNChannels];
113 memset(fSumRmsV[iCr], 0, sizeof(Double_t) * fNChannels);
114 fNoisyChipChannels[iCr] = new Bool_t*[fNChannels];
115 fAdcProfiles[iCr] = new UInt_t*[fNChannels];
116 // fAdc[iCr] = new Double_t*[fNChannels];
117 fAdc[iCr] = new Float_t*[fNChannels];
118 fCMode[iCr] = new Float_t[fNChannels];
119 fCMode0[iCr] = new Float_t[fNChannels];
120 fSMode[iCr] = new Float_t[fNChannels];
121 fSMode0[iCr] = new Float_t[fNChannels];
122 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
123 // fPedVal[iCr][iCh] = new Double_t[fNSamples];
124 fPedVal[iCr][iCh] = new Float_t[fNSamples];
125 fPedValTemp[iCr][iCh] = new Double_t[fNSamples];
126 fNvalsADC[iCr][iCh] = new UInt_t[fNSamples];
127 memset(fNvalsADC[iCr][iCh], 0, sizeof(UInt_t) * fNSamples);
128 fNvalsCMod[iCr][iCh] = new UInt_t[fNSamples];
129 memset(fNvalsCMod[iCr][iCh], 0, sizeof(UInt_t) * fNSamples);
130 fPedRms[iCr][iCh] = new Double_t[fNSamples];
131 fPedCMod[iCr][iCh] = new Double_t[fNSamples];
132 fPedCMod2[iCr][iCh] = new Double_t[fNSamples];
133 fNoisyChipChannels[iCr][iCh] = new Bool_t[fNSamples];
134 fAdcProfiles[iCr][iCh] = new UInt_t[fNSamples];
135 // fAdc[iCr][iCh] = new Double_t[fNSamples];
136 fAdc[iCr][iCh] = new Float_t[fNSamples];
137 fCMode[iCr][iCh] = 0.0;
138 fCMode0[iCr][iCh] = 0.0;
139 fSMode[iCr][iCh] = 0.0;
140 fSMode0[iCr][iCh] = 0.0;
141 for (Int_t iSmpl = 0; iSmpl < fNSamples; ++iSmpl) {
142 fPedVal[iCr][iCh][iSmpl] = 0.0;
143 fPedValTemp[iCr][iCh][iSmpl] = 0.0;
144 fPedRms[iCr][iCh][iSmpl] = 0.0;
145 fPedCMod[iCr][iCh][iSmpl] = 0.0;
146 fPedCMod2[iCr][iCh][iSmpl] = 0.0;
147 fNoisyChipChannels[iCr][iCh][iSmpl] = kFALSE;
148 fAdcProfiles[iCr][iCh][iSmpl] = 0;
149 fAdc[iCr][iCh][iSmpl] = 0.0;
150 }
151 }
152 }
153
154 fPedDat = new Double_t***[fNSerials];
155 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
156 fPedDat[iCr] = new Double_t**[N_EV_FOR_PEDESTALS];
157 for (UInt_t iEv = 0; iEv < N_EV_FOR_PEDESTALS; ++iEv) {
158 fPedDat[iCr][iEv] = new Double_t*[fNChannels];
159 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
160 fPedDat[iCr][iEv][iCh] = new Double_t[fNSamples];
161 for (Int_t iSmpl = 0; iSmpl < fNSamples; ++iSmpl)
162 fPedDat[iCr][iEv][iCh][iSmpl] = 0.0;
163 }
164 }
165 }
166 // Int_t nevents = N_EV_FOR_PEDESTALS - 2;
167 // for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
168 // vector<TH1*> hv;
169 // vector<TH1*> hcm;
170 // vector<TH1*> hsm;
171 // for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
172 // TString hname = TString(Form("%08X:%02d pedestals SM ev", fAdcSerials[iCr], iCh));
173 // TH1* h = new TH2F(hname, hname,
174 // N_EV_FOR_PEDESTALS, 0, N_EV_FOR_PEDESTALS,
175 // fNSamples, 0, fNSamples);
176 // h->GetXaxis()->SetTitle("Event #");
177 // h->GetYaxis()->SetTitle("Sample(channel) #");
178 // hv.push_back(h);
179 //
180 // hname = TString(Form("%08X:%02d cmods SM ev", fAdcSerials[iCr], iCh));
181 // TH1* hc = new TH1F(hname, hname, nevents, 0, nevents);
182 // hc->GetXaxis()->SetTitle("Event #");
183 // hcm.push_back(hc);
184 //
185 // hname = TString(Form("%08X:%02d smods SM ev", fAdcSerials[iCr], iCh));
186 // TH1* hs = new TH1F(hname, hname, nevents, 0, nevents);
187 // hs->GetXaxis()->SetTitle("Event #");
188 // hsm.push_back(hs);
189 // }
190 // hPedLine.push_back(hv);
191 // hCMode.push_back(hcm);
192 // hSMode.push_back(hsm);
193 // }
194}
195
197{
198 for (Int_t iCr = 0; iCr < fNSerials; iCr++) {
199 for (Int_t iEv = 0; iEv < N_EV_FOR_PEDESTALS; iEv++) {
200 for (Int_t iCh = 0; iCh < fNChannels; iCh++)
201 delete[] fPedDat[iCr][iEv][iCh];
202 delete[] fPedDat[iCr][iEv];
203 }
204 delete[] fPedDat[iCr];
205 }
206 delete[] fPedDat;
207 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
208 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
209 delete[] fPedVal[iCr][iCh];
210 delete[] fPedValTemp[iCr][iCh];
211 delete[] fNvalsCMod[iCr][iCh];
212 delete[] fNvalsADC[iCr][iCh];
213 delete[] fPedRms[iCr][iCh];
214 delete[] fPedCMod[iCr][iCh];
215 delete[] fPedCMod2[iCr][iCh];
216 delete[] fAdcProfiles[iCr][iCh];
217 delete[] fAdc[iCr][iCh];
218 delete[] fNoisyChipChannels[iCr][iCh];
219 }
220 delete[] fPedVal[iCr];
221 delete[] fPedValTemp[iCr];
222 delete[] fNvalsADC[iCr];
223 delete[] fNvals[iCr];
224 delete[] fNvalsCMod[iCr];
225 delete[] fPedRms[iCr];
226 delete[] fPedCMod[iCr];
227 delete[] fPedCMod2[iCr];
228 delete[] fSumRmsV[iCr];
229 delete[] fAdcProfiles[iCr];
230 delete[] fAdc[iCr];
231 delete[] fCMode[iCr];
232 delete[] fCMode0[iCr];
233 delete[] fSMode[iCr];
234 delete[] fSMode0[iCr];
235 delete[] fNoisyChipChannels[iCr];
236 }
237 delete[] fPedVal;
238 delete[] fPedValTemp;
239 delete[] fNvalsADC;
240 delete[] fNvals;
241 delete[] fNvalsCMod;
242 delete[] fPedRms;
243 delete[] fPedCMod;
244 delete[] fPedCMod2;
245 delete[] fSumRmsV;
246 delete[] fAdcProfiles;
247 delete[] fAdc;
248 delete[] fCMode;
249 delete[] fCMode0;
250 delete[] fSMode;
251 delete[] fSMode0;
252 delete[] fNoisyChipChannels;
253 /*canStrip->cd(1);
254 h->Draw("colz");
255 canStrip->cd(2);
256 hp->Draw("colz");
257 canStrip->cd(3);
258 hcms->Draw("colz");
259 canStrip->cd(4);
260 hscms_adc->Draw("colz");
261 canStrip->cd(5);
262 hcms1p->Draw("colz");
263 canStrip->cd(6);
264 hscms1p_adc->Draw("colz");
265 canStrip->cd(7);
266 hcms1->Draw("colz");
267 canStrip->cd(8);
268 hscms1_adc->Draw("colz");
269 canStrip->SaveAs("can-cms.png");*/
270 // for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
271 // delete hSModeSi[iCr];
272 // delete hCModeSi[iCr];
273 // for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
274 // delete hPedLineSi[iCr][iCh];
275 // }
276 // }
277 LOGF(info, "%10s: Real time %4.3f s, CPU time %4.3f s ADC processor. Precalc", fDetName.Data(), workTime_real_p,
279 LOGF(info, "%10s: Real time %4.3f s, CPU time %4.3f s ADC processor", fDetName.Data(), workTime_real, workTime_cpu);
280}
281
282void BmnAdcProcessor::CopyData2PedMap(TClonesArray* adc, UInt_t ev)
283{
284 for (Int_t iAdc = 0; iAdc < adc->GetEntriesFast(); ++iAdc) {
285 BmnADCDigit* adcDig = (BmnADCDigit*)adc->At(iAdc);
286 // printf("GEM ser 0x%08X, ev %d\n", adcDig->GetSerial(), ev);
287 auto serIter = fSerMap.find(adcDig->GetSerial());
288 if (serIter == fSerMap.end()) {
289 // printf("Serial %08X not found in the map\n", ser);
290 continue;
291 }
292 Int_t iSer = serIter->second;
293 // printf("GEM ser = 0x%08x, iSer = %02d, ev= %05d, ch = %d\n", adcDig->GetSerial(), iSer, ev,
294 // adcDig->GetChannel());
295 for (Int_t iSmpl = 0; iSmpl < fNSamples; ++iSmpl) {
296 if ((fRun > GetBoundaryRun(ADC32_N_SAMPLES)) || (fPeriod >= 7)) {
297 fPedDat[iSer][ev][adcDig->GetChannel()][iSmpl] = (Double_t)(adcDig->GetShortValue())[iSmpl] / 16.0;
298 } else
299 fPedDat[iSer][ev][adcDig->GetChannel()][iSmpl] = (Double_t)(adcDig->GetUShortValue())[iSmpl] / 16.0;
300 // printf("adc = %i pedData = %f\n", (adcDig->GetShortValue())[iSmpl],
301 // pedData[iSer][ev][adcDig->GetChannel()][iSmpl]);
302 }
303 }
304}
305
307{
308 TString docName = TString(Form("%s-mods-ev-%d.pdf", fDetName.Data(), drawCnt));
309 const UInt_t PAD_WIDTH_SIL = 1920; // 8192;
310 const UInt_t PAD_HEIGHT_SIL = 1080;
311 // TString docName = "sil-ped-cms.pdf";
312 TCanvas* c = new TCanvas("can cms", "can", PAD_WIDTH_SIL, 2 * PAD_HEIGHT_SIL);
313 c->Divide(1, 2);
314 c->Print(docName + "[");
315 for (int iCr = 0; iCr < fNSerials; iCr++)
316 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
317 // auto *p = c->cd(iCr * fNChannels + iCh + 1);
318 c->Clear("D");
319 // hPedLine[iCr][iCh]->Draw("colz");
320 c->cd(1);
321 hSMode[iCr][iCh]->Draw("colz");
322 c->cd(2);
323 hCMode[iCr][iCh]->Draw("colz");
324 c->Print(docName);
325 // c->SaveAs(Form("%s.png", hPedLine[iCr][iCh]->GetName()));
326 }
327 c->Print(docName + "]");
328 // c->SaveAs();
329 drawCnt++;
330}
331
333{
334 TString docName = TString(Form("%s-run-%d-mods-%d.pdf", fDetName.Data(), fRun, drawCnt2d));
335 const UInt_t PAD_WIDTH_SIL = 1920; // 8192;
336 const UInt_t PAD_HEIGHT_SIL = 1080;
337 gStyle->SetOptStat(0);
338 // TColor::InvertPalette();
339 gStyle->SetPalette(kDeepSea);
340 TCanvas* c = new TCanvas("can cms", "can", PAD_WIDTH_SIL, 2 * PAD_HEIGHT_SIL);
341 c->Divide(1, 3);
342 c->Print(docName + "[");
343 for (int iCr = 0; iCr < fNSerials; iCr++) {
344 c->Clear("D");
345 c->cd(1);
346 hPedSi[iCr]->Draw("box");
347 c->cd(2);
348 hCModeSi[iCr]->Draw("colz");
349 c->cd(3);
350 hSModeSi[iCr]->Draw("colz");
351 c->Print(docName);
352 }
353 // c->Print(docName + "]");
354 delete c;
355 c = new TCanvas("can pedestals", "can", PAD_WIDTH_SIL, PAD_HEIGHT_SIL);
356 // c->Print(docName + "[");
357 c->Clear("D");
358 // c->Divide(1, 2);
359 for (int iCr = 0; iCr < fNSerials; iCr++)
360 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
361 c->Clear("D");
362 c->cd(1);
363 hPedLineSi[iCr][iCh]->Draw("colz");
364 c->Print(docName);
365 // c->SaveAs(Form("%s.png", hPedLine[iCr][iCh]->GetName()));
366 }
367 c->Print(docName + "]");
368 delete c;
369 drawCnt2d++;
370}
371
373{
374 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
375 hPedSi[iCr]->Reset();
376 hSModeSi[iCr]->Reset();
377 hCModeSi[iCr]->Reset();
378 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
379 hPedLine[iCr][iCh]->Reset();
380 hSMode[iCr][iCh]->Reset();
381 hCMode[iCr][iCh]->Reset();
382 }
383 }
384}
385
387{
388 const UShort_t nSmpl = fNSamples;
389
390 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
391 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
392 for (Int_t iSmpl = 0; iSmpl < fNSamples; ++iSmpl) {
393 fPedVal[iCr][iCh][iSmpl] = 0.0;
394 fPedRms[iCr][iCh][iSmpl] = 0.0;
395 fAdcProfiles[iCr][iCh][iSmpl] = 0;
396 }
397 }
398 // cout << fDetName << " pedestals calculation..." << endl;
399 for (Int_t iEv = 0; iEv < N_EV_FOR_PEDESTALS; ++iEv) {
400 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
401 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
402 Double_t signals[nSmpl];
403 for (Int_t i = 0; i < nSmpl; ++i)
404 signals[i] = 0.0;
405 Int_t nOk = 0;
406 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl) {
407 if (fPedDat[iCr][iEv][iCh][iSmpl] == 0)
408 continue;
409 signals[iSmpl] = fPedDat[iCr][iEv][iCh][iSmpl];
410 nOk++;
411 }
412 Double_t CMS = CalcCMS(signals, nOk);
413 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl) {
414 fPedVal[iCr][iCh][iSmpl] += ((signals[iSmpl] - CMS) / N_EV_FOR_PEDESTALS);
415 }
416 }
417 }
418
419 // cout << fDetName << " RMS calculation..." << endl;
420 for (Int_t iEv = 0; iEv < N_EV_FOR_PEDESTALS; ++iEv) {
421 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
422 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
423 Double_t signals[nSmpl];
424 for (Int_t i = 0; i < nSmpl; ++i)
425 signals[i] = 0.0;
426 Int_t nOk = 0;
427 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl) {
428 if (fPedDat[iCr][iEv][iCh][iSmpl] == 0)
429 continue;
430 signals[iSmpl] = fPedDat[iCr][iEv][iCh][iSmpl];
431 nOk++;
432 }
433 Double_t CMS = CalcCMS(signals, nOk);
434 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl) {
435 Float_t ped = fPedVal[iCr][iCh][iSmpl];
436 fPedRms[iCr][iCh][iSmpl] += ((signals[iSmpl] - CMS - ped) * (signals[iSmpl] - CMS - ped));
437 }
438 }
439 }
440
441 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
442 for (Int_t iCh = 0; iCh < fNChannels; ++iCh)
443 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl)
444 fPedRms[iCr][iCh][iSmpl] = Sqrt(fPedRms[iCr][iCh][iSmpl] / N_EV_FOR_PEDESTALS);
445
446 ofstream pedFile(Form("%s/input/%s_pedestals_%d.txt", getenv("VMCWORKDIR"), fDetName.Data(), fRun));
447 pedFile << "Serial\tCh_id\tPed\tRMS" << endl;
448 pedFile << "============================================" << endl;
449 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
450 for (Int_t iCh = 0; iCh < fNChannels; ++iCh)
451 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl)
452 pedFile << hex << fAdcSerials[iCr] << dec << "\t" << iCh * nSmpl + iSmpl << "\t"
453 << fPedVal[iCr][iCh][iSmpl] << "\t" << fPedRms[iCr][iCh][iSmpl] << endl;
454 pedFile.close();
455
456 return kBMNSUCCESS;
457}
458
460{
461 // if (hSModeSi.size() == 0) {
462 // for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
463 // vector<TH1*> hv;
464 // vector<TH1*> hcm;
465 // vector<TH1*> hsm;
466 // for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
467 // // printf("Creating icr %2d ich %2d %s\n", iCr, iCh, Form("%08X:%02d
468 // pedestal line MK", fSerials[iCr], iCh)); TString hname = TString(Form("%08X:%02d pedestals SM",
469 // fAdcSerials[iCr], iCh)); TH1* h = new TH2F(hname, hname,
470 // 500, 0, 500,
471 // fNSamples, 0, fNSamples);
472 // h->GetXaxis()->SetTitle("Event #");
473 // h->GetYaxis()->SetTitle("Sample(channel) #");
474 // h->SetDirectory(0);
475 // hv.push_back(h);
476 // }
477 // hPedLineSi.push_back(hv);
478 // }
479 // const Int_t maxAdc = 8192;
480 // const Int_t MaxSig = 2300;
481 // const Int_t RngSig = 400;
482 // const Int_t StripSi = fNChannels * fNSamples;
483 // const Int_t StripSiLo = 9 * fNSamples;
484 // const Int_t StripSiHi = 14 * fNSamples;
485 // const Int_t StripSiBins = (StripSiHi - StripSiLo);
486 // for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
487 // TString hname = TString(Form("%08X pedestals SM", fAdcSerials[iCr]));
488 // // TH1* h = new TH2F(hname, hname, maxAdc, 0, maxAdc, MaxSig, 0, MaxSig);
489 // TH1* h = new TH2F(hname, hname, StripSiBins, StripSiLo, StripSiHi, 2 * RngSig, -RngSig, RngSig);
490 // // printf("maxAdc %04d max %04f peds\n", maxAdc, h->GetXaxis()->GetXmax());
491 // h->GetXaxis()->SetTitle("Channel #");
492 // h->GetYaxis()->SetTitle("Signal");
493 // h->SetDirectory(0);
494 // hPedSi.push_back(h);
495 // hname = TString(Form("%08X cmods SM", fAdcSerials[iCr]));
496 // TH1* hc = new TH2F(hname, hname, StripSiBins, StripSiLo, StripSiHi, 2 * RngSig, -RngSig, RngSig);
497 // // printf("maxAdc %04d max %04f cmode\n", maxAdc, hc->GetXaxis()->GetXmax());
498 // hc->GetXaxis()->SetTitle("Channel #");
499 // hc->GetYaxis()->SetTitle("Signal");
500 // hc->SetDirectory(0);
501 // hCModeSi.push_back(hc);
502 // hname = TString(Form("%08X smods SM", fAdcSerials[iCr]));
503 // TH1* hs = new TH2F(hname, hname, StripSiBins, StripSiLo, StripSiHi, 2 * RngSig, -RngSig, RngSig);
504 // hs->GetXaxis()->SetTitle("Channel #");
505 // hs->GetYaxis()->SetTitle("Signal");
506 // hs->SetDirectory(0);
507 // hSModeSi.push_back(hs);
508 // }
509 // }
510 if (fVerbose)
511 printf("%s %s started niter %d thrMax %4.2f\n", fDetName.Data(), __func__, fNIters, thrMax);
512 const UShort_t nSmpl = fNSamples;
513 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
514 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
515 // memset(fNvalsADC[iCr][iCh], 0, sizeof (UInt_t) * fNSamples);
516 for (Int_t iSmpl = 0; iSmpl < fNSamples; ++iSmpl) {
517 fPedVal[iCr][iCh][iSmpl] = 0.0;
518 fPedValTemp[iCr][iCh][iSmpl] = 0.0;
519 fPedRms[iCr][iCh][iSmpl] = 0.0;
520 fAdcProfiles[iCr][iCh][iSmpl] = 0;
521 fNvalsADC[iCr][iCh][iSmpl] = 0;
522 // fNoisyChipChannels[iCr][iCh][iSmpl] = kFALSE;
523 }
524 }
525 for (Int_t iEv = 0; iEv < N_EV_FOR_PEDESTALS; ++iEv) {
526 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
527 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
528 Int_t nOk = 0;
529 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl) {
530 if (fPedDat[iCr][iEv][iCh][iSmpl] == 0.0 || fNoisyChipChannels[iCr][iCh][iSmpl] == kTRUE)
531 continue;
532 fPedVal[iCr][iCh][iSmpl] += fPedDat[iCr][iEv][iCh][iSmpl]; // / N_EV_FOR_PEDESTALS);
533 fNvalsADC[iCr][iCh][iSmpl]++;
534 nOk++;
535 // static_cast<TH2*> (hPedLineSi[iCr][iCh])->Fill(iEv, iSmpl, fPedDat[iCr][iEv][iCh][iSmpl]);
536 }
537 }
538 }
539 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
540 for (Int_t iCh = 0; iCh < fNChannels; ++iCh)
541 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl)
542 if (fNvalsADC[iCr][iCh][iSmpl])
543 fPedVal[iCr][iCh][iSmpl] /= fNvalsADC[iCr][iCh][iSmpl];
544 // iteratively calculate pedestals and CMSs
545 // Double_t rmsthr = 200.0;
546 // Double_t rmsthrf = 200.0;
547 Double_t sumRms = thrMax / 3.5;
548 for (Int_t iter = 0; iter < fNIters; iter++) {
549 Double_t thr = thrMax - thrDif * iter;
550 if (thr < 0)
551 thr = 0;
552 if (fVerbose)
553 printf("iter %d thr %4.2f\n", iter, thr);
554 UInt_t nFiltered = 0;
555 // clear
556 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
557 memset(fSumRmsV[iCr], 0, sizeof(Double_t) * fNChannels);
558 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
559 memset(fNvalsCMod[iCr][iCh], 0, sizeof(UInt_t) * fNSamples);
560 memset(fNvalsADC[iCr][iCh], 0, sizeof(UInt_t) * fNSamples);
561 memset(fPedValTemp[iCr][iCh], 0, sizeof(Double_t) * fNSamples);
562 memset(fPedCMod[iCr][iCh], 0, sizeof(Double_t) * fNSamples);
563 memset(fPedCMod2[iCr][iCh], 0, sizeof(Double_t) * fNSamples);
564 }
565 }
566 for (Int_t iEv = 0; iEv < N_EV_FOR_PEDESTALS - 1; ++iEv) {
567 // clear
568 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
569 // memset(fNvals[iCr], 0, sizeof (UInt_t) * fNChannels);
570 memset(fNvals[iCr], 0, sizeof(Float_t) * fNChannels);
571 memset(fCMode[iCr], 0, sizeof(Float_t) * fNChannels);
572 memset(fSMode[iCr], 0, sizeof(Float_t) * fNChannels);
573 }
574 // Pedestals pre filtering
575 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
576 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
577 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl) {
578 // if (iCr == 3 && iCh == 8) {
579 // printf("iter %i iEv %i fpedDat %f noise %i\n", iter, iEv,
580 // fPedDat[iCr][iEv][iCh][iSmpl],
581 // fNoisyChipChannels[iCr][iCh][iSmpl]);
582 // }
583 if (fPedDat[iCr][iEv][iCh][iSmpl] == 0 || fNoisyChipChannels[iCr][iCh][iSmpl] == kTRUE)
584 continue;
585 Double_t Asig = TMath::Abs(fPedDat[iCr][iEv][iCh][iSmpl] - fPedVal[iCr][iCh][iSmpl]);
586 if (Asig < thr) {
587 fSMode[iCr][iCh] += fPedDat[iCr][iEv][iCh][iSmpl]; // CMS from current event
588 fCMode[iCr][iCh] += fPedVal[iCr][iCh][iSmpl]; // CMS over all pedestals
589 // fPedValTemp[iCr][iCh][iSmpl] += fPedDat[iCr][iEv][iCh][iSmpl]; //
590 // CMS from current event fNvalsADC[iCr][iCh][iSmpl]++;
591 fNvals[iCr][iCh]++;
592 }
593 }
594 }
595 // normalize
596 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
597 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
598 if (fNvals[iCr][iCh]) {
599 fSMode[iCr][iCh] /= fNvals[iCr][iCh];
600 fCMode[iCr][iCh] /= fNvals[iCr][iCh];
601 // hSMode[iCr][iCh]->SetBinContent(iEv, fSMode[iCr][iCh]);
602 // hCMode[iCr][iCh]->SetBinContent(iEv, fCMode[iCr][iCh]);
603 } else {
604 fSMode[iCr][iCh] = 0;
605 fCMode[iCr][iCh] = 0;
606 }
607 }
608 // Pedestals filtering
609 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
610 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
611 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl) {
612 Double_t Adc = fPedDat[iCr][iEv][iCh][iSmpl];
613 if ((Adc == 0) || (fNoisyChipChannels[iCr][iCh][iSmpl] == kTRUE))
614 continue;
615 Double_t sig = fPedDat[iCr][iEv][iCh][iSmpl] - fPedVal[iCr][iCh][iSmpl] + fCMode[iCr][iCh]
616 - fSMode[iCr][iCh];
617 Double_t Asig = TMath::Abs(sig);
618 if (Asig < thr) {
619 fPedValTemp[iCr][iCh][iSmpl] += fPedDat[iCr][iEv][iCh][iSmpl];
620 fNvalsADC[iCr][iCh][iSmpl]++;
621
622 Adc = fPedDat[iCr][iEv][iCh][iSmpl] - fSMode[iCr][iCh];
623 fPedCMod[iCr][iCh][iSmpl] += Adc;
624 fPedCMod2[iCr][iCh][iSmpl] += Adc * Adc;
625 fNvalsCMod[iCr][iCh][iSmpl]++;
626 nFiltered++;
627 }
628 // if (iter == fEndIter - 1) {
629 // Int_t ic = iCh * nSmpl + iSmpl;
630 // printf("iCh %4d iSmpl %4d ic %4d cmod %5f smod %5f\n", iCh,
631 // iSmpl, ic, fCMode[iCr][iCh], fSMode[iCr][iCh]);
632 // hCModeSi[iCr]->Fill(ic, fCMode[iCr][iCh]);
633 // hSModeSi[iCr]->Fill(ic, fSMode[iCr][iCh]);
634 // hPedSi[iCr]->Fill(ic, fPedVal[iCr][iCh][iSmpl]);
635 //}
636 }
637 }
638 } // event loop
639
640 sumRms = 0.0;
641 Int_t nrms = 0;
642
643 // hists fill
644 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
645 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
646 Int_t nvrms = 0;
647 fSumRmsV[iCr][iCh] = 0.0;
648 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl) {
649 if (fNoisyChipChannels[iCr][iCh][iSmpl] == kTRUE)
650 continue;
651 if (fNvalsCMod[iCr][iCh][iSmpl]) {
652 fPedCMod[iCr][iCh][iSmpl] /= fNvalsCMod[iCr][iCh][iSmpl];
653 fPedCMod2[iCr][iCh][iSmpl] = Sqrt(Abs(fPedCMod2[iCr][iCh][iSmpl] / fNvalsCMod[iCr][iCh][iSmpl]
654 - Sq(fPedCMod[iCr][iCh][iSmpl])));
655 sumRms += fPedCMod2[iCr][iCh][iSmpl];
656 fSumRmsV[iCr][iCh] += fPedCMod2[iCr][iCh][iSmpl];
657 nrms++;
658 nvrms++;
659 }
660 if (fNvalsADC[iCr][iCh][iSmpl])
661 fPedVal[iCr][iCh][iSmpl] = fPedValTemp[iCr][iCh][iSmpl] / fNvalsADC[iCr][iCh][iSmpl];
662 else
663 fPedVal[iCr][iCh][iSmpl] = 0.0;
664 // printf("fPedVal[%3d][%3d][%3d] = %4.2f\n", iCr, iCh,
665 // iSmpl,fPedVal[iCr][iCh][iSmpl]);
666 fNvalsADC[iCr][iCh][iSmpl] = 0;
667 // fCMode[iCr][iCh] += fPedVal[iCr][iCh][iSmpl]; // CMS over all pedestals
668 }
669 if (nvrms)
670 fSumRmsV[iCr][iCh] /= nvrms;
671 }
672 if (nrms > 0)
673 sumRms /= nrms;
674
675 // noise ch detection
676 if (fDetName.Contains("SiBT") || fDetName.Contains("SiProf"))
677 continue;
678 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
679 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
680 if ((iter == fNIters - 1) && (fVerbose > 1))
681 printf("cr %8X ich %2d fSumRmsV %4f sumRms %4f\n", fAdcSerials[iCr], iCh, fSumRmsV[iCr][iCh],
682 sumRms);
683 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl) {
684 if (fNoisyChipChannels[iCr][iCh][iSmpl] == kTRUE)
685 continue;
686 if ((iter == fNIters - 1) && (fVerbose > 1))
687 printf("cr %8X ich %2d iSmpl %2d fPedCMod2 %4f\n", fAdcSerials[iCr], iCh, iSmpl,
688 fPedCMod2[iCr][iCh][iSmpl]);
689 if ((fPedCMod2[iCr][iCh][iSmpl]
690 > 4 * sumRms) /* || (fPedCMod2[iCr][iCh][iSmpl] > 5 * fSumRmsV[iCr][iCh])*/)
691 {
692 fNoisyChipChannels[iCr][iCh][iSmpl] = kTRUE;
693 if (fVerbose)
694 printf("new noisy ch on cr %i %08X ch %i smpl %i by signal %4.2f\n", iCr,
695 fAdcSerials[iCr], iCh, iSmpl, fPedCMod2[iCr][iCh][iSmpl]);
696 }
697 }
698 }
699 } // iter loop
700 return kBMNSUCCESS;
701}
702
704{
705 TStopwatch timer;
706 // Double_t rtime;
707 // Double_t ctime;
708 timer.Start();
709 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
710 memset(fNvals[iCr], 0, sizeof(Float_t) * fNChannels);
711 memset(fCMode[iCr], 0, sizeof(Float_t) * fNChannels);
712 memset(fSMode[iCr], 0, sizeof(Float_t) * fNChannels);
713 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
714 memset(fAdc[iCr][iCh], 0, sizeof(Float_t) * fNSamples);
715 }
716 }
717 timer.Stop();
718 // rtime = timer.RealTime();
719 // ctime = timer.CpuTime();
720 // printf("\nReal time %f s, CPU time %f s clear\n", rtime, ctime);
721 timer.Start();
722 for (Int_t iAdc = 0; iAdc < adc->GetEntriesFast(); ++iAdc) {
723 BmnADCDigit* adcDig = (BmnADCDigit*)adc->At(iAdc);
724 UInt_t iCh = adcDig->GetChannel();
725 UInt_t ser = adcDig->GetSerial();
726 // printf("Serial %08X \n", ser);
727 auto serIter = fSerMap.find(ser);
728 // printf("iter %08X end %08X\n", serIter, fSerMap.end());
729 if (serIter == fSerMap.end()) {
730 // printf("Serial %08X not found in the map\n", ser);
731 continue;
732 }
733 Int_t iCr = serIter->second;
734 for (Int_t iSmpl = 0; iSmpl < fNSamples; iSmpl++) {
735 if ((fNoisyChipChannels[iCr][iCh][iSmpl] == kTRUE))
736 continue;
737 Double_t val = static_cast<Double_t>(adcDig->GetUShortValue()[iSmpl] / 16);
738 fAdc[iCr][iCh][iSmpl] = val;
739
740 Double_t Sig = fAdc[iCr][iCh][iSmpl] - fPedVal[iCr][iCh][iSmpl];
741 Double_t Asig = TMath::Abs(Sig);
742 // printf("adc %6f ped %6f\n", fAdc[iCr][iCh][iSmpl], fPedVal[iCr][iCh][iSmpl]);
743
744 if ((Asig < thrMax)) {
745 // printf("adc %6f < thrMax %6f\n", fAdc[iCr][iCh][iSmpl], thrMax);
746 fSMode[iCr][iCh] += fAdc[iCr][iCh][iSmpl];
747 fCMode[iCr][iCh] += fPedVal[iCr][iCh][iSmpl];
748 fNvals[iCr][iCh]++;
749 }
750 }
751 }
752 timer.Stop();
753 // rtime = timer.RealTime();
754 // ctime = timer.CpuTime();
755 // printf("\nReal time %f s, CPU time %f s fill array\n", rtime, ctime);
756}
757
759{
760 TStopwatch timer;
761 // Double_t rtime;
762 // Double_t ctime;
763 timer.Start();
764 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
765 // memset(fNvals[iCr], 0, sizeof (UInt_t) * fNChannels);
766 memset(fNvals[iCr], 0, sizeof(Float_t) * fNChannels);
767 memset(fCMode[iCr], 0, sizeof(Float_t) * fNChannels);
768 memset(fSMode[iCr], 0, sizeof(Float_t) * fNChannels);
769 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
770 memset(fAdc[iCr][iCh], 0, sizeof(Float_t) * fNSamples);
771 }
772 }
773 timer.Stop();
774 // rtime = timer.RealTime();
775 // ctime = timer.CpuTime();
776 // printf("\nReal time %f s, CPU time %f s clear\n", rtime, ctime);
777 timer.Start();
778 for (Int_t iAdc = 0; iAdc < adc->GetEntriesFast(); ++iAdc) {
779 BmnADCDigit* adcDig = (BmnADCDigit*)adc->At(iAdc);
780 UInt_t iCh = adcDig->GetChannel();
781 UInt_t ser = adcDig->GetSerial();
782 // printf("Serial %08X channel %d\n", ser, iCh);
783 auto serIter = fSerMap.find(ser);
784 if (serIter == fSerMap.end()) {
785 // printf("Serial %08X not found in the map\n", ser);
786 continue;
787 }
788 Int_t iCr = serIter->second;
789 for (Int_t iSmpl = 0; iSmpl < fNSamples; iSmpl++) {
790 if ((fNoisyChipChannels[iCr][iCh][iSmpl] == kTRUE))
791 continue;
792 Double_t val = static_cast<Double_t>(adcDig->GetShortValue()[iSmpl] / 16);
793 fAdc[iCr][iCh][iSmpl] = val;
794
795 Double_t Sig = fAdc[iCr][iCh][iSmpl] - fPedVal[iCr][iCh][iSmpl];
796 Double_t Asig = TMath::Abs(Sig);
797 // printf("adc %6f ped %6f\n", fAdc[iCr][iCh][iSmpl], fPedVal[iCr][iCh][iSmpl]);
798
799 if ((Asig < thrMax)) {
800 // printf("adc %6f < thrMax %6f\n", fAdc[iCr][iCh][iSmpl], thrMax);
801 fSMode[iCr][iCh] += fAdc[iCr][iCh][iSmpl];
802 fCMode[iCr][iCh] += fPedVal[iCr][iCh][iSmpl];
803 fNvals[iCr][iCh]++;
804 }
805 }
806 }
807 timer.Stop();
808 workTime_cpu_p += timer.RealTime();
809 workTime_real_p += timer.CpuTime();
810}
811
812// simd function
813
815{
816 const Int_t VecLenSIMD = 4;
817 TStopwatch timer;
818 timer.Start();
819 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
820 // memset(fNvals[iCr], 0, sizeof (UInt_t) * fNChannels);
821 memset(fNvals[iCr], 0, sizeof(Float_t) * fNChannels);
822 memset(fCMode[iCr], 0, sizeof(Float_t) * fNChannels);
823 memset(fSMode[iCr], 0, sizeof(Float_t) * fNChannels);
824 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
825 memset(fAdc[iCr][iCh], 0, sizeof(Float_t) * fNSamples);
826 }
827 }
828 // timer.Stop();
829 // rtime = timer.RealTime();
830 // ctime = timer.CpuTime();
831 // // printf("\nReal time %f s, CPU time %f s clear\n", rtime, ctime);
832 // timer.Start();
833
834 // fvec * fNvals_vec, * fSMode_vec, * fCMode_vec;
835 fvec *fPedVal_vec, *fAdc_vec, *fNoisyChipChannels_vec;
836 Float_t*** fNoisy_float;
837
838 fNoisy_float = new Float_t**[fNSerials];
839 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
840 fNoisy_float[iCr] = new Float_t*[fNChannels];
841 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
842 fNoisy_float[iCr][iCh] = new Float_t[fNSamples];
843 for (Int_t iSmpl = 0; iSmpl < fNSamples; ++iSmpl) {
844 // fNoisy_float[iCr][iCh][iSmpl] = (Float_t)fNoisyChipChannels[iCr][iCh][iSmpl];
845 if (fNoisyChipChannels[iCr][iCh][iSmpl])
846 fNoisy_float[iCr][iCh][iSmpl] = 1.0;
847 else
848 fNoisy_float[iCr][iCh][iSmpl] = 0.0;
849 }
850 }
851 }
852 __m128 outThresh = _mm_set1_ps(10000); // some number definitely higher than threshold
853 for (Int_t iAdc = 0; iAdc < adc->GetEntriesFast(); ++iAdc) {
854 BmnADCDigit* adcDig = (BmnADCDigit*)adc->At(iAdc);
855 UInt_t iCh = adcDig->GetChannel();
856 UInt_t ser = adcDig->GetSerial();
857 // printf("Serial %08X \n", ser);
858 auto serIter = fSerMap.find(ser);
859 if (serIter == fSerMap.end()) {
860 // printf("Serial %08X not found in the map\n", ser);
861 continue;
862 }
863
864 Int_t iCr = serIter->second;
865 // printf("ser %08X iCr %d\n", ser, iCr);
866 for (Int_t iSmpl = 0; iSmpl < fNSamples; iSmpl++) {
867 if ((fNoisyChipChannels[iCr][iCh][iSmpl] == kTRUE))
868 continue;
869 Float_t val = static_cast<Float_t>(adcDig->GetShortValue()[iSmpl] / 16);
870 fAdc[iCr][iCh][iSmpl] = val;
871 }
872
873 fPedVal_vec = (fvec*)fPedVal[iCr][iCh];
874 fAdc_vec = (fvec*)fAdc[iCr][iCh];
875 fNoisyChipChannels_vec = (fvec*)fNoisy_float[iCr][iCh];
876 fvec fSMode_sum = 0.0;
877 fvec fCMode_sum = 0.0;
878 fvec Nval = 0.0;
879 Float_t sum1;
880 Float_t sum2;
881 Float_t sum3;
882
883 for (Int_t iSmpl = 0; iSmpl < fNSamples / VecLenSIMD; iSmpl++) {
884 // if ((fNoisyChipChannels[iCr][iCh][iSmpl] == kTRUE)) continue;
885 // fvec sig = if3(fvec(fNoisyChipChannels_vec[iSmpl] == 1.0), 10000, fAdc_vec[iSmpl]
886 // - fPedVal_vec[iSmpl]);
887 fvec sig = _mm_blendv_ps(fAdc_vec[iSmpl] - fPedVal_vec[iSmpl], outThresh,
888 fvec(fNoisyChipChannels_vec[iSmpl] == 1.0));
889 // cout << "result of the comparison: " << sig << endl;
890 fvec Asig = fabs(sig);
891 // cout << "absolute value: " << Asig << endl;
892 fvec thrMax_vec = thrMax;
893 // Nval = if3(fvec(Asig < thrMax_vec), Nval + 1.0, Nval + 0.0);
894 Nval = _mm_blendv_ps(Nval, Nval + 1.0, fvec(Asig < thrMax_vec));
895 fSMode_sum = _mm_blendv_ps(fSMode_sum, fSMode_sum + fAdc_vec[iSmpl], fvec(Asig < thrMax_vec));
896 fCMode_sum = _mm_blendv_ps(fCMode_sum, fCMode_sum + fPedVal_vec[iSmpl], fvec(Asig < thrMax_vec));
897 }
898 fSMode_sum = _mm_dp_ps(fSMode_sum, _f32vec4_one, 0xFF);
899 _mm_store_ss(&sum1, fSMode_sum);
900 fSMode[iCr][iCh] += sum1;
901 fCMode_sum = _mm_dp_ps(fCMode_sum, _f32vec4_one, 0xFF);
902 _mm_store_ss(&sum2, fCMode_sum);
903 fCMode[iCr][iCh] += sum2;
904 Nval = _mm_dp_ps(Nval, _f32vec4_one, 0xFF);
905 _mm_store_ss(&sum3, Nval);
906 fNvals[iCr][iCh] += sum3;
907 }
908 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
909 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
910 delete[] fNoisy_float[iCr][iCh];
911 }
912 delete[] fNoisy_float[iCr];
913 }
914 delete[] fNoisy_float;
915
916 timer.Stop();
917 workTime_cpu_p += timer.RealTime();
918 workTime_real_p += timer.CpuTime();
919}
920
922{
923 TStopwatch timer;
924 timer.Start();
925 // normalize
926 Int_t fNvalsMin = 0;
927 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
928 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
929 if (fNvals[iCr][iCh] > fNvalsMin) {
930 fSMode[iCr][iCh] /= fNvals[iCr][iCh];
931 fCMode[iCr][iCh] /= fNvals[iCr][iCh];
932 } else {
933 fSMode[iCr][iCh] = 0.0;
934 fCMode[iCr][iCh] = 0.0;
935 }
936 }
937 // filter out sigs to get mods
938 for (Int_t iter = 0; iter < fNIters; ++iter) {
939 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
940 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
941 fSMode0[iCr][iCh] = 0.0;
942 fCMode0[iCr][iCh] = 0.0;
943 fNvals[iCr][iCh] = 0;
944 }
945 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
946 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
947 Double_t cs = fCMode[iCr][iCh] - fSMode[iCr][iCh];
948 for (Int_t iSmpl = 0; iSmpl < fNSamples; ++iSmpl) {
949 if (fPedVal[iCr][iCh][iSmpl] == 0 || fNoisyChipChannels[iCr][iCh][iSmpl] == kTRUE)
950 continue;
951 Double_t sig = fAdc[iCr][iCh][iSmpl] - fPedVal[iCr][iCh][iSmpl] + cs;
952 Double_t Asig = TMath::Abs(sig);
953 Double_t thr = thrMax - iter * thrDif;
954 if (Asig < thr) {
955 fSMode0[iCr][iCh] += fAdc[iCr][iCh][iSmpl];
956 fCMode0[iCr][iCh] += fPedVal[iCr][iCh][iSmpl];
957 fNvals[iCr][iCh]++;
958 // rmsthrf += Sq(fPedDat[iCr][iEv][iCh][iSmpl] - fSigCMS[iCr][iCh]);
959 }
960 }
961 }
962
963 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
964 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
965 if (fNvals[iCr][iCh] > fNvalsMin) {
966 fSMode[iCr][iCh] = fSMode0[iCr][iCh] / fNvals[iCr][iCh];
967 fCMode[iCr][iCh] = fCMode0[iCr][iCh] / fNvals[iCr][iCh];
968 } else {
969 fSMode[iCr][iCh] = 0.0;
970 fCMode[iCr][iCh] = 0.0;
971 }
972 }
973 }
974 timer.Stop();
975 workTime_cpu += (Double_t)timer.CpuTime();
976 workTime_real += (Double_t)timer.RealTime();
977}
978
979// simd function
980
982{
983 const Int_t VecLenSIMD = 4;
984 TStopwatch timer;
985 timer.Start();
986
987 // normalize
988 fvec fNvalsMin_vec = 0;
989 __m128 outThresh = _mm_set1_ps(10000);
990 fvec *fNvals_vec, *fSMode_vec, *fCMode_vec;
991 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
992 fNvals_vec = (fvec*)fNvals[iCr];
993 fSMode_vec = (fvec*)fSMode[iCr];
994 fCMode_vec = (fvec*)fCMode[iCr];
995 for (Int_t iCh = 0; iCh < fNChannels / VecLenSIMD; iCh++) {
996 fSMode_vec[iCh] =
997 _mm_blendv_ps(_f32vec4_zero, fSMode_vec[iCh] / fNvals_vec[iCh], fvec(fNvals_vec[iCh] > fNvalsMin_vec));
998 fCMode_vec[iCh] =
999 _mm_blendv_ps(_f32vec4_zero, fCMode_vec[iCh] / fNvals_vec[iCh], fvec(fNvals_vec[iCh] > fNvalsMin_vec));
1000 }
1001 fNvals[iCr] = (Float_t*)fNvals_vec;
1002 fSMode[iCr] = (Float_t*)fSMode_vec;
1003 fCMode[iCr] = (Float_t*)fCMode_vec;
1004 }
1005 // filter out sigs to get mods
1006
1007 Float_t*** fNoisy_float;
1008 fNoisy_float = new Float_t**[fNSerials];
1009 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
1010 fNoisy_float[iCr] = new Float_t*[fNChannels];
1011 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
1012 fNoisy_float[iCr][iCh] = new Float_t[fNSamples];
1013 for (Int_t iSmpl = 0; iSmpl < fNSamples; ++iSmpl) {
1014 // fNoisy_float[iCr][iCh][iSmpl] = (Float_t)fNoisyChipChannels[iCr][iCh][iSmpl];
1015 if (fNoisyChipChannels[iCr][iCh][iSmpl] == kTRUE)
1016 fNoisy_float[iCr][iCh][iSmpl] = 1.0;
1017 else
1018 fNoisy_float[iCr][iCh][iSmpl] = 0.0;
1019 }
1020 }
1021 }
1022
1023 fvec *fSMode0_vec, *fCMode0_vec;
1024 fvec *fPedVal_vec, *fAdc_vec, *fNoisyChipChannels_vec;
1025 for (Int_t iter = 0; iter < fNIters; ++iter) {
1026 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
1027 fSMode0_vec = (fvec*)fSMode0[iCr];
1028 fCMode0_vec = (fvec*)fCMode0[iCr];
1029 fNvals_vec = (fvec*)fNvals[iCr];
1030 for (Int_t iCh = 0; iCh < fNChannels / VecLenSIMD; iCh++) {
1031 fSMode0_vec[iCh] = 0.0;
1032 fCMode0_vec[iCh] = 0.0;
1033 fNvals_vec[iCh] = 0;
1034 }
1035 fSMode0[iCr] = (Float_t*)fSMode0_vec;
1036 fCMode0[iCr] = (Float_t*)fCMode0_vec;
1037 fNvals[iCr] = (Float_t*)fNvals_vec;
1038 }
1039 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
1040 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
1041 fPedVal_vec = (fvec*)fPedVal[iCr][iCh];
1042 fAdc_vec = (fvec*)fAdc[iCr][iCh];
1043 fNoisyChipChannels_vec = (fvec*)fNoisy_float[iCr][iCh];
1044 Float_t cs = fCMode[iCr][iCh] - fSMode[iCr][iCh];
1045 fvec cs_vec = cs;
1046 fvec fSMode0_sum = 0.0;
1047 fvec fCMode0_sum = 0.0;
1048 fvec Nval = 0.0;
1049 Float_t sum1;
1050 Float_t sum2;
1051 Float_t sum3;
1052
1053 for (Int_t iSmpl = 0; iSmpl < fNSamples / VecLenSIMD; iSmpl++) {
1054 // if (fPedVal[iCr][iCh][iSmpl] == 0 || fNoisyChipChannels[iCr][iCh][iSmpl] == kTRUE) continue;
1055 fvec sig =
1056 _mm_blendv_ps(fAdc_vec[iSmpl] - fPedVal_vec[iSmpl] + cs_vec, outThresh,
1057 fvec(fPedVal_vec[iSmpl] == 0.0) | fvec(fNoisyChipChannels_vec[iSmpl] == 1.0));
1058 fvec Asig = fabs(sig);
1059 Float_t thr = thrMax - iter * thrDif;
1060 fvec thr_vec = thr;
1061 Nval = _mm_blendv_ps(Nval, Nval + 1.0, fvec(Asig < thr_vec));
1062 fSMode0_sum = _mm_blendv_ps(fSMode0_sum, fSMode0_sum + fAdc_vec[iSmpl], fvec(Asig < thr_vec));
1063 fCMode0_sum = _mm_blendv_ps(fCMode0_sum, fCMode0_sum + fPedVal_vec[iSmpl], fvec(Asig < thr_vec));
1064 }
1065 fSMode0_sum = _mm_dp_ps(fSMode0_sum, _f32vec4_one, 0xFF);
1066 _mm_store_ss(&sum1, fSMode0_sum);
1067 fSMode0[iCr][iCh] += sum1;
1068 fCMode0_sum = _mm_dp_ps(fCMode0_sum, _f32vec4_one, 0xFF);
1069 _mm_store_ss(&sum2, fCMode0_sum);
1070 fCMode0[iCr][iCh] += sum2;
1071 Nval = _mm_dp_ps(Nval, _f32vec4_one, 0xFF);
1072 _mm_store_ss(&sum3, Nval);
1073 fNvals[iCr][iCh] += sum3;
1074 }
1075
1076 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
1077 fNvals_vec = (fvec*)fNvals[iCr];
1078 fSMode_vec = (fvec*)fSMode[iCr];
1079 fCMode_vec = (fvec*)fCMode[iCr];
1080 fSMode0_vec = (fvec*)fSMode0[iCr];
1081 fCMode0_vec = (fvec*)fCMode0[iCr];
1082 for (Int_t iCh = 0; iCh < fNChannels / VecLenSIMD; iCh++) {
1083 fSMode_vec[iCh] = _mm_blendv_ps(_f32vec4_zero, fSMode0_vec[iCh] / fNvals_vec[iCh],
1084 fvec(fNvals_vec[iCh] > fNvalsMin_vec));
1085 fCMode_vec[iCh] = _mm_blendv_ps(_f32vec4_zero, fCMode0_vec[iCh] / fNvals_vec[iCh],
1086 fvec(fNvals_vec[iCh] > fNvalsMin_vec));
1087 }
1088 fNvals[iCr] = (Float_t*)fNvals_vec;
1089 fSMode[iCr] = (Float_t*)fSMode_vec;
1090 fCMode[iCr] = (Float_t*)fCMode_vec;
1091 // if (iter == fEndIter - 1) {
1092 // for (Int_t iCh = 0; iCh < fNChannels; iCh++) {
1093 // Int_t ic = iCh * fNSamples + fNSamples / 2;
1094 // // printf("iCh %4d iSmpl %4d ic %4d cmod %5f smod %5f\n",
1095 // iCh, iSmpl, ic, fCMode[iCr][iCh], fSMode[iCr][iCh]); hCModeSi[iCr]->Fill(ic,
1096 // fCMode[iCr][iCh]); hSModeSi[iCr]->Fill(ic, fSMode[iCr][iCh]);
1097 // }
1098 // }
1099 }
1100 }
1101
1102 for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
1103 for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
1104 delete[] fNoisy_float[iCr][iCh];
1105 }
1106 delete[] fNoisy_float[iCr];
1107 }
1108 delete[] fNoisy_float;
1109
1110 timer.Stop();
1111 workTime_cpu += timer.CpuTime();
1112 workTime_real += timer.RealTime();
1113}
1114
1115Double_t BmnAdcProcessor::CalcCMS(Double_t* samples, Int_t size)
1116{
1117
1118 const UShort_t kNITER = 10;
1119 Double_t CMS = 0.0;
1120 UInt_t nStr = size;
1121 Double_t upd[size];
1122 for (Int_t iSmpl = 0; iSmpl < size; ++iSmpl)
1123 upd[iSmpl] = samples[iSmpl];
1124
1125 for (Int_t itr = 0; itr < kNITER; ++itr) {
1126 if (nStr == 0)
1127 break;
1128 Double_t cms = 0.0; // common mode shift
1129 for (UInt_t iSmpl = 0; iSmpl < nStr; ++iSmpl)
1130 cms += upd[iSmpl];
1131 cms /= nStr;
1132 Double_t rms = 0.0; // chip noise
1133 for (UInt_t iSmpl = 0; iSmpl < nStr; ++iSmpl)
1134 rms += (upd[iSmpl] - cms) * (upd[iSmpl] - cms);
1135 rms = Sqrt(rms / nStr);
1136
1137 UInt_t nOk = 0;
1138 for (UInt_t iSmpl = 0; iSmpl < nStr; ++iSmpl)
1139 if (Abs(upd[iSmpl] - cms) < 3 * rms)
1140 upd[nOk++] = upd[iSmpl];
1141 nStr = nOk;
1142 CMS = cms;
1143 }
1144 return CMS;
1145}
1146
1148{
1149 BmnCalibData data;
1150 vector<vector<vector<bool>>>& noise_obj(data.GetNoise());
1151 vector<vector<vector<Float_t>>>& calib_obj(data.GetCalibration());
1152 noise_obj.resize(fNSerials, vector<vector<bool>>(fNChannels, vector<bool>(fNSamples, false)));
1153 calib_obj.resize(fNSerials, vector<vector<Float_t>>(fNChannels, vector<Float_t>(fNSamples, 0.0)));
1154 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
1155 for (Int_t iCh = 0; iCh < fNChannels; ++iCh)
1156 for (Int_t iSmpl = 0; iSmpl < fNSamples; ++iSmpl) {
1157 noise_obj[iCr][iCh][iSmpl] = fNoisyChipChannels[iCr][iCh][iSmpl];
1158 calib_obj[iCr][iCh][iSmpl] = fPedVal[iCr][iCh][iSmpl];
1159 }
1160 calibFile->cd(); // ??? global variables
1161 calibFile->WriteObject(&data, (fDetName + ".Calibration").Data());
1162 return kBMNSUCCESS;
1163}
1164
1166{
1167 BmnCalibData* data(nullptr);
1168 calibFile->GetObject((fDetName + ".Calibration").Data(), data);
1169 if (data) {
1170 for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
1171 for (Int_t iCh = 0; iCh < fNChannels; ++iCh)
1172 for (Int_t iSmpl = 0; iSmpl < fNSamples; ++iSmpl) {
1173 // printf("iCr %3d iCh %3d iSmpl %4d\n", iCr, iCh,iSmpl);
1174 fNoisyChipChannels[iCr][iCh][iSmpl] = data->GetNoise()[iCr][iCh][iSmpl];
1175 fPedVal[iCr][iCh][iSmpl] = data->GetCalibration()[iCr][iCh][iSmpl];
1176 }
1177 delete data;
1178 } else
1179 return kBMNERROR;
1180 // for (Int_t iCr = 0; iCr < fNSerials; ++iCr) {
1181 // printf("\n");
1182 // for (Int_t iCh = 0; iCh < fNChannels; ++iCh) {
1183 // printf("\n\t");
1184 // for (Int_t iSmpl = 0; iSmpl < fNSamples; ++iSmpl) {
1185 // printf(" %3.2f ", fPedVal[iCr][iCh][iSmpl]);
1186 // }
1187 // }
1188 // }
1189 return kBMNSUCCESS;
1190}
1191
1193{
1194 // if (!data)
1195 // return kBMNERROR;
1196 // for (auto )
1197 // for (Int_t iCr = 0; iCr < fNSerials; ++iCr)
1198 // for (Int_t iCh = 0; iCh < fNChannels; ++iCh)
1199 // for (Int_t iSmpl = 0; iSmpl < fNSamples; ++iSmpl) {
1200 // // printf("iCr %3d iCh %3d iSmpl %4d\n", iCr, iCh,iSmpl);
1201 // fNoisyChipChannels[iCr][iCh][iSmpl] = data->GetNoise()[iCr][iCh][iSmpl];
1202 // fPedVal[iCr][iCh][iSmpl] = data->GetCalibration()[iCr][iCh][iSmpl];
1203 // }
1204 return kBMNSUCCESS;
1205}
1206
1207void BmnAdcProcessor::SetThreshold(Double_t final_thr, Double_t thr_dif, Double_t n_iters, Double_t cmod_cut)
1208{
1209 FinalThr = final_thr;
1210 if (thr_dif >= 0)
1211 thrDif = thr_dif;
1212 if (n_iters >= 0)
1213 fNIters = n_iters;
1214 if (cmod_cut >= 0)
1215 cmodcut = cmod_cut;
1216 thrMax = FinalThr + (fNIters - 1) * thrDif;
1217 return;
1218}
const float thr
void memset(T *dest, T i, size_t num)
uses binary expansion of copied volume for speed up
Definition L1Grid.h:25
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
#define _f32vec4_zero
Definition P4_F32vec4.h:38
int i
Definition P4_F32vec4.h:22
#define _f32vec4_one
Definition P4_F32vec4.h:39
F32vec4 fvec
Definition P4_F32vec4.h:231
BmnStatus
Definition BmnEnums.h:24
@ kBMNERROR
Definition BmnEnums.h:26
@ kBMNSUCCESS
Definition BmnEnums.h:25
@ kBMNSETUP
Definition BmnEnums.h:90
UShort_t GetChannel() const
Definition BmnADCDigit.h:33
UInt_t GetSerial() const
Definition BmnADCDigit.h:31
Short_t * GetShortValue() const
Definition BmnADCDigit.h:39
UShort_t * GetUShortValue() const
Definition BmnADCDigit.h:37
void PrecalcEventMods(TClonesArray *adc)
Double_t *** fPedRms
Double_t **** fPedDat
data set to calculate pedestals
BmnAdcProcessor(uint32_t period, uint32_t run, TString det, Int_t nCh, const Int_t nSmpl)
BmnStatus RecalculatePedestalsAugmented()
UInt_t *** fNvalsCMod
Bool_t *** fNoisyChipChannels
Double_t *** fPedCMod
Double_t *** fPedValTemp
Double_t *** fPedCMod2
static BmnStatus DrawFilterInfo(BmnCalibData *data)
void SetThreshold(Double_t final_thr, Double_t thr_dif=-1, Double_t n_iters=-1, Double_t cmod_cut=-1)
BmnStatus SaveFilterInfo(TFile *cTalibFile)
Float_t *** fPedVal
vector< TH1 * > hSModeSi
vector< TH1 * > hPedSi
vector< TH1 * > hCModeSi
vector< vector< TH1 * > > hCMode
void CopyData2PedMap(TClonesArray *adc, UInt_t ev)
Double_t ** fSumRmsV
void GrabNewSerial(UInt_t serial)
const Int_t fNSamples
vector< UInt_t > fAdcSerials
list of serial id for ADC-detector
void PrecalcEventMods_simd(TClonesArray *adc)
void PrecalcEventModsOld(TClonesArray *adc)
vector< vector< TH1 * > > hPedLine
vector< vector< TH1 * > > hSMode
unordered_map< UInt_t, Int_t > fSerMap
ADC serials map.
vector< vector< TH1 * > > hPedLineSi
void(BmnAdcProcessor::* PrecalcEventModsImp)(TClonesArray *adc)
Double_t CalcCMS(Double_t *samples, Int_t size)
BmnStatus LoadFilterInfo(TFile *cTalibFile)
BmnStatus RecalculatePedestals()
UInt_t GetBoundaryRun(UInt_t nSmpl)
FloatVecADC & GetCalibration()
BitVecADC & GetNoise()
#define N_EV_FOR_PEDESTALS
#define ADC32_N_SAMPLES
#define N_EV_FOR_PEDESTALS
#define PAD_WIDTH_SIL
#define PAD_HEIGHT_SIL
-clang-format