From 92d831678e7f34c04f2937c035f9548420bc41e4 Mon Sep 17 00:00:00 2001 From: James Szalkie Date: Thu, 4 Jun 2026 17:26:34 -0400 Subject: [PATCH] QQQ --- Armory/ClassAnasen.h | 6 +- Armory/ClassQQQ.h | 288 +++++++++++++++++++++++++++++++++++++++++++ Armory/anasenMS.cpp | 114 +++++++++++++++-- 3 files changed, 399 insertions(+), 9 deletions(-) create mode 100644 Armory/ClassQQQ.h diff --git a/Armory/ClassAnasen.h b/Armory/ClassAnasen.h index 549392e..aa86676 100644 --- a/Armory/ClassAnasen.h +++ b/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 wireID = pw->GetNearestID(); diff --git a/Armory/ClassQQQ.h b/Armory/ClassQQQ.h new file mode 100644 index 0000000..1ad1d88 --- /dev/null +++ b/Armory/ClassQQQ.h @@ -0,0 +1,288 @@ +#ifndef ClassQQQ_h +#define ClassQQQ_h + +#include +#include +#include +#include +#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 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(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 \ No newline at end of file diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index 96a344b..a9402cd 100644 --- a/Armory/anasenMS.cpp +++ b/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 #include #include @@ -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; + double thetab, phib, Tb, qqqTb; + double thetaB, phiB, TB, qqqTB; std::array 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; + double thetab2, phib2, Tb2, qqqTb2; + double thetaB2, phiB2, TB2, qqqTB2; std::array 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();