new file: GainMatchSX3Front.C using a 2 step fit for the same method, going to implement a inverse fit now wherein the fronts are fit first and then the sum is fit against the backs.
new file: GainMatchSX3Front.h new file: GainMatchSX3Front1.C
This commit is contained in:
parent
dd2ec66db1
commit
4fc05ea338
4
.vscode/settings.json
vendored
4
.vscode/settings.json
vendored
|
|
@ -111,7 +111,9 @@
|
|||
"RelBack_Fix_new.C": "cpp",
|
||||
"SiRelativeGains_Step1_new.C": "cpp",
|
||||
"charconv": "cpp",
|
||||
"format": "cpp"
|
||||
"format": "cpp",
|
||||
"GainMatchSX3Front.C": "cpp",
|
||||
"GainMatchSX3Front1.C": "cpp"
|
||||
},
|
||||
"github-enterprise.uri": "https://fsunuc.physics.fsu.edu"
|
||||
}
|
||||
226
GainMatchSX3.C
226
GainMatchSX3.C
|
|
@ -12,8 +12,6 @@
|
|||
#include <algorithm>
|
||||
#include <TProfile.h>
|
||||
#include "Armory/ClassSX3.h"
|
||||
#include "TGraphErrors.h"
|
||||
#include "TMultiDimFit.h"
|
||||
|
||||
#include "TVector3.h"
|
||||
|
||||
|
|
@ -28,10 +26,8 @@ int padID = 0;
|
|||
|
||||
SX3 sx3_contr;
|
||||
TCutG *cut;
|
||||
TCutG *cut1;
|
||||
std::map<std::tuple<int, int, int, int>, std::vector<std::tuple<double, double, double>>> dataPoints;
|
||||
|
||||
|
||||
void GainMatchSX3::Begin(TTree * /*tree*/)
|
||||
{
|
||||
TString option = GetOption();
|
||||
|
|
@ -60,21 +56,6 @@ void GainMatchSX3::Begin(TTree * /*tree*/)
|
|||
return;
|
||||
}
|
||||
cut->SetName("sx3cut"); // Ensure the cut has the correct name
|
||||
|
||||
// Load the TCutG object
|
||||
TFile *cutFile1 = TFile::Open("UvD.root");
|
||||
if (!cutFile1 || cutFile1->IsZombie())
|
||||
{
|
||||
std::cerr << "Error: Could not open UvD.root" << std::endl;
|
||||
return;
|
||||
}
|
||||
cut1 = dynamic_cast<TCutG *>(cutFile1->Get("UvD"));
|
||||
if (!cut1)
|
||||
{
|
||||
std::cerr << "Error: Could not find TCutG named 'UvD' in UvD.root" << std::endl;
|
||||
return;
|
||||
}
|
||||
cut1->SetName("UvD");
|
||||
}
|
||||
|
||||
Bool_t GainMatchSX3::Process(Long64_t entry)
|
||||
|
|
@ -104,11 +85,11 @@ Bool_t GainMatchSX3::Process(Long64_t entry)
|
|||
for (int i = 0; i < sx3.multi; i++)
|
||||
{
|
||||
|
||||
for (int j = i + 1; j < sx3.multi; j++)
|
||||
{
|
||||
if (sx3.id[i] == 3)
|
||||
hsx3Coin->Fill(sx3.index[i], sx3.index[j]);
|
||||
}
|
||||
// for (int j = i + 1; j < sx3.multi; j++)
|
||||
// {
|
||||
// if (sx3.id[i] == 3)
|
||||
// hsx3Coin->Fill(sx3.index[i], sx3.index[j]);
|
||||
// }
|
||||
if (sx3.e[i] > 100)
|
||||
{
|
||||
ID.push_back(std::pair<int, int>(sx3.id[i], i));
|
||||
|
|
@ -186,7 +167,7 @@ Bool_t GainMatchSX3::Process(Long64_t entry)
|
|||
|
||||
for (int i = 0; i < sx3.multi; i++)
|
||||
{
|
||||
if (sx3.id[i] == 3 && sx3.e[i] > 100)
|
||||
if (sx3.id[i] == 3)
|
||||
{
|
||||
// Fill the histogram for the front vs back with gain correction
|
||||
hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
|
|
@ -208,7 +189,7 @@ Bool_t GainMatchSX3::Process(Long64_t entry)
|
|||
|
||||
hist2d->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
|
||||
if (cut && cut->IsInside(sx3EUp + sx3EDn, sx3EBk))
|
||||
// if (cut && cut->IsInside(sx3EUp + sx3EDn, sx3EBk))
|
||||
// if (sx3.id[i] < 24 && sx3ChUp < 4 && sx3ChBk < 4 && std::isfinite(sx3EUp) && std::isfinite(sx3EDn) && std::isfinite(sx3EBk))
|
||||
{
|
||||
// Accumulate data for gain matching
|
||||
|
|
@ -224,9 +205,6 @@ Bool_t GainMatchSX3::Process(Long64_t entry)
|
|||
|
||||
void GainMatchSX3::Terminate()
|
||||
{
|
||||
|
||||
// --- Store fit coefficients in memory ---
|
||||
std::map<std::tuple<int, int, int, int>, TVectorD> fitCoefficients;
|
||||
const int MAX_DET = 24;
|
||||
const int MAX_UP = 4;
|
||||
const int MAX_DOWN = 4;
|
||||
|
|
@ -234,124 +212,126 @@ void GainMatchSX3::Terminate()
|
|||
|
||||
double gainArray[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
|
||||
bool gainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
|
||||
double fbgain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
|
||||
bool fbgainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
|
||||
|
||||
std::ofstream outFile("sx3_MultiDimFit_results.txt");
|
||||
// std::map<int, TH2F *> updn2DHistos;
|
||||
std::map<int, double> upCorrFactor;
|
||||
|
||||
// === Gain matching ===
|
||||
|
||||
std::ofstream outFile("sx3_GainMatchback.txt");
|
||||
if (!outFile.is_open())
|
||||
{
|
||||
std::cerr << "Error opening output file!" << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
// === Loop over all (id, bk, up, dn) combinations ===
|
||||
for (const auto &kv : dataPoints) {
|
||||
// Gain fit using up+dn vs bk
|
||||
for (const auto &kv : dataPoints)
|
||||
|
||||
{
|
||||
// kv.first is a tuple of (id, up, bk)
|
||||
// kv.second is a vector of tuples (bkE, upE, dnE)
|
||||
auto [id, bk, u, d] = kv.first;
|
||||
const auto &pts = kv.second;
|
||||
// Check if we have enough points for fitting
|
||||
if (pts.size() < 5)
|
||||
continue;
|
||||
|
||||
if (pts.size() < 20) continue;
|
||||
std::vector<double> bkE, udE;
|
||||
|
||||
std::vector<double> x_bk, x_up, y_fsum;
|
||||
for (const auto &pr : pts) {
|
||||
double eBk, eUp, eDn;
|
||||
for (const auto &pr : pts)
|
||||
{
|
||||
double eUp, eDn, eBk;
|
||||
std::tie(eBk, eUp, eDn) = pr;
|
||||
if (eBk > 0 && eUp > 0 && eDn > 0) {
|
||||
x_bk.push_back(eBk);
|
||||
x_up.push_back(eUp);
|
||||
y_fsum.push_back(eUp + eDn);
|
||||
bkE.push_back(eBk);
|
||||
udE.push_back(eUp + eDn);
|
||||
}
|
||||
|
||||
// Fill the TGraph with bkE and udE
|
||||
TGraph g(bkE.size(), bkE.data(), udE.data());
|
||||
// Fit the graph to a linear function
|
||||
TF1 f("f", "[0]*x", 0, 16000);
|
||||
g.Fit(&f, "QNR");
|
||||
if (TMath::Abs(f.GetParameter(0) - 1) > 1)
|
||||
{
|
||||
continue; // Skip this fit if the slope is too far from 1
|
||||
}
|
||||
gainArray[id][bk][u][d] = f.GetParameter(0);
|
||||
gainValid[id][bk][u][d] = true;
|
||||
}
|
||||
|
||||
// Output results
|
||||
for (int id = 0; id < MAX_DET; ++id)
|
||||
{
|
||||
for (int bk = 0; bk < MAX_BK; ++bk)
|
||||
{
|
||||
for (int u = 0; u < MAX_UP; ++u)
|
||||
{
|
||||
for (int d = 0; d < MAX_DOWN; ++d)
|
||||
{
|
||||
// Check if the gain is valid for this detector, back, up, and down
|
||||
if (gainValid[id][bk][u][d])
|
||||
{
|
||||
outFile << id << " " << bk << " " << u << " " << d << " " << gainArray[id][u][d][bk] << std::endl;
|
||||
if (TMath::Abs(gainArray[id][u][d][bk] - 1) < 0.3)
|
||||
{
|
||||
printf("Gain match Det%d Up%dDn%d Back%d → %.4f \n", id, u, d, bk, gainArray[id][u][d][bk]);
|
||||
}
|
||||
else if (gainArray[id][u][d][bk]!=0)
|
||||
{
|
||||
std::cerr << "Warning: Gain value out of range for Det " << id << " Up " << u << " Dn " << d << " Back " << bk << ": "
|
||||
<< gainArray[id][u][d][bk] << std::endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
int nPoints = y_fsum.size();
|
||||
if (nPoints < 20) continue;
|
||||
|
||||
TMultiDimFit *mdf = new TMultiDimFit(2, TMultiDimFit::kMonomials);
|
||||
mdf->SetMaxPowers(new Int_t[2]{1, 1});
|
||||
mdf->SetMinAngle(10);
|
||||
mdf->SetMinRelativeError(1e-4);
|
||||
|
||||
double *x_row = new double[2];
|
||||
for (int i = 0; i < nPoints; ++i) {
|
||||
x_row[0] = x_bk[i];
|
||||
x_row[1] = x_up[i];
|
||||
mdf->AddRow(x_row, y_fsum[i]);
|
||||
}
|
||||
delete[] x_row;
|
||||
|
||||
mdf->Fit();
|
||||
|
||||
const TVectorD *coeffs = mdf->GetCoefficients();
|
||||
if (!coeffs || coeffs->GetNoElements() == 0 || !TMath::Finite((*coeffs)(0))) {
|
||||
std::cerr << "Fit failed for Det" << id << " B" << bk << " U" << u << " D" << d << std::endl;
|
||||
delete mdf;
|
||||
continue;
|
||||
}
|
||||
|
||||
// Store coefficients in the map and write to file
|
||||
fitCoefficients[kv.first] = *coeffs;
|
||||
|
||||
int nCoeffs = mdf->GetNCoefficients();
|
||||
outFile << id << " " << bk << " " << u << " " << d;
|
||||
printf("Fit for Det%d B%d U%d D%d -> ", id, bk, u, d);
|
||||
for (int i = 0; i < nCoeffs; ++i) {
|
||||
outFile << " " << (*coeffs)(i);
|
||||
printf("p%d: %.4f ", i, (*coeffs)(i));
|
||||
}
|
||||
outFile << std::endl;
|
||||
printf("\n");
|
||||
|
||||
delete mdf;
|
||||
}
|
||||
|
||||
// for (int bk = 0; bk < MAX_BK; ++bk)
|
||||
// {
|
||||
// TString name = Form("hUpDnVsBk_%d", bk);
|
||||
// TString title = Form("Up/Bk vs Dn/Bk for Back %d;Dn/Bk;Up/Bk", bk);
|
||||
// updn2DHistos[bk] = new TH2F(name, title, 400, 0, 1, 400, 0, 1);
|
||||
// }
|
||||
|
||||
outFile.close();
|
||||
std::cout << "Multi-dimensional gain matching complete. Results saved." << std::endl;
|
||||
std::cout << "Gain matching complete." << std::endl;
|
||||
|
||||
// --- Stage 2: Apply corrections and create new histograms ---
|
||||
std::cout << "--- Stage 2: Applying Corrections and Visualizing Results ---" << std::endl;
|
||||
TH2F *hCorrectedFvB = new TH2F("hCorrectedFvB", "Gain Corrected Data;Predicted Front Sum (from fit);Measured Front Sum", 400, 0, 16000, 400, 0, 16000);
|
||||
// === Create histograms ===
|
||||
TH2F *hFVB = new TH2F("hFVB", "Corrected Up+Dn vs Corrected Back;Corrected Back E;Up+Dn E",
|
||||
400, 0, 16000, 400, 0, 16000);
|
||||
TH2F *hAsym = new TH2F("hAsym", "Up vs Dn dvide corrected back;Up/Back E;Dn/Back E",
|
||||
400, 0.0, 1.0, 400, 0.0, 1.0);
|
||||
|
||||
for (const auto &kv : dataPoints) {
|
||||
// Find the coefficients for this segment
|
||||
if (fitCoefficients.find(kv.first) == fitCoefficients.end()) {
|
||||
continue; // Skip if no valid fit was found
|
||||
}
|
||||
const TVectorD& coeffs = fitCoefficients[kv.first];
|
||||
double p0 = coeffs(0);
|
||||
double p1 = coeffs(1);
|
||||
double p2 = coeffs(2);
|
||||
// Fill histograms
|
||||
for (const auto &kv : dataPoints)
|
||||
{
|
||||
auto [id, u, d, bk] = kv.first;
|
||||
if (!gainValid[id][u][d][bk])
|
||||
continue;
|
||||
double gain = gainArray[id][u][d][bk];
|
||||
|
||||
// Loop over the data points for this segment and apply the correction
|
||||
const auto &pts = kv.second;
|
||||
for (const auto &pr : pts) {
|
||||
// Prepare vectors to hold the points for TGraph
|
||||
std::vector<double> xVals;
|
||||
std::vector<double> yVals;
|
||||
|
||||
for (const auto &pr : kv.second)
|
||||
{
|
||||
double eBk, eUp, eDn;
|
||||
std::tie(eBk, eUp, eDn) = pr;
|
||||
|
||||
// Calculate the predicted front sum using the fit parameters
|
||||
double predicted_front_sum = p0 + p1 * eBk + p2 * eUp;
|
||||
|
||||
// The measured front sum is just the raw sum
|
||||
double measured_front_sum = eUp + eDn;
|
||||
double updn = eUp + eDn;
|
||||
if (updn == 0 || eBk == 0)
|
||||
continue;
|
||||
|
||||
// Fill the corrected histogram
|
||||
hCorrectedFvB->Fill(predicted_front_sum, measured_front_sum);
|
||||
double asym = (eUp - eDn) / updn;
|
||||
double correctedBack = eBk * gain;
|
||||
|
||||
hFVB->Fill(correctedBack, updn);
|
||||
hAsym->Fill(eUp / correctedBack, eDn / correctedBack);
|
||||
}
|
||||
}
|
||||
|
||||
// --- Stage 3: Draw the comparison canvases ---
|
||||
gStyle->SetOptStat(1110);
|
||||
TCanvas *c1 = new TCanvas("c1", "Gain Correction Results", 1200, 600);
|
||||
c1->Divide(2, 1);
|
||||
|
||||
c1->cd(1);
|
||||
hSX3FvsB_g->SetTitle("Before Correction (Gated)");
|
||||
hSX3FvsB_g->GetXaxis()->SetTitle("Measured Front Sum (E_Up + E_Dn)");
|
||||
hSX3FvsB_g->GetYaxis()->SetTitle("Measured Back E");
|
||||
hSX3FvsB_g->Draw("colz");
|
||||
|
||||
c1->cd(2);
|
||||
hCorrectedFvB->SetTitle("After Correction");
|
||||
hCorrectedFvB->Draw("colz");
|
||||
// Draw a perfect y=x line for comparison
|
||||
TF1 *diag = new TF1("diag", "x", 0, 16000);
|
||||
diag->SetLineColor(kRed);
|
||||
diag->SetLineWidth(2);
|
||||
diag->Draw("same");
|
||||
}
|
||||
338
GainMatchSX3Front.C
Normal file
338
GainMatchSX3Front.C
Normal file
|
|
@ -0,0 +1,338 @@
|
|||
#define GainMatchSX3Front_cxx
|
||||
|
||||
#include "GainMatchSX3Front.h"
|
||||
#include <TH2.h>
|
||||
#include <TF1.h>
|
||||
#include <TStyle.h>
|
||||
#include <TCanvas.h>
|
||||
#include <TMath.h>
|
||||
#include <TCutG.h>
|
||||
#include <fstream>
|
||||
#include <utility>
|
||||
#include <algorithm>
|
||||
#include <TProfile.h>
|
||||
#include "Armory/ClassSX3.h"
|
||||
#include "TGraphErrors.h"
|
||||
#include "TMultiDimFit.h"
|
||||
|
||||
#include "TVector3.h"
|
||||
|
||||
TH2F *hSX3FvsB;
|
||||
TH2F *hSX3FvsB_g;
|
||||
TH2F *hsx3IndexVE;
|
||||
TH2F *hsx3IndexVE_g;
|
||||
TH2F *hSX3;
|
||||
TH2F *hsx3Coin;
|
||||
|
||||
int padID = 0;
|
||||
|
||||
SX3 sx3_contr;
|
||||
TCutG *cut;
|
||||
TCutG *cut1;
|
||||
std::map<std::tuple<int, int, int, int>, std::vector<std::tuple<double, double, double>>> dataPoints;
|
||||
|
||||
// Gain arrays
|
||||
|
||||
const int MAX_DET = 24;
|
||||
const int MAX_UP = 4;
|
||||
const int MAX_DOWN = 4;
|
||||
const int MAX_BK = 4;
|
||||
double backGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
|
||||
bool backGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
|
||||
double frontGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
|
||||
bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
|
||||
|
||||
void GainMatchSX3Front::Begin(TTree * /*tree*/)
|
||||
{
|
||||
TString option = GetOption();
|
||||
|
||||
hSX3FvsB = new TH2F("hSX3FvsB", "SX3 Front vs Back; Front E; Back E", 400, 0, 16000, 400, 0, 16000);
|
||||
hSX3FvsB_g = new TH2F("hSX3FvsB_g", "SX3 Front vs Back; Front E; Back E", 400, 0, 16000, 400, 0, 16000);
|
||||
hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000);
|
||||
hsx3IndexVE_g = new TH2F("hsx3IndexVE_g", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000);
|
||||
hSX3 = new TH2F("hSX3", "SX3 Front v Back; Fronts; Backs", 8, 0, 8, 4, 0, 4);
|
||||
|
||||
hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24 * 12, 0, 24 * 12, 24 * 12, 0, 24 * 12);
|
||||
|
||||
sx3_contr.ConstructGeo();
|
||||
|
||||
// Load the TCutG object
|
||||
TFile *cutFile = TFile::Open("sx3cut.root");
|
||||
bool cutLoaded = (cut != nullptr);
|
||||
cut = dynamic_cast<TCutG *>(cutFile->Get("sx3cut"));
|
||||
if (!cut)
|
||||
{
|
||||
std::cerr << "Error: Could not find TCutG named 'sx3cut' in sx3cut.root" << std::endl;
|
||||
return;
|
||||
}
|
||||
cut->SetName("sx3cut"); // Ensure the cut has the correct name
|
||||
|
||||
// Load the TCutG object
|
||||
TFile *cutFile1 = TFile::Open("UvD.root");
|
||||
bool cut1Loaded = (cut1 != nullptr);
|
||||
cut1 = dynamic_cast<TCutG *>(cutFile1->Get("UvD"));
|
||||
if (!cut1)
|
||||
{
|
||||
std::cerr << "Error: Could not find TCutG named 'UvD' in UvD.root" << std::endl;
|
||||
return;
|
||||
}
|
||||
cut1->SetName("UvD");
|
||||
std::string filename = "sx3_GainMatchback.txt";
|
||||
|
||||
std::ifstream infile(filename);
|
||||
if (!infile.is_open())
|
||||
{
|
||||
std::cerr << "Error opening " << filename << "!" << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
int id, bk, u, d;
|
||||
double gain;
|
||||
while (infile >> id >> bk >> u >> d >> gain)
|
||||
{
|
||||
backGain[id][bk][u][d] = gain;
|
||||
backGainValid[id][bk][u][d] = true;
|
||||
}
|
||||
|
||||
infile.close();
|
||||
std::cout << "Loaded back gains from " << filename << std::endl;
|
||||
SX3 sx3_contr;
|
||||
}
|
||||
|
||||
Bool_t GainMatchSX3Front::Process(Long64_t entry)
|
||||
{
|
||||
|
||||
b_sx3Multi->GetEntry(entry);
|
||||
b_sx3ID->GetEntry(entry);
|
||||
b_sx3Ch->GetEntry(entry);
|
||||
b_sx3E->GetEntry(entry);
|
||||
b_sx3T->GetEntry(entry);
|
||||
|
||||
sx3.CalIndex();
|
||||
|
||||
std::vector<std::pair<int, int>> ID;
|
||||
for (int i = 0; i < sx3.multi; i++)
|
||||
{
|
||||
|
||||
for (int j = i + 1; j < sx3.multi; j++)
|
||||
{
|
||||
if (sx3.id[i] == 3)
|
||||
hsx3Coin->Fill(sx3.index[i], sx3.index[j]);
|
||||
}
|
||||
if (sx3.e[i] > 100)
|
||||
{
|
||||
ID.push_back(std::pair<int, int>(sx3.id[i], i));
|
||||
hsx3IndexVE->Fill(sx3.index[i], sx3.e[i]);
|
||||
}
|
||||
}
|
||||
|
||||
if (ID.size() > 0)
|
||||
{
|
||||
std::sort(ID.begin(), ID.end(), [](const std::pair<int, int> &a, const std::pair<int, int> &b)
|
||||
{ return a.first < b.first; });
|
||||
|
||||
// start with the first entry in the sorted array: channels that belong to the same detector are together in sequenmce
|
||||
std::vector<std::pair<int, int>> sx3ID;
|
||||
sx3ID.push_back(ID[0]);
|
||||
bool found = false;
|
||||
|
||||
for (size_t i = 1; i < ID.size(); i++)
|
||||
{ // Check if id of i belongs to the same detector and then add it to the detector ID vector
|
||||
if (ID[i].first == sx3ID.back().first)
|
||||
{ // count the nunmber of hits that belong to the same detector
|
||||
sx3ID.push_back(ID[i]);
|
||||
|
||||
if (sx3ID.size() >= 3)
|
||||
{
|
||||
found = true;
|
||||
}
|
||||
}
|
||||
else
|
||||
{ // the next event does not belong to the same detector, abandon the first event and continue with the next one
|
||||
if (!found)
|
||||
{
|
||||
sx3ID.clear();
|
||||
sx3ID.push_back(ID[i]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (found)
|
||||
{
|
||||
int sx3ChUp = -1, sx3ChDn = -1, sx3ChBk = -1;
|
||||
float sx3EUp = 0.0, sx3EDn = 0.0, sx3EBk = 0.0;
|
||||
|
||||
for (size_t i = 0; i < sx3ID.size(); i++)
|
||||
{
|
||||
int index = sx3ID[i].second;
|
||||
// Check the channel number and assign it to the appropriate channel type
|
||||
if (sx3.ch[index] < 8)
|
||||
{
|
||||
if (sx3.ch[index] % 2 == 0)
|
||||
{
|
||||
sx3ChDn = sx3.ch[index] / 2;
|
||||
sx3EDn = sx3.e[index];
|
||||
}
|
||||
else
|
||||
{
|
||||
sx3ChUp = sx3.ch[index] / 2;
|
||||
sx3EUp = sx3.e[index];
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
sx3ChBk = sx3.ch[index] - 8;
|
||||
// if (sx3ChBk == 2)
|
||||
// printf("Found back channel Det %d Back %d \n", sx3.id[index], sx3ChBk);
|
||||
sx3EBk = sx3.e[index];
|
||||
}
|
||||
}
|
||||
// If we have a valid front and back channel, fill the histograms
|
||||
hSX3->Fill(sx3ChDn + 4, sx3ChBk);
|
||||
hSX3->Fill(sx3ChUp, sx3ChBk);
|
||||
|
||||
// Fill the histogram for the front vs back
|
||||
hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
|
||||
for (int i = 0; i < sx3.multi; i++)
|
||||
{
|
||||
if (sx3.id[i] == 3 && sx3.e[i] > 100)
|
||||
{
|
||||
// back gain correction
|
||||
|
||||
// Fill the histogram for the front vs back with gain correction
|
||||
hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
// Fill the index vs energy histogram
|
||||
hsx3IndexVE_g->Fill(sx3.index[i], sx3.e[i]);
|
||||
// }
|
||||
// {
|
||||
TString histName = Form("hSX3FVB_id%d_U%d_D%d_B%d", sx3.id[i], sx3ChUp, sx3ChDn, sx3ChBk);
|
||||
TH2F *hist2d = (TH2F *)gDirectory->Get(histName);
|
||||
if (!hist2d)
|
||||
{
|
||||
hist2d = new TH2F(histName, Form("hSX3FVB_id%d_U%d_D%d_B%d", sx3.id[i], sx3ChUp, sx3ChDn, sx3ChBk), 400, 0, 16000, 400, 0, 16000);
|
||||
}
|
||||
|
||||
// if (sx3ChBk == 2)
|
||||
// printf("Found back channel Det %d Back %d \n", sx3.id[i], sx3ChBk);
|
||||
// hsx3IndexVE_g->Fill(sx3.index[i], sx3.e[i]);
|
||||
// hSX3FvsB_g->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
|
||||
hist2d->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
|
||||
if (cut && cut->IsInside(sx3EUp + sx3EDn, sx3EBk) &&
|
||||
cut1 && cut1->IsInside(sx3EUp / sx3EBk, sx3EDn / sx3EBk))
|
||||
{
|
||||
|
||||
if (backGainValid[sx3.id[i]][sx3ChBk][sx3ChUp][sx3ChDn])
|
||||
{
|
||||
sx3EBk *= backGain[sx3.id[i]][sx3ChBk][sx3ChUp][sx3ChDn];
|
||||
}
|
||||
// Accumulate data for gain matching
|
||||
dataPoints[{sx3.id[i], sx3ChBk, sx3ChUp, sx3ChDn}].emplace_back(sx3EBk, sx3EUp, sx3EDn);
|
||||
}
|
||||
// if (sx3.id[i] < 24 && sx3ChUp < 4 && sx3ChBk < 4 && std::isfinite(sx3EUp) && std::isfinite(sx3EDn) && std::isfinite(sx3EBk))
|
||||
{
|
||||
// Accumulate data for gain matching
|
||||
dataPoints[{sx3.id[i], sx3ChBk, sx3ChUp, sx3ChDn}].emplace_back(sx3EBk, sx3EUp, sx3EDn);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return kTRUE;
|
||||
}
|
||||
|
||||
void GainMatchSX3Front::Terminate()
|
||||
{
|
||||
|
||||
std::map<std::tuple<int, int, int, int>, TVectorD> fitCoefficients;
|
||||
|
||||
// === Gain matching ===
|
||||
|
||||
std::ofstream outFile("sx3_GainMatchfront.txt");
|
||||
if (!outFile.is_open())
|
||||
{
|
||||
std::cerr << "Error opening output file!" << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
TH2F *hUvD = new TH2F("hUvD", " UvD; Up/CorrBack; Down/CorrBack", 600, 0, 1, 600, 0, 1);
|
||||
|
||||
for (const auto &kv : dataPoints)
|
||||
{
|
||||
auto [id, bk, u, d] = kv.first;
|
||||
const auto &pts = kv.second;
|
||||
|
||||
if (pts.size() < 5)
|
||||
continue;
|
||||
|
||||
std::vector<double> udE, corrBkE;
|
||||
|
||||
for (const auto &pr : pts)
|
||||
{
|
||||
double eBkCorr, eUp, eDn;
|
||||
std::tie(eBkCorr, eUp, eDn) = pr;
|
||||
udE.push_back(eUp + eDn);
|
||||
corrBkE.push_back(eBkCorr);
|
||||
hUvD->Fill(eUp / eBkCorr, eDn / eBkCorr);
|
||||
}
|
||||
|
||||
TGraph g(udE.size(), udE.data(), corrBkE.data());
|
||||
TF1 f("f", "[0]*x", 0, 40000);
|
||||
g.Fit(&f, "QNR");
|
||||
|
||||
frontGain[id][bk][u][d] = f.GetParameter(0);
|
||||
frontGainValid[id][bk][u][d] = true;
|
||||
|
||||
outFile << id << " " << bk << " " << u << " " << d << " " << frontGain[id][bk][u][d] << std::endl;
|
||||
printf("Front gain Det%d Back%d Up%dDn%d → %.4f\n", id, bk, u, d, frontGain[id][bk][u][d]);
|
||||
}
|
||||
|
||||
outFile.close();
|
||||
std::cout << "Gain matching complete." << std::endl;
|
||||
|
||||
// === Stage 3: Create corrected histogram ===
|
||||
TH2F *hCorrectedFvB = new TH2F("hCorrectedFvB", "Corrected;Corrected Front Sum;Corrected Back", 800, 0, 16000, 800, 0, 16000);
|
||||
TH2F *hCorrectedUvD = new TH2F("hCorrectedUvD", "Corrected UvD; UvD Up; UvD Down", 600, 0, 1, 600, 0, 1);
|
||||
|
||||
for (const auto &kv : dataPoints)
|
||||
{
|
||||
|
||||
auto [id, bk, u, d] = kv.first;
|
||||
double front = frontGain[id][bk][u][d];
|
||||
|
||||
for (const auto &pr : kv.second)
|
||||
{
|
||||
double eBk, eUp, eDn;
|
||||
std::tie(eBk, eUp, eDn) = pr;
|
||||
double corrUp = eUp * front;
|
||||
double corrDn = eDn * front;
|
||||
|
||||
hCorrectedFvB->Fill(corrUp + corrDn, eBk);
|
||||
hCorrectedUvD->Fill(corrUp / eBk, corrDn / eBk);
|
||||
}
|
||||
}
|
||||
|
||||
// === Final canvas ===
|
||||
gStyle->SetOptStat(1110);
|
||||
TCanvas *c1 = new TCanvas("c1", "Gain Correction Results", 1200, 600);
|
||||
c1->Divide(2, 1);
|
||||
|
||||
c1->cd(1);
|
||||
hSX3FvsB_g->SetTitle("Before Correction (Gated)");
|
||||
hSX3FvsB_g->GetXaxis()->SetTitle("Measured Front Sum (E_Up + E_Dn)");
|
||||
hSX3FvsB_g->GetYaxis()->SetTitle("Measured Back E");
|
||||
hSX3FvsB_g->Draw("colz");
|
||||
|
||||
c1->cd(2);
|
||||
hCorrectedFvB->SetTitle("After Correction");
|
||||
hCorrectedFvB->Draw("colz");
|
||||
TF1 *diag = new TF1("diag", "x", 0, 40000);
|
||||
diag->SetLineColor(kRed);
|
||||
diag->SetLineWidth(2);
|
||||
diag->Draw("same");
|
||||
|
||||
std::cout << "Terminate() completed successfully." << std::endl;
|
||||
}
|
||||
104
GainMatchSX3Front.h
Normal file
104
GainMatchSX3Front.h
Normal file
|
|
@ -0,0 +1,104 @@
|
|||
#ifndef GainMatchSX3Front_h
|
||||
#define GainMatchSX3Front_h
|
||||
|
||||
#include <TROOT.h>
|
||||
#include <TChain.h>
|
||||
#include <TFile.h>
|
||||
#include <TSelector.h>
|
||||
|
||||
#include "Armory/ClassDet.h"
|
||||
|
||||
class GainMatchSX3Front : public TSelector {
|
||||
public :
|
||||
TTree *fChain; //!pointer to the analyzed TTree or TChain
|
||||
|
||||
// Fixed size dimensions of array or collections stored in the TTree if any.
|
||||
|
||||
// Declaration of leaf types
|
||||
Det sx3;
|
||||
Det qqq;
|
||||
Det pc ;
|
||||
|
||||
ULong64_t evID;
|
||||
UInt_t run;
|
||||
|
||||
// List of branches
|
||||
TBranch *b_eventID; //!
|
||||
TBranch *b_run; //!
|
||||
TBranch *b_sx3Multi; //!
|
||||
TBranch *b_sx3ID; //!
|
||||
TBranch *b_sx3Ch; //!
|
||||
TBranch *b_sx3E; //!
|
||||
TBranch *b_sx3T; //!
|
||||
TBranch *b_qqqMulti; //!
|
||||
TBranch *b_qqqID; //!
|
||||
TBranch *b_qqqCh; //!
|
||||
TBranch *b_qqqE; //!
|
||||
TBranch *b_qqqT; //!
|
||||
TBranch *b_pcMulti; //!
|
||||
TBranch *b_pcID; //!
|
||||
TBranch *b_pcCh; //!
|
||||
TBranch *b_pcE; //!
|
||||
TBranch *b_pcT; //!
|
||||
|
||||
GainMatchSX3Front(TTree * /*tree*/ =0) : fChain(0) { }
|
||||
virtual ~GainMatchSX3Front() { }
|
||||
virtual Int_t Version() const { return 2; }
|
||||
virtual void Begin(TTree *tree);
|
||||
virtual void SlaveBegin(TTree *tree);
|
||||
virtual void Init(TTree *tree);
|
||||
virtual Bool_t Notify();
|
||||
virtual Bool_t Process(Long64_t entry);
|
||||
virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
|
||||
virtual void SetOption(const char *option) { fOption = option; }
|
||||
virtual void SetObject(TObject *obj) { fObject = obj; }
|
||||
virtual void SetInputList(TList *input) { fInput = input; }
|
||||
virtual TList *GetOutputList() const { return fOutput; }
|
||||
virtual void SlaveTerminate();
|
||||
virtual void Terminate();
|
||||
|
||||
ClassDef(GainMatchSX3Front,0);
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
#ifdef GainMatchSX3Front_cxx
|
||||
void GainMatchSX3Front::Init(TTree *tree){
|
||||
|
||||
// Set branch addresses and branch pointers
|
||||
if (!tree) return;
|
||||
fChain = tree;
|
||||
fChain->SetMakeClass(1);
|
||||
|
||||
fChain->SetBranchAddress("evID", &evID, &b_eventID);
|
||||
fChain->SetBranchAddress("run", &run, &b_run);
|
||||
|
||||
sx3.SetDetDimension(24,12);
|
||||
qqq.SetDetDimension(4,32);
|
||||
pc.SetDetDimension(2,24);
|
||||
|
||||
fChain->SetBranchAddress("sx3Multi", &sx3.multi, &b_sx3Multi);
|
||||
fChain->SetBranchAddress("sx3ID", &sx3.id, &b_sx3ID);
|
||||
fChain->SetBranchAddress("sx3Ch", &sx3.ch, &b_sx3Ch);
|
||||
fChain->SetBranchAddress("sx3E", &sx3.e, &b_sx3E);
|
||||
fChain->SetBranchAddress("sx3T", &sx3.t, &b_sx3T);
|
||||
|
||||
}
|
||||
|
||||
Bool_t GainMatchSX3Front::Notify(){
|
||||
|
||||
return kTRUE;
|
||||
}
|
||||
|
||||
void GainMatchSX3Front::SlaveBegin(TTree * /*tree*/){
|
||||
|
||||
TString option = GetOption();
|
||||
|
||||
}
|
||||
|
||||
void GainMatchSX3Front::SlaveTerminate(){
|
||||
|
||||
}
|
||||
|
||||
|
||||
#endif // #ifdef GainMatchSX3Front_cxx
|
||||
245
GainMatchSX3Front1.C
Normal file
245
GainMatchSX3Front1.C
Normal file
|
|
@ -0,0 +1,245 @@
|
|||
#define GainMatchSX3_cxx
|
||||
|
||||
#include "GainMatchSX3.h"
|
||||
#include "Armory/ClassSX3.h"
|
||||
#include <TFile.h>
|
||||
#include <TTree.h>
|
||||
#include <TGraph.h>
|
||||
#include <TF1.h>
|
||||
#include <TH2F.h>
|
||||
#include <TCanvas.h>
|
||||
#include <TStyle.h>
|
||||
#include <TApplication.h>
|
||||
#include <map>
|
||||
#include <vector>
|
||||
#include <tuple>
|
||||
#include <fstream>
|
||||
#include <iostream>
|
||||
#include <algorithm>
|
||||
|
||||
// Constants
|
||||
const int MAX_DET = 24;
|
||||
const int MAX_BK = 4;
|
||||
const int MAX_UP = 4;
|
||||
const int MAX_DOWN = 4;
|
||||
|
||||
// Gain arrays
|
||||
double backGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
|
||||
bool backGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
|
||||
|
||||
double frontGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
|
||||
bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
|
||||
|
||||
// Data container
|
||||
std::map<std::tuple<int, int, int, int>, std::vector<std::tuple<double, double, double>>> dataPoints;
|
||||
|
||||
// Load back gains
|
||||
void LoadBackGains(const std::string &filename)
|
||||
{
|
||||
std::ifstream infile(filename);
|
||||
if (!infile.is_open())
|
||||
{
|
||||
std::cerr << "Error opening " << filename << "!" << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
int id, bk, u, d;
|
||||
double gain;
|
||||
while (infile >> id >> bk >> u >> d >> gain)
|
||||
{
|
||||
backGain[id][bk][u][d] = gain;
|
||||
backGainValid[id][bk][u][d] = true;
|
||||
}
|
||||
|
||||
infile.close();
|
||||
std::cout << "Loaded back gains from " << filename << std::endl;
|
||||
SX3 sx3_contr;
|
||||
}
|
||||
|
||||
// Front gain matching function
|
||||
Bool_t GainMatchSX3::Process(Long64_t entry)
|
||||
{
|
||||
// Link SX3 branches
|
||||
|
||||
b_sx3Multi->GetEntry(entry);
|
||||
b_sx3ID->GetEntry(entry);
|
||||
b_sx3Ch->GetEntry(entry);
|
||||
b_sx3E->GetEntry(entry);
|
||||
b_sx3T->GetEntry(entry);
|
||||
|
||||
sx3.CalIndex();
|
||||
|
||||
Long64_t nentries = tree->GetEntries(Long64_t entry);
|
||||
std::cout << "Total entries: " << nentries << std::endl;
|
||||
|
||||
TH2F *hBefore = new TH2F("hBefore", "Before Correction;E_Up+E_Dn;Back Energy", 400, 0, 40000, 400, 0, 40000);
|
||||
TH2F *hAfter = new TH2F("hAfter", "After Correction;E_Up+E_Dn;Corrected Back Energy", 400, 0, 40000, 400, 0, 40000);
|
||||
|
||||
for (Long64_t entry = 0; entry < nentries; ++entry)
|
||||
{
|
||||
tree->GetEntry(entry);
|
||||
sx3.CalIndex();
|
||||
|
||||
std::vector<std::pair<int, int>> ID;
|
||||
|
||||
for (int i = 0; i < sx3.multi; i++)
|
||||
{
|
||||
if (sx3.e[i] > 100)
|
||||
{
|
||||
ID.push_back({sx3.id[i], i});
|
||||
}
|
||||
}
|
||||
|
||||
if (ID.empty())
|
||||
continue;
|
||||
|
||||
// Sort by id
|
||||
std::sort(ID.begin(), ID.end(), [](auto &a, auto &b) { return a.first < b.first; });
|
||||
|
||||
std::vector<std::pair<int, int>> sx3ID;
|
||||
sx3ID.push_back(ID[0]);
|
||||
bool found = false;
|
||||
|
||||
for (size_t i = 1; i < ID.size(); i++)
|
||||
{
|
||||
if (ID[i].first == sx3ID.back().first)
|
||||
{
|
||||
sx3ID.push_back(ID[i]);
|
||||
if (sx3ID.size() >= 3)
|
||||
found = true;
|
||||
}
|
||||
else if (!found)
|
||||
{
|
||||
sx3ID.clear();
|
||||
sx3ID.push_back(ID[i]);
|
||||
}
|
||||
}
|
||||
|
||||
if (!found)
|
||||
continue;
|
||||
|
||||
int sx3ChUp = -1, sx3ChDn = -1, sx3ChBk = -1;
|
||||
float sx3EUp = 0.0, sx3EDn = 0.0, sx3EBk = 0.0;
|
||||
int sx3id = sx3ID[0].first;
|
||||
|
||||
for (auto &[id, idx] : sx3ID)
|
||||
{
|
||||
if (sx3.ch[idx] < 8)
|
||||
{
|
||||
if (sx3.ch[idx] % 2 == 0)
|
||||
{
|
||||
sx3ChDn = sx3.ch[idx] / 2;
|
||||
sx3EDn = sx3.e[idx];
|
||||
}
|
||||
else
|
||||
{
|
||||
sx3ChUp = sx3.ch[idx] / 2;
|
||||
sx3EUp = sx3.e[idx];
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
sx3ChBk = sx3.ch[idx] - 8;
|
||||
sx3EBk = sx3.e[idx];
|
||||
}
|
||||
}
|
||||
|
||||
if (sx3ChUp < 0 || sx3ChDn < 0 || sx3ChBk < 0)
|
||||
continue;
|
||||
|
||||
if (!backGainValid[sx3id][sx3ChBk][sx3ChUp][sx3ChDn])
|
||||
continue;
|
||||
|
||||
double corrBk = sx3EBk * backGain[sx3id][sx3ChBk][sx3ChUp][sx3ChDn];
|
||||
|
||||
hBefore->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
hAfter->Fill(sx3EUp + sx3EDn, corrBk);
|
||||
|
||||
dataPoints[{sx3id, sx3ChBk, sx3ChUp, sx3ChDn}].emplace_back(corrBk, sx3EUp, sx3EDn);
|
||||
}
|
||||
|
||||
// === Fit front gains ===
|
||||
std::ofstream outFile("sx3_GainMatchfront.txt");
|
||||
if (!outFile.is_open())
|
||||
{
|
||||
std::cerr << "Error opening sx3_GainMatchfront.txt!" << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
for (const auto &kv : dataPoints)
|
||||
{
|
||||
auto [id, bk, u, d] = kv.first;
|
||||
const auto &pts = kv.second;
|
||||
|
||||
if (pts.size() < 5)
|
||||
continue;
|
||||
|
||||
std::vector<double> udE, corrBkE;
|
||||
|
||||
for (const auto &pr : pts)
|
||||
{
|
||||
double eBkCorr, eUp, eDn;
|
||||
std::tie(eBkCorr, eUp, eDn) = pr;
|
||||
udE.push_back(eUp + eDn);
|
||||
corrBkE.push_back(eBkCorr);
|
||||
}
|
||||
|
||||
TGraph g(udE.size(), udE.data(), corrBkE.data());
|
||||
TF1 f("f", "[0]*x", 0, 40000);
|
||||
g.Fit(&f, "QNR");
|
||||
|
||||
frontGain[id][bk][u][d] = f.GetParameter(0);
|
||||
frontGainValid[id][bk][u][d] = true;
|
||||
|
||||
outFile << id << " " << bk << " " << u << " " << d << " " << frontGain[id][bk][u][d] << std::endl;
|
||||
printf("Front gain Det%d Back%d Up%dDn%d → %.4f\n", id, bk, u, d, frontGain[id][bk][u][d]);
|
||||
}
|
||||
|
||||
outFile.close();
|
||||
std::cout << "Front gain matching complete." << std::endl;
|
||||
|
||||
// === Draw diagnostic plots ===
|
||||
gStyle->SetOptStat(1110);
|
||||
TCanvas *c = new TCanvas("c", "Gain Matching Diagnostics", 1200, 600);
|
||||
c->Divide(2, 1);
|
||||
|
||||
c->cd(1);
|
||||
hBefore->Draw("colz");
|
||||
TF1 *diag1 = new TF1("diag1", "x", 0, 40000);
|
||||
diag1->SetLineColor(kRed);
|
||||
diag1->Draw("same");
|
||||
|
||||
c->cd(2);
|
||||
hAfter->Draw("colz");
|
||||
TF1 *diag2 = new TF1("diag2", "x", 0, 40000);
|
||||
diag2->SetLineColor(kRed);
|
||||
diag2->Draw("same");
|
||||
}
|
||||
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
TApplication app("app", &argc, argv);
|
||||
|
||||
// Load back gains
|
||||
LoadBackGains("sx3_GainMatchback.txt");
|
||||
|
||||
// Open tree
|
||||
TFile *f = TFile::Open("input_tree.root"); // <<< Change file name
|
||||
if (!f || f->IsZombie())
|
||||
{
|
||||
std::cerr << "Cannot open input_tree.root!" << std::endl;
|
||||
return 1;
|
||||
}
|
||||
TTree *tree = (TTree *)f->Get("tree");
|
||||
if (!tree)
|
||||
{
|
||||
std::cerr << "Tree not found!" << std::endl;
|
||||
return 1;
|
||||
}
|
||||
|
||||
// Run front gain matching
|
||||
GainMatchSX3(tree);
|
||||
|
||||
app.Run();
|
||||
return 0;
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user