diff --git a/bin/crop_pdf b/bin/crop_pdf new file mode 100755 index 0000000000000000000000000000000000000000..488871761f598c6c661499204a23c37479a3990c --- /dev/null +++ b/bin/crop_pdf @@ -0,0 +1,28 @@ +#!/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 + diff --git a/bin/crop_views b/bin/crop_views new file mode 100644 index 0000000000000000000000000000000000000000..8d5191233723b37be2ab8e96d12384e9690c6424 --- /dev/null +++ b/bin/crop_views @@ -0,0 +1,35 @@ +#!/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 + diff --git a/bin/download_views b/bin/download_views index ca90dec83b13b4c0f4c886b95299ee44c5dde862..4dc90790f83a8f2defb35300efb3413689bb4cf7 100644 --- a/bin/download_views +++ b/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' diff --git a/ecal/emcal_electrons.sh b/ecal/emcal_electrons.sh index ecaa1f652b2646d1957b73d5e45487e71f729432..9cfc1f3095e2986930f4dc91dcd5b4b1e978a5b9 100644 --- a/ecal/emcal_electrons.sh +++ b/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\")" diff --git a/ecal/emcal_pi0s.sh b/ecal/emcal_pi0s.sh index 2387fbf4811194eaedc8156f81f9644006f5deaf..60ef1a1b65e79403945fa5a7c11ca3389714fba2 100644 --- a/ecal/emcal_pi0s.sh +++ b/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 \ diff --git a/ecal/scripts/crystal_cal_electrons.cxx b/ecal/scripts/crystal_cal_electrons.cxx new file mode 100644 index 0000000000000000000000000000000000000000..f5d781eade75ad4304a1cd2c10e125cde2e3e1df --- /dev/null +++ b/ecal/scripts/crystal_cal_electrons.cxx @@ -0,0 +1,130 @@ +#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; +} diff --git a/ecal/scripts/rdf_test.cxx b/ecal/scripts/rdf_test.cxx new file mode 100644 index 0000000000000000000000000000000000000000..7dc3b14e6138ba526b859a00981c6541d0bba311 --- /dev/null +++ b/ecal/scripts/rdf_test.cxx @@ -0,0 +1,66 @@ +#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; +}