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

Compare commits

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

No commits in common. "master" and "devel_ruchi" have entirely different histories.
master ... devel_ruch

5 changed files with 285 additions and 430 deletions
Show Stats Expand all files Collapse all files

1
.gitignore vendored
View File

@ -10,4 +10,3 @@ AnasenMS
data/ data/
data_proton/ data_proton/
root_data/

15
.vscode/settings.json vendored
View File

@ -90,17 +90,6 @@
"processRun.C": "cpp", "processRun.C": "cpp",
"TrackRecon.C": "cpp", "TrackRecon.C": "cpp",
"processRuns.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"
},
"github-enterprise.uri": "https://fsunuc.physics.fsu.edu"
} }

254
Analyzer.C
View File

@ -32,22 +32,18 @@ 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;
void Analyzer::Begin(TTree * /*tree*/) 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->SetNdivisions( -1204, "x");
hqqqIndexVE = new TH2F("hqqqIndexVE", "QQQ index vs Energy; QQQ index ; Energy", 4 * 2 * 16, 0, 4 * 2 * 16, 400, 0, 5000); hpcIndexVE = new TH2F("hpcIndexVE", "PC index vs Energy; PC index ; Energy", 2*24, 0, 2*24, 400, 0, 4000); hpcIndexVE->SetNdivisions( -1204, "x");
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); 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); hqqqCoin = new TH2F("hqqqCoin", "QQQ Coincident", 4*2*16, 0, 4*2*16, 4*2*16, 0, 4*2*16);
@ -75,10 +71,13 @@ void Analyzer::Begin(TTree * /*tree*/)
hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 10000, 400, 0 , 16000); hanVScatsum = new TH2F("hanVScatsum", "Anode vs Cathode Sum; Anode E; Cathode E", 400,0 , 10000, 400, 0 , 16000);
sx3_contr.ConstructGeo(); sx3_contr.ConstructGeo();
pw_contr.ConstructGeo(); pw_contr.ConstructGeo();
} }
Bool_t Analyzer::Process(Long64_t entry)
{
Bool_t Analyzer::Process(Long64_t entry){
// if ( entry > 100 ) return kTRUE; // if ( entry > 100 ) return kTRUE;
@ -114,19 +113,16 @@ 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 ++){
{
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] );
@ -134,30 +130,25 @@ Bool_t Analyzer::Process(Long64_t entry)
} }
} }
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) std::sort(ID.begin(), ID.end(), [](const std::pair<int, int> & a, const std::pair<int, int> & b) {
{ return a.first < b.first; }); 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{
else if( !found ){
{
if (!found)
{
sx3ID.clear(); sx3ID.clear();
sx3ID.push_back(ID[i]); sx3ID.push_back(ID[i]);
} }
@ -166,38 +157,29 @@ 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[index] > 8) if( sx3.ch[index] > 8 ){
{
hsx3VpcE->Fill( sx3.e[i], pc.e[j] ); hsx3VpcE->Fill( sx3.e[i], pc.e[j] );
// hpcIndexVE->Fill( pc.index[i], pc.e[i] ); // hpcIndexVE->Fill( pc.index[i], pc.e[i] );
} }
@ -209,50 +191,44 @@ Bool_t Analyzer::Process(Long64_t entry)
HitNonZero = true; HitNonZero = true;
// hitPos.Print(); // hitPos.Print();
} }
} }
// //======================= QQQ // //======================= QQQ
for (int i = 0; i < qqq.multi; i++) for( int i = 0; i < qqq.multi; i ++){
{
// for( int j = 0; j < pc.multi; j++){ // for( int j = 0; j < pc.multi; j++){
// if(pc.index[j]==4){ // if(pc.index[j]==4){
hqqqIndexVE->Fill( qqq.index[i], qqq.e[i] ); hqqqIndexVE->Fill( qqq.index[i], qqq.e[i] );
// } // }
// } // }
for (int j = 0; j < qqq.multi; j++) for( int j = 0; j < qqq.multi; j++){
{ if ( j == i ) continue;
if (j == i)
continue;
hqqqCoin->Fill( qqq.index[i], qqq.index[j]); hqqqCoin->Fill( qqq.index[i], qqq.index[j]);
} }
for (int j = i + 1; j < qqq.multi; j++)
{ for( int j = i + 1; j < qqq.multi; j++){
for (int k = 0; k < pc.multi; k++) for( int k = 0; k < pc.multi; k++){
{ if(pc.index[k]<24 && pc.e[k]>50 ){
if (pc.index[k] < 24 && pc.e[k] > 50)
{
hqqqVpcE->Fill( qqq.e[i], pc.e[k] ); hqqqVpcE->Fill( qqq.e[i], pc.e[k] );
// 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;
//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]) if( qqq.id[i] == qqq.id[j] ){ // must be same detector
{ // 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; chRing = qqq.ch[j] - 16;
chWedge = qqq.ch[i]; chWedge = qqq.ch[i];
} }else{
else
{
chRing = qqq.ch[i]; chRing = qqq.ch[i];
chWedge = qqq.ch[j] - 16; chWedge = qqq.ch[j] - 16;
} }
@ -267,8 +243,7 @@ Bool_t Analyzer::Process(Long64_t entry)
// 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 x = rho * TMath::Cos(theta);
double y = rho * TMath::Sin(theta); double y = rho * TMath::Sin(theta);
hitPos.SetXYZ(x, y, 23 + 75 + 30); hitPos.SetXYZ(x, y, 23 + 75 + 30);
@ -276,6 +251,8 @@ Bool_t Analyzer::Process(Long64_t entry)
} }
} }
} }
} }
// //======================= PC // //======================= PC
@ -283,62 +260,22 @@ Bool_t Analyzer::Process(Long64_t entry)
int counter=0; int counter=0;
std::vector<std::pair<int, double>> E; std::vector<std::pair<int, double>> E;
E.clear(); E.clear();
for (int i = 0; i < pc.multi; i++) for( int i = 0; i < pc.multi; i ++){
{
if (pc.e[i] > 100) if( pc.e[i] > 100 ) ID.push_back(std::pair<int, int>(pc.id[i], i));
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]));
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] ); hpcIndexVE->Fill( pc.index[i], pc.e[i] );
for (int j = i + 1; j < pc.multi; j++) 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( size_t i = 0; i < E.size(); i++) printf("%zu | %d %d \n", i, E[i].first, E[i].second );
pwinstance.ConstructGeo(); if( E.size()>=3 ){
Coord Crossover[24][24];
TVector3 a, c, diff;
double a2, ac, c2, adiff, cdiff, denom, alpha, beta;
int index = 0;
for (int i = 0; i < pwinstance.An.size(); i++)
{
a = pwinstance.An[i].first - pwinstance.An[i].second;
for (int j = 0; j < pwinstance.Ca.size(); j++)
{
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);
}
}
}
}
}
if (E.size() >= 3)
{
int aID = 0; int aID = 0;
int cID = 0; int cID = 0;
@ -354,65 +291,58 @@ Bool_t Analyzer::Process(Long64_t entry)
// } // }
// printf("anode= %d, cathode = %d\n", aID, cID); // printf("anode= %d, cathode = %d\n", aID, cID);
for (int k = 0; k < qqq.multi; k++) for( int k = 0; k < qqq.multi; k++){
{ if(qqq.index[k]==75 && pc.index[k]==2 && pc.e[k]>100){
if (qqq.index[k] == 75 && pc.index[k] == 2 && pc.e[k] > 100)
{
int multi_an =0; int multi_an =0;
for (int l = 0; l < E.size(); l++) for(int l=0;l<E.size();l++){
{ if(E[l].first<24){
if (E[l].first < 24)
{
multi_an++; multi_an++;
} }
} }
if (multi_an >= 1) if(multi_an>=1){
{ for(int l=0;l<E.size();l++){
for (int l = 0; l < E.size(); l++) if(E[l].first<24 && E[l].first!=19 && E[l].first!=12){
{
if (E[l].first < 24 && E[l].first != 19 && E[l].first != 12)
{
aE = E[l].second; aE = E[l].second;
} }else if(E[l].first>24){
else if (E[l].first > 24)
{
cE = E[l].second; cE = E[l].second;
} }
} }
} }
} }
} }
hanVScatsum->Fill(aE,cE); hanVScatsum->Fill(aE,cE);
if (ID[0].first < 1) if( ID[0].first < 1 ) {
{
aID = pc.ch[ID[0].second]; aID = pc.ch[ID[0].second];
cID = pc.ch[ID[1].second]; cID = pc.ch[ID[1].second];
} }else{
else
{
cID = pc.ch[ID[0].second]; cID = pc.ch[ID[0].second];
aID = pc.ch[ID[1].second]; aID = pc.ch[ID[1].second];
} }
if (HitNonZero) if( HitNonZero){
{
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
//############################## DO THE KINEMATICS //############################## DO THE KINEMATICS
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);
@ -421,61 +351,45 @@ void Analyzer::Terminate()
//hsx3VpcIndex->Draw("colz"); //hsx3VpcIndex->Draw("colz");
//=============================================== pad-1 //=============================================== pad-1
padID++; padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
canvas->cd(padID);
canvas->cd(padID)->SetGrid(1);
hsx3IndexVE->Draw("colz"); hsx3IndexVE->Draw("colz");
//=============================================== pad-2 //=============================================== pad-2
padID++; padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
canvas->cd(padID);
canvas->cd(padID)->SetGrid(1);
hqqqIndexVE->Draw("colz"); hqqqIndexVE->Draw("colz");
//=============================================== pad-3 //=============================================== pad-3
padID++; padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
canvas->cd(padID);
canvas->cd(padID)->SetGrid(1);
hpcIndexVE->Draw("colz"); hpcIndexVE->Draw("colz");
//=============================================== pad-4 //=============================================== pad-4
padID++; padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
canvas->cd(padID);
canvas->cd(padID)->SetGrid(1);
hsx3Coin->Draw("colz"); hsx3Coin->Draw("colz");
//=============================================== pad-5 //=============================================== pad-5
padID++; padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
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 //=============================================== pad-6
padID++; padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
canvas->cd(padID);
canvas->cd(padID)->SetGrid(1);
hpcCoin->Draw("colz"); hpcCoin->Draw("colz");
//=============================================== pad-7 //=============================================== pad-7
padID++; padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
canvas->cd(padID);
canvas->cd(padID)->SetGrid(1);
// hsx3VpcIndex ->Draw("colz"); // hsx3VpcIndex ->Draw("colz");
hsx3VpcE->Draw("colz") ; hsx3VpcE->Draw("colz") ;
//=============================================== pad-8 //=============================================== pad-8
padID++; padID ++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
canvas->cd(padID);
canvas->cd(padID)->SetGrid(1);
// hqqqVpcIndex ->Draw("colz"); // hqqqVpcIndex ->Draw("colz");
@ -486,8 +400,8 @@ void Analyzer::Terminate()
// canvas->cd(padID)->DrawFrame(-50, -50, 50, 50); // canvas->cd(padID)->DrawFrame(-50, -50, 50, 50);
// hqqqPolar->Draw("same colz pol"); // hqqqPolar->Draw("same colz pol");
canvas->cd(padID); canvas->cd(padID); canvas->cd(padID)->SetGrid(1);
canvas->cd(padID)->SetGrid(1);
// hZProj->Draw(); // hZProj->Draw();
hanVScatsum->Draw("colz"); hanVScatsum->Draw("colz");
} }

121
Armory/ClassPW.h
View File

@ -6,16 +6,14 @@
#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;
@ -27,21 +25,8 @@ struct PWHitInfo
} }
}; };
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
}
};
//!######################################################## //!########################################################
class PW class PW{ // proportional wire
{ // proportional wire
public: public:
PW(){ ClearHitInfo();}; PW(){ ClearHitInfo();};
~PW(){}; ~PW(){};
@ -52,9 +37,6 @@ public:
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;}
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 GetTrackPos() const {return trackPos;}
TVector3 GetTrackVec() const {return trackVec;} TVector3 GetTrackVec() const {return trackVec;}
double GetTrackTheta() const {return trackVec.Theta();} double GetTrackTheta() const {return trackVec.Theta();}
@ -84,8 +66,7 @@ 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,
@ -98,6 +79,7 @@ public:
} }
private: private:
PWHitInfo hitInfo; PWHitInfo hitInfo;
TVector3 trackPos; TVector3 trackPos;
@ -113,25 +95,23 @@ 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();
@ -143,8 +123,7 @@ inline void PW::ConstructGeo()
//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()),
@ -169,29 +148,23 @@ inline void PW::ConstructGeo()
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) if( phi > 0 && phiS > phiL ) phiL = phiL + TMath::TwoPi();
phiL = phiL + TMath::TwoPi(); if( phi < 0 && phiS > phiL ) phiS = phiS - 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); disA = Distance( pos, pos + direction, An[i].first, An[i].second);
if (disA < hitInfo.nearestDist.first) if( disA < hitInfo.nearestDist.first ){
{
hitInfo.nearestDist.first = disA; hitInfo.nearestDist.first = disA;
hitInfo.nearestWire.first = i; hitInfo.nearestWire.first = i;
} }
@ -201,70 +174,53 @@ inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose)
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) if( phi > 0 && phiS > phiL ) phiL = phiL + TMath::TwoPi();
phiL = phiL + TMath::TwoPi(); if( phi < 0 && phiS > phiL ) phiS = phiS - 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); disC = Distance( pos, pos + direction, Ca[i].first, Ca[i].second);
if (disC < hitInfo.nearestDist.second) if( disC < hitInfo.nearestDist.second ){
{
hitInfo.nearestDist.second = disC; hitInfo.nearestDist.second = disC;
hitInfo.nearestWire.second = i; hitInfo.nearestWire.second = i;
} }
} }
if (verbose) if(verbose) printf(" %2d | %8.2f, %8.2f\n", i, disA, disC);
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; short aaa1 = anode1 - 1; if( aaa1 < 0 ) aaa1 += nWire;
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; short ccc1 = cathode1 - 1; if( ccc1 < 0 ) ccc1 += nWire;
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) if( verbose ) Print();
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;
@ -274,12 +230,11 @@ 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) if( verbose ) printf("Theta, Phi = %f, %f \n", trackVec.Theta() *TMath::RadToDeg(), trackVec.Phi()*TMath::RadToDeg());
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;
@ -311,12 +266,11 @@ 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) if( verbose ) printf("Theta, Phi = %f, %f \n", trackVec.Theta() *TMath::RadToDeg(), trackVec.Phi()*TMath::RadToDeg());
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();
@ -324,6 +278,7 @@ inline double PW::GetZ0()
double theta = trackVec.Theta(); double theta = trackVec.Theta();
return trackPos.Z() - rho / TMath::Tan(theta); return trackPos.Z() - rho / TMath::Tan(theta);
} }
#endif #endif

12
ProcessRun.sh
View File

@ -10,20 +10,18 @@ runID=$1
timeWindow=$2 timeWindow=$2
rawFolder=/home/tandem/Desktop/analysis/data rawFolder=/home/tandem/Desktop/analysis/data
rootFolder=../root_data rootFolder=/home/tandem/Desktop/analysis/data/root_data
if [ $option -eq 0 ]; then rsync -a splitpole@128.186.111.223:/media/nvmeData/ANASEN27Alap/*.fsu /home/tandem/Desktop/analysis/data
rsync -auh --info=progress2 splitpole@128.186.111.223:/media/nvmeData/2024_09_17Fap/*.fsu /home/tandem/data1/2024_09_17Fap/data fileList=`\ls -1 ${rawFolder}/Run_${runID}_*.fsu`
fileList=`\ls -1 ${rawFolder}/*Run_${runID}_*.fsu` ./EventBuilder ${timeWindow} 0 0 10000000 ${fileList}
./EventBuilder ${timeWindow} 0 0 100000000 ${fileList}
outFile=${rawFolder}/*${runID}*${timeWindow}.root outFile=${rawFolder}/*${runID}*${timeWindow}.root
mv -vf ${outFile} ${rootFolder}/. mv -vf ${outFile} ${rootFolder}/.
./Mapper ${rootFolder}/*${runID}*${timeWindow}.root ./Mapper ${rootFolder}/*${runID}*${timeWindow}.root
fi
root "processRun.C(\"${rootFolder}/Run_${runID}_mapped.root\")" root "processRun.C(\"${rootFolder}/Run_${runID}_mapped.root\")"