#include <NavPropagator.h>

Inheritance diagram for NavPropagator:

Public Member Functions
virtual NavPropagator *	clone () const
virtual const MagneticField *	magneticField () const
	NavPropagator (const MagneticField *field, PropagationDirection dir=alongMomentum)
TrajectoryStateOnSurface	propagate (const FreeTrajectoryState &ts, const Plane &plane) const
	propagation of FTS to plane
TrajectoryStateOnSurface	propagate (const FreeTrajectoryState &fts, const Cylinder &cylinder) const
	propagation to cylinder
TrajectoryStateOnSurface	propagate (const TrajectoryStateOnSurface &ts, const Plane &plane) const
	propagation of TSOS to plane
std::pair < TrajectoryStateOnSurface, double >	propagateWithPath (const TrajectoryStateOnSurface &, const Plane &plane) const
	propagation of TSOS to plane with path length
std::pair < TrajectoryStateOnSurface, double >	propagateWithPath (const FreeTrajectoryState &fts, const Cylinder &cylinder) const
	propagation to cylinder with path length
std::pair < TrajectoryStateOnSurface, double >	propagateWithPath (const FreeTrajectoryState &, const Plane &plane) const
	propagation of FTS to plane with path length
	~NavPropagator ()
Private Types
typedef std::map< const MagVolume , NavVolume >	MagVolumeMap
typedef RKPropagatorInS	PropagatorType
typedef TrajectoryStateOnSurface	TSOS
typedef std::pair < TrajectoryStateOnSurface, double >	TsosWP
Private Member Functions
bool	destinationCrossed (const TSOS &startState, const TSOS &endState, const Plane &plane) const
const NavVolume *	findVolume (const TrajectoryStateOnSurface &inputState) const
const NavVolume *	navVolume (const MagVolume *magVolume) const
TrajectoryStateOnSurface	noNextVolume (TrajectoryStateOnSurface startingState) const
std::pair < TrajectoryStateOnSurface, double >	propagateInVolume (const NavVolume *currentVolume, const TrajectoryStateOnSurface &startingState, const Plane &targetPlane) const
Private Attributes
bool	isIronVolume [272]
const VolumeBasedMagneticField *	theField
MagVolumeMap	theNavVolumeMap

Detailed Description

Definition at line 16 of file NavPropagator.h.

Member Typedef Documentation

typedef std::map<const MagVolume*, NavVolume*> NavPropagator::MagVolumeMap [private]

Definition at line 72 of file NavPropagator.h.

typedef RKPropagatorInS NavPropagator::PropagatorType [private]

Definition at line 71 of file NavPropagator.h.

typedef TrajectoryStateOnSurface NavPropagator::TSOS [private]

Definition at line 73 of file NavPropagator.h.

typedef std::pair<TrajectoryStateOnSurface,double> NavPropagator::TsosWP [private]

Definition at line 70 of file NavPropagator.h.

Constructor & Destructor Documentation

NavPropagator::NavPropagator	(	const MagneticField *	field,
		PropagationDirection	dir = `alongMomentum`
	)

Definition at line 14 of file NavPropagator.cc.

References i, isIronVolume, and theField.

Referenced by clone().

                                                        :
  Propagator(dir)
{
  theField = dynamic_cast<const VolumeBasedMagneticField*>(field);
  //FIXME: iron volumes are hard-coded below... will change with new .xml geometry MM 28/6/07
  int myIronVolumes[126] = {6,9,11,14,15,18,22,23,26,29,30,33,36,37,43,46,49,53,56,57,60,62,
                            63,65,71,77,105,106,107,111,112,113,114,115,116,117,118,122,123,
                            125,126,128,129,130,131,132,133,135,137,141,145,147,148,149,153,
                            154,155,156,157,158,159,160,164,165,167,168,170,171,172,173,174,
                            175,177,178,179,180,184,185,186,187,188,189,190,191,195,196,198,
                            200,204,208,210,211,212,216,217,218,219,220,221,222,223,227,228,
                            230,231,233,234,235,236,237,238,240,241,242,243,247,248,249,250,
                            251,252,253,254,258,259,261}; 
  for(int i=0; i<272; i++) isIronVolume[i]=false;
  for(int i=0; i<126 ; i++) isIronVolume[myIronVolumes[i]]=true; 
  //  for(int i=1; i<272 ; i++) {
  //  if(isIronVolume[i]) std::cout << "Volume no." << i << " is made of iron " << std::endl;
  //  if(!isIronVolume[i]) std::cout << "Volume no." << i << " is made of air " << std::endl;
  //}
}

NavPropagator::~NavPropagator ( )

Definition at line 36 of file NavPropagator.cc.

References i, and theNavVolumeMap.

                              {
    for (MagVolumeMap::iterator i = theNavVolumeMap.begin(); i != theNavVolumeMap.end(); ++i) {
      delete i->second;
    }
}

Member Function Documentation

virtual NavPropagator* NavPropagator::clone ( void ) const [inline, virtual]

Implements Propagator.

Definition at line 64 of file NavPropagator.h.

References NavPropagator().

{return new NavPropagator(*this);}

bool NavPropagator::destinationCrossed	(	const TSOS &	startState,
		const TSOS &	endState,
		const Plane &	plane
	)		const `[private]`

Definition at line 313 of file NavPropagator.cc.

References TrajectoryStateOnSurface::globalPosition(), and Plane::side().

Referenced by propagateWithPath().

{
  bool res = ( plane.side( startState.globalPosition(), 1.e-6) != 
               plane.side( endState.globalPosition(), 1.e-6));
  /*
   LogDebug("NavPropagator") <<  "NavPropagator::destinationCrossed called with startState "
       << startState << std::endl
       << " endState " << endState 
       << std::endl
       << " plane at " << plane.position()
       << std::endl;

   LogDebug("NavPropagator") <<  "plane.side( startState) " << plane.side( startState.globalPosition(), 1.e-6);
   LogDebug("NavPropagator") <<  "plane.side( endState)   " << plane.side( endState.globalPosition(), 1.e-6);
  */

  return res;
}

const NavVolume * NavPropagator::findVolume ( const TrajectoryStateOnSurface & inputState ) const [private]

Definition at line 176 of file NavPropagator.cc.

References gather_cfg::cout, VolumeBasedMagneticField::findVolume(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalPosition(), LogDebug, navVolume(), theField, Vector3DBase< T, FrameTag >::unit(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by propagateWithPath().

{
   LogDebug("NavPropagator") <<  "NavPropagator: calling theField->findVolume ";

  GlobalPoint gp = inputState.globalPosition() + 0.1*inputState.globalMomentum().unit();
  // Next protection only needed when 1 mm linear move crosses the z==0 plane:
  GlobalPoint gpSym(gp.x(), gp.y(), (gp.z()<0? gp.z() : -gp.z()));

  const MagVolume* magVolume = theField->findVolume( gpSym);

   LogDebug("NavPropagator") <<  "NavPropagator: got MagVolume* " << magVolume 
       << " when searching with pos " << gpSym ;


   if (!magVolume) {
     cout << "Got invalid volume pointer " << magVolume << " at position " << gpSym;
     return 0;
   } 
  return navVolume(magVolume);
}

const MagneticField * NavPropagator::magneticField ( ) const [virtual]

Implements Propagator.

Definition at line 42 of file NavPropagator.cc.

References theField.

{return theField;}

const NavVolume * NavPropagator::navVolume ( const MagVolume * magVolume ) const [private]

Definition at line 198 of file NavPropagator.cc.

References gather_cfg::cout, i, isIronVolume, n, MagVolume6Faces::name, query::result, and theNavVolumeMap.

Referenced by findVolume().

{
  NavVolume* result;
  MagVolumeMap::iterator i = theNavVolumeMap.find( magVolume);
  if (i == theNavVolumeMap.end()) {
    // Create a NavVolume from a MagVolume if the NavVolume doesn't exist yet
    // FIXME: hardcoded iron/air classification should be made into something more general
    // will break with new magnetic field volume .xml file MM 28/6/07
    const MagVolume6Faces* pVol6 = dynamic_cast<const MagVolume6Faces*> ( magVolume );
    int n=0;
    bool isIron=false;
    if(pVol6) { 
      std::stringstream ss(pVol6->name.substr(5));
      ss >> n;
      //      std::cout << " Found volume with number " << n << std::endl;
    } else {
      std::cout << "Error (NavVolume6Faces) failed to get MagVolume6Faces pointer" << std::endl;
    }
    if(n<1 || n>271) {
      std::cout << "Error (NavVolume6Faces) unexpected Volume number!" << std::endl;
    } else {
      isIron = isIronVolume[n];
    }

    result= new NavVolume6Faces( *magVolume, isIron);
    theNavVolumeMap[magVolume] = result;
  }
  else result = i->second;
  return result;
}

TrajectoryStateOnSurface NavPropagator::noNextVolume ( TrajectoryStateOnSurface startingState ) const [private]

Definition at line 334 of file NavPropagator.cc.

References LogDebug.

Referenced by propagateWithPath().

{
   LogDebug("NavPropagator") <<  std::endl
       << "Propagation reached end of volume geometry without crossing target surface" 
       << std::endl
       << "starting state of last propagation " << startingState;

  return TrajectoryStateOnSurface();
}

TrajectoryStateOnSurface NavPropagator::propagate	(	const FreeTrajectoryState &	fts,
		const Cylinder &	cylinder
	)		const `[inline, virtual]`

propagation to cylinder

Implements Propagator.

Definition at line 55 of file NavPropagator.h.

References propagateWithPath().

                                                                     {
    return propagateWithPath(fts,cylinder).first;
  }

TrajectoryStateOnSurface NavPropagator::propagate	(	const FreeTrajectoryState &	ts,
		const Plane &	plane
	)		const `[inline, virtual]`

propagation of FTS to plane

Implements Propagator.

Definition at line 45 of file NavPropagator.h.

References propagateWithPath().

                                                               {
    return propagateWithPath(ts,plane).first;
  }

TrajectoryStateOnSurface NavPropagator::propagate	(	const TrajectoryStateOnSurface &	ts,
		const Plane &	plane
	)		const `[inline, virtual]`

propagation of TSOS to plane

Reimplemented from Propagator.

Definition at line 35 of file NavPropagator.h.

References propagateWithPath().

                                                               {
    return propagateWithPath(ts,plane).first;
  }

std::pair< TrajectoryStateOnSurface, double > NavPropagator::propagateInVolume	(	const NavVolume *	currentVolume,
		const TrajectoryStateOnSurface &	startingState,
		const Plane &	targetPlane
	)		const `[private]`

Definition at line 230 of file NavPropagator.cc.

References GlobalTrajectoryParameters::charge(), Vector3DBase< T, FrameTag >::cross(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalParameters(), TrajectoryStateOnSurface::globalPosition(), MagVolume::inside(), RKPropagatorInS::magneticField(), Plane::normalVector(), GloballyPositioned< T >::position(), RKPropagatorInS::propagateWithPath(), TrajectoryStateOnSurface::surface(), Vector3DBase< T, FrameTag >::unit(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by propagateWithPath().

{
  PropagatorType prop( *currentVolume);

  // OK, fix is rather ugly for now: need to reflect z>0 states to z<0 MM 25/6/07
  bool isReflected = false;
  TSOS okState = startingState;
  PlaneBuilder::ReturnType ReflectedPlane;
  Plane okPlane = targetPlane;
  
  if (startingState.globalPosition().z()>0 || 
      (fabs(startingState.globalPosition().z())<1.e-4 && startingState.globalMomentum().z()>1.e-4)) {
    GlobalTrajectoryParameters gtp( GlobalPoint(startingState.globalPosition().x(),
                                                startingState.globalPosition().y(),
                                                -startingState.globalPosition().z()),
                                    GlobalVector(startingState.globalMomentum().x(),
                                                 startingState.globalMomentum().y(),
                                                 -startingState.globalMomentum().z()),
                                    startingState.globalParameters().charge(), prop.magneticField() );
    
    FreeTrajectoryState fts(gtp);
    TSOS ReflectedState( fts, startingState.surface());
    okState = ReflectedState;
      
    GlobalVector TempVec = targetPlane.normalVector();
    GlobalVector ReflectedVector = GlobalVector ( TempVec.x(), TempVec.y(), -TempVec.z());
    GlobalVector zAxis = ReflectedVector.unit();
    GlobalVector yAxis( zAxis.y(), -zAxis.x(), 0);
    GlobalVector xAxis = yAxis.cross( zAxis);
    Surface::RotationType rot = Surface::RotationType( xAxis, yAxis, zAxis);

    GlobalPoint gp = targetPlane.position();
    GlobalPoint gpSym(gp.x(), gp.y(), -gp.z());
    PlaneBuilder pb;
    ReflectedPlane = pb.plane( gpSym, rot);
    
    okPlane =  *ReflectedPlane;
    isReflected = true;    
  }
  
  // Done Reflecting


  TsosWP res = prop.propagateWithPath( okState, okPlane);

  // Now reflect back the result if necessary...
  if (isReflected && res.first.isValid()) { // reflect back... nobody should know we secretely z-reflected the tsos 
   

    TSOS TempState = res.first; 
    GlobalTrajectoryParameters gtp( GlobalPoint(TempState.globalPosition().x(),
                                                TempState.globalPosition().y(),
                                                -TempState.globalPosition().z()),
                                    GlobalVector(TempState.globalMomentum().x(),
                                                 TempState.globalMomentum().y(),
                                                 -TempState.globalMomentum().z()),
                                    TempState.globalParameters().charge(), prop.magneticField() );
    
    FreeTrajectoryState fts(gtp);
    TSOS ReflectedState( fts, TempState.surface());
    res.first = ReflectedState;
    isReflected = false;
  }
  // Done Reflecting back !

  if (res.first.isValid()) {

    GlobalPoint gp = res.first.globalPosition();
    GlobalPoint gpSym(gp.x(), gp.y(), (gp.z()<0? gp.z() : -gp.z()));
 
    if (currentVolume->inside( gpSym )) {
      return res;
    }
    // sometimes fails when propagation on plane and surface are less than 0.1 mm apart... 
  } 

  //  return TsosWP( TrajectoryStateOnSurface(), 0); // Gives 'double free' errors... return res even when Outside volume... 
  return res;

}

std::pair< TrajectoryStateOnSurface, double > NavPropagator::propagateWithPath	(	const TrajectoryStateOnSurface &	inputState,
		const Plane &	plane
	)		const `[virtual]`

propagation of TSOS to plane with path length

Reimplemented from Propagator.

Definition at line 45 of file NavPropagator.cc.

References GlobalTrajectoryParameters::charge(), prof2calltree::count, gather_cfg::cout, NavVolume::crossToNextVolume(), destinationCrossed(), findVolume(), TrajectoryStateOnSurface::globalMomentum(), TrajectoryStateOnSurface::globalParameters(), TrajectoryStateOnSurface::globalPosition(), TrajectoryStateOnSurface::isValid(), LogDebug, noNextVolume(), VolumeCrossReturnType::path(), GloballyPositioned< T >::position(), propagateInVolume(), TrajectoryStateOnSurface::surface(), theField, VolumeCrossReturnType::tsos(), Vector3DBase< T, FrameTag >::unit(), VolumeCrossReturnType::volume(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by propagate().

{

   LogDebug("NavPropagator") <<  "NavPropagator::propagateWithPath(TrajectoryStateOnSurface, Plane) called with "
       << inputState;

  VolumeCrossReturnType exitState( 0, inputState, 0.0);
  TSOS startingState = inputState;
  TSOS TempState = inputState;
  const NavVolume* currentVolume = 0; // explicit initialization to get rid of compiler warning

  int count = 0;
  int maxCount = 100;
  while (!destinationCrossed( startingState, TempState, targetPlane)) {

     LogDebug("NavPropagator") <<  "NavPropagator:: at beginning of while loop at iteration " << count ;

    startingState = TempState;
    bool isReflected = false;

    if (exitState.volume() != 0) { // next volume connected
      currentVolume = exitState.volume();
    }
    else {
      // 1 mm shifted peek is necessary only because geometry is not yet glued and exists only for z<0
      GlobalPoint gp = startingState.globalPosition() + 0.1*startingState.globalMomentum().unit();

      if (gp.z()>0) { // Reflect in Z (as long as MagVolumes exist only for z<0
        GlobalTrajectoryParameters gtp( GlobalPoint(startingState.globalPosition().x(),
                                                       startingState.globalPosition().y(),
                                                    -startingState.globalPosition().z()),
                                        GlobalVector(startingState.globalMomentum().x(),
                                                       startingState.globalMomentum().y(),
                                                       -startingState.globalMomentum().z()),
                                        startingState.globalParameters().charge(), theField );
        FreeTrajectoryState fts(gtp);
        TempState = TSOS( fts, startingState.surface());
        isReflected = true;
      }
      currentVolume = findVolume( TempState );
      if (currentVolume == 0) {
         std::cout <<  "NavPropagator: findVolume failed to find volume containing point "
             << startingState.globalPosition() ;
        return std::pair<TrajectoryStateOnSurface,double>(noNextVolume( startingState), 0);
      }
    }

    PropagatorType currentPropagator( *currentVolume);

     LogDebug("NavPropagator") <<  "NavPropagator: calling crossToNextVolume" ;
    exitState = currentVolume->crossToNextVolume( TempState, currentPropagator);
     LogDebug("NavPropagator") <<  "NavPropagator: crossToNextVolume returned" ;
     LogDebug("NavPropagator") <<  "Volume pointer: " << exitState.volume() << " and new ";
     LogDebug("NavPropagator") <<  exitState.tsos() ;
     LogDebug("NavPropagator") <<  "So that was a path length " << exitState.path() ;


    if ( !exitState.tsos().isValid()) {
      // return propagateInVolume( currentVolume, startingState, targetPlane);
       std::cout <<  "NavPropagator: failed to crossToNextVolume in volume at pos "
           << currentVolume->position();
      break;
    }
    else {
       LogDebug("NavPropagator") <<  "NavPropagator: crossToNextVolume reached volume ";
      if (exitState.volume() != 0) {       
         LogDebug("NavPropagator") <<  " at pos "
             << exitState.volume()->position() 
             << "with state " << exitState.tsos() ;
      }
      else  LogDebug("NavPropagator") <<  " unknown";
    }

    TempState = exitState.tsos();

    if (fabs(exitState.path())<0.01) {
      LogDebug("NavPropagator") <<  "failed to move at all!! at position: " << exitState.tsos().globalPosition();

      GlobalTrajectoryParameters gtp( exitState.tsos().globalPosition()+0.01*exitState.tsos().globalMomentum().unit(),
                                      exitState.tsos().globalMomentum(),
                                      exitState.tsos().globalParameters().charge(), theField );

      FreeTrajectoryState fts(gtp);
      TSOS ShiftedState( fts, exitState.tsos().surface());
      //      exitState.tsos() = ShiftedState;
      TempState = ShiftedState;

      LogDebug("NavPropagator") <<  "Shifted to new position " << TempState.globalPosition();

    }


    // std::cout << "Just moved " << exitState.path() << " cm through " << ( currentVolume->isIron()? "IRON":"AIR" ); 
    // std::cout << " at radius: " << TempState.globalPosition().perp() << std::endl;
    // reflect back to normal universe if necessary:
    if (isReflected) { // reflect back... nobody should know we secretely z-reflected the tsos 
      

      GlobalTrajectoryParameters gtp( GlobalPoint(TempState.globalPosition().x(),
                                                     TempState.globalPosition().y(),
                                                  -TempState.globalPosition().z()),
                                      GlobalVector(TempState.globalMomentum().x(),
                                                   TempState.globalMomentum().y(),
                                                   -TempState.globalMomentum().z()),
                                      TempState.globalParameters().charge(), theField );

      FreeTrajectoryState fts(gtp);
      TSOS ReflectedState( fts, TempState.surface());
      TempState = ReflectedState;
      isReflected = false;
    }

    ++count;
    if (count > maxCount) {
       LogDebug("NavPropagator") <<  "Ohoh, NavPropagator in infinite loop, count = " << count;
      return TsosWP();
    }

  }

  

  // Arriving here only if destinationCrossed or if crossToNextVolume returned invalid state,
  // in which case we should check if the targetPlane is not crossed in the current volume

  return propagateInVolume( currentVolume, startingState, targetPlane);  

}

std::pair< TrajectoryStateOnSurface, double > NavPropagator::propagateWithPath	(	const FreeTrajectoryState &	fts,
		const Cylinder &	cylinder
	)		const `[virtual]`

propagation to cylinder with path length

Implements Propagator.

Definition at line 347 of file NavPropagator.cc.

References LogDebug.

{
   LogDebug("NavPropagator") <<  "propagateWithPath(const FreeTrajectoryState&) not implemented yet in NavPropagator" ;
  return TsosWP();
}

std::pair< TrajectoryStateOnSurface, double > NavPropagator::propagateWithPath	(	const FreeTrajectoryState &	fts,
		const Plane &	plane
	)		const `[virtual]`

propagation of FTS to plane with path length

Implements Propagator.

Definition at line 354 of file NavPropagator.cc.

References LogDebug.

{
   LogDebug("NavPropagator") <<  "propagateWithPath(const FreeTrajectoryState&) not implemented yet in NavPropagator" ;
  return TsosWP();
}

Member Data Documentation

bool NavPropagator::isIronVolume[272] [private]

Definition at line 78 of file NavPropagator.h.

Referenced by NavPropagator(), and navVolume().

const VolumeBasedMagneticField* NavPropagator::theField [private]

Definition at line 76 of file NavPropagator.h.

Referenced by findVolume(), magneticField(), NavPropagator(), and propagateWithPath().

MagVolumeMap NavPropagator::theNavVolumeMap [mutable, private]

Definition at line 77 of file NavPropagator.h.

Referenced by navVolume(), and ~NavPropagator().

NavPropagator Class Reference

Public Member Functions

Private Types

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation