diff --git a/TrackPerf/EfficiencyHists.hxx b/TrackPerf/EfficiencyHists.hxx index 990f5d3..1a18332 100644 --- a/TrackPerf/EfficiencyHists.hxx +++ b/TrackPerf/EfficiencyHists.hxx @@ -5,6 +5,10 @@ #include #include +// DD4hep +#include +#include + namespace EVENT { class Track; class MCParticle; @@ -22,12 +26,10 @@ class EfficiencyHists { // Fill histograms with a single track void fillMC(const EVENT::MCParticle* track, bool passed); - void fillTrack(const EVENT::Track* track, bool passed); + void fillTrack(const EVENT::Track* track, bool passed, + dd4hep::Detector* lcdd); private: - //! magnetic field to use for curvature -> pT conversion - float _Bz = 3.57; - //! Efficiency plots TEfficiency* h_effpt; TEfficiency* h_efftheta; diff --git a/TrackPerf/ResoHists.hxx b/TrackPerf/ResoHists.hxx index 8db26a1..d95cde3 100644 --- a/TrackPerf/ResoHists.hxx +++ b/TrackPerf/ResoHists.hxx @@ -4,6 +4,12 @@ #include #include +// DD4hep +#include +#include + + + namespace EVENT { class Track; class MCParticle; @@ -20,12 +26,11 @@ class ResoHists { ResoHists(); // Fill histograms with a single track - void fill(const EVENT::Track* track, const EVENT::MCParticle* particle); + void fill(const EVENT::Track* track, + const EVENT::MCParticle* particle, + dd4hep::Detector* lcdd); private: - //! magnetic field to use for curvature -> pT conversion - float _Bz = 3.57; - //! Reconstructed track pT TH2* h_track_truth_pt; TH1* h_reso_pt_rel; diff --git a/TrackPerf/TrackHists.hxx b/TrackPerf/TrackHists.hxx index 5a74dd4..79aa099 100644 --- a/TrackPerf/TrackHists.hxx +++ b/TrackPerf/TrackHists.hxx @@ -3,6 +3,10 @@ #include #include +// DD4hep +#include +#include + namespace EVENT { class Track; } @@ -18,12 +22,10 @@ class TrackHists { TrackHists(); // Fill histograms with a single track - void fill(const EVENT::Track* track); + void fill(const EVENT::Track* track, + dd4hep::Detector* lcdd); private: - //! magnetic field to use for curvature -> pT conversion - float _Bz = 3.57; - //! Reconstructed track pT TH1* h_pt; TH1* h_lambda; diff --git a/TrackPerf/TrackPerfHistProc.hxx b/TrackPerf/TrackPerfHistProc.hxx index 3f85ce1..6328f9e 100644 --- a/TrackPerf/TrackPerfHistProc.hxx +++ b/TrackPerf/TrackPerfHistProc.hxx @@ -4,6 +4,9 @@ #include +// DD4hep +#include + namespace TrackPerf { class TrackHists; class TruthHists; @@ -49,6 +52,10 @@ class TrackPerfHistProc : public marlin::Processor { //! Track to MC truth match collection std::string _trkMatchColName{}; + //! Magnetic field to use for curvature -> pT conversion + dd4hep::Detector* _lcdd; + + //! Determination of good vs bad match float _matchProb = 0.5; diff --git a/src/EfficiencyHists.cxx b/src/EfficiencyHists.cxx index 6af1f9c..1344baa 100644 --- a/src/EfficiencyHists.cxx +++ b/src/EfficiencyHists.cxx @@ -17,8 +17,16 @@ EfficiencyHists::EfficiencyHists(bool effi) { } } -void EfficiencyHists::fillTrack(const EVENT::Track* track, bool passed) { - float pt = fabs(0.3 * _Bz / track->getOmega() / 1000); +void EfficiencyHists::fillTrack(const EVENT::Track* track, bool passed, + dd4hep::Detector* lcdd) { + //TODO: This assumes uniform magnetic field + const double position[3] = {0, 0, 0}; // position to calculate magnetic field (here, the origin) + double magneticFieldVector[3] = {0, 0, 0}; // initialise object to hold magnetic field + lcdd->field().magneticField( + position, magneticFieldVector); // get the magnetic field vector from DD4hep + float Bz = magneticFieldVector[2] / dd4hep::tesla; + + float pt = fabs(0.3 * Bz / track->getOmega() / 1000); float theta = TMath::Pi() - std::atan(track->getTanLambda()); h_effpt->Fill(passed, pt); diff --git a/src/ResoHists.cxx b/src/ResoHists.cxx index b9bfad8..fd8a9ea 100644 --- a/src/ResoHists.cxx +++ b/src/ResoHists.cxx @@ -23,8 +23,16 @@ ResoHists::ResoHists() { } void ResoHists::fill(const EVENT::Track* track, - const EVENT::MCParticle* particle) { - float track_pt = fabs(0.3 * _Bz / track->getOmega() / 1000); + const EVENT::MCParticle* particle, + dd4hep::Detector* lcdd) { + //TODO: This assumes uniform magnetic field + const double position[3] = {0, 0, 0}; // position to calculate magnetic field (here, the origin) + double magneticFieldVector[3] = {0, 0, 0}; // initialise object to hold magnetic field + lcdd->field().magneticField( + position, magneticFieldVector); // get the magnetic field vector from DD4hep + float Bz = magneticFieldVector[2] / dd4hep::tesla; + + float track_pt = fabs(0.3 * Bz / track->getOmega() / 1000); float track_lambda = std::atan(track->getTanLambda()); const double* mom = particle->getMomentum(); diff --git a/src/TrackHists.cxx b/src/TrackHists.cxx index 4332fae..4320593 100644 --- a/src/TrackHists.cxx +++ b/src/TrackHists.cxx @@ -32,8 +32,16 @@ TrackHists::TrackHists() { 20, -0.5, 19.5); } -void TrackHists::fill(const EVENT::Track* track) { - float pt = fabs(0.3 * _Bz / track->getOmega() / 1000); +void TrackHists::fill(const EVENT::Track* track, + dd4hep::Detector* lcdd) { + //TODO: This assumes uniform magnetic field + const double position[3] = {0, 0, 0}; // position to calculate magnetic field (here, the origin) + double magneticFieldVector[3] = {0, 0, 0}; // initialise object to hold magnetic field + lcdd->field().magneticField( + position, magneticFieldVector); // get the magnetic field vector from DD4hep + float Bz = magneticFieldVector[2] / dd4hep::tesla; + + float pt = fabs(0.3 * Bz / track->getOmega() / 1000); h_pt->Fill(pt); float lambda = std::atan(track->getTanLambda()); diff --git a/src/TrackPerfHistProc.cxx b/src/TrackPerfHistProc.cxx index cbcd6c5..95cfc4c 100644 --- a/src/TrackPerfHistProc.cxx +++ b/src/TrackPerfHistProc.cxx @@ -6,6 +6,10 @@ #include #include +// DD4hep +#include +#include + #include #include #include @@ -76,6 +80,8 @@ void TrackPerfHistProc::init() { tree->cd("../efficiency"); _effiPlots = std::make_shared(true); _fakePlots = std::make_shared(false); + + _lcdd = &dd4hep::Detector::getInstance(); } void TrackPerfHistProc::processRunHeader(LCRunHeader* /*run*/) {} @@ -135,7 +141,7 @@ void TrackPerfHistProc::processEvent(LCEvent* evt) { static_cast(trkCol->getElementAt(i)); trkSet.insert(trk); - _allTracks->fill(trk); + _allTracks->fill(trk, _lcdd); } h_number_of_tracks->Fill(trkSet.size()); @@ -160,11 +166,11 @@ void TrackPerfHistProc::processEvent(LCEvent* evt) { if (rel->getWeight() > _matchProb) { if (trkSet.find(trk) != trkSet.end()) { - _realTracks->fill(trk); + _realTracks->fill(trk, _lcdd); _realTruths->fill(mcp); _effiPlots->fillMC(mcp, true); - _realReso->fill(trk, mcp); - _fakePlots->fillTrack(trk, false); + _realReso->fill(trk, mcp, _lcdd); + _fakePlots->fillTrack(trk, false, _lcdd); mcpSet.erase(mcp); trkSet.erase(trk); @@ -179,8 +185,8 @@ void TrackPerfHistProc::processEvent(LCEvent* evt) { _effiPlots->fillMC(mcp, false); } for (const EVENT::Track* trk : trkSet) { - _fakeTracks->fill(trk); - _fakePlots->fillTrack(trk, true); + _fakeTracks->fill(trk, _lcdd); + _fakePlots->fillTrack(trk, true, _lcdd); } h_number_of_fakes->Fill(trkSet.size()); }