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;
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)
406 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl) {
407 if (
fPedDat[iCr][iEv][iCh][iSmpl] == 0)
409 signals[iSmpl] =
fPedDat[iCr][iEv][iCh][iSmpl];
412 Double_t CMS =
CalcCMS(signals, nOk);
413 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl) {
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)
427 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl) {
428 if (
fPedDat[iCr][iEv][iCh][iSmpl] == 0)
430 signals[iSmpl] =
fPedDat[iCr][iEv][iCh][iSmpl];
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));
441 for (Int_t iCr = 0; iCr <
fNSerials; ++iCr)
443 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl)
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)
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;
513 for (Int_t iCr = 0; iCr <
fNSerials; ++iCr)
514 for (Int_t iCh = 0; iCh <
fNChannels; ++iCh) {
516 for (Int_t iSmpl = 0; iSmpl <
fNSamples; ++iSmpl) {
517 fPedVal[iCr][iCh][iSmpl] = 0.0;
519 fPedRms[iCr][iCh][iSmpl] = 0.0;
520 fAdcProfiles[iCr][iCh][iSmpl] = 0;
526 for (Int_t iCr = 0; iCr <
fNSerials; ++iCr)
527 for (Int_t iCh = 0; iCh <
fNChannels; ++iCh) {
529 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl) {
539 for (Int_t iCr = 0; iCr <
fNSerials; ++iCr)
541 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl)
547 Double_t sumRms =
thrMax / 3.5;
548 for (Int_t iter = 0; iter <
fNIters; iter++) {
553 printf(
"iter %d thr %4.2f\n", iter,
thr);
554 UInt_t nFiltered = 0;
556 for (Int_t iCr = 0; iCr <
fNSerials; ++iCr) {
558 for (Int_t iCh = 0; iCh <
fNChannels; ++iCh) {
568 for (Int_t iCr = 0; iCr <
fNSerials; ++iCr) {
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) {
585 Double_t Asig = TMath::Abs(
fPedDat[iCr][iEv][iCh][iSmpl] -
fPedVal[iCr][iCh][iSmpl]);
596 for (Int_t iCr = 0; iCr <
fNSerials; ++iCr)
597 for (Int_t iCh = 0; iCh <
fNChannels; ++iCh) {
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];
617 Double_t Asig = TMath::Abs(sig);
644 for (Int_t iCr = 0; iCr <
fNSerials; ++iCr)
645 for (Int_t iCh = 0; iCh <
fNChannels; ++iCh) {
648 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl) {
663 fPedVal[iCr][iCh][iSmpl] = 0.0;
678 for (Int_t iCr = 0; iCr <
fNSerials; ++iCr)
679 for (Int_t iCh = 0; iCh <
fNChannels; ++iCh) {
681 printf(
"cr %8X ich %2d fSumRmsV %4f sumRms %4f\n",
fAdcSerials[iCr], iCh,
fSumRmsV[iCr][iCh],
683 for (Int_t iSmpl = 0; iSmpl < nSmpl; ++iSmpl) {
687 printf(
"cr %8X ich %2d iSmpl %2d fPedCMod2 %4f\n",
fAdcSerials[iCr], iCh, iSmpl,
694 printf(
"new noisy ch on cr %i %08X ch %i smpl %i by signal %4.2f\n", iCr,
816 const Int_t VecLenSIMD = 4;
819 for (Int_t iCr = 0; iCr <
fNSerials; ++iCr) {
824 for (Int_t iCh = 0; iCh <
fNChannels; ++iCh) {
835 fvec *fPedVal_vec, *fAdc_vec, *fNoisyChipChannels_vec;
836 Float_t*** fNoisy_float;
839 for (Int_t iCr = 0; iCr <
fNSerials; ++iCr) {
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) {
846 fNoisy_float[iCr][iCh][iSmpl] = 1.0;
848 fNoisy_float[iCr][iCh][iSmpl] = 0.0;
852 __m128 outThresh = _mm_set1_ps(10000);
853 for (Int_t iAdc = 0; iAdc < adc->GetEntriesFast(); ++iAdc) {
858 auto serIter =
fSerMap.find(ser);
859 if (serIter ==
fSerMap.end()) {
864 Int_t iCr = serIter->second;
866 for (Int_t iSmpl = 0; iSmpl <
fNSamples; iSmpl++) {
869 Float_t val =
static_cast<Float_t
>(adcDig->
GetShortValue()[iSmpl] / 16);
870 fAdc[iCr][iCh][iSmpl] = val;
875 fNoisyChipChannels_vec = (
fvec*)fNoisy_float[iCr][iCh];
876 fvec fSMode_sum = 0.0;
877 fvec fCMode_sum = 0.0;
883 for (Int_t iSmpl = 0; iSmpl <
fNSamples / VecLenSIMD; iSmpl++) {
887 fvec sig = _mm_blendv_ps(fAdc_vec[iSmpl] - fPedVal_vec[iSmpl], outThresh,
888 fvec(fNoisyChipChannels_vec[iSmpl] == 1.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));
899 _mm_store_ss(&sum1, fSMode_sum);
902 _mm_store_ss(&sum2, fCMode_sum);
905 _mm_store_ss(&sum3, Nval);
908 for (Int_t iCr = 0; iCr <
fNSerials; ++iCr) {
909 for (Int_t iCh = 0; iCh <
fNChannels; ++iCh) {
910 delete[] fNoisy_float[iCr][iCh];
912 delete[] fNoisy_float[iCr];
914 delete[] fNoisy_float;
983 const Int_t VecLenSIMD = 4;
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) {
995 for (Int_t iCh = 0; iCh <
fNChannels / VecLenSIMD; iCh++) {
997 _mm_blendv_ps(
_f32vec4_zero, fSMode_vec[iCh] / fNvals_vec[iCh],
fvec(fNvals_vec[iCh] > fNvalsMin_vec));
999 _mm_blendv_ps(
_f32vec4_zero, fCMode_vec[iCh] / fNvals_vec[iCh],
fvec(fNvals_vec[iCh] > fNvalsMin_vec));
1001 fNvals[iCr] = (Float_t*)fNvals_vec;
1002 fSMode[iCr] = (Float_t*)fSMode_vec;
1003 fCMode[iCr] = (Float_t*)fCMode_vec;
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) {
1016 fNoisy_float[iCr][iCh][iSmpl] = 1.0;
1018 fNoisy_float[iCr][iCh][iSmpl] = 0.0;
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) {
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;
1035 fSMode0[iCr] = (Float_t*)fSMode0_vec;
1036 fCMode0[iCr] = (Float_t*)fCMode0_vec;
1037 fNvals[iCr] = (Float_t*)fNvals_vec;
1039 for (Int_t iCr = 0; iCr <
fNSerials; ++iCr)
1040 for (Int_t iCh = 0; iCh <
fNChannels; ++iCh) {
1043 fNoisyChipChannels_vec = (
fvec*)fNoisy_float[iCr][iCh];
1046 fvec fSMode0_sum = 0.0;
1047 fvec fCMode0_sum = 0.0;
1053 for (Int_t iSmpl = 0; iSmpl <
fNSamples / VecLenSIMD; iSmpl++) {
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));
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));
1065 fSMode0_sum = _mm_dp_ps(fSMode0_sum,
_f32vec4_one, 0xFF);
1066 _mm_store_ss(&sum1, fSMode0_sum);
1068 fCMode0_sum = _mm_dp_ps(fCMode0_sum,
_f32vec4_one, 0xFF);
1069 _mm_store_ss(&sum2, fCMode0_sum);
1072 _mm_store_ss(&sum3, Nval);
1073 fNvals[iCr][iCh] += sum3;
1076 for (Int_t iCr = 0; iCr <
fNSerials; ++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));
1088 fNvals[iCr] = (Float_t*)fNvals_vec;
1089 fSMode[iCr] = (Float_t*)fSMode_vec;
1090 fCMode[iCr] = (Float_t*)fCMode_vec;
1102 for (Int_t iCr = 0; iCr <
fNSerials; ++iCr) {
1103 for (Int_t iCh = 0; iCh <
fNChannels; ++iCh) {
1104 delete[] fNoisy_float[iCr][iCh];
1106 delete[] fNoisy_float[iCr];
1108 delete[] fNoisy_float;
1118 const UShort_t kNITER = 10;
1122 for (Int_t iSmpl = 0; iSmpl < size; ++iSmpl)
1123 upd[iSmpl] = samples[iSmpl];
1125 for (Int_t itr = 0; itr < kNITER; ++itr) {
1129 for (UInt_t iSmpl = 0; iSmpl < nStr; ++iSmpl)
1133 for (UInt_t iSmpl = 0; iSmpl < nStr; ++iSmpl)
1134 rms += (upd[iSmpl] - cms) * (upd[iSmpl] - cms);
1135 rms = Sqrt(rms / nStr);
1138 for (UInt_t iSmpl = 0; iSmpl < nStr; ++iSmpl)
1139 if (Abs(upd[iSmpl] - cms) < 3 * rms)
1140 upd[nOk++] = upd[iSmpl];