Belle II Software development
TrgEclDigitizer.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <framework/datastore/StoreArray.h>
10#include <framework/gearbox/Unit.h>
11
12#include <ecl/geometry/ECLGeometryPar.h>
13#include <ecl/dataobjects/ECLSimHit.h>
14#include <ecl/dataobjects/ECLHit.h>
15
16#include "trg/ecl/TrgEclDigitizer.h"
17
18
19#include "trg//ecl/dataobjects/TRGECLDigi0.h"
20#include "trg/ecl/dataobjects/TRGECLWaveform.h" // by shebalin
21
22#include <iostream>
23#include <math.h>
24#include <TRandom.h>
25
26using namespace std;
27using namespace Belle2;
28//using namespace TRG;
29//
30//
31//
32TrgEclDigitizer::TrgEclDigitizer():
33 m_TimeRange(0), m_SaveTCWaveForm(0), m_FADC(1), m_BeambkgTag(0)
34{
35 m_MatrixParallel.clear();
36 m_MatrixSerial.clear();
37
38 m_TCMap = new TrgEclMapping();
39 m_DataBase = new TrgEclDataBase();
40
41 for (int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
42 m_TCEnergy_tot[iTCIdm] = 0;
43 m_TCTiming_tot[iTCIdm] = 0;
44 for (int iTime = 0; iTime < 80; iTime++) {
45 m_TCEnergy[iTCIdm][iTime] = 0;
46 m_TCTiming[iTCIdm][iTime] = 0;
47 m_TCBkgContribution[iTCIdm][iTime] = 0;
48 m_TCSigContribution[iTCIdm][iTime] = 0;
49 m_TCBeambkgTag[iTCIdm][iTime] = 0;
50 }
51
52 for (int iii = 0; iii < 60 ; iii++) {
53 m_TCRawEnergy[iTCIdm][iii] = 0;
54 m_TCRawTiming[iTCIdm][iii] = 0;
55 m_TCRawBkgTag[iTCIdm][iii] = 0;
56 }
57 for (int iii = 0; iii < 64 ; iii++) {
58 m_WaveForm[iTCIdm][iii] = 0;
59 }
60 }
61}
62//
63//
64//
70//
71//
72//
73
74void
75TrgEclDigitizer::setup(int SourceOfTC)
76{
77
78 // prepare Matrix for Noise generation
79 m_DataBase->readNoiseLMatrix(m_MatrixParallel, m_MatrixSerial);
80
81 // Set TC data
82 // SourceOfTC => 1=ECLHit, 2=ECLSimHit, 3=ECLHit+TRGECLBGTCHit
83 // ("1:=ECLHit" is used for signal w/o bkg, and real time background monitor)
84 getTCHit(SourceOfTC);
85
86 return;
87}
88//
89//
90//
91void
92TrgEclDigitizer::getTCHit(int SourceOfTC)
93{
94
95 std::vector< std::vector<float> > E_cell(8736, std::vector<float>(80, 0.0));
96 std::vector< std::vector<float> > T_ave(8736, std::vector<float>(80, 0.0));
97 std::vector< std::vector<float> > Tof_ave(8736, std::vector<float>(80, 0.0));
98 std::vector< std::vector<float> > beambkg_tag(8736, std::vector<float>(80, 0.0));
99
100 int nBinTime = 80;
101 m_TimeRange = 4000; // -4us ~ 4us
102 //-------------------------------------------------------------------
103 // read Xtal data
104 //---------------------------------------------------------------------
105 if (SourceOfTC == 1) { // read ECLHit table
106
107 StoreArray<ECLHit> eclHitArray("ECLHits");
108
109 int nHits_hit = eclHitArray.getEntries();
110 //
111 for (int iHits = 0; iHits < nHits_hit; iHits++) {
112 // Get a hit
113 ECLHit* aECLHit = eclHitArray[iHits];
114 int beambkg = aECLHit->getBackgroundTag();
115
116 // Hit geom. info
117 int hitCellId = aECLHit->getCellId() - 1;
118 float hitE = aECLHit->getEnergyDep() / Unit::GeV;
119 float aveT = aECLHit->getTimeAve(); // ns :time from IP to PD
120 if (aveT < - m_TimeRange || aveT > m_TimeRange) {continue;} //Choose - TimeRange ~ TimeTange
121 int TimeIndex = (int)((aveT + m_TimeRange) / 100); //Binning : -4000 = 1st bin ~ 4000 80th bin.
122
123 int iTCIdm = m_TCMap->getTCIdFromXtalId(hitCellId + 1) - 1;
124 m_TCEnergy[iTCIdm][TimeIndex] += hitE;
125 m_TCTiming[iTCIdm][TimeIndex] += hitE * aveT;
126 if (beambkg > 0) {m_TCBkgContribution[iTCIdm][TimeIndex] += hitE ;}
127 else if (beambkg == 0) {m_TCSigContribution[iTCIdm][TimeIndex] += hitE;}
128
129 }
130 for (int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
131 for (int iTime = 0; iTime < nBinTime; iTime++) {
132
133 if (m_TCEnergy[iTCIdm][iTime] < 1e-9) {continue;}// 0.01MeV cut
134 /* cppcheck-suppress variableScope */
135 double maxbkgE = 0;
136 /* cppcheck-suppress variableScope */
137 int maxbkgtag = 0;
138 m_TCTiming[iTCIdm][iTime] /= m_TCEnergy[iTCIdm][iTime];
139 if (m_BeambkgTag == 1) {
140 if (m_TCBkgContribution[iTCIdm][iTime] < m_TCSigContribution[iTCIdm][iTime]) {
141 m_TCBeambkgTag[iTCIdm][iTime] = 0; //signal Tag : 0
142 } else {
143 for (int iXtalIdm = 0; iXtalIdm < 8736; iXtalIdm++) {
144 int iTCId = m_TCMap->getTCIdFromXtalId(iXtalIdm + 1) - 1;
145 if (iTCIdm != iTCId) {continue;}
146 if (maxbkgE < E_cell[iXtalIdm][iTime]) {
147 maxbkgE = E_cell[iXtalIdm][iTime];
148 maxbkgtag = beambkg_tag[iXtalIdm][iTime];
149 }
150 }
151 m_TCBeambkgTag[iTCIdm][iTime] = maxbkgtag;
152 }
153 }
154 }
155 }
156 }
157
158
159 if (SourceOfTC == 2) { // read ECLSimHit
160 ECL::ECLGeometryPar* eclp = ECL::ECLGeometryPar::Instance();
161 //=====================
162 // Loop over all hits of steps
163 //=====================
164 StoreArray<ECLSimHit> eclArray("ECLSimHits");
165 int nHits = eclArray.getEntries();
166 //
167 for (int iHits = 0; iHits < nHits; iHits++) {
168 // Get a hit
169 ECLSimHit* aECLSimHit = eclArray[iHits];
170
171 int hitCellId = aECLSimHit->getCellId() - 1;
172 float hitE = aECLSimHit->getEnergyDep() / Unit::GeV;
173 float hitTOF = aECLSimHit->getFlightTime() / Unit::ns;
174
175 G4ThreeVector t = aECLSimHit->getPosIn(); // [cm], Hit position in Xtal (based on from IP)
176 ROOT::Math::XYZVector HitInPos(t.x(), t.y(), t.z()); // = aECLSimHit->getPosIn(); // [cm], Hit position in Xtal (based on from IP)
177 ROOT::Math::XYZVector PosCell = eclp->GetCrystalPos(hitCellId);// [cm], Xtal position (based on from IP)
178 ROOT::Math::XYZVector VecCell = eclp->GetCrystalVec(hitCellId);
179 float local_pos_r = 15.0 - (HitInPos - PosCell).Dot(VecCell);
180 if (hitTOF < - m_TimeRange || hitTOF > m_TimeRange) {continue;}
181 int TimeIndex = (int)((hitTOF + m_TimeRange) / 100);
182 E_cell[hitCellId][TimeIndex] = E_cell[hitCellId][TimeIndex] + hitE;
183 T_ave[hitCellId][TimeIndex] = T_ave[hitCellId][TimeIndex] + hitE * local_pos_r;
184 Tof_ave[hitCellId][TimeIndex] = Tof_ave[hitCellId][TimeIndex] + hitE * hitTOF;
185 }
186 //
187 //
188 //
189 //===============
190 // Xtal energy and timing (0-8us, 0.2us interval, 40 bins)
191 //===============
192 for (int iECLCell = 0; iECLCell < 8736; iECLCell++) {
193 for (int TimeIndex = 0; TimeIndex < nBinTime; TimeIndex++) {
194
195 if (E_cell[iECLCell][TimeIndex] < 1e-9) {continue;} // 0.01MeV cut
196
197 T_ave[iECLCell][TimeIndex] = T_ave[iECLCell][TimeIndex] / E_cell[iECLCell][TimeIndex];
198 T_ave[iECLCell][TimeIndex] =
199 6.05 +
200 0.0749 * T_ave[iECLCell][TimeIndex] -
201 0.00112 * T_ave[iECLCell][TimeIndex] * T_ave[iECLCell][TimeIndex];
202 Tof_ave[iECLCell][TimeIndex] = Tof_ave[iECLCell][TimeIndex] / E_cell[iECLCell][TimeIndex];
203
204 } // 40bins, Time interval = 200ns
205 }
206 //
207 //
208 //
209 for (int iXtalIdm = 0; iXtalIdm < 8736; iXtalIdm++) {
210 int iTCIdm = m_TCMap->getTCIdFromXtalId(iXtalIdm + 1) - 1;
211 for (int iTime = 0; iTime < nBinTime; iTime++) {
212
213 if (E_cell[iXtalIdm][iTime] < 1e-9) {continue;} // 0.01MeV cut
214 m_TCEnergy[iTCIdm][iTime] += E_cell[iXtalIdm][iTime];
215 m_TCTiming[iTCIdm][iTime] += E_cell[iXtalIdm][iTime] * (T_ave[iXtalIdm][iTime]);
216 if (beambkg_tag[iXtalIdm][iTime] > 0) {m_TCBkgContribution[iTCIdm][iTime] += E_cell[iXtalIdm][iTime];}
217 if (beambkg_tag[iXtalIdm][iTime] == 0) {m_TCSigContribution[iTCIdm][iTime] += E_cell[iXtalIdm][iTime];}
218
219 }
220 }
221 for (int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
222 for (int iTime = 0; iTime < nBinTime; iTime++) {
223 /* cppcheck-suppress variableScope */
224 double maxbkgE = 0;
225 /* cppcheck-suppress variableScope */
226 int maxbkgtag = 0;
227
228 if (m_TCEnergy[iTCIdm][iTime] < 1e-9) {continue;} // 0.01MeV cut
229 m_TCTiming[iTCIdm][iTime] /= m_TCEnergy[iTCIdm][iTime];
230 if (m_BeambkgTag == 1) {
231 if (m_TCBkgContribution[iTCIdm][iTime] < m_TCSigContribution[iTCIdm][iTime]) {
232 m_TCBeambkgTag[iTCIdm][iTime] = 0; //signal Tag : 0
233 } else {
234 for (int iXtalIdm = 0; iXtalIdm < 8736; iXtalIdm++) {
235 int iTCId = m_TCMap->getTCIdFromXtalId(iXtalIdm + 1) - 1;
236 if (iTCIdm != iTCId) {continue;}
237 if (maxbkgE < E_cell[iXtalIdm][iTime]) {
238 maxbkgE = E_cell[iXtalIdm][iTime];
239 maxbkgtag = beambkg_tag[iXtalIdm][iTime];
240 }
241 }
242 m_TCBeambkgTag[iTCIdm][iTime] = maxbkgtag;
243 }
244 }
245 }
246
247 }
248
249 }
250 //--------------------------------------------------------
251 //
252 //--------------------------------------------------------
253 if (SourceOfTC == 3) {
254
255 StoreArray<ECLHit> eclHitArray("ECLHits");
256
257 int nHits_hit = eclHitArray.getEntries();
258 // signal hit
259 for (int iHits = 0; iHits < nHits_hit; iHits++) {
260 // Get a hit
261 ECLHit* aECLHit = eclHitArray[iHits];
262 // Hit geom. info
263 int hitCellId = aECLHit->getCellId() - 1;
264 float hitE = aECLHit->getEnergyDep() / Unit::GeV;
265 float aveT = aECLHit->getTimeAve(); // ns :time from IP to PD
266 // Choose - TimeRange ~ TimeTange
267 if (aveT < - m_TimeRange || aveT > m_TimeRange) {continue;}
268 // Binning : -4000 = 1st bin ~ 4000 80th bin.
269 int TimeIndex = (int)((aveT + m_TimeRange) / 100);
270
271 int iTCIdm = m_TCMap->getTCIdFromXtalId(hitCellId + 1) - 1;
272 m_TCEnergy[iTCIdm][TimeIndex] += hitE;
273 m_TCTiming[iTCIdm][TimeIndex] += hitE * aveT;
274 }
275 // background hit
276 StoreArray<TRGECLBGTCHit> m_trgeclBGTCHits("TRGECLBGTCHits_beamBG");
277 for (const TRGECLBGTCHit& ttt : m_trgeclBGTCHits) {
278 // TC timing
279 double tcbg_t = ttt.getTimeAve();
280 // Timing cut
281 if (abs(tcbg_t) > m_TimeRange) continue;
282 //Binning : -4000 = 1st bin ~ 4000 80th bin.
283 int TimeIndex = (int)((tcbg_t + m_TimeRange) / 100);
284 // TC energy
285 double tcbg_e = ttt.getEnergyDep();
286 // TC ID index
287 int iTCIdm = ttt.getTCId() - 1;
288 // TC energy and timing
289 double tc_e = m_TCEnergy[iTCIdm][TimeIndex];
290 double tc_t = m_TCTiming[iTCIdm][TimeIndex];
291 m_TCEnergy[iTCIdm][TimeIndex] += tcbg_e;
292 m_TCTiming[iTCIdm][TimeIndex] += (tcbg_e * tcbg_t + tc_e * tc_t);
293 }
294 for (int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
295 for (int iTime = 0; iTime < nBinTime; iTime++) {
296 m_TCTiming[iTCIdm][iTime] /= m_TCEnergy[iTCIdm][iTime];
297 }
298 }
299 }
300 //--------------------------
301 // TC energy and timing in t=0-1us as true values.
302 //--------------------------
303 for (int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
304 for (int iTime = 40; iTime < 50; iTime++) {
305 m_TCEnergy_tot[iTCIdm] += m_TCEnergy[iTCIdm][iTime];
306 m_TCTiming_tot[iTCIdm] += m_TCTiming[iTCIdm][iTime] * m_TCEnergy[iTCIdm][iTime];
307 }
308 m_TCTiming_tot[iTCIdm] /= m_TCEnergy_tot[iTCIdm];
309 }
310
311 return;
312}
313//
314//
315//
316void
317TrgEclDigitizer::digitization01(std::vector<std::vector<double>>& TCDigiE,
318 std::vector<std::vector<double>>& TCDigiT)
319{
320 TCDigiE.clear();
321 TCDigiE.resize(576, vector<double>(64, 0));
322 TCDigiT.clear();
323 TCDigiT.resize(576, vector<double>(64, 0));
324
325 //
326 double cut_energy_tot = 0.03; // [GeV]
327 int nbin_pedestal = 4; // = nbin_pedestal*fam_sampling_interval [ns] in total
328 double fam_sampling_interval = 125; // [ns]
329 int NSampling = 64; // # of sampling array
330
331
332 std::vector< std::vector<float> > noise_pileup(576, std::vector<float>(64, 0.0)); // [GeV]
333 std::vector< std::vector<float> > noise_parallel(576, std::vector<float>(64, 0.0)); // [GeV]
334 std::vector< std::vector<float> > noise_serial(576, std::vector<float>(64, 0.0)); // [GeV]
335 std::vector<float> X_pr(64, 0.0);
336 std::vector<float> X_sr(64, 0.0);
337
338
339 // (Make sampling time random between FAM sampling intervals)
340 double random_sampling_correction = 0; // [ns]
341 random_sampling_correction = gRandom->Rndm() * fam_sampling_interval;
342 //==================
343 // (01)Signal digitization
344 //==================
345 for (int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
346 for (int iTimeBin = 0; iTimeBin < 80; iTimeBin++) {
347
348 if (m_TCEnergy[iTCIdm][iTimeBin] < 0.0001) {continue;} // 0.1MeV cut on TC bin_energy
349 for (int iSampling = 0; iSampling < NSampling; iSampling++) {
350 // inputTiming is in [us] <-- Be careful, here is NOT [ns]
351 double inputTiming
352 = (-m_TCTiming[iTCIdm][iTimeBin] - m_TimeRange + (-nbin_pedestal + iSampling) * fam_sampling_interval) * 0.001;
353 inputTiming += random_sampling_correction * 0.001;
354 if (inputTiming < 0 || inputTiming > 2.0) {continue;} // Shaping in t0 ~t0+ 2.0 us
355 if (m_FADC == 1) {
356
357 TCDigiE[iTCIdm][iSampling] += interFADC(inputTiming) * m_TCEnergy[iTCIdm][iTimeBin];
358 } else {
359 TCDigiE[iTCIdm][iSampling] += SimplifiedFADC(0, inputTiming) * m_TCEnergy[iTCIdm][iTimeBin];
360 }
361 }
362
363 for (int iSampling = 0; iSampling < NSampling; iSampling++) {
364 TCDigiT[iTCIdm][iSampling] = (-nbin_pedestal + iSampling - m_TimeRange / fam_sampling_interval) * fam_sampling_interval;
365 TCDigiT[iTCIdm][iSampling] += random_sampling_correction;
366 }
367 }
368 }
369 //
370 //
371 //
372 if (0) {
373 FILE* f_out_dat = fopen("ztsim.no_noise.dat", "w");
374 for (int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
375 if (m_TCEnergy_tot[iTCIdm] < cut_energy_tot) { continue; } // TC energy_tot cut
376 fprintf(f_out_dat, "%5i %8.5f %8.5f %8.1f ",
377 iTCIdm + 1, m_TCEnergy_tot[iTCIdm], m_TCEnergy[iTCIdm][0], TCDigiT[iTCIdm][0]);
378 for (int iSampling = 0; iSampling < NSampling; iSampling++) {
379 fprintf(f_out_dat, "%7.4f ", TCDigiE[iTCIdm][iSampling]);
380 }
381 fprintf(f_out_dat, "\n");
382 }
383 fclose(f_out_dat);
384 }
385 //==================
386 // (01)noise embedding
387 //==================
388 double tmin_noise = -4; // original
389 double tgen = 10.3; // original
390 int bkg_level = 1030;
391 double ttt0 = 0; // [us]
392 /* cppcheck-suppress variableScope */
393 double ttt1 = 0; // [us]
394 /* cppcheck-suppress variableScope */
395 double ttt2 = 0; // [us]
396 //
397 double frac_pileup = 0.035; // pileup noise fraction?
398 double frac_parallel = 0.023; // parallel noise fraction?
399 double frac_serial = 0.055; // serial noise fraction?
400 double times_pileup = 1; // noise scale based on Belle noise.
401 double times_parallel = 3.15; // noise scale
402 double times_serial = 3.15; // noise scale
403
404 double corr_pileup = times_pileup * frac_pileup * sqrt(fam_sampling_interval * 0.001);
405 double corr_parallel = times_parallel * frac_parallel * sqrt(fam_sampling_interval * 0.001);
406 double corr_serial = times_serial * frac_serial * sqrt(fam_sampling_interval * 0.001);
407
408
409 if (0) {
410 for (int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
411 if (m_TCEnergy_tot[iTCIdm] < cut_energy_tot) { continue; } // TC energy_tot cut
412 for (int jjj = 0; jjj < bkg_level; jjj++) {
413 ttt0 = -(tmin_noise + tgen * gRandom->Rndm()); // [us]
414 ttt1 = -(tmin_noise + tgen * gRandom->Rndm()); // [us]
415 ttt2 = -(tmin_noise + tgen * gRandom->Rndm()); // [us]
416 for (int iSampling = 0; iSampling < NSampling; iSampling++) {
417 // (pile-up noise)
418 if (ttt0 > 0) { TCDigiE[iTCIdm][iSampling] += SimplifiedFADC(0, ttt0) * corr_pileup * 0.001; }
419 // (parallel noise)
420 if (ttt1 > 0) { TCDigiE[iTCIdm][iSampling] += SimplifiedFADC(1, ttt1) * corr_parallel; }
421 // (serial noise)
422 if (ttt2 > 0) { TCDigiE[iTCIdm][iSampling] += SimplifiedFADC(2, ttt2) * corr_serial; }
423 ttt0 += fam_sampling_interval * 0.001;
424 ttt1 += fam_sampling_interval * 0.001;
425 ttt2 += fam_sampling_interval * 0.001;
426 }
427 }
428 }
429 }
430 if (1) { //use L Matrix
431 for (int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
432
433 for (int iSampling = 0; iSampling < NSampling; iSampling++) {
434 X_pr[iSampling] = gRandom ->Gaus(0, 1);
435 X_sr[iSampling] = gRandom ->Gaus(0, 1);
436 }
437
438 for (int iSampling = 0; iSampling < NSampling; iSampling++) {
439 for (int jSampling = 0; jSampling < NSampling; jSampling++) {
440 noise_parallel[iTCIdm][iSampling] += 10 * corr_parallel * m_MatrixParallel[iSampling][jSampling] * X_pr[jSampling];
441 noise_serial[iTCIdm][iSampling] += 10 * corr_serial * m_MatrixSerial[iSampling][jSampling] * X_sr[jSampling];
442 }
443 }
444
445 for (int iSampling = 0; iSampling < NSampling; iSampling++) {
446
447 TCDigiE[iTCIdm][iSampling] += noise_pileup[iTCIdm][iSampling] + noise_parallel[iTCIdm][iSampling] +
448 noise_serial[iTCIdm][iSampling];
449 }
450 }
451 if (0) {
452 for (int iTCIdm = 0; iTCIdm < 576; iTCIdm++) { //Only Pile-up noise use old method.
453 if (m_TCEnergy_tot[iTCIdm] < cut_energy_tot) { continue; } // TC energy_tot cut
454 for (int jjj = 0; jjj < bkg_level; jjj++) {
455 ttt0 = -(tmin_noise + tgen * gRandom->Rndm()); // [us]
456 for (int iSampling = 0; iSampling < NSampling; iSampling++) {
457 // (pile-up noise)
458 if (ttt0 > 0) { TCDigiE[iTCIdm][iSampling] += SimplifiedFADC(0, ttt0) * corr_pileup * 0.001; }
459 ttt0 += fam_sampling_interval * 0.001;
460 }
461 }
462 }
463 }
464 }
465
466 if (m_SaveTCWaveForm == 1) {
467
468 for (int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
469 for (int iSampling = 0; iSampling < NSampling; iSampling++) {
470 m_WaveForm[iTCIdm][iSampling] = TCDigiE[iTCIdm][iSampling];
471 }
472
473 }
474 }
475
476}
477void
478TrgEclDigitizer::digitization02(std::vector<std::vector<double>>& TCDigiE, std::vector<std::vector<double>>& TCDigiT)
479{
480 //===============
481 // (03)Signal digitization (w/ 12ns interval for method-0)
482 //===============
483 TCDigiE.clear();
484 TCDigiE.resize(576, vector<double>(666, 0));
485 TCDigiT.clear();
486 TCDigiT.resize(576, vector<double>(666, 0));
487 double cut_energy_tot = 0.03; // [GeV]
488 int nbin_pedestal = 100;
489 float fam_sampling_interval = 12; // [ns]
490
491 std::vector< std::vector<float> > TCDigiEnergy(576, std::vector<float>(666, 0.0)); // [GeV]
492 std::vector< std::vector<float> > TCDigiTiming(576, std::vector<float>(666, 0.0)); // [ns]
493 int NSampling = 666;
494
495 // Make sampling time random between FAM sampling intervals
496 float random_sampling_correction = 0; // [ns]
497 random_sampling_correction = gRandom->Rndm() * fam_sampling_interval;
498 for (int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
499 if (m_TCEnergy_tot[iTCIdm] < cut_energy_tot) {continue;} // TC energy_tot cut
500 for (int iTimeBin = 0; iTimeBin < 80; iTimeBin++) {
501 if (m_TCEnergy[iTCIdm][iTimeBin] < 0.0001) { continue; } // 0.1 MeV cut on TC bin_energy
502 for (int iSampling = 0; iSampling < NSampling; iSampling++) {
503 // inputTiming is in [us] <-- Be careful, here is NOT [ns]
504 float inputTiming
505 = (-m_TCTiming[iTCIdm][iTimeBin] - m_TimeRange + (-nbin_pedestal + iSampling) * fam_sampling_interval) * 0.001;
506
507 inputTiming += random_sampling_correction * 0.001;
508 if (inputTiming < 0 || inputTiming > 2.0) {continue;} // Shaping in t0 ~t0+ 2.0 us
509
510 TCDigiEnergy[iTCIdm][iSampling] += FADC(0, inputTiming) * m_TCEnergy[iTCIdm][iTimeBin];
511 }
512 }
513 for (int iSampling = 0; iSampling < NSampling; iSampling++) {
514 TCDigiTiming[iTCIdm][iSampling] = (-nbin_pedestal + iSampling - m_TimeRange / fam_sampling_interval) * fam_sampling_interval;
515 TCDigiTiming[iTCIdm][iSampling] += random_sampling_correction;
516 }
517 }
518 //==================
519 // (03)noise embedding
520 //==================
521
522 double tmin_noise = -4; // original
523 double tgen = 10.3; //
524 int bkg_level = 1030;
525 double ttt0 = 0; // [us]
526 double ttt1 = 0; // [us]
527 double ttt2 = 0; // [us]
528 //double frac_pileup = 0.035; // pileup noise fraction?
529 //double frac_parallel = 0.023; // parallel noise fraction?
530 //double frac_serial = 0.055; // serial noise fraction?
531 //double times_pileup = 1; // noise scale based on Belle noise.
532 //double times_parallel = 1; // noise scale
533 //double times_serial = 1; // noise scale
534 //double corr_pileup = times_pileup * frac_pileup * sqrt(fam_sampling_interval * 0.001);
535 //double corr_parallel = times_parallel * frac_parallel * sqrt(fam_sampling_interval * 0.001);
536 //double corr_serial = times_serial * frac_serial * sqrt(fam_sampling_interval * 0.001);
537 double corr_pileup = 0.011068;
538 double corr_parallel = 0.00727324;
539 double corr_serial = 0.0173925;
540
541
542 for (int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
543 if (m_TCEnergy_tot[iTCIdm] < cut_energy_tot) { continue; } // 1 MeV TC energy cut
544 for (int jjj = 0; jjj < bkg_level; jjj++) {
545 ttt0 = -(tmin_noise + tgen * gRandom->Rndm()); // [us]
546 ttt1 = -(tmin_noise + tgen * gRandom->Rndm()); // [us]
547 ttt2 = -(tmin_noise + tgen * gRandom->Rndm()); // [us]
548 for (int iSampling = 0; iSampling < NSampling; iSampling++) {
549 // (pile-up noise)
550 if (ttt0 > 0) { TCDigiEnergy[iTCIdm][iSampling] += FADC(0, ttt0) * corr_pileup * 0.001; }
551 // (parallel noise)
552 if (ttt1 > 0) { TCDigiEnergy[iTCIdm][iSampling] += FADC(1, ttt1) * corr_parallel; }
553 // (serial noise)
554 if (ttt2 > 0) { TCDigiEnergy[iTCIdm][iSampling] += FADC(2, ttt2) * corr_serial; }
555 ttt0 += fam_sampling_interval * 0.001;
556 ttt1 += fam_sampling_interval * 0.001;
557 ttt2 += fam_sampling_interval * 0.001;
558 }
559 }
560 }
561
562}
563//
564//
565//
566void
567TrgEclDigitizer::save(int m_nEvent)
568{
569 //---------------
570 // Root Output
571 //---------------
572 int m_hitNum = 0;
573 for (int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
574 for (int iBinTime = 0; iBinTime < 80; iBinTime++) {
575 if (m_TCEnergy[iTCIdm][iBinTime] < 0.001) {continue;}
576 StoreArray<TRGECLDigi0> TCDigiArray;
577 TCDigiArray.appendNew();
578 m_hitNum = TCDigiArray.getEntries() - 1;
579
580 TCDigiArray[m_hitNum]->setEventId(m_nEvent);
581 TCDigiArray[m_hitNum]->setTCId(iTCIdm + 1);
582 TCDigiArray[m_hitNum]->setiBinTime(iBinTime);
583 TCDigiArray[m_hitNum]->setRawEnergy(m_TCEnergy[iTCIdm][iBinTime]);
584 TCDigiArray[m_hitNum]->setRawTiming(m_TCTiming[iTCIdm][iBinTime]);
585 TCDigiArray[m_hitNum]->setBeamBkgTag(m_TCBeambkgTag[iTCIdm][iBinTime]);
586 }
587 }
588
589 if (m_SaveTCWaveForm == 1) {
590 StoreArray<TRGECLWaveform> TCWaveformArray;
591 for (int iTCIdm = 0; iTCIdm < 576; iTCIdm++) {
592 if (iTCIdm == 80) iTCIdm = 512; // skip barrel
593 int tc_phi_id = m_TCMap->getTCPhiIdFromTCId(iTCIdm + 1);
594 int tc_theta_id = m_TCMap->getTCThetaIdFromTCId(iTCIdm + 1);
595 TRGECLWaveform* newWf =
596 TCWaveformArray.appendNew(iTCIdm + 1, m_WaveForm[iTCIdm]);
597 newWf->setThetaPhiIDs(tc_theta_id, tc_phi_id);
598 }
599
600 }
601
602
603 return;
604}
605//
606//
607//
608double
609TrgEclDigitizer::interFADC(double timing)
610{
611
612 std::vector<double> SignalAmp;
613 std::vector<double> SignalTime;
614
615 SignalAmp.clear();
616 SignalTime.clear();
617
618 timing = timing * 1000;
619 int startbin = (int)(timing) / 12;
620 int endbin = startbin + 1;
621
622 SignalAmp = { 0, 3.47505e-06, 0.000133173, 0.000939532, 0.00337321, 0.00845977, 0.0170396, 0.0296601, 0.0465713, 0.0677693, 0.0930545, 0.122086, 0.154429, 0.189595, 0.227068, 0.266331, 0.306882, 0.348243, 0.389971, 0.431664, 0.472959, 0.513539, 0.553127, 0.59149, 0.628431, 0.663791, 0.697445, 0.729297, 0.759281, 0.787354, 0.813494, 0.837698, 0.859982, 0.880372, 0.898908, 0.915639, 0.930622, 0.943919, 0.955598, 0.965731, 0.97439, 0.98165, 0.987587, 0.992275, 0.995788, 0.9982, 0.999581, 1, 0.999525, 0.998219, 0.996143, 0.993356, 0.989846, 0.985364, 0.979439, 0.971558, 0.961304, 0.948421, 0.932809, 0.914509, 0.893664, 0.870494, 0.845267, 0.818279, 0.789837, 0.76025, 0.729815, 0.698815, 0.667511, 0.636141, 0.604919, 0.574035, 0.543654, 0.513915, 0.484939, 0.456822, 0.429643, 0.403462, 0.378324, 0.354259, 0.331286, 0.309412, 0.288633, 0.26894, 0.250316, 0.232736, 0.216174, 0.200597, 0.185972, 0.172262, 0.159428, 0.147431, 0.136232, 0.12579, 0.116067, 0.107022, 0.0986191, 0.0908193, 0.0835871, 0.0768874, 0.0706867, 0.0649528, 0.0596551, 0.0547641, 0.0502523, 0.0460932, 0.042262, 0.0387353, 0.0354909, 0.0325082, 0.0297677, 0.0272511, 0.0249416, 0.0228232, 0.020881, 0.0191014, 0.0174714, 0.0159792, 0.0146138, 0.0133649, 0.012223, 0.0111793, 0.0102258, 0.00935504, 0.00856, 0.00783438, 0.00717232, 0.00656842, 0.00601773, 0.00551567, 0.00505807, 0.00464106, 0.00426114, 0.00391506, 0.00359985, 0.0033128, 0.00305143, 0.00281346, 0.00259681, 0.00239957, 0.00222002, 0.00205655, 0.00190773, 0.00177223, 0.00164885, 0.00153648, 0.00143413, 0.00134087, 0.00125589, 0.00117842, 0.00110777, 0.00104332, 0.000984488, 0.000930766, 0.000881678, 0.000836798, 0.000795734, 0.000758134, 0.000723677, 0.00069207, 0.000663051, 0.000636377, 0.000611832, 0.000589219, 0.000568358, 0.000549086, 0.000531258, 0.000514738, 0.000499407, 0.000485156, 0.000471884, 0.000459502, 0.00044793, 0.000437092, 0.000426923, 0.000417363, 0.000408356, 0.000399853, 0.000391811, 0.000384187, 0.000376946, 0.000370055, 0.000363483, 0.000357203, 0.000351192, 0.000345426, 0.000339886, 0.000334554, 0.000329413, 0.000324448, 0.000319645, 0.000314993, 0.000310481, 0.000306098, 0.000301836, 0.000297686, 0.00029364, 0.000289693, 0.000285837, 0.000282068, 0.00027838, 0.000274768, 0.000271229, 0.000267759, 0.000264354, 0.000261011, 0.000257727, 0.000254499, 0.000251326, 0.000248204, 0.000245132, 0.000242108, 0.000239131, 0.000236198, 0.000233308, 0.00023046, 0.000227653, 0.000224885, 0.000222155, 0.000219463, 0.000216807, 0.000214187, 0.000211602, 0.00020905, 0.000206532, 0.000204046, 0.000201593, 0.00019917, 0.000196779, 0.000194417, 0.000192085, 0.000189782, 0.000187508, 0.000185262, 0.000183044, 0.000180853, 0.000178689, 0.000176552, 0.000174441, 0.000172355, 0.000170295, 0.00016826, 0.000166249, 0.000164263, 0.000162301, };
623
624 SignalTime = { 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 156, 168, 180, 192, 204, 216, 228, 240, 252, 264, 276, 288, 300, 312, 324, 336, 348, 360, 372, 384, 396, 408, 420, 432, 444, 456, 468, 480, 492, 504, 516, 528, 540, 552, 564, 576, 588, 600, 612, 624, 636, 648, 660, 672, 684, 696, 708, 720, 732, 744, 756, 768, 780, 792, 804, 816, 828, 840, 852, 864, 876, 888, 900, 912, 924, 936, 948, 960, 972, 984, 996, 1008, 1020, 1032, 1044, 1056, 1068, 1080, 1092, 1104, 1116, 1128, 1140, 1152, 1164, 1176, 1188, 1200, 1212, 1224, 1236, 1248, 1260, 1272, 1284, 1296, 1308, 1320, 1332, 1344, 1356, 1368, 1380, 1392, 1404, 1416, 1428, 1440, 1452, 1464, 1476, 1488, 1500, 1512, 1524, 1536, 1548, 1560, 1572, 1584, 1596, 1608, 1620, 1632, 1644, 1656, 1668, 1680, 1692, 1704, 1716, 1728, 1740, 1752, 1764, 1776, 1788, 1800, 1812, 1824, 1836, 1848, 1860, 1872, 1884, 1896, 1908, 1920, 1932, 1944, 1956, 1968, 1980, 1992, 2004, 2016, 2028, 2040, 2052, 2064, 2076, 2088, 2100, 2112, 2124, 2136, 2148, 2160, 2172, 2184, 2196, 2208, 2220, 2232, 2244, 2256, 2268, 2280, 2292, 2304, 2316, 2328, 2340, 2352, 2364, 2376, 2388, 2400, 2412, 2424, 2436, 2448, 2460, 2472, 2484, 2496, 2508, 2520, 2532, 2544, 2556, 2568, 2580, 2592, 2604, 2616, 2628, 2640, 2652, 2664, 2676, 2688, 2700, 2712, 2724, 2736, 2748, 2760, 2772, 2784, 2796, 2808, 2820, 2832, 2844, 2856, 2868, 2880, 2892, 2904, 2916, 2928, };
625
626
627 double E_out = 0;
628 if (timing < 0 || timing > 2928) {return E_out;}
629
630 if (timing > SignalTime[startbin] && timing < SignalTime[endbin]) {
631 E_out = (((SignalAmp[endbin] - SignalAmp[startbin])) / (SignalTime[endbin] - SignalTime[startbin])) *
632 (timing - SignalTime[startbin]) + SignalAmp[startbin];
633 return E_out;
634 } else {
635 E_out = 0;
636 }
637
638
639 return E_out;
640
641}
642//
643//
644//
645double
646TrgEclDigitizer::FADC(int flag_gen,
647 double timing)
648{
649
650 //--------------------------------------
651 //
652 // o "timing" unit is [us]
653 // o flag_gen = 0(=signal), 1(=parallel), 2(=serial)
654 // o return value(PDF) is [GeV]
655 //
656 //--------------------------------------
657 double tsh, dd;
658 static double tc, tc2, tsc, tris, b1, b2;
659 static double amp, dft, as;
660 static double td, t1, t2;
661
662 static int ifir = 0;
663
664 if (ifir == 0) {
665
666 td = 0.10; // diff time ( 0.10)
667 t1 = 0.10; // integ1 real ( 0.10)
668 b1 = 30.90; // integ1 imag ( 30.90)
669 t2 = 0.01; // integ2 real ( 0.01)
670 b2 = 30.01; // integ2 imag ( 30.01)
671 double ts = 1.00; // scint decay ( 1.00)
672 dft = 0.600; // diff delay ( 0.600)
673 as = 0.548; // diff frac ( 0.548)
674 //
675 amp = 1.0;
676 tris = 0.01;
677 tsc = ts;
678 //
679 int im = 0;
680 int ij = 0;
681 int fm = 0;
682 tc = 0;
683 double tt = u_max(u_max(td, t1), u_max(t2, ts)) * 2;
684 int flag_once = 0;
685 while (flag_once == 0) {
686 double dt = tt / 1000;
687 double tm = 0;
688 for (int j = 1; j <= 1000; j++) {
689 tc2 = tc - dft;
690 double fff =
691 (ShapeF(tc, t1, b1, t2, b2, td, tsc) -
692 ShapeF(tc, t1, b1, t2, b2, td, tris) * 0.01) -
693 (ShapeF(tc2, t1, b1, t2, b2, td, tsc) -
694 ShapeF(tc2, t1, b1, t2, b2, td, tris) * 0.01) * as;
695 if (fff > fm) {
696 fm = fff;
697 tm = tc;
698 im = j;
699 }
700 tc = tc + dt;
701 }
702 if (im >= 1000) {
703 tt = 2 * tt;
704 flag_once = 0;
705 continue;
706 }
707 if (ij == 0) {
708 ij = 1;
709 tc = 0.99 * tm;
710 dt = tm * 0.02 / 1000;
711 flag_once = 0;
712 continue;
713 }
714 flag_once = 1;
715 }
716 amp = 1.0 / fm;
717 ifir = 1;
718 }
719 //
720 //
721 double pdf = 0;
722 if (flag_gen == 0) {
723 //-----<signal>
724 tc2 = timing - dft;
725 pdf = amp * (
726 (ShapeF(timing, t1, b1, t2, b2, td, tsc) -
727 ShapeF(timing, t1, b1, t2, b2, td, tris) * 0.01) -
728 (ShapeF(tc2, t1, b1, t2, b2, td, tsc) -
729 ShapeF(tc2, t1, b1, t2, b2, td, tris) * 0.01) * as);
730 } else if (flag_gen == 1) {
731 //-----<parallel>
732 tc2 = timing - dft;
733 tsh = 0.001;
734 pdf = amp * (
735 ShapeF(timing, t1, b1, t2, b2, td, tsh) -
736 ShapeF(tc2, t1, b1, t2, b2, td, tsh) * as);
737 pdf = pdf * 0.001; // GeV
738 } else {
739 //-----<serial>
740 tc2 = timing - dft;
741 tsh = 0.001;
742 pdf = amp * (
743 ShapeF(timing, t1, b1, t2, b2, td, tsh) -
744 ShapeF(tc2, t1, b1, t2, b2, td, tsh) * as);
745 //
746 tc = timing - 0.01;
747 if (tc < 0) { tc = 0; }
748 dd = timing - tc;
749 tc2 = tc - dft;
750 pdf = (amp * (
751 ShapeF(tc, t1, b1, t2, b2, td, tsh) -
752 ShapeF(tc2, t1, b1, t2, b2, td, tsh) * as) - pdf) / dd;
753 pdf = pdf * 0.001; // GeV
754 }
755
756 return pdf;
757}
758//
759//
760//
761double
762TrgEclDigitizer::ShapeF(double t00,
763 double t01,
764 double tb1,
765 double t02,
766 double tb2,
767 double td1,
768 double ts1)
769{
770 double dr;
771 double dzna;
772 double das1, das0, dac0;
773 double dcs0s, dsn0s, dcs0d, dsn0d;
774 double dcs1s, dsn1s, dcs1d, dsn1d;
775
776 double sv123 = 0.0;
777 if (t00 < 0) return 0;
778
779 dr = (ts1 - td1) / td1;
780 if (fabs(dr) <= 1.0e-5) {
781 if (ts1 > td1) { ts1 = td1 * 1.00001; }
782 else { ts1 = td1 * 0.99999; }
783 }
784 //
785 dr = (pow((t01 - t02), 2) + pow((tb1 - tb2), 2)) / (pow((t01), 2) + pow((tb1), 2));
786 if (dr <= 1.e-10) {
787 if (t01 < t02) { t01 = t02 / 1.00001; }
788 else { t01 = t02 / 0.99999; }
789 }
790 if (t00 <= 0.0) return 0;
791 //
792 //
793 //
794 double a1 = 1 / t01;
795 double a2 = 1 / tb1;
796 double b1 = 1 / t02;
797 double b2 = 1 / tb2;
798 double c1 = 1 / td1;
799 double c2 = 1 / ts1;
800
801 das0 = b2 * (pow((b1 - a1), 2) + (b2 + a2) * (b2 - a2));
802 dac0 = -2 * (b1 - a1) * a2 * b2;
803 das1 = a2 * (pow((b1 - a1), 2) - (b2 + a2) * (b2 - a2));
804
805 dsn0s = ((c2 - a1) * das0 - (-a2) * dac0) / (pow(a2, 2) + pow((c2 - a1), 2));
806 dcs0s = ((-a2) * das0 + (c2 - a1) * dac0) / (pow(a2, 2) + pow((c2 - a1), 2));
807
808 dsn1s = ((c2 - b1) * das1 - (-b2) * (-dac0)) / (pow(b2, 2) + pow((c2 - b1), 2));
809 dcs1s = ((-b2) * das1 + (c2 - b1) * (-dac0)) / (pow(b2, 2) + pow((c2 - b1), 2));
810
811 dsn0d = ((c1 - a1) * das0 - (-a2) * dac0) / (pow(a2, 2) + pow((c1 - a1), 2));
812 dcs0d = ((-a2) * das0 + (c1 - a1) * dac0) / (pow(a2, 2) + pow((c1 - a1), 2));
813
814 dsn1d = ((c1 - b1) * das1 - (-b2) * (-dac0)) / (pow(b2, 2) + pow((c1 - b1), 2));
815 dcs1d = ((-b2) * das1 + (c1 - b1) * (-dac0)) / (pow(b2, 2) + pow((c1 - b1), 2));
816 //
817 //
818 //
819 dzna = (pow((b1 - a1), 2) + pow((b2 - a2), 2)) * (pow((b1 - a1), 2) + pow((a2 + b2), 2));
820
821 sv123 = (
822 (dcs0s + dcs1s) * exp(-c2 * t00) * (-1) +
823 (dcs0d + dcs1d) * exp(-c1 * t00) +
824 ((dsn0s - dsn0d) * sin(a2 * t00) + (dcs0s - dcs0d) * cos(a2 * t00)) * exp(-a1 * t00) +
825 ((dsn1s - dsn1d) * sin(b2 * t00) + (dcs1s - dcs1d) * cos(b2 * t00)) * exp(-b1 * t00)
826 )
827 / dzna / (1 / c2 - 1 / c1);
828
829 return sv123;
830}
831
832//
833//
834//
835
836
837double
838TrgEclDigitizer::SimplifiedFADC(int flag_gen,
839 double timing)
840{
841
842 //--------------------------------------
843 //
844 // o "timing" unit is [us]
845 // o flag_gen = 0(=signal), 1(=parallel), 2(=serial)
846 // o return value(PDF) is [GeV]
847 // o Generate signal shape using a simplified function.
848 //
849 //
850 //--------------------------------------
851 double tsh, dd;
852 static double tc, tc2, tsc, tris;
853 static double amp, dft, as;
854
855
856 // int im, ij;
857
858 static int ifir = 0;
859
860 if (ifir == 0) {
861
862 const double td = 0.10; // diff time ( 0.10)
863 const double t1 = 0.10; // integ1 real ( 0.10)
864 // b1 = 30.90; // integ1 imag ( 30.90)
865 const double t2 = 0.01; // integ2 real ( 0.01)
866 // b2 = 30.01; // integ2 imag ( 30.01)
867 double ts = 1.00; // scint decay ( 1.00)
868 dft = 0.600; // diff delay ( 0.600)
869 as = 0.548; // diff frac ( 0.548)
870 //
871 amp = 1.0;
872 tris = 0.01;
873 // ts0 = 0.0;
874 tsc = ts;
875 double fm = 0;
876 //
877 int im = 0;
878 int ij = 0;
879 tc = 0;
880 double tt = u_max(u_max(td, t1), u_max(t2, ts)) * 2;
881 int flag_once = 0;
882 while (flag_once == 0) {
883 double dt = tt / 1000;
884
885 double tm = 0;
886 for (int j = 1; j <= 1000; j++) {
887 tc2 = tc - dft;
888 double fff =
889 (ShapeF(tc, tsc) -
890 ShapeF(tc, tris) * 0.01) -
891 (ShapeF(tc2, tsc) -
892 ShapeF(tc2, tris) * 0.01) * as;
893 if (fff > fm) {
894 fm = fff;
895 tm = tc;
896 im = j;
897 }
898 tc = tc + dt;
899 }
900 if (im >= 1000) {
901 tt = 2 * tt;
902 flag_once = 0;
903 continue;
904 }
905 if (ij == 0) {
906 ij = 1;
907 tc = 0.99 * tm;
908 dt = tm * 0.02 / 1000;
909 flag_once = 0;
910 continue;
911 }
912 flag_once = 1;
913 }
914 amp = 1.0 / fm;
915 ifir = 1;
916 }
917 //
918 //
919 double pdf = 0;
920 if (flag_gen == 0) {
921 //-----<signal>
922 tc2 = timing - dft;
923 pdf = amp * (
924 (ShapeF(timing, tsc) -
925 ShapeF(timing, tris) * 0.01) -
926 (ShapeF(tc2, tsc) -
927 ShapeF(tc2, tris) * 0.01) * as);
928 } else if (flag_gen == 1) {
929 //-----<parallel>
930 tc2 = timing - dft;
931 tsh = 0.001;
932 pdf = amp * (
933 ShapeF(timing, tsh) -
934 ShapeF(tc2, tsh) * as);
935 pdf = pdf * 0.001; // GeV
936 } else {
937 //-----<serial>
938 tc2 = timing - dft;
939 tsh = 0.001;
940 pdf = amp * (
941 ShapeF(timing, tsh) -
942 ShapeF(tc2, tsh) * as);
943 //
944 tc = timing - 0.01;
945 if (tc < 0) { tc = 0; }
946 dd = timing - tc;
947 tc2 = tc - dft;
948 pdf = (amp * (
949 ShapeF(tc, tsh) -
950 ShapeF(tc2, tsh) * as) - pdf) / dd;
951 pdf = pdf * 0.001; // GeV
952 }
953
954 return pdf;
955}
956
957
958double TrgEclDigitizer::ShapeF(double t00, double ts1)
959{
960 static const double realNaN = std::numeric_limits<double>::quiet_NaN();
961
962 double dzna;
963 // double das1,das0,dac0;
964 double dcs0s = realNaN, dsn0s = realNaN, dcs0d, dsn0d;
965 double dcs1s = realNaN, dsn1s = realNaN, dcs1d, dsn1d;
966 double a1, a2, b1, b2, c1, c2 = realNaN;
967 double sv123 = 0.0;
968
969 if (t00 <= 0.0) return 0;
970
971
972
973 a1 = 10 ;
974 a2 = 0.0323625 ;
975 b1 = 100;
976 b2 = 0.0333222;
977 c1 = 10;
978
979 if (ts1 == 1) {
980 c2 = 1;
981 dsn0s = -29.9897;
982 dcs0s = -0.08627;
983
984 dsn1s = -2.64784;
985 dcs1s = -0.00285194;
986 }
987 if (ts1 == 0.01) {
988 c2 = 100;
989 dsn0s = 2.999 ;
990 dcs0s = -0.00323517;
991
992 dsn1s = 5.82524;
993 dcs1s = -7866.7;
994 }
995 if (ts1 == 0.001) {
996 c2 = 1000;
997 dsn0s = 0.272636;
998 dcs0s = -0.000204983;
999
1000 dsn1s = 0.291262;
1001 dcs1s = 0.000204894;
1002 }
1003
1004 dsn0d = -5.998;
1005 dcs0d = -8340.22;
1006
1007 dsn1d = -2.91262;
1008 dcs1d = -0.00323517;
1009 //
1010 //
1011 //
1012 dzna = 6.561e+07;
1013
1014 sv123 = (
1015 (dcs0s + dcs1s) * exp(-c2 * t00) * (-1) +
1016 (dcs0d + dcs1d) * exp(-c1 * t00) +
1017 ((dsn0s - dsn0d) * sin(a2 * t00) + (dcs0s - dcs0d) * cos(a2 * t00)) * exp(-a1 * t00) +
1018 ((dsn1s - dsn1d) * sin(b2 * t00) + (dcs1s - dcs1d) * cos(b2 * t00)) * exp(-b1 * t00)
1019 )
1020 / dzna / (1 / c2 - 1 / c1);
1021
1022 return sv123;
1023}
1024//
1025//
1026//
1027double
1028TrgEclDigitizer::u_max(double aaa, double bbb)
1029{
1030 if (aaa > bbb) { return aaa; }
1031 else { return bbb; }
1032}
1033//
1034//
1035//
Belle2::ECLHit
Class to store simulated hits which equate to average of ECLSImHit on crystals input for digitization...
Definition ECLHit.h:25
Belle2::ECLHit::getTimeAve
double getTimeAve() const
Get average time.
Definition ECLHit.h:75
Belle2::ECLHit::getCellId
int getCellId() const
Get Cell ID.
Definition ECLHit.h:65
Belle2::ECLHit::getEnergyDep
double getEnergyDep() const
Get deposit energy.
Definition ECLHit.h:70
Belle2::ECLSimHit
ClassECLSimHit - Geant4 simulated hit for the ECL.
Definition ECLSimHit.h:29
Belle2::ECLSimHit::getCellId
int getCellId() const
Get Cell ID.
Definition ECLSimHit.h:86
Belle2::ECLSimHit::getFlightTime
double getFlightTime() const
Get Flight time from IP.
Definition ECLSimHit.h:101
Belle2::ECLSimHit::getEnergyDep
double getEnergyDep() const
Get Deposit energy.
Definition ECLSimHit.h:106
Belle2::ECLSimHit::getPosIn
G4ThreeVector getPosIn() const
Get Position.
Definition ECLSimHit.h:121
Belle2::ECL::ECLGeometryPar
The Class for ECL Geometry Parameters.
Definition ECLGeometryPar.h:44
Belle2::ECL::ECLGeometryPar::Instance
static ECLGeometryPar * Instance()
Static method to get a reference to the ECLGeometryPar instance.
Definition ECLGeometryPar.cc:168
Belle2::ECL::ECLGeometryPar::GetCrystalPos
ROOT::Math::XYZVector GetCrystalPos(int cid)
The Position of crystal.
Definition ECLGeometryPar.h:101
Belle2::ECL::ECLGeometryPar::GetCrystalVec
ROOT::Math::XYZVector GetCrystalVec(int cid)
The direction of crystal.
Definition ECLGeometryPar.h:109
Belle2::SimHitBase::getBackgroundTag
virtual unsigned short getBackgroundTag() const
Get background tag.
Definition SimHitBase.h:46
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition StoreArray.h:113
Belle2::StoreArray::appendNew
T * appendNew()
Construct a new T object at the end of the array.
Definition StoreArray.h:246
Belle2::StoreArray::getEntries
int getEntries() const
Get the number of objects in the array.
Definition StoreArray.h:216
Belle2::TRGECLBGTCHit
Example Detector.
Definition TRGECLBGTCHit.h:24
Belle2::TRGECLWaveform::setThetaPhiIDs
void setThetaPhiIDs(int thid, int phid)
Set Theta and Phi Id of TC.
Definition TRGECLWaveform.h:80
Belle2::TrgEclDataBase
class TrgEclDataBase;
Definition TrgEclDataBase.h:20
Belle2::TrgEclDigitizer::~TrgEclDigitizer
virtual ~TrgEclDigitizer()
Destructor.
Definition TrgEclDigitizer.cc:65
Belle2::TrgEclDigitizer::save
void save(int)
save fitting result into tables
Definition TrgEclDigitizer.cc:567
Belle2::TrgEclDigitizer::m_MatrixParallel
std::vector< std::vector< double > > m_MatrixParallel
Noise Matrix of Parallel and Serial Noise.
Definition TrgEclDigitizer.h:97
Belle2::TrgEclDigitizer::setup
void setup(int)
setup fam module
Definition TrgEclDigitizer.cc:75
Belle2::TrgEclDigitizer::m_TCTiming
double m_TCTiming[576][80]
TC Timing converted from Xtarl Timing [GeV].
Definition TrgEclDigitizer.h:72
Belle2::TrgEclDigitizer::m_TCRawBkgTag
double m_TCRawBkgTag[576][60]
Input Beambackground tag.
Definition TrgEclDigitizer.h:88
Belle2::TrgEclDigitizer::digitization01
void digitization01(std::vector< std::vector< double > > &, std::vector< std::vector< double > > &)
fit method, digi with 125ns interval
Definition TrgEclDigitizer.cc:317
Belle2::TrgEclDigitizer::m_SaveTCWaveForm
int m_SaveTCWaveForm
Flag of waveform table.
Definition TrgEclDigitizer.h:107
Belle2::TrgEclDigitizer::m_TCRawEnergy
double m_TCRawEnergy[576][60]
Input TC energy[GeV].
Definition TrgEclDigitizer.h:84
Belle2::TrgEclDigitizer::m_TCBkgContribution
double m_TCBkgContribution[576][80]
Beambackground contribution.
Definition TrgEclDigitizer.h:101
Belle2::TrgEclDigitizer::m_TCEnergy_tot
double m_TCEnergy_tot[576]
TC Energy converted from Xtarl Energy [GeV].
Definition TrgEclDigitizer.h:74
Belle2::TrgEclDigitizer::m_TCRawTiming
double m_TCRawTiming[576][60]
Input TC timing[ns].
Definition TrgEclDigitizer.h:86
Belle2::TrgEclDigitizer::interFADC
double interFADC(double)
Faster FADC using interpolation.
Definition TrgEclDigitizer.cc:609
Belle2::TrgEclDigitizer::m_TCSigContribution
double m_TCSigContribution[576][80]
Signal contribution.
Definition TrgEclDigitizer.h:103
Belle2::TrgEclDigitizer::m_DataBase
TrgEclDataBase * m_DataBase
Object of DataBase.
Definition TrgEclDigitizer.h:93
Belle2::TrgEclDigitizer::m_TCTiming_tot
double m_TCTiming_tot[576]
TC Timing converted from Xtarl Timing [GeV].
Definition TrgEclDigitizer.h:76
Belle2::TrgEclDigitizer::FADC
double FADC(int, double)
FADC.
Definition TrgEclDigitizer.cc:646
Belle2::TrgEclDigitizer::SimplifiedFADC
double SimplifiedFADC(int, double)
FADC.
Definition TrgEclDigitizer.cc:838
Belle2::TrgEclDigitizer::m_MatrixSerial
std::vector< std::vector< double > > m_MatrixSerial
Noise Low triangle Matrix of Serial noise.
Definition TrgEclDigitizer.h:99
Belle2::TrgEclDigitizer::m_WaveForm
double m_WaveForm[576][64]
TC Energy converted from Xtarl Energy [GeV].
Definition TrgEclDigitizer.h:109
Belle2::TrgEclDigitizer::u_max
double u_max(double, double)
Find max value between 2 vals;.
Definition TrgEclDigitizer.cc:1028
Belle2::TrgEclDigitizer::m_TCEnergy
double m_TCEnergy[576][80]
TC Energy converted from Xtarl Energy [GeV].
Definition TrgEclDigitizer.h:70
Belle2::TrgEclDigitizer::ShapeF
double ShapeF(double, double, double, double, double, double, double)
return shape using FADC function
Definition TrgEclDigitizer.cc:762
Belle2::TrgEclDigitizer::digitization02
void digitization02(std::vector< std::vector< double > > &, std::vector< std::vector< double > > &)
original no fit method, digi with 12ns interval
Definition TrgEclDigitizer.cc:478
Belle2::TrgEclDigitizer::m_TCBeambkgTag
int m_TCBeambkgTag[576][80]
Beambackground tag.
Definition TrgEclDigitizer.h:105
Belle2::TrgEclDigitizer::getTCHit
void getTCHit(int)
get TC Hits from Xtal hits
Definition TrgEclDigitizer.cc:92
Belle2::TrgEclDigitizer::m_FADC
int m_FADC
Flag of choosing the method of waveform generation function 0: use simplifiedFADC,...
Definition TrgEclDigitizer.h:111
Belle2::TrgEclDigitizer::m_TCMap
TrgEclMapping * m_TCMap
Object of TC Mapping.
Definition TrgEclDigitizer.h:91
Belle2::TrgEclDigitizer::m_BeambkgTag
int m_BeambkgTag
Flag of saving beam background tag or not.
Definition TrgEclDigitizer.h:113
Belle2::TrgEclDigitizer::m_TimeRange
double m_TimeRange
time range(default : -4000 ~ 4000 ns )
Definition TrgEclDigitizer.h:68
Belle2::TrgEclMapping
A class of TC Mapping.
Definition TrgEclMapping.h:26
Belle2::Unit::ns
static const double ns
Standard of [time].
Definition Unit.h:48
Belle2::Unit::GeV
static const double GeV
Standard of [energy, momentum, mass].
Definition Unit.h:51
Belle2::realNaN
static const double realNaN
This collects the B-meson properties in the hadronic B-decays It is used for the Ecms calibration in ...
Definition PhysicsObjectsMiraBelleEcmsBBModule.h:28
Belle2::sqrt
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Belle2
Abstract base class for different kinds of events.
Definition MillepedeAlgorithm.h:17
std
STL namespace.