add benchmarks/u_rho

.gitlab-ci.yml

@@ -117,6 +117,7 @@ include:
   - local: 'benchmarks/dvcs/config.yml'
   - local: 'benchmarks/tcs/config.yml'
   - local: 'benchmarks/u_omega/config.yml'
+  - local: 'benchmarks/u_rho/config.yml'
   - local: 'benchmarks/single/config.yml'
   - local: 'benchmarks/backgrounds/config.yml'
 

Snakefile

@@ -18,3 +18,4 @@ root -l -b -q -e '.L {input}+'
 """
 
 include: "benchmarks/diffractive_vm/Snakefile"
+include: "benchmarks/u_rho/Snakefile"

benchmarks/u_rho/Snakefile

@@ -0,0 +1,91 @@
+import os
+
+from snakemake.remote.S3 import RemoteProvider as S3RemoteProvider
+
+
+S3 = S3RemoteProvider(
+    endpoint_url="https://dtn01.sdcc.bnl.gov:9000",
+    access_key_id=os.environ["S3_ACCESS_KEY"],
+    secret_access_key=os.environ["S3_SECRET_KEY"],
+)
+
+
+rule uchannelrho_compile:
+    input:
+        ROOT_BUILD_DIR_PREFIX + "benchmarks/u_rho/analysis/uchannelrho_cxx.so",
+
+
+rule uchannelrho_campaign_reco_get:
+    input:
+        lambda wildcards: S3.remote(f"eictest/EPIC/RECO/23.12.0/epic_craterlake/EXCLUSIVE/UCHANNEL_RHO/10x100/rho_10x100_uChannel_Q2of0to10_hiDiv.{wildcards.INDEX}.eicrecon.tree.edm4eic.root"),
+    output:
+        "benchmark_output/campaign_23.12.0_rho_10x100_uChannel_Q2of0to10_hiDiv_{INDEX}_eicrecon.edm4eic.root",
+    shell:
+        """
+echo "Getting file for INDEX {wildcards.INDEX}"
+ln {input} {output}
+"""
+
+
+rule uchannelrho_analysis:
+    input:
+        script="benchmarks/u_rho/analysis/uchannelrho.cxx",
+        #script_compiled=ROOT_BUILD_DIR_PREFIX + "benchmarks/u_rho/analysis/uchannelrho_cxx.so",
+        data="benchmark_output/campaign_23.12.0_rho_10x100_uChannel_Q2of0to10_hiDiv_{INDEX}_eicrecon.edm4eic.root",
+    output:
+        plots="benchmark_output/campaign_23.12.0_{INDEX}_eicrecon.edm4eic/plots.root",
+    shell:
+        """
+mkdir -p $(dirname "{output.plots}")
+root -l -b -q '{input.script}+("{input.data}","{output.plots}")'
+"""
+
+
+rule uchannelrho_combine:
+    input:
+        #lambda wildcards: [f"benchmark_output/campaign_23.12.0_{ix:04d}_eicrecon.edm4eic/plots.root" for ix in range(int(wildcards.NUM_FILES))],
+        lambda wildcards: expand(
+           "benchmark_output/campaign_23.12.0_{INDEX:04d}_eicrecon.edm4eic/plots.root",
+           INDEX=range(int(wildcards.N)),
+        ),	
+    wildcard_constraints:
+        N="\d+",
+    output:
+        "benchmark_output/campaign_23.12.0_combined_{N}files_eicrecon.edm4eic.plots.root",
+    shell:
+        """
+hadd {output} {input}
+"""
+
+
+ruleorder: uchannelrho_combine > uchannelrho_analysis
+
+
+rule uchannelrho_plots:
+    input:
+        script="benchmarks/u_rho/macros/plot_rho_physics_benchmark.C",
+        plots="benchmark_output/campaign_23.12.0_combined_{N}files_eicrecon.edm4eic.plots.root",
+    output:
+        "benchmark_output/campaign_23.12.0_combined_{N}files_eicrecon.edm4eic.plots_figures/benchmark_rho_mass.pdf",
+    shell:
+        """
+if [ ! -d "{input.plots}_figures" ]; then
+    mkdir "{input.plots}_figures"
+    echo "{input.plots}_figures directory created successfully."
+else
+    echo "{input.plots}_figures directory already exists."
+fi
+root -l -b -q '{input.script}("{input.plots}")'
+cat *.json
+#cp {input.plots}_figures/*.pdf figures/
+"""
+
+
+# Couple examples of invocation:
+
+rule uchannelrho_run_over_a_campaign:
+    input:
+        "benchmark_output/figures/plots_benchmark_rho_dsigmadt.pdf",
+    message:
+        "See output in {input[0]}"
+

benchmarks/u_rho/analysis/pleaseIncludeMe.h

@@ -0,0 +1,116 @@
+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <vector>
+#include <algorithm>
+#include <utility>
+
+#include "ROOT/RDataFrame.hxx"
+#include <TH1D.h>
+#include <TFitResult.h>
+#include <TRandom3.h>
+#include <TCanvas.h>
+#include <TSystem.h>
+
+#include "TFile.h"
+#include "TLorentzVector.h"
+
+#include "TLorentzRotation.h"
+#include "TVector2.h"
+#include "TVector3.h"
+
+#include "fmt/color.h"
+#include "fmt/core.h"
+
+// #include "nlohmann/json.hpp"
+#include "edm4eic/InclusiveKinematicsData.h"
+#include "edm4eic/ReconstructedParticleData.h"
+#include "edm4hep/MCParticleData.h"
+
+#define PI            3.1415926
+#define MASS_ELECTRON 0.00051
+#define MASS_PROTON   0.93827
+#define MASS_PION     0.13957
+#define MASS_KAON     0.493667
+#define MASS_AU197    183.45406466643374
+
+auto getNtrk(const std::vector<edm4eic::ReconstructedParticleData>& parts)
+{
+  std::vector<int> mult;
+  int n=0;
+  for(auto& i1 : parts){
+    if(i1.charge!=0) n++;
+  }
+  mult.push_back( n );
+  return mult;
+}
+auto getNtrkMC(const std::vector<edm4hep::MCParticleData>& parts)
+{
+  std::vector<int> mult;
+  int n=0;
+  for(auto& i1 : parts){
+    if(i1.charge!=0 && i1.generatorStatus==1) n++;
+  }
+  mult.push_back( n );
+  return mult;
+}
+auto momenta_from_chargedparticles(const std::vector<edm4eic::ReconstructedParticleData>& parts) {
+  std::vector<TVector3> momenta;
+  for(auto& i1 : parts){
+		TVector3 trk(i1.momentum.x,i1.momentum.y,i1.momentum.z);
+		if(i1.charge!=0) momenta.push_back(trk);
+	}
+  return momenta;
+}
+auto momenta_from_mcparticles(const std::vector<edm4hep::MCParticleData>& parts) {
+  std::vector<TVector3> momenta;
+  for(auto& i1 : parts){
+		TVector3 trk(i1.momentum.x,i1.momentum.y,i1.momentum.z);
+		if(i1.charge!=0 && i1.generatorStatus==1) momenta.push_back(trk);
+	}
+  return momenta;
+}
+auto pt_resolution(const std::vector<edm4hep::MCParticleData>& mcs,
+							const std::vector<edm4eic::ReconstructedParticleData>& recos){
+
+	std::vector<double> resolution;
+	for(auto& i1: recos){
+		TVector3 trk(i1.momentum.x,i1.momentum.y,i1.momentum.z);
+		if(i1.charge==0) continue;
+		double minR=99;
+		TVector3 matchMCtrk(-99,-99,-99);
+		for(auto& i2 : mcs){
+			TVector3 trkMC(i2.momentum.x,i2.momentum.y,i2.momentum.z);
+			if(i2.charge!=0 && i2.generatorStatus==1){
+				if(trk.DeltaR(trkMC) < minR ){
+					minR=trk.DeltaR(trkMC);
+					matchMCtrk=trkMC;
+				}
+			}
+		}
+		double res= (matchMCtrk.Perp()-trk.Perp()) / matchMCtrk.Perp();
+		resolution.push_back( res );
+
+	}
+
+	return resolution;
+}
+auto getPt(const std::vector<TVector3>& tracks)
+{
+	std::vector<double> Pt;
+	for(auto& i1 : tracks){Pt.push_back(i1.Perp());}
+	return Pt;
+}
+auto getEta(const std::vector<TVector3>& tracks)
+{
+	std::vector<double> eta;
+	for(auto& i1 : tracks){eta.push_back(i1.Eta());}
+	return eta;
+}
+auto getPhi(const std::vector<TVector3>& tracks)
+{
+	std::vector<double> Phi;
+	for(auto& i1 : tracks){Phi.push_back(i1.Phi());}
+	return Phi;
+}