Skip to content
Snippets Groups Projects
Commit b8703f9b authored by Matteo Tenti's avatar Matteo Tenti
Browse files

kalman fitter (genfit) implemented in analysis.cpp; Makefile changed...

kalman fitter (genfit) implemented in analysis.cpp; Makefile changed accordingly; reconstruction.cpp modified to be more efficient
parent 43ca9e0f
No related branches found
No related tags found
No related merge requests found
......@@ -7,9 +7,13 @@ ROOTINCDIR = $(shell root-config --incdir)
EDEPSIM = /wd/sw/EDEPSIM/edep-sim.binary
EDEPGLIBS = -L$(EDEPSIM)/lib/ -ledepsim -ledepsim_io
EDEPGLIBS = -L$(EDEPSIM)/lib/ -ledepsim_io
EDEPINCDIR = $(EDEPSIM)/include/EDepSim
GENFIT=/wd/sw/GENFIT/GenFit.binary
EIGEN3=/usr/include/eigen3
all: Digitize Reconstruct Analyze
struct.cxx: include/struct.h include/Linkdef.h
......@@ -27,5 +31,5 @@ Reconstruct: libStruct.so
$(EDEPGLIBS) -Llib -lStruct src/reconstruction.cpp
Analyze: libStruct.so
g++ -o bin/$@ $(CFLAGS) $(LDFLAGS) -I$(EDEPINCDIR) -Iinclude $(ROOTGLIBS) -lGeom -lEG \
g++ -o bin/$@ $(CFLAGS) $(LDFLAGS) -I$(EDEPINCDIR) -Iinclude -I${GENFIT}/include -I${EIGEN3} -L${GENFIT}/lib64 -lgenfit2 $(ROOTGLIBS) -lGeom -lEG \
$(EDEPGLIBS) -Llib -lStruct src/analysis.cpp
......@@ -310,7 +310,7 @@ int main(int argc, char* argv[]) {
for(unsigned int j = 0; j < n; j++)
{
hitCoords[0] = evt->particles[k].tr.digits[j].x*mm2cm;
hitCoords[1] = evt->particles[k].tr.digits[j].y*mm2cm;;
hitCoords[1] = evt->particles[k].tr.digits[j].y*mm2cm;
genfit::PlanarMeasurement* measurement = new genfit::PlanarMeasurement(hitCoords, hitCov, detId, ++hitId, nullptr);
measurement->setPlane(genfit::SharedPlanePtr(new genfit::DetPlane(TVector3(0,0,evt->particles[k].tr.digits[j].z*mm2cm), TVector3(1,0,0), TVector3(0,1,0))), ++planeId);
fitTrack.insertPoint(new genfit::TrackPoint(measurement, &fitTrack));
......
......@@ -10,3 +10,6 @@ fi
export PATH=${PATH}:${DIR}/../bin
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:${DIR}/../lib
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/wd/sw/RAVE/rave-0.6.25.binary/lib
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/wd/sw/CLHEP/2.4.1.2/binary/lib/
export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:/wd/sw/GENFIT/GenFit.binary/lib64
......@@ -16,10 +16,17 @@ do
done
ifile=${1}
gfile=${2}
if [ "${ifile}" == "" ]
then
echo -e "\e[5m\e[91mERROR\e[0m: source ${BASH_SOURCE[0]} <input file>"
echo -e "\e[5m\e[91mERROR\e[0m: source ${BASH_SOURCE[0]} <input file> <geometry file>"
return 1
fi
if [ "${gfile}" == "" ]
then
echo -e "\e[5m\e[91mERROR\e[0m: source ${BASH_SOURCE[0]} <input file> <geometry file>"
return 1
fi
......@@ -57,4 +64,4 @@ Reconstruct "${rfile}"
echo "===============================" &&
echo "Analysis =" &&
echo "===============================" &&
Analyze "${rfile}"
Analyze "${rfile} ${gfile}"
......@@ -7,6 +7,18 @@
#include "TG4Event.h"
#include "TG4HitSegment.h"
#include <ConstField.h>
#include <FieldManager.h>
#include <KalmanFitterRefTrack.h>
#include <StateOnPlane.h>
#include <Track.h>
#include <TrackPoint.h>
#include <MaterialEffects.h>
#include <RKTrackRep.h>
#include <TGeoMaterialInterface.h>
#include <PlanarMeasurement.h>
#include <Exception.h>
#include "struct.h"
#include "utils.h"
......@@ -193,7 +205,20 @@ void CellXYZTE(cell c, double& x, double& y, double& z, double& t, double& e)
void RecoFromTrack(particle& p)
{
if(p.tr.ret_ln == 0 && p.tr.ret_cr == 0)
if(p.pdg > 10000 || !(p.tr.ret_ln == 0 && p.tr.ret_cr == 0))
{
p.charge_reco = 0.;
p.pxreco = 0.;
p.pyreco = 0.;
p.pzreco = 0.;
p.Ereco = 0.;
p.xreco = 0.;
p.yreco = 0.;
p.zreco = 0.;
p.treco = 0.;
return;
}
else
{
std::sort(p.tr.digits.begin(), p.tr.digits.end(), isDigBefore);
......@@ -219,6 +244,11 @@ void RecoFromTrack(particle& p)
l /= TMath::Abs(l);
const double mm2cm = 0.1;
const double cm2mm = 1./mm2cm;
const double GeV2MeV = 1000.;
const double MeV2GeV = 1./GeV2MeV;
double r0z = p.tr.z0 - p.tr.zc;
double r0y = p.tr.y0 - p.tr.yc;
......@@ -226,27 +256,105 @@ void RecoFromTrack(particle& p)
double ang_yz = TMath::ATan2(r0z, -r0y);
double ang_x = 0.5 * TMath::Pi() - TMath::ATan(1./p.tr.b);
p.charge_reco = -l;
p.pxreco = mom_yz * TMath::Tan(ang_x);
p.pyreco = p.charge_reco * mom_yz * TMath::Sin(ang_yz);
p.pzreco = p.charge_reco * mom_yz * TMath::Cos(ang_yz);
p.Ereco = TMath::Sqrt(p.pxreco*p.pxreco + p.pyreco*p.pyreco + p.pzreco*p.pzreco + p.mass*p.mass);
p.xreco = p.tr.x0;
p.yreco = p.tr.y0;
p.zreco = p.tr.z0;
p.treco = p.tr.t0;
}
else
{
p.charge_reco = 0.;
p.pxreco = 0.;
p.pyreco = 0.;
p.pzreco = 0.;
p.Ereco = 0.;
p.xreco = 0.;
p.yreco = 0.;
p.zreco = 0.;
p.treco = 0.;
double px_guess = mom_yz * TMath::Tan(ang_x) * MeV2GeV;
double py_guess = p.charge_reco * mom_yz * TMath::Sin(ang_yz) * MeV2GeV;
double pz_guess = p.charge_reco * mom_yz * TMath::Cos(ang_yz) * MeV2GeV;
double x_guess = p.tr.x0*mm2cm;
double y_guess = p.tr.y0*mm2cm;
double z_guess = p.tr.z0*mm2cm;
// init fitter
genfit::AbsKalmanFitter* fitter = new genfit::KalmanFitterRefTrack();
TVector3 pos(x_guess, y_guess, z_guess);
TVector3 mom(px_guess, py_guess, pz_guess);
// trackrep
genfit::AbsTrackRep* rep = new genfit::RKTrackRep(p.pdg);
// create track
genfit::Track fitTrack(rep, pos, mom);
const int detId(0); // detector ID
int planeId(0); // detector plane ID
int hitId(0); // hit ID
double detectorResolutionT(0.02); // resolution of planar detectors (tranverse) cm
double detectorResolutionL(2.); // resolution of planar detectors (longitudinal) cm
TMatrixDSym hitCovH(2);
hitCovH.UnitMatrix();
hitCovH[0][0] = 0.;
hitCovH[0][1] = detectorResolutionL*detectorResolutionL;
hitCovH[1][0] = 0.;
hitCovH[1][1] = detectorResolutionT*detectorResolutionT;
TMatrixDSym hitCovV(2);
hitCovV.UnitMatrix();
hitCovV[0][0] = 0.;
hitCovV[0][1] = detectorResolutionT*detectorResolutionT;
hitCovV[1][0] = 0.;
hitCovV[1][1] = detectorResolutionL*detectorResolutionL;
// add some planar hits to track with coordinates I just made up
TVectorD hitCoords(2);
TMatrixDSym* hitCov = 0;
const int n = p.tr.digits.size();
for(unsigned int j = 0; j < n; j++)
{
hitCoords[0] = p.tr.digits[j].x*mm2cm;
hitCoords[1] = p.tr.digits[j].y*mm2cm;
if(p.tr.digits[j].hor)
{
hitCov = &hitCovH;
}
else
{
hitCov = &hitCovV;
}
genfit::PlanarMeasurement* measurement = new genfit::PlanarMeasurement(hitCoords, *hitCov, detId, ++hitId, nullptr);
measurement->setPlane(genfit::SharedPlanePtr(new genfit::DetPlane(TVector3(0,0,p.tr.digits[j].z*mm2cm), TVector3(1,0,0), TVector3(0,1,0))), ++planeId);
fitTrack.insertPoint(new genfit::TrackPoint(measurement, &fitTrack));
}
//check
fitTrack.checkConsistency();
//fit
fitter->processTrack(&fitTrack);
if(fitTrack.getPointWithFitterInfo(0) != 0)
{
p.charge_reco = l;
p.pxreco = fitTrack.getFittedState().getMom().X()*GeV2MeV;
p.pyreco = fitTrack.getFittedState().getMom().Y()*GeV2MeV;
p.pzreco = fitTrack.getFittedState().getMom().Z()*GeV2MeV;
p.Ereco = TMath::Sqrt(p.pxreco*p.pxreco+p.pyreco*p.pyreco+p.pzreco*p.pzreco+p.mass*p.mass);
p.xreco = fitTrack.getFittedState().getPos().X()*cm2mm;
p.yreco = fitTrack.getFittedState().getPos().Y()*cm2mm;
p.zreco = fitTrack.getFittedState().getPos().Z()*cm2mm;
p.treco = p.tr.t0;
return;
}
else
{
p.charge_reco = l;
p.pxreco = px_guess;
p.pyreco = py_guess;
p.pzreco = pz_guess;
p.Ereco = TMath::Sqrt(p.pxreco*p.pxreco+p.pyreco*p.pyreco+p.pzreco*p.pzreco+p.mass*p.mass);
p.xreco = x_guess;
p.yreco = y_guess;
p.zreco = z_guess;
p.treco = p.tr.t0;
return;
}
}
}
......@@ -798,13 +906,13 @@ void EvalNuEnergy(event& ev)
}
}
void Analyze(const char* fIn)
void Analyze(const char* fIn, const char* fgeo)
{
TFile f(fIn,"UPDATE");
TTree* tReco = (TTree*) f.Get("tReco");
TTree* tTrueMC = (TTree*) f.Get("EDepSimEvents");
TTree* gRooTracker = (TTree*) f.Get("gRooTracker");
TGeoManager* geo = (TGeoManager*) f.Get("EDepSimGeometry");
//TGeoManager* geo = (TGeoManager*) f.Get("EDepSimGeometry");
tReco->AddFriend(tTrueMC);
tReco->AddFriend(gRooTracker);
......@@ -828,6 +936,13 @@ void Analyze(const char* fIn)
std::map<int, particle> map_part;
TGeoManager::Import(fgeo);
const double Bx = 6.; // kGauss
genfit::MaterialEffects::getInstance()->init(new genfit::TGeoMaterialInterface());
genfit::FieldManager::getInstance()->init(new genfit::ConstField(Bx, 0., 0.)); // kGauss
event evt;
TTree tout("tEvent","tEvent");
......@@ -916,15 +1031,15 @@ void Analyze(const char* fIn)
void help_ana()
{
std::cout << "Analyze <input file>" << std::endl;
std::cout << "Analyze <input file> <geometry file>" << std::endl;
}
int main(int argc, char* argv[])
{
if(argc != 2)
if(argc != 3)
help_ana();
else
Analyze(argv[1]);
Analyze(argv[1],argv[2]);
}
......@@ -66,47 +66,45 @@ double angle(double x1, double y1, double z1,
return TMath::ACos(prod/(mag1*mag2));
}
bool ishitok(TG4Event* ev, int trackid, TG4HitSegment hit,
bool ishitok_trj(TG4Trajectory trj, TG4HitSegment hit,
double postol = 5., double angtol = 0.3)
{
double x = 0.5*(hit.Start.X()+hit.Stop.X());
double y = 0.5*(hit.Start.Y()+hit.Stop.Y());
double z = 0.5*(hit.Start.Z()+hit.Stop.Z());
for(unsigned int kk = 0; kk < trj.Points.size()-1; kk++)
{
double dpos = mindist(trj.Points[kk].Position.X(),
trj.Points[kk].Position.Y(),
trj.Points[kk].Position.Z(),
trj.Points[kk+1].Position.X(),
trj.Points[kk+1].Position.Y(),
trj.Points[kk+1].Position.Z(),
x,y,z);
double dang = angle(trj.Points[kk+1].Position.X()-trj.Points[kk].Position.X(),
trj.Points[kk+1].Position.Y()-trj.Points[kk].Position.Y(),
trj.Points[kk+1].Position.Z()-trj.Points[kk].Position.Z(),
hit.Stop.X()-hit.Start.X(),
hit.Stop.Y()-hit.Start.Y(),
hit.Stop.Z()-hit.Start.Z());
if(dpos < postol && dang < angtol)
return true;
}
return false;
}
bool ishitok_tid(TG4Event* ev, int trackid, TG4HitSegment hit,
double postol = 5., double angtol = 0.3)
{
for(unsigned int jj = 0; jj < ev->Trajectories.size(); jj++)
{
//if(ev->Trajectories[jj].TrackId == hit.PrimaryId)
if(ev->Trajectories[jj].TrackId == trackid)
{
for(unsigned int kk = 0; kk < ev->Trajectories[jj].Points.size()-1; kk++)
{
double dpos = mindist(ev->Trajectories[jj].Points[kk].Position.X(),
ev->Trajectories[jj].Points[kk].Position.Y(),
ev->Trajectories[jj].Points[kk].Position.Z(),
ev->Trajectories[jj].Points[kk+1].Position.X(),
ev->Trajectories[jj].Points[kk+1].Position.Y(),
ev->Trajectories[jj].Points[kk+1].Position.Z(),
x,y,z);
double dang = angle(ev->Trajectories[jj].Points[kk+1].Position.X()-ev->Trajectories[jj].Points[kk].Position.X(),
ev->Trajectories[jj].Points[kk+1].Position.Y()-ev->Trajectories[jj].Points[kk].Position.Y(),
ev->Trajectories[jj].Points[kk+1].Position.Z()-ev->Trajectories[jj].Points[kk].Position.Z(),
hit.Stop.X()-hit.Start.X(),
hit.Stop.Y()-hit.Start.Y(),
hit.Stop.Z()-hit.Start.Z());
/*
std::cout << ev->Trajectories[jj].Points[kk].Position.X() << " " <<
ev->Trajectories[jj].Points[kk].Position.Y() << " " <<
ev->Trajectories[jj].Points[kk].Position.Z() << " " <<
ev->Trajectories[jj].Points[kk+1].Position.X() << " " <<
ev->Trajectories[jj].Points[kk+1].Position.Y() << " " <<
ev->Trajectories[jj].Points[kk+1].Position.Z() << " " <<
x << " " << y << " " << z << " " << dpos.back() << " " <<
ev->Trajectories[jj].TrackId << " " << ev->Trajectories[jj].PDGCode << std::endl;
*/
if(dpos < postol && dang < angtol)
return true;
}
return ishitok_trj(ev->Trajectories[jj], hit, postol, angtol);
}
}
return false;
......@@ -291,7 +289,7 @@ void TrackFind(TG4Event* ev, std::vector<digit>* vec_digi, std::vector<track>& v
//if(hseg.PrimaryId == tr.tid)
//{
//if(ishitok(ev, hseg->PrimaryId, hseg))
if(ishitok(ev, tr.tid, hseg))
if(ishitok_trj(ev->Trajectories.at(j), hseg))
{
tr.digits.push_back(vec_digi->at(k));
break;
......@@ -300,9 +298,12 @@ void TrackFind(TG4Event* ev, std::vector<digit>* vec_digi, std::vector<track>& v
}
}
std::sort(tr.digits.begin(), tr.digits.end(), isDigBefore);
vec_tr.push_back(tr);
if(tr.digits.size() > 0)
{
std::sort(tr.digits.begin(), tr.digits.end(), isDigBefore);
vec_tr.push_back(tr);
}
}
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment