diff --git a/dis/config.yml b/dis/config.yml index ac6797a17ffcaef8df24f028a2966e30cedb351f..a300e5f034632f6d966711e166770c19696d9f07 100644 --- a/dis/config.yml +++ b/dis/config.yml @@ -1,24 +1,11 @@ -dis:dummy_test: +dis:run_test: stage: analyze timeout: 1 hours script: - - dis/dummy_test.sh - -dis:dummy_test2: - stage: analyze - timeout: 1 hours - script: - - dis/dummy_test2.sh - -dis:dummy_fail_test: - stage: analyze - timeout: 1 hours - allow_failure: true - script: - - dis/dummy_fail_test.sh + - bash dis/dis.sh dis:results: stage: collect - needs: ["dis:dummy_test", "dis:dummy_test2", "dis:dummy_fail_test"] + needs: ["dis:run_test"] script: - echo "All DIS benchmarks successful" diff --git a/dis/dis.sh b/dis/dis.sh new file mode 100644 index 0000000000000000000000000000000000000000..657c9b80c04d5c169a1adb86748533d783aa8ed7 --- /dev/null +++ b/dis/dis.sh @@ -0,0 +1,83 @@ +#!/bin/bash + +if [[ ! -n "${JUGGLER_DETECTOR}" ]] ; then + export JUGGLER_DETECTOR="topside" +fi + +if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then + export JUGGLER_N_EVENTS=100 +fi + +if [[ ! -n "${JUGGLER_INSTALL_PREFIX}" ]] ; then + export JUGGLER_INSTALL_PREFIX="/usr/local" +fi + +export JUGGLER_FILE_NAME_TAG="dis" +export JUGGLER_GEN_FILE="${JUGGLER_FILE_NAME_TAG}.hepmc" + +export JUGGLER_SIM_FILE="sim_${JUGGLER_FILE_NAME_TAG}.root" +export JUGGLER_REC_FILE="rec_${JUGGLER_FILE_NAME_TAG}.root" + +echo "JUGGLER_N_EVENTS = ${JUGGLER_N_EVENTS}" +echo "JUGGLER_DETECTOR = ${JUGGLER_DETECTOR}" + + +### Build the detector constructors. +git clone https://eicweb.phy.anl.gov/EIC/detectors/${JUGGLER_DETECTOR}.git +git clone https://eicweb.phy.anl.gov/EIC/detectors/accelerator.git +pushd ${JUGGLER_DETECTOR} +ln -s ../accelerator/eic +popd +mkdir ${JUGGLER_DETECTOR}/build +pushd ${JUGGLER_DETECTOR}/build +cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j30 install +popd + +# generate the input events +# temporary standin until hepmc output from pythia is generated. +root -b -q "dis/scripts/gen_central_electrons.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")" +if [[ "$?" -ne "0" ]] ; then + echo "ERROR running script" + exit 1 +fi + +# +pushd ${JUGGLER_DETECTOR} + +## run geant4 simulations +npsim --runType batch \ + --part.minimalKineticEnergy 1000*GeV \ + -v WARNING \ + --numberOfEvents ${JUGGLER_N_EVENTS} \ + --compactFile ${JUGGLER_DETECTOR}.xml \ + --inputFiles ../${JUGGLER_FILE_NAME_TAG}.hepmc \ + --outputFile ${JUGGLER_SIM_FILE} +if [[ "$?" -ne "0" ]] ; then + echo "ERROR running script" + exit 1 +fi + +# Need to figure out how to pass file name to juggler from the commandline +xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \ + gaudirun.py ../options/tracker_reconstruction.py + +if [[ "$?" -ne "0" ]] ; then + echo "ERROR running juggler" + exit 1 +fi +ls -l +popd + +pwd +mkdir -p results/dis + +root -b -q "dis/scripts/rec_dis_electrons.cxx(\"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\")" +if [[ "$?" -ne "0" ]] ; then + echo "ERROR running root script" + exit 1 +fi + +if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then +cp ${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE} results/. +fi + diff --git a/dis/dummy_fail_test.sh b/dis/dummy_fail_test.sh deleted file mode 100755 index f8c506b3ea3d71d5acb606b58de5e31bbac6b62d..0000000000000000000000000000000000000000 --- a/dis/dummy_fail_test.sh +++ /dev/null @@ -1,7 +0,0 @@ -#!/bin/bash - -echo "Dummy Test..." -echo "..." -echo "Fails!" - -exit 1 diff --git a/dis/dummy_test.sh b/dis/dummy_test.sh deleted file mode 100755 index 5af1dc34f2f597a6327a3a5ef7d8ca734e413c1b..0000000000000000000000000000000000000000 --- a/dis/dummy_test.sh +++ /dev/null @@ -1,7 +0,0 @@ -#!/bin/bash - -echo "Dummy Test..." -echo "..." -echo "Passes!" - -#exit 1 diff --git a/dis/dummy_test2.sh b/dis/dummy_test2.sh deleted file mode 100755 index dfede272a5390cb72e5417698f6390b0e7c55349..0000000000000000000000000000000000000000 --- a/dis/dummy_test2.sh +++ /dev/null @@ -1,7 +0,0 @@ -#!/bin/bash - -echo "Dummy Test number 2..." -echo "..." -echo "Passes!" - -#exit 1 diff --git a/dis/scripts/gen_central_electrons.cxx b/dis/scripts/gen_central_electrons.cxx new file mode 100644 index 0000000000000000000000000000000000000000..5b06ada728ef1ee90d12a993fe5edf730f05a304 --- /dev/null +++ b/dis/scripts/gen_central_electrons.cxx @@ -0,0 +1,83 @@ +#include "HepMC3/GenEvent.h" +#include "HepMC3/ReaderAscii.h" +#include "HepMC3/WriterAscii.h" +#include "HepMC3/Print.h" + +#include <iostream> +#include<random> +#include<cmath> +#include <math.h> +#include <TMath.h> + +using namespace HepMC3; + +/** Generate electrons in the central region. + * This is for testing detectors in the "barrel" region. + */ +void gen_central_electrons(int n_events = 100, + const char* out_fname = "central_electrons.hepmc") +{ + double cos_theta_min = std::cos( 50.0*(M_PI/180.0)); + double cos_theta_max = std::cos(130.0*(M_PI/180.0)); + + WriterAscii hepmc_output(out_fname); + int events_parsed = 0; + GenEvent evt(Units::GEV, Units::MM); + + // Random number generator + TRandom *r1 = new TRandom(); + + for (events_parsed = 0; events_parsed < n_events; events_parsed++) { + // FourVector(px,py,pz,e,pdgid,status) + // type 4 is beam + // pdgid 11 - electron + // pdgid 111 - pi0 + // pdgid 2212 - proton + GenParticlePtr p1 = + std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4); + GenParticlePtr p2 = std::make_shared<GenParticle>( + FourVector(0.0, 0.0, 0.0, 0.938), 2212, 4); + + // Define momentum + Double_t p = r1->Uniform(1.0, 10.0); + Double_t phi = r1->Uniform(0.0, 2.0 * M_PI); + Double_t costh = r1->Uniform(cos_theta_min, cos_theta_max); + Double_t th = std::acos(costh); + Double_t px = p * std::cos(phi) * std::sin(th); + Double_t py = p * std::sin(phi) * std::sin(th); + Double_t pz = p * std::cos(th); + // Generates random vectors, uniformly distributed over the surface of a + // sphere of given radius, in this case momentum. + // r1->Sphere(px, py, pz, p); + + //std::cout << std::sqrt(px*px + py*py + pz*pz) - p << " is zero? \n"; + + // type 1 is final state + // pdgid 11 - electron 0.510 MeV/c^2 + GenParticlePtr p3 = std::make_shared<GenParticle>( + FourVector( + px, py, pz, + sqrt(p*p + (0.000511 * 0.000511))), + 11, 1); + + GenVertexPtr v1 = std::make_shared<GenVertex>(); + v1->add_particle_in(p1); + v1->add_particle_in(p2); + + v1->add_particle_out(p3); + evt.add_vertex(v1); + + if (events_parsed == 0) { + std::cout << "First event: " << std::endl; + Print::listing(evt); + } + + hepmc_output.write_event(evt); + if (events_parsed % 10000 == 0) { + std::cout << "Event: " << events_parsed << std::endl; + } + evt.clear(); + } + hepmc_output.close(); + std::cout << "Events parsed and written: " << events_parsed << std::endl; +} diff --git a/dis/scripts/rec_dis_electrons.cxx b/dis/scripts/rec_dis_electrons.cxx new file mode 100644 index 0000000000000000000000000000000000000000..b7227a1ddb33c2564fd160e2124da774a9711161 --- /dev/null +++ b/dis/scripts/rec_dis_electrons.cxx @@ -0,0 +1,113 @@ +#include "ROOT/RDataFrame.hxx" +#include <iostream> + +#include "dd4pod/Geant4ParticleCollection.h" +#include "eicd/TrackParametersCollection.h" +#include "eicd/ClusterCollection.h" +#include "eicd/ClusterData.h" + +using ROOT::RDataFrame; +using namespace ROOT::VecOps; + +auto p_track = [](std::vector<eic::TrackParametersData> const& in) { + std::vector<double> result; + for (size_t i = 0; i < in.size(); ++i) { + result.push_back(std::abs(1.0/(in[i].qOverP))); + } + return result; +}; + +std::vector<float> pt (std::vector<dd4pod::Geant4ParticleData> const& in){ + std::vector<float> result; + for (size_t i = 0; i < in.size(); ++i) { + result.push_back(std::sqrt(in[i].psx * in[i].psx + in[i].psy * in[i].psy)); + } + return result; +} + +auto momentum = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) { + std::vector<double> result; + for (size_t i = 0; i < in.size(); ++i) { + result.push_back(in[i].E()); + } + return result; +}; +auto theta = [](std::vector<ROOT::Math::PxPyPzMVector> const& in) { + std::vector<double> result; + for (size_t i = 0; i < in.size(); ++i) { + result.push_back(in[i].Theta()*180/M_PI); + } + return result; +}; +auto fourvec = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) { + std::vector<ROOT::Math::PxPyPzMVector> result; + ROOT::Math::PxPyPzMVector lv; + for (size_t i = 0; i < in.size(); ++i) { + lv.SetCoordinates(in[i].psx, in[i].psy, in[i].psz, in[i].mass); + result.push_back(lv); + } + return result; +}; + +auto delta_p = [](const std::vector<double>& tracks, const std::vector<double>& thrown) { + std::vector<double> res; + for (const auto& p1 : thrown) { + for (const auto& p2 : tracks) { + res.push_back(p1 - p2); + } + } + return res; +}; + +int rec_dis_electrons(const char* fname = "topside/rec_central_electrons.root") +{ + + ROOT::EnableImplicitMT(); + ROOT::RDataFrame df("events", fname); + + auto df0 = df.Define("isThrown", "mcparticles2.genStatus == 1") + .Define("thrownParticles", "mcparticles2[isThrown]") + .Define("thrownP", fourvec, {"thrownParticles"}) + .Define("p_thrown", momentum, {"thrownP"}) + .Define("nTracks", "outputTrackParameters.size()") + .Define("p_track", p_track, {"outputTrackParameters"}) + .Define("p_track1", p_track, {"outputTrackParameters1"}) + .Define("p_track2", p_track, {"outputTrackParameters2"}) + .Define("delta_p",delta_p, {"p_track", "p_thrown"}) + .Define("delta_p1",delta_p, {"p_track1", "p_thrown"}) + .Define("delta_p2",delta_p, {"p_track2", "p_thrown"}); + + auto h_nTracks = df0.Histo1D({"h_nTracks", "; N tracks ", 10, 0, 10}, "nTracks"); + auto h_pTracks = df0.Histo1D({"h_pTracks", "; GeV/c ", 100, 0, 10}, "p_track"); + + auto h_delta_p = df0.Histo1D({"h_delta_p", "; GeV/c ", 100, -10, 10}, "delta_p"); + auto h_delta_p1 = df0.Histo1D({"h_delta_p1", "; GeV/c ", 100, -10, 10}, "delta_p1"); + auto h_delta_p2 = df0.Histo1D({"h_delta_p2", "; GeV/c ", 100, -10, 10}, "delta_p2"); + + auto c = new TCanvas(); + + h_nTracks->DrawCopy(); + c->SaveAs("results/dis/rec_central_electrons_nTracks.png"); + c->SaveAs("results/dis/rec_central_electrons_nTracks.pdf"); + + h_pTracks->DrawCopy(); + c->SaveAs("results/dis/rec_central_electrons_pTracks.png"); + c->SaveAs("results/dis/rec_central_electrons_pTracks.pdf"); + + THStack * hs = new THStack("hs_delta_p","; GeV/c "); + TH1D* h1 = (TH1D*) h_delta_p->Clone(); + hs->Add(h1); + h1 = (TH1D*) h_delta_p1->Clone(); + h1->SetLineColor(2); + hs->Add(h1); + h1 = (TH1D*) h_delta_p2->Clone(); + h1->SetLineColor(4); + h1->SetFillStyle(3001); + h1->SetFillColor(4); + hs->Add(h1); + hs->Draw("nostack"); + c->SaveAs("results/dis/rec_central_electrons_delta_p.png"); + c->SaveAs("results/dis/rec_central_electrons_delta_p.pdf"); + + return 0; +}