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

modified: .gitignore

Browse Source 
modified:   .vscode/settings.json
	modified:   Analyzer.C
	modified:   Calibration.C
	modified:   GainMatchQQQ.C
	modified:   QQQ_Calcheck.C
This commit is contained in:
Vignesh Sitaraman 2025-12-16 15:41:21 -05:00
parent aee3a2467d
commit 97880940be
6 changed files with 657 additions and 606 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
.gitignore vendored
View File

@ -10,6 +10,7 @@ AnasenMS
data/ data/
data_proton/ data_proton/
Sudarshan/
Analyzer_C_ACLiC_dict0713aaa966_dictContent.h Analyzer_C_ACLiC_dict0713aaa966_dictContent.h
.gitignore .gitignore
Analyzer_C_ACLiC_dict5411fecd5c_dictUmbrella.h Analyzer_C_ACLiC_dict5411fecd5c_dictUmbrella.h
@ -19,3 +20,7 @@ MakePlotsQQQ.C
MakePlotsQQQ.h MakePlotsQQQ.h
MakePlotsSX3.C MakePlotsSX3.C
MakePlotsSX3.h MakePlotsSX3.h
qqq_gains_det3.dat
qqq_relative_gains.dat
Armory/CorrelateQQQ.h
QQQStage2.C

5
.vscode/settings.json vendored
View File

@ -121,7 +121,10 @@
"MakePlotsSX3.C": "cpp", "MakePlotsSX3.C": "cpp",
"QQQ_Calibcheck.C": "cpp", "QQQ_Calibcheck.C": "cpp",
"QQQ_Calcheck.C": "cpp", "QQQ_Calcheck.C": "cpp",
"makeplots.C": "cpp" "makeplots.C": "cpp",
"GlobalMinimizeQQQ.C": "cpp",
"QQQStage2.C": "cpp",
"inspect.C": "cpp"
}, },
"github-enterprise.uri": "https://fsunuc.physics.fsu.edu" "github-enterprise.uri": "https://fsunuc.physics.fsu.edu"
} }

921
Analyzer.C
View File

@ -48,7 +48,7 @@ PW pwinstance;
TVector3 hitPos; TVector3 hitPos;
// TVector3 anodeIntersection; // TVector3 anodeIntersection;
std::map<int, std::pair<double, double>> slopeInterceptMap; std::map<int, std::pair<double, double>> slopeInterceptMap;
// SX3 Calibration Arrays
const int MAX_DET = 24; const int MAX_DET = 24;
const int MAX_UP = 4; const int MAX_UP = 4;
const int MAX_DOWN = 4; const int MAX_DOWN = 4;
@ -58,6 +58,15 @@ bool backGainValid[MAX_DET][MAX_BK] = {{false}};
double frontGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}}; double frontGain[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{0}}}};
bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}}; bool frontGainValid[MAX_DET][MAX_BK][MAX_UP][MAX_DOWN] = {{{{false}}}};
// QQQ Calibration Arrays
const int MAX_QQQ = 4;
const int MAX_RING = 16;
const int MAX_WEDGE = 16;
double qqqGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
bool qqqGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
double qqqCalib[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
bool qqqCalibValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
bool HitNonZero; bool HitNonZero;
bool sx3ecut; bool sx3ecut;
bool qqqEcut; bool qqqEcut;
@ -187,6 +196,52 @@ void Analyzer::Begin(TTree * /*tree*/)
frontGain[idf][bkf][uf][df] = fgain; frontGain[idf][bkf][uf][df] = fgain;
frontGainValid[idf][bkf][uf][df] = true; frontGainValid[idf][bkf][uf][df] = true;
} }
// QQQ Gain Matching and Calibration
// ----------------------- Load QQQ Gains
{
std::string filename = "qqq_GainMatch.txt";
std::ifstream infile(filename);
if (!infile.is_open())
{
std::cerr << "Error opening " << filename << "!" << std::endl;
}
else
{
int det, ring, wedge;
double gainw, gainr;
while (infile >> det >> ring >> wedge >> gainw >> gainr)
{
qqqGain[det][ring][wedge] = gainw;
// qqqrGain[det][ring][wedge] = gainr;
qqqGainValid[det][ring][wedge] = (gainw > 0);
// qqqrGainValid[det][ring][wedge] = (gainr > 0);
}
infile.close();
std::cout << "Loaded QQQ gains from " << filename << std::endl;
}
}
// ----------------------- Load QQQ Calibrations
{
std::string filename = "qqq_Calib.txt";
std::ifstream infile(filename);
if (!infile.is_open())
{
std::cerr << "Error opening " << filename << "!" << std::endl;
}
else
{
int det, ring, wedge;
double slope;
while (infile >> det >> ring >> wedge >> slope)
{
qqqCalib[det][ring][wedge] = slope;
qqqCalibValid[det][ring][wedge] = (slope > 0);
}
infile.close();
std::cout << "Loaded QQQ calibrations from " << filename << std::endl;
}
}
} }
Bool_t Analyzer::Process(Long64_t entry) Bool_t Analyzer::Process(Long64_t entry)
@ -323,9 +378,9 @@ Bool_t Analyzer::Process(Long64_t entry)
} }
} }
sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, sx3EUp, sx3EDn); // sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, sx3EUp, sx3EDn);
hitPos = sx3_contr.GetHitPos(); // hitPos = sx3_contr.GetHitPos();
HitNonZero = true; // HitNonZero = true;
// hitPos.Print(); // hitPos.Print();
} }
} }
@ -372,431 +427,463 @@ Bool_t Analyzer::Process(Long64_t entry)
if (qqq.id[i] == qqq.id[j]) if (qqq.id[i] == qqq.id[j])
{ // must be same detector { // must be same detector
int chWedge = -1; if (qqq.e[i] > 100)
int chRing = -1; qqqEcut = true;
if (qqq.ch[i] < qqq.ch[j]) if (qqq.id[i] == qqq.id[j])
{ {
chRing = qqq.ch[j] - 16; int chWedge = -1;
chWedge = qqq.ch[i]; int chRing = -1;
} float eWedgeRaw = 0.0;
else float eWedge = 0.0;
{ float eWedgeMeV = 0.0;
chRing = qqq.ch[i]; float eRingRaw = 0.0;
chWedge = qqq.ch[j] - 16; float eRing = 0.0;
} float eRingMeV = 0.0;
// printf(" ID : %d , chWedge : %d, chRing : %d \n", qqq.id[i], chWedge, chRing); // plug in gains
if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && /*qqqrGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16] &&*/ qqqGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16])
double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
double rho = 50. + 40. / 16. * (chRing + 0.5);
// if(qqq.e[i]>50){
hqqqPolar->Fill(theta, rho);
// }
// qqq.used[i] = true;
// qqq.used[j] = true;
if (!HitNonZero)
{
double x = rho * TMath::Cos(theta);
double y = rho * TMath::Sin(theta);
hitPos.SetXYZ(x, y, 23 + 75 + 30);
HitNonZero = true;
}
}
}
}
// //======================= PC
// Calculate the crossover points and put them into an array
pwinstance.ConstructGeo();
Coord Crossover[24][24][2];
TVector3 a, c, diff;
double a2, ac, c2, adiff, cdiff, denom, alpha, beta;
int index = 0;
for (int i = 0; i < pwinstance.An.size(); i++)
{
a = pwinstance.An[i].first - pwinstance.An[i].second;
for (int j = 0; j < pwinstance.Ca.size(); j++)
{
// Ok so this method uses what is essentially the solution of 2 equations to find the point of intersection between the anode and cathode wires
// here a and c are the vectors of the anode and cathode wires respectively
// diff is the perpendicular vector between the anode and cathode wires
// The idea behind this is to then find the scalars alpha and beta that give a ratio between 0 and -1,
c = pwinstance.Ca[j].first - pwinstance.Ca[j].second;
diff = pwinstance.An[i].first - pwinstance.Ca[j].first;
a2 = a.Dot(a);
c2 = c.Dot(c);
ac = a.Dot(c);
adiff = a.Dot(diff);
cdiff = c.Dot(diff);
denom = a2 * c2 - ac * ac;
alpha = (ac * cdiff - c2 * adiff) / denom;
beta = (a2 * cdiff - ac * adiff) / denom;
Crossover[i][j][0].x = pwinstance.An[i].first.X() + alpha * a.X();
Crossover[i][j][0].y = pwinstance.An[i].first.Y() + alpha * a.Y();
Crossover[i][j][0].z = pwinstance.An[i].first.Z() + alpha * a.Z();
if (Crossover[i][j][0].z < -190 || Crossover[i][j][0].z > 190)
{
Crossover[i][j][0].z = 9999999;
}
// placeholder variable Crossover[i][j][1].x has nothing to do with the geometry of the crossover and is being used to store the alpha value-
//-so that it can be used to sort "good" hits later
Crossover[i][j][1].x = alpha;
Crossover[i][j][1].y = 0;
// if(i==0){
// printf("CID, Crossover z and alpha are : %d %f %f \n", j, Crossover[i][j][0].z, Crossover[i][j][1].x /*this is alpha*/);
// }
// }
// }
}
}
// printf("Anode and cathode indices, alpha, denom, andiff, cndiff : %d %d %f %f %f %f\n", i, j, alpha, denom, adiff, cdiff);
// anodeIntersection.Clear();
for (int i = 0; i < pc.multi; i++)
{
if (pc.e[i] > 100)
{
hpcIndexVE->Fill(pc.index[i], pc.e[i]); // non gain matched energy
}
// Gain Matching of PC wires
if (pc.index[i] >= 0 && pc.index[i] < 48)
{
// printf("index: %d, Old cathode energy: %d \n", pc.index[i],pc.e[i]);
auto it = slopeInterceptMap.find(pc.index[i]);
if (it != slopeInterceptMap.end())
{
double slope = it->second.first;
double intercept = it->second.second;
// printf("slope: %f, intercept:%f\n" ,slope, intercept);
pc.e[i] = slope * pc.e[i] + intercept;
// printf("index: %d, New cathode energy: %d \n",pc.index[i], pc.e[i]);
}
hpcIndexVE_GM->Fill(pc.index[i], pc.e[i]);
// hPC_E[pc.index[i]]->Fill(pc.e[i]); // gain matched energy per channel
}
}
std::vector<std::pair<int, double>> anodeHits = {};
std::vector<std::pair<int, double>> cathodeHits = {};
std::vector<std::pair<int, double>> corrcatMax = {};
std::vector<std::pair<int, double>> corrcatnextMax = {};
std::vector<std::pair<int, double>> commcat = {};
int aID = 0;
int cID = 0;
float aE = 0;
float cE = 0;
float aESum = 0;
float cESum = 0;
float aEMax = 0;
float cEMax = 0;
float aEnextMax = 0;
float cEnextMax = 0;
int aIDMax = 0;
int cIDMax = 0;
int aIDnextMax = 0;
int cIDnextMax = 0;
// Define the excluded SX3 and QQQ channels
// std::unordered_set<int> excludeSX3 = {34, 35, 36, 37, 61, 62, 67, 73, 74, 75, 76, 77, 78, 79, 80, 93, 97, 100, 103, 108, 109, 110, 111, 112};
// std::unordered_set<int> excludeQQQ = {0, 17, 109, 110, 111, 112, 113, 119, 127, 128};
// inCuth=false;
// inCutl=false;
// inPCCut=false;
for (int i = 0; i < pc.multi; i++)
{
if (pc.e[i] > 100 /*&& pc.multi < 7*/)
{
// creating a vector of pairs of anode and cathode hits
if (pc.index[i] < 24)
{
anodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
}
else if (pc.index[i] >= 24)
{
cathodeHits.push_back(std::pair<int, double>(pc.index[i] - 24, pc.e[i]));
}
for (int j = i + 1; j < pc.multi; j++)
{
// if(PCCoinc_cut1->IsInside(pc.index[i], pc.index[j]) || PCCoinc_cut2->IsInside(pc.index[i], pc.index[j])){
// // hpcCoin->Fill(pc.index[i], pc.index[j]);
// inPCCut = true;
// }
hpcCoin->Fill(pc.index[i], pc.index[j]);
}
}
}
// sorting the anode and cathode hits in descending order of energy
std::sort(anodeHits.begin(), anodeHits.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b)
{ return a.second > b.second; });
std::sort(cathodeHits.begin(), cathodeHits.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b)
{ return a.second > b.second; });
bool SiPCflag;
corrcatMax.clear();
if (anodeHits.size() >= 1 && cathodeHits.size() > 1)
{
if (((TMath::TanH(hitPos.Y() / hitPos.X())) > (TMath::TanH(a.Y() / a.X()) - TMath::PiOver4())) || ((TMath::TanH(hitPos.Y() / hitPos.X())) < (TMath::TanH(a.Y() / a.X()) + TMath::PiOver4())))
{
for (const auto &anode : anodeHits)
{
aID = anode.first;
aE = anode.second;
aESum += aE;
if (aE > aEMax)
{
aEMax = aE;
aIDMax = aID;
}
if (aE > aEnextMax && aE < aEMax)
{
aEnextMax = aE;
aIDnextMax = aID;
}
// for(const auto &cat : cathodeHits){
// hAVCcoin->Fill(aID, cat.first);
// }
}
// std::cout << " Anode iD : " << aIDMax << " Energy : " << aEMax << std::endl;
// printf("aID : %d, aE : %f, cE : %f\n", aID, aE, cE);
for (const auto &cathode : cathodeHits)
{
cID = cathode.first;
cE = cathode.second;
// std::cout << "Cathode ID : " << cID << " Energy : " << cE << std::endl;
hAVCcoin->Fill(aIDMax, cID);
// This section of code is used to find the cathodes are correlated with the max and next max anodes, as well as to figure out if there are any common cathodes
// the anodes are correlated with the cathodes +/-3 from the anode number in the reverse order
for (int j = -4; j < 3; j++)
{
if ((aIDMax + 24 + j) % 24 == 23 - cID)
/* the 23-cID is used to accomodate for the fact that the order of the cathodes was reversed relative top the physical geometry */
// if (Crossover[aIDMax][cID][0].z != 9999999)
{ {
corrcatMax.push_back(std::pair<int, double>(cID, cE)); chWedge = qqq.ch[i];
cESum += cE; eWedgeRaw = qqq.e[i];
// printf("Max Anode : %d Correlated Cathode : %d Anode Energy : %f z value : %f \n", aIDMax, cID, cESum, Crossover[aIDMax][cID][1].z /*prints alpha*/); eWedge = qqq.e[i] * qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
// std::cout << " Cathode iD : " << cID << " Energy : " << cE << std::endl; // printf("Wedge E: %.2f Gain: %.4f \n", eWedge, qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]);
chRing = qqq.ch[j] - 16;
eRingRaw = qqq.e[j];
eRing = qqq.e[j]; //* qqqrGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i]-16];
}
else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 /* && qqqrGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16] */ && qqqGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16])
{
chWedge = qqq.ch[j];
eWedge = qqq.e[j] * qqqGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
eWedgeRaw = qqq.e[j];
chRing = qqq.ch[i] - 16;
eRing = qqq.e[i]; // * qqqrGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
eRingRaw = qqq.e[i];
}
else
continue;
// plug in calibrations
if (qqqCalibValid[qqq.id[i]][chRing][chWedge])
{
eWedgeMeV = eWedge * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000;
eRingMeV = eRing * qqqCalib[qqq.id[i]][chRing][chWedge] / 1000;
}
else
continue;
// printf(" ID : %d , chWedge : %d, chRing : %d \n", qqq.id[i], chWedge, chRing);
double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
double rho = 50. + 40. / 16. * (chRing + 0.5);
// if(qqq.e[i]>50){
hqqqPolar->Fill(theta, rho);
// }
// qqq.used[i] = true;
// qqq.used[j] = true;
if (!HitNonZero)
{
double x = rho * TMath::Cos(theta);
double y = rho * TMath::Sin(theta);
hitPos.SetXYZ(x, y, 23 + 75 + 30);
HitNonZero = true;
} }
} }
} }
} }
} }
// //======================= PC
TVector3 anodeIntersection; // Calculate the crossover points and put them into an array
anodeIntersection.Clear();
// Implementing a method for PC reconstruction using a single Anode event pwinstance.ConstructGeo();
// if (anodeHits.size() == 1) Coord Crossover[24][24][2];
{ TVector3 a, c, diff;
float x, y, z = 0; double a2, ac, c2, adiff, cdiff, denom, alpha, beta;
for (const auto &corr : corrcatMax) int index = 0;
for (int i = 0; i < pwinstance.An.size(); i++)
{ {
if (cESum > 0) a = pwinstance.An[i].first - pwinstance.An[i].second;
{
x += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].x;
y += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].y;
z += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].z;
// printf("Max Anode : %d Correlated Cathode : %d cathode Energy : %f cESum Energy : %f z value : %f \n", aIDMax, corr.first, corr.second, cESum, Crossover[aIDMax][corr.first][1].z /*prints alpha*/);
}
else
{
printf("Warning: No valid cathode hits to correlate with anode %d! \n", aIDMax);
}
// printf("EventID : %llu, Max Anode : %d Cathode: %d PC X and Y : (%f, %f) \n", entry, aIDMax, cID, Crossover[aIDMax][cID][0].x, Crossover[aIDMax][cID][0].y);
// for (int i = 0; i < sx3.multi; i++)
// {
// printf("EventID : %llu, HitPos X, Y, Z: %f %f %f SX3ID : %d %d \n", entry, hitPos.X(), hitPos.Y(), hitPos.Z(), sx3.id[i], sx3.ch[i]);
// }
// for (int i = 0; i < qqq.multi; i++) for (int j = 0; j < pwinstance.Ca.size(); j++)
// { {
// printf("Max Anode : %d Cathode: %d PC X and Y : %f %f \n", aIDMax, cID, Crossover[aIDMax][cID][0].x, Crossover[aIDMax][cID][0].y); // Ok so this method uses what is essentially the solution of 2 equations to find the point of intersection between the anode and cathode wires
// printf("HitPos X, Y, Z, QQQID : %f %f %f %d \n", hitPos.X(), hitPos.Y(), hitPos.Z(), qqq.id[i]); // here a and c are the vectors of the anode and cathode wires respectively
// } // diff is the perpendicular vector between the anode and cathode wires
// The idea behind this is to then find the scalars alpha and beta that give a ratio between 0 and -1,
c = pwinstance.Ca[j].first - pwinstance.Ca[j].second;
diff = pwinstance.An[i].first - pwinstance.Ca[j].first;
a2 = a.Dot(a);
c2 = c.Dot(c);
ac = a.Dot(c);
adiff = a.Dot(diff);
cdiff = c.Dot(diff);
denom = a2 * c2 - ac * ac;
alpha = (ac * cdiff - c2 * adiff) / denom;
beta = (a2 * cdiff - ac * adiff) / denom;
Crossover[i][j][0].x = pwinstance.An[i].first.X() + alpha * a.X();
Crossover[i][j][0].y = pwinstance.An[i].first.Y() + alpha * a.Y();
Crossover[i][j][0].z = pwinstance.An[i].first.Z() + alpha * a.Z();
if (Crossover[i][j][0].z < -190 || Crossover[i][j][0].z > 190)
{
Crossover[i][j][0].z = 9999999;
}
// placeholder variable Crossover[i][j][1].x has nothing to do with the geometry of the crossover and is being used to store the alpha value-
//-so that it can be used to sort "good" hits later
Crossover[i][j][1].x = alpha;
Crossover[i][j][1].y = 0;
// if(i==0){
// printf("CID, Crossover z and alpha are : %d %f %f \n", j, Crossover[i][j][0].z, Crossover[i][j][1].x /*this is alpha*/);
// }
// }
// }
}
} }
anodeIntersection = TVector3(x, y, z); // printf("Anode and cathode indices, alpha, denom, andiff, cndiff : %d %d %f %f %f %f\n", i, j, alpha, denom, adiff, cdiff);
// std::cout << "Anode Intersection " << anodeIntersection.Z() << " " << x << " " << y << " " << z << std::endl;
// anodeIntersection.Clear();
for (int i = 0; i < pc.multi; i++)
{
if (pc.e[i] > 100)
{
hpcIndexVE->Fill(pc.index[i], pc.e[i]); // non gain matched energy
}
// Gain Matching of PC wires
if (pc.index[i] >= 0 && pc.index[i] < 48)
{
// printf("index: %d, Old cathode energy: %d \n", pc.index[i],pc.e[i]);
auto it = slopeInterceptMap.find(pc.index[i]);
if (it != slopeInterceptMap.end())
{
double slope = it->second.first;
double intercept = it->second.second;
// printf("slope: %f, intercept:%f\n" ,slope, intercept);
pc.e[i] = slope * pc.e[i] + intercept;
// printf("index: %d, New cathode energy: %d \n",pc.index[i], pc.e[i]);
}
hpcIndexVE_GM->Fill(pc.index[i], pc.e[i]);
// hPC_E[pc.index[i]]->Fill(pc.e[i]); // gain matched energy per channel
}
}
std::vector<std::pair<int, double>> anodeHits = {};
std::vector<std::pair<int, double>> cathodeHits = {};
std::vector<std::pair<int, double>> corrcatMax = {};
std::vector<std::pair<int, double>> corrcatnextMax = {};
std::vector<std::pair<int, double>> commcat = {};
int aID = 0;
int cID = 0;
float aE = 0;
float cE = 0;
float aESum = 0;
float cESum = 0;
float aEMax = 0;
float cEMax = 0;
float aEnextMax = 0;
float cEnextMax = 0;
int aIDMax = 0;
int cIDMax = 0;
int aIDnextMax = 0;
int cIDnextMax = 0;
// Define the excluded SX3 and QQQ channels
// std::unordered_set<int> excludeSX3 = {34, 35, 36, 37, 61, 62, 67, 73, 74, 75, 76, 77, 78, 79, 80, 93, 97, 100, 103, 108, 109, 110, 111, 112};
// std::unordered_set<int> excludeQQQ = {0, 17, 109, 110, 111, 112, 113, 119, 127, 128};
// inCuth=false;
// inCutl=false;
// inPCCut=false;
for (int i = 0; i < pc.multi; i++)
{
if (pc.e[i] > 100 /*&& pc.multi < 7*/)
{
// creating a vector of pairs of anode and cathode hits
if (pc.index[i] < 24)
{
anodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
}
else if (pc.index[i] >= 24)
{
cathodeHits.push_back(std::pair<int, double>(pc.index[i] - 24, pc.e[i]));
}
for (int j = i + 1; j < pc.multi; j++)
{
// if(PCCoinc_cut1->IsInside(pc.index[i], pc.index[j]) || PCCoinc_cut2->IsInside(pc.index[i], pc.index[j])){
// // hpcCoin->Fill(pc.index[i], pc.index[j]);
// inPCCut = true;
// }
hpcCoin->Fill(pc.index[i], pc.index[j]);
}
}
}
// sorting the anode and cathode hits in descending order of energy
std::sort(anodeHits.begin(), anodeHits.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b)
{ return a.second > b.second; });
std::sort(cathodeHits.begin(), cathodeHits.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b)
{ return a.second > b.second; });
bool SiPCflag;
corrcatMax.clear();
if (anodeHits.size() >= 1 && cathodeHits.size() > 1)
{
if (((TMath::TanH(hitPos.Y() / hitPos.X())) > (TMath::TanH(a.Y() / a.X()) - TMath::PiOver4())) || ((TMath::TanH(hitPos.Y() / hitPos.X())) < (TMath::TanH(a.Y() / a.X()) + TMath::PiOver4())))
{
for (const auto &anode : anodeHits)
{
aID = anode.first;
aE = anode.second;
aESum += aE;
if (aE > aEMax)
{
aEMax = aE;
aIDMax = aID;
}
if (aE > aEnextMax && aE < aEMax)
{
aEnextMax = aE;
aIDnextMax = aID;
}
// for(const auto &cat : cathodeHits){
// hAVCcoin->Fill(aID, cat.first);
// }
}
// std::cout << " Anode iD : " << aIDMax << " Energy : " << aEMax << std::endl;
// printf("aID : %d, aE : %f, cE : %f\n", aID, aE, cE);
for (const auto &cathode : cathodeHits)
{
cID = cathode.first;
cE = cathode.second;
// std::cout << "Cathode ID : " << cID << " Energy : " << cE << std::endl;
hAVCcoin->Fill(aIDMax, cID);
// This section of code is used to find the cathodes are correlated with the max and next max anodes, as well as to figure out if there are any common cathodes
// the anodes are correlated with the cathodes +/-3 from the anode number in the reverse order
for (int j = -4; j < 3; j++)
{
if ((aIDMax + 24 + j) % 24 == 23 - cID)
/* the 23-cID is used to accomodate for the fact that the order of the cathodes was reversed relative top the physical geometry */
// if (Crossover[aIDMax][cID][0].z != 9999999)
{
corrcatMax.push_back(std::pair<int, double>(cID, cE));
cESum += cE;
// printf("Max Anode : %d Correlated Cathode : %d Anode Energy : %f z value : %f \n", aIDMax, cID, cESum, Crossover[aIDMax][cID][1].z /*prints alpha*/);
// std::cout << " Cathode iD : " << cID << " Energy : " << cE << std::endl;
}
}
}
}
}
TVector3 anodeIntersection;
anodeIntersection.Clear();
// Implementing a method for PC reconstruction using a single Anode event
// if (anodeHits.size() == 1)
{
float x, y, z = 0;
for (const auto &corr : corrcatMax)
{
if (cESum > 0)
{
x += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].x;
y += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].y;
z += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].z;
// printf("Max Anode : %d Correlated Cathode : %d cathode Energy : %f cESum Energy : %f z value : %f \n", aIDMax, corr.first, corr.second, cESum, Crossover[aIDMax][corr.first][1].z /*prints alpha*/);
}
else
{
printf("Warning: No valid cathode hits to correlate with anode %d! \n", aIDMax);
}
// printf("EventID : %llu, Max Anode : %d Cathode: %d PC X and Y : (%f, %f) \n", entry, aIDMax, cID, Crossover[aIDMax][cID][0].x, Crossover[aIDMax][cID][0].y);
// for (int i = 0; i < sx3.multi; i++)
// {
// printf("EventID : %llu, HitPos X, Y, Z: %f %f %f SX3ID : %d %d \n", entry, hitPos.X(), hitPos.Y(), hitPos.Z(), sx3.id[i], sx3.ch[i]);
// }
// for (int i = 0; i < qqq.multi; i++)
// {
// printf("Max Anode : %d Cathode: %d PC X and Y : %f %f \n", aIDMax, cID, Crossover[aIDMax][cID][0].x, Crossover[aIDMax][cID][0].y);
// printf("HitPos X, Y, Z, QQQID : %f %f %f %d \n", hitPos.X(), hitPos.Y(), hitPos.Z(), qqq.id[i]);
// }
}
anodeIntersection = TVector3(x, y, z);
// std::cout << "Anode Intersection " << anodeIntersection.Z() << " " << x << " " << y << " " << z << std::endl;
}
if (anodeIntersection.Z() != 0)
{
hPCZProj->Fill(anodeIntersection.Z());
}
// Filling the PC Z projection histogram
// std::cout << anodeIntersection.Z() << std::endl;
// hPCZProj->Fill(anodeIntersection.Z());
// }
// inCuth = false;
// inCutl = false;
// inPCCut = false;
// for(int j=i+1;j<pc.multi;j++){
// if(PCCoinc_cut1->IsInside(pc.index[i], pc.index[j]) || PCCoinc_cut2->IsInside(pc.index[i], pc.index[j])){
// // hpcCoin->Fill(pc.index[i], pc.index[j]);
// inPCCut = true;
// }
// hpcCoin->Fill(pc.index[i], pc.index[j]);
// }
// Check if the accumulated energies are within the defined ranges
// if (AnCatSum_high && AnCatSum_high->IsInside(aESum, cESum)) {
// inCuth = true;
// }
// if (AnCatSum_low && AnCatSum_low->IsInside(aESum, cESum)) {
// inCutl = true;
// }
// Fill histograms based on the cut conditions
// if (inCuth && inPCCut) {
// hanVScatsum_hcut->Fill(aESum, cESum);
// }
// if (inCutl && inPCCut) {
// hanVScatsum_lcut->Fill(aESum, cESum);
// }
// for(auto anode : anodeHits){
// float aE = anode.second;
// aESum += aE;
// if(inPCCut){
hanVScatsum->Fill(aEMax, cESum);
// }
// if (sx3ecut)
// {
hCat4An->Fill(corrcatMax.size());
hNosvAe->Fill(corrcatMax.size(), aEMax);
hAnodehits->Fill(anodeHits.size());
// }
// }
if (anodeHits.size() < 1)
{
hCat0An->Fill(cathodeHits.size());
}
if (HitNonZero && anodeIntersection.Z() != 0)
{
pw_contr.CalTrack2(hitPos, anodeIntersection);
hZProj->Fill(pw_contr.GetZ0());
}
// ########################################################### Track constrcution
// ############################## DO THE KINEMATICS
return kTRUE;
} }
if (anodeIntersection.Z() != 0) void Analyzer::Terminate()
{ {
hPCZProj->Fill(anodeIntersection.Z());
// gStyle->SetOptStat("neiou");
// TCanvas *canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000);
// canvas->Divide(3, 3);
// // hsx3VpcIndex->Draw("colz");
// //=============================================== pad-1
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// hsx3IndexVE->Draw("colz");
// //=============================================== pad-2
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// hqqqIndexVE->Draw("colz");
// //=============================================== pad-3
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// hpcIndexVE->Draw("colz");
// //=============================================== pad-4
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// hsx3Coin->Draw("colz");
// //=============================================== pad-5
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// canvas->cd(padID)->SetLogz(true);
// hqqqCoin->Draw("colz");
// //=============================================== pad-6
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// hpcCoin->Draw("colz");
// //=============================================== pad-7
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// // hsx3VpcIndex ->Draw("colz");
// hsx3VpcE->Draw("colz");
// //=============================================== pad-8
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// // hqqqVpcIndex ->Draw("colz");
// hqqqVpcE->Draw("colz");
// //=============================================== pad-9
// padID++;
// // canvas->cd(padID)->DrawFrame(-50, -50, 50, 50);
// // hqqqPolar->Draw("same colz pol");
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// // hZProj->Draw();
// hanVScatsum->Draw("colz");
// // TFile *outRoot = new TFile("Histograms.root", "RECREATE");
// // if (!outRoot->IsOpen())
// // {
// // std::cerr << "Error opening file for writing!" << std::endl;
// // return;
// // }
// // // Loop through histograms and write them to the ROOT file
// // for (int i = 0; i < 48; i++)
// // {
// // if (hPC_E[i] != nullptr)
// // {
// // hPC_E[i]->Write(); // Write histogram to file
// // }
// // }
// // outRoot->Close();
} }
// Filling the PC Z projection histogram
// std::cout << anodeIntersection.Z() << std::endl;
// hPCZProj->Fill(anodeIntersection.Z());
// }
// inCuth = false;
// inCutl = false;
// inPCCut = false;
// for(int j=i+1;j<pc.multi;j++){
// if(PCCoinc_cut1->IsInside(pc.index[i], pc.index[j]) || PCCoinc_cut2->IsInside(pc.index[i], pc.index[j])){
// // hpcCoin->Fill(pc.index[i], pc.index[j]);
// inPCCut = true;
// }
// hpcCoin->Fill(pc.index[i], pc.index[j]);
// }
// Check if the accumulated energies are within the defined ranges
// if (AnCatSum_high && AnCatSum_high->IsInside(aESum, cESum)) {
// inCuth = true;
// }
// if (AnCatSum_low && AnCatSum_low->IsInside(aESum, cESum)) {
// inCutl = true;
// }
// Fill histograms based on the cut conditions
// if (inCuth && inPCCut) {
// hanVScatsum_hcut->Fill(aESum, cESum);
// }
// if (inCutl && inPCCut) {
// hanVScatsum_lcut->Fill(aESum, cESum);
// }
// for(auto anode : anodeHits){
// float aE = anode.second;
// aESum += aE;
// if(inPCCut){
hanVScatsum->Fill(aEMax, cESum);
// }
// if (sx3ecut)
// {
hCat4An->Fill(corrcatMax.size());
hNosvAe->Fill(corrcatMax.size(), aEMax);
hAnodehits->Fill(anodeHits.size());
// }
// }
if (anodeHits.size() < 1)
{
hCat0An->Fill(cathodeHits.size());
}
if (HitNonZero && anodeIntersection.Z() != 0)
{
pw_contr.CalTrack2(hitPos, anodeIntersection);
hZProj->Fill(pw_contr.GetZ0());
}
// ########################################################### Track constrcution
// ############################## DO THE KINEMATICS
return kTRUE;
}
void Analyzer::Terminate()
{
// gStyle->SetOptStat("neiou");
// TCanvas *canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000);
// canvas->Divide(3, 3);
// // hsx3VpcIndex->Draw("colz");
// //=============================================== pad-1
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// hsx3IndexVE->Draw("colz");
// //=============================================== pad-2
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// hqqqIndexVE->Draw("colz");
// //=============================================== pad-3
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// hpcIndexVE->Draw("colz");
// //=============================================== pad-4
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// hsx3Coin->Draw("colz");
// //=============================================== pad-5
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// canvas->cd(padID)->SetLogz(true);
// hqqqCoin->Draw("colz");
// //=============================================== pad-6
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// hpcCoin->Draw("colz");
// //=============================================== pad-7
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// // hsx3VpcIndex ->Draw("colz");
// hsx3VpcE->Draw("colz");
// //=============================================== pad-8
// padID++;
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// // hqqqVpcIndex ->Draw("colz");
// hqqqVpcE->Draw("colz");
// //=============================================== pad-9
// padID++;
// // canvas->cd(padID)->DrawFrame(-50, -50, 50, 50);
// // hqqqPolar->Draw("same colz pol");
// canvas->cd(padID);
// canvas->cd(padID)->SetGrid(1);
// // hZProj->Draw();
// hanVScatsum->Draw("colz");
// // TFile *outRoot = new TFile("Histograms.root", "RECREATE");
// // if (!outRoot->IsOpen())
// // {
// // std::cerr << "Error opening file for writing!" << std::endl;
// // return;
// // }
// // // Loop through histograms and write them to the ROOT file
// // for (int i = 0; i < 48; i++)
// // {
// // if (hPC_E[i] != nullptr)
// // {
// // hPC_E[i]->Write(); // Write histogram to file
// // }
// // }
// // outRoot->Close();
}

26
Calibration.C
View File

@ -28,9 +28,9 @@ const int MAX_QQQ = 4;
const int MAX_RING = 16; const int MAX_RING = 16;
const int MAX_WEDGE = 16; const int MAX_WEDGE = 16;
double qqqwGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}}; double qqqwGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
double qqqrGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}}; // double qqqrGain[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{0}}};
bool qqqwGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}}; bool qqqwGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
bool qqqrGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}}; // bool qqqrGainValid[MAX_QQQ][MAX_RING][MAX_WEDGE] = {{{false}}};
void Calibration::Begin(TTree * /*tree*/) void Calibration::Begin(TTree * /*tree*/)
{ {
@ -50,9 +50,9 @@ void Calibration::Begin(TTree * /*tree*/)
while (infile >> det >> ring >> wedge >> gainw >> gainr) while (infile >> det >> ring >> wedge >> gainw >> gainr)
{ {
qqqwGain[det][ring][wedge] = gainw; qqqwGain[det][ring][wedge] = gainw;
qqqrGain[det][ring][wedge] = gainr; // qqqrGain[det][ring][wedge] = gainr;
qqqwGainValid[det][ring][wedge] = (gainw > 0); qqqwGainValid[det][ring][wedge] = (gainw > 0);
qqqrGainValid[det][ring][wedge] = (gainr > 0); // qqqrGainValid[det][ring][wedge] = (gainr > 0);
} }
infile.close(); infile.close();
std::cout << "Loaded QQQ gains from " << filename << std::endl; std::cout << "Loaded QQQ gains from " << filename << std::endl;
@ -96,7 +96,7 @@ Bool_t Calibration::Process(Long64_t entry)
float eWedge = 0.0; float eWedge = 0.0;
float eRingRaw = 0.0; float eRingRaw = 0.0;
float eRing = 0.0; float eRing = 0.0;
if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqrGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16] && qqqwGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]) if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && /*qqqrGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16] &&*/ qqqwGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16])
{ {
chWedge = qqq.ch[i]; chWedge = qqq.ch[i];
eWedgeRaw = qqq.e[i]; eWedgeRaw = qqq.e[i];
@ -104,16 +104,16 @@ Bool_t Calibration::Process(Long64_t entry)
// printf("Wedge E: %.2f Gain: %.4f \n", eWedge, qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]); // printf("Wedge E: %.2f Gain: %.4f \n", eWedge, qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]);
chRing = qqq.ch[j] - 16; chRing = qqq.ch[j] - 16;
eRingRaw = qqq.e[j]; eRingRaw = qqq.e[j];
eRing = qqq.e[j] * qqqrGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16]; eRing = qqq.e[j];// * qqqrGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
} }
else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqrGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16] && qqqwGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16]) else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && /*qqqrGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16] &&*/ qqqwGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16])
{ {
chWedge = qqq.ch[j]; chWedge = qqq.ch[j];
eWedge = qqq.e[j] * qqqwGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16]; eWedge = qqq.e[j] * qqqwGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
eWedgeRaw = qqq.e[j]; eWedgeRaw = qqq.e[j];
chRing = qqq.ch[i] - 16; chRing = qqq.ch[i] - 16;
eRing = qqq.e[i] * qqqrGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]; eRing = qqq.e[i];// * qqqrGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
eRingRaw = qqq.e[i]; eRingRaw = qqq.e[i];
} }
else else
@ -137,8 +137,14 @@ Bool_t Calibration::Process(Long64_t entry)
const double MIN_ADC = 1500.0; const double MIN_ADC = 1500.0;
const double MAX_ADC = 3000.0; const double MAX_ADC = 3000.0;
if (eWedge >= MIN_ADC && eWedge <= MAX_ADC && // if (eWedge >= MIN_ADC && eWedge <= MAX_ADC &&
eRing >= MIN_ADC && eRing <= MAX_ADC) // eRing >= MIN_ADC && eRing <= MAX_ADC)
double ratio = (eWedge > 0.0) ? (eRing / eWedge) : 0.0;
double maxslope = 1.5;
bool validPoint = false;
if (ratio < maxslope && ratio > 1. / maxslope)
{ {
// Accumulate data for gain matching // Accumulate data for gain matching
dataPoints[{qqq.id[i], chRing, chWedge}].emplace_back(eWedge, eRing); dataPoints[{qqq.id[i], chRing, chWedge}].emplace_back(eWedge, eRing);

278
GainMatchQQQ.C
View File

@ -14,6 +14,9 @@
#include <numeric> #include <numeric>
#include "Armory/HistPlotter.h" #include "Armory/HistPlotter.h"
#include "TVector3.h" #include "TVector3.h"
#include "TGraphErrors.h"
#include "TF1.h"
#include <cmath>
TH2F *hQQQFVB; TH2F *hQQQFVB;
HistPlotter *plotter; HistPlotter *plotter;
@ -123,11 +126,11 @@ Bool_t GainMatchQQQ::Process(Long64_t entry)
bool validPoint = false; bool validPoint = false;
if (ratio < maxslope && ratio > 1. / maxslope) if (ratio < maxslope && ratio > 1. / maxslope)
{ {
// Accumulate data for gain matching // Accumulate data for gain matching
dataPoints[{id, chr, chw}].emplace_back(ew, er); dataPoints[{id, chr, chw}].emplace_back(ew, er);
plotter->Fill2D("hAll_in", 4000, 0, 16000, 4000, 0, 16000, ew, er); plotter->Fill2D("hAll_in", 4000, 0, 16000, 4000, 0, 16000, ew, er);
validPoint = true; validPoint = true;
} }
if (!validPoint) if (!validPoint)
{ {
plotter->Fill2D("hAll_out", 4000, 0, 16000, 4000, 0, 16000, ew, er); plotter->Fill2D("hAll_out", 4000, 0, 16000, 4000, 0, 16000, ew, er);
@ -137,223 +140,170 @@ Bool_t GainMatchQQQ::Process(Long64_t entry)
return kTRUE; return kTRUE;
} }
void GainMatchQQQ::Terminate() void GainMatchQQQ::Terminate()
{ {
const int MAX_DET = 4; const int MAX_DET = 4;
const int MAX_RING = 16; const int MAX_RING = 16;
const int MAX_WEDGE = 16; const int MAX_WEDGE = 16;
// We store gains locally just for the "corrected" plot,
// but the file will output Slopes for the global minimizer.
double gainW[MAX_DET][MAX_RING][MAX_WEDGE] = {{{0}}}; double gainW[MAX_DET][MAX_RING][MAX_WEDGE] = {{{0}}};
double gainR[MAX_DET][MAX_RING][MAX_WEDGE] = {{{0}}}; double gainR[MAX_DET][MAX_RING][MAX_WEDGE] = {{{0}}};
bool gainValid[MAX_DET][MAX_RING][MAX_WEDGE] = {{{false}}}; bool gainValid[MAX_DET][MAX_RING][MAX_WEDGE] = {{{false}}};
// Output file for the Minimizer
std::ofstream outFile("qqq_GainMatch.txt"); std::ofstream outFile("qqq_GainMatch.txt");
if (!outFile.is_open())
{
std::cerr << "Error opening output file!" << std::endl;
return;
}
// Parameters for sigma-clipping // Benchmark/Debug file
const int MIN_POINTS = 5; std::ofstream benchFile("benchmark_diff.txt");
const int MAX_ITER = 5; benchFile << "ID Wedge Ring Chi2NDF Slope SlopeErr" << std::endl;
const double CLIP_SIGMA = 3.0;
if (!outFile.is_open()) { std::cerr << "Error opening output file!" << std::endl; return; }
const int MIN_POINTS = 50;
const int MAX_ITER = 3; // Outlier rejection passes
const double CLIP_SIGMA = 2.5; // Sigma threshold for outliers
for (const auto &kv : dataPoints) for (const auto &kv : dataPoints)
{ {
auto key = kv.first; auto key = kv.first;
auto [id, ring, wedge] = key; auto [id, ring, wedge] = key;
const auto &pts_in = kv.second; const auto &pts = kv.second;
if (pts_in.size() < (size_t)MIN_POINTS)
continue;
// Make a working copy of the points for clipping if (pts.size() < (size_t)MIN_POINTS) continue;
std::vector<std::pair<double, double>> pts = pts_in;
bool solved = false; std::vector<std::pair<double, double>> current_pts = pts;
double final_gW = 0.0;
double final_gR = 0.0;
// make a vector of ring wedge as I go through all points and fit a tgraph on the fly
// Iterative sigma-clipping loop double finalSlope = 0.0;
double finalSlopeErr = 0.0;
bool converged = false;
// --- Iterative Fitting ---
for (int iter = 0; iter < MAX_ITER; ++iter) for (int iter = 0; iter < MAX_ITER; ++iter)
{ {
// Compute sums if (current_pts.size() < (size_t)MIN_POINTS) break;
double sum_w2 = 0.0;
double sum_wr = 0.0; std::vector<double> x, y, ex, ey;
double sum_r2 = 0.0;
for (const auto &p : pts) for (const auto &p : current_pts)
{ {
double w = p.first; x.push_back(p.first); // Wedge E
double r = p.second; y.push_back(p.second); // Ring E
sum_w2 += w * w; ex.push_back(std::sqrt(std::abs(p.first))); // Error in X (Poisson)
sum_wr += w * r; ey.push_back(std::sqrt(std::abs(p.second))); // Error in Y (Poisson)
sum_r2 += r * r;
// Sanity check to avoid 0 error
if(ex.back() < 1.0) ex.back() = 1.0;
if(ey.back() < 1.0) ey.back() = 1.0;
} }
// Guard against degenerate cases // 2. Create Graph
if (sum_w2 <= 0.0 || sum_wr <= 0.0) TGraphErrors *gr = new TGraphErrors(current_pts.size(), x.data(), y.data(), ex.data(), ey.data());
{
// // fallback to single-parameter linear fit (original method) // 3. Fit Linear Function through Origin
// // Use ROOT TGraph fitting as fallback TF1 *f1= new TF1("calibFit", "[0]*x", 0, 16000);
// if (pts.size() >= 2) f1->SetParameter(0, 1.0);
// {
// std::vector<double> wE, rE; // "Q"=Quiet, "N"=NoDraw, "S"=ResultPtr
// wE.reserve(pts.size()); // We do NOT use "W" (Ignore weights), we want to use the errors we set.
// rE.reserve(pts.size()); int fitStatus = gr->Fit(f1, "QNS");
// for (const auto &pr : pts)
// { if (fitStatus != 0) {
// wE.push_back(pr.first); delete gr; delete f1;
// rE.push_back(pr.second);
// }
// TGraph g(static_cast<int>(wE.size()), wE.data(), rE.data());
// TF1 f("f_fallback", "[0]*x", 0, 16000);
// g.Fit(&f, "QNR");
// double slope = f.GetParameter(0);
// if (slope > 0)
// {
// // distribute correction symmetrically:
// double alpha = slope; // r ≈ slope * w => alpha = slope
// double gW_try = std::sqrt(alpha);
// double gR_try = 1.0 / gW_try;
// final_gW = gW_try;
// final_gR = gR_try;
// solved = true;
// }
// }
break; break;
} }
// alpha = sum(w*r) / sum(w^2) finalSlope = f1->GetParameter(0);
double alpha = sum_wr / sum_w2; double chi2 = f1->GetChisquare();
double ndf = f1->GetNDF();
if (!(alpha > 0.0) || !std::isfinite(alpha)) // Get the statistical error on the slope
{ double rawErr = f1->GetParError(0);
// // degenerate; fallback to TF1 fit as above
// if (pts.size() >= 2) // SCALING ERROR:
// { // If Chi2/NDF > 1, the data scatters more than Poisson stats predict.
// std::vector<double> wE, rE; // // We inflate the error by sqrt(Chi2/NDF) to be conservative for the Minimizer.
// wE.reserve(pts.size()); // double redChi2 = (ndf > 0) ? (chi2 / ndf) : 1.0;
// rE.reserve(pts.size()); // double inflation = (redChi2 > 1.0) ? std::sqrt(redChi2) : 1.0;
// for (const auto &pr : pts)
// { // finalSlopeErr = rawErr * inflation;
// wE.push_back(pr.first);
// rE.push_back(pr.second); // 4. Outlier Rejection
// } if (iter == MAX_ITER - 1) {
// TGraph g(static_cast<int>(wE.size()), wE.data(), rE.data()); converged = true;
// TF1 f("f_fallback2", "[0]*x", 0, 16000); delete gr; delete f1;
// g.Fit(&f, "QNR");
// double slope = f.GetParameter(0);
// if (slope > 0)
// {
// double gW_try = std::sqrt(slope);
// double gR_try = 1.0 / gW_try;
// final_gW = gW_try;
// final_gR = gR_try;
// solved = true;
// }
// }
break; break;
} }
// distribute correction between W and R // Calculate Residuals
double gW = std::sqrt(alpha);
double gR = 1.0 / gW;
// compute residuals and sigma
std::vector<double> residuals; std::vector<double> residuals;
residuals.reserve(pts.size()); double sumSqResid = 0.0;
for (const auto &p : pts) for(size_t k=0; k<current_pts.size(); ++k) {
{ double val = f1->Eval(current_pts[k].first);
double w = p.first; double res = current_pts[k].second - val;
double r = p.second;
double res = gW * w - gR * r;
residuals.push_back(res); residuals.push_back(res);
sumSqResid += res*res;
} }
// double sigma = std::sqrt(sumSqResid / current_pts.size());
// mean and stddev (use population sigma here) // // Filter
double mean = 0.0; // std::vector<std::pair<double, double>> next_pts;
for (double v : residuals) // for(size_t k=0; k<current_pts.size(); ++k) {
mean += v; // if(std::abs(residuals[k]) < CLIP_SIGMA * sigma) {
mean /= residuals.size(); // next_pts.push_back(current_pts[k]);
// }
// }
double var = 0.0; // if (next_pts.size() == current_pts.size()) {
for (double v : residuals) // converged = true;
var += (v - mean) * (v - mean); // delete gr; delete f1;
var /= residuals.size(); // break;
double sigma = std::sqrt(var); // }
// current_pts = next_pts;
// delete gr; delete f1;
}
// If sigma is NaN or zero, accept and break if (!converged || finalSlope <= 0) continue;
if (!std::isfinite(sigma) || sigma == 0.0)
{
final_gW = gW;
final_gR = gR;
solved = true;
break;
}
// Clip > CLIP_SIGMA and build new pts // --- Store/Output ---
size_t before = pts.size();
std::vector<std::pair<double, double>> new_pts;
new_pts.reserve(pts.size());
for (size_t k = 0; k < pts.size(); ++k)
{
if (std::fabs(residuals[k] - mean) <= CLIP_SIGMA * sigma)
new_pts.push_back(pts[k]);
}
pts.swap(new_pts);
size_t after = pts.size();
// If no points removed or too few remain, accept current solution // 1. Save locally for the verification plot (hAll)
final_gW = gW; // Approximate local gain for plotting purposes only
final_gR = gR; double gW_local = std::sqrt(finalSlope);
solved = true; double gR_local = 1.0 / gW_local;
if (before == after || pts.size() < (size_t)MIN_POINTS) gainW[id][ring][wedge] = gW_local;
break; gainR[id][ring][wedge] = gR_local;
// otherwise iterate again with clipped pts
} // end iter loop
if (!solved)
continue;
// sanity checks: avoid ridiculous gains
if (!(final_gW > 0.0) || !(final_gR > 0.0) || !std::isfinite(final_gW) || !std::isfinite(final_gR))
continue;
// store gains
gainW[id][ring][wedge] = final_gW;
gainR[id][ring][wedge] = final_gR;
gainValid[id][ring][wedge] = true; gainValid[id][ring][wedge] = true;
// write out both gains: id wedge ring gW gR // 2. Write to File for Minimizer
outFile << id << " " << wedge << " " << ring << " " << final_gW << " " << final_gR << std::endl; // Format: ID Wedge Ring Slope Error
printf("Gain match Det%d Ring%d Wedge%d → gW=%.6f gR=%.6f\n", id, ring, wedge, final_gW, final_gR); outFile << id << " " << wedge << " " << ring << " " << finalSlope << " " << finalSlopeErr << std::endl;
// 3. Benchmark Info
benchFile << id << " " << wedge << " " << ring << " "
<< finalSlope << " " << finalSlopeErr << std::endl;
} }
outFile.close(); outFile.close();
benchFile.close();
std::cout << "Gain matching with Errors complete." << std::endl;
std::cout << "Gain matching complete." << std::endl; // Plotting the corrected data (Visual check using local approx gains)
// === Plot all gain-matched QQQ points together with a 2D histogram ===
// Fill the combined TH2F with corrected data
for (auto &kv : dataPoints) for (auto &kv : dataPoints)
{ {
int id, ring, wedge; int id, ring, wedge;
std::tie(id, ring, wedge) = kv.first; std::tie(id, ring, wedge) = kv.first;
if (!gainValid[id][ring][wedge]) if (!gainValid[id][ring][wedge]) continue;
continue;
auto &pts = kv.second; auto &pts = kv.second;
for (auto &pr : pts) for (auto &pr : pts)
{ {
double corrWedge = pr.first * gainW[id][ring][wedge]; double corrWedge = pr.first * gainW[id][ring][wedge];
double corrRing = pr.second * gainR[id][ring][wedge]; double corrRing = pr.second * gainR[id][ring][wedge];
// hAll->Fill(corrWedge, corrRing);
plotter->Fill2D("hAll", 4000, 0, 16000, 4000, 0, 16000, corrWedge, corrRing); plotter->Fill2D("hAll", 4000, 0, 16000, 4000, 0, 16000, corrWedge, corrRing);
plotter->Fill2D(Form("hGMQQQ_id%d_ring%d_wedge%d", id, ring, wedge), 400, 0, 16000, 400, 0, 16000, corrWedge, corrRing, "GainMatched");
} }
} }
// Optionally keep previous global histos too
plotter->FlushToDisk(); plotter->FlushToDisk();
} }

14
QQQ_Calcheck.C
View File

@ -52,9 +52,9 @@ void QQQ_Calcheck::Begin(TTree * /*tree*/)
while (infile >> det >> ring >> wedge >> gainw>> gainr) while (infile >> det >> ring >> wedge >> gainw>> gainr)
{ {
qqqwGain[det][ring][wedge] = gainw; qqqwGain[det][ring][wedge] = gainw;
qqqrGain[det][ring][wedge] = gainr; // qqqrGain[det][ring][wedge] = gainr;
qqqwGainValid[det][ring][wedge] = (gainw > 0); qqqwGainValid[det][ring][wedge] = (gainw > 0);
qqqrGainValid[det][ring][wedge] = (gainr > 0); // qqqrGainValid[det][ring][wedge] = (gainr > 0);
} }
infile.close(); infile.close();
std::cout << "Loaded QQQ gains from " << filename << std::endl; std::cout << "Loaded QQQ gains from " << filename << std::endl;
@ -110,7 +110,7 @@ Bool_t QQQ_Calcheck::Process(Long64_t entry)
float eRing = 0.0; float eRing = 0.0;
float eRingMeV = 0.0; float eRingMeV = 0.0;
// plug in gains // plug in gains
if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && qqqrGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16] && qqqwGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]) if (qqq.ch[i] < 16 && qqq.ch[j] >= 16 && /*qqqrGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16] &&*/ qqqwGainValid[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16])
{ {
chWedge = qqq.ch[i]; chWedge = qqq.ch[i];
eWedgeRaw = qqq.e[i]; eWedgeRaw = qqq.e[i];
@ -118,16 +118,16 @@ Bool_t QQQ_Calcheck::Process(Long64_t entry)
// printf("Wedge E: %.2f Gain: %.4f \n", eWedge, qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]); // printf("Wedge E: %.2f Gain: %.4f \n", eWedge, qqqGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]);
chRing = qqq.ch[j] - 16; chRing = qqq.ch[j] - 16;
eRingRaw = qqq.e[j]; eRingRaw = qqq.e[j];
eRing = qqq.e[j] * qqqrGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i]-16]; eRing = qqq.e[j];//* qqqrGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i]-16];
} }
else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16 && qqqrGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16] && qqqwGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16]) else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16/* && qqqrGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16] */&& qqqwGainValid[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16])
{ {
chWedge = qqq.ch[j]; chWedge = qqq.ch[j];
eWedge = qqq.e[j] * qqqwGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16]; eWedge = qqq.e[j] * qqqwGain[qqq.id[j]][qqq.ch[j]][qqq.ch[i] - 16];
eWedgeRaw = qqq.e[j]; eWedgeRaw = qqq.e[j];
chRing = qqq.ch[i] - 16; chRing = qqq.ch[i] - 16;
eRing = qqq.e[i] * qqqrGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16]; eRing = qqq.e[i];// * qqqrGain[qqq.id[i]][qqq.ch[i]][qqq.ch[j] - 16];
eRingRaw = qqq.e[i]; eRingRaw = qqq.e[i];
} }
else else
@ -151,7 +151,7 @@ Bool_t QQQ_Calcheck::Process(Long64_t entry)
plotter->Fill2D(Form("hRCal_qqq%d", qqq.id[i]), 16, 0, 15, 1000, 0, 30, chRing, eRingMeV, "RingCal"); plotter->Fill2D(Form("hRCal_qqq%d", qqq.id[i]), 16, 0, 15, 1000, 0, 30, chRing, eRingMeV, "RingCal");
plotter->Fill2D(Form("hWCal_qqq%d", qqq.id[i]), 16, 0, 15, 1000, 0, 30, chWedge, eWedgeMeV, "WedgeCal"); plotter->Fill2D(Form("hWCal_qqq%d", qqq.id[i]), 16, 0, 15, 1000, 0, 30, chWedge, eWedgeMeV, "WedgeCal");
plotter->Fill2D("hRawQQQ", 4000, 0, 8000, 4000, 0, 8000, eWedgeRaw, eRingRaw); plotter->Fill2D("hRawQQQ", 4000, 0, 8000, 4000, 0, 8000, eWedgeRaw, eRingRaw);
plotter->Fill2D("hGMQQQ", 4000, 0, 16000, 4000, 0, 16000, eWedge, eRing); plotter->Fill2D("hGMQQQ", 4000, 0, 8000, 4000, 0, 8000, eWedge, eRing);
plotter->Fill2D("hCalQQQ", 4000, 0, 10, 4000, 0, 10, eWedgeMeV, eRingMeV); plotter->Fill2D("hCalQQQ", 4000, 0, 10, 4000, 0, 10, eWedgeMeV, eRingMeV);
} }
} }