From cf728adf66ebe85a36bc21675ea367eda5fb30ca Mon Sep 17 00:00:00 2001
From: Whitney Armstrong <warmstrong@anl.gov>
Date: Thu, 22 Oct 2020 10:35:30 -0500
Subject: [PATCH] Adding clustering benchmarks.
- Added barrel clustering test.
- Running with existing clustering algorithms to unknown efficacy.
---
.gitlab-ci.yml | 1 +
clustering/barrel_clusters.sh | 56 +++++++++++++
clustering/clustering_config.yml | 14 ++++
clustering/options/calorimeter_clustering.py | 77 ++++++++++++++++++
clustering/scripts/gen_central_electrons.cxx | 83 ++++++++++++++++++++
5 files changed, 231 insertions(+)
create mode 100644 clustering/barrel_clusters.sh
create mode 100644 clustering/clustering_config.yml
create mode 100644 clustering/options/calorimeter_clustering.py
create mode 100644 clustering/scripts/gen_central_electrons.cxx
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index d2d2f164..d0396a27 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -27,6 +27,7 @@ get_data:
include:
- local: 'ecal/ecal_config.yml'
- local: 'tracking/tracking_config.yml'
+ - local: 'clustering/clustering_config.yml'
final_report:
diff --git a/clustering/barrel_clusters.sh b/clustering/barrel_clusters.sh
new file mode 100644
index 00000000..c04b3ae3
--- /dev/null
+++ b/clustering/barrel_clusters.sh
@@ -0,0 +1,56 @@
+#!/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="barrel_clusters"
+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
+mkdir ${JUGGLER_DETECTOR}/build
+pushd ${JUGGLER_DETECTOR}/build
+cmake ../. -DCMAKE_INSTALL_PREFIX=/usr/local && make -j30 install
+popd
+#
+## generate the input events
+## note datasets is now only used to develop datasets.
+##git clone https://eicweb.phy.anl.gov/EIC/datasets.git datasets
+root -b -q "clustering/scripts/gen_central_electrons.cxx(${JUGGLER_N_EVENTS}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+#
+pushd ${JUGGLER_DETECTOR}
+#ls -l
+## 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}
+
+# Need to figure out how to pass file name to juggler from the commandline
+xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv gaudirun.py ../clustering/options/calorimeter_clustering.py
+ls -l
+popd
+
+pwd
+mkdir -p results
+cp topside/${JUGGLER_REC_FILE} results/.
+
diff --git a/clustering/clustering_config.yml b/clustering/clustering_config.yml
new file mode 100644
index 00000000..89bfc9ca
--- /dev/null
+++ b/clustering/clustering_config.yml
@@ -0,0 +1,14 @@
+clustering_barrel_electrons:
+ image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:latest
+ tags:
+ - silicon
+ needs: ["get_data"]
+ timeout: 24 hours
+ artifacts:
+ expire_in: 20 weeks
+ paths:
+ - results/
+ stage: run
+ script:
+ - bash clustering/barrel_clusters.sh
+
diff --git a/clustering/options/calorimeter_clustering.py b/clustering/options/calorimeter_clustering.py
new file mode 100644
index 00000000..1dce32dd
--- /dev/null
+++ b/clustering/options/calorimeter_clustering.py
@@ -0,0 +1,77 @@
+from Gaudi.Configuration import *
+import os
+
+from GaudiKernel.DataObjectHandleBase import DataObjectHandleBase
+from Configurables import ApplicationMgr, EICDataSvc, PodioOutput, GeoSvc
+from GaudiKernel import SystemOfUnits as units
+
+detector_name = "topside"
+if "JUGGLER_DETECTOR" in os.environ :
+ detector_name = str(os.environ["JUGGLER_DETECTOR"])
+
+# todo add checks
+input_sim_file = str(os.environ["JUGGLER_SIM_FILE"])
+output_rec_file = str(os.environ["JUGGLER_REC_FILE"])
+n_events = str(os.environ["JUGGLER_N_EVENTS"])
+
+geo_service = GeoSvc("GeoSvc", detectors=["{}.xml".format(detector_name)])
+podioevent = EICDataSvc("EventDataSvc", inputs=[input_sim_file])
+
+from Configurables import PodioInput
+from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollection_dd4pod__Geant4ParticleCollection_ as MCCopier
+from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
+
+from Configurables import Jug__Digi__HadronicCalDigi as HadCalorimeterDigi
+from Configurables import Jug__Digi__EMCalorimeterDigi as EMCalorimeterDigi
+from Configurables import Jug__Digi__CrystalEndcapsDigi as CrystalEndcapsDigi
+
+from Configurables import Jug__Reco__CrystalEndcapsReco as CrystalEndcapsReco
+from Configurables import Jug__Reco__EMCalReconstruction as EMCalReconstruction
+
+from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
+from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
+
+podioinput = PodioInput("PodioReader", collections=["mcparticles","CrystalEcalHits","EcalBarrelHits"], OutputLevel=DEBUG)
+
+## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
+copier = MCCopier("MCCopier", inputCollection="mcparticles", outputCollection="mcparticles2",OutputLevel=DEBUG)
+calcopier = CalCopier("CalCopier", inputCollection="CrystalEcalHits", outputCollection="CrystalEcalHits2",OutputLevel=DEBUG)
+
+emcaldigi = CrystalEndcapsDigi("ecal_digi", inputHitCollection="CrystalEcalHits", outputHitCollection="RawDigiEcalHits")
+ecdigi = EMCalorimeterDigi("ec_barrel_digi", inputHitCollection="EcalBarrelHits", outputHitCollection="RawEcalBarrelHits")
+
+crystal_ec_reco = CrystalEndcapsReco("crystal_ec_reco", inputHitCollection="RawDigiEcalHits", outputHitCollection="RecoEcalHits",
+ minModuleEdep=1.0*units.MeV,OutputLevel=DEBUG)
+ecal_reco = EMCalReconstruction("ecal_reco", inputHitCollection="RawEcalBarrelHits", outputHitCollection="RecEcalBarrelHits",
+ minModuleEdep=5.0*units.MeV,OutputLevel=DEBUG)
+
+ec_barrel_cluster = IslandCluster("ec_barrel_cluster",
+ inputHitCollection="RecEcalBarrelHits",
+ outputClusterCollection="EcalBarrelClusters",
+ splitHitCollection="splitEcalBarrelHitCollection",
+ minClusterCenterEdep=1*units.MeV,
+ groupRange=2.0,
+ OutputLevel=DEBUG)
+crystal_ec_cluster = IslandCluster("crystal_ec_cluster",
+ inputHitCollection="RecoEcalHits",
+ outputClusterCollection="EcalClusters",
+ minClusterCenterEdep=30*units.MeV,
+ groupRange=2.0,
+ OutputLevel=DEBUG)
+
+clusterreco = RecoCoG("cluster_reco", clusterCollection="EcalClusters", logWeightBase=4.2, moduleDimZName="CrystalBox_z_length",OutputLevel=DEBUG)
+
+out = PodioOutput("out", filename=output_rec_file)
+
+out.outputCommands = ["keep *"]
+
+ApplicationMgr(
+ TopAlg = [podioinput, copier, calcopier,
+ ecdigi, emcaldigi,
+ crystal_ec_reco, ecal_reco,
+ ec_barrel_cluster, crystal_ec_cluster, clusterreco, out],
+ EvtSel = 'NONE',
+ EvtMax = n_events,
+ ExtSvc = [podioevent],
+ OutputLevel=DEBUG
+ )
diff --git a/clustering/scripts/gen_central_electrons.cxx b/clustering/scripts/gen_central_electrons.cxx
new file mode 100644
index 00000000..06867dfa
--- /dev/null
+++ b/clustering/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( 70.0*(M_PI/180.0));
+ double cos_theta_max = std::cos(110.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;
+}
--
GitLab