Commit 02f6cfc7 authored by Johnston's avatar Johnston
Browse files

GenFindGenFit has genfit work for a single track found by genfind

parent 02af0235
// this script tests combined function of GenFind and GenFit
//
// First, simulated data events are artificially combined
// to create multitrack events.
//
// Then, GenFind functions are called to separate the hits
// in the event into individual tracks.
//
// Finally, the individual tracks found by GenFind are fed
// into the GenFit framework, and displayed in a visualized
// version of the simulated detector.
#include "ConstField.h"
#include "Exception.h"
#include "FieldManager.h"
......@@ -18,7 +6,7 @@
#include "StateOnPlane.h"
#include "Track.h"
#include "TrackPoint.h"
#include <MaterialEffects.h>
#include <RKTrackRep.h>
#include <TGeoMaterialInterface.h>
......@@ -36,65 +24,63 @@
#include <TRandom.h>
#include <TVector3.h>
#include <vector>
#include "TDatabasePDG.h"
#include <TMath.h>
void GenFindGenFit(
int i_event = 20,
void simplefit_test(
int i_event = 0,
int Ntracks = 1)
{
// read the simulated data and set up GenFind
gSystem->Load("libGenFind.so");
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; }
int nentries=t->GetEntries();
if(!t) { std::cout << " Tree not found " << std::endl; return;}
if(i_event+Ntracks>nentries){
std::cout << "Only " << nentries << " in file" << 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);
HoughTransform* ht = new HoughTransform();
// set up the visualization for the final tracks
// Load the geometry
DD4hep::Geometry::LCDD& lcdd = DD4hep::Geometry::LCDD::getInstance();
lcdd.fromCompact("JLEIC.xml");
const double position[3] = {0,0,0};
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)); //genfit uses kgauss
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 geometry to eve display
TGeoNode* node1 = gGeoManager->GetTopNode();
TEveGeoTopNode* its = new TEveGeoTopNode(gGeoManager, node1);
// Add the geometry to eve display
TGeoNode* node1 = gGeoManager->GetTopNode();
TEveGeoTopNode* its = new TEveGeoTopNode(gGeoManager, node1);
gEve->AddGlobalElement(its);
// initialize genfit
//------------------------------------------------
//
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);
......@@ -105,7 +91,8 @@ void GenFindGenFit(
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
......@@ -116,14 +103,96 @@ void GenFindGenFit(
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(i_event);
i_event++;
t->GetEntry(increment_event);
//for( auto thit : (*g4hits) ) {
// hits.push_back( {thit->position.X(), thit->position.Y(), thit->position.Z(), 0.0} );
......@@ -138,78 +207,95 @@ void GenFindGenFit(
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 ) {
// here I can call to genfit with the atrack? or something else?
// maybe make a new object with only the position information of the
// hits in atrack, and then call genfit with that?
auto fhit = atrack[0];
t->GetEntry(i_event); // might have problems with multiple tracks this way?
//start values for the fit, e.g. from pattern recognition
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(39.5163, -3.46122, -16.6242);
TVector3 mom(first_hit_pos.X(), first_hit_pos.Y(), first_hit_pos.Z());
mom.SetMag(2.0);
std::cout << "first hit pos x " << hits[std::get<0>(fhit)].X() << std::endl;
std::cout << "first hit pos y " << hits[std::get<0>(fhit)].Y() << std::endl;
std::cout << "first hit pos z " << hits[std::get<0>(fhit)].Z() << std::endl;
stateSmeared.setPosMomCov(pos, mom, covM);
// create track
TVectorD seedState(6);
TVectorD seedState(6);
TMatrixDSym seedCov(6);
stateSmeared.get6DStateCov(seedState, seedCov);
genfit::Track fitTrack(rep, seedState, seedCov);
std::cout << "about to go through hits in the track" << std::endl;
for( auto thit : atrack ) {
TVector3 point = { hits[std::get<0>(thit)].X()*.1, hits[std::get<0>(thit)].Y()*.1, hits[std::get<0>(thit)].Z()*.1 };
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()*.1, hits[std::get<0>(thit)].Y()*.1, 0.0 };
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) );
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;
hitCoords[0] = 0.0;
hitCoords[1] = 0.0;
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));
std::cout << "did one of the hits" << 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;
std::cout << "cellID is " << cellID_hits[std::get<0>(thit)] << std::endl;
}
// check
std::cout << "about to check consistency" << std::endl;
}
//check
fitTrack.checkConsistency();
// do the fit
std::cout << "about to do the fit" << std::endl;
fitter->processTrack(&fitTrack);
// print fit result
fitTrack.getFittedState().Print();
// check
std::cout << "again check consistency" << std::endl;
//check
fitTrack.checkConsistency();
std::cout << "adding event track to display" << std::endl;
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();
}
......@@ -104,7 +104,7 @@ void fit_test1_nextevent(
for(int i_track = 0; i_track < Ntracks; i_track++) {
t->GetEntry(i_event);
i_event++;
//i_event++;
if( (*g4hits).size() == 0 ) {
continue;
......@@ -256,6 +256,7 @@ void fit_test1_nextevent(
display->addEvent(&fitTrack);
i_event++;
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment