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

Compare commits

base: rtang:master
Branches Tags
rtang:master
rtang:devel_vignesh
rtang:devel_ruchi
..
compare: rtang:devel_vignesh
Branches Tags
rtang:devel_vignesh
rtang:master
rtang:devel_ruchi

27 Commits
master ... devel_vign

Author SHA1 Message Date
Vignesh Sitaraman 2225b3a942 modified: .gitignore 2025-02-17 17:06:33 -05:00
Vignesh Sitaraman 6cfb38b564 modified: Analyzer.C 2025-02-17 17:03:58 -05:00
Vignesh Sitaraman cb72c14ca4 modified: Analyzer.C modfied the algorithm for selection of valid cathodes, the section is now based on geometry 
modified:   Armory/ClassPW.h the rotation of the cathodes was reversed
 2025-02-17 16:59:08 -05:00
Vignesh Sitaraman 5995081396 modified: Analyzer.C 
modified:   Armory/ClassPW.h
    adjusted the Cathode geometry to  incorporate the offsett of 3 wires
 2025-02-05 14:59:02 -05:00
Vignesh Sitaraman d623e0cd17 modified: Analyzer.C 
made changes to correct for geomtry and eliminating redundant variables to a make the  code more readable
	modified:   Armory/ClassPW.h
    made changes to correct for the fact that the physical geometry of the detector is not the same as the geometry in the simulation, it is reversed and has an offset of 3
 2025-02-05 08:23:23 -05:00
Vignesh Sitaraman 8d7322cf5a modified: Analyzer.C 2025-02-03 10:32:07 -05:00
Vignesh Sitaraman 699b0f8701 modified: Analyzer.C 
HPZProjection implemented, testing seems to result in a plot with 0s. Need to debug
 2025-01-31 09:25:37 -05:00
Vignesh Sitaraman 02213caaee modified: Analyzer.C 
implementation of function to fgiure out the cqreelated cathodes in events
 2025-01-28 10:19:14 -05:00
Vignesh Sitaraman 56a6389b4f modified: Analyzer.C 
modified the analyser to include gain matching for the anodes and the cathodes
	modified:   FitHistogramsWithTSpectrum_Sequential_Improved.C
	modified:   MatchAndPlotCentroids.C
	modified:   centroids.txt
	modified:   slope_intercept_results.txt
 2025-01-27 16:34:39 -05:00
Vignesh Sitaraman b99ad4e4d7 new file: FitHistogramsWithTSpectrum_Sequential_Improved.C 
new file:   MatchAndPlotCentroids.C
	new file:   centroids.txt
	new file:   centroids_edited.txt
	new file:   slope_intercept_cathode.txt
	new file:   slope_intercept_results.txt
	new file:   slope_intercept_results_anode.txt
 2025-01-27 15:11:27 -05:00
Vignesh Sitaraman a5dfa2ecd3 modified: .vscode/settings.json 
modified:   Analyzer.C
	modified:   Analyzer.h
	modified:   ProcessRun.sh
 2025-01-27 10:55:22 -05:00
Vignesh Sitaraman 26e943adc8 modified: .vscode/settings.json 
modified:   Analyzer.C
 2025-01-27 09:51:25 -05:00
vs19g 42e093b104 modified: .vscode/c_cpp_properties.json 
modified:   .vscode/settings.json
	modified:   Analyzer.C
	new file:   Armory/#ClassPW.h#
	modified:   Armory/ClassDet.h
	modified:   Armory/ClassPW.h
	modified:   Armory/Mapper.cpp
 2025-01-17 11:18:07 -05:00
dirac 3129339647 modified: Analyzer.C 
modified:   ProcessRun.sh
	modified:   mapping.h
 2024-11-05 08:48:21 -05:00
dirac 7a70340b18 modified: Analyzer.C 
modified:   mapping.h
 2024-10-30 09:28:01 -04:00
dirac a10081ea81 modified: .vscode/settings.json 
modified:   Analyzer.C
	modified:   Analyzer.h
	modified:   Armory/ClassDet.h
	modified:   Armory/ClassPW.h
	modified:   Armory/Mapper.cpp
	modified:   PCGainMatch.C
	modified:   ProcessRun.sh
	modified:   mapping.h
 2024-10-25 15:02:59 -04:00
vs19g 4cb8a2c48c new file: Analyzer1.C 
new file:   Analyzer1.h
 2024-10-04 12:43:24 -04:00
vs19g fe6dbee171 Refactor PCGainMatch.C and Analyzer.C for improved efficiency and readability 2024-10-01 14:13:53 -04:00
vs19g 7805481ead Refactor PCGainMatch.C and Analyzer.C for improved efficiency and readability 2024-09-16 10:25:13 -04:00
vs19g 511b4aa808 modified: PCGainMatch.C 2024-09-06 15:06:32 -04:00
vs19g 43233ceb02 modified: PCGainMatch.C 2024-09-06 14:32:22 -04:00
vs19g 68fc36a8f6 modified: PCGainMatch.C 2024-09-03 16:41:17 -04:00
vs19g a6e754b958 modified: Analyzer.C 
new file:   PCGainMatch.C
	new file:   PCGainMatch.h
 2024-08-28 10:49:36 -04:00
vs19g 48ede97992 modified: Analyzer.C implemented cuts and histograms for the 2D histograms 2024-08-26 11:41:06 -04:00
vs19g f0a393abe2 modified: Analyzer.C 2024-08-26 10:34:13 -04:00
vs19g 4ba9c73b98 modified: Analyzer.C 2024-08-26 10:33:29 -04:00
vs19g 238ec8961e modified: Analyzer.C 2024-08-26 10:33:29 -04:00
22 changed files with 2053 additions and 1114 deletions
Show Stats Expand all files Collapse all files

3
.gitignore vendored
View File

@ -4,10 +4,9 @@ EventBuilder*
*.pcm
*.root
*.exe
*.txt
Mapper
AnasenMS
data/
data_proton/
root_data/

7
.vscode/settings.json vendored
View File

@ -100,7 +100,12 @@
"PCPulser_All_new.C": "cpp",
"PosCal_2.C": "cpp",
"AutoFit.C": "cpp",
"Fitting.C": "cpp"
"Fitting.C": "cpp",
"PCGainMatch.C": "cpp",
"Analyzer1.C": "cpp",
"FitHistogramsWithTSpectrum_Sequential_Improved.C": "cpp",
"PlotAndFitCentroids.C": "cpp",
"MatchAndPlotCentroids.C": "cpp"
},
"github-enterprise.uri": "https://fsunuc.physics.fsu.edu"
}

2
Analysis.C
View File

@ -16,5 +16,5 @@ void Analysis(int start, int end) {
// Define a macro with the same name as the script
void Analysis() {
Analysis(150, 194); // Adjust the range if needed
Analysis(72, 194); // Adjust the range if needed
}

476
Analyzer.C
View File

@ -1,26 +1,31 @@
#define Analyzer_cxx
#include "Analyzer.h"
#include "Armory/ClassSX3.h"
#include "Armory/ClassPW.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include "TVector3.h"
#include <fstream>
#include <iostream>
#include <sstream>
#include <map>
#include <utility>
#include <algorithm>
#include "Armory/ClassSX3.h"
#include "Armory/ClassPW.h"
#include "TVector3.h"
TH2F *hsx3IndexVE;
TH2F *hqqqIndexVE;
TH2F *hpcIndexVE;
TH2F *hpcIndexVE_GM;
TH2F *hsx3Coin;
TH2F *hqqqCoin;
TH2F *hpcCoin;
TH2F *hAVCcoin;
TH2F *hqqqPolar;
TH2F *hsx3VpcIndex;
@ -28,15 +33,29 @@ TH2F *hqqqVpcIndex;
TH2F *hqqqVpcE;
TH2F *hsx3VpcE;
TH2F *hanVScatsum;
TH2F *hanVScatsum_a[24];
TH1F *hPC_E[48];
TH1F *hCat4An;
TH1F *hCat0An;
TH1F *hAnodehits;
TH2F *hNosvAe;
TH2F *hUncorrCAN;
int padID = 0;
SX3 sx3_contr;
PW pw_contr;
PW pwinstance;
TVector3 hitPos;
// TVector3 anodeIntersection;
std::map<int, std::pair<double, double>> slopeInterceptMap;
bool HitNonZero;
bool sx3ecut;
bool qqqEcut;
TH1F *hZProj;
TH1F *hPCZProj;
void Analyzer::Begin(TTree * /*tree*/)
{
@ -46,12 +65,15 @@ void Analyzer::Begin(TTree * /*tree*/)
hsx3IndexVE->SetNdivisions(-612, "x");
hqqqIndexVE = new TH2F("hqqqIndexVE", "QQQ index vs Energy; QQQ index ; Energy", 4 * 2 * 16, 0, 4 * 2 * 16, 400, 0, 5000);
hqqqIndexVE->SetNdivisions(-1204, "x");
hpcIndexVE = new TH2F("hpcIndexVE", "PC index vs Energy; PC index ; Energy", 2 * 24, 0, 2 * 24, 400, 0, 4000);
hpcIndexVE = new TH2F("hpcIndexVE", "PC index vs Energy; PC index ; Energy", 2 * 24, 0, 2 * 24, 400, 0, 16000);
hpcIndexVE->SetNdivisions(-1204, "x");
hpcIndexVE_GM = new TH2F("hpcIndexVE_GM", "PC index vs Energy; PC index ; Energy", 2 * 24, 0, 2 * 24, 400, 0, 16000);
hpcIndexVE_GM->SetNdivisions(-1204, "x");
hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24 * 12, 0, 24 * 12, 24 * 12, 0, 24 * 12);
hqqqCoin = new TH2F("hqqqCoin", "QQQ Coincident", 4 * 2 * 16, 0, 4 * 2 * 16, 4 * 2 * 16, 0, 4 * 2 * 16);
hpcCoin = new TH2F("hpcCoin", "PC Coincident", 2 * 24, 0, 2 * 24, 2 * 24, 0, 2 * 24);
hAVCcoin = new TH2F("hAVCcoin", "Anode vs Cathode Coincident", 24, 0, 24, 24, 0, 24);
hqqqPolar = new TH2F("hqqqPolar", "QQQ Polar ID", 16 * 4, -TMath::Pi(), TMath::Pi(), 16, 10, 50);
@ -70,11 +92,54 @@ void Analyzer::Begin(TTree * /*tree*/)
hsx3VpcE->SetNdivisions(-612, "x");
hsx3VpcE->SetNdivisions(-12, "y");
hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600);
// hZProj = new TH1F("hZProj", "Z Projection", 1200, -600, 600);
hPCZProj = new TH1F("hPCZProj", "PC Z Projection", 600, -300, 300);
hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400, 0, 10000, 400, 0, 16000);
hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400, 0, 16000, 400, 0, 20000);
hCat4An = new TH1F("hCat4An", "Number of Cathodes/Anode", 24, 0, 24);
hCat0An = new TH1F("hCat0An", "Number of Cathodes without Anode", 24, 0, 24);
hAnodehits = new TH1F("hAnodehits", "Number of Anode hits", 24, 0, 24);
hNosvAe = new TH2F("hnosvAe", "Number of Cathodes/Anode vs Anode Energy", 20, 0, 20, 400, 0, 16000);
hUncorrCAN = new TH2F("hUncorrCAn", "Uncorrelated Cathodes/Anode", 24, 0, 24, 24, 0, 24);
// for (int i = 0; i < 24; i++)
// {
// TString histName = Form("hAnodeVsCathode_%d", i);
// TString histTitle = Form("Anode %d vs Cathode Sum; Anode E; Cathode Sum E", i);
// hanVScatsum_a[i] = new TH2F(histName, histTitle, 400, 0, 16000, 400, 0, 20000);
// }
// for (int i = 0; i < 48; i++)
// {
// TString histName = Form("hCathode_%d", i);
// TString histTitle = Form("Cathode_E_%d;", i);
// hPC_E[i] = new TH1F(histName, histTitle, 3200, 0, 32000);
// }
sx3_contr.ConstructGeo();
pw_contr.ConstructGeo();
std::ifstream inputFile("slope_intercept_results.txt");
if (inputFile.is_open())
{
std::string line;
int index;
double slope, intercept;
while (std::getline(inputFile, line))
{
std::stringstream ss(line);
ss >> index >> slope >> intercept;
// wires 37, 39, 44 have fit data that is incorrect or not present, they have thus been set to 1,0 (slope, intercept) for convenience
// wire 19 the 4th point was genereated using the slope of the line produced uising the other 3 points from the wire 1 vs wire 19 plot
if (index >= 0 && index <= 47)
{
slopeInterceptMap[index] = std::make_pair(slope, intercept);
}
}
inputFile.close();
}
else
{
std::cerr << "Error opening slope_intercept.txt" << std::endl;
}
}
Bool_t Analyzer::Process(Long64_t entry)
@ -112,13 +177,17 @@ Bool_t Analyzer::Process(Long64_t entry)
// ########################################################### Raw data
// //======================= SX3
sx3ecut = false;
std::vector<std::pair<int, int>> ID; // first = id, 2nd = index
for (int i = 0; i < sx3.multi; i++)
{
ID.push_back(std::pair<int, int>(sx3.id[i], i));
hsx3IndexVE->Fill(sx3.index[i], sx3.e[i]);
if (sx3.e[i] > 1000 && sx3.e[i] < 2200)
{
sx3ecut = true;
}
for (int j = i + 1; j < sx3.multi; j++)
{
@ -210,7 +279,7 @@ Bool_t Analyzer::Process(Long64_t entry)
// hitPos.Print();
}
}
qqqEcut = false;
// //======================= QQQ
for (int i = 0; i < qqq.multi; i++)
{
@ -218,6 +287,10 @@ Bool_t Analyzer::Process(Long64_t entry)
// if(pc.index[j]==4){
hqqqIndexVE->Fill(qqq.index[i], qqq.e[i]);
// }
if (qqq.e[i] > 1600 && qqq.e[i] < 2600)
{
qqqEcut = true;
}
// }
for (int j = 0; j < qqq.multi; j++)
{
@ -279,129 +352,315 @@ Bool_t Analyzer::Process(Long64_t entry)
}
// //======================= PC
ID.clear();
int counter = 0;
std::vector<std::pair<int, double>> E;
E.clear();
for (int i = 0; i < pc.multi; i++)
{
if (pc.e[i] > 100)
ID.push_back(std::pair<int, int>(pc.id[i], i));
if (pc.e[i] > 100)
E.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
hpcIndexVE->Fill(pc.index[i], pc.e[i]);
for (int j = i + 1; j < pc.multi; j++)
{
hpcCoin->Fill(pc.index[i], pc.index[j]);
}
}
// for( size_t i = 0; i < E.size(); i++) printf("%zu | %d %d \n", i, E[i].first, E[i].second );
// Calculate the crossover points and put them into an array
pwinstance.ConstructGeo();
Coord Crossover[24][24];
TVector3 a, c, diff;
double a2, ac, c2, adiff, cdiff, denom, alpha, beta;
Coord Crossover[24][24][2];
TVector3 a, c, diff, an, cn;
double a2, an2, ac, c2, cn2, 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);
ac = a.Dot(c);
c2 = c.Dot(c);
ac = a.Dot(c);
adiff = a.Dot(diff);
cdiff = c.Dot(diff);
an2 = an.Dot(an);
cn2 = cn.Dot(cn);
denom = a2 * c2 - ac * ac;
alpha = (ac * cdiff - c2 * adiff) / denom;
beta = (a2 * cdiff - ac * adiff) / denom;
Crossover[i][j].x = pwinstance.An[i].first.X() + alpha * a.X();
Crossover[i][j].y = pwinstance.An[i].first.Y() + alpha * a.Y();
Crossover[i][j].z = pwinstance.An[i].first.Z() + alpha * a.Z();
if (i == 23)
{
if (abs(i - j) < 7 || abs(i - j) > 17)
{
if (alpha < 0 && alpha > -1)
{
printf("Anode and cathode indices and coord : %d %d %f %f %f %f\n", i, j, pwinstance.Ca[j].first.X(), pwinstance.Ca[j].first.Y(), pwinstance.Ca[j].first.Z(), alpha);
printf("Crossover wires, points and alpha are : %f %f %f %f \n", Crossover[i][j].x, Crossover[i][j].y, Crossover[i][j].z, alpha);
}
}
}
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][2].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;
// printf("AID, CID, Crossover z and alpha are : %d %d %f %f \n", i, j, Crossover[i][j][0].z, Crossover[i][j][1].x /*this is alpha*/);
}
}
if (E.size() >= 3)
// if (i == 4)
{
// for (int k = -4; k < 3; k++)
// {
// if ((i + 24 + k) % 24 == 23 - j) // the 23-j is used to accomodate for the fact that the order of the cathodes was reversed
// {
// Crossover[i][j][1].y = 1;
// }
// }
// for (int k = -4; k < 3; k++)
// {
// if (Crossover[i][j][k].y != 1)
// // if (alpha < 1 && alpha >= -1)
// {
// // printf("i an: %d %f %f %f \n", i, an.X(), an.Y(), an.Z());
// Crossover[i][j][0].z = 9999999; // this is a placeholder value to indicate that the anode and cathode wires do not intersect
// }
// }
}
// 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++)
{
int aID = 0;
int cID = 0;
float aE = 0;
float cE = 0;
// if( ID[0].first < 1 ) {
// aID = pc.ch[ID[0].second];
// cID = pc.ch[ID[1].second];
// }else{
// cID = pc.ch[ID[0].second];
// aID = pc.ch[ID[1].second];
// }
// printf("anode= %d, cathode = %d\n", aID, cID);
for (int k = 0; k < qqq.multi; k++)
if (pc.e[i] > 100)
{
if (qqq.index[k] == 75 && pc.index[k] == 2 && pc.e[k] > 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
}
}
int multi_an = 0;
for (int l = 0; l < E.size(); l++)
{
if (E[l].first < 24)
{
multi_an++;
}
}
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;
if (multi_an >= 1)
{
for (int l = 0; l < E.size(); l++)
{
if (E[l].first < 24 && E[l].first != 19 && E[l].first != 12)
{
aE = E[l].second;
}
else if (E[l].first > 24)
{
cE = E[l].second;
}
}
}
// 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]);
}
}
hanVScatsum->Fill(aE, cE);
}
// 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; });
if (ID[0].first < 1)
if (anodeHits.size() >= 1 && cathodeHits.size() >=1)
{
for (const auto &anode : anodeHits)
{
aID = pc.ch[ID[0].second];
cID = pc.ch[ID[1].second];
}
else
{
cID = pc.ch[ID[0].second];
aID = pc.ch[ID[1].second];
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);
// }
}
if (HitNonZero)
// std::cout << " Anode iD : " << aIDMax << " Energy : " << aEMax << std::endl;
// printf("aID : %d, aE : %f, cE : %f\n", aID, aE, cE);
corrcatMax.clear();
for (const auto &cathode : cathodeHits)
{
pw_contr.CalTrack(hitPos, aID, cID);
hZProj->Fill(pw_contr.GetZ0());
cID = cathode.first;
cE = cathode.second;
// std::cout << "Cathode ID : " << cID << " Energy : " << cE << std::endl;
// if (cE > cEMax)
// {
// cIDnextMax = cIDMax;
// cEnextMax = cEMax;
// cEMax = cE;
// cIDMax = cID;
// }
// if (cE > cEnextMax && cE < cEMax)
// {
// cEnextMax = cE;
// cIDnextMax = 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
// bool corr = false;
hAVCcoin->Fill(aIDMax, cID);
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 */
{
corrcatMax.push_back(std::pair<int, double>(cID, 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;
cESum += cE;
// corr = true;
}
// if ((aIDnextMax + 24 + j) % 24 == cID)
// {
// corrcatnextMax.push_back(std::pair<int, double>(cathode.first, cathode.second));
// std::sort(corrcatMax.begin(), corrcatMax.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b)S
// { return a.second > b.second; });
// }
// for (int k = 0; k < 5; k++)
// {
// if ((aIDMax + 24 + j) % 24 == (aIDnextMax + 24 + k) % 24)
// {
// commcat.push_back(std::pair<int, double>(cathode.first, cathode.second));
// }
// }
}
// if (!corr)
// {
// hUncorrCAN->Fill(aIDMax, cID);
// }
}
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);
}
}
anodeIntersection = TVector3(x, y, z);
// std::cout << "Anode Intersection " << anodeIntersection.Z() << " " << x << " " << y << " " << z << std::endl;
}
// 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)
// {
// pw_contr.CalTrack(hitPos, aID, cID);
// hZProj->Fill(pw_contr.GetZ0());
// }
}
// ########################################################### Track constrcution
@ -490,4 +749,23 @@ void Analyzer::Terminate()
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();
}

15
Analyzer.h
View File

@ -18,6 +18,7 @@ public :
Det sx3;
Det qqq;
Det pc ;
Det misc;
ULong64_t evID;
UInt_t run;
@ -40,6 +41,13 @@ public :
TBranch *b_pcCh; //!
TBranch *b_pcE; //!
TBranch *b_pcT; //!
TBranch *b_miscMulti; //!
TBranch *b_miscID; //!
TBranch *b_miscCh; //!
TBranch *b_miscE; //!
TBranch *b_miscT; //!
TBranch *b_miscTf; //!
Analyzer(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~Analyzer() { }
@ -92,6 +100,13 @@ void Analyzer::Init(TTree *tree){
fChain->SetBranchAddress("pcCh", &pc.ch, &b_pcCh);
fChain->SetBranchAddress("pcE", &pc.e, &b_pcE);
fChain->SetBranchAddress("pcT", &pc.t, &b_pcT);
fChain->SetBranchAddress("miscMulti", &misc.multi, &b_miscMulti);
fChain->SetBranchAddress("miscID", &misc.id, &b_miscID);
fChain->SetBranchAddress("miscCh", &misc.ch, &b_miscCh);
fChain->SetBranchAddress("miscE", &misc.e, &b_miscE);
fChain->SetBranchAddress("miscT", &misc.t, &b_miscT);
// fChain->SetBranchAddress("miscF", &misc.tf, &b_miscTf);
}

402
Analyzer1.C Normal file
View File

@ -0,0 +1,402 @@
#define Analyzer1_cxx
#include "Analyzer1.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <utility>
#include <algorithm>
#include "Armory/ClassSX3.h"
#include "Armory/ClassPW.h"
#include "TVector3.h"
TH2F * hsx3IndexVE;
TH2F * hqqqIndexVE;
TH2F * hpcIndexVE;
TH2F * hsx3Coin;
TH2F * hqqqCoin;
TH2F * hpcCoin;
TH2F * hqqqPolar;
TH2F * hsx3VpcIndex;
TH2F * hqqqVpcIndex;
TH2F * hqqqVpcE;
TH2F * hsx3VpcE;
TH2F * hanVScatsum;
int padID = 0;
SX3 sx3_contr;
PW pw_contr;
TVector3 hitPos;
bool HitNonZero;
TH1F * hZProj;
void Analyzer1::Begin(TTree * /*tree*/){
TString option = GetOption();
hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24*12, 0, 24*12, 400, 0, 5000); hsx3IndexVE->SetNdivisions( -612, "x");
hqqqIndexVE = new TH2F("hqqqIndexVE", "QQQ index vs Energy; QQQ index ; Energy", 4*2*16, 0, 4*2*16, 400, 0, 5000); hqqqIndexVE->SetNdivisions( -1204, "x");
hpcIndexVE = new TH2F("hpcIndexVE", "PC index vs Energy; PC index ; Energy", 2*24, 0, 2*24, 400, 0, 4000); hpcIndexVE->SetNdivisions( -1204, "x");
hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24*12, 0, 24*12, 24*12, 0, 24*12);
hqqqCoin = new TH2F("hqqqCoin", "QQQ Coincident", 4*2*16, 0, 4*2*16, 4*2*16, 0, 4*2*16);
hpcCoin = new TH2F("hpcCoin", "PC Coincident", 2*24, 0, 2*24, 2*24, 0, 2*24);
hqqqPolar = new TH2F("hqqqPolar", "QQQ Polar ID", 16*4, -TMath::Pi(), TMath::Pi(),16, 10, 50);
hsx3VpcIndex = new TH2F("hsx3Vpcindex", "sx3 vs pc; sx3 index; pc index", 24*12, 0, 24*12, 48, 0, 48);
hsx3VpcIndex->SetNdivisions( -612, "x");
hsx3VpcIndex->SetNdivisions( -12, "y");
hqqqVpcIndex = new TH2F("hqqqVpcindex", "qqq vs pc; qqq index; pc index", 4*2*16, 0, 4*2*16, 48, 0, 48);
hqqqVpcIndex->SetNdivisions( -612, "x");
hqqqVpcIndex->SetNdivisions( -12, "y");
hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 400, 0, 5000, 400, 0, 5000);
hqqqVpcE->SetNdivisions( -612, "x");
hqqqVpcE->SetNdivisions( -12, "y");
hsx3VpcE = new TH2F("hsx3VpcEnergy", "sx3 vs pc; sx3 energy; pc energy", 400, 0, 5000, 400, 0, 5000);
hsx3VpcE->SetNdivisions( -612, "x");
hsx3VpcE->SetNdivisions( -12, "y");
hZProj = new TH1F("hZProj", "Z Projection", 1200, -600, 600);
hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 10000, 400, 0 , 16000);
sx3_contr.ConstructGeo();
pw_contr.ConstructGeo();
}
Bool_t Analyzer1::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);
b_sx3E->GetEntry(entry);
b_sx3T->GetEntry(entry);
b_qqqMulti->GetEntry(entry);
b_qqqID->GetEntry(entry);
b_qqqCh->GetEntry(entry);
b_qqqE->GetEntry(entry);
b_qqqT->GetEntry(entry);
b_pcMulti->GetEntry(entry);
b_pcID->GetEntry(entry);
b_pcCh->GetEntry(entry);
b_pcE->GetEntry(entry);
b_pcT->GetEntry(entry);
sx3.CalIndex();
qqq.CalIndex();
pc.CalIndex();
// sx3.Print();
//########################################################### Raw data
// //======================= SX3
std::vector<std::pair<int, int>> ID; // first = id, 2nd = index
for( int i = 0; i < sx3.multi; i ++){
ID.push_back(std::pair<int, int>(sx3.id[i], i));
hsx3IndexVE->Fill( sx3.index[i], sx3.e[i] );
for( int j = i+1; j < sx3.multi; j++){
hsx3Coin->Fill( sx3.index[i], sx3.index[j]);
}
for( int j = 0; j < pc.multi; j++){
hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
// if( sx3.ch[index] > 8 ){
// hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
// }
}
}
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;
} );
// printf("##############################\n");
// for( size_t i = 0; i < ID.size(); i++) printf("%zu | %d %d \n", i, ID[i].first, ID[i].second );
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]);
}
}
}
// printf("---------- sx3ID Multi : %zu \n", sx3ID.size());
if( found ){
int sx3ChUp, sx3ChDn, sx3ChBk;
float sx3EUp, sx3EDn;
// printf("------ sx3 ID : %d, multi: %zu\n", sx3ID[0].first, sx3ID.size());
for( size_t i = 0; i < sx3ID.size(); i++ ){
int index = sx3ID[i].second;
// printf(" %zu | index %d | ch : %d, energy : %d \n", i, index, sx3.ch[index], sx3.e[index]);
if( sx3.ch[index] < 8 ){
if( sx3.ch[index] % 2 == 0) {
sx3ChDn = sx3.ch[index];
sx3EDn = sx3.e[index];
}else{
sx3ChUp = sx3.ch[index];
sx3EUp = sx3.e[index];
}
}else{
sx3ChBk = sx3.ch[index];
}
for( int j = 0; j < pc.multi; j++){
// hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
if( sx3.ch[index] > 8 ){
hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
// hpcIndexVE->Fill( pc.index[i], pc.e[i] );
}
}
}
sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, sx3EUp, sx3EDn);
hitPos = sx3_contr.GetHitPos();
HitNonZero = true;
// hitPos.Print();
}
}
// //======================= QQQ
for( int i = 0; i < qqq.multi; i ++){
// for( int j = 0; j < pc.multi; j++){
// if(pc.index[j]==4){
hqqqIndexVE->Fill( qqq.index[i], qqq.e[i] );
// }
// }
for( int j = 0; j < qqq.multi; j++){
if ( j == i ) continue;
hqqqCoin->Fill( qqq.index[i], qqq.index[j]);
}
for( int j = i + 1; j < qqq.multi; j++){
for( int k = 0; k < pc.multi; k++){
if(pc.index[k]<24 && pc.e[k]>50 ){
hqqqVpcE->Fill( qqq.e[i], pc.e[k] );
// hpcIndexVE->Fill( pc.index[i], pc.e[i] );
hqqqVpcIndex->Fill( qqq.index[i], pc.index[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;
if( qqq.ch[i] < qqq.ch[j]){
chRing = qqq.ch[j] - 16;
chWedge = qqq.ch[i];
}else{
chRing = qqq.ch[i];
chWedge = qqq.ch[j] - 16;
}
// 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 = 10.+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
ID.clear();
int counter=0;
std::vector<std::pair<int, double>> E;
E.clear();
for( int i = 0; i < pc.multi; i ++){
if( pc.e[i] > 100 ) ID.push_back(std::pair<int, int>(pc.id[i], i));
if( pc.e[i] > 100 ) E.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
hpcIndexVE->Fill( pc.index[i], pc.e[i] );
for( int j = i+1; j < pc.multi; j++){
hpcCoin->Fill( pc.index[i], pc.index[j]);
}
}
// for( size_t i = 0; i < E.size(); i++) printf("%zu | %d %d \n", i, E[i].first, E[i].second );
if( E.size()>=3 ){
int aID = 0;
int cID = 0;
float aE = 0;
float cE = 0;
bool multi_an =false;
// if( ID[0].first < 1 ) {
// aID = pc.ch[ID[0].second];
// cID = pc.ch[ID[1].second];
// }else{
// cID = pc.ch[ID[0].second];
// aID = pc.ch[ID[1].second];
// }
// printf("anode= %d, cathode = %d\n", aID, cID);
// for( int k = 0; k < qqq.multi; k++){
// if(qqq.index[k]==75 && pc.index[k]==2 && pc.e[k]>100){
for(int l=0;l<E.size();l++){
if(E[l].first<24 ){
if(!multi_an){
aE = E[l].second;
}
multi_an=true;
}
else if (E[l].first>=24){
cE = E[l].second + cE;
}
}
// }
// }
hanVScatsum->Fill(aE,cE);
if( ID[0].first < 1 ) {
aID = pc.ch[ID[0].second];
cID = pc.ch[ID[1].second];
}else{
cID = pc.ch[ID[0].second];
aID = pc.ch[ID[1].second];
}
if( HitNonZero){
pw_contr.CalTrack( hitPos, aID, cID);
hZProj->Fill(pw_contr.GetZ0());
}
}
//########################################################### Track constrcution
//############################## DO THE KINEMATICS
return kTRUE;
}
void Analyzer1::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");
}

24
TrackRecon.h → Analyzer1.h
View File

@ -1,5 +1,5 @@
#ifndef TrackRecon_h
#define TrackRecon_h
#ifndef Analyzer1_h
#define Analyzer1_h
#include <TROOT.h>
#include <TChain.h>
@ -8,7 +8,7 @@
#include "Armory/ClassDet.h"
class TrackRecon : public TSelector {
class Analyzer1 : public TSelector {
public :
TTree *fChain; //!pointer to the analyzed TTree or TChain
@ -41,8 +41,8 @@ public :
TBranch *b_pcE; //!
TBranch *b_pcT; //!
TrackRecon(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~TrackRecon() { }
Analyzer1(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~Analyzer1() { }
virtual Int_t Version() const { return 2; }
virtual void Begin(TTree *tree);
virtual void SlaveBegin(TTree *tree);
@ -57,13 +57,13 @@ public :
virtual void SlaveTerminate();
virtual void Terminate();
ClassDef(TrackRecon,0);
ClassDef(Analyzer1,0);
};
#endif
#ifdef TrackRecon_cxx
void TrackRecon::Init(TTree *tree){
#ifdef Analyzer1_cxx
void Analyzer1::Init(TTree *tree){
// Set branch addresses and branch pointers
if (!tree) return;
@ -95,20 +95,20 @@ void TrackRecon::Init(TTree *tree){
}
Bool_t TrackRecon::Notify(){
Bool_t Analyzer1::Notify(){
return kTRUE;
}
void TrackRecon::SlaveBegin(TTree * /*tree*/){
void Analyzer1::SlaveBegin(TTree * /*tree*/){
TString option = GetOption();
}
void TrackRecon::SlaveTerminate(){
void Analyzer1::SlaveTerminate(){
}
#endif // #ifdef TrackRecon_cxx
#endif // #ifdef Analyzer_cxx

283
Armory/#ClassPW.h#
View File

@ -1,283 +0,0 @@
#ifndef ClassPW_h
#define ClassPW_h
#include <cstdio>
#include <TMath.h>
#include <TVector3.h>
struct PWHitInfo{
std::pair<short, short> nearestWire; // anode, cathode
std::pair<double, double> nearestDist; // anode, cathode
std::pair<short, short> nextNearestWire; // anode, cathode
std::pair<double, double> nextNearestDist; // anode, cathode
void Clear(){
nearestWire.first = -1;
nearestWire.second = -1;
nearestDist.first = 999999999;
nearestDist.second = 999999999;
nextNearestWire.first = -1;
nextNearestWire.second = -1;
nextNearestDist.first = 999999999;
nextNearestDist.second = 999999999;
}
};
//!########################################################
class PW{ // proportional wire
public:
PW(){ ClearHitInfo();};
~PW(){};
PWHitInfo GetHitInfo() const {return hitInfo;}
std::pair<short, short> GetNearestID() const {return hitInfo.nearestWire;}
std::pair<double, double> GetNearestDistance() const {return hitInfo.nearestDist;}
std::pair<short, short> Get2ndNearestID() const {return hitInfo.nextNearestWire;}
std::pair<double, double> Get2ndNearestDistance() const {return hitInfo.nextNearestDist;}
TVector3 GetTrackPos() const {return trackPos;}
TVector3 GetTrackVec() const {return trackVec;}
double GetTrackTheta() const {return trackVec.Theta();}
double GetTrackPhi() const {return trackVec.Phi();}
double GetZ0();
int GetNumWire() const {return nWire;}
double GetDeltaAngle() const {return dAngle;}
double GetAnodeLength() const {return anodeLength;}
double GetCathodeLength() const {return cathodeLength;}
TVector3 GetAnodeDn(short id) const {return An[id].first;}
TVector3 GetAnodeUp(short id) const {return An[id].second;}
TVector3 GetCathodeDn(short id) const {return Ca[id].first;}
TVector3 GetCathodeUp(short id) const {return Ca[id].second;}
TVector3 GetAnodneMid(short id) const {return (An[id].first + An[id].second) * 0.5; }
double GetAnodeTheta(short id) const {return (An[id].first - An[id].second).Theta();}
double GetAnodePhi(short id) const {return (An[id].first - An[id].second).Phi();}
TVector3 GetCathodneMid(short id) const {return (Ca[id].first + Ca[id].second) * 0.5; }
double GetCathodeTheta(short id) const {return (Ca[id].first - Ca[id].second).Theta();}
double GetCathodePhi(short id) const {return (Ca[id].first - Ca[id].second).Phi();}
void ClearHitInfo();
void ConstructGeo();
void FindWireID(TVector3 pos, TVector3 direction, bool verbose = false);
void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose = false);
void CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA = 0, double sigmaC = 0, bool verbose = false);
void Print(){
printf(" The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nearestWire.first,
hitInfo.nearestDist.first,
hitInfo.nearestWire.second,
hitInfo.nearestDist.second);
printf(" The 2nd nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nextNearestWire.first,
hitInfo.nextNearestDist.first,
hitInfo.nextNearestWire.second,
hitInfo.nextNearestDist.second);
}
private:
PWHitInfo hitInfo;
TVector3 trackPos;
TVector3 trackVec;
const int nWire = 24;
const int wireShift = 3;
const float zLen = 380; //mm
const float radiusA = 37;
const float radiusC = 43;
double dAngle;
double anodeLength;
double cathodeLength;
std::vector<std::pair<TVector3,TVector3>> An; // the anode wire position vector in space
std::vector<std::pair<TVector3,TVector3>> Ca; // the cathode wire position vector in space
double Distance(TVector3 a1, TVector3 a2, TVector3 b1, TVector3 b2){
TVector3 na = a1 - a2;
TVector3 nb = b1 - b2;
TVector3 nd = (na.Cross(nb)).Unit();
return TMath::Abs(nd.Dot(a1-b2));
}
};
inline void PW::ClearHitInfo(){
hitInfo.Clear();
}
inline void PW::ConstructGeo(){
An.clear();
Ca.clear();
std::pair<TVector3, TVector3> p1; // anode
std::pair<TVector3, TVector3> q1; // cathode
//anode and cathode start at pos-Y axis and count in right-Hand
//anode wire shift is right-hand.
//cathode wire shift is left-hand.
for(int i = 0; i < nWire; i++ ){
// Anode rotate right-hand
p1.first.SetXYZ( radiusA * TMath::Cos( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
radiusA * TMath::Sin( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
zLen/2);
p1.second.SetXYZ( radiusA * TMath::Cos( TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()),
radiusA * TMath::Sin( TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()),
-zLen/2);
An.push_back(p1);
// Cathod rotate left-hand
q1.first.SetXYZ( radiusC * TMath::Cos( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
radiusC * TMath::Sin( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
zLen/2);
q1.second.SetXYZ( radiusC * TMath::Cos( TMath::TwoPi() / nWire * (i - wireShift) + TMath::PiOver2()),
radiusC * TMath::Sin( TMath::TwoPi() / nWire * (i - wireShift) + TMath::PiOver2()),
-zLen/2);
Ca.push_back(q1);
}
dAngle = wireShift * TMath::TwoPi() / nWire;
anodeLength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusA * TMath::Sin(dAngle/2),2) );
cathodeLength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusC * TMath::Sin(dAngle/2),2) );
}
inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose ){
hitInfo.Clear();
double phi = direction.Phi();
for( int i = 0; i < nWire; i++){
double disA = 99999999;
double phiS = An[i].first.Phi() - TMath::PiOver4();
double phiL = An[i].second.Phi() + TMath::PiOver4();
// printf("A%2d: %f %f | %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg(), phi * TMath::RadToDeg());
if( phi > 0 && phiS > phiL ) phiL = phiL + TMath::TwoPi();
if( phi < 0 && phiS > phiL ) phiS = phiS - TMath::TwoPi();
if( phiS < phi && phi < phiL) {
disA = Distance( pos, pos + direction, An[i].first, An[i].second);
if( disA < hitInfo.nearestDist.first ){
hitInfo.nearestDist.first = disA;
hitInfo.nearestWire.first = i;
}
}
double disC = 99999999;
phiS = Ca[i].second.Phi()- TMath::PiOver4();
phiL = Ca[i].first.Phi() + TMath::PiOver4();
// printf("C%2d: %f %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg());
if( phi > 0 && phiS > phiL ) phiL = phiL + TMath::TwoPi();
if( phi < 0 && phiS > phiL ) phiS = phiS - TMath::TwoPi();
if(phiS < phi && phi < phiL) {
disC = Distance( pos, pos + direction, Ca[i].first, Ca[i].second);
if( disC < hitInfo.nearestDist.second ){
hitInfo.nearestDist.second = disC;
hitInfo.nearestWire.second = i;
}
}
if(verbose) printf(" %2d | %8.2f, %8.2f\n", i, disA, disC);
}
//==== find the 2nd nearest wire
short anode1 = hitInfo.nearestWire.first;
short aaa1 = anode1 - 1; if( aaa1 < 0 ) aaa1 += nWire;
short aaa2 = (anode1 + 1) % nWire;
double haha1 = Distance( pos, pos + direction, An[aaa1].first, An[aaa1].second);
double haha2 = Distance( pos, pos + direction, An[aaa2].first, An[aaa2].second);
if( haha1 < haha2){
hitInfo.nextNearestWire.first = aaa1;
hitInfo.nextNearestDist.first = haha1;
}else{
hitInfo.nextNearestWire.first = aaa2;
hitInfo.nextNearestDist.first = haha2;
}
short cathode1 = hitInfo.nearestWire.second;
short ccc1 = cathode1 - 1; if( ccc1 < 0 ) ccc1 += nWire;
short ccc2 = (cathode1 + 1) % nWire;
haha1 = Distance( pos, pos + direction, Ca[ccc1].first, Ca[ccc1].second);
haha2 = Distance( pos, pos + direction, Ca[ccc2].first, Ca[ccc2].second);
if( haha1 < haha2){
hitInfo.nextNearestWire.second = ccc1;
hitInfo.nextNearestDist.second = haha1;
}else{
hitInfo.nextNearestWire.second = ccc2;
hitInfo.nextNearestDist.second = haha2;
}
if( verbose ) Print();
}
inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose){
trackPos = sx3Pos;
TVector3 n1 = (An[anodeID].first - An[anodeID].second).Cross((sx3Pos - An[anodeID].second)).Unit();
TVector3 n2 = (Ca[cathodeID].first - Ca[cathodeID].second).Cross((sx3Pos - Ca[cathodeID].second)).Unit();
// if the handiness of anode and cathode revered, it should be n2 cross n1
trackVec = (n2.Cross(n1)).Unit();
if( verbose ) printf("Theta, Phi = %f, %f \n", trackVec.Theta() *TMath::RadToDeg(), trackVec.Phi()*TMath::RadToDeg());
}
inline void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA, double sigmaC, bool verbose){
trackPos = sx3Pos;
double p1 = TMath::Abs(hitInfo.nearestDist.first + gRandom->Gaus(0, sigmaA));
double p2 = TMath::Abs(hitInfo.nextNearestDist.first + gRandom->Gaus(0, sigmaA));
double fracA = p1 / (p1 + p2);
short anodeID1 = hitInfo.nearestWire.first;
short anodeID2 = hitInfo.nextNearestWire.first;
TVector3 shiftA1 = (An[anodeID2].first - An[anodeID1].first) * fracA;
TVector3 shiftA2 = (An[anodeID2].second - An[anodeID1].second) * fracA;
double q1 = TMath::Abs(hitInfo.nearestDist.second + gRandom->Gaus(0, sigmaC));
double q2 = TMath::Abs(hitInfo.nextNearestDist.second + gRandom->Gaus(0, sigmaC));
double fracC = q1 / (q1 + q2);
short cathodeID1 = hitInfo.nearestWire.second;
short cathodeID2 = hitInfo.nextNearestWire.second;
TVector3 shiftC1 = (Ca[cathodeID2].first - Ca[cathodeID1].first) * fracC;
TVector3 shiftC2 = (Ca[cathodeID2].second - Ca[cathodeID1].second) * fracC;
TVector3 a1 = An[anodeID1].first + shiftA1;
TVector3 a2 = An[anodeID1].second + shiftA2;
TVector3 c1 = Ca[cathodeID1].first + shiftC1;
TVector3 c2 = Ca[cathodeID1].second + shiftC2;
TVector3 n1 = (a1 - a2).Cross((sx3Pos - a2)).Unit();
TVector3 n2 = (c1 - c2).Cross((sx3Pos - c2)).Unit();
// if the handiness of anode and cathode revered, it should be n2 cross n1
trackVec = (n2.Cross(n1)).Unit();
if( verbose ) printf("Theta, Phi = %f, %f \n", trackVec.Theta() *TMath::RadToDeg(), trackVec.Phi()*TMath::RadToDeg());
}
inline double PW::GetZ0(){
double x = trackPos.X();
double y = trackPos.Y();
double rho = TMath::Sqrt(x*x + y*y);
double theta = trackVec.Theta();
return trackPos.Z() - rho / TMath::Tan(theta);
}
#endif

4
Armory/ClassPW.h
View File

@ -163,6 +163,10 @@ inline void PW::ConstructGeo()
-zLen / 2);
Ca.push_back(q1);
}
// correcting for the fact that the order of the cathode wires is reversed
std::reverse(Ca.begin(), Ca.end());
// adjusting for the 3 wire offset, the rbegin and rend are used as the rotation of the wires is done in the opposite direction i.e. 1,2,3 -> 3,1,2
std::rotate(Ca.rbegin(), Ca.rbegin() + 3, Ca.rend());
dAngle = wireShift * TMath::TwoPi() / nWire;
anodeLength = TMath::Sqrt(zLen * zLen + TMath::Power(2 * radiusA * TMath::Sin(dAngle / 2), 2));

124
FitHistogramsWithTSpectrum_Sequential_Improved.C Normal file
View File

@ -0,0 +1,124 @@
#include <TFile.h>
#include <TH1.h>
#include <TSpectrum.h>
#include <TF1.h>
#include <TCanvas.h>
#include <vector>
#include <iostream>
#include <algorithm>
#include <fstream>
#include <TText.h>
void FitHistogramsWithTSpectrum_Sequential_Improved() {
TFile *inputFile = new TFile("Histograms_anodes.root", "READ");
if (!inputFile || inputFile->IsZombie()) {
std::cerr << "Error opening the input file!" << std::endl;
return;
}
TCanvas *c1 = new TCanvas("c1", "Histogram Viewer", 800, 600);
// Open the output ASCII file to save the centroids
std::ofstream outFile("centroids.txt");
if (!outFile.is_open()) {
std::cerr << "Error opening output file!" << std::endl;
return;
}
outFile << "HistogramIndex\tPeakNumber\tCentroid\tAmplitude\tSigma" << std::endl;
for (int i = 0; i < 24; ++i) {
TH1 *histogram = dynamic_cast<TH1*>(inputFile->Get(Form("hCathode_%d", i)));
if (!histogram) {
std::cerr << "Failed to retrieve histogram_" << i << " from the file." << std::endl;
continue;
}
// Set range for peak search
double minX = 700;
double maxX = 25000;
histogram->GetXaxis()->SetRangeUser(minX, maxX);
// Draw the histogram
c1->cd();
histogram->Draw();
// Peak search using TSpectrum
const int maxPeaks = 5;
TSpectrum spectrumFinder(maxPeaks);
int nFound = spectrumFinder.Search(histogram, 2, "", 0.01);
if (nFound <= 0) {
std::cerr << "No peaks found for histogram " << i << std::endl;
continue;
}
Double_t *xPositions = spectrumFinder.GetPositionX();
Double_t *yPositions = spectrumFinder.GetPositionY();
std::vector<std::pair<Double_t, Double_t>> peaks;
// Collect and sort peaks by X position
for (int j = 0; j < nFound; ++j) {
peaks.emplace_back(xPositions[j], yPositions[j]);
}
std::sort(peaks.begin(), peaks.end());
// Fit each peak with a Gaussian
for (int j = 0; j < peaks.size(); ++j) {
Double_t peakX = peaks[j].first;
Double_t peakY = peaks[j].second;
Double_t initialAmplitude = peakY; // Better initial guess
Double_t initialCentroid = peakX; // Centroid based on peak position
Double_t initialSigma = 60.0;
// Define Gaussian with initial parameters
TF1 *gaussFit = new TF1(Form("gauss_%d", j), "gaus", peakX - 200, peakX + 200);
//gaussFit->SetParameters(peakY, peakX, 25.0); // Initial guesses for amplitude, mean, sigma
gaussFit->SetParameters(initialAmplitude, initialCentroid, initialSigma);
// Perform fit
int fitStatus = histogram->Fit(gaussFit, "RQ+");
if (fitStatus != 0) {
std::cerr << "Fit failed for peak " << j + 1 << " in histogram " << i << std::endl;
delete gaussFit;
continue;
}
// Retrieve fit parameters
double amplitude = gaussFit->GetParameter(0);
double centroid = gaussFit->GetParameter(1);
double sigma = gaussFit->GetParameter(2);
double amplitudeError = gaussFit->GetParError(0);
double centroidError = gaussFit->GetParError(1);
double sigmaError = gaussFit->GetParError(2);
// Chi-squared value
double chi2 = gaussFit->GetChisquare();
int ndf = gaussFit->GetNDF();
outFile << i << "\t" << j + 1 << "\t" << centroid << std::endl;
gaussFit->SetLineColor(kRed);
gaussFit->Draw("SAME");
TText *text = new TText();
text->SetNDC();
text->SetTextSize(0.03);
text->SetTextColor(kRed);
//text->DrawText(0.15, 0.8 - j * 0.05, Form("Peak %d: Amp=%.2f, Mean=%.2f, Sigma=%.2f", j + 1, amplitude, centroid, sigma));
text->DrawText(0.15, 0.8 - j * 0.05,
Form("Peak %d: Amp=%.2f±%.2f, Mean=%.2f±%.2f, Sigma=%.2f±%.2f, Chi2/NDF=%.2f",
j + 1, amplitude, amplitudeError, centroid, centroidError, sigma, sigmaError, chi2 / ndf));
// Save results
// Clean up
delete gaussFit;
}
// Update canvas for visualization
c1->Update();
std::cout << "Press Enter to view the next histogram..." << std::endl;
c1->WaitPrimitive(); // Wait until Enter is pressed in the ROOT console
}
// Close resources
inputFile->Close();
outFile.close();
delete c1;
}

133
MatchAndPlotCentroids.C Normal file
View File

@ -0,0 +1,133 @@
#include <fstream>
#include <sstream>
#include <vector>
#include <map>
#include <iostream>
#include <TGraph.h>
#include <TF1.h>
#include <TCanvas.h>
#include <TH1.h>
void MatchAndPlotCentroids() {
// Open the centroid data file
std::ifstream inputFile("centroids.txt");
if (!inputFile.is_open()) {
std::cerr << "Error: Could not open Centroids.txt" << std::endl;
return;
}
// Data structure to store centroids by histogram and peak number
std::map<int, std::map<int, double>> centroidData;
// Read data from the file
std::string line;
while (std::getline(inputFile, line)) {
std::istringstream iss(line);
int histogramIndex, peakNumber;
double centroid;
if (iss >> histogramIndex >> peakNumber >> centroid) {
centroidData[histogramIndex][peakNumber] = centroid;
}
}
inputFile.close();
// Ensure histogram 24 exists and has data
if (centroidData.find(1) == centroidData.end()) {
std::cerr << "Error: Histogram 0 not found in the data!" << std::endl;
return;
}
// Reference centroids from histogram 24
const auto& referenceCentroids = centroidData[1];
std::ofstream outputFile("slope_intercept_results.txt");
if (!outputFile.is_open()) {
std::cerr << "Error: Could not open the output file for writing!" << std::endl;
return;
}
outputFile << "Histogram Number\tSlope\tIntercept\n";
// Loop through histograms 25 to 47
for (int targetHist = 0; targetHist <= 23; targetHist++) {
// Ensure the target histogram exists and matches in peak numbers
if (centroidData.find(targetHist) == centroidData.end() || centroidData[targetHist].size() != referenceCentroids.size()) {
//4th cnetroid data point for 19 was generated using the 3 datqa points for the slope of wires 0 and 19
std::cout << "Skipping Histogram " << targetHist << " due to mismatched or missing data." << std::endl;
continue;
}
// Prepare x and y values for TGraph
std::vector<double> xValues, yValues;
for (const auto& [peakNumber, refCentroid] : referenceCentroids) {
if (centroidData[targetHist].find(peakNumber) != centroidData[targetHist].end()) {
yValues.push_back(refCentroid);
xValues.push_back(centroidData[targetHist][peakNumber]);
} else {
std::cerr << "Warning: Peak " << peakNumber << " missing in histogram " << targetHist << std::endl;
}
}
if (xValues.size() < 3) {
std::cout << "Skipping Histogram " << targetHist << " as it has less than 3 matching centroids." << std::endl;
continue;
}
// Create a TGraph
TCanvas *c1 = new TCanvas(Form("c_centroid_1_vs_%d", targetHist), Form("Centroid 1 vs %d", targetHist), 800, 600);
TGraph *graph = new TGraph(xValues.size(), &xValues[0], &yValues[0]);
graph->SetTitle(Form("Centroid of Histogram %d vs 1", targetHist));
graph->GetYaxis()->SetTitle("Centroid of Histogram 1");
graph->GetXaxis()->SetTitle(Form("Centroid of Histogram %d", targetHist));
graph->SetMarkerStyle(20); // Full circle marker
graph->SetMarkerSize(1.0);
graph->SetMarkerColor(kBlue);
// Draw the graph
graph->Draw("AP");
double minX = *std::min_element(xValues.begin(), xValues.end());
double maxX = *std::max_element(xValues.begin(), xValues.end());
// Fit the data with a linear function
TF1 *fitLine = new TF1("fitLine", "pol1", minX, maxX); // Adjust range as needed
fitLine->SetLineColor(kRed); // Set the line color to distinguish it
fitLine->SetLineWidth(2); // Thicker line for visibility
graph->Fit(fitLine, "M");
fitLine->Draw("same");
fitLine->SetParLimits(0, -10, 10); // Limit intercept between -10 and 10
fitLine->SetParLimits(1, 0, 2);
// Extract slope and intercept
double slope = fitLine->GetParameter(1);
double intercept = fitLine->GetParameter(0);
outputFile << targetHist << "\t" << slope << "\t" << intercept << "\n";
std::cout << "Histogram 24 vs " << targetHist << ": Slope = " << slope << ", Intercept = " << intercept << std::endl;
std::vector<double> residuals;
for (size_t i = 0; i < xValues.size(); ++i) {
double fittedY = fitLine->Eval(xValues[i]); // Evaluate fitted function at x
double residual = yValues[i] - fittedY; // Residual = observed - fitted
residuals.push_back(residual);
}
// Create a graph for the residuals
/*TGraph *residualGraph = new TGraph(residuals.size(), &xValues[0], &residuals[0]);
residualGraph->SetTitle(Form("Residuals for Histogram 24 vs %d", targetHist));
residualGraph->GetYaxis()->SetTitle("Residuals");
residualGraph->GetXaxis()->SetTitle(Form("Centroid of Histogram %d", targetHist));
residualGraph->SetMarkerStyle(20);
residualGraph->SetMarkerSize(1.0);
residualGraph->SetMarkerColor(kGreen);
// Draw the residuals plot below the original plot (can be on a new canvas if preferred)
TCanvas *c2 = new TCanvas(Form("c_residuals_24_vs_%d", targetHist), Form("Residuals for Centroid 24 vs %d", targetHist), 800, 400);
residualGraph->Draw("AP");*/
c1->Update();
//c2->Update();
std::cout << "Press Enter to continue..." << std::endl;
//std::cin.get();
c1->WaitPrimitive();
//c2->WaitPrimitive();
//std::cin.get();
//std::cin.get();
}
outputFile.close();
std::cout << "Results written to slope_intercept_results.txt" << std::endl;
}

454
PCGainMatch.C Normal file
View File

@ -0,0 +1,454 @@
#define PCGainMatch_cxx
#include "PCGainMatch.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TCutG.h>
#include <utility>
#include <algorithm>
#include "Armory/ClassSX3.h"
#include "Armory/ClassPW.h"
#include "TVector3.h"
TH2F * hsx3IndexVE;
TH2F * hqqqIndexVE;
TH2F * hpcIndexVE;
TH2F * hsx3Coin;
TH2F * hqqqCoin;
TH2F * hpcCoin;
TH2F * hqqqPolar;
TH2F * hsx3VpcIndex;
TH2F * hqqqVpcIndex;
TH2F * hqqqVpcE;
TH2F * hsx3VpcE;
TH2F * hanVScatsum;
TH2F * hanVScatsum_a[24];
TH2F * hanVScatsum_hcut;
TH2F * hanVScatsum_lcut;
TH2F * hAnodeHits;
TH1F * hAnodeMultiplicity;
int padID = 0;
SX3 sx3_contr;
PW pw_contr;
TVector3 hitPos;
bool HitNonZero;
TH1F * hZProj;
TCutG *AnCatSum_high;
TCutG *AnCatSum_low;
TCutG *PCCoinc_cut1;
TCutG *PCCoinc_cut2;
bool inCuth;
bool inCutl;
bool inPCCut;
void PCGainMatch::Begin(TTree * /*tree*/){
TString option = GetOption();
hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24*12, 0, 24*12, 400, 0, 5000); hsx3IndexVE->SetNdivisions( -612, "x");
hqqqIndexVE = new TH2F("hqqqIndexVE", "QQQ index vs Energy; QQQ index ; Energy", 4*2*16, 0, 4*2*16, 400, 0, 5000); hqqqIndexVE->SetNdivisions( -1204, "x");
hpcIndexVE = new TH2F("hpcIndexVE", "PC index vs Energy; PC index ; Energy", 2*24, 0, 2*24, 800, 0, 16000); hpcIndexVE->SetNdivisions( -1204, "x");
hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24*12, 0, 24*12, 24*12, 0, 24*12);
hqqqCoin = new TH2F("hqqqCoin", "QQQ Coincident", 4*2*16, 0, 4*2*16, 4*2*16, 0, 4*2*16);
hpcCoin = new TH2F("hpcCoin", "PC Coincident", 2*24, 0, 2*24, 2*24, 0, 2*24);
hqqqPolar = new TH2F("hqqqPolar", "QQQ Polar ID", 16*4, -TMath::Pi(), TMath::Pi(),16, 10, 50);
hsx3VpcIndex = new TH2F("hsx3Vpcindex", "sx3 vs pc; sx3 index; pc index", 24*12, 0, 24*12, 48, 0, 48);
hsx3VpcIndex->SetNdivisions( -612, "x");
hsx3VpcIndex->SetNdivisions( -12, "y");
hqqqVpcIndex = new TH2F("hqqqVpcindex", "qqq vs pc; qqq index; pc index", 4*2*16, 0, 4*2*16, 48, 0, 48);
hqqqVpcIndex->SetNdivisions( -612, "x");
hqqqVpcIndex->SetNdivisions( -12, "y");
hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 400, 0, 5000, 400, 0, 5000);
hqqqVpcE->SetNdivisions( -612, "x");
hqqqVpcE->SetNdivisions( -12, "y");
hsx3VpcE = new TH2F("hsx3VpcEnergy", "sx3 vs pc; sx3 energy; pc energy", 400, 0, 5000, 400, 0, 5000);
hsx3VpcE->SetNdivisions( -612, "x");
hsx3VpcE->SetNdivisions( -12, "y");
hZProj = new TH1F("hZProj", "Nos of anodes", 20, 0, 19);
hAnodeHits = new TH2F("hAnodeHits", "Anode vs Anode Energy, Anode ID; Anode E", 24,0 , 23, 400, 0 , 20000);
hAnodeMultiplicity = new TH1F("hAnodeMultiplicity", "Number of Anodes/Event", 40, 0, 40);
hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 10000, 800, 0 , 16000);
for (int i = 0; i < 24; i++) {
TString histName = Form("hAnodeVsCathode_%d", i);
TString histTitle = Form("Anode %d vs Cathode Sum; Anode E; Cathode Sum E", i);
hanVScatsum_a[i] = new TH2F(histName, histTitle, 400, 0, 10000, 400, 0, 16000);
}
hanVScatsum_lcut = new TH2F("hanVScatsum_LCUT", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 16000, 400, 0 , 16000);
hanVScatsum_hcut = new TH2F("hanVScatsum_HCUT", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 16000, 400, 0 , 16000);
sx3_contr.ConstructGeo();
pw_contr.ConstructGeo();
// TFile *f1 = new TFile("AnCatSum_high.root");
// TFile *f2 = new TFile("AnCatSum_low.root");
// TFile *f3 = new TFile("PCCoinc_cut1.root");
// TFile *f4 = new TFile("PCCoinc_cut2.root");
// AnCatSum_high= (TCutG*)f1->Get("AnCatSum_high");
// AnCatSum_low= (TCutG*)f2->Get("AnCatSum_low");
// PCCoinc_cut1= (TCutG*)f3->Get("PCCoinc_cut1");
// PCCoinc_cut2= (TCutG*)f4->Get("PCCoinc_cut2");
}
Bool_t PCGainMatch::Process(Long64_t entry){
// if (entry % 1000000 == 0) {
// std::cout << "Processing entry: " << entry << std::endl;
// }
// 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);
b_sx3E->GetEntry(entry);
b_sx3T->GetEntry(entry);
b_qqqMulti->GetEntry(entry);
b_qqqID->GetEntry(entry);
b_qqqCh->GetEntry(entry);
b_qqqE->GetEntry(entry);
b_qqqT->GetEntry(entry);
b_pcMulti->GetEntry(entry);
b_pcID->GetEntry(entry);
b_pcCh->GetEntry(entry);
b_pcE->GetEntry(entry);
b_pcT->GetEntry(entry);
sx3.CalIndex();
qqq.CalIndex();
pc.CalIndex();
// sx3.Print();
//########################################################### Raw data
// //======================= PC
std::vector<std::pair<int, double>> anodeHits={};
std::vector<std::pair<int, double>> cathodeHits={};
int aID = 0;
int cID = 0;
float aE = 0;
float cE = 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]>50 && pc.multi<7){
float aESum = 0;
float cESum = 0;
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], 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]);
}
if (anodeHits.size()==1 && cathodeHits.size() >= 1) {
for (const auto& anode : anodeHits) {
// for(int l=0; l<sx3.multi; l++){
// if (sx3.index[l]==80){
aID = anode.first;
aE = anode.second;
aESum += aE;
printf("aID : %d, aE : %f\n", aID, aE);
}
printf("aID : %d, aE : %f, cE : %f\n", aID, aE, cE);
for (const auto& cathode : cathodeHits) {
cID = cathode.first;
cE = cathode.second;
// if(cE>cEMax){
// cEMax = cE;
// cIDMax = cID;
// }
// if(cE>cEnextMax && cE<cEMax){
// cEnextMax = cE;
// cIDnextMax = cID;
// }
cESum += cE;
}
// }
// 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(aESum, cESum);
// }
if (aID < 24 && aE > 50) {
hanVScatsum_a[aID]->Fill(aE, cESum);
}
// }
// Fill histograms for the `pc` data
hpcIndexVE->Fill(pc.index[i], pc.e[i]);
// if(inPCCut){
hAnodeMultiplicity->Fill(anodeHits.size());
// }
}
}
}
// //======================= SX3
std::vector<std::pair<int, int>> ID; // first = id, 2nd = index
for( int i = 0; i < sx3.multi; i ++){
if(sx3.e[i]>50){
ID.push_back(std::pair<int, int>(sx3.id[i], i));
hsx3IndexVE->Fill( sx3.index[i], sx3.e[i] );
for( int j = i+1; j < sx3.multi; j++){
hsx3Coin->Fill( sx3.index[i], sx3.index[j]);
}
for( int j = 0; j < pc.multi; j++){
hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
// if( sx3.ch[index] > 8 ){
// hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
// }
}
}
}
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;
} );
// printf("##############################\n");
// for( size_t i = 0; i < ID.size(); i++) printf("%zu | %d %d \n", i, ID[i].first, ID[i].second );
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]);
}
}
}
// printf("---------- sx3ID Multi : %zu \n", sx3ID.size());
if( found ){
int sx3ChUp, sx3ChDn, sx3ChBk;
float sx3EUp, sx3EDn;
// printf("------ sx3 ID : %d, multi: %zu\n", sx3ID[0].first, sx3ID.size());
for( size_t i = 0; i < sx3ID.size(); i++ ){
int index = sx3ID[i].second;
// printf(" %zu | index %d | ch : %d, energy : %d \n", i, index, sx3.ch[index], sx3.e[index]);
if( sx3.ch[index] < 8 ){
if( sx3.ch[index] % 2 == 0) {
sx3ChDn = sx3.ch[index];
sx3EDn = sx3.e[index];
}else{
sx3ChUp = sx3.ch[index];
sx3EUp = sx3.e[index];
}
}else{
sx3ChBk = sx3.ch[index];
}
for( int j = 0; j < pc.multi; j++){
// hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
if( sx3.ch[index] > 8 && pc.index[j]<24 && pc.e[j]>50 ){
hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
// hpcIndexVE->Fill( pc.index[i], pc.e[i] );
}
}
}
sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, sx3EUp, sx3EDn);
hitPos = sx3_contr.GetHitPos();
HitNonZero = true;
// hitPos.Print();
}
}
// //======================= QQQ
for( int i = 0; i < qqq.multi; i ++){
// for( int j = 0; j < pc.multi; j++){
if(qqq.e[i]>50 ){
hqqqIndexVE->Fill( qqq.index[i], qqq.e[i] );
for( int j = 0; j < qqq.multi; j++){
if ( j == i ) continue;
hqqqCoin->Fill( qqq.index[i], qqq.index[j]);
}
for( int j = i + 1; j < qqq.multi; j++){
for( int k = 0; k < pc.multi; k++){
// if(qqq.e[i>50]){
hqqqVpcE->Fill( qqq.e[i], pc.e[k] );
hqqqVpcIndex->Fill( qqq.index[i], pc.index[j] );
}
// }
}
}
// }
}
// hanVScatsum->Fill(aE,cE);
if( HitNonZero){
pw_contr.CalTrack( hitPos, aID, cID);
hZProj->Fill(pw_contr.GetZ0());
}
//########################################################### Track constrcution
//############################## DO THE KINEMATICS
return kTRUE;
}
void PCGainMatch::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);
hanVScatsum->Draw("colz");
// hAnodeHits->Draw("colz");
// hAnodeMultiplicity->Draw();
}

228
gainmatch.h → PCGainMatch.h
View File

@ -1,114 +1,114 @@
#ifndef gainmatch_h
#define gainmatch_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TSelector.h>
#include "Armory/ClassDet.h"
class gainmatch : 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; //!
gainmatch(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~gainmatch() { }
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(gainmatch,0);
};
#endif
#ifdef gainmatch_cxx
void gainmatch::Init(TTree *tree){
// Set branch addresses and branch pointers
if (!tree) return;
fChain = tree;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("evID", &evID, &b_eventID);
fChain->SetBranchAddress("run", &run, &b_run);
sx3.SetDetDimension(24,12);
qqq.SetDetDimension(4,32);
pc.SetDetDimension(2,24);
fChain->SetBranchAddress("sx3Multi", &sx3.multi, &b_sx3Multi);
fChain->SetBranchAddress("sx3ID", &sx3.id, &b_sx3ID);
fChain->SetBranchAddress("sx3Ch", &sx3.ch, &b_sx3Ch);
fChain->SetBranchAddress("sx3E", &sx3.e, &b_sx3E);
fChain->SetBranchAddress("sx3T", &sx3.t, &b_sx3T);
fChain->SetBranchAddress("qqqMulti", &qqq.multi, &b_qqqMulti);
fChain->SetBranchAddress("qqqID", &qqq.id, &b_qqqID);
fChain->SetBranchAddress("qqqCh", &qqq.ch, &b_qqqCh);
fChain->SetBranchAddress("qqqE", &qqq.e, &b_qqqE);
fChain->SetBranchAddress("qqqT", &qqq.t, &b_qqqT);
fChain->SetBranchAddress("pcMulti", &pc.multi, &b_pcMulti);
fChain->SetBranchAddress("pcID", &pc.id, &b_pcID);
fChain->SetBranchAddress("pcCh", &pc.ch, &b_pcCh);
fChain->SetBranchAddress("pcE", &pc.e, &b_pcE);
fChain->SetBranchAddress("pcT", &pc.t, &b_pcT);
}
Bool_t gainmatch::Notify(){
return kTRUE;
}
void gainmatch::SlaveBegin(TTree * /*tree*/){
TString option = GetOption();
}
void gainmatch::SlaveTerminate(){
}
#endif // #ifdef gainmatch_cxx
#ifndef PCGainMatch_h
#define PCGainMatch_h
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include <TSelector.h>
#include "Armory/ClassDet.h"
class PCGainMatch : 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; //!
PCGainMatch(TTree * /*tree*/ =0) : fChain(0) { }
virtual ~PCGainMatch() { }
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(PCGainMatch,0);
};
#endif
#ifdef PCGainMatch_cxx
void PCGainMatch::Init(TTree *tree){
// Set branch addresses and branch pointers
if (!tree) return;
fChain = tree;
fChain->SetMakeClass(1);
fChain->SetBranchAddress("evID", &evID, &b_eventID);
fChain->SetBranchAddress("run", &run, &b_run);
sx3.SetDetDimension(24,12);
qqq.SetDetDimension(4,32);
pc.SetDetDimension(2,24);
fChain->SetBranchAddress("sx3Multi", &sx3.multi, &b_sx3Multi);
fChain->SetBranchAddress("sx3ID", &sx3.id, &b_sx3ID);
fChain->SetBranchAddress("sx3Ch", &sx3.ch, &b_sx3Ch);
fChain->SetBranchAddress("sx3E", &sx3.e, &b_sx3E);
fChain->SetBranchAddress("sx3T", &sx3.t, &b_sx3T);
fChain->SetBranchAddress("qqqMulti", &qqq.multi, &b_qqqMulti);
fChain->SetBranchAddress("qqqID", &qqq.id, &b_qqqID);
fChain->SetBranchAddress("qqqCh", &qqq.ch, &b_qqqCh);
fChain->SetBranchAddress("qqqE", &qqq.e, &b_qqqE);
fChain->SetBranchAddress("qqqT", &qqq.t, &b_qqqT);
fChain->SetBranchAddress("pcMulti", &pc.multi, &b_pcMulti);
fChain->SetBranchAddress("pcID", &pc.id, &b_pcID);
fChain->SetBranchAddress("pcCh", &pc.ch, &b_pcCh);
fChain->SetBranchAddress("pcE", &pc.e, &b_pcE);
fChain->SetBranchAddress("pcT", &pc.t, &b_pcT);
}
Bool_t PCGainMatch::Notify(){
return kTRUE;
}
void PCGainMatch::SlaveBegin(TTree * /*tree*/){
TString option = GetOption();
}
void PCGainMatch::SlaveTerminate(){
}
#endif // #ifdef Analyzer_cxx

14
ProcessRun.sh
View File

@ -1,7 +1,8 @@
#!/bin/bash
if [ "$#" -ne 2 ]; then
echo "Usage: $0 runID timeWindow_ns"
if [ "$#" -ne 3 ]; then
echo "Usage: $0 runID timeWindow_ns option"
echo "option: 0 - process raw data, 1 - process mapped data"
echo "Exiting..."
exit 1
fi
@ -9,7 +10,9 @@ fi
runID=$1
timeWindow=$2
rawFolder=/home/tandem/Desktop/analysis/data
option=$3
rawFolder=/home/tandem/data1/2024_09_17Fap/data
rootFolder=../root_data
if [ $option -eq 0 ]; then
@ -18,7 +21,7 @@ if [ $option -eq 0 ]; then
fileList=`\ls -1 ${rawFolder}/*Run_${runID}_*.fsu`
./EventBuilder ${timeWindow} 0 0 100000000 ${fileList}
./EventBuilder ${timeWindow} 0 0 100000000 ${fileList}
outFile=${rawFolder}/*${runID}*${timeWindow}.root
@ -26,4 +29,5 @@ if [ $option -eq 0 ]; then
./Mapper ${rootFolder}/*${runID}*${timeWindow}.root
fi
root "processRun.C(\"${rootFolder}/Run_${runID}_mapped.root\")"
root "processRun.C(\"${rootFolder}/ProtonRun_${runID}_mapped.root\")"

152
TrackRecon.C
View File

@ -1,6 +1,6 @@
#define TrackRecon_cxx
#define Analyzer_cxx
#include "TrackRecon.h"
#include "Analyzer.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
@ -14,6 +14,20 @@
#include "TVector3.h"
TH2F * hsx3IndexVE;
TH2F * hqqqIndexVE;
TH2F * hpcIndexVE;
TH2F * hsx3Coin;
TH2F * hqqqCoin;
TH2F * hpcCoin;
TH2F * hqqqPolar;
TH2F * hsx3VpcIndex;
TH2F * hqqqVpcIndex;
TH2F * hqqqVpcE;
TH2F * hsx3VpcE;
TH2F * hanVScatsum;
int padID = 0;
SX3 sx3_contr;
@ -23,7 +37,7 @@ bool HitNonZero;
TH1F * hZProj;
void TrackRecon::Begin(TTree * /*tree*/){
void Analyzer::Begin(TTree * /*tree*/){
TString option = GetOption();
hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600);
@ -33,14 +47,17 @@ void TrackRecon::Begin(TTree * /*tree*/){
}
Bool_t TrackRecon::Process(Long64_t entry){
Bool_t Analyzer::Process(Long64_t entry){
// if ( entry > 100 ) return kTRUE;
hitPos.Clear();
HitNonZero = false;
if( entry > 1) return kTRUE;
// if( entry > 1) return kTRUE;
// printf("################### ev : %llu \n", entry);
b_sx3Multi->GetEntry(entry);
@ -116,6 +133,13 @@ Bool_t TrackRecon::Process(Long64_t entry){
}else{
sx3ChBk = sx3.ch[index];
}
for( int j = 0; j < pc.multi; j++){
// hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
if( sx3.ch[index] > 8 ){
hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
// hpcIndexVE->Fill( pc.index[i], pc.e[i] );
}
}
}
sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, sx3EUp, sx3EDn);
@ -146,7 +170,7 @@ Bool_t TrackRecon::Process(Long64_t entry){
double theta = -TMath::Pi()/2 + 2*TMath::Pi()/16/4.*(qqq.id[i]*16 + chWedge +0.5);
double rho = 10.+40./16.*(chRing+0.5);
// if(qqq.e[i]>50){
// hqqqPolar->Fill( theta, rho);
hqqqPolar->Fill( theta, rho);
// }
// qqq.used[i] = true;
// qqq.used[j] = true;
@ -163,6 +187,7 @@ Bool_t TrackRecon::Process(Long64_t entry){
}
// //======================= PC
PCHit_1An hitInfo;
ID.clear();
int counter=0;
@ -172,51 +197,33 @@ Bool_t TrackRecon::Process(Long64_t entry){
if( E.size()==3 ){
float aE = 0;
float cE = 0;
int multi_an =0;
for(int l=0;l<E.size();l++){
if(E[l].first<24 && E[l].first!=20 && E[l].first!=12){
multi_an++;
}
}
if(multi_an==1){
bool multi_an =false;
for(int l=0;l<E.size();l++){
if(E[l].first<24 && E[l].first!=20 && E[l].first!=12){
aE = E[l].second;
}else if(E[l].first>24){
cE = E[l].second;
if(!multi_an){
aE = E[l].second;
}
multi_an=true;
}
else {
cE = E[l].second + cE;
}
}
// printf("anode= %d, cathode = %d\n", aID, cID);
// }
if( ID[0].first < 1 ) {
aID = pc.ch[ID[0].second];
cID = pc.ch[ID[1].second];
}else{
cID = pc.ch[ID[0].second];
aID = pc.ch[ID[1].second];
}
//using CalTrack3 to get the track position and direction
// hanVScatsum->Fill(aE,cE);
hanVScatsum->Fill(aE,cE);
if( HitNonZero){
if (ID.size() == 3) {
int aID = -1;
int cID1 = -1;
int cID2 = -1;
for (int i = 0; i < ID.size(); i++) {
if (pc.ch[ID[i].second] < 24 && pc.ch[ID[i].second] != 20 && pc.ch[ID[i].second] != 12) {
aID = pc.ch[ID[i].second];
} else if (pc.ch[ID[i].second] > 24) {
if (cID1 == -1) {
cID1 = pc.ch[ID[i].second];
} else {
cID2 = pc.ch[ID[i].second];
}
}
}
if (aID != -1 && cID1 != -1 && cID2 != -1) {
pw_contr.CalTrack3(hitPos, aID, cID1, cID2);
pw_contr.Print();
printf("###################\n");
hZProj->Fill(pw_contr.GetZ0());
}
}
pw_contr.CalTrack3( hitPos, hitinfo, cID);
hZProj->Fill(pw_contr.GetZ0());
}
// }
@ -233,13 +240,64 @@ Bool_t TrackRecon::Process(Long64_t entry){
return kTRUE;
}
void TrackRecon::Terminate(){
void Analyzer::Terminate(){
gStyle->SetOptStat("neiou");
TCanvas * canvas = new TCanvas("cANASEN", "ANASEN", 200, 200);
padID=1;
canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
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);
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");
}

97
centroids.txt Normal file
View File

@ -0,0 +1,97 @@
HistogramIndex PeakNumber Centroid Amplitude Sigma
0 1 991.118
0 2 2026.83
0 3 3060.26
0 4 4092.45
1 1 922.213
1 2 1885.55
1 3 2845.53
1 4 3810.32
2 1 955.591
2 2 1953.17
2 3 2949.37
2 4 3950.79
3 1 995.787
3 2 2036.58
3 3 3076.91
3 4 4112.05
4 1 1017.48
4 2 2080.19
4 3 3142.24
4 4 4206.1
5 1 1022.78
5 2 2091.21
5 3 3158.28
5 4 4226.97
6 1 1076.22
6 2 2203.37
6 3 3329.53
6 4 4457.69
7 1 977.46
7 2 1998.02
7 3 3017.36
7 4 4040.47
8 1 1049.74
8 2 2144.38
8 3 3238.2
8 4 4335.25
9 1 1000.59
9 2 2046.42
9 3 3090.29
9 4 4129.63
10 1 1014.92
10 2 2076.16
10 3 3134.59
10 4 4213.42
11 1 1004.85
11 2 2052.88
11 3 3100.3
11 4 4164.75
12 1 945.861
12 2 1932.49
12 3 2917.95
12 4 3955.15
13 1 998.307
13 2 2040.38
13 3 3078.76
13 4 4135.51
14 1 966.429
14 2 1972.15
14 3 2974.84
14 4 4056.41
15 1 958.352
15 2 1958.64
15 3 2957.7
15 4 3970.41
16 1 970.732
16 2 1977.63
16 3 2984.97
16 4 4002.56
17 1 1013.65
17 2 2064.9
17 3 3114.19
17 4 4190.98
18 1 975.538
18 2 1990.64
18 3 3005.46
18 4 4048.99
19 1 1082.91
19 2 2194.08
19 3 3303.65
19 4 4411.32
20 1 912.778
20 2 1866.83
20 3 2819.21
20 4 3781.63
21 1 1002.36
21 2 1989.95
21 3 2975.53
21 4 3986.71
22 1 1075.38
22 2 2144.25
22 3 3210.17
22 4 4312.84
23 1 988.828
23 2 2016.35
23 3 3044.19
23 4 4082.41

89
centroids_edited.txt Normal file
View File

@ -0,0 +1,89 @@
HistogramIndex PeakNumber Centroid Amplitude Sigma
1 1 922.213
1 2 1885.55
1 3 2845.53
1 4 3810.32
2 1 955.591
2 2 1953.17
2 3 2949.37
2 4 3950.79
3 1 995.787
3 2 2036.58
3 3 3076.91
3 4 4112.05
4 1 1017.48
4 2 2080.19
4 3 3142.24
4 4 4206.1
5 1 1022.78
5 2 2091.21
5 3 3158.28
5 4 4226.97
6 1 1076.22
6 2 2203.37
6 3 3329.53
6 4 4457.69
7 1 977.46
7 2 1998.02
7 3 3017.36
7 4 4040.47
8 1 1049.74
8 2 2144.38
8 3 3238.2
8 4 4335.25
9 1 1000.59
9 2 2046.42
9 3 3090.29
9 4 4129.63
10 1 1014.92
10 2 2076.16
10 3 3134.59
10 4 4213.42
11 1 1004.85
11 2 2052.88
11 3 3100.3
11 4 4164.75
12 1 945.861
12 2 1932.49
12 3 2917.95
12 4 3955.15
13 1 998.307
13 2 2040.38
13 3 3078.76
13 4 4135.51
14 1 966.429
14 2 1972.15
14 3 2974.84
14 4 4056.41
15 1 958.352
15 2 1958.64
15 3 2957.7
15 4 3970.41
16 1 970.732
16 2 1977.63
16 3 2984.97
16 4 4002.56
17 1 1013.65
17 2 2064.9
17 3 3114.19
17 4 4190.98
18 1 975.538
18 2 1990.64
18 3 3005.46
18 4 4048.99
20 1 912.778
20 2 1866.83
20 3 2819.21
20 4 3781.63
21 1 1002.36
21 2 1989.95
21 3 2975.53
21 4 3986.71
22 1 1075.38
22 2 2144.25
22 3 3210.17
22 4 4312.84
23 1 988.828
23 2 2016.35
23 3 3044.19
23 4 4082.41

541
gainmatch.C
View File

@ -1,541 +0,0 @@
#define gainmatch_cxx
#include "gainmatch.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TCutG.h>
#include <utility>
#include <algorithm>
#include "Armory/ClassSX3.h"
#include "Armory/ClassPW.h"
#include "TVector3.h"
TH2F * hsx3IndexVE;
TH2F * hqqqIndexVE;
TH2F * hqqqIndexVE_cut;
TH2F * hpcIndexVE;
TH2F * hsx3Coin;
TH2F * hqqqCoin;
TH2F * hpcCoin;
TH2F * hpcCoin_cut;
TH2F * hGoodQQQ;
TH2F * hGoodQQQRingVWedge;
TH2F * hqqqPolar;
TH2F * hsx3VpcIndex;
TH2F * hqqqVpcIndex;
TH2F * hqqqVpcIndex_cut;
TH2F * hqqqVpcE;
TH2F * hqqqVpcE_cut;
TH2F * hqqqVpcE_cut1;
TH2F * hqqqVpcE_cut2;
TH2F * hqqqVpcE_cutCoinc;
TH2F * hsx3VpcE;
TH2F * hanVScatsum;
TH2F * hanVScatsum_cut;
TH2F * hanVScatsum_cut1;
TH2F * hanVScatsum_cut2;
TH2F * hsx3Vsx3;
TH2F * hsx3uVsx3d_01;
TH2F * hsx3uVsx3d_23;
TH2F * hsx3uVsx3d_45;
TH2F * hsx3uVsx3d_67;
TH2F * hVCID;
TH1F *hsx3bk_9_shifted ;
TH1F *hsx3bk_10_shifted ;
TH1F *hsx3bk_11_shifted ;
int padID = 0;
TCutG *Coinc_cut_set1;
//TCutG *crap_cut;
TCutG *AnCathCoinc_cut;
TCutG *AnCathCoinc_cut1;
TCutG *AnCathCoinc_cut2;
SX3 sx3_contr;
PW pw_contr;
TVector3 hitPos;
bool HitNonZero;
bool inCut;
bool inCut1;
bool inCut2;
bool inCutCoinc;
TH1F *hZd_01_1;
TH1F *hZd_01_2;
TH1F *hZd_01_3;
TH1F *hZd_01_4;
TH1F * hZProj;
TH1F * hsx3bk_11;
TH1F * hsx3bk_10;
TH1F * hsx3bk_9;
TH1F * hsx3bk_8;
void gainmatch::Begin(TTree * /*tree*/){
TString option = GetOption();
hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24*12, 0, 24*12, 400, 0, 5000); hsx3IndexVE->SetNdivisions( -612, "x");
hqqqIndexVE = new TH2F("hqqqIndexVE", "QQQ index vs Energy; QQQ index ; Energy", 4*2*16, 0, 4*2*16, 400, 0, 5000); hqqqIndexVE->SetNdivisions( -1204, "x");
hqqqIndexVE_cut = new TH2F("hqqqIndexVE_cut", "QQQ index vs Energy gated; QQQ index ; Energy", 4*2*16, 0, 4*2*16, 400, 0, 5000); hqqqIndexVE->SetNdivisions( -1204, "x");
hpcIndexVE = new TH2F("hpcIndexVE", "PC index vs Energy; PC index ; Energy", 2*24, 0, 2*24, 400, 0, 4000); hpcIndexVE->SetNdivisions( -1204, "x");
hGoodQQQ = new TH2F("hGoodQQQ", "number of good QQQ vs QQQ id", 10, 0, 10, 4, 0, 4);
hGoodQQQRingVWedge = new TH2F("hGoodQQQRingVWedge", "Ring index, Wedge index", 16*4, 0, 16*4, 16*4, 0, 16*4);
hZd_01_1 =new TH1F("hZd_01_1", "Z position", 100, -1, 1);
hZd_01_2 =new TH1F("hZd_01_2", "Z position", 100, -1, 1);
hZd_01_3 =new TH1F("hZd_01_3", "Z position", 100, -1, 1);
hZd_01_4 =new TH1F("hZd_01_4", "Z position", 100, -1, 1);
hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24*12, 0, 24*12, 24*12, 0, 24*12);
hqqqCoin = new TH2F("hqqqCoin", "QQQ Coincident", 4*2*16, 0, 4*2*16, 4*2*16, 0, 4*2*16);
hpcCoin = new TH2F("hpcCoin", "PC Coincident", 2*24, 0, 2*24, 2*24, 0, 2*24);
hpcCoin_cut = new TH2F("hpcCoin_cut", "PC Coincident gated", 2*24, 0, 2*24, 2*24, 0, 2*24);
hqqqPolar = new TH2F("hqqqPolar", "QQQ Polar ID", 16*4, -TMath::Pi(), TMath::Pi(),16, 10, 50);
hsx3VpcIndex = new TH2F("hsx3Vpcindex", "sx3 vs pc; sx3 index; pc index", 24*12, 0, 24*12, 48, 0, 48);
hsx3Vsx3 = new TH2F("hsx3Vsx3", "sx3 vs sx3; sx3 E; sx3 E", 8000, 0, 16000, 8000, 0, 16000);
hsx3uVsx3d_01 = new TH2F("hsx3uVsx3d_01", "sx3u vs sx3d; sx3u E; sx3d E", 100, 0, 1, 100, 0, 1);
hsx3uVsx3d_23 = new TH2F("hsx3uVsx3d_23", "sx3u vs sx3d; sx3u E; sx3d E", 100, 0, 1, 100, 0, 1);
hsx3uVsx3d_45 = new TH2F("hsx3uVsx3d_45", "sx3u vs sx3d; sx3u E; sx3d E", 1000, 0, 5000, 1000, 0, 5000);
hsx3uVsx3d_67 = new TH2F("hsx3uVsx3d_67", "sx3u vs sx3d; sx3u E; sx3d E", 1000, 0, 5000, 1000, 0, 5000);
hsx3VpcIndex->SetNdivisions( -612, "x");
hsx3VpcIndex->SetNdivisions( -12, "y");
hqqqVpcIndex = new TH2F("hqqqVpcindex", "qqq vs pc; qqq index; pc index", 4*2*16, 0, 4*2*16, 48, 0, 48);
hqqqVpcIndex->SetNdivisions( -612, "x");
hqqqVpcIndex->SetNdivisions( -12, "y");
hqqqVpcIndex_cut = new TH2F("hqqqVpcindex_cut", "qqq vs pc gated; qqq index; pc index", 4*2*16, 0, 4*2*16, 48, 0, 48);
hqqqVpcIndex_cut->SetNdivisions( -612, "x");
hqqqVpcIndex_cut->SetNdivisions( -12, "y");
hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 8000, 0, 16000, 8000, 0, 16000);
hqqqVpcE->SetNdivisions( -612, "x");
hqqqVpcE->SetNdivisions( -12, "y");
hqqqVpcE_cut = new TH2F("hqqqVpcEnergy_cut", "qqq vs pc gated; qqq energy; pc energy", 8000, 0, 16000, 8000, 0, 16000);
hqqqVpcE_cut->SetNdivisions( -612, "x");
hqqqVpcE_cut->SetNdivisions( -12, "y");
hqqqVpcE_cut1 = new TH2F("hqqqVpcEnergy_cut1", "qqq vs pc gated; qqq energy; pc energy", 8000, 0, 16000, 8000, 0, 16000);
hqqqVpcE_cut1->SetNdivisions( -612, "x");
hqqqVpcE_cut1->SetNdivisions( -12, "y");
hqqqVpcE_cut2 = new TH2F("hqqqVpcEnergy_cut2", "qqq vs pc gated; qqq energy; pc energy", 8000, 0, 16000, 8000, 0, 16000);
hqqqVpcE_cut2->SetNdivisions( -612, "x");
hqqqVpcE_cut2->SetNdivisions( -12, "y");
hqqqVpcE_cutCoinc = new TH2F("hqqqVpcEnergy_cutCoinc", "qqq vs pc gated; qqq energy; pc energy", 8000, 0, 16000, 8000, 0, 16000);
hqqqVpcE_cutCoinc->SetNdivisions( -612, "x");
hqqqVpcE_cutCoinc->SetNdivisions( -12, "y");
hsx3bk_8=new TH1F("hsx3bk_8", "hsx3bk_8",1000, 0,5000);
hsx3bk_9=new TH1F("hsx3bk_9", "hsx3bk_9",1000, 0,5000);
hsx3bk_10=new TH1F("hsx3bk_10", "hsx3bk_10",1000, 0,5000);
hsx3bk_11=new TH1F("hsx3bk_11", "hsx3bk_11",1000, 0,5000);
hsx3VpcE = new TH2F("hsx3VpcEnergy", "sx3 vs pc; sx3 energy; pc energy", 400, 0, 5000, 400, 0, 5000);
hsx3VpcE->SetNdivisions( -612, "x");
hsx3VpcE->SetNdivisions( -12, "y");
hsx3bk_9_shifted = new TH1F("hsx3bk_9_shifted", "hsx3bk_9",1000, 0,5000);
hsx3bk_10_shifted = new TH1F("hsx3bk_10_shifted", "hsx3bk_9",1000, 0,5000);
hsx3bk_11_shifted = new TH1F("hsx3bk_11_shifted", "hsx3bk_9",1000, 0,5000);
hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600);
hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 8000,0 , 16000, 8000, 0 , 16000);
hanVScatsum_cut = new TH2F("hanVScatsum_cut", "Anode vs Cathode Sum gated; Anode E; Cathode E", 1600,0 , 16000, 1600, 0 , 16000);
hanVScatsum_cut1 = new TH2F("hanVScatsum_cut1", "Anode vs Cathode Sum gated; Anode E; Cathode E", 1600,0 , 16000, 1600, 0 , 16000);
hanVScatsum_cut2 = new TH2F("hanVScatsum_cut2", "Anode vs Cathode Sum gated; Anode E; Cathode E", 1600,0 , 16000, 1600, 0 , 16000);
hVCID = new TH2F("hVCID", "Virtual Cathod ID vs total Cath. Energy", 200, 0, 24, 200, 0, 10000);
sx3_contr.ConstructGeo();
pw_contr.ConstructGeo();
TFile *f3 = new TFile("Coinc_cut_set1.root");
//TFile *f4 = new TFile("crap_cut.root");
TFile *f = new TFile("AnCathCoinc_cut.root");
TFile *f1 = new TFile("AnCathCoinc_cut1.root");
TFile *f2 = new TFile("AnCathCoinc_cut2.root");
Coinc_cut_set1 = (TCutG*)f3->Get("Coinc_cut_set1");
//crap_cut = (TCutG*)f4->Get("crap_cut");
AnCathCoinc_cut = (TCutG*)f->Get("AnCathCoinc_cut");
AnCathCoinc_cut1 = (TCutG*)f1->Get("AnCathCoinc_cut1");
AnCathCoinc_cut2 = (TCutG*)f2->Get("AnCathCoinc_cut2");
}
Bool_t gainmatch::Process(Long64_t entry){
// if ( entry > 100 ) return kTRUE;
hitPos.Clear();
HitNonZero = false;
inCut = false;
// if( entry > 1) return kTRUE;
// printf("################### ev : %llu \n", entry);
b_sx3Multi->GetEntry(entry);
b_sx3ID->GetEntry(entry);
b_sx3Ch->GetEntry(entry);
b_sx3E->GetEntry(entry);
b_sx3T->GetEntry(entry);
b_qqqMulti->GetEntry(entry);
b_qqqID->GetEntry(entry);
b_qqqCh->GetEntry(entry);
b_qqqE->GetEntry(entry);
b_qqqT->GetEntry(entry);
b_pcMulti->GetEntry(entry);
b_pcID->GetEntry(entry);
b_pcCh->GetEntry(entry);
b_pcE->GetEntry(entry);
b_pcT->GetEntry(entry);
sx3.CalIndex();
qqq.CalIndex();
pc.CalIndex();
// sx3.Print();
//########################################################### Raw data
// //======================= SX3
std::vector<std::pair<int, int>> ID; // first = id, 2nd = index
for (int i = 0; i < sx3.multi; i++) {
ID.push_back(std::pair<int, int>(sx3.id[i], i));
hsx3IndexVE->Fill(sx3.index[i], sx3.e[i]);
for (int j = i + 1; j < sx3.multi; j++) {
hsx3Coin->Fill(sx3.index[i], sx3.index[j]);
}
for (int j = 0; j < pc.multi; j++) {
hsx3VpcIndex->Fill(sx3.index[i], pc.index[j]);
}
}
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;
});
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) {
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;
if (sx3.ch[index] < 8) {
if (sx3.ch[index] % 2 == 0) {
sx3ChDn = sx3.ch[index];
sx3EDn = sx3.e[index];
} else {
sx3ChUp = sx3.ch[index];
sx3EUp = sx3.e[index];
}
} else {
sx3ChBk = sx3.ch[index];
sx3EBk = sx3.e[index];
}
int ch = sx3.ch[index];
float energy = sx3.e[index];
if (sx3ID[0].first == 9) {
float peak8 = 0.0;
float peak9 = 0.0;
int peak10 = 0.0;
float peak11 = 0.0;
float shift9 =0.0;
float shift10 =0.0;
float shift11 =0.0;
int minBin_8 = hsx3bk_8->FindBin(1);
int maxBin_8 = hsx3bk_8->FindBin(5000);
int maxRangeBinContent_8 = -1;
double maxBinCenter_8 = 0.0;
int minBin_9 = hsx3bk_9->FindBin(1);
int maxBin_9 = hsx3bk_9->FindBin(5000);
int maxRangeBinContent_9 = -1;
double maxBinCenter_9 = 0.0;
int minBin_10 = hsx3bk_10->FindBin(1);
int maxBin_10 = hsx3bk_10->FindBin(5000);
int maxRangeBinContent_10 = -1;
double maxBinCenter_10 = 0.0;
int minBin_11 = hsx3bk_11->FindBin(1);
int maxBin_11 = hsx3bk_11->FindBin(5000);
int maxRangeBinContent_11 = -1;
double maxBinCenter_11 = 0.0;
if (sx3ChBk == 8) {
hsx3bk_8->Fill(sx3EBk);
for (int bin = minBin_8; bin <= maxBin_8; ++bin) {
if (hsx3bk_8->GetBinContent(bin) > maxRangeBinContent_8) {
maxRangeBinContent_8 = hsx3bk_8->GetBinContent(bin);
maxBinCenter_8 = hsx3bk_8->GetBinCenter(bin);
}
}
//peak8 = hsx3bk_8->GetMaximumBin();
//peak8 = hsx3bk_8->GetMaximumBin();
//printf("peak8: %f\n", maxBinCenter_8);
}
//printf("peak8_mm: %f\n", maxBinCenter);
else if (sx3ChBk == 9) {
hsx3bk_9->Fill(sx3EBk);
for (int bin = minBin_9; bin <= maxBin_9; ++bin) {
if (hsx3bk_9->GetBinContent(bin) > maxRangeBinContent_9) {
maxRangeBinContent_9 = hsx3bk_9->GetBinContent(bin);
maxBinCenter_9 = hsx3bk_9->GetBinCenter(bin);
}
}
//peak8 = hsx3bk_8->GetMaximumBin();
//peak8 = hsx3bk_8->GetMaximumBin();
//printf("peak9: %f\n", maxBinCenter_9);
//hsx3bk_9_shifted->Fill(sx3EBk*0.76);
peak9 = 2097.5/maxBinCenter_9;
//printf("peak9_shift: %f\n", peak9);
hsx3bk_9_shifted->Fill(sx3EBk*(2097.5/maxBinCenter_9));
//printf("peak9 %d\n", peak9);
}
else if(sx3ChBk == 10) {
hsx3bk_10->Fill(sx3EBk);
for (int bin = minBin_10; bin <= maxBin_10; ++bin) {
if (hsx3bk_10->GetBinContent(bin) > maxRangeBinContent_10) {
maxRangeBinContent_10 = hsx3bk_10->GetBinContent(bin);
maxBinCenter_10 = hsx3bk_10->GetBinCenter(bin);
}
}
//peak8 = hsx3bk_8->GetMaximumBin();
//peak8 = hsx3bk_8->GetMaximumBin();
//printf("peak10: %f\n", maxBinCenter_10);
//hsx3bk_9_shifted->Fill(sx3EBk*0.76);
peak10= 2097.5/maxBinCenter_10;
//printf("peak10_shift: %f\n", 1787.5/maxBinCenter_10);
hsx3bk_10_shifted->Fill(sx3EBk*(2097.5/maxBinCenter_10));
//printf("peak9 %d\n", peak9);
}
//peak10 = hsx3bk_10->GetMaximumBin();
// printf("peak10 %d\n" ,peak10);
else if(sx3ChBk == 11) {
hsx3bk_11->Fill(sx3EBk);
for (int bin = minBin_11; bin <= maxBin_11; ++bin) {
if (hsx3bk_11->GetBinContent(bin) > maxRangeBinContent_11) {
maxRangeBinContent_11 = hsx3bk_11->GetBinContent(bin);
maxBinCenter_11 = hsx3bk_11->GetBinCenter(bin);
}
}
//peak8 = hsx3bk_8->GetMaximumBin();
//peak8 = hsx3bk_8->GetMaximumBin();
//printf("peak9: %f\n", maxBinCenter_11);
//hsx3bk_9_shifted->Fill(sx3EBk*0.76);
peak11 = 2097.5/maxBinCenter_11;
//printf("peak11_shift: %f\n", peak11);
hsx3bk_11_shifted->Fill(sx3EBk*(2097.5/maxBinCenter_11));
//printf("peak9 %d\n", peak9);
}
float sx3EBk_shifted = 0.0;
float sx3E_u_matched_01 = 0.0;
float sx3E_d_matched_01 = 0.0;
float sx3E_fb_matched_01 = 0.0;
float sx3E_fbu_matched_01 = 0.0;
float sx3E_fbd_matched_01 = 0.0;
float diff =0.0;
float ratio = 0.0;
float coeff = 0.0;
if (sx3ChBk == 9) {
sx3EBk_shifted = (sx3EBk *(2097.5/maxBinCenter_9));
} else if (sx3ChBk == 10) {
sx3EBk_shifted = (sx3EBk * (2097.5/maxBinCenter_10));
} else if (sx3ChBk == 11) {
sx3EBk_shifted = (sx3EBk * (2097.5/maxBinCenter_11)) ;
} else {
sx3EBk_shifted = sx3EBk; // Use unshifted value for sx3ChBk == 8
}
if ((sx3ChUp == 1 && sx3ChDn == 0)) {
sx3E_u_matched_01= (sx3EUp-0.898729)/0.836243;
//sx3E_u_matched_01= (0.836243*sx3EDn)+0.898729;
sx3E_d_matched_01= (sx3EDn-0.898729)/0.836243;
sx3E_fb_matched_01=(sx3EBk_shifted+9.2423)/0.924773 ;
sx3E_fbu_matched_01=(sx3E_u_matched_01+9.2423)/0.924773 ;
sx3E_fbd_matched_01=(sx3E_d_matched_01+9.2423)/0.924773 ;
diff = sx3E_fb_matched_01 - (sx3EUp+sx3E_fbd_matched_01);
ratio = sx3EUp/sx3E_fbd_matched_01;
coeff = ((sx3EUp+diff) - (sx3E_fbd_matched_01*ratio))/(diff*(1+ratio));
}
//TH2F *hsx3uVsx3d_01 = nullptr;
if (sx3ChBk >=8) {
//if (sx3ChBk == 9) {
if ((sx3ChUp == 1 && sx3ChDn == 0)) {
if (sx3ChUp != -1 && sx3ChDn != -1 && sx3ChBk !=-1) {
if (sx3EBk_shifted > 50 && sx3EUp > 50 && sx3EDn > 50) {
printf("sx3EUp: %f, sx3EDn: %f, sx3E_u_matched_01: %f,sx3E_d_matched_01: %f\n", sx3EUp, sx3EDn, sx3E_u_matched_01,sx3E_d_matched_01);
//printf("Filling hsx3uVsx3d_01_shifted: %f\n", sx3EBk_ud_matched_01 / sx3EBk_shifted);
// hsx3uVsx3d_01->Fill(sx3E_u_matched_01 / sx3EBk_shifted, sx3E_d_matched_01 / sx3EBk_shifted);
hsx3uVsx3d_01->Fill(sx3EUp / sx3EBk_shifted, sx3E_d_matched_01 / sx3EBk_shifted);
hsx3uVsx3d_23->Fill(sx3EUp / sx3EBk_shifted, sx3EDn/ sx3EBk_shifted);
}
}
//}
}
else if ((sx3ChUp == 3 && sx3ChDn == 2)) {
if (sx3ChUp != -1 && sx3ChDn != -1 && sx3ChBk !=-1) {
if (sx3EBk_shifted != 0 && sx3EBk_shifted > 50 && sx3EUp > 50 && sx3EDn > 50) {
printf("sx3EUp: %f, sx3EDn: %f, sx3EBk_shifted: %f\n", sx3EUp, sx3EDn, sx3EBk_shifted);
printf("Filling hsx3uVsx3d_23_shifted: %f\n", sx3EUp / sx3EBk_shifted);
// hsx3uVsx3d_23->Fill(sx3EUp / sx3EBk_shifted, (-0.924754*sx3EUp+0.916671) / sx3EBk_shifted);
}
}
}
/* if ((sx3ChUp == 1 && sx3ChDn == 0)) {
if (sx3ChUp != -1 && sx3ChDn != -1 && sx3ChBk !=-1) {
if (sx3EBk != 0 && sx3EBk > 50 && sx3EUp > 50 && sx3EDn > 50) {
printf("sx3EUp: %f, sx3EDn: %f, sx3EBk: %f\n", sx3EUp, sx3EDn, sx3EBk);
printf("Filling hsx3uVsx3d_01: %f\n", sx3EUp / sx3EBk);
hsx3uVsx3d_45->Fill(sx3EUp / sx3EBk, sx3EDn / sx3EBk);
}
}
}
else if ((sx3ChUp == 3 && sx3ChDn == 2)) {
if (sx3ChUp != -1 && sx3ChDn != -1 && sx3ChBk !=-1) {
if (sx3EBk != 0 && sx3EBk > 50 && sx3EUp > 50 && sx3EDn > 50) {
printf("sx3EUp: %f, sx3EDn: %f, sx3EBk: %f\n", sx3EUp, sx3EDn, sx3EBk);
printf("Filling hsx3uVsx3d_23: %f\n", sx3EUp / sx3EBk);
hsx3uVsx3d_67->Fill(sx3EUp / sx3EBk, sx3EDn / sx3EBk);
}
}
}*/
if (sx3ChUp == 1 && sx3ChDn == 0){
//if (sx3ChUp == 1 || sx3ChDn == 0 || sx3ChUp == 3 || sx3ChDn == 2 || sx3ChUp == 5 || sx3ChDn == 4 || sx3ChUp == 7 || sx3ChDn == 6) {
if (sx3ChUp != -1 && sx3ChBk !=-1 && sx3ChDn !=-1) {
if (sx3EBk_shifted > 50 && sx3EUp > 50 && sx3EDn>50 &&sx3E_u_matched_01>50 && sx3E_u_matched_01>50) {
//printf("sx3EUp: %f, sx3EDn: %f, sx3E_u_matched_01: %f,sx3E_d_matched_01: %f\n", sx3EUp, sx3EDn, sx3E_u_matched_01,sx3E_d_matched_01);
printf("Filling hsx3uVsx3d_nn: %f, gggggg: %f \n", (sx3EUp+sx3EDn),(sx3E_u_matched_01+sx3E_d_matched_01) );
hsx3uVsx3d_45->Fill((sx3EUp+sx3E_d_matched_01),sx3EBk_shifted);
hsx3uVsx3d_67->Fill((sx3EUp+sx3E_d_matched_01),sx3E_fb_matched_01);
}
}
}
/*if (sx3ChBk > 8) {
if ((sx3ChUp == 7 && sx3ChDn == 6) ||
(sx3ChUp == 5 && sx3ChDn == 4) ||
(sx3ChUp == 3 && sx3ChDn == 2) ||
(sx3ChUp == 1 && sx3ChDn == 0)) {
if (sx3ChUp != -1 && sx3ChDn != -1 && sx3ChBk !=-1) {
if (sx3EBk != 0 && sx3EBk > 50 && sx3EUp > 50 && sx3EDn > 50) {
hsx3uVsx3d->Fill(sx3EUp / sx3EBk, sx3EDn / sx3EBk);
hsx3Vsx3->Fill(sx3EUp ,sx3EDn);
printf("sx3EUp: %f | sx3EDn: %f | sx3EBk: %f | sx3ChUp: %d | sx3ChDn: %d | sx3ChBk: %d\n", sx3EUp, sx3EDn, sx3EBk, sx3ChUp, sx3ChDn, sx3ChBk);
}
}
}
}*/
//else {
//printf("sx3EUp\n");
//}
if (sx3ChUp == 1 && sx3ChDn == 0){
if (sx3ChUp != -1 && sx3ChBk !=-1 && sx3ChDn !=-1) {
if (sx3E_d_matched_01> sx3EUp ) {
//printf("hZd_01_1_dn: %f\n", sx3E_d_matched_01);
//printf("hZd_01_1_b: %f\n", sx3E_fb_matched_01);
hZd_01_1->Fill((2*(sx3E_d_matched_01+(coeff*diff))/sx3E_fb_matched_01)-1);
}
else if(sx3EUp> sx3E_d_matched_01) {
//printf("hZd_01_2_sx3EUp: %f\n",sx3EUp );
//printf("hZd_01_2_sx3EDn: %f\n",sx3E_fb_matched_01);
hZd_01_2->Fill(1-(2*(sx3EUp+(1-coeff)*diff))/sx3E_fb_matched_01);
}
else if(sx3EUp>0.0 && sx3E_d_matched_01>0.0 && sx3E_d_matched_01>=sx3EUp ) {
hZd_01_3->Fill((2*(sx3E_d_matched_01+ coeff*diff)/sx3E_fb_matched_01)-1);
}
else if(sx3EUp>0.0 && sx3E_d_matched_01>0.0 && sx3E_d_matched_01<sx3EUp ) {
hZd_01_4->Fill(1-(2*(sx3EUp+ (1-coeff)*diff)/sx3E_fb_matched_01));
}
}
}
}
for (int j = 0; j < pc.multi; j++) {
if (sx3.ch[index] > 8) {
hsx3VpcE->Fill(sx3.e[i], pc.e[j]);
}
}
}
sx3_contr.CalSX3Pos(sx3ID[0].first, sx3ChUp, sx3ChDn, sx3ChBk, sx3EUp, sx3EDn);
hitPos = sx3_contr.GetHitPos();
HitNonZero = true;
// hitPos.Print();
}
}
// //======================= PC
//########################################################### Track constrcution
//############################## DO THE KINEMATICS
}
return kTRUE;
}
void gainmatch::Terminate(){
}

26
mapping.h
View File

@ -19,21 +19,23 @@ const std::map<int, unsigned short> board = {
{4, 22129},
{5, 15529},
{6, 15528},
// {7,89},
{7, 334},
{8, 379},
{9, 325},
{10, 405}
{10, 405}
};
const int nBd = board.size();
const int nV1740 = 7;
const int nV1725 = 3;
const int nV1725 = 4;
//+++++++++++++++++++ detID;
// The detectors are seperated into 2 type: SuperX3, QQQ, and PC
// the SuperX3 has 24 detectors for each kind, wach detector has 12 channels
// the QQQ has 4 detectors for each kind, each detector has 32 channels
// the PC has 2 types, anode and cathode, each has 24 channels
// the MISC has 6 channels, the lollipop IC and siliscon followed by the hotneedle IC, as well as the Rf and MCP
// The detID = Type * 10000 + index * 100 + channel
// fro example, detID(superX3-8, ch-7) = 00807
@ -76,17 +78,23 @@ const std::vector<int> mapping = {
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
//================== 89
// 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
// 30004, -1, 30003, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1,
//================== 334
20116, 20117, 20118, 20119, 20120, 20121, 20122, 20123, 20016, 20017, 20018, 20019, 20020, 20021, 20022, 20023,
20116, 20117, 20118, 20119, -1, 20121, 20122, 20123, 20016, 20017, 20018, -1, 20020, 20021, 20022, 20023,
//================== 379
20000, 20001, 20002, 20003, 20004, 20005, -1, 20007, 20008, -1, 20010, 20011, 20012, 20013, 20014, 20015,
-1 , 20001, 20002, 20003, 20004, 20005, -1, 20007, 20008, -1, 20010, 20011, 20012, 20013, 20014, 20015,
//================== 325
20100, 20101, 20102, 20103, 20104, 20105, 20106, 20107, 20108, 20109, 20110, 20111, 20112, 20113, 20114, 20115,
20100, 20101, 20102, 20103, 20104, 20105, 20106, 20107, 20108, 20109, 20110, 20111, 20112, -1, 20114, 20115,
//================== 405
20006, -1, -1, 20009, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1
// 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
20006, -1, 30005, 20009, -1, 20120, 20000, 20019, 20113, 30000, 30004, 30001, 30002, -1, 30003, -1
};
//MCP moved from channel 1 to 2 after Run number 322
//MCP and Rf moved to ch 0 and 1 after Run number after Run282
//moved back to ch
void PrintMapping(){
int digiID = 0;
@ -141,8 +149,9 @@ void PrintMapping(){
printf("\033[35m%3d(%2d)\033[0m|", detID, ch);
}else{
}else if( typeID == 3){ // MISC
printf("\033[33m%3d(%2d)\033[0m|", detID, ch);
}
}
@ -215,7 +224,6 @@ void GenMapping(std::string mapFile){
detID += 20000;
if( words[3] == "ANODE") detID += atoi(words[4].c_str());
if( words[3] == "CATHODE") detID += 100 + atoi(words[4].c_str());
}
if( words[2] == "blank") {

23
slope_intercept_cathode.txt Normal file
View File

@ -0,0 +1,23 @@
Histogram Number Slope Intercept
24 1 -2.89219e-10
25 0.942098 -0.105169
26 0.980862 -0.732032
27 0.982975 -2.22704
28 0.978815 -1.51477
29 0.965245 -2.19515
30 0.945384 -0.892599
31 0.977408 -0.908592
32 0.919546 3.25464
33 0.972194 2.44956
34 0.92852 5.44745
35 0.947098 1.40531
36 0.875491 -1.13145
37 1.95496 -1735.58
38 0.970862 2.86019
40 0.91793 -3.80615
41 0.913897 -2.12964
42 0.954014 -0.760604
43 0.993616 -1.40278
45 0.926169 -21.2016
46 1.00577 -2.14281
47 0.943312 -1.26464

49
slope_intercept_results.txt Normal file
View File

@ -0,0 +1,49 @@
Histogram Number Slope Intercept
0 0.931015 -1.35431
1 1 -1.87356e-10
2 0.964185 1.49989
3 0.92638 -1.30621
4 0.905569 1.00834
5 0.901182 0.470903
6 0.853932 3.32687
7 0.942785 1.08887
8 0.878904 -0.0107433
9 0.922662 -2.32259
10 0.903343 8.38332
11 0.914227 6.56108
12 0.961008 23.0982
13 0.920976 5.22104
14 0.936584 31.5073
15 0.959044 5.43267
16 0.95263 -0.404053
17 0.90953 4.82833
18 0.940277 10.3629
19 0.86746 -17.8678
20 1.00683 4.76371
21 0.968342 -43.9496
22 0.892882 -32.0742
23 0.933615 1.10704
24 1 -2.89219e-10
25 0.942098 -0.105169
26 0.980862 -0.732032
27 0.982975 -2.22704
28 0.978815 -1.51477
29 0.965245 -2.19515
30 0.945384 -0.892599
31 0.977408 -0.908592
32 0.919546 3.25464
33 0.972194 2.44956
34 0.92852 5.44745
35 0.947098 1.40531
36 0.875491 -1.13145
37 1 0
38 0.970862 2.86019
39 1 0
40 0.91793 -3.80615
41 0.913897 -2.12964
42 0.954014 -0.760604
43 0.993616 -1.40278
44 1 0
45 0.926169 -21.2016
46 1.00577 -2.14281
47 0.943312 -1.26464

21
slope_intercept_results_anode.txt Normal file
View File

@ -0,0 +1,21 @@
Histogram Number Slope Intercept
1 1 -1.87356e-10
2 0.964185 1.49989
3 0.92638 -1.30621
4 0.905569 1.00834
5 0.901182 0.470903
7 0.942785 1.08887
8 0.878904 -0.0107433
10 0.903343 8.38332
11 0.914227 6.56108
12 0.961008 23.0982
13 0.920976 5.22104
14 0.936584 31.5073
15 0.959044 5.43267
16 0.95263 -0.404053
17 0.90953 4.82833
18 0.940277 10.3629
20 1.00683 4.76371
21 0.968342 -43.9496
22 0.892882 -32.0742
23 0.933615 1.10704