From 4f87dcc8257ec49b234e356e6586ef264bcba0be Mon Sep 17 00:00:00 2001 From: fengyvoid <67118022+fengyvoid@users.noreply.github.com> Date: Mon, 8 Sep 2025 10:51:05 -0500 Subject: [PATCH] LAPPD offset Fitting script for EBV2 --- offsetFit_MultipleLAPPD.cpp | 1315 +++++++++++++++++++++++++++++++++++ 1 file changed, 1315 insertions(+) create mode 100644 offsetFit_MultipleLAPPD.cpp diff --git a/offsetFit_MultipleLAPPD.cpp b/offsetFit_MultipleLAPPD.cpp new file mode 100644 index 000000000..0b7d64004 --- /dev/null +++ b/offsetFit_MultipleLAPPD.cpp @@ -0,0 +1,1315 @@ +// This is a script for offset fitting, use it on the output LAPPDTree.root +// This script require all LAPPD events appear in sequence, and all CTC events appear in sequence +/* procedure: + 1. load LAPPDTree.root (TODO: or use LAPPDTree_runNumber.root) + 2. load TimeStamp tree, get run number and part file number + 3. for each unique run number and part file number, do fit: + 4. in this file, for each LAPPD ID: + 5. loop TimeStamp tree + load data timestamp + load data beamgate + load PPS for board 0 and board 1 + 6. loop Trig (or GTrig) tree + load CTC PPS (word = 32) + load target trigger word + + 7. Find the number of resets in LAPPD PPS: + There must be at least one event in each reset. + Based on the event index order, only fit before the reset which doesn't have data event. + Save the data event and PPS index order. + + 8. Use the target trigger word, fit the offset + 9. save the offset for this LAPPD ID, run number, part file number, index, reset number. + 10. Print info to txt file + +Laser trigger word: 47 +Undelayed beam trigger word: 14 + +root -l -q 'offsetFit_MultipleLAPPD.cpp("LAPPDTree.root", 14, 1, 10, 0)' +root -l -q 'offsetFit_MultipleLAPPD.cpp("LAPPDTree.root", 47, 0, 10, 0)' +*/ +#include +#include +#include +#include +#include +#include +#include "TFile.h" +#include "TTree.h" +#include "TH2D.h" +#include "TH1D.h" +#include "TString.h" + +vector> fitInThisReset( + const std::vector &LAPPDDataTimeStampUL, + const std::vector &LAPPDDataBeamgateUL, + const std::vector &LAPPD_PPS, + const int fitTargetTriggerWord, + const std::vector &CTCTrigger, + const std::vector &CTCPPS, + const ULong64_t PPSDeltaT, + const int partFileNumber, + const int ACDCNumber) +{ + cout << "***************************************" << endl; + cout << "Fitting part file number: " << partFileNumber << ", ACDC number: " << ACDCNumber << endl; + cout << "Fitting in this reset with:" << endl; + cout << "LAPPDDataTimeStampUL size: " << LAPPDDataTimeStampUL.size() << endl; + cout << "LAPPDDataBeamgateUL size: " << LAPPDDataBeamgateUL.size() << endl; + cout << "LAPPD_PPS size: " << LAPPD_PPS.size() << endl; + cout << "CTCTrigger size: " << CTCTrigger.size() << endl; + cout << "CTCPPS size: " << CTCPPS.size() << endl; + cout << "PPSDeltaT: " << PPSDeltaT << endl; // in ps + // for this input PPS, fit an offset + // return the offset and other information in order + // precedure: + // 1. check drift + // 2. shift the timestmap and beamgate based on drift + // 3. fit the offset + + std::vector PPSInterval_ACDC; + + for (int i = 1; i < LAPPD_PPS.size(); i++) + { + + double diff = static_cast(LAPPD_PPS[i] - LAPPD_PPS[i - 1]); + cout< 2 microseconds + cout<<"PPSInterval_ACDC["< 2) + { + h->Fill(PPSInterval_ACDC[i]); // fill in microseconds + cout << "Fill histogram: PPSInterval_ACDC[" << i << "]: " << PPSInterval_ACDC[i] << endl; + } + } + + TF1 *gausf = new TF1("gausf", "gaus", 0, 1E3); + h->Fit(gausf, "Q"); // Q for quiet mode + ULong64_t drift = static_cast(gausf->GetParameter(1) * 1E3 * 1E3); // drift in ps + ULong64_t trueInterval = PPSDeltaT - drift; + std::cout << "Gaussian Drift in ps is " << drift << std::endl; + std::cout << "True PPS interval in ps is " << trueInterval << std::endl; + delete gausf; + delete h; + // initialize variables need for fitting + int orphanCount = 0; // count the number that how many data event timestamp doesn't matched to a target trigger within an interval. + std::map>> DerivationMap; + cout << "LAPPD_PPS.size() " << LAPPD_PPS.size() << " CTCPPS.size() " << CTCPPS.size() << endl; + for (int i = 0; i < LAPPD_PPS.size(); i++) + { + if (i > 5) + { + if (i % static_cast(LAPPD_PPS.size() / 5) == 0) + cout << "Fitting PPS " << i << " of " << LAPPD_PPS.size() << endl; + } + ULong64_t LAPPD_PPS_ns = LAPPD_PPS.at(i) / 1000; + ULong64_t LAPPD_PPS_truncated_ps = LAPPD_PPS.at(i) % 1000; + + for (int j = 0; j < CTCPPS.size(); j++) + { + + vector diffSum; + ULong64_t offsetNow_ns = 0; + if (drift == 0) + { + offsetNow_ns = CTCPPS.at(j) - LAPPD_PPS_ns; + } + else + { + double LAPPD_PPS_ns_dd = static_cast(LAPPD_PPS_ns); + double driftScaling = LAPPD_PPS_ns_dd / (trueInterval / 1000); + ULong64_t totalDriftedClock = static_cast(drift * driftScaling / 1000); + offsetNow_ns = CTCPPS.at(j) - (LAPPD_PPS_ns + totalDriftedClock); + if (i < 3 && j < 3) + cout<<"totalDriftedClock in ns = "< notOrphanIndex; + vector ctcPairedIndex; + vector ctcOrphanPairedIndex; + + for (int lappdb = 0; lappdb < LAPPDDataBeamgateUL.size(); lappdb++) + { + ULong64_t TS_ns = LAPPDDataTimeStampUL.at(lappdb) / 1000; + ULong64_t TS_truncated_ps = LAPPDDataTimeStampUL.at(lappdb) % 1000; + // if fit for the undelayed beam trigger, use BG + if (fitTargetTriggerWord == 14) + { + TS_ns = LAPPDDataBeamgateUL.at(lappdb) / 1000; + TS_truncated_ps = LAPPDDataBeamgateUL.at(lappdb) % 1000; + } + + ULong64_t driftCorrectionForTS = 0; + if (drift != 0) + { + double TS_ns_dd = static_cast(TS_ns); + double driftScaling = TS_ns_dd / (trueInterval / 1000); + driftCorrectionForTS = static_cast(drift * driftScaling / 1000); + if(lappdb<3 && i<3 && j<3) + cout<::max(); + bool useThisValue = true; + int minPairIndex = 0; + string reason = "none"; + Long64_t first_derivation = 0; + int minIndex = 0; + double useValue = true; + + // find the best match of target trigger to this TS, in ns level. + for (int ctcb = 0; ctcb < CTCTrigger.size(); ctcb++) + { + ULong64_t CTCTrigger_ns = CTCTrigger.at(ctcb); + Long64_t diff = CTCTrigger_ns - DriftCorrectedTS_ns; + if (diff < 0) + diff = -diff; + + if (diff < minMatchDiff) + { + minMatchDiff = diff; + minPairIndex = ctcb; + } + + Long64_t LastBound = diff - minMatchDiff; + if (LastBound < 0) + LastBound = -LastBound; + Long64_t FirstBound = diff - first_derivation; + if (FirstBound < 0) + FirstBound = -FirstBound; + + // save the dt for the first matching to determin the matching position in range. + if (ctcb == 0) + first_derivation = diff; + + if (ctcb == CTCTrigger.size() - 1 && LastBound / 1E9 < 0.01) + { + useValue = false; + reason = "When matching the last TS, this diff is too close to (or even is) the minMatchDiff at index " + std::to_string(minPairIndex) + " in " + std::to_string(CTCTrigger.size()) + ". All LAPPD TS is out of the end of CTC trigger range"; + } + if (ctcb == CTCTrigger.size() - 1 && FirstBound / 1E9 < 0.01) + { + useValue = false; + reason = "When matching the last TS, this diff is too close to (or even is) the first_derivation at index " + std::to_string(minPairIndex) + " in " + std::to_string(CTCTrigger.size()) + ". All LAPPD TS is out of the beginning of CTC trigger range"; + } + } + + if (useValue) + { + diffSum.push_back(minMatchDiff); + //cout<<"LAPPD PPS "< maxAllowedDiff || minMatchDiff < minAllowedDiff) + { + // TODO: adjust the limit for laser. + orphanCount += 1; + ctcOrphanPairedIndex.push_back(minPairIndex); + } + else + { + notOrphanIndex.push_back(lappdb); + ctcPairedIndex.push_back(minPairIndex); + } + } + else + { + orphanCount += 1; + } + } + + double mean_dev = 0; + for (int k = 0; k < diffSum.size(); k++) + { + mean_dev += diffSum[k]; + } + if (diffSum.size() > 0) + mean_dev = mean_dev / diffSum.size(); + + + double mean_dev_noOrphan = 0; + double pairedCount = 0; + for (int k = 0; k < notOrphanIndex.size(); k++) + { + if (fitTargetTriggerWord == 14) + { + if (diffSum.at(k) > 322E3 && diffSum.at(k) < 326E3) + { + mean_dev_noOrphan += diffSum[notOrphanIndex[k]]; + pairedCount += 1; + } + } + else + { + mean_dev_noOrphan += diffSum[notOrphanIndex[k]]; + } + } + if ((diffSum.size() - orphanCount) != 0) + { + mean_dev_noOrphan = mean_dev_noOrphan / (diffSum.size() - orphanCount); + } + else + { + mean_dev_noOrphan = -1; // if all timestamps are out of the range, set it to -1 + } + + bool increMean_dev = false; + int maxAttempts = 1000; + int attemptCount = 0; + + // the mean_dev can't be larger than 1 s because the beam spill is 15Hz. + // so, a good match will give a large number of events in desired range on integer level, minus the mean _dev. + // by selecting the largest quality Number, for matchs with the same matched beamgate, smaller mean_dev is better. + double qualityNumber = pairedCount*1e9 - mean_dev; + bool debug = false; + if (debug) + { + cout<<"LAPPD PPS "< 0) + { + int increament_dev = 0; + while (true) + { + // TODO + auto iter = DerivationMap.find(mean_dev); + if (iter == DerivationMap.end() || iter->second.empty()) + { + vector Info = {i, j, orphanCount, static_cast(mean_dev_noOrphan * 1000), increament_dev}; + /* + DerivationMap[mean_dev].push_back(Info); + DerivationMap[mean_dev].push_back(notOrphanIndex); + DerivationMap[mean_dev].push_back(ctcPairedIndex); + DerivationMap[mean_dev].push_back(ctcOrphanPairedIndex); + */ + + DerivationMap[qualityNumber].push_back(Info); + DerivationMap[qualityNumber].push_back(notOrphanIndex); + DerivationMap[qualityNumber].push_back(ctcPairedIndex); + DerivationMap[qualityNumber].push_back(ctcOrphanPairedIndex); + + break; + } + else + { + increament_dev += 1; + mean_dev += 0.001; // if the mean_dev is already in the map, increase it by 1ps + qualityNumber += 1e-6; + attemptCount += 1; + if (attemptCount > maxAttempts) + break; + } + } + } + } + } + // finish matching, found the minimum mean_dev in the map, extract the matching information + double min_mean_dev = std::numeric_limits::max(); + double max_qualityNumber = 0; + int final_i = 0; + int final_j = 0; + int gotOrphanCount = 0; + double gotMin_mean_dev_noOrphan = 0; + double increament_times = 0; + vector final_notOrphanIndex; + vector final_ctcPairedIndex; + vector final_ctcOrphanIndex; + for (const auto &minIter : DerivationMap) + { + //if (minIter.first > 10 && minIter.first < min_mean_dev) + if (minIter.first > -1 && minIter.first > max_qualityNumber) + { + //min_mean_dev = minIter.first; + max_qualityNumber = minIter.first; + double qnum_ns = minIter.first / 1e9; + min_mean_dev = (std::ceil(qnum_ns) - qnum_ns)*1e9; + final_i = minIter.second[0][0]; + final_j = minIter.second[0][1]; + gotOrphanCount = minIter.second[0][2]; + gotMin_mean_dev_noOrphan = static_cast(minIter.second[0][3]) / 1000; + increament_times = minIter.second[0][4]; + final_notOrphanIndex = minIter.second[1]; + final_ctcPairedIndex = minIter.second[2]; + final_ctcOrphanIndex = minIter.second[3]; + } + } + + ULong64_t final_offset_ns = 0; + ULong64_t final_offset_ps_negative = 0; + if (drift == 0) + { + cout<<"Drift is 0 above ns level." << " Using CTC PPS at index " << final_j <<" = "< TimeStampRaw; + vector BeamGateRaw; + vector TimeStamp_ns; + vector BeamGate_ns; + vector TimeStamp_ps; + vector BeamGate_ps; + vector EventIndex; + vector EventDeviation_ns; + vector CTCTriggerIndex; + vector CTCTriggerTimeStamp_ns; + vector BeamGate_correction_tick; + vector TimeStamp_correction_tick; + vector PPS_tick_correction; // if timestamp falls in between of PPS i and i+1, corrected timestamp = timestamp + PPS_tick_correction[i] + vector LAPPD_PPS_missing_ticks; + vector LAPPD_PPS_interval_ticks; + + vector BG_PPSBefore; + vector BG_PPSAfter; + vector BG_PPSDiff; + vector BG_PPSMiss; + vector TS_PPSBefore; + vector TS_PPSAfter; + vector TS_PPSDiff; + vector TS_PPSMiss; + + vector TS_driftCorrection_ns; + vector BG_driftCorrection_ns; + + // calculate the missing ticks for each LAPPD PPS + PPS_tick_correction.push_back(0); + LAPPD_PPS_missing_ticks.push_back(0); + LAPPD_PPS_interval_ticks.push_back(0); + ULong64_t intervalTicks = 320000000 * (PPSDeltaT / 1000 / 1E9); + for (int i = 1; i < LAPPD_PPS.size(); i++) + { + ULong64_t thisPPS = LAPPD_PPS.at(i) / 3125; + ULong64_t prevPPS = LAPPD_PPS.at(i - 1) / 3125; + long long thisInterval = thisPPS - prevPPS; + long long thisMissingTicks = intervalTicks - thisInterval; + LAPPD_PPS_interval_ticks.push_back(thisInterval); + + // if the difference between this PPS and previous PPS is > intervalTicks-20 and < intervalTicks+5, + // the missed value is thisInterval - intervalTicks + // push this value plus the sum of previous missing ticks to the vector + // else just push the sum of previous missing ticks + long long sumOfPreviousMissingTicks = 0; + for (int j = 0; j < LAPPD_PPS_missing_ticks.size(); j++) + { + sumOfPreviousMissingTicks += LAPPD_PPS_missing_ticks.at(j); + } + + if (thisMissingTicks > -20 && thisMissingTicks < 20) + { + LAPPD_PPS_missing_ticks.push_back(thisMissingTicks); + PPS_tick_correction.push_back(thisMissingTicks + sumOfPreviousMissingTicks); + } + else + { + // some time one pps might be wired, but the combination of two is ok. + bool combined = false; + // if there are missing PPS, interval is like 22399999990 % 3200000000 = 3199999990 + long long pInterval = thisInterval % intervalTicks; + long long missingPTicks = 0; + if (intervalTicks - pInterval < 30) + missingPTicks = intervalTicks - pInterval; + else if (pInterval < 30) + missingPTicks = -pInterval; + ////// + if (missingPTicks == 1 || missingPTicks == -1) + { + cout << "Found missing tick is " << missingPTicks << ",continue." << endl; + LAPPD_PPS_missing_ticks.push_back(0); + PPS_tick_correction.push_back(sumOfPreviousMissingTicks); + continue; + } + + if (missingPTicks != 0) + { + LAPPD_PPS_missing_ticks.push_back(missingPTicks); + PPS_tick_correction.push_back(missingPTicks + sumOfPreviousMissingTicks); + combined = true; + cout << "Pushing PPS correction " << i << ", this PPS interval tick is " << thisInterval << ", missing ticks: " << thisMissingTicks << ", push missing " << LAPPD_PPS_missing_ticks.at(i) << ", push correction " << PPS_tick_correction.at(i) << endl; + continue; + } + + // if one PPS is not recoreded correctly, like one interval is 1574262436, followed by a 4825737559 + // then 1574262436 + 4825737559 = 6399999995 = 3200000000 * 2 - 5 + long long nextInterval = 0; + if (i < LAPPD_PPS.size() - 1) + { + nextInterval = LAPPD_PPS.at(i + 1) / 3125 - thisPPS; + long long combinedMissingTicks = intervalTicks - (nextInterval + thisInterval) % intervalTicks; + if (combinedMissingTicks > -30 && combinedMissingTicks < 30) + { + LAPPD_PPS_missing_ticks.push_back(combinedMissingTicks); + PPS_tick_correction.push_back(combinedMissingTicks + sumOfPreviousMissingTicks); + combined = true; + cout << "Pushing PPS correction " << i << ", this PPS interval tick is " << thisInterval << ", missing ticks: " << thisMissingTicks << ", push missing " << LAPPD_PPS_missing_ticks.at(i) << ", push correction " << PPS_tick_correction.at(i) << endl; + continue; + } + } + + if (!combined) + { + LAPPD_PPS_missing_ticks.push_back(0); + PPS_tick_correction.push_back(sumOfPreviousMissingTicks); + } + } + cout << "Pushing PPS correction " << i << ", this PPS interval tick is " << thisInterval << ", missing ticks: " << thisMissingTicks << ", push missing " << LAPPD_PPS_missing_ticks.at(i) << ", push correction " << PPS_tick_correction.at(i) << endl; + } + + // loop all data events, plus the offset, save the event time and beamgate time + // fing the closest CTC trigger, also save all information + cout << "Start saving results. PPS size: " << LAPPD_PPS.size() << ", beamgate size: " << LAPPDDataBeamgateUL.size() << endl; + cout << "First PPS: " << LAPPD_PPS.at(0) / 3125 << ", Last PPS: " << LAPPD_PPS.at(LAPPD_PPS.size() - 1) / 3125 << endl; + cout << "First BG: " << LAPPDDataBeamgateUL.at(0) / 3125 << ", Last BG: " << LAPPDDataBeamgateUL.at(LAPPDDataBeamgateUL.size() - 1) / 3125 << endl; + cout << "First TS: " << LAPPDDataTimeStampUL.at(0) / 3125 << ", Last TS: " << LAPPDDataTimeStampUL.at(LAPPDDataTimeStampUL.size() - 1) / 3125 << endl; + + for (int l = 0; l < LAPPDDataTimeStampUL.size(); l++) + { + ULong64_t TS_ns = LAPPDDataTimeStampUL.at(l) / 1000; + ULong64_t TS_truncated_ps = LAPPDDataTimeStampUL.at(l) % 1000; + ULong64_t BG_ns = LAPPDDataBeamgateUL.at(l) / 1000; + ULong64_t BG_truncated_ps = LAPPDDataBeamgateUL.at(l) % 1000; + ULong64_t driftCorrectionForTS = 0; + ULong64_t driftCorrectionForBG = 0; + if (drift != 0) + { + + double driftScaling = static_cast(TS_ns) / (trueInterval / 1000); + driftCorrectionForTS = static_cast(drift/ 1000 * driftScaling ); + double driftScalingBG = static_cast(BG_ns) / (trueInterval / 1000); + driftCorrectionForBG = static_cast(drift/ 1000 * driftScalingBG ); + cout<<"drift = "< LAPPD_PPS.at(LAPPD_PPS.size() - 1) / 3125) + { + BeamGate_correction_tick.push_back(PPS_tick_correction.at(LAPPD_PPS.size() - 1) + 1000); + cout << "BGraw = " << LAPPDDataBeamgateUL.at(l) / 3125 << ", pps = " << LAPPD_PPS.at(LAPPD_PPS.size() - 1) << ", pps/3125 = " << LAPPD_PPS.at(LAPPD_PPS.size() - 1) / 3125 << endl; + BGFound = true; + } + if (!BGFound && LAPPDDataBeamgateUL.at(l) / 3125 < LAPPD_PPS.at(0) / 3125) + { + BeamGate_correction_tick.push_back(0 + 1000); + cout << "BGraw less than pps0" << endl; + BGFound = true; + } + if (!TSFound || !BGFound) + { + cout << "Error: PPS not found for event " << l << ", TSFound: " << TSFound << ", BGFound: " << BGFound << endl; + } + /* + cout << "******Found result:" << endl; + cout << "LAPPDDataTimeStampUL.at(" << l << "): " << LAPPDDataTimeStampUL.at(l)/3125 << endl; + cout << "LAPPDDataBeamgateUL.at(" << l << "): " << LAPPDDataBeamgateUL.at(l)/3125 << endl; + cout << "TS_ns: " << TS_ns << endl; + cout << "BG_ns: " << BG_ns << endl; + cout << "final_offset_ns: " << final_offset_ns << endl; + cout << "drift: " << drift << endl; + cout << "TS driftscaling: " << TS_ns / trueInterval / 1000 << endl; + cout << "DriftCorrectedTS_ns: " << DriftCorrectedTS_ns << endl; + cout << "TS_truncated_ps: " << TS_truncated_ps << endl; + cout << "DriftCorrectedBG_ns: " << DriftCorrectedBG_ns << endl; + cout << "BG_truncated_ps: " << BG_truncated_ps << endl; + cout << "Found min_mean_dev_match: " << min_mean_dev_match << endl; + cout << "MatchedIndex: " << matchedIndex << endl; + cout << "CTCTrigger.at(" << matchedIndex << "): " << CTCTrigger.at(matchedIndex) << endl; + cout << "BeamGate_correction_tick at " << l << ": " << BeamGate_correction_tick.at(l) << endl; + cout << "TimeStamp_correction_tick at " << l << ": " << TimeStamp_correction_tick.at(l) << endl; + */ + + // for this LAPPDDataBeamgateUL.at(l), in the LAPPD_PPS vector, find it's closest PPS before and after, and also calculate the time difference between them, and the missing tick between them. + // save as BG_PPSBefore_tick, BG_PPSAfter_tick, BG_PPSDiff_tick, BG_PPSMissing_tick + // Do the same thing for LAPPDDataTimeStampUL.at(l), save as TS_PPSBefore_tick, TS_PPSAfter_tick, TS_PPSDiff_tick, TS_PPSMissing_tick + + ULong64_t BG_PPSBefore_tick = 0; + ULong64_t BG_PPSAfter_tick = 0; + ULong64_t BG_PPSDiff_tick = 0; + ULong64_t BG_PPSMissing_tick = 0; + + ULong64_t TS_PPSBefore_tick = 0; + ULong64_t TS_PPSAfter_tick = 0; + ULong64_t TS_PPSDiff_tick = 0; + ULong64_t TS_PPSMissing_tick = 0; + // if the first PPS is later than beamgate, then set before = 0, after is the first, diff is the first, missing is the first + + cout << "Start finding PPS before and after the beamgate" << endl; + + if ((LAPPD_PPS.at(0) / 3125 > LAPPDDataBeamgateUL.at(l) / 3125) || (LAPPD_PPS.at(LAPPD_PPS.size() - 1) / 3125 < LAPPDDataBeamgateUL.at(l) / 3125)) + { + if (LAPPD_PPS.at(0) / 3125 > LAPPDDataBeamgateUL.at(l) / 3125) + { + BG_PPSBefore_tick = 0; + BG_PPSAfter_tick = LAPPD_PPS.at(0) / 3125; + BG_PPSDiff_tick = LAPPD_PPS.at(0) / 3125; + BG_PPSMissing_tick = LAPPD_PPS.at(0) / 3125; + cout << "First PPS is later than beamgate, before is 0, after is the first, diff is the first, missing is the first" << endl; + } + // if the last PPS is earlier than beamgate, then set before is the last, after = 0, diff is the last, missing is the last + if (LAPPD_PPS.at(LAPPD_PPS.size() - 1) / 3125 < LAPPDDataBeamgateUL.at(l) / 3125) + { + BG_PPSBefore_tick = LAPPD_PPS.at(LAPPD_PPS.size() - 1) / 3125; + BG_PPSAfter_tick = 0; + BG_PPSDiff_tick = LAPPD_PPS.at(LAPPD_PPS.size() - 1) / 3125; + BG_PPSMissing_tick = LAPPD_PPS.at(LAPPD_PPS.size() - 1) / 3125; + cout << "Last PPS is earlier than beamgate, before is the last, after is 0, diff is the last, missing is the last" << endl; + } + } + else + { + // if the first PPS is earlier than beamgate, and the last PPS is later than beamgate, then find the closest PPS before and after + for (int i = 0; i < LAPPD_PPS.size() - 1; i++) + { + if (LAPPD_PPS.at(i) / 3125 < LAPPDDataBeamgateUL.at(l) / 3125 && LAPPD_PPS.at(i + 1) / 3125 > LAPPDDataBeamgateUL.at(l) / 3125) + { + BG_PPSBefore_tick = LAPPD_PPS.at(i) / 3125; + BG_PPSAfter_tick = LAPPD_PPS.at(i + 1) / 3125; + BG_PPSDiff_tick = LAPPD_PPS.at(i + 1) / 3125 - LAPPD_PPS.at(i) / 3125; + ULong64_t DiffTick = (LAPPD_PPS.at(i + 1) - LAPPD_PPS.at(i)) / 3125 - 1000; + BG_PPSMissing_tick = 3200000000 - DiffTick; + cout << "Found PPS before and after the beamgate, before: " << BG_PPSBefore_tick << ", after: " << BG_PPSAfter_tick << ", diff: " << BG_PPSDiff_tick << ", missing: " << BG_PPSMissing_tick << endl; + break; + } + } + } + + // do the samething for timestamp + cout << "Start finding PPS before and after the timestamp" << endl; + + if ((LAPPD_PPS.at(0) / 3125 > LAPPDDataTimeStampUL.at(l) / 3125) || (LAPPD_PPS.at(LAPPD_PPS.size() - 1) / 3125 < LAPPDDataTimeStampUL.at(l) / 3125)) + { + if (LAPPD_PPS.at(0) / 3125 > LAPPDDataTimeStampUL.at(l) / 3125) + { + TS_PPSBefore_tick = 0; + TS_PPSAfter_tick = LAPPD_PPS.at(0) / 3125; + TS_PPSDiff_tick = LAPPD_PPS.at(0) / 3125; + TS_PPSMissing_tick = LAPPD_PPS.at(0) / 3125; + cout << "First PPS is later than timestamp, before is 0, after is the first, diff is the first, missing is the first" << endl; + } + if (LAPPD_PPS.at(LAPPD_PPS.size() - 1) / 3125 < LAPPDDataTimeStampUL.at(l) / 3125) + { + TS_PPSBefore_tick = LAPPD_PPS.at(LAPPD_PPS.size() - 1) / 3125; + TS_PPSAfter_tick = 0; + TS_PPSDiff_tick = LAPPD_PPS.at(LAPPD_PPS.size() - 1) / 3125; + TS_PPSMissing_tick = LAPPD_PPS.at(LAPPD_PPS.size() - 1) / 3125; + cout << "Last PPS is earlier than timestamp, before is the last, after is 0, diff is the last, missing is the last" << endl; + } + } + else + { + for (int i = 0; i < LAPPD_PPS.size() - 1; i++) + { + if (LAPPD_PPS.at(i) / 3125 < LAPPDDataTimeStampUL.at(l) / 3125 && LAPPD_PPS.at(i + 1) / 3125 > LAPPDDataTimeStampUL.at(l) / 3125) + { + TS_PPSBefore_tick = LAPPD_PPS.at(i) / 3125; + TS_PPSAfter_tick = LAPPD_PPS.at(i + 1) / 3125; + TS_PPSDiff_tick = LAPPD_PPS.at(i + 1) / 3125 - LAPPD_PPS.at(i) / 3125; + ULong64_t DiffTick = (LAPPD_PPS.at(i + 1) - LAPPD_PPS.at(i)) / 3125 - 1000; + TS_PPSMissing_tick = 3200000000 - DiffTick; + cout << "Found PPS before and after the timestamp, before: " << TS_PPSBefore_tick << ", after: " << TS_PPSAfter_tick << ", diff: " << TS_PPSDiff_tick << ", missing: " << TS_PPSMissing_tick << endl; + break; + } + } + } + + TimeStampRaw.push_back(LAPPDDataTimeStampUL.at(l) / 3125); + BeamGateRaw.push_back(LAPPDDataBeamgateUL.at(l) / 3125); + TimeStamp_ns.push_back(DriftCorrectedTS_ns); + BeamGate_ns.push_back(DriftCorrectedBG_ns); + TimeStamp_ps.push_back(TS_truncated_ps); + BeamGate_ps.push_back(BG_truncated_ps); + EventIndex.push_back(l); + EventDeviation_ns.push_back(min_mean_dev_match); + // to get ps, should minus the TS_truncated_ps + CTCTriggerIndex.push_back(matchedIndex); + CTCTriggerTimeStamp_ns.push_back(CTCTrigger.at(matchedIndex)); + + BG_PPSBefore.push_back(BG_PPSBefore_tick); + BG_PPSAfter.push_back(BG_PPSAfter_tick); + BG_PPSDiff.push_back(BG_PPSDiff_tick); + BG_PPSMiss.push_back(BG_PPSMissing_tick); + TS_PPSBefore.push_back(TS_PPSBefore_tick); + TS_PPSAfter.push_back(TS_PPSAfter_tick); + TS_PPSDiff.push_back(TS_PPSDiff_tick); + TS_PPSMiss.push_back(TS_PPSMissing_tick); + + TS_driftCorrection_ns.push_back(driftCorrectionForTS); + BG_driftCorrection_ns.push_back(driftCorrectionForBG); + } + ULong64_t totalEventNumber = LAPPDDataTimeStampUL.size(); + ULong64_t gotOrphanCount_out = gotOrphanCount; + ULong64_t gotMin_mean_dev_noOrphan_out = gotMin_mean_dev_noOrphan; + ULong64_t increament_times_out = increament_times; + ULong64_t min_mean_dev_out = min_mean_dev; + ULong64_t final_i_out = final_i; + ULong64_t final_j_out = final_j; + ULong64_t drift_out = drift; + + vector FitInfo = {final_offset_ns, final_offset_ps_negative, gotOrphanCount_out, gotMin_mean_dev_noOrphan_out, increament_times_out, min_mean_dev_out, final_i_out, final_j_out, totalEventNumber, drift_out}; + + vector> Result = {FitInfo, TimeStampRaw, BeamGateRaw, TimeStamp_ns, BeamGate_ns, TimeStamp_ps, BeamGate_ps, EventIndex, EventDeviation_ns, CTCTriggerIndex, CTCTriggerTimeStamp_ns, BeamGate_correction_tick, TimeStamp_correction_tick, LAPPD_PPS_interval_ticks, BG_PPSBefore, BG_PPSAfter, BG_PPSDiff, BG_PPSMiss, TS_PPSBefore, TS_PPSAfter, TS_PPSDiff, TS_PPSMiss, TS_driftCorrection_ns, BG_driftCorrection_ns}; + return Result; +} + +vector> fitInPartFile(TTree *lappdTree, TTree *triggerTree, int runNumber, int subRunNumber, int partFileNumber, int LAPPD_ID, int fitTargetTriggerWord, bool triggerGrouped, int intervalInSecond) +{ + cout << "***************************************" << endl; + cout << "Fitting in run " << runNumber << ", sub run " << subRunNumber << ", part file " << partFileNumber << " for LAPPD ID " << LAPPD_ID << endl; + + vector LAPPD_PPS0; + vector LAPPD_PPS1; + vector LAPPDDataTimeStampUL; + vector LAPPDDataBeamgateUL; + + vector CTCTargetTimeStamp; + vector CTCPPSTimeStamp; + + int LAPPD_ID_inTree; + int runNumber_inTree; + int subRunNumber_inTree; + int partFileNumber_inTree; + + ULong64_t ppsTime0; + ULong64_t ppsTime1; + ULong64_t LAPPDTimeStampUL; + ULong64_t LAPPDBeamgateUL; + ULong64_t LAPPDDataTimestamp; + ULong64_t LAPPDDataBeamgate; + double LAPPDDataTimestampFloat; + double LAPPDDataBeamgateFloat; + + vector *CTCTriggerWord = nullptr; + ULong64_t CTCTimeStamp; + vector *groupedTriggerWords = nullptr; + vector *groupedTriggerTimestamps = nullptr; + + lappdTree->SetBranchAddress("RunNumber", &runNumber_inTree); + lappdTree->SetBranchAddress("SubRunNumber", &subRunNumber_inTree); + lappdTree->SetBranchAddress("PartFileNumber", &partFileNumber_inTree); + lappdTree->SetBranchAddress("LAPPD_ID", &LAPPD_ID_inTree); + lappdTree->SetBranchAddress("LAPPDDataTimeStampUL", &LAPPDTimeStampUL); + lappdTree->SetBranchAddress("LAPPDDataBeamgateUL", &LAPPDBeamgateUL); + lappdTree->SetBranchAddress("LAPPDDataTimestamp", &LAPPDDataTimestamp); + lappdTree->SetBranchAddress("LAPPDDataBeamgate", &LAPPDDataBeamgate); + lappdTree->SetBranchAddress("LAPPDDataTimestampFloat", &LAPPDDataTimestampFloat); + lappdTree->SetBranchAddress("LAPPDDataBeamgateFloat", &LAPPDDataBeamgateFloat); + lappdTree->SetBranchAddress("ppsTime0", &ppsTime0); + lappdTree->SetBranchAddress("ppsTime1", &ppsTime1); + + // triggerTree->Print(); + triggerTree->SetBranchAddress("RunNumber", &runNumber_inTree); + triggerTree->SetBranchAddress("SubRunNumber", &subRunNumber_inTree); + triggerTree->SetBranchAddress("PartFileNumber", &partFileNumber_inTree); + if (triggerGrouped) + { + triggerTree->SetBranchAddress("gTrigWord", &groupedTriggerWords); + triggerTree->SetBranchAddress("gTrigTime", &groupedTriggerTimestamps); + } + else + { + triggerTree->SetBranchAddress("CTCTriggerWord", &CTCTriggerWord); + triggerTree->SetBranchAddress("CTCTimeStamp", &CTCTimeStamp); + } + + int l_nEntries = lappdTree->GetEntries(); + int t_nEntries = triggerTree->GetEntries(); + int repeatedPPSNumber0 = -1; + int repeatedPPSNumber1 = -1; + + // 5. loop TimeStamp tree + for (int i = 0; i < l_nEntries; i++) + { + lappdTree->GetEntry(i); + // if(i<100) + // cout<<"LAPPD_ID_inTree: "<(LAPPDTimeStampUL * 1000 / 8 * 25) / 1E9 / 1000 << ", Beamgate: " << static_cast(LAPPDBeamgateUL * 1000 / 8 * 25) / 1E9 / 1000 << endl; + LAPPDDataTimeStampUL.push_back(LAPPDTimeStampUL * 1000 / 8 * 25); + LAPPDDataBeamgateUL.push_back(LAPPDBeamgateUL * 1000 / 8 * 25); + /* + cout << "LAPPDDataTimestamp: " << LAPPDDataTimestamp << ", LAPPDDataBeamgate: " << LAPPDDataBeamgate << endl; + cout << "LAPPDDataTimestampFloat: " << LAPPDDataTimestampFloat << ", LAPPDDataBeamgateFloat: " << LAPPDDataBeamgateFloat << endl; + ULong64_t ULTS = LAPPDTimeStampUL*1000/8*25; + double DTS = static_cast(LAPPDDataTimestamp); + double plusDTS = DTS + LAPPDDataTimestampFloat; + ULong64_t ULBG = LAPPDBeamgateUL*1000/8*25; + double DBG = static_cast(LAPPDDataBeamgate); + double plusDBG = DBG + LAPPDDataBeamgateFloat; + cout << std::fixed << " ULTS: " << ULTS << ", DTS: " << DTS << ", plusDTS: " << plusDTS << endl; + cout << std::fixed << " ULBG: " << ULBG << ", DBG: " << DBG << ", plusDBG: " << plusDBG << endl; + cout << endl; + */ + } + else if (LAPPDTimeStampUL == 0) + { + // cout << "ppsTime0: " << ppsTime0 << ", ppsTime1: " << ppsTime1 << endl; + if (LAPPD_PPS0.size() == 0) + LAPPD_PPS0.push_back(ppsTime0 * 1000 / 8 * 25); + if (LAPPD_PPS1.size() == 0) + LAPPD_PPS1.push_back(ppsTime1 * 1000 / 8 * 25); + if (LAPPD_PPS0.size() > 0 && ppsTime0 * 1000 / 8 * 25 != LAPPD_PPS0.at(LAPPD_PPS0.size() - 1)) + LAPPD_PPS0.push_back(ppsTime0 * 1000 / 8 * 25); + else + repeatedPPSNumber0 += 1; + if (LAPPD_PPS1.size() > 0 && ppsTime1 * 1000 / 8 * 25 != LAPPD_PPS1.at(LAPPD_PPS1.size() - 1)) + LAPPD_PPS1.push_back(ppsTime1 * 1000 / 8 * 25); + else + repeatedPPSNumber1 += 1; + } + } + } + cout << "repeatedPPSNumber0: " << repeatedPPSNumber0 << ", loaded PPS0 size: " << LAPPD_PPS0.size() << endl; + cout << "repeatedPPSNumber1: " << repeatedPPSNumber1 << ", loaded PPS1 size: " << LAPPD_PPS1.size() << endl; + + // 6. loop Trig (or GTrig) tree + for (int i = 0; i < t_nEntries; i++) + { + triggerTree->GetEntry(i); + if (runNumber_inTree == runNumber && subRunNumber_inTree == subRunNumber && partFileNumber_inTree == partFileNumber) + { + // cout<<"triggerGrouped: "<size(); j++) + { + // cout<<"At j = "<at(j)<<", fitTargetTriggerWord: "<at(j) == fitTargetTriggerWord) + CTCTargetTimeStamp.push_back(groupedTriggerTimestamps->at(j)); + if (groupedTriggerWords->at(j) == 32) + CTCPPSTimeStamp.push_back(groupedTriggerTimestamps->at(j)); + } + } + else + { + for (int j = 0; j < CTCTriggerWord->size(); j++) + { + if (CTCTriggerWord->at(j) == fitTargetTriggerWord) + CTCTargetTimeStamp.push_back(CTCTimeStamp); + if (CTCTriggerWord->at(j) == 32) + CTCPPSTimeStamp.push_back(CTCTimeStamp); + } + } + } + } + cout << "Vector for partfile " << partFileNumber << " for LAPPD ID " << LAPPD_ID << " loaded." << endl; + cout << "LAPPDDataTimeStampUL in ps size: " << LAPPDDataTimeStampUL.size() << endl; + cout << "LAPPDDataBeamgateUL in ps size: " << LAPPDDataBeamgateUL.size() << endl; + cout << "LAPPD_PPS0 size: " << LAPPD_PPS0.size() << endl; + cout << "LAPPD_PPS1 size: " << LAPPD_PPS1.size() << endl; + cout << "CTCTargetTimeStamp size: " << CTCTargetTimeStamp.size() << endl; + cout << "CTCPPSTimeStamp size: " << CTCPPSTimeStamp.size() << endl; + // 7. Find the number of resets in LAPPD PPS: + int LAPPDDataFitStopIndex = 0; + int resetNumber = 0; + // first check is there a reset, if no, set the LAPPDDataFitStopIndex to the size of LAPPDDataTimeStampUL + // if all PPS in this part file was increamented, then there is no reset + for (int i = 1; i < LAPPD_PPS0.size(); i++) + { + if (LAPPD_PPS0[i] < LAPPD_PPS0[i - 1]) + { + resetNumber += 1; + cout << "For LAPPD ID " << LAPPD_ID << ", run number " << runNumber << ", sub run number " << subRunNumber << ", part file number " << partFileNumber << ", reset " + << " found at PPS_ACDC0 index " << i << endl; + break; + } + } + for (int i = 1; i < LAPPD_PPS1.size(); i++) + { + if (LAPPD_PPS1[i] < LAPPD_PPS1[i - 1]) + { + resetNumber += 1; + cout << "For LAPPD ID " << LAPPD_ID << ", run number " << runNumber << ", sub run number " << subRunNumber << ", part file number " << partFileNumber << ", reset " + << " found at PPS_ACDC1 index " << i << endl; + break; + } + } + if (resetNumber == 0) + cout << "For LAPPD ID " << LAPPD_ID << ", run number " << runNumber << ", sub run number " << subRunNumber << ", part file number " << partFileNumber << ", no reset found." << endl; + + // if reset found, loop the timestamp in order to find the LAPPDDataFitStopIndex + // + if (resetNumber != 0) + { + for (int i = 1; i < LAPPDDataTimeStampUL.size(); i++) + { + if (LAPPDDataTimeStampUL[i] < LAPPDDataTimeStampUL[i - 1]) + { + LAPPDDataFitStopIndex = i - 1; + break; + } + } + // TODO: extend this to later offsets + } + // 8. Use the target trigger word, fit the offset + // TODO: fit for each reset + vector> ResultTotal; + if (LAPPDDataTimeStampUL.size() == LAPPDDataBeamgateUL.size()) + { + if (LAPPD_PPS0.size() == 0 || LAPPD_PPS1.size() == 0) + { + cout << "Error: PPS0 or PPS1 is empty, return empty result." << endl; + return ResultTotal; + } + vector> ResultACDC0 = fitInThisReset(LAPPDDataTimeStampUL, LAPPDDataBeamgateUL, LAPPD_PPS0, fitTargetTriggerWord, CTCTargetTimeStamp, CTCPPSTimeStamp, intervalInSecond * 1E9 * 1000, partFileNumber, 0); + vector> ResultACDC1 = fitInThisReset(LAPPDDataTimeStampUL, LAPPDDataBeamgateUL, LAPPD_PPS1, fitTargetTriggerWord, CTCTargetTimeStamp, CTCPPSTimeStamp, intervalInSecond * 1E9 * 1000, partFileNumber, 1); + + // 9. save the offset for this LAPPD ID, run number, part file number, index, reset number. + + cout << "Fitting in part file " << partFileNumber << " for LAPPD ID " << LAPPD_ID << " done." << endl; + + // Combine ResultACDC0 and ResultACDC1 to ResultTotal + for (int i = 0; i < ResultACDC0.size(); i++) + { + ResultTotal.push_back(ResultACDC0[i]); + } + for (int i = 0; i < ResultACDC1.size(); i++) + { + ResultTotal.push_back(ResultACDC1[i]); + } + } + return ResultTotal; +} + +void offsetFit_MultipleLAPPD(string fileName, int fitTargetTriggerWord, bool triggerGrouped, int intervalInSecond, int processPFNumber) +{ + // 1. load LAPPDTree.root + const string file = fileName; + TFile *f = new TFile(file.c_str(), "READ"); + if (!f->IsOpen()) + { + std::cerr << "Error: cannot open file " << file << std::endl; + return; + } + + std::cout << "Opened file " << file << std::endl; + + std::ofstream outputOffset("offset.txt"); + + outputOffset << "runNumber" + << "\t" + << "subRunNumber" + << "\t" + << "partFileNumber" + << "\t" + << "resetNumber" + << "\t" + << "LAPPD_ID" + << "\t" + << "offset_ACDC0_ns" + << "\t" + << "offset_ACDC1_ns" + << "\t" + << "offset_ACDC0_ps_negative" + << "\t" + << "offset_ACDC1_ps_negative" + << "\t" + << "gotOrphanCount_ACDC0" + << "\t" + << "gotOrphanCount_ACDC1" + << "\t" + << "EventNumInThisPartFile" + << "\t" + << "min_mean_dev_noOrphan_ACDC0" + << "\t" + << "min_mean_dev_noOrphan_ACDC1" + << "\t" + << "increament_times_ACDC0" + << "\t" + << "increament_times_ACDC1" + << "\t" + << "min_mean_dev_ACDC0" + << "\t" + << "min_mean_dev_ACDC1" + << std::endl; + + // 2. load TimeStamp tree, get run number and part file number + TTree *TSTree = (TTree *)f->Get("TimeStamp"); + TTree *CTCTree = (TTree *)f->Get("Trig"); + TTree *GCTCTree = (TTree *)f->Get("GTrig"); + int runNumber; + int subRunNumber; + int partFileNumber; + int LAPPD_ID; + ULong64_t ppsCount0; + TSTree->SetBranchAddress("RunNumber", &runNumber); + TSTree->SetBranchAddress("SubRunNumber", &subRunNumber); + TSTree->SetBranchAddress("PartFileNumber", &partFileNumber); + TSTree->SetBranchAddress("LAPPD_ID", &LAPPD_ID); + TSTree->SetBranchAddress("ppsCount0", &ppsCount0); + std::vector> loopInfo; + int nEntries = TSTree->GetEntries(); + for (int i = 0; i < nEntries; i++) + { + TSTree->GetEntry(i); + // if this entry is PPS event, continue. + if (ppsCount0 != 0) + continue; + vector info = {runNumber, subRunNumber, partFileNumber, LAPPD_ID}; + // find if this info is already in loopInfo + bool found = std::find_if(loopInfo.begin(), loopInfo.end(), [&info](const std::vector &vec) + { + return vec == info; // compare vec and info + }) != loopInfo.end(); // if find_if doesn't return end，found + + if (!found) + loopInfo.push_back(info); + } + + // 3. for each unique run number and part file number + // 4. in this file, for each LAPPD ID + // do fit: + std::map>> ResultMap; + std::vector>::iterator it; + int pfNumber = 0; + for (it = loopInfo.begin(); it != loopInfo.end(); it++) + { + int runNumber = (*it)[0]; + int subRunNumber = (*it)[1]; + int partFileNumber = (*it)[2]; + int LAPPD_ID = (*it)[3]; + if (processPFNumber != 0) + { + if (pfNumber >= processPFNumber) + break; + } + + vector> Result; + if (!triggerGrouped) + { // cout<<"trigger not grouped"<> Result = {FitInfo, TimeStampRaw, BeamGateRaw, TimeStamp_ns, BeamGate_ns, TimeStamp_ps, BeamGate_ps, EventIndex, EventDeviation_ns, CTCTriggerIndex, CTCTriggerTimeStamp_ns,BeamGate_correction_tick,TimeStamp_correction_tick}; + vector FitInfo = {final_offset_ns, final_offset_ps_negative, gotOrphanCount, gotMin_mean_dev_noOrphan, increament_times, min_mean_dev, final_i, final_j, totalEventNumber}; + vector TimeStampRaw; + vector BeamGateRaw; + vector TimeStamp_ns; + vector BeamGate_ns; + vector TimeStamp_ps; + vector BeamGate_ps; + vector EventIndex; + vector EventDeviation_ns; + vector CTCTriggerIndex; + vector CTCTriggerTimeStamp_ns; + vector BeamGate_correction_tick; + vector TimeStamp_correction_tick; + */ + // Loop the ResultMap, save the result to a root tree in a root file + cout << "Start saving the result to root file and txt file..." << endl; + TFile *fOut = new TFile("offsetFitResult.root", "RECREATE"); + TTree *tOut = new TTree("Events", "Events"); + int runNumber_out; + int subRunNumber_out; + int partFileNumber_out; + int LAPPD_ID_out; + ULong64_t EventIndex; + ULong64_t EventNumberInThisPartFile; + ULong64_t final_offset_ns_0; + ULong64_t final_offset_ns_1; + ULong64_t final_offset_ps_negative_0; + ULong64_t final_offset_ps_negative_1; + double drift0; + double drift1; + ULong64_t gotOrphanCount_0; + ULong64_t gotOrphanCount_1; + ULong64_t gotMin_mean_dev_noOrphan_0; + ULong64_t gotMin_mean_dev_noOrphan_1; + ULong64_t increament_times_0; + ULong64_t increament_times_1; + ULong64_t min_mean_dev_0; + ULong64_t min_mean_dev_1; + ULong64_t TimeStampRaw; + ULong64_t BeamGateRaw; + ULong64_t TimeStamp_ns; + ULong64_t BeamGate_ns; + ULong64_t TimeStamp_ps; + ULong64_t BeamGate_ps; + ULong64_t EventDeviation_ns_0; + ULong64_t EventDeviation_ns_1; + ULong64_t CTCTriggerIndex; + ULong64_t CTCTriggerTimeStamp_ns; + long long BGMinusTrigger_ns; + long long BGCorrection_tick; + long long TSCorrection_tick; + ULong64_t LAPPD_PPS_interval_ticks; + ULong64_t BG_PPSBefore_tick; + ULong64_t BG_PPSAfter_tick; + ULong64_t BG_PPSDiff_tick; + ULong64_t BG_PPSMissing_tick; + ULong64_t TS_PPSBefore_tick; + ULong64_t TS_PPSAfter_tick; + ULong64_t TS_PPSDiff_tick; + ULong64_t TS_PPSMissing_tick; + + ULong64_t TS_driftCorrection_ns; + ULong64_t BG_driftCorrection_ns; + + tOut->Branch("runNumber", &runNumber_out, "runNumber/I"); + tOut->Branch("subRunNumber", &subRunNumber_out, "subRunNumber/I"); + tOut->Branch("partFileNumber", &partFileNumber_out, "partFileNumber/I"); + tOut->Branch("LAPPD_ID", &LAPPD_ID_out, "LAPPD_ID/I"); + tOut->Branch("EventIndex", &EventIndex, "EventIndex/l"); + tOut->Branch("EventNumInThisPF", &EventNumberInThisPartFile, "EventNumInThisPF/l"); + tOut->Branch("final_offset_ns_0", &final_offset_ns_0, "final_offset_ns_0/l"); + tOut->Branch("final_offset_ns_1", &final_offset_ns_1, "final_offset_ns_1/l"); + tOut->Branch("final_offset_ps_negative_0", &final_offset_ps_negative_0, "final_offset_ps_negative_0/l"); + tOut->Branch("final_offset_ps_negative_1", &final_offset_ps_negative_1, "final_offset_ps_negative_1/l"); + tOut->Branch("drift0", &drift0, "drift0/D"); + tOut->Branch("drift1", &drift1, "drift1/D"); + tOut->Branch("gotOrphanCount_0", &gotOrphanCount_0, "gotOrphanCount_0/l"); + tOut->Branch("gotOrphanCount_1", &gotOrphanCount_1, "gotOrphanCount_1/l"); + tOut->Branch("gotMin_mean_dev_noOrphan_0", &gotMin_mean_dev_noOrphan_0, "gotMin_mean_dev_noOrphan_0/l"); + tOut->Branch("gotMin_mean_dev_noOrphan_1", &gotMin_mean_dev_noOrphan_1, "gotMin_mean_dev_noOrphan_1/l"); + tOut->Branch("increament_times_0", &increament_times_0, "increament_times_0/l"); + tOut->Branch("increament_times_1", &increament_times_1, "increament_times_1/l"); + tOut->Branch("min_mean_dev_0", &min_mean_dev_0, "min_mean_dev_0/l"); + tOut->Branch("min_mean_dev_1", &min_mean_dev_1, "min_mean_dev_1/l"); + tOut->Branch("TimeStampRaw", &TimeStampRaw, "TimeStampRaw/l"); + tOut->Branch("BeamGateRaw", &BeamGateRaw, "BeamGateRaw/l"); + tOut->Branch("TimeStamp_ns", &TimeStamp_ns, "TimeStamp_ns/l"); + tOut->Branch("BeamGate_ns", &BeamGate_ns, "BeamGate_ns/l"); + tOut->Branch("TimeStamp_ps", &TimeStamp_ps, "TimeStamp_ps/l"); + tOut->Branch("BeamGate_ps", &BeamGate_ps, "BeamGate_ps/l"); + tOut->Branch("EventDeviation_ns_0", &EventDeviation_ns_0, "EventDeviation_ns_0/l"); + tOut->Branch("EventDeviation_ns_1", &EventDeviation_ns_1, "EventDeviation_ns_1/l"); + tOut->Branch("CTCTriggerIndex", &CTCTriggerIndex, "CTCTriggerIndex/l"); + tOut->Branch("CTCTriggerTimeStamp_ns", &CTCTriggerTimeStamp_ns, "CTCTriggerTimeStamp_ns/l"); + tOut->Branch("BGMinusTrigger_ns", &BGMinusTrigger_ns, "BGMinusTrigger_ns/L"); + tOut->Branch("BGCorrection_tick", &BGCorrection_tick, "BGCorrection_tick/l"); + tOut->Branch("TSCorrection_tick", &TSCorrection_tick, "TSCorrection_tick/l"); + tOut->Branch("LAPPD_PPS_interval_ticks", &LAPPD_PPS_interval_ticks, "LAPPD_PPS_interval_ticks/l"); + tOut->Branch("BG_PPSBefore_tick", &BG_PPSBefore_tick, "BG_PPSBefore_tick/l"); + tOut->Branch("BG_PPSAfter_tick", &BG_PPSAfter_tick, "BG_PPSAfter_tick/l"); + tOut->Branch("BG_PPSDiff_tick", &BG_PPSDiff_tick, "BG_PPSDiff_tick/l"); + tOut->Branch("BG_PPSMissing_tick", &BG_PPSMissing_tick, "BG_PPSMissing_tick/l"); + tOut->Branch("TS_PPSBefore_tick", &TS_PPSBefore_tick, "TS_PPSBefore_tick/l"); + tOut->Branch("TS_PPSAfter_tick", &TS_PPSAfter_tick, "TS_PPSAfter_tick/l"); + tOut->Branch("TS_PPSDiff_tick", &TS_PPSDiff_tick, "TS_PPSDiff_tick/l"); + tOut->Branch("TS_PPSMissing_tick", &TS_PPSMissing_tick, "TS_PPSMissing_tick/l"); + tOut->Branch("TS_driftCorrection_ns", &TS_driftCorrection_ns, "TS_driftCorrection_ns/l"); + tOut->Branch("BG_driftCorrection_ns", &BG_driftCorrection_ns, "BG_driftCorrection_ns/l"); + + std::ofstream outputEvents("outputEvents.txt"); + for (auto it = ResultMap.begin(); it != ResultMap.end(); it++) + { + string key = it->first; + vector> Result = it->second; + + if (Result.size() == 0) + continue; + + runNumber_out = std::stoi(key.substr(0, key.find("_"))); + subRunNumber_out = std::stoi(key.substr(key.find("_") + 1, key.find("_", key.find("_") + 1) - key.find("_") - 1)); + partFileNumber_out = std::stoi(key.substr(key.find("_", key.find("_") + 1) + 1, key.find("_", key.find("_", key.find("_") + 1) + 1) - key.find("_", key.find("_") + 1) - 1)); + LAPPD_ID_out = std::stoi(key.substr(key.find("_", key.find("_", key.find("_") + 1) + 1) + 1, key.size() - key.find("_", key.find("_", key.find("_") + 1) + 1) - 1)); + final_offset_ns_0 = Result[0][0]; + final_offset_ns_1 = Result[24][0]; + final_offset_ps_negative_0 = Result[0][1]; + final_offset_ps_negative_1 = Result[24][1]; + gotOrphanCount_0 = Result[0][2]; + gotOrphanCount_1 = Result[24][2]; + gotMin_mean_dev_noOrphan_0 = Result[0][3]; + gotMin_mean_dev_noOrphan_1 = Result[24][3]; + increament_times_0 = Result[0][4]; + increament_times_1 = Result[24][4]; + min_mean_dev_0 = Result[0][5]; + min_mean_dev_1 = Result[24][5]; + EventNumberInThisPartFile = Result[0][8]; + drift0 = Result[0][9]; + drift1 = Result[24][9]; + + // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 + // vector> Result = {FitInfo, TimeStampRaw, BeamGateRaw, TimeStamp_ns, BeamGate_ns, TimeStamp_ps, BeamGate_ps, EventIndex, EventDeviation_ns, CTCTriggerIndex, CTCTriggerTimeStamp_ns, BeamGate_correction_tick, TimeStamp_correction_tick, LAPPD_PPS_interval_ticks, BG_PPSBefore, BG_PPSAfter, BG_PPSDiff, BG_PPSMiss, TS_PPSBefore, TS_PPSAfter, TS_PPSDiff, TS_PPSMiss}; + // any Result[x] , if x>13, x = x + 8 + for (int j = 0; j < Result[1].size(); j++) + { + long long BGTdiff = Result[4][j] - Result[10][j] - 325250; + // cout<<"BGTDiff: "<Fill(); + } + outputOffset << runNumber_out << "\t" << subRunNumber_out << "\t" << partFileNumber_out << "\t" << 0 << "\t" << LAPPD_ID_out << "\t" << final_offset_ns_0 << "\t" << final_offset_ns_1 << "\t" << final_offset_ps_negative_0 << "\t" << final_offset_ps_negative_1 << "\t" << gotOrphanCount_0 << "\t" << gotOrphanCount_1 << "\t" << EventNumberInThisPartFile << "\t" << gotMin_mean_dev_noOrphan_0 << "\t" << gotMin_mean_dev_noOrphan_1 << "\t" << increament_times_0 << "\t" << increament_times_1 << "\t" << min_mean_dev_0 << "\t" << min_mean_dev_1 << std::endl; + } + outputOffset.close(); + fOut->cd(); + tOut->Write(); + fOut->Close(); + cout << "Result saved." << endl; +}