FEAT: add AmpTools Zlm and TF3 plot scripts

.github/workflows/amptools.yml

@@ -76,7 +76,20 @@ jobs:
           fi
           pushd extern/AmpTools; make -j$(nproc); popd
           make -j$(nproc)
+      - uses: actions/setup-python@v4
+        with:
+          python-version: "3.10"
+      - run: pip install .
       - name: Run executables
         run: |
-          cd script
-          ./print_amplitudes -c zlm.cfg
+          ./print_amplitudes -c zlm.cfg > .print_amplitudes.log
+          python ./makeComparison.py
+          ./draw_zlm
+        working-directory: script
+      - uses: actions/upload-artifact@v3
+        if: always()
+        with:
+          path: |
+            script/**.log
+            script/**.root
+            !script/data.root

.pre-commit-config.yaml

@@ -17,6 +17,11 @@ repos:
       - id: mixed-line-ending
       - id: trailing-whitespace
 
+  - repo: https://github.com/psf/black
+    rev: 23.7.0
+    hooks:
+      - id: black
+
   - repo: https://github.com/pre-commit/mirrors-clang-format
     rev: v16.0.6
     hooks:
@@ -26,6 +31,11 @@ repos:
           - c
           - cuda
 
+  - repo: https://github.com/PyCQA/isort
+    rev: 5.12.0
+    hooks:
+      - id: isort
+
   - repo: https://github.com/nbQA-dev/nbQA
     rev: 1.7.0
     hooks:
@@ -57,6 +67,6 @@ repos:
             metadata.vscode
 
   - repo: https://github.com/pre-commit/mirrors-prettier
-    rev: v3.0.0
+    rev: v3.0.1
     hooks:
       - id: prettier

.vscode/extensions.json

@@ -1,5 +1,6 @@
 {
   "recommendations": [
+    "albertopdrf.root-file-viewer",
     "eamodio.gitlens",
     "editorconfig.editorconfig",
     "esbenp.prettier-vscode",

pyproject.toml

@@ -10,7 +10,9 @@ dependencies = [
     "ipywidgets",
     "ipywidgets",
     "matplotlib",
+    "pandas",
     "sympy",
+    "uproot",
 ]
 name = "gluex-amplitude"
 requires-python = ">=3.7"

script/.gitignore

@@ -1,5 +1,8 @@
 *
+!.gitignore
 !*.cc
 !*.cfg
-!.gitignore
+!*.md
+!*.py
+!data.root
 !Makefile

script/Makefile

@@ -4,7 +4,7 @@ LIB_DIR := -L$(GLUEXAMPLITUDE)/lib
 
 # add extra library linking flags here if they are needed -- be sure the locations
 # are specified in the LIB_DIR arugments above
-EXTRA_LIBS := -lGlueXAmplitude
+EXTRA_LIBS := -lGlueXAmplitude -lAmpTools$(ATSUFFIX)
 
 ######  most things below probably don't need adjustment
 
@@ -25,7 +25,7 @@ CXX_FLAGS += $(shell root-config --cflags)
 # need some linking information available
 LIB_DIR += -L$(shell root-config --libdir) -L$(AMPTOOLS)/lib -L$(AMPPLOTTER)/lib
 
-LIBS :=  -lAmpTools$(ATSUFFIX) $(EXTRA_LIBS) $(shell root-config --glibs) -lstdc++
+LIBS := $(EXTRA_LIBS) $(shell root-config --glibs) -lstdc++
 
 # any source file falls into one of these two categories -- we'll
 # assume that anything that has MPI in the name is an MPI executable

script/README.md

@@ -0,0 +1,8 @@
+# Usage
+
+The data comes from the spring 2017 run, thrown monte-carlo simulations used in the etapi analysis
+
+```shell
+./print_amplitudes -c zlm.cfg > .print_amplitudes.log
+python ./makeComparison.py
+```

script/draw_zlm.cc

@@ -0,0 +1,94 @@
+#include "AMPTOOLS_AMPS/wignerD.h"
+#include "TAxis.h"
+#include "TCanvas.h"
+#include "TF3.h"
+#include "TH3F.h"
+#include "TLorentzRotation.h"
+#include "TLorentzVector.h"
+#include "TMath.h"
+#include "TVector3.h"
+
+// complex< GDouble > zlm(Double_t *xs, Int_t *pars)
+Double_t zlm(Double_t *xs, Double_t *pars) {
+  // Parameters
+  Int_t m_j = (Int_t)pars[0];        // spin
+  Int_t m_m = (Int_t)pars[1];        // projection
+  Int_t m_r = (Int_t)pars[2];        // reflectivity
+  Int_t m_s = (Int_t)pars[3];        // sign
+  Int_t pAngle = (Int_t)pars[4];     // polarization angle
+  GDouble pGamma = (GDouble)pars[5]; // polarization fraction
+  GDouble factor = sqrt(1 + m_s * pGamma);
+  // cout << "j: " << m_j << endl;
+  // cout << "m: " << m_m << endl;
+  // cout << "reflectivity: " << m_r << endl;
+  // cout << "sign: " << m_s << endl;
+  // cout << "pAngle: " << pAngle << endl;
+  // cout << "pGamma: " << pGamma << endl;
+  // cout << "factor: " << factor << endl;
+
+  // Variables
+  GDouble cosTheta = xs[0];
+  GDouble phi = xs[1];
+  GDouble bigPhi = xs[2];
+
+  // Calculate Intensity
+  GDouble zlm = 0;
+  complex<GDouble> rotateY = polar((GDouble)1., (GDouble)(-1. * bigPhi));
+  if (m_r == 1)
+    zlm = real(Y(m_j, m_m, cosTheta, phi) * rotateY);
+  if (m_r == -1)
+    zlm = imag(Y(m_j, m_m, cosTheta, phi) * rotateY);
+
+  complex<GDouble> amplitude(factor * zlm);
+  Double_t intensity = std::norm(amplitude);
+  return intensity;
+}
+
+int main(int argc, char *argv[]) {
+  cout << "Drawing Zlm Amplitude" << endl;
+  // Double_t xs[] = {0,0,0};
+  // Int_t pars[] = {2,0,1,1,0};
+  // cout << zlm(xs, pars) << endl; // Test
+
+  // Define Function
+  TF3 *function = new TF3("zlm", zlm, -1, 1, -3.14, 3.14, -3.14, 3.14, 6);
+  Int_t m_j = 2;
+  Int_t m_m = 0;
+  Int_t m_r = 1;
+  Int_t m_s = 1;
+  Int_t pAngle = 0;
+  GDouble pGamma = 1;
+  function->SetParameters(m_j, m_m, m_r, m_s, pAngle, pGamma);
+
+  // Draw Function
+  TCanvas canvas = new TCanvas("canvas", "canvas", 600, 600);
+  int nBins = 50;
+  string par_string(
+      Form("l:%i m:%i e:%i s:%i a:%i g:%f;cos#theta;#phi [deg];#Phi [deg]", m_j,
+           m_m, m_r, m_s, pAngle, pGamma));
+  TH3F *fcnHist = new TH3F("zlm", par_string.c_str(), nBins, -1, +1, nBins,
+                           -3.14, 3.14, nBins, -3.14, 3.14);
+  TAxis *xAxis = fcnHist->GetXaxis();
+  TAxis *yAxis = fcnHist->GetYaxis();
+  TAxis *zAxis = fcnHist->GetZaxis();
+  double x, y, z;
+  for (int xBin = 1; xBin < nBins + 1; ++xBin) {
+    for (int yBin = 1; yBin < nBins + 1; ++yBin) {
+      for (int zBin = 1; zBin < nBins + 1; ++zBin) {
+        x = xAxis->GetBinCenter(xBin);
+        y = yAxis->GetBinCenter(yBin);
+        z = zAxis->GetBinCenter(zBin);
+        fcnHist->SetBinContent(xBin, yBin, zBin, function->Eval(x, y, z));
+        // cout << function->Eval(x, y, z) << endl;
+      }
+    }
+  }
+  fcnHist->SetMinimum(0);
+  fcnHist->GetXaxis()->SetTitleOffset(1.5);
+  fcnHist->GetYaxis()->SetTitleOffset(2);
+  fcnHist->GetZaxis()->SetTitleOffset(1.5);
+  fcnHist->Draw("BOX2");
+  // canvas.SaveAs("zlm_3d.pdf");
+  canvas.SaveAs("zlm_3d.root");
+  return 0;
+}

script/makeComparison.py

@@ -0,0 +1,38 @@
+import re
+
+import numpy as np
+import uproot
+
+# Load the data from the root file
+data = uproot.open("data.root")["kin"].arrays(
+    ["cosTheta_eta_hel_thrown", "phi_eta_hel_thrown", "Phi_thrown"], library="pd"
+)
+
+# Load the intensities extracted by print_amplitudes.cc
+intensities = []
+with open(".print_amplitudes.log", "r") as file:
+    for line in file:
+        if "INTENSITY =" in line:
+            line = re.sub(" +", " ", line.rstrip().lstrip())
+            intensities.append(float(line.split()[2]))
+intensities = np.array(intensities)
+
+assert len(intensities) == len(data["cosTheta_eta_hel_thrown"])
+
+# Rename dictionary keys
+data["intensity"] = intensities
+data = data.rename(
+    {
+        "cosTheta_eta_hel_thrown": "theta",
+        "phi_eta_hel_thrown": "phi",
+        "Phi_thrown": "Phi",
+    },
+    axis=1,
+)
+
+# Standardize the data to be angles in radians
+DEGREE_TO_RADIAN = np.pi / 180.0
+data["theta"] = np.arccos(data["theta"])
+data["Phi"] = data["Phi"] * DEGREE_TO_RADIAN
+
+data.to_csv("comparison.csv", index=False, sep=",")

script/print_amplitudes.cc

@@ -9,6 +9,7 @@
 
 #include "TSystem.h"
 
+#include "AMPTOOLS_AMPS/ROOTDataReader.h"
 #include "AMPTOOLS_AMPS/Zlm.h"
 
 #include "IUAmpTools/AmpToolsInterface.h"
@@ -23,7 +24,6 @@ int main(int argc, char *argv[]) {
 
   for (int i = 1; i < argc; i++) {
     string arg(argv[i]);
-
     if (arg == "-c") {
       configfile = argv[++i];
     }
@@ -36,25 +36,45 @@ int main(int argc, char *argv[]) {
 
   ConfigFileParser parser(configfile);
   ConfigurationInfo *cfgInfo = parser.getConfigurationInfo();
-  cfgInfo->display();
+  string reactionName = cfgInfo->reactionList()[0]->reactionName();
+
+  // cfgInfo->display();
 
   AmpToolsInterface::registerAmplitude(Zlm());
+  AmpToolsInterface::registerDataReader(ROOTDataReader());
 
   AmpToolsInterface ati(cfgInfo);
-  vector<TLorentzVector> p4List;
 
-  // list lab frame p4s, in the order specified in the AmpTools config file
-  p4List.push_back(TLorentzVector(0, 0, 8.209470, 8.209470)); // beam photon
-  p4List.push_back(TLorentzVector(-0.291810, -0.315238, 0.695296,
-                                  1.244320)); // recoil proton
-  p4List.push_back(TLorentzVector(0.183147, 0.002392, 0.384435,
-                                  0.448137)); // resonance daughter 1
-  p4List.push_back(TLorentzVector(0.151104, 0.333728, 1.914300,
-                                  1.953710)); // resonance daughter 2
+  int iDataSet = 0; // if use used {data, accmc, genmc datasets}, we can choose
+                    // one of them here
 
-  Kinematics kin(p4List);
+  //////////////////////////////////////////////////////////////////////////
+  ///       TESTING PURPOSES
+  // vector< TLorentzVector > p4List;
+  ////// list lab frame p4s, in the order specified in the AmpTools config file
+  // p4List.push_back( TLorentzVector(0.000000,0.000000,8.419198,8.419198) );
+  // p4List.push_back( TLorentzVector(-0.105035,0.014094,1.106842,1.454883) );
+  // p4List.push_back( TLorentzVector(-0.330669,-0.777849,0.458704,0.971092) );
+  // p4List.push_back( TLorentzVector(0.435704,0.763756,6.855184,6.933026) );
+  // Kinematics kin_test( p4List );
+  // ati.printIntensity( reactionName, &kin_test );
+  //////////////////////////////////////////////////////////////////////////
 
-  ati.printEventDetails(cfgInfo->reactionList()[0]->reactionName(), &kin);
+  // Have to loadEvents before using kinematics (and numEvents) but cannot run
+  // processEvents.
+  //   To use the intensity method, we need to processEvents first
+  Kinematics kin;
+  ati.loadEvents(ati.dataReader(reactionName), iDataSet);
+  int nentries = ati.numEvents(iDataSet);
+  printf("Number of entries: %i\n", nentries);
+  // cout << ati.processEvents(reactionName,iDataSet) << endl;
+  for (int i = 0; i < nentries; ++i) {
+    kin = *ati.kinematics(i, iDataSet);
+    // for (auto four_vector: kin.particleList()){
+    //     four_vector.Print();
+    // }
+    ati.printIntensity(reactionName, &kin);
+  }
 
   return 0;
 }