Added dataset, reconstruction, test analysis scripts

ecal/ecal_config.yml

@@ -25,3 +25,18 @@ ecal_1_emcal_pi0s:
   stage: run
   script:
     - bash ecal/emcal_pi0s.sh
+
+emcal_barrel_electrons:
+  image: eicweb.phy.anl.gov:4567/eic/juggler/juggler:$JUGGLER_TAG
+  tags:
+    - silicon
+  needs: ["configure"]
+  timeout: 48 hours
+  artifacts:
+    expire_in: 20 weeks
+    paths:
+      - results/
+  stage: run
+  script:
+    - bash ecal/full_emcal_electrons.sh
+

ecal/full_emcal_electrons.sh

+#!/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
+
+if [[ ! -n  "${E_start}" ]] ; then
+  export E_start=0.5
+fi
+
+if [[ ! -n  "${E_end}" ]] ; then
+  export E_end=5.0
+fi
+
+export JUGGLER_FILE_NAME_TAG="emcal_uniform_electrons"
+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
+# note datasets is now only used to develop datasets.
+#git clone https://eicweb.phy.anl.gov/EIC/datasets.git datasets
+root -b -q "ecal/scripts/full_emcal_electrons.cxx(${JUGGLER_N_EVENTS}, ${E_start}, ${E_end}, \"${JUGGLER_FILE_NAME_TAG}.hepmc\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running script"
+  exit 1
+fi
+
+pushd ${JUGGLER_DETECTOR}
+ls -l
+# run geant4 simulations
+npsim --runType batch \
+      -v WARNING \
+      --part.minimalKineticEnergy 1*GeV  \
+      --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
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running npdet"
+  exit 1
+fi
+xenv -x ${JUGGLER_INSTALL_PREFIX}/Juggler.xenv \
+  gaudirun.py ../ecal/options/full_em_calorimeter_reco.py
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running juggler"
+  exit 1
+fi
+ls -l
+popd
+
+pwd
+mkdir -p results
+
+root -b -q "ecal/scripts/full_emcal_electrons_analysis.cxx(\"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\")"
+if [[ "$?" -ne "0" ]] ; then
+  echo "ERROR running analysis script"
+  exit 1
+fi
+
+root_filesize=$(stat --format=%s "${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}")
+if [[ "${JUGGLER_N_EVENTS}" -lt "500" ]] ; then 
+  # file must be less than 10 MB to upload
+  if [[ "${root_filesize}" -lt "10000000" ]] ; then 
+    cp ${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE} results/.
+  fi
+fi
+

ecal/options/full_em_calorimeter_reco.py

@@ -33,15 +33,13 @@ from Configurables import Jug__Base__InputCopier_dd4pod__Geant4ParticleCollectio
 from Configurables import Jug__Base__InputCopier_dd4pod__CalorimeterHitCollection_dd4pod__CalorimeterHitCollection_ as CalCopier
 
 from Configurables import Jug__Digi__CrystalEndcapsDigi as CrystalEndcapsDigi
-from Configurables import Jug__Digi__EMCalorimeterDigi as EMCalorimeterDigi
+from Configurables import Jug__Digi__EcalTungstenSamplingDigi as EcalTungstenSamplingDigi
 
 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
-from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
 
-podioinput = PodioInput("PodioReader", collections=["mcparticles","CrystalEcalHits","EcalBarrelHits"], OutputLevel=DEBUG)
+podioinput = PodioInput("PodioReader", collections=["mcparticles","CrystalEcalHits","EcalBarrelHits","EcalEndcapHits"], OutputLevel=DEBUG)
 
 ## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
 copier = MCCopier("MCCopier", 
@@ -54,9 +52,16 @@ calcopier = CalCopier("CalCopier",
         outputCollection="CrystalEcalHits2",
         OutputLevel=DEBUG)
 
-embarreldigi = EMCalorimeterDigi("ec_barrel_digi", 
+embarreldigi = EcalTungstenSamplingDigi("ec_barrel_digi", 
         inputHitCollection="EcalBarrelHits", 
         outputHitCollection="RawEcalBarrelHits",
+        energyResolution=0.11,
+        OutputLevel=DEBUG)
+
+emendcapdigi = EcalTungstenSamplingDigi("ec_endcap_digi", 
+        inputHitCollection="EcalEndcapHits", 
+        outputHitCollection="RawEcalEndcapHits",
+        energyResolution=0.07,
         OutputLevel=DEBUG)
 
 emcaldigi = CrystalEndcapsDigi("ecal_digi", 
@@ -68,11 +73,6 @@ emcalreco = CrystalEndcapsReco("ecal_reco",
         outputHitCollection="RecoEcalHits",
         minModuleEdep=0.00001*units.MeV)
 
-embarrelreco = EMCalReconstruction("ecal_reco", 
-        inputHitCollection="RawEcalBarrelHits", 
-        outputHitCollection="RecEcalBarrelHits",
-        minModuleEdep=0.00001*units.MeV)
-
 emcalcluster = IslandCluster("emcal_cluster", 
         inputHitCollection="RecoEcalHits", 
         outputClusterCollection="EcalClusters",
@@ -90,7 +90,7 @@ out.outputCommands = ["keep *"]
 
 ApplicationMgr(
     TopAlg = [podioinput, copier, calcopier,
-              embarreldigi, emcaldigi, emcalreco, embarrelreco, emcalcluster, clusterreco, out],
+              embarreldigi, emendcapdigi, emcaldigi, emcalreco, emcalcluster, clusterreco, out],
     EvtSel = 'NONE',
     EvtMax   = n_events,
     ExtSvc = [podioevent],

ecal/scripts/full_emcal_electrons.cxx

+//////////////////////////////////////////////////////////////
+// Tungsten Sampling EMCAL detector
+// Barrel and Endcap 
+// Single Electron dataset
+// J.KIM 02/12/2021
+//
+//////////////////////////////////////////////////////////////
+#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;
+
+void full_emcal_electrons(int n_events = 1e6, double e_start = 0.0, double e_end = 30.0, const char* out_fname = "./data/emcal_electron_0GeVto30GeV_1000kEvt.hepmc")
+{
+  WriterAscii hepmc_output(out_fname);
+  int events_parsed = 0;
+  GenEvent evt(Units::GEV, Units::MM);
+
+  // Random number generator
+  TRandom *r1 = new TRandom();
+
+  // Constraining the solid angle, but larger than that subtended by the detector
+  // https://indico.bnl.gov/event/7449/contributions/35966/attachments/27177/41430/EIC-DWG-Calo-03192020.pdf
+  // See a figure on slide 26
+  double cos_theta_min = std::cos(M_PI * (45.0 / 180.0));
+  double cos_theta_max = std::cos(M_PI * (135.0 / 180.0));
+
+  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
+    // Momentum starting at 0 GeV/c to 30 GeV/c
+    Double_t p        = r1->Uniform(e_start, e_end);
+    Double_t phi      = r1->Uniform(0.0, 2.0 * M_PI);
+    Double_t costheta = r1->Uniform(cos_theta_min, cos_theta_max);
+    Double_t theta    = std::acos(costheta);
+    Double_t px       = p * std::cos(phi) * std::sin(theta);
+    Double_t py       = p * std::sin(phi) * std::sin(theta);
+    Double_t pz       = p * std::cos(theta);
+
+    // Generates random vectors, uniformly distributed over the surface of a
+    // sphere of given radius, in this case momentum.
+    //r1->Sphere(px, py, pz, p);
+
+    // 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;
+}

ecal/scripts/full_emcal_electrons_analysis.cxx

+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/ClusterData.h"
+
+#include "TCanvas.h"
+#include "TStyle.h"
+#include "TMath.h"
+#include "TH1.h"
+#include "TH1D.h"
+
+// If root cannot find a dd4pod header make sure to add the include dir to ROOT_INCLUDE_PATH:
+//  export ROOT_INCLUDE_PATH=$HOME/include:$ROOT_INCLUDE_PATH
+//  export ROOT_INCLUDE_PATH=/home/jihee/stow/development/include:$ROOT_INCLUDE_PATH
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+void full_emcal_electrons_analysis(const char* input_fname = "rec_emcal_uniform_electrons.root")
+{
+  // Setting for graphs
+  gROOT->SetStyle("Plain");
+  gStyle->SetOptFit(1);
+  gStyle->SetLineWidth(2);
+  gStyle->SetPadTickX(1);
+  gStyle->SetPadTickY(1);
+  gStyle->SetPadGridX(1);
+  gStyle->SetPadGridY(1);
+  gStyle->SetPadLeftMargin(0.14);
+  gStyle->SetPadRightMargin(0.14);
+
+  ROOT::EnableImplicitMT();
+  ROOT::RDataFrame d0("events", input_fname);
+
+  // Thrown Energy [GeV]
+  auto E_thr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+  	std::vector<float> result;
+    	result.push_back(TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass));
+  return result;
+  };
+
+  // Define variables
+  auto d1 = d0.Define("E_thr",    E_thr,    {"mcparticles2"})
+	      ;
+
+  // Define Histograms
+  auto hEthr     = d1.Histo1D({"hEthr",     "Thrown Energy; E_{thr}[GeV]; Events",          100, -0.5, 30.5},   "E_thr");
+
+  // Event Counts
+  auto nevents_thrown    = d1.Count();
+
+  // Print out number of events
+  std::cout << "Number of Events: " << (*nevents_thrown) << " = all thrown events \n";
+
+  // Draw Histogram
+  TCanvas *c1 = new TCanvas("c1", "c1", 500, 500);
+  hEthr->GetYaxis()->SetTitleOffset(1.4);
+  hEthr->SetLineWidth(2);
+  hEthr->SetLineColor(kBlue);
+  hEthr->DrawClone();
+  c1->SaveAs("results/emcal_electrons_Ethr.png");
+  c1->SaveAs("results/emcal_electrons_Ethr.pdf");
+}