new file: bin/crop_pdf

	new file:   bin/crop_views
	modified:   bin/download_views
	new file:   ecal/scripts/crystal_cal_electrons.cxx
	new file:   ecal/scripts/rdf_test.cxx

bin/crop_pdf

+#!/bin/bash
+
+if [ $# -lt 5 ]
+then
+  echo "Usage: `basename $0` <pdf-file> <x_min> <x_max> <y_min> <y_max>"
+  echo "Notes:"
+  echo " - all coordinates are absolute; no calculation of width/height necessary"
+  echo " - use 'gv' to determine the coordinates"
+  exit 65
+fi
+
+file="$1"
+xmin="$2"
+xmax="$3"
+ymin="$4"
+ymax="$5"
+
+base="${file%.*}"
+outfile="${base}_cropped.pdf"
+
+echo "writing to: $outfile"
+
+gs \
+  -o $outfile \
+  -sDEVICE=pdfwrite \
+  -c "[/CropBox [$xmin $ymin $xmax $ymax] /PAGES pdfmark" \
+  -f $file
+

bin/crop_views

+#!/bin/bash
+
+mkdir -p cropped
+#convert view1.png -crop 1100x760+100+400 cropped/view1.png
+#
+#convert view2a.png -crop 1100x960+100+300 cropped/view2a.png
+#convert view2b.png -crop 1100x960+100+300 cropped/view2b.png
+#convert view2c.png -crop 1100x960+100+300 cropped/view2c.png
+#convert view2d.png -crop 1100x960+100+300 cropped/view2d.png
+#convert view2e.png -crop 1100x960+100+300 cropped/view2e.png
+
+gs -o cropped/view1.pdf -sDEVICE=pdfwrite \
+  -c "[/CropBox [51 250 550 590] /PAGES pdfmark" \
+  -f view1.pdf
+ 
+gs -o cropped/view2a.pdf -sDEVICE=pdfwrite \
+  -c "[/CropBox [50 175 550 675] /PAGES pdfmark" \
+  -f view2a.pdf
+
+gs -o cropped/view2b.pdf -sDEVICE=pdfwrite \
+  -c "[/CropBox [50 175 550 675] /PAGES pdfmark" \
+  -f view2b.pdf
+
+gs -o cropped/view2c.pdf -sDEVICE=pdfwrite \
+  -c "[/CropBox [50 175 550 675] /PAGES pdfmark" \
+  -f view2c.pdf
+
+gs -o cropped/view2d.pdf -sDEVICE=pdfwrite \
+  -c "[/CropBox [50 175 550 675] /PAGES pdfmark" \
+  -f view2d.pdf
+
+gs -o cropped/view2e.pdf -sDEVICE=pdfwrite \
+  -c "[/CropBox [50 175 550 675] /PAGES pdfmark" \
+  -f view2e.pdf
+

bin/download_views

@@ -7,3 +7,13 @@ curl -o view2c.png 'https://eicweb.phy.anl.gov/api/v4/projects/390/jobs/artifact
 curl -o view2d.png 'https://eicweb.phy.anl.gov/api/v4/projects/390/jobs/artifacts/master/raw/images/view2d.png?job=compile'
 curl -o view2e.png 'https://eicweb.phy.anl.gov/api/v4/projects/390/jobs/artifacts/master/raw/images/view2e.png?job=compile'
 curl -o view3.png  'https://eicweb.phy.anl.gov/api/v4/projects/390/jobs/artifacts/master/raw/images/view3.png?job=compile'
+
+
+
+curl -o view1.pdf  'https://eicweb.phy.anl.gov/api/v4/projects/390/jobs/artifacts/master/raw/images/view1.pdf?job=compile'
+curl -o view2a.pdf 'https://eicweb.phy.anl.gov/api/v4/projects/390/jobs/artifacts/master/raw/images/view2a.pdf?job=compile'
+curl -o view2b.pdf 'https://eicweb.phy.anl.gov/api/v4/projects/390/jobs/artifacts/master/raw/images/view2b.pdf?job=compile'
+curl -o view2c.pdf 'https://eicweb.phy.anl.gov/api/v4/projects/390/jobs/artifacts/master/raw/images/view2c.pdf?job=compile'
+curl -o view2d.pdf 'https://eicweb.phy.anl.gov/api/v4/projects/390/jobs/artifacts/master/raw/images/view2d.pdf?job=compile'
+curl -o view2e.pdf 'https://eicweb.phy.anl.gov/api/v4/projects/390/jobs/artifacts/master/raw/images/view2e.pdf?job=compile'
+curl -o view3.pdf  'https://eicweb.phy.anl.gov/api/v4/projects/390/jobs/artifacts/master/raw/images/view3.pdf?job=compile'

ecal/emcal_electrons.sh

@@ -48,6 +48,7 @@ ls -l
 # run geant4 simulations
 npsim --runType batch \
       -v WARNING \
+      --part.minimalKineticEnergy 1000*GeV  \
       --numberOfEvents ${JUGGLER_N_EVENTS} \
       --compactFile ${JUGGLER_DETECTOR}.xml \
       --inputFiles ../${JUGGLER_FILE_NAME_TAG}.hepmc \
@@ -64,6 +65,8 @@ root -b -q "ecal/scripts/rec_emcal_electrons_reader.C(${E_start}, ${E_end}, \"${
 #root -b -q "ecal/scripts/makeplot.C(${E_start}, ${E_end}, \"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\", \"results/rec_${JUGGLER_FILE_NAME_TAG}.txt\")"
 #root -b -q "ecal/scripts/makeplot_input.C(\"${JUGGLER_DETECTOR}/${JUGGLER_SIM_FILE}\", \"results/sim_${JUGGLER_FILE_NAME_TAG}.txt\")"
 
+root -b -q "ecal/scripts/crystal_cal_electrons.cxx(\"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\")"
+
 #paste results/sim_${JUGGLER_FILE_NAME_TAG}.txt results/rec_${JUGGLER_FILE_NAME_TAG}.txt > results/eng_${JUGGLER_FILE_NAME_TAG}.txt
 #root -b -q "ecal/scripts/read_eng.C(\"results/eng_${JUGGLER_FILE_NAME_TAG}.root\", \"results/eng_${JUGGLER_FILE_NAME_TAG}.txt\")"
 #root -b -q "ecal/scripts/cal_eng_res.C(\"results/eng_${JUGGLER_FILE_NAME_TAG}.root\")"

ecal/emcal_pi0s.sh

@@ -41,6 +41,7 @@ echo "JUGGLER_SIM_FILE = ${JUGGLER_SIM_FILE}"
 # run geant4 simulations
 npsim --runType batch \
       -v WARNING \
+      --part.minimalKineticEnergy 1000*GeV  \
       --numberOfEvents ${JUGGLER_N_EVENTS} \
       --compactFile ${JUGGLER_DETECTOR}.xml \
       --inputFiles ../${JUGGLER_FILE_NAME_TAG}.hepmc \

ecal/scripts/crystal_cal_electrons.cxx

+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#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;
+};
+
+void crystal_cal_electrons(const char* in_fname = "topside/rec_emcal_uniform_electrons.root")
+{
+  ROOT::EnableImplicitMT();
+  ROOT::RDataFrame df("events", in_fname);
+
+  //TH1D* h1 = new TH1D("h1", "Scattering Angle(#theta)", 100, 130.0, 180.0);
+  //TH1D* h2 = new TH1D("h2", "Pseudo-rapidity(#eta)", 100, -5.0, 0.0);
+  //TH2D* h3 = new TH2D("h3", "Cluster E vs Pseudo-rapidity", 100, e_start - 0.5, e_end + 0.5, 100, -5.0, 0.0);
+  //TH1D* h4 = new TH1D("h4", "Reconstructed energy per event", 100, e_start - 0.5, e_end + 0.5);
+  //TH1D* h5 = new TH1D("h5", "Number of Clusters per event", 5, -0.5, 4.5);
+  //TH1D* h6 = new TH1D("h6", "Scattering Angle(#theta) with CUT", 100, 130.0, 180.0);
+  //TH1D* h7 = new TH1D("h7", "Pseudo-rapidity(#eta) with CUT", 100, -5.0, 0.0);
+  //TH2D *h8  = new TH2D("h8","Cluster Hit Position",                        62,-62.0,62.0,62,-62.0,62.0);
+  //TH2D *h9  = new TH2D("h9","All Hit Position",                            62,-62.0,62.0,62,-62.0,62.0);
+  //TH1D *h10 = new TH1D("hEnergyRes","Energy Resolution",                   100,-0.3,0.3);
+  //TH1D *h11 = new TH1D("hEnergyResCUT","Energy Resolution with CUT",       100,-0.3,0.3);
+  //TH1D *h12 = new TH1D("h12","Thrown momentum",                            61,e_start-0.5,e_end+0.5);
+  //TH1D *h13 = new TH1D("h13","Thrown momentum for reconstructed particle", 61,e_start-0.5,e_end+0.5);
+  //TH1D *h14 = new TH1D("h14","Ratio p_{rec}/p_{thr}",                      61,e_start-0.5,e_end+0.5);
+
+  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("Ecluster", [](const std::vector<eic::ClusterData>& in) {
+                  std::vector<float> res;
+                  for (const auto& i : in)
+                    res.push_back(i.energy);
+                  return res;
+                },{"EcalClusters"})
+                .Define("nclusters","EcalClusters.size()") ;
+
+  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,30},"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_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/crystal_cal_electrons_etaThrown.png");
+
+  h_theta_thrown->DrawCopy();
+  c->SaveAs("results/crystal_cal_electrons_thetaThrown.png");
+
+  h_nclusters->DrawCopy();
+  c->SaveAs("results/crystal_cal_electrons_nclusters.png");
+
+  h_Ecluster->DrawCopy();
+  h_Ecluster1->SetLineColor(2);
+  h_Ecluster1->DrawCopy("same");
+  //h_Ecluster2->SetLineColor(4);
+  //h_Ecluster2->DrawCopy("same");
+  c->SaveAs("results/crystal_cal_electrons_Ecluster.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;
+}

ecal/scripts/rdf_test.cxx

+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+
+//namespace edm4hep {
+//
+//std::vector<float> pt (std::vector<MCParticleData> const& in){
+//  std::vector<float> result;
+//  for (size_t i = 0; i < in.size(); ++i) {
+//    result.push_back(std::sqrt(in[i].momentum.x * in[i].momentum.x + in[i].momentum.y * in[i].momentum.y));
+//  }
+//  return result;
+//}
+//
+//std::vector<float> eta(std::vector<MCParticleData> const& in){
+//  std::vector<float> result;
+//  ROOT::Math::PxPyPzMVector lv;
+//  for (size_t i = 0; i < in.size(); ++i) {
+//    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
+//    result.push_back(lv.Eta());
+//  }
+//  return result;
+//}
+//
+//std::vector<float> cos_theta(std::vector<MCParticleData> const& in){
+//  std::vector<float> result;
+//  ROOT::Math::PxPyPzMVector lv;
+//  for (size_t i = 0; i < in.size(); ++i) {
+//    lv.SetCoordinates(in[i].momentum.x, in[i].momentum.y, in[i].momentum.z, in[i].mass);
+//    result.push_back(cos(lv.Theta()));
+//  }
+//  return result;
+//}
+//
+//}
+
+
+std::vector<float> 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;
+}
+
+int rdf_test() {
+
+  ROOT::EnableImplicitMT();
+
+  ROOT::RDataFrame df("events", "sim_barrel_clusters.root");
+
+  auto df2 =
+      df.Define("MCParticles_pt", pt, {"mcparticles"});
+          //.Define("MCParticles_eta", edm4hep::eta, {"mcparticles"})
+          //.Define("MCParticles_cosTheta", edm4hep::cos_theta, {"mcparticles"});
+
+  std::string outfilename = "rdf_test.root";
+
+  df2.Snapshot("events", outfilename,
+               {
+                   "MCParticles_pt",
+                   "mcparticles"
+               });
+  return 0;
+}