diff --git a/.gitignore b/.gitignore
index a743d854509495f714935dff92e7324677b2aa96..831c69456261a98ab4d5a87b60908cdb83f9c6e2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -2,3 +2,4 @@ results/*
data/*
datasets/*
data
+calorimeters/test/
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index c8127fd89dc164e6adf67dc388022136369b18de..a2502a9efcb27b5705d011796db01f5b6b0667ce 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -58,6 +58,15 @@ roman_pot_simu:
- cp NPDet/src/GenericDetectors/trackers/compact/materials.xml ./.
- bash trackers/roman_pot_simu.sh
+zdc_simulation:
+ stage: simulate
+ tags:
+ - sodium
+ script:
+ - cp NPDet/src/GenericDetectors/calorimeters/compact/elements.xml ./.
+ - cp NPDet/src/GenericDetectors/calorimeters/compact/materials.xml ./.
+ - bash calorimeters/run_simulation_zdc.sh
+
cal_test_1_dummy_test2:
stage: benchmarks
tags:
@@ -106,6 +115,14 @@ roman_pot_benchmark:
- root -b -q trackers/simple_tracking.cxx+
allow_failure: true
+zdc_benchmark:
+ stage: benchmarks
+ tags:
+ - sodium
+ script:
+ - ls -lrth sim_output
+ - root -b -q calorimeters/simple_checking.cxx+
+ allow_failure: true
deploy_results:
stage: deploy
diff --git a/calorimeters/ZDC_example.xml b/calorimeters/ZDC_example.xml
new file mode 100644
index 0000000000000000000000000000000000000000..fce8c083eb7d8030191fb1f81a93fb1ef511f9be
--- /dev/null
+++ b/calorimeters/ZDC_example.xml
@@ -0,0 +1,174 @@
+<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0"
+ xmlns:xs="http://www.w3.org/2001/XMLSchema"
+ xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd">
+
+ <!-- Some information about detector -->
+ <info name="ZDC_example" title="Zero Degree Calorimeter detector example"
+ author="Jihee Kim"
+ url="https://eicweb.phy.anl.gov/EIC/NPDet"
+ status="development"
+ version="v2.0 2020-07-22">
+ <comment>Zero Degree Calorimeter detector</comment>
+ </info>
+
+ <!-- Use DD4hep elements and materials definitions -->
+ <includes>
+ <gdmlFile ref="elements.xml"/>
+ <gdmlFile ref="materials.xml"/>
+ </includes>
+
+ <!-- Define the dimensions of the world volume -->
+ <define>
+ <constant name="world_side" value="50*m"/>
+ <constant name="world_x" value="world_side"/>
+ <constant name="world_y" value="world_side"/>
+ <constant name="world_z" value="world_side"/>
+
+ <constant name="CrossingAngle" value="0.020*rad"/>
+ <constant name="tracker_region_zmax" value="5*m"/>
+ <constant name="tracker_region_rmax" value="5*m"/>
+
+ <constant name="offset_ZDC" value="5.0*mm"/>
+ <constant name="st_length" value="20.0*mm"/>
+ <constant name="lt_length" value="40.0*mm"/>
+ <constant name="st_ZDC_x_pos" value="0.60*m"/>
+ <constant name="st_ZDC_y_pos" value="0.0*m"/>
+ <constant name="st_ZDC_z_pos" value="34.0*m"/>
+ <constant name="lt_ZDC_x_pos" value="0.60*m"/>
+ <constant name="lt_ZDC_y_pos" value="offset_ZDC + (st_length+lt_length)/sqrt(2)"/>
+ <constant name="lt_ZDC_z_pos" value="34.0*m"/>
+ </define>
+
+ <limits>
+ </limits>
+
+ <regions>
+ </regions>
+
+ <!-- Common Generic visualization attributes -->
+ <comment>Common Generic visualization attributes</comment>
+ <display>
+ <vis name="InvisibleNoDaughters" showDaughters="false" visible="false"/>
+ <vis name="InvisibleWithDaughters" showDaughters="true" visible="false"/>
+ <vis name="GreenVis" alpha="0.5" r= "0.0" g="1.0" b="0.0" showDaughters="true" visible="true"/>
+ <vis name="RedVis" alpha="0.5" r= "1.0" g="0.0" b="0.0" showDaughters="true" visible="true"/>
+ <vis name="BlueVis" alpha="0.5" r= "0.0" g="0.0" b="1.0" showDaughters="true" visible="true"/>
+ <vis name="OrangeVis" alpha="0.5" r= "1.0" g="0.45" b="0.0" showDaughters="true" visible="true"/>
+ <vis name="RedGreenVis" alpha="0.5" r= "1.0" g="1.0" b="0.0" showDaughters="true" visible="true"/>
+ <vis name="BlueGreenVis" alpha="0.5" r= "0.0" g="1.0" b="1.0" showDaughters="true" visible="true"/>
+ <vis name="PurpleVis" alpha="0.5" r= "1.0" g="0.0" b="1.0" showDaughters="true" visible="true"/>
+ <vis name="DoubleRedG" alpha="0.5" r= "2.0" g=".10" b="0.0" showDaughters="true" visible="true"/>
+ <vis name="RBG015" alpha="0.5" r= "0.0" g=".2" b="1.0" showDaughters="true" visible="true"/>
+ <vis name="RBG510" alpha="0.5" r= "1.0" g=".2" b="0.0" showDaughters="true" visible="true"/>
+ <vis name="RBG" alpha="0.5" r= "1.0" g="1.0" b="1.0" showDaughters="true" visible="true"/>
+ <vis name="GrayVis" alpha="0.5" r= "0.75" g="0.75" b="0.75" showDaughters="true" visible="true"/>
+ </display>
+
+ <!-- Define detector -->
+ <detectors>
+ <detector id="1" name="smallZDC" type="ZDC" readout="ZDCHits" vis="RedVis">
+ <position x="st_ZDC_x_pos" y="st_ZDC_y_pos" z="st_ZDC_z_pos"/>
+ <dimensions x = "st_length" y = "st_length"/>
+ <layer repeat="2">
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ <slice name="Scint_slice" material="PlasticScint" thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
+ </layer>
+ <layer repeat="1">
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ <slice name="SciFi_belt" material="PlasticScint" thickness="1*mm" vis = "RedVis" sensitive = "true"/>
+ <slice name="SciFi_belt" material="PlasticScint" thickness="1*mm" vis = "RedVis" sensitive = "true"/>
+ </layer>
+ <layer repeat="2">
+ <slice name="Scint_slice" material="PlasticScint" thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ </layer>
+ <layer repeat="2">
+ <slice name="SciFi_belt" material="PlasticScint" thickness="1*mm" vis = "RedVis" sensitive = "true"/>
+ </layer>
+ <layer repeat="7">
+ <slice name="Scint_slice" material="PlasticScint" thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ </layer>
+ <layer repeat="1">
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ <slice name="Scint_slice" material="PlasticScint" thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ </layer>
+ <layer repeat="2">
+ <slice name="SciFi_belt" material="PlasticScint" thickness="1*mm" vis = "RedVis" sensitive = "true"/>
+ </layer>
+ <layer repeat="3">
+ <slice name="Scint_slice" material="PlasticScint" thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ </layer>
+ <layer repeat="2">
+ <slice name="SciFi_belt" material="PlasticScint" thickness="1*mm" vis = "RedVis" sensitive = "true"/>
+ </layer>
+ <layer repeat="1">
+ <slice name="Scint_slice" material="PlasticScint" thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ </layer>
+ </detector>
+ <detector id="2" name="largeZDC" type="ZDC" readout="ZDCHits" vis="RedVis">
+ <position x="lt_ZDC_x_pos" y="lt_ZDC_y_pos" z="lt_ZDC_z_pos"/>
+ <dimensions x = "lt_length" y = "lt_length"/>
+ <layer repeat="2">
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ <slice name="Scint_slice" material="PlasticScint" thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
+ </layer>
+ <layer repeat="1">
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ <slice name="SciFi_belt" material="PlasticScint" thickness="1*mm" vis = "RedVis" sensitive = "true"/>
+ <slice name="SciFi_belt" material="PlasticScint" thickness="1*mm" vis = "RedVis" sensitive = "true"/>
+ </layer>
+ <layer repeat="2">
+ <slice name="Scint_slice" material="PlasticScint" thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ </layer>
+ <layer repeat="2">
+ <slice name="SciFi_belt" material="PlasticScint" thickness="1*mm" vis = "RedVis" sensitive = "true"/>
+ </layer>
+ <layer repeat="7">
+ <slice name="Scint_slice" material="PlasticScint" thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ </layer>
+ <layer repeat="1">
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ <slice name="Scint_slice" material="PlasticScint" thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ </layer>
+ <layer repeat="2">
+ <slice name="SciFi_belt" material="PlasticScint" thickness="1*mm" vis = "RedVis" sensitive = "true"/>
+ </layer>
+ <layer repeat="3">
+ <slice name="Scint_slice" material="PlasticScint" thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ </layer>
+ <layer repeat="2">
+ <slice name="SciFi_belt" material="PlasticScint" thickness="1*mm" vis = "RedVis" sensitive = "true"/>
+ </layer>
+ <layer repeat="1">
+ <slice name="Scint_slice" material="PlasticScint" thickness="3*mm" vis = "BlueVis" sensitive = "true"/>
+ <slice name="Tungsten_slice" material="TungstenDens24" thickness="7*mm" vis = "GrayVis"/>
+ </layer>
+ </detector>
+ </detectors>
+
+ <!-- Definition of the readout segmentation/definition -->
+ <readouts>
+ <readout name="ZDCHits">
+ <segmentation type="CartesianGridXY" grid_size_x="1.0*mm" grid_size_y="1.0*mm" />
+ <id>system:8,layer:12,slice:12,x:48:-8,y:-8</id>
+ </readout>
+ </readouts>
+
+ <plugins>
+ </plugins>
+
+ <fields>
+ </fields>
+</lccdd>
diff --git a/calorimeters/run_simulation_zdc.sh b/calorimeters/run_simulation_zdc.sh
new file mode 100755
index 0000000000000000000000000000000000000000..76216c4b313af59c29eec22e294b8d4c0cc4615a
--- /dev/null
+++ b/calorimeters/run_simulation_zdc.sh
@@ -0,0 +1,8 @@
+#!/bin/bash
+
+source /usr/local/bin/thisdd4hep.sh
+
+ddsim --runType batch --numberOfEvents 100 \
+ --compactFile ./calorimeters/ZDC_example.xml \
+ --inputFiles ./data/zdc_photons.hepmc \
+ --outputFile ./sim_output/output_zdc_photons.root
diff --git a/calorimeters/simple_checking.cxx b/calorimeters/simple_checking.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..774317ffbacf23cc093a55b6d4afaf9fcf09727a
--- /dev/null
+++ b/calorimeters/simple_checking.cxx
@@ -0,0 +1,124 @@
+//R__LOAD_LIBRARY(libfmt.so)
+//#include "fmt/core.h"
+R__LOAD_LIBRARY(libDDG4IO.so)
+//
+//#include "DD4hep/Detector.h"
+#include "DDG4/Geant4Data.h"
+//#include "DDRec/CellIDPositionConverter.h"
+//#include "DDRec/SurfaceManager.h"
+//#include "DDRec/Surface.h"
+#include "ROOT/RDataFrame.hxx"
+//
+//#include "lcio2/MCParticleData.h"
+//#include "lcio2/ReconstructedParticleData.h"
+
+//#include "Math/Vector3D.h"
+//#include "Math/Vector4D.h"
+//#include "Math/VectorUtil.h"
+#include "TCanvas.h"
+//#include "TLegend.h"
+//#include "TMath.h"
+//#include "TRandom3.h"
+//#include "TFile.h"
+//#include "TH1F.h"
+//#include "TH1D.h"
+//#include "TTree.h"
+#include "TChain.h"
+//#include "TF1.h"
+#include <random>
+//#include "lcio2/TrackerRawDataData.h"
+//#include "lcio2/TrackerRawData.h"
+
+void simple_checking(const char* fname = "sim_output/output_zdc_photons.root"){
+std::cout << "testing 1\n";
+ ROOT::EnableImplicitMT(); // Tell ROOT you want to go parallel
+ //using namespace lcio2;
+ double degree = TMath::Pi()/180.0;
+
+ TChain* t = new TChain("EVENT");
+ t->Add(fname);
+
+ ROOT::RDataFrame d0(*t);//, {"GEMTrackerHintits","MCParticles"});
+std::cout << "testing 2\n";
+ //std::cout << t->GetBranch("GEMTrackerHits")->GetClassName() << std::endl;
+ //std::vector<dd4hep::sim::Geant4Tracker::Hit*>
+
+ // -------------------------
+ // Get the DD4hep instance
+ // Load the compact XML file
+ // Initialize the position converter tool
+ //dd4hep::Detector& detector = dd4hep::Detector::getInstance();
+ //detector.fromCompact("gem_tracker_disc.xml");
+ //dd4hep::rec::CellIDPositionConverter cellid_converter(detector);
+
+ //// -------------------------
+ //// Get the surfaces map
+ //dd4hep::rec::SurfaceManager& surfMan = *detector.extension<dd4hep::rec::SurfaceManager>() ;
+ //auto surfMap = surfMan.map( "world" ) ;
+
+ auto nhits = [] (std::vector<dd4hep::sim::Geant4Calorimeter::Hit*>& hits){ return (int) hits.size(); };
+std::cout << "testing 3\n";
+ //for(const auto& h: hits){
+ // //std::cout << (h->position/10.0) << std::endl;
+ // //std::cout << cellid_converter.position(h->cellID) << std::endl;
+ // //dd4hep::rec::SurfaceMap::const_iterator
+ // const auto si = _surfMap.find( cellid_converter.findContext(h->cellID)->identifier ); //identifier=volumeID
+ // dd4hep::rec::ISurface* surf = (si != _surfMap.end() ? si->second : 0);
+ // dd4hep::rec::Vector3D pos = surf->origin();//fit_global(pivot[0],pivot[1],pivot[2]);
+ // //std::cout << pos.x() << ", " << pos.y() << ", " << pos.z()<< std::endl;
+ // // transform lcio units to dd4hep units, see documentation for other functions
+ // //DDSurfaces::Vector2D fit_local = surf->globalToLocal( dd4hep::mm * fit_global );
+ //}
+ // return hits.size(); };
+
+ //auto digitize_gem_hits =
+ // [&](const std::vector<dd4hep::sim::Geant4Tracker::Hit*>& hits) {
+ // std::vector<lcio2::TrackerRawDataData> digi_hits;
+ // std::normal_distribution<> time_dist(0,2.0);
+ // std::normal_distribution<> adc_dist(5.0,3.0);
+
+ // std::map<int64_t,lcio2::TrackerRawDataData> hits_by_id;
+ // for(const auto& h: hits) {
+ // //lcio2::TrackerRawDataData ahit;
+ // auto& ahit = hits_by_id[(int64_t)h->cellID];
+ // auto pos = h->position/10.0; //cm
+
+ // ahit.cellID0 = h->cellID;
+ // ahit.cellID1 = h->cellID;
+ // ahit.channelID = h->cellID;
+ // //fmt::print("{} vs {} vs {}\n", id1, id2,id3);
+ // fmt::print("{} vs {}\n", h->cellID, ahit.cellID0);
+ // // time is not kept from dd4hep hit, instead using z position as crude substitute
+ // ahit.time = pos.z() + time_dist(gen);
+ // ahit.adc = adc_dist(gen);
+ // //digi_hits.push_back(ahit);
+ // }
+ // for (auto& [cell, hit] : hits_by_id) {
+ // //fmt::print("{} vs {}\n", cell, hit.cellID0);
+ // digi_hits.push_back(hit);
+ // }
+ // return digi_hits;
+ // };
+
+ auto d1 = d0.Define("nhits", nhits, {"ZDCHits"})
+ //.Filter([](int n){ return (n>4); },{"nhits"})
+ //.Define("delta",hit_position, {"GEMTrackerHits"})
+ //.Define("RawTrackerHits", digitize_gem_hits, {"GEMTrackerHits"})
+ ;
+std::cout << "testing 4\n";
+ auto h0 = d1.Histo1D(TH1D("h0", "nhits; ", 20, 0,20), "nhits");
+
+ auto n0 = d1.Filter([](int n){ return (n>=0); },{"nhits"}).Count();
+
+ TCanvas* c = new TCanvas();
+
+ //d1.Snapshot("digitized_EVENT","test_gem_tracker_digi.root");
+std::cout << "testing 5\n";
+ std::cout << *n0 << " events with nonzero hits\n";
+std::cout << "testing 6\n";
+ if(*n0<5) {
+ std::quick_exit(1);
+ }
+
+}
+