Compare commits
37 Commits
master
...
devel_vign
Author | SHA1 | Date | |
---|---|---|---|
|
e86ab5ed4d | ||
|
9cadfdd191 | ||
|
d4582d80ff | ||
|
18870ed82b | ||
|
5c1c5348f4 | ||
|
65ab69ebe6 | ||
|
1df7470ca1 | ||
|
a7a765c059 | ||
|
39e8f41ab1 | ||
|
9225620426 | ||
|
2225b3a942 | ||
|
6cfb38b564 | ||
|
cb72c14ca4 | ||
|
5995081396 | ||
|
d623e0cd17 | ||
|
8d7322cf5a | ||
|
699b0f8701 | ||
|
02213caaee | ||
|
56a6389b4f | ||
|
b99ad4e4d7 | ||
|
a5dfa2ecd3 | ||
|
26e943adc8 | ||
|
42e093b104 | ||
|
3129339647 | ||
|
7a70340b18 | ||
|
a10081ea81 | ||
|
4cb8a2c48c | ||
|
fe6dbee171 | ||
|
7805481ead | ||
|
511b4aa808 | ||
|
43233ceb02 | ||
|
68fc36a8f6 | ||
|
a6e754b958 | ||
|
48ede97992 | ||
|
f0a393abe2 | ||
|
4ba9c73b98 | ||
|
238ec8961e |
5
.gitignore
vendored
5
.gitignore
vendored
|
@ -4,9 +4,12 @@ EventBuilder*
|
|||
*.pcm
|
||||
*.root
|
||||
*.exe
|
||||
|
||||
*.txt
|
||||
Mapper
|
||||
AnasenMS
|
||||
|
||||
data/
|
||||
data_proton/
|
||||
Analyzer_C_ACLiC_dict0713aaa966_dictContent.h
|
||||
.gitignore
|
||||
Analyzer_C_ACLiC_dict5411fecd5c_dictUmbrella.h
|
||||
|
|
46
.vscode/c_cpp_properties.json
vendored
46
.vscode/c_cpp_properties.json
vendored
|
@ -1,23 +1,11 @@
|
|||
{
|
||||
"configurations": [
|
||||
{
|
||||
"name": "splitpole",
|
||||
"name": "Hades",
|
||||
"includePath": [
|
||||
"${workspaceFolder}/**",
|
||||
"/home/splitpole/cern/root/**"
|
||||
],
|
||||
"defines": [],
|
||||
"compilerPath": "/usr/bin/gcc",
|
||||
"cStandard": "c17",
|
||||
"cppStandard": "gnu++17",
|
||||
"intelliSenseMode": "linux-gcc-x64"
|
||||
},
|
||||
{
|
||||
"name": "Ryan",
|
||||
"includePath": [
|
||||
"${workspaceFolder}/**",
|
||||
"/home/ryan/Downloads/root_build/**",
|
||||
"/home/ryan/Downloads/root_build/include"
|
||||
"/usr/include/x86_64-linux-gnu/qt6/**",
|
||||
"/usr/local/cern/root_v6.26.06/include/**"
|
||||
],
|
||||
"defines": [],
|
||||
"compilerPath": "/usr/bin/gcc",
|
||||
|
@ -29,7 +17,8 @@
|
|||
"name": "RyanUbuntu",
|
||||
"includePath": [
|
||||
"${workspaceFolder}/**",
|
||||
"/opt/root/**"
|
||||
"/usr/include/x86_64-linux-gnu/qt6/**",
|
||||
"/opt/root/include/**"
|
||||
],
|
||||
"defines": [],
|
||||
"compilerPath": "/usr/bin/gcc",
|
||||
|
@ -38,10 +27,11 @@
|
|||
"intelliSenseMode": "linux-gcc-x64"
|
||||
},
|
||||
{
|
||||
"name": "RyanHome",
|
||||
"name": "Anasen",
|
||||
"includePath": [
|
||||
"${workspaceFolder}/**",
|
||||
"/home/ryan/root_v6.30.06/**"
|
||||
"/usr/include/x86_64-linux-gnu/qt6/**",
|
||||
"/opt/root/include/**"
|
||||
],
|
||||
"defines": [],
|
||||
"compilerPath": "/usr/bin/gcc",
|
||||
|
@ -50,10 +40,26 @@
|
|||
"intelliSenseMode": "linux-gcc-x64"
|
||||
},
|
||||
{
|
||||
"name": "Dirac",
|
||||
"name": "Splitpole",
|
||||
"includePath": [
|
||||
"${workspaceFolder}/**",
|
||||
"/usr/opt/root/**"
|
||||
"/usr/include/x86_64-linux-gnu/qt6/**",
|
||||
"/home/splitpole/cern/root/include/**",
|
||||
"/usr/include/x86_64-linux-gnu/qt6/QtWidgets",
|
||||
"/usr/include/x86_64-linux-gnu/qt6/QtCore"
|
||||
],
|
||||
"defines": [],
|
||||
"compilerPath": "/usr/bin/gcc",
|
||||
"cStandard": "c17",
|
||||
"cppStandard": "gnu++17",
|
||||
"intelliSenseMode": "linux-gcc-x64"
|
||||
},
|
||||
{
|
||||
"name": "Penguin",
|
||||
"includePath": [
|
||||
"${workspaceFolder}/**",
|
||||
"/usr/include/x86_64-linux-gnu/qt6/**",
|
||||
"/usr/local/cern/root/include/**",
|
||||
],
|
||||
"defines": [],
|
||||
"compilerPath": "/usr/bin/gcc",
|
||||
|
|
24
.vscode/settings.json
vendored
24
.vscode/settings.json
vendored
|
@ -90,6 +90,26 @@
|
|||
"processRun.C": "cpp",
|
||||
"TrackRecon.C": "cpp",
|
||||
"processRuns.C": "cpp",
|
||||
"Analysis.C": "cpp"
|
||||
}
|
||||
"Analysis.C": "cpp",
|
||||
"datastructs.h": "c",
|
||||
"ANASENPlotEdit.C": "cpp",
|
||||
"GetMean_Q3_new.C": "cpp",
|
||||
"AlphaCal_new.C": "cpp",
|
||||
"f1.C": "cpp",
|
||||
"GeoCal_Maria_new.C": "cpp",
|
||||
"PCPulser_All_new.C": "cpp",
|
||||
"PosCal_2.C": "cpp",
|
||||
"AutoFit.C": "cpp",
|
||||
"Fitting.C": "cpp",
|
||||
"PCGainMatch.C": "cpp",
|
||||
"Analyzer1.C": "cpp",
|
||||
"FitHistogramsWithTSpectrum_Sequential_Improved.C": "cpp",
|
||||
"PlotAndFitCentroids.C": "cpp",
|
||||
"MatchAndPlotCentroids.C": "cpp",
|
||||
"GainMatch.C": "cpp",
|
||||
"GainMatchSX3.C": "cpp",
|
||||
"RelBack_Fix_new.C": "cpp",
|
||||
"SiRelativeGains_Step1_new.C": "cpp"
|
||||
},
|
||||
"github-enterprise.uri": "https://fsunuc.physics.fsu.edu"
|
||||
}
|
774
Analyzer.C
774
Analyzer.C
|
@ -1,84 +1,147 @@
|
|||
#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"
|
||||
TH2F *hsx3IndexVE;
|
||||
TH2F *hqqqIndexVE;
|
||||
TH2F *hpcIndexVE;
|
||||
|
||||
#include "TVector3.h"
|
||||
TH2F *hpcIndexVE_GM;
|
||||
TH2F *hsx3Coin;
|
||||
TH2F *hqqqCoin;
|
||||
TH2F *hpcCoin;
|
||||
TH2F *hAVCcoin;
|
||||
|
||||
TH2F * hsx3IndexVE;
|
||||
TH2F * hqqqIndexVE;
|
||||
TH2F * hpcIndexVE;
|
||||
TH2F *hqqqPolar;
|
||||
TH2F *hsx3VpcIndex;
|
||||
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 * 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;
|
||||
PW pwinstance;
|
||||
TVector3 hitPos;
|
||||
// TVector3 anodeIntersection;
|
||||
std::map<int, std::pair<double, double>> slopeInterceptMap;
|
||||
|
||||
bool HitNonZero;
|
||||
bool sx3ecut;
|
||||
bool qqqEcut;
|
||||
|
||||
TH1F * hZProj;
|
||||
TH1F *hZProj;
|
||||
TH1F *hPCZProj;
|
||||
|
||||
void Analyzer::Begin(TTree * /*tree*/){
|
||||
void Analyzer::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");
|
||||
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, 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);
|
||||
|
||||
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);
|
||||
|
||||
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");
|
||||
|
||||
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, 800, 0, 16000);
|
||||
hqqqVpcE->SetNdivisions(-612, "x");
|
||||
hqqqVpcE->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");
|
||||
hsx3VpcE = new TH2F("hsx3VpcEnergy", "sx3 vs pc; sx3 energy; pc energy", 400, 0, 5000, 800, 0, 16000);
|
||||
hsx3VpcE->SetNdivisions(-612, "x");
|
||||
hsx3VpcE->SetNdivisions(-12, "y");
|
||||
|
||||
hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600);
|
||||
|
||||
hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 10000, 400, 0 , 16000);
|
||||
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, 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);
|
||||
// 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){
|
||||
Bool_t Analyzer::Process(Long64_t entry)
|
||||
{
|
||||
|
||||
// if ( entry > 100 ) return kTRUE;
|
||||
|
||||
|
@ -110,46 +173,61 @@ Bool_t Analyzer::Process(Long64_t entry){
|
|||
|
||||
// sx3.Print();
|
||||
|
||||
//########################################################### Raw data
|
||||
// //======================= SX3
|
||||
|
||||
// ########################################################### Raw data
|
||||
// //======================= SX3
|
||||
sx3ecut = false;
|
||||
std::vector<std::pair<int, int>> ID; // first = id, 2nd = index
|
||||
for( int i = 0; i < sx3.multi; i ++){
|
||||
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]);
|
||||
|
||||
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]);
|
||||
if (sx3.e[i] > 100)
|
||||
{
|
||||
sx3ecut = true;
|
||||
}
|
||||
|
||||
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] );
|
||||
// }
|
||||
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++)
|
||||
{
|
||||
if (pc.index[j] < 24 && pc.e[j] > 100)
|
||||
{
|
||||
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;
|
||||
} );
|
||||
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;
|
||||
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) {
|
||||
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) {
|
||||
if (sx3ID.size() >= 3)
|
||||
{
|
||||
found = true;
|
||||
}
|
||||
}else{
|
||||
if( !found ){
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
if (!found)
|
||||
{
|
||||
sx3ID.clear();
|
||||
sx3ID.push_back(ID[i]);
|
||||
}
|
||||
|
@ -158,33 +236,42 @@ Bool_t Analyzer::Process(Long64_t entry){
|
|||
|
||||
// printf("---------- sx3ID Multi : %zu \n", sx3ID.size());
|
||||
|
||||
if( found ){
|
||||
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++ ){
|
||||
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) {
|
||||
if (sx3.ch[index] < 8)
|
||||
{
|
||||
if (sx3.ch[index] % 2 == 0)
|
||||
{
|
||||
sx3ChDn = sx3.ch[index];
|
||||
sx3EDn = sx3.e[index];
|
||||
}else{
|
||||
}
|
||||
else
|
||||
{
|
||||
sx3ChUp = sx3.ch[index];
|
||||
sx3EUp = sx3.e[index];
|
||||
}
|
||||
}else{
|
||||
}
|
||||
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] );
|
||||
}
|
||||
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);
|
||||
|
@ -192,211 +279,474 @@ Bool_t Analyzer::Process(Long64_t entry){
|
|||
HitNonZero = true;
|
||||
// hitPos.Print();
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
// //======================= QQQ
|
||||
for( int i = 0; i < qqq.multi; i ++){
|
||||
|
||||
qqqEcut = false;
|
||||
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] );
|
||||
// }
|
||||
// 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]);
|
||||
// printf("QQQ ID : %d, ch : %d, e : %d \n", qqq.id[i], qqq.ch[i], qqq.e[i]);
|
||||
if (qqq.e[i] > 100)
|
||||
{
|
||||
qqqEcut = true;
|
||||
}
|
||||
// }
|
||||
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] );
|
||||
|
||||
}
|
||||
// }
|
||||
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[k]);
|
||||
}
|
||||
}
|
||||
|
||||
// }
|
||||
|
||||
for (int j = i + 1; j < qqq.multi; j++)
|
||||
{
|
||||
// if( qqq.used[i] == true ) continue;
|
||||
|
||||
//if( qqq.id[i] == qqq.id[j] && (16 - qqq.ch[i]) * (16 - qqq.ch[j]) < 0 ){ // must be same detector and wedge and ring
|
||||
if( qqq.id[i] == qqq.id[j] ){ // must be same detector
|
||||
|
||||
|
||||
// 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]){
|
||||
int chRing = -1;
|
||||
if (qqq.ch[i] < qqq.ch[j])
|
||||
{
|
||||
chRing = qqq.ch[j] - 16;
|
||||
chWedge = qqq.ch[i];
|
||||
}else{
|
||||
}
|
||||
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);
|
||||
double theta = -TMath::Pi() / 2 + 2 * TMath::Pi() / 16 / 4. * (qqq.id[i] * 16 + chWedge + 0.5);
|
||||
double rho = 50. + 40. / 16. * (chRing + 0.5);
|
||||
// if(qqq.e[i]>50){
|
||||
hqqqPolar->Fill( theta, rho);
|
||||
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;
|
||||
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 ++){
|
||||
// Calculate the crossover points and put them into an array
|
||||
|
||||
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]));
|
||||
pwinstance.ConstructGeo();
|
||||
Coord Crossover[24][24][2];
|
||||
TVector3 a, c, diff;
|
||||
double a2, ac, c2, adiff, cdiff, denom, alpha, beta;
|
||||
int index = 0;
|
||||
for (int i = 0; i < pwinstance.An.size(); i++)
|
||||
{
|
||||
a = pwinstance.An[i].first - pwinstance.An[i].second;
|
||||
|
||||
hpcIndexVE->Fill( pc.index[i], pc.e[i] );
|
||||
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,
|
||||
|
||||
for( int j = i+1; j < pc.multi; j++){
|
||||
hpcCoin->Fill( pc.index[i], pc.index[j]);
|
||||
|
||||
c = pwinstance.Ca[j].first - pwinstance.Ca[j].second;
|
||||
diff = pwinstance.An[i].first - pwinstance.Ca[j].first;
|
||||
a2 = a.Dot(a);
|
||||
c2 = c.Dot(c);
|
||||
ac = a.Dot(c);
|
||||
adiff = a.Dot(diff);
|
||||
cdiff = c.Dot(diff);
|
||||
denom = a2 * c2 - ac * ac;
|
||||
alpha = (ac * cdiff - c2 * adiff) / denom;
|
||||
beta = (a2 * cdiff - ac * adiff) / denom;
|
||||
|
||||
Crossover[i][j][0].x = pwinstance.An[i].first.X() + alpha * a.X();
|
||||
Crossover[i][j][0].y = pwinstance.An[i].first.Y() + alpha * a.Y();
|
||||
Crossover[i][j][0].z = pwinstance.An[i].first.Z() + alpha * a.Z();
|
||||
if (Crossover[i][j][0].z < -190 || Crossover[i][j][0].z > 190)
|
||||
{
|
||||
Crossover[i][j][0].z = 9999999;
|
||||
}
|
||||
// placeholder variable Crossover[i][j][1].x has nothing to do with the geometry of the crossover and is being used to store the alpha value-
|
||||
//-so that it can be used to sort "good" hits later
|
||||
Crossover[i][j][1].x = alpha;
|
||||
Crossover[i][j][1].y = 0;
|
||||
// if(i==0){
|
||||
// printf("CID, Crossover z and alpha are : %d %f %f \n", j, Crossover[i][j][0].z, Crossover[i][j][1].x /*this is alpha*/);
|
||||
// }
|
||||
// }
|
||||
// }
|
||||
}
|
||||
}
|
||||
// printf("Anode and cathode indices, alpha, denom, andiff, cndiff : %d %d %f %f %f %f\n", i, j, alpha, denom, adiff, cdiff);
|
||||
|
||||
// anodeIntersection.Clear();
|
||||
for (int i = 0; i < pc.multi; i++)
|
||||
{
|
||||
|
||||
if (pc.e[i] > 100)
|
||||
{
|
||||
hpcIndexVE->Fill(pc.index[i], pc.e[i]); // non gain matched energy
|
||||
}
|
||||
|
||||
// Gain Matching of PC wires
|
||||
if (pc.index[i] >= 0 && pc.index[i] < 48)
|
||||
{
|
||||
// printf("index: %d, Old cathode energy: %d \n", pc.index[i],pc.e[i]);
|
||||
auto it = slopeInterceptMap.find(pc.index[i]);
|
||||
if (it != slopeInterceptMap.end())
|
||||
{
|
||||
double slope = it->second.first;
|
||||
double intercept = it->second.second;
|
||||
// printf("slope: %f, intercept:%f\n" ,slope, intercept);
|
||||
pc.e[i] = slope * pc.e[i] + intercept;
|
||||
// printf("index: %d, New cathode energy: %d \n",pc.index[i], pc.e[i]);
|
||||
}
|
||||
hpcIndexVE_GM->Fill(pc.index[i], pc.e[i]);
|
||||
// hPC_E[pc.index[i]]->Fill(pc.e[i]); // gain matched energy per channel
|
||||
}
|
||||
}
|
||||
// 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 ){
|
||||
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;
|
||||
|
||||
int aID = 0;
|
||||
int cID = 0;
|
||||
// Define the excluded SX3 and QQQ channels
|
||||
// std::unordered_set<int> excludeSX3 = {34, 35, 36, 37, 61, 62, 67, 73, 74, 75, 76, 77, 78, 79, 80, 93, 97, 100, 103, 108, 109, 110, 111, 112};
|
||||
// std::unordered_set<int> excludeQQQ = {0, 17, 109, 110, 111, 112, 113, 119, 127, 128};
|
||||
// inCuth=false;
|
||||
// inCutl=false;
|
||||
// inPCCut=false;
|
||||
for (int i = 0; i < pc.multi; i++)
|
||||
{
|
||||
if (pc.e[i] > 100 /*&& pc.multi < 7*/)
|
||||
{
|
||||
// creating a vector of pairs of anode and cathode hits
|
||||
if (pc.index[i] < 24)
|
||||
{
|
||||
anodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
|
||||
}
|
||||
else if (pc.index[i] >= 24)
|
||||
{
|
||||
cathodeHits.push_back(std::pair<int, double>(pc.index[i] - 24, pc.e[i]));
|
||||
}
|
||||
|
||||
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 j = i + 1; j < pc.multi; j++)
|
||||
{
|
||||
// if(PCCoinc_cut1->IsInside(pc.index[i], pc.index[j]) || PCCoinc_cut2->IsInside(pc.index[i], pc.index[j])){
|
||||
// // hpcCoin->Fill(pc.index[i], pc.index[j]);
|
||||
// inPCCut = true;
|
||||
// }
|
||||
hpcCoin->Fill(pc.index[i], pc.index[j]);
|
||||
}
|
||||
}
|
||||
}
|
||||
// sorting the anode and cathode hits in descending order of energy
|
||||
std::sort(anodeHits.begin(), anodeHits.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b)
|
||||
{ return a.second > b.second; });
|
||||
std::sort(cathodeHits.begin(), cathodeHits.end(), [](const std::pair<int, double> &a, const std::pair<int, double> &b)
|
||||
{ return a.second > b.second; });
|
||||
|
||||
// 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 && E[l].first!=20 && E[l].first!=12){
|
||||
if(!multi_an){
|
||||
aE = E[l].second;
|
||||
}
|
||||
multi_an=true;
|
||||
bool SiPCflag;
|
||||
|
||||
corrcatMax.clear();
|
||||
if (anodeHits.size() >= 1 && cathodeHits.size() > 1)
|
||||
{
|
||||
if (((TMath::TanH(hitPos.Y() / hitPos.X())) > (TMath::TanH(a.Y() / a.X()) - TMath::PiOver4())) || ((TMath::TanH(hitPos.Y() / hitPos.X())) < (TMath::TanH(a.Y() / a.X()) + TMath::PiOver4())))
|
||||
{
|
||||
|
||||
for (const auto &anode : anodeHits)
|
||||
{
|
||||
aID = anode.first;
|
||||
aE = anode.second;
|
||||
aESum += aE;
|
||||
if (aE > aEMax)
|
||||
{
|
||||
aEMax = aE;
|
||||
aIDMax = aID;
|
||||
}
|
||||
else {
|
||||
cE = E[l].second + cE;
|
||||
if (aE > aEnextMax && aE < aEMax)
|
||||
{
|
||||
aEnextMax = aE;
|
||||
aIDnextMax = aID;
|
||||
}
|
||||
// for(const auto &cat : cathodeHits){
|
||||
// hAVCcoin->Fill(aID, cat.first);
|
||||
// }
|
||||
}
|
||||
|
||||
// std::cout << " Anode iD : " << aIDMax << " Energy : " << aEMax << std::endl;
|
||||
|
||||
// printf("aID : %d, aE : %f, cE : %f\n", aID, aE, cE);
|
||||
for (const auto &cathode : cathodeHits)
|
||||
{
|
||||
cID = cathode.first;
|
||||
cE = cathode.second;
|
||||
// std::cout << "Cathode ID : " << cID << " Energy : " << cE << std::endl;
|
||||
|
||||
hAVCcoin->Fill(aIDMax, cID);
|
||||
|
||||
// This section of code is used to find the cathodes are correlated with the max and next max anodes, as well as to figure out if there are any common cathodes
|
||||
// the anodes are correlated with the cathodes +/-3 from the anode number in the reverse order
|
||||
|
||||
for (int j = -4; j < 3; j++)
|
||||
{
|
||||
if ((aIDMax + 24 + j) % 24 == 23 - cID)
|
||||
/* the 23-cID is used to accomodate for the fact that the order of the cathodes was reversed relative top the physical geometry */
|
||||
// if (Crossover[aIDMax][cID][0].z != 9999999)
|
||||
{
|
||||
corrcatMax.push_back(std::pair<int, double>(cID, cE));
|
||||
cESum += cE;
|
||||
// printf("Max Anode : %d Correlated Cathode : %d Anode Energy : %f z value : %f \n", aIDMax, cID, cESum, Crossover[aIDMax][cID][1].z /*prints alpha*/);
|
||||
// std::cout << " Cathode iD : " << cID << " Energy : " << cE << std::endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
// }
|
||||
// }
|
||||
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());
|
||||
TVector3 anodeIntersection;
|
||||
anodeIntersection.Clear();
|
||||
// Implementing a method for PC reconstruction using a single Anode event
|
||||
// if (anodeHits.size() == 1)
|
||||
{
|
||||
float x, y, z = 0;
|
||||
for (const auto &corr : corrcatMax)
|
||||
{
|
||||
if (cESum > 0)
|
||||
{
|
||||
x += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].x;
|
||||
y += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].y;
|
||||
z += (corr.second) / cESum * Crossover[aIDMax][corr.first][0].z;
|
||||
// printf("Max Anode : %d Correlated Cathode : %d cathode Energy : %f cESum Energy : %f z value : %f \n", aIDMax, corr.first, corr.second, cESum, Crossover[aIDMax][corr.first][1].z /*prints alpha*/);
|
||||
}
|
||||
else
|
||||
{
|
||||
printf("Warning: No valid cathode hits to correlate with anode %d! \n", aIDMax);
|
||||
}
|
||||
// printf("EventID : %llu, Max Anode : %d Cathode: %d PC X and Y : (%f, %f) \n", entry, aIDMax, cID, Crossover[aIDMax][cID][0].x, Crossover[aIDMax][cID][0].y);
|
||||
// for (int i = 0; i < sx3.multi; i++)
|
||||
// {
|
||||
// printf("EventID : %llu, HitPos X, Y, Z: %f %f %f SX3ID : %d %d \n", entry, hitPos.X(), hitPos.Y(), hitPos.Z(), sx3.id[i], sx3.ch[i]);
|
||||
// }
|
||||
|
||||
// for (int i = 0; i < qqq.multi; i++)
|
||||
// {
|
||||
// printf("Max Anode : %d Cathode: %d PC X and Y : %f %f \n", aIDMax, cID, Crossover[aIDMax][cID][0].x, Crossover[aIDMax][cID][0].y);
|
||||
// printf("HitPos X, Y, Z, QQQID : %f %f %f %d \n", hitPos.X(), hitPos.Y(), hitPos.Z(), qqq.id[i]);
|
||||
// }
|
||||
}
|
||||
|
||||
|
||||
anodeIntersection = TVector3(x, y, z);
|
||||
// std::cout << "Anode Intersection " << anodeIntersection.Z() << " " << x << " " << y << " " << z << std::endl;
|
||||
}
|
||||
|
||||
if(anodeIntersection.Z() != 0){
|
||||
hPCZProj->Fill(anodeIntersection.Z());
|
||||
}
|
||||
// Filling the PC Z projection histogram
|
||||
// std::cout << anodeIntersection.Z() << std::endl;
|
||||
// hPCZProj->Fill(anodeIntersection.Z());
|
||||
|
||||
// }
|
||||
|
||||
//########################################################### Track constrcution
|
||||
// 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;
|
||||
// }
|
||||
|
||||
//############################## DO THE KINEMATICS
|
||||
// Fill histograms based on the cut conditions
|
||||
// if (inCuth && inPCCut) {
|
||||
// hanVScatsum_hcut->Fill(aESum, cESum);
|
||||
// }
|
||||
// if (inCutl && inPCCut) {
|
||||
// hanVScatsum_lcut->Fill(aESum, cESum);
|
||||
// }
|
||||
// for(auto anode : anodeHits){
|
||||
|
||||
// float aE = anode.second;
|
||||
// aESum += aE;
|
||||
// if(inPCCut){
|
||||
hanVScatsum->Fill(aEMax, cESum);
|
||||
// }
|
||||
|
||||
// if (sx3ecut)
|
||||
// {
|
||||
hCat4An->Fill(corrcatMax.size());
|
||||
hNosvAe->Fill(corrcatMax.size(), aEMax);
|
||||
hAnodehits->Fill(anodeHits.size());
|
||||
// }
|
||||
|
||||
// }
|
||||
if (anodeHits.size() < 1)
|
||||
{
|
||||
hCat0An->Fill(cathodeHits.size());
|
||||
}
|
||||
|
||||
if (HitNonZero && anodeIntersection.Z() != 0)
|
||||
{
|
||||
pw_contr.CalTrack2(hitPos, anodeIntersection);
|
||||
hZProj->Fill(pw_contr.GetZ0());
|
||||
}
|
||||
|
||||
// ########################################################### Track constrcution
|
||||
|
||||
// ############################## DO THE KINEMATICS
|
||||
|
||||
return kTRUE;
|
||||
}
|
||||
|
||||
void Analyzer::Terminate(){
|
||||
void Analyzer::Terminate()
|
||||
{
|
||||
|
||||
gStyle->SetOptStat("neiou");
|
||||
TCanvas * canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000);
|
||||
canvas->Divide(3,3);
|
||||
// gStyle->SetOptStat("neiou");
|
||||
// TCanvas *canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000);
|
||||
// canvas->Divide(3, 3);
|
||||
|
||||
//hsx3VpcIndex->Draw("colz");
|
||||
// // hsx3VpcIndex->Draw("colz");
|
||||
|
||||
//=============================================== pad-1
|
||||
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
|
||||
// //=============================================== pad-1
|
||||
// padID++;
|
||||
// canvas->cd(padID);
|
||||
// canvas->cd(padID)->SetGrid(1);
|
||||
|
||||
hsx3IndexVE->Draw("colz");
|
||||
// hsx3IndexVE->Draw("colz");
|
||||
|
||||
//=============================================== pad-2
|
||||
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
|
||||
// //=============================================== pad-2
|
||||
// padID++;
|
||||
// canvas->cd(padID);
|
||||
// canvas->cd(padID)->SetGrid(1);
|
||||
|
||||
hqqqIndexVE->Draw("colz");
|
||||
// hqqqIndexVE->Draw("colz");
|
||||
|
||||
//=============================================== pad-3
|
||||
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
|
||||
// //=============================================== pad-3
|
||||
// padID++;
|
||||
// canvas->cd(padID);
|
||||
// canvas->cd(padID)->SetGrid(1);
|
||||
|
||||
hpcIndexVE->Draw("colz");
|
||||
// hpcIndexVE->Draw("colz");
|
||||
|
||||
//=============================================== pad-4
|
||||
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
|
||||
// //=============================================== pad-4
|
||||
// padID++;
|
||||
// canvas->cd(padID);
|
||||
// canvas->cd(padID)->SetGrid(1);
|
||||
|
||||
hsx3Coin->Draw("colz");
|
||||
// hsx3Coin->Draw("colz");
|
||||
|
||||
//=============================================== pad-5
|
||||
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
|
||||
// //=============================================== pad-5
|
||||
// padID++;
|
||||
// canvas->cd(padID);
|
||||
// canvas->cd(padID)->SetGrid(1);
|
||||
|
||||
canvas->cd(padID)->SetLogz(true);
|
||||
// canvas->cd(padID)->SetLogz(true);
|
||||
|
||||
hqqqCoin->Draw("colz");
|
||||
// hqqqCoin->Draw("colz");
|
||||
|
||||
//=============================================== pad-6
|
||||
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
|
||||
// //=============================================== pad-6
|
||||
// padID++;
|
||||
// canvas->cd(padID);
|
||||
// canvas->cd(padID)->SetGrid(1);
|
||||
|
||||
hpcCoin->Draw("colz");
|
||||
// hpcCoin->Draw("colz");
|
||||
|
||||
//=============================================== pad-7
|
||||
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
|
||||
// //=============================================== pad-7
|
||||
// padID++;
|
||||
// canvas->cd(padID);
|
||||
// canvas->cd(padID)->SetGrid(1);
|
||||
|
||||
// hsx3VpcIndex ->Draw("colz");
|
||||
hsx3VpcE->Draw("colz") ;
|
||||
// // hsx3VpcIndex ->Draw("colz");
|
||||
// hsx3VpcE->Draw("colz");
|
||||
|
||||
//=============================================== pad-8
|
||||
padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
|
||||
// //=============================================== pad-8
|
||||
// padID++;
|
||||
// canvas->cd(padID);
|
||||
// canvas->cd(padID)->SetGrid(1);
|
||||
|
||||
// hqqqVpcIndex ->Draw("colz");
|
||||
// // hqqqVpcIndex ->Draw("colz");
|
||||
|
||||
hqqqVpcE ->Draw("colz");
|
||||
//=============================================== pad-9
|
||||
padID ++;
|
||||
// hqqqVpcE->Draw("colz");
|
||||
// //=============================================== pad-9
|
||||
// padID++;
|
||||
|
||||
// canvas->cd(padID)->DrawFrame(-50, -50, 50, 50);
|
||||
// hqqqPolar->Draw("same colz pol");
|
||||
// // 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");
|
||||
// canvas->cd(padID);
|
||||
// canvas->cd(padID)->SetGrid(1);
|
||||
// // hZProj->Draw();
|
||||
// hanVScatsum->Draw("colz");
|
||||
|
||||
// // TFile *outRoot = new TFile("Histograms.root", "RECREATE");
|
||||
|
||||
// // if (!outRoot->IsOpen())
|
||||
// // {
|
||||
// // std::cerr << "Error opening file for writing!" << std::endl;
|
||||
// // return;
|
||||
// // }
|
||||
|
||||
// // // Loop through histograms and write them to the ROOT file
|
||||
// // for (int i = 0; i < 48; i++)
|
||||
// // {
|
||||
// // if (hPC_E[i] != nullptr)
|
||||
// // {
|
||||
// // hPC_E[i]->Write(); // Write histogram to file
|
||||
// // }
|
||||
// // }
|
||||
|
||||
// // outRoot->Close();
|
||||
}
|
||||
|
|
15
Analyzer.h
15
Analyzer.h
|
@ -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
402
Analyzer1.C
Normal 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");
|
||||
|
||||
}
|
114
Analyzer1.h
Normal file
114
Analyzer1.h
Normal file
|
@ -0,0 +1,114 @@
|
|||
#ifndef Analyzer1_h
|
||||
#define Analyzer1_h
|
||||
|
||||
#include <TROOT.h>
|
||||
#include <TChain.h>
|
||||
#include <TFile.h>
|
||||
#include <TSelector.h>
|
||||
|
||||
#include "Armory/ClassDet.h"
|
||||
|
||||
class Analyzer1 : 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; //!
|
||||
|
||||
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);
|
||||
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(Analyzer1,0);
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
#ifdef Analyzer1_cxx
|
||||
void Analyzer1::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 Analyzer1::Notify(){
|
||||
|
||||
return kTRUE;
|
||||
}
|
||||
|
||||
void Analyzer1::SlaveBegin(TTree * /*tree*/){
|
||||
|
||||
TString option = GetOption();
|
||||
|
||||
}
|
||||
|
||||
void Analyzer1::SlaveTerminate(){
|
||||
|
||||
}
|
||||
|
||||
|
||||
#endif // #ifdef Analyzer_cxx
|
|
@ -103,8 +103,8 @@ void ANASEN_model(int anodeID1 = -1, int anodeID2 = -1, int cathodeID1 = -1, int
|
|||
new TGeoRotation("rot1", 360/nSX3 * (i + 0.5), 0., 0.)));
|
||||
}
|
||||
|
||||
const int qqqR1 = 10;
|
||||
const int qqqR2 = 50;
|
||||
const int qqqR1 = 50;
|
||||
const int qqqR2 = 100;
|
||||
TGeoVolume *qqq = geom->MakeTubs("qqq", Al, qqqR1, qqqR2, 0.5, 5, 85);
|
||||
qqq->SetLineColor(7);
|
||||
for( int i = 0; i < 4; i++){
|
||||
|
|
327
Armory/ClassPW.h
327
Armory/ClassPW.h
|
@ -4,15 +4,18 @@
|
|||
#include <cstdio>
|
||||
#include <TMath.h>
|
||||
#include <TVector3.h>
|
||||
#include <TRandom.h>
|
||||
|
||||
struct PWHitInfo{
|
||||
std::pair<short, short> nearestWire; // anode, cathode
|
||||
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<short, short> nextNearestWire; // anode, cathode
|
||||
std::pair<double, double> nextNearestDist; // anode, cathode
|
||||
|
||||
void Clear(){
|
||||
void Clear()
|
||||
{
|
||||
nearestWire.first = -1;
|
||||
nearestWire.second = -1;
|
||||
nearestDist.first = 999999999;
|
||||
|
@ -24,63 +27,77 @@ struct PWHitInfo{
|
|||
}
|
||||
};
|
||||
|
||||
//!########################################################
|
||||
class PW{ // proportional wire
|
||||
public:
|
||||
PW(){ ClearHitInfo();};
|
||||
~PW(){};
|
||||
struct Coord
|
||||
{
|
||||
float x, y, z;
|
||||
Coord() : x(0), y(0), z(0) {}
|
||||
Coord(const TVector3 &vec)
|
||||
{
|
||||
x = vec.X(); // TVector3's X() returns the x-coordinate
|
||||
y = vec.Y(); // TVector3's Y() returns the y-coordinate
|
||||
z = vec.Z(); // TVector3's Z() returns the z-coordinate
|
||||
}
|
||||
};
|
||||
|
||||
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;}
|
||||
//! ########################################################
|
||||
class PW
|
||||
{ // proportional wire
|
||||
public:
|
||||
PW() { ClearHitInfo(); };
|
||||
~PW() {};
|
||||
|
||||
TVector3 GetTrackPos() const {return trackPos;}
|
||||
TVector3 GetTrackVec() const {return trackVec;}
|
||||
double GetTrackTheta() const {return trackVec.Theta();}
|
||||
double GetTrackPhi() const {return trackVec.Phi();}
|
||||
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; }
|
||||
|
||||
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
|
||||
|
||||
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;}
|
||||
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 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();}
|
||||
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 CalTrack2(TVector3 sx3Pos, TVector3 anodeInt, bool verbose = false);
|
||||
|
||||
double CircularMean(std::vector<std::pair<int, double>> wireList);
|
||||
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);
|
||||
|
||||
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);
|
||||
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;
|
||||
|
@ -88,7 +105,7 @@ private:
|
|||
|
||||
const int nWire = 24;
|
||||
const int wireShift = 3;
|
||||
const float zLen = 380; //mm
|
||||
const float zLen = 380; // mm
|
||||
const float radiusA = 37;
|
||||
const float radiusC = 43;
|
||||
|
||||
|
@ -96,23 +113,25 @@ private:
|
|||
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
|
||||
// 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){
|
||||
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));
|
||||
}
|
||||
|
||||
return TMath::Abs(nd.Dot(a1 - b2));
|
||||
}
|
||||
};
|
||||
|
||||
inline void PW::ClearHitInfo(){
|
||||
inline void PW::ClearHitInfo()
|
||||
{
|
||||
hitInfo.Clear();
|
||||
}
|
||||
|
||||
inline void PW::ConstructGeo(){
|
||||
inline void PW::ConstructGeo()
|
||||
{
|
||||
|
||||
An.clear();
|
||||
Ca.clear();
|
||||
|
@ -120,108 +139,137 @@ inline void PW::ConstructGeo(){
|
|||
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.
|
||||
// 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++ ){
|
||||
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);
|
||||
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);
|
||||
// Cathod rotate left-hand with the 3 wire offset accounted for (+1 from the calculated offset from the PC coincidence spectrum)
|
||||
q1.first.SetXYZ(radiusC * TMath::Cos(TMath::TwoPi() / nWire * (i + wireShift + 1) + TMath::PiOver2()),
|
||||
radiusC * TMath::Sin(TMath::TwoPi() / nWire * (i + wireShift + 1) + TMath::PiOver2()),
|
||||
zLen / 2);
|
||||
q1.second.SetXYZ(radiusC * TMath::Cos(TMath::TwoPi() / nWire * (i + 1) + TMath::PiOver2()),
|
||||
radiusC * TMath::Sin(TMath::TwoPi() / nWire * (i + 1) + TMath::PiOver2()),
|
||||
-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
|
||||
// NOT NECESSARY ANY MORE, HAS BEEN IMCORPORATED INTO THE WIREOFFSET IN THE BEGINNING
|
||||
// std::rotate(Ca.rbegin(), Ca.rbegin() + 4, Ca.rend());
|
||||
|
||||
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) );
|
||||
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 ){
|
||||
inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose)
|
||||
{
|
||||
|
||||
hitInfo.Clear();
|
||||
double phi = direction.Phi();
|
||||
|
||||
for( int i = 0; i < nWire; i++){
|
||||
for (int i = 0; i < nWire; i++)
|
||||
{
|
||||
|
||||
double disA = 99999999;
|
||||
double phiS = An[i].first.Phi() - TMath::PiOver4();
|
||||
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 (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 ){
|
||||
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();
|
||||
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 (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 ){
|
||||
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);
|
||||
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;
|
||||
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){
|
||||
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{
|
||||
}
|
||||
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;
|
||||
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){
|
||||
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{
|
||||
}
|
||||
else
|
||||
{
|
||||
hitInfo.nextNearestWire.second = ccc2;
|
||||
hitInfo.nextNearestDist.second = haha2;
|
||||
}
|
||||
|
||||
if( verbose ) Print();
|
||||
if (verbose)
|
||||
Print();
|
||||
}
|
||||
|
||||
inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose){
|
||||
inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbose)
|
||||
{
|
||||
|
||||
trackPos = sx3Pos;
|
||||
|
||||
|
@ -231,72 +279,37 @@ inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbo
|
|||
// 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());
|
||||
|
||||
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;
|
||||
inline void PW::CalTrack2(TVector3 siPos, TVector3 anodeInt, bool verbose)
|
||||
{
|
||||
|
||||
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());
|
||||
float mx, my;
|
||||
double z;
|
||||
mx = siPos.X() / (siPos.X() - anodeInt.X());
|
||||
my = siPos.Y() / (siPos.Y() - anodeInt.Y());
|
||||
z=siPos.Z() + mx * (anodeInt.Z() - siPos.Z());
|
||||
// if (mx == my)
|
||||
{
|
||||
trackVec=TVector3(0,0,z);
|
||||
}
|
||||
|
||||
if (verbose)
|
||||
printf("X slope = %f and Y slope = %f \n", mx, my);
|
||||
}
|
||||
|
||||
inline double PW::GetZ0(){
|
||||
inline double PW::GetZ0()
|
||||
{
|
||||
|
||||
double x = trackPos.X();
|
||||
double y = trackPos.Y();
|
||||
double rho = TMath::Sqrt(x*x + y*y);
|
||||
double rho = TMath::Sqrt(x * x + y * y);
|
||||
double theta = trackVec.Theta();
|
||||
|
||||
return trackPos.Z() - rho / TMath::Tan(theta);
|
||||
|
||||
return trackVec.Z();
|
||||
}
|
||||
|
||||
inline double PW::CircularMean(std::vector<std::pair<int, double>> wireList){
|
||||
|
||||
//use unit vector, wireID start from Zero
|
||||
double xCom = 0, yCom = 0;
|
||||
for( size_t i = 0; i < wireList.size() ; i++){
|
||||
xCom += TMath::Cos(TMath::TwoPi() * wireList[i].first / nWire) * wireList[i].second;
|
||||
yCom += TMath::Sin(TMath::TwoPi() * wireList[i].first / nWire) * wireList[i].second;
|
||||
}
|
||||
|
||||
//calculate the angle of the summed unit vectors
|
||||
double angle = TMath::ATan2(yCom, xCom);
|
||||
if( angle < 0 ) angle += TMath::TwoPi(); // convert the angle from 0 to 2 pi
|
||||
|
||||
return angle/ TMath::TwoPi() * nWire;
|
||||
|
||||
}
|
||||
|
||||
#endif
|
||||
#endif
|
124
FitHistogramsWithTSpectrum_Sequential_Improved.C
Normal file
124
FitHistogramsWithTSpectrum_Sequential_Improved.C
Normal 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;
|
||||
}
|
||||
|
264
GainMatch.C
Normal file
264
GainMatch.C
Normal file
|
@ -0,0 +1,264 @@
|
|||
#define GainMatch_cxx
|
||||
|
||||
#include "GainMatch.h"
|
||||
#include <TH2.h>
|
||||
#include <TF1.h>
|
||||
#include <TStyle.h>
|
||||
#include <TCanvas.h>
|
||||
#include <TMath.h>
|
||||
#include <TCutG.h>
|
||||
#include <fstream>
|
||||
#include <utility>
|
||||
#include <algorithm>
|
||||
|
||||
#include "Armory/ClassSX3.h"
|
||||
|
||||
#include "TVector3.h"
|
||||
|
||||
TH2F *hSX3FvsB;
|
||||
TH2F *hQQQFVB;
|
||||
|
||||
int padID = 0;
|
||||
|
||||
SX3 sx3_contr;
|
||||
TCutG *cut;
|
||||
std::map<std::tuple<int, int, int>, std::vector<std::pair<double, double>>> dataPoints;
|
||||
|
||||
void GainMatch::Begin(TTree * /*tree*/)
|
||||
{
|
||||
TString option = GetOption();
|
||||
|
||||
hSX3FvsB = new TH2F("hSX3FvsB", "SX3 Front vs Back; Front E; Back E", 400, 0, 16000, 400, 0, 16000);
|
||||
hQQQFVB = new TH2F("hQQQFVB", "number of good QQQ vs QQQ id", 400, 0, 16000, 400, 0, 16000);
|
||||
|
||||
sx3_contr.ConstructGeo();
|
||||
|
||||
// Load the TCutG object
|
||||
TFile *cutFile = TFile::Open("qqqcorr.root");
|
||||
if (!cutFile || cutFile->IsZombie())
|
||||
{
|
||||
std::cerr << "Error: Could not open qqqcorr.root" << std::endl;
|
||||
return;
|
||||
}
|
||||
cut = dynamic_cast<TCutG *>(cutFile->Get("qqqcorr"));
|
||||
if (!cut)
|
||||
{
|
||||
std::cerr << "Error: Could not find TCutG named 'qqqcorr' in qqqcorr.root" << std::endl;
|
||||
return;
|
||||
}
|
||||
cut->SetName("qqqcorr"); // Ensure the cut has the correct name
|
||||
}
|
||||
|
||||
Bool_t GainMatch::Process(Long64_t 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();
|
||||
|
||||
std::vector<std::pair<int, int>> ID;
|
||||
for (int i = 0; i < sx3.multi; i++)
|
||||
{
|
||||
ID.push_back(std::pair<int, int>(sx3.id[i], i));
|
||||
}
|
||||
|
||||
if (ID.size() > 0)
|
||||
{
|
||||
std::sort(ID.begin(), ID.end(), [](const std::pair<int, int> &a, const std::pair<int, int> &b)
|
||||
{ return a.first < b.first; });
|
||||
|
||||
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];
|
||||
}
|
||||
}
|
||||
hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
}
|
||||
}
|
||||
|
||||
for (int i = 0; i < qqq.multi; i++)
|
||||
{
|
||||
for (int j = i + 1; j < qqq.multi; j++)
|
||||
{
|
||||
if (qqq.id[i] == qqq.id[j])
|
||||
{
|
||||
int chWedge = -1;
|
||||
int chRing = -1;
|
||||
float eWedge = 0.0;
|
||||
float eRing = 0.0;
|
||||
if (qqq.ch[i] < 16 && qqq.ch[j] >= 16)
|
||||
{
|
||||
chWedge = qqq.ch[i];
|
||||
eWedge = qqq.e[i];
|
||||
chRing = qqq.ch[j] - 16;
|
||||
eRing = qqq.e[j];
|
||||
}
|
||||
else if (qqq.ch[j] < 16 && qqq.ch[i] >= 16)
|
||||
{
|
||||
chWedge = qqq.ch[j];
|
||||
eWedge = qqq.e[j];
|
||||
chRing = qqq.ch[i] - 16;
|
||||
eRing = qqq.e[i];
|
||||
}
|
||||
else
|
||||
continue;
|
||||
|
||||
hQQQFVB->Fill(eWedge, eRing);
|
||||
|
||||
TString histName = Form("hQQQFVB_id%d_r%d_w%d", qqq.id[i], chRing, chWedge);
|
||||
TH2F *hist2d = (TH2F *)gDirectory->Get(histName);
|
||||
if (!hist2d)
|
||||
{
|
||||
hist2d = new TH2F(histName, Form("QQQ Det%d R%d W%d;Wedge E;Ring E", qqq.id[i], chRing, chWedge), 400, 0, 16000, 400, 0, 16000);
|
||||
}
|
||||
|
||||
hist2d->Fill(eWedge, eRing);
|
||||
if (cut && cut->IsInside(eWedge, eRing))
|
||||
{
|
||||
// Accumulate data for gain matching
|
||||
dataPoints[{qqq.id[i], chRing, chWedge}].emplace_back(eWedge, eRing);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return kTRUE;
|
||||
}
|
||||
|
||||
void GainMatch::Terminate()
|
||||
{
|
||||
const int MAX_DET = 4;
|
||||
const int MAX_RING = 16;
|
||||
const int MAX_WEDGE = 16;
|
||||
|
||||
double gainArray[MAX_DET][MAX_RING][MAX_WEDGE] = {{{0}}};
|
||||
bool gainValid[MAX_DET][MAX_RING][MAX_WEDGE] = {{{false}}};
|
||||
|
||||
std::ofstream outFile("qqq_gainmatch.txt");
|
||||
if (!outFile.is_open())
|
||||
{
|
||||
std::cerr << "Error opening output file!" << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
for (const auto &kv : dataPoints)
|
||||
{
|
||||
auto [id, ring, wedge] = kv.first;
|
||||
const auto &pts = kv.second;
|
||||
if (pts.size() < 5)
|
||||
continue;
|
||||
|
||||
std::vector<double> wE, rE;
|
||||
for (const auto &pr : pts)
|
||||
{
|
||||
wE.push_back(pr.first);
|
||||
rE.push_back(pr.second);
|
||||
}
|
||||
|
||||
TGraph g(wE.size(), wE.data(), rE.data());
|
||||
TF1 f("f", "[0]*x", 0, 16000);
|
||||
g.Fit(&f, "QNR");
|
||||
gainArray[id][ring][wedge] = f.GetParameter(0);
|
||||
gainValid[id][ring][wedge] = true;
|
||||
}
|
||||
|
||||
for (int id = 0; id < MAX_DET; ++id)
|
||||
{
|
||||
for (int ring = 0; ring < MAX_RING; ++ring)
|
||||
{
|
||||
for (int wedge = 0; wedge < MAX_WEDGE; ++wedge)
|
||||
{
|
||||
if (gainValid[id][ring][wedge])
|
||||
{
|
||||
outFile << id << " " << wedge << " " << ring << " " << gainArray[id][ring][wedge] << std::endl;
|
||||
printf("Gain match Det%d Ring%d Wedge%d → %.4f \n", id, ring, wedge, gainArray[id][ring][wedge]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
outFile.close();
|
||||
std::cout << "Gain matching complete." << std::endl;
|
||||
|
||||
// === Plot all gain-matched QQQ points together with a 2D histogram ===
|
||||
TH2F *hAll = new TH2F("hAll", "All QQQ Gain-Matched;Corrected Wedge E;Ring E",
|
||||
400, 0, 16000, 400, 0, 16000);
|
||||
|
||||
// Fill the combined TH2F with corrected data
|
||||
for (auto &kv : dataPoints)
|
||||
{
|
||||
int id, ring, wedge;
|
||||
std::tie(id, ring, wedge) = kv.first;
|
||||
if (!gainValid[id][ring][wedge])
|
||||
continue;
|
||||
auto &pts = kv.second;
|
||||
for (auto &pr : pts)
|
||||
{
|
||||
double corrWedge = pr.first * gainArray[id][ring][wedge];
|
||||
double ringE = pr.second;
|
||||
hAll->Fill(corrWedge, ringE);
|
||||
}
|
||||
}
|
||||
}
|
114
GainMatch.h
Normal file
114
GainMatch.h
Normal file
|
@ -0,0 +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
|
270
GainMatchSX3.C
Normal file
270
GainMatchSX3.C
Normal file
|
@ -0,0 +1,270 @@
|
|||
#define GainMatchSX3_cxx
|
||||
|
||||
#include "GainMatchSX3.h"
|
||||
#include <TH2.h>
|
||||
#include <TF1.h>
|
||||
#include <TStyle.h>
|
||||
#include <TCanvas.h>
|
||||
#include <TMath.h>
|
||||
#include <TCutG.h>
|
||||
#include <fstream>
|
||||
#include <utility>
|
||||
#include <algorithm>
|
||||
|
||||
#include "Armory/ClassSX3.h"
|
||||
|
||||
#include "TVector3.h"
|
||||
|
||||
TH2F *hSX3FvsB;
|
||||
TH2F *hsx3IndexVE;
|
||||
int padID = 0;
|
||||
|
||||
SX3 sx3_contr;
|
||||
TCutG *cut;
|
||||
std::map<std::tuple<int, int, int>, std::vector<std::tuple<double, double, double>>> dataPoints;
|
||||
|
||||
void GainMatchSX3::Begin(TTree * /*tree*/)
|
||||
{
|
||||
TString option = GetOption();
|
||||
|
||||
hSX3FvsB = new TH2F("hSX3FvsB", "SX3 Front vs Back; Front E; Back E", 400, 0, 16000, 400, 0, 16000);
|
||||
hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000);
|
||||
|
||||
sx3_contr.ConstructGeo();
|
||||
|
||||
// Load the TCutG object
|
||||
TFile *cutFile = TFile::Open("sx3cut.root");
|
||||
if (!cutFile || cutFile->IsZombie())
|
||||
{
|
||||
std::cerr << "Error: Could not open sx3cut.root" << std::endl;
|
||||
return;
|
||||
}
|
||||
cut = dynamic_cast<TCutG *>(cutFile->Get("sx3cut"));
|
||||
if (!cut)
|
||||
{
|
||||
std::cerr << "Error: Could not find TCutG named 'sx3cut' in sx3cut.root" << std::endl;
|
||||
return;
|
||||
}
|
||||
cut->SetName("sx3cut"); // Ensure the cut has the correct name
|
||||
}
|
||||
|
||||
Bool_t GainMatchSX3::Process(Long64_t 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();
|
||||
|
||||
std::vector<std::pair<int, int>> ID;
|
||||
for (int i = 0; i < sx3.multi; i++)
|
||||
{
|
||||
ID.push_back(std::pair<int, int>(sx3.id[i], i));
|
||||
}
|
||||
|
||||
if (ID.size() > 0)
|
||||
{
|
||||
std::sort(ID.begin(), ID.end(), [](const std::pair<int, int> &a, const std::pair<int, int> &b)
|
||||
{ return a.first < b.first; });
|
||||
|
||||
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] / 2;
|
||||
sx3EDn = sx3.e[index];
|
||||
}
|
||||
else
|
||||
{
|
||||
sx3ChUp = sx3.ch[index] / 2;
|
||||
sx3EUp = sx3.e[index];
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
sx3ChBk = sx3.ch[index] - 8;
|
||||
// if (sx3ChBk == 2)
|
||||
// printf("Found back channel Det %d Back %d \n", sx3.id[index], sx3ChBk);
|
||||
sx3EBk = sx3.e[index];
|
||||
}
|
||||
}
|
||||
hSX3FvsB->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
|
||||
// Fill the histogram for the front vs back
|
||||
std::array<int, 24> detectorIDs;
|
||||
|
||||
for (int i = 0; i < sx3.multi; i++)
|
||||
{
|
||||
if (sx3.id[i]==3)
|
||||
{
|
||||
TString histName = Form("hSX3FVB_id%d_U%d_D%d_B%d", sx3.id[i], sx3ChUp, sx3ChDn, sx3ChBk);
|
||||
TH2F *hist2d = (TH2F *)gDirectory->Get(histName);
|
||||
if (!hist2d)
|
||||
{
|
||||
hist2d = new TH2F(histName, Form("hSX3FVB_id%d_U%d_D%d_B%d", sx3.id[i], sx3ChUp, sx3ChDn, sx3ChBk), 400, 0, 16000, 400, 0, 16000);
|
||||
}
|
||||
|
||||
// if (sx3ChBk == 2)
|
||||
// printf("Found back channel Det %d Back %d \n", sx3.id[i], sx3ChBk);
|
||||
hsx3IndexVE->Fill(sx3.index[i], sx3.e[i]);
|
||||
|
||||
hist2d->Fill(sx3EUp + sx3EDn, sx3EBk);
|
||||
if (cut && cut->IsInside(sx3EUp + sx3EDn, sx3EBk))
|
||||
{
|
||||
// Accumulate data for gain matching
|
||||
dataPoints[{sx3.id[i], sx3ChUp, sx3ChBk}].emplace_back(sx3EBk, sx3EUp, sx3EDn);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return kTRUE;
|
||||
}
|
||||
void GainMatchSX3::Terminate()
|
||||
{
|
||||
const int MAX_DET = 24;
|
||||
const int MAX_UP = 4;
|
||||
const int MAX_DOWN = 4;
|
||||
const int MAX_BK = 4;
|
||||
|
||||
double gainArray[MAX_DET][MAX_UP][MAX_BK] = {{{0}}};
|
||||
bool gainValid[MAX_DET][MAX_UP][MAX_BK] = {{{false}}};
|
||||
|
||||
std::ofstream outFile("sx3_GainMatch.txt");
|
||||
if (!outFile.is_open())
|
||||
{
|
||||
std::cerr << "Error opening output file!" << std::endl;
|
||||
return;
|
||||
}
|
||||
|
||||
// === Updated dataPoints type ===
|
||||
// std::map<std::tuple<int, int, int>, std::vector<std::tuple<double, double, double>>> dataPoints;
|
||||
|
||||
// Gain fit using up+dn vs bk
|
||||
for (const auto &kv : dataPoints)
|
||||
{
|
||||
auto [id, ud, bk] = kv.first;
|
||||
const auto &pts = kv.second;
|
||||
if (pts.size() < 5)
|
||||
continue;
|
||||
|
||||
std::vector<double> bkE, udE;
|
||||
for (const auto &pr : pts)
|
||||
{
|
||||
double eUp, eDn, eBk;
|
||||
std::tie(eBk, eUp, eDn) = pr;
|
||||
bkE.push_back(eBk);
|
||||
udE.push_back(eUp + eDn);
|
||||
}
|
||||
|
||||
TGraph g(bkE.size(), bkE.data(), udE.data());
|
||||
TF1 f("f", "[0]*x", 0, 16000);
|
||||
g.Fit(&f, "QNR");
|
||||
gainArray[id][ud][bk] = f.GetParameter(0);
|
||||
gainValid[id][ud][bk] = true;
|
||||
}
|
||||
|
||||
// Output results
|
||||
for (int id = 0; id < MAX_DET; ++id)
|
||||
{
|
||||
for (int bk = 0; bk < MAX_BK; ++bk)
|
||||
{
|
||||
for (int ud = 0; ud < MAX_UP; ++ud)
|
||||
{
|
||||
if (gainValid[id][ud][bk])
|
||||
{
|
||||
outFile << id << " " << bk << " " << ud << " " << gainArray[id][ud][bk] << std::endl;
|
||||
printf("Gain match Det%d Up+Dn%d Back%d → %.4f \n", id, ud, bk, gainArray[id][ud][bk]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
outFile.close();
|
||||
std::cout << "Gain matching complete." << std::endl;
|
||||
|
||||
// === Create histograms ===
|
||||
TH2F *hFVB = new TH2F("hFVB", "Corrected Up+Dn vs Corrected Back;Corrected Back E;Up+Dn E",
|
||||
400, 0, 16000, 400, 0, 16000);
|
||||
TH2F *hAsym = new TH2F("hAsym", "Up vs Dn dvide corrected back;Up/Back E;Dn/Back E",
|
||||
400, 0.0, 1.0, 400, 0.0, 1.0);
|
||||
|
||||
// Fill histograms
|
||||
for (const auto &kv : dataPoints)
|
||||
{
|
||||
auto [id, ud, bk] = kv.first;
|
||||
if (!gainValid[id][ud][bk])
|
||||
continue;
|
||||
double gain = gainArray[id][ud][bk];
|
||||
|
||||
for (const auto &pr : kv.second)
|
||||
{
|
||||
double eBk, eUp, eDn;
|
||||
std::tie(eBk, eUp, eDn) = pr;
|
||||
|
||||
double updn = eUp + eDn;
|
||||
if (updn == 0)
|
||||
continue;
|
||||
|
||||
double asym = (eUp - eDn) / updn;
|
||||
double correctedBack = eBk * gain;
|
||||
|
||||
hFVB->Fill(correctedBack, updn);
|
||||
hAsym->Fill(eUp / correctedBack, eDn / correctedBack);
|
||||
}
|
||||
}
|
||||
|
||||
// Optional: save histograms to a file
|
||||
// TFile *outHist = new TFile("sx3_gainmatch_hists.root", "RECREATE");
|
||||
// hFVB->Write();
|
||||
// hAsym->Write();
|
||||
// outHist->Close();
|
||||
}
|
114
GainMatchSX3.h
Normal file
114
GainMatchSX3.h
Normal file
|
@ -0,0 +1,114 @@
|
|||
#ifndef GainMatchSX3_h
|
||||
#define GainMatchSX3_h
|
||||
|
||||
#include <TROOT.h>
|
||||
#include <TChain.h>
|
||||
#include <TFile.h>
|
||||
#include <TSelector.h>
|
||||
|
||||
#include "Armory/ClassDet.h"
|
||||
|
||||
class GainMatchSX3 : 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; //!
|
||||
|
||||
GainMatchSX3(TTree * /*tree*/ =0) : fChain(0) { }
|
||||
virtual ~GainMatchSX3() { }
|
||||
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(GainMatchSX3,0);
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
#ifdef GainMatchSX3_cxx
|
||||
void GainMatchSX3::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 GainMatchSX3::Notify(){
|
||||
|
||||
return kTRUE;
|
||||
}
|
||||
|
||||
void GainMatchSX3::SlaveBegin(TTree * /*tree*/){
|
||||
|
||||
TString option = GetOption();
|
||||
|
||||
}
|
||||
|
||||
void GainMatchSX3::SlaveTerminate(){
|
||||
|
||||
}
|
||||
|
||||
|
||||
#endif // #ifdef GainMatchSX3_cxx
|
133
MatchAndPlotCentroids.C
Normal file
133
MatchAndPlotCentroids.C
Normal 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
454
PCGainMatch.C
Normal 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();
|
||||
}
|
114
PCGainMatch.h
Normal file
114
PCGainMatch.h
Normal file
|
@ -0,0 +1,114 @@
|
|||
#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
|
|
@ -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,19 +10,24 @@ fi
|
|||
runID=$1
|
||||
timeWindow=$2
|
||||
|
||||
rawFolder=/home/tandem/Desktop/analysis/data
|
||||
rootFolder=/home/tandem/Desktop/analysis/data/root_data
|
||||
option=$3
|
||||
|
||||
rsync -a splitpole@128.186.111.223:/media/nvmeData/ANASEN27Alap/*.fsu /home/tandem/Desktop/analysis/data
|
||||
rawFolder=/home/tandem/data1/2024_09_17Fap/data
|
||||
rootFolder=../root_data
|
||||
|
||||
fileList=`\ls -1 ${rawFolder}/Run_${runID}_*.fsu`
|
||||
if [ $option -eq 0 ]; then
|
||||
|
||||
./EventBuilder ${timeWindow} 0 0 10000000 ${fileList}
|
||||
rsync -auh --info=progress2 splitpole@128.186.111.223:/media/nvmeData/2024_09_17Fap/*.fsu /home/tandem/data1/2024_09_17Fap/data
|
||||
|
||||
outFile=${rawFolder}/*${runID}*${timeWindow}.root
|
||||
fileList=`\ls -1 ${rawFolder}/*Run_${runID}_*.fsu`
|
||||
|
||||
mv -vf ${outFile} ${rootFolder}/.
|
||||
./EventBuilder ${timeWindow} 0 0 100000000 ${fileList}
|
||||
|
||||
./Mapper ${rootFolder}/*${runID}*${timeWindow}.root
|
||||
outFile=${rawFolder}/*${runID}*${timeWindow}.root
|
||||
|
||||
root "processRun.C(\"${rootFolder}/Run_${runID}_mapped.root\")"
|
||||
mv -vf ${outFile} ${rootFolder}/.
|
||||
|
||||
./Mapper ${rootFolder}/*${runID}*${timeWindow}.root
|
||||
fi
|
||||
|
||||
root "processRun.C(\"${rootFolder}/ProtonRun_${runID}_mapped.root\")"
|
||||
|
|
97
centroids.txt
Normal file
97
centroids.txt
Normal 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
89
centroids_edited.txt
Normal 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
|
26
mapping.h
26
mapping.h
|
@ -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
23
slope_intercept_cathode.txt
Normal 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
49
slope_intercept_results.txt
Normal 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
21
slope_intercept_results_anode.txt
Normal 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
|
Loading…
Reference in New Issue
Block a user