add diffractive_vmp

.gitlab-ci.yml

@@ -108,6 +108,7 @@ include:
   - local: 'benchmarks/dis/config.yml'
     #- local: 'benchmarks/dvmp/config.yml'
   - local: 'benchmarks/dvcs/config.yml'
+  - local: 'benchmarks/diffractive_vmp/config.yml'
   - local: 'benchmarks/tcs/config.yml'
   - local: 'benchmarks/u_omega/config.yml'
   - local: 'benchmarks/single/config.yml'

benchmarks/diffractive_vmp/analysis/simple_analysis.cxx

+#include <cmath>
+#include <fstream>
+#include <iostream>
+#include <string>
+#include <vector>
+#include <algorithm>
+#include <utility>
+
+#include <TH1D.h>
+#include <TH2D.h>
+#include <TFitResult.h>
+#include <TRandom3.h>
+#include <TCanvas.h>
+#include <TTreeReader.h>
+#include <TTreeReaderArray.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TLorentzVector.h>
+
+#include <TLorentzRotation.h>
+#include <TVector2.h>
+#include <TVector3.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 giveme_t_method_L(TLorentzVector eIn, 
+					   TLorentzVector eOut, 
+					   TLorentzVector pIn, 
+					   TLorentzVector vmOut)
+{
+	TLorentzVector aInVec(pIn.Px()*197,pIn.Py()*197,pIn.Pz()*197,sqrt(pIn.Px()*197*pIn.Px()*197 + pIn.Py()*197*pIn.Py()*197 + pIn.Pz()*197*pIn.Pz()*197 + MASS_AU197*MASS_AU197) );
+	double method_L = 0;
+	TLorentzVector a_beam_scattered = aInVec-(vmOut+eOut-eIn);
+	double p_Aplus = a_beam_scattered.E()+a_beam_scattered.Pz();
+	double p_TAsquared = TMath::Power(a_beam_scattered.Pt(),2);
+	double p_Aminus = (MASS_AU197*MASS_AU197 + p_TAsquared) / p_Aplus;
+	TLorentzVector a_beam_scattered_corr; 
+	a_beam_scattered_corr.SetPxPyPzE(a_beam_scattered.Px(),a_beam_scattered.Py(),(p_Aplus-p_Aminus)/2., (p_Aplus+p_Aminus)/2. );
+	method_L = -(a_beam_scattered_corr-aInVec).Mag2();
+
+	return method_L;
+}
+
+
+int diffractive_vm_simple_analysis(TString rec_file, TString outputfile)
+{	
+// read our configuration	
+TString name_of_input = (TString) rec_file;
+std::cout << "Input file = " << name_of_input << endl;
+auto tree = new TChain("events");
+tree->Add(name_of_input);
+TTreeReader tree_reader(tree);       // !the tree reader
+
+TTreeReaderArray<int> mc_genStatus_array = {tree_reader, "MCParticles.generatorStatus"};
+// MC particle pz array for each MC particle
+TTreeReaderArray<float> mc_px_array = {tree_reader, "MCParticles.momentum.x"};
+TTreeReaderArray<float> mc_py_array = {tree_reader, "MCParticles.momentum.y"};
+TTreeReaderArray<float> mc_pz_array = {tree_reader, "MCParticles.momentum.z"};
+TTreeReaderArray<double> mc_mass_array = {tree_reader, "MCParticles.mass"};
+TTreeReaderArray<int> mc_pdg_array = {tree_reader, "MCParticles.PDG"};
+
+//Reconstructed EcalEndcapNClusters
+TTreeReaderArray<float> em_energy_array = {tree_reader, "EcalEndcapNClusters.energy"};
+TTreeReaderArray<float> em_x_array = {tree_reader, "EcalEndcapNClusters.position.x"};
+TTreeReaderArray<float> em_y_array = {tree_reader, "EcalEndcapNClusters.position.y"};
+TTreeReaderArray<float> emhits_x_array = {tree_reader, "EcalEndcapNRecHits.position.x"};
+TTreeReaderArray<float> emhits_y_array = {tree_reader, "EcalEndcapNRecHits.position.y"};
+TTreeReaderArray<float> emhits_energy_array = {tree_reader, "EcalEndcapNRecHits.energy"};
+
+TTreeReaderArray<unsigned int> em_rec_id_array = {tree_reader, "EcalEndcapNClustersAssociations.recID"};
+TTreeReaderArray<unsigned int> em_sim_id_array = {tree_reader, "EcalEndcapNClustersAssociations.simID"};
+
+// Reconstructed particles pz array for each reconstructed particle
+TTreeReaderArray<float> reco_px_array = {tree_reader, "ReconstructedChargedParticles.momentum.x"};
+TTreeReaderArray<float> reco_py_array = {tree_reader, "ReconstructedChargedParticles.momentum.y"};
+TTreeReaderArray<float> reco_pz_array = {tree_reader, "ReconstructedChargedParticles.momentum.z"};
+TTreeReaderArray<float> reco_charge_array = {tree_reader, "ReconstructedChargedParticles.charge"};
+
+TTreeReaderArray<unsigned int> rec_id = {tree_reader, "ReconstructedChargedParticlesAssociations.recID"};
+TTreeReaderArray<unsigned int> sim_id = {tree_reader, "ReconstructedChargedParticlesAssociations.simID"};
+
+TString output_name_dir = outputfile+"_output.root";
+cout << "Output file = " << output_name_dir << endl;
+TFile* output = new TFile(output_name_dir,"RECREATE");
+
+//events
+TH1D* h_Q2_e = new TH1D("h_Q2_e",";Q^{2}_{e,MC}",100,0,20);
+TH1D* h_y_e = new TH1D("h_y_e",";y_{e,MC}",100,0,1);
+TH1D* h_energy_MC = new TH1D("h_energy_MC",";E_{MC} (GeV)",100,0,20);
+TH1D* h_t_MC = new TH1D("h_t_MC",";t_{MC}; counts",100,0,0.2);
+
+TH1D* h_Q2REC_e = new TH1D("h_Q2REC_e",";Q^{2}_{e,REC}",100,0,20);
+TH1D* h_yREC_e = new TH1D("h_yREC_e",";y_{e,REC}",100,0,1);
+TH1D* h_energy_REC = new TH1D("h_energy_REC",";E_{REC} (GeV)",100,0,20);
+TH1D* h_trk_energy_REC = new TH1D("h_trk_energy_REC",";E_{REC} (GeV)",100,0,20);
+TH1D* h_trk_Epz_REC = new TH1D("h_trk_Epz_REC",";E - p_{z} (GeV)",200,0,50);
+
+//track
+TH1D* h_eta = new TH1D("h_eta",";#eta",100,-5,5);
+TH2D* h_trk_energy_res = new TH2D("h_trk_energy_res",";E_{MC} (GeV); E_{MC}-E_{REC}/E_{MC} track-base ",100,0,20,1000,-1,1);
+TH2D* h_trk_Pt_res = new TH2D("h_trk_Pt_res",";p_{T,MC} (GeV); P_{T,MC}-P_{T,REC}/P_{T,MC} track-base ",100,0,15,1000,-1,1);
+TH1D* h_Epz_REC = new TH1D("h_Epz_REC",";E - p_{z} (GeV)",200,0,50);
+
+//VM & t
+TH1D* h_VM_mass_REC = new TH1D("h_VM_mass_REC",";mass (GeV)",200,0,4);
+TH1D* h_VM_pt_REC = new TH1D("h_VM_pt_REC",";p_{T} (GeV/c)",200,0,2);
+TH2D* h_VM_res = new TH2D("h_VM_res",";p_{T,MC} (GeV); p_{T,MC}-E_{T,REC}/p_{T,MC}",100,0,2,1000,-1,1);
+TH1D* h_t_REC = new TH1D("h_t_REC",";t_{REC} (GeV^{2}); counts",100,0,0.2);
+TH1D* h_t_trk_REC = new TH1D("h_t_trk_REC",";t_{REC}(GeV^{2}) track-base; counts",100,0,0.2);
+TH1D* h_t_combo_REC = new TH1D("h_t_combo_REC",";t_{combo,REC}(GeV^{2}); counts",100,0,0.2);
+TH2D* h_t_res = new TH2D("h_t_res",";t_{MC} (GeV^{2}); t_{MC}-t_{REC}/t_{MC}",100,0,0.2,1000,-10,10);
+TH2D* h_trk_t_res = new TH2D("h_trk_t_res",";t_{MC} (GeV^{2}); t_{MC}-t_{REC}/t_{MC} track-base",100,0,0.2,1000,-10,10);
+TH2D* h_t_2D = new TH2D("h_t_2D",";t_{MC} (GeV^{2}); t_{REC} (GeV^{2}) track-base",100,0,0.2,100,0,0.2);
+TH2D* h_t_REC_2D = new TH2D("h_t_REC_2D",";t_{trk,REC} (GeV^{2}); t_{EEMC,REC} (GeV^{2})",100,0,0.2,100,0,0.2);
+TH2D* h_t_RECMC_2D = new TH2D("h_t_RECMC_2D",";t_{MC} (GeV^{2}); t_{trk,REC} / t_{EEMC,REC} ",100,0,0.2,200,-10,10);
+
+//energy clus
+TH2D* h_emClus_position_REC = new TH2D("h_emClus_position_REC",";x (mm);y (mm)",80,-800,800,80,-800,800);
+TH2D* h_emHits_position_REC = new TH2D("h_emHits_position_REC",";x (mm);y (mm)",80,-800,800,80,-800,800);
+TH2D* h_energy_res = new TH2D("h_energy_res",";E_{MC} (GeV); E_{MC}-E_{REC}/E_{MC} emcal",100,0,20,1000,-1,1);
+TH1D* h_energy_calibration_REC = new TH1D("h_energy_calibration_REC",";E (GeV)",200,0,2);
+TH1D* h_EoverP_REC = new TH1D("h_EoverP_REC",";E/p",200,0,2);
+TH1D* h_ClusOverHit_REC = new TH1D("h_ClusOverHit_REC",";cluster energy / new cluster energy",200,0,2);
+
+tree_reader.SetEntriesRange(0, tree->GetEntries());
+while (tree_reader.Next()) {
+
+	/*
+	Beam particles
+	*/
+	TLorentzVector ebeam(0,0,0,0);
+	TLorentzVector pbeam(0,0,0,0);
+
+	TLorentzVector vmMC(0,0,0,0);
+	TLorentzVector kplusMC(0,0,0,0);
+	TLorentzVector kminusMC(0,0,0,0);
+
+	//MC level
+	TLorentzVector scatMC(0,0,0,0);
+	int mc_elect_index=-1;
+	double maxPt=-99.;
+	for(int imc=0;imc<mc_px_array.GetSize();imc++){
+		TVector3 mctrk(mc_px_array[imc],mc_py_array[imc],mc_pz_array[imc]);	
+		if(mc_genStatus_array[imc]==4){//4 is Sartre.
+			if(mc_pdg_array[imc]==11) ebeam.SetVectM(mctrk, MASS_ELECTRON);
+				if(mc_pdg_array[imc]==2212) pbeam.SetVectM(mctrk, MASS_PROTON);
+		}
+		if(mc_genStatus_array[imc]!=1) continue;
+		if(mc_pdg_array[imc]==11 	
+			&& mctrk.Perp()>maxPt){
+			maxPt=mctrk.Perp();
+			mc_elect_index=imc;
+			scatMC.SetVectM(mctrk,mc_mass_array[imc]);
+		}
+		if(mc_pdg_array[imc]==321
+			&& mc_genStatus_array[imc]==1) kplusMC.SetVectM(mctrk,MASS_KAON);
+		if(mc_pdg_array[imc]==-321
+			&& mc_genStatus_array[imc]==1) kminusMC.SetVectM(mctrk,MASS_KAON);
+
+	}
+	vmMC=kplusMC+kminusMC;
+	//protection.
+	if(ebeam.E()==pbeam.E() && ebeam.E()==0) {
+		std::cout << "problem with MC incoming beams" << std::endl;
+		continue;
+	}
+	TLorentzVector qbeam=ebeam-scatMC;
+	double Q2=-(qbeam).Mag2();  
+	double pq=pbeam.Dot(qbeam);
+	double y= pq/pbeam.Dot(ebeam);
+	
+	//MC level phase space cut
+	if(Q2<1.||Q2>10.) continue;
+	if(y<0.01||y>0.85) continue;
+
+	h_Q2_e->Fill(Q2);
+	h_y_e->Fill(y);
+	h_energy_MC->Fill(scatMC.E());
+
+	double t_MC=0.;
+	if(vmMC.E()!=0 
+		&& fabs(vmMC.Rapidity())<3.5)
+	{
+		double method_E = -(qbeam-vmMC).Mag2();
+		t_MC=method_E;
+		h_t_MC->Fill( method_E );
+	}
+
+	//rec level
+	//leading cluster
+	double maxEnergy=-99.;
+	double xpos=-999.;
+	double ypos=-999.;
+	for(int iclus=0;iclus<em_energy_array.GetSize();iclus++){
+		if(em_energy_array[iclus]>maxEnergy){
+			maxEnergy=em_energy_array[iclus];
+			xpos=em_x_array[iclus];
+			ypos=em_y_array[iclus];
+		}
+	}
+	//leading hit energy
+	double maxHitEnergy=0.01;//threshold 10 MeV
+	double xhitpos=-999.;
+	double yhitpos=-999.;
+	int hit_index=-1;
+	for(int ihit=0;ihit<emhits_energy_array.GetSize();ihit++){	
+		if(emhits_energy_array[ihit]>maxHitEnergy){
+			maxHitEnergy=emhits_energy_array[ihit];
+			        xhitpos=emhits_x_array[ihit];
+			        yhitpos=emhits_y_array[ihit];
+			        hit_index=ihit;
+		}
+	}
+	//sum over all 3x3 towers around the leading tower
+	double xClus=xhitpos*maxHitEnergy;
+	double yClus=yhitpos*maxHitEnergy;
+	for(int ihit=0;ihit<emhits_energy_array.GetSize();ihit++){
+		double hitenergy=emhits_energy_array[ihit];
+		double x=emhits_x_array[ihit];
+		double y=emhits_y_array[ihit];
+		double d=sqrt( (x-xhitpos)*(x-xhitpos) + (y-yhitpos)*(y-yhitpos));
+		if(d<70. && ihit!=hit_index && hitenergy>0.01)  {
+			maxHitEnergy+=hitenergy;//clustering around leading tower 3 crystal = 60mm.
+			xClus+=x*hitenergy;
+			yClus+=y*hitenergy;
+		}
+	}
+	
+	h_ClusOverHit_REC->Fill( maxEnergy / maxHitEnergy );
+	//weighted average cluster position.
+	xClus = xClus/maxHitEnergy;
+	yClus = yClus/maxHitEnergy;
+	double radius=sqrt(xClus*xClus+yClus*yClus);
+	if( radius>550. ) continue; //geometric acceptance cut
+	//4.4% energy calibration.
+	double clusEnergy=1.044*maxHitEnergy; 
+
+	h_energy_REC->Fill(clusEnergy);
+	//ratio of reco / truth Energy
+	h_energy_calibration_REC->Fill( clusEnergy / scatMC.E() );
+	//energy resolution
+	double res= (scatMC.E()-clusEnergy)/scatMC.E();
+	h_energy_res->Fill(scatMC.E(), res);
+	h_emClus_position_REC->Fill(xpos,ypos);//default clustering position
+	h_emHits_position_REC->Fill(xClus,yClus); //self clustering position
+	
+	//association of rec level scat' e
+	int rec_elect_index=-1;
+	for(int i=0;i<sim_id.GetSize();i++){
+		if(sim_id[i]==mc_elect_index){
+			//find the rec_id
+			rec_elect_index = rec_id[i];
+		}
+	}
+    
+    TLorentzVector scatMCmatchREC(0,0,0,0);
+    TLorentzVector scatREC(0,0,0,0);
+    TLorentzVector scatClusEREC(0,0,0,0);
+    TLorentzVector hfs(0,0,0,0);
+    TLorentzVector particle(0,0,0,0);
+    TLorentzVector kplusREC(0,0,0,0);
+    TLorentzVector kminusREC(0,0,0,0);
+    TLorentzVector vmREC(0,0,0,0);
+
+    double maxP=-1.;
+    //track loop
+	for(int itrk=0;itrk<reco_pz_array.GetSize();itrk++){
+		TVector3 trk(reco_px_array[itrk],reco_py_array[itrk],reco_pz_array[itrk]);
+		if(rec_elect_index!=-1
+			&& itrk==rec_elect_index){
+			scatMCmatchREC.SetVectM(trk,MASS_ELECTRON);//Reserved to calculate t.
+		}
+		if(trk.Mag()>maxP){
+			//track-base 4 vector
+			maxP=trk.Mag();
+			scatREC.SetVectM(trk,MASS_ELECTRON);
+
+			//use emcal energy to define 4 vector
+			double p = sqrt(clusEnergy*clusEnergy- MASS_ELECTRON*MASS_ELECTRON );
+			double eta=scatREC.Eta();
+			double phi=scatREC.Phi();
+			double pt = TMath::Sin(scatREC.Theta())*p;
+			scatClusEREC.SetPtEtaPhiM(pt,eta,phi,MASS_ELECTRON);
+		}
+	}
+	//loop over track again;
+	for(int itrk=0;itrk<reco_pz_array.GetSize();itrk++){
+		TVector3 trk(reco_px_array[itrk],reco_py_array[itrk],reco_pz_array[itrk]);
+		particle.SetVectM(trk,MASS_PION);//assume pions;
+		if(itrk!=rec_elect_index) {
+    		hfs += particle; //hfs 4vector sum.
+    		//selecting phi->kk daughters;
+    		h_eta->Fill(trk.Eta());
+    		if(fabs(trk.Eta())<3.0){
+    			if(reco_charge_array[itrk]>0) kplusREC.SetVectM(trk,MASS_KAON);
+    			if(reco_charge_array[itrk]<0) kminusREC.SetVectM(trk,MASS_KAON);
+    		}
+		}
+	}
+	//4vector of VM;
+	if(kplusREC.E()!=0. && kminusREC.E()!=0.){
+		vmREC=kplusREC+kminusREC;
+	}
+
+	//track-base e' energy REC;
+	h_trk_energy_REC->Fill(scatMCmatchREC.E());
+	
+	//track-base e' energy resolution;
+	res= (scatMC.E()-scatMCmatchREC.E())/scatMC.E();
+	h_trk_energy_res->Fill(scatMC.E(), res);
+	
+	//track-base e' pt resolution;
+	res= (scatMC.Pt()-scatMCmatchREC.Pt())/scatMC.Pt();
+	h_trk_Pt_res->Fill(scatMC.Pt(), res);
+
+	//track-base Epz scat' e
+	double EpzREC= (scatMCmatchREC+hfs).E() - (scatMCmatchREC+hfs).Pz();
+	h_trk_Epz_REC->Fill( EpzREC );
+
+	//EEMC cluster Epz scat' e
+	EpzREC= (scatClusEREC+hfs).E() - (scatClusEREC+hfs).Pz();
+	h_Epz_REC->Fill( EpzREC );
+
+	//E over p
+	double EoverP=scatClusEREC.E() / scatMCmatchREC.P();
+	h_EoverP_REC->Fill( EoverP );
+
+	//cluster-base DIS kine;
+	TLorentzVector qbeamREC=ebeam-scatClusEREC;
+	double Q2REC=-(qbeamREC).Mag2();  
+	double pqREC=pbeam.Dot(qbeamREC);
+	double yREC= pqREC/pbeam.Dot(ebeam);
+	h_Q2REC_e->Fill(Q2REC);
+	h_yREC_e->Fill(yREC);
+
+	//Event selection:
+	if( EpzREC<27||EpzREC>40 ) continue;
+	if( EoverP<0.8||EoverP>1.18 ) continue;		
+
+	//MC level phase space cut
+	if(Q2REC<1.||Q2REC>10.) continue;
+	if(yREC<0.01||yREC>0.85) continue;
+
+	//VM rec
+	if(vmREC.E()==0) continue;
+	double phi_mass = vmREC.M();
+	h_VM_mass_REC->Fill(phi_mass);
+	h_VM_pt_REC->Fill(vmREC.Pt());
+
+	//select phi mass and rapidity window 
+	if( fabs(phi_mass-1.02)<0.02
+    		&& fabs(vmREC.Rapidity())<3.5 ){
+    	//2 versions: track and energy cluster:
+	double t_trk_REC = giveme_t_method_L(ebeam,scatMCmatchREC,pbeam,vmREC);
+    	double t_REC = giveme_t_method_L(ebeam,scatClusEREC,pbeam,vmREC);
+    	h_t_trk_REC->Fill( t_trk_REC );
+    	h_t_REC->Fill( t_REC );
+    	h_t_REC_2D->Fill(t_trk_REC,t_REC);
+    	if( (t_trk_REC/t_REC) > 0.5 && (t_trk_REC/t_REC) < 1.5 ){
+    		h_t_combo_REC->Fill( (t_trk_REC+t_REC)/2. );//w=1./(fabs(1.0-(t_trk_REC/t_REC)))
+    	} 
+    	h_t_RECMC_2D->Fill(t_MC,t_trk_REC/t_REC);
+
+    	//t track resolution 
+	res= (t_MC-t_trk_REC)/t_MC;
+	h_trk_t_res->Fill(t_MC, res);
+	
+    	//t EEMC resolution;
+	res= (t_MC-t_REC)/t_MC;
+	h_t_res->Fill(t_MC, res);
+	
+	//2D t
+	h_t_2D->Fill(t_MC,t_trk_REC);
+
+    	//VM pt resolution;
+    	res= (vmMC.Pt()-vmREC.Pt())/vmMC.Pt();
+	h_VM_res->Fill(vmMC.Pt(), res);
+    	}
+
+}
+output->Write();
+output->Close();
+
+return 0;
+}

benchmarks/diffractive_vmp/config.yml

+diffractive_vmp:compile:
+  stage: compile
+  extends: .compile_benchmark
+  script:
+    - compile_analyses.py diffractive_vmp
+
+diffractive_vmp:generate:
+  stage: generate
+  extends: .phy_benchmark
+  needs: ["common:detector", "diffractive_vmp:compile"]
+  parallel:
+    matrix:
+      - VM: jpsi
+  timeout: 1 hours
+  script:
+    - bash benchmarks/diffractive_vmp/get.sh --config diffractive_${vm}
+
+diffractive_vmp:simulate:
+  stage: simulate
+  extends: .phy_benchmark
+  needs: ["diffractive_vmp:generate"]
+  parallel:
+    matrix:
+      - VM: jpsi
+  timeout: 96 hour
+  script:
+    - bash benchmarks/diffractive_vmp/diffractive_vmp.sh --config diffractive_${vm}
+  retry:
+    max: 2
+    when:
+      - runner_system_failure
+
+diffractive_vmp:results:
+  stage: collect
+  needs: ["diffractive_vmp:simulate"]
+  script:
+    - collect_tests.py diffractive_vmp

benchmarks/diffractive_vmp/diffractive_vmp.sh

+#!/bin/bash
+source strict-mode.sh
+
+## =============================================================================
+## Run the Diffractive VMP benchmarks in 5 steps:
+## 1. Parse the command line and setup environment
+## 2. Detector simulation through ddsim
+## 3. Digitization and reconstruction through Juggler
+## 4. Root-based Physics analyses
+## 5. Finalize
+## =============================================================================
+
+## make sure we launch this script from the project root directory
+PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/../..
+pushd ${PROJECT_ROOT}
+
+echo "Running the Diffractive VMP benchmarks"
+
+
+## =============================================================================
+## Step 1: Setup the environment variables
+##
+## First parse the command line flags.
+## This sets the following environment variables:
+## - CONFIG:   The specific generator configuration
+## - EBEAM:    The electron beam energy
+## - PBEAM:    The ion beam energy
+## - LEADING:  Leading particle of interest (J/psi)
+#export REQUIRE_LEADING=1
+source parse_cmd.sh $@
+
+## We also need the following benchmark-specific variables:
+##
+## - BENCHMARK_TAG: Unique identified for this benchmark process.
+## - BEAM_TAG:      Identifier for the chosen beam configuration
+## - INPUT_PATH:    Path for generator-level input to the benchmarks
+## - TMP_PATH:      Path for temporary data (not exported as artifacts)
+## - RESULTS_PATH:  Path for benchmark output figures and files
+##
+## You can read dvmp/env.sh for more in-depth explanations of the variables.
+source benchmarks/diffractive_vmp/env.sh
+
+## Get a unique file names based on the configuration options
+GEN_FILE=${INPUT_PATH}/gen-${CONFIG}_${JUGGLER_N_EVENTS}.hepmc
+
+SIM_FILE=${TMP_PATH}/sim-${CONFIG}.edm4hep.root
+SIM_LOG=${TMP_PATH}/sim-${CONFIG}.log
+
+
+REC_FILE=${TMP_PATH}/rec-${CONFIG}.root
+REC_LOG=${TMP_PATH}/sim-${CONFIG}.log
+
+PLOT_TAG=${CONFIG}
+
+## =============================================================================
+## Step 2: Run the simulation
+echo "Running Geant4 simulation"
+ls -lrth 
+ls -lrth input
+echo ${TMP_PATH}
+ls -lrth ${TMP_PATH}
+ddsim --runType batch \
+      --part.minimalKineticEnergy 100*GeV  \
+      --filter.tracker edep0 \
+      -v WARNING \
+      --numberOfEvents ${JUGGLER_N_EVENTS} \
+      --compactFile ${DETECTOR_PATH}/${DETECTOR_CONFIG}.xml \
+      --inputFiles ${GEN_FILE} \
+      --outputFile ${SIM_FILE}
+if [ "$?" -ne "0" ] ; then
+  echo "ERROR running ddsim"
+  exit 1
+fi
+
+## =============================================================================
+## Step 3: Run digitization & reconstruction
+echo "Running the digitization and reconstruction"
+## FIXME Need to figure out how to pass file name to juggler from the commandline
+## the tracker_reconstruction.py options file uses the following environment
+## variables:
+## - JUGGLER_SIM_FILE:    input detector simulation
+## - JUGGLER_REC_FILE:    output reconstructed data
+## - JUGGLER_N_EVENTS:    number of events to process (part of global environment)
+## - DETECTOR:    detector package (part of global environment)
+export JUGGLER_SIM_FILE=${SIM_FILE}
+export JUGGLER_REC_FILE=${REC_FILE}
+if [ ${RECO} == "eicrecon" ] ; then
+  eicrecon ${JUGGLER_SIM_FILE} -Ppodio:output_file=${JUGGLER_REC_FILE}
+  if [[ "$?" -ne "0" ]] ; then
+    echo "ERROR running eicrecon"
+    exit 1
+  fi
+fi
+
+if [[ ${RECO} == "juggler" ]] ; then
+  gaudirun.py options/reconstruction.py || [ $? -eq 4 ]
+  if [ "$?" -ne "0" ] ; then
+    echo "ERROR running juggler"
+    exit 1
+  fi
+fi
+## =============================================================================
+## Step 4: Analysis
+
+## write a temporary configuration file for the analysis script
+CONFIG="${TMP_PATH}/${PLOT_TAG}.json"
+cat << EOF > ${CONFIG}
+{
+  "rec_file": "${REC_FILE}",
+  "vm_name": "${LEADING}",
+  "detector": "${DETECTOR_CONFIG}",
+  "output_prefix": "${RESULTS_PATH}/${PLOT_TAG}",
+  "test_tag": "${LEADING}_${BEAM_TAG}"
+}
+EOF
+#cat ${CONFIG}
+
+## run the analysis script with this configuration
+root -b -q "benchmarks/exclusive/diffractive_vmp/analysis/simple_analysis.cxx+(\"${CONFIG}\")"
+if [ "$?" -ne "0" ] ; then
+  echo "ERROR running vm_mass script"
+  exit 1
+fi
+
+## =============================================================================
+## Step 5: finalize
+echo "Finalizing Diffractive VMP benchmark"
+
+## Copy over reconsturction artifacts as long as we don't have
+## too many events
+if [ "${JUGGLER_N_EVENTS}" -lt "500" ] ; then 
+  cp ${REC_FILE} ${RESULTS_PATH}
+fi
+
+## Always move over log files to the results path
+mv ${REC_LOG} ${RESULTS_PATH}
+
+
+## cleanup output files
+#rm -f ${REC_FILE} ${SIM_FILE} ## --> not needed for CI
+
+## =============================================================================
+## All done!
+echo "Diffractive VMP benchmarks complete"

benchmarks/diffractive_vmp/env.sh

+#!/bin/bash
+source strict-mode.sh
+
+## =============================================================================
+## Local configuration variables for this particular benchmark 
+## It defines the following additional variables: 
+##
+##  - BENCHMARK_TAG:          Tag to identify this particular benchmark
+##  - BEAM_TAG                Tag to identify the beam configuration
+##  - INPUT_PATH:             Path for generator-level input to the benchmarks
+##  - TMP_PATH:               Path for temporary data (not exported as artifacts)
+##  - RESULTS_PATH:           Path for benchmark output figures and files
+##
+## This script assumes that EBEAM and PBEAM are set as part of the
+## calling script (usually as command line argument).
+##
+## =============================================================================
+
+## Tag for the local benchmark. Should be the same as the directory name for
+## this particular benchmark set (for clarity). 
+## This tag is used for the output artifacts directory (results/${JUGGLER_TAG}) 
+## and a tag in some of the output files.
+export BENCHMARK_TAG="diffractive_vmp"
+echo "Setting up the local environment for the ${BENCHMARK_TAG^^} benchmarks"
+
+## Extra beam tag to identify the desired beam configuration
+export BEAM_TAG="${EBEAM}on${PBEAM}"
+
+if [[ ! -d "input" ]] ; then
+  echo " making local link to input "
+  mkdir_local_data_link  input
+fi
+
+## Data path for input data (generator-level hepmc file)
+INPUT_PATH="input/${BENCHMARK_TAG}/${BEAM_TAG}"
+#export INPUT_PATH=`realpath ${INPUT_PATH}`
+mkdir -p ${INPUT_PATH}
+echo "INPUT_PATH:             ${INPUT_PATH}"
+
+
+## Data path for temporary data (not exported as artifacts)
+TMP_PATH=${LOCAL_DATA_PATH}/tmp/${BEAM_TAG}
+#export TMP_PATH=`realpath ${TMP_PATH}`
+mkdir -p ${TMP_PATH}
+echo "TMP_PATH:               ${TMP_PATH}"
+
+## Data path for benchmark output (plots and reconstructed files
+## if not too big).
+RESULTS_PATH="results/${BENCHMARK_TAG}/${BEAM_TAG}"
+mkdir -p ${RESULTS_PATH}
+export RESULTS_PATH=`realpath ${RESULTS_PATH}`
+echo "RESULTS_PATH:           ${RESULTS_PATH}"
+
+## =============================================================================
+## That's all!
+echo "Local environment setup complete."

benchmarks/diffractive_vmp/get.sh

+#!/bin/bash
+source strict-mode.sh
+
+## =============================================================================
+## Standin for a proper pythia generation process, similar to how we
+## generate events for DVMP
+## Runs in 5 steps:
+##   1. Parse the command line and setup the environment
+##   2. Check if we can download the file
+##   3. Finalize
+## =============================================================================
+## =============================================================================
+
+## make sure we launch this script from the project root directory
+PROJECT_ROOT="$( cd "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"/../..
+pushd ${PROJECT_ROOT}
+
+## =============================================================================
+## Step 1: Setup the environment variables
+## First parse the command line flags.
+## This sets the following environment variables:
+## - CONFIG:   The specific generator configuration --> not currenlty used FIXME
+## - EBEAM:    The electron beam energy --> not currently used FIXME
+## - PBEAM:    The ion beam energy --> not currently used FIXME
+source parse_cmd.sh $@
+
+## To run the generator, we need the following global variables:
+##
+## - LOCAL_PREFIX:      Place to cache local packages and data
+## - JUGGLER_N_EVENTS:  Number of events to process
+## - JUGGLER_RNG_SEED:  Random seed for event generation.
+##
+## defined in common_bench repo
+## You can ready bin/env.sh for more in-depth explanations of the variables
+## and how they can be controlled.
+
+
+## We also need the following benchmark-specific variables:
+##
+## - BENCHMARK_TAG: Unique identified for this benchmark process.
+## - INPUT_PATH:    Path for generator-level input to the benchmarks
+## - TMP_PATH:      Path for temporary data (not exported as artifacts)
+##
+## You can read dvmp/env.sh for more in-depth explanations of the variables.
+source benchmarks/diffractive_vmp/env.sh
+
+## Get a unique file name prefix based on the configuration options
+GEN_TAG=gen-${CONFIG}_${JUGGLER_N_EVENTS} ## Generic file prefix
+
+## =============================================================================
+## Step 2: Check if we can find the file
+if [ -f "${INPUT_PATH}/${GEN_TAG}.hepmc" ]; then
+  echo "Found cached generator output for $GEN_TAG, no need to rerun"
+  exit 0
+fi
+
+## =============================================================================
+## Step 3: Copy the file (about 180 lines per event in DIS NC files)
+nlines=$((190*${JUGGLER_N_EVENTS}))
+DATA_URL=S3/eictest/EPIC/EVGEN/EXCLUSIVE/DIFFRACTIVE_JPSI_ABCONV/Sartre/Coherent/sartre_bnonsat_Au_jpsi_ab_eAu_1_000.hepmc.gz
+mc config host add S3 https://dtn01.sdcc.bnl.gov:9000 ${S3_ACCESS_KEY} ${S3_SECRET_KEY}
+mc head -n ${nlines} ${DATA_URL} | sanitize_hepmc3 > ${TMP_PATH}/${GEN_TAG}.hepmc
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR downloading file"
+  exit 1
+fi
+
+## =============================================================================
+## Step 4: Finally, move relevant output into the artifacts directory and clean up
+## =============================================================================
+echo "Moving generator output to ${INPUT_PATH}/${GEN_TAG}.hepmc"
+mv ${TMP_PATH}/${GEN_TAG}.hepmc ${INPUT_PATH}/${GEN_TAG}.hepmc
+## this step only matters for local execution
+echo "Cleaning up"
+## does nothing
+
+## =============================================================================
+## All done!
+echo "$BENCHMARK_TAG event generation complete"