| File: | trg/ecl/src/TrgEclDigitizer.cc |
| Warning: | line 710, column 9 Value stored to 'dt' is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 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 | |
| 26 | using namespace std; |
| 27 | using namespace Belle2; |
| 28 | //using namespace TRG; |
| 29 | // |
| 30 | // |
| 31 | // |
| 32 | TrgEclDigitizer::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 | // |
| 65 | TrgEclDigitizer::~TrgEclDigitizer() |
| 66 | { |
| 67 | delete m_TCMap; |
| 68 | delete m_DataBase; |
| 69 | } |
| 70 | // |
| 71 | // |
| 72 | // |
| 73 | |
| 74 | void |
| 75 | TrgEclDigitizer::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 | // |
| 91 | void |
| 92 | TrgEclDigitizer::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 | // |
| 316 | void |
| 317 | TrgEclDigitizer::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 | } |
| 477 | void |
| 478 | TrgEclDigitizer::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 | // |
| 566 | void |
| 567 | TrgEclDigitizer::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 | // |
| 608 | double |
| 609 | TrgEclDigitizer::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 | // |
| 645 | double |
| 646 | TrgEclDigitizer::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; |
Value stored to 'dt' is never read | |
| 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 | // |
| 761 | double |
| 762 | TrgEclDigitizer::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 | |
| 837 | double |
| 838 | TrgEclDigitizer::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 | |
| 958 | double 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 | // |
| 1027 | double |
| 1028 | TrgEclDigitizer::u_max(double aaa, double bbb) |
| 1029 | { |
| 1030 | if (aaa > bbb) { return aaa; } |
| 1031 | else { return bbb; } |
| 1032 | } |
| 1033 | // |
| 1034 | // |
| 1035 | // |