14 #include <framework/datastore/StoreArray.h>
16 #include <framework/logging/Logger.h>
17 #include "trg/ecl/TrgEclFAMFit.h"
19 #include "trg/ecl/dataobjects/TRGECLHit.h"
20 #include "trg/ecl/dataobjects/TRGECLDigi0.h"
21 #include "trg/ecl/dataobjects/TRGECLFAMAna.h"
33 TrgEclFAMFit::TrgEclFAMFit(): _BeamBkgTag(0), _AnaTag(0), EventId(0)
100 double* TCFitSample =
new double [12];
101 double* preped =
new double [4];
103 int nbin_pedestal = 4;
104 int fam_sampling_interval = 125;
115 double IntervaldT = 125 * 0.001 / Nsmalldt;
117 int FitSleepCounter = 100;
118 int FitSleepThreshold = 2;
124 for (
int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
128 for (
int iShift = 20; iShift < (NSampling - 12); iShift++) {
131 if (FitSleepCounter <= FitSleepThreshold) {
continue;}
132 for (
int iFitSample = 0; iFitSample < 12; iFitSample++) {
133 int iReplace = iFitSample + iShift;
134 TCFitSample[iFitSample] = digiEnergy[iTCIdm][iReplace] * 1000.0;
136 if (pedFlag == 1 && iFitSample < 4) {
137 TCFitSample[iFitSample] = preped[iFitSample];
146 dTBin = (int)(ShiftdTBin + Nsmalldt);
147 if (dTBin < 1) {dTBin = 1;}
148 if (dTBin > 20) {dTBin = 20;}
153 for (
int iFitSample = 0; iFitSample < 12; iFitSample++) {
154 CoeffAAA +=
CoeffNoise31[dTBin - 1][iFitSample] * TCFitSample[iFitSample];
155 CoeffBBB +=
CoeffNoise32[dTBin - 1][iFitSample] * TCFitSample[iFitSample];
156 CoeffPPP +=
CoeffNoise33[dTBin - 1][iFitSample] * TCFitSample[iFitSample];
158 double deltaT = CoeffBBB / CoeffAAA;
160 ShiftdTBin = int(deltaT / IntervaldT + dTBin);
167 double condition_t = -(deltaT + dTBin * IntervaldT - fam_sampling_interval * 0.001);
169 if (fabs(condition_t) < 0.8 * (fam_sampling_interval * 0.001) && FitE >
Threshold[iTCIdm]) {
170 FitT = condition_t + (SmallOffset + iShift + nbin_pedestal - 5) * (fam_sampling_interval * 0.001);
174 double rand_sampling_correction = digiTiming[iTCIdm][iShift] + (nbin_pedestal - iShift + 32) * fam_sampling_interval;
190 delete [] TCFitSample;
207 int ta_id[20] = {1000};
208 double ttt_a[20] = {0};
209 double ttt_b[20] = {0};
210 for (
int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
214 int maxId[500] = {0};
215 for (
int iii = 0 ; iii < 500 ; iii++) {
225 for (
int iSampling = 1; iSampling < NSampling; iSampling++) {
230 if (TCDigiE[iTCIdm][iSampling] >= max) {
232 max = TCDigiE[iTCIdm][iSampling];
233 maxId[noutput] = iSampling;
238 if (count_down >= flag_down) {
239 if (count_up >= flag_up) {
245 double NoiseLevel = 0;
246 double NoiseCount = 0;
247 for (
int iNoise = 0; iNoise < 5; iNoise++) {
248 int iNoiseReplace = (maxId[noutput] - 10) + iNoise;
249 if (iNoiseReplace >= 0) {
250 NoiseLevel += TCDigiE[iTCIdm][iNoiseReplace];
254 if (NoiseCount != 0) { NoiseLevel /= NoiseCount; }
256 TCFitEnergy[iTCIdm].push_back(TCDigiE[iTCIdm][maxId[noutput]] - NoiseLevel);
257 if (!(maxId[noutput] - 1)) {
258 for (
int jSampling = 1; jSampling < maxId[noutput] + 3; jSampling++) {
259 TCDigiE[iTCIdm][jSampling] -= NoiseLevel;
262 for (
int jSampling = maxId[noutput] - 1; jSampling < maxId[noutput] + 3; jSampling++) {
263 TCDigiE[iTCIdm][jSampling] -= NoiseLevel;
267 for (
int iSearch = 0; iSearch < 5; iSearch++) {
269 if (TCDigiE[iTCIdm][maxId[noutput] - iSearch] > 0.6 *
TCFitEnergy[iTCIdm][noutput] &&
270 TCDigiE[iTCIdm][maxId[noutput] - iSearch - 1] < 0.6 *
TCFitEnergy[iTCIdm][noutput]) {
271 ta_id[noutput] = maxId[noutput] - iSearch - 1;
276 if (ta_id[noutput] == 1000) {
277 printf(
"TrgEclFAMFit::digi02> Cannot find TC Timing (TCId=%5i, E=%8.5f)!!!\n", iTCIdm - 1,
TCFitEnergy[iTCIdm][0]);
278 B2ERROR(
"TrgEclFAMFit::digi02> Cannot find TC Timing");
280 ttt_a[noutput] = TCDigiT[iTCIdm][ta_id[noutput]];
281 ttt_b[noutput] = TCDigiT[iTCIdm][ta_id[noutput] + 1];
283 (0.6 *
TCFitEnergy[iTCIdm][noutput] - TCDigiE[iTCIdm][ta_id[noutput]]) * (ttt_b[noutput] -
285 / (TCDigiE[iTCIdm][ta_id[noutput] + 1] - TCDigiE[iTCIdm][ta_id[noutput]])) - (278.7 + 2) + (
_DataBase->
GetTCFLatency(iTCIdm + 1)));
325 float max_shape_time = 563.48;
327 for (
int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
329 int maxId[500] = {0};
335 for (
int iSampling = 1; iSampling < NSampling; iSampling++) {
337 if (TCDigiEnergy[iTCIdm][iSampling] >= max) {
339 max = TCDigiEnergy[iTCIdm][iSampling];
340 maxId[noutput] = iSampling;
345 if (count_down >= flag_down) {
346 if (count_up >= flag_up) {
352 float NoiseLevel = 0;
353 float NoiseCount = 0;
354 for (
int iNoise = 0; iNoise < 42; iNoise++) {
355 int iNoiseReplace = (maxId[noutput] - 88) + iNoise;
356 if (iNoiseReplace >= 0) {
357 NoiseLevel += TCDigiEnergy[iTCIdm][iNoiseReplace];
361 if (NoiseCount != 0) { NoiseLevel /= NoiseCount; }
362 TCFitEnergy[iTCIdm].push_back(TCDigiEnergy[iTCIdm][maxId[noutput]] - NoiseLevel);
385 std::vector<int> TCId;
386 std::vector<double> RawTCTiming;
387 std::vector<double> RawTCEnergy;
388 std::vector<double> RawBeamBkgTag;
389 RawBeamBkgTag.clear();
403 for (
int ii = 0; ii < trgeclDigiArray.
getEntries(); ii++) {
406 if (
EventId != eventid) {
continue;}
407 TCId.push_back(aTRGECLDigi->
getTCId());
408 RawTCTiming.push_back(aTRGECLDigi -> getRawTiming());
409 RawTCEnergy.push_back(aTRGECLDigi -> getRawEnergy());
410 RawBeamBkgTag.push_back(aTRGECLDigi -> getBeamBkgTag());
413 for (
int iTCId = 0; iTCId < 576 ; iTCId ++) {
415 for (
int iHit = 0; iHit < hitsize; iHit++) {
416 const int rawsize = TCId.size();
417 for (
int iDigi = 0; iDigi < rawsize; iDigi++) {
418 if (TCId[iDigi] != (iTCId + 1)) {
continue;}
419 if (abs(
TCFitEnergy[iTCId][iHit] - RawTCEnergy[iDigi]) < 12) {
420 if (abs(
TCFitTiming[iTCId][iHit] - RawTCTiming[iDigi]) < 30) {
421 BeamBkgTag[iTCId].push_back(RawBeamBkgTag[iDigi]);
445 for (
int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
447 for (
int iHit = 0; iHit < hitsize; iHit++) {
451 TrgEclHitArray[m_hitNum]->setEventId(m_nEvent);
452 TrgEclHitArray[m_hitNum]->setTCId(iTCIdm + 1);
453 TrgEclHitArray[m_hitNum]->setEnergyDep(
TCFitEnergy[iTCIdm][iHit]);
454 TrgEclHitArray[m_hitNum] ->setTimeAve(
TCFitTiming[iTCIdm][iHit]);
456 TrgEclHitArray[m_hitNum]->setBeamBkgTag(
BeamBkgTag[iTCIdm][iHit]);
464 for (
int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
467 for (
int iHit = 0; iHit < hitsize; iHit++) {
471 TrgEclAnaArray[m_hitNum]->setEventId(m_nEvent);
472 TrgEclAnaArray[m_hitNum]->setTCId(iTCIdm + 1);
475 TrgEclAnaArray[m_hitNum]->setRawEnergy(
TCRawEnergy[iTCIdm][iHit]);
476 TrgEclAnaArray[m_hitNum]->setRawTiming(
TCRawTiming[iTCIdm][iHit]);
477 TrgEclAnaArray[m_hitNum]->setFitEnergy(
TCFitEnergy[iTCIdm][iHit]);
478 TrgEclAnaArray[m_hitNum]->setFitTiming(
TCFitTiming[iTCIdm][iHit]);
479 TrgEclAnaArray[m_hitNum]->setBeamBkgTag(
BeamBkgTag[iTCIdm][iHit]);