Added clustering script

	modified:   clustering/barrel_clusters.sh
	modified:   clustering/options/calorimeter_clustering.py
	new file:   clustering/scripts/barrel_clusters.cxx

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/.
 

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],

clustering/scripts/barrel_clusters.cxx

+#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;
+}

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);