modified: .vscode/c_cpp_properties.json

modified:   .vscode/settings.json
	modified:   Analyzer.C
	new file:   Armory/#ClassPW.h#
	modified:   Armory/ClassDet.h
	modified:   Armory/ClassPW.h
	modified:   Armory/Mapper.cpp
This commit is contained in:
vs19g 2025-01-17 11:18:07 -05:00
parent 3129339647
commit 42e093b104
7 changed files with 692 additions and 487 deletions

View File

@ -1,23 +1,11 @@
{ {
"configurations": [ "configurations": [
{ {
"name": "splitpole", "name": "Hades",
"includePath": [ "includePath": [
"${workspaceFolder}/**", "${workspaceFolder}/**",
"/home/splitpole/cern/root/**" "/usr/include/x86_64-linux-gnu/qt6/**",
], "/usr/local/cern/root_v6.26.06/include/**"
"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"
], ],
"defines": [], "defines": [],
"compilerPath": "/usr/bin/gcc", "compilerPath": "/usr/bin/gcc",
@ -29,7 +17,8 @@
"name": "RyanUbuntu", "name": "RyanUbuntu",
"includePath": [ "includePath": [
"${workspaceFolder}/**", "${workspaceFolder}/**",
"/opt/root/**" "/usr/include/x86_64-linux-gnu/qt6/**",
"/opt/root/include/**"
], ],
"defines": [], "defines": [],
"compilerPath": "/usr/bin/gcc", "compilerPath": "/usr/bin/gcc",
@ -38,10 +27,11 @@
"intelliSenseMode": "linux-gcc-x64" "intelliSenseMode": "linux-gcc-x64"
}, },
{ {
"name": "RyanHome", "name": "Anasen",
"includePath": [ "includePath": [
"${workspaceFolder}/**", "${workspaceFolder}/**",
"/home/ryan/root_v6.30.06/**" "/usr/include/x86_64-linux-gnu/qt6/**",
"/opt/root/include/**"
], ],
"defines": [], "defines": [],
"compilerPath": "/usr/bin/gcc", "compilerPath": "/usr/bin/gcc",
@ -50,10 +40,26 @@
"intelliSenseMode": "linux-gcc-x64" "intelliSenseMode": "linux-gcc-x64"
}, },
{ {
"name": "Dirac", "name": "Splitpole",
"includePath": [ "includePath": [
"${workspaceFolder}/**", "${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": [], "defines": [],
"compilerPath": "/usr/bin/gcc", "compilerPath": "/usr/bin/gcc",

16
.vscode/settings.json vendored
View File

@ -1,7 +1,5 @@
{ {
"files.associations": { "files.associations": {
"NiceMatStyle.C": "cpp",
"*.rmd": "markdown",
"ryanScript.C": "cpp", "ryanScript.C": "cpp",
"ryanSelector.C": "cpp", "ryanSelector.C": "cpp",
"array": "cpp", "array": "cpp",
@ -93,8 +91,16 @@
"TrackRecon.C": "cpp", "TrackRecon.C": "cpp",
"processRuns.C": "cpp", "processRuns.C": "cpp",
"Analysis.C": "cpp", "Analysis.C": "cpp",
"Analyzer1.C": "cpp", "datastructs.h": "c",
"PCGainMatch.C": "cpp" "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"
}, },
"C_Cpp.clang_format_fallbackStyle": "{BasedonStyle: Google, IndentWidth: 2, ColumnLimit: 0}" "github-enterprise.uri": "https://fsunuc.physics.fsu.edu"
} }

View File

@ -1,16 +1,16 @@
#define Analyzer_cxx #define Analyzer_cxx
#include "Analyzer.h" #include "Analyzer.h"
#include <TCanvas.h>
#include <TCutG.h>
#include <TH2.h> #include <TH2.h>
#include <TMath.h>
#include <TStyle.h> #include <TStyle.h>
#include <algorithm> #include <TCanvas.h>
#include <utility> #include <TMath.h>
#include <utility>
#include <algorithm>
#include "Armory/ClassPW.h"
#include "Armory/ClassSX3.h" #include "Armory/ClassSX3.h"
#include "Armory/ClassPW.h"
#include "TVector3.h" #include "TVector3.h"
@ -28,43 +28,25 @@ TH2F *hqqqVpcIndex;
TH2F *hqqqVpcE; TH2F *hqqqVpcE;
TH2F *hsx3VpcE; TH2F *hsx3VpcE;
TH2F *hanVScatsum; TH2F *hanVScatsum;
TH2F *hICvsSi;
TH2F *hAnodeHits;
TH2F *hSiEvsMCPt;
TH2F *hRfvsMCPt;
TH1F *hAnodeHits1d;
TH1F *hPCMultiplicity;
TH1F *hRFtime;
TH1F *hSi;
TH1F *hSi_gated;
TH1F *hSiMCPt;
int padID = 0; int padID = 0;
SX3 sx3_contr; SX3 sx3_contr;
PW pw_contr; PW pw_contr;
PW pwinstance;
TVector3 hitPos; TVector3 hitPos;
bool HitNonZero; bool HitNonZero;
TH1F *hZProj; TH1F *hZProj;
TCutG *PCCoinc;
TCutG *alpha_cut_up; void Analyzer::Begin(TTree * /*tree*/)
TCutG *alpha_cut_down; {
TCutG *cutg;
bool inCut;
bool inCutUp;
bool inCutDown;
bool inCutG;
void Analyzer::Begin(TTree * /*tree*/) {
TString option = GetOption(); TString option = GetOption();
hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000); hsx3IndexVE = new TH2F("hsx3IndexVE", "SX3 index vs Energy; sx3 index ; Energy", 24 * 12, 0, 24 * 12, 400, 0, 5000);
hsx3IndexVE->SetNdivisions(-612, "x"); 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 = new TH2F("hqqqIndexVE", "QQQ index vs Energy; QQQ index ; Energy", 4 * 2 * 16, 0, 4 * 2 * 16, 400, 0, 5000);
hqqqIndexVE->SetNdivisions(-1204, "x"); hqqqIndexVE->SetNdivisions(-1204, "x");
hpcIndexVE = new TH2F("hpcIndexVE", "PC index vs Energy; PC index ; Energy", 2 * 24, 0, 2 * 24, 6400, 0, 30000); hpcIndexVE = new TH2F("hpcIndexVE", "PC index vs Energy; PC index ; Energy", 2 * 24, 0, 2 * 24, 400, 0, 4000);
hpcIndexVE->SetNdivisions(-1204, "x"); hpcIndexVE->SetNdivisions(-1204, "x");
hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24 * 12, 0, 24 * 12, 24 * 12, 0, 24 * 12); hsx3Coin = new TH2F("hsx3Coin", "SX3 Coincident", 24 * 12, 0, 24 * 12, 24 * 12, 0, 24 * 12);
@ -80,43 +62,23 @@ void Analyzer::Begin(TTree * /*tree*/) {
hqqqVpcIndex->SetNdivisions(-612, "x"); hqqqVpcIndex->SetNdivisions(-612, "x");
hqqqVpcIndex->SetNdivisions(-12, "y"); hqqqVpcIndex->SetNdivisions(-12, "y");
hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 400, 0, 5000, 6400, 0, 30000); hqqqVpcE = new TH2F("hqqqVpcEnergy", "qqq vs pc; qqq energy; pc energy", 400, 0, 5000, 800, 0, 16000);
hqqqVpcE->SetNdivisions(-612, "x"); hqqqVpcE->SetNdivisions(-612, "x");
hqqqVpcE->SetNdivisions(-12, "y"); hqqqVpcE->SetNdivisions(-12, "y");
hsx3VpcE = new TH2F("hsx3VpcEnergy", "sx3 vs pc; sx3 energy; pc energy", 400, 0, 5000, 6400, 0, 30000); hsx3VpcE = new TH2F("hsx3VpcEnergy", "sx3 vs pc; sx3 energy; pc energy", 400, 0, 5000, 800, 0, 16000);
hsx3VpcE->SetNdivisions(-612, "x"); hsx3VpcE->SetNdivisions(-612, "x");
hsx3VpcE->SetNdivisions(-12, "y"); hsx3VpcE->SetNdivisions(-12, "y");
hZProj = new TH1F("hZProj", "ZProjection", 600, -600, 600); hZProj = new TH1F("hZProj", "Z Projection", 200, -600, 600);
hAnodeHits1d = new TH1F("hAnodeHits1d", "Anode Hits", 24, 0, 24);
hAnodeHits = new TH2F("hAnodeHits", "Anode vs Anode Energy, Anode ID; Anode E", 24, 0, 23, 400, 0, 30000);
hPCMultiplicity = new TH1F("hPCMultiplicity", "Number of PC/Event", 40, 0, 40);
hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 6400, 0, 30000, 6400, 0, 30000);
hICvsSi = new TH2F("hICvsSi", "IC vs Si; Si E; IC E", 800, 0, 20000, 400, 0, 8000);
hSi = new TH1F("hSi", "Si E", 800, 0, 20000);
hSi_gated = new TH1F("hSi_gated", "Si E", 800, 0, 20000);
hRFtime = new TH1F("hRFtime", "Rf-MCP time(ns)", 3000, -3000, 3000);
hSiEvsMCPt = new TH2F("hSiEsMCPt", "Si E vs MCP time; SI E; MCP time", 800, 0, 20000, 3000, -3000, 3000);
hSiMCPt = new TH1F("hSiMCPt", "Si vs MCP time", 1500, -3000, 3000);
hRfvsMCPt = new TH2F("hRfvsMCPt", "RF vs MCP time; RF(ns) ; MCP time(ns)", 1000, -2000, 2000, 1000, -2000, 2000);
hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400, 0, 10000, 400, 0, 16000);
sx3_contr.ConstructGeo(); sx3_contr.ConstructGeo();
pw_contr.ConstructGeo(); pw_contr.ConstructGeo();
TFile *f1 = new TFile("PCCoinc.root");
PCCoinc = (TCutG *)f1->Get("PCCoinc");
TFile *f2 = new TFile("alpha_cut_up.root");
alpha_cut_up = (TCutG *)f2->Get("alpha_cut_up");
TFile *f3 = new TFile("alpha_cut_down.root");
alpha_cut_down = (TCutG *)f3->Get("alpha_cut_down");
TFile *f4 = new TFile("CUTG.root");
cutg = (TCutG *)f4->Get("CUTG");
// TFile *f1 = new TFile("AnCatSum.root");
// AnCatSum= (TCutG*)f1->Get("AnCatSum");
} }
Bool_t Analyzer::Process(Long64_t entry) { Bool_t Analyzer::Process(Long64_t entry)
{
// if ( entry > 100 ) return kTRUE; // if ( entry > 100 ) return kTRUE;
@ -141,12 +103,6 @@ Bool_t Analyzer::Process(Long64_t entry) {
b_pcCh->GetEntry(entry); b_pcCh->GetEntry(entry);
b_pcE->GetEntry(entry); b_pcE->GetEntry(entry);
b_pcT->GetEntry(entry); b_pcT->GetEntry(entry);
b_miscCh->GetEntry(entry);
b_miscE->GetEntry(entry);
b_miscID->GetEntry(entry);
b_miscMulti->GetEntry(entry);
b_miscT->GetEntry(entry);
b_miscTf->GetEntry(entry);
sx3.CalIndex(); sx3.CalIndex();
qqq.CalIndex(); qqq.CalIndex();
@ -158,41 +114,50 @@ Bool_t Analyzer::Process(Long64_t entry) {
// //======================= SX3 // //======================= SX3
std::vector<std::pair<int, int>> ID; // first = id, 2nd = index 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++)
if (sx3.e[i] > 50) { {
ID.push_back(std::pair<int, int>(sx3.id[i], 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++) { for (int j = i + 1; j < sx3.multi; j++)
hsx3Coin->Fill(sx3.index[i], sx3.index[j]); {
} hsx3Coin->Fill(sx3.index[i], sx3.index[j]);
}
// for( int j = 0; j < pc.multi; j++){ for (int j = 0; j < pc.multi; j++)
// hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] ); {
hsx3VpcIndex->Fill(sx3.index[i], pc.index[j]);
// if( sx3.ch[index] > 8 ){ // if( sx3.ch[index] > 8 ){
// hsx3VpcE->Fill( sx3.e[i], pc.e[j] ); // hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
// } // }
// }
} }
} }
if (ID.size() > 0) { 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::sort(ID.begin(), ID.end(), [](const std::pair<int, int> &a, const std::pair<int, int> &b)
{ return a.first < b.first; });
// printf("##############################\n"); // printf("##############################\n");
// for( size_t i = 0; i < ID.size(); i++) printf("%zu | %d %d \n", i, ID[i].first, ID[i].second ); // 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]); sx3ID.push_back(ID[0]);
bool found = false; bool found = false;
for (size_t i = 1; i < ID.size(); i++) { for (size_t i = 1; i < ID.size(); i++)
if (ID[i].first == sx3ID.back().first) { {
if (ID[i].first == sx3ID.back().first)
{
sx3ID.push_back(ID[i]); sx3ID.push_back(ID[i]);
if (sx3ID.size() >= 3) { if (sx3ID.size() >= 3)
{
found = true; found = true;
} }
} else { }
if (!found) { else
{
if (!found)
{
sx3ID.clear(); sx3ID.clear();
sx3ID.push_back(ID[i]); sx3ID.push_back(ID[i]);
} }
@ -201,30 +166,39 @@ Bool_t Analyzer::Process(Long64_t entry) {
// printf("---------- sx3ID Multi : %zu \n", sx3ID.size()); // printf("---------- sx3ID Multi : %zu \n", sx3ID.size());
if (found) { if (found)
{
int sx3ChUp, sx3ChDn, sx3ChBk; int sx3ChUp, sx3ChDn, sx3ChBk;
float sx3EUp, sx3EDn; float sx3EUp, sx3EDn;
// printf("------ sx3 ID : %d, multi: %zu\n", sx3ID[0].first, sx3ID.size()); // 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; int index = sx3ID[i].second;
// printf(" %zu | index %d | ch : %d, energy : %d \n", i, index, sx3.ch[index], sx3.e[index]); // 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] < 8)
if (sx3.ch[index] % 2 == 0) { {
if (sx3.ch[index] % 2 == 0)
{
sx3ChDn = sx3.ch[index]; sx3ChDn = sx3.ch[index];
sx3EDn = sx3.e[index]; sx3EDn = sx3.e[index];
} else { }
else
{
sx3ChUp = sx3.ch[index]; sx3ChUp = sx3.ch[index];
sx3EUp = sx3.e[index]; sx3EUp = sx3.e[index];
} }
} else { }
else
{
sx3ChBk = sx3.ch[index]; sx3ChBk = sx3.ch[index];
} }
for (int j = 0; j < pc.multi; j++) { for (int j = 0; j < pc.multi; j++)
{
// hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] ); // hsx3VpcIndex->Fill( sx3.index[i], pc.index[j] );
if (sx3.ch[i] > 8 && pc.index[j] < 24 && pc.e[j] > 50) { if (sx3.ch[index] > 8)
{
hsx3VpcE->Fill(sx3.e[i], pc.e[j]); hsx3VpcE->Fill(sx3.e[i], pc.e[j]);
// printf(" sx3 Ch: %d, pc Ch: %d , : %d\n", sx3.index[i], pc.index[j], sx3.t[i] - pc.t[j]);
// hpcIndexVE->Fill( pc.index[i], pc.e[i] ); // hpcIndexVE->Fill( pc.index[i], pc.e[i] );
} }
} }
@ -238,276 +212,196 @@ Bool_t Analyzer::Process(Long64_t entry) {
} }
// //======================= QQQ // //======================= QQQ
for (int i = 0; i < qqq.multi; i++) { for (int i = 0; i < qqq.multi; i++)
if (qqq.e[i] > 50) { {
// 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 = 0; j < pc.multi; j++){ for (int j = i + 1; j < qqq.multi; j++)
// if(pc.index[j]==4){ {
hqqqIndexVE->Fill(qqq.index[i], qqq.e[i]); for (int k = 0; k < pc.multi; k++)
// } {
// } if (pc.index[k] < 24 && pc.e[k] > 50)
for (int j = 0; j < qqq.multi; j++) { {
if (j == i) hqqqVpcE->Fill(qqq.e[i], pc.e[k]);
continue;
hqqqCoin->Fill(qqq.index[i], qqq.index[j]);
}
for (int j = 0; j < pc.multi; j++) {
if (pc.index[j] < 24 && pc.e[j] > 50) {
hqqqVpcE->Fill(qqq.e[i], pc.e[j]);
// hpcIndexVE->Fill( pc.index[i], pc.e[i] ); // hpcIndexVE->Fill( pc.index[i], pc.e[i] );
hqqqVpcIndex->Fill(qqq.index[i], pc.index[j]); hqqqVpcIndex->Fill(qqq.index[i], pc.index[j]);
} }
// }
} }
// }
// if( qqq.used[i] == true ) continue; // if( qqq.used[i] == true ) continue;
for (int j = i + 1; j < qqq.multi; j++) {
// 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] && (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])
{ // must be same detector
int chWedge = -1; int chWedge = -1;
int chRing = -1; int chRing = -1;
if (qqq.ch[i] < qqq.ch[j]) { if (qqq.ch[i] < qqq.ch[j])
chRing = qqq.ch[j] - 16; {
chWedge = qqq.ch[i]; chRing = qqq.ch[j] - 16;
} else { chWedge = qqq.ch[i];
chRing = qqq.ch[i]; }
chWedge = qqq.ch[j] - 16; else
} {
chRing = qqq.ch[i];
chWedge = qqq.ch[j] - 16;
}
// printf(" ID : %d , chWedge : %d, chRing : %d \n", qqq.id[i], chWedge, chRing); // 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 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 rho = 10. + 40. / 16. * (chRing + 0.5);
// if(qqq.e[i]>50){ // if(qqq.e[i]>50){
hqqqPolar->Fill(theta, rho); hqqqPolar->Fill(theta, rho);
// } // }
// qqq.used[i] = true; // qqq.used[i] = true;
// qqq.used[j] = true; // qqq.used[j] = true;
if (!HitNonZero) { if (!HitNonZero)
double x = rho * TMath::Cos(theta); {
double y = rho * TMath::Sin(theta); double x = rho * TMath::Cos(theta);
hitPos.SetXYZ(x, y, 23 + 75 + 30); double y = rho * TMath::Sin(theta);
HitNonZero = true; hitPos.SetXYZ(x, y, 23 + 75 + 30);
} HitNonZero = true;
} }
} }
} }
} }
// //======================= PC // //======================= PC
std::vector<std::pair<int, double>> anodeHits; ID.clear();
std::vector<std::pair<int, double>> cathodeHits; int counter = 0;
int aID = 0; std::vector<std::pair<int, double>> E;
int cID = 0; E.clear();
float cEMax = 0; for (int i = 0; i < pc.multi; i++)
int cIDMax = 0; {
float cEnextMax = 0;
int cIDnextMax = 0;
float aE = 0;
float cE = 0;
// Define the excluded SX3 and QQQ channels if (pc.e[i] > 100)
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}; ID.push_back(std::pair<int, int>(pc.id[i], i));
std::unordered_set<int> excludeQQQ = {0, 17, 109, 110, 111, 112, 113, 119, 127, 128}; if (pc.e[i] > 100)
E.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
for (int i = 0; i < pc.multi; i++) { hpcIndexVE->Fill(pc.index[i], pc.e[i]);
// for(int j=0; j<pc.multi;j++){
// if(pc.id[j]==0){
// anodeCount++;
// }
// }
if (pc.e[i] > 100 & pc.multi < 7) { for (int j = i + 1; j < pc.multi; j++)
// 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]);
// 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 );
// for (int j=0;j<sx3.multi;j++){ pwinstance.ConstructGeo();
// if(excludeSX3.find(sx3.index[j]) == excludeSX3.end()){ Coord Crossover[24][24];
TVector3 a, c, diff;
// hpcIndexVE->Fill( pc.index[i], pc.e[i] ); double a2, ac, c2, adiff, cdiff, denom, alpha, beta;
int index = 0;
for (int j = i + 1; j < pc.multi; j++) { for (int i = 0; i < pwinstance.An.size(); i++)
inCut = false; {
if (PCCoinc->IsInside(pc.index[i], pc.index[j])) { a = pwinstance.An[i].first - pwinstance.An[i].second;
inCut = true; for (int j = 0; j < pwinstance.Ca.size(); j++)
{
c = pwinstance.Ca[j].first - pwinstance.Ca[j].second;
diff = pwinstance.An[i].first - pwinstance.Ca[j].first;
a2 = a.Dot(a);
ac = a.Dot(c);
c2 = c.Dot(c);
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].x = pwinstance.An[i].first.X() + alpha * a.X();
Crossover[i][j].y = pwinstance.An[i].first.Y() + alpha * a.Y();
Crossover[i][j].z = pwinstance.An[i].first.Z() + alpha * a.Z();
if (i == 23)
{
if (abs(i - j) < 7 || abs(i - j) > 17)
{
if (alpha < 0 && alpha > -1)
{
printf("Anode and cathode indices and coord : %d %d %f %f %f %f\n", i, j, pwinstance.Ca[j].first.X(), pwinstance.Ca[j].first.Y(), pwinstance.Ca[j].first.Z(), alpha);
printf("Crossover wires, points and alpha are : %f %f %f %f \n", Crossover[i][j].x, Crossover[i][j].y, Crossover[i][j].z, alpha);
}
} }
// hpcCoin->Fill( pc.index[i], pc.index[j]);
} }
// if(pc.e[i]>100){
if (pc.index[i] < 24) {
anodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
// anodeCount++;
} else if (pc.index[i] >= 24) {
cathodeHits.push_back(std::pair<int, double>(pc.index[i], pc.e[i]));
}
// }
// }
// }
// hpcIndexVE->Fill( pc.index[i], pc.e[i] );
} }
} }
hPCMultiplicity->Fill(pc.multi);
float aESum = 0; if (E.size() >= 3)
float cESum = 0; {
if (anodeHits.size() == 1 && cathodeHits.size() >= 1) {
inCutDown = false; int aID = 0;
inCutUp = false; int cID = 0;
for (const auto &anode : anodeHits) {
// for(int l=0; l<sx3.multi; l++){ float aE = 0;
// if (sx3.index[l]==80){ float cE = 0;
// if( ID[0].first < 1 ) {
aID = anode.first; // aID = pc.ch[ID[0].second];
aE = anode.second; // cID = pc.ch[ID[1].second];
aESum += aE; // }else{
// printf("aID : %d, aE : %f, cE : %f\n", aID, aE, cE); // cID = pc.ch[ID[0].second];
} // aID = pc.ch[ID[1].second];
for (const auto &cathode : cathodeHits) {
cID = cathode.first;
cE = cathode.second;
cESum += cE;
if (cE > cEMax) {
cEMax = cE;
cIDMax = cID;
}
if (cE > cEnextMax && cE < cEMax) {
cEnextMax = cE;
cIDnextMax = cID;
}
}
if (alpha_cut_down->IsInside(aE, cESum)) {
inCutDown = true;
}
if (alpha_cut_up->IsInside(aE, cESum)) {
inCutUp = true;
}
// if (inCutUp)
// {
for (int i = 0; i < pc.multi; i++) {
for (int j = i + 1; j < pc.multi; j++) {
hpcCoin->Fill(pc.index[i], pc.index[j]);
hpcIndexVE->Fill(pc.index[i], pc.e[i]);
}
}
// } // }
// if (inCut) { // printf("anode= %d, cathode = %d\n", aID, cID);
hanVScatsum->Fill(aE, cESum);
hAnodeHits->Fill(aID, aE);
hAnodeHits1d->Fill(anodeHits.size());
// }
// }
}
// Miscellaneous channels including the Lollipop IC and Si detectors and hot needle IC for (int k = 0; k < qqq.multi; k++)
// Misc ch 0,1, 2, 3, 4 in order are the LIC, LSi, HNIC-difference, MCP, and Rf {
bool timing = false; if (qqq.index[k] == 75 && pc.index[k] == 2 && pc.e[k] > 100)
inCutG = false; {
double SiE = 0;
double SiT = 0;
double MCPt = 0;
double MCPE = 0;
double Rft = 0;
double ICt = 0;
double ICe = 0;
double SiCFDt = 0;
for (int i = 0; i < misc.multi; i++) {
// if (misc.ch[i] == 1 && misc.e[i] > 10000 && misc.e[i] < 15000) {
// if(misc.e[i] > 7500 && misc.e[i]<15000) hSi->Fill(misc.e[i]);
if (misc.ch[i] == 1) { int multi_an = 0;
// hSi->Fill(misc.e[i]); for (int l = 0; l < E.size(); l++)
SiE = misc.e[i]; {
SiT = misc.t[i] + misc.tf[i] * 4. / 1000; if (E[l].first < 24)
// hSi->Fill(misc.e[i]); {
} multi_an++;
if (misc.ch[i] == 2) { }
ICt = misc.t[i] + misc.tf[i] * 4. / 1000; }
ICe = misc.e[i];
hSi->Fill(misc.e[i]); if (multi_an >= 1)
} {
if (misc.ch[i] == 3) { for (int l = 0; l < E.size(); l++)
// only analyze the first MCP in any event {
if (MCPt == 0) { if (E[l].first < 24 && E[l].first != 19 && E[l].first != 12)
MCPt = misc.t[i] + misc.tf[i] * 4. / 1000; {
MCPE = misc.e[i]; aE = E[l].second;
} }
} else if (E[l].first > 24)
if (misc.ch[i] == 4) { {
// only analyze the first RF in any event cE = E[l].second;
if (Rft == 0) { }
Rft = misc.t[i] + misc.tf[i] * 4. / 1000; }
}
}
if (misc.ch[i] == 5) {
if (SiCFDt == 0) {
SiCFDt = misc.t[i] + misc.tf[i] * 4. / 1000;
}
}
// hSiEvsMCPt1->Fill(SiE, Rft-MCPt);
// hSiEvsMCPt->Fill(ICe, MCPt - Rft);
if (MCPt != 0 && Rft != 0) {
// if (SiE > 10200 && SiE < 12200) {
// hRfvsMCPt->Fill(Rft - ICt, MCPt - ICt);
hSiMCPt->Fill(MCPt - ICt);
// if(misc.ch[i] == 2 && misc.e[i] > 1000 && misc.e[i]<2000)
hRFtime->Fill(Rft - ICt);
// }
// printf("RF time : %lld %lld %lld\n", Rft, MCPt, (MCPt - Rft));
// }
}
// inCutG = true;
// if (misc.ch[i] == 1) hSi->Fill(misc.e[i]);
// for (int j = 0; j < qqq.multi; j++) {
// if (pc.id[j] == 0) {
hRfvsMCPt->Fill(Rft-ICt, MCPt -ICt);
hSiEvsMCPt->Fill(ICe, MCPt - ICt);
// }
// }
for (int j = i + 1; j < misc.multi; j++) {
// if (cutg->IsInside(misc.e[i], misc.e[j])) {
// inCutG = true;
// })
if (misc.ch[j] == 4 && misc.ch[i] == 3) {
// hRFtime->Fill(misc.t[j]*1. + misc.tf[j] * 4. / 1000 - (misc.t[i]*1. + misc.tf[i] * 4. / 1000));
if (misc.t[j] + misc.tf[j] * 4. / 1000 - (misc.t[i] + misc.tf[i] * 4. / 1000) > 20 && misc.t[j] + misc.tf[j] * 4. / 1000 - (misc.t[i] + misc.tf[i] * 4. / 1000) < 100) {
timing = true;
} }
// printf("RF time : %lld %lld %lld %lld %lld\n", misc.t[i], misc.t[j], misc.tf[i], misc.tf[j], (misc.t[j]*1000 + misc.tf[j]*4 - (misc.t[i]*1000 + misc.tf[i]*4)));
} }
} }
hanVScatsum->Fill(aE, cE);
// for (int j = i + 1; j < misc.multi; j++) { if (ID[0].first < 1)
if (timing == true) { {
// hICvsSi->Fill(misc.e[i], misc.e[j]); aID = pc.ch[ID[0].second];
if (misc.ch[i] == 1) { cID = pc.ch[ID[1].second];
hSi_gated->Fill(misc.e[i]); }
// } else
} {
// } cID = pc.ch[ID[0].second];
aID = pc.ch[ID[1].second];
} }
}
if (HitNonZero) { if (HitNonZero)
// pw_contr.CalTrack1( hitPos, aID, cIDMax, cIDnextMax, cEMax, cEnextMax,1); {
pw_contr.CalTrack(hitPos, aID, cID);
// pw_contr.CalTrack(hitPos, aID, cID); hZProj->Fill(pw_contr.GetZ0());
hZProj->Fill(pw_contr.GetZ0()); }
} }
// ########################################################### Track constrcution // ########################################################### Track constrcution
@ -517,18 +411,16 @@ Bool_t Analyzer::Process(Long64_t entry) {
return kTRUE; return kTRUE;
} }
void Analyzer::Terminate() { void Analyzer::Terminate()
{
gStyle->SetOptStat("neiou"); gStyle->SetOptStat("neiou");
TCanvas *canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000); TCanvas *canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000);
// TCanvas *a = new TCanvas("aANASEN", "ANASEN", 800, 600);
canvas->Divide(3, 3); canvas->Divide(3, 3);
// hRFtime->Draw();
// TCanvas *b = new TCanvas("bANASEN", "ANASEN", 800, 600);
// // hICvsSi->Draw("colz");
// hSi->Draw();
// =============================================== pad-1 // hsx3VpcIndex->Draw("colz");
//=============================================== pad-1
padID++; padID++;
canvas->cd(padID); canvas->cd(padID);
canvas->cd(padID)->SetGrid(1); canvas->cd(padID)->SetGrid(1);
@ -598,6 +490,4 @@ void Analyzer::Terminate() {
canvas->cd(padID)->SetGrid(1); canvas->cd(padID)->SetGrid(1);
// hZProj->Draw(); // hZProj->Draw();
hanVScatsum->Draw("colz"); hanVScatsum->Draw("colz");
// hAnodeHits->Draw("colz");
// // hAnodeMultiplicity->Draw();
} }

283
Armory/#ClassPW.h# Normal file
View File

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

View File

@ -14,7 +14,7 @@ public:
unsigned short ch[MAXMULTI]; unsigned short ch[MAXMULTI];
unsigned short e[MAXMULTI]; unsigned short e[MAXMULTI];
unsigned long long t[MAXMULTI]; unsigned long long t[MAXMULTI];
unsigned long long tf[MAXMULTI];
unsigned short sn[MAXMULTI]; unsigned short sn[MAXMULTI];
unsigned short digiCh[MAXMULTI]; unsigned short digiCh[MAXMULTI];
@ -28,7 +28,6 @@ public:
ch[i] = 0; ch[i] = 0;
e[i] = 0; e[i] = 0;
t[i] = 0; t[i] = 0;
tf[i] = 0;
index[i] = 0; index[i] = 0;
sn[i] = 0; sn[i] = 0;
digiCh[i] = 0; digiCh[i] = 0;

View File

@ -6,14 +6,16 @@
#include <TVector3.h> #include <TVector3.h>
#include <TRandom.h> #include <TRandom.h>
struct PWHitInfo{ struct PWHitInfo
std::pair<short, short> nearestWire; // anode, cathode {
std::pair<short, short> nearestWire; // anode, cathode
std::pair<double, double> nearestDist; // 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 std::pair<double, double> nextNearestDist; // anode, cathode
void Clear(){ void Clear()
{
nearestWire.first = -1; nearestWire.first = -1;
nearestWire.second = -1; nearestWire.second = -1;
nearestDist.first = 999999999; nearestDist.first = 999999999;
@ -25,40 +27,56 @@ struct PWHitInfo{
} }
}; };
//!######################################################## struct Coord
class PW{ // proportional wire {
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
}
};
//! ########################################################
class PW
{ // proportional wire
public: public:
PW(){ ClearHitInfo();}; PW() { ClearHitInfo(); };
~PW(){}; ~PW() {};
PWHitInfo GetHitInfo() const {return hitInfo;} PWHitInfo GetHitInfo() const { return hitInfo; }
std::pair<short, short> GetNearestID() const {return hitInfo.nearestWire;} std::pair<short, short> GetNearestID() const { return hitInfo.nearestWire; }
std::pair<double, double> GetNearestDistance() const {return hitInfo.nearestDist;} std::pair<double, double> GetNearestDistance() const { return hitInfo.nearestDist; }
std::pair<short, short> Get2ndNearestID() const {return hitInfo.nextNearestWire;} std::pair<short, short> Get2ndNearestID() const { return hitInfo.nextNearestWire; }
std::pair<double, double> Get2ndNearestDistance() const {return hitInfo.nextNearestDist;} std::pair<double, double> Get2ndNearestDistance() const { return hitInfo.nextNearestDist; }
TVector3 GetTrackPos() const {return trackPos;} std::vector<std::pair<TVector3, TVector3>> An; // the anode wire position vector in space
TVector3 GetTrackVec() const {return trackVec;} std::vector<std::pair<TVector3, TVector3>> Ca; // the cathode wire position vector in space
double GetTrackTheta() const {return trackVec.Theta();}
double GetTrackPhi() const {return trackVec.Phi();} TVector3 GetTrackPos() const { return trackPos; }
TVector3 GetTrackVec() const { return trackVec; }
double GetTrackTheta() const { return trackVec.Theta(); }
double GetTrackPhi() const { return trackVec.Phi(); }
double GetZ0(); double GetZ0();
int GetNumWire() const {return nWire;} int GetNumWire() const { return nWire; }
double GetDeltaAngle() const {return dAngle;} double GetDeltaAngle() const { return dAngle; }
double GetAnodeLength() const {return anodeLength;} double GetAnodeLength() const { return anodeLength; }
double GetCathodeLength() const {return cathodeLength;} double GetCathodeLength() const { return cathodeLength; }
TVector3 GetAnodeDn(short id) const {return An[id].first;} TVector3 GetAnodeDn(short id) const { return An[id].first; }
TVector3 GetAnodeUp(short id) const {return An[id].second;} TVector3 GetAnodeUp(short id) const { return An[id].second; }
TVector3 GetCathodeDn(short id) const {return Ca[id].first;} TVector3 GetCathodeDn(short id) const { return Ca[id].first; }
TVector3 GetCathodeUp(short id) const {return Ca[id].second;} TVector3 GetCathodeUp(short id) const { return Ca[id].second; }
TVector3 GetAnodneMid(short id) const {return (An[id].first + An[id].second) * 0.5; } 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 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();} 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; } 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 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();} double GetCathodePhi(short id) const { return (Ca[id].first - Ca[id].second).Phi(); }
void ClearHitInfo(); void ClearHitInfo();
void ConstructGeo(); void ConstructGeo();
@ -66,20 +84,20 @@ public:
void CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, 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, PWHitInfo hitInfo, double sigmaA = 0, double sigmaC = 0, bool verbose = false);
void Print(){ void Print()
{
printf(" The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nearestWire.first, printf(" The nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nearestWire.first,
hitInfo.nearestDist.first, hitInfo.nearestDist.first,
hitInfo.nearestWire.second, hitInfo.nearestWire.second,
hitInfo.nearestDist.second); hitInfo.nearestDist.second);
printf(" The 2nd nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nextNearestWire.first, printf(" The 2nd nearest | Anode: %2d(%5.2f) Cathode: %2d(%5.2f)\n", hitInfo.nextNearestWire.first,
hitInfo.nextNearestDist.first, hitInfo.nextNearestDist.first,
hitInfo.nextNearestWire.second, hitInfo.nextNearestWire.second,
hitInfo.nextNearestDist.second); hitInfo.nextNearestDist.second);
} }
private: private:
PWHitInfo hitInfo; PWHitInfo hitInfo;
TVector3 trackPos; TVector3 trackPos;
@ -87,7 +105,7 @@ private:
const int nWire = 24; const int nWire = 24;
const int wireShift = 3; const int wireShift = 3;
const float zLen = 380; //mm const float zLen = 380; // mm
const float radiusA = 37; const float radiusA = 37;
const float radiusC = 43; const float radiusC = 43;
@ -95,23 +113,25 @@ private:
double anodeLength; double anodeLength;
double cathodeLength; double cathodeLength;
std::vector<std::pair<TVector3,TVector3>> An; // the anode 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 // 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 na = a1 - a2;
TVector3 nb = b1 - b2; TVector3 nb = b1 - b2;
TVector3 nd = (na.Cross(nb)).Unit(); 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(); hitInfo.Clear();
} }
inline void PW::ConstructGeo(){ inline void PW::ConstructGeo()
{
An.clear(); An.clear();
Ca.clear(); Ca.clear();
@ -119,108 +139,132 @@ inline void PW::ConstructGeo(){
std::pair<TVector3, TVector3> p1; // anode std::pair<TVector3, TVector3> p1; // anode
std::pair<TVector3, TVector3> q1; // cathode std::pair<TVector3, TVector3> q1; // cathode
//anode and cathode start at pos-Y axis and count in right-Hand // anode and cathode start at pos-Y axis and count in right-Hand
//anode wire shift is right-hand. // anode wire shift is right-hand.
//cathode wire shift is left-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 // Anode rotate right-hand
p1.first.SetXYZ( radiusA * TMath::Cos( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()), p1.first.SetXYZ(radiusA * TMath::Cos(TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
radiusA * TMath::Sin( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()), radiusA * TMath::Sin(TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
zLen/2); zLen / 2);
p1.second.SetXYZ( radiusA * TMath::Cos( TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()), p1.second.SetXYZ(radiusA * TMath::Cos(TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()),
radiusA * TMath::Sin( TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()), radiusA * TMath::Sin(TMath::TwoPi() / nWire * (i + wireShift) + TMath::PiOver2()),
-zLen/2); -zLen / 2);
An.push_back(p1); An.push_back(p1);
// Cathod rotate left-hand // Cathod rotate left-hand
q1.first.SetXYZ( radiusC * TMath::Cos( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()), q1.first.SetXYZ(radiusC * TMath::Cos(TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
radiusC * TMath::Sin( TMath::TwoPi() / nWire * (i) + TMath::PiOver2()), radiusC * TMath::Sin(TMath::TwoPi() / nWire * (i) + TMath::PiOver2()),
zLen/2); zLen / 2);
q1.second.SetXYZ( radiusC * TMath::Cos( TMath::TwoPi() / nWire * (i - wireShift) + TMath::PiOver2()), q1.second.SetXYZ(radiusC * TMath::Cos(TMath::TwoPi() / nWire * (i - wireShift) + TMath::PiOver2()),
radiusC * TMath::Sin( TMath::TwoPi() / nWire * (i - wireShift) + TMath::PiOver2()), radiusC * TMath::Sin(TMath::TwoPi() / nWire * (i - wireShift) + TMath::PiOver2()),
-zLen/2); -zLen / 2);
Ca.push_back(q1); Ca.push_back(q1);
} }
dAngle = wireShift * TMath::TwoPi() / nWire; dAngle = wireShift * TMath::TwoPi() / nWire;
anodeLength = TMath::Sqrt( zLen*zLen + TMath::Power(2* radiusA * 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) ); 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(); hitInfo.Clear();
double phi = direction.Phi(); double phi = direction.Phi();
for( int i = 0; i < nWire; i++){ for (int i = 0; i < nWire; i++)
{
double disA = 99999999; 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(); 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()); // 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)
if( phi < 0 && phiS > phiL ) phiS = phiS - TMath::TwoPi(); phiL = phiL + TMath::TwoPi();
if (phi < 0 && phiS > phiL)
phiS = phiS - TMath::TwoPi();
if( phiS < phi && phi < phiL) { if (phiS < phi && phi < phiL)
disA = Distance( pos, pos + direction, An[i].first, An[i].second); {
if( disA < hitInfo.nearestDist.first ){ disA = Distance(pos, pos + direction, An[i].first, An[i].second);
if (disA < hitInfo.nearestDist.first)
{
hitInfo.nearestDist.first = disA; hitInfo.nearestDist.first = disA;
hitInfo.nearestWire.first = i; hitInfo.nearestWire.first = i;
} }
} }
double disC = 99999999; double disC = 99999999;
phiS = Ca[i].second.Phi()- TMath::PiOver4(); phiS = Ca[i].second.Phi() - TMath::PiOver4();
phiL = Ca[i].first.Phi() + TMath::PiOver4(); phiL = Ca[i].first.Phi() + TMath::PiOver4();
// printf("C%2d: %f %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg()); // 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)
if( phi < 0 && phiS > phiL ) phiS = phiS - TMath::TwoPi(); phiL = phiL + TMath::TwoPi();
if (phi < 0 && phiS > phiL)
phiS = phiS - TMath::TwoPi();
if(phiS < phi && phi < phiL) { if (phiS < phi && phi < phiL)
disC = Distance( pos, pos + direction, Ca[i].first, Ca[i].second); {
if( disC < hitInfo.nearestDist.second ){ disC = Distance(pos, pos + direction, Ca[i].first, Ca[i].second);
if (disC < hitInfo.nearestDist.second)
{
hitInfo.nearestDist.second = disC; hitInfo.nearestDist.second = disC;
hitInfo.nearestWire.second = i; 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 //==== find the 2nd nearest wire
short anode1 = hitInfo.nearestWire.first; short anode1 = hitInfo.nearestWire.first;
short aaa1 = anode1 - 1; if( aaa1 < 0 ) aaa1 += nWire; short aaa1 = anode1 - 1;
if (aaa1 < 0)
aaa1 += nWire;
short aaa2 = (anode1 + 1) % nWire; short aaa2 = (anode1 + 1) % nWire;
double haha1 = Distance( pos, pos + direction, An[aaa1].first, An[aaa1].second); double haha1 = Distance(pos, pos + direction, An[aaa1].first, An[aaa1].second);
double haha2 = Distance( pos, pos + direction, An[aaa2].first, An[aaa2].second); double haha2 = Distance(pos, pos + direction, An[aaa2].first, An[aaa2].second);
if( haha1 < haha2){ if (haha1 < haha2)
{
hitInfo.nextNearestWire.first = aaa1; hitInfo.nextNearestWire.first = aaa1;
hitInfo.nextNearestDist.first = haha1; hitInfo.nextNearestDist.first = haha1;
}else{ }
else
{
hitInfo.nextNearestWire.first = aaa2; hitInfo.nextNearestWire.first = aaa2;
hitInfo.nextNearestDist.first = haha2; hitInfo.nextNearestDist.first = haha2;
} }
short cathode1 = hitInfo.nearestWire.second; short cathode1 = hitInfo.nearestWire.second;
short ccc1 = cathode1 - 1; if( ccc1 < 0 ) ccc1 += nWire; short ccc1 = cathode1 - 1;
if (ccc1 < 0)
ccc1 += nWire;
short ccc2 = (cathode1 + 1) % nWire; short ccc2 = (cathode1 + 1) % nWire;
haha1 = Distance( pos, pos + direction, Ca[ccc1].first, Ca[ccc1].second); haha1 = Distance(pos, pos + direction, Ca[ccc1].first, Ca[ccc1].second);
haha2 = Distance( pos, pos + direction, Ca[ccc2].first, Ca[ccc2].second); haha2 = Distance(pos, pos + direction, Ca[ccc2].first, Ca[ccc2].second);
if( haha1 < haha2){ if (haha1 < haha2)
{
hitInfo.nextNearestWire.second = ccc1; hitInfo.nextNearestWire.second = ccc1;
hitInfo.nextNearestDist.second = haha1; hitInfo.nextNearestDist.second = haha1;
}else{ }
else
{
hitInfo.nextNearestWire.second = ccc2; hitInfo.nextNearestWire.second = ccc2;
hitInfo.nextNearestDist.second = haha2; 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; trackPos = sx3Pos;
@ -230,11 +274,12 @@ 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 // if the handiness of anode and cathode revered, it should be n2 cross n1
trackVec = (n2.Cross(n1)).Unit(); 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){ inline void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA, double sigmaC, bool verbose)
{
trackPos = sx3Pos; trackPos = sx3Pos;
@ -266,19 +311,19 @@ inline void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA, dou
// if the handiness of anode and cathode revered, it should be n2 cross n1 // if the handiness of anode and cathode revered, it should be n2 cross n1
trackVec = (n2.Cross(n1)).Unit(); 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 double PW::GetZ0(){ inline double PW::GetZ0()
{
double x = trackPos.X(); double x = trackPos.X();
double y = trackPos.Y(); double y = trackPos.Y();
double rho = TMath::Sqrt(x*x + y*y); double rho = TMath::Sqrt(x * x + y * y);
double theta = trackVec.Theta(); double theta = trackVec.Theta();
return trackPos.Z() - rho / TMath::Tan(theta); return trackPos.Z() - rho / TMath::Tan(theta);
} }
#endif #endif

View File

@ -67,7 +67,6 @@ int main(int argc, char **argv){
Det sx3; Det sx3;
Det qqq; Det qqq;
Det pc ; Det pc ;
Det misc;
printf(" Raw root file : %s\n", inFileName.c_str()); printf(" Raw root file : %s\n", inFileName.c_str());
printf(" Run : %03d\n", run); printf(" Run : %03d\n", run);
@ -100,14 +99,6 @@ int main(int argc, char **argv){
newTree->Branch("pcE", &pc.e, "pcEnergy[pcMulti]/s"); newTree->Branch("pcE", &pc.e, "pcEnergy[pcMulti]/s");
newTree->Branch("pcT", &pc.t, "pcTime[pcMulti]/l"); newTree->Branch("pcT", &pc.t, "pcTime[pcMulti]/l");
newTree->Branch("miscMulti", &misc.multi, "miscMulti/s");
newTree->Branch("miscID", &misc.id, "miscID[miscMulti]/s");
newTree->Branch("miscCh", &misc.ch, "miscCh[miscMulti]/s");
newTree->Branch("miscE", &misc.e, "miscEnergy[miscMulti]/s");
newTree->Branch("miscT", &misc.t, "miscTime[miscMulti]/l");
newTree->Branch("miscF", &misc.tf, "miscFineTime[miscMulti]/l");
///================== looping old tree and apply mapping ///================== looping old tree and apply mapping
//clock //clock
@ -121,12 +112,8 @@ int main(int argc, char **argv){
sx3.multi = 0; sx3.multi = 0;
qqq.multi = 0; qqq.multi = 0;
pc.multi = 0; pc.multi = 0;
misc.multi=0;
sx3.Clear();
qqq.Clear(); qqq.Clear();
pc.Clear();
misc.Clear();
for( unsigned int i = 0; i < multi; i++){ for( unsigned int i = 0; i < multi; i++){
@ -173,17 +160,6 @@ int main(int argc, char **argv){
pc.t[pc.multi] = e_t[i]; pc.t[pc.multi] = e_t[i];
pc.multi ++; pc.multi ++;
} }
//=================================== misc
if( 30000 <= ID && ID < 40000 ) {
misc.id[misc.multi] = (ID - 30000) / 100;
misc.ch[misc.multi] = (ID - 30000) % 100;
misc.e[misc.multi] = e[i];
misc.t[misc.multi] = e_t[i];
misc.tf[misc.multi] = e_f[i];
// if( ID == 30002 || ID == 30004 ) printf("sn : %d ch: %2d | gID %3d | ID %6d | e_f : %d\n", sn[i], ch[i], globalCh, ID, e_f[i]);
misc.multi ++;
}
} }
saveFile->cd(); //set focus on this file saveFile->cd(); //set focus on this file