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

Adding comparison for conformal finding

parent 04ffa426
Branches
No related tags found
1 merge request!179Adding comparison for conformal finding
...@@ -38,6 +38,7 @@ from Configurables import Jug__Reco__TrackParamTruthInit as TrackParamTruthInit ...@@ -38,6 +38,7 @@ from Configurables import Jug__Reco__TrackParamTruthInit as TrackParamTruthInit
from Configurables import Jug__Reco__TrackParamClusterInit as TrackParamClusterInit from Configurables import Jug__Reco__TrackParamClusterInit as TrackParamClusterInit
from Configurables import Jug__Reco__TrackParamVertexClusterInit as TrackParamVertexClusterInit from Configurables import Jug__Reco__TrackParamVertexClusterInit as TrackParamVertexClusterInit
from Configurables import Jug__Reco__ConformalXYPeakProtoTracks as ConformalXYPeakProtoTracks
from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm from Configurables import Jug__Reco__TrackFindingAlgorithm as TrackFindingAlgorithm
from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit from Configurables import Jug__Reco__ParticlesFromTrackFit as ParticlesFromTrackFit
from Configurables import Jug__Reco__EMCalReconstruction as EMCalReconstruction from Configurables import Jug__Reco__EMCalReconstruction as EMCalReconstruction
...@@ -48,7 +49,7 @@ from Configurables import Jug__Reco__SimpleClustering as SimpleClustering ...@@ -48,7 +49,7 @@ from Configurables import Jug__Reco__SimpleClustering as SimpleClustering
algorithms = [ ] algorithms = [ ]
podioinput = PodioInput("PodioReader", podioinput = PodioInput("PodioReader",
collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","GEMTrackerEndcapHits"])#, OutputLevel=DEBUG) collections=["mcparticles","TrackerEndcapHits","TrackerBarrelHits","VertexBarrelHits","VertexEndcapHits","GEMTrackerEndcapHits"])
algorithms.append( podioinput ) algorithms.append( podioinput )
## copiers to get around input --> output copy bug. Note the "2" appended to the output collection. ## copiers to get around input --> output copy bug. Note the "2" appended to the output collection.
...@@ -128,6 +129,14 @@ trk_hit_col = TrackingHitsCollector("trk_hit_col", ...@@ -128,6 +129,14 @@ trk_hit_col = TrackingHitsCollector("trk_hit_col",
OutputLevel=DEBUG) OutputLevel=DEBUG)
algorithms.append( trk_hit_col ) algorithms.append( trk_hit_col )
# Hit Source linker
conformal_find = ConformalXYPeakProtoTracks("conformal_find",
inputTrackerHits=trk_hit_col.trackingHits,
outputProtoTracks="outputProtoTracks",
nProtoTracks="nProtoTracks",
OutputLevel=DEBUG)
algorithms.append( conformal_find )
# Hit Source linker # Hit Source linker
sourcelinker = TrackerSourceLinker("sourcelinker", sourcelinker = TrackerSourceLinker("sourcelinker",
inputHitCollection=trk_hit_col.trackingHits, inputHitCollection=trk_hit_col.trackingHits,
......
#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;
}
...@@ -19,7 +19,7 @@ using namespace HepMC3; ...@@ -19,7 +19,7 @@ using namespace HepMC3;
*/ */
void gen_multiple_tracks(int n_events = 100, void gen_multiple_tracks(int n_events = 100,
const char* out_fname = "multiple_tracks.hepmc", const char* out_fname = "multiple_tracks.hepmc",
int n_parts = 2) int n_parts = 5)
{ {
double cos_theta_min = std::cos( 10.0*(M_PI/180.0)); double cos_theta_min = std::cos( 10.0*(M_PI/180.0));
double cos_theta_max = std::cos(170.0*(M_PI/180.0)); double cos_theta_max = std::cos(170.0*(M_PI/180.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_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;
}
...@@ -87,6 +87,7 @@ int rec_multiple_tracks(const char* fname = "topside/rec_multiple_tracks.root") ...@@ -87,6 +87,7 @@ int rec_multiple_tracks(const char* fname = "topside/rec_multiple_tracks.root")
auto df0 = df.Define("isThrown", "mcparticles2.genStatus == 1") auto df0 = df.Define("isThrown", "mcparticles2.genStatus == 1")
.Define("thrownParticles", "mcparticles2[isThrown]") .Define("thrownParticles", "mcparticles2[isThrown]")
.Define("thrownP", fourvec, {"thrownParticles"}) .Define("thrownP", fourvec, {"thrownParticles"})
.Define("nThrown", "thrownParticles.size()")
.Define("p_thrown", momentum, {"thrownP"}) .Define("p_thrown", momentum, {"thrownP"})
.Define("theta_thrown", theta, {"thrownP"}) .Define("theta_thrown", theta, {"thrownP"})
.Define("theta0", "theta_thrown[0]") .Define("theta0", "theta_thrown[0]")
...@@ -106,6 +107,9 @@ int rec_multiple_tracks(const char* fname = "topside/rec_multiple_tracks.root") ...@@ -106,6 +107,9 @@ int rec_multiple_tracks(const char* fname = "topside/rec_multiple_tracks.root")
auto h_delta_p0 = df0.Histo1D({"h_delta_p0", "Truth Track Init; GeV/c ", 100, -10, 10}, "delta_p0"); auto h_delta_p0 = df0.Histo1D({"h_delta_p0", "Truth Track Init; GeV/c ", 100, -10, 10}, "delta_p0");
auto h_nProtoTracks = df0.Histo1D({"h_nProtoTracks", "; n ", 10, 0, 10}, "nProtoTracks");
auto h_nThrown = df0.Histo1D({"h_nThrown", "; n ", 10, 0, 10}, "nThrown");
auto h_delta_p0_over_p = df0.Histo1D({"h_delta_p0_over_p", "Truth Track Init; delta p/p ", 100, -0.1, 0.1}, "delta_p_over_p0"); auto h_delta_p0_over_p = df0.Histo1D({"h_delta_p0_over_p", "Truth Track Init; delta p/p ", 100, -0.1, 0.1}, "delta_p_over_p0");
auto hNhits_vs_theta = df0.Histo1D({"hNhits_vs_theta", "; #theta [deg.]", 40, 0, 180 }, "theta0", "N_Hits"); auto hNhits_vs_theta = df0.Histo1D({"hNhits_vs_theta", "; #theta [deg.]", 40, 0, 180 }, "theta0", "N_Hits");
...@@ -123,6 +127,12 @@ int rec_multiple_tracks(const char* fname = "topside/rec_multiple_tracks.root") ...@@ -123,6 +127,12 @@ int rec_multiple_tracks(const char* fname = "topside/rec_multiple_tracks.root")
// ----------------------------------------------- // -----------------------------------------------
auto c = new TCanvas(); auto c = new TCanvas();
h_nThrown->SetLineColor(2);
h_nThrown->DrawCopy();
h_nProtoTracks->DrawCopy("same");
c->SaveAs("results/track_finding/rec_multiple_tracks_nProtoTracks.png");
c->SaveAs("results/track_finding/rec_multiple_tracks_nProtoTracks.pdf");
h_nTracks->DrawCopy(); h_nTracks->DrawCopy();
c->SaveAs("results/track_finding/rec_multiple_tracks_nTracks.png"); c->SaveAs("results/track_finding/rec_multiple_tracks_nTracks.png");
c->SaveAs("results/track_finding/rec_multiple_tracks_nTracks.pdf"); c->SaveAs("results/track_finding/rec_multiple_tracks_nTracks.pdf");
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment