From a668f64ab7dcda333cb45a99a7d0fb4a9043c2cd Mon Sep 17 00:00:00 2001 From: Ryan Tang Date: Fri, 2 Feb 2024 19:13:28 -0500 Subject: [PATCH] added uncertainties in calculating th etrack --- Armory/ClassAnasen.h | 17 +++++++++++++---- Armory/ClassPW.h | 17 +++++++++-------- Armory/ClassSX3.h | 18 ++++++++++++++++++ Armory/anasenMS.cpp | 25 +++++++++++++++++-------- script.C | 17 ++++++++++------- 5 files changed, 67 insertions(+), 27 deletions(-) diff --git a/Armory/ClassAnasen.h b/Armory/ClassAnasen.h index a08b804..90585b7 100644 --- a/Armory/ClassAnasen.h +++ b/Armory/ClassAnasen.h @@ -21,6 +21,13 @@ public: ANASEN(); ~ANASEN(); + void SetUncertainties(double sx3W, double sx3L, double anode, double cathode){ + sigmaA = anode; + sigmaC = cathode; + sigmaW = sx3W; + sigmaL = sx3L; + } + void DrawTrack(TVector3 pos, TVector3 direction, bool drawEstimatedTrack = false); void DrawDeducedTrack(TVector3 sx3Pos, int anodeID, int cathodeID); void DrawAnasen(int anodeID1 = -1, @@ -40,6 +47,9 @@ private: PW * pw; SX3 * sx3; + double sigmaA, sigmaC; // pw + double sigmaW, sigmaL; // sx3 + const float qqqR1 = 50; const float qqqR2 = 100; const float qqqZPos = 23 + 75 + 30; @@ -113,11 +123,9 @@ inline void ANASEN::Construct3DModel(int anodeID1, int anodeID2, int cathodeID1, TGeoVolume *axisX = geom->MakeTube("axisX", Al, 0, 0.1, 5.); axisX->SetLineColor(1); worldBox->AddNode(axisX, 1, new TGeoCombiTrans(5, 0, 0., new TGeoRotation("rotA", 90., 90., 0.))); - TGeoVolume *axisY = geom->MakeTube("axisY", Al, 0, 0.1, 5.); axisY->SetLineColor(1); worldBox->AddNode(axisY, 1, new TGeoCombiTrans(0, 5, 0., new TGeoRotation("rotB", 0., 90., 0.))); - TGeoVolume *axisZ = geom->MakeTube("axisZ", Al, 0, 0.1, 5.); axisZ->SetLineColor(1); worldBox->AddNode(axisZ, 1, new TGeoTranslation(0, 0, 5)); @@ -238,7 +246,8 @@ inline void ANASEN::DrawTrack(TVector3 pos, TVector3 direction, bool drawEstima worldBox->AddNode(startPos, 3, new TGeoCombiTrans( pos.X(), pos.Y(), pos.Z(), new TGeoRotation("rotA", 0, 0, 0.))); if( sx3->GetID() >= 0 ){ - TVector3 hitPos = sx3->GetHitPos(); + //TVector3 hitPos = sx3->GetHitPos(); + TVector3 hitPos = sx3->GetHitPosWithSigma(sigmaW, sigmaL); TGeoVolume * hit = geom->MakeSphere("hitpos", 0, 0, 3); hit->SetLineColor(kRed); @@ -259,7 +268,7 @@ inline void ANASEN::DrawTrack(TVector3 pos, TVector3 direction, bool drawEstima {//===== complicated PWHitInfo hitInfo = pw->GetHitInfo(); - pw->CalTrack2(hitPos, hitInfo, true); + pw->CalTrack2(hitPos, hitInfo, sigmaA, sigmaC, true); double thetaDeduce = pw->GetTrackTheta() * TMath::RadToDeg(); double phiDeduce = pw->GetTrackPhi() * TMath::RadToDeg(); diff --git a/Armory/ClassPW.h b/Armory/ClassPW.h index 947f28b..d1f5ef4 100644 --- a/Armory/ClassPW.h +++ b/Armory/ClassPW.h @@ -27,7 +27,7 @@ struct PWHitInfo{ //!######################################################## class PW{ // proportional wire public: - PW(){ ClearHitInfo(); }; + PW(){ ClearHitInfo();}; ~PW(){}; PWHitInfo GetHitInfo() const {return hitInfo;} @@ -62,7 +62,7 @@ public: 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, 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, @@ -232,20 +232,21 @@ inline void PW::CalTrack(TVector3 sx3Pos, int anodeID, int cathodeID, bool verbo } -inline void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, bool verbose){ +inline void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA, double sigmaC, bool verbose){ trackPos = sx3Pos; - // fraction between the nearest wire and and the 2nd nearest wire by distance - double totDistA = hitInfo.nearestDist.first + hitInfo.nextNearestDist.first; - double fracA = hitInfo.nearestDist.first / totDistA; + 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 totDistC = hitInfo.nearestDist.second + hitInfo.nextNearestDist.second; - double fracC = hitInfo.nearestDist.second / totDistC; + 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; diff --git a/Armory/ClassSX3.h b/Armory/ClassSX3.h index 3f29e59..0a57557 100644 --- a/Armory/ClassSX3.h +++ b/Armory/ClassSX3.h @@ -16,6 +16,7 @@ public: short GetChBk() const {return chBk;} TVector3 GetHitPos() const {return hitPos;} + TVector3 GetHitPosWithSigma(double sigmaY_mm, double sigmaZ_mm); double GetZFrac() const {return zFrac;} // range from -0.5 to 0.5 @@ -187,6 +188,23 @@ inline void SX3::FindSX3Pos(TVector3 pos, TVector3 direction, bool verbose){ if( verbose) Print(); } +inline TVector3 SX3::GetHitPosWithSigma(double sigmaY_mm, double sigmaZ_mm){ + + double phi = SNorml[id].Phi(); + + TVector3 haha = hitPos; + haha.RotateZ(-phi); + double y = haha.Y() + gRandom->Gaus(0, sigmaY_mm); + double z = haha.Z() + gRandom->Gaus(0, sigmaZ_mm); + + haha.SetY(y); + haha.SetZ(z); + haha.RotateZ(phi); + + return haha; + +} + inline void SX3::CalSX3Pos(unsigned short ID, unsigned short chUp, unsigned short chDown, unsigned short chBack, float eUp, float eDown){ diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index 238968d..778cd49 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -35,9 +35,14 @@ int main(int argc, char **argv){ std::vector ExAList = {0}; std::vector ExList = {0, 1, 2}; - double vertexXRange[2] = { 0, 0}; - double vertexYRange[2] = { 0, 0}; - double vertexZRange[2] = { 0, 0}; + double vertexXRange[2] = { -10, 10}; // mm + double vertexYRange[2] = { -10, 10}; + double vertexZRange[2] = { -70, 70}; + + double sigmaSX3_W = 10; // mm + double sigmaSX3_L = 0; // mm + double sigmaPW_A = 0; // from 0 to 1. + double sigmaPW_C = 0; // from 0 to 1. //################################################### transfer.CalReactionConstant(); @@ -183,15 +188,19 @@ int main(int argc, char **argv){ sx3Bk = sx3->GetChBk(); sx3ZFrac = sx3->GetZFrac(); - sx3X = sx3->GetHitPos().X(); - sx3Y = sx3->GetHitPos().Y(); - sx3Z = sx3->GetHitPos().Z(); + //Introduce uncertaity + // TVector3 hitPos = sx3->GetHitPos(); + TVector3 hitPos = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L); + + sx3X = hitPos.X(); + sx3Y = hitPos.Y(); + sx3Z = hitPos.Z(); - pw->CalTrack(sx3->GetHitPos(), anodeID[0], cathodeID[0], false); + pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false); reTheta = pw->GetTrackTheta() * TMath::RadToDeg(); rePhi = pw->GetTrackPhi() * TMath::RadToDeg(); - pw->CalTrack2(sx3->GetHitPos(), hitInfo, false); + pw->CalTrack2(hitPos, hitInfo, sigmaPW_A, sigmaPW_C, false); reTheta1 = pw->GetTrackTheta() * TMath::RadToDeg(); rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg(); diff --git a/script.C b/script.C index f24508a..24da7f0 100644 --- a/script.C +++ b/script.C @@ -246,8 +246,15 @@ void script(TString fileName = "", int maxEvent = -1){ TCanvas * c2 = new TCanvas("c2", "ANASEN Simulation", 800, 800); c2->cd(); + + ANASEN * haha = new ANASEN(); + haha->SetUncertainties(10, 10, 1, 1); + + PW * pw = haha->GetPW(); + SX3 * sx3 = haha->GetSX3(); + // haha->DrawAnasen(0, 1, 0, 1, -1, false); gRandom->SetSeed(0); @@ -273,13 +280,9 @@ void script(TString fileName = "", int maxEvent = -1){ haha->DrawTrack(pos, dir, true); - PW pc; - pc.ConstructGeo(); - pc.FindWireID(pos, dir, true); - - SX3 sx3; - sx3.ConstructGeo(); - sx3.FindSX3Pos(pos, dir, true); + pw->FindWireID(pos, dir, true); + sx3->FindSX3Pos(pos, dir, true); + // haha->CalTrack(sx3.hitPos, wireID.first, wireID.second, true);