From 6dc71040a198d44d4f6863aae23b3f90d722f64f Mon Sep 17 00:00:00 2001 From: Jarred Roberts Date: Mon, 1 Jun 2026 18:15:00 -0700 Subject: [PATCH 1/2] Add StripEnergyThresholdFinder app --- apps/StripEnergyThresholdFinder.cxx | 1892 +++++++++++++++++++ apps/StripEnergyThresholdFinder_config.yaml | 18 + 2 files changed, 1910 insertions(+) create mode 100644 apps/StripEnergyThresholdFinder.cxx create mode 100644 apps/StripEnergyThresholdFinder_config.yaml diff --git a/apps/StripEnergyThresholdFinder.cxx b/apps/StripEnergyThresholdFinder.cxx new file mode 100644 index 00000000..642715ac --- /dev/null +++ b/apps/StripEnergyThresholdFinder.cxx @@ -0,0 +1,1892 @@ +/* + * StripEnergyThresholdFinder.cpp + * + * Author: Jarred Roberts + * Affiliation: UC San Diego, Department of Astronomy & Astrophysics + * + * Description: + * Application for determining per-strip energy thresholds for germanium detectors. + * Includes both slow (energy-based) and fast (timing-based) threshold estimation, + * along with diagnostic outputs and ROOT-based visualization tools. + * + * Notes: + * - Built on MEGAlib and Nuclearizer frameworks + * - Uses ROOT for analysis and visualization + * + * Copyright (C) 2026 Jarred Roberts + * + * This software is provided for research and academic use. + * Redistribution and modification should follow the licensing + * terms of the parent frameworks (MEGAlib/Nuclearizer) where applicable. + */ + +// NOTES: +// Compile Command +/* + +h5c++ -O2 StripEnergyThresholdFinder.cxx -o StripEnergyThresholdFinder \ +$(root-config --cflags --libs) \ +-I$MEGALIB/include \ +-I~/COSItools/nuclearizer/include \ +-L$MEGALIB/lib \ +-lMEGAlib -lNuclearizer -lyaml-cpp + +*/ + +#include +#include +#include +#include +#include +#include +#include + +#include +namespace fs = std::filesystem; + + +using namespace std; + +/* ROOT */ +#include +#include +#include +#include +#include +#include "TLegend.h" +#include +#include +#include "TStyle.h" + +/* YAML */ +#include + +/* MEGAlib */ +#include "MGlobal.h" +#include "MSupervisor.h" +#include "MModuleLoaderMeasurementsHDF.h" +#include "MReadOutAssembly.h" +#include "MStripHit.h" +#include "MString.h" + + +/* ------------------------------------------------------------- */ +/* Strip identifier structure */ +/* ------------------------------------------------------------- */ + +struct StripKey +{ + int det; + char side; + int strip; + + bool operator<(const StripKey& o) const + { + if(det!=o.det) return det>tag; + + if(tag!="CM") continue; + + string unused; + int det=-1; + int strip=-1; + string side; + string order; + + ss>>unused>>det>>strip>>side>>order; + + if(det<0||det>=m_NDet) continue; + if(strip<0||strip>=m_NStrip) continue; + + int sideInt=(side=="l")?0:1; + + Coeff c; + + if(order=="poly1zero") + ss>>c.c1; + else if(order=="poly1") + ss>>c.c0>>c.c1; + else if(order=="poly2") + ss>>c.c0>>c.c1>>c.c2; + else + ss>>c.c0>>c.c1>>c.c2>>c.c3; + + m_Coeffs[det][sideInt][strip]=c; + } + + return true; + } + + double ADCToEnergy(int det,char side,int strip,double adc) const + { + int sideInt=(side=='l')?0:1; + const Coeff& c=m_Coeffs[det][sideInt][strip]; + + return c.c3*pow(adc,3)+c.c2*pow(adc,2)+c.c1*adc+c.c0; + } + +private: + + int m_NDet; + int m_NStrip; + + vector>> m_Coeffs; +}; + +/* ------------------------------------------------------------- */ +/* TAC calibration helper */ +/* ------------------------------------------------------------- */ + +class TACCalHelper +{ +public: + + struct Coeff + { + double slope = 0; + double offset = 0; + }; + + bool Load(const string& fileName) + { + ifstream in(fileName); + if(!in) + { + cout<<"Failed to open TAC calibration file: "<> strip_id >> det >> side >> strip + >> slope >> slope_err >> offset >> offset_err; + + char sideChar = (side==0) ? 'l' : 'h'; + + StripKey key{det, sideChar, strip}; + + m_Coeffs[key] = {slope, offset}; + } + + return true; + } + + double TACToEnergy(int det, char side, int strip, double tac) const + { + StripKey key{det, side, strip}; + + auto it = m_Coeffs.find(key); + if(it == m_Coeffs.end()) return 0; + + return it->second.slope * tac + it->second.offset; + } + +private: + map m_Coeffs; +}; + + +/* ------------------------------------------------------------- */ +/* Hit selection */ +/* ------------------------------------------------------------- */ + +bool PassHitSelection(MStripHit* SH) +{ + if(SH==nullptr) return false; + + if(SH->IsGuardRing()) return true; + if(SH->IsNearestNeighbor()) return false; + + return true; +} + + + + +int main(int argc,char** argv) +{ + + // Force unbuffered stdout so progress bar updates live + setvbuf(stdout, NULL, _IONBF, 0); + + /* ------------------------------------------------------------- */ + /* We Include this help menu so that users can define their */ + /* Variables,, without needing a yaml config file */ + /* ------------------------------------------------------------- */ + + for(int i=1;i(); + + if(config["analysis"]["fallback_threshold_keV"]) + fallbackThreshold=config["analysis"]["fallback_threshold_keV"].as(); + + if(config["analysis"]["histogram_bins"]) + histogramBins=config["analysis"]["histogram_bins"].as(); + + if(config["analysis"]["histogram_max_adc"]) + histogramMaxADC=config["analysis"]["histogram_max_adc"].as(); + + if(config["analysis"]["noise_search_max_adc"]) + noise_search_max_adc=config["analysis"]["noise_search_max_adc"].as(); + } + + + vector inputFiles=config["input"]["data_files"].as>(); + + if (inputFiles.empty()) { + cerr << "Error: No input files provided." << endl; + return 1; + } + + string calibrationFile=config["input"]["calibration_file"].as(); + string stripMapFileStr=config["input"]["strip_map"].as(); + string outfile=config["output"]["prefix"].as(); + + + EnergyCalHelper helperCal; + TACCalHelper helperCal_TAC; + + if(!helperCal.Load(calibrationFile)) + { + cout<<"Failed to load calibration file: "<(); + + + + /* ------------------------------------------------------------- */ + /* Modern command line parser using getopt_long */ + /* ------------------------------------------------------------- */ + + static struct option long_options[] = + { + {"min_entries", required_argument, 0, 'm'}, + {"noise_search_max_adc", required_argument, 0, 'n'}, + {"fallback_threshold_keV", required_argument, 0, 'f'}, + + /* NEW overrides */ + {"data_file", required_argument, 0, 'd'}, + {"calibration_file", required_argument, 0, 'c'}, + {"strip_map", required_argument, 0, 's'}, + {"output_prefix", required_argument, 0, 'o'}, + {"max_events", required_argument, 0, 'e'}, + + {"help", no_argument, 0, 'h'}, + {0,0,0,0} + }; + + optind = 2; // skip program name and YAML file + + int opt; + while((opt = getopt_long(argc, argv, "", long_options, NULL)) != -1) + { + switch(opt) + { + case 'm': + cmd_min_entries = atoi(optarg); + break; + + case 'n': + cmd_noise_search_max_adc = atof(optarg); + break; + + case 'f': + cmd_fallback_threshold_keV = atof(optarg); + break; + + //case 'd': + //cmd_data_file = optarg; + //break; + case 'd': + inputFiles.push_back(optarg); + break; + + case 'e': + max_events = atol(optarg); + break; + + case 'c': + cmd_calibration_file = optarg; + break; + + case 's': + cmd_strip_map = optarg; + break; + + case 'o': + cmd_output_prefix = optarg; + break; + + case 'h': + cout<= 0) + minEntries = cmd_min_entries; + + if(cmd_noise_search_max_adc > 0) + noise_search_max_adc = cmd_noise_search_max_adc; + + if(cmd_fallback_threshold_keV > 0) + fallbackThreshold = cmd_fallback_threshold_keV; + + + MSupervisor* S=MSupervisor::GetSupervisor(); + + map histograms; + map histograms_TAC; + + map>> timingCounts; + + /* ------------------------------------------------------------- */ + /* Save configuration to log file */ + /* ------------------------------------------------------------- */ + + cout<SetFileName(inputFile.c_str()); + Loader->SetFileNameStripMap(stripMapFile); + + S->SetModule(Loader,0); + + if(!Loader->Initialize()) return -1; + + MReadOutAssembly* Event=new MReadOutAssembly(); + + long event_counter = 0; + while(Loader->IsFinished()==false) + { + if(max_events > 0 && event_counter >= max_events) + break; + Event->Clear(); + + if(Loader->IsReady()) + { + Loader->AnalyzeEvent(Event); + event_counter++; + + for(unsigned int sh=0;shGetNStripHits();sh++) + { + MStripHit* SH=Event->GetStripHit(sh); + + if(!PassHitSelection(SH)) continue; + + double adc=SH->GetADCUnits(); + + int det=SH->GetDetectorID(); + int strip=SH->GetStripID(); + char side=SH->IsLowVoltageStrip()?'l':'h'; + + StripKey key{det,side,strip}; + + /* ------------------------------------------------------------- */ + /* Fast threshold data accumulation (dt0 vs dt1) */ + /* ------------------------------------------------------------- */ + + int adc_int = (int)adc; + + // Define timing type + int timing_type = SH->HasFastTiming() ? 1 : 0; + + // Accumulate counts + auto& entry = timingCounts[key][adc_int]; + + if(timing_type == 0) + entry.first++; + else + entry.second++; + + /* ------------------------------------------------------------- */ + /* TAC (fast shaper) histogram filling */ + /* ------------------------------------------------------------- */ + + if(SH->HasFastTiming() && SH->GetTAC() > 0) + { + double tac = SH->GetTAC(); + + //StripKey key{det,side,strip}; + + if(histograms_TAC[key]==nullptr) + { + string name="h_TAC_"+to_string(det)+"_"+side+"_"+to_string(strip); + + histograms_TAC[key]=new TH1D( + name.c_str(), + name.c_str(), + histogramBins,0,histogramMaxADC); + + histograms_TAC[key]->GetXaxis()->SetTitle("TAC ADC"); + histograms_TAC[key]->GetYaxis()->SetTitle("Counts"); + } + + histograms_TAC[key]->Fill(tac); + } + + + + static int debugCounter = 0; + + if(histograms[key]==nullptr) + { + string name="h_"+to_string(det)+"_"+side+"_"+to_string(strip); + + histograms[key]=new TH1D(name.c_str(),name.c_str(), + histogramBins,0,histogramMaxADC); + + histograms[key]->GetXaxis()->SetTitle("ADC"); + histograms[key]->GetYaxis()->SetTitle("Counts"); + } + + histograms[key]->Fill(adc); + } + } + } + + delete Event; + } + + + /* ------------------------------------------------------------- */ + /* Threshold finding algorithm */ + /* ------------------------------------------------------------- */ + + map thresholds; + map thresholdsADC; + + // Progress tracking (slow thresholds) + int totalStrips = histograms.size(); + int processedStrips = 0; + + map thresholdLV; + map thresholdHV; + + map thresholdLV_ADC; + map thresholdHV_ADC; + + /* ------------------------------------------------------------- */ + /* Fast shaper (TAC) histograms */ + /* ------------------------------------------------------------- */ + + //map histograms_TAC; + + /* ------------------------------------------------------------- */ + /* TAC threshold storage */ + /* ------------------------------------------------------------- */ + + map thresholds_TAC; + map thresholds_TAC_ADC; + //map>> timingCounts; + + + /* HV/LV split */ + + map thresholdLV_TAC; + map thresholdHV_TAC; + + map thresholdLV_TAC_ADC; + map thresholdHV_TAC_ADC; + + + + /* ------------------------------------------------------------- */ + /* Diagnostic storage vectors */ + /* These vectors allow us to build ROOT diagnostic plots later */ + /* ------------------------------------------------------------- */ + + vector stripIndex; + vector thresholdValues; + vector noisePeakADC; + + /* ------------------------------------------------------------- */ + /* Separate vectors for HV and LV strip diagnostics */ + /* ------------------------------------------------------------- */ + + vector stripIndex_LV; + vector stripIndex_HV; + + vector thresholdValues_LV; + vector thresholdValues_HV; + + vector noisePeakADC_LV; + vector noisePeakADC_HV; + + /* ------------------------------------------------------------- */ + /* TAC diagnostic vectors */ + /* ------------------------------------------------------------- */ + + vector stripIndex_TAC_LV; + vector stripIndex_TAC_HV; + + vector thresholdValues_TAC_LV; + vector thresholdValues_TAC_HV; + + vector tacPeakADC_LV; + vector tacPeakADC_HV; + + vector slowEnergyVec; + vector tacEnergyVec; + + for(auto& kv:histograms) + { + + StripKey key=kv.first; + TH1D* hist=kv.second; + + if(hist->GetEntries()Smooth(3); + + int maxSearchBin=hist->FindBin(noise_search_max_adc); + + int startBin=-1; + + for(int b=1;b<=maxSearchBin;b++) + if(hist->GetBinContent(b)>5){ startBin=b; break; } + + if(startBin<0) + { + thresholds[key]=fallbackThreshold; + continue; + } + + int peakBin=startBin; + double peakCounts=hist->GetBinContent(startBin); + + for(int b=startBin+1;b<=maxSearchBin;b++) + { + double c=hist->GetBinContent(b); + + if(c>peakCounts){ peakCounts=c; peakBin=b; } + else if(b>peakBin && cGetBinContent(b)<=peakCounts*0.5) + { thresholdBin=b; break; } + } + else + { + for(int b=peakBin+1;b<=maxSearchBin;b++) + { + double c=hist->GetBinContent(b); + + if(cGetBinContent(thresholdBin)) + thresholdBin=b; + + if(b>peakBin && c>peakCounts*0.5) + break; + } + } + + + /* ------------------------------------------------------------- */ + /* Shift threshold slightly to the right of the noise trough */ + /* This prevents thresholds from sitting inside the noise tail */ + /* ------------------------------------------------------------- */ + + int shiftBins = 2; // move threshold up by a couple of bins + thresholdBin = min(thresholdBin + shiftBins, hist->GetNbinsX()); + + double thresholdADC = hist->GetBinCenter(thresholdBin); + + /* Convert ADC → keV using SLOW calibration */ + double thresholdKeV = + helperCal.ADCToEnergy(key.det,key.side,key.strip,thresholdADC); + + /* Store SLOW thresholds */ + thresholds[key] = thresholdKeV; + thresholdsADC[key] = thresholdADC; + + // Progress bar update + processedStrips++; + + int barWidth = 40; + float progress = (float)processedStrips / totalStrips; + + cout << "\r["; + int pos = barWidth * progress; + + for(int i = 0; i < barWidth; ++i) + { + if(i < pos) cout << "="; + else if(i == pos) cout << ">"; + else cout << " "; + } + + cout << "] " << int(progress * 100.0) << "%" << flush; + + /* ------------------------------------------------------------- */ + /* Store values for diagnostic plots */ + /* ------------------------------------------------------------- */ + + stripIndex.push_back(key.strip); + thresholdValues.push_back(thresholdKeV); + noisePeakADC.push_back(hist->GetBinCenter(peakBin)); + + /* -------------------------------------------------------------- */ + /* Store HV and LV diagnostics separately so both appear in plots */ + /* -------------------------------------------------------------- */ + + if(key.side=='l') + { + stripIndex_LV.push_back(key.strip); + thresholdValues_LV.push_back(thresholdKeV); + noisePeakADC_LV.push_back(hist->GetBinCenter(peakBin)); + } + else + { + stripIndex_HV.push_back(key.strip); + thresholdValues_HV.push_back(thresholdKeV); + noisePeakADC_HV.push_back(hist->GetBinCenter(peakBin)); + } + + if(key.side=='l') + { + thresholdLV[key.strip]=thresholdKeV; + thresholdLV_ADC[key.strip]=thresholdADC; + } + else + { + thresholdHV[key.strip]=thresholdKeV; + thresholdHV_ADC[key.strip]=thresholdADC; + } + + } + cout << endl; + + /* ------------------------------------------------------------- */ + /* Fast threshold finder (dt0 vs dt1 crossover) */ + /* ------------------------------------------------------------- */ + + for(auto& kv : timingCounts) + { + StripKey key = kv.first; + auto& adcMap = kv.second; + + // --- FIX: skip guard ring for FAST thresholds --- + if(key.strip == 64) + { + continue; + } + + int totalCounts = 0; + for(auto& a : adcMap) + totalCounts += a.second.first + a.second.second; + + //if(totalCounts < minEntries) + + // Proper minEntries guard + if(totalCounts < minEntries) + { + thresholds_TAC[key] = fallbackThreshold; + continue; + } + + // Find first nonzero ADC + int first_nonzero = -1; + + for(auto& a : adcMap) + { + if(a.second.first + a.second.second > 0) + { + first_nonzero = a.first + 10; + break; + } + } + + if(first_nonzero < 0) + { + thresholds_TAC[key] = fallbackThreshold; + continue; + } + + + //cout << "FAST threshold RMS (keV): " << hFastThreshDist->GetRMS() << endl; + + int bestADC = -1; + int minDiff = 1e9; + + // ------------------------------------------------------------- + // Direct crossover detection (no smoothing bias) + // ------------------------------------------------------------- + int crossoverADC = -1; + + for(auto& a : adcMap) + { + int adc_val = a.first; + int n0 = a.second.first; + int n1 = a.second.second; + + // Require minimum statistics to avoid noise triggers + if(n0 + n1 < 10) continue; + + if(n1 > n0) + { + crossoverADC = adc_val; + break; + } + } + + int nbins = 21; + + // Extend search window for stability + int searchMax = first_nonzero + 800; + + for(int adc = first_nonzero; adc < searchMax; adc++) + { + int n_dt0 = 0; + int n_dt1 = 0; + + for(int a = adc; a < adc + nbins; a++) + { + if(adcMap.find(a) == adcMap.end()) continue; + + n_dt0 += adcMap[a].first; + n_dt1 += adcMap[a].second; + } + + int diff; + + if(n_dt0 == 0 && n_dt1 == 0) + diff = 1000000; + else + diff = abs(n_dt1 - n_dt0); + + if(diff < minDiff) + { + minDiff = diff; + bestADC = adc; + } + } + + if(bestADC < 0) + { + thresholds_TAC[key] = fallbackThreshold; + continue; + } + + + // prefer physical crossover if available + int fast_thresh_adc; + + if(crossoverADC > 0) + { + fast_thresh_adc = crossoverADC; + } + else + { + // fallback to old method + fast_thresh_adc = bestADC + nbins/2; + } + + /*cout << "[FAST DEBUG] strip " << key.strip + << " bestADC=" << bestADC + << " crossoverADC=" << crossoverADC + << endl;*/ + + thresholds_TAC_ADC[key] = fast_thresh_adc; + + double fast_thresh_keV = + helperCal.ADCToEnergy(key.det,key.side,key.strip,fast_thresh_adc); + + thresholds_TAC[key] = fast_thresh_keV; + + /* Store diagnostics */ + + if(key.side=='l') + { + stripIndex_TAC_LV.push_back(key.strip); + thresholdValues_TAC_LV.push_back(fast_thresh_keV); + } + else + { + stripIndex_TAC_HV.push_back(key.strip); + thresholdValues_TAC_HV.push_back(fast_thresh_keV); + } + + /* Debug print */ + + /*cout << "[FAST] DET " << key.det + << " " << key.side + << " strip " << key.strip + << " ADC: " << fast_thresh_adc + << " keV: " << fast_thresh_keV + << endl;*/ + + } + + + + /* ------------------------------------------------------------- */ + /* Write CSV threshold tables (HV and LV) */ + /* ------------------------------------------------------------- */ + + ofstream csv_HV(outfile + "_Slow_HV_thresholds.csv"); + ofstream csv_LV(outfile + "_Slow_LV_thresholds.csv"); + + /* CSV headers */ + + csv_HV << "detector_side,strip,threshold_adc,threshold_keV\n"; + csv_LV << "detector_side,strip,threshold_adc,threshold_keV\n"; + + /* Write rows */ + + for(const auto& kv : thresholds) + { + char side = kv.first.side; + int strip = kv.first.strip; + + double thr_keV = kv.second; + double thr_adc = thresholdsADC[kv.first]; + + + if(side == 'h') + { + csv_HV << "h," + << strip << "," + << thr_adc << "," + << thr_keV << "\n"; + } + else if(side == 'l') + { + csv_LV << "l," + << strip << "," + << thr_adc << "," + << thr_keV << "\n"; + } + } + + + /* ------------------------------------------------------------- */ + /* Pixel threshold maps */ + /* ------------------------------------------------------------- */ + + TH2D pixelThresholdMap( + "PixelThresholdMap", + "Pixel Threshold Map (keV);LV Strip;HV Strip", + 64,0,64, + 64,64,0); + + TH2D pixelThresholdADCMap( + "PixelThresholdADCMap", + "Pixel Threshold Map (ADC);LV Strip;HV Strip", + 64,0,64, + 64,64,0); + + /* force ROOT to display full strip axes */ + + pixelThresholdMap.GetXaxis()->SetNdivisions(64,false); + pixelThresholdMap.GetYaxis()->SetNdivisions(64,false); + + pixelThresholdADCMap.GetXaxis()->SetNdivisions(64,false); + pixelThresholdADCMap.GetYaxis()->SetNdivisions(64,false); + + + for(int lv=0;lv<64;lv++) + { + if(thresholdLV.find(lv)==thresholdLV.end()) continue; + + for(int hv=0;hv<64;hv++) + { + if(thresholdHV.find(hv)==thresholdHV.end()) continue; + + double pixelThr=max(thresholdLV[lv],thresholdHV[hv]); + pixelThresholdMap.Fill(lv,hv,pixelThr); + + double pixelADC=max(thresholdLV_ADC[lv],thresholdHV_ADC[hv]); + pixelThresholdADCMap.Fill(lv,hv,pixelADC); + } + } + + + /* ------------------------------------------------------------- */ + /* Write CSV threshold tables (SLOW + FAST) */ + /* ------------------------------------------------------------- */ + + // --- FAST CSV --- + ofstream csv_TAC_HV(outfile + "_Fast_HV_thresholds.csv"); + ofstream csv_TAC_LV(outfile + "_Fast_LV_thresholds.csv"); + + csv_TAC_HV << "detector_side,strip,threshold_adc,threshold_keV\n"; + csv_TAC_LV << "detector_side,strip,threshold_adc,threshold_keV\n"; + + cout << "Writing FAST CSV entries: " << thresholds_TAC.size() << endl; + + for(const auto& kv : thresholds_TAC) + { + char side = kv.first.side; + int strip = kv.first.strip; + + double thr_adc = thresholds_TAC_ADC[kv.first]; + double thr_keV = kv.second; + + if(side == 'h') + { + csv_TAC_HV << "h," << strip << "," << thr_adc << "," << thr_keV << "\n"; + } + else if(side == 'l') + { + csv_TAC_LV << "l," << strip << "," << thr_adc << "," << thr_keV << "\n"; + } + } + + csv_TAC_HV.close(); + csv_TAC_LV.close(); + + + + /* ------------------------------------------------------------- */ + /* ROOT output */ + /* ------------------------------------------------------------- */ + //cout << "DEBUG: outfile being used = " << outfile << endl; + TFile f((outfile+"_diagnostics.root").c_str(),"RECREATE"); + + + + /* ------------------------------------------------------------- */ + /* Slow and Fast Threshold diagnostic plots */ + /* ------------------------------------------------------------- */ + + + // ------------------------------------------------------------- + // SLOW threshold per strip (HV vs LV scatter) + // ------------------------------------------------------------- + TGraph* gSlowThresh_LV = new TGraph(); + TGraph* gSlowThresh_HV = new TGraph(); + + gSlowThresh_LV->SetName("SlowThresh_LV"); + gSlowThresh_HV->SetName("SlowThresh_HV"); + + // Axis titles + gSlowThresh_LV->SetTitle("Slow Threshold per Strip;Strip;Threshold (keV)"); + + // Styling + gSlowThresh_LV->SetMarkerStyle(20); // circle + gSlowThresh_LV->SetMarkerSize(1.0); + gSlowThresh_LV->SetMarkerColor(kBlue); + + gSlowThresh_HV->SetMarkerStyle(20); + gSlowThresh_HV->SetMarkerSize(1.0); + gSlowThresh_HV->SetMarkerColor(kRed); + + // Fill graphs + for(const auto& kv : thresholds) + { + int strip = kv.first.strip; + char side = kv.first.side; + double thr = kv.second; + + + if(side == 'l') + gSlowThresh_LV->SetPoint(gSlowThresh_LV->GetN(), strip, thr); + + if(side == 'h') + gSlowThresh_HV->SetPoint(gSlowThresh_HV->GetN(), strip, thr); + } + + // ------------------------------------------------------------- + // Canvas with legend + // ------------------------------------------------------------- + TCanvas* cSlowThresh = new TCanvas("cSlowThresh","Slow Threshold per Strip",800,600); + cSlowThresh->cd(); + + // Draw LV first (sets axes) + gSlowThresh_LV->Draw("AP"); + + // Optional: auto-scale Y axis + double ymin = 1e9, ymax = -1e9; + for(const auto& kv : thresholds) + { + if(kv.first.strip == 64) continue; + double v = kv.second; + if(v < ymin) ymin = v; + if(v > ymax) ymax = v; + } + double pad = 0.1 * (ymax - ymin); + gSlowThresh_LV->GetYaxis()->SetRangeUser(ymin - pad, ymax + pad); + + // Overlay HV + gSlowThresh_HV->Draw("P SAME"); + + // Legend + TLegend* legSlowScatter = new TLegend(0.7,0.72,0.8,0.8); + legSlowScatter->AddEntry(gSlowThresh_LV, "LV", "p"); + legSlowScatter->AddEntry(gSlowThresh_HV, "HV", "p"); + + legSlowScatter->SetTextSize(0.02); + legSlowScatter->SetBorderSize(1); + legSlowScatter->SetFillStyle(0); + + legSlowScatter->Draw(); + + // Finalize + cSlowThresh->Modified(); + cSlowThresh->Update(); + + // Save + cSlowThresh->Write(); + gSlowThresh_LV->Write(); + gSlowThresh_HV->Write(); + + // ------------------------------------------------------------- + // SLOW threshold distribution (HV vs LV separated) + // ------------------------------------------------------------- + TH1D* hSlowDist_LV = new TH1D( + "SlowThresholdDistribution_LV", + "Slow Threshold Distribution;Threshold (keV);Counts", + 100, 0, 50 + ); + + TH1D* hSlowDist_HV = new TH1D( + "SlowThresholdDistribution_HV", + "Slow Threshold Distribution;Threshold (keV);Counts", + 100, 0, 50 + ); + + // Styling + hSlowDist_LV->SetLineColor(kBlue); + hSlowDist_LV->SetLineWidth(2); + + hSlowDist_HV->SetLineColor(kRed); + hSlowDist_HV->SetLineWidth(2); + + // Fill histograms + for(const auto& kv : thresholds) + { + int strip = kv.first.strip; + char side = kv.first.side; + double thr = kv.second; + + if(strip == 64) continue; + + if(side == 'l') hSlowDist_LV->Fill(thr); + if(side == 'h') hSlowDist_HV->Fill(thr); + } + + // ------------------------------------------------------------- + // Canvas with legend + // ------------------------------------------------------------- + TCanvas* cSlowDist = new TCanvas("cSlowDist","Slow Threshold Distribution",800,600); + cSlowDist->cd(); + + // Draw LV first + hSlowDist_LV->Draw("HIST"); + + // Overlay HV + hSlowDist_HV->Draw("HIST SAME"); + + // Legend + TLegend* legSlow = new TLegend(0.7,0.72,0.8,0.8); + legSlow->AddEntry(hSlowDist_LV, "LV", "l"); + legSlow->AddEntry(hSlowDist_HV, "HV", "l"); + legSlow->Draw(); + + // Finalize + cSlowDist->Modified(); + cSlowDist->Update(); + + // Write everything + cSlowDist->Write(); + hSlowDist_LV->Write(); + hSlowDist_HV->Write(); + + + + + // ------------------------------------------------------------- + // FAST threshold per strip histogram + // ------------------------------------------------------------- + TH1D* hFastThreshPerStrip = new TH1D( + "FastThresholds_per_strip", + "Fast Threshold per Strip;Strip;Threshold (keV)", + 65, 0, 65 + ); + + for(const auto& kv : thresholds_TAC) + { + int strip = kv.first.strip; + double thr_keV = kv.second; + + if(strip == 64) continue; // skip guard ring + + hFastThreshPerStrip->SetBinContent(strip + 1, thr_keV); + } + + hFastThreshPerStrip->Write(); + + + + + /* ------------------------------------------------------------- */ + /* Threshold vs strip for HV and LV sides */ + /* ------------------------------------------------------------- */ + + TGraph Threshold_vs_Strip_LV( + stripIndex_LV.size(), + stripIndex_LV.data(), + thresholdValues_LV.data()); + + Threshold_vs_Strip_LV.SetName("Threshold_vs_Strip_LV"); + Threshold_vs_Strip_LV.SetTitle("Threshold vs Strip;Strip;Threshold (keV)"); + + Threshold_vs_Strip_LV.SetMarkerStyle(20); + Threshold_vs_Strip_LV.SetMarkerColor(kBlue); + Threshold_vs_Strip_LV.SetMarkerSize(1); + + Threshold_vs_Strip_LV.Write(); + + + TGraph Threshold_vs_Strip_HV( + stripIndex_HV.size(), + stripIndex_HV.data(), + thresholdValues_HV.data()); + + Threshold_vs_Strip_HV.SetName("Threshold_vs_Strip_HV"); + + Threshold_vs_Strip_HV.SetMarkerStyle(20); + Threshold_vs_Strip_HV.SetMarkerColor(kRed); + Threshold_vs_Strip_HV.SetMarkerSize(1); + + Threshold_vs_Strip_HV.Write(); + + + + /* ------------------------------------------------------------- */ + /* Noise peak ADC vs strip for HV and LV sides */ + /* ------------------------------------------------------------- */ + + TGraph NoisePeakADC_vs_Strip_LV( + stripIndex_LV.size(), + stripIndex_LV.data(), + noisePeakADC_LV.data()); + + NoisePeakADC_vs_Strip_LV.SetName("NoisePeakADC_vs_Strip_LV"); + NoisePeakADC_vs_Strip_LV.SetTitle("Noise Peak ADC vs Strip;Strip;Noise Peak (ADC)"); + + NoisePeakADC_vs_Strip_LV.SetMarkerStyle(20); + NoisePeakADC_vs_Strip_LV.SetMarkerColor(kBlue); + NoisePeakADC_vs_Strip_LV.SetMarkerSize(1); + + NoisePeakADC_vs_Strip_LV.Write(); + + + TGraph NoisePeakADC_vs_Strip_HV( + stripIndex_HV.size(), + stripIndex_HV.data(), + noisePeakADC_HV.data()); + + NoisePeakADC_vs_Strip_HV.SetName("NoisePeakADC_vs_Strip_HV"); + + NoisePeakADC_vs_Strip_HV.SetMarkerStyle(20); + NoisePeakADC_vs_Strip_HV.SetMarkerColor(kRed); + NoisePeakADC_vs_Strip_HV.SetMarkerSize(1); + + NoisePeakADC_vs_Strip_HV.Write(); + + + + + /* ------------------------------------------------------------- */ + /* TAC Threshold vs Strip */ + /* ------------------------------------------------------------- */ + + TGraph Threshold_TAC_vs_Strip_LV( + stripIndex_TAC_LV.size(), + stripIndex_TAC_LV.data(), + thresholdValues_TAC_LV.data()); + + Threshold_TAC_vs_Strip_LV.SetName("Threshold_TAC_vs_Strip_LV"); + Threshold_TAC_vs_Strip_LV.SetTitle("TAC Threshold vs Strip;Strip;Threshold (keV)"); + Threshold_TAC_vs_Strip_LV.SetMarkerStyle(20); + Threshold_TAC_vs_Strip_LV.SetMarkerColor(kBlue); + Threshold_TAC_vs_Strip_LV.SetMarkerSize(1); + + Threshold_TAC_vs_Strip_LV.Write(); + + + TGraph Threshold_TAC_vs_Strip_HV( + stripIndex_TAC_HV.size(), + stripIndex_TAC_HV.data(), + thresholdValues_TAC_HV.data()); + + Threshold_TAC_vs_Strip_HV.SetName("Threshold_TAC_vs_Strip_HV"); + Threshold_TAC_vs_Strip_HV.SetMarkerStyle(20); + Threshold_TAC_vs_Strip_HV.SetMarkerColor(kRed); + Threshold_TAC_vs_Strip_HV.SetMarkerSize(1); + + Threshold_TAC_vs_Strip_HV.Write(); + + + + + + /* ------------------------------------------------------------- */ + /* Write ADC spectra and create energy spectra with thresholds */ + /* ------------------------------------------------------------- */ + + for(auto& kv:histograms) + { + StripKey key=kv.first; + TH1D* adcHist=kv.second; + + adcHist->Write(); + + string name="Energy_"+to_string(key.det)+"_"+key.side+"_"+to_string(key.strip); + + TH1D* energyHist=new TH1D( + name.c_str(), + name.c_str(), + adcHist->GetNbinsX(), + 0, + helperCal.ADCToEnergy(key.det,key.side,key.strip,histogramMaxADC) + ); + + for(int b=1;b<=adcHist->GetNbinsX();b++) + { + double adc=adcHist->GetBinCenter(b); + double energy=helperCal.ADCToEnergy(key.det,key.side,key.strip,adc); + double counts=adcHist->GetBinContent(b); + + int ebin=energyHist->FindBin(energy); + energyHist->AddBinContent(ebin,counts); + } + + double thr=thresholds[key]; + + TLine* line=new TLine(thr,0,thr,energyHist->GetMaximum()); + line->SetLineColor(kRed); + line->SetLineWidth(2); + + energyHist->GetListOfFunctions()->Add(line); + + energyHist->Write(); + } + + + // ------------------------------------------------------------- + // FAST threshold distribution (HV vs LV separated) + // ------------------------------------------------------------- + TH1D* hFastDist_LV = new TH1D( + "FastThresholdDistribution_LV", + "Fast Threshold Distribution;Threshold (keV);Counts", + 100, 0, 100 + ); + + TH1D* hFastDist_HV = new TH1D( + "FastThresholdDistribution_HV", + "Fast Threshold Distribution;Threshold (keV);Counts", + 100, 0, 100 + ); + + // Styling + hFastDist_LV->SetLineColor(kBlue); + hFastDist_LV->SetLineWidth(2); + + hFastDist_HV->SetLineColor(kRed); + hFastDist_HV->SetLineWidth(2); + + // Fill histograms + for(const auto& kv : thresholds_TAC) + { + int strip = kv.first.strip; + char side = kv.first.side; + double thr = kv.second; + + if(strip == 64) continue; + + if(side == 'l') hFastDist_LV->Fill(thr); + if(side == 'h') hFastDist_HV->Fill(thr); + } + + // ------------------------------------------------------------- + // Canvas with legend + // ------------------------------------------------------------- + TCanvas* cFastDist = new TCanvas("cFastDist","Fast Threshold Distribution",800,600); + cFastDist->cd(); + + // Draw LV first + hFastDist_LV->Draw("HIST"); + + // Overlay HV + hFastDist_HV->Draw("HIST SAME"); + + // Legend + TLegend* leg2 = new TLegend(0.70,0.72,0.8,0.8); + leg2->AddEntry(hFastDist_LV, "LV", "l"); + leg2->AddEntry(hFastDist_HV, "HV", "l"); + leg2->Draw(); + + // Finalize + cFastDist->Modified(); + cFastDist->Update(); + + // Write everything + cFastDist->Write(); + hFastDist_LV->Write(); + hFastDist_HV->Write(); + + + + // Distribution in ADC + TH1D* FastThresholdDistribution_ADC = new TH1D( + "FastThresholdDistribution_ADC", + "Fast Threshold Distribution (ADC);ADC;Counts", + 200, 0, 2000 + ); + + for(const auto& kv : thresholds_TAC_ADC) + { + FastThresholdDistribution_ADC->Fill(kv.second); + } + + FastThresholdDistribution_ADC->Write(); + + + // ------------------------------------------------------------- + // FAST threshold per strip (use TGraph for proper axes) + // ------------------------------------------------------------- + TGraph* gFastThresh_LV = new TGraph(); + TGraph* gFastThresh_HV = new TGraph(); + + gFastThresh_LV->SetName("FastThresh_LV"); + gFastThresh_HV->SetName("FastThresh_HV"); + + // Axis titles + gFastThresh_LV->SetTitle("Fast Threshold per Strip;Strip;Threshold (keV)"); + + // Styling + gFastThresh_LV->SetMarkerStyle(20); // circle + gFastThresh_LV->SetMarkerSize(1.0); + gFastThresh_LV->SetMarkerColor(kBlue); + + gFastThresh_HV->SetMarkerStyle(20); + gFastThresh_HV->SetMarkerSize(1.0); + gFastThresh_HV->SetMarkerColor(kRed); + + // Fill graphs + for(const auto& kv : thresholds_TAC) + { + int strip = kv.first.strip; + char side = kv.first.side; + double thr = kv.second; + + if(strip == 64) continue; + + if(side == 'l') + gFastThresh_LV->SetPoint(gFastThresh_LV->GetN(), strip, thr); + + if(side == 'h') + gFastThresh_HV->SetPoint(gFastThresh_HV->GetN(), strip, thr); + } + + // Make sure ROOT file is the active directory + f.cd(); + + // ------------------------------------------------------------- + // Canvas with legend (FIXED VERSION) + // ------------------------------------------------------------- + TCanvas* cFast = new TCanvas("cFastThresh","Fast Threshold per Strip",800,600); + cFast->cd(); + + // Draw LV first + gFastThresh_LV->Draw("AP"); + + // Force axis AFTER draw + gFastThresh_LV->GetYaxis()->SetRangeUser(0, 60); + + // Draw HV + gFastThresh_HV->Draw("P SAME"); + + // Create legend AFTER graphs are drawn + TLegend* leg = new TLegend(0.7,0.72,0.8,0.8); + leg->AddEntry(gFastThresh_LV, "LV", "p"); + leg->AddEntry(gFastThresh_HV, "HV", "p"); + + // Make legend clearly visible + leg->SetBorderSize(1); + leg->SetFillColor(0); + leg->SetTextSize(0.025); + + leg->Draw(); + cFast->GetListOfPrimitives()->Add(leg); + + // force rendering + cFast->Modified(); + cFast->Update(); + + // Write AFTER everything is finalized + cFast->Write(); + + gFastThresh_LV->Write(); + gFastThresh_HV->Write(); + + + /* ------------------------------------------------------------- */ + /* Build Slow vs TAC comparison */ + /* ------------------------------------------------------------- */ + + for(const auto& kv : thresholds) + { + StripKey key = kv.first; + + if(thresholds_TAC.find(key) == thresholds_TAC.end()) + continue; + + double slowE = kv.second; + double tacE = thresholds_TAC[key]; + + if(slowE > 0 && tacE > 0) + { + slowEnergyVec.push_back(slowE); + tacEnergyVec.push_back(tacE); + } + } + + TGraph Slow_vs_TAC( + slowEnergyVec.size(), + slowEnergyVec.data(), + tacEnergyVec.data()); + + Slow_vs_TAC.SetName("Slow_vs_TAC"); + Slow_vs_TAC.SetTitle("Slow vs TAC Energy;Slow Energy (keV);TAC Energy (keV)"); + Slow_vs_TAC.SetMarkerStyle(20); + Slow_vs_TAC.SetMarkerSize(1); + + Slow_vs_TAC.Write(); + + + /* ------------------------------------------------------------- */ + /* dt0 vs dt1 diagnostic histogram */ + /* ------------------------------------------------------------- */ + + + + for(auto& kv : timingCounts) + { + StripKey key = kv.first; + auto& adcMap = kv.second; + + if(key.strip == 64) + { + continue; + } + + string name0 = "dt0_" + to_string(key.det) + "_" + key.side + "_" + to_string(key.strip); + string name1 = "dt1_" + to_string(key.det) + "_" + key.side + "_" + to_string(key.strip); + + TH1D* h0 = new TH1D(name0.c_str(), name0.c_str(), 500, 0, 100); + TH1D* h1 = new TH1D(name1.c_str(), name1.c_str(), 500, 0, 100); + + for(auto& a : adcMap) + { + int adc = a.first; + int n0 = a.second.first; + int n1 = a.second.second; + + //double energy = calib[key.det][key.side][key.strip].Evaluate(adc); + double energy = helperCal.ADCToEnergy(key.det, key.side, key.strip, adc); + + //h0->Fill(energy, n0); + double e0 = helperCal.ADCToEnergy(key.det, key.side, key.strip, adc); + double e1 = helperCal.ADCToEnergy(key.det, key.side, key.strip, adc + 1); + + // distribute counts across the interval + int nSub = 5; // small subdivision + for(int i = 0; i < nSub; ++i) + { + double e = e0 + (e1 - e0)*(i + 0.5)/nSub; + h0->Fill(e, n0 / (double)nSub); + h1->Fill(e, n1 / (double)nSub); + } + + // h1->Fill(energy, n1); + } + + h0->SetLineColor(kRed); + h1->SetLineColor(kBlue); + + // --- FAST threshold line (energy space) --- + if(thresholds_TAC_ADC.find(key) != thresholds_TAC_ADC.end()) + { + double thr_adc = thresholds_TAC_ADC[key]; + //double thr_keV = calib[key.det][key.side][key.strip].Evaluate(thr_adc); + double thr_keV = helperCal.ADCToEnergy(key.det, key.side, key.strip, thr_adc); + + double maxY = max(h0->GetMaximum(), h1->GetMaximum()); + + TLine* line = new TLine(thr_keV, 0, thr_keV, maxY); + line->SetLineColor(kGreen+2); + line->SetLineWidth(2); + line->SetLineStyle(2); // optional (dashed) + + h0->GetListOfFunctions()->Add(line); + } + + h0->Write(); + h1->Write(); + } + + // ------------------------------------------------------------- + // SLOW threshold pixel map (LV vs HV) + // ------------------------------------------------------------- + TH2D* hSlowPixelMap = new TH2D( + "SlowThresholdPixelMap", + "Slow Threshold Pixel Map;LV Strip;HV Strip;Threshold Value [keV]", + 64, -0.5, 63.5, // LV: 0–63 + 61, -0.5, 63.5 // HV: 0–60 + ); + + // Temporary storage + map slowLV; + map slowHV; + + // Separate thresholds by side + for(const auto& kv : thresholds) + { + int strip = kv.first.strip; + char side = kv.first.side; + double thr = kv.second; + + if(strip == 64) continue; // skip guard ring + + if(side == 'l') slowLV[strip] = thr; + if(side == 'h') slowHV[strip] = thr; + } + + // Fill pixel map (HV vs LV) + for(const auto& lv : slowLV) + { + int lv_strip = lv.first; + double lv_thr = lv.second; + + for(const auto& hv : slowHV) + { + int hv_strip = hv.first; + double hv_thr = hv.second; + + //double value = 0.5 * (lv_thr + hv_thr); // average + double value = std::max(lv_thr, hv_thr); // select the higher threshold value + + //hSlowPixelMap->Fill(lv_strip, hv_strip, value); + int binX = hSlowPixelMap->GetXaxis()->FindBin(lv_strip); + int binY = hSlowPixelMap->GetYaxis()->FindBin(hv_strip); + + hSlowPixelMap->SetBinContent(binX, binY, value); + } + } + + // ------------------------------------------------------------- + // Draw heatmap + // ------------------------------------------------------------- + TCanvas* cSlowPixel = new TCanvas("cSlowPixel","Slow Threshold Pixel Map",800,700); + cSlowPixel->cd(); + + gStyle->SetPalette(112); // safe Viridis + + hSlowPixelMap->SetStats(0); + + // Flip Y-axis so 0 is at top (your requirement) + hSlowPixelMap->GetYaxis()->SetRangeUser(60.5, -0.5); + + hSlowPixelMap->Draw("COLZ"); + + cSlowPixel->Write(); + hSlowPixelMap->Write(); + + + + + + + f.Close(); + /* ------------------------------------------------------------- */ + /* Helpful ROOT instructions */ + /* ------------------------------------------------------------- */ + + cout<Draw()"<Draw()"<Draw()"<GetXaxis()->SetRangeUser(0,30)"<Draw()"<Draw()"<Draw()"<Draw()"<SetLineColor(kBlue)"<Draw(\"SAME\")"< Date: Mon, 22 Jun 2026 22:18:21 -0700 Subject: [PATCH 2/2] Add StripEnergyThresholdFinder for per-strip slow and fast threshold extraction --- apps/StripEnergyThresholdFinder | Bin 0 -> 198944 bytes apps/StripEnergyThresholdFinder.cxx | 1353 +++++++++++++++------------ 2 files changed, 753 insertions(+), 600 deletions(-) create mode 100755 apps/StripEnergyThresholdFinder diff --git a/apps/StripEnergyThresholdFinder b/apps/StripEnergyThresholdFinder new file mode 100755 index 0000000000000000000000000000000000000000..2fce71ea5d302231dd204358144b00de2c6f3c0b GIT binary patch literal 198944 zcmeFa33wF6);`{0GblO(Dn>;exPt~2$OHn3h)zfbdSn8T1dZz?Boh+JW-@`Gg24o3 z?ikHQuNU`-S3j;Bxu}SU7(o{Ih#OuPA}(|r1Vlw#F8tn8RXts4W)RW;_k927v9(iu z-a1uv>eQ*KQ&m09^&a2&ZgFvnadlVDQ;5ock|S;r(V=t~Pm1DJQkA3ecc9`>dI9Z$ z+bv?P;pT2U!U*F9ALT~r5=6GvaDd=5LPJ6dDP629?*<-Ygo>Thl&&W|Hgnn^O#}Ywj!L5pohYrL)iG>8xR^5W@%!eUglH zrwV;q!z9GfWrWsrB_f?Q%;d=wBcw9u%K0K#bU5W@i*yS{3cXmuBGyQXVuVI{7a|?$ z;qM>D_6}5Fk>rzJz9W`=TP35p;mVo4P z;3jb=Oq!~oHf@9nH7eceIc?-eamR^x;=g_L;8%R_ef$2TrI%^S&~3L}u{C}V@sQk< zhIr_r{6v#8dL90Ta651xh5KaZ<=tG}6j#6cc%_HxRu;D_@s93loZZo*TaP~J8S(Ya z%}Nh<3ncA3$O);ymwJ_El-rZ5(eq4PYw-V=WPf9?icK6PQqz?C*ZuVmw z@er!#qV)J4;2rE%s(Lu=eg8TwVe*l7d%>a}iE+*w93$i7dZ#N#ZY(KM9(y_7-c3~o zLMwI$?v{f1{CT}xr=%+hedg!e&+5@5zW#JaoDxP{PbJY&8rQD{d`kVm9zsz^#@k&> z;%{>(Zl%8dZ3O|;0hX3!E$IV#+^F=4Z*=#GQ~GpM63%o0Ctds%w{hKKJczJagnvVL z8}2)B|3hGRBD_mH8^PV;;T{p*hwy&!{D254cu+h)gz#bU{67f)h5K>bPY7%&!dCH2 z!N0}xG6|!h_G#S9CCqqy1`jK6KP!Oe5U#@gyr3y~K|HsKkb)Pb=QVhK3HJuv8wLI{ z!dG#>CTI%Y!2KrfcHD2_eh2saxHsee5AF|fe~kMR+*@(e^(F3YxOLo}xas;D_fFiq zaDRh)H}3Cne~viJj zJEpg;Uw+i*r;PpEm*c1IA3b%*3wyIa&6;(>;OG4xxf5D%epq?-?N_HKzw|@D<@sL? zyne|1MWZiFYdi0!lYdUzv2F0NjQ`SXxzFXv zTKAPRe#+e!zWeH~Vdu&^9`eJ>9C%*WeO?hel2Pb$_ns)a{JO0_zZ`k?G zIWy+9`5#`;E%f;7Cx5l?{sC1h`Yv93-UZLA;n#cDjT?H|obr2*={RB5X@QRy&6s=a zgD;d#yyec!__})vh9AGEPsU$9Khra$?HPyKKl`GwCtk4VtYber*)f08Js%F=^2Q1C zA3S~V#Pv%F;lO8?Z~vpb9}*fZ-4&8Bdb3;KJUW6x%>Z^(=_CltLp#Kvts2* zPkwe^)l)S&v&VFwtzTPR{z{*7{&Ca|-<~ih|Jd`IdpG{$$BH$}kNW3{WB;~n&cdm$ ze~`a<^f3=_A9d{D=Vy;}Jv{B&e=dFdPIv#-)t7X;{_MA>j`O^`{ifE_7C*W8hkoDB z-`2DC@~;N&d+e2c-@LN?s0{;u`FQWfnz6;pyN;VX{f2)` zD%tb*Jr6p|`(*r3^wq!t=bZP0Ybd(Ub^ed;;P!Pn3ANMuE*)~p%>(q}(*}L; zKWFYc>)kys%&fce*xh$+8}{$7i*hr*UhsJGv=6@?aNR3262H6q-D{uf@lE()rStut zihi8jxwCp~_TptV)6aSG$+E&d+im?<&#g&c@zBF5*$wZH*<3T~s&9Pos)n?T z>T|#A*70b!e|_ZM_VxC&HeB`hr-u5=KYcKK!dLr_pL^eXk5Bxp?d;LHz2CZb?RiVo zrLM*2b*}yB{>NUpbH<&gobqY$r`>bZTORla8~& zd)4)KUAOYHh1YGlJ1fI8`L0*C%yfTy`7i1FY(MUKYwJ;G*Ij(Z`%@mOJ6-!?ZP~={ zUV7%=C)VyQTUxDGmw&qF(Mv)Pp8n3G|C%-7`ct=#+xpGcXHG3HP?WcmT-WN=cU8@t z)ckjM|M%mjo#PBF{oB~YxA~-M?J%HM9EmX_3f_E7bo@tKbi4$PRkVDnFzH6alaG&% zKW~qY*Y}E!PdOzz{!*Xl`28{D`Lu6z{*z%JKIpBqD-busXUBchiZ7o%Oy=^vfH zE=K;hogAHi65g_J4XM#D2ATjj!|yjNzu!_ zJcb>fjh^K=pm!y7Q;UG#F$5JkI^s4ZH(3~-Wd8^5yRfjj=_Iy3_UN!dOKQwSs6obTVsq@ zkHi?Cq0XrF`XUDZvKadLM+`d*#OS{xW3*#+4E?N&k^i(9^TOL$>_^kvxEOYIMGX51 z#po9qG1_-T3_s_K81scAhMs?nQLhtYwEH_T>|{#}{m+b%Pg4v%e;h*|FZ2x4J>bfU z(Y`0gn2)ZBQEo?!@isGtAN6gFcE3MHzj!xByZkFgJ02Iq9>&MWXEQVt%?|g*$mh)% z_1zm|KB7&nX!-ZX_!g~QddFy&Nipm&F2+2*I)?pBiD7SPG4g*k#ymJTMn2cXu=5=; z+UwF7^}QQ<9)#V)U9_oy3u*PT^pANxq)OVhcCy6#ra8k!ZCE9?i&#I9?GalBpV0o{aM?Jw-dVTiWr#BwlWd#ZG<_^6Z+=^AUmK747H{a$-1Q z*SRn$8pqn^@p{>LGUZm$kKEUB;8ej+`_4=ccXRw|A^$C+Ud=*pp9y|N^q02DoFBsl zyXG7H?k0{K{ZU6O%U|%<3jRh>-|Xi&zhMs_oy+;_@8kGF!S5F1Rh7_>=%&h2jN359 z1Fn9eee*5+w_)O-adnqvT+M_(K=!jh*pEZx*a>?k`SUk(c?|o!PV^T==qDif!@^Fo zE#s{Yg;PFl`8=N*QGqH^u0zaUbt0dDX!mB(zGzl<^+*%?Igz(-JFhkasx&$o<+tA(93f6U7@>}NPO2uL1#53U~rzZxBcwWH{#7&h58 z_b6@;N;2oaP~hjosEEJ$ah|`??gNA!CRyf}=TTp>D~GTv48QCeBj(e3VgE+C-@{I5 z{4bfu%N6X(FG3H^Vjj#95qAjxq($(fS=m*Bza(dpWgPen<2K>!KnOg8IlEFI^69vN zm&>1yG_h@*?YR`VLAE{@wYAOs{d2^a=WUV$@6!Ja$AM|t^Qk%uur?i zKEHy%l)t+l&qs7Kj0GQEr=p=reyhI}6nd+At*hRK z!p_NlT5$MImzVujlr1U3t^{5r@LNSYZWirm*u!k$KWJh+Nfi9~#<*dzTQ7>Eayy3e zaxuKKD@E8()e#&w>|Ytn>(zcamw&0?cVV4EdbSHa*9-g?*f|?#Ec5sfj5Fk)81DsP zm}S=<=#%*^J2|})?M>Gfh(Y?yx7azoT0r@?2z$%s>`DeB zqkQr#^K!kD^Ev93wqf z(V;f3YGEkv3VGU(=j~p?qm^d(EtG#!0>_Q{rA&+)?hTyLkmn}!N3yGYF&<7~rQi%) zv{y+D&!=ACGl%f}JN!JKw*>wy`YDy$Cj7lDf!~adN$s^r*tyZ}wP=^VC|7at{Ds<- ze+$2*1#f!LWy}jNi+*AC53hp%K>F!;ofr7M=wO4yc%CoTxuP2==U~8R`nPx0&wmB( z_>$MlDEC~U=Z>{q@greB%$~)(?BK;KYrsZ$vt=Cnu0Jo=>JR6OalkF~Y>We=;Q&)U zyDZ~#68fWETK|e-l^-Dh@z-1W#UmJhNIxx)a?UnkaEk`>bJzN6 zVjSrBk#jacKDu5Lc5bzYN0BX!+wGQdFU8K~Tq63rVGoJ02a-Qg$e+NoSAG`jj%JHL zGavEPj;-QF2%q535c8nb{|sW_B7W<7zXbMAxbkk-`hG0%Hc_sz{vC_{&g|_yUSJ;V zo38%CeyslRZlRz0Ll);7`EbHj~Fg%jS>MxxQu?e?ZPIPlZ*G>q zC{P`kURo0jROe=8lvP#)ato%F1r)!3dU<7qzb05v9rXJ-Yjnm-C%KYx(}TgRg6bK8 zYS+k|fOyQPEUT>c)F}Q-C#4d@_|md6S2CjML5j=?x<>hPs{DnOwH3iYk$+l2O~C7$ zlm-zhYW;Cf_42g(9q^_LaY1pT#mD6A|i^30_2 zM(0kRc2%G-SXAc8X24;T;Yu>g;JM}mN6kd0mDS$ToC0sD$KxU~icuSXX+=$-I_Qyw z9_{y6R92T4AQH+gpk&#lIWE7)Q$A&4DWuHx9)?s%H*!3(&7N3zjgg1TpH&+S)E$cA zNI4U8f+>Z;*;N5=p(nMd(m$)TC{W=kq^ME3Ilsj7 zVo|KYIWn`lV3r43kinFTtD$Q`OG&}843lMM7-kL{reJ!&<>D5TURG9F=qYEZUAY&P zRtIYf%CZZJM&PzuK@qH=$m7X$@d&1%+&~>P0rk;CPU$rP4^y6LI+JpOune+lmn#{j zURsDwR$W>#-CIaHEG#Lg_6MsAN`p0M9u!dM@lcSH0>g2QW+qucG76LPUch2JGk?wW zz?L%_d1O{i&P9K|G$$CiI(CuMrc8`2Xs$mu$)8J>oIMfQlCBOkj*Cl+f+bwZUB-tI zk?|p#b_|s!ODgR$ZkTPyZ0@%X93<;+Z-9(sZ_Wr--NGpo@kF_#ldnvHq0sk+I9q7#hFRb82p%9j=uQEeuU%AHVMP*svTArQQzpsuurHzMP8QAZzF zT9HwSF}@<`DP<`{0%XKS=@GeVvWpBvX6(W_8j+I<$|DM#h$OBwls9P_hTuG&jmXtd zt64~C_QWDI6v=9KNMUkO`_`0+=-|%S%JCU(OEfcRApKa z&%7!?%4(+MR8E;F@#J145r=3ThhrK>DVa$^F@|W#G>*u^3hT|yF_LtN7T2We45)h+9a%vWj@>LeVa2%<)DEqXpAwZw<-7TW*r;0Et{KDjcm@rIrRgSCJmN2$zN~#fyU@A@XEgBQHtM z@{KAA1Pe;b(kd#80{)_cU_r_>_`AX`8AmEfNu!(RV3S_#*%<=mdMB~0TqDL66k&*) z$sHIEMO`GKk`a|tJFTV=6$!E!a+)wJo1Sk~5I!Pj7h=0cln2Tyt7rR5;M)eOJqLG- zNj1?cJpz`(TP=FLYh)plI_Sr&RbEg9RY16;8h9h{$VskAuo$Z}REww;&*b(Qsk_UI zp`@&|Xc8JhR!Z`SDS@i$Kn*OTAPA4oQ&$+MVh@p=g(60AJ-|V#4ip3faIKLOn>1on zD22P^R2hF&L20#T3e{jFOUI{hL%OssrAC_a0GSV=xrCI5+Dxu4HevL5!4nyNm) zpi^rCL4QSUS=oX5urydod7yea%$(1GM*X80Z44`dL7=Fb!h*5_=<^^R;z2u$uozjt z##9}nEtAW!w!})c7>f=vlau`ZU`ci5EPq*`VtTN|AE>Ub ztQO!5w-`jIh#jFY~E9IirNf5%K#| zWWRq}ZE>{xxa)5y<}hWis-k|6`m+U1PHAU<4N{%o(!zPO{!No8& z4;q^mC48qut_fV6i0)}Z=&I+Dvq}nrA{8bJDxqdJc9O~sLmydH5G+aYlcQ)%`jR0Z zdM}`5Z_VtALQi&KIW;*mK{GkJTCAFZxr>TQD=>eMBN9V%hRt&GKA6r9zh{b#rK+-m z!holoc5)(h65XouSGpCHI6V-=uVJx;gX%@yG>APMjrLqrb2qSXjd9h$0TL z-MZ@XH|rxOII<8E8n;oSkGT%LNsL7YbL|e=TU{gk-tww4Ou5(t3{)3l8+B4nvOk5P zq*v5=lrhy#D=30CYJy^;H8QSa7zV^t#&)d8Z3`F<9Vtw#l<9lu)!gs zp-OJ&##|Pqx5|El;5PjR?5~aRmyy)+rY+Qv#dmF4ReyB^X4}*?RnURIu)HcL7V1*3 zKE!=UmQQtn4P5?$+BytgRm2rMRrNrX(Wdp?D6F#IZ)H zD0<|>1_lblQ-ia>Mq2E{Ru>coX{UBdPU#e}`B_#&yYH5fGg?~|=K4Jt5xEapG0PaL zf5m;`hW0B=Bgs%bG!vP9I5L-EM~P>0DE1a*c0I&4B}a=dsXS9?Ze`uII>4N0VRtbH zP*-Y6L5;t5 zjs#Gsbx33l2(LdfL3L4jQQ(mLpIB37#Fs-SOC|?1Cs-YI{lI(U@ApAi5iT4Bid9|Z zbPO@SZ&8P#3MOA!fswhkFzBHXjOnbHyliQHnp!cd6l$CjsHrUvD3#cjs3=mf{+$W8 z2T8CT!jWWQ$qaunwwIJZU1?Ck9z}WO%z%QdgO%9erWpyxJ8*Z(0#pO+as`k_uo|mV z#A0WyFnI(W5(F#D02h%xp`e25>5;VT{I$HGl=#X0!4jBllChskdziT}gBse+ z!LlrVYYV_yG23xd7#-Tcms1xOlAlC@v0G*$3qsg!pIWhubXR$Wm@kHAX zWZQi%>-G^XAjKxsL1e%Vg|V{!Lz;qSTB#{YN--lz?SW}wgi#(Jp3L_K;C~qHbx`e@ zhgyfzS;%PUgqSx0YXMhu{TN0qx01QR&TR+7?BN;@W-l^=gXW%e*tyfb;jgr8PH+Ua zjq#R;%kM2MMQ7tZ3r|oHu!*Hu89p5 zK>tX3IZ61_zuNYemwXOrGOlv!v`aKH&n36nntg z3Af=+{#}Cpiscs_VCi`0#cC$>8q6Vh6;j7`jq+5`@h@|- z@KT7m|M15s?tqmUH`vH{S`^jUGiR8nb6inm1H+Z#slch(?2LjkHBg4_q>*%7%pAwx z?vL`<(xQ%3v~fYcXk?_~T+t@rKWK&!zVSf^oIj{28hQ^}4u%}+r2o4!^Y0n3F|`}h zcZ}&Yat8kEBZTFc_E3+H*}Pd==S|1vH+(HR_2Dk|L7dye8MPyQftjRPD)YgtrE-@) zxr@us<%uNZE`ymcH|LKQSQyi<4KTk=?_JiCF=mqAF+d;Gvp=k6V_JwkjFd)J$A7ei zvW;T$dBynwhPIqwTI67c_w+xuari6{z+OtFc#-I!73bUe(CYseQb@am!m#Cc9e&s8 zbk`mjwcE$AOsp;scRj}6RR2S~9Q=@AGB&)5N-OA%eRC2v#|9>*IW`l8CxYS(lJ9Sg-Q+0_X6 zLc5E`*y=Un>YH=8_}PNhr;}WN_L?)V_h|M4fu$PQpE;l!##fF0T1<@H z6A#dg%axv)kz2`^Xr9tS{$qiz10>tUjx=n5-T|iNi1C_bv=tT~X!&bK2Y*5L&?X=9 zcn&%|vRU9ThsW54dN}6UbtsBz|7egy_MZ+-KSq8NV+~ZJ*fB|FX~Fc0%9>zlVU&@I z?e&ieT#esr)rhmvLprFi?u?-amSTovCoTs#2YmB4nVbU*Kgnwu)G{EE7A{A@DsvpdQaq_z70S5VZVv| z6;JcGt=b}8&zx0VNz-WB^Rw3%rCpc9-h4cG-~K%ts zC6(oY;a3$@ z#jzfYm*Th#JP%5B;D6E41NjaWG34sX9}VRtn=yv;T~o%Qcjm%t#zbnP{W#nVV z9ZGqe$4q17VuT0A9Ecda){xt00lvY7FHIi#^|GzHw1lLX(_(=aD?)K zpleSm#iu5^E1!X;D^7`5W{C2Lp1JoQ_;SQC%1*&Qj?sOTeS!{MKLqi;m7~w&@)WT6 zc;#e4U(M)V%9({ucwkG=(m^7z$YzwDB}ek9I_B!e(9!M2%0YA-`x_QgHha1 zalBc4uf@RY#dlE_S#U-0FR|dp_o-SfxbeNL6&5^Se8+5~1vkEP(r&?x@0M(~;4R`i z9vv3k_^zF1!CM9YE(>maM@13e?J)FUe3vrbf*aotv0HHCyC?}(Tzo&xVZqylJc$;( zUf@X<+|chL@f{sQPUHJmODwqYJ-1d1ZhVKQ&4L@>JKAi)jqj*v7Toy0t0KNnW5{EC zpD@9KJH&Uuk}SCK-LO7QC1EK3j_gKZAxlTx}M-pTIjTxLxF<%DnvwL-=OgZrtpW_@4t{@EiSx`0?qg2yTpjR@|^h!skf* zhCNzw!@jJzVJ}wPXn!l7DAh|Y*Bt}jDD08)DZ_1CyJFz+Qa&=jBL-e4;ggKa5w^s@ zH^#ts#lTlw&+P%ufPzm^MXojp&xoM-ci$v@qJ+0gc)NsemhfB&?~w3P3D+e2MhVXs z_DOm!k#OTXXY@TP7e$K8_|6>RWVgm;d|!(2F8}74_?{f$%?2|<};(x-(1flVrLc;q>`550lCH!-V-}p`~;crX0@%=Nx#isyStnpoB z!l^$Rm+_rW!uJ>;LgRaigd5++q-f*2afByJa@vKR6MmG0*NgTd{5IUi)gs~IQ;#gF zO~Th3%m|I|y-_~KzjsH`#`iA??a{@)@NqAq0 zKVQP}RmaFxBH`jwc`U9=O=tmN6e8ZWdmq<9iLL0eSCH#~Kitn3C_^A@!CgH}v5kS!!CEOwLw@dgy3EwQ? zgCx8|!trU=$fZfRGlF7#knq70zE{FemvCi{Y5(}dV&saK@Wcp;u}H#)NO*#T<5Q)P z%OT-sMNrtkgyWNpkt<2UhelA?zl0ByaJPix6Qhw!m2iBzFmh!}I6fsBx$-1@LZ;UyA|PXtGj)o{7MOL zmT

FOu+l32%|`0tsIt;nO6%Rl*A;e1(J;NqC!t2PAx>gcnPAyM#}d@XZo_m4tUl z_zVfxB)nY0cS(4KgzuH`N(onlpHDtcm4wGj_|+0_m+)!{Pmu5$33o_%P{I=>yjH@K zBz&fXr%L!N33p3)orJ3rK3l@GCHxu*&y(yT#9u-&`qusX_ZK%aI*#Z~CQVyp`Wln&PV_R9?m_folkQ3MLX)P0 z75xU2rp*Su+N9%&o@UZ?fTB+^X*xL3$C)%8nCQtSO$R0VAd{v861|^E)4_<|&7|o- zMBlwHQeQd<(YKm39f0VYOqvcp^fe|;2Ojz|lcs|XeX&W?0f)ZOr27+ngGtjthhA;c zbikocGif^5(5IL*9cbv|Oqvcd^kkE!0}OqTNz=iF-p{1zz(Vh4(sWRv@BVM3{tlwI zn)E=TH<|PxqSu)8Ux;32(oUimoAh9!7n(F3IOsQ+^ch4~n{*=4(@dHU7W645O$Q44 zIFqJ>1U=cL=>S0=WYTnSp!YLrIxx_?nKT^~=(~T3)IW*nttL$e2KpwGrh@`~jY-o1 zfxgV7>0m%#Y|?Zfpf5COItb8jFljmf(5p?FHvaW#CQTdu`V^Bshv;!8O&kAuvPsj1 zzdp#MX`^57XVSF6uXi(P+Su23{~W2mo9L}3olf*7lO9L(8k5c-dYMUQ61~`@Jwz`w zY1+8gZ!qZzL|2=%O7t|7rVV<1ib>PPygtsPClZ}((zHRZ4>D=mnAiK6^dzFYnKW(C z>$`u7)Souy^{pmN8}j-llctS$eT_-e2E4w^q;rX0Y|^w5uP-!d+JM(@FlpL&*Q-sM zHr(}TCQTdd`V^C<4R(E;Nz=xI>o=J6bfT+GnkFE9nn{-uJ;kJHLtGzc z(zFq-C!2H`(SuByHpKOQCQTdRdN-4%4RC$;-bnpv<6GZq(zM~NZ!&4x=+@VmG;MI} z%S@U!w)Mp(O&i+!LX)PAZ2bn4rVVVp+N5dYTAyapv|+7JF=^VU*2kGNZBXmU5qkgN z_qN3;#kSMX4{841KlqK$y3*?{^lA&e+(JKLp&z!;_gUyWEc7iF`bG;q*Fw*<&=nTC z*g{`vp)az~lPvUj3w@r29%Z42S?JR(^r;rQzlA==Lie)J`?hv%pFI|OhlT#qLVseR z-?PwfSm^Z@dbNdKZlRyB&<|Va`z-Vw7Wx(oeWQh*YoTXa=n4y6Y@x5T&=*K{7J9XXUT&eEu+R@%==&`69Txf) z3w@)7o@=3JTIdQ3U2LJRw9pq>=t&lOyoEl`LXWc0!z}da7Wz~R-QPkVW1)Ll=zUu( z?H{4FFs5~1Xs5Ooi%9Jgtl8L#KD1U_hG#XjMY|M!oEkpYsXYJ=HJs^8C{oYrF__j+ zweiG|1_tdNWTTzW^4g}oj|HE$P|B-V3GVzs zaQi~rv>hO{G=kre;r9uq^{65H9<046_a_iR>GoiSsl`bQ?}DN8V}pHzV6TMBW?)YW z>@o?vUtkk}HES1AsW$`CZCJpm;p?1nw4l* zJwYY!)HVTDn;M-kJLkIg+9|%aQH&|wN~N)n5T`RG0h%QH z%Klv_VO4rU=XsmP+tsuU!K1xl zJTz6sX&-*Of4{e3y`%FZ>bc0GvrP>@{eyN~oa?RK*R;@OVHCW~<@jzy1$(!;~Fw@EUvq;dw6Kk$sVZUeca@hRl5j zIx=<0dI%D_oz+8Ix{@hr3F_|+^>7|o#K}<7nRhZJJ-LaB*rr{JBD9&1Ub`B9cpcYK z#9}l_R%p9^IvPn-T|Y@Y3uRsp7QGiR^sCRX2}uPmfjnm4U;}s^6>%(bZ&sUDIa|?x zwW)8B5UZTcM5uzOB4UgnvWXZbh$JEg3c^mrae~-|;zq~Edw7D(m5>=h2@eW*a48RRc;MrKhX>~&P+>iGw59gVF`S?QN&B1;FiPwtTdF^s zlzZH6v|l*rtRus`=^G+eIol}Z%9EK84TE@Y`GnmevS#(Zq#ut_D8?rUR6xi$np} z?EG_iJ0eExKn!?pB7i!^fi(OlVWL`(+J*^jdHFdk@Hn<%(x;LwJl zY02pK_LTKJn)KYj6}tUp>d>p4o6&gGUE8S2`LejesQqqM8JJ&)qs&g#}{VbLS_D{sA zq!5&ur+ok^B5miRxIT<=n|3}C`=LgW%W$H1GPUm11`)B`&~p_Dq@<5$P#Jl ziEsm|`Mm$QBj5i0D3o(@VsM*Yc>;)qp zI}b*zKZmYVZ1T{|kbbf2<(25qm~?(fzsR*GeQJ8{RJCd7uJv(BPO)mMcsKMj1s~I& zFX+$L>b&iXP!Bac_S^U3Fhqs9!%sbgz zTiJUAY{bJB+mKvm62(w=W(k3QY8!I7(}9N_wjo7MyKTtXPHh<~f~j`&JU-Q~y$9kn zrKB=)oAwk)KFMk)LR+j^eDF$GW4JXBMo~j~zNUXVTTqym1hCjvdA76IwzDnt363mR z-$d)8N))MXXmn~|qRSWh2_O@w62h+MveInKCUjP##G@D9NFzrhn=c8c$%4X-6aw#{ z5V_KHNEfmpWiTwJcf%QrDntb(L<2#@5C3kZtmoGo}pkKYfDR)TR+K`W$Z8ph}X zOmfg=>&qv93gN;rdLs)Tb8e<^$UmJ5cx`>}aVAnY4(jPy9J#s!aRXj3A$c{3j|IRB<0g8!W9B<7j%aRHG$*f5v=+OXu42k zyLcp#6Qt)Plf^B;s`Q+1vQ$ZMo%Gynvb0F>66v|kWNDY+&C>HOlZ9$P7l|8hJVOI^ zsu{urQA3A_qZUlWGli)lR7IF4!V(eIA=F0W0JO2KHuX7bJ{o)Opw46+du;P3fmxmR z(jo+*&6iz~etG&8=~t%v*HV&iaY)wK7FwhpjR)TIJi#yMvQF){>vH;0cF8JsQPMuMTJ>v(iG5q-3X4sf`afcLPv zY|ti&RO5N{8KSKBMf7kEe|as}TbqCr0J}caZ^LGSw(5(>mXiJ^yla||)zDh~C+Yc9 z>G>Y%d7pT8eeDgsb)_2Gsy2MSD?8WYYID7e>3$#E z4jH^o>kecPzrTZ@`x)t}KJ{aJ(tJ~U)>O4&>`HK=yMGbv1LxsO5ZcS1K@UZcb$!o| z-;7h@u0YyI`IKItg-9XK)X?`i{Cy(7U7i17?U9`xcFqW$g>xKgI3pqLE7dkR!*09( zMF>=jDfyGemusD_wzPLu+oWBH+e~jCsiDnkX!m7Q7OSt?l@OomS(R zo8k@UIuq2i{{#cnWEnK_hSp@JRXOHdqc+qfD7Hpyr=a2f^8!t4otex|yrGZh1ifi+ z-fWH7!b0>GZSUuiq;t<#8?H@In5$%K{3{r|Y3~PnU}Brga`AikjUGR z30J7o1-*D1u)V9)x-&z2Xy%(rOM&@q*etUBy-l6JsdmL`sw%BOF8CBtFeP6|lcrFk zeW+57%(QC9oTpuFD7&_?^G13v4jOaq_l&7^_(ET~ZTDVzQ5uj=z#)(nTG&n4*4qA3ntpSK}VY<>Qnc5b8V^QSF;s6KPmA-lby@20Au|H35S zhDq9`aZa`2C)=!XY8dwNQ(Udfw!GJZBh>H!?AgshD_;@6z6bhB80_Zu z%{a8Wotm7DC-XN2-|sxBSq**Z+N51XHbL8>@M*&2OL*+CTVLo)?KSF9BHzx3wI)fZ zgyRyt*gmJ41^bJ6ZjBmtL%O!w?|4*r9Atv}X_KC*LAT+_7*DuxcXTk}E_2X;)He^p z@utwLP&1zI2DaW{do3j20rHF_9fW^_Jl&Ci%kv`!L+;mX&>urLmu8e^m~O-N#4bH7 z0X;132V3K6YCoE&u19D19wI}vmpkpWC$PcTE_l&41l?u+eleJb-B98LEQ65}D(guL zhj_FvZ|5+2RTEk_RD(`Z;{eWP7vjt&?+GsL(QHd+jpQ$d-crL;py5rm->IRANVrC{ z1L0@~!oyP&$>~QtM9hG+If>eBbZmeDKspUQ5&bMt^fPq>BqKTXE6`fn(;pftTB|Qb zmygWX)c4q6q&94@#~zR{=*-_!+mpnZtA;QPHT-Owm4Ux;wL{Zw%X>H%(1;RsSh`9Y zmeB6=KMcOtIT({a-$TNT#fPXTKSF(?4pBosp_=lcuk$X;pw|ZuqVg!2T`Ixh}tE&&sSH)8i|4i@`iaD+FsMIVe-#?%|38(ga0#TOsjy;$p7t0qw~%`awHFWMpYimJPMmw)c2@AJCo=j&(A}LEDq<1=FG{^fU9;nU6jm z{g@_f9B0zF?8RVahYB!}XF46>j6?$1RF)1P5feT(N8-~%Y^V+%O@7Tk$=DIe=62HS zJt;J+b1F9qJNepEq#+b6Nz?3~BGG+In`Fa8}M>=D^uKgc!aD1T- z?r*O(=AtqTNi@yRpzUsQ_b>*aBWhS|p(9@Ibo2TC!{88UKA>gO32c_w!j?^>Qf=nD zloOUstMfsUB@{8_>)J<|!)K;Hz`rurL(IdX2@j5(FFY+C>Jpub@+_uS_TgVE133Gs-$AzD4PT$2U5&RD zuvwS||A8&OwDz~Dcki$Ll-X5NO{%ahI09}~_3#`w>}*2@(+yk%RJJA?u@S|1wC$u} zD8`&0yZ8rOJ(&=xi+-l&!8d@9k>Q4KkQc*{5AZ#|B3tGMR}ovbHR2r!v>W_SuABLv zLwhge1JqinQ;S}67`3xC7Q*ngBY6#teG#esQHQAY(@#b^gJ14xY;#!5dmB8YVnxq` zxY7JwI@ja=x3D;XjkqJ#_A7xCsb(YL-^H~GYI|0>gQ*Wx|v^nX!Sl8{L~F}YA0hP z+lcAecz#=Y{zQ7#BhPDf1>@`gALeT}EJi!AoixwRw?FkkbiE!t;@iaB2s z6OGcV1yk4gn#bQ)Skl-W%0jw)E4bnN)^lhA|`6j9j|kyuI4n57w(OrJX=$(#5vn6CAdM zZEMZ(&V{{Qc|9G=|2`VaducsI6ELRHUJM06TMC%G`Q*Bk&9U8uQWgA3C#Jlv{J*@B(g0Q zzKD7E2FeqCTx^rRul?mUW(u|Iz2Uiu$cC1M*uZ)WjEo0IXxjc)8IP^;4k|dbqq7-< zK9=G$1oQOF}IH$S)F1G86379$>kD9(>dn`a(}Y0(cd8 z$B@n)NRmMYFC9qMvdljWj-#a-y)o`iY#y%i(MFfDdt7&Ka+hLd(i)7%lLWQ7hq#U+R@>!GGkMx;RSjp zK&d?jPubn@WNUQe1ut*Yy{#bVH4J(w8cq+ovJYnhRa8DRtA(w13hYA$T$7F_xbX?rox$OOr^`@pmF8UVzq(@@u-s zSp`YD{(T|kzhC=?Mx>^c#URPPwu1DA@1@r%vYMW!1Fx*6 zrxj%GOZ%B0MUdaT+h(7bb|I=<; zp6v`*C61|c)E<$Z+Fsi&?R&TF{`K%A65Ne@R9hx?4sDqmY47JPU)nb~wm>en<>}7w zg^6ji9k#~29-!}UYizU4uS66?xG)K2<~r@^^?OTrJf6yIjaTA5II7DXWMQGT-n0$2 z`4_-8>7bh$q8bPDG|sNaF9MnFkaINMIke1|DHUuVXY5O`RqP9`;;+2FN?PODVa0+x z$0JXzmNEpYik97`b*BvIWh!c8w{|;Mc1ku^b}fycO(~~?BzcM%dCKH}$0;jMf@vV? z2xLoD`4Elh4L?oxM&$O5NQjNT+fhTbmK!iyD~Sb9=ObW;Oti@?{MPXN5}-y?u~E|t zErn~)N>{4Zno3Q?8C30;ys2KFyT$NdxB6&Tvkwy6mQQeoFHTIG zg~oXmxr@g62HHg9_|Z7_LTaG_PW&k{nkfy-J>E>)XyJ#->?9qDX3C*vYRaS+oHo!< zZr5BJfY2!aGqE={8gDRDY^rvw9P4YsF9c9Jne|WWMY8U!anXd@HiDXt&wdChZ94h` zG$-x9>hxzfe<|y1Zol84B)l4htx(?dDetnTf67^cer_8gOkfeqZk4kG#^`Gr{ON-* zfFcqHQnV5l`X}qlzNVqiqc3am)VkqQ(!r;OhmNBFoRlTjx|$k7A4V$l_Z}4b(C;Zu zdl5|xVMgL*40L12K$;)&l(P!;(vCA^AcEHqg)|pu7&c56;r;uUs7EtxU5;^;s)UqyXQf9K(-#QsZ^)fq z`swn+6Wjc2NRoLk(V!mXyfQtt4Fl4s<5h}YLrEZq19H5Fyo|I<5h>)jA1siA7qJEkHh97&&$XOK$a*ie zf+EL#60*gd9nbV-n=hJ25+B@M(lvFw;T%q@7d1N+Vq#aI@bI zXoHA%0_*iDUtU0cV(jF{kl&?Zgj+_%w(O-!;zt-v<7)RhEj-K7@J0`YHnM z`CX&w$Y3Bvw<21*8~GSZ+ke4A`SV)+G3i;8o_9r_H|R+)dL!)!>3LM-d96;)gOTSF z={X=hA20D=47}O&zbUEYbI;qp5e){91CAEv*DhWVy{PDFYuFcOyJark+TBQ`Ue@ic-IdSsdKID`~WY0WP62hNh~A zm|{o=Y8caGqINF6P2g_CLSrq(Kt(iFlc#|FTx7UYdlkPOfsb)6lVAraFZ7c({y%0Z z!7I9y@M>>LA1Ogukm)=_dkMtXj={hFMN!WO9@aZA-tJIC1~3|AoeZV_4o1c zA0qj!6iTuU=dwM9f2)e8B-zf{8dJ7e(PTSEWS$&LwwKZJrfeU6%4N#9pNBu=;}2A| zZz+^yyNyycjXk4_Y{R~iWD9iB+qW415766vBJ(>C8hSfVE@fyJ+1}K-OcU;fCxaJP z$xca+7mI8rTV&k%CX;a?{y`Q~e<$shWSrhb#?zw7_!GW{K?>N7(2&tB%NU|kB5nKJ zH(aLfXxrV99z#<8TV(Vv3MEZ?xs>~|WGQ?66jRCq=D1MBnT>M`%oUmYBq?zqVcPj; zZ^%+U^dpz@3TTo`xmW1&D3MWjgq?k{f4{}(HQFwGCxm4QUoOC;-KlNIDjnqI5#hC~ zj*OFf=w6yl8dk9$iViy*9w+TY;Oz{`8DmW|%B9OTq|sTAP&?_kI94S!d^?N$);9NM z_VYG8a26~?8rNa;qGhpz4mXnEM7~Jp21&GgkO-~g7>SA@S2a5eN6=9UqM!WrYsPTs zZP?(@f|w?<()QKvq}>MWF5pN5hq%5WJGB2ihZ>;yc2|?ic4}WCq(T$DX(+|kIF3zS zP&Cvq*+*@-Q+p3Dx`QlcBCip>1d$*r^)&j))48i?U6D)6 z7TvYS8*1}n4dd44x3RMb?e!PAFaaUV6%pntPm{5p@?<495T6jAh7S93jVrBD=IEsEYTjMe)}wxnf6e?$l` z`od^M`$W-gmZDcs(YL<9YyCb%peCpkMNdcQMbU1xDVJ9HFj4$&1S!zW6woC%Ilbgg zb)&FF#4SNc7tOQT+FsRA-3+5G1g#yvtWBr4$G5OS{z6FMX~f19y#sI5{0j1% zpu*4h0c$$Gijl15w>aFNqjK%=q!RF!E|&Qya-0{5j)&Pa?8O@lJ>Y3IXoU@X z?X%wT7G};@d~w=&?4?>G{a7!ubnVfmpgNtcc!@R*^h(T*1kOQobl%O*SD+R(oH78B zr1~CAwdp^n)VR*&@1; z!+#g>@d7UA@WUdreP|{s<01}k5b!nupU+`U!0iG~dm9i6w!By z=-n$>^>PHx^Oyj*EDX!^rwkP~^nYwkzBucL0fr$cLP zR8fQpT$p)ZfP^ijYR^#VcuLY97f*@WgW^fm?i5e)+VFDfxNJNHrsvV(2UeRvFe+(B zsE^-63oT~zioo}l0dK`hNZ=Z1;_o1Q_YV=gUEmuDr_Ue={uKgmCHz#AzlCV}@g?y! zrNrT%f`mSgBz?Y%DpFO4!mRZZAZMii3l~VVOT!xMV6e#O6q2PW#V=AkffSMM@c;WJD>#b3*VvTu zeFh(=`hwpZz9`-s9z|j;L|apMDJ91T6?|!!$!gyu^-SOe<&EVDH#3EA(2_X(j>zvt z!0Ax-X09HsfTEFMZ_c(+u-yeVxU_aOwO~}Kfx*~85;Uc(6pU9#7-C&@-lX9<@(oU0j zCP~l9;+gm-OXX%r&!fe&PVR+KkwFqef9ZJ+1_p!YTIqSW^!%mte6IBThV*=e^n980 zJOveYJiDyLOkf{`QL?HwU6JZd`^eV#AZ9tNLp-jp!xIv_PjIHktw>Ztgw&yvl^ z-*)T2X-;gIlc&_CQ>QLujO$b^Qj=L zB4yzfBX=C{rLA~Rg+$P8%>ZAaP`oWP|FU+hGDQ=e`PKlYv*pSH$!^EgCkA1#f@ z`!T&W2v=r*dei*#5eL$)ymHqVPgCt_P(T9dI}hLh;vrqd!!Js28s7;8cwCn+MICw%OLW z6RndzMs?V3{sK=to6e)fd-(@X?@L>2YkV`3ejT1X-SeGcpVNm$UT3_DjUc>ti}LZa zitsfFX&0b;HS`|2I1p4*(>@9I95cx==M1&UlSHS~NJnSyxpZhxKm4Una1^5rc$VuM zUX07^UXh5`BYl}ToV~d_l#?0y3Zu6@bI8_E`wyFK>Kae)?JM2mY>F*&P1>$GeV`;e zY~whD0}v*pH>{6KPg`%B{{)sm$alc}mu(CD5FIb6rh;Jr_DVMSFnu`42~0$@;{Er@ zYrRb`4U z;9p{U8ID~KH|{kn6Tqw5GV!{j;(hV+p1$r*J5AF?&IEQy8L$2O33ggCz=?7bdAZ17 z<8Ip}K4(0h*5YX`l3aq??Y8B;NU0G6^i@IT^>E@$aqe2J7o=1}3HrB~r1*M;vqY;-be3xqwK+IzH%6`tzOL{p;~ZS(2dAZTwrlR{|F z5pJbkIDgOeLs=iRHP*5Dk;aFnptDGAT0vD(!;3n=rWOAMU(&dY%$(BFR19Ce3r8k# zq=6$m!GUH>)i+RnO)0q}*&^Ygh49{6p)iQ4#XpX^Ou{6qJ@GGdWG#VPWV~OYHvMA( ziswJZ!cQe|%zzF^zm~wyOQdgatN@B(kap74@J*CN9dlDNg5VKqh&`xKBzo%KT3Bu7 zu>kCDbUOBW5s8lEZKMIqSYi5+x;2puAMs-RU7W(&MT&~+lL7x*GUQ`IQHR)KF6 z^kxxu2z)Pvkd2lSOg23470-B#N3;RkrRPNP43<<0c1zEBCQFF~S4qzcOqN9w+#)@% zFj+QAaJ%%ZnJjxHIPwh%Xd#~2P9y(*1hkzf=u{E9DFn{;$-qx%&!oL7fiDpBA`!Mw zh`J{rqzjLfC0zt3NzbauLb)3(`OM$NEsCe1)U9d&DQuHq)ZR(@T6_O{!+0o)-^w5&(hoRyKcFV zBK!`~HtZ#jZNhgX;c+Opr;hy}NQjD{{Ek_@g^Y6^Q^UvNwVCP(q`@0DQ=6o1w9W4a zal;LtLyef;yCHMT@aeWm>qG1JFYG`Zn~B-j%BEgk4ZMd`VD^0TP`;eXZtf>+0`pk^i)o7TVWunA16nu?!ms3FXVtDFkO z14;n|A>bq-AA9GS(*82_59K;;kjpf3&oi2U|i(qn^sc+2hjj11Bk(oxTZtTUp zf=Kdxkjqm{1h^j?A&K??y^8^iWB^9vxYh)0Hx zMz-pU`?_l{EansRbQFOE_lX31A7%-b@B}nIlSfYpZV?GaApy2q-(1atjXZdf2P=8- zBo7|r!9zT_4}lt9#lKH@p&&X)@JvQvPZpLTEdub3DOv<<1Zpk;==bF?>A8J6$2J=s zWyFCU<>=INnZYveZ1~3NcRZml^=v4Yx(Rg=^bFk%n*4Hh8b>19f~?x^`%T~NoyQEnHS(*!mq35`oJms-mT5B*wTB16RIw~c^M*FPfF{R2 zC;j?|{cagYT_5;Do3trpdkvn%IQp?0J_4W`-cl zScIU9lS*!hs|}~YJ#F|d5;+Bz7>9Q|HsDm4Z3pMmJV@P4yOMeVcs3%> z7wF@Vm^ta0E}8-Hru8!lb%_xCONNPEM5JLe7`6uK!y zWS)q`T7rMnjzXHBDRO~eUL?X6fv=zt;*r}x7d>nIk)$SR=pvXTHl7p2GhmX*fRm(W z)np;b43>Q9xz1!E`3#n3>3NCCLZTWhZN@XIxDn4()piOY&}I>*O=uxgpF4K)ENuAV zwQUV;LqlL_ZL!x@@g3T2pKS=fBD(Zm8Ym|?*ejvFa3Vf&N31L}>XmfSdy?e8(z!rm ze>?4%?4&mqsF!Lx?x9YiyEd_-blMg<_W`ni=cp#w)7P5N>97MkekKj&=t&Tk()~a3 z{sg}2;`;x`6HPQC@dhOt_oz`5mo})WiHJ=kf%nB5K?Q}n(V~bI1@%U72a+h)>s8vS zrD`o!TbGYoKdn_ng#-wTD{8eaRZvmp8Wa>|ktM(9Yv%pF_h!+yem>vNzmLa159Pk! zvz$3|=A1KU&YYR|Sjj@7V$+f0uK|_KEVc%!BG>nV4jUh=)yqwb<~0?KEtU)TeRKPr zykO3uJ%<&u+dVt`1Eweg$nD0;D%nS+|_z_sV< zJhD9xQ5zp-4@BfC%6|L0s#j}5zvH0ju%f!A*RXxf2gi0FOvA9(5u$f+#LzbKQ z8cdzG^?9=SQfnF9?dBFzt1YKa)nrH;5_~M+BjZPDFw+dKX{#ulxN)Lw6i*Sx%=nVh z!2fMebH|6~xitXf|Fy^upa_!S^Sn34b%D4_@%MOd%tfVG^z(rvN04*iwc2$->=8>< zPPs!Rx$?Z$k69rcp%uc}S|OZ84-GxD#yUMdp4tDC)ItGGF0iBYf1MTzBZ7qjDI@Eg zT~=-$7$C`G&$yf9Z7^~2YG$(t3MB8(rA~x?1}jV2iv#XnV_umFi8~30lfdgjz)u`} zH8rgu!qqxQUM_Wt4srUz%5VY*JqqInbKtK5-dPyW?be~LGLvclaaouY;^ zq%eoR<796?Jn)Bd!;j>!CJ`;krqM4?L=#;Z_Pw+osI8peuy1YE(rfmYBGbJ)Sfl#{=U*llN8b2~F4=-_2I7E^Pt|zKGQh1#GU1N??x|PCG`-Bwk&x4b~2A=)}DZGwi zl0q?Y6O01?{+n|6AWY7bI4cmxEojQ`ni6{j^8WvFUYtos5Wtzf0P58or2FwR1@Ox| zl)!zv8dWENR@9nX5x{VKo~NG0Zyaem`0Gjw>o)6M)zZ)2=(QGedXZi1iE_HoDNBZM zms`D9x4N6EwJF+S!97*2pT#$-4gG~?zya;u}Z%8PZZ0OQ+1D6w1_GGe_9= ze4j@+cB`N$J9Vde=X*VFvuh&C*=Az4efe($JskVBS~YjFbFMUDY7)HQ7&@KUO{Lia z7ve0?OSuXhp*>q(vh=dQ|FR5lW>@>-J{EEMg)$=0kj$0lN(m^W*vLzUYe92W&g_dW_16l6Pd7)tOoY{#T zP*;Cy^;(zYc3k;mu3<@o%-$?XcK2<1p%OnYoUYQHp-?W9)gY!WrO9 ziAp=?91D6#*Eq}zd&phRMnc0wjr&PU$E#Vbv6DAtazS?7>s`0!XSMOsBso`mP2BoED)o)PNCKCg|fVrr}|VehxdG(M~dwP>^k(_edR zm`#NEh_kp$o-GCQ=V~OG3`$QVlPUwJwj-VtYO@4;p?1P zOywgvr9b^$o<&PLh%&;b%l;y_wEYLMS_6LJY8JJ#r> zqK8$|&1=2NA*@nj{MMsO^z zGWX^q9q#BUUddFcgBN7+Ttm!Wf79nBA0Vrg6(I$}j_ z)r#wWUe~;G^7)jdFBT$mBYOfBe~2z`4QPMU$B*vr}Vr4&nv^8vSCG zx+RZN9R5;_C8h!)08Nte`Y|g{Y|NRJO+OoRO0wfQC3wcig_?gy-;Ic%d!aGuORm_3!Ib!Rp1%0v>Q^UUK)NO*`|bpaNrtt5`Q zdz|lQ_i5k%m;06rDFU zJE}hs#lUe&x74z7jNZgNjMt9nRf_T`+h!h)5kh})&~fujFKT5B)C&p*UM!@tl8>y5I}0mi9qelCRgbke$WH*swmmh#1M;YEcB6(i4{92>ltzEVUfXqiPOjqfy=;k3 zE~Ug53`F}s0{bIT@i&+Eyz<^qUJDgkUE*Jr*g~REhMem{Kgr>6`>5P!iIadW*P7FA z78B`8|6fqW+(ubYxmNWyl9jxTj)$F|gB~(L^rv8srq5*ICvNtSzpOFSZ({akrIY-d zGGJ@POZNk6__>ot-B9^v+r`HK`anR3H+<~Pe(6Q`#}&)r*3WiU8c|8DCf7^kOwrC- zO<9}C0Uud*CF%YyE3&kBXCf$m42mgg?FF?O-lxs^j-InsVgV)eMZStBT;hF7d`z82 z#Mr)FY+B=VGp|D_-6r$@QeP7Tez-f8o}S)gk?F2TZxy+try>m(0^zelvM2SJJ%sF+ zl|8b@Y@OQuld`LN%+{&h8Oq*^Z1eFA>CKQfjPf8{@I-8t;JH5Id058tsEp_HGoCNa zc)mE}IXB}uFXMS+#`C@z&wVnU6B*CR3_SBQ(pxjqe~^)0oAF$n@w`vQ^P3pMuD@-` z!!y!f&VX5(k!~{5f12@pZN~GR8P9j{-0A+GZt;?N^1F2l6w(K;dcKk6;)5uU%-Bh2 znWl+E;huZI>A)E2B_4kqHZzMVl+L`tWF%+ATa>GreZJ+u{7X3iNTu~Kq6C`hMXYIF zgC+1~-*t|1yL?dVNmTZ{7M1Fa5ZGD^Ouwceaf39vK((2DNR*>HQSmzqpMC!ob}qwD z%~@lcwkAkZu`wZ`L9w;A*uc^7>Ec6W3yc}f6fn}O`XY}<7IyxKdu4$w$4r&@rI~8LZQRv;|pSaDNE@7X3eS#1*ke)wJ>@t z0pGPHB{f*9rE+($2snOx!_$yn8a#eE5SFJQv_BizZ}T3+WH*DDLu+P>3{Mx!gtwnf zgeP`HK`BWL8cQL~-X@PHH>SOn&tO3g*Q4|BqQ+;dTal_&p{D-O9TR21J^Na_ zZD03`wSq4a*Ja+=8I9MMR!0w@%wXJKr5nEtook~F&*t(ZYdWhS&*lm75g~$5=tWA~ zLfF|98kIZJLf08oa4|nepHU(r8z`^YEQV zaTQZH#KbH1bH!ATPT12!#b7eAL}2lO!}Z9`vQYD*lL75n4tm|4FPN`0S-!h0BL#6@^{3M>DIDMNO5z!Ihswd0ZFfDM~S*)zBt_I*cUF zuSP`iW+0mn?8^mSR={Ho9)uW8p<|^pMY*fVd>S%ljg8euaQz5EFAHHz8p1#d&9R~z zJ1v0lHH>W`>_O&i2Vo!((qbK{LKYehYd_tsi;(Buk{J<7{3kA0G1bX1i-GFs4{G#9 zifSjul3>iB=-bFJD+#x!XX~{ua`2qo-hvkBDz^(6`MEGE^n9Hk^XYX?Zuy=hPd?u- zPyx&LPLOBYcL}nSg9-pc=zVjs!bEWW0{0LH_s6?9IoKf6$m$9!2P>|(a*(&1E&D4w zII#N~v!T&utWip$aJtPvK-QS|$q-+GYkR(RNm@S4De?JmRltWi?3D51_%Np=;6swW zh3mSrremnH&?A-)Yf^lu{td83Ug60*3B!sfkO}iMJ9F6qJs(b=X?H1d%NaZAYPrgy zM~|Js{w4JXmN#yzj_yj>H$<)JxUEFHXPhmd=yJ13WWmA4FnbNJ1Oj7=IpoihE2GiP z&A)O|JwJMu7mtgvys~8;VdEpMi-wiaa#IG|fMh91>a~}9{h|MQ5v>JnEJt@#0%PQM zeYqb=qc8n-+vmj0iP9grv-B_LHVD;VTP{-dQK;z=6i?d!WNUfx84s$U77vri zt^rhO4SnE46;F;YDUnHD3IKHdvUW}})C%a5(34&XSM2p4iO;CI@+EvmOeFGfSf7y` zBDJg}L$&NXF4evxFhWg5fVcY{Czn>8Q{1p~b){+8f#&k6T78$udL~Gn#A6?Wm2vvs zxA_XaiCPYszUbPhK?+p0eyWUmuW_RgLE=_CBvHmzu0MQZaY>kD1WCPKlb?uwk%o zHyojr+v;H$L^mT0v!rL&t;6ggmYCh){_M`V{G*|!J~Ym;$p!d6=8&N3`pm~8%xP>dtcXVEEG?7_8w>epg&u2KiI zN80fQEl-xxD#t3J>B-V(ZdL&Xn&A^>RTY#~HeRh)ylstRH;oJ%Z59hyywZ#B=v1oH z`8D((u=`d4R^a!vMo-~Us+*`d3D*{>xd$m#OvVSNuzglT32dLfXm5JC!M3c(TeH1H zKasy#ezQd3hXUd*u&nHM6M>LW*I>{29rxiFuCY6#0$Y%bM_ICPIfoSd|0UY9C6NJRMfl40rAM@y@@M&qAPbwa){7#H=I8u|$Dp9EDJkMyRpoeGNYzhk z%XqE^OHvR=!+43lZFm=*C35*>nc8kj`nhevs?MbLM+& zc9eSC_b#X6l|^3cCA*+0Fh?^2BvH#6CIFF7n)QhT_pKotY`$IBW^NbaTK;X%olr|a znDu}n)g!WMC+D(V5ML>TQN+@a+Oc~{YJyHH2IvN=qilwuWAflAJPU~d^RqB|e%oM2 zPP94k3h6jFSyX6arq@y;L)X+i2jHyaHvYkn>_^5|wS{g)70@tY{&t6IhK-E}4l;s9 zbBj7noyAg4rKnBTa0U4#G9}_U1;n=Ky`hMLOnIiaSSA2w=_LX=QYj!_kuCcFTcca@ z+(84ab|AoOmV9&%P>2x#o_9wiyfJItf~fE&L;z|Ib2d?^S?+(W)zlBsu?-$s%Tf~V z=)+NfHay4MXHdWC%WWbdQPDLJl>$A=WvU`y0-8$+nBwO-K+m>)|bRh;+=HQfv*>B{|GyGqUM}q9Sk}Dx3EWeYL23c6yn=z^KOH<{?1?870ASnK) z*%Y7W!9b>sVOEUhd-{jNo6Ik4Dzb>rQ^F!Mg8Z7flNpoQ{Et9Z+`$xyaO1XVC>k59 z+s7OTzcS(6B$s@mu)%Hy*^7a$L>ndhexq%sYBz$6f5$TWWpY_z)yxEaQdU@P|4WM- zdR1a%VZb5XJ{|6-ghkbhUs!H-7gs6_Tu$+V+&%2tMT-$c5T1HVeDVhDwiHINJY#MH zi7@-bIu`Yu`cdXwHC}xsfCZ` zrc$X}RM9qd-j94r!yK9~qNgvQNb0|SC&t?UMu3AfWjlDFJ zv=PA?{00kUOm%7%|5gCyY0(ZSQ3r(&TnNg$f?_!XDCu>$mTX#wdku#rf93R4PwViJ zmpNss1zPlye{UV$@mlgcL_9seOY&o_$!R{nS=R&1w%>Zke|BNZ`^vA~yE8Z!p>}&( z$7gMC$EA=TbKrzv9Y6gNjW^a@qGEG41v1w0Ywde_9e*)^;m~~tU<&Y?uV19&92HAc zOhS2!e}Aad^g6z1BJdOW)wc1%7q(NTUB~}WGM&gj)Md2T47-lsoeX~+e-JRSz5idY z<8Pxy%g;-s4g!6Vf4g147hb|`e)Mk@-86SvMdWO>iDHPvtkL{%;y+qKd@6K=1OwU8 zksDgfh4DH_6wpP@3)vT)&E=FK%;ub4`QeWl(5;Zp@Dz5cM zCDZji$)>4INxChzY@=%8mM zpPsLC0JL%KcpGHCnG(LkYy@psQ-4-k7W=nfmR{da`36U~1~TR3ZV|!Lj_zk*cHiLk zqE&9?A9`|lNr~^})-?`)pYqN4hQC@UfukEJ9DbCyy!~nEeZ}a=KWg4(Xr0iVg4lNz zIwtqR_6`xYVkx9%TwBDi$h4+1ssY;dfOS?-|EK%_xyi02Col3P)+{`Gio5u09lgOWCkt_U>ess;> zDs!ik)*D^M;WooMnWvKBY?=YPGv6{_r#!~vSfk3ZEQ-K>Of7Prr(9L5_@k4#U@>sJ^>wBy#cA-z2be)r#hT5z?pYnKM5Nu;5 z9F9L%D*21=NI#?bs_{wYOt^|N!A$u7##5aKyO`r3n?~g)SP)NFtGj`ADI7{t&Fle= z1jnr>i{WC_7@)WLf4-S}p$3TqQWQ~KrQsJc>s5c6npxg8*)vf3wbZoMA*EVu8Fl*! zclaqj)9MDY-tx0}6iEFk%3!T~v^6x2N@L|_|2)2#T-~r=(&HUFQLH*qEbJVS*gQvR zm2gt2zG=W$7t4vgx|@oPo2sL)A3r`iiHY(BFMgYU<*<#--?{yVGt&2MTvnj_Hu?`K zzMlIwma(c(+?CPIR6FUakHoz|Q{x_%=%j zsSSe7UTfH;{b(z{vEZJMp4*lEL$vZnF8lC`ujMRjy5Yyu4bRdeS^6~x&~PL+x-_|h z2|V4tuC=@4=THg<0)r*@8msu`!{H65H}*s$)y#etRt+}GMASs#=bHt*B5Mn<8!nLE z5YgWqs~u%-$7)jB@c|c|!flYR&7CL=&K;*_m|iH!Oi<3+o$$P2jF=ek1I@3A=ls{) zcX;ijR(`Mcg1JAKp}u1urbSJjrAf&$w3!NKZqgd;5>c~eJc<5##c`iyqhbxv6zXT8 zucIfaUk8cReP)V|iVrLW<&(;#jXwxId0L4y^6THG*W2f&>?IYc!ABOgYqG;yq_BLw z{TrvZt({Z2&h^#W#ZIYJRFNj7V$)n{k7)ml2pZq;OLrzLt&E9rQwzqZ~S*#8zq(Ne6_q7wWY#U~Rpt39BjD@{L`v`MIKs+4YLhCpGpa42QLzJ?LHgg@B7s@iZ5ed8oyq`pzh#u=<|}BcIbSN$ zY_R@ZLVES`sdWYO{TIvWMDCymGs18;X-Xmj7pGlD1a69773;GEFrBqS(}$bUC5Zu- zAi<7q%9qRM%u;*c?QW~l?b5g)Qnii)aa4fz%ru8>^b&QGrDzZ#_PUGPk3(Na;ty&M zsozq2TLQ{q8F&E?r-gJVyGL(G3O`dfV++u_p{5!@Qeq5X&OJ{b;;oxp8^5U7Rw|Z$ z6#Vl_aEN3<;Tn05tgaCY1bT+>(>9kM<9HTKAD>D44Ytd5Rl zEYYuPOCG&e&on`S%ifWyy~|!r*Xwm;x=;O6-hAJjjP4>17-|id2yH6h-k2Nuz&LMe`J89pLR*S zRGLT|&)DBh1z%7i-R<6b?L-*h&djWEaOVmK=ryAs14m;lJ!LDPb8r|mZ6SU;HQ}Us z1K3L)Si53w+DNI4&VQ-b46s@k3>cjP`%S26r?0w$wX27{Qec;ansOd&%dhrp&s z?sgWvc0~WogNE>)Rdxv9@%5(=oJm-Kq zWxvv6_QPZ!qU^_e%#M@2qq2X|Q+BAC`3RmHFi%BJrbrN;_84YJ5Sq6%{c%tKS@ZOPKq^Md3WjGBNBX_yN-!IV!WD3ghENk;m&wfkNM z{G&5y63VDska})Q_DMZ2NET*1@1A;YOYV_+UXXknTPDEm%8a(QPnBs)=A@n%Bs1y1 zQ!2eJ*(>$DAo&+yrq|w*7qt*u%U+(pKV@5FUvSnGloa!j-p9r!<#8#+6trjX-n+|t z*Zf}pPClN4!M%k-t%D!GaV0GQ!GHJSH@XXUIQX4*y{Gd>^joUWU2{!`Jg9C4V&_!z zXXd=g6&mcnERQWZlZ!OY?pTq6AA4$gr}nzygQ~evo^4vR{nxu(*TFmZY2nJNqi5*J z8+E`w^Bt=s{ho1shy5O3*LS`V*=anwzRxWJTJq07VJvNZ_x}4X@4xKwJ}vc5ROcU# zlXGFmX`a7xQ}9fo7oWuYVs6RFP8L$j@>j;#ruYr~>!|Oy=f$sQJC_eeS1pd_V;U8C ziHX_fZzCO%N-fbB@s^g4n$u%n)Qo6N?jMfv2}ib&zZjtpzE54fg2Lx{8TSjyZ_ux0 z!+Rm`^VrLcTXzcGupUNj+&U(7!(teC>Q;T3E_8#gD4V)9bf}(lc#gHt#waG~cQ@z& z%RC$;z@j%~Q8tf5%LcudqcrR(zcsy5UNmy{Vf!mfqPvbi;64 z9?uQ#qPhZ3Hy@MD*O0U397&2oG;S7oLlbgWA~E53J=cq9pU%BUvcJ&nZfy9ucfAzr zWfg_>Qm&Uz{g-OJeCEFl)5}`_Wwc(_`7aaovfjQll}rJP%=Z~x-e+}rmn$JvzE4xh z)4aFuAi0mVdl$aGe{=gOG2^v#Uu_Id1$+SrHLEs4q5p(k^4LQfuleD(C> z`_ym~5vLDUZ#O-s&(v)q`=pks+w4%>NVr5+V_SLk^fK8Q8o!LS$uu%qo2(!;65It; za}K_T!3A7NW@!pTno7$p6LM#;_{pTH-ZN>c_Y_Tk<-b^(-r>Jkn*PRru{6Egf3Y$G+m3&}pN5||u(v?>8ENe01BSvIpFG?oJaa77qQEAsRv>>e)MUpoG zv_>ACEZx_|V75P4^Q@j`W8CWXH#N^nB9g}f$jlxg6qNdQc-{*&?PUowR{0;f^amj*R6hGU;8nx#4t+gu-S$VBN;Xk^|#O&=M! z4p^5Ig{#(tnnsYS3KQu@B=)sXn%t!?yQtzUTk#SqoGFCCzACf2`l|>QIwDnyd<(3?%BRo+WQ)$6rNy%n{V#{zd3F3Hq@RznWCl zDfVJ-so7GC^k>8r6?^eMUW_lj#*bqhL(S!qK9=cuX4R>%Q(QT`MdDZT(bNlz>ecnR zkjsl5K@ihsLURO<@1AjWu%bXPPA&qb7S&>mhKtPafD?%~P@NAlSA7!QlW!*%n}^^r z9dvJn&BC!o+Jss~sbY?iMfayfk(o5vmZ*f_l zuhhGI1(yO3UqP6#gyu^-)C2G&pau*#>h=#j1{TRt#rpE_kzOpPBpiE%FEkgK9cT~` zUi?hI!3t5J`CS@}D>QeLx45QQ?L=Z1((Z+{TjRIO_YuTvx()e3qzqni7azwyR@>|; zoNJrg-?np)Hj*>jdFd3to$%tpwt@QCq*DGb=iN5I0@ut`1GqEPu^nkc4V0@#l}Jhk z^y;qAC#0bdve3gZxREO=!ixv_(ASvbEc7<}ZQU#E&}>%ZiYN4ksg;pxZ#190W`y{v zl>Kt{8w4XTz3I!@ZiI4ZitL6D0PQz*4o#ht(G>R&iaFWxsahCneqjqe$B8sCLSr1H zWJ~Ac`uJGdSQ$MaOp8TkqHSbVJ&5?N+%rtAtZZ1GZl!yNF0`%Ko64(0&4k4*s09I4)gX+95O=&5rEA?V^<;gPOrR_27xrQ?A!Nr3MsD!Ms zfi)PtJ~(gD;G!%KvF$g)OH3dxYy&pmI(wbqLbVlpJL4*`UEwp`@^|LhT30Oamwpi1 zcvf(in_19V$rkN+KBu>WUA_3&&PQ5`4`ZIfsRMR5(g2ieA zXLGnsI(u|DaqT#>pL?m~rP}n-OZ;M{ShQZ%d(p$3f3!bhR}d^Vo%~#&RvQUfKI*5CyQn3hLwM=gCkxE9IXEqrUUP@qJ3v!a<6s>91kptkdNM8rj|s zw2aJ#k8C+W6P0LfA`AP|FlgeM0dB*8+Xk57_Mq{iGf`S+N6J2_IjiREnsbt^6bwWe z71>1Zo&610?TG(}OAn_E*!vg=94(6v?420a5j(ASX6)rVN z7}@QgqF<~Nin+k5J)OeMv$rwwYu5zau0?0n_7vft>`(jY<&O2ISe+ifzFdUNH9w~S z?3nTq;AgnJ=uGQn_>-5Cd)vj@6g|%6iiCHg2ikfTN{OK0bdjPxGkFCS-TXzplrtyj zf6ed_r$yoy6v)&V8pZ&yz{({P+*519V%j)+>{x>i#l;w1bIXCOIb|`fnPW*}267h8O)i|rPR>17zl4{yKS z#0KlFQ!DLE01v;6^d=~oIRJtG?Xv?i{X-riJ2DVi;@_DGpBMztMDBgs3)A(@DVjxLk#P|6%;ael@cnJh-r%|ghQO1j3192f)k}~Eu2;>U ze@J#`r<0_I?`EIS5!yZfQ8Ghiy-T>yMT^3OyPjW`0jh=d5ofO0HoZ|;u(B#J??sph z;|Cu?f=nZAO^zjCWCkF;c;Ox-hj0ouPL z2tO99&;)c|{6xbpcoxK_J`E*?4D$%2k*@XK1%2-TcrYwDvM_=-IHJsT`6^^S}pn|ZUd&q=!8tE0J)S?${uDlF zf2uJSuVHNCxwGuMMu>J{7Ue{8SRGjo?8a>YQyL@mK}CX*(lnvw8->}_g_5JrMym91 z<8?V5ZV}jcT^>)mluaHi>u^5Gi=Ue3Wq)b5{fyYD`|Z4*@Yidorcy*V`64=5dTM9s zDY{*1?qxu_qtPVpF(028*p_YWOPy`m*1jY#GDiguI|UDW$qF|^H~X+*hl`5Y44-S8 zfI_LhukNW)zK;cFuQDK(b`X95wS9qw*7_t4BHR9i_~>DZrB$R5s2Nc+9377?3YI^W zohF6Sel20#3tHcRp!@W-kv-u6C9O6&4`^kiDq{g)2AZCwBfA&xpcY&O3&91X0?(ZEDf#ucz=1})qj;LMSUM{KO2zJE0gwY)pU zXi|PJbW-mUbG~hR-^ZW;u{G?o<*wt=sbBP&STvqLlO(BK@v(j#=LS)X1C71-h&=uI zhk7G9!=!yj+NAi^Xx9$2Q=N#Cdn(YRd`HtO|x)H zC{ephfZ7l(MsqZHZ+$lcC8G3&w}L{f!6S76GN(m^rFBW@P>6N~A=2s84l)67D5| z=VCqo!982&{X_1#TIr9v=Xx54oc*nujdpaL4FOdxfRvGp%mzdd|29I<#tf4RKI4Kub#>Q zu#fu7##5wfgOmJfi~hv}?6{v^yl<#ITF435zaW@i$ZS?5SDP?M8i1n&($25;FNe7d>XXAoKzCt7Y$S1NhY zJQ*Dzq^9sdvusyVPiSS(>T?Y6dezY4C8X7eCFw20yVQV2}DQVBYL& zKZdHS(-OHil4WJiAJ?aPMC!N6v5aZ}b_YGQ-zV?MlCQ*rc4n_w$wdz4Sn9VwklJr$ z-C&~wO0^wxh8G{KTP$&Kb3>p08B4guSCUkciWMfoTjbx`a|k1)N3W#>38P0 z0OveyJv5Dhd%VW$E3={@(85KiC9z$f%qvi{W0Asn_g>9{K(&8+!0td*B#`NX;l;)NZb*uT z$-6UH_82#Yo2*Vqjp0oj0U9q{fQ7;sF0>R`W1gTZ=3GV_k+6$ubD$ar-*3nXm~>b7 zu(*BC;NeFa@u$Qk6n#ZN-l>$vkbA%`)URH{T3UlDS$mm0vnQp^6Fjl#E-)`#l(E)+ zO!=D_GQYv(g&nD9(`xx(xw(*5QtMY0OmkOxhOL0sl-;yap56i!(|%dr?4SMgyy7tv zF1f;<{1e+~@|%6h73W`e`GpfM(d+KgRe8IJb$&|3o<9{lW~sD&%6nU2m2H;#kjUN1 z-h;%M*-oa;0WfyC7+VpTW`uC>R3N&GsiwZ)(e`&lzgr=B+B8hbdte_q8FDLCFLn$c zy~+rUkP*|b#Y)$2as2vPA=I?D)~RuJB1-WUoSoJ5Zm8++cnfP9uPxXv)V!5PFA=O_ z*=DkWE`d163Qom%I_Hkx=l(a3ofokAr;{YhX(o;9_d&_EZ%lBjgSDWhlcqm66WC>C zFdeY@XDz73{Y9N43)B|IQyo+0cormyYjUcwpt4Dp`E;)>mHV2;BA&mFA5HYT>25d- zC#LlRdrh2Ox8@#iFXfTRk-40rr&B;hbAL5GRc|a0C?RtRJ4iB`$UTybCHhw&MYwTGcBnZa zLctG3p(nGZJEM*W2yJrD1eU_VPh6B!Y^PFBhu+nKkA-xfSn5Q?K(ho&>=V<>3sXnq z{Wh3;t7*34rS98ln}Qyc-Zb^6MY_M%+@Su&3s`oU^|b}XxuGVgY00@z#7^XHEf7_^ ztV$vAU+q`~A0>o}+PUbwF8%ZvhZi!EdxJj;o>}MyuW&cPwZ6l`e)c*4K>OqDCD089 zKVG=WR<;Y$L6F9}Smgak9#C28XUJQ?X-1*3c0&Y#C^U%g0MNO7(lq^!^cLkpRmQOq!B; zz0VRBc!`U29(mHSnncdiZ2_)iayy`DP+wY2b*?#u@znAOwq=nkk$aRShaMb}^i!yiVOxCpMGvion#+K|bk$eZR-&SiOMf)u zuwH3I8wS z?)ew+Kb^!dYhhOqDr5Cz?zUraz3T9O3q#-6}+D-IFuAxpU97$5JD3` zHCVl(w{z#(+oh#seeK&K1kk>LSim90>R^V_4c4$`sTxVJ$cB6x%#rNy>hR&?z)*kq zaXRpcl34;v+JpS+k439WSq`i;6Dx2C?yc8E#iR3_){Gz^bMjo;U;VU735{FfHl1ZD zh{TyM;ay2A`jx)VH@m7P`drYO{F>X2k|4WD<=fnXhgpcaAlkVFCnG;SFLqW~p~{`g zWdtz9d6*_>;8|)Ba&;o5T!O>2^s}m0UA6|vgdYlJZ8IlwcZSJSZ=?ykLAz54@%(2j z^GT#0vI%bWTFp5B!WEdNA1X5RPAT&cDCJJv6EK;1%WC{4>pw}YG#f5Q={m8<+_08H z1mI-x9pPY&FK3%8e}Lfq9>ELzey+>QO*zy}wv!=^vm1k{zLlGYd>togN;Ma1x&tbh z&wv%End;!-XYQe`eO{2_uK?TH|6uzTFHBg5mC?bfz1k@b&+?xWB=c2IjXJ;I+dkDR zO=y*mBP?<_5Y2VeYM(8o^8|%5_M;S8$Ceap3_L^Y%nhIok72J?{Vm$tT8;jGktCXB z^UCO_Sj%V)Zj5~2JOyi>d3JM4%~_$_h^@}p@Ji_``4n*|T4$OgRL{2J@wxDhTS11J zh6c8xuC6+wq`K*|P*Wahx&IuftYJzKuz zILUH!=GI64J|z7wSdNK)qIm9`T=;;}vSK2x^)1K6_8y2G%>}O1OufeQpXNuFL_+vPd_|CcOByp8a)^(E<_iKaI^ThAV{hX`)u#Jhr-AQuW(WD;llW*kpe z%&hZm1}_Y9!%q4NVeCa>(m+?GxUkFuyz}`^M96 zXUg6IUS6UBb&i{CfjQL916_C}I@S!ukK1JmTF}0aZ@czFPP*ByzW?Dzk|BRxZ7yW| zyV0rP_nIV6##&}(UnMI!QDHb?uhh3v&9>=vCrkb=*m3ut3~gqrUcd;DJS2FyD0n!* ze@GesS6gh{Hsr5!)S!(2HM&P8AG;RNOHwS}b&V2gVoBPOd|RxQNfDNowx_8Bjrx)B|Y+e-H9(j z9Okx@Tg8F-6SCCbOn0I_mBN&0eTIt%mczPP<~JIYME+G2j>IOyBV><8?c-964`c~g zX%0O^D9MjVT=>OVX7Q)A!oso`Ph+JfTAYvs7)k3l-HLu?31?2Tbu|6Dh>*~Y+0wj3 zg@l?3h?Bw@ku@E=JQCwu#&{aEQ4)^ll=#uwSDG(Mv~EI=)Z-x|c>FUtDz7azn_dw% zYfK;52H?b3H>^^*;&jODq9L|OI%M{1V7ZW4TByf+QyNBt7OVeyqoUW! zfT}7;?OzmI3-(dLLFcgndU(#krIZde&*WdPaiN#FU6C?*aJvxf59_cu(XMzc$i@|p z;~E<84U`C`R#v^H-$EwxU{wJ?H(4lA8j_l4=KRU9t0i)mztw3MKXIXkzQl2W?gD- zP$_v!e|vg{w{(h4&8))xc96n6l@j5$_S^l2llzhrV6$)V(A$5oTM;7KBVtzh!S(gC zF#W<&{wkZ`GI4>%7=*?0t|*G?b0$Jq-kim8sThFUv;i;{OTauwgSU630>?SH2H`6`j^l~WVglR$hf$SoUX%mo&vpQV%a6h|! zxeB$vttwNbpgI8#HG2%GCg~gI8EMf-s811o+9r-GIcpGs!!Ap@38r+Pllvn6E*_5;lXp2Z((-0Zx(GwcK*{*@(5Q@QszodgaKABdnh?*as;4nys z%QQ@a!-0lN1Wjb^*ntHc7r3Nh{CF*w2rJ8>$wdSEX0^9jZnR!SQ-PeiZ}ZH=e{@542_x|v zS#^nuJ2wjOa&yyu!AP{V?*pU%e1Wj3IEibblzow(-3~F9o8M}z8G|$h{IgmN>AjS6 z@v=6?_zd*k7m>snd7ZQf0kPbg6&)7eArgDDeJ9or^YzZQ_7J}9c4AT0iN)7&U)P&n zjBDt*uPe+{aCrYg#n&9-B?iI&!)-9EuU-_+ROILmG)I4hTM#gGOn(rKEoxtC9^NNY zC5ox047eR;%!=;ZI5;m`p%RSSTKDYsxA`Upw>~+K%h&$q{KR#k3B5XXWt`%c*o%~J zvoh3YyHKAu?>$C%gb(WBiVf{Td9MdT47u z9G?9cV;H9VguQ8X_Po2A5y%Pb|#{V0getu)*zJ+vAy zX@}d`8I1NuIq@1-7c8)xCvyj z2H)nQF7*8y9r#}O9@h=uhcd)Z^`$?_a-efJMr%(oFkq1 zclt~J4t;;ptu*)^+!fy_1pi+2eFx*PGrr-j_`dbxzk~0?-AaRR!|q+^yO0iiFMJQ{ zhVO4Sb)j!*Hw?auAHI2^Ywy6ft*{HeeZMchS79A?rf+g%7knpo!{EE{9n!5d_!f4> z_dE#xz3BTq7Gr08tGeR5{c3l&+prXzE9AB?}hJP-SFj0D4pf4s2c{~ zh3}_JyY>!z|GHZjd~?1pzCXwQ>`dRc*LT5pTsI893*Y0ql?LD3uJ{gz;G(aLfSNPH z_#Eb7L+CH0SOo) zUoX)+%d2X=>_R@l5^83pN4!YCrMj1LM-IdG$*V8>YG~qkt|aV}(}(%}`U|lmuG`&) z@NzsDUoe*CO;FN{_b$2em zwN@{CLHjgkCvgl1H&!H;Q(}J8hX!J;G<|&yvHX}+ERlPcM9Mxl=O~M|+iwi^BiNA7 zzU)r#Y1c?chzqX0*b820&|6;XN87jGWUV{d#CA6=&%6j-b(kfu{r;3khDxnoXz&s* zv5!tAXC+So8XkqYEJe*-xFpgBlC$qQ3^&P}^Is3UKZUb}W82b$;>FL!5il5g{>NS! zcrYrVKh5{5)?OC&IP+_xv=TY~ddWXq&GF8nP@jQ@&z3c`bMHW27itFmrPdVT7>@m; zzHAHQ?t@?E%^urLuRtP4$dUH~Bs+K{4`H{pC z7}Zgtqgy}TUS8Qp;Y+)$64K#V3+}FBzP54&|1WdnySTbQ_%IIO&sJsCvW9kDSP`2a zu38?&F~_<3I0w(k=WI7+h8Qn4WTX1|Q9b4%cWPY05M#NrREU&qGgpN$W3j$)Kh<80 zw{?ebBB!N3+k|6-^6F!Qy0RJLRQ^ZCX$HhEB^n<_~gT#|HtF=PX50;J~QrT zd`@FnGsfqFOsTPl7VVl5&*=*eUBFVsu-sa>x?wH1f^-)d z*0q7Wq|_(Ji{FlvZQ}cLNDi`sTz(mvwgMIoKj$kEzS&cRAp6VGNU9B5UGFDpnGhrrr;z}GaaMLv)ciO4}#jwC1ikH42*At>~0 z$cj*(Q|aA;h7Wrblv8^90}+XRfSBlmMv?4RD=8;`i=^yuA0x4p0R*WgiLF1s_)`kw z;EnmxF8HUeHaYTMU!Qqv1&0bq?AH<-n&ZWX((#S4J-m5=9`wpYQt3J_^+5mh+ZnaM zy)x4tXA3QA*3X{yBp zpIb_!JWj__kA@qyY8uT{YZmRImxe4GnI@ma8&)=M^YPJIfBUP9-RnAt6RvDs)b7;L zv7Q_AAxJJ)X_4sB2rm7M7&sSiqODeozPjYlh9Yz^{zy&>)!DWSHT?+y1ZNbsK{f%J zMM2IjUt&9i51GKkKa-`9z3SwzsTJ^JXzqxCQmFVlGz?p9uh@VI0!5nFS%!SpH63a? zK`qii3yiIxK4*H_kt{HOa%A*qwI3JjYCOks658IC5zkM4;>PMhC~1e*Qq+NR=Pv1? zF3ixMLa80>n{(_hdD=kIpH3$bf40F4wW|5ZnE>v$o}g}{aGGHfxx}+M_BWo}NTxy` z527rbGhDCn{04q3dZ7vSoHgWxOKd}l+?J12rOCD5*fLDdEfy;A8v1pdTNx_5*|ruc zYbqH{x_&46zf<@%^VDV{_wf%k7%mi0T_j6?0y~mtJMj9AKVqpYsX+1>L=w@24LfOw zn;+UyRURpsqN+BqC`$sUz@p|)V0AW(cbpCjSx6QF;CN9~(V`rJiV}~w=HkR%lI~+t zWHDJ|{(_*22MXjwT(1|6GMTWm6P;{c2uEF zHSbIv>9L8c%fTUQO9Y*y-4|HuRVBmlv_Gt(x%tX6dT!lHfmg5N`9HPKX-*JBk$fV_ zDZYY9gn@rdIDdDW*M(ELPe_n57*>PXqE;`lM5e>rso5?v+|F=LxweuW^9$th!Wx}m zv9Dw3@>f1c^eXd;s?pCSp{79;nrl;ZQ*~(Wvg9$U#*E~D31E%)bZV*!{i)R}TV!7T zTI4v*&J4*9@Jq~)@!Yr=)xOA_!n!nBi*rT#+bet7RNI9rbIF*2ImvaD#-15F^~VG2 zL(NBR0}u?d7K~sJ4^6N_$kAz3to>7`SMa1_1sW~*uUNmNXW92A`^iYX+jPp=bcGDV zO#+>T-+Z#{5Pq-GuNugsvi>L|xVF*KbXvJhgEfoI{wP@qE8geHLc@VqPMrqA;x#@S?%42#TLsy*1um-Wa+XC?tIo13>|D?q5i|GO^ zK;Kvdy=v$B?b-X_BS{nPx44G-j4m+8h|~yDu{q3t8fl9Cr(tGyo-{2!wj(UI88%JX zV`B4-!ot@5+KPBWIGF>N&rYJ(%1pR794N;nR~krd^8H@7kUobun~~H{gTN2 zJjfCVL?N(aOd*6-M`f?3HgQ7D$4Sp{SeVrGUXG3nr{ix}pkwa0+S?)re+?A2p|GYu z_vy_Qe;gg&&wicE>V+=%DLtD-71r_;(0xk#>6#yIYPy1H#NMWKH~H9tSf?8r{-VqK zu%a(2i0_ts_baPU&`LIZ=C(D+)NaW(n_7!3DYiVv>`V2Onbc%GQHcugR(h2fa4!ap zdBhGQcO=+bJ{FrJwpSMxo3Hro_inP*>Zz1qnCDbna_2!E#HBNp_%Vs* zmM+DwOTSzBxg`CrC?|7{YW|nNbsvkz#Dj?e18=0{RAHOPb2H2*4UYqsTP-pz_SS+qU(9 zO4IMuopU?3`MYup+7ENva|@CQa)RfdqLTgR_C)aTfZ85sPLv45bDt5S@%(4`sZ(5j z#^mSQ8$OQ9X=)4AB} z>pzV&yZKMU%uYOs&F)1lBrg_&QsuTMd8yv1c|UzgoZWv0f``TSFE9_y7@Eq_vHRv>sj~jcGjQqPH37{bTO8vpX%dS({7-L-aH0P zf7g47V~_zY!dRU*{!)3c8cS7OFWdUjc0#?t{26Aq!G!W@6zN0OiJXG?PE;;hzW}^r z^tUAk5Yax^LG~>#=0S}`4PQD+OcP`|bs14XY}p(@%znMT?2}0L=5XUSD5i^s7#mBp z<>oNjjbyK5VzXD=p*YugxPHgj%B<$=p-?z>c3wDpF^_TwX3|IaM`B-iWp8_pn}L0x z#ev%o-q8rlU1O?%ge}XM@aAQ|6)F25jCVj07mb^!g4l>1gihHjW@k!AvOhwfhngl@ z=A?oKt}r!lp6+=bP9Ai~96!<<)+GCPVq$Nu8-$MiX&E;OUe`s~5YU{i*@K$?d>go% z_J^PpS_}363Stmc1y8cAF)xIg|H$@rO+2TBuLa!2W2ku^c{XAD=aUixdjA4t^E>+9gZxQ9q3O}+HuJ$-%*5%%`Yg{>dJ3+45#g)Sq?C#&^ zcvx_^-f0ws-aX!rheKZ8_;BoRX}6yZiegoNuJ%jtCc6C+JU8r@bXW`I%W#a#8?|3z zFT1`;D%kCpgr=c_=_`9CPf>JsF>n1o>$lX)e!dh>78TJAz{Bt3Kzv001F=IMG0{=Y?AKlJ@u?lqYj`!BxKbgTp~fhD@4_}`g9oG3f;qiQ z-Sj*gve~R9k?OWT6xfl^t#FKaz2q&%KFyyQd@t_r)9~HR>f$78#FX9H*isy>dOtcT zQnoIfIK4nyIr}jl3Q>(8<;5m~Z?!Z8_=0oQThS3-{4%(BvcjAQF{#AwDtsqqH*U+j zd~)^rMTD4s4c0jC7lLcA+UX$@w~Ia%ThD#3k)yWZ;EEi(&A&Iqi=heOe|DjC8RaCJ zOYIU>GM#CcC>#rM%P~Zs&~?)w&D4>2uShIZ6Z_)i`0!pEPmT=_5u56v?p_OXX>6qI zmF#WdvTqtc*$9@QX{8W==YAbEvCWN(^M=Mt_8J;1*|TcHWuIH(jGYmQ?@}MjtBGv_ z-Yy%5#)jtsZwTty7SXDorb*n%1+lmKK+)yB(KCr6R(>jU<= zLeoy?J+;@xUf&Y7bC&SUdikG|1Ie5bNes@eY2281`3W`aTc9Xk!pf>sO!>9whB|$8 zuAp}MsRd*bmUN-gVl9IR-Q8i^*uk;k+1rvMq(k^_xx&aJu~Q4`V)IXqpS3$(IZI)? zr_z;m;l>Z?%BywROX|yBs3DA3ErzD;MOCYwloXD=S{utLp=Y&x|D2u?tG`<~mb0t$ zPF2fgNe{um_F>Ue- zYQ#Fk6+|2>d7DjshvYipO`mE(eNcS(4qMc#Ey?T2N_+5w{XRVA5HHzFcqazDg1KoP zS)Zx$-1l4lv-V$mZFJoIefMAAvC{D!`>%U`MFTzWzrMfDwf+6~Un914CfDNA=$Wqt2r>UY6FY`O*{SNCQTt6YWi<<;FJMRzAHL3&*EjRuuk5| zlhdR2-4ZooX_&pLYL0f?BQ}U0%WXz*j^nHh&%%ERuc-?iD`p-L z59U{KTqHi4x}~Y$%n5%J8O2cv?j4m|u3T;c*s0vSzLY!|GX7VUiOS8g#(a)2s00i% z*h&FO>f*TI)s0> zkCp}r*$sq$+Z9~fKb-rG_QB??k|N(uR7F?fKEK54Dgn{_c{ z5(D})E+|j_n)i9tklE4$H%X4^EoDUaC8WAM>-gsYFJqpC?J9H2qr4?@m%l(2N^y+ScywEaTHI%S~VCnHr}X|DaxX8S{o`G$N0P z&XNSjv}gBN@)v?Qjgo4p{bu{`4MC$^K3}-Yo9c7L{?E${TxB+-kZTJrcWY}$<19Dt zK@OEDR4=zTchw_^VF+avGeteW6dKx;`!7^MtJ{`-L){kCQ%pTqSQW-9v$pQD3&sU zj!qT#i@cvHMtO}X<;9)k%i<(Ls7XSU{vb{hH^x~r&z9YCpO0W@*JWHca@`x&0eXj| z&qr=A+xWnuMcWf=E)OnD(B&fr1lp9EFPEtkV=6U@e*~&kq!|*uQZO4AdiHQbYW-ag zFBa?0fFcj>;+rn+p4@PJXm$Gz*iCGOa{d*e7*@@1kyrIfZUBP*omLyZSBRICE_&yG|bQ^I`&q|Ym*_of17iNq!shn~EqWPdgzbm$R_);tGi zgd5TL+4kVt*coOs7fSH)9zL=&*rLRjT%vo~@6=@f+WA%XR2_Dws&YyiKGAo6(jTgM zxsNT9{!qA&kSYJ_oe{#U5GrqTlhr z9+?Kj?`e-5Vq5snw?~%0pg#9(j|^-j4~E?Qh|2i($T^=e4F&ed2lr%5nZtPzsK6fS zYg97B9$B#Je`1d`zy8nckuiRWF7`;oPyW94$TcX7Onc;!Cq$OO9@)lBhVe0UfuAUQa9LN#qo*%o1PWDLCZ!?Iy;cgMPM|)%+(K&68ywqb!Yma>Cx7d?C zavgE7*dw>LI(y`T6mnhdk$M!O?2%@EyJ7UdV2>mo{#WdgL$Q^8dt@ITyV)Z<%9j34 zd*lyyW!NJ#dHGlEk@1wZ2Gh6NBR4HkM}DyK|H>X&*akh9uLF?(+#Wf}f$UC81f=bxq&_MGV(y8^N*YYfl^EOg||CUfj!c+TqQH?k)uiY z{{E3Vjw~rW>D&Dy2l*ws*dx9DbtVUU6+9UhOJWJUlmuUJ=H|@3H4-;e05s}aWCMd+wAZBBR7lAX?tWqk0q@=@-s?G?e%1je7?}yBOSAy zJ@WSya$W6_=a?^Kk9^HXd*qvWYW$zuBd<{|ut#nH>VK3yQV6s3WRLt7&gg89yvtQcJ=-H6G&_4_ z6OQ0-@sB)~)x{opPtg9u?2+F z9(mU;-F*MZ3IAXk3ha^AM8=eP3@-u|*dt%PsgfD?$bBSye|sbsZ_RhxBOfsbsYDlh z~+#dZSTSVuyJ@ROeC9OU38YREY9y#Y(YmZDI1WMM%V81BWadp|RIr9Hw?_1!T zDzg74Z7waZwxF&cs1d3bg_`t53n(pZpoO+{X+=~t4Q)!Zrb$dvD2iI7qFB(XsI1TW z#y9H&QCX?)$D*<-%BriZx>^=@Wqs_Ab(Q?TXJ&44liakWpz{BGyq`~dZqA)KbLPyM zGiT1soeP6VDI9T?V_Uw;ba=FGIX~B=a`Xl8? z+jD&6cf9c4{gG)%r}`uRSxF77jpP5q@sUfQ0eb$(P~ao|ky&rW^GAMrf$&G1Q1v+e z$WhpkfF+S|e58V?9c2EKXS-r@%)j${5fn2H1Uy7E>-=J?R3ZekrY_B@Lr3b za=qqG9tZxodp^nz`6HXD^?%;@$Q0ru!XJ5UMW6kV3o$8-^+$F-8clMsO0su<u49RqrWbgjS_9yk?BgHC7)gO5a?tMS~kwv_OgWn%H=@DTD_4bckgWVzIkIbSW zSB;N6g#HlnM=TIfC_Zv$n;@n7BlT}kN{v5q2oU}EM_z>S?%5x?L#BxGM^;I6fBlgy zFcw;Wiu zBa0sn`6Dje#>7WX1F~m-!h82e9ziT;BS`y)p{E?^~Wd6t!be#QANSZ0*BPVf5{H*bj2kxixwegWj;HZ!B zk$?R`&mTEGhU6FUoH@zf{gE35$tZthKn%$>D#_mcky1f2Dn9b)eX`=JKQavz`{|E- zPOsYhf3TnR)_a8+)bmFYZ^bkCBTr*O81hGs`V)pBH9m6NQq7QgD1K<;Bk#OODK-Ac zav=Kek6Z%d-LpSZCsRcEBjpm^Uw>o@jD^-8aZzMF5FCC~8M3tW?Wrjt+o7Qh`}y@BYXa*z7CgBOl^6#vgeZ$e#TX#}bV{G7mrc z>W_>^+MfN9nY{4c{gE-e5@J8=yDikf2T^>)0S(aeNB#-5jr2!O#O_Z0^%{#y_#hYZ3W*+_TXANd!=^Yg?< z+ISrJj|hL{zUDsrBWGh$7#kmX_4a6zr>P`+_eTZ`l2QK1 zq8O4!m1OV!$i{Vg@sY7AN!1^D1nzx5{gJu6go8gml6spkgL?jm|2jN_KXM`sxoUjm zisvv4ss6|?2q+XEY4ZzGYJ8;RkCam5k9-4% z^+!IW$a=^h$)w?bwQ>=5v_mZ9e{`QoBEmTuK=tI08D5$^H$5i|OW z_((a)IqZ-8c&=Vb?vIQ?%6R_BL$`B(WHY|0N#4fwA?BjvBNxCRQi{g`V>mhh{>UFO z73|#~ISreqr9V=R+ZcajB9J}%BipWLtNBSc-Zj9%T+*T+2ktd&y>yNmx9VNmaxtXXPWd6wM=s5fDkJP#2 z`6J^{XVnyF{E;kV5dKI8-En{9>?eONf8<*p2mZMGfhN*S86QdIlK5HukyXS;ggwh44p`WYXC9$Uh+t@<)!PT1Inp0Q`{&i~Hh_yo}A$(jVD`+Zcc39w2-6N34w+ zeZlvSA@L@Ie$Gu|fm%{E^q8wvqnG;PvtRkzr>F ze`Flx(T|UOgsmkJ{>V_Gc98iauR<34?~lCgi06;|{^77G(D);_BZKfq*3upKN8W&V z^y4G)YY+I0-m)+EVDyPNb)*}J_5t+i+P3-x9GroZ6=qnvo*V%*%+@!h_lNvvycJDH zV6801QCR%2te1)l4KFQkMM}KAzWg%4ohBol%`D&N$H)6{fReeqy(Jm#mU?X$4n$79 z_BqqmqLnWMZbVPPyXPzMF8T0odUgMw{MG#`5J2~(ck|P_|92Jd!2G2Od@~sZw4S^I zlQ}v*wi5uSk>ucod2zNjovBG{4|L%F+wq8A%>O5U zF`sx>U(8=7O2+H>conau_2e{ODo#8{O>1zz1>PR_QUQ1&zJ-@UwNGt30VJpp`Y0Zs zcJPC$aY$w&saPHr+^mc(|yThb%i8Wdey%iM;=mH-5$O!?h$8Qs^4%= z9p-NohZ5}Q5N4m!mf}=sigJ{$7RQ>kui=6K$(Ga&LQ?o5;$^(0)HBaSdvy;KrHB)s zsg$jtlRE5aPL3aKvRD?hofkl-Jg!{65%=A_Yl~1JhpwmjJA<9?Q zr_bviinU^BLrr=q4U$!Z{ue0owb6?Zqz7ru0iT0b_3db z*$aDsjXc7FpPNQ~FrPR0Mx;XuI<^QsETZ-4R~4S3!EnAq>CQoL>eTaEWJ{@5!}8Cb zYz8#~^#Zy=Mk4iH(Hh(wQ1eMUAh_^ST3}oF!M+{z3l|s>ZrRS=|KsPB?G78Do#z3y zo0|$A&x5Xdagba+x|iYWon3ecT`laMMB-YFPS*`VHBqTb;AOz^OQsG+cj2y$NTq&u~+!#P?q zb~%2~DHhA9*Z+7ndOhNj+$l8QC(0JbpsF8}$N3Xz(TMTrG9cvFiw$%}e)|kMF)Fp~ zD^d}Ba=Lvu&fB}JJUz7yAF{q3$mQpyC8l1o6?f(BdF6EE@yhSYSAJJf@GU+e3v75y zB-^q4MN@alL+~u{#f`68C*zwv&B>{k90j|r_C)a|i~nf*BI#BJH7S{x7V4Ro;Rj(& zgtc05wmD8&p?kVX<47EuJho_g5GVYn;1n-@J|euhw9*11v8IdWBO58h>YLKL0O`yv zZz;}M@1!4eChN!DllPG5KOzIW6^nyq?aj$$lbh30+b$;iNGBPt#%DB;zI5epsKF^K zKVSYnZa)VQjz4Z+CBII9PZF)HOfG}snaErhLO$Kw&Z36h#7|i7euaLZ=q|vb<$K{3 z(0Q!5!)at|aAz&}ztl_cjJu>JI?D?d_2uqilnvE`1?>u?5+*o5-2cD^HLIq>&n1z- zHmA}#RCuI0H8$K{=uxR`*RvXbViu0iItP8foCf!kq^ahMd;vkG7K z?Ed5)nu36x67;Srj4SSK_Gk52*wo*Yr`I;5x^U=}Fue3x*l_!z zDiGJp-2s$9-)A3>DnVsZ+Y)HPNJWxj9!aup2WQ|40+C!5Msfwd&D~Sb2e>Csr&5_p zFz6r?sRTGJq`MI~;w`o9xgS9tP4FX{oZ#c|m3KZ!NF1iqk}P!Uef-{eKmO4u?(TSKYaMLyE8frK_dJ6}+6pNUkkjc}82HWB}bdpbea-z}J72{!wYe`@jH1xTy(9oK75iW(NEyAl#$x9iu=pkIY-O+n z47Zo$|A1=!u(JvU^E1Y}|3C^I{%-d)I8@Z?k-e&I1$gRixIKH@`sIb}8D349?LQEbzM=pYMHa9a1B|N#-|okKy?@ac6~mrboNe zfZ6qcTHG@@C+wq@)S|f`UqFiTRri}3wo)k_@pmuX?D~$Lf3=P}M@*Fvvk30Y;Eno| zonN%(64Vb>FP}$Jk@3}sA@QGiKm2&&1ML$&Pw}bK`W#<%x1;f6&M*24-zu#oxlkq9 z`+DyvK{6`7x*~?;p(@GVxfb!aL^2;=F)kAfQkPaLuWLlp0^Hxs|%ECcc^oME~Qf-$5;V z-hcMCOcAx-yGf$^8(-aYJDG&=d5T{ms2q;3o`TUB!PG0yBl{xOd(R;etyZ=G)l+=6 zAY8Fy*Q>SyD;xldB zPQKYfZr>N(_k@^>UhiEE3rZAmNOXl{G!(`qQZNRuilPy>U!^7P^FF-K4j%jNLcJi?4=1sC!_}Q#4)xjYNr?W|L-0f8=N`iJx`7_ok&( zzSem-0UY&lz4xPw_56{t7?ST*_IbUxU673ONB(nuG?#5E$=>~ud_gkGAGt#%srn;X zpxDp#-glRg91ebeAiYR|%sYLhJ9~lp0q4h`RQ_wrK-uqBBeyvu5=Z5`} z4`4EYn&Ru(A9?KBuo0T7WdUwvqnGsH@}o zBgdR9{E-66qrcw!uP37Xkt2!PLFSMA1+v(Gf8^5{@%)i{uL_$2jX&}TG6;WU9o=z% z!`WSt-x<&PYql2rYX+d;9P{>XG*!oly41o`f1 zeSf5Gp-{+CG~}xO$c5Kn7*hQa0|XR`j|566S=b-Rxq?z^{E<&r?~gyyjMnMdAE}lp zqWqC565U^at)M;v?x0hdL>ms^gE` zS*#K6;rJ0V`iuC;6q0k;ANl7ry_DP^8IF|k{E=Jy+#h*(q3}m84KWuTA8CR?B!A>C zz!;7WfIsperh>iuBQxMDOMhevZe#qBTdyMGBddwp zLFSK?q2uhoKT)fMJtJBDrsEH&{I+U*y@r3jw2eOXRMLI|2>XnsE59jeKQB3T!}Rp_xfvys z-KKzTX>(ZtR)$mCzNPZO;2Dr<$;$6b=zDWtwU=VexD;PCd5&}0JrvPyj63`rLHNGo z%J)l`pFe_?rCMLZ-&^<`OftUnrCC_-D9V4S`KF?R_m-wybbaZ{e_P(#^}g6&v+^&< zO+}9g@TUBx9p9_mP?b(XC}=b_Czs@Z-JDob@ZSr@qM+q1*t(W#{S#!o9r<=$tFE_} zF25*)HIEb}V|e`{T=>r2ST@B6Elur%O4=)uQ30#5y*eXxgFn6fd{g%{gtXhb0x$Av z?1aI=8gBS|B}J!VOGi;n`Rawqkmrv=p7>H4Uw)_Ww72IKrEa*uh_y;9zLgqk?*gx> zc~D_~OL{Y!u`q+HR~1*UlBzVp{Z@RRxA{H#LiX?izp43x!u<2on-`Z96lYx6(y=)o z%}|DBs75caV)F&QJ&u(as0TjDim#{o@$so@6ZV3cu#?Y>w5NE2McPm*47Qb69$$ z!q6{h#52TF`O010?2^5~UE2hHGeiL7-^x+y%8$COmtp){Zp|?4dR#n*Ty~lvQ|idr zd4rD%SEDHo!v1(Zt_rH{5=j{#b_={aO3tBM3J4vp-HeK>cwH`eW)P zS3$F>Po`dSDekCWre1;dGHj9<~E&vDz-XM73FwxF!$F#LGA~2e?4M(5B>FZ z1h9Xa{`x3d_Gj&{3A`%(8b4I26p|*jBT6>(K_9KL=sS8gJ`3HQwO})f{p;50d zV?Wu)&(xwQ>c<@ypdbH5cf8%!s`le2!{cq*{}1EsV6?+Q-G9%xsE7VL^!Gni|NRCn z`m^@mX=|g$TR!Qcc>~;dI(hxq>sJZ<|0V$g!`5hKnugc5rVe`Vjw!Ei*)nC&dt;|8 zu=|~LHh;k9bT6D_s`I#U+gulLdfcW4pQp*>Zg#m$PQS_R378hyUCw%0@zd-+H_uY( zF$EfZ4u7M^Rqr?17ulUI`vRB46!4f9`ix!&=i z{5E@io#xj$j=6(WIiKpcFLX>YWu8<`S*&h{Z{d;}wcP1WcfG?$+&4HEW-hTextP`G z^Y|v2z@5|WZ4Q_koGypoWNF2$6jx$L)8}e?4sUyGRY(hHubqk z)0rMG<@cXOSxu%16YA{&yNwDqO)o3Anl#t!MGl`2jF?U@swg*6I?kZXr(q^4dnlhK ziO5OiacSC`?B4k0@gycBk=5gTdYS_ekaHUPI-HB)ljHk~-hZ!U_bH;WH{O;H>kj`Kn!DJk&8vUUUoIVGDnG`j;NJY3VM2H^y~mc!FmH`*y-I9mqH!O6$g zn~tqBam~d2xXGe(n>d4^GNbrZGm;fKx%b>L{Pwq=+I*%AXy&-d_!r_bMxIs97Kh7J zVfRihpDP=pRJ{p{VXocP>@b~)7U8yomV9+H*c@i@nJU zOEqWq>>4O_Qq0Z}_WWsoW4s8#ZTu!&c#KxlyT+K~x69YLI7mx4Evxun6 zgn`*^^5fR-V^kB)!kBUF#7w8V0sqUy|H|oX0>@0?82kq&(`|RN%u0JD{uNa+8!3ye zbZ*#w`ZMeNizL207Ehg6Cq0IfNt5gW);OLGVy^Mb#?!^-KMno>8N}&D)ioyg1nO9T z!B(}WOf!}+ef!QrV~K5n)9u&%ltTy88&A-X7~@!INcJ;saA>q-VsK*td!0=)uJWOm z3BQBQ7DM|aRz_q_n$1T{lSdBYQ0F-gzsI!*9)%n!#UP2AnY6mVX--l{89WOy5kks` zcqS;xVOjee&%*VfuAvc{HXRkNFGJ2(C5@-meALA-3o|D3B3Gn{!!#ovNpY*WVBowS zCk6-*hm8x}m&3oJRwqXYk1Cs+7GMnUG-w?&;KL3XT!w%X{@8gA?wdJ$$FQMivXj{< z>?{*I7XKWJ;6zX{jW&k7AfOO}me`LGejhf{+iz0?>iSRL=hBw&(n)PT>VHRk9L z1#NAJkqEvc9->@&XbNRDX&;NSxg_F~fWm4Dh%639Dw@8gVZ$^d8MU#%rR_oEcbd` zODyL(d>+Q{`2Q?zdNh|yCBw`(A*LI<8n)uE3xBNyfg`B>Q=fN4J*&lU)rTTj9k5tz z70m%hiw&NJqpZBLz-C+MZno96v{$*m}&ld36odJJY zoweGIe{@L}Cst4&abv?#pIwk+&avcV!X+7k1Gq>kBzF-zFeSF2o8Dqd#%&zDPpZ z*Z0WZhZZ7!hwIC!;`t5|cje8BW$vDD?;8(f;f*gO8c!DPfZ8W+K)0`39q z1b)Qgy}|MLz5PPqQQj+9+1QcH*gPyX96AuoS+s}*n2hz}`GBQ>Er6E+t^wSRrJ5ap z8CY5z4?ed8)&k~Z3HAZNwU>ep_}+-E0=0l^v6W&CU}HP-0p???(|FWp1U8JV0$dBY ziGE{Cs2TO^UIRLS8?mi>N-BN>j!pxASY#C&L0*6vfFs6$ zAA)01f56g=U@&bs(gT(Q&N~rw0Q1KOgWCZ&TF`!nfNmb@3D^j@5itK`&^;9876ya0 zfa|c?w+(PR;99_Z?7@Boa4j}w?*SY!1LZ@nHUe$}+zz;vU>WNB8}K(1@&L@oE@l^C zH{cq;wRjh56JROcG5ZW~58#l)Q6ApV%K@A>I~be~n2gQ3ZGamA9{_CpZ7}!_;2yl~ zYZ?JMHQ);{A1~6Ujf8wp3kEL(Y-xZz0rxnemq#Go?}EVqU^jNeb^*3{QQjz&*Nl1s zCNF_}0rM|}9*#!(Whfu;l@(}rz>LeG|3`waE1`#g8?OdGfZf*zgX}2O<3^MN*wP6- z1l)cH+7U4MUbOeo;QImSKj7LwKyHA09z(mBAfKnvE`YB*i+%vO=Orv+AA@pU3kGij zyzEW1H(>W$kR$r-i0!B^V9N*4Q`oP@4^bZA#=oLI=*P()BR}le9>7(ArGGS{- zgHw(}eLqEe1Mc|@atCbr9R2Kg(D?%L0i+Af$}hutbC{te%`j?6QgVku;2{cDwPbG) zSz#=bpoxs3SeApoJ*UE!ppUVXwCO1$W~L5aoZQN$9#L>Yb~=xA5#9Ot>pW|3kY~g_ zU4#$dZ!Pd$lwV$yW)=RthG0-d6Rs0~F5uAJ_(izK@Ye_&Oe?<#_X_^%qi`PrR~v=f z4cz=FoDub%7lj)ETul^iJaARO!F2PB%B6w095|Tk5UvKeQsAIs{34k(N_vwCV@*k0 zkx-m6VtL|>6mxCLi0LV5B`L{_=#lKHzFK;B0!Q`LmP@CqQGK;?*#aEZSBu*L9Mu=Y z2EV9`Z-JxwYH`U4qP|+(XyB;6iYBiaII6D}R{|W>SKDs$fGdf@xqvGK4knUc#7ir1 zQ=)KdfXk1<4M#gtGuC}qTql(eFh z1!U6@n^`h7I|Pw`9*Tt1Kgd!SqKCDTb@vtl04&iM$Zm4 zqCe?8E5$Q2o_&dD70@9zEhX*BglQ=wu1K7gVp?sqrev&A%2LcL224-MU!F80r7&UV zpp^WQ6mwBZ22vm;a!gA}E*s2<)>6>g-Ln6*NG`8{*0-=#ci^7*T9q&(WyFd^$a1+6 zoh+ezP>QK2Wdx**^5HCIAlVPd<2`S0@VAr)GFXv-6x1h8=kt`1zV-Q+7$ZAKm*Kps za&9G)ZY-F_*eKOj7UEet{@OsR20CWkg~0$kTtTu~O|n^KEJ?{&p+FML2P9553`)rW zVLWg^TE6=HzECSmF7yKS-3oLMFJz4{eVtfg*{Sl;aKT}Rbo+!X$8rDxl#^Acw7h{5)-;uG?=&VXULysld``E`5}8{g>hBl>I5$L3>ddI z&_5US8G7a_(iN#+33FA&n)UokTVQsOBnJ)l4$_|o9WxdSeq#CxmsN@ED~W$h(CsWU z2Ae>y4fL)%dvEY%+#63M{tA~Ttr)OMS#7)`@yY})&n@6_2j&+0D^Ic+Ba&eUew6ls ze&cUEaD~9}B{6PWo!qukTbANk4(2oaTQ<@1!ZwXSIhO$%{pfw@n^z=2H&>Hxt|B|L zLMf)sF(aiTA@K|JCo|d??HuY-XjP$OG+&#~g2C^h>qJM}24Q`bZLr$_mDkZ%P|HE~ zWVB`4fMC#zd*jC#-5befluDa%rrPsLo_`6f@Lq=gg|vmpSBiWi$u>NSe7*6}VMrW_ zO6ge6X;d$=M;qb%9Ey3+rMM@5VtGO(h7C9G`qX~g@jQU`7)#G_kxqUIoELd$Y!))* z_Q^^zMGsv;eHZpA6M9JY$*Oh{>tH<9N*CEp_?rTUq;(K{TigpC?1NI0XAYi$|1j~~ z*Q5FF14z?MY0!oy3^T-k3-E_x9>sNu^IwNtRVXQ$*KaMJ=iqr3<>Sw5@q8A4PCBv? z&tJgK{wAsmda{1IhV0PM)VE>{QPE&(S0i{$Mp+KrlU!uqrbaX-l)$2vq@>NDDxu4Z z0Yvx{Q;@F^`7rh1S3dq%0ap&(r-T8m)d_Hqr451 zro6m=d3gV#zPAp~(($Y!^lT}AM*7tWJ30@~R*?G?RED*k=HCQ^4--!VK$H8?kRpxQh3AHbzS6;GhwjM`5I-*< z@14kd8uh^sAW_~n;qigD4IG7 zS8B&lDI1DC%o`aD`bjo-iF`fiFK@40k_~Jg_|$q*Ge|a6zdJ#*?I`p`(t-8}nq}mh zRYs1@GgCT@iHGa!0Evv(Z!{eBR>T^r@mtd`(EiJfX9x+ynJOK82JaR$gAu%#VQ1DM z-|NsJV^#z&+)s?LyU@ov{rZs^;q_|+&Gd1>;Qcfv934wDtmC}T&q!%CBz{9SG*+Kx z5ls{RexUk|4+bB_*e`5M!i_L8rGw}3|9sA1LVRX)W-xdRlF&6Ce`&Bi`M}LbKH~@Q zsw4C;VY*iND19N)H=Yy>o=NFnjY)5fPVYtfwdP>(B1->6MEVk_)Hz5Giv&wKUHEt} zfX8mcc&;M`F+wK%4*cU8>BT0{8iKgvcW4~rrJxnpHb}X34NKgn=*o?^g$Xp%Ai%fg z?hPgof15!wrE8eh8A(H$dHDA&c!I~tuTuOkCBbJQ_Cfis7rb0S_EGxC&~51>&q(Ps zBqpNy%*Bx-59c8rWGGJhST^D^xHmR{77X+9gqd*gFb1U+rzAs*DcvTdn}QdIdAjoW z={`fc7u0lF@zV_%0@Iu;>5Yh=E(ht}QPb^3&?^pq^O0^o;&Gfm3Xa7|*M@YTsp(#e zpY8#qb4`@=*2hox4$^(Arn@eFx*w2kX}+YlG=92t1Rab8GTj3Gbfg1b#3NTB-6hZg zp^LLa^T)7VkmGh$LWiy(_c`$e`>^YkqstX?< zxGm?tfwW1r34hkBOOaZaT5a1U4%IjNMT5!j+>LT?LcIBE+)JB(5v*l7IhLfKX^3Z! zM*O)1<7l+465GV2{#%C;9vuKf-O)Z$6iT z%;+Y`=Mxb$`_1Q8klAytR`apAQVq&6xAcBvJWIt)$bPLuMc|ZCB2KoTBdlzV~yGzRR z-3Xffmgk{xed!wmVmykV+3$951kH{IC7+i^(Cjy#n?UoX^^(uo5j6YFXBwuTEssk+ zPl=$}Z$9fmvmUYVd1OZuU5Fry=vNEn%L$e1QX?2Sh+7n>+E+viS_kR8Hn9z*PX1PSQ3@-F-QRuMGEtVD)74% zp`&>?hixy#O3v+rtalp{)5!0M^=Em!%mkY0*k|y&^1Z?1ac`Uenkl;tT3J)MrAV^_ zhiN{AScY+!T5egW9$`OP`n|2mi7Ck_H@1ywyAw0v?M5EZ5yvb8CBRLL2Kn2yT|(!pA$s3@&CnZ<3eQL4qAat`&Z}D z$3Vt^{AugliZX^ggMAu(s&nEYq=vDIT|AiTyaluhH}7Aa)5gF*0*#-lFPDSQ0Loa4 z(=WE+UaZq$sbLkZ)vhp>i5R&#!G@qKR$?_Vath5NO$ZYjGdr;ymT~jm;Bxe9Js2`Uk@+M#GDVn46Ls(-pjNm?3$nGMpU~Mh1h10aqB5_5`fQ8H}UxK%^UZ z7;aBn@Q%SCa_i>3G1oBeI)n0VvSE3G@~>pW!->i>$%am&8TVHy%4^AnZ395&t|aA? zWTJBcQ78Hq!*Rov`HAC{3la;JTMg_qrQN_<6a;;5ppVr(ME{rtmQ2EMM1H0Uyylh}!8=xvE=Ol1)+d#I7`8I)Fq-;+# z++$GQOf-CGpl`Nam8iUw$nK$*dN$F_*wsemc_Vwos61h0Z&7VT3vgZrouPCY*cK8E zILrIYpxii!U7nz98OVNbRK7}PcN>*Yli3F2TM6$cvxk%JFg!nyT}NE)9H_i9h+Q{G zKJu@0$XnwwJL%A&7k}tfn6$}9yT9O zA2p2HmcZH*lno$}P$}E!GK1j)C!Q`hC{GS!Pa2e+!`KG~<*&oo_lDMlKM!Nq8i^JGy5JxL2DQ`ZCmDoC8)Q1}%N~di2pm_U{Dc z!9=z$k?y}Dg~6ji)07(x3<`B{{o;DxW-u(h)~Gx+8he41TSv2P#)*G8g55Ge**l7D zNmAZA0{fDcr;fmGUgf1D*j1}SfjW)BP~+%7Uupj;lY>A{#%@YfHV$F887#Cn(KL)X&qFg?Ny^A4wFne{7vTO*ubMRTe>UEsHzAll5 zh$TwbaqMw}^4M|g9fPvvIQES}c{rV2ouKSFmOY#>_se71wTa5Qf{hTUW&%pZ)(j$_%y3f*@o%6((mT1B~KEZfL2uPe$s$5Jce=>r3l=f=`VxM3{& zkYoNmKxsRUJ(omJw*_F=T8bbKbhA97+&b}I=EI*dr zHdMLtSoZo*!hbwec{-i_W2o|SI=d-Fc{`o0NL8*+XE&xQccinI`Tbu~mG*S@?^NaE zbav0MHsejlv6l|Rkax>(jzG%mMhpd=KaEfx8N;3*so>%Ek;-dh*p89P-DB84M=Cdr zVecNHTrq}ybcAy87!n%D=|2*GDOzjA4HrrMx?awU7Q!A|8Kb zB0|fLQNA9_UOQ$a?qQh6-A{w!cZs83HL}e{=)m`Bo0HS0yyzA)592ENaj@ve;KfWwV*> zH7bvq**c{j_%{Y9ms;3g2PpqDv;PjT0JAFOA^LvBxRwdhPP(H6>uNjp0tEjKw9;5IZ zRqWq}-zVb7y^8W&6?;_KVZe{)la%|a*oR3S3HY&QkaBw^yJZk%d0~*Uw}O2=xC-}w z8=`z&!FCNXmUMpws{6_htoUJ)rdAFQBezyA=xsbH@jrF>t`Hkg!um9y_m%4@UP700}u zfFJ9QtwZL!k5is0XRjW26z=yNHx2kL$D2Nx&F;z=iTf`yh{rYKly_&dTgLIqj#D;N zvF}ef1^BKLiO$<6CV=SQPaF;kmyD-0_l{ToQOO=1PeQn6f^zX}_P_)#gui7f|EXlF zPRhoE7f({It7LDUL>b;cNx6GA`|_k~Na(j@D>uw$>#|9l>vEJUX0vB=wvo(VpQwCP z$(BvpV7R)9y*pW1QN?ylUPC{=pQ3y}o3)*+;Qq>!m4D4-{o={HpqX8#rZDuG`m@T&w4xCC~b6lrHn^dc^<4rvF~0&68+{jI`{QTVM< zztwcR7(%*WxB1W2z$C)LY1BhcgS=*q@87P((RGy%<_3v zoOEf@uD8nPopI8s<*QI|iFhv71L;zKe--t6Z`AJ;tet7o9~JdGGwS!0sNa=QzwJ@K z&x!hdUexc)7D#b54w}mbV*cFmL+SfoAwyy9?U0anu>2+;> zb)7HxRkgBHX2cdZe&KZ=eqs9$zcwxsP{r?XiQg*U^u8Qj5`@Z^Vhq*qe?d_{=mF2V zSN`riN94yFQ@sAPeVDG@^7r0G@%z8>_fIt~Ok1K#M3krgH_P@el(0&|S_!=pwo2F` zVW))aCG3)LtAslw+$|xyN|cx;p-DosgoP4TNmwhPSHe~aJ0$FsaJ__G5^j}nhlIN& zWLL}bB{WHBmatI5DhX>P^h($&VTXjB60Vo9OTw)Z?vQY|gzOqwzJw+T%@P($SS4Yt zgkA|-CG3!}Q^NHUc1gHZ!W|OsmXNKHf4cZ83M14gKCLQ^5)Ap&)xF*V{x80$mav3pA${LUOqp5 zoPPOr^7-&M&;9awS{(Y^Xd}&0ah^M5{&7)$m_+ll8`S&0z!}9jMmOWs1vo{q*<{Ji z%+54V$ZO_5vo6fa$u#F=TE@u-J<%gM8QH(-HQk7d#xVn)tdaw!Nj^OUP%nKJ@CkV1 zvFb6A@N${{28rJhg-=Gol-~QeNO-tRpDyto5^quINc`3)yb;Ywbka76g#26)+3j~lYT=kr`mS5W@8bkuBjKPB z-n2;|=$s9@7D@c)QR%l!eBmEOy0-*i$H39xbe<6iel7xTuaWrOnL&LFZ4<^&A{6N}tuqf3lYaV7I_;;&dcJ z<#QYyMZ$Yu6$w@RohtrM0?*&S$L;ZOBq{y(odQqqhSF6n@wFEU!~+7&wn}{3MUnVQ zykkamx}xw~B;MN^nSLg=rBQm*vPk^P62CPHZ%q^FTQ81G|GdPTmq+4@hKck$qVSs~ zzN0NN{mH{c`a=HF7tTJA^?eF>0~^ImtwI3{MK<<>O#j1eA|dV7rYjpAmFji(UE=rg z0x+klK9>CiQg9sz-9rjXJ_Wh^c|1Ngc8h?_^v4YITFtg{nBJv20ZbZ z7RBdkjz65K^Ak18W0H;bxZ-+%J-qgd@S zksjTdUju(b@*Kw0d6z`9RNxKl2&T?&RQ$0VKaxfJ_n90&ibaptlO?`UmaEd40ek{O zCszG+`Je)L#tvbZaHhr^a%B3oz(aLou4m%lyMZ4A{ZZ%p$I6^jF>nkA zzVm*8=z@Kus|0vT-#SHTAj}lMPM7$uT7joC6zIBGP5)bgFXUK!Z<5oOKQtf-6HKD7 z#(5+!iV4*q#br|Vonht74N>#uRrpMXFFmAhNoRSDP%i8n>b;SS)B z!7N^lBSGx^dP$~VUn$Bhm-vwgWDxzt@cbMh;LihIuf7}O;IYk6 zH~n|O6Q8XzK9#{sMf?MSn_<8g`h`BIcnk29zER>^W%{MSbNNX9;rsXT><|Pn^!TLs zlODb(4*qoDsa!RVrpoi#IO#tCp7>GYcNm`dMF&Cr zRZ;ERA@NME%QQw9=)zg}{Tz?-BT;12QY3W-nKE%-TE>fv7C_4pq;M%bm!WrDsmSIi{wT@s%r$gw<* zAHlk%{;2IZ8+fu0gi8q}m@X}pl$Bq^JY`r*v=U4d| z4?NN7lJ>cRKV{4byk5Pokm+~W1pPnD^pC|!|CvnRdZ;K@)&CXgLT~3q^{)?s=j}3F zq_2|lIUEC%9zUnV!8?E_I_mo4B*|y9#CM)12-7)8biJVB=L&pKpxG!GWTJz^FgZrm zhiu?!+}R!_&*{MH@%bwJH{OoYKZjWOHAm9#kQMz@;@3)i?RHVF+P=FbUR^(x?6X2R zqI%^v#KCvP!G9PB{|)eD_gFx*i>jYPVE{=!RsV_Pvq9onnnq?+2t`os$(T$*U)^_L zQRz#3l}tEK>d$Y1r~2-Z>m~Uz{oTOp@$;0V(;@8u+!cOplX&xDQLkep{$CQ``aglE zvpnb;dXnICt6b+&+r3`m( z=|tYtnlJh$_b&Jan*3cOyuZj(Mzg4n6=pjYr+G zT`+v`>-TZeKL9-O-#JxiIGtfi*GqBIC*+CrJ7j$i5y{zcz?0sp>+59_@0aPjq~D8X z<<~Zzegsq3+tv0RI#JMHFZFY>OrIg~(f-L|6+a}>K3ogDp4{G#gU^7WBe|tD3p%RY zP63|eW7;6(qw4?lz?1%~`?gg6KjU;zvPsZU=iT1`ug8y}Kp%f{9Q+x;bA5YS@T1zp zC2`WfA=7t8_3NYH_z@k3_v7f|=UM@}CBE}~fxnJF#e1s2lYH`{ z#P62vI9TGpl=vObNAf=y4Nd%XwFvx~0cB zMACUk;(*z>+!jc zr$@h#_h5t~mHZPS&S?0`MwLNt%P{anfH1Jk=|@eLH}sda3(#21;T#^YmKE z;(nu~v;G+&7_}b_!+1gTI}Q=@|A)*<`zz9P#*rI=*OS|Crt0IT0#Ea6H4dh7HY-kg zJMhGhSI#>|$_yLgr2jysPm_AC%3)8O^ht$+erMEt>TKZ2o?IsN$|TFZ4tTx#Zi$2c z6nGPhy}kxsk4{n%;o&F$SLC5H3Fyk@cyX`&KNEPa54RXN`fB;mE%Aj*1^y_B-w3>3 zxgX2)=BRO#P18@G4m_2+Th5zQ|8_&1^smLiA6=|Ze<{Z!o*)&F&RU`C!8qw(2A=pa zeJ<*&+TqV7Ufs8)j*qiT^vm@C&*jMvNo4G1L5{r&JoPViA6va(hgq!pbnXD2^k<#a zAJrfGAWr)6)Ai|`4?OYHb-7?fvc)z_yt-fQSjo?)YWi6+y~GcmA?TRzi+eYvt)d{}6$v z^H1n%jzj0AIQUQF;P(Pg<2Y+iE{Q5Zf=i=Z8!=RCTqT>}O0MEy- zcZI;{JTkiG15fr--8X!+K(mX0*VEfaIUSiu1KtLn%8ibnd@Av+QR{amICy&UoCrMG z2lago)y}Mv_&w5Yx5~0#cSn_bc#V|*t%6T} z&Mhe2C-LiLe@T(_ji-wAg;9EbjKp`y^r{_L1-xFpK9uQeWqL)@&zURecl=fGqmKJ? zB)&DuFIoh=9{sCi`mOT?orfj;dw?gorOEhGtHggK(-$^~^mJ}MU5A_|`H}wWrvlAR z2VSpSuS{=RDd^lR0NX6_-a7R;N4N_~vRRKGw@@XW)am`L;O=@KIccqs&#C3leVoLL8hA8mg3k8yi78PtNGMk-OnHSF(e|0vF~R?Vtw0;}2aSdIMV zc@9?*q^=GW*vjf1?tn9}gj?}S40T-T4y57VC%e_N`9~Ls)2;ezCu2I4#bPURxjc2e zV;taqbWk#>(Y%mVDAVmY)q!k$l^wp5g?loQJ`!-VoYgt7QL!d;8cb-etxPp!_}Pp3 zM2a}={F5^lr4Tdr&B9Vm%MX-=J&^MLd8*ZSU8=NX)f5#M+g+s&m)GITE%rDX8e%mk zH+Jv^su%4+iha?Zs+^iqOG(VQe}Ic7x8L7r_c`iq-hi*HGUxQ7it@aAN5JlM34?4! z?}@a?=uK)Tqn<@m5$a$)+eoC))P7c$rGifoq9#T8r}^N67FW-;E>auqj0m0zXppeg z+XMD&%zDa-CfeY4&8b8m7TnI61s^P`PIzd^ZgMnvd`oPN7>*r27=y}Lr@GWRM-^4E z1XZxrxvX$S!i@G@IOb=+Wa+4bCC6zAsmnKbFONaXR&TBZLm$(XWO??cfVsb6vC^X(Mvdk5BWWWGw?h*oDr`n za;G&rUG=3-f55ZQXKxA%bVm8?X+`BWo4kaIVAWN4nxhN8LGAN)m67x8(6rzQk!?==S&n&bkmAO?IbS z=mi3){y=j>LuMVrc{#^f^TpR-3i5 z#Aah9r)Bya4VgCE!j={rd_}(pa=?t$wg_cYG0Y}~I(z!`YHN+Hrf6EZ)y8Zkr&ktL zlob<$evhpYoyz6FkD7{NNtu@n%#6xYZPrp*XlcnD$h4*!&P|@J+Kb2&-a3ZYwgA!3 zvhvCT@ZoN@A*PJDRJo__9E5;f&blSmEGusd3;eBSaJk6o_xK8`9f1l@eX|R3(}@-0 z?svRVHQry?)3RhW*?s5O>K&++qu#nGl+u!4v&8GDpJuNw_P7^`<}X7)dtXBZR%?kF zWU2MJ-ekk`Etug5R4#D3{qr;hPDGdz`7np@{`*nQ-c zXW7a;e)@~=gqTaXXXedEV|IH|bWGN{+5weoavzP~&KUFri6w zQ|&y@VQmmK?`68aj~2G4^^Gj_CaApKwg3jgZL6~f>KY5t3Psb(Y?jPyVu!Oj-RW{! zvcM|c6nk7ApOprzdRNt=CQbUx0?LtBGs9>1HWmwsR|mW0CeI?Nf#LSnFUyjP9)%8ML&U^JD^1$D6GCqx zE_shCwv-`+hxrvIRrEyDwIXY)(^G_iAVmZKb8K1-%FMC|VtI?OwB_-Y)gx`a)rxhB z3LEs0hq(lswB~^))r5hyw&OEG-rLFk+*t{ zoaE{nFGSh}Z~JUy`9hozt97u)HW$1W8->3q`YU^`wm~e}D(k|2!-ipwniE1Uvb&lc zRx8X~*w|_(<$4wwCMjE}V1v`;h-2-@etEE57Sp9w^MD3JlS7(LZX~o!=ZPsa3@92Y z;1H2W;;&||(Ci*08cJd*gAsB$>$$-ec{N(qcY|%!)u;n&dJ3SAnx$3ES{z?r#X3*? z^)vWna}c$LeZu&F(YUVh99sioDj~C+Lmk)e2!s&QP9y4*RT6H6`g&(QtsnHLwc*vG zesobDznKU_B^`ylia8E@{p{vIk>BrVTHsnDEn$ystd38rH!-oyPgIt~>Tk5&!#LW% zDIN7wifsB=eK$9+SeWwoa!(kbM= zH*8v#ja(wk+Q^GRgXqmaCh|yxd!eceTbqT^q*{QO*x_%*+upBn0438Uhq43zp(A~A$t&V-gBnIzl9riUQ{pXV9@zY1N~jS5mN8P0`++pOz$HrZD+Dc#|Nrv^EB9ecK1RDv`Z_${;4Xu zyJi!WRkulI@&3q19!s9pP5QuR4=w0Th#aC|EsOA!76ux{kf6)rfzem)B-T6q*k~yy zX%lP89Sh+rLPzUyp=`{eNGhEDkP}j`8@b|Jge)S{%Z?o{_;#CEW^ z%HP}sAFj9<0Vlq@lR{?s(AgSifTDV{>QFK!U7C)PuZ>T{5BP|AyWaFlSXsT0Up=?d zR98!`*y|cvLC>Ln+)rB^I|U7QqFf&G(Rg{5s8L=F0MclOnn8Dp8Xm5-^6I0p4Nsrc z(MYxJxjf2EZNV!NngZ)d@Bp#Et&XP6`&b7Wo5iHwI;W{g+pC3D?2&bh^d#PD?FbKC zoDp6fy5WThC9RxTF@uNkE#@)=+?@yxAe5*kgD1h&%2G4e=?gU5UGVX<@h2Mo=beM> zk8g$zIMO#meW^pkd66OKI7>D1am`%#L+1XN7gjpk|@@O(OUKrnE^rxeRj7(Iid|s_76$0x+AdFlI zt)_=egou1*MMOTmt~s;Rc2mU_nhsIb#ndXqgw5q#fOQA7vdck3Cqf~Pa6~>k3x(B` zsRZ!@1DFLX4m=o0d3|uplBOh8NT$K_ zt!ePGqJwJ`^g)Osx@pu{DF!}2wmLZdjdF1nMDuJmkMPmdNwdxp3x<@CsgtS2x|6j! z&j#MIZPE~kK$BSBI^eRSXlcCz%zAjVx2%TdSz^5&OWL$rUt?b&0zM=zGtKM$OR${Z zq^{HAm8U?X#p1;-M7-@nKO#+YW>y|)Mzu4?>u{VC`Hak)Tz<&Pp5yTP9DWe7)1H*j z1l~%a_PQo-wDN_`P(i>WSBd%BDA)|muONT5wIk5+8oLKm3k4QqK3~P}Wl-#*azthJTB%#OJzADrYLv%s-Ju0z?ime;ab!vB}75-x_GXQ0$?!pFt1q=65%42P7|W$ zv^?z$FL-40EftYopO8k#;fe?*@58JIOB~LHZel3T8ufvZn(7^dM`IZn<|7X*BG$6G zN?Chb=%%+EtX8~LBlmE~S07~L6)$qCFTw>kpJg6OC03rX^0qiuWV6J#(?F zvdAHNvB!sy9QUEZi^HK#2T4VG|C+cBtr%cVNa24N^2nAH`&q*+*`v2ACY#VP>~|7t zaW2n71dI6!0Q_NXztt&ot% z*m#5?)KcoQyY9|6u{B3%c0XP+Ac6KFMxjow_jR$T*R4_Xngs1X@X?M^w%iW4R8C?pS9)tF(q&}yIdmL6uYTN-DlwS zEDq}vU&!JF!-IR|nq$A0kZQ0v(Wga!J)y<%yX)uD=x2!?jzc@)Cc=TM#3or0u(K^} z^0-AKVG}GAJiK4dSv8CHVG%M;P|Fg#|Cmy%5!0L!pM9|!7u2>%^wubK5({zct6z%F zTw$}53b^6uS)g=2Y@!W;ay%oK1*YE^_dd2ySBd1~%Fto$iYSq<7Vx+y+;(KWe&3Poi8D{=$#OZ{Ux0gA`#8BEGqdUn$fL z2MEAfQx#m>3BrpX*%+|Z2)H_URO3ZfTUqP^jC_TWDD#6M0BLQQeq)cC*mQ}V-S%b| zHoJSVJE-2;1PS!$g&m+uk|R_d@ab-}7OlMhU+pbYZjwsL%B^3-$Dl}u33DKsgv~U`}-y(DP#O0pg$6dH2AtB z3pyZ=9_qYr6y;wC*0ixTIhn5xIL?jl{4eOKal`L3Ah4XveZ1;=VkpL_IeDkUJ0U)M zU4Ul&bE=(lbT)g1v^sir>}YYu?23+feaVtdyCOnwGpQqCzPf^@%hmykf2jkry4xe4 zwsX;oWN=$mHwwRpvl))~>B}-@v~+m>i3?*E8a8D~18 z23F%dj`O-TIP1|R7cVX$;s_C-dBiq;+PLR;_Pmpyw`BtkAB;O)RU>FB%|~d) z_Guy+IG9Q50KS#M@6bTBmBCCQA{suFy;=?iv+eknwx* zBDW1I#MXlSC$cjJdA%EEd`=Po7=Dy^z;HzJnvRBNqSXc8|l(9p_FV;^J_t;W2Zd@H?f zQ{i8jl{KxX#8!RkwCdtHWmPq{)5>Z}ZL_SWPs^oa2zU<`y1%auU2#_I3k=mxAC9!u zSq@M;*1kl);>JJ1{j{f-k};|2y~Al=NygMW;0bR4q%|j&>0i}Y6Yne|KDC_=!1Pdt%_ zG(`Nc;Tc|cv{MFo=XC{GCPlMYCT=tTUuoy>8buVu@r8i+qeuiRO>y9*iF#ZmraJ`x z2P+{Y!3!ELDv06_VwXa}&MO2@5K`D&sg>7Q3YJ2eH1=0%m-9P2^Vz%%f`tz@-+8n1 zotgJ`Zv%VB%s1Sj^5a4T35K}+#l5MC^{ScD4o$?mp5h#v!?H6>(}6_fNH-2 zs=azu4fla;ek={StW#dy+Egs}uCIkEl?B<)kohWWUHzKP z2IM~uKG=&w;3vXQ&VN<bq}BXCdKlFKeO!TA!P8RpML^5n_Ky7Bs;XXv9)AaS{S5fyl+TF|pY^{9 z{&hF2K0R~+9NFe-82`ME%zrPPrHKEg5BO^-HSGQ_N51}Vb+Y8Qf(L$&>yM9x5Wu~} zpW1>nK56*-d+i_waQ~qFnfTpX&IjA` zSUTqTzooX`0`b_WcB8)wnD?9b;G7n3A>a>#4`S!hkM4i;^Wxk7!KW`6xZQj(+(QdT zfr}!W3H*h}zJ72V7lzxAxUB1=xMlpw6DNS(v|;eUrIdfu7jnn_^v!uOT7_)lC7m~4 z literal 0 HcmV?d00001 diff --git a/apps/StripEnergyThresholdFinder.cxx b/apps/StripEnergyThresholdFinder.cxx index 642715ac..60661f3c 100644 --- a/apps/StripEnergyThresholdFinder.cxx +++ b/apps/StripEnergyThresholdFinder.cxx @@ -257,10 +257,14 @@ class TACCalHelper bool PassHitSelection(MStripHit* SH) { - if(SH==nullptr) return false; + if (SH == nullptr) return false; - if(SH->IsGuardRing()) return true; - if(SH->IsNearestNeighbor()) return false; + // Basic sanity checks only + if (SH->GetStripID() < 0) return false; + if (SH->GetDetectorID() < 0) return false; + + // Require signal + if (SH->GetADCUnits() <= 0) return false; return true; } @@ -268,65 +272,171 @@ bool PassHitSelection(MStripHit* SH) -int main(int argc,char** argv) -{ +//! Finds per-strip slow and fast energy thresholds for germanium detectors +class MStripThresholdFinder +{ +public: + MStripThresholdFinder() {} + + bool ParseCommandLine(int Argc, char** Argv); + bool BuildHistograms(); + + void FindSlowThresholds(); + void FindFastThresholds(); + + void WriteCSV() const; + void WriteDiagnostics() const; + + MString GetOutputPrefix() const { return m_OutputPrefix; } + +private: + + // --- Configuration --- + EnergyCalHelper m_EnergyCal; + vector m_InputFiles; + MString m_CalibrationFile; + MString m_TACCalibrationFile; + MString m_StripMapFile; + MString m_OutputPrefix; + + int m_MinEntries = 10; + double m_FallbackThreshold = 20.0; + + int m_HistogramBins = 2048; + double m_HistogramMaxADC = 4096.0; + double m_NoiseSearchMaxADC = 2200.0; + + long m_MaxEvents = -1; + + // --- Data --- + map m_ADCHistograms; + map>> m_TimingCounts; + + map dt0_hists; + map dt1_hists; - // Force unbuffered stdout so progress bar updates live + // --- Results --- + map m_SlowThresholds; + map m_SlowThresholdsADC; + + map m_FastThresholds; + map m_FastThresholdsADC; + + // --- Helpers --- + int FindNoisePeakBin(TH1D* Histogram) const; + int FindThresholdBin(TH1D* Histogram, int PeakBin) const; + int FindCrossoverADC(const map>& ADCMap, int FirstNonzero) const; + + // --- Diagnostics (temporary restore for plotting) --- + vector m_stripIndex_LV; + vector m_stripIndex_HV; + + vector m_thresholdValues_LV; + vector m_thresholdValues_HV; + + +}; + + + + + +int main(int Argc, char** Argv) +{ setvbuf(stdout, NULL, _IONBF, 0); + MGlobal::Initialize("Standalone", "ThresholdFinder"); + + MStripThresholdFinder Finder; + + if (!Finder.ParseCommandLine(Argc, Argv)) return 1; + if (!Finder.BuildHistograms()) return 1; + + Finder.FindSlowThresholds(); + Finder.FindFastThresholds(); + + Finder.WriteCSV(); + Finder.WriteDiagnostics(); + /* ------------------------------------------------------------- */ - /* We Include this help menu so that users can define their */ - /* Variables,, without needing a yaml config file */ + /* Helpful ROOT instructions */ /* ------------------------------------------------------------- */ - for(int i=1;iDraw()"<Draw()"<Draw()"<GetXaxis()->SetRangeUser(0,30)"<Draw()"<Draw()"<Draw()"<Draw();" << endl; + + + return 0; +} + + +bool MStripThresholdFinder::ParseCommandLine(int argc, char** argv) +{ + if(argc < 2) { - cout<<"Usage: ./StripEnergyThresholdFinder config.yaml"<(); if(config["analysis"]["histogram_bins"]) - histogramBins=config["analysis"]["histogram_bins"].as(); + m_HistogramBins = config["analysis"]["histogram_bins"].as(); if(config["analysis"]["histogram_max_adc"]) - histogramMaxADC=config["analysis"]["histogram_max_adc"].as(); + m_HistogramMaxADC = config["analysis"]["histogram_max_adc"].as(); if(config["analysis"]["noise_search_max_adc"]) - noise_search_max_adc=config["analysis"]["noise_search_max_adc"].as(); + m_NoiseSearchMaxADC = config["analysis"]["noise_search_max_adc"].as(); } - vector inputFiles=config["input"]["data_files"].as>(); - if (inputFiles.empty()) { + m_InputFiles = config["input"]["data_files"].as>(); + + if (m_InputFiles.empty()) + { cerr << "Error: No input files provided." << endl; return 1; } - string calibrationFile=config["input"]["calibration_file"].as(); - string stripMapFileStr=config["input"]["strip_map"].as(); - string outfile=config["output"]["prefix"].as(); + + m_CalibrationFile = config["input"]["calibration_file"].as().c_str(); + m_StripMapFile = config["input"]["strip_map"].as().c_str(); + m_OutputPrefix = config["output"]["prefix"].as().c_str(); - EnergyCalHelper helperCal; - TACCalHelper helperCal_TAC; + //EnergyCalHelper m_EnergyCal; + TACCalHelper m_EnergyCal_TAC; - if(!helperCal.Load(calibrationFile)) + if(!m_EnergyCal.Load(m_CalibrationFile.Data())) { - cout<<"Failed to load calibration file: "<(); @@ -430,11 +543,8 @@ int main(int argc,char** argv) cmd_fallback_threshold_keV = atof(optarg); break; - //case 'd': - //cmd_data_file = optarg; - //break; case 'd': - inputFiles.push_back(optarg); + m_InputFiles.push_back(optarg); break; case 'e': @@ -484,191 +594,363 @@ int main(int argc,char** argv) if(cmd_calibration_file != "") { - calibrationFile = cmd_calibration_file; + m_CalibrationFile = cmd_calibration_file; } if(cmd_strip_map != "") { - stripMapFileStr = cmd_strip_map; + m_StripMapFile = cmd_strip_map; } if(cmd_output_prefix != "") { - outfile = cmd_output_prefix; + m_OutputPrefix = cmd_output_prefix; } // Determine output directory from first data file - fs::path dataPath(inputFiles.back()); + fs::path dataPath(m_InputFiles.back()); fs::path outputDir = dataPath.parent_path(); - // Resolve outfile path - fs::path outPath(outfile); + // Resolve m_OutputPrefix path + fs::path outPath(m_OutputPrefix.Data()); if (!outPath.is_absolute()) { - outfile = (outputDir / outPath).string(); + m_OutputPrefix = (outputDir / outPath).string(); } - //outfile = fs::weakly_canonical(outfile).string(); // Debug print - cout << "Resolved output prefix: " << outfile << endl; + cout << "Resolved output prefix: " << m_OutputPrefix << endl; - MString stripMapFile = stripMapFileStr.c_str(); - - /* ------------------------------------------------------------- */ - /* Apply command line overrides if provided */ - /* ------------------------------------------------------------- */ - - if(cmd_min_entries >= 0) - minEntries = cmd_min_entries; - - if(cmd_noise_search_max_adc > 0) - noise_search_max_adc = cmd_noise_search_max_adc; - - if(cmd_fallback_threshold_keV > 0) - fallbackThreshold = cmd_fallback_threshold_keV; - - - MSupervisor* S=MSupervisor::GetSupervisor(); - - map histograms; - map histograms_TAC; - - map>> timingCounts; /* ------------------------------------------------------------- */ /* Save configuration to log file */ /* ------------------------------------------------------------- */ cout< histograms; + map> hist_counts; + map> dt0_counts; + map> dt1_counts; + map adc_to_keV_scale; + map histograms_TAC; + + map>> timingCounts; /* ------------------------------------------------------------- */ /* Build ADC histograms from data */ /* ------------------------------------------------------------- */ + + MModuleLoaderMeasurementsHDF* Loader = new MModuleLoaderMeasurementsHDF(); - for(string inputFile:inputFiles) - { - - MModuleLoaderMeasurementsHDF* Loader=new MModuleLoaderMeasurementsHDF(); - Loader->SetFileName(inputFile.c_str()); - Loader->SetFileNameStripMap(stripMapFile); + Loader->SetFileName(m_InputFiles[0].c_str()); - S->SetModule(Loader,0); + cout << "Loading file: " << m_InputFiles[0] << endl; + cout << "Number of input files: " << m_InputFiles.size() << endl; + + Loader->SetFileNameStripMap(m_StripMapFile.Data()); - if(!Loader->Initialize()) return -1; + MSupervisor* S = MSupervisor::GetSupervisor(); + S->SetModule(Loader, 0); + if (!Loader->Initialize()) + { + cerr << "Failed to initialize loader!" << endl; + return false; + } + + + MReadOutAssembly* Event=new MReadOutAssembly(); long event_counter = 0; - while(Loader->IsFinished()==false) + + while (Loader->IsFinished() == false) { - if(max_events > 0 && event_counter >= max_events) + + + + if (max_events > 0 && event_counter >= max_events) break; - Event->Clear(); - if(Loader->IsReady()) + Event->Clear(); + + if (Loader->IsReady()) { Loader->AnalyzeEvent(Event); - event_counter++; - - for(unsigned int sh=0;shGetNStripHits();sh++) + event_counter++; + + // ------------------------------------------------------------- + // GLOBAL progress bar (event processing) + // ------------------------------------------------------------- + if (event_counter % 10000 == 0) { - MStripHit* SH=Event->GetStripHit(sh); + int barWidth = 40; - if(!PassHitSelection(SH)) continue; + // If we know total → real progress + if (max_events > 0) + { + float progress = (float)event_counter / max_events; - double adc=SH->GetADCUnits(); - - int det=SH->GetDetectorID(); - int strip=SH->GetStripID(); - char side=SH->IsLowVoltageStrip()?'l':'h'; - - StripKey key{det,side,strip}; - - /* ------------------------------------------------------------- */ - /* Fast threshold data accumulation (dt0 vs dt1) */ - /* ------------------------------------------------------------- */ + cout << "\r["; + int pos = barWidth * progress; + + for (int i = 0; i < barWidth; ++i) + { + if (i < pos) cout << "="; + else if (i == pos) cout << ">"; + else cout << " "; + } - int adc_int = (int)adc; + cout << "] " << int(progress * 100.0) << "%"; + } + else + { + // ------------------------------------------------------------- + // UNKNOWN TOTAL → animated progress bar + // ------------------------------------------------------------- + int pos = (event_counter / 10000) % barWidth; - // Define timing type - int timing_type = SH->HasFastTiming() ? 1 : 0; + cout << "\r["; - // Accumulate counts - auto& entry = timingCounts[key][adc_int]; + for (int i = 0; i < barWidth; ++i) + { + if (i == pos) cout << ">"; + else cout << " "; + } - if(timing_type == 0) - entry.first++; - else - entry.second++; - - /* ------------------------------------------------------------- */ - /* TAC (fast shaper) histogram filling */ - /* ------------------------------------------------------------- */ + cout << "] events: " << event_counter; + } + + cout << flush; + } + + + + int NStrips = Event->GetNStripHits(); + + for (int i = 0; i < NStrips; ++i) + { + MStripHit* SH = Event->GetStripHit(i); + + if (!PassHitSelection(SH)) continue; - if(SH->HasFastTiming() && SH->GetTAC() > 0) + double adc = SH->GetADCUnits(); + int det = SH->GetDetectorID(); + int strip = SH->GetStripID(); + char side = SH->IsLowVoltageStrip() ? 'l' : 'h'; + + StripKey key{det, side, strip}; + + // ------------------------------------------------------------- + // Precompute ADC→keV scale (once per strip) + // ------------------------------------------------------------- + if (adc_to_keV_scale.find(key) == adc_to_keV_scale.end()) { - double tac = SH->GetTAC(); + // approximate slope using ADC=0 → ADC=1000 + double e0 = m_EnergyCal.ADCToEnergy(det, side, strip, 0); + double e1 = m_EnergyCal.ADCToEnergy(det, side, strip, 1000); - //StripKey key{det,side,strip}; + adc_to_keV_scale[key] = (e1 - e0) / 1000.0; + } + - if(histograms_TAC[key]==nullptr) - { - string name="h_TAC_"+to_string(det)+"_"+side+"_"+to_string(strip); - histograms_TAC[key]=new TH1D( - name.c_str(), - name.c_str(), - histogramBins,0,histogramMaxADC); + //it_hist->second->Fill(adc); + int bin = (int)(adc / m_HistogramMaxADC * m_HistogramBins); - histograms_TAC[key]->GetXaxis()->SetTitle("TAC ADC"); - histograms_TAC[key]->GetYaxis()->SetTitle("Counts"); + if (bin >= 0 && bin < m_HistogramBins) + { + if (hist_counts.find(key) == hist_counts.end()) + { + hist_counts[key] = vector(m_HistogramBins, 0); } - histograms_TAC[key]->Fill(tac); + hist_counts[key][bin]++; } - - - - static int debugCounter = 0; + + // ------------------------------------------------------------- + // FAST timing accumulation (dt0 vs dt1) + // ------------------------------------------------------------- - if(histograms[key]==nullptr) - { - string name="h_"+to_string(det)+"_"+side+"_"+to_string(strip); + int adc_bin = static_cast(adc); + + double tac = SH->GetTAC(); - histograms[key]=new TH1D(name.c_str(),name.c_str(), - histogramBins,0,histogramMaxADC); + // --- dt0 vs dt1 separation --- + bool is_dt1 = (tac > 8000); // initial threshold - histograms[key]->GetXaxis()->SetTitle("ADC"); - histograms[key]->GetYaxis()->SetTitle("Counts"); + + // Initialize bin if needed + if (timingCounts[key].find(adc_bin) == timingCounts[key].end()) + { + timingCounts[key][adc_bin] = {0, 0}; } - histograms[key]->Fill(adc); + + double energy = adc * adc_to_keV_scale[key]; + + + int ebin = (int)(energy / (adc_to_keV_scale[key] * m_HistogramMaxADC) * m_HistogramBins); + + if (ebin >= 0 && ebin < m_HistogramBins) + { + if (is_dt1) + { + timingCounts[key][adc_bin].second++; + + auto& counts = dt1_counts[key]; + if (counts.empty()) counts.resize(m_HistogramBins, 0); + counts[ebin]++; + } + else + { + timingCounts[key][adc_bin].first++; + + auto& counts = dt0_counts[key]; + if (counts.empty()) counts.resize(m_HistogramBins, 0); + counts[ebin]++; + } + } } } } + cout << endl; + + // Save into class + + + for (auto& kv : hist_counts) + { + StripKey key = kv.first; + vector& counts = kv.second; + + string name = "h_" + to_string(key.det) + "_" + key.side + "_" + to_string(key.strip); + + TH1D* h = new TH1D( + name.c_str(), + name.c_str(), + m_HistogramBins, + 0, + m_HistogramMaxADC + ); + + for (int i = 0; i < m_HistogramBins; ++i) + { + h->SetBinContent(i+1, counts[i]); + } + + m_ADCHistograms[key] = h; + } + + // ------------------------------------------------------------- + // Rebuild dt0 histograms + // ------------------------------------------------------------- + for (auto& kv : dt0_counts) + { + StripKey key = kv.first; + vector& counts = kv.second; + + double maxE = m_EnergyCal.ADCToEnergy(key.det, key.side, key.strip, m_HistogramMaxADC); + + string name = "dt0_" + to_string(key.det) + "_" + key.side + "_" + to_string(key.strip); + + TH1D* h = new TH1D( + name.c_str(), + name.c_str(), + m_HistogramBins, + 0, + maxE + ); + + for (int i = 0; i < m_HistogramBins; ++i) + { + h->SetBinContent(i+1, counts[i]); + } - delete Event; + dt0_hists[key] = h; } + // ------------------------------------------------------------- + // Rebuild dt1 histograms + // ------------------------------------------------------------- + for (auto& kv : dt1_counts) + { + StripKey key = kv.first; + vector& counts = kv.second; + + double maxE = m_EnergyCal.ADCToEnergy(key.det, key.side, key.strip, m_HistogramMaxADC); + + string name = "dt1_" + to_string(key.det) + "_" + key.side + "_" + to_string(key.strip); + + TH1D* h = new TH1D( + name.c_str(), + name.c_str(), + m_HistogramBins, + 0, + maxE + ); + + for (int i = 0; i < m_HistogramBins; ++i) + { + h->SetBinContent(i+1, counts[i]); + } + + dt1_hists[key] = h; + } + + + m_TimingCounts = timingCounts; + + + + return true; + +} + + + +void MStripThresholdFinder::FindSlowThresholds() +{ + + m_stripIndex_LV.clear(); + m_stripIndex_HV.clear(); + m_thresholdValues_LV.clear(); + m_thresholdValues_HV.clear(); + + /* ------------------------------------------------------------- */ /* Threshold finding algorithm */ /* ------------------------------------------------------------- */ @@ -677,7 +959,7 @@ int main(int argc,char** argv) map thresholdsADC; // Progress tracking (slow thresholds) - int totalStrips = histograms.size(); + int totalStrips = m_ADCHistograms.size(); int processedStrips = 0; map thresholdLV; @@ -686,22 +968,12 @@ int main(int argc,char** argv) map thresholdLV_ADC; map thresholdHV_ADC; - /* ------------------------------------------------------------- */ - /* Fast shaper (TAC) histograms */ - /* ------------------------------------------------------------- */ - - //map histograms_TAC; /* ------------------------------------------------------------- */ /* TAC threshold storage */ /* ------------------------------------------------------------- */ - - map thresholds_TAC; - map thresholds_TAC_ADC; - //map>> timingCounts; - - /* HV/LV split */ + // HV/LV split map thresholdLV_TAC; map thresholdHV_TAC; @@ -730,40 +1002,34 @@ int main(int argc,char** argv) vector thresholdValues_LV; vector thresholdValues_HV; - vector noisePeakADC_LV; - vector noisePeakADC_HV; - + /* ------------------------------------------------------------- */ /* TAC diagnostic vectors */ /* ------------------------------------------------------------- */ - vector stripIndex_TAC_LV; - vector stripIndex_TAC_HV; - vector thresholdValues_TAC_LV; vector thresholdValues_TAC_HV; vector tacPeakADC_LV; vector tacPeakADC_HV; - vector slowEnergyVec; - vector tacEnergyVec; - for(auto& kv:histograms) - { + for(auto& kv : m_ADCHistograms) + { + StripKey key=kv.first; TH1D* hist=kv.second; - if(hist->GetEntries()GetEntries()Smooth(3); - int maxSearchBin=hist->FindBin(noise_search_max_adc); + int maxSearchBin = hist->FindBin(m_NoiseSearchMaxADC); int startBin=-1; @@ -772,7 +1038,7 @@ int main(int argc,char** argv) if(startBin<0) { - thresholds[key]=fallbackThreshold; + thresholds[key]=m_FallbackThreshold; continue; } @@ -822,74 +1088,78 @@ int main(int argc,char** argv) /* Convert ADC → keV using SLOW calibration */ double thresholdKeV = - helperCal.ADCToEnergy(key.det,key.side,key.strip,thresholdADC); - - /* Store SLOW thresholds */ + m_EnergyCal.ADCToEnergy(key.det,key.side,key.strip,thresholdADC); + + if(key.side == 'l') + { + m_stripIndex_LV.push_back(key.strip); + m_thresholdValues_LV.push_back(thresholdKeV); + } + else if(key.side == 'h') + { + m_stripIndex_HV.push_back(key.strip); + m_thresholdValues_HV.push_back(thresholdKeV); + } + + + + // Store SLOW thresholds thresholds[key] = thresholdKeV; thresholdsADC[key] = thresholdADC; // Progress bar update - processedStrips++; - - int barWidth = 40; - float progress = (float)processedStrips / totalStrips; + //processedStrips++; - cout << "\r["; - int pos = barWidth * progress; + //int barWidth = 40; + //if (totalStrips == 0) totalStrips = 1; + //float progress = (float)processedStrips / totalStrips; + + + + //cout << "\r["; + //int pos = barWidth * progress; - for(int i = 0; i < barWidth; ++i) - { - if(i < pos) cout << "="; - else if(i == pos) cout << ">"; - else cout << " "; - } + //for(int i = 0; i < barWidth; ++i) + //{ + // if(i < pos) cout << "="; + // else if(i == pos) cout << ">"; + // else cout << " "; + //} - cout << "] " << int(progress * 100.0) << "%" << flush; + //cout << "] " << int(progress * 100.0) << "%" << flush; + + } + cout << endl; - /* ------------------------------------------------------------- */ - /* Store values for diagnostic plots */ - /* ------------------------------------------------------------- */ + // Save results into class members + m_SlowThresholds = thresholds; + m_SlowThresholdsADC = thresholdsADC; + + + cout << "LV points: " << m_stripIndex_LV.size() << endl; + cout << "HV points: " << m_stripIndex_HV.size() << endl; +} - stripIndex.push_back(key.strip); - thresholdValues.push_back(thresholdKeV); - noisePeakADC.push_back(hist->GetBinCenter(peakBin)); - - /* -------------------------------------------------------------- */ - /* Store HV and LV diagnostics separately so both appear in plots */ - /* -------------------------------------------------------------- */ - if(key.side=='l') - { - stripIndex_LV.push_back(key.strip); - thresholdValues_LV.push_back(thresholdKeV); - noisePeakADC_LV.push_back(hist->GetBinCenter(peakBin)); - } - else - { - stripIndex_HV.push_back(key.strip); - thresholdValues_HV.push_back(thresholdKeV); - noisePeakADC_HV.push_back(hist->GetBinCenter(peakBin)); - } - if(key.side=='l') - { - thresholdLV[key.strip]=thresholdKeV; - thresholdLV_ADC[key.strip]=thresholdADC; - } - else - { - thresholdHV[key.strip]=thresholdKeV; - thresholdHV_ADC[key.strip]=thresholdADC; - } - +void MStripThresholdFinder::FindFastThresholds() +{ + + if (m_TimingCounts.empty()) + { + cout << "Warning: No timing data available for fast threshold calculation." << endl; } - cout << endl; + + //map m_FastThresholds; + map thresholds_TAC_ADC; + //map m_FastThresholdsADC; + /* ------------------------------------------------------------- */ /* Fast threshold finder (dt0 vs dt1 crossover) */ /* ------------------------------------------------------------- */ - for(auto& kv : timingCounts) + for(auto& kv : m_TimingCounts) { StripKey key = kv.first; auto& adcMap = kv.second; @@ -904,17 +1174,18 @@ int main(int argc,char** argv) for(auto& a : adcMap) totalCounts += a.second.first + a.second.second; - //if(totalCounts < minEntries) + //if(totalCounts < m_MinEntries) - // Proper minEntries guard - if(totalCounts < minEntries) + // Proper m_MinEntries guard + if(totalCounts < m_MinEntries) { - thresholds_TAC[key] = fallbackThreshold; + m_FastThresholds[key] = m_FallbackThreshold; continue; } // Find first nonzero ADC int first_nonzero = -1; + if (kv.second.empty()) continue; for(auto& a : adcMap) { @@ -927,12 +1198,10 @@ int main(int argc,char** argv) if(first_nonzero < 0) { - thresholds_TAC[key] = fallbackThreshold; + m_FastThresholds[key] = m_FallbackThreshold; continue; } - - //cout << "FAST threshold RMS (keV): " << hFastThreshDist->GetRMS() << endl; int bestADC = -1; int minDiff = 1e9; @@ -992,7 +1261,7 @@ int main(int argc,char** argv) if(bestADC < 0) { - thresholds_TAC[key] = fallbackThreshold; + m_FastThresholds[key] = m_FallbackThreshold; continue; } @@ -1009,51 +1278,34 @@ int main(int argc,char** argv) // fallback to old method fast_thresh_adc = bestADC + nbins/2; } - - /*cout << "[FAST DEBUG] strip " << key.strip - << " bestADC=" << bestADC - << " crossoverADC=" << crossoverADC - << endl;*/ - - thresholds_TAC_ADC[key] = fast_thresh_adc; - - double fast_thresh_keV = - helperCal.ADCToEnergy(key.det,key.side,key.strip,fast_thresh_adc); - thresholds_TAC[key] = fast_thresh_keV; - /* Store diagnostics */ + m_FastThresholdsADC[key] = fast_thresh_adc; - if(key.side=='l') - { - stripIndex_TAC_LV.push_back(key.strip); - thresholdValues_TAC_LV.push_back(fast_thresh_keV); - } - else - { - stripIndex_TAC_HV.push_back(key.strip); - thresholdValues_TAC_HV.push_back(fast_thresh_keV); - } + double fast_thresh_keV = + m_EnergyCal.ADCToEnergy(key.det,key.side,key.strip,fast_thresh_adc); - /* Debug print */ + m_FastThresholds[key] = fast_thresh_keV; - /*cout << "[FAST] DET " << key.det - << " " << key.side - << " strip " << key.strip - << " ADC: " << fast_thresh_adc - << " keV: " << fast_thresh_keV - << endl;*/ - } +} + +void MStripThresholdFinder::WriteCSV() const +{ + + if (m_SlowThresholds.empty() && m_FastThresholds.empty()) + { + cout << "Warning: No thresholds available to write to CSV." << endl; + } /* ------------------------------------------------------------- */ /* Write CSV threshold tables (HV and LV) */ /* ------------------------------------------------------------- */ - ofstream csv_HV(outfile + "_Slow_HV_thresholds.csv"); - ofstream csv_LV(outfile + "_Slow_LV_thresholds.csv"); + ofstream csv_HV(m_OutputPrefix + "_Slow_HV_thresholds.csv"); + ofstream csv_LV(m_OutputPrefix + "_Slow_LV_thresholds.csv"); /* CSV headers */ @@ -1062,13 +1314,13 @@ int main(int argc,char** argv) /* Write rows */ - for(const auto& kv : thresholds) + for(const auto& kv : m_SlowThresholds) { char side = kv.first.side; int strip = kv.first.strip; double thr_keV = kv.second; - double thr_adc = thresholdsADC[kv.first]; + double thr_adc = m_SlowThresholdsADC.at(kv.first); if(side == 'h') @@ -1088,46 +1340,7 @@ int main(int argc,char** argv) } - /* ------------------------------------------------------------- */ - /* Pixel threshold maps */ - /* ------------------------------------------------------------- */ - - TH2D pixelThresholdMap( - "PixelThresholdMap", - "Pixel Threshold Map (keV);LV Strip;HV Strip", - 64,0,64, - 64,64,0); - - TH2D pixelThresholdADCMap( - "PixelThresholdADCMap", - "Pixel Threshold Map (ADC);LV Strip;HV Strip", - 64,0,64, - 64,64,0); - - /* force ROOT to display full strip axes */ - - pixelThresholdMap.GetXaxis()->SetNdivisions(64,false); - pixelThresholdMap.GetYaxis()->SetNdivisions(64,false); - - pixelThresholdADCMap.GetXaxis()->SetNdivisions(64,false); - pixelThresholdADCMap.GetYaxis()->SetNdivisions(64,false); - - - for(int lv=0;lv<64;lv++) - { - if(thresholdLV.find(lv)==thresholdLV.end()) continue; - - for(int hv=0;hv<64;hv++) - { - if(thresholdHV.find(hv)==thresholdHV.end()) continue; - - double pixelThr=max(thresholdLV[lv],thresholdHV[hv]); - pixelThresholdMap.Fill(lv,hv,pixelThr); - double pixelADC=max(thresholdLV_ADC[lv],thresholdHV_ADC[hv]); - pixelThresholdADCMap.Fill(lv,hv,pixelADC); - } - } /* ------------------------------------------------------------- */ @@ -1135,20 +1348,20 @@ int main(int argc,char** argv) /* ------------------------------------------------------------- */ // --- FAST CSV --- - ofstream csv_TAC_HV(outfile + "_Fast_HV_thresholds.csv"); - ofstream csv_TAC_LV(outfile + "_Fast_LV_thresholds.csv"); + ofstream csv_TAC_HV(m_OutputPrefix + "_Fast_HV_thresholds.csv"); + ofstream csv_TAC_LV(m_OutputPrefix + "_Fast_LV_thresholds.csv"); csv_TAC_HV << "detector_side,strip,threshold_adc,threshold_keV\n"; csv_TAC_LV << "detector_side,strip,threshold_adc,threshold_keV\n"; - cout << "Writing FAST CSV entries: " << thresholds_TAC.size() << endl; + cout << "Writing FAST CSV entries: " << m_FastThresholds.size() << endl; - for(const auto& kv : thresholds_TAC) + for(const auto& kv : m_FastThresholds) { char side = kv.first.side; int strip = kv.first.strip; - double thr_adc = thresholds_TAC_ADC[kv.first]; + double thr_adc = m_FastThresholdsADC.at(kv.first); double thr_keV = kv.second; if(side == 'h') @@ -1160,25 +1373,31 @@ int main(int argc,char** argv) csv_TAC_LV << "l," << strip << "," << thr_adc << "," << thr_keV << "\n"; } } - + + + csv_HV.close(); + csv_LV.close(); csv_TAC_HV.close(); csv_TAC_LV.close(); +} + + +void MStripThresholdFinder::WriteDiagnostics() const +{ - + if (m_SlowThresholds.empty() && m_FastThresholds.empty()) + { + cout << "Warning: No thresholds available for diagnostics." << endl; + } /* ------------------------------------------------------------- */ /* ROOT output */ /* ------------------------------------------------------------- */ - //cout << "DEBUG: outfile being used = " << outfile << endl; - TFile f((outfile+"_diagnostics.root").c_str(),"RECREATE"); + //cout << "DEBUG: m_OutputPrefix being used = " << m_OutputPrefix << endl; + TFile f((m_OutputPrefix + "_diagnostics.root").Data(), "RECREATE"); - /* ------------------------------------------------------------- */ - /* Slow and Fast Threshold diagnostic plots */ - /* ------------------------------------------------------------- */ - - // ------------------------------------------------------------- // SLOW threshold per strip (HV vs LV scatter) // ------------------------------------------------------------- @@ -1201,7 +1420,7 @@ int main(int argc,char** argv) gSlowThresh_HV->SetMarkerColor(kRed); // Fill graphs - for(const auto& kv : thresholds) + for(const auto& kv : m_SlowThresholds) { int strip = kv.first.strip; char side = kv.first.side; @@ -1226,7 +1445,7 @@ int main(int argc,char** argv) // Optional: auto-scale Y axis double ymin = 1e9, ymax = -1e9; - for(const auto& kv : thresholds) + for(const auto& kv : m_SlowThresholds) { if(kv.first.strip == 64) continue; double v = kv.second; @@ -1282,7 +1501,7 @@ int main(int argc,char** argv) hSlowDist_HV->SetLineWidth(2); // Fill histograms - for(const auto& kv : thresholds) + for(const auto& kv : m_SlowThresholds) { int strip = kv.first.strip; char side = kv.first.side; @@ -1321,8 +1540,7 @@ int main(int argc,char** argv) hSlowDist_LV->Write(); hSlowDist_HV->Write(); - - + // ------------------------------------------------------------- // FAST threshold per strip histogram @@ -1333,7 +1551,7 @@ int main(int argc,char** argv) 65, 0, 65 ); - for(const auto& kv : thresholds_TAC) + for(const auto& kv : m_FastThresholds) { int strip = kv.first.strip; double thr_keV = kv.second; @@ -1344,145 +1562,171 @@ int main(int argc,char** argv) } hFastThreshPerStrip->Write(); - - - + - /* ------------------------------------------------------------- */ - /* Threshold vs strip for HV and LV sides */ - /* ------------------------------------------------------------- */ - - TGraph Threshold_vs_Strip_LV( - stripIndex_LV.size(), - stripIndex_LV.data(), - thresholdValues_LV.data()); - - Threshold_vs_Strip_LV.SetName("Threshold_vs_Strip_LV"); - Threshold_vs_Strip_LV.SetTitle("Threshold vs Strip;Strip;Threshold (keV)"); - - Threshold_vs_Strip_LV.SetMarkerStyle(20); - Threshold_vs_Strip_LV.SetMarkerColor(kBlue); - Threshold_vs_Strip_LV.SetMarkerSize(1); + // ------------------------------------------------------------- + // Threshold vs strip for HV and LV sides + // ------------------------------------------------------------- - Threshold_vs_Strip_LV.Write(); - - - TGraph Threshold_vs_Strip_HV( - stripIndex_HV.size(), - stripIndex_HV.data(), - thresholdValues_HV.data()); + if (!m_stripIndex_LV.empty()) + { + TGraph Threshold_vs_Strip_LV( + m_stripIndex_LV.size(), + m_stripIndex_LV.data(), + m_thresholdValues_LV.data()); - Threshold_vs_Strip_HV.SetName("Threshold_vs_Strip_HV"); + Threshold_vs_Strip_LV.SetName("Threshold_vs_Strip_LV"); + Threshold_vs_Strip_LV.SetTitle("Threshold vs Strip (LV)"); + Threshold_vs_Strip_LV.SetMarkerStyle(20); + Threshold_vs_Strip_LV.SetMarkerColor(kBlue); - Threshold_vs_Strip_HV.SetMarkerStyle(20); - Threshold_vs_Strip_HV.SetMarkerColor(kRed); - Threshold_vs_Strip_HV.SetMarkerSize(1); + Threshold_vs_Strip_LV.Write(); + } - Threshold_vs_Strip_HV.Write(); + if (!m_stripIndex_HV.empty()) + { + TGraph Threshold_vs_Strip_HV( + m_stripIndex_HV.size(), + m_stripIndex_HV.data(), + m_thresholdValues_HV.data()); + Threshold_vs_Strip_HV.SetName("Threshold_vs_Strip_HV"); + Threshold_vs_Strip_HV.SetTitle("Threshold vs Strip (HV)"); + Threshold_vs_Strip_HV.SetMarkerStyle(20); + Threshold_vs_Strip_HV.SetMarkerColor(kRed); + Threshold_vs_Strip_HV.Write(); + } - /* ------------------------------------------------------------- */ - /* Noise peak ADC vs strip for HV and LV sides */ - /* ------------------------------------------------------------- */ - TGraph NoisePeakADC_vs_Strip_LV( - stripIndex_LV.size(), - stripIndex_LV.data(), - noisePeakADC_LV.data()); + + for (auto& kv : dt0_hists) + { + const StripKey& key = kv.first; + + if (dt1_hists.find(key) == dt1_hists.end()) continue; + + TH1D* dt0 = kv.second; + TH1D* dt1 = dt1_hists.at(key); + + int det = key.det; + char side = key.side; + int strip = key.strip; + + // ------------------------------- + // Create canvas + // ------------------------------- + string cname = "cFast_" + to_string(det) + "_" + side + "_" + to_string(strip); + TCanvas* c = new TCanvas(cname.c_str(), cname.c_str(), 800, 600); + + // ------------------------------- + // Styling + // ------------------------------- + dt0->SetLineColor(kBlack); + dt1->SetLineColor(kBlue); + + dt0->SetTitle(Form( + "Fast Shaper Threshold (det=%d %c strip=%d);Energy (keV);Counts", + det, side, strip + )); + dt0->GetXaxis()->SetRangeUser(0, 100); + + - NoisePeakADC_vs_Strip_LV.SetName("NoisePeakADC_vs_Strip_LV"); - NoisePeakADC_vs_Strip_LV.SetTitle("Noise Peak ADC vs Strip;Strip;Noise Peak (ADC)"); + // ------------------------------- + // Draw both + // ------------------------------- + dt0->Draw("HIST"); + dt1->Draw("HIST SAME"); + + gPad->Update(); - NoisePeakADC_vs_Strip_LV.SetMarkerStyle(20); - NoisePeakADC_vs_Strip_LV.SetMarkerColor(kBlue); - NoisePeakADC_vs_Strip_LV.SetMarkerSize(1); - NoisePeakADC_vs_Strip_LV.Write(); + TLegend* leg = new TLegend(0.65,0.7,0.88,0.88); + leg->AddEntry(dt0, "dt0", "l"); + leg->AddEntry(dt1, "dt1", "l"); + // ------------------------------- + // Threshold line (NOW IN keV) + // ------------------------------- + if (m_FastThresholds.find(key) != m_FastThresholds.end()) + { - TGraph NoisePeakADC_vs_Strip_HV( - stripIndex_HV.size(), - stripIndex_HV.data(), - noisePeakADC_HV.data()); + double thr_keV = m_FastThresholds.at(key); - NoisePeakADC_vs_Strip_HV.SetName("NoisePeakADC_vs_Strip_HV"); + double ymin = gPad->GetUymin(); + double ymax = gPad->GetUymax(); - NoisePeakADC_vs_Strip_HV.SetMarkerStyle(20); - NoisePeakADC_vs_Strip_HV.SetMarkerColor(kRed); - NoisePeakADC_vs_Strip_HV.SetMarkerSize(1); + TLine* line = new TLine(thr_keV, ymin, thr_keV, ymax); + line->SetLineColor(kRed); + line->SetLineWidth(2); + line->Draw("SAME"); - NoisePeakADC_vs_Strip_HV.Write(); + // Legend entry (fix from earlier) + leg->AddEntry(line, "Fast Threshold", "l"); + } - + //leg->AddEntry((TObject*)0, "Fast Threshold", ""); // label only + leg->Draw(); + // ------------------------------- + // Save + // ------------------------------- + c->Write(); + } - /* ------------------------------------------------------------- */ - /* TAC Threshold vs Strip */ - /* ------------------------------------------------------------- */ - - TGraph Threshold_TAC_vs_Strip_LV( - stripIndex_TAC_LV.size(), - stripIndex_TAC_LV.data(), - thresholdValues_TAC_LV.data()); - - Threshold_TAC_vs_Strip_LV.SetName("Threshold_TAC_vs_Strip_LV"); - Threshold_TAC_vs_Strip_LV.SetTitle("TAC Threshold vs Strip;Strip;Threshold (keV)"); - Threshold_TAC_vs_Strip_LV.SetMarkerStyle(20); - Threshold_TAC_vs_Strip_LV.SetMarkerColor(kBlue); - Threshold_TAC_vs_Strip_LV.SetMarkerSize(1); - - Threshold_TAC_vs_Strip_LV.Write(); - - - TGraph Threshold_TAC_vs_Strip_HV( - stripIndex_TAC_HV.size(), - stripIndex_TAC_HV.data(), - thresholdValues_TAC_HV.data()); - - Threshold_TAC_vs_Strip_HV.SetName("Threshold_TAC_vs_Strip_HV"); - Threshold_TAC_vs_Strip_HV.SetMarkerStyle(20); - Threshold_TAC_vs_Strip_HV.SetMarkerColor(kRed); - Threshold_TAC_vs_Strip_HV.SetMarkerSize(1); - Threshold_TAC_vs_Strip_HV.Write(); - - - + /* ------------------------------------------------------------- */ /* Write ADC spectra and create energy spectra with thresholds */ /* ------------------------------------------------------------- */ - for(auto& kv:histograms) + for(auto& kv:m_ADCHistograms) { StripKey key=kv.first; TH1D* adcHist=kv.second; adcHist->Write(); - - string name="Energy_"+to_string(key.det)+"_"+key.side+"_"+to_string(key.strip); + + int det = key.det; + char side = key.side; + int strip = key.strip; + + + //string name="Slow Threshold Det "+to_string(key.det)+", Side"+key.side+", Strip"+to_string(key.strip); + string name = "Energy_" + to_string(det) + "_" + side + "_" + to_string(strip); + + double maxE = m_EnergyCal.ADCToEnergy(det, side, strip, m_HistogramMaxADC); + TH1D* energyHist=new TH1D( name.c_str(), name.c_str(), adcHist->GetNbinsX(), 0, - helperCal.ADCToEnergy(key.det,key.side,key.strip,histogramMaxADC) + maxE ); - + + energyHist->SetTitle(Form( + "Slow Threshold (det=%d side=%c strip=%d);Energy (keV);Counts", + det, side, strip + )); + + energyHist->GetXaxis()->SetRangeUser(0, 100); + for(int b=1;b<=adcHist->GetNbinsX();b++) { double adc=adcHist->GetBinCenter(b); - double energy=helperCal.ADCToEnergy(key.det,key.side,key.strip,adc); + double energy=m_EnergyCal.ADCToEnergy(key.det,key.side,key.strip,adc); double counts=adcHist->GetBinContent(b); int ebin=energyHist->FindBin(energy); energyHist->AddBinContent(ebin,counts); } - double thr=thresholds[key]; + double thr = m_SlowThresholds.at(key); TLine* line=new TLine(thr,0,thr,energyHist->GetMaximum()); line->SetLineColor(kRed); @@ -1517,7 +1761,7 @@ int main(int argc,char** argv) hFastDist_HV->SetLineWidth(2); // Fill histograms - for(const auto& kv : thresholds_TAC) + for(const auto& kv : m_FastThresholds) { int strip = kv.first.strip; char side = kv.first.side; @@ -1565,7 +1809,7 @@ int main(int argc,char** argv) 200, 0, 2000 ); - for(const auto& kv : thresholds_TAC_ADC) + for(const auto& kv : m_FastThresholdsADC) { FastThresholdDistribution_ADC->Fill(kv.second); } @@ -1595,7 +1839,7 @@ int main(int argc,char** argv) gFastThresh_HV->SetMarkerColor(kRed); // Fill graphs - for(const auto& kv : thresholds_TAC) + for(const auto& kv : m_FastThresholds) { int strip = kv.first.strip; char side = kv.first.side; @@ -1652,111 +1896,6 @@ int main(int argc,char** argv) gFastThresh_HV->Write(); - /* ------------------------------------------------------------- */ - /* Build Slow vs TAC comparison */ - /* ------------------------------------------------------------- */ - - for(const auto& kv : thresholds) - { - StripKey key = kv.first; - - if(thresholds_TAC.find(key) == thresholds_TAC.end()) - continue; - - double slowE = kv.second; - double tacE = thresholds_TAC[key]; - - if(slowE > 0 && tacE > 0) - { - slowEnergyVec.push_back(slowE); - tacEnergyVec.push_back(tacE); - } - } - - TGraph Slow_vs_TAC( - slowEnergyVec.size(), - slowEnergyVec.data(), - tacEnergyVec.data()); - - Slow_vs_TAC.SetName("Slow_vs_TAC"); - Slow_vs_TAC.SetTitle("Slow vs TAC Energy;Slow Energy (keV);TAC Energy (keV)"); - Slow_vs_TAC.SetMarkerStyle(20); - Slow_vs_TAC.SetMarkerSize(1); - - Slow_vs_TAC.Write(); - - - /* ------------------------------------------------------------- */ - /* dt0 vs dt1 diagnostic histogram */ - /* ------------------------------------------------------------- */ - - - - for(auto& kv : timingCounts) - { - StripKey key = kv.first; - auto& adcMap = kv.second; - - if(key.strip == 64) - { - continue; - } - - string name0 = "dt0_" + to_string(key.det) + "_" + key.side + "_" + to_string(key.strip); - string name1 = "dt1_" + to_string(key.det) + "_" + key.side + "_" + to_string(key.strip); - - TH1D* h0 = new TH1D(name0.c_str(), name0.c_str(), 500, 0, 100); - TH1D* h1 = new TH1D(name1.c_str(), name1.c_str(), 500, 0, 100); - - for(auto& a : adcMap) - { - int adc = a.first; - int n0 = a.second.first; - int n1 = a.second.second; - - //double energy = calib[key.det][key.side][key.strip].Evaluate(adc); - double energy = helperCal.ADCToEnergy(key.det, key.side, key.strip, adc); - - //h0->Fill(energy, n0); - double e0 = helperCal.ADCToEnergy(key.det, key.side, key.strip, adc); - double e1 = helperCal.ADCToEnergy(key.det, key.side, key.strip, adc + 1); - - // distribute counts across the interval - int nSub = 5; // small subdivision - for(int i = 0; i < nSub; ++i) - { - double e = e0 + (e1 - e0)*(i + 0.5)/nSub; - h0->Fill(e, n0 / (double)nSub); - h1->Fill(e, n1 / (double)nSub); - } - - // h1->Fill(energy, n1); - } - - h0->SetLineColor(kRed); - h1->SetLineColor(kBlue); - - // --- FAST threshold line (energy space) --- - if(thresholds_TAC_ADC.find(key) != thresholds_TAC_ADC.end()) - { - double thr_adc = thresholds_TAC_ADC[key]; - //double thr_keV = calib[key.det][key.side][key.strip].Evaluate(thr_adc); - double thr_keV = helperCal.ADCToEnergy(key.det, key.side, key.strip, thr_adc); - - double maxY = max(h0->GetMaximum(), h1->GetMaximum()); - - TLine* line = new TLine(thr_keV, 0, thr_keV, maxY); - line->SetLineColor(kGreen+2); - line->SetLineWidth(2); - line->SetLineStyle(2); // optional (dashed) - - h0->GetListOfFunctions()->Add(line); - } - - h0->Write(); - h1->Write(); - } - // ------------------------------------------------------------- // SLOW threshold pixel map (LV vs HV) // ------------------------------------------------------------- @@ -1771,8 +1910,8 @@ int main(int argc,char** argv) map slowLV; map slowHV; - // Separate thresholds by side - for(const auto& kv : thresholds) + // Separate m_SlowThresholds by side + for(const auto& kv : m_SlowThresholds) { int strip = kv.first.strip; char side = kv.first.side; @@ -1824,69 +1963,83 @@ int main(int argc,char** argv) cSlowPixel->Write(); hSlowPixelMap->Write(); + + + //------------------------------------------------------------ + // New Plots + //------------------------------------------------------------ + // Plot Slow Thresholds per Strip - + TGraph* gSlowLV = new TGraph(); + TGraph* gSlowHV = new TGraph(); - f.Close(); - /* ------------------------------------------------------------- */ - /* Helpful ROOT instructions */ - /* ------------------------------------------------------------- */ + for (const auto& kv : m_SlowThresholds) + { + const StripKey& key = kv.first; + double thr = kv.second; - cout<SetPoint(gSlowLV->GetN(), key.strip, thr); + } + else if (key.side == 'h') + { + gSlowHV->SetPoint(gSlowHV->GetN(), key.strip, thr); + } + } - cout<<"Open file:"<SetTitle("Slow Threshold vs Strip (LV);Strip;Threshold [keV]"); + gSlowHV->SetTitle("Slow Threshold vs Strip (HV);Strip;Threshold [keV]"); - cout<<"Slow threshold diagnostic plots:"<Divide(2,1); - cout<<"Slow threshold distribution:"<Draw()"<cd(1); + gSlowLV->Draw("AP"); - cout<Draw()"<Draw()"<GetXaxis()->SetRangeUser(0,30)"<cd(2); + gSlowHV->Draw("AP"); - cout<<"Pixel threshold heat maps:"<Draw()"<Write(); + // PLot Fast Thresholds per strip - - cout<Draw()"<Draw()"<Draw()"<SetLineColor(kBlue)"<Draw(\"SAME\")"<SetPoint(gFastLV->GetN(), key.strip, thr); + } + else if (key.side == 'h') + { + gFastHV->SetPoint(gFastHV->GetN(), key.strip, thr); + } + } + + TCanvas* cFastStrip = new TCanvas("cFastStrip","Fast Threshold vs Strip",1200,500); + cFastStrip->Divide(2,1); + + cFastStrip->cd(1); + gFastLV->Draw("AP"); + + cFastStrip->cd(2); + gFastHV->Draw("AP"); + + cFastStrip->Write(); + + f.Close(); }