40 if ( p==
nullptr )
return nullptr;
42 if ( p->absPdgId()==pdg_id ) {
45 auto vertex = p->prodVtx();
46 if ( vertex ==
nullptr ) {
49 if ( vertex->nIncomingParticles() < 1 ) {
52 for(
unsigned ip = 0; ip < vertex->nIncomingParticles(); ip++ ) {
53 auto* in = vertex->incomingParticle(ip);
55 if ( parent!=
nullptr ) {
56 if (parent->absPdgId()==pdg_id)
return parent;
74 if ( p==
nullptr )
return nullptr;
76 for (
size_t i=ids.size() ; i-- ; ) {
77 if ( p->absPdgId()==ids[i] )
return p;
80 auto vertex = p->prodVtx();
81 if ( vertex ==
nullptr )
return nullptr;
83 if ( vertex->nIncomingParticles()<1 )
return nullptr;
85 for(
unsigned ip = 0; ip < vertex->nIncomingParticles(); ip++ ) {
86 auto* in = vertex->incomingParticle(ip);
88 if ( parent!=
nullptr ) {
89 for (
size_t i=ids.size() ; i-- ; ) {
90 if ( parent->absPdgId()==ids[i] )
return parent;
123 static const int hpmap[20] = { 0, 1, 2, 7, 8, 9, 3, 4, 5, 6, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 };
127#ifdef TRKPARAMETERS_MEASUREDPERIGEE_H
134 double pT = measPer->
pT();
135 double eta = measPer->
eta();
138 double d0 = measPer->parameters()[
Trk::d0];
144 if ( measPer->parameters()[
Trk::qOverP]<0 && pT>0 ) pT *= -1;
146#ifdef TRKPARAMETERS_MEASUREDPERIGEE_H
158 double deta = 0.5*dtheta/(std::cos(0.5*
theta)*std::cos(0.5*
theta)*std::tan(0.5*
theta));
160#ifdef TRKPARAMETERS_MEASUREDPERIGEE_H
162 double dz0 = err.error(
Trk::z0);
163 double dd0 = err.error(
Trk::d0);
166 double dz0 = std::sqrt((*measPer->covariance())(
Trk::z0,
Trk::z0));
167 double dd0 = std::sqrt((*measPer->covariance())(
Trk::d0,
Trk::d0));
173 double sintheta = std::sin(
theta);
174 double costheta = std::cos(
theta);
175 double dpt2 = (p*p*sintheta)*(p*p*sintheta)*dqovp*dqovp + (p*costheta)*(p*costheta)*dtheta*dtheta - 2*(p*p*sintheta)*(p*costheta)*covthetaOvP;
177 if ( dpt2>0 ) dpT = std::sqrt( dpt2 );
190 int nSiHits = nPixelHits + nSctHits;
191 bool expectBL =
false;
200 unsigned long id = (
unsigned long)track;
202 for (
int ih=0 ; ih<20 ; ih++ ) {
212 int fitter = track->info().trackFitter();
213 std::string dumpinfo = track->info().dumpInfo();
215 int trackAuthor = -1;
217 if ( dumpinfo.find(
"TRTStandalone")!=std::string::npos) trackAuthor = 2;
218 else if ( dumpinfo.find(
"TRTSeededTrackFinder")!=std::string::npos) trackAuthor = 1;
219 else trackAuthor = 0;
223 std::cout <<
"\t\t\tSUTT TP track"
229 <<
"\tNsi=" << nSiHits
230 <<
"\tNtrt=" << nTrHits
231 <<
"\tNstr=" << nStrawHits
232 <<
"\tauthor=" << trackAuthor
239 deta, dphi, dz0, dd0, dpT,
240 nBlayerHits, nPixelHits, nSctHits, nSiHits,
241 nStrawHits, nTrHits, bitmap, 0,
242 trackAuthor,
false, -1, -1,
266 while ( trackitr!=trackend ) {
285 while ( trackitr!=trackend ) {
311 std::vector<double> xpos(Nx,0);
312 std::vector<double> ypos(Ny,0);
315 std::vector<int> xn(Nx,0);
316 std::vector<int> yn(Ny,0);
321 double deltax = 3.0/Nx;
322 double deltay = 3.0/Ny;
329 for ( ; trackitr!=trackend ; ++trackitr ) {
333 if ( !
MC::isStable(track) || !track->hasProdVtx() )
continue;
337 double xp[3] = { track->prodVtx()->x(), track->prodVtx()->y(), track->prodVtx()->z() };
341 int ix = xp[0]/deltax + xoffset;
342 int iy = xp[1]/deltay + yoffset;
344 if ( ix<0 || ix>=Nx || iy<0 || iy>=Nx )
continue;
359 for (
size_t i=0 ; i<xpos.size() ; i++ ) {
360 if ( xn[i]>xn[imx] ) imx = i;
361 if ( yn[i]>yn[imy] ) imy = i;
367 if ( xn[imx]>1 ) x0 = xpos[imx]/xn[imx];
368 if ( yn[imy]>1 ) y0 = ypos[imy]/yn[imy];
386 for ( ; trackitr!=trackend; ++trackitr) {
390 double q = (*trackitr)->charge();
394 if ( q==-999 ) q = ptype.
charge( (*trackitr)->pdgId() );
401 bool gotPdgId =
true;
405 bool gotParentPdgId =
true;
410 if ( gotParentPdgId && gotPdgId )
selectTrack( *trackitr, x0, y0);
426 m_id = (
unsigned long)(track.get());
447 if ( t == 0 )
return false;
467 double pT = measPer->
pt();
468 double eta = measPer->
eta();
469 double phi = measPer->
phi();
472 if ( measPer->
charge()<0 && pT>0 ) pT *= -1;
473 double q = track->charge();
478 if ( q==-999 ) q = ptype.
charge( track->pdgId() );
481 if ( q==0 )
return 0;
490 if ( !track->hasProdVtx() )
return false;
495 double xb[3] = { xp[0]-x0, xp[1]-y0, measPer->
prodVtx()->
z() };
496 double xd[3] = { 0, 0, 0 };
498 if ( track->hasDecayVtx() ) {
499 xd[0] = track->decayVtx()->x();
500 xd[1] = track->decayVtx()->y();
501 xd[2] = track->decayVtx()->z();
504 double rp = std::sqrt( xp[0]*xp[0] + xp[1]*xp[1] );
505 double rd = std::sqrt( xd[0]*xd[0] + xd[1]*xd[1] );
509 double theta = 2*std::atan( std::exp( -
eta ) );
510 double z0 = xb[2] - (xb[0]*std::cos(
phi) + xb[1]*std::sin(
phi))/std::tan(
theta);
511 double d0 = xb[1]*std::cos(
phi) - xb[0]*std::sin(
phi);
513 bool final_state =
false;
526 const double inner_radius =
m_radius;
527 const double outer_radius =
m_radius;
529 if ( ( track->hasProdVtx() &&
rp<=inner_radius ) &&
530 ( !track->hasDecayVtx() || rd>outer_radius ) ) final_state =
true;
532 if ( !final_state )
return false;
550 bool expectBL =
false;
558 unsigned long id = (
unsigned long)track;
562 int trackAuthor = track->pdgId();
566 std::cout <<
"\t\t\tSUTT TP track"
572 <<
"\tauthor=" << trackAuthor
573 <<
"\tVTX x " << xp[0]<<
"\ty " << xp[1] <<
"\tz " << xp[2]
580 deta, dphi, dz0, dd0, dpT,
581 nBlayerHits, nPixelHits, nSctHits, nSiHits,
582 nStrawHits, nTrtHits, bitmap, 0,
583 trackAuthor,
false, uniqueID, -1,
605 unsigned long id = (
unsigned long)(track.get());
613 if ( track==0 )
return 0;
617 double phi = track->phi();
618 double eta = track->eta();
623 double xp[3] = { 0, 0, 0 };
625 if ( track->genParticle()->production_vertex() ) {
626 xp[0] = track->genParticle()->production_vertex()->position().x();
627 xp[1] = track->genParticle()->production_vertex()->position().y();
628 xp[2] = track->genParticle()->production_vertex()->position().z();
633 double theta = 2*std::atan( exp( (-1)*
eta ) );
634 double z0 = xp[2] - (xp[0]*std::cos(
phi) + xp[1]*std::sin(
phi))/std::tan(
theta);
636 double xd[3] = { 0, 0, 0 };
638 if ( track->genParticle()->end_vertex() ) {
639 xd[0] = track->genParticle()->end_vertex()->position().x();
640 xd[1] = track->genParticle()->end_vertex()->position().y();
641 xd[2] = track->genParticle()->end_vertex()->position().z();
644 double rp = std::sqrt( xp[0]*xp[0] + xp[1]*xp[1] );
645 double rd = std::sqrt( xd[0]*xd[0] + xd[1]*xd[1] );
648 bool final_state =
false;
661 const double inner_radius =
m_radius;
662 const double outer_radius =
m_radius;
663 if ( ( track->genParticle()->production_vertex() &&
rp<=inner_radius ) &&
664 ( track->genParticle()->end_vertex()==0 || rd>outer_radius ) ) final_state =
true;
667 if ( !final_state )
return 0;
680 double q = track->charge();
685 if ( q==-999 ) q = ptype.
charge( track->pdgId() );
688 if ( q==0 )
return 0;
690 double pT = q*track->pt();
705 d0 = xp[1]*std::cos(
phi) - xp[0]*std::sin(
phi);
720 int author = track->pdgId();
724 unsigned long id = (
unsigned long)track;
725 if ( tid!=0 )
id = tid;
749 author,
false, uniqueID, -1,
764 static const int hpmap[20] = { 0, 1, 2, 7, 8, 9, 3, 4, 5, 6, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 };
771#ifdef TRKPARAMETERS_MEASUREDPERIGEE_H
775 const Trk::Perigee* startPerigee = track->perigeeParameters();
793 double p = 1/startPerigee->parameters()[
Trk::qOverP];
796 double eta = startPerigee->
eta();
798 double z0 = startPerigee->parameters()[
Trk::z0];
799 double d0 = startPerigee->parameters()[
Trk::d0];
801 double pT = (1./qOverPt);
808#ifdef TRKPARAMETERS_MEASUREDPERIGEE_H
815 double dz0 = err.error(
Trk::z0);
816 double dd0 = err.error(
Trk::d0);
823 double dz0 = std::sqrt((*measPer->covariance())(
Trk::z0,
Trk::z0));
824 double dd0 = std::sqrt((*measPer->covariance())(
Trk::d0,
Trk::d0));
827 double deta = 0.5*dtheta/(std::cos(0.5*
theta)*std::cos(0.5*
theta)*std::tan(0.5*
theta));
835 double sintheta = std::sin(
theta);
836 double costheta = std::cos(
theta);
837 double dpT2 = (p*p*sintheta)*(p*p*sintheta)*dqovp*dqovp + (p*costheta)*(p*costheta)*dtheta*dtheta - 2*(p*p*sintheta)*(p*costheta)*covthetaOvP;
839 if ( dpT2>0 ) dpT = std::sqrt( dpT2 );
858 bool expectBL =
false;
862 std::cerr <<
"Could not create TrackSummary - Track will likely fail hits requirements" << std::endl;
870 nSiHits = nPixelHits + nSctHits;
872 for (
int ih=0 ; ih<20 ; ih++ ) {
877 unsigned long id = (
unsigned long)track;
882 if(quality==0) std::cerr <<
"Could not create FitQuality - Track will likely fail hits requirements" << std::endl;
888 int trackAuthor = -1;
896 int fitter = track->info().trackFitter();
900 if ((track->info().dumpInfo()).find(
"TRTStandalone") != std::string::npos) trackAuthor = 2;
901 else if ((track->info().dumpInfo()).find(
"TRTSeededTrackFinder") != std::string::npos) trackAuthor = 1;
902 else trackAuthor = 0;
906 std::cout <<
"\t\t\tSUTT TP track"
912 <<
"\tNsi=" << nSiHits
913 <<
"\tNtrt=" << nTrHits
914 <<
"\tNstr=" << nStrawHits
915 <<
"\tauthor=" << trackAuthor
920 deta, dphi, dz0, dd0, dpT,
921 nBlayerHits, nPixelHits, nSctHits, nSiHits,
922 nStrawHits, nTrHits, bitmap, 0,
923 trackAuthor,
false, -1, -1,
947 while ( trackitr!=trackend ) {
968 double pT = measPer->
pt();
969 double eta = measPer->
eta();
971 double z0 = measPer->
z0() + measPer->
vz();
972 double d0 = measPer->
d0();
981 if ( measPer->
qOverP()==0 )
throw std::runtime_error(
"probable corrupted track - this should never happen" );
982 double p = 1/measPer->
qOverP();
985 if ( measPer->
qOverP()<0 && pT>0 ) pT *= -1;
991 double deta = 0.5*dtheta/(std::cos(0.5*
theta)*std::cos(0.5*
theta)*std::tan(0.5*
theta));
1004 double sintheta = std::sin(
theta);
1005 double costheta = std::cos(
theta);
1006 double dpt2 = (p*p*sintheta)*(p*p*sintheta)*dqovp*dqovp + (p*costheta)*(p*costheta)*dtheta*dtheta - 2*(p*p*sintheta)*(p*costheta)*covthetaOvP;
1008 if ( dpt2>0 ) dpT = std::sqrt( dpt2 );
1015 uint8_t sum_nBlayerHits = 0;
1017 int nBlayerHits = 2*sum_nBlayerHits;
1019 uint8_t sum_nPixelHits = 0;
1021 int nPixelHits = 2*sum_nPixelHits;
1023 uint8_t sum_nSctHits = 0;
1025 int nSctHits = sum_nSctHits;
1027 uint8_t sum_nStrawHits = 0;
1029 int nStrawHits = sum_nStrawHits;
1031 uint8_t sum_nTrtHits = 0;
1033 int nTrtHits = sum_nTrtHits;
1036 uint8_t sum_expectBL = 0;
1038 bool expectBL = ( sum_expectBL ? true : false );
1042 uint8_t sum_sctholes = 0;
1045 uint8_t sum_pixholes = 0;
1052 nSctHits += 1000*sum_sctholes;
1053 nPixelHits += 1000*sum_pixholes;
1056 int nSiHits = nPixelHits + nSctHits;
1060 double chi2 = track->chiSquared();
1061 double dof = track->numberDoF();
1063 unsigned long id = (
unsigned long)track;
1065 unsigned bitmap = track->hitPattern();
1069 double xbeam = track->vx();
1070 double ybeam = track->vy();
1071 double zbeam = track->vz();
1075 int trackAuthor = 0;
1077 int fitter = track->trackFitter();
1078 std::bitset<xAOD::NumberOfTrackRecoInfo> patternrec = track->patternRecoInfo();
1081 for (
unsigned ipr=patternrec.size() ; ipr-- ; ) {
1082 if ( patternrec[ipr] ) {
1084 trackAuthor |= (ipr >> 16);
1093 trackAuthor |= fitter;
1102 std::cout <<
"\t\t\tSUTT TP track"
1108 <<
"\tNsi=" << nSiHits
1111 <<
"\tfitter=" << fitter
1112 <<
"\tauthor=" << trackAuthor
1113 <<
"\tVTX x " << track->vx() <<
"\ty " << track->vy() <<
"\tz " << track->vz()
1120 deta, dphi, dz0, dd0, dpT,
1121 nBlayerHits, nPixelHits, nSctHits, nSiHits,
1122 nStrawHits, nTrtHits, bitmap, 0,
1123 trackAuthor,
false, -1, -1,
1145 while ( trackitr!=trackend ) {
1158 while ( trackitr!=trackend ) {
1174 double& d0,
double& dd0,
Scalar eta() const
pseudorapidity method
Scalar phi() const
phi method
Scalar theta() const
theta method
double charge(const T &p)
ATLAS-specific HepMC functions.
DataVector< Trk::Track > TrackCollection
This typedef represents a collection of Trk::Track objects.
DataModel_detail::const_iterator< DataVector > const_iterator
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
virtual bool addTrack(TIDA::Track *t, bool(*f)(const TIDA::Track *)=0)
TrackSelector(TrackFilter *selector=0)
const std::vector< TIDA::Track * > & tracks() const
void setBeamline(double x, double y, double z=0)
const xAOD::TruthParticle * fromAncestor(const int pdg_id, const xAOD::TruthParticle *p) const
recursive functions to identify whether a particle comes from some particle of a specific PDG ID,...
void truthBeamline(const xAOD::TruthParticleContainer *truthtracks, double &x0, double &y0)
extract all the tracks from a xAOD::TruthParticle collection and histogram the x and y production coo...
TrigTrackSelector(TrackFilter *selector)
use a radius of 47 mm corresponding to the Run 1 pixel inner radius For the IBL it should be 32 mm,...
TIDA::Track * makeTrack(HepMC::ConstGenParticlePtr track)
void correctToBeamline(double &z0, double &dz0, double &d0, double &dd0, double theta, double phi)
bool selectTrack(const Rec::TrackParticle *track)
neater code to make use of vector function also for a single ancestor pdgid, instead of the full code...
static const double s_default_radius
NB: This was 47 for Run 2, but with the addition of the IBL it should be 32 It was kept at 47 for all...
void selectTracks(const Rec::TrackParticleContainer *trigtracks)
std::vector< ElementLink< xAOD::TrackParticleContainer > > TrackParticleLinks_t
Class to represent and store fit qualities from track reconstruction in terms of and number of degre...
int numberDoF() const
returns the number of degrees of freedom of the overall track or vertex fit as integer
double chiSquared() const
returns the of the overall track fit
double eta() const
Access method for pseudorapidity - from momentum.
double charge() const
Returns the charge.
double pT() const
Access method for transverse momentum.
@ NumberOfTrackFitters
maximum number of enums
A summary of the information contained by a track.
(HepMC) Monte Carlo particle.
double charge(int id) const
float z0() const
Returns the parameter.
float theta() const
Returns the parameter, which has range 0 to .
const ParametersCovMatrix_t definingParametersCovMatrix() const
Returns the 5x5 symmetric matrix containing the defining parameters covariance matrix.
float vz() const
The z origin for the parameters.
float d0() const
Returns the parameter.
float qOverP() const
Returns the parameter.
virtual double pt() const override final
The transverse momentum ( ) of the particle.
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
float phi0() const
Returns the parameter, which has range to .
virtual double pt() const override final
The transverse momentum ( ) of the particle.
const TruthVertex_v1 * prodVtx() const
The production vertex of this particle.
double charge() const
Physical charge.
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
virtual double phi() const override final
The azimuthal angle ( ) of the particle.
float z() const
Vertex longitudinal distance along the beam line form the origin.
float y() const
Vertex y displacement.
float x() const
Vertex x displacement.
double chi2(TH1 *h0, TH1 *h1)
HepMC3::ConstGenParticlePtr ConstGenParticlePtr
bool isElectron(const T &p)
bool isStable(const T &p)
Identify if the particle is stable, i.e. has not decayed.
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
DetectorType
enumerates the various detector types currently accessible from the isHit() method.
@ numberOfSCTHits
number of SCT holes
@ numberOfPixelHits
number of pixel layers on track with absence of hits
@ numberOfBLayerHits
these are the hits in the 0th pixel layer?
@ numberOfTRTHits
number of TRT outliers
@ numberOfTRTHighThresholdHits
total number of TRT hits which pass the high threshold
TrackParticle_v1 TrackParticle
Reference the current persistent version:
TruthParticle_v1 TruthParticle
Typedef to implementation.
TrackParticleContainer_v1 TrackParticleContainer
Definition of the current "TrackParticle container version".
@ expectInnermostPixelLayerHit
Do we expect a 0th-layer barrel hit for this track?
@ numberOfPixelHoles
number of pixel layers on track with absence of hits [unit8_t].
@ numberOfTRTHits
number of TRT hits [unit8_t].
@ numberOfSCTHits
number of hits in SCT [unit8_t].
@ numberOfInnermostPixelLayerHits
these are the hits in the 0th pixel barrel layer
@ numberOfPixelHits
these are the pixel hits, including the b-layer [unit8_t].
@ numberOfTRTHighThresholdHits
number of TRT hits which pass the high threshold (only xenon counted) [unit8_t].
@ numberOfSCTHoles
number of SCT holes [unit8_t].
TruthParticleContainer_v1 TruthParticleContainer
Declare the latest version of the truth particle container.