rtang/ANASEN_analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

modified: GainMatch.C

Browse Source
This commit is contained in:
Vignesh Sitaraman 2025-04-24 10:48:14 -04:00
parent 1df7470ca1
commit 65ab69ebe6
1 changed files with 143 additions and 89 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

232
GainMatch.C
View File

@ -7,7 +7,7 @@
#include <TCanvas.h>
#include <TMath.h>
#include <TCutG.h>
#include <fstream>
#include <utility>
#include <algorithm>
@ -26,28 +26,30 @@ PW pw_contr;
TVector3 hitPos;
bool HitNonZero;
const int MAX_DET = 4;
const int MAX_WEDGE = 16;
const int MAX_RING = 16;
bool gainValid[MAX_DET][MAX_RING][MAX_WEDGE] = {{{false}}};
double gainArray[MAX_DET][MAX_RING][MAX_WEDGE] = {{{0}}};
void GainMatch::Begin(TTree * /*tree*/)
{
TString option = GetOption();
hSX3FvsB = new TH2F("hSX3FvsB", "SX3 Front vs Back; Front E; Back E", 400, 0, 16000, 400, 0, 16000);
hQQQFVB = new TH2F("hQQQFVB", "number of good QQQ vs QQQ id", 400, 0, 16000, 400, 0, 16000);
sx3_contr.ConstructGeo();
pw_contr.ConstructGeo();
}
Bool_t GainMatch::Process(Long64_t entry)
{
// if ( entry > 100 ) return kTRUE;
hitPos.Clear();
HitNonZero = false;
// if( entry > 1) return kTRUE;
// printf("################### ev : %llu \n", entry);
b_sx3Multi->GetEntry(entry);
b_sx3ID->GetEntry(entry);
b_sx3Ch->GetEntry(entry);
@ -68,41 +70,26 @@ Bool_t GainMatch::Process(Long64_t entry)
qqq.CalIndex();
pc.CalIndex();
// sx3.Print();
// ########################################################### Raw data
// //======================= SX3
// Initialize vector to store pairs of ID and index
std::vector<std::pair<int, int>> ID; // first = id, second = index
// Loop through each entry in sx3.multi
std::vector<std::pair<int, int>> ID;
for (int i = 0; i < sx3.multi; i++)
{
// Store ID and index pair in ID vector
ID.push_back(std::pair<int, int>(sx3.id[i], i));
}
// Check if the ID vector is not empty
if (ID.size() > 0)
{
// Sort the ID vector by the first element (ID)
std::sort(ID.begin(), ID.end(), [](const std::pair<int, int> &a, const std::pair<int, int> &b)
{ return a.first < b.first; });
// Create a new vector to store IDs that are the same
std::vector<std::pair<int, int>> sx3ID;
sx3ID.push_back(ID[0]); // Start with the first ID
sx3ID.push_back(ID[0]);
bool found = false;
// Loop through the sorted IDs and group by same IDs
for (size_t i = 1; i < ID.size(); i++)
{
if (ID[i].first == sx3ID.back().first)
{
// If the ID matches the last one in sx3ID, add to the group
sx3ID.push_back(ID[i]);
// If 3 or more IDs are grouped, set found to true
if (sx3ID.size() >= 3)
{
found = true;
@ -110,7 +97,6 @@ Bool_t GainMatch::Process(Long64_t entry)
}
else
{
// If a new ID is encountered and no group is found, reset the group
if (!found)
{
sx3ID.clear();
@ -119,21 +105,16 @@ Bool_t GainMatch::Process(Long64_t entry)
}
}
// If a group of 3 or more IDs is found, process the channels and energies
if (found)
{
int sx3ChUp = -1, sx3ChDn = -1, sx3ChBk = -1;
float sx3EUp = 0.0, sx3EDn = 0.0, sx3EBk = 0.0;
// Loop through the grouped IDs
for (size_t i = 0; i < sx3ID.size(); i++)
{
int index = sx3ID[i].second; // Get the index from the pair
// Categorize channel and energy
int index = sx3ID[i].second;
if (sx3.ch[index] < 8)
{
// If channel is less than 8, assign it to up or down
if (sx3.ch[index] % 2 == 0)
{
sx3ChDn = sx3.ch[index];
@ -147,88 +128,99 @@ Bool_t GainMatch::Process(Long64_t entry)
}
else
{
// If channel is >= 8, assign it to back
sx3ChBk = sx3.ch[index];
sx3EBk = sx3.e[index];
}
}
// make a histogram for Sx3 back energy vs sum of sx3 front energies
hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk);
// Further energy matching calculations can be added here...
}
}
// //======================= QQQ
TH2F *hist1 = NULL;
// Loop through each entry in qqq.multi
for (int i = 0; i < qqq.multi; i++)
{
for (int j = i + 1; j < qqq.multi; j++)
{
// if( qqq.used[i] == true ) continue;
// if( qqq.id[i] == qqq.id[j] && (16 - qqq.ch[i]) * (16 - qqq.ch[j]) < 0 ){ // must be same detector and wedge and ring
if (qqq.id[i] == qqq.id[j])
{ // must be same detector
{
int chWedge = -1;
int chRing = -1;
float eWedge = 0.0;
float eRing = 0.0;
if (qqq.ch[i] < qqq.ch[j])
if (qqq.ch[i] < 16 && qqq.ch[j] >= 16)
{
chRing = qqq.ch[j] - 16;
eRing = qqq.e[j];
chWedge = qqq.ch[i];
eWedge = qqq.e[i];
chRing = qqq.ch[j] - 16;
eRing = qqq.e[j];
}
else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16)
{
chWedge = qqq.ch[j];
eWedge = qqq.e[j];
chRing = qqq.ch[i] - 16;
eRing = qqq.e[i];
}
else
{
chRing = qqq.ch[i];
eRing = qqq.e[i];
chWedge = qqq.ch[j] - 16;
eWedge = qqq.e[j];
}
// printf(" ID : %d , chWedge : %d, chRing : %d \n", qqq.id[i], chWedge, chRing);
hQQQFVB->Fill(qqq.e[i], qqq.e[j]);
// 1. Create/get individual 2D histogram
TString histName = Form("hQQQFVB_r%d_w%d_id%d", chRing, chWedge, qqq.id[i]);
continue;
hQQQFVB->Fill(eWedge, eRing);
TString histName = Form("hQQQFVB_id%d_r%d_w%d", qqq.id[i], chRing, chWedge);
TH2F *hist2d = (TH2F *)gDirectory->Get(histName);
if (!hist2d)
{
hist2d = new TH2F(histName, Form("QQQ GainMatch ID%d R%d W%d;Wedge E;Ring E", qqq.id[i], chRing, chWedge),
400, 0, 16000, 400, 0, 16000);
hist2d = new TH2F(histName, Form("QQQ GainMatch Det%d R%d W%d;Wedge E;Ring E", qqq.id[i], chRing, chWedge), 400, 0, 16000, 400, 0, 16000);
}
hist2d->Fill(eWedge, eRing);
}
}
}
return kTRUE;
}
void GainMatch::Terminate()
{
TFile *cutFile = TFile::Open("qqqcorr.root");
TCutG *cut = (TCutG *)cutFile->Get("qqqcorr");
const int MAX_DET = 4;
const int MAX_RING = 16;
const int MAX_WEDGE = 16;
// Store gains and validity
static double gainArray[MAX_DET][MAX_RING][MAX_WEDGE] = {{{0}}};
static bool gainValid[MAX_DET][MAX_RING][MAX_WEDGE] = {{{false}}};
std::ofstream outFile("qqq_gainmatch.txt");
if (!outFile.is_open())
{
printf("Error opening output file!");
return;
}
// Collect all (wedgeE, ringE) pairs per detector/ring/wedge
std::map<std::tuple<int, int, int>, std::vector<std::pair<double, double>>> dataPoints;
TIter next(gDirectory->GetList());
TObject *obj;
while ((obj = next()))
{
if (!obj->InheritsFrom("TH2F"))
continue;
TH2F *hist2d = (TH2F *)obj;
TH2F *hist2d = static_cast<TH2F *>(obj);
TString name = hist2d->GetName();
if (!name.BeginsWith("hQQQFVB_r"))
if (!name.BeginsWith("hQQQFVB_id"))
continue;
if (hist2d->GetEntries() < 100)
continue;
if (hist2d->GetEntries() < 1000)
continue;
std::vector<double> wedge_vals;
std::vector<double> peak_vals;
int id, ring, wedge;
sscanf(name.Data(), "hQQQFVB_id%d_r%d_w%d", &id, &ring, &wedge);
for (int binX = 1; binX <= hist2d->GetNbinsX(); ++binX)
{
TH1D *projY = hist2d->ProjectionY("_py", binX, binX);
TH1D *projY = hist2d->ProjectionY("_py", binX, binX + 1);
if (projY->GetEntries() < 30)
{
delete projY;
@ -240,34 +232,96 @@ void GainMatch::Terminate()
projY->Fit(fit, "QNR");
double wedgeE = hist2d->GetXaxis()->GetBinCenter(binX);
double ringPeak = fit->GetParameter(1);
double ringE = fit->GetParameter(1);
dataPoints[{id, ring, wedge}].emplace_back(wedgeE, ringE);
wedge_vals.push_back(wedgeE);
peak_vals.push_back(ringPeak);
delete projY;
delete fit;
delete projY;
}
}
// Fit slopes with sigma-clipping
for (auto &kv : dataPoints)
{
auto [id, ring, wedge] = kv.first;
auto &pts = kv.second;
if (pts.size() < 5)
continue;
// Build vectors
std::vector<double> wE, rE;
for (auto &pr : pts)
{
wE.push_back(pr.first);
rE.push_back(pr.second);
}
if (wedge_vals.size() >= 5)
// Initial fit
TGraph *g0 = new TGraph(wE.size(), wE.data(), rE.data());
TF1 *f0 = new TF1("f0", "[0]*x + [1]", 0, 16000);
g0->Fit(f0, "QNR");
double m0 = f0->GetParameter(0), b0 = f0->GetParameter(1);
// Clip to cut
std::vector<double> wC, rC;
for (size_t i = 0; i < wE.size(); ++i)
{
TGraph *gr = new TGraph(wedge_vals.size(), &wedge_vals[0], &peak_vals[0]);
gr->SetName(name + "_fit");
TF1 *line = new TF1("line", "pol1", 0, 16000);
gr->Fit(line, "Q");
if (cut->IsInside(wE[i], rE[i]))
{
wC.push_back(wE[i]);
rC.push_back(rE[i]);
}
}
double gain = line->GetParameter(1);
double offset = line->GetParameter(0);
printf("Gain match %s → Gain = %.4f, Offset = %.2f\n", name.Data(), gain, offset);
// Final fit
if (wC.size() >= 5)
{
TGraph *g1 = new TGraph(wC.size(), wC.data(), rC.data());
TF1 *f1 = new TF1("f1", "[0]*x + [1]", 0, 16000);
g1->Fit(f1, "QNR");
gainArray[id][ring][wedge] = f1->GetParameter(0);
gainValid[id][ring][wedge] = true;
delete g1;
delete f1;
}
delete g0;
delete f0;
}
TCanvas *c1 = new TCanvas(name + "_c", name + "_c");
gr->SetTitle(Form("Gain Match: %s", name.Data()));
gr->GetXaxis()->SetTitle("Wedge Energy");
gr->GetYaxis()->SetTitle("Ring Energy");
gr->Draw("AP");
line->Draw("same");
c1->SaveAs(Form("%s_fit.png", name.Data()));
delete c1;
// Write out valid entries
for (int id = 0; id < MAX_DET; ++id)
{
for (int ring = 0; ring < MAX_RING; ++ring)
{
for (int wedge = 0; wedge < MAX_WEDGE; ++wedge)
{
if (!gainValid[id][ring][wedge])
continue;
outFile << id << " " << wedge << " " << ring << " " << gainArray[id][ring][wedge] << std::endl;
// printf("Gain match Det%d Ring%d Wedge%d → %.4f \n", id, ring, wedge, gainArray[id][ring][wedge]);
}
}
}
outFile.close();
printf("Gain matching complete (with sigma-clipping)");
// === Plot all gain-matched QQQ points together with a 2D histogram ===
TH2F *hAll = new TH2F("hAll", "All QQQ Gain-Matched;Corrected Wedge E;Ring E",
400, 0, 16000, 400, 0, 16000);
// Fill the combined TH2F with corrected data
for (auto &kv : dataPoints)
{
int id, ring, wedge;
std::tie(id, ring, wedge) = kv.first;
if (!gainValid[id][ring][wedge])
continue;
auto &pts = kv.second;
for (auto &pr : pts)
{
double corrWedge = pr.first * gainArray[id][ring][wedge];
double ringE = pr.second;
hAll->Fill(corrWedge, ringE);
}
}
}