TrigTrackSelector.cxx

Go to the documentation of this file.

 
 
#include "TrigInDetAnalysisUtils/TrigTrackSelector.h"
#include "TruthUtils/HepMCHelpers.h"
 
#include "xAODTruth/TruthVertexContainer.h"
 
#include <stdexcept>

const double TrigTrackSelector::s_default_radius = 47;
 
 
TrigTrackSelector::TrigTrackSelector( TrackFilter* selector, double radius, int selectPdgId, int selectParentPdgId ) : 
    TrackSelector(selector), m_id(0), m_xBeam(0), m_yBeam(0), m_zBeam(0),
    m_correctTrkTracks(false), 
    m_radius(radius), m_selectPdgId(selectPdgId), m_selectParentPdgId(selectParentPdgId) { 
} 
 
 

 
const xAOD::TruthParticle* TrigTrackSelector::fromAncestor(const int pdg_id,  const xAOD::TruthParticle *p) const { 
  if ( p==nullptr ) return nullptr;
  if ( MC::isElectron(p) || MC::isMuon(p) ) return nullptr; //don't want light leptons from eg tau decays - they are found directly
  if ( p->absPdgId()==pdg_id ) {
    return p;   // recursive stopping conditions
  }
  auto vertex = p->prodVtx();
  if ( vertex == nullptr ) {
    return nullptr; // has no production vertex !!!
  }
  if ( vertex->nIncomingParticles() < 1 ) {
    return nullptr;  // recursive stopping conditions
  }
  for( unsigned ip = 0; ip < vertex->nIncomingParticles(); ip++ ) {
    auto* in = vertex->incomingParticle(ip);
    auto parent = fromAncestor( pdg_id, in);
    if ( parent!=nullptr ) { 
      if (parent->absPdgId()==pdg_id) return parent;
    }
  }
  
  return nullptr;
}
 

 
const xAOD::TruthParticle* TrigTrackSelector::fromAncestor( const std::vector<int>& ids,  const xAOD::TruthParticle *p ) const { 
  if ( p==nullptr ) return nullptr;
  if (MC::isElectron(p) || MC::isMuon(p)) return nullptr; //don't want light leptons from eg tau decays - they are found directly
  for ( size_t i=ids.size() ; i-- ; ) { 
    if ( p->absPdgId()==ids[i] ) return p; // recursive stopping conditions
  }
 
  auto vertex = p->prodVtx();
  if ( vertex == nullptr ) return nullptr; // has no production vertex !!!
 
  if ( vertex->nIncomingParticles()<1 ) return nullptr;  // recursive stopping conditions
 
  for( unsigned ip = 0; ip < vertex->nIncomingParticles(); ip++ ) {
    auto* in = vertex->incomingParticle(ip);
    auto parent = fromAncestor( ids, in);
    if ( parent!=nullptr ) { 
      for ( size_t i=ids.size() ; i-- ; ) { 
    if ( parent->absPdgId()==ids[i] ) return parent; 
      }
    }
  }  
 
  return nullptr;
}
 
 

 
// const xAOD::TruthParticle* TrigTrackSelector::fromAncestor(const int pdg_id,  const xAOD::TruthParticle *p) const { 
//   return fromAncestor( std::vector<int>(1,pdg_id), p );
// }
 
 
// extract all the tracks from a TruthParticle collection and add them
void TrigTrackSelector::selectTracks( const TruthParticleContainer* truthtracks ) { 
    //    std::cout << "\t\t\tSUTT \tTrackParticleContainer->size() = " << trigtracks->size() << std::endl;
    
    TruthParticleContainer::const_iterator  trackitr = truthtracks->begin();
    TruthParticleContainer::const_iterator  trackend = truthtracks->end();
   
    while ( trackitr!=trackend ) { 
    
      selectTrack( *trackitr );
      
      ++trackitr;
 
    } // loop over tracks
    
}
 

 
void TrigTrackSelector::truthBeamline( const xAOD::TruthParticleContainer* truthtracks, double& x0, double& y0 ) { 
 
  x0 = 0;
  y0 = 0;
    
 
  int Nx = 300;
  int Ny = 300;

  std::vector<double> xpos(Nx,0);
  std::vector<double> ypos(Ny,0);

  std::vector<int> xn(Nx,0);
  std::vector<int> yn(Ny,0);
 
  int xoffset = Nx/2;; 
  int yoffset = Ny/2; 
 
  double deltax = 3.0/Nx;
  double deltay = 3.0/Ny;

 
  xAOD::TruthParticleContainer::const_iterator  trackitr = truthtracks->begin();
  xAOD::TruthParticleContainer::const_iterator  trackend = truthtracks->end();
    
  for ( ; trackitr!=trackend ; ++trackitr ) { 
    
    const xAOD::TruthParticle* track = (*trackitr);
    
    if ( !MC::isStable(track) || !track->hasProdVtx() ) continue; 

 
    double xp[3] = { track->prodVtx()->x(), track->prodVtx()->y(), track->prodVtx()->z() };

 
    int ix = xp[0]/deltax + xoffset;
    int iy = xp[1]/deltay + yoffset;
 
    if ( ix<0 || ix>=Nx || iy<0 || iy>=Nx ) continue;
 
    xpos[ix] += xp[0];
    ypos[iy] += xp[1];
 
    xn[ix]++;
    yn[iy]++;
 
  } // loop over tracks

 
  int imx = 0; 
  int imy = 0; 
 
  for ( size_t i=0 ; i<xpos.size() ; i++ ) {
    if ( xn[i]>xn[imx] ) imx = i;
    if ( yn[i]>yn[imy] ) imy = i;
  }

  if ( xn[imx]>1 ) x0 = xpos[imx]/xn[imx];
  if ( yn[imy]>1 ) y0 = ypos[imy]/yn[imy];
 
}
 
 
// extract all the tracks from a xAOD::TruthParticle collection and add them
void TrigTrackSelector::selectTracks( const xAOD::TruthParticleContainer* truthtracks) { 
 
  xAOD::TruthParticleContainer::const_iterator  trackitr = truthtracks->begin();
  xAOD::TruthParticleContainer::const_iterator  trackend = truthtracks->end();
  
  double x0 = 0;
  double y0 = 0;
 
  truthBeamline( truthtracks, x0, y0 );
 
 
  for ( ; trackitr!=trackend; ++trackitr) {
 
 
    // Only select charged final state particles
    double q = (*trackitr)->charge();

    static const particleType ptype;
    if ( q==-999 ) q = ptype.charge( (*trackitr)->pdgId() );
 
    if (q == 0 || !MC::isStable(*trackitr) ) continue;
 
    // If looking for tau parents, don't select mu or e children
 
    // select based on the pdg of final state particle 
    bool gotPdgId = true;
    if (m_selectPdgId!=0) gotPdgId = (*trackitr)->absPdgId()==m_selectPdgId;
 
    // select based on the pdg of the parent or ancestor
    bool gotParentPdgId = true;
    if   ( gotPdgId && m_selectParentPdgId!=0 )         gotParentPdgId = fromAncestor(m_selectParentPdgId, (*trackitr))!=nullptr;
 
    if ( gotParentPdgId && gotPdgId ) selectTrack( *trackitr, x0, y0);
 
  } // loop over tracks
    
}
 
 
 
// add a TruthParticle from a GenParticle - easy, bet it doesn't work 
bool TrigTrackSelector::selectTrack( HepMC::ConstGenParticlePtr track ) {
  
    if ( !MC::isStable(track) ) return false;

//AV Using memory to get some value is not a good idea. This is not a repruducible/portable way, but I leave it as is.
    m_id = (unsigned long)(track.get());
    bool sel;
    sel = selectTrack( TruthParticle(track) );
    m_id = 0;
    
    return sel; 
    
}
 
 
// add a TruthParticle 
bool TrigTrackSelector::selectTrack( const TruthParticle& track ) {
 
  return selectTrack( &track );
  
}
 
 
// add a TruthParticle 
bool TrigTrackSelector::selectTrack( const TruthParticle* track ) { 
    TIDA::Track* t = makeTrack( track, m_id );
    if ( t == 0 ) return false;
    if ( !addTrack(t) ) {
      delete t;
      return false;
    }
    return true;
}
 
 
 
 
// add an xAOD::TruthParticle 
bool TrigTrackSelector::selectTrack( const xAOD::TruthParticle* track, double x0, double y0) { 
  if ( track ) { 
        
    if (!MC::isStable(track) ) return false;

    const xAOD::TruthParticle* measPer = track;
 
    double pT    = measPer->pt(); 
    double eta   = measPer->eta();
    double phi   = measPer->phi();
  
    // AAARCH!!!!! the TrackParticle pT is NOT SIGNED!!!! ( I ask you! ) 
    if ( measPer->charge()<0 && pT>0 ) pT *= -1;
    double q = track->charge();
 
    static const particleType ptype;

    if ( q==-999 ) q = ptype.charge( track->pdgId() );
    
    if ( q==0 ) return 0;

    //      double xbeam = getBeamX(); // track->vx(); 
    //      double ybeam = getBeamY(); // track->vy(); 
    //      double zbeam = getBeamZ(); // track->vz(); 
      
    if ( !track->hasProdVtx() ) return false; 

 
    double xp[3] = { measPer->prodVtx()->x(), measPer->prodVtx()->y(), measPer->prodVtx()->z() };
    double xb[3] = { xp[0]-x0, xp[1]-y0, measPer->prodVtx()->z() };
    double xd[3] = { 0, 0, 0 };
      
    if ( track->hasDecayVtx() ) { 
      xd[0] = track->decayVtx()->x();
      xd[1] = track->decayVtx()->y();
      xd[2] = track->decayVtx()->z();
    }
      
    double rp = std::sqrt( xp[0]*xp[0] + xp[1]*xp[1] ); 
    double rd = std::sqrt( xd[0]*xd[0] + xd[1]*xd[1] ); 
      
    double theta = 2*std::atan( std::exp( -eta ) );
    double z0 = xb[2] - (xb[0]*std::cos(phi) + xb[1]*std::sin(phi))/std::tan(theta);
    double d0 = xb[1]*std::cos(phi) -  xb[0]*std::sin(phi);
 
    bool final_state = false; 
      
    const double inner_radius = m_radius; 
    const double outer_radius = m_radius;
 
    if ( (  track->hasProdVtx() && rp<=inner_radius ) && 
     ( !track->hasDecayVtx() || rd>outer_radius ) ) final_state = true; 
      
    if ( !final_state ) return false; 
    
    double deta = 0;
    double dphi = 0;
    double dz0  = 0;
    double dd0  = 0;
 
    double dpT  = 0; 
 
    int nBlayerHits = 0;
    int nPixelHits  = 0;
    int nSctHits    = 0;
    int nStrawHits  = 0;
    int nTrtHits    = 0;
      
    double chi2 = 0;
    double dof  = 0;
 
    bool expectBL = false;
 
    nSctHits   += 0;
    nPixelHits += 0;

    int nSiHits     = 0;
 
    unsigned long id = (unsigned long)track;
 
    unsigned bitmap = 0;
 
    int trackAuthor = track->pdgId();
    int uniqueID     = HepMC::uniqueID(track);
 
#if 0
    std::cout << "\t\t\tSUTT TP track" 
          << "\teta=" << eta  
          << "\tphi=" << phi  
          << "\tz0="  << z0 
          << "\tpT="  << pT 
          << "\td0="  << d0
          << "\tauthor=" << trackAuthor
          << "\tVTX x " << xp[0]<< "\ty " << xp[1] << "\tz " << xp[2] 
          << std::endl;
#endif
 
    // Create and save Track
      
    TIDA::Track* t = new TIDA::Track( eta,   phi,  z0,  d0,  pT, chi2, dof,
                      deta,  dphi, dz0, dd0, dpT,
                      nBlayerHits, nPixelHits, nSctHits, nSiHits,
                      nStrawHits,  nTrtHits,   bitmap, 0,
                      trackAuthor,  false, uniqueID, -1,  
                      expectBL, id) ;  

    //      std::cout << "SUTT TP track " << *t << "\t0x" << std::hex << bitmap << std::dec << std::endl; 
      
    // addTrack applies additional cuts using the Filter
    if ( !addTrack( t ) ){
      delete t;
      return false;
    }
  }
  return false;
     
}
 
 
 
 
// make a TIDA::Track from a GenParticle 
TIDA::Track* TrigTrackSelector::makeTrack(HepMC::ConstGenParticlePtr track ) { 
//AV Using memory to get some value is not a good idea. This is not a repruducible/portable way, but I leave it as is.
    unsigned long id = (unsigned long)(track.get());
    TruthParticle t = TruthParticle(track); 
    return  makeTrack( &t, id );
}
 
// make a TIDA::Track from a TruthParticle 
TIDA::Track* TrigTrackSelector::makeTrack( const TruthParticle* track, unsigned long tid ) { 
  
    if ( track==0 ) return 0; 
    if ( !MC::isStable(track)) return 0;   
 
 
    double phi = track->phi();
    double eta = track->eta();

    double xp[3] = { 0, 0, 0 };
 
    if ( track->genParticle()->production_vertex() ) { 
      xp[0] = track->genParticle()->production_vertex()->position().x();
      xp[1] = track->genParticle()->production_vertex()->position().y();
      xp[2] = track->genParticle()->production_vertex()->position().z();
    }
 
    // CHANGED BY JK - z0 with respect to (0,0)
    //    double z0 = xp[2];
    double theta = 2*std::atan( exp( (-1)*eta ) );
    double z0 = xp[2] - (xp[0]*std::cos(phi) + xp[1]*std::sin(phi))/std::tan(theta);
    
    double xd[3] = { 0, 0, 0 };
 
    if ( track->genParticle()->end_vertex() ) { 
      xd[0] = track->genParticle()->end_vertex()->position().x();
      xd[1] = track->genParticle()->end_vertex()->position().y();
      xd[2] = track->genParticle()->end_vertex()->position().z();
    }
 
    double rp = std::sqrt( xp[0]*xp[0] + xp[1]*xp[1] ); 
    double rd = std::sqrt( xd[0]*xd[0] + xd[1]*xd[1] ); 
 
    
    bool final_state = false; 

    const double inner_radius = m_radius;
    const double outer_radius = m_radius;
    if ( ( track->genParticle()->production_vertex() && rp<=inner_radius ) && 
     ( track->genParticle()->end_vertex()==0 || rd>outer_radius ) ) final_state = true; 
     
    
    if ( !final_state ) return 0; 
    
 
    double q = track->charge();
    
    static const particleType ptype;

    if ( q==-999 ) q = ptype.charge( track->pdgId() );

    if ( q==0 ) return 0;
    
    double pT  = q*track->pt();
 
 
    double d0  = 0;


 
 
    // CHANGED BY JK - d0 with respect to (0,0)
    // d0 = q*rp*std::sin(phi);
    d0  =  xp[1]*std::cos(phi) -  xp[0]*std::sin(phi);
 

 
    double dz0 = 0;
    double dd0 = 0;
   
    correctToBeamline( z0, dz0, d0, dd0, theta, phi );
 
 


    int author  = track->pdgId();   
    int uniqueID = HepMC::uniqueID(track); 
 
 
    unsigned long id = (unsigned long)track;
    if ( tid!=0 ) id = tid;

 

        
    //    std::cout << "\t\t\tSUTT Truth track" 
    //              << "\teta=" << eta  // << " +- " << (*trackitr)->params()->deta()
    //              << "\tphi=" << phi  // << " +- " << (*trackitr)->params()->dphi()
    //              << "\tz0="  << z0 
    //              << "\tpT="  << pT // << "\t( " << 1/qoverp << ")"
    //              << "\td0="  << d0
    //              << "\tauthor=" << author
    //              << std::endl;
    
 
 
    TIDA::Track* t = new TIDA::Track(eta, phi, z0, d0, pT, 0, 0,
                                     0, 0, 0, 0, 0,
                                     0, 0, 0, 0,
                                     0, 0, 0, 0,
                                     author, false, uniqueID, -1,
                                     false, 
                                     id ) ;  
 
    return t;
 
}
 
 
 
// add a Trk::Track
bool TrigTrackSelector::selectTrack( const Trk::Track* track ) { 
        
    // do the track extraction stuff here....
 
    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 };
    //    std::cout << __FILE__<< " " <<__LINE__<<std::endl;
 
    if ( track ) { 
        
      // const Trk::Perigee* startPerigee = track->perigeeParameters();
 
#ifdef TRKPARAMETERS_MEASUREDPERIGEE_H
      const Trk::MeasuredPerigee* startPerigee = dynamic_cast<const Trk::MeasuredPerigee*>(track->perigeeParameters());
      //      const Trk::MeasuredPerigee* measPer = startPerigee; // just out of laziness
#else
      const Trk::Perigee* startPerigee = track->perigeeParameters();
      const Trk::Perigee* measPer = startPerigee; // just out of laziness
#endif
 
 
      //      CLHEP::HepVector perigeeParams = measPer->parameters();      
      //      double pT    = measPer->pT(); 
      //      double eta   = measPer->eta();
      //      double phi   = perigeeParams[Trk::phi0];
      //      double z0    = perigeeParams[Trk::z0];
      //      double d0    = perigeeParams[Trk::d0];
      //      // AAARCH!!!!! the TrackParticle pT is NOT SIGNED!!!! ( I ask you! ) 
      //      if ( perigeeParams[Trk::qOverP]<0 ) pT *= -1;
      //      std::cout <<pT1<<" pt1vspT "<<pT<<std::endl;
      
      if (startPerigee){
     
    double theta   = startPerigee->parameters()[Trk::theta];  
    double p       = 1/startPerigee->parameters()[Trk::qOverP];
    double qOverPt = startPerigee->parameters()[Trk::qOverP]/std::sin(theta);
    double charge  = startPerigee->charge();
    double eta     = startPerigee->eta();
    double phi     = startPerigee->parameters()[Trk::phi0];
    double z0      = startPerigee->parameters()[Trk::z0];
    double d0      = startPerigee->parameters()[Trk::d0];
    //  double pT = (1./qOverPt)*(charge);   
    double pT      = (1./qOverPt); // always use signed PT
    
    if ( charge<0 && pT>0 ) pT *= -1;
    if ( charge<0 && p>0 )  p  *= -1;
    
 
 
#ifdef TRKPARAMETERS_MEASUREDPERIGEE_H
    const Trk::ErrorMatrix err = startPerigee->localErrorMatrix();
    double dtheta = err.error(Trk::theta);
    double dqovp  = err.error(Trk::qOverP); 
    double covthetaOvP  = err.covValue(Trk::qOverP,Trk::theta);  
    
    double dphi = err.error(Trk::phi0);
    double dz0  = err.error(Trk::z0);
    double dd0  = err.error(Trk::d0);
#else
    double dtheta      = std::sqrt((*measPer->covariance())(Trk::theta,Trk::theta));
    double dqovp       = std::sqrt((*measPer->covariance())(Trk::qOverP,Trk::qOverP)); 
    double covthetaOvP = (*measPer->covariance())(Trk::qOverP,Trk::theta);    
 
    double dphi = std::sqrt((*measPer->covariance())(Trk::phi0,Trk::phi0));
    double dz0  = std::sqrt((*measPer->covariance())(Trk::z0,Trk::z0));
    double dd0  = std::sqrt((*measPer->covariance())(Trk::d0,Trk::d0));
#endif  
    
    double deta = 0.5*dtheta/(std::cos(0.5*theta)*std::cos(0.5*theta)*std::tan(0.5*theta));  
 
 
    if ( m_correctTrkTracks ) correctToBeamline( z0, dz0, d0, dd0, theta, phi );
 
    double dpT  = 0; 
    
    
    double sintheta = std::sin(theta);
    double costheta = std::cos(theta);
    double dpT2 = (p*p*sintheta)*(p*p*sintheta)*dqovp*dqovp + (p*costheta)*(p*costheta)*dtheta*dtheta - 2*(p*p*sintheta)*(p*costheta)*covthetaOvP;
    
    if ( dpT2>0 ) dpT = std::sqrt( dpT2 );
 
    // Check number of hits
    // NB: a spacepoint is two offline "hits", so a pixel spacepoint is really 
    // 2 "hits" and an offline SCT "hit" is really a 1D cluster, so two intersetcting
    // stereo clusters making a spacepoint are two "hits"
    // const Trk::TrackSummary *summary = dynamic_cast<const Trk::TrackSummary*>(track->trackSummary());
    //ToolHandle< Trk::ITrackSummaryTool > m_trackSumTool;
    //m_trackSumTool  = ToolHandle<Trk::ITrackSummaryTool>("Trk::TrackSummaryTool/InDetTrackSummaryTool");
    //const Trk::TrackSummary* summary = NULL;       
    //summary = m_trackSumTool->createSummary(*track); 
 
    const Trk::TrackSummary * summary = track->trackSummary();
    int nBlayerHits = 0;
    int nPixelHits  = 0;
    int nSctHits    = 0;
    int nStrawHits  = 0;
    int nTrHits     = 0;
    int nSiHits     = 0;
    bool expectBL = false; // Not stored for Trk::Track
    unsigned bitmap = 0;
 
    if(summary==0){
            std::cerr << "Could not create TrackSummary  - Track will likely fail hits requirements" << std::endl;
    }    
    else{      
            nBlayerHits = 2*summary->get(Trk::numberOfBLayerHits); 
            nPixelHits  = 2*summary->get(Trk::numberOfPixelHits);  
        nSctHits    = summary->get(Trk::numberOfSCTHits); 
            nStrawHits  = summary->get(Trk::numberOfTRTHits);
            nTrHits     = summary->get(Trk::numberOfTRTHighThresholdHits);
        nSiHits     = nPixelHits + nSctHits;
 
        for ( int ih=0 ; ih<20 ; ih++ ) {
          if ( summary->isHit(Trk::DetectorType(ih)) ) bitmap |= ( 1<<hpmap[ih] );  
        }
    }
 
    unsigned long id = (unsigned long)track;
    double chi2 = 0;
    double dof  = 0; 
    //const Trk::FitQuality *quality   = dynamic_cast<const Trk::FitQuality*>(track->fitQuality());
    const Trk::FitQuality *quality   = (track->fitQuality());
    if(quality==0) std::cerr << "Could not create FitQuality  - Track will likely fail hits requirements" << std::endl;
    else{
            chi2 = quality->chiSquared();          
        dof  = quality->numberDoF();
    }
 
    int trackAuthor = -1;

 
    //      std::cout << "fetching author info :" << track->info().trackFitter() << ":" 
    //      << track->info().dumpInfo() << ": bm 0x" << std::hex << bitmap << std::dec << std::endl;
 
    int fitter = track->info().trackFitter();
    //      std::string dumpinfo = track->info().dumpInfo();
 
    if ( fitter>0 && fitter<Trk::TrackInfo::NumberOfTrackFitters ) {
      if     ((track->info().dumpInfo()).find("TRTStandalone")        != std::string::npos) trackAuthor = 2;
      else if   ((track->info().dumpInfo()).find("TRTSeededTrackFinder") != std::string::npos) trackAuthor = 1;
      else                                                                                    trackAuthor = 0;
    }
 
  #if 0
    std::cout << "\t\t\tSUTT TP track" 
          << "\teta=" << eta  // << " +- " << (*trackitr)->params()->deta()
          << "\tphi=" << phi  // << " +- " << (*trackitr)->params()->dphi()
          << "\tz0="  << z0 
          << "\tpT="  << pT // << "\t( " << 1/qoverp << ")"
          << "\td0="  << d0
          << "\tNsi=" << nSiHits 
          << "\tNtrt=" << nTrHits 
          << "\tNstr=" << nStrawHits
          << "\tauthor=" << trackAuthor
          << std::endl;
  #endif    
    // Create and save Track      
    TIDA::Track* t = new TIDA::Track(eta, phi, z0, d0, pT, chi2, dof,
                                         deta, dphi, dz0, dd0, dpT,
                                         nBlayerHits, nPixelHits, nSctHits, nSiHits,
                                         nStrawHits, nTrHits, bitmap, 0,
                                         trackAuthor,  false, -1, -1,  
                                         expectBL, id) ;  
 
    if ( !addTrack( t ) ){
      delete t;
      return false;
    }
    return true;
    
    //std::cout << "SUTT TP track " << *t << "\t0x" << std::hex << bitmap << std::dec << std::endl; 
      }
    }
    
    return false;
}
 
// extract all the tracks from a TrackCollection and add them
void  TrigTrackSelector::selectTracks( const TrackCollection* trigtracks ) { 
    
    //    std::cout << "\t\t\tSUTT \tTrackContainer->size() = " << trigtracks->size() << std::endl;
    
    TrackCollection::const_iterator  trackitr = trigtracks->begin();
    TrackCollection::const_iterator  trackend = trigtracks->end();
   
    while ( trackitr!=trackend ) { 
      selectTrack( *trackitr );
      ++trackitr;
    } // loop over tracks
    
}
 
 
 
bool TrigTrackSelector::selectTrack( const xAOD::TrackParticle* track, void* ) {
    // do the track extraction stuff here....
 
    //    static int hpmap[20] = { 0, 1, 2,  7, 8, 9,  3, 4, 5, 6, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 };
 
    if ( track ) { 
        
      const xAOD::TrackParticle* measPer = track;
 
      //     CLHEP::HepVector perigeeParams = measPer->parameters();
      
      double pT    = measPer->pt(); 
      double eta   = measPer->eta();
      double phi   = measPer->phi0();
      double z0    = measPer->z0() + measPer->vz();  
      double d0    = measPer->d0();
 
      double theta = measPer->theta();

      // if ( measPer->qOverP()==0 ) return false;
      if ( measPer->qOverP()==0 ) throw std::runtime_error( "probable corrupted track - this should never happen" );
      double  p    = 1/measPer->qOverP();
  
      // AAARCH!!!!! the TrackParticle pT is NOT SIGNED!!!! ( I ask you! ) 
      if ( measPer->qOverP()<0 && pT>0 ) pT *= -1;
 
      double dtheta      = std::sqrt(measPer->definingParametersCovMatrix()(Trk::theta,Trk::theta));
      double dqovp       = std::sqrt(measPer->definingParametersCovMatrix()(Trk::qOverP,Trk::qOverP)); 
      double covthetaOvP = measPer->definingParametersCovMatrix()(Trk::qOverP,Trk::theta);  
 
      double deta = 0.5*dtheta/(std::cos(0.5*theta)*std::cos(0.5*theta)*std::tan(0.5*theta));  // ???? CHECK THIS <<--
 
      double dphi = std::sqrt(measPer->definingParametersCovMatrix()(Trk::phi0,Trk::phi0));
      double dz0  = std::sqrt(measPer->definingParametersCovMatrix()(Trk::z0,Trk::z0));
      double dd0  = std::sqrt(measPer->definingParametersCovMatrix()(Trk::d0,Trk::d0));
 
      double dpT  = 0; 
 

      //      if ( m_xBeam!=0 || m_yBeam!=0 ) correctToBeamline( z0, dz0, d0, dd0, theta, phi );
      
 
      double sintheta = std::sin(theta);
      double costheta = std::cos(theta);
      double dpt2 = (p*p*sintheta)*(p*p*sintheta)*dqovp*dqovp + (p*costheta)*(p*costheta)*dtheta*dtheta - 2*(p*p*sintheta)*(p*costheta)*covthetaOvP;   
 
      if ( dpt2>0 ) dpT = std::sqrt( dpt2 );
 
      // Check number of hits
      // NB: a spacepoint is two offline "hits", so a pixel spacepoint is really 
      // 2 "hits" and an offline SCT "hit" is really a 1D cluster, so two intersetcting
      // stereo clusters making a spacepoint are two "hits"
      
      uint8_t sum_nBlayerHits = 0;
      track->summaryValue( sum_nBlayerHits, xAOD::numberOfInnermostPixelLayerHits);
      int nBlayerHits = 2*sum_nBlayerHits;
 
      uint8_t  sum_nPixelHits = 0;
      track->summaryValue( sum_nPixelHits, xAOD::numberOfPixelHits);
      int nPixelHits = 2*sum_nPixelHits;
 
      uint8_t  sum_nSctHits = 0;
      track->summaryValue( sum_nSctHits, xAOD::numberOfSCTHits);
      int nSctHits = sum_nSctHits;
 
      uint8_t  sum_nStrawHits = 0;
      track->summaryValue( sum_nStrawHits, xAOD::numberOfTRTHits);
      int nStrawHits  = sum_nStrawHits;
 
      uint8_t  sum_nTrtHits = 0;
      track->summaryValue( sum_nTrtHits, xAOD::numberOfTRTHighThresholdHits);
      int nTrtHits  = sum_nTrtHits;
      
 
      uint8_t sum_expectBL  = 0;
      track->summaryValue( sum_expectBL, xAOD::expectInnermostPixelLayerHit);
      bool expectBL = ( sum_expectBL ? true : false );

 
      uint8_t sum_sctholes  = 0;
      track->summaryValue( sum_sctholes, xAOD::numberOfSCTHoles);
 
      uint8_t sum_pixholes  = 0;
      track->summaryValue( sum_pixholes, xAOD::numberOfPixelHoles);

 
      nSctHits   += 1000*sum_sctholes;
      nPixelHits += 1000*sum_pixholes;

      int nSiHits     = nPixelHits + nSctHits;

 
      double chi2 = track->chiSquared();
      double dof  = track->numberDoF();
      
      unsigned long id = (unsigned long)track;
 
      unsigned bitmap = track->hitPattern();
 
 
 
      double xbeam = track->vx(); 
      double ybeam = track->vy(); 
      double zbeam = track->vz(); 
      
      if ( xbeam!=getBeamX() || ybeam!=getBeamY() || zbeam!=getBeamZ() ) setBeamline( xbeam, ybeam, zbeam );
    
      int trackAuthor = 0;
 
      int fitter = track->trackFitter();
      std::bitset<xAOD::NumberOfTrackRecoInfo>  patternrec = track->patternRecoInfo();
 
      //int icount = 0;<- never used if section below is commented
      for ( unsigned ipr=patternrec.size() ; ipr-- ; ) { 
    if ( patternrec[ipr] ) {
      //icount++; <- never used if section below is commented
      trackAuthor |= (ipr >> 16);
      // static bool first = true;
      // if ( first && icount>1 ) { 
      //   std::cerr << "more than one pattern rec strategy " << ipr << "\t(suppressing further output)" << std::endl;
      //   first = false;
      // }
    }
      }
      
      trackAuthor |= fitter;
 
      //   if ( fitter>0 && fitter<Trk::TrackInfo::NumberOfTrackFitters ) {
      //    if      ( dumpinfo.find("TRTStandalone")!=std::string::npos)        trackAuthor = 2;
      //    else if ( dumpinfo.find("TRTSeededTrackFinder")!=std::string::npos) trackAuthor = 1;
      //    else                                                                trackAuthor = 0;
      //   }
      
#if 0
      std::cout << "\t\t\tSUTT TP track" 
        << "\teta=" << eta  // << " +- " << (*trackitr)->params()->deta()
        << "\tphi=" << phi  // << " +- " << (*trackitr)->params()->dphi()
        << "\tz0="  << z0 
        << "\tpT="  << pT // << "\t( " << 1/qoverp << ")"
        << "\td0="  << d0
        << "\tNsi=" << nSiHits 
    //      << "\tNtrt=" << nTrtHits 
    //      << "\tNstr=" << nStrawHits
        << "\tfitter=" << fitter
        << "\tauthor=" << trackAuthor
        << "\tVTX x " << track->vx() << "\ty " << track->vy() << "\tz " << track->vz() 
        << std::endl;
#endif  
 
      // Create and save Track
      
      TIDA::Track* t = new TIDA::Track( eta,   phi,  z0,  d0,  pT, chi2, dof,
                        deta,  dphi, dz0, dd0, dpT,
                        nBlayerHits, nPixelHits, nSctHits, nSiHits,
                        nStrawHits,  nTrtHits,   bitmap, 0,
                        trackAuthor,  false, -1, -1,  
                        expectBL, id) ;  
 
      //      std::cout << "SUTT TP track " << *t << "\t0x" << std::hex << bitmap << std::dec << std::endl; 
      
      if ( !addTrack( t ) ){
    delete t;
    return false;
      }
      return true;
      
    }
    return false;
     
 
}
 
 
void TrigTrackSelector::selectTracks( const xAOD::TrackParticleContainer* tracks, void* ) {
    //    std::cout << "\t\t\tSUTT \tTrackContainer->size() = " << trigtracks->size() << std::endl;
    xAOD::TrackParticleContainer::const_iterator  trackitr = tracks->begin();
    xAOD::TrackParticleContainer::const_iterator  trackend = tracks->end();
    while ( trackitr!=trackend ) { 
      selectTrack( *trackitr );
      ++trackitr;
    } // loop over tracks     
}
 
 
 
void TrigTrackSelector::selectTracks( xAOD::TrackParticleContainer::const_iterator trackitr, 
                      xAOD::TrackParticleContainer::const_iterator trackend, void* ) {
    //    std::cout << "\t\t\tSUTT \tTrackContainer->size() = " << trigtracks->size() << std::endl;
    while ( trackitr!=trackend ) { 
      selectTrack( *trackitr );
      ++trackitr;
    } // loop over tracks     
}
 
 
void TrigTrackSelector::selectTracks( const TrackParticleLinks_t& tracks ) {
  for( const auto& track : tracks ) selectTrack( *track );
}
 
 
 
 
 
void TrigTrackSelector::correctToBeamline( double& z0,    double& dz0, 
                       double& d0,    double& dd0, 
                       double  theta, double  phi ) {
    
    
    //   if ( m_first ) { 
    //     if ( m_xBeam==0 && m_yBeam==0 ) { 
    //      std::cerr << "TrigTrackSelector::correctToBeamline() WARNING -- Beamline set to (0,0) -- WARNING" << std::endl;  
    //     }
    //     else { 
    //     std::cout << "TrigTrackSelector::correctToBeamline() Beamline set to " << m_xBeam << " " << m_yBeam << std::endl;  
    //     }
    //     m_first = false;
    //   }
    
    
    //      double theta = 2*std::atan( exp( (-1)*eta ) );
    double z0t = z0 + ((std::cos(phi)*m_xBeam + std::sin(phi)*m_yBeam)/std::tan(theta));    
    double a0t = d0 +   std::sin(phi)*m_xBeam - std::cos(phi)*m_yBeam; 
    
    double dz0t = dz0 + ((std::cos(phi)*m_xBeam + std::sin(phi)*m_yBeam)/std::tan(theta));
    double da0t = dd0 +   std::sin(phi)*m_xBeam - std::cos(phi)*m_yBeam;
    
    z0  = z0t;
    d0  = a0t;
    
    dz0 = dz0t;
    dd0 = da0t;
}