diff --git a/Armory/ClassPC1An.h b/Armory/ClassPC1An.h deleted file mode 100644 index 8e06560..0000000 --- a/Armory/ClassPC1An.h +++ /dev/null @@ -1,246 +0,0 @@ -#ifndef ClassPC_h -#define ClassPC_h - -#include -#include -#include -#include - -struct PCHit_1An{ - std::pair nearestWire; // anode, cathode - std::pair nearestDist; // anode, cathode - - short nextNearestWire; // cathode - double nextNearestDist; // cathode - - void Clear(){ - nearestWire.first = -1; - nearestWire.second = -1; - nearestDist.first = 999999999; - nearestDist.second = 999999999; - nextNearestWire= -1; - nextNearestDist = 999999999; - } -}; - -//!######################################################## -class PC{ // proportional wire -public: - PC(){ ClearHitInfo();}; - ~PC(){}; - - PCHit_1An GetHitInfo() const {return hitInfo;} - std::pair GetNearestID() const {return hitInfo.nearestWire;} - std::pair GetNearestDistance() const {return hitInfo.nearestDist;} - short Get2ndNearestID() const {return hitInfo.nextNearestWire;} - 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 CalTrack3(TVector3 sx3Pos, PCHit_1An 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 Cathode: %2d(%5.2f)\n", hitInfo.nextNearestWire, - hitInfo.nextNearestDist); - } - -private: - - // PCHitInfo hitInfo; - PCHit_1An 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 PC::ClearHitInfo(){ - hitInfo.Clear(); -} - -inline void PC::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 PC::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); - } - - short cathode1 = hitInfo.nearestWire.second; - short ccc1 = cathode1 - 1; if( ccc1 < 0 ) ccc1 += nWire; - short ccc2 = (cathode1 + 1) % nWire; - - double haha1 = Distance( pos, pos + direction, Ca[ccc1].first, Ca[ccc1].second); - double haha2 = Distance( pos, pos + direction, Ca[ccc2].first, Ca[ccc2].second); - if( haha1 < haha2){ - hitInfo.nextNearestWire = ccc1; - hitInfo.nextNearestDist = haha1; - }else{ - hitInfo.nextNearestWire = ccc2; - hitInfo.nextNearestDist= haha2; - } - - if( verbose ) Print(); -} - - -inline void PC::CalTrack3(TVector3 sx3Pos, PCHit_1An hitInfo, double sigmaA, double sigmaC, bool verbose){ - - trackPos = sx3Pos; - - double p1 = TMath::Abs(hitInfo.nearestDist.first + gRandom->Gaus(0, sigmaA)); - short anodeID1 = hitInfo.nearestWire.first; - - double q1 = TMath::Abs(hitInfo.nearestDist.second + gRandom->Gaus(0, sigmaC)); - double q2 = TMath::Abs(hitInfo.nextNearestDist+ gRandom->Gaus(0, sigmaC)); - double fracC = q1 / (q1 + q2); - short cathodeID1 = hitInfo.nearestWire.second; - short cathodeID2 = hitInfo.nextNearestWire; - TVector3 shiftC1 = (Ca[cathodeID2].first - Ca[cathodeID1].first) * fracC; - TVector3 shiftC2 = (Ca[cathodeID2].second - Ca[cathodeID1].second) * fracC; - - TVector3 a1 = An[anodeID1].first; - - TVector3 c1 = Ca[cathodeID1].first + shiftC1; - TVector3 c2 = Ca[cathodeID1].second + shiftC2; - - TVector3 n1 = (sx3Pos - a1).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 PC::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/ClassPW.h b/Armory/ClassPW.h index 1db9d69..2ed454b 100755 --- a/Armory/ClassPW.h +++ b/Armory/ClassPW.h @@ -154,8 +154,8 @@ inline void PW::ConstructGeo() std::pair q1; // cathode double k = TMath::TwoPi() / 24.; // 48 solder thru holes, wires in every other one - double offset_a1 = -6 * k - 3 * k; - double offset_c1 = -4 * k - 2 * k - TMath::TwoPi() / 48; // correct for a half-turn + double offset_a1 = -5 * k - 5 * k; + double offset_c1 = -5 * k - 2 * k - TMath::TwoPi() / 48; // correct for a half-turn // std::cerr << "Here!" << std::endl; // #include "../scratch/testing.h" double offset_a2 = offset_a1 + wireShift * k; diff --git a/Armory/PC_StepLadder_Correction.h b/Armory/PC_StepLadder_Correction.h index 6618596..439fac3 100644 --- a/Armory/PC_StepLadder_Correction.h +++ b/Armory/PC_StepLadder_Correction.h @@ -32,7 +32,7 @@ double model_invert(double* y, double* p) { break; } } - return result+80; + return result; } double model_a1c1(double* x, double* p) { diff --git a/anasen_fem/paraview_plotter.py b/anasen_fem/paraview_plotter.py index 5270785..e2762b4 100755 --- a/anasen_fem/paraview_plotter.py +++ b/anasen_fem/paraview_plotter.py @@ -1,5 +1,5 @@ -#!/home/vsitaraman/ParaView-6.1.0-RC1-MPI-Linux-Python3.12-x86_64/bin/pvbatch -######## !/home/vs19g/ParaView-6.1.0-MPI-Linux-Python3.12-x86_64/bin/pvbatch +#!/home/vsitaraman/ParaView-6.1.0-MPI-Linux-Python3.12-x86_64/bin/pvbatch +########### !/home/vsitaraman/ParaView-6.1.0-RC1-MPI-Linux-Python3.12-x86_64/bin/pvbatch import numpy as np import sys from paraview.simple import * diff --git a/run_17F.sh b/run_17F.sh index 906fb88..bd672dc 100644 --- a/run_17F.sh +++ b/run_17F.sh @@ -2,7 +2,8 @@ rm results_run*.root export DATASET="17F" export flip180="0" export flipa=0 -export anode_offset=0 +export anode_offset=2 +export cathode_offset=0 #root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Source_005_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run05.root; #root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Source_006_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run06.root; #root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Source_007_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run07.root; @@ -26,18 +27,19 @@ export anode_offset=0 #root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_021_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run21.root; #17F reaction data -export flip180="0" -declare -i run=231 #49 -while [[ $run -lt 258 ]]; do #392 +# export flip180="0" +declare -i run=322 #49 +while [[ $run -lt 399 ]]; do #392 wrun=$(printf "%03d" $run) file_exists=$(test -f ../ANASEN_analysis/data/17F_Data/Run_"$wrun"_mapped.root) if [[ $file_exists -ne 0 ]]; then continue; fi - root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Run_"$wrun"_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run$wrun.root; + root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Run_"$wrun"_mapped.root -e 'tree->Process("MakeVertex.C+O","Analyzer_17F.root")'; mv Analyzer_17F.root 17F_output/results_run$wrun.root; run=run+1 done -rm output_17F.root -hadd -j 4 -k output_17F.root results_run2*.root +rm output.root +hadd -k -j 4 output.root 17F_output/results_run*.root +mv output.root output_17F.root unset souce_vertex unset DATASET diff --git a/run_27Al.sh b/run_27Al.sh index b0e3011..c86e467 100644 --- a/run_27Al.sh +++ b/run_27Al.sh @@ -2,7 +2,9 @@ rm results_run*.root export DATASET="27Al" export flip180="0" export flipa=0 +export flipc=0 export anode_offset=0 +export cathode_offset=0 #declare -i run=28 #while [[ $run -lt 34 ]]; do #runs 1 to 84 # wrun=$(printf "%03d" $run) diff --git a/run_sx3.sh b/run_sx3.sh index 51409d7..fe6e156 100755 --- a/run_sx3.sh +++ b/run_sx3.sh @@ -2,6 +2,8 @@ # rm results_run*.root export flipa=0 export anode_offset=0 +export cathode_offset=0 + export DATASET="27Al" if [[ 1 -eq 0 ]]; then #root -b -q -l -x ../ANASEN_analysis/data/27Al_Data/Run_009_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run09.root; @@ -18,9 +20,8 @@ fi #exit #alpha+gas 27Al export DATASET="27Al" -export flip180="0" #root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_009_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run09.root; -if [[ 1 -eq 0 ]]; then +if [[ 1 -eq 1 ]]; then #export timecut_low=230.0; export timecut_low=400.0; #export timecut_high=400.0; @@ -29,7 +30,7 @@ export timecut_low=400.0; #export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_010_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run10.root; #export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_011_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run11.root; export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_012_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run12.root; -exit +# exit #export source_vertex=53.44; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_013_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run13.root; #exit fi @@ -37,12 +38,10 @@ fi #protons+gas, 27Al #export flip180="1" #export flip180="0" -if [[ 1 -eq 1 ]]; then -export flipa=0 -export anode_offset=0 +if [[ 1 -eq 0 ]]; then export source_vertex=-200.0; #put the 'source' on the entrance window root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_018_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run18.root; -exit +# exit root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_015_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run15.root; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_017_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run17.root; root -q -b -x ../ANASEN_analysis/data/27Al_Data/Run_019_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run19.root; @@ -60,7 +59,6 @@ fi #root -l -x results_run19.root results_run12.root -e "new TBrowser" #exit export DATASET="17F" -export flip180="0" if [[ 1 -eq 0 ]]; then root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Source_005_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run05.root; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/Source_006_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run06.root; @@ -86,8 +84,7 @@ export source_vertex=-24.96; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/S export source_vertex=-73.96; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/SourceRun_021_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run21.root; fi #17F reaction data -#export flip180="0" -if [[ 1 -eq 0 ]]; then +if [[ 1 -eq 1 ]]; then export source_vertex=-57.28; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/ProtonRun_035_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run35.root; #export source_vertex=-8.28; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/ProtonRun_036_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root resulrs_run36.root; #export source_vertex=-27.88; root -q -l -b -x ../ANASEN_analysis/data/17F_Data/ProtonRun_037_mapped.root -e 'tree->Process("MakeVertex.C+O")'; mv Analyzer_SX3.root results_run37.root;