Commit bd0c4c2a authored by Johnston's avatar Johnston
Browse files

update documentation, remove most confusing, broken, scripts

parent 20e04efe
GenFind
=======
GenFit: existing package with good general fitting framework for a number of
detectors.
Generic "Track Finder" package to provide appropriate "seed" tracks for GenFit
without any hard-coded algorithmic dependencies to specific detector
geometries.
Requires already found seed tracks.
Basic, general finding algorithm needed, this basic functionality to be
provided by GenFind, this package.
Motivation for development:
---------------------------
Track hit pruning, optimization algorithms in GenFit make even imperfect
Finding algorithm sufficient.
GenFit exists as a "generic" package with general fitting methods available to
reconstruct tracks in a variety of detectors. However, GenFit requires "seed"
tracks to start with. Once given seeds, GenFit can use hit pruning and
optimization algorithms to produce good reconstruction with imperfect
beginnings.
Proposal:
GenFind is a generic package that determines suitable "seed" tracks from whole
events. The input events can contain an unspecified number of tracks, and
GenFind will separate out the hits corresponding to each component track as an
isolated "seed," and pass these along to GenFit.
Simple transform (Hough or otherwise) to feature space given event hits,
binning, and maximum bin selection process to identify candidate hits.
GenFind contains methods which create a conformal map of input hits, which will
turn helical tracks into straight lines in conformal space. GenFind then uses a
modified Hough Transform method to find lines via locating intersections in
rho-theta space.
Need:
In addition to the library of transform and intersection finding methods,
GenFind contains sample scripts to illustrate how to use these methods, and how
to interface the track output in GenFind with the fitters from GenFit.
Transforms suitable to geometry, B field. (Hough transform first)
Data structure to hold candidate hits. (vector of slcio hit type?)
API of functions to talk to GenFit, dd4hep (TGeo geometry, as necessary)
Wrap this, and Genfit, in Marlin Processors to actually use.
Further development:
--------------------
Currently, thresholds necessary for locating the intersections corresponding to
track signatures are set by hand. These need to be dynamically set by reading
geometries in from DD4hep.
Currently, relatively simple hit and track data structures are used, add
something more sophisticated.
Project Template
================
Further develop standardized communication methods between GenFind, GenFit,
DD4hep, and TGeo.
GenFind built on top of template cmake build including Root-compatible library
creation and access to root data types. Template documentation follows:
Test performance numerically for efficiency and reconstruction, in more than
one concept detector, to maximize true generality.
A project template for building shared libraries and binaries with cmake. It
also installs headers and the proper cmake files so that you can use
`find_package(this_project)` in subsequent cmake builds.
Getting started
---------------
Example scripts:
----------------
These have been developed and tested on simulation output generated by DD4hep
using the JLEIC tracker barrel and vertex detector geometries.
Here we install to your home directory (instead of /usr/local):
```bash
git clone https://github.com/whit2333/project_template.git
mkdir project_build
cd project_build
cmake ../project_template/. -DCMAKE_INSTALL_PREFIX=$HOME
make && make install
```
GenFind/src/GenFind/tests/HT_test.cxx
```
export PATH=$HOME/bin:$PATH # if it is not included already
export LD_LIBRARY_PATH=$HOME/lib:$LD_LIBRARY_PATH # if it is not included already
```
Most up to date script illustrating basic tests:
1) Combines single track simulated events to "fake multi-track events."
2) Conformally maps all points in the "multi-track event."
3) Performs modified hough transform and line-finding (one call in script.)
4) Illustrates plotting of up to three separated tracks, as well as the hough
transform functions for each conformal hit.
5) Includes limited performance testing for accuracy of track separation.
Make sure the last two lines are in your bashrc if not a standard location.
Comments
--------
GenFind/src/GenFind/tests/DisaggregatedHits.cxx
This template is not meant to be the best example and it is probably not using
the most current features of cmake. So please don't use it as a sole source for
learning cmake because cmake has wonderful documentation.
Script combines the output of the GenFind track separation algorithms, with the
GenFit reconstruction of track kinematics.
It is meant to provide a quick starting point.
1) Initializes geometry from DD4hep and initializes GenFit parameters.
2) Performs steps 1-3 from HT_test.cxx script.
3) Feeds the "seed tracks" separated by GenFind into GenFit track fitter.
4) prints the kinematic results
GenFind/src/GenFind/tests/IncludeVertexHits.cxx
Please feel free to send questions or comments to Whitney Armstrong
(whit@jlab.org). Also pull requests are greatly appreciated!
Does the same as DisaggregatedHits, but for both vertex and tracker barrel
hits.
Todo
----
Caveats:
For weird events, strong curlers, lots of extra hits, etc. these algorithms can
fail and do so with segmentation faults, limiting ability to automatically
check performance numerically.
1. Add renaming script
2. Clean up cmake files
3. Document features/convention
#include "ConstField.h"
#include "Exception.h"
#include "FieldManager.h"
#include "KalmanFitterRefTrack.h"
#include "DAF.h"
#include "StateOnPlane.h"
#include "Track.h"
#include "TrackPoint.h"
#include <MaterialEffects.h>
#include <RKTrackRep.h>
#include <TGeoMaterialInterface.h>
#include <EventDisplay.h>
#include <HelixTrackModel.h>
#include <MeasurementCreator.h>
#include <WireMeasurement.h>
#include "PlanarMeasurement.h"
#include <TDatabasePDG.h>
#include <TEveManager.h>
#include <TGeoManager.h>
#include <TRandom.h>
#include <TVector3.h>
#include <vector>
#include "TDatabasePDG.h"
#include <TMath.h>
void GenFindGenFit(
int i_event = 0,
int Ntracks = 1)
{
using namespace genfind;
using namespace ROOT::Math;
int increment_event = i_event;
TFile * f = new TFile("simple_example_out.root","READ");
TTree * t = (TTree*)gROOT->FindObject("EVENT");
if(!t) { std::cout << " Tree not found " << std::endl; return;}
std::vector<DD4hep::Simulation::Geant4Particle*> * mc_particles = nullptr;
std::vector<DD4hep::Simulation::Geant4Tracker::Hit*> * g4hits = nullptr;
std::vector<DD4hep::Simulation::Geant4Tracker::Hit*> * g4hits2 = nullptr;
t->SetBranchAddress("SiVertexBarrelHits", &g4hits);
t->SetBranchAddress("SiTrackerBarrelHits", &g4hits2);
t->SetBranchAddress("MCParticles", &mc_particles);
// Load the geometry
DD4hep::Geometry::LCDD& lcdd = DD4hep::Geometry::LCDD::getInstance();
lcdd.fromCompact("JLEIC.xml");
const double position[3] = {0,0,0}; // position to calculate magnetic field at (the origin in this case)
double bField[3] = {0,0,0};
lcdd.field().magneticField(position,bField);
double Bz = bField[2]/dd4hep::tesla;
std::cout << " Magnetic Field Bz = " << Bz << std::endl;
auto vol_man = lcdd.volumeManager();
genfit::FieldManager::getInstance()->init(new genfit::ConstField(0.,0., Bz*10.0)); // gentfit uses kilo-Gauss
genfit::MaterialEffects::getInstance()->init(new genfit::TGeoMaterialInterface());
genfit::EventDisplay* display = genfit::EventDisplay::getInstance();
display->setOptions("X");
// Add the geometry to eve display
TGeoNode* node1 = gGeoManager->GetTopNode();
TEveGeoTopNode* its = new TEveGeoTopNode(gGeoManager, node1);
gEve->AddGlobalElement(its);
//------------------------------------------------
//
double degree = TMath::Pi()/180.0;
const int pdg = 11; // particle pdg code
genfit::AbsKalmanFitter* fitter = new genfit::DAF();
// trackrep
genfit::AbsTrackRep* rep = new genfit::RKTrackRep(pdg);
// smeared start state
genfit::MeasuredStateOnPlane stateSmeared(rep);
// approximate covariance
TMatrixDSym covM(6);
double res1 = 0.01;
for (int i = 0; i < 3; ++i)
covM(i,i) = res1*res1;
for (int i = 3; i < 6; ++i)
covM(i,i) = pow(res1/9.0/sqrt(3.0), 2);
const int detId(0); // detector ID
int planeId(0); // detector plane ID
int hitId(0); // hit ID
double detectorResolution(0.0001); // resolution of planar detectors
TMatrixDSym hitCov(2);
hitCov.UnitMatrix();
hitCov *= detectorResolution*detectorResolution;
TVectorD hitCoords(2);
for(int i_track = 0; i_track < Ntracks; i_track++) {
/*
t->GetEntry(increment_event);
//i_event++;
if( (*g4hits).size() == 0 ) {
continue;
}
auto first_hit_pos = (*((*g4hits)[0])).position;
std::cout << "first hit pos x " << first_hit_pos.X() << std::endl;
std::cout << "first hit pos y " << first_hit_pos.Y() << std::endl;
std::cout << "first hit pos z " << first_hit_pos.Z() << std::endl;
// start values for the fit, e.g. from pattern recognition
TVector3 pos(0.0,0.0,0.0);
pos = 0.1*pos;
TVector3 mom(first_hit_pos.X(), first_hit_pos.Y(), first_hit_pos.Z());
mom.SetMag(2.0);
stateSmeared.setPosMomCov(pos, mom, covM);
// create track
TVectorD seedState(6);
TMatrixDSym seedCov(6);
stateSmeared.get6DStateCov(seedState, seedCov);
genfit::Track fitTrack(rep, seedState, seedCov);
*/
/*
for( auto thit : (*g4hits2) ) {
TVector3 point = {thit->position.X()*0.1, thit->position.Y()*0.1, thit->position.Z()*0.1};
//TVector3 point = {401.223*0.1, -8.89552*0.1, -167.461*0.1};
TVector3 zdir = {0,0,1};
TVector3 planeNorm = {thit->position.X()*0.1, thit->position.Y()*0.1,0.0};
//TVector3 planeNorm = {401.223*0.1, -8.89552*0.1,0.0};
planeNorm = planeNorm.Unit();
genfit::SharedPlanePtr plane(new genfit::DetPlane(point, zdir.Cross(planeNorm) , zdir) );
// add some planar hits to track with coordinates
hitCoords[0] = 0.0;//thit->position.X()*0.1;
hitCoords[1] = 0.0;//thit->position.Y()*0.1;
genfit::PlanarMeasurement* measurement = new genfit::PlanarMeasurement(hitCoords, hitCov, detId, ++hitId, nullptr);
//measurement->setPlane(plane, 545917228);
measurement->setPlane(plane, thit->cellID);
fitTrack.insertPoint(new genfit::TrackPoint(measurement, &fitTrack));
std::cout << " x position is " << thit->position.X() << std::endl;
std::cout << " y position is " << thit->position.Y() << std::endl;
std::cout << " z position is " << thit->position.Z() << std::endl;
std::cout << "celID is " << thit->cellID << std::endl;
}
//check
fitTrack.checkConsistency();
// do the fit
fitter->processTrack(&fitTrack);
// print fit result
fitTrack.getFittedState().Print();
//check
fitTrack.checkConsistency();
display->addEvent(&fitTrack);
*/
increment_event++;
}
//std::cout << " |P| = " << (*mc_particles)[0]->psx << ", " << (*mc_particles)[0]->psy << ", " << (*mc_particles)[0]->psz << std::endl ;
// open event display
//display->open();
// finished doing it the cheating way. Now, check what GenFind would have done.
// combine hits to create composit event
std::vector<XYZTVector> hits ;
std::vector<int> cellID_hits ;
std::vector<int> tracktruth ;
increment_event = i_event;
for(int i_track = 0; i_track < Ntracks; i_track++){ // made tracks by combining events
t->GetEntry(increment_event);
//for( auto thit : (*g4hits) ) {
// hits.push_back( {thit->position.X(), thit->position.Y(), thit->position.Z(), 0.0} );
// cellID_hits.push_back( thit->cellID );
// tracktruth.push_back(i_track);
// std::cout << thit->position.X() << " , " << thit->position.Y() << std::endl;
//}
for( auto thit : (*g4hits2) ) {
hits.push_back( {thit->position.X(), thit->position.Y(), thit->position.Z(), 0.0} );
cellID_hits.push_back( thit->cellID );
tracktruth.push_back(i_track);
std::cout << thit->position.X() << " , " << thit->position.Y() << std::endl;
}
increment_event++;
}
std::cout << "Created 'event' by combining " << Ntracks << " real events from file\n";
HoughTransform* ht = new HoughTransform();
// run the standard GenFind method over composit event to disaggregate via algorithm
std::vector<ROOT::Math::XYZTVector> conf_hits = ht->UseConformalMap(hits);
auto mg = ht->FillHoughTransforms(conf_hits);
std::vector<std::vector<std::tuple<int,double,double>>> FoundTracks = ht->FindTracks(conf_hits);
std::cout << "FindTracks found " << FoundTracks.size() << " tracks " << endl;
for( std::vector<std::tuple<int,double,double>> atrack : FoundTracks ) {
t->GetEntry(i_event); // might have problems with multiple tracks this way?
if( (*g4hits).size() == 0 ) {
continue;
}
auto first_hit_pos = (*((*g4hits)[0])).position;
std::cout << "first hit pos x " << first_hit_pos.X() << std::endl;
std::cout << "first hit pos y " << first_hit_pos.Y() << std::endl;
std::cout << "first hit pos z " << first_hit_pos.Z() << std::endl;
// start values for the fit, e.g. from pattern recognition
TVector3 pos(0.0,0.0,0.0);
pos = 0.1*pos;
TVector3 mom(first_hit_pos.X(), first_hit_pos.Y(), first_hit_pos.Z());
mom.SetMag(2.0);
stateSmeared.setPosMomCov(pos, mom, covM);
// create track
TVectorD seedState(6);
TMatrixDSym seedCov(6);
stateSmeared.get6DStateCov(seedState, seedCov);
genfit::Track fitTrack(rep, seedState, seedCov);
for( auto thit : atrack ) {
std::cout << " did a hit! " << std::endl;
std::cout << "x position is " << hits[std::get<0>(thit)].X() << std::endl;
std::cout << "y position is " << hits[std::get<0>(thit)].Y() << std::endl;
std::cout << "z position is " << hits[std::get<0>(thit)].Z() << std::endl;
TVector3 point = {hits[std::get<0>(thit)].X()*0.1, hits[std::get<0>(thit)].Y()*0.1, hits[std::get<0>(thit)].Z()*0.1};
TVector3 zdir = {0,0,1};
TVector3 planeNorm = {hits[std::get<0>(thit)].X()*0.1, hits[std::get<0>(thit)].Y()*0.1,0.0};
planeNorm = planeNorm.Unit();
genfit::SharedPlanePtr plane(new genfit::DetPlane(point, zdir.Cross(planeNorm) , zdir) );
// add some planar hits to track with coordinates
hitCoords[0] = 0.0;//thit->position.X()*0.1;
hitCoords[1] = 0.0;//thit->position.Y()*0.1;
genfit::PlanarMeasurement* measurement = new genfit::PlanarMeasurement(hitCoords, hitCov, detId, ++hitId, nullptr);
measurement->setPlane(plane, cellID_hits[std::get<0>(thit)]);
std::cout << "cellID is " << cellID_hits[std::get<0>(thit)] << std::endl;
fitTrack.insertPoint(new genfit::TrackPoint(measurement, &fitTrack));
}
//check
fitTrack.checkConsistency();
// do the fit
fitter->processTrack(&fitTrack);
// print fit result
std::cout << " printing fit results: " << std::endl;
fitTrack.getFittedState().Print();
//check
std::cout << " checking fit results: " << std::endl;
fitTrack.checkConsistency();
display->addEvent(&fitTrack);
}
std::cout << " |P| = " << (*mc_particles)[0]->psx << ", " << (*mc_particles)[0]->psy << ", " << (*mc_particles)[0]->psz << std::endl ;
// open event display
display->open();
}
#include "GenFind/GenFindHits.h"
void conformal_transform_test()
{
using namespace genfind;
std::vector<std::shared_ptr<Hit>> hits;
std::vector<ConformalHit> chits;
std::vector<ROOT::Math::XYZTVector> hits2;
std::vector<ConformalHit> chits2;
double degree = TMath::Pi()/180.0;
double r0 = 10.0;
double x0 = 0.0;
double y0 = r0;
ROOT::Math::XYZTVector ref_hit = {0.0,0.0,0.0,0.0};
for(int ith=0; ith < 10; ith++) {
double x = r0*TMath::Sin((5.0 + ith*3.0)*degree)- x0 ;
double y = TMath::Sqrt(r0*r0 - (x+x0)*(x+x0)) - y0 ;
auto ahit = std::make_shared<Hit>(x,y,0);
hits.push_back(ahit);
chits.push_back(compute_conformal_hit(ahit));
std::cout << " x " << x << ", y " << y << std::endl;
ROOT::Math::XYZTVector pos_hit = {x,y,0.0,0.0};
hits2.push_back(pos_hit);
//chits2.push_back(compute_conformal_hit(pos_hit, ref_hit));
}
chits2 = compute_conformal_hits(hits2,ref_hit);
TGraph* gr_hits = new TGraph();
int i_point = 0;
for(auto ahit : hits){
gr_hits->SetPoint(i_point, ahit->X(), ahit->Y());
i_point++;
}
TGraph* gr_chits = new TGraph();
i_point = 0;
for(auto ahit : chits){
gr_chits->SetPoint(i_point, ahit.X(), ahit.Y());
std::cout << " x " << ahit.X() << ", y " << ahit.Y() << std::endl;
i_point++;
}
TGraph* gr_hits2 = new TGraph();
i_point = 0;
for(auto ahit : hits2){
gr_hits2->SetPoint(i_point, ahit.X(), ahit.Y());
i_point++;
}
TGraph* gr_chits2 = new TGraph();
i_point = 0;
for(auto ahit : chits2){
gr_chits2->SetPoint(i_point, ahit.X(), ahit.Y());
std::cout << " x " << ahit.X() << ", y " << ahit.Y() << std::endl;
i_point++;
}
TCanvas* c = new TCanvas();
c->Divide(2,2);
c->cd(1);
gr_hits->SetMarkerStyle(20);
gr_hits->Draw("ap");
c->cd(2);
gr_chits->SetMarkerStyle(20);
gr_chits->Draw("ap");
c->cd(3);
gr_hits2->SetMarkerStyle(20);
gr_hits2->Draw("ap");
c->cd(4);
gr_chits2->SetMarkerStyle(20);
gr_chits2->Draw("ap");
}
#include "ConstField.h"
#include "Exception.h"
#include "FieldManager.h"
#include "KalmanFitterRefTrack.h"
#include "DAF.h"
#include "StateOnPlane.h"
#include "Track.h"
#include "TrackPoint.h"
#include <MaterialEffects.h>
#include <RKTrackRep.h>
#include <TGeoMaterialInterface.h>
#include <EventDisplay.h>
#include <HelixTrackModel.h>
#include <MeasurementCreator.h>
#include <WireMeasurement.h>
#include "PlanarMeasurement.h"
#include <TDatabasePDG.h>
#include <TEveManager.h>
#include <TGeoManager.h>
#include <TRandom.h>
#include <TVector3.h>
#include <vector>
#include "TDatabasePDG.h"
#include <TMath.h>
void fit_test1_nextevent(
int i_event = 20,
int Ntracks = 1)
{
using namespace genfind;
using namespace ROOT::Math;
TFile * f = new TFile("simple_example_out.root","READ");
TTree * t = (TTree*)gROOT->FindObject("EVENT");
if(!t) { std::cout << " Tree not found " << std::endl; return;}
std::vector<DD4hep::Simulation::Geant4Particle*> * mc_particles = nullptr;
std::vector<DD4hep::Simulation::Geant4Tracker::Hit*> * g4hits = nullptr;
std::vector<DD4hep::Simulation::Geant4Tracker::Hit*> * g4hits2 = nullptr;
t->SetBranchAddress("SiVertexBarrelHits", &g4hits);
t->SetBranchAddress("SiTrackerBarrelHits", &g4hits2);
t->SetBranchAddress("MCParticles", &mc_particles);
// Load the geometry
DD4hep::Geometry::LCDD& lcdd = DD4hep::Geometry::LCDD::getInstance();
lcdd.fromCompact("JLEIC.xml");
const double position[3] = {0,0,0}; // position to calculate magnetic field at (the origin in this case)
double bField[3] = {0,0,0};
lcdd.field().magneticField(position,bField);
double Bz = bField[2]/dd4hep::tesla;
std::cout << " Magnetic Field Bz = " << Bz << std::endl;
auto vol_man = lcdd.volumeManager();
genfit::FieldManager::getInstance()->init(new genfit::ConstField(0.,0., Bz*10.0)); // gentfit uses kilo-Gauss
genfit::MaterialEffects::getInstance()->init(new genfit::TGeoMaterialInterface());
genfit::EventDisplay* display = genfit::EventDisplay::getInstance();
display->setOptions("X");
// Add the geometry to eve display
TGeoNode* node1 = gGeoManager->GetTopNode();
TEveGeoTopNode* its = new TEveGeoTopNode(gGeoManager, node1);
gEve->AddGlobalElement(its);
//------------------------------------------------
//
double degree = TMath::Pi()/180.0;
const int pdg = 11; // particle pdg code
genfit::AbsKalmanFitter* fitter = new genfit::DAF();
// trackrep
genfit::AbsTrackRep* rep = new genfit::RKTrackRep(pdg);
// smeared start state
genfit::MeasuredStateOnPlane stateSmeared(rep);
// approximate covariance