diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json index e523710..341687e 100644 --- a/.vscode/c_cpp_properties.json +++ b/.vscode/c_cpp_properties.json @@ -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", diff --git a/.vscode/settings.json b/.vscode/settings.json index 59115c2..fbd381c 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -1,7 +1,5 @@ { "files.associations": { - "NiceMatStyle.C": "cpp", - "*.rmd": "markdown", "ryanScript.C": "cpp", "ryanSelector.C": "cpp", "array": "cpp", @@ -93,8 +91,16 @@ "TrackRecon.C": "cpp", "processRuns.C": "cpp", "Analysis.C": "cpp", - "Analyzer1.C": "cpp", - "PCGainMatch.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" }, - "C_Cpp.clang_format_fallbackStyle": "{BasedonStyle: Google, IndentWidth: 2, ColumnLimit: 0}" + "github-enterprise.uri": "https://fsunuc.physics.fsu.edu" } \ No newline at end of file diff --git a/Analyzer.C b/Analyzer.C index d500897..3047704 100644 --- a/Analyzer.C +++ b/Analyzer.C @@ -1,16 +1,16 @@ #define Analyzer_cxx #include "Analyzer.h" -#include -#include #include -#include #include -#include -#include +#include +#include + +#include +#include -#include "Armory/ClassPW.h" #include "Armory/ClassSX3.h" +#include "Armory/ClassPW.h" #include "TVector3.h" @@ -28,43 +28,25 @@ TH2F *hqqqVpcIndex; TH2F *hqqqVpcE; TH2F *hsx3VpcE; 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; SX3 sx3_contr; PW pw_contr; +PW pwinstance; TVector3 hitPos; bool HitNonZero; TH1F *hZProj; -TCutG *PCCoinc; -TCutG *alpha_cut_up; -TCutG *alpha_cut_down; -TCutG *cutg; -bool inCut; -bool inCutUp; -bool inCutDown; -bool inCutG; - -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, 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"); 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(-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(-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(-12, "y"); - hZProj = new TH1F("hZProj", "ZProjection", 600, -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); + 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); sx3_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; @@ -141,12 +103,6 @@ Bool_t Analyzer::Process(Long64_t entry) { b_pcCh->GetEntry(entry); b_pcE->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(); qqq.CalIndex(); @@ -158,41 +114,50 @@ Bool_t Analyzer::Process(Long64_t entry) { // //======================= SX3 std::vector> ID; // first = id, 2nd = index - for (int i = 0; i < sx3.multi; i++) { - if (sx3.e[i] > 50) { - ID.push_back(std::pair(sx3.id[i], i)); + for (int i = 0; i < sx3.multi; i++) + { + ID.push_back(std::pair(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]); - } + 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] ); + 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 &a, const std::pair &b) { return a.first < b.first; }); + if (ID.size() > 0) + { + std::sort(ID.begin(), ID.end(), [](const std::pair &a, const std::pair &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> 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]); } @@ -201,30 +166,39 @@ 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++) { + for (int j = 0; j < pc.multi; 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]); - // 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] ); } } @@ -238,276 +212,196 @@ Bool_t Analyzer::Process(Long64_t entry) { } // //======================= QQQ - for (int i = 0; i < qqq.multi; i++) { - if (qqq.e[i] > 50) { + 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 = 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++) { - if (pc.index[j] < 24 && pc.e[j] > 50) { - hqqqVpcE->Fill(qqq.e[i], pc.e[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; - 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]) { // 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]) { - chRing = qqq.ch[j] - 16; - chWedge = qqq.ch[i]; - } else { - chRing = qqq.ch[i]; - chWedge = qqq.ch[j] - 16; - } + 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); + // 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; + 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; - } + if (!HitNonZero) + { + double x = rho * TMath::Cos(theta); + double y = rho * TMath::Sin(theta); + hitPos.SetXYZ(x, y, 23 + 75 + 30); + HitNonZero = true; } } } } // //======================= PC - std::vector> anodeHits; - std::vector> cathodeHits; - int aID = 0; - int cID = 0; - float cEMax = 0; - int cIDMax = 0; - float cEnextMax = 0; - int cIDnextMax = 0; - float aE = 0; - float cE = 0; + ID.clear(); + int counter = 0; + std::vector> E; + E.clear(); + for (int i = 0; i < pc.multi; i++) + { - // Define the excluded SX3 and QQQ channels - std::unordered_set 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 excludeQQQ = {0, 17, 109, 110, 111, 112, 113, 119, 127, 128}; + if (pc.e[i] > 100) + ID.push_back(std::pair(pc.id[i], i)); + if (pc.e[i] > 100) + E.push_back(std::pair(pc.index[i], pc.e[i])); - for (int i = 0; i < pc.multi; i++) { - // for(int j=0; jFill(pc.index[i], pc.e[i]); - if (pc.e[i] > 100 & pc.multi < 7) { - // 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 (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 ); - // for (int j=0;jFill( pc.index[i], pc.e[i] ); - - for (int j = i + 1; j < pc.multi; j++) { - inCut = false; - if (PCCoinc->IsInside(pc.index[i], pc.index[j])) { - inCut = true; + pwinstance.ConstructGeo(); + 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); + } } - // hpcCoin->Fill( pc.index[i], pc.index[j]); } - // if(pc.e[i]>100){ - if (pc.index[i] < 24) { - anodeHits.push_back(std::pair(pc.index[i], pc.e[i])); - // anodeCount++; - } else if (pc.index[i] >= 24) { - cathodeHits.push_back(std::pair(pc.index[i], pc.e[i])); - } - // } - // } - // } - // hpcIndexVE->Fill( pc.index[i], pc.e[i] ); } } - hPCMultiplicity->Fill(pc.multi); - float aESum = 0; - float cESum = 0; - if (anodeHits.size() == 1 && cathodeHits.size() >= 1) { + if (E.size() >= 3) + { - inCutDown = false; - inCutUp = false; - for (const auto &anode : anodeHits) { + int aID = 0; + int cID = 0; - // for(int l=0; l 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]); - } - } + float aE = 0; + float cE = 0; + // if( ID[0].first < 1 ) { + // aID = pc.ch[ID[0].second]; + // cID = pc.ch[ID[1].second]; + // }else{ + // cID = pc.ch[ID[0].second]; + // aID = pc.ch[ID[1].second]; // } - // if (inCut) { - hanVScatsum->Fill(aE, cESum); - hAnodeHits->Fill(aID, aE); - hAnodeHits1d->Fill(anodeHits.size()); - // } - // } - } + // printf("anode= %d, cathode = %d\n", aID, cID); - // Miscellaneous channels including the Lollipop IC and Si detectors and hot needle IC - // Misc ch 0,1, 2, 3, 4 in order are the LIC, LSi, HNIC-difference, MCP, and Rf - bool timing = false; - 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]); + for (int k = 0; k < qqq.multi; k++) + { + if (qqq.index[k] == 75 && pc.index[k] == 2 && pc.e[k] > 100) + { - if (misc.ch[i] == 1) { - // hSi->Fill(misc.e[i]); - SiE = misc.e[i]; - SiT = misc.t[i] + misc.tf[i] * 4. / 1000; - // hSi->Fill(misc.e[i]); - } - if (misc.ch[i] == 2) { - ICt = misc.t[i] + misc.tf[i] * 4. / 1000; - ICe = misc.e[i]; - hSi->Fill(misc.e[i]); - } - if (misc.ch[i] == 3) { - // only analyze the first MCP in any event - if (MCPt == 0) { - MCPt = misc.t[i] + misc.tf[i] * 4. / 1000; - MCPE = misc.e[i]; - } - } - if (misc.ch[i] == 4) { - // only analyze the first RF in any event - 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; + int multi_an = 0; + for (int l = 0; l < E.size(); l++) + { + if (E[l].first < 24) + { + multi_an++; + } + } + + if (multi_an >= 1) + { + for (int l = 0; l < E.size(); l++) + { + if (E[l].first < 24 && E[l].first != 19 && E[l].first != 12) + { + aE = E[l].second; + } + else if (E[l].first > 24) + { + cE = E[l].second; + } + } } - // 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 (timing == true) { - // hICvsSi->Fill(misc.e[i], misc.e[j]); - if (misc.ch[i] == 1) { - hSi_gated->Fill(misc.e[i]); - // } - } - // } + 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.CalTrack1( hitPos, aID, cIDMax, cIDnextMax, cEMax, cEnextMax,1); - - // pw_contr.CalTrack(hitPos, aID, cID); - hZProj->Fill(pw_contr.GetZ0()); + if (HitNonZero) + { + pw_contr.CalTrack(hitPos, aID, cID); + hZProj->Fill(pw_contr.GetZ0()); + } } // ########################################################### Track constrcution @@ -517,18 +411,16 @@ Bool_t Analyzer::Process(Long64_t entry) { return kTRUE; } -void Analyzer::Terminate() { +void Analyzer::Terminate() +{ gStyle->SetOptStat("neiou"); TCanvas *canvas = new TCanvas("cANASEN", "ANASEN", 2000, 2000); - // TCanvas *a = new TCanvas("aANASEN", "ANASEN", 800, 600); 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++; canvas->cd(padID); canvas->cd(padID)->SetGrid(1); @@ -598,6 +490,4 @@ void Analyzer::Terminate() { canvas->cd(padID)->SetGrid(1); // hZProj->Draw(); hanVScatsum->Draw("colz"); - // hAnodeHits->Draw("colz"); - // // hAnodeMultiplicity->Draw(); } diff --git a/Armory/#ClassPW.h# b/Armory/#ClassPW.h# new file mode 100644 index 0000000..ddfa498 --- /dev/null +++ b/Armory/#ClassPW.h# @@ -0,0 +1,283 @@ +#ifndef ClassPW_h +#define ClassPW_h + +#include +#include +#include + +struct PWHitInfo{ + std::pair nearestWire; // anode, cathode + std::pair nearestDist; // anode, cathode + + std::pair nextNearestWire; // anode, cathode + std::pair 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 GetNearestID() const {return hitInfo.nearestWire;} + std::pair GetNearestDistance() const {return hitInfo.nearestDist;} + std::pair Get2ndNearestID() const {return hitInfo.nextNearestWire;} + std::pair 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> An; // the anode wire position vector in space + std::vector> 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 p1; // anode + std::pair 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 \ No newline at end of file diff --git a/Armory/ClassDet.h b/Armory/ClassDet.h index e6d4b9c..510cbb4 100644 --- a/Armory/ClassDet.h +++ b/Armory/ClassDet.h @@ -14,7 +14,7 @@ public: unsigned short ch[MAXMULTI]; unsigned short e[MAXMULTI]; unsigned long long t[MAXMULTI]; - unsigned long long tf[MAXMULTI]; + unsigned short sn[MAXMULTI]; unsigned short digiCh[MAXMULTI]; @@ -28,7 +28,6 @@ public: ch[i] = 0; e[i] = 0; t[i] = 0; - tf[i] = 0; index[i] = 0; sn[i] = 0; digiCh[i] = 0; diff --git a/Armory/ClassPW.h b/Armory/ClassPW.h index de4fc58..d65e559 100644 --- a/Armory/ClassPW.h +++ b/Armory/ClassPW.h @@ -6,14 +6,16 @@ #include #include -struct PWHitInfo{ - std::pair nearestWire; // anode, cathode +struct PWHitInfo +{ + std::pair nearestWire; // anode, cathode std::pair nearestDist; // anode, cathode - std::pair nextNearestWire; // anode, cathode + std::pair nextNearestWire; // anode, cathode std::pair nextNearestDist; // anode, cathode - void Clear(){ + void Clear() + { nearestWire.first = -1; nearestWire.second = -1; nearestDist.first = 999999999; @@ -25,40 +27,56 @@ 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 GetNearestID() const {return hitInfo.nearestWire;} - std::pair GetNearestDistance() const {return hitInfo.nearestDist;} - std::pair Get2ndNearestID() const {return hitInfo.nextNearestWire;} - std::pair 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 GetNearestID() const { return hitInfo.nearestWire; } + std::pair GetNearestDistance() const { return hitInfo.nearestDist; } + std::pair Get2ndNearestID() const { return hitInfo.nextNearestWire; } + std::pair Get2ndNearestDistance() const { return hitInfo.nextNearestDist; } + + std::vector> An; // the anode wire position vector in space + std::vector> 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(); @@ -66,20 +84,20 @@ public: 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); + 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; @@ -87,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; @@ -95,23 +113,25 @@ private: double anodeLength; double cathodeLength; - std::vector> An; // the anode wire position vector in space - std::vector> Ca; // the cathode wire position vector in space + // std::vector> An; // the anode wire position vector in space + // std::vector> 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(); @@ -119,108 +139,132 @@ inline void PW::ConstructGeo(){ std::pair p1; // anode std::pair 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); + 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) ); + 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; @@ -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 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; @@ -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 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 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 trackPos.Z() - rho / TMath::Tan(theta); } -#endif \ No newline at end of file +#endif \ No newline at end of file diff --git a/Armory/Mapper.cpp b/Armory/Mapper.cpp index 0ebb949..8030de1 100644 --- a/Armory/Mapper.cpp +++ b/Armory/Mapper.cpp @@ -67,7 +67,6 @@ int main(int argc, char **argv){ Det sx3; Det qqq; Det pc ; - Det misc; printf(" Raw root file : %s\n", inFileName.c_str()); 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("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 //clock @@ -121,12 +112,8 @@ int main(int argc, char **argv){ sx3.multi = 0; qqq.multi = 0; pc.multi = 0; - misc.multi=0; - sx3.Clear(); qqq.Clear(); - pc.Clear(); - misc.Clear(); 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.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