Commit e936f42a authored by Whitney Armstrong's avatar Whitney Armstrong
Browse files

modified: fit_test1.cxx

	new file:   fit_test2.cxx
parent d37766db
......@@ -50,11 +50,11 @@ void fit_test1(
// 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;
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& id_decoder = DD4hep::DDRec::IDDecoder::getInstance();
auto vol_man = lcdd.volumeManager();
......@@ -65,40 +65,29 @@ void fit_test1(
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::KalmanFitterRefTrack();
genfit::AbsKalmanFitter* fitter = new genfit::DAF();
//HoughTransform* ht = new HoughTransform();
//TH2F * hxy = new TH2F("hxy", "hxy", 100, -1000, 1000, 100, -1000, 1000);
//TH2F * hxy2 = new TH2F("hxy2", "hxy2", 50, -100 , 100, 50, -100, 100);
//TH2F * hpeaky = new TH2F("hpeakxy", "hxy", 100, -1000, 1000, 100, -1000, 1000);
//TH2F * hpeaky2 = new TH2F("hpeakxy2", "hxy2", 50, -100 , 100, 50, -100, 100);
//TH2F * huv = new TH2F("huv", "huv", 100, -0.12, 0.12 , 100, -0.12, 0.12 );
//TH2F * hrphi = new TH2F("hrphi", "hrphi", 90, 0, 180, 100, -0.2, 0.2);
//TH1F * hphi = new TH1F("hphi", "hphi", 50, -4, 4 );
//TH1F * htheta = new TH1F("htheta", "htheta", 180, 0, 180 );
//std::vector<XYZTVector> hits ;
//
// helix track model
const int pdg = 11;//1000020040;//2212; // particle pdg code
//const double charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge()/(3.);
//genfit::HelixTrackModel* helix = new genfit::HelixTrackModel(pos, mom, charge);
double degree = TMath::Pi()/180.0;
// trackrep
genfit::AbsTrackRep* rep = new genfit::RKTrackRep(pdg);
// smeared start state
genfit::MeasuredStateOnPlane stateSmeared(rep);
// approximate covariance
TMatrixDSym covM(6);
double res1 = 0.01;
......
#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_test2(
int i_event = 26,
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& id_decoder = DD4hep::DDRec::IDDecoder::getInstance();
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::KalmanFitterRefTrack();
genfit::AbsKalmanFitter* fitter = new genfit::DAF();
// helix track model
//const double charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge()/(3.);
//genfit::HelixTrackModel* helix = new genfit::HelixTrackModel(pos, mom, charge);
//// trackrep
//genfit::AbsTrackRep* rep = new genfit::RKTrackRep(pdg);
//// smeared start state
//genfit::MeasuredStateOnPlane stateSmeared(rep);
// approximate covariance
TMatrixDSym covM(6);
double res1 = 0.001;
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);
//-----------------------------------------------
std::vector<genfit::Track*> all_tracks;
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;
// trackrep
genfit::AbsTrackRep* rep = new genfit::RKTrackRep(pdg);
// smeared start state
genfit::MeasuredStateOnPlane stateSmeared(rep);
// 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(1.0);
stateSmeared.setPosMomCov(pos, mom, covM);
// Create track.
TVectorD seedState(6);
TMatrixDSym seedCov(6);
stateSmeared.get6DStateCov(seedState, seedCov);
auto aTrack = new genfit::Track(rep, seedState, seedCov);
auto& fitTrack = *aTrack;
all_tracks.push_back(aTrack);
//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};//point.Unit();//clas12::geo::Convert(svt_geometry.GetChannelNorm(chan)/CLHEP::cm);
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);
static_cast<genfit::PlanarMeasurement*>(measurement)->setPlane(plane, thit->cellID);
fitTrack.insertPoint(new genfit::TrackPoint(measurement, &fitTrack));
}
for( auto thit : (*g4hits2) ) {
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};//point.Unit();//clas12::geo::Convert(svt_geometry.GetChannelNorm(chan)/CLHEP::cm);
planeNorm = planeNorm.Unit();
genfit::SharedPlanePtr plane(new genfit::DetPlane(point, zdir.Cross(planeNorm) , zdir) );
//genfit::SharedPlanePtr plane(new genfit::DetPlane(point, planeNorm.Cross(zdir) , (planeNorm.Cross(zdir)).Cross(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, thit->cellID);
fitTrack.insertPoint(new genfit::TrackPoint(measurement, &fitTrack));
}
//check
fitTrack.checkConsistency();
// do the fit
fitter->processTrack(&fitTrack);
// print fit result
fitTrack.getFittedState().Print();
//check
fitTrack.checkConsistency();
//genfit::SharedPlanePtr test_plane(new genfit::DetPlane({0,0,0}, {1,0,0} , {0,1,0}) );
//std::vector<genfit::MeasurementOnPlane*> constructMeasurementsOnPlane ...
//delete fitter;
}
display->addEvent(all_tracks);
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