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

GenFind use scripts added/modified

parent a91020b0
// 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"
#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 = 20,
int Ntracks = 1)
{
// read the simulated data and set up GenFind
gSystem->Load("libGenFind.so");
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; }
int nentries=t->GetEntries();
if(i_event+Ntracks>nentries){
std::cout << "Only " << nentries << " in file" << std::endl;
return;
}
std::vector<DD4hep::Simulation::Geant4Tracker::Hit*> * g4hits = nullptr;
std::vector<DD4hep::Simulation::Geant4Tracker::Hit*> * g4hits2 = nullptr;
t->SetBranchAddress("SiVertexBarrelHits", &g4hits);
t->SetBranchAddress("SiTrackerBarrelHits", &g4hits2);
HoughTransform* ht = new HoughTransform();
// set up the visualization for the final tracks
DD4hep::Geometry::LCDD& lcdd = DD4hep::Geometry::LCDD::getInstance();
lcdd.fromCompact("JLEIC.xml");
const double position[3] = {0,0,0};
double bField[3] = {0,0,0};
lcdd.field().magneticField(position,bField);
double Bz = bField[2]/dd4hep::tesla;
auto vol_man = lcdd.volumeManager();
genfit::FieldManager::getInstance()->init(new genfit::ConstField(0.,0.,Bz*10)); //genfit uses kgauss
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);
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);
// 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);
// combine hits to create composit event
std::vector<XYZTVector> hits ;
std::vector<int> cellID_hits ;
std::vector<int> tracktruth ;
for(int i_track = 0; i_track < Ntracks; i_track++){ // made tracks by combining events
t->GetEntry(i_event);
i_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;
}
}
std::cout << "Created 'event' by combining " << Ntracks << " real events from file\n";
// 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];
//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);
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);
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 };
TVector3 zdir = {0,0,1};
TVector3 planeNorm = { hits[std::get<0>(thit)].X()*.1, hits[std::get<0>(thit)].Y()*.1, 0.0 };
planeNorm = planeNorm.Unit();
genfit::SharedPlanePtr plane(new genfit::DetPlane(point, zdir.Cross(planeNorm), zdir) );
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)]);
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;
fitTrack.checkConsistency();
// do the fit
std::cout << "about to do the fit" << std::endl;
fitTrack.getFittedState().Print();
// check
std::cout << "again check consistency" << std::endl;
fitTrack.checkConsistency();
std::cout << "adding event track to display" << std::endl;
display->addEvent(&fitTrack);
}
display->open();
}
......@@ -44,7 +44,7 @@ void HT_test(
for(int i_track = 0; i_track < Ntracks; i_track++){ // made tracks by combining events
t->GetEntry(i_event);
i_event++;
i_event++; // this causes an error, because it increments before you do anything to the track
for( auto thit : (*g4hits) ) {
hits.push_back( {thit->position.X(), thit->position.Y(), thit->position.Z(), 0.0} );
......@@ -124,7 +124,7 @@ void HT_test(
int ref_truth_track = tracktruth[std::get<0>(atrack[0])];
int ref_diff = tracknumber - ref_truth_track;
//std::cout << "ref diff is " << ref_diff << std::endl;
double grouped_test = 0;
int hit_diff=0;
int numhit=0;
for( std::tuple<int,double,double> thit : atrack ) {
......@@ -136,7 +136,11 @@ void HT_test(
//singletrack->Fill( std::get<1>(thit), std::get<2>(thit) );
singletrack->Fill( hits[std::get<0>(thit)].X(), hits[std::get<0>(thit)].Y() );
}
double grouped_test = hit_diff/numhit;
if (numhit>0) {
grouped_test = hit_diff/numhit;
} else {
grouped_test = -1;
}
std::cout << "grouped test " << grouped_test << std::endl;
FoundTrackHists.push_back(singletrack);
}
......
#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
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(i_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 : (*g4hits) ) {
TVector3 point = {thit->position.X()*0.1, thit->position.Y()*0.1, thit->position.Z()*0.1};
TVector3 zdir = {0,0,1};
TVector3 planeNorm = {thit->position.X()*0.1, thit->position.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;
hitCoords[1] = 0.0;
genfit::PlanarMeasurement* measurement = new genfit::PlanarMeasurement(hitCoords, hitCov, detId, ++hitId, nullptr);
static_cast<genfit::PlanarMeasurement*>(measurement)->setPlane(plane, thit->cellID);
fitTrack.insertPoint(new genfit::TrackPoint(measurement, &fitTrack));
}*/
for( auto thit : (*g4hits2) ) {
std::cout << "hits coords: " << thit->position.X() << " , " << thit->position.Y() << std::endl;
std::cout << "cellID value " << thit->cellID << std::endl;
}
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;
}
for( auto thit : (*g4hits2) ) {
std::cout << "hits coords: " << thit->position.X() << " , " << thit->position.Y() << std::endl;
std::cout << "cellID value " << thit->cellID << std::endl;
}
// combine hits to create composit event
std::vector<XYZTVector> hits ;
std::vector<int> cellID_hits ;
std::vector<int> tracktruth ;
for(int i_track = 0; i_track < Ntracks; i_track++){ // made tracks by combining events
t->GetEntry(i_event);
i_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 << "hits coords: " << thit->position.X() << " , " << thit->position.Y() << std::endl;
std::cout << "cellID value " << thit->cellID << std::endl;
}
}
std::cout << "take two! " << 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 << "hits coords: " << thit->position.X() << " , " << thit->position.Y() << std::endl;
std::cout << "cellID value " << thit->cellID << std::endl;
}
/* for( int ihit=0; ihit< hits.size(); ihit++ ) {
TVector3 point = {401.223*0.1, -8.89552*0.1, -167.461*0.1};
//TVector3 point = { hits[1].X()*.1, hits[1].Y()*.1, hits[1].Z()*.1 };
TVector3 zdir = {0,0,1};
TVector3 planeNorm = {401.223*0.1, -8.89552*0.1,0.0};
//TVector3 planeNorm = { hits[1].X()*.1, hits[1].Y()*.1, 0.0 };
planeNorm = planeNorm.Unit();
genfit::SharedPlanePtr plane(new genfit::DetPlane(point, zdir.Cross(planeNorm), zdir) );
hitCoords[0] = 0.0;
hitCoords[1] = 0.0;
genfit::PlanarMeasurement* measurement = new genfit::PlanarMeasurement(hitCoords, hitCov, detId, ++hitId, nullptr);
measurement->setPlane(plane, 545917228);
// measurement->setPlane(plane, cellID_hits[1]);
fitTrack.insertPoint(new genfit::TrackPoint(measurement, &fitTrack));
std::cout << " x position is " << hits[0].X() << std::endl;
std::cout << " y position is " << hits[0].Y() << std::endl;
std::cout << " z position is " << hits[0].Z() << std::endl;
std::cout << "cellID is " << cellID_hits[0] << std::endl;
// }*/
//check
fitTrack.checkConsistency();
// do the fit
fitter->processTrack(&fitTrack);
// print fit result
fitTrack.getFittedState().Print();
//check
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();
}
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