QQQ tracking (set to alpha scatter)

Armory/anasenMS.cpp

@@ -158,19 +158,52 @@ int main(int argc, char **argv){
   std::signal(SIGINT, handler);
 
   //To set beam energy loss, use energy loss app, and create table with target isotope, set Initial beam energy as max energy 
-  transfer.SetA(27, 13, 0);            // 22Mg projectile
+  transfer.SetA(1, 1, 0);            // 22Mg projectile
   transfer.Seta(4, 2);                 // 4He target
-  transfer.Setb(1, 1);                 // outgoing proton from the primary transfer
+  transfer.Setb(4, 2);                 // outgoing proton from the primary transfer
-  transfer.SetB(30, 14);               // 30Si* heavy product
+  transfer.SetB(1, 1);               // 30Si* heavy product
-  double beamE = 56.1; //56.1
+  double beamE = 7; //56.1
   const ReactionConfig reactionConfig = transfer.GetRectionConfig();
   const double beamA = reactionConfig.beamA;                   // mass number of 14N beam
 
   // Excited state lists (projectile and heavy-product excitation states)
   std::vector<float> ExAList = {0};              // Beam excited energy
-  std::vector<float> ExList = {0, 2.235, 3.498, 3.774, 4.831, 4.809, 5.279, 5.614, 5.951, 6.550};  // a/p aluminum
+  std::vector<float> ExList = {0};//{0, 2.235, 3.498, 3.774, 4.831, 4.809, 5.279, 5.614, 5.951, 6.550};  // a/p aluminum
   //std::vector<float> ExList = {0.84, 1.01, 2.21, 2.73, 2.98, 3}; 
 
+  TGraph* elossBeam = LoadELoss("../ELoss/HeLoss/E_vs_x_proton.dat"); // x = path length (cm), y = beam energy (MeV)
+  TGraph* sigmaXBeam = LoadSigmaXVsEnergy("../ELoss/HeLoss/E_vs_x_proton.dat"); // x = beam energy (MeV), y = distance straggle sigma_x (cm)
+  //TGraph* sigmaABeam = LoadSigmaAVsEnergy("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = beam energy (MeV), y = distance straggle sigma_a (cm)
+  // Build a temporary inverse (energy -> path) to locate the path at beamE.
+  TGraph* elossBeamInverseRaw = new TGraph(elossBeam->GetN());
+  for( int p = 0; p < elossBeam->GetN(); p++ ){
+    double x, y;
+    elossBeam->GetPoint(p, x, y);
+    elossBeamInverseRaw->SetPoint(p, y, x);
+  }
+  elossBeamInverseRaw->Sort(); // TGraph::Eval requires ascending x
+
+  const double pathAtBeamE = elossBeamInverseRaw->Eval(beamE);
+
+  // Keep only energies <= beamE and shift path so beamE corresponds to x = 0.
+  TGraph* elossBeamFiltered = new TGraph();
+  for( int p = 0; p < elossBeam->GetN(); p++ ){
+    double x, y;
+    elossBeam->GetPoint(p, x, y);
+    if( y <= beamE ){
+      const int n = elossBeamFiltered->GetN();
+      elossBeamFiltered->SetPoint(n, x - pathAtBeamE, y);
+    }
+  }
+
+  TGraph* elossBeamInverse = new TGraph(elossBeamFiltered->GetN());
+  for( int p = 0; p < elossBeamFiltered->GetN(); p++ ){
+    double x, y;
+    elossBeamFiltered->GetPoint(p, x, y);
+    elossBeamInverse->SetPoint(p, y, x);
+  }
+  elossBeamInverse->Sort(); // TGraph::Eval requires ascending x
+
   const int kMBeam = reactionConfig.beamA; // mass number of beam
   const int kMTarget = reactionConfig.targetA; // mass number of target
   const int kMLight = reactionConfig.recoilLightA; // mass number of light ejectile
@@ -208,38 +241,7 @@ int main(int argc, char **argv){
   double vertexZRange[2] = { -174.3, 174.3}; // -174.3, 174.3 (full length of gas volume, centered at 0)
   const double beamEntranceZ = -280 - 174.3; //vertexZRange[0]; // mm, assumed beam entrance into the gas
   
-  TGraph* elossBeam = LoadELoss("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = path length (cm), y = beam energy (MeV)
-  TGraph* sigmaXBeam = LoadSigmaXVsEnergy("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = beam energy (MeV), y = distance straggle sigma_x (cm)
-  //TGraph* sigmaABeam = LoadSigmaAVsEnergy("../ELoss/HeLoss/E_vs_x_Al-27.dat"); // x = beam energy (MeV), y = distance straggle sigma_a (cm)
-  // Build a temporary inverse (energy -> path) to locate the path at beamE.
-  TGraph* elossBeamInverseRaw = new TGraph(elossBeam->GetN());
-  for( int p = 0; p < elossBeam->GetN(); p++ ){
-    double x, y;
-    elossBeam->GetPoint(p, x, y);
-    elossBeamInverseRaw->SetPoint(p, y, x);
-  }
-  elossBeamInverseRaw->Sort(); // TGraph::Eval requires ascending x
 
-  const double pathAtBeamE = elossBeamInverseRaw->Eval(beamE);
-
-  // Keep only energies <= beamE and shift path so beamE corresponds to x = 0.
-  TGraph* elossBeamFiltered = new TGraph();
-  for( int p = 0; p < elossBeam->GetN(); p++ ){
-    double x, y;
-    elossBeam->GetPoint(p, x, y);
-    if( y <= beamE ){
-      const int n = elossBeamFiltered->GetN();
-      elossBeamFiltered->SetPoint(n, x - pathAtBeamE, y);
-    }
-  }
-
-  TGraph* elossBeamInverse = new TGraph(elossBeamFiltered->GetN());
-  for( int p = 0; p < elossBeamFiltered->GetN(); p++ ){
-    double x, y;
-    elossBeamFiltered->GetPoint(p, x, y);
-    elossBeamInverse->SetPoint(p, y, x);
-  }
-  elossBeamInverse->Sort(); // TGraph::Eval requires ascending x
 
   // detector resolution / uncertainty parameters
   double sigmaSX3_W = -1; // mm, if < 0 use mid-point (no spread in SX3 horizontal dimension)
@@ -517,6 +519,11 @@ int main(int argc, char **argv){
     Ex_kin = TMath::QuietNaN();
     EBeam_track = TMath::QuietNaN();
     EBeam_kin = TMath::QuietNaN();
+    vX_track = TMath::QuietNaN();
+    vY_track = TMath::QuietNaN();
+    vZ_track = TMath::QuietNaN();
+    qqqUp = -1;
+    qqqBk = -1;
     sx3ID = -1;
     sx3Up = -1;
     sx3Dn = -1;
@@ -535,7 +542,7 @@ int main(int argc, char **argv){
     transfer.CalReactionConstant();
 
     // vertex position in target volume
-    beamEnergy = gRandom->Uniform(0, beamE); // MeV, sample beam energy at vertex from uniform distribution between 0 and initial beam energy
+    beamEnergy = gRandom->Uniform(6.56, beamE); // MeV, sample beam energy at vertex from uniform distribution between 0 and initial beam energy
     KEA = beamEnergy / beamA;
     beamPath_cm = elossBeamInverse->Eval(beamEnergy); // beamE maps to x=0 after path shift
     vertexZ = beamEntranceZ + beamPath_cm * 10.0; // cm -> mm
@@ -801,7 +808,6 @@ int main(int argc, char **argv){
 
         }
 
-
       } else {
         theta_recon = TMath::QuietNaN();
         vX_kin = TMath::QuietNaN();
@@ -826,6 +832,8 @@ int main(int argc, char **argv){
       qqqX = hitPos.X();
       qqqY = hitPos.Y();
       qqqZ = hitPos.Z();
+      qqqUp = qqq->GetChUp();
+      qqqBk = qqq->GetChBk();
 
         auto [cathodeIntersection, anodeIntersection, pcPathLengthCm] =
           find_PC_PathLength(vertex, hitPos);
@@ -864,10 +872,92 @@ int main(int argc, char **argv){
 
       if (Eqqq <= 0 || Eanode <= 0 || Ecathode <= 0) {
           Eqqq = NAN;
+          beamEnergy = NAN;
           Ex = NAN;
+          qqqID = -1;
           //continue;
       }
 
+      // PW track fits from QQQ hit
+      pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false);
+      reTheta = pw->GetTrackTheta() * TMath::RadToDeg();
+      rePhi   = pw->GetTrackPhi() * TMath::RadToDeg();
+      pw->CalTrack2(hitPos, hitInfo, sigmaPW_A, sigmaPW_C, false);
+      reTheta1 = pw->GetTrackTheta() * TMath::RadToDeg();
+      rePhi1   = pw->GetTrackPhi() * TMath::RadToDeg();
+      z0 = pw->GetZ0();
+
+      // QQQ sits at fixed z: smear in-plane (radial, azimuthal), z stays fixed
+      const double qHitPhi = hitPos.Phi();
+      const double qSigR = 2.0, qSigPhi = 2.0; // mm
+      const double dR   = gRandom->Gaus(0, qSigR);
+      const double dTan = gRandom->Gaus(0, qSigPhi);
+      const double qHitX = hitPos.X() + dR * std::cos(qHitPhi) - dTan * std::sin(qHitPhi);
+      const double qHitY = hitPos.Y() + dR * std::sin(qHitPhi) + dTan * std::cos(qHitPhi);
+      const double qHitZ = hitPos.Z();
+
+      const double qAX = gRandom->Gaus(aX, 2.3);
+      const double qAY = gRandom->Gaus(aY, 2.3);
+      const double qAZ = gRandom->Gaus(aZ, 2.3);
+
+      if (qqqID >= 0 && !std::isnan(aX) && !std::isnan(aY) && !std::isnan(aZ) && aZ != 54321 && (aX != 0 || aY != 0 || aZ != 0)) {
+        // line through QQQ hit and anode point; vertex = closest approach to the beam (z) axis
+        TVector3 x1(qHitX, qHitY, qHitZ);
+        TVector3 v = TVector3(qAX, qAY, qAZ) - x1;
+        const double vT2 = v.X()*v.X() + v.Y()*v.Y();
+
+        if (vT2 > 0.0) {
+          const double tMin = -(x1.X()*v.X() + x1.Y()*v.Y()) / vT2;
+          TVector3 rVtx = x1 + tMin * v;
+
+          vX_track = rVtx.X();
+          vY_track = rVtx.Y();
+          vZ_track = rVtx.Z();
+
+          EBeam_track = elossBeam->Eval((vZ_track + 454.3) / 10);
+
+          Tb_recon = CalculateOriginalEnergy(qHitX, qHitY, qHitZ,
+                                              vX_track, vY_track, vZ_track,
+                                              b, "He", Eqqq,
+                                              distance_qqq);
+
+          // polar angle of the track relative to the beam (z) axis
+          const TVector3 projectileDirection = x1 - TVector3(qAX, qAY, qAZ);
+          theta_recon = projectileDirection.Theta() * 180.0 / TMath::Pi();
+
+          if (reactionConfig.recoilLightA == 1 && !TMath::IsNaN(EBeam_track)) {
+            Kinematics apkin_27Al(26.981538408,4.00260325413,1.00782503224,29.973770136,EBeam_track/26.981538408);
+            EBeam_kin = apkin_27Al.getEbeam_givenQ(Tb_recon, Ex, theta_recon);
+            vZ_kin = elossBeamInverse->Eval(EBeam_kin) * 10 - 454.3;
+            Ex_kin = apkin_27Al.getExc(Tb_recon, theta_recon);
+          } else if (reactionConfig.recoilLightA == 4) {
+            double beamMomentum = std::sqrt((EBeam_track + al27Iso.Mass) * (EBeam_track + al27Iso.Mass) - al27Iso.Mass * al27Iso.Mass);
+            TLorentzVector beam4(0.0, 0.0, beamMomentum, EBeam_track + al27Iso.Mass);
+            TLorentzVector helium4(0, 0, 0, alphaIso.Mass);
+
+            double alphaEnergy = Tb_recon + alphaIso.Mass;
+            double alphaMomentum = std::sqrt(alphaEnergy * alphaEnergy - alphaIso.Mass * alphaIso.Mass);
+            TLorentzVector alpha4(alphaMomentum * std::sin(theta_recon * TMath::Pi() / 180.0),
+                                  0,
+                                  alphaMomentum * std::cos(theta_recon * TMath::Pi() / 180.0),
+                                  alphaEnergy);
+            TLorentzVector recoil4 = beam4 + helium4 - alpha4;
+            Ex_kin = recoil4.M() - al27Iso.Mass;
+
+            if (std::isfinite(Ex) && std::isfinite(Tb_recon) && std::isfinite(theta_recon)) {
+              EBeam_kin = InvertToBeamEnergy(Ex, Tb_recon, theta_recon, al27Iso, alphaIso);
+            } else {
+              EBeam_kin = TMath::QuietNaN();
+            }
+
+            vZ_kin = std::isfinite(EBeam_kin) && EBeam_kin >= 0.0
+                         ? elossBeamInverse->Eval(EBeam_kin) * 10 - 454.3
+                         : TMath::QuietNaN();
+          }
+        }
+      }
+      hitTheta = hitPos.Theta() * 180.0 / TMath::Pi();
+
       qqqTb = Tb; // for simplicity, using the same kinetic energy for QQQ hit events, can be modified to simulate energy loss if desired
       qqqTB = TB;
       Edet = Eqqq;
@@ -875,7 +965,6 @@ int main(int argc, char **argv){
       AutoHist2D::Fill("beamEnergy_vs_vZ", vertexZ / 10, beamEnergy, "vZ (cm)", "beamEnergy (MeV)");
       AutoHist2D::Fill("EPC x sin(theta) vs Eqqq", Eqqq, EPC * sin(thetab * TMath::DegToRad()), "Eqqq (MeV)", "EPC x sin(theta) (MeV)");
 
-      beamEnergy = TMath::QuietNaN(); // mark beam energy as invalid for QQQ hit case
       //tree1->Fill();
 
     }else{

Armory/testMacro.C

@@ -32,7 +32,7 @@
 	TH2F *h5 = new TH2F("Excitation Energy vs vZ", "Excitation Energy vs vZ;vZ (mm);Excitation energy (MeV)", 300, 0, 0, 300, 0, 0); //arguments are (name, title, nbinsX, xlow, xup, nbinsY, ylow, yup)
 	tree1->Draw("Ex_recon:vZ >> hExReconVsVZ", "sx3ID >= 0", "colz");
 	tree1->Draw("Ex:vZ >> h5", "sx3ID >= 0 && EPC >= 0.03", "same");
-*/
+
 
 	new TCanvas("c1", "Canvas", 900, 600);
 	tree1->Draw("Ex_kin:vZ_track>>h2(800,-500,300,110,-1,10)", "sx3ID >=0 && EPC > .02", "colz");
@@ -76,6 +76,11 @@
 	tree1->GetHistogram()->SetLineStyle(2);
 	tree1->Draw("Ex_kin", "sx3ID >= 0 && EPC > .02", "same");
 	c7->SaveAs("Ex_kin.png");
+*/
 
-
+	new TCanvas("c3", "Canvas", 900, 600);
+	tree1->Draw("vZ_track:vZ", "qqqID >=0 && EPC > 0", "");
+	TF1 *line_track = new TF1("line_track", "x", gPad->GetUxmin(), gPad->GetUxmax());
+	line_track->SetLineColor(kRed); line_track->Draw("SAME");
+	c3->SaveAs("vZ_track_vs_vZ.png");
 }