Skip to content
Snippets Groups Projects
Commit f3299de6 authored by Whitney Armstrong's avatar Whitney Armstrong
Browse files

deleted: scripts/gen_central_particles.cxx

	deleted:    scripts/hits_central_electrons.cxx
	deleted:    scripts/hits_central_pions.cxx
parent 7fe0ac30
No related branches found
No related tags found
No related merge requests found
#include "HepMC3/GenEvent.h"
#include "HepMC3/ReaderAscii.h"
#include "HepMC3/WriterAscii.h"
#include "HepMC3/Print.h"
#include <iostream>
#include <random>
#include <cmath>
#include "TMath.h"
#include "common_bench/particles.h"
using namespace HepMC3;
/** Generate particles in the central region.
*/
void gen_central_particles(int n_events = 100,
const char* out_fname = "central_electrons.hepmc")
{
double cos_theta_min = std::cos( 10.0*(M_PI/180.0));
double cos_theta_max = std::cos(170.0*(M_PI/180.0));
const auto& partmap = common_bench::particleMap;
WriterAscii hepmc_output(out_fname);
GenEvent evt(Units::GEV, Units::MM);
int events_parsed = 0;
// Random number generator
TRandom *r1 = new TRandom();
for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
// FourVector(px,py,pz,e,pdgid,status)
// type 4 is beam
// pdgid 11 - electron
// pdgid 111 - pi0
// pdgid 2212 - proton
GenParticlePtr p1 =
std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 11, 4);
GenParticlePtr p2 = std::make_shared<GenParticle>(
FourVector(0.0, 0.0, 0.0, partmap.at(2212).mass), 2212, 4);
// Define momentum
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 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);
// Generates random vectors, uniformly distributed over the surface of a
// sphere of given radius, in this case momentum.
// r1->Sphere(px, py, pz, p);
//std::cout << std::sqrt(px*px + py*py + pz*pz) - p << " is zero? \n";
// type 1 is final state
// pdgid 11 - electron 0.510 MeV/c^2
GenParticlePtr p3 = std::make_shared<GenParticle>(
FourVector(
px, py, pz,
sqrt(p*p + std::pow(partmap.at(11).mass,2))),
11, 1);
GenVertexPtr v1 = std::make_shared<GenVertex>();
v1->add_particle_in(p1);
v1->add_particle_in(p2);
v1->add_particle_out(p3);
evt.add_vertex(v1);
if (events_parsed == 0) {
std::cout << "First event: " << std::endl;
Print::listing(evt);
}
hepmc_output.write_event(evt);
if (events_parsed % 10000 == 0) {
std::cout << "Event: " << events_parsed << std::endl;
}
evt.clear();
}
hepmc_output.close();
std::cout << "Events parsed and written: " << events_parsed << std::endl;
}
#include "ROOT/RDataFrame.hxx"
#include "TCanvas.h"
#include "TLegend.h"
#include "TH1D.h"
#include "TProfile.h"
#include <iostream>
R__LOAD_LIBRARY(libeicd.so)
R__LOAD_LIBRARY(libDD4pod.so)
#include "dd4pod/TrackerHitCollection.h"
#include "dd4pod/TrackerHitData.h"
#include "dd4pod/Geant4ParticleCollection.h"
#include "eicd/TrackParametersCollection.h"
#include "eicd/ClusterCollection.h"
#include "eicd/ClusterData.h"
#include "eicd/TrackerHitCollection.h"
using ROOT::RDataFrame;
using namespace ROOT::VecOps;
auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
std::vector<double> result;
for (size_t i = 0; i < in.size(); ++i) {
result.push_back(std::abs(1.0/(in[i].qOverP)));
}
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].ps.x * in[i].ps.x + in[i].ps.y * in[i].ps.y));
}
return result;
}
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].ps.x, in[i].ps.y, in[i].ps.z, in[i].mass);
result.push_back(lv);
}
return result;
};
auto delta_p = [](const std::vector<double>& tracks, const std::vector<double>& thrown) {
std::vector<double> res;
for (const auto& p1 : thrown) {
for (const auto& p2 : tracks) {
res.push_back(p1 - p2);
}
}
return res;
};
auto delta_p_over_p = [](const std::vector<double>& tracks, const std::vector<double>& thrown) {
std::vector<double> res;
for (const auto& p1 : thrown) {
for (const auto& p2 : tracks) {
res.push_back((p1 - p2)/p1);
}
}
return res;
};
int hits_central_electrons(const char* fname = "sim_central_electrons.root")
{
ROOT::EnableImplicitMT();
ROOT::RDataFrame df("events", fname);
auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
.Define("thrownParticles", "mcparticles[isThrown]")
.Define("thrownP", fourvec, {"thrownParticles"})
.Define("p_thrown", momentum, {"thrownP"})
.Define("theta_thrown", theta, {"thrownP"})
.Define("theta0", "theta_thrown[0]")
//.Define("nTracks", "outputTrackParameters.size()")
//.Define("p_track", p_track, {"outputTrackParameters"})
//.Define("p_track1", p_track, {"outputTrackParameters1"})
//.Define("p_track2", p_track, {"outputTrackParameters2"})
//.Define("delta_p0",delta_p, {"p_track", "p_thrown"})
//.Define("delta_p1",delta_p, {"p_track1", "p_thrown"})
////.Define("delta_p2",delta_p, {"p_track2", "p_thrown"})
//.Define("delta_p_over_p0",delta_p_over_p, {"p_track", "p_thrown"})
//.Define("delta_p_over_p1",delta_p_over_p, {"p_track1", "p_thrown"})
//.Define("delta_p_over_p2",delta_p_over_p, {"p_track2", "p_thrown"})
//.Define("N_VtxBarrelHits",[](std::vector<eic::TrackerHitData> hits) { return hits.size();},{"VertexBarrelRecHits"})
//.Define("N_Hits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"trackingHits"})
.Define("N_BarrelHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerBarrelHits"})
.Define("N_EndcapHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerEndcapHits"})
;
auto hBarrel_x_vs_y = df0.Histo2D({"hBarrel_x_vs_y", "; x ; y ", 100, -900, 900,100, -900, 900 }, "TrackerBarrelHits.position.x", "TrackerBarrelHits.position.y");
auto hEndcap_x_vs_y = df0.Histo2D({"hEndcap_x_vs_y", "; x ; y ", 100, -900, 900,100, -900, 900 }, "TrackerEndcaplHits.position.x", "TrackerEndcapHits.position.y");
//auto hvtxBarrel_x_vs_y = df0.Histo2D({"hvtxBarrel_x_vs_y", "; x ; y ", 100, -900, 900,100, -900, 900 }, "VertexBarrelHits.position.x", "VertexBarrelHits.position.y");
//auto hvtxEndcap_x_vs_y = df0.Histo2D({"hvtxEndcap_x_vs_y", "; x ; y ", 100, -900, 900,100, -900, 900 }, "VertexEndcaplHits.position.x", "VertexEndcapHits.position.y");
//auto hAllHits_x_vs_y = df0.Histo2D({"hAllHitsx_vs_y", "; x ; y ", 100, -900, 900,100, -900, 900 }, "allHits.position.x", "allHits.position.y");
auto hBarrel_x_vs_z = df0.Histo2D({"hBarrel_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "TrackerBarrelHits.position.z" , "TrackerBarrelHits.position.y");
auto hEndcap_x_vs_z = df0.Histo2D({"hEndcap_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "TrackerEndcaplHits.position.z", "TrackerEndcapHits.position.y");
//auto hvtxBarrel_x_vs_z = df0.Histo2D({"hvtxBarrel_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "VertexBarrelHits.position.z" , "VertexBarrelHits.position.y" );
//auto hvtxEndcap_x_vs_z = df0.Histo2D({"hvtxEndcap_x_vs_z", "; z ; x ", 100, -900, 900,100, -900, 900 }, "VertexEndcaplHits.position.z" , "VertexEndcapHits.position.y" );
//auto hAllHits_x_vs_z = df0.Histo2D({"hAllHitsx_vs_z","; z ; x ", 100, -900, 900,100, -900, 900 }, "allHits.position.z" , "allHits.position.y" );
auto hBarrel_N_vs_theta = df0.Histo1D({"hBarrel_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_BarrelHits");
auto hEndcap_N_vs_theta = df0.Histo1D({"hEndcap_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_EndcapHits");
//auto hVtxBarrel_N_vs_theta = df0.Histo1D({"hVtxBarrel_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_VtxBarrelHits");
auto hBarrel_Nhits = df0.Histo1D({"hBarrel_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_BarrelHits");
auto hEndcap_Nhits = df0.Histo1D({"hEndcap_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_EndcapHits");
//auto hVtxBarrel_Nhits = df0.Histo1D({"hVtxBarrel_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_VtxBarrelHits");
auto hBarrel_Ntheta = df0.Histo1D({"hBarrel_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
auto hEndcap_Ntheta = df0.Histo1D({"hEndcap_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
//auto hVtxBarrel_Ntheta = df0.Histo1D({"hVtxBarrel_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
auto c = new TCanvas();
auto hs = new THStack("n_hits","; #theta ");
auto h1 = (TH1D*) hBarrel_N_vs_theta->Clone();
auto h2 = (TH1D*) hBarrel_Ntheta->Clone();
h1->Divide(h2);
hs->Add(h1);
h1 = (TH1D*) hEndcap_N_vs_theta->Clone();
h2 = (TH1D*) hEndcap_Ntheta->Clone();
h1->Divide(h2);
h1->SetLineColor(2);
hs->Add(h1);
//h1 = (TH1D*) hVtxBarrel_vs_theta->Clone();
//h1->SetLineColor(4);
//h1->SetFillStyle(3001);
//h1->SetFillColor(4);
//hs->Add(h1);
hs->Draw("nostack, hist");
c->BuildLegend();
c->SaveAs("results/tracking/hits_central_electrons_n_hits_vs_theta.png");
c->SaveAs("results/tracking/hits_central_electrons_n_hits_vs_theta.pdf");
c = new TCanvas();
hs = new THStack("theta","; #theta ");
h1 = (TH1D*) hBarrel_N_vs_theta->Clone();
h2 = (TH1D*) hBarrel_Ntheta->Clone();
//h1->Divide(h2);
hs->Add(h2);
h1 = (TH1D*) hEndcap_N_vs_theta->Clone();
h2 = (TH1D*) hEndcap_Ntheta->Clone();
//h1->Divide(h2);
h1->SetLineColor(2);
h2->SetLineColor(2);
hs->Add(h2);
//h1 = (TH1D*) hVtxBarrel_vs_theta->Clone();
//h1->SetLineColor(4);
//h1->SetFillStyle(3001);
//h1->SetFillColor(4);
//hs->Add(h1);
hs->Draw("nostack hist");
c->BuildLegend();
c->SaveAs("results/tracking/hits_central_electrons_theta.png");
c->SaveAs("results/tracking/hits_central_electrons_theta.pdf");
c = new TCanvas();
hs = new THStack("hits","; hits ");
h1 = (TH1D*) hBarrel_Nhits->Clone();
hs->Add(h1);
h1 = (TH1D*) hEndcap_Nhits->Clone();
h1->SetLineColor(2);
h2->SetLineColor(2);
hs->Add(h2);
//h1 = (TH1D*) hVtxBarrel_Nhits->Clone();
//h1->SetLineColor(4);
//h1->SetFillStyle(3001);
//h1->SetFillColor(4);
//hs->Add(h1);
//hs->Draw("nostack hist");
c->BuildLegend();
c->SaveAs("results/tracking/hits_central_electrons_nhits.png");
c->SaveAs("results/tracking/hits_central_electrons_nhits.pdf");
c = new TCanvas();
hBarrel_x_vs_y->DrawCopy("colz");
c->SaveAs("results/tracking/hits_central_electrons_trkBarrel_xy.png");
c->SaveAs("results/tracking/hits_central_electrons_trkBarrel_xy.pdf");
c = new TCanvas();
hBarrel_x_vs_y->DrawCopy("colz");
hEndcap_x_vs_y->DrawCopy("colz same");
//hvtxBarrel_x_vs_y->DrawCopy("colz same");
//hvtxEndcap_x_vs_y->DrawCopy("colz same");
c->SaveAs("results/tracking/hits_central_electrons_Hits_xy.png");
c->SaveAs("results/tracking/hits_central_electrons_Hits_xy.pdf");
//hAllHits_x_vs_z->DrawCopy("colz");
hBarrel_x_vs_z->DrawCopy("colz");
hEndcap_x_vs_z->DrawCopy("colz same");
c->SaveAs("results/tracking/hits_central_electrons_Hits_xz.png");
c->SaveAs("results/tracking/hits_central_electrons_Hits_xz.pdf");
return 0;
}
#include "ROOT/RDataFrame.hxx"
#include "TCanvas.h"
#include "TLegend.h"
#include "TH1D.h"
#include "TProfile.h"
#include <iostream>
R__LOAD_LIBRARY(libeicd.so)
R__LOAD_LIBRARY(libDD4pod.so)
#include "dd4pod/TrackerHitCollection.h"
#include "dd4pod/TrackerHitData.h"
#include "dd4pod/Geant4ParticleCollection.h"
#include "eicd/TrackParametersCollection.h"
#include "eicd/ClusterCollection.h"
#include "eicd/ClusterData.h"
#include "eicd/TrackerHitCollection.h"
using ROOT::RDataFrame;
using namespace ROOT::VecOps;
auto p_track = [](std::vector<eic::TrackParametersData> const& in) {
std::vector<double> result;
for (size_t i = 0; i < in.size(); ++i) {
result.push_back(std::abs(1.0/(in[i].qOverP)));
}
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].ps.x * in[i].ps.x + in[i].ps.y * in[i].ps.y));
}
return result;
}
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].ps.x, in[i].ps.y, in[i].ps.z, in[i].mass);
result.push_back(lv);
}
return result;
};
auto delta_p = [](const std::vector<double>& tracks, const std::vector<double>& thrown) {
std::vector<double> res;
for (const auto& p1 : thrown) {
for (const auto& p2 : tracks) {
res.push_back(p1 - p2);
}
}
return res;
};
auto delta_p_over_p = [](const std::vector<double>& tracks, const std::vector<double>& thrown) {
std::vector<double> res;
for (const auto& p1 : thrown) {
for (const auto& p2 : tracks) {
res.push_back((p1 - p2)/p1);
}
}
return res;
};
int hits_central_pions(const char* fname = "sim_central_pions.root")
{
ROOT::EnableImplicitMT();
ROOT::RDataFrame df("events", fname);
auto df0 = df.Define("isThrown", "mcparticles.genStatus == 1")
.Define("thrownParticles", "mcparticles[isThrown]")
.Define("thrownP", fourvec, {"thrownParticles"})
.Define("p_thrown", momentum, {"thrownP"})
.Define("theta_thrown", theta, {"thrownP"})
.Define("theta0", "theta_thrown[0]")
//.Define("nTracks", "outputTrackParameters.size()")
//.Define("p_track", p_track, {"outputTrackParameters"})
//.Define("p_track1", p_track, {"outputTrackParameters1"})
//.Define("p_track2", p_track, {"outputTrackParameters2"})
//.Define("delta_p0",delta_p, {"p_track", "p_thrown"})
//.Define("delta_p1",delta_p, {"p_track1", "p_thrown"})
////.Define("delta_p2",delta_p, {"p_track2", "p_thrown"})
//.Define("delta_p_over_p0",delta_p_over_p, {"p_track", "p_thrown"})
//.Define("delta_p_over_p1",delta_p_over_p, {"p_track1", "p_thrown"})
//.Define("delta_p_over_p2",delta_p_over_p, {"p_track2", "p_thrown"})
//.Define("N_VtxBarrelHits",[](std::vector<eic::TrackerHitData> hits) { return hits.size();},{"VertexBarrelRecHits"})
.Define("N_BarrelHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerBarrelHits"})
.Define("N_EndcapHits", [](std::vector<dd4pod::TrackerHitData> hits) { return hits.size();}, {"TrackerEndcapHits"})
;
auto hBarrel_x_vs_y = df0.Histo2D({"hBarrel_x_vs_y", "; x ; y ", 100, -900, 900,100, -900, 900 }, "TrackerBarrelHits.position.x", "TrackerBarrelHits.position.y");
auto hBarrel_N_vs_theta = df0.Histo1D({"hBarrel_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_BarrelHits");
auto hEndcap_N_vs_theta = df0.Histo1D({"hEndcap_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_EndcapHits");
//auto hVtxBarrel_N_vs_theta = df0.Histo1D({"hVtxBarrel_N_vs_theta", "; #theta [deg.]", 20, 0, 180 }, "theta0", "N_VtxBarrelHits");
auto hBarrel_Nhits = df0.Histo1D({"hBarrel_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_BarrelHits");
auto hEndcap_Nhits = df0.Histo1D({"hEndcap_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_EndcapHits");
//auto hVtxBarrel_Nhits = df0.Histo1D({"hVtxBarrel_Nhits", "; #theta [deg.]", 20, 0, 20 }, "N_VtxBarrelHits");
auto hBarrel_Ntheta = df0.Histo1D({"hBarrel_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
auto hEndcap_Ntheta = df0.Histo1D({"hEndcap_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
//auto hVtxBarrel_Ntheta = df0.Histo1D({"hVtxBarrel_Ntheta", "; #theta [deg.]", 20, 0, 180 }, "theta0");
auto c = new TCanvas();
auto hs = new THStack("n_hits","; #theta ");
auto h1 = (TH1D*) hBarrel_N_vs_theta->Clone();
auto h2 = (TH1D*) hBarrel_Ntheta->Clone();
h1->Divide(h2);
hs->Add(h1);
h1 = (TH1D*) hEndcap_N_vs_theta->Clone();
h2 = (TH1D*) hEndcap_Ntheta->Clone();
h1->Divide(h2);
h1->SetLineColor(2);
hs->Add(h1);
//h1 = (TH1D*) hVtxBarrel_vs_theta->Clone();
//h1->SetLineColor(4);
//h1->SetFillStyle(3001);
//h1->SetFillColor(4);
//hs->Add(h1);
hs->Draw("nostack, hist");
c->BuildLegend();
c->SaveAs("results/tracking/hits_central_pions_n_hits_vs_theta.png");
c->SaveAs("results/tracking/hits_central_pions_n_hits_vs_theta.pdf");
c = new TCanvas();
hs = new THStack("theta","; #theta ");
h1 = (TH1D*) hBarrel_N_vs_theta->Clone();
h2 = (TH1D*) hBarrel_Ntheta->Clone();
//h1->Divide(h2);
hs->Add(h2);
h1 = (TH1D*) hEndcap_N_vs_theta->Clone();
h2 = (TH1D*) hEndcap_Ntheta->Clone();
//h1->Divide(h2);
h1->SetLineColor(2);
h2->SetLineColor(2);
hs->Add(h2);
//h1 = (TH1D*) hVtxBarrel_vs_theta->Clone();
//h1->SetLineColor(4);
//h1->SetFillStyle(3001);
//h1->SetFillColor(4);
//hs->Add(h1);
hs->Draw("nostack hist");
c->BuildLegend();
c->SaveAs("results/tracking/hits_central_pions_theta.png");
c->SaveAs("results/tracking/hits_central_pions_theta.pdf");
c = new TCanvas();
hs = new THStack("hits","; hits ");
h1 = (TH1D*) hBarrel_Nhits->Clone();
hs->Add(h1);
h1 = (TH1D*) hEndcap_Nhits->Clone();
h1->SetLineColor(2);
h2->SetLineColor(2);
hs->Add(h2);
//h1 = (TH1D*) hVtxBarrel_Nhits->Clone();
//h1->SetLineColor(4);
//h1->SetFillStyle(3001);
//h1->SetFillColor(4);
//hs->Add(h1);
//hs->Draw("nostack hist");
c->BuildLegend();
c->SaveAs("results/tracking/hits_central_pions_nhits.png");
c->SaveAs("results/tracking/hits_central_pions_nhits.pdf");
c = new TCanvas();
hBarrel_x_vs_y->DrawCopy("colz");
c->SaveAs("results/tracking/hits_central_pions_xy.png");
c->SaveAs("results/tracking/hits_central_pions_xy.pdf");
return 0;
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment