diff --git a/clustering/barrel_clusters.sh b/clustering/barrel_clusters.sh
index c04b3ae3f154a77c4d5e657566baad9384d61875..ecee95847ffacd5531826698aca1758028442886 100644
--- a/clustering/barrel_clusters.sh
+++ b/clustering/barrel_clusters.sh
@@ -28,15 +28,12 @@ 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
+### run geant4 simulations
npsim --runType batch \
--part.minimalKineticEnergy 1000*GeV \
-v WARNING \
@@ -51,6 +48,9 @@ ls -l
popd
pwd
-mkdir -p results
-cp topside/${JUGGLER_REC_FILE} results/.
+mkdir -p results/clustering
+
+root -b -q "clustering/scripts/barrel_clusters.cxx(\"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\")"
+
+cp topside/${JUGGLER_REC_FILE} results/clustering/.
diff --git a/clustering/options/calorimeter_clustering.py b/clustering/options/calorimeter_clustering.py
index 1dce32ddf09a6ad140a85db24be78200c574c31a..281a201ab329a279b7ce9fe17e1cbc3babec8d75 100644
--- a/clustering/options/calorimeter_clustering.py
+++ b/clustering/options/calorimeter_clustering.py
@@ -28,6 +28,7 @@ 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__SimpleClustering as SimpleClustering
from Configurables import Jug__Reco__CalorimeterIslandCluster as IslandCluster
from Configurables import Jug__Reco__ClusterRecoCoG as RecoCoG
@@ -43,7 +44,7 @@ ecdigi = EMCalorimeterDigi("ec_barrel_digi", inputHitCollection="EcalBarrelHits"
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)
+ minModuleEdep=0.0*units.MeV,OutputLevel=DEBUG)
ec_barrel_cluster = IslandCluster("ec_barrel_cluster",
inputHitCollection="RecEcalBarrelHits",
@@ -59,6 +60,12 @@ crystal_ec_cluster = IslandCluster("crystal_ec_cluster",
groupRange=2.0,
OutputLevel=DEBUG)
+simple_cluster = SimpleClustering("simple_cluster",
+ inputHitCollection="RecEcalBarrelHits",
+ outputClusters="SimpleClusters",
+ OutputLevel=DEBUG)
+
+
clusterreco = RecoCoG("cluster_reco", clusterCollection="EcalClusters", logWeightBase=4.2, moduleDimZName="CrystalBox_z_length",OutputLevel=DEBUG)
out = PodioOutput("out", filename=output_rec_file)
@@ -69,7 +76,9 @@ ApplicationMgr(
TopAlg = [podioinput, copier, calcopier,
ecdigi, emcaldigi,
crystal_ec_reco, ecal_reco,
- ec_barrel_cluster, crystal_ec_cluster, clusterreco, out],
+ ec_barrel_cluster, crystal_ec_cluster, clusterreco,
+ simple_cluster,
+ out],
EvtSel = 'NONE',
EvtMax = n_events,
ExtSvc = [podioevent],
diff --git a/clustering/scripts/barrel_clusters.cxx b/clustering/scripts/barrel_clusters.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..2cfb8783cff975d8d369e3904a79bbf9ac7dd578
--- /dev/null
+++ b/clustering/scripts/barrel_clusters.cxx
@@ -0,0 +1,137 @@
+#include <iostream>
+#include "ROOT/RDataFrame.hxx"
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "eicd/ClusterCollection.h"
+#include "eicd/ClusterData.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
+
+using ROOT::RDataFrame;
+using namespace ROOT::VecOps;
+
+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 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 eta = [](ROOT::VecOps::RVec<dd4pod::Geant4ParticleData> const& in) {
+ std::vector<float> 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.Eta());
+ }
+ return result;
+};
+
+auto delta_E = [](std::vector<ROOT::Math::PxPyPzMVector> const& thrown, const std::vector<eic::ClusterData>& clusters) {
+ std::vector<double> result;
+ double best = 1000000.0;
+ for (const auto& p : thrown) {
+ for (const auto& c : clusters) {
+ double dE = (p.E() - c.energy);
+ result.push_back(dE );
+ if( dE < best) {
+ best = dE;
+ }
+ }
+ }
+ return best;
+};
+
+
+void barrel_clusters(const char* in_fname = "topside/rec_barrel_clusters.root")
+{
+ ROOT::EnableImplicitMT();
+ ROOT::RDataFrame df("events", in_fname);
+
+ auto d0 = df.Define("isThrown", "mcparticles2.genStatus == 1")
+ .Define("thrownParticles", "mcparticles2[isThrown]")
+ .Define("thrownP", fourvec, {"thrownParticles"})
+ .Define("thrownEta", eta, {"thrownParticles"})
+ .Define("thrownTheta", theta, {"thrownP"})
+ .Define("thrownMomentum", momentum, {"thrownP"})
+ .Define("delta_E", delta_E, {"thrownP","SimpleClusters"})
+ .Define("nclusters", "SimpleClusters.size()")
+ .Define("Ecluster",
+ [](const std::vector<eic::ClusterData>& in) {
+ std::vector<double> res;
+ for (const auto& i : in)
+ res.push_back(i.energy);
+ return res;
+ },
+ {"SimpleClusters"});
+
+ auto d1 = d0.Filter("nclusters==1");
+ auto c_nclusters1 = d1.Count();
+ auto c_thrown = d0.Count();
+
+ auto h_eta_thrown = d0.Histo1D({"h_eta_thrown", " ; #eta ", 100, -5.0, 0.0}, "thrownEta");
+ auto h_theta_thrown = d0.Histo1D({"h_theta_thrown", "; #theta", 100, 30.0, 180.0}, "thrownTheta");
+ auto h_momentum_thrown = d0.Histo1D({"h_momentum_thrown", "; E [GeV]", 100, 0.0, 30.0}, "thrownMomentum");
+ auto h_nclusters = d0.Histo1D({"h_nclusters", "; N clusters", 6, 0, 6}, "nclusters");
+ auto h_Ecluster = d0.Histo1D({"h_Ecluster", "; cluster E [GeV]", 100, 0, 1}, "Ecluster");
+ auto h_Ecluster1 = d1.Histo1D({"h_Ecluster1", "One cluster events; cluster E [GeV]", 100, 0, 30}, "Ecluster");
+ auto h_momentum_thrown1 = d1.Histo1D({"h_momentum_thrown", "; E [GeV]", 100, 0.0, 30.0},"thrownMomentum");
+
+ auto h_delta_E = d0.Histo1D({"h_delta_E", "; #Delta E [GeV]", 100, -10, 10}, "delta_E");
+
+ //auto h_Ecluster2 = d1.Filter("thrownMomentum > 4").Histo1D({"h_Ecluster1", "One cluster events; cluster E [GeV]",100, 0,30},"Ecluster");
+ //auto h_momentum_thrown2 = d1.Filter("thrownMomentum > 4").Histo1D({"h_momentum_thrown", "; E [GeV]", 100, 0.0, 30.0},"thrownMomentum");
+
+ auto c = new TCanvas();
+
+ h_eta_thrown->DrawCopy();
+ c->SaveAs("results/clustering/barrel_clusters_etaThrown.png");
+
+ h_theta_thrown->DrawCopy();
+ c->SaveAs("results/clustering/barrel_clusters_thetaThrown.png");
+
+ h_nclusters->DrawCopy();
+ c->SaveAs("results/clustering/barrel_clusters_nclusters.png");
+
+ h_Ecluster->DrawCopy();
+ h_Ecluster1->SetLineColor(2);
+ h_Ecluster1->DrawCopy("same");
+ //h_Ecluster2->SetLineColor(4);
+ //h_Ecluster2->DrawCopy("same");
+ c->SaveAs("results/clustering/barrel_clusters_Ecluster.png");
+
+ h_momentum_thrown->DrawCopy();
+ c->SaveAs("results/clustering/barrel_clusters_momentum_thrown.png");
+
+ h_delta_E->DrawCopy();
+ c->SaveAs("results/clustering/barrel_clusters_delta_E.png");
+
+ std::cout << (*c_nclusters1) << " / " << (*c_thrown) << " = single cluster events/all \n";
+
+ //c->SaveAs("results/crystal_cal_electrons_Ecluster.png");
+ //std::string outfilename = "rdf_test.root";
+ //df2.Snapshot("events", outfilename, {"MCParticles_pt", "mcparticles"});
+ return 0;
+}
diff --git a/clustering/scripts/gen_central_electrons.cxx b/clustering/scripts/gen_central_electrons.cxx
index 06867dfaea7507522bf41a0ffc3ca87d413fb301..e3e08bd44007fd907bbe0188d11ac2686e0897a4 100644
--- a/clustering/scripts/gen_central_electrons.cxx
+++ b/clustering/scripts/gen_central_electrons.cxx
@@ -42,7 +42,7 @@ void gen_central_electrons(int n_events = 100,
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 th = M_PI/2.1;//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);