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

bug fix wireID, FindWireID

Browse Source
This commit is contained in:
Ryan Tang 2024-02-02 16:54:58 -05:00
parent 6790ef41df
commit d710df842d
5 changed files with 50 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.gitignore vendored
View File

@ -5,4 +5,4 @@ EventBuilder*
*.root
Mapper
anasenMS
AnasenMS

26
Armory/ClassPW.h
View File

@ -188,26 +188,30 @@ inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose ){
//==== 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[anode1-1].first, An[anode1-1].second);
double haha2 = Distance( pos, pos + direction, An[anode1+1].first, An[anode1+1].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);
if( haha1 < haha2){
hitInfo.nextNearestWire.first = anode1-1;
hitInfo.nextNearestWire.first = aaa1;
hitInfo.nextNearestDist.first = haha1;
}else{
hitInfo.nextNearestWire.first = anode1+1;
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[cathode1-1].first, Ca[cathode1-1].second);
haha2 = Distance( pos, pos + direction, Ca[cathode1+1].first, Ca[cathode1+1].second);
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 = cathode1-1;
hitInfo.nextNearestWire.second = ccc1;
hitInfo.nextNearestDist.second = haha1;
}else{
hitInfo.nextNearestWire.second = cathode1+1;
hitInfo.nextNearestWire.second = ccc2;
hitInfo.nextNearestDist.second = haha2;
}
@ -241,11 +245,11 @@ inline void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, bool verbose){
TVector3 shiftA2 = (An[anodeID2].second - An[anodeID1].second) * fracA;
double totDistC = hitInfo.nearestDist.second + hitInfo.nextNearestDist.second;
double fracC = hitInfo.nearestDist.second / totDistA;
double fracC = hitInfo.nearestDist.second / totDistC;
short cathodeID1 = hitInfo.nearestWire.second;
short cathodeID2 = hitInfo.nextNearestWire.second;
TVector3 shiftC1 = (Ca[anodeID2].first - Ca[anodeID1].first) * fracC;
TVector3 shiftC2 = (Ca[anodeID2].second - Ca[anodeID1].second) * fracC;
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;

4
Armory/Makefile
View File

@ -23,6 +23,6 @@ Mapper : Mapper.cpp ../mapping.h ClassDet.h
@echo "--------- making Mapper"
$(CC) $(COPTS) -o Mapper Mapper.cpp $(ROOTLIBS)
AnasenMS : ClassAnasen.h anasenMS.cpp Isotope.h constant.h
AnasenMS : constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassPW.h ClassAnasen.h anasenMS.cpp
@echo "--------- making ANASEN Monte Carlo"
$(CC) $(COPTS) -o anasenMS anasenMS.cpp $(ROOTLIBS)
$(CC) $(COPTS) -o AnasenMS anasenMS.cpp $(ROOTLIBS)

40
Armory/anasenMS.cpp
View File

@ -90,9 +90,13 @@ int main(int argc, char **argv){
tree->Branch("sx3Y", &sx3Y, "sx3Y/D");
tree->Branch("sx3Z", &sx3Z, "sx3Z/D");
int anodeID, cathodeID;
tree->Branch("aID", &anodeID, "anodeID/I");
tree->Branch("cID", &cathodeID, "cathodeID/I");
int anodeID[2], cathodeID[2];
tree->Branch("aID", anodeID, "anodeID/I");
tree->Branch("cID", cathodeID, "cathodeID/I");
double anodeDist[2], cathodeDist[2];
tree->Branch("aDist", anodeDist, "anodeDist/D");
tree->Branch("cDist", cathodeDist, "cathodeDist/D");
int sx3ID, sx3Up, sx3Dn, sx3Bk;
double sx3ZFrac;
@ -106,6 +110,10 @@ int main(int argc, char **argv){
tree->Branch("reTheta", &reTheta, "reconstucted_theta/D");
tree->Branch("rePhi", &rePhi, "reconstucted_phi/D");
double reTheta1, rePhi1;
tree->Branch("reTheta1", &reTheta1, "reconstucted_theta1/D");
tree->Branch("rePhi1", &rePhi1, "reconstucted_phi1/D");
//========timer
TBenchmark clock;
@ -156,10 +164,17 @@ int main(int argc, char **argv){
pw->FindWireID(vertex, dir, false);
sx3->FindSX3Pos(vertex, dir, false);
std::pair<int, int> wireID = pw->GetNearestID();
PWHitInfo hitInfo = pw->GetHitInfo();
anodeID = wireID.first;
cathodeID = wireID.second;
anodeID[0] = hitInfo.nearestWire.first;
cathodeID[0] = hitInfo.nearestWire.second;
anodeID[1] = hitInfo.nextNearestWire.first;
cathodeID[1] = hitInfo.nextNearestWire.second;
anodeDist[0] = hitInfo.nearestDist.first;
cathodeDist[0] = hitInfo.nearestDist.second;
anodeDist[1] = hitInfo.nextNearestDist.first;
cathodeDist[1] = hitInfo.nextNearestDist.second;
sx3ID = sx3->GetID();
if( sx3ID >= 0 ){
@ -172,11 +187,14 @@ int main(int argc, char **argv){
sx3Y = sx3->GetHitPos().Y();
sx3Z = sx3->GetHitPos().Z();
pw->CalTrack(sx3->GetHitPos(), wireID.first, wireID.second, false);
pw->CalTrack(sx3->GetHitPos(), anodeID[0], cathodeID[0], false);
reTheta = pw->GetTrackTheta() * TMath::RadToDeg();
rePhi = pw->GetTrackPhi() * TMath::RadToDeg();
pw->CalTrack2(sx3->GetHitPos(), hitInfo, false);
reTheta1 = pw->GetTrackTheta() * TMath::RadToDeg();
rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg();
}else{
sx3Up = -1;
sx3Dn = -1;
@ -192,7 +210,11 @@ int main(int argc, char **argv){
// }
reTheta = TMath::QuietNaN();
rePhi = TMath::QuietNaN();
rePhi = TMath::QuietNaN();
reTheta1 = TMath::QuietNaN();
rePhi1 = TMath::QuietNaN();
}
tree->Fill();

3
Armory/constant.h
View File

@ -7,7 +7,6 @@
* email: goluckyryan@gmail.com
* ********************************************************************/
#ifndef constant
#define constant
#include <cmath>
@ -69,7 +68,7 @@ const double hbar = 197.326979;
//======================================================================
std::vector<std::string> SplitStr(std::string tempLine, std::string splitter, int shift = 0){
inline std::vector<std::string> SplitStr(std::string tempLine, std::string splitter, int shift = 0){
std::vector<std::string> output;