QQQ

Armory/ClassAnasen.h

@@ -15,6 +15,7 @@
 
 #include "ClassSX3.h"
 #include "ClassPW.h"
+#include "ClassQQQ.h"
 
 //to not include certain wires in the simulation, pass the anode and cathode IDs to the constructor, e.g. for anode wires 5-10, pass anodeID1 = 5, anodeID2 = 10, and for cathode wires 20-25, pass cathodeID1 = 20, cathodeID2 = 25. To include all wires, pass -1 for all IDs.
 
@@ -43,6 +44,7 @@ public:
 
   PW * GetPW() {return pw;}
   SX3 * GetSX3() {return sx3;}
+  QQQ * GetQQQ() {return qqq;}
   TGeoManager * GetGeoManager() {return geom;}
   TGeoVolume * GetWorldBox() {return worldBox;}
 
@@ -50,6 +52,7 @@ private:
 
   PW * pw;
   SX3 * sx3;
+  QQQ * qqq;
 
   double sigmaA, sigmaC; // pw
   double sigmaW, sigmaL; // sx3
@@ -77,7 +80,7 @@ inline ANASEN::ANASEN(){
 
   pw = new PW();
   sx3 = new SX3();
-
+  qqq = new QQQ();
   CalGeometry();
 
   geom = nullptr;
@@ -236,6 +239,7 @@ inline  void ANASEN::DrawTrack(TVector3 pos, TVector3 direction, bool drawEstima
 
   pw->FindWireID(pos, direction);
   sx3->FindSX3Pos(pos, direction);
+  qqq->FindQQQPos(pos, direction);
 
   std::pair<short, short> wireID = pw->GetNearestID();
   

Armory/ClassQQQ.h

@@ -0,0 +1,288 @@
+#ifndef ClassQQQ_h
+#define ClassQQQ_h
+
+#include <cstdio>
+#include <TMath.h>
+#include <TVector3.h>
+#include <TRandom.h>
+#include "TGeoManager.h"
+#include "TGeoVolume.h"
+#include "TGeoBBox.h"
+
+class QQQ{
+public:
+  QQQ(){Clear();};
+  ~QQQ(){}
+  
+  short GetID() const {return id;}
+  short GetChUp() const {return chUp;}
+  short GetChDn() const {return chDn;}
+  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
+
+  void Clear();
+  void ConstructGeo();
+  void FindQQQPos(TVector3 pos, TVector3 direction, bool verbose = false);
+  void CalQQQPos(unsigned short ID, unsigned short chUp, unsigned short chDown, unsigned short chBack, float eUp, float eDown);
+
+  double GetNumDet() const {return numDet;}
+
+  void Print(){ 
+    if( id == -1 ){
+      printf("Did not hit any QQQ.\n"); 
+    }else{
+      printf("ID: %d, U,D,B: %d %d %d| zFrac : %.2f\n", id, chUp, chDn, chBk, zFrac); 
+      printf("Hit Pos: %.2f, %.2f, %.2f\n", hitPos.X(), hitPos.Y(), hitPos.Z()); 
+    }
+  }
+
+  // void CalZFrac(){
+  //   zFrac = (eUp - eDn)/(eUp + eDn); 
+  // }
+
+private:
+
+  const int numDet = 4;
+  const float qqqR1 = 50;
+  const float qqqR2 = 100;
+  const float qqqZPos = 23 + 75 + 30;
+
+  short id; // -1 when no hit
+  short chUp;
+  short chDn;
+  short chBk;
+
+  double zFrac; // from +1 (downstream) to -1 (upstream)
+
+  double eUp;
+  double eDn;
+  double eBk;
+
+  TVector3 hitPos;
+
+  TGeoManager *geom;
+  TGeoVolume *worldBox;
+  TGeoMedium *Al;
+
+  // helper function to calculate intersection between line segments, return pair of (fraction along line1, fraction along line2) where the intersection occurs. If no intersection, return (0, -1).
+  std::pair<double, double> Intersect(TVector3 p1, TVector3 p2, TVector3 q1, TVector3 q2, bool verbose){ 
+
+    //see https://nukephysik101.wordpress.com/2023/12/30/intersect-between-2-line-segments/
+    //zero all z-component
+    TVector3 a0 = p1; a0.SetZ(0);
+    TVector3 a1 = p2; a1.SetZ(0);
+
+    TVector3 b0 = q1; b0.SetZ(0);
+    TVector3 b1 = q2; b1.SetZ(0);
+
+    double h = 0, k = 0; // placeholder values, implementation of intersection logic
+    if( verbose ) printf(" ----h, k : %f, %f\n", h, k);
+
+    return std::pair<double,double>(h,k);
+  }
+
+};
+
+inline void QQQ::Clear(){
+  id = -1;
+  chUp = -1;
+  chDn = -1;
+  chBk = -1;
+  zFrac = TMath::QuietNaN();
+
+  eUp = TMath::QuietNaN();
+  eDn = TMath::QuietNaN();
+  eBk = TMath::QuietNaN();
+}
+
+inline void QQQ::ConstructGeo(){
+  TGeoVolume *qqq = geom->MakeTubs("qqq", Al, qqqR1, qqqR2, 0.5, 5, 85);
+  qqq->SetLineColor(7);
+  for( int i = 0; i < 4; i++){
+    worldBox->AddNode(qqq, i+1, new TGeoCombiTrans( 0, 
+                                                    0, 
+                                                    qqqZPos, 
+                                                    new TGeoRotation("rot1", 360/4 * (i), 0., 0.)));
+  }
+}
+
+
+inline void QQQ::FindQQQPos(TVector3 pos,
+                            TVector3 direction,
+                            bool verbose){
+
+  id   = -1;
+  chUp = -1;
+  chDn = -1;
+  chBk = -1;
+
+  //--------------------------------------------
+  // Intersect trajectory with QQQ plane
+  //--------------------------------------------
+
+  if( TMath::Abs(direction.Z()) < 1e-10 ) return;
+
+  double t = (qqqZPos - pos.Z()) / direction.Z();
+
+  if( t <= 0 ) return;
+
+  hitPos = pos + t * direction;
+
+  //--------------------------------------------
+  // Cylindrical coordinates
+  //--------------------------------------------
+
+  double x = hitPos.X();
+  double y = hitPos.Y();
+
+  double r = TMath::Sqrt(x*x + y*y);
+
+  if( r < qqqR1 || r > qqqR2 ) return;
+
+  double phi = hitPos.Phi() * TMath::RadToDeg();
+
+  if( phi < 0 ) phi += 360.0;
+
+  //--------------------------------------------
+  // Determine detector ID
+  //--------------------------------------------
+
+  id = -1;
+
+  for(int det = 0; det < 4; det++){
+
+    double phiMin = det*90.0 + 5.0;
+    double phiMax = phiMin + 85.0;
+
+    if( phi >= phiMin && phi <= phiMax ){
+      id = det;
+      break;
+    }
+  }
+
+  if( id < 0 ) return;
+
+  //--------------------------------------------
+  // Ring number (32 strips)
+  //--------------------------------------------
+
+  const double ringWidth =
+      (qqqR2 - qqqR1)/32.0;
+
+  int ring =
+      (int)((r - qqqR1)/ringWidth);
+
+  if( ring < 0 ) ring = 0;
+  if( ring > 31 ) ring = 31;
+
+  //--------------------------------------------
+  // Sector number (4 strips)
+  //--------------------------------------------
+
+  double localPhi =
+      phi - (id*90.0 + 5.0);
+
+  int sector =
+      (int)(localPhi/(85.0/4.0));
+
+  if( sector < 0 ) sector = 0;
+  if( sector > 3 ) sector = 3;
+
+  chBk = ring;
+  chDn = sector;
+  chUp = sector;
+
+  zFrac = 0.0;
+
+  if(verbose){
+
+    printf("\nQQQ Hit\n");
+    printf(" ID     = %d\n", id);
+    printf(" Ring   = %d\n", ring);
+    printf(" Sector = %d\n", sector);
+    printf(" r      = %.2f mm\n", r);
+    printf(" phi    = %.2f deg\n", phi);
+
+    hitPos.Print();
+  }
+}
+/*s
+inline TVector3 QQQ::GetHitPosWithSigma(double sigmaY_mm, double sigmaZ_mm){
+
+  double phi = SNorml[id%numDet].Phi();
+
+  TVector3 haha = hitPos;
+  haha.RotateZ(-phi);
+
+  double y = haha.Y() + gRandom->Gaus(0, sigmaY_mm);
+  if( sigmaY_mm < 0 ){
+    double deltaW = width/4;
+    y = TMath::Floor((haha.Y()-deltaW)/deltaW)*deltaW + deltaW*1.5; // when ever land on each strip, set the position to be center of the strip.
+    if( y >= 25 ) y = 15;
+  }
+
+  double z = haha.Z() + gRandom->Gaus(0, sigmaZ_mm);
+  if( sigmaZ_mm < 0 ){
+    haha.Z();
+    double delta = length/4;
+    int sign = z > 0 ? 1 : -1;
+    z = TMath::Floor( (abs(z)-gap/2)/delta )*delta + 0.5 * delta + gap/2;
+    if( z >= 107.375 ) z = 88.625;
+    z = sign * z;
+  }
+
+  haha.SetY(y);
+  haha.SetZ(z);
+  haha.RotateZ(phi);
+
+  return haha;
+
+}*/
+
+
+inline void QQQ::CalQQQPos(unsigned short ID,
+                           unsigned short chUp,
+                           unsigned short chDown,
+                           unsigned short chBack,
+                           float eUp,
+                           float eDown){
+
+  hitPos.Clear();
+
+  if( ID > 3 ) return;
+  if( chBack > 31 ) return;
+  if( chDown > 3 ) return;
+
+  const double ringWidth =
+      (qqqR2 - qqqR1)/32.0;
+
+  double r =
+      qqqR1 + (chBack + 0.5)*ringWidth;
+
+  const double sectorWidth =
+      85.0/4.0;
+
+  double phiDeg =
+      ID*90.0 + 5.0 +
+      (chDown + 0.5)*sectorWidth;
+
+  double phi =
+      phiDeg * TMath::DegToRad();
+
+  hitPos.SetXYZ(
+      r*TMath::Cos(phi),
+      r*TMath::Sin(phi),
+      qqqZPos
+  );
+
+  id   = ID;
+  chBk = chBack;
+  chDn = chDown;
+  chUp = chUp;
+}
+
+#endif 

Armory/anasenMS.cpp

@@ -11,6 +11,7 @@
 #include "TApplication.h" // ROOT app loop
 #include "ClassTransfer.h" // Reaction kinematics and MC event generation
 #include "ClassAnasen.h"   // ANASEN detector model classes (SX3, PW, etc.)
+#include "ClassQQQ.h"      // QQQ detector model class
 #include <stdio.h>
 #include <stdlib.h>
 #include <set>
@@ -160,6 +161,7 @@ int main(int argc, char **argv){
   ANASEN * anasen = new ANASEN();       // top-level detector object
   SX3 * sx3 = anasen->GetSX3();         // silicon array part
   PW * pw = anasen->GetPW();            // proportional wire chamber part
+  QQQ * qqq = anasen->GetQQQ();         // optional QQQ detector part, not used in this simulation but can be enabled for visualization
 
   // output file + trees
   TString saveFileName = "SimAnasen1.root";
@@ -192,8 +194,8 @@ int main(int argc, char **argv){
   tree2->Branch("phiCM",     &phiCM2, "phiCM/D");
 
   // outgoing particles in lab frame (light/heavy)
-  double thetab, phib, Tb;
+  double thetab, phib, Tb, qqqTb;
-  double thetaB, phiB, TB;
+  double thetaB, phiB, TB, qqqTB;
   std::array<double, 2> T;
   tree1->Branch("thetab", &thetab, "thetab/D"); // polar angle of light particle in lab frame
   tree1->Branch("phib",     &phib, "phib/D"); // azimuthal angle of light particle in lab frame
@@ -202,9 +204,11 @@ int main(int argc, char **argv){
   tree1->Branch("phiB",     &phiB, "phiB/D");
   tree1->Branch("TB",         &TB, "TB/D"); // kinetic energy of heavy particle at vertex (before energy loss)
   tree1->Branch("T",          &T, "T/D");   // placeholder for true Q-value, currently set to 0 for simplicity
+  tree1->Branch("qqqTb", &qqqTb, "qqqTb/D"); // kinetic energy of light particle at vertex (before energy loss) for events where the light particle hits the QQQ, currently set to 0 for simplicity
+  tree1->Branch("qqqTB", &qqqTB, "qqqTB/D"); // kinetic energy of heavy particle at vertex (before energy loss) for events where the light
 
-  double thetab2, phib2, Tb2;
+  double thetab2, phib2, Tb2, qqqTb2;
-  double thetaB2, phiB2, TB2;
+  double thetaB2, phiB2, TB2, qqqTB2;
   std::array<double, 2> T2;
   tree2->Branch("thetab", &thetab2, "thetab/D");
   tree2->Branch("phib",     &phib2, "phib/D");
@@ -213,6 +217,8 @@ int main(int argc, char **argv){
   tree2->Branch("phiB",     &phiB2, "phiB/D");
   tree2->Branch("TB",         &TB2, "TB/D");
   tree2->Branch("T",          &T2, "T/D");
+  tree2->Branch("qqqTb", &qqqTb2, "qqqTb/D");
+  tree2->Branch("qqqTB", &qqqTB2, "qqqTB/D");
 
   // excitation state identifiers
   int ExAID;
@@ -257,6 +263,16 @@ int main(int argc, char **argv){
   tree2->Branch("sx3Y",     &sx3Y2, "sx3Y/D");
   tree2->Branch("sx3Z",     &sx3Z2, "sx3Z/D");
 
+  double qqqX, qqqY, qqqZ;
+  tree1->Branch("qqqX",     &qqqX, "qqqX/D"); // reconstructed X position from QQQ (with optional smearing) in mm
+  tree1->Branch("qqqY",     &qqqY, "qqqY/D"); // reconstructed Y position from QQQ (with optional smearing)
+  tree1->Branch("qqqZ",     &qqqZ, "qqqZ/D"); // reconstructed Z position from QQQ (with optional smearing)
+
+  double qqqX2, qqqY2, qqqZ2;
+  tree2->Branch("qqqX",     &qqqX2, "qqqX/D");
+  tree2->Branch("qqqY",     &qqqY2, "qqqY/D");
+  tree2->Branch("qqqZ",     &qqqZ2, "qqqZ/D");
+
   // PW nearest and next nearest wires
   int anodeID[2], cathodeID[2];
   int anodeID2[2], cathodeID2[2];
@@ -274,16 +290,18 @@ int main(int argc, char **argv){
   tree2->Branch("cDist",  cathodeDist2, "cathodeDist/D");
 
   // SX3 channel assignment and Z fraction (depth) information
-  int sx3ID, sx3Up, sx3Dn, sx3Bk;
+  int sx3ID, sx3Up, sx3Dn, sx3Bk, qqqID;
   double sx3ZFrac;
-  int sx3ID2, sx3Up2, sx3Dn2, sx3Bk2;
+  int sx3ID2, sx3Up2, sx3Dn2, sx3Bk2, qqqID2;
   double sx3ZFrac2;
   tree1->Branch("sx3ID",   &sx3ID, "sx3ID/I");
+  tree1->Branch("qqqID",   &qqqID, "qqqID/I");
   tree1->Branch("sx3Up",   &sx3Up, "sx3Up/I");
   tree1->Branch("sx3Dn",   &sx3Dn, "sx3Dn/I");
   tree1->Branch("sx3Bk",   &sx3Bk, "sx3Bk/I");
   tree1->Branch("sx3ZFrac", &sx3ZFrac, "sx3ZFrac/D");
   tree2->Branch("sx3ID",   &sx3ID2, "sx3ID/I");
+  tree2->Branch("qqqID",   &qqqID2, "qqqID/I");
   tree2->Branch("sx3Up",   &sx3Up2, "sx3Up/I");
   tree2->Branch("sx3Dn",   &sx3Dn2, "sx3Dn/I");
   tree2->Branch("sx3Bk",   &sx3Bk2, "sx3Bk/I");
@@ -401,6 +419,7 @@ int main(int argc, char **argv){
     // run detector response models for PW and SX3
     pw->FindWireID(vertex, dir, false);
     sx3->FindSX3Pos(vertex, dir, false);
+    qqq->FindQQQPos(vertex, dir, false);
 
     PWHitInfo hitInfo = pw->GetHitInfo();
 
@@ -417,12 +436,14 @@ int main(int argc, char **argv){
 
     // SX3 hit channel info and depth fraction
     sx3ID = sx3->GetID();
+    qqqID = qqq->GetID();
 
     if(IsDeadSX3(sx3ID)) continue;
 
     anodeDist[0]   = hitInfo.nearestDist.first; // distance to nearest anode wire
     cathodeDist[0] = hitInfo.nearestDist.second; // distance to nearest cathode wire
     
+    //start HERE
     if( sx3ID >= 0 ){
       sx3Up    = sx3->GetChUp();
       sx3Dn    = sx3->GetChDn();
@@ -529,6 +550,82 @@ int main(int argc, char **argv){
 
       tree2->Fill();
 
+    }else if (qqqID >= 0){
+      // handle QQQ hit case
+      sx3Up = -1;
+      sx3Dn = -1;
+      sx3Bk = -1;
+      sx3ZFrac = TMath::QuietNaN();
+
+      sx3X = TMath::QuietNaN();
+      sx3Y = TMath::QuietNaN();
+      sx3Z = TMath::QuietNaN();
+
+      reTheta  = TMath::QuietNaN();
+      rePhi    = TMath::QuietNaN();
+      reTheta1 = TMath::QuietNaN();
+      rePhi1   = TMath::QuietNaN();
+      z0       = TMath::QuietNaN();
+
+      qqqTb = Tb; // for simplicity, using the same kinetic energy for QQQ hit events, can be modified to simulate energy loss if desired
+      qqqTB = TB;
+
+      Tb = TMath::QuietNaN(); // mark kinetic energy as invalid for SX3 hit case
+      TB = TMath::QuietNaN();
+
+      TVector3 hitPos = qqq->GetHitPos();
+
+      qqqX = hitPos.X();
+      qqqY = hitPos.Y();
+      qqqZ = hitPos.Z();
+
+      if(enableVis){
+        visTrackVertex.push_back(vertex);
+        visTrackDir.push_back(dir);
+        visTrackHitPos.push_back(hitPos);
+        //visTrackWires.push_back({anodeID[0], cathodeID[0]});
+      }
+
+      tree1->Fill();
+
+      TVector3 dir2(1, 0, 0);
+      dir2.SetTheta(thetab2 * TMath::DegToRad());
+      dir2.SetPhi(phib2 * TMath::DegToRad());
+
+      qqq->FindQQQPos(vertex, dir2, false);
+
+      if(qqqID2 < 0){
+        qqqID2 = -1;
+        qqqX2 = TMath::QuietNaN();
+        qqqY2 = TMath::QuietNaN();
+        qqqZ2 = TMath::QuietNaN();
+        anodeDist2[0] = TMath::QuietNaN();
+        cathodeDist2[0] = TMath::QuietNaN();
+        reTheta2 = TMath::QuietNaN();
+        rePhi2 = TMath::QuietNaN();
+        reTheta12 = TMath::QuietNaN();
+        rePhi12 = TMath::QuietNaN();
+        z02 = TMath::QuietNaN();
+        anodeID2[0]   = TMath::QuietNaN(); // no valid anode wire for QQQ hit case
+        cathodeID2[0] = TMath::QuietNaN(); // no valid cathode wire for QQQ hit case
+        anodeID2[1]   = TMath::QuietNaN(); // no valid next nearest anode wire for QQQ hit case
+        cathodeID2[1] = TMath::QuietNaN(); // no valid next nearest cathode wire for QQQ hit case
+        anodeDist2[1]   = TMath::QuietNaN();
+        cathodeDist2[1] = TMath::QuietNaN();
+      }
+
+      KEA2 = KEA;
+      thetaCM2 = thetaCM;
+      phiCM2 = phiCM;
+      ExAID2 = ExAID;
+      ExA2 = ExA;
+      ExID2 = ExID;
+      Ex2 = Ex;
+      vertexX2 = vertexX;
+      vertexY2 = vertexY;
+      vertexZ2 = vertexZ;
+
+      tree2->Fill();
     }else{
       // no valid SX3 hit: mark clearly invalid
       sx3Up = -1;
@@ -545,7 +642,8 @@ int main(int argc, char **argv){
       reTheta1 = TMath::QuietNaN();
       rePhi1   = TMath::QuietNaN();
       z0       = TMath::QuietNaN();
-      //Tb = -12354567; // mark kinetic energy as invalid for no hit case
+      Tb       = TMath::QuietNaN(); // mark kinetic energy as invalid for no hit case
+      TB       = TMath::QuietNaN();
       // fill tree with original data (no energy loss for these events)
       //comment out tree fill for no hit case  
       //tree1->Fill();