Skip to content
Snippets Groups Projects
Commit fd00c3b3 authored by Andy Haas's avatar Andy Haas
Browse files

Closer to real

parent 7ccbdbf3
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Trigger/TrigFTK/FTK_RecTools/FTK_RecTools/FTK_DuplicateTrackRemovalTool.h
+ 2
− 87
View file @ fd00c3b3
......@@ -5,118 +5,33 @@
////////////////////////////////////////////////////////////////////////////////
// FTK_DuplicateTrackRemovalTool tool
// -------------------------------
// ATLAS Collaboration
//
// Remove duplicate (overlapping) tracks
//
// June 2017: Tool created
//
// Author: Andy Haas, NYU
// e-mail: ahaas@cern.ch
//
////////////////////////////////////////////////////////////////////////////////
#ifndef __TRIG_FTK_DUPLICATETRACKREMOVAL_TOOL_H__
#define __TRIG_FTK_DUPLICATETRACKREMOVAL_TOOL_H__
#include "GaudiKernel/ToolHandle.h"
#include "AthenaBaseComps/AthAlgTool.h"
#include "FTK_RecToolInterfaces/IFTK_DuplicateTrackRemovalTool.h"
#include "FTK_DataProviderInterfaces/IFTK_UncertaintyTool.h"
#include "TrkVxEdmCnv/IVxCandidateXAODVertex.h"
#include "TrigFTK_RawData/FTK_RawTrack.h"
#include "TrigFTK_RawData/FTK_RawTrackContainer.h"
#include "TrkTrack/TrackCollection.h"
#include "VxVertex/VxContainer.h"
//#include "xAODTracking/Vertex.h"
//#include "xAODTracking/TrackParticle.h"
//#include "xAODTracking/VertexContainer.h"
//#include "xAODTracking/TrackParticleContainer.h"
#include "xAODTracking/VertexFwd.h"
//#include "xAODTracking/TrackParticleFwd.h"
#include "xAODTracking/VertexContainerFwd.h"
#include "xAODTracking/VertexAuxContainer.h"
//#include "xAODTracking/TrackParticleContainerFwd.h"
//#include "Tracking/TrkVertexFitter/TrkVxEdmCnv/TrkVxEdmCnv/IVxCandidateXAODVertex.h"
class VxContainer;
using std::vector;
namespace Trk {
class Track;
//class Trk::IVxCandidateXAODVertex;
}
class FTK_DuplicateTrackRemovalTool : public AthAlgTool, virtual public IFTK_DuplicateTrackRemovalTool
{
//struct to hold track parameters
struct MyTrack{
MyTrack(int idx,float pt,float theta,float phi,float d0,float z0,float pterr,float thetaerr,float phierr,float d0err,float z0err){
m_idx=idx;
m_pt=pt;
m_theta=theta;
m_phi=phi;
m_d0=d0;
m_z0=z0;
m_pterr=pterr;
// m_qoverperr=qoverperr;
m_thetaerr=thetaerr;
m_phierr=phierr;
m_d0err=d0err;
m_z0err=z0err;
}
int m_idx;
double m_pt;
double m_theta;
double m_phi;
double m_d0;
double m_z0;
double m_pterr;
//double m_qoverperr=qoverperr;
double m_thetaerr;
double m_phierr;
double m_d0err;
double m_z0err;
//Track pT sort
bool operator<(const MyTrack &a)const{
return fabs(m_pt) > fabs(a.m_pt);
}
};
public:
FTK_DuplicateTrackRemovalTool( const std::string&, const std::string&, const IInterface* );
virtual ~FTK_DuplicateTrackRemovalTool(){};
virtual StatusCode initialize();
virtual StatusCode finalize ();
FTK_RawTrackContainer* removeDuplicates(const FTK_RawTrackContainer* trks);
private:
bool m_barrelOnly;
bool m_hasIBL;
double m_cluster_size;
double m_chi2cut;
double m_constTrkPt;
double m_constTrkEta;
double m_z0errfactor;
FTK_RawTrackContainer* m_trks_nodups;
std::pair<xAOD::VertexContainer*, xAOD::VertexAuxContainer*> findVertex(vector<MyTrack> trks);
double ctheta2eta(double cot);
vector<MyTrack> getTracks(const FTK_RawTrackContainer* trks);
vector<MyTrack> getTracks(const TrackCollection* trks);
ToolHandle<IFTK_UncertaintyTool> m_uncertaintyTool;
bool match(const FTK_RawTrack* track, const FTK_RawTrack* oldtrack) const;
};
#endif
Trigger/TrigFTK/FTK_RecTools/src/FTK_DuplicateTrackRemovalTool.cxx
+ 22
− 374
View file @ fd00c3b3
......@@ -2,63 +2,20 @@
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#include <cmath>
#include <iostream>
#include <vector>
#include "TVector3.h"
#include "GaudiKernel/MsgStream.h"
#include "FTK_RecTools/FTK_DuplicateTrackRemovalTool.h"
#include "FTK_DataProviderInterfaces/IFTK_UncertaintyTool.h"
#include "TrigFTK_RawData/FTK_RawTrack.h"
#include "TrkTrack/TrackCollection.h"
#include "VxVertex/VxContainer.h"
#include "VxVertex/ExtendedVxCandidate.h"
#include "VxVertex/RecVertex.h"
#include "VxVertex/Vertex.h"
#include "VxVertex/VxTrackAtVertex.h"
#include "TrigFTK_RawData/FTK_RawTrackContainer.h"
#include "TrkParameters/TrackParameters.h"
#include "TrkTrack/Track.h"
#include "xAODTracking/TrackParticleAuxContainer.h"
#include "xAODTracking/TrackParticleContainer.h"
#include "xAODTracking/VertexContainer.h"
#include "xAODTracking/VertexAuxContainer.h"
#include "VxVertex/VxCandidate.h"
#include "TrkVxEdmCnv/IVxCandidateXAODVertex.h"
#include <map>
#include <vector>
#include <utility>
//#include "TrkSurfaces/PlaneSurface.h"
using std::map; using std::string;
using std::cout; using std::endl;
using std::vector; using std::iterator;
using Trk::RecVertex;using Trk::VxTrackAtVertex;
FTK_DuplicateTrackRemovalTool::FTK_DuplicateTrackRemovalTool(const std::string& t,
const std::string& n,
const IInterface* p ):
AthAlgTool(t,n,p),
m_barrelOnly(false),
m_hasIBL(false),
m_cluster_size(1.87),
m_chi2cut(5.),
m_constTrkPt(1.),
m_constTrkEta(1.1),
m_z0errfactor(1.0),
m_trks_nodups(NULL),
m_uncertaintyTool("FTK_UncertaintyTool",this)
{
declareInterface< IFTK_DuplicateTrackRemovalTool >( this );
declareProperty( "HasIBL", m_hasIBL);
declareProperty( "BarrelOnly", m_barrelOnly);
declareProperty( "Cluster_size", m_cluster_size);
declareProperty( "Chi2cut", m_chi2cut);
declareProperty( "ConstTrkPt", m_constTrkPt);
declareProperty( "ConstTrkEta", m_constTrkEta);
declareProperty( "UncertaintyTool", m_uncertaintyTool);
declareProperty( "z0ErrFactor", m_z0errfactor);
}
StatusCode FTK_DuplicateTrackRemovalTool::initialize() {
......@@ -66,8 +23,7 @@ StatusCode FTK_DuplicateTrackRemovalTool::initialize() {
StatusCode sc = AlgTool::initialize();
MsgStream athenaLog(msgSvc(), name());
m_trks_nodups = new FTK_RawTrackContainer; //we DO own the tracks
//(SG::VIEW_ELEMENTS);//we don't own the tracks, we're just going to hold them
m_trks_nodups = new FTK_RawTrackContainer(SG::VIEW_ELEMENTS);//we don't own the tracks, we're just going to hold them
athenaLog << MSG::INFO << "FTK_DuplicateTrackRemovalTool initialized "<< endmsg;
return sc;
......@@ -75,345 +31,37 @@ StatusCode FTK_DuplicateTrackRemovalTool::initialize() {
StatusCode FTK_DuplicateTrackRemovalTool::finalize() {
StatusCode sc = AlgTool::finalize();
m_trks_nodups->clear();
delete m_trks_nodups;
return sc;
}
bool FTK_DuplicateTrackRemovalTool::match(const FTK_RawTrack* track, const FTK_RawTrack* oldtrack) const {
return true; // TODO
}
FTK_RawTrackContainer* FTK_DuplicateTrackRemovalTool::removeDuplicates(const FTK_RawTrackContainer* trks){
ATH_MSG_INFO("ACH99 - I'm in removeDuplicates!");
m_trks_nodups->clear();
m_trks_nodups->reserve(trks->size());
for (unsigned int i = 0; i!=trks->size(); i++) m_trks_nodups->push_back(new FTK_RawTrack(*(trks->at(i))));
return m_trks_nodups;
}
std::pair<xAOD::VertexContainer*, xAOD::VertexAuxContainer*> FTK_DuplicateTrackRemovalTool::findVertex(vector<MyTrack> mytrk)
{
MsgStream athenaLog(msgSvc(), name());
double mGkx=0.0001;
double mGky=0.0001;
double mGkz=16.;
athenaLog << MSG::DEBUG << "FTK_DuplicateTrackRemovalTool:: barrel " << m_barrelOnly<< endmsg;
xAOD::VertexContainer* theVertexContainer = new xAOD::VertexContainer;
xAOD::VertexAuxContainer* theVertexAuxContainer = new xAOD::VertexAuxContainer;
theVertexContainer->setStore( theVertexAuxContainer );
vector<MyTrack>::iterator trkbegin;
vector<MyTrack>::iterator trkend;
vector<MyTrack> trkatvtx;
athenaLog << MSG::DEBUG << "FTK_DuplicateTrackRemovalTool:: total " << mytrk.size()<< " tracks "<< endmsg;
sort(mytrk.begin(),mytrk.end());//sorting mytrk by pt(Track7.h)
int numdo=0;
do{
athenaLog << MSG::DEBUG << "findVertex: finding vertex "<< endmsg;
//////////////////////
//clustering
//////////////////////
numdo++;
vector<MyTrack> vxtrk;
//high pt track selecting
trkbegin=mytrk.begin();
trkend=mytrk.end();
float seed_z=(*trkbegin).m_z0;
//selecting tracks near high pt track
for(vector<MyTrack>::iterator i=trkbegin ; i<trkend;){
if(fabs((*i).m_z0 - seed_z) < m_cluster_size){
vxtrk.push_back(*i);
i=mytrk.erase(i);
}
else{
i++;
}
trkbegin=mytrk.begin();
trkend=mytrk.end();
}
//calculating seed z position(z average weighted by z error)
double oldz=0;
double z_weight=0;
trkbegin=vxtrk.begin();
trkend=vxtrk.end();
for(vector<MyTrack>::iterator i=trkbegin ; i<trkend;i++){
oldz += (*i).m_z0*(*i).m_z0err;
z_weight += (*i).m_z0err;
}
oldz /= z_weight;
double oldx=0;
double oldy=0;
double m_Gk[3][3]={{mGkx,0.,0.},{0.,mGky,0.},{0.,0.,mGkz}};//error
Amg::MatrixX C22_mat(3,3);
double chi2=0;
trkbegin=vxtrk.begin();
trkend=vxtrk.end();
int trknumber=0;
athenaLog << MSG::DEBUG << "FTK_DuplicateTrackRemovalTool:: fit VTX, using "<<vxtrk.size()<<" trks."<< endmsg;
for( vector<MyTrack>::iterator i=trkbegin ; i<trkend ;){
trknumber++;
double tmpchi2 = chi2;
athenaLog << MSG::VERBOSE << "getTracks: "<<trknumber<<", pt "<<(*i).m_pt<<", eta "<<(*i).m_theta<<", phi "<<(*i).m_phi<<", d0 "<<(*i).m_d0<<", z0 "<<(*i).m_z0<<", pt err "<<(*i).m_pterr<<", theta err "<<(*i).m_thetaerr<<", phi err "<< (*i).m_phierr<<", d0 err"<< (*i).m_d0err<<", z0 err "<<(*i).m_z0err<< endmsg;
//track parameter reading
double xv,yv,zv,P0,phi0,theta0;
TVector3 trk;
xv=oldx;
yv=oldy;
zv=oldz;
P0=(*i).m_pt;
phi0=(*i).m_phi;
theta0=(*i).m_theta;
//parameter of angle
double cosPhi0,sinPhi0,sinPsi,cosPsi;
double psi,ctt,sint;
double alpha=0.02997*20.84;///1000.0;
cosPhi0=cos(phi0);sinPhi0=sin(phi0);
sinPsi=-alpha*(xv*cosPhi0+yv*sinPhi0)/P0;
cosPsi=sqrt(1.0-sinPsi*sinPsi);
psi=asin(sinPsi);
sint=sin(theta0);
ctt=cos(theta0)/sint;
//difference of values expected and observed
double m_A[2][3],m_B[2][3],m_u[2],m_h[2],m_resid[2];
m_A[0][0]=-sin(phi0+psi)/cosPsi;
m_A[0][1]= cos(phi0+psi)/cosPsi;
m_A[0][2]=0.0;
m_A[1][0]=-ctt*cosPhi0/cosPsi;
m_A[1][1]=-ctt*sinPhi0/cosPsi;
m_A[1][2]=1.0;
m_B[0][0]=-xv*m_A[0][1]+yv*m_A[0][0];
m_B[0][1]=0.0;
m_B[0][2]=(1.0-1.0/cosPsi)/alpha;
m_B[1][0]=-xv*m_A[1][1]+yv*m_A[1][0];
m_B[1][1]=-P0*psi/(alpha*sint*sint);
m_B[1][2]=ctt*(psi-sinPsi/cosPsi)/alpha;
m_u[0]=(*i).m_d0;
m_u[1]=(*i).m_z0;
m_h[0]=yv*cosPhi0-xv*sinPhi0+P0*(1-cosPsi)/alpha;
m_h[1]=zv+P0*ctt*psi/alpha;
for (unsigned int i = 0; i!=trks->size(); i++) {
const FTK_RawTrack *track = trks->at(i);
m_resid[0]=m_u[0]-m_h[0];
m_resid[1]=m_u[1]-m_h[1];
//now we should see whether this track overlaps with one (or more?) tracks already in the nodups container
std::vector<int> matching_oldtracks;
for (unsigned int e = 0; e!=m_trks_nodups->size(); e++) {
const FTK_RawTrack *oldtrack = m_trks_nodups->at(e);
if (match(track,oldtrack)) {
matching_oldtracks.push_back(e);
}
}
//error of difference
//if it does, either replace the (worst?) matching track with this new track, or ignore this new track, depending on which track we like better
double m_Vqq[3][3]={{(*i).m_phierr*(*i).m_phierr,0.,0.},{0.,(*i).m_thetaerr*(*i).m_thetaerr,0.},{0.,0.,(*i).m_pterr*(*i).m_pterr}};
double m_Vuu[2][2]={{(*i).m_d0err*(*i).m_d0err,0.},{0.,(*i).m_z0err*(*i).m_z0err}};
double AC[2][3],BV[2][3],Sk[2][2];
for(int i1=0;i1<2;i1++) for(int j=0;j<2;j++) Sk[i1][j]=m_Vuu[i1][j];
for(int i1=0;i1<2;i1++) for(int j=0;j<3;j++)
{
AC[i1][j]=0.0;
for(int k=0;k<3;k++) AC[i1][j]+=m_A[i1][k]*m_Gk[j][k];
}
for(int i1=0;i1<2;i1++) for(int j=0;j<2;j++)
{
for(int k=0;k<3;k++) Sk[i1][j]+=AC[i1][k]*m_A[j][k];
}
for(int i1=0;i1<2;i1++)
for(int j=0;j<3;j++)
{
BV[i1][j]=0.0;
for(int k=0;k<3;k++) BV[i1][j]+=m_B[i1][k]*m_Vqq[k][j];
}
for(int i1=0;i1<2;i1++)
for(int j=0;j<2;j++)
{
for(int k=0;k<3;k++) Sk[i1][j]+=BV[i1][k]*m_B[j][k];
}
//error determinant
double detr,m_V[2][2];
detr=1.0/(Sk[0][0]*Sk[1][1]-Sk[0][1]*Sk[1][0]);
m_V[0][0]=Sk[1][1]*detr;
m_V[1][1]=Sk[0][0]*detr;
m_V[0][1]=m_V[1][0]=-Sk[0][1]*detr;
//chi2
chi2=m_V[0][0]*m_resid[0]*m_resid[0]+m_V[1][1]*m_resid[1]*m_resid[1]+2.0*m_V[0][1]*m_resid[1]*m_resid[0];
if(chi2>m_chi2cut || chi2<0){
i=vxtrk.erase(i);
chi2 = tmpchi2;
}
else{
trkatvtx.push_back(*i);
i++;
}
trkbegin=vxtrk.begin();
trkend=vxtrk.end();
//vertex position & covariance update
double K[2][3],dV[3];
for(int i2=0;i2<2;i2++){
for(int j=0;j<3;j++){
for(int k=0;k<2;k++)K[i2][j]=AC[k][j]*m_V[k][i2];
}
}
for(int i2=0;i2<3;i2++){
dV[i2]=K[0][i2]*m_resid[0]+K[1][i2]*m_resid[1];
for(int j2=i2;j2<3;j2++){
m_Gk[i2][j2]-=K[0][i2]*AC[0][j2]+K[1][i2]*AC[1][j2];
m_Gk[j2][i2]=m_Gk[i2][j2];
}
}
oldx+=dV[0];
oldy+=dV[1];
oldz+=dV[2];
C22_mat(0,0)=m_Gk[0][0];C22_mat(0,1)=m_Gk[0][1];C22_mat(0,2)=m_Gk[0][0];
C22_mat(1,0)=m_Gk[1][0];C22_mat(1,1)=m_Gk[1][1];C22_mat(1,2)=m_Gk[1][2];
C22_mat(2,0)=m_Gk[2][0];C22_mat(2,1)=m_Gk[2][1];C22_mat(2,2)=m_Gk[2][2];
}//track loop end
if (vxtrk.size()==0)continue;
athenaLog << MSG::DEBUG << "findVertex: find vertex at "<< oldx<<";"<<oldy<<";"<<oldz<< endmsg;
Amg::Vector3D frameOrigin(oldx,oldy,oldz);
int ndf = 2*vxtrk.size()-3;
xAOD::Vertex* vx = new xAOD::Vertex;
vx->makePrivateStore();
vx->setPosition (frameOrigin);
vx->setCovariancePosition (C22_mat);
vx->setFitQuality(chi2,static_cast<float>(ndf));
vx->setVertexType(xAOD::VxType::PriVtx);
std::vector<VxTrackAtVertex> & tracksAtVertex = vx->vxTrackAtVertex(); tracksAtVertex.clear();
// Store the tracks at vertex
Amg::Vector3D Vertex(frameOrigin[0],frameOrigin[1],frameOrigin[2]);
const Trk::PerigeeSurface Surface(Vertex);
Trk::Perigee * refittedPerigee(0);
vector<MyTrack>::iterator i;
for (i = trkatvtx.begin(); i != trkatvtx.end() ; ++i)
{
AmgSymMatrix(5) * CovMtxP = new AmgSymMatrix(5);
double eta = ctheta2eta((*i).m_theta);
double etaerr=(*i).m_thetaerr*cosh(eta);
double qoverperr=cosh(eta)? 2./cosh(eta) * sqrt(4*(*i).m_pterr*(*i).m_pterr + 1/4/(*i).m_pt/(*i).m_pt*tanh(eta)*tanh(eta)*etaerr*etaerr) : 10;
(*CovMtxP)(0,0)=(*i).m_d0err*(*i).m_d0err;
(*CovMtxP)(0,0)=(*i).m_z0err*(*i).m_z0err;
(*CovMtxP)(0,0)=(*i).m_phierr*(*i).m_phierr;
(*CovMtxP)(0,0)=(*i).m_thetaerr*(*i).m_thetaerr;
(*CovMtxP)(0,0)=qoverperr*qoverperr;
double qoverp = 2/(*i).m_pt/cosh(eta);
refittedPerigee = new Trk::Perigee ((*i).m_d0,(*i).m_z0,(*i).m_phi,(*i).m_theta,qoverp, Surface, CovMtxP);
tracksAtVertex.emplace_back(-1, refittedPerigee);
}
theVertexContainer->push_back(vx);
trkatvtx.clear();
//delete vx;
}while(mytrk.size()>0);//vertex loop end
for (auto pv = theVertexContainer->begin(); pv != theVertexContainer->end(); ++pv) {
athenaLog << MSG::DEBUG << "vertex x = " <<(*pv)->position().x()<< endmsg;
m_trks_nodups->push_back((FTK_RawTrack*)track);
}
return std::make_pair(theVertexContainer, theVertexAuxContainer);
//return theVertexContainer;
}
//maybe at the end we should do a check that no further matches exist?
vector<FTK_DuplicateTrackRemovalTool::MyTrack> FTK_DuplicateTrackRemovalTool::getTracks(const FTK_RawTrackContainer* trks){
MsgStream athenaLog(msgSvc(), name());
vector<MyTrack> mytrk;
athenaLog << MSG::DEBUG << "getRawTracks: find "<< trks->size()<< " tracks "<< endmsg;
for (unsigned int ftk_track_index = 0; ftk_track_index != trks->size(); ++ftk_track_index){
const FTK_RawTrack* ftk_track= trks->at(ftk_track_index);
float trk_theta = std::atan2(1.0,ftk_track->getCotTh());
double invpt= ftk_track->getInvPt();
float trk_eta = -std::log(std::tan(trk_theta/2));
float trk_phi = ftk_track->getPhi();
float trk_d0 = ftk_track->getD0();
float trk_z0 = ftk_track->getZ0();
float trk_pt = 1/invpt;
/* double trk_phierr=sqrt(3.446e-7+29.24*invpt*invpt);
double trk_thetaerr=sqrt(7.093e-6+38.4*invpt*invpt);
double trk_pterr=1;
double trk_d0err=sqrt(1.662e-3+6.947e4*invpt*invpt);
double trk_z0err=sqrt(6.472e-2+1.587e5*invpt*invpt);
*/
double trk_d0err= sqrt(m_uncertaintyTool->getParamCovMtx(*ftk_track, m_hasIBL, FTKTrackParam::d0, FTKTrackParam::d0));
double trk_z0err= sqrt(m_uncertaintyTool->getParamCovMtx(*ftk_track, m_hasIBL, FTKTrackParam::z0, FTKTrackParam::z0))*m_z0errfactor;
double trk_phierr= sqrt(m_uncertaintyTool->getParamCovMtx(*ftk_track, m_hasIBL, FTKTrackParam::phi, FTKTrackParam::phi));
double trk_thetaerr= sqrt(m_uncertaintyTool->getParamCovMtx(*ftk_track, m_hasIBL, FTKTrackParam::theta, FTKTrackParam::theta));
double trk_pterr= sqrt(m_uncertaintyTool->getParamCovMtx(*ftk_track, m_hasIBL, FTKTrackParam::pt, FTKTrackParam::pt));
if(fabs(trk_pt)<m_constTrkPt)continue;
if(m_barrelOnly&&fabs(trk_eta)>m_constTrkEta)continue;
if(ftk_track_index<10)athenaLog << MSG::DEBUG << "getRawTracks: get pt "<< trk_pt<<"+/-"<<trk_pterr << ", eta "<<trk_eta<<"+/-"<<trk_thetaerr<<", phi"<<trk_phi<<"+/-"<<trk_phierr<<", d0"<<trk_d0<<"+/-"<<trk_d0err<<", z0"<<trk_z0<<"+/-"<<trk_z0err <<endmsg;
MyTrack Trk(ftk_track_index,trk_pt,trk_theta,trk_phi,trk_d0,trk_z0,trk_pterr,trk_thetaerr,trk_phierr,trk_d0err,trk_z0err);
mytrk.push_back(Trk);
}
athenaLog << MSG::DEBUG << "getRawTracks: total "<< mytrk.size() << " track "<< endmsg;
return mytrk;
}
vector<FTK_DuplicateTrackRemovalTool::MyTrack> FTK_DuplicateTrackRemovalTool::getTracks(const TrackCollection* trks){
MsgStream athenaLog(msgSvc(), name());
vector<MyTrack> mytrk;
athenaLog << MSG::DEBUG << "getRefitTracks: find "<< trks->size()<< " tracks "<< endmsg;
TrackCollection::const_iterator track_it = trks->begin();
TrackCollection::const_iterator last_track = trks->end();
for (int iTrack = 0 ; track_it!= last_track; track_it++, iTrack++) {
// double invpt= sin(trk_Theta)*ftk_track->getCurv()/2.;
double trk_qOverp= (*track_it)->perigeeParameters()->parameters()[Trk::qOverP];
float trk_theta = (*track_it)->perigeeParameters()->parameters()[Trk::theta];
float trk_eta = (*track_it)->perigeeParameters()->eta();//ctheta2eta(trk_theta);
float trk_phi = (*track_it)->perigeeParameters()->parameters()[Trk::phi0];
float trk_d0 = (*track_it)->perigeeParameters()->parameters()[Trk::d0];
float trk_z0 = (*track_it)->perigeeParameters()->parameters()[Trk::z0];
float trk_pt = 1/(trk_qOverp/sin(trk_theta));
/*
float invpt = 1/(trk_pt);
double trk_phierr=sqrt(3.446e-7+29.24*invpt*invpt);
double trk_thetaerr=sqrt(7.093e-6+38.4*invpt*invpt);
double trk_pterr=1;
double trk_d0err=sqrt(1.662e-3+6.947e4*invpt*invpt);
double trk_z0err=sqrt(6.472e-2+1.587e5*invpt*invpt)*m_z0errfactor;
*/
const AmgSymMatrix(5) *cov=(*track_it)->perigeeParameters()->covariance();
double trk_d0err=sqrt((*cov)(0,0));
double trk_z0err=sqrt((*cov)(1,1));
double trk_phierr=sqrt((*cov)(2,2));
double trk_thetaerr=sqrt((*cov)(3,3));
double trk_qOverperr=sqrt((*cov)(4,4));
/* double trk_d0err=(*track_it)->measurementsOnTrack()->localCovariance().error(Trk::d0);
double trk_z0err=(*track_it)->measurementsOnTrack()->localErrorMatrix().error(Trk::z0);
double trk_phi0err=(*track_it)->measurementsOnTrack()->localErrorMatrix().error(Trk::phi0);
double trk_thetaerr=(*track_it)->measurementsOnTrack()->localErrorMatrix().error(Trk::theta);
double trk_qOverperr=(*track_it)->measurementsOnTrack()->localErrorMatrix().error(Trk::qOverP);
*/ double trk_pterr=-sin(trk_theta)*trk_qOverperr/trk_qOverp/trk_qOverp+cos(trk_theta)*trk_thetaerr/trk_qOverp;
if(fabs(trk_pt)<m_constTrkPt)continue;
if(m_barrelOnly&&fabs(trk_eta)>m_constTrkEta)continue;
if(iTrack<10)athenaLog << MSG::DEBUG << "getRefitTracks: get pt "<< trk_pt<<"+/-"<<trk_qOverperr << ", eta "<<trk_eta<<"+/-"<<trk_thetaerr<<", phi"<<trk_phi<<"+/-"<<trk_phierr<<", d0"<<trk_d0<<"+/-"<<trk_d0err<<", z0"<<trk_z0<<"+/-"<<trk_z0err <<endmsg;
mytrk.push_back(MyTrack(iTrack,trk_pt,trk_theta,trk_phi,trk_d0,trk_z0,trk_pterr,trk_thetaerr,trk_phierr,trk_d0err,trk_z0err));
}
return mytrk;
}
double FTK_DuplicateTrackRemovalTool::ctheta2eta(double cot){
double eta=0.0;
double tempX=0.0;
if(atan(1.0/cot) >= 0) tempX = atan(1.0/cot);
if(atan(1.0/cot) < 0) tempX = atan(1.0/cot)+M_PI;
eta = -log(tan(tempX/2.0));
return eta;
return m_trks_nodups;
}
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