FWPFTrackUtils.cc

Go to the documentation of this file.

#include "Fireworks/ParticleFlow/interface/FWPFTrackUtils.h"

FWPFTrackSingleton *FWPFTrackSingleton::pInstance = NULL;
bool FWPFTrackSingleton::instanceFlag = false;

//______________________________________________________________________________
FWPFTrackSingleton *
FWPFTrackSingleton::Instance()
{
   if( !instanceFlag )  // Instance doesn't exist yet
   {
      pInstance = new FWPFTrackSingleton();
      instanceFlag = true;
   }

   return pInstance;    // Pointer to sole instance
}

//______________________________________________________________________________
void
FWPFTrackSingleton::initPropagator()
{
   m_magField = new FWMagField();

   // Common propagator, helix stepper
   m_trackPropagator = new TEveTrackPropagator();
   m_trackPropagator->SetMagFieldObj( m_magField, false );
   m_trackPropagator->SetMaxR( FWPFGeom::caloR3() );
   m_trackPropagator->SetMaxZ( FWPFGeom::caloZ2() );
   m_trackPropagator->SetDelta( 0.01 );
   m_trackPropagator->SetProjTrackBreaking( TEveTrackPropagator::kPTB_UseLastPointPos );
   m_trackPropagator->SetRnrPTBMarkers( kTRUE );
   m_trackPropagator->IncDenyDestroy();

   // Tracker propagator
   m_trackerTrackPropagator = new TEveTrackPropagator();
   m_trackerTrackPropagator->SetStepper( TEveTrackPropagator::kRungeKutta );
   m_trackerTrackPropagator->SetMagFieldObj( m_magField, false );
   m_trackerTrackPropagator->SetDelta( 0.01 );
   m_trackerTrackPropagator->SetMaxR( FWPFGeom::caloR3() );
   m_trackerTrackPropagator->SetMaxZ( FWPFGeom::caloZ2() );
   m_trackerTrackPropagator->SetProjTrackBreaking( TEveTrackPropagator::kPTB_UseLastPointPos );
   m_trackerTrackPropagator->SetRnrPTBMarkers( kTRUE );
   m_trackerTrackPropagator->IncDenyDestroy();
}

//______________________________________________________________________________
FWPFTrackUtils::FWPFTrackUtils()
{
   m_singleton = FWPFTrackSingleton::Instance();
}


//______________________________________________________________________________
TEveTrack *
FWPFTrackUtils::getTrack( const reco::Track &iData )
{
   TEveTrackPropagator *propagator = ( !iData.extra().isAvailable() ) ? 
                           m_singleton->getTrackerTrackPropagator() : m_singleton->getTrackPropagator();

   TEveRecTrack t;
   t.fBeta = 1;
   t.fP = TEveVector( iData.px(), iData.py(), iData.pz() );
   t.fV = TEveVector( iData.vertex().x(), iData.vertex().y(), iData.vertex().z() );
   t.fSign = iData.charge();
   TEveTrack *trk = new TEveTrack( &t, propagator );
   trk->MakeTrack();

   return trk;
}

//______________________________________________________________________________
TEveStraightLineSet *
FWPFTrackUtils::setupLegoTrack( const reco::Track &iData )
{
   using namespace FWPFGeom;

   // Declarations
   int wraps[3] = { -1, -1, -1 };
   bool ECAL = false;
   TEveTrack *trk = getTrack( iData );
   std::vector<TEveVector> trackPoints( trk->GetN() - 1 );
   const Float_t *points = trk->GetP();
   TEveStraightLineSet *legoTrack = new TEveStraightLineSet();

   if( m_singleton->getField()->getSource() == FWMagField::kNone )
   {
      if( fabs( iData.eta() ) < 2.0 && iData.pt() > 0.5 && iData.pt() < 30 )
      {
         double estimate = fw::estimate_field( iData, true );
         if( estimate >= 0 ) m_singleton->getField()->guessField( estimate );
      }
   } 

   // Convert to Eta/Phi and store in vector
   for( Int_t i = 1; i < trk->GetN(); ++i )
   {
      int j = i * 3;
      TEveVector temp = TEveVector( points[j], points[j+1], points[j+2] );
      TEveVector vec = TEveVector( temp.Eta(), temp.Phi(), 0.001 );

      trackPoints[i-1] = vec;
   }

   // Add first point to ps if necessary
   for( Int_t i = 1; i < trk->GetN(); ++i )
   {
      int j = i * 3;
      TEveVector v1 = TEveVector( points[j], points[j+1], points[j+2] );

      if( !ECAL )
      {
         if( FWPFMaths::checkIntersect( v1, caloR1() ) )
         {        
            TEveVector v2 = TEveVector( points[j-3], points[j-2], points[j-1] );
            TEveVector xyPoint = FWPFMaths::lineCircleIntersect( v1, v2, caloR1() );
            TEveVector zPoint;
            if( v1.fZ < 0 )
               zPoint = TEveVector( xyPoint.fX, xyPoint.fY, v1.fZ - 50.f );
            else
               zPoint = TEveVector( xyPoint.fX, xyPoint.fY, v1.fZ + 50.f );

            TEveVector vec = FWPFMaths::lineLineIntersect( v1, v2, xyPoint, zPoint );
            legoTrack->AddMarker( vec.Eta(), vec.Phi(), 0.001, 0 );

            wraps[0] = i;        // There is now a chance that the track will also reach the HCAL radius
            ECAL = true;
         }
         else if( fabs( v1.fZ ) >= caloZ1() )
         {
            TEveVector p1, p2;
            TEveVector vec, v2 = TEveVector( points[j-3], points[j-2], points[j-1] );
            float z, y = FWPFMaths::linearInterpolation( v2, v1, caloZ1() );

            if( v2.fZ > 0 )
               z = caloZ1();
            else
               z = caloZ1() * -1;

            p1 = TEveVector( v2.fX + 50, y, z );
            p2 = TEveVector( v2.fX - 50, y, z );
            vec = FWPFMaths::lineLineIntersect( v1, v2, p1, p2 );

            legoTrack->AddMarker( vec.Eta(), vec.Phi(), 0.001, 0 );
            wraps[0] = i;
            ECAL = true;
         }
      }
      else if( FWPFMaths::checkIntersect( v1, caloR2() ) )
      {
         TEveVector v2 = TEveVector( points[j-3], points[j-2], points[j-1] );
         TEveVector xyPoint = FWPFMaths::lineCircleIntersect( v1, v2, caloR2() );
         TEveVector zPoint;
         if( v1.fZ < 0 )
            zPoint = TEveVector( xyPoint.fX, xyPoint.fY, v1.fZ - 50.f );
         else
            zPoint = TEveVector( xyPoint.fX, xyPoint.fY, v1.fZ + 50.f );

         TEveVector vec = FWPFMaths::lineLineIntersect( v1, v2, xyPoint, zPoint );
         legoTrack->AddMarker( vec.Eta(), vec.Phi(), 0.001, 1 );

         wraps[1] = i;        // There is now a chance that the track will also reach the HCAL radius
         break;
      }
   }

   if( wraps[0] != -1 ) //if( ECAL )
   {
      int i = ( trk->GetN() - 1 ) * 3;
      int j = trk->GetN() - 2;   // This is equal to the last index in trackPoints vector
      TEveVector vec = TEveVector( points[i], points[i+1], points[i+2] );

      if( FWPFMaths::checkIntersect( vec, caloR3() - 1 ) )
      {
         legoTrack->AddMarker( vec.Eta(), vec.Phi(), 0.001, 2 );
         wraps[2] = j;
      }
      else if( fabs( vec.fZ ) >= caloZ2() )
      {
         legoTrack->AddMarker( vec.Eta(), vec.Phi(), 0.001, 2 );
         wraps[2] = j;
      }
   }

   /* Handle phi wrapping */
   for( int i = 0; i < static_cast<int>( trackPoints.size() - 1 ); ++i )
   {
      if( ( trackPoints[i+1].fY - trackPoints[i].fY ) > 1 )
      {
         trackPoints[i+1].fY -= TMath::TwoPi();
         if( i == wraps[0] )
         {
            TEveChunkManager::iterator mi( legoTrack->GetMarkerPlex() );
            mi.next();  // First point
            TEveStraightLineSet::Marker_t &m = * ( TEveStraightLineSet::Marker_t* ) mi();
            m.fV[0] = trackPoints[i+1].fX; m.fV[1] = trackPoints[i+1].fY; m.fV[2] = 0.001;
         }
         else if( i == wraps[1] )
         {
            TEveChunkManager::iterator mi( legoTrack->GetMarkerPlex() );
            mi.next(); mi.next();  // Second point
            TEveStraightLineSet::Marker_t &m = * ( TEveStraightLineSet::Marker_t* ) mi();
            m.fV[0] = trackPoints[i+1].fX; m.fV[1] = trackPoints[i+1].fY; m.fV[2] = 0.001;
         }
      }

      if( ( trackPoints[i].fY - trackPoints[i+1].fY ) > 1 )
      {
         trackPoints[i+1].fY += TMath::TwoPi();
         if( i == wraps[0] )
         {
            TEveChunkManager::iterator mi( legoTrack->GetMarkerPlex() );
            mi.next();  // First point
            TEveStraightLineSet::Marker_t &m = * ( TEveStraightLineSet::Marker_t* ) mi();
            m.fV[0] = trackPoints[i+1].fX; m.fV[1] = trackPoints[i+1].fY; m.fV[2] = 0.001;
         }
         else if( i == wraps[1] )
         {
            TEveChunkManager::iterator mi( legoTrack->GetMarkerPlex() );
            mi.next();  mi.next(); // Second point
            TEveStraightLineSet::Marker_t &m = * ( TEveStraightLineSet::Marker_t* ) mi();
            m.fV[0] = trackPoints[i+1].fX; m.fV[1] = trackPoints[i+1].fY; m.fV[2] = 0.001;
         }
      }
   }

   int end = static_cast<int>( trackPoints.size() - 1 );
   if( wraps[2] == end )
   {
      TEveChunkManager::iterator mi( legoTrack->GetMarkerPlex() );
      mi.next(); mi.next(); mi.next();   // Third point
      TEveStraightLineSet::Marker_t &m = * ( TEveStraightLineSet::Marker_t* ) mi();
      m.fV[0] = trackPoints[end].fX; m.fV[1] = trackPoints[end].fY; m.fV[2] = 0.001;
   }

   // Set points on TEveLineSet object ready for displaying
   for( unsigned int i = 0;i < trackPoints.size() - 1; ++i )
      legoTrack->AddLine( trackPoints[i], trackPoints[i+1] );

   legoTrack->SetDepthTest( false );
   legoTrack->SetMarkerStyle( 4 );
   legoTrack->SetMarkerSize( 1 );
   legoTrack->SetRnrMarkers( true );

   delete trk;    // Release memory that is no longer required

   return legoTrack;
}

//______________________________________________________________________________
TEveTrack *
FWPFTrackUtils::setupTrack( const reco::Track &iData )
{
    if( m_singleton->getField()->getSource() == FWMagField::kNone )
   {
      if( fabs( iData.eta() ) < 2.0 && iData.pt() > 0.5 && iData.pt() < 30 )
      {
         double estimate = fw::estimate_field( iData, true );
         if( estimate >= 0 ) m_singleton->getField()->guessField( estimate );
      }
   } 

   TEveTrack *trk = getTrack( iData );
   
   return trk;
}


//______________________________________________________________________________
TEvePointSet *
FWPFTrackUtils::getCollisionMarkers( const TEveTrack *trk )
{
   using namespace FWPFGeom;

   bool ECAL = false;
   const Float_t *points = trk->GetP();
   TEvePointSet *ps = new TEvePointSet();

   for( Int_t i = 1; i < trk->GetN(); ++i )
   {
      int j = i * 3;
      TEveVector v1 = TEveVector( points[j], points[j+1], points[j+2] );

      if( !ECAL )
      {
         if( FWPFMaths::checkIntersect( v1, caloR1() ) )
         {
            TEveVector v2 = TEveVector( points[j-3], points[j-2], points[j-1] );
            TEveVector xyPoint = FWPFMaths::lineCircleIntersect( v1, v2, caloR1() );
            TEveVector zPoint;
            if( v1.fZ < 0 )
               zPoint = TEveVector( xyPoint.fX, xyPoint.fY, v1.fZ - 50.f );
            else
               zPoint = TEveVector( xyPoint.fX, xyPoint.fY, v1.fZ + 50.f );

            TEveVector vec = FWPFMaths::lineLineIntersect( v1, v2, xyPoint, zPoint );
            ps->SetNextPoint( vec.fX, vec.fY, vec.fZ );

            ECAL = true;
         }
         else if( fabs( v1.fZ ) >= caloZ1() )
         {
            TEveVector p1, p2;
            TEveVector vec, v2 = TEveVector( points[j-3], points[j-2], points[j-1] );
            float z, y = FWPFMaths::linearInterpolation( v2, v1, caloZ1() );

            if( v2.fZ > 0 )
               z = caloZ1();
            else
               z = caloZ1() * -1;

            p1 = TEveVector( v2.fX + 50, y, z );
            p2 = TEveVector( v2.fX - 50, y, z );
            vec = FWPFMaths::lineLineIntersect( v1, v2, p1, p2 );

            ps->SetNextPoint( vec.fX, vec.fY, vec.fZ );
            ECAL = true;
         }
      }
      else if ( FWPFMaths::checkIntersect( v1, caloR2() ) )
      {
         TEveVector v2 = TEveVector( points[j-3], points[j-2], points[j-1] );
         TEveVector xyPoint = FWPFMaths::lineCircleIntersect( v1, v2,  caloR2() );
         TEveVector zPoint;
         if( v1.fZ < 0 )
            zPoint = TEveVector( xyPoint.fX, xyPoint.fY, v1.fZ - 50.f );
         else
            zPoint = TEveVector( xyPoint.fX, xyPoint.fY, v1.fZ + 50.f );

         TEveVector vec = FWPFMaths::lineLineIntersect( v1, v2, xyPoint, zPoint );
         ps->SetNextPoint( vec.fX, vec.fY, vec.fZ );
         break;   // ECAL and HCAL collisions found so stop looping
      }
   }

   // HCAL collisions (outer radius and endcap)
   int i = ( trk->GetN() - 1 ) * 3;
   TEveVector vec = TEveVector( points[i], points[i+1], points[i+2] );

   if( FWPFMaths::checkIntersect( vec, caloR3() - 1 ) )
      ps->SetNextPoint( vec.fX, vec.fY, vec.fZ );
   else if( fabs( vec.fZ ) >= caloZ2() )
      ps->SetNextPoint( vec.fX, vec.fY, vec.fZ );

   return ps;
}