diff --git a/.vscode/settings.json b/.vscode/settings.json index 8dd812a..4f0ff34 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -118,7 +118,10 @@ "GainMatchQQQ.C": "cpp", "UTF-8gainmatch.C": "cpp", "MakePlotsQQQ.C": "cpp", - "MakePlotsSX3.C": "cpp" + "MakePlotsSX3.C": "cpp", + "QQQ_Calibcheck.C": "cpp", + "QQQ_Calcheck.C": "cpp", + "makeplots.C": "cpp" }, "github-enterprise.uri": "https://fsunuc.physics.fsu.edu" } \ No newline at end of file diff --git a/Calibration.C b/Calibration.C index ff659fe..a5ec1f7 100644 --- a/Calibration.C +++ b/Calibration.C @@ -27,8 +27,10 @@ bool qqqEcut = false; const int MAX_QQQ = 4; const int MAX_RING = 16; const int MAX_WEDGE = 16; -double qqqGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}}; -bool qqqGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}}; +double qqqwGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}}; +double qqqrGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}}; +bool qqqwGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}}; +bool qqqrGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}}; void Calibration::Begin(TTree * /*tree*/) { @@ -44,11 +46,13 @@ void Calibration::Begin(TTree * /*tree*/) else { int det, ring, wedge; - double gain; - while (infile >> det >> ring >> wedge >> gain) + double gainw,gainr; + while (infile >> det >> ring >> wedge >> gainw>>gainr) { - qqqGain[det][ring][wedge] = gain; - qqqGainValid[det][ring][wedge] = (gain > 0); + qqqwGain[det][ring][wedge] = gainw; + qqqrGain[det][ring][wedge] = gainr; + qqqwGainValid[det][ring][wedge] = (gainw > 0); + qqqrGainValid[det][ring][wedge] = (gainr > 0); } infile.close(); std::cout << "Loaded QQQ gains from " << filename << std::endl; @@ -92,32 +96,32 @@ Bool_t Calibration::Process(Long64_t entry) float eWedge = 0.0; float eRingRaw = 0.0; float eRing = 0.0; - if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]) + if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqrGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16] && qqqwGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]) { chWedge = qqq.ch[i]; eWedgeRaw = qqq.e[i]; - eWedge = qqq.e[i] * qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]; - printf("Wedge E: %.2f Gain: %.4f \n", eWedge, qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]); + eWedge = qqq.e[i] * qqqwGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]; + // printf("Wedge E: %.2f Gain: %.4f \n", eWedge, qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]); chRing = qqq.ch[j] - 16; eRingRaw = qqq.e[j]; - eRing = qqq.e[j]; //*qqqGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i]-16]; + eRing = qqq.e[j] * qqqrGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i]-16]; } - else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16]) + else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqrGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16] && qqqwGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16]) { chWedge = qqq.ch[j]; - eWedge = qqq.e[j] * qqqGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16]; + eWedge = qqq.e[j] * qqqwGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16]; eWedgeRaw = qqq.e[j]; chRing = qqq.ch[i] - 16; - eRing = qqq.e[i];// * qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]; + eRing = qqq.e[i] * qqqrGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]; eRingRaw = qqq.e[i]; } else continue; // hQQQFVB->Fill(eWedge, eRing); - plotter->Fill2D(Form("hRaw_qqq%d_ring%d_wedge%d", qqq.id[i], chRing, chWedge),400,0,16000,400,0,16000, eWedgeRaw, eRingRaw,"ERaw"); - plotter->Fill2D(Form("hGM_qqq%d_ring%d_wedge%d", qqq.id[i], chRing, chWedge),400,0,16000,400,0,16000, eWedge, eRing,"EGM"); + plotter->Fill2D(Form("hRaw_qqq%d_ring%d_wedge%d", qqq.id[i], chRing, chWedge), 400, 0, 16000, 400, 0, 16000, eWedgeRaw, eRingRaw, "ERaw"); + plotter->Fill2D(Form("hGM_qqq%d_ring%d_wedge%d", qqq.id[i], chRing, chWedge), 400, 0, 16000, 400, 0, 16000, eWedge, eRing, "EGM"); plotter->Fill2D("hRawQQQ", 4000, 0, 16000, 4000, 0, 16000, eWedgeRaw, eRingRaw); plotter->Fill2D("hGMQQQ", 4000, 0, 16000, 4000, 0, 16000, eWedge, eRing); @@ -143,89 +147,113 @@ Bool_t Calibration::Process(Long64_t entry) void Calibration::Terminate() { - const double AM241_PEAK = 5485.56; // keV + const double AM241_PEAK = 5485.56; + const double P_PEAK=7000; // keV + double calibArray[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}}; bool calibValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}}; std::ofstream outFile("qqq_Calib.txt"); if (!outFile.is_open()) { - std::cerr << "Error opening output file!" << std::endl; + std::cerr << "Error opening qqq_Calib.txt!" << std::endl; return; } - // ======================= Loop over all channels ======================= -for (auto &kv : dataPoints) { + //---------------------------------------------------------------------- + // 1. Create per–channel 1D spectra in ADC from stored gain-matched data + //---------------------------------------------------------------------- - int det, ring, wedge; - std::tie(det, ring, wedge) = kv.first; + std::map, TH1F *> spectra; - const std::vector> &pts = kv.second; - - if (pts.size() < 5) - continue; - - // Build TGraph from stored (wedgeGM, ringE) - std::vector wedgeGM, ringE; - wedgeGM.reserve(pts.size()); - ringE.reserve(pts.size()); - - for (auto &p : pts) { - wedgeGM.push_back(p.first); // gain-matched wedge energy (ADC) - ringE.push_back(p.second); // ring energy (ADC) - } - - TGraph g(pts.size(), wedgeGM.data(), ringE.data()); - g.SetTitle(Form("QQQ Det%d Ring%d Wedge%d", det, ring, wedge)); - - // Fit a line through origin: E_ring = a * E_wedge - TF1 f("f","[0]*x",0,16000); - g.Fit(&f,"QNR"); // Quiet, No draw, use Range - - double slope_raw = f.GetParameter(0); - - if (slope_raw <= 0) - continue; - - // Convert raw slope into keV calibration: - // Use the Am241 peak expected position: - // E_keV = ADC * slope_keV - double slope_keV = AM241_PEAK / (AM241_PEAK / slope_raw); - // Simplifies to: - // slope_keV = slope_raw; // slope now directly converts ADC → keV - - calibArray[det][ring][wedge] = slope_keV; - calibValid[det][ring][wedge] = true; - - outFile << det << " " << ring << " " << wedge << " " - << slope_keV << "\n"; - - printf("Calib Det=%d Ring=%d Wedge=%d slope=%.5f\n", - det, ring, wedge, slope_keV); -} - - outFile.close(); - std::cout << "Gain matching complete." << std::endl; - - // === Plot all gain-matched QQQ points together with a 2D histogram === - TH2F *hAll = new TH2F("hAll", "All QQQ Gain-Matched;Corrected Wedge E;Ring E", - 800, 0, 16000, 800, 0, 16000); - - // Fill the combined TH2F with corrected data for (auto &kv : dataPoints) { - int id, ring, wedge; - std::tie(id, ring, wedge) = kv.first; - if (!calibValid[id][ring][wedge]) - continue; - auto &pts = kv.second; - for (auto &pr : pts) + int det, ring, wedge; + std::tie(det, ring, wedge) = kv.first; + + TString hname = Form("hSpec_d%d_r%d_w%d", det, ring, wedge); + TH1F *h = new TH1F(hname, hname, 4000, 0, 16000); + + for (auto &p : kv.second) { - double corrWedge = pr.first * calibArray[id][ring][wedge]; - double corrRing = pr.second * calibArray[id][ring][wedge]; - hAll->Fill(corrWedge, corrRing); - plotter->Fill2D("hAll", 4000, 0, 16000, 4000, 0, 16000, corrWedge, corrRing); // Create the histogram in the plotter + double eWedge = p.first; // already gain-matched ADC + double eRing = p.second; + + // Use ring ADC for calibration (cleaner alpha peak) + h->Fill(eRing); + } + + spectra[kv.first] = h; + } + + //---------------------------------------------------------------------- + // 2. Fit Am-241 peak and extract keV/ADC calibration slope + //---------------------------------------------------------------------- + + for (auto &kv : spectra) + { + int det, ring, wedge; + std::tie(det, ring, wedge) = kv.first; + TH1F *h = kv.second; + + if (!h || h->GetEntries() < 50) + continue; + + int binMax = h->GetMaximumBin(); + double adcPeak = h->GetXaxis()->GetBinCenter(binMax); + + if (adcPeak <= 0) + continue; + + // double slope_keV = AM241_PEAK / adcPeak; // keV per ADC + double slope_keV = P_PEAK / adcPeak; // keV per ADC + + calibArray[det][ring][wedge] = slope_keV; + calibValid[det][ring][wedge] = true; + + outFile << det << " " << ring << " " << wedge << " " + << slope_keV << "\n"; + + // printf("QQQ DET=%d R=%d W=%d ADCpeak=%.1f slope_keV=%.6f\n",det, ring, wedge, adcPeak, slope_keV); + } + + outFile.close(); + std::cout << "Wrote QQQ calibration file qqq_Calib.txt\n"; + + //---------------------------------------------------------------------- + // 3. Build fully calibrated 2D combined histogram + //---------------------------------------------------------------------- + + TH2F *hCal = new TH2F("hCal", + "All QQQ Calibrated;Wedge Energy (keV);Ring Energy (keV)", + 800, 0, 7000, + 800, 0, 7000); + + for (auto &kv : dataPoints) + { + int det, ring, wedge; + std::tie(det, ring, wedge) = kv.first; + + if (!calibValid[det][ring][wedge]) + continue; + + double slope = calibArray[det][ring][wedge]; + + for (auto &p : kv.second) + { + double eWGM = p.first; // gain matched ADC + double eRGM = p.second; + + double eWkeV = eWGM * slope / 1000; + double eRkeV = eRGM * slope / 1000; + + hCal->Fill(eWkeV, eRkeV); + plotter->Fill2D("hCalQQQ", 4000, 0, 100, 4000, 0, 100, eWkeV, eRkeV); + plotter->Fill2D(Form("hRCal_qqq%d", det ), 16,0,15, 400, 0, 24, ring, eRkeV, "RingCal"); + plotter->Fill2D(Form("hWCal_qqq%d", det ), 16,0,15, 400, 0, 24, wedge, eWkeV, "WedgeCal"); } } + plotter->FlushToDisk(); + std::cout << "Calibrated 2D QQQ histogram saved.\n"; } diff --git a/GainMatchQQQ.C b/GainMatchQQQ.C index 3a8c73e..e72975e 100644 --- a/GainMatchQQQ.C +++ b/GainMatchQQQ.C @@ -10,6 +10,8 @@ #include #include #include +#include +#include #include "Armory/HistPlotter.h" #include "TVector3.h" @@ -28,7 +30,6 @@ void GainMatchQQQ::Begin(TTree * /*tree*/) hQQQFVB = new TH2F("hQQQFVB", "QQQ Front vs Back; Front E; Back E", 800, 0, 16000, 800, 0, 16000); - // Load the TCutG object TFile *cutFile = TFile::Open("qqqcorr.root"); if (!cutFile || cutFile->IsZombie()) @@ -47,6 +48,11 @@ void GainMatchQQQ::Begin(TTree * /*tree*/) Bool_t GainMatchQQQ::Process(Long64_t entry) { + + int ringMults[16] = {0}; + int wedgeMults[16] = {0}; + std::vector> events; + b_qqqMulti->GetEntry(entry); b_qqqID->GetEntry(entry); b_qqqCh->GetEntry(entry); @@ -81,40 +87,54 @@ Bool_t GainMatchQQQ::Process(Long64_t entry) } else continue; - + ringMults[chRing]++; + wedgeMults[chWedge]++; hQQQFVB->Fill(eWedge, eRing); + events.emplace_back(qqq.id[i], chRing, chWedge, eRing, eWedge); + plotter->Fill2D(Form("hRaw_qqq%d_ring%d_wedge%d", qqq.id[i], chRing, chWedge), 800, 0, 3000, 800, 0, 3000, eWedge, eRing, "hRawQQQ"); + // double ratio = (eWedge > 0.0) ? (eRing / eWedge) : 0.0; + // double maxslope = 1.5; - // TString histName = Form("hQQQFVB_id%d_r%d_w%d", qqq.id[i], chRing, chWedge); - // TH2F *hist2d = (TH2F *)gDirectory->Get(histName); - // if (!hist2d) + // bool validPoint = false; + // if (ratio < maxslope && ratio > 1. / maxslope) // { - // hist2d = new TH2F(histName, Form("QQQ Det%d R%d W%d;Wedge E;Ring E", qqq.id[i], chRing, chWedge), 400, 0, 16000, 400, 0, 16000); + // // Accumulate data for gain matching + // dataPoints[{qqq.id[i], chRing, chWedge}].emplace_back(eWedge, eRing); + // plotter->Fill2D("hAll_in", 4000, 0, 16000, 4000, 0, 16000, eWedge, eRing); + // validPoint = true; // } - // hist2d->Fill(eWedge, eRing); - // if (cut && cut->IsInside(eWedge, eRing)) - double ratio = eRing / eWedge; - double maxslope=1.5; - //gate gets rid of noisy off diagonal events forming a 'V' about the center - //TODO: These are very likely nearest-neighbor charge-sharing events, that will go away if appropriately summed - - // if(ratio < maxslope && ratio > 1./maxslope || eWedge<200 || eRing<200) //method adopted from Sudarshan's approach - bool validPoint = false; - if(ratio < maxslope && ratio > 1./maxslope)// || eWedge<200 || eRing<200) //method adopted from Sudarshan's approach - { - // Accumulate data for gain matching - dataPoints[{qqq.id[i], chRing, chWedge}].emplace_back(eWedge, eRing); - plotter->Fill2D("hAll_in", 4000, 0, 16000, 4000, 0, 16000, eWedge, eRing); - validPoint = true; - } - - if(!validPoint){ - plotter->Fill2D("hAll_out", 4000, 0, 16000, 4000, 0, 16000, eWedge, eRing); - } + // if (!validPoint) + // { + // plotter->Fill2D("hAll_out", 4000, 0, 16000, 4000, 0, 16000, eWedge, eRing); + // } } } } + for (auto tuple : events) + { + auto [id, chr, chw, er, ew] = tuple; + if (ringMults[chr] > 1 || wedgeMults[chw] > 1) + continue; // ignore multiplicity > 1 events + double ratio = (ew > 0.0) ? (er / ew) : 0.0; + double maxslope = 1.5; + + bool validPoint = false; + if (ratio < maxslope && ratio > 1. / maxslope) + { + // Accumulate data for gain matching + dataPoints[{id, chr, chw}].emplace_back(ew, er); + plotter->Fill2D("hAll_in", 4000, 0, 16000, 4000, 0, 16000, ew, er); + validPoint = true; + } + + if (!validPoint) + { + plotter->Fill2D("hAll_out", 4000, 0, 16000, 4000, 0, 16000, ew, er); + } + } + return kTRUE; } @@ -124,7 +144,8 @@ void GainMatchQQQ::Terminate() const int MAX_RING = 16; const int MAX_WEDGE = 16; - double gainArray[MAX_DET][MAX_RING][MAX_WEDGE] = {{{0}}}; + double gainW[MAX_DET][MAX_RING][MAX_WEDGE] = {{{0}}}; + double gainR[MAX_DET][MAX_RING][MAX_WEDGE] = {{{0}}}; bool gainValid[MAX_DET][MAX_RING][MAX_WEDGE] = {{{false}}}; std::ofstream outFile("qqq_GainMatch.txt"); @@ -134,48 +155,186 @@ void GainMatchQQQ::Terminate() return; } + // Parameters for sigma-clipping + const int MIN_POINTS = 5; + const int MAX_ITER = 5; + const double CLIP_SIGMA = 3.0; + for (const auto &kv : dataPoints) { - auto [id, ring, wedge] = kv.first; - const auto &pts = kv.second; - if (pts.size() < 5) + auto key = kv.first; + auto [id, ring, wedge] = key; + const auto &pts_in = kv.second; + if (pts_in.size() < (size_t)MIN_POINTS) continue; - std::vector wE, rE; - for (const auto &pr : pts) - { - wE.push_back(pr.first); - rE.push_back(pr.second); - } + // Make a working copy of the points for clipping + std::vector> pts = pts_in; - TGraph g(wE.size(), wE.data(), rE.data()); - TF1 f("f", "[0]*x", 0, 16000); - g.Fit(&f, "QNR"); - gainArray[id][ring][wedge] = f.GetParameter(0); - gainValid[id][ring][wedge] = true; - } + bool solved = false; + double final_gW = 0.0; + double final_gR = 0.0; - for (int id = 0; id < MAX_DET; ++id) - { - for (int ring = 0; ring < MAX_RING; ++ring) + // Iterative sigma-clipping loop + for (int iter = 0; iter < MAX_ITER; ++iter) { - for (int wedge = 0; wedge < MAX_WEDGE; ++wedge) + // Compute sums + double sum_w2 = 0.0; + double sum_wr = 0.0; + double sum_r2 = 0.0; + for (const auto &p : pts) { - if (gainValid[id][ring][wedge]) - { - outFile << id << " " << wedge << " " << ring << " " << gainArray[id][ring][wedge] << std::endl; - printf("Gain match Det%d Ring%d Wedge%d → %.4f \n", id, ring, wedge, gainArray[id][ring][wedge]); - } + double w = p.first; + double r = p.second; + sum_w2 += w * w; + sum_wr += w * r; + sum_r2 += r * r; } - } + + // Guard against degenerate cases + if (sum_w2 <= 0.0 || sum_wr <= 0.0) + { + // // fallback to single-parameter linear fit (original method) + // // Use ROOT TGraph fitting as fallback + // if (pts.size() >= 2) + // { + // std::vector wE, rE; + // wE.reserve(pts.size()); + // rE.reserve(pts.size()); + // for (const auto &pr : pts) + // { + // wE.push_back(pr.first); + // rE.push_back(pr.second); + // } + // TGraph g(static_cast(wE.size()), wE.data(), rE.data()); + // TF1 f("f_fallback", "[0]*x", 0, 16000); + // g.Fit(&f, "QNR"); + // double slope = f.GetParameter(0); + // if (slope > 0) + // { + // // distribute correction symmetrically: + // double alpha = slope; // r ≈ slope * w => alpha = slope + // double gW_try = std::sqrt(alpha); + // double gR_try = 1.0 / gW_try; + // final_gW = gW_try; + // final_gR = gR_try; + // solved = true; + // } + // } + break; + } + + // alpha = sum(w*r) / sum(w^2) + double alpha = sum_wr / sum_w2; + + if (!(alpha > 0.0) || !std::isfinite(alpha)) + { + // // degenerate; fallback to TF1 fit as above + // if (pts.size() >= 2) + // { + // std::vector wE, rE; + // wE.reserve(pts.size()); + // rE.reserve(pts.size()); + // for (const auto &pr : pts) + // { + // wE.push_back(pr.first); + // rE.push_back(pr.second); + // } + // TGraph g(static_cast(wE.size()), wE.data(), rE.data()); + // TF1 f("f_fallback2", "[0]*x", 0, 16000); + // g.Fit(&f, "QNR"); + // double slope = f.GetParameter(0); + // if (slope > 0) + // { + // double gW_try = std::sqrt(slope); + // double gR_try = 1.0 / gW_try; + // final_gW = gW_try; + // final_gR = gR_try; + // solved = true; + // } + // } + break; + } + + // distribute correction between W and R + double gW = std::sqrt(alpha); + double gR = 1.0 / gW; + + // compute residuals and sigma + std::vector residuals; + residuals.reserve(pts.size()); + for (const auto &p : pts) + { + double w = p.first; + double r = p.second; + double res = gW * w - gR * r; + residuals.push_back(res); + } + + // mean and stddev (use population sigma here) + double mean = 0.0; + for (double v : residuals) + mean += v; + mean /= residuals.size(); + + double var = 0.0; + for (double v : residuals) + var += (v - mean) * (v - mean); + var /= residuals.size(); + double sigma = std::sqrt(var); + + // If sigma is NaN or zero, accept and break + if (!std::isfinite(sigma) || sigma == 0.0) + { + final_gW = gW; + final_gR = gR; + solved = true; + break; + } + + // Clip > CLIP_SIGMA and build new pts + size_t before = pts.size(); + std::vector> new_pts; + new_pts.reserve(pts.size()); + for (size_t k = 0; k < pts.size(); ++k) + { + if (std::fabs(residuals[k] - mean) <= CLIP_SIGMA * sigma) + new_pts.push_back(pts[k]); + } + pts.swap(new_pts); + size_t after = pts.size(); + + // If no points removed or too few remain, accept current solution + final_gW = gW; + final_gR = gR; + solved = true; + if (before == after || pts.size() < (size_t)MIN_POINTS) + break; + // otherwise iterate again with clipped pts + } // end iter loop + + if (!solved) + continue; + + // sanity checks: avoid ridiculous gains + if (!(final_gW > 0.0) || !(final_gR > 0.0) || !std::isfinite(final_gW) || !std::isfinite(final_gR)) + continue; + + // store gains + gainW[id][ring][wedge] = final_gW; + gainR[id][ring][wedge] = final_gR; + gainValid[id][ring][wedge] = true; + + // write out both gains: id wedge ring gW gR + outFile << id << " " << wedge << " " << ring << " " << final_gW << " " << final_gR << std::endl; + printf("Gain match Det%d Ring%d Wedge%d → gW=%.6f gR=%.6f\n", id, ring, wedge, final_gW, final_gR); } outFile.close(); + std::cout << "Gain matching complete." << std::endl; // === Plot all gain-matched QQQ points together with a 2D histogram === - TH2F *hAll = new TH2F("hAll", "All QQQ Gain-Matched;Corrected Wedge E;Ring E", - 4000, 0, 16000, 4000, 0, 16000); // Fill the combined TH2F with corrected data for (auto &kv : dataPoints) @@ -187,11 +346,14 @@ void GainMatchQQQ::Terminate() auto &pts = kv.second; for (auto &pr : pts) { - double corrWedge = pr.first * gainArray[id][ring][wedge]; - double ringE = pr.second; - hAll->Fill(corrWedge, ringE); + double corrWedge = pr.first * gainW[id][ring][wedge]; + double corrRing = pr.second * gainR[id][ring][wedge]; + // hAll->Fill(corrWedge, corrRing); + plotter->Fill2D("hAll", 4000, 0, 16000, 4000, 0, 16000, corrWedge, corrRing); + plotter->Fill2D(Form("hGMQQQ_id%d_ring%d_wedge%d", id, ring, wedge), 400, 0, 16000, 400, 0, 16000, corrWedge, corrRing, "GainMatched"); } } - + + // Optionally keep previous global histos too plotter->FlushToDisk(); } diff --git a/QQQ_Calcheck.C b/QQQ_Calcheck.C new file mode 100644 index 0000000..60dd60a --- /dev/null +++ b/QQQ_Calcheck.C @@ -0,0 +1,164 @@ + +#define QQQ_Calcheck_cxx + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "Armory/HistPlotter.h" +#include "TVector3.h" +#include "QQQ_Calcheck.h" + +TH2F *hQQQFVB; +HistPlotter *plotter; +int padID = 0; + +TCutG *cut; +std::map, std::vector>> dataPoints; + +bool qqqEcut = false; + +// Gain Arrays +const int MAX_QQQ = 4; +const int MAX_RING = 16; +const int MAX_WEDGE = 16; +double qqqwGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}}; +double qqqrGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}}; +bool qqqwGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}}; +bool qqqrGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}}; +double qqqCalib[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}}; +bool qqqCalibValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}}; + +void QQQ_Calcheck::Begin(TTree * /*tree*/) +{ + plotter = new HistPlotter("Cal_checkQQQ.root", "TFILE"); + // ----------------------- Load QQQ Gains + { + std::string filename = "qqq_GainMatch.txt"; + std::ifstream infile(filename); + if (!infile.is_open()) + { + std::cerr << "Error opening " << filename << "!" << std::endl; + } + else + { + int det, ring, wedge; + double gainw,gainr; + while (infile >> det >> ring >> wedge >> gainw>> gainr) + { + qqqwGain[det][ring][wedge] = gainw; + qqqrGain[det][ring][wedge] = gainr; + qqqwGainValid[det][ring][wedge] = (gainw > 0); + qqqrGainValid[det][ring][wedge] = (gainr > 0); + } + infile.close(); + std::cout << "Loaded QQQ gains from " << filename << std::endl; + } + } + // ----------------------- Load QQQ Calibrations + { + std::string filename = "qqq_Calib.txt"; + std::ifstream infile(filename); + if (!infile.is_open()) + { + std::cerr << "Error opening " << filename << "!" << std::endl; + } + else + { + int det, ring, wedge; + double slope; + while (infile >> det >> ring >> wedge >> slope) + { + qqqCalib[det][ring][wedge] = slope; + qqqCalibValid[det][ring][wedge] = (slope > 0); + } + infile.close(); + std::cout << "Loaded QQQ calibrations from " << filename << std::endl; + } + } +} + +Bool_t QQQ_Calcheck::Process(Long64_t entry) +{ + b_qqqMulti->GetEntry(entry); + b_qqqID->GetEntry(entry); + b_qqqCh->GetEntry(entry); + b_qqqE->GetEntry(entry); + b_qqqT->GetEntry(entry); + + qqq.CalIndex(); + + for (int i = 0; i < qqq.multi; i++) + { + for (int j = i + 1; j < qqq.multi; j++) + { + if (qqq.e[i] > 100) + qqqEcut = true; + if (qqq.id[i] == qqq.id[j]) + { + int chWedge = -1; + int chRing = -1; + float eWedgeRaw = 0.0; + float eWedge = 0.0; + float eWedgeMeV = 0.0; + float eRingRaw = 0.0; + float eRing = 0.0; + float eRingMeV = 0.0; + // plug in gains + if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqrGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16] && qqqwGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]) + { + chWedge = qqq.ch[i]; + eWedgeRaw = qqq.e[i]; + eWedge = qqq.e[i] * qqqwGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]; + // printf("Wedge E: %.2f Gain: %.4f \n", eWedge, qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]); + chRing = qqq.ch[j] - 16; + eRingRaw = qqq.e[j]; + eRing = qqq.e[j] * qqqrGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i]-16]; + } + else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqrGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16] && qqqwGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16]) + { + chWedge = qqq.ch[j]; + eWedge = qqq.e[j] * qqqwGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16]; + eWedgeRaw = qqq.e[j]; + + chRing = qqq.ch[i] - 16; + eRing = qqq.e[i] * qqqrGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]; + eRingRaw = qqq.e[i]; + } + else + continue; + // plug in calibrations + if (qqqCalibValid[qqq.id[i]][chRing][chWedge]) + { + eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000; + eRingMeV = eRing * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000; + } + else + continue; + + // hQQQFVB->Fill(eWedge, eRing); + plotter->Fill2D(Form("hRaw_qqq%d_ring%d_wedge%d", qqq.id[i], chRing, chWedge), 400, 0, 16000, 400, 0, 16000, eWedgeRaw, eRingRaw, "ERaw"); + plotter->Fill2D(Form("hGM_qqq%d_ring%d_wedge%d", qqq.id[i], chRing, chWedge), 400, 0, 16000, 400, 0, 16000, eWedge, eRing, "EGM"); + plotter->Fill2D(Form("hCal_qqq%d_ring%d_wedge%d", qqq.id[i], chRing, chWedge), 400, 0, 24, 400, 0, 24, eWedgeMeV, eRingMeV, "ECal"); + plotter->Fill2D(Form("hRCal_qqq%d", qqq.id[i]), 16, 0, 15, 1000, 0, 30, chRing, eRingMeV, "RingCal"); + plotter->Fill2D(Form("hWCal_qqq%d", qqq.id[i]), 16, 0, 15, 1000, 0, 30, chWedge, eWedgeMeV, "WedgeCal"); + plotter->Fill2D("hRawQQQ", 4000, 0, 16000, 4000, 0, 16000, eWedgeRaw, eRingRaw); + plotter->Fill2D("hGMQQQ", 4000, 0, 16000, 4000, 0, 16000, eWedge, eRing); + plotter->Fill2D("hCalQQQ", 4000, 0, 24, 4000, 0, 24, eWedgeMeV, eRingMeV); + } + } + } + + return kTRUE; +} + +void QQQ_Calcheck::Terminate() +{ + plotter->FlushToDisk(); + std::cout << "Calibration check file for 2D QQQ histogram saved.\n"; +} diff --git a/QQQ_Calcheck.h b/QQQ_Calcheck.h new file mode 100644 index 0000000..90fdaca --- /dev/null +++ b/QQQ_Calcheck.h @@ -0,0 +1,114 @@ +#ifndef QQQ_Calcheck_h +#define QQQ_Calcheck_h + +#include +#include +#include +#include + +#include "Armory/ClassDet.h" + +class QQQ_Calcheck : public TSelector { +public : + TTree *fChain; //!pointer to the analyzed TTree or TChain + +// Fixed size dimensions of array or collections stored in the TTree if any. + + // Declaration of leaf types + Det sx3; + Det qqq; + Det pc ; + + ULong64_t evID; + UInt_t run; + + // List of branches + TBranch *b_eventID; //! + TBranch *b_run; //! + TBranch *b_sx3Multi; //! + TBranch *b_sx3ID; //! + TBranch *b_sx3Ch; //! + TBranch *b_sx3E; //! + TBranch *b_sx3T; //! + TBranch *b_qqqMulti; //! + TBranch *b_qqqID; //! + TBranch *b_qqqCh; //! + TBranch *b_qqqE; //! + TBranch *b_qqqT; //! + TBranch *b_pcMulti; //! + TBranch *b_pcID; //! + TBranch *b_pcCh; //! + TBranch *b_pcE; //! + TBranch *b_pcT; //! + + QQQ_Calcheck(TTree * /*tree*/ =0) : fChain(0) { } + virtual ~QQQ_Calcheck() { } + virtual Int_t Version() const { return 2; } + virtual void Begin(TTree *tree); + virtual void SlaveBegin(TTree *tree); + virtual void Init(TTree *tree); + virtual Bool_t Notify(); + virtual Bool_t Process(Long64_t entry); + virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; } + virtual void SetOption(const char *option) { fOption = option; } + virtual void SetObject(TObject *obj) { fObject = obj; } + virtual void SetInputList(TList *input) { fInput = input; } + virtual TList *GetOutputList() const { return fOutput; } + virtual void SlaveTerminate(); + virtual void Terminate(); + + ClassDef(QQQ_Calcheck,0); +}; + +#endif + +#ifdef QQQ_Calcheck_cxx +void QQQ_Calcheck::Init(TTree *tree){ + + // Set branch addresses and branch pointers + if (!tree) return; + fChain = tree; + fChain->SetMakeClass(1); + + fChain->SetBranchAddress("evID", &evID, &b_eventID); + fChain->SetBranchAddress("run", &run, &b_run); + + sx3.SetDetDimension(24,12); + qqq.SetDetDimension(4,32); + pc.SetDetDimension(2,24); + + fChain->SetBranchAddress("sx3Multi", &sx3.multi, &b_sx3Multi); + fChain->SetBranchAddress("sx3ID", &sx3.id, &b_sx3ID); + fChain->SetBranchAddress("sx3Ch", &sx3.ch, &b_sx3Ch); + fChain->SetBranchAddress("sx3E", &sx3.e, &b_sx3E); + fChain->SetBranchAddress("sx3T", &sx3.t, &b_sx3T); + fChain->SetBranchAddress("qqqMulti", &qqq.multi, &b_qqqMulti); + fChain->SetBranchAddress("qqqID", &qqq.id, &b_qqqID); + fChain->SetBranchAddress("qqqCh", &qqq.ch, &b_qqqCh); + fChain->SetBranchAddress("qqqE", &qqq.e, &b_qqqE); + fChain->SetBranchAddress("qqqT", &qqq.t, &b_qqqT); + fChain->SetBranchAddress("pcMulti", &pc.multi, &b_pcMulti); + fChain->SetBranchAddress("pcID", &pc.id, &b_pcID); + fChain->SetBranchAddress("pcCh", &pc.ch, &b_pcCh); + fChain->SetBranchAddress("pcE", &pc.e, &b_pcE); + fChain->SetBranchAddress("pcT", &pc.t, &b_pcT); + +} + +Bool_t QQQ_Calcheck::Notify(){ + + return kTRUE; +} + +void QQQ_Calcheck::SlaveBegin(TTree * /*tree*/){ + + TString option = GetOption(); + +} + +void QQQ_Calcheck::SlaveTerminate(){ + +} + + +#endif // #ifdef QQQ_Calcheck_cxx \ No newline at end of file diff --git a/makeplots.C b/makeplots.C new file mode 100644 index 0000000..aebed56 --- /dev/null +++ b/makeplots.C @@ -0,0 +1,25 @@ +#include "GainMatchQQQ.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "Armory/HistPlotter.h" +#include "TVector3.h" + +void make_plots(); + +TH2F *h1, *h2, *h3, *h4 =nullptr; + +int main(){ + TFile* inFile = TFile::Open("Cal_checkQQQ.root"); + + make_plots(); + return 0; +}