added uncertainties in calculating th etrack

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();

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 p1 = TMath::Abs(hitInfo.nearestDist.first + gRandom->Gaus(0, sigmaA));
-  double totDistA = hitInfo.nearestDist.first + hitInfo.nextNearestDist.first;
+  double p2 = TMath::Abs(hitInfo.nextNearestDist.first + gRandom->Gaus(0, sigmaA));
-  double fracA = hitInfo.nearestDist.first / totDistA;
+  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 q1 = TMath::Abs(hitInfo.nearestDist.second + gRandom->Gaus(0, sigmaC));
-  double fracC = hitInfo.nearestDist.second / totDistC;
+  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;

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){
 

Armory/anasenMS.cpp

@@ -35,9 +35,14 @@ int main(int argc, char **argv){
   std::vector<float> ExAList = {0};
   std::vector<float> ExList = {0, 1, 2};
 
-  double vertexXRange[2] = { 0, 0}; 
+  double vertexXRange[2] = { -10, 10}; // mm
-  double vertexYRange[2] = { 0, 0}; 
+  double vertexYRange[2] = { -10, 10}; 
-  double vertexZRange[2] = { 0, 0}; 
+  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();
+      //Introduce uncertaity
-      sx3Y = sx3->GetHitPos().Y();
+      // TVector3 hitPos = sx3->GetHitPos();
-      sx3Z = sx3->GetHitPos().Z();
+      TVector3 hitPos = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L);
 
-      pw->CalTrack(sx3->GetHitPos(), anodeID[0], cathodeID[0], false);
+      sx3X = hitPos.X();
+      sx3Y = hitPos.Y();
+      sx3Z = hitPos.Z();
+      
+      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();
 

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;
+  pw->FindWireID(pos, dir, true);
-  pc.ConstructGeo();
+  sx3->FindSX3Pos(pos, dir, true);
-  pc.FindWireID(pos, dir, true);
 
-  SX3 sx3;
-  sx3.ConstructGeo();
-  sx3.FindSX3Pos(pos, dir, true);
 
   // haha->CalTrack(sx3.hitPos, wireID.first, wireID.second, true);