DecayTreeFitter::InternalParticle Class Reference

#include <InternalParticle.h>

Inheritance diagram for DecayTreeFitter::InternalParticle:

Public Types
enum	ParticleType {   kInteractionPoint, kRecoComposite, kRecoResonance, kInternalParticle,   kRecoTrack, kResonance, kRecoPhoton, kMissingParticle }
typedef std::vector < ParticleBase * >	ParticleContainer
typedef std::vector < ParticleBase * >	daucontainer
typedef daucontainer::const_iterator	const_iterator
typedef std::vector< std::pair < const ParticleBase *, int > >	indexmap
typedef std::vector < DecayTreeFitter::Constraint >	constraintlist

Public Member Functions
	InternalParticle (RhoCandidate *bc, const ParticleBase *mother, const Configuration &config)
virtual int	dim () const
virtual ErrCode	initPar1 (FitParams *)
virtual ErrCode	initPar2 (FitParams *)
virtual int	type () const
virtual int	posIndex () const
virtual int	lenIndex () const
virtual int	momIndex () const
virtual bool	hasEnergy () const
virtual bool	hasPosition () const
virtual std::string	parname (int index) const
ErrCode	projectKineConstraint (const FitParams *, Projection &) const
ErrCode	projectLifeTimeConstraint (const FitParams *, Projection &) const
ErrCode	projectConversionConstraint (const FitParams *, Projection &p) const
ErrCode	projectMassConstraintTwoBody (const FitParams *fitparams, Projection &p) const
virtual ErrCode	projectConstraint (Constraint::Type type, const FitParams *fitparams, Projection &p) const
virtual void	addToConstraintList (constraintlist &alist, int depth) const
virtual void	updateIndex (int &offset)
virtual ErrCode	initCov (FitParams *) const
virtual void	print (const FitParams *) const
const ParticleBase *	locate (RhoCandidate *bc) const
RhoCandidate *	particle () const
virtual int	index () const
const ParticleBase *	mother () const
const std::string &	name () const
virtual ErrCode	projectGeoConstraint (const FitParams *, Projection &) const
virtual ErrCode	projectMassConstraint (const FitParams *, Projection &) const
int	eneIndex () const
virtual double	chiSquareD (const FitParams *) const
double	pdtMass () const
double	pdtWidth () const
double	pdtCLifeTime () const
double	pdtTau () const
int	charge () const
const daucontainer &	daughters () const
const_iterator	begin () const
const_iterator	end () const
ParticleBase *	addDaughter (RhoCandidate *, const Configuration &config)
void	removeDaughter (const ParticleBase *pb)
virtual void	retrieveIndexMap (indexmap &anindexmap) const
void	setMother (const ParticleBase *m)
virtual int	nFinalChargedCandidates () const
void	setParticle (RhoCandidate *bc)
void	collectVertexDaughters (daucontainer &particles, int posindex)
bool	setMassConstraint (bool add)
void	setMassConstraint (double mass)
ChiSquare	chiSquare (const FitParams *fitparams) const

Static Public Member Functions
static ParticleBase *	createParticle (RhoCandidate *bc, const ParticleBase *mother, const Configuration &config)

Protected Member Functions
ErrCode	initMom (FitParams *fitparams) const
daucontainer &	daughters ()
ErrCode	initTau (FitParams *par) const
void	makeName (RhoCandidate *bc)
bool	hasMassConstraint () const
void	setIndex (int i)
void	setName (const std::string &n)

Static Protected Member Functions
static double	pdtCLifeTime (RhoCandidate *bc)
static bool	isAResonance (const TParticlePDG *bc)
static double	bFieldOverC ()

Private Attributes
bool	m_lifetimeconstraint
bool	m_isconversion

Detailed Description

Definition at line 18 of file InternalParticle.h.

Member Typedef Documentation

typedef daucontainer::const_iterator DecayTreeFitter::ParticleBase::const_iterator

inherited

Definition at line 97 of file ParticleBase.h.

typedef std::vector<DecayTreeFitter::Constraint> DecayTreeFitter::ParticleBase::constraintlist

inherited

Definition at line 110 of file ParticleBase.h.

typedef std::vector<ParticleBase*> DecayTreeFitter::ParticleBase::daucontainer

inherited

Definition at line 96 of file ParticleBase.h.

typedef std::vector< std::pair<const ParticleBase*,int> > DecayTreeFitter::ParticleBase::indexmap

inherited

Definition at line 105 of file ParticleBase.h.

typedef std::vector<ParticleBase*> DecayTreeFitter::ParticleBase::ParticleContainer

inherited

Definition at line 33 of file ParticleBase.h.

Member Enumeration Documentation

enum DecayTreeFitter::ParticleBase::ParticleType

inherited

Enumerator
kInteractionPoint
kRecoComposite
kRecoResonance
kInternalParticle
kRecoTrack
kResonance
kRecoPhoton
kMissingParticle

Definition at line 28 of file ParticleBase.h.

                      {kInteractionPoint,
      kRecoComposite,kRecoResonance,
      kInternalParticle,kRecoTrack,
      kResonance,kRecoPhoton,
      kMissingParticle} ;

Constructor & Destructor Documentation

DecayTreeFitter::InternalParticle::InternalParticle	(	RhoCandidate *	bc,
		const ParticleBase *	mother,
		const Configuration &	config
	)

Definition at line 30 of file InternalParticle.cxx.

References DecayTreeFitter::ParticleBase::addDaughter(), RhoCandidate::Daughter(), DecayTreeFitter::ParticleBase::daughters(), m_isconversion, m_lifetimeconstraint, RhoCandidate::NDaughters(), and RhoCandidate::PdgCode().

: ParticleBase(bc,aMother),m_lifetimeconstraint(false)
{
  //    BOOST_FOREACH( RhoCandidate* daughter, bc.daughters() )
  //      addDaughter(*daughter,config) ;
  for( int iDau=0;iDau<bc->NDaughters();iDau++ ){
    addDaughter(const_cast<RhoCandidate*>(bc->Daughter(iDau)),config) ;
  }
  // copy constraints
  m_lifetimeconstraint = false ; //bc && bc->constraint(BtaConstraint::Life) ;
  m_isconversion = daughters().size()==2 &&  bc->PdgCode() == 22 ;
}

Member Function Documentation

ParticleBase * DecayTreeFitter::ParticleBase::addDaughter ( RhoCandidate *  cand, const Configuration &  config  )

inherited

Definition at line 94 of file ParticleBase.cxx.

References DecayTreeFitter::ParticleBase::createParticle().

Referenced by InternalParticle().

{
  m_daughters.push_back( DecayTreeFitter::ParticleBase::createParticle(cand,this,config) ) ;
  return m_daughters.back() ;
}

void DecayTreeFitter::InternalParticle::addToConstraintList ( constraintlist &  alist, int  depth  ) const

virtual

Implements DecayTreeFitter::ParticleBase.

Reimplemented in DecayTreeFitter::InteractionPoint.

Definition at line 509 of file InternalParticle.cxx.

References DecayTreeFitter::Constraint::conversion, DecayTreeFitter::Constraint::geometric, DecayTreeFitter::Constraint::kinematic, DecayTreeFitter::Constraint::lifetime, and DecayTreeFitter::Constraint::mass.

Referenced by DecayTreeFitter::InteractionPoint::addToConstraintList().

{
  // first the daughters
  for(daucontainer::const_iterator it = daughters().begin() ;
      it != daughters().end() ; ++it)
    (*it)->addToConstraintList(alist,depth-1) ;

  //double geoprecision  = 1e-5 ; // 1mu
  //double massprecision = 4*pdtMass()*pdtMass()*1e-5 ; // 0.01 MeV

  // the lifetime constraint
  if(lenIndex()>=0 && m_lifetimeconstraint)
    alist.push_back(Constraint(this,Constraint::lifetime,depth,1)) ;
  // the kinematic constraint
  if(momIndex()>=0)
    alist.push_back(Constraint(this,Constraint::kinematic,depth,4)) ;
  // the geometric constraint
  if(mother() && lenIndex()>=0)
    alist.push_back(Constraint(this,Constraint::geometric,depth,3,3)) ;
  // the mass constraint. FIXME: move to ParticleBase
  if(hasMassConstraint()) {
    if( !m_isconversion )
      alist.push_back(Constraint(this,Constraint::mass,depth,1,10)) ;
    else
      alist.push_back(Constraint(this,Constraint::conversion,depth,1,3)) ;
  }
}

const_iterator DecayTreeFitter::ParticleBase::begin ( ) const

inlineinherited

Definition at line 100 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_daughters.

{ return m_daughters.begin() ; }

double DecayTreeFitter::ParticleBase::bFieldOverC ( )

staticprotectedinherited

Definition at line 582 of file ParticleBase.cxx.

References C(), RhoCalculationTools::GetBz(), and pos.

{ 
  TVector3 pos(0.,0.,0.); //TODO get a sensible position, but let's assume zero first...
  return RhoCalculationTools::GetBz ( pos ) / TMath::C() ; 
}

int DecayTreeFitter::ParticleBase::charge ( ) const

inlineinherited

Definition at line 90 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_charge.

{ return m_charge ; }

ChiSquare DecayTreeFitter::ParticleBase::chiSquare ( const FitParams *  fitparams ) const

inherited

Definition at line 568 of file ParticleBase.cxx.

References DecayTreeFitter::FitParams::chiSquare().

Referenced by DecayTreeFitter::DecayChain::chiSquare().

{
  ChiSquare chi2 ;
  // add contribution from daughters
  for(daucontainer::const_iterator it = m_daughters.begin() ;
      it != m_daughters.end() ; ++it) {
    chi2 += (*it)->chiSquare(fitparams) ;
  }
  // add own chisquare, adjust for number of parameters
  chi2 += fitparams->chiSquare( *this ) ;
  chi2 += ChiSquare( 0, -dim() ) ;
  return chi2 ;
}

double DecayTreeFitter::ParticleBase::chiSquareD ( const FitParams *  fitparams ) const

virtualinherited

Definition at line 551 of file ParticleBase.cxx.

{
  double rc = 0;
  for(daucontainer::const_iterator it = m_daughters.begin() ;
      it != m_daughters.end(); ++it)
    rc += (*it)->chiSquareD(fitparams) ;
  return rc ;
}

void DecayTreeFitter::ParticleBase::collectVertexDaughters ( daucontainer &  particles, int  posindex  )

inherited

Definition at line 227 of file ParticleBase.cxx.

References name, and vtxverbose.

{
  // collect all particles emitted from vertex with position posindex
  if(vtxverbose>=3) {
    std::cout << "DecayTreeFitter::ParticleBase::collectVertexDaughters " << posindex << std::endl ;
  }
  //skip: head of tree, particles from different vertex, resonances
  if( mother() && mother()->posIndex() == posindex && type()!=kRecoResonance && type()!=kResonance )
  {
    particles.push_back( this ) ;
    if(vtxverbose>=3) {
      std::cout << "DecayTreeFitter::ParticleBase::collectVertexDaughters - added a particle "<<name()<<" to vertex " << posindex << std::endl ;
    }
  }
  for( daucontainer::const_iterator idau = daughters().begin() ;
      idau != daughters().end() ; ++idau )
    (*idau)->collectVertexDaughters(particles,posindex ) ; //FIXME: RK Caution here!!
    //collectVertexDaughters(particles,posindex ) ;
}

ParticleBase * DecayTreeFitter::ParticleBase::createParticle ( RhoCandidate *  bc, const ParticleBase *  mother, const Configuration &  config  )

staticinherited

Definition at line 113 of file ParticleBase.cxx.

References RhoCandidate::Charge(), fabs(), RhoCandidate::GetRecoCandidate(), DecayTreeFitter::ParticleBase::index(), RhoCandidate::IsLocked(), DecayTreeFitter::ParticleBase::name(), RhoCandidate::NDaughters(), RhoCandidate::PdgCode(), RhoCandidate::PdtEntry(), DecayTreeFitter::ParticleBase::type(), RhoCandidate::Uid(), and vtxverbose.

Referenced by DecayTreeFitter::ParticleBase::addDaughter().

{
  // This routine interpretes a beta candidate as one of the
  // 'Particles' used by the fitter.

  const TParticlePDG* prop = particle->PdtEntry();

  if(vtxverbose>=5)
    std::cout << "DecayTreeFitter::ParticleBase::createParticle from     " <<particle->PdgCode() << " | " << particle->Uid() << std::endl ;
  ParticleBase* rc=0 ;
  //bool bsconstraint = false ;

  // We refit invalid fits, kinematic fits and composites with beamspot
  // constraint if not at head of tree.
  bool validfit = particle->IsLocked();//TODO needed? && particle.endVertex() != 0 ;
  bool iscomposite = (particle->NDaughters()>0) ;
  bool isresonance = iscomposite && prop && isAResonance(prop) ;

  if(!mother) { // 'head of tree' particles
                //if ( bsconstraint )
                //rc = new InteractionPoint(particle) ;
                //else

    if( iscomposite )
    { if(vtxverbose>=2) std::cout<<" H "; if(vtxverbose>=5) std::cout<<std::endl;
      rc = new InternalParticle(particle,0,config) ; // still need proper head-of-tree class
    }
    else {
      std::cout << "DecayTreeFitter::ParticleBase::createParticle: You are fitting a decay tree that exists of "
        << "a single, non-composite particle and which does not have a beamconstraint."
        << "I do not understand what you want me to do with this." << std::endl ;
      std::cout<<" X ";
      rc = new InternalParticle(particle,0,config) ; // still need proper head-of-tree class
    }
  } else if( !iscomposite ) { // external particles (i.e. tracks/neutrals)
    bool hasrecocand = ( NULL !=particle->GetRecoCandidate() );
    bool ischarged = ( fabs(particle->Charge())>0 );
    //bool hastrack = proto && (proto->GetTrackIndex() >= 0) ;
    //bool hascalo  = proto && (proto->GetEmcNumberOfCrystals() > 0) ;
    //std::cout<<" particle="<<particle<<" proto="<<proto<<" ntrk="<<proto->GetTrackIndex()<<" ncry"<<proto->GetEmcNumberOfCrystals()<<std::endl;
    if( hasrecocand )
    {
      if( ischarged )
      { if(vtxverbose>=2) std::cout<<" T "; if(vtxverbose>=5) std::cout<<std::endl;
        rc = new RecoTrack(particle,mother,config) ;  // reconstructed track
      } else
      { if(vtxverbose>=2) std::cout<<" P "; if(vtxverbose>=5) std::cout<<std::endl;
        rc = new RecoPhoton(particle,mother) ; // reconstructed photon
      }
    } else if( validfit ) {  // fitted composites w/o daughters?
      if( isresonance )
      { if(vtxverbose>=2) std::cout<<" RF "; if(vtxverbose>=5) std::cout<<std::endl;
        rc = new RecoResonance(particle,mother) ;
      } else
      { if(vtxverbose>=2) std::cout<<" CF "; if(vtxverbose>=5) std::cout<<std::endl;
        rc = new RecoComposite(particle,mother) ;
      }
    } else { // missing particle!
      if(vtxverbose>=2) std::cout<<" M "; if(vtxverbose>=5) std::cout<<std::endl;
      rc = new MissingParticle(particle,mother) ;
    }
  } else { // 'internal' particles
    if( validfit /*|| isconversion*/ ) {  // fitted composites
      if( isresonance ) { 
        if(vtxverbose>=2) std::cout<<" Rf "; if(vtxverbose>=5) std::cout<<std::endl;
        rc = new RecoResonance(particle,mother) ;
      } else { 
        if(vtxverbose>=2) std::cout<<" Cf "; if(vtxverbose>=5) std::cout<<std::endl;
        rc = new RecoComposite(particle,mother) ;
      }
    } else {         // unfited composites
      if( isresonance ) { 
        if(vtxverbose>=2) std::cout<<" R "; if(vtxverbose>=5) std::cout<<std::endl;
        rc = new Resonance(particle,mother,config) ;
      } else { 
        if(vtxverbose>=2) std::cout<<" I "; if(vtxverbose>=5) std::cout<<std::endl;
        rc = new InternalParticle(particle,mother,config) ;
      }
    }
  }
  if(vtxverbose>=5)
    std::cout << "DecayTreeFitter::ParticleBase::createParticle finished " <<particle->PdgCode() << " | " << particle->Uid() << std::endl ;

  if(vtxverbose>=2)
    std::cout << "DecayTreeFitter::ParticleBase::createParticle returns type=" << rc->type()
    << " index=" << rc->index() << " with name \""<< rc->name() << "\""<<std::endl ;
  return rc ;
}

const daucontainer& DecayTreeFitter::ParticleBase::daughters ( ) const

inlineinherited

Definition at line 99 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_daughters.

Referenced by DecayTreeFitter::InteractionPoint::InteractionPoint(), and InternalParticle().

{ return m_daughters ; }

daucontainer& DecayTreeFitter::ParticleBase::daughters ( )

inlineprotectedinherited

Definition at line 138 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_daughters.

{ return m_daughters ; }

virtual int DecayTreeFitter::InternalParticle::dim ( ) const

inlinevirtual

Implements DecayTreeFitter::ParticleBase.

Reimplemented in DecayTreeFitter::InteractionPoint, and DecayTreeFitter::Resonance.

Definition at line 24 of file InternalParticle.h.

References DecayTreeFitter::ParticleBase::mother().

{ return mother() ? 8 : 7 ; }

const_iterator DecayTreeFitter::ParticleBase::end ( ) const

inlineinherited

Definition at line 101 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_daughters.

Referenced by DecayTreeFitter::ParticleBase::locate().

{ return m_daughters.end() ; }

int DecayTreeFitter::ParticleBase::eneIndex ( ) const

inlineinherited

Definition at line 80 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::hasEnergy(), and DecayTreeFitter::ParticleBase::momIndex().

{ return hasEnergy() ? momIndex()+3 : -1 ; }

virtual bool DecayTreeFitter::InternalParticle::hasEnergy ( ) const

inlinevirtual

Reimplemented from DecayTreeFitter::ParticleBase.

Definition at line 34 of file InternalParticle.h.

{ return true ; }

bool DecayTreeFitter::ParticleBase::hasMassConstraint ( ) const

inlineprotectedinherited

Definition at line 139 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_hasMassConstraint.

Referenced by DecayTreeFitter::MissingParticle::dim(), and DecayTreeFitter::MissingParticle::hasEnergy().

{ return m_hasMassConstraint ; }

virtual bool DecayTreeFitter::InternalParticle::hasPosition ( ) const

inlinevirtual

Reimplemented from DecayTreeFitter::ParticleBase.

Reimplemented in DecayTreeFitter::Resonance.

Definition at line 35 of file InternalParticle.h.

{ return true ; }

virtual int DecayTreeFitter::ParticleBase::index ( ) const

inlinevirtualinherited

Definition at line 59 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_index.

Referenced by DecayTreeFitter::ParticleBase::createParticle(), DecayTreeFitter::DecayChain::index(), lenIndex(), DecayTreeFitter::RecoComposite::lenIndex(), DecayTreeFitter::MissingParticle::momIndex(), DecayTreeFitter::RecoParticle::momIndex(), DecayTreeFitter::RecoResonance::momIndex(), DecayTreeFitter::Resonance::momIndex(), momIndex(), DecayTreeFitter::RecoComposite::momIndex(), posIndex(), and DecayTreeFitter::RecoComposite::posIndex().

{ return m_index ; }

ErrCode DecayTreeFitter::ParticleBase::initCov ( FitParams *  fitparams ) const

virtualinherited

Reimplemented in DecayTreeFitter::InteractionPoint, and DecayTreeFitter::RecoPhoton.

Definition at line 248 of file ParticleBase.cxx.

References DecayTreeFitter::FitParams::cov(), name, row, status, and vtxverbose.

Referenced by DecayTreeFitter::Fitter::add(), and DecayTreeFitter::InteractionPoint::initCov().

{
  ErrCode status ;

  if(vtxverbose>=2) {
    std::cout << "DecayTreeFitter::ParticleBase::initCov for " << name() << std::endl ;
  }

  // position
  int posindex = posIndex() ;
  if( posindex>=0 ) {
    const double sigx   = 10; // * Gaudi::Units::cm ;
    const double sigy   = 10; // * Gaudi::Units::cm ;
    const double sigz   = 50; // * Gaudi::Units::cm ;
    fitparams->cov()(posindex+0,posindex+0) = sigx*sigx ;
    fitparams->cov()(posindex+1,posindex+1) = sigy*sigy ;
    fitparams->cov()(posindex+2,posindex+2) = sigz*sigz ;
  }

  // momentum
  int momindex = momIndex() ;
  if(momindex>=0) {
    // TODO: calo at high energy?
    const double sigmom = 10; // * Gaudi::Units::GeV ; // GeV
    int maxrow = hasEnergy() ? 4 : 3 ;
    for(int row=momindex; row<momindex+maxrow; row++)
      fitparams->cov()(row,row) = sigmom*sigmom ;
  }

  // lifetime
  int lenindex = lenIndex() ;
  if(lenindex>=0) {
    const double sigz   = 50; // * Gaudi::Units::cm ;
    fitparams->cov()(lenindex,lenindex) = sigz*sigz ;
  }

  for(daucontainer::const_iterator it = m_daughters.begin() ;
      it != m_daughters.end() ; ++it)
    status |= (*it)->initCov(fitparams) ;
  return status ;
}

ErrCode DecayTreeFitter::InternalParticle::initMom ( FitParams *  fitparams ) const

protected

Definition at line 275 of file InternalParticle.cxx.

References DecayTreeFitter::FitParams::par(), sqrt(), DecayTreeFitter::ErrCode::success, and vtxverbose.

{
  int momindex = momIndex() ;
  // reset
  for( int irow=0; irow<4; irow++)
    fitparams->par(momindex+irow) = 0 ;

  // now add daughter momenta
  for(daucontainer::const_iterator it = begin() ; it != end() ; ++it) {
    int daumomindex = (*it)->momIndex() ;
    double e2(0) ;
    int maxrow = (*it)->hasEnergy() && !hasMassConstraint() ? 4 : 3 ;
    for(int irow=0; irow<maxrow; irow++) {
      double px = fitparams->par()(daumomindex+irow) ;
      e2 += px*px ;
      fitparams->par(momindex+irow) += px ;
    }
    if(maxrow==3) {
      double mass = (*it)->pdtMass() ;
      fitparams->par(momindex+3) += sqrt(e2+mass*mass) ;
      if(vtxverbose>6){std::cout<<"InternalParticle::initMom: mass="<<mass<<", p^2="<<e2<<", E="<<sqrt(e2+mass*mass)<<" from daughter momindex "<<daumomindex<<" p=("<<fitparams->par()(daumomindex+0)<<","<<fitparams->par()(daumomindex+1)<<","<<fitparams->par()(daumomindex+2)<<")"<<std::endl;}
    }
  }
  return ErrCode::success ;
}

ErrCode DecayTreeFitter::InternalParticle::initPar1 ( FitParams *  fitparams )

virtual

Implements DecayTreeFitter::ParticleBase.

Reimplemented in DecayTreeFitter::InteractionPoint, and DecayTreeFitter::Resonance.

Definition at line 44 of file InternalParticle.cxx.

References DecayTreeFitter::ErrCode::badsetup, DecayTreeFitter::closestPointParams(), DecayTreeFitter::compTrkTransverseMomentum(), fabs(), DecayTreeFitter::ParticleBase::kRecoTrack, mu1, name, p1, p2, DecayTreeFitter::FitParams::par(), particle, point, DecayTreeFitter::State::position(), row, DecayTreeFitter::RecoTrack::setFlightLength(), DecayTreeFitter::State::slopes(), DecayTreeFitter::RecoTrack::state(), status, and vtxverbose.

Referenced by DecayTreeFitter::InteractionPoint::initPar1().

{
  // This is the most complicated part of the vertexer: an
  // initialization that always works.

  // There are two ways out: If the RhoCandidate was vertexed
  // before, we can rely on the existing vertex (case A). If not, we
  // need to estimate the vertex position from the daughters; that
  // is very complicated (case B). The momentum is always
  // initialized from the sum of the daughter four-vectors. In the
  // end, it doesn't really matter.

  // FIXME: Currently, this scheme does not work for B->K*D0, with
  // D0->pi0Ks, because the D0 is initialized before there is a B
  // vertex.

  if(vtxverbose>5){std::cout<<"InternalParticle::initPar1(): - "<<std::endl;}
  if(vtxverbose>=3)
    std::cout << "InternalParticle::initPar1(): "
    << particle() << " "
    << hasPosition() << " " << posIndex() << std::endl ;

  ErrCode status ;
  int posindex = posIndex() ;

  // logic check: we do not want to call this routine for resonances.
  assert( hasPosition() ) ;
  if(vtxverbose>=3)
    std::cout << "InternalParticle::initPar1(): assertion passed"<<std::endl;

  // Start with origin
  for(int row=0; row<3; row++) fitparams->par()(row+posindex) = 0 ;

  // Step 1: pre-initialization of all daughters
  for(daucontainer::const_iterator it = begin() ; it != end() ; ++it)
    status |= (*it)->initPar1(fitparams) ;

  // Step 2: initialize the vertex. if we are lucky, we had a
  // 'resonant' daughter, and we are already done.
  if( fitparams->par()(posindex+0)==0 && fitparams->par()(posindex+1)==0 &&
     fitparams->par()(posindex+2)==0 ) {

    //const RhoVector3Err* vtx(0) ;
    RhoVector3Err vtx = particle()->DecayVtx();
    if(vtx.Mag()>0){
      if(vtxverbose>=3)    std::cout << "InternalParticle::initPar1(): A"<<std::endl;
      // we found an existing valid vertex. that's fine as well ...
      //TVector3 point = vtx->position() ;
      fitparams->par(posindex+0) = vtx.x() ;
      fitparams->par(posindex+1) = vtx.y() ;
      fitparams->par(posindex+2) = vtx.z() ;
      if(vtxverbose>=2)
        std::cout << "using existing vertex: " << vtx << std::endl ;

    } else {
      if(vtxverbose>=3)    std::cout << "InternalParticle::initPar1: B"<<std::endl;
      // Case B: the hard way ... use the daughters to estimate the
      // vertex. First we check if there are sufficient tracks
      // attached to this vertex. If so, estimate the poca of the
      // two tracks with the highest momentum. This will work for
      // the majority of the cases. If there are not sufficient
      // tracks, add the composites and take the two with the best
      // doca.

      // create a vector with all daughters that constitute a
      // 'trajectory' (ie tracks, composites and daughters of
      // resonances.)
      daucontainer alldaughters ;
      collectVertexDaughters( alldaughters, posindex ) ;
      if(vtxverbose>=3)    std::cout << "InternalParticle::initPar1(): number of daughters for initializing vertex: "
      << name() << " " << alldaughters.size() << std::endl ;

      if(vtxverbose>=3)    std::cout << "InternalParticle::initPar1(): B - r?"<<std::endl;
      // select daughters that are either charged, or have an initialized vertex
      daucontainer vtxdaughters ;
      std::vector<RecoTrack*> trkdaughters ;
      for(daucontainer::const_iterator it = alldaughters.begin() ;
        it != alldaughters.end() ; ++it) {
        if( (*it)->type()==ParticleBase::kRecoTrack ) {
          trkdaughters.push_back( static_cast<RecoTrack*>(*it) )  ;
        } else if( (*it)->hasPosition() && fitparams->par((*it)->posIndex()+0)!=0 ) {
          vtxdaughters.push_back( *it ) ;
        }
      }
      if(vtxverbose>=3)    std::cout << "InternalParticle::initPar1(): B -daughters?"<<std::endl;

      if( trkdaughters.size() >=2 ) {
        if(vtxverbose>=3)    std::cout << "InternalParticle::initPar1(): B -a?"<<std::endl;
        // sort in pT. not very efficient, but it works.
        if( trkdaughters.size()>2 )
          std::sort(trkdaughters.begin(),trkdaughters.end(),compTrkTransverseMomentum) ;
        // now, just take the first two ...
        RecoTrack* dau1 = trkdaughters[0] ;
        RecoTrack* dau2 = trkdaughters[1] ;

        //FIXME [R.K.03/2017] We at Panda don't have so many hard tracks.
        // We 
        // get the poca of the two statevectors
        const DecayTreeFitter::State& state1 = dau1->state() ;
        const DecayTreeFitter::State& state2 = dau2->state() ;
        DecayTreeFitter::Line line1(state1.position(),state1.slopes()) ;
        DecayTreeFitter::Line line2(state2.position(),state2.slopes()) ;
        double mu1(0),mu2(0) ;
        DecayTreeFitter::closestPointParams(line1,line2,mu1,mu2) ;
        TVector3 p1 = line1.position(mu1) ;
        TVector3 p2 = line2.position(mu2) ;
        fitparams->par()(posindex+0) = 0.5*(p1.x()+p2.x()) ;
        fitparams->par()(posindex+1) = 0.5*(p1.y()+p2.y()) ;
        fitparams->par()(posindex+2) = 0.5*(p1.z()+p2.z()) ;
        dau1->setFlightLength( mu1 ) ;
        dau2->setFlightLength( mu2 ) ;

      } else if(trkdaughters.size()+vtxdaughters.size()>=2)  {
        if(vtxverbose>=3)    std::cout << "InternalParticle::initPar1(): B -b?"<<std::endl;
        std::cout << "InternalParticle::initPar1(): InternalParticle: Trying new code!"<< std::endl ;
        // ###################

        // that's unfortunate: no enough charged tracks from this
        // vertex. need all daughters. create trajectories and use
        // normal TrkPoca.

        //std::vector<LHCb::StateZTraj> trajectories ;
        std::vector<DecayTreeFitter::Line> trajectories ;
        for(std::vector<RecoTrack*>::const_iterator it = trkdaughters.begin() ; it != trkdaughters.end() ; ++it)
        {
          const DecayTreeFitter::State& state1 = (*it)->state() ;
          trajectories.push_back( DecayTreeFitter::Line(state1.position(),state1.slopes()) );
          //trajectories.push_back( (*it)->traj() ) ;
        }

        //trajectories.push_back(&((*it)->particle().trkAbsFit()->traj())) ;

        //std::vector<DecayTreeFitter::Line> linetrajectories ; // store trajectories of composites
        //linetrajectories.reserve(  vtxdaughters.size() ) ;
        for(daucontainer::const_iterator it = vtxdaughters.begin() ; it != vtxdaughters.end() ; ++it) 
        {
          //std::cout << (*it)->particle().pdtEntry()->name() << std::endl ;
          int dauposindex = (*it)->posIndex() ;
          int daumomindex = (*it)->momIndex() ;
          TVector3 point(fitparams->par()(dauposindex+0), fitparams->par()(dauposindex+1), fitparams->par()(dauposindex+2)) ;
          TVector3 direction(fitparams->par()(daumomindex+0), fitparams->par()(daumomindex+1), fitparams->par()(daumomindex+2)) ;
          trajectories.push_back(DecayTreeFitter::Line(point,direction) ) ;
          //linetrajectories.push_back(DecayTreeFitter::Line(point,direction,1) ) ;
          //trajectories.push_back(&(linetrajectories.back())) ;
          //daupoint = point ;
        }

        // we select the two trajectories with the best poca
        double docabest(99999);
        for( std::vector<DecayTreeFitter::Line>::iterator it1 = trajectories.begin() ; it1 != trajectories.end(); ++it1 )
        {
          for( std::vector<DecayTreeFitter::Line>::iterator it2 = trajectories.begin() ; it2 != it1; ++it2 ) 
          {
            double mu1(0),mu2(0) ;
            DecayTreeFitter::closestPointParams(*it1,*it2,mu1,mu2) ;
            TVector3 p1 = (*it1).position(mu1) ;
            TVector3 p2 = (*it2).position(mu2) ;
            double doca = (p1-p2).Mag();
            if(fabs(doca)<fabs(docabest)) 
            {
              fitparams->par()(posindex+0) = 0.5*(p1.x()+p2.x()) ;
              fitparams->par()(posindex+1) = 0.5*(p1.y()+p2.y()) ;
              fitparams->par()(posindex+2) = 0.5*(p1.z()+p2.z()) ;
              docabest = doca ;
            }
          }
        }
  // ###########################33
      } else if( mother() && mother()->posIndex()>=0 ) {
        if(vtxverbose>=3)    std::cout << "InternalParticle::initPar1(): B -c?"<<std::endl;

        // let's hope the mother was initialized
        int posindexmother = mother()->posIndex() ;
        for(int ipos=0; ipos<3; ipos++) {
          fitparams->par()(posindex+ipos) =
          fitparams->par()(posindexmother+ipos) ;
        }
      } else {
        if(vtxverbose>=3)    std::cout << "InternalParticle::initPar1(): B -d?"<<std::endl;
       // something is wrong!
        std::cout << "InternalParticle::initPar1(): There are not sufficient geometric constraints to fit "
          << "this decay tree. Perhaps you should add a beam constraint. "
          << std::endl ;
        //<< particle().constraint(BtaConstraint::Beam)
        //   << std::endl
        //   << treeprinter.print(*bc()) << endmsg ;
        status |= ErrCode::badsetup ;
      }
    if(vtxverbose>=3)    std::cout << "InternalParticle::initPar1(): C ?"<<std::endl;
    }
  }
  if(vtxverbose>=3)
    std::cout << "InternalParticle::initPar1(): big code chunk passed"<<std::endl;

  // step 3: do the post initialization step of all daughters
  if(vtxverbose>5){std::cout<<"InternalParticle::initPar1(): - calling all initPar2()"<<std::endl;}
  for(daucontainer::const_iterator it = daughters().begin() ;
      it != daughters().end() ; ++it)
    (*it)->initPar2(fitparams) ;

  // step 4: initialize the momentum by adding up the daughter 4-vectors
  initMom(fitparams) ;

  if(vtxverbose>=3)
    std::cout << "InternalParticle::initPar1(): End of initpar: "
    << name() << " fitpar pos("
    << fitparams->par()(posindex+0) << ","
    << fitparams->par()(posindex+1) << ","
    << fitparams->par()(posindex+2) << ")" << std::endl ;

  return status ;
}

ErrCode DecayTreeFitter::InternalParticle::initPar2 ( FitParams *  fitparams )

virtual

Implements DecayTreeFitter::ParticleBase.

Reimplemented in DecayTreeFitter::Resonance.

Definition at line 258 of file InternalParticle.cxx.

References DecayTreeFitter::FitParams::par(), and vtxverbose.

{
  if(vtxverbose>5){std::cout<<"InternalParticle::initPar2: - "<<std::endl;}
  // FIXME: in the unfortunate case (the B-->D0K*- above) that our
  // vertex is still the origin, we copy the mother vertex.
  int posindex = posIndex() ;
  if( hasPosition() && mother() && fitparams->par(posindex+0)==0 &&
     fitparams->par(posindex+1)==0 && fitparams->par(posindex+2)==0 ) {
    int posindexmom = mother()->posIndex() ;
    for(int irow=0; irow<3; irow++)
      fitparams->par(posindex+irow) = fitparams->par(posindexmom+irow) ;
  }
  // step 5: initialize the lifetime
  return initTau(fitparams) ;
}

ErrCode DecayTreeFitter::ParticleBase::initTau ( FitParams *  par ) const

protectedinherited

Definition at line 522 of file ParticleBase.cxx.

References mom, DecayTreeFitter::FitParams::par(), DecayTreeFitter::ParticleBase::posIndex(), and DecayTreeFitter::ErrCode::success.

{
  int lenindex = lenIndex() ;
  if(lenindex>=0 && hasPosition() ) {
    const ParticleBase* amother = mother() ;
    int momposindex = amother ? amother->posIndex() : -1 ;
    int posindex = posIndex() ;
    int momindex = momIndex() ;
    assert(momposindex>=0) ; // check code logic: no mother -> no tau
                            //assert(fitparams->par(momposindex+0)!=0 ||fitparams->par(momposindex+1)!=0
                            //      ||fitparams->par(momposindex+2)!=0) ; // mother must be initialized

    TVector3 dX,mom ;
    double mom2(0) ;
    for(int irow=0; irow<3; irow++) {
      dX(irow)  = fitparams->par(posindex+irow) - fitparams->par(momposindex+irow) ;
      double px = fitparams->par(momindex+irow) ;
      mom(irow) = px ;
      mom2 += px*px ;
    }
    double tau = dX.Dot(mom)/mom2 ;
    // we don't like 0 and we don't like very negative values
    if( tau==0 ) tau=pdtTau() ;
    //tau = tau==0 ? pdtTau() : std::max(tau,-pdtTau()) ;
    fitparams->par(lenindex) = tau ;
  }
  return ErrCode::success ;
}

bool DecayTreeFitter::ParticleBase::isAResonance ( const TParticlePDG *  bc )

staticprotectedinherited

Definition at line 205 of file ParticleBase.cxx.

References C(), and vtxverbose.

                                                                  {
  bool rc = false ;
  switch(prop->PdgCode()) {
    case 11: // bremstrahlung is treated as a resonance
    case -11:
      rc = true ;
      break ;
    //case 22: // conversions are not treated as a resonance
    //case 310: // K shorts count as "stable" in PDT Table
      //rc = false;
      //break ;
    default: // this should take care of the pi0
      rc = (prop->Lifetime()>0) && (100.*TMath::C()*prop->Lifetime() < 0.0001); //[cm]
      if(rc) break;
      if(vtxverbose>4)std::cout<<"ParticleBase::isAResonance l."<<__LINE__<<": ctau="<<100.*TMath::C()*prop->Lifetime()<<" cm"<<std::endl;
      if(prop->Stable()) return false;
      //rc = prop.ctau() < 0.001*Gaudi::Units::mm ;
      break;
  }
  return rc ;
}

virtual int DecayTreeFitter::InternalParticle::lenIndex ( ) const

inlinevirtual

Reimplemented from DecayTreeFitter::ParticleBase.

Reimplemented in DecayTreeFitter::Resonance.

Definition at line 32 of file InternalParticle.h.

References DecayTreeFitter::ParticleBase::index(), and DecayTreeFitter::ParticleBase::mother().

{ return mother() ? index()+3 : -1 ; }

const ParticleBase * DecayTreeFitter::ParticleBase::locate ( RhoCandidate *  bc ) const

inherited

Definition at line 353 of file ParticleBase.cxx.

References DecayTreeFitter::ParticleBase::end(), and DecayTreeFitter::ParticleBase::locate().

Referenced by DecayTreeFitter::Fitter::decayLengthSum(), DecayTreeFitter::DecayChain::locate(), DecayTreeFitter::ParticleBase::locate(), DecayTreeFitter::Fitter::name(), and DecayTreeFitter::Fitter::updateCand().

{
  // does there exist an 'iscloneof' in lhcb?
  const ParticleBase* rc = m_particle == abc ? this : 0 ;
  for(daucontainer::const_iterator it = m_daughters.begin() ;
      !rc && it != m_daughters.end(); ++it)
    rc = (*it)->locate(abc) ;
  return rc ;
}

void DecayTreeFitter::ParticleBase::makeName ( RhoCandidate *  bc )

protectedinherited

virtual int DecayTreeFitter::InternalParticle::momIndex ( ) const

inlinevirtual

Reimplemented from DecayTreeFitter::ParticleBase.

Reimplemented in DecayTreeFitter::Resonance.

Definition at line 33 of file InternalParticle.h.

References DecayTreeFitter::ParticleBase::index(), and DecayTreeFitter::ParticleBase::mother().

{ return mother() ? index()+4 : index() + 3 ; }

const ParticleBase* DecayTreeFitter::ParticleBase::mother ( ) const

inlineinherited

Definition at line 60 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_mother.

Referenced by DecayTreeFitter::DecayChain::chiSquare(), DecayTreeFitter::Fitter::decayLengthSum(), dim(), lenIndex(), momIndex(), DecayTreeFitter::RecoResonance::posIndex(), DecayTreeFitter::Resonance::posIndex(), DecayTreeFitter::Fitter::remove(), and DecayTreeFitter::Fitter::updateCand().

{ return m_mother ; }

const std::string& DecayTreeFitter::ParticleBase::name ( ) const

inlineinherited

Definition at line 61 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_name.

Referenced by DecayTreeFitter::ParticleBase::createParticle(), and DecayTreeFitter::Fitter::name().

{ return m_name ; }

int DecayTreeFitter::ParticleBase::nFinalChargedCandidates ( ) const

virtualinherited

Reimplemented in DecayTreeFitter::RecoTrack.

Definition at line 560 of file ParticleBase.cxx.

Referenced by DecayTreeFitter::sortByType().

                                                               {
  int rc=0;
  for(daucontainer::const_iterator it = m_daughters.begin() ;
      it != m_daughters.end() ; ++it)
    rc +=(*it)->nFinalChargedCandidates() ;
  return rc ;
}

std::string DecayTreeFitter::InternalParticle::parname ( int  index ) const

virtual

Reimplemented from DecayTreeFitter::ParticleBase.

Reimplemented in DecayTreeFitter::Resonance.

Definition at line 539 of file InternalParticle.cxx.

References DecayTreeFitter::ParticleBase::parname().

{
  int id = thisindex ;
  // skip the lifetime parameter if there is no mother
  if(!mother() && id>=3) ++id ;
  return ParticleBase::parname(id) ;
}

RhoCandidate* DecayTreeFitter::ParticleBase::particle ( ) const

inlineinherited

Definition at line 57 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_particle.

Referenced by DecayTreeFitter::compTrkTransverseMomentum(), DecayTreeFitter::InteractionPoint::InteractionPoint(), DecayTreeFitter::DecayChain::locate(), and DecayTreeFitter::sortByType().

{ return m_particle ; }

double DecayTreeFitter::ParticleBase::pdtCLifeTime ( ) const

inlineinherited

Definition at line 88 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_pdtCLifeTime.

{ return m_pdtCLifeTime ; }

static double DecayTreeFitter::ParticleBase::pdtCLifeTime ( RhoCandidate *  bc )

staticprotectedinherited

double DecayTreeFitter::ParticleBase::pdtMass ( ) const

inlineinherited

Definition at line 86 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_pdtMass.

Referenced by projectMassConstraintTwoBody(), and DecayTreeFitter::Fitter::updateCand().

{ return m_pdtMass ; }

double DecayTreeFitter::ParticleBase::pdtTau ( ) const

inlineinherited

Definition at line 89 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_pdtCLifeTime, and DecayTreeFitter::ParticleBase::m_pdtMass.

{ return m_pdtMass >0 ? m_pdtCLifeTime/m_pdtMass : 0 ; }

double DecayTreeFitter::ParticleBase::pdtWidth ( ) const

inlineinherited

Definition at line 87 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_pdtWidth.

{ return m_pdtWidth ; }

virtual int DecayTreeFitter::InternalParticle::posIndex ( ) const

inlinevirtual

Reimplemented from DecayTreeFitter::ParticleBase.

Reimplemented in DecayTreeFitter::Resonance.

Definition at line 31 of file InternalParticle.h.

References DecayTreeFitter::ParticleBase::index().

{ return index()   ; }

void DecayTreeFitter::ParticleBase::print ( const FitParams *  fitpar ) const

virtualinherited

Definition at line 309 of file ParticleBase.cxx.

References DecayTreeFitter::FitParams::cov(), i, name, DecayTreeFitter::FitParams::par(), pz, and sqrt().

{
  std::cout << std::setw(5) << "[" << type() << "]" << std::setw(15) << std::flush << name().c_str()
  << std::setw(15)<< " val"  << std::setw(15)<< "err" << std::setw(15) << "sigma^2"<< std::endl ;
  std::cout << std::setprecision(5) ;
  for(int i=0; i<dim(); i++) {
    int theindex = index()+i ;
    std::cout << std::setw(2) << theindex << " "
    << std::setw(20) << parname(i).c_str()
    << std::setw(15) << fitpar->par()(theindex)
    << std::setw(15) << sqrt(fitpar->cov()(theindex,theindex))
    << std::setw(15) << fitpar->cov()(theindex,theindex) <<std::endl ;
  }
  if( hasEnergy() ) {
    int momindex = momIndex() ;
    double E  = fitpar->par()(momindex+3) ;
    double px = fitpar->par()(momindex+0) ;
    double py = fitpar->par()(momindex+1) ;
    double pz = fitpar->par()(momindex+2) ;
    double mass2 = E*E-px*px-py*py-pz*pz ;
    double mass = mass2>0 ? sqrt(mass2) : -sqrt(-mass2) ;

    // TODO this does not work? fitpar->cov()).GetSub(momindex+0,momindex+3)
    TMatrixDSym cov = fitpar->cov().GetSub(momindex+0,momindex+3,momindex+0,momindex+3);
    //HepSymMatrix cov = fitpar->cov().sub(momindex+0,momindex+3) ;
    TVectorD G(4) ; //was at G(4,0) .. ??why
    G(0) = -px/mass ;
    G(1) = -py/mass ;
    G(2) = -pz/mass ;
    G(3) =   E/mass ;
    double massvar = cov.Similarity(G) ;
    std::cout << std::setw(2)<< "-"<<" " << std::setw(20) << "mass: "
    << std::setw(15) << mass
    << std::setw(15) << sqrt(massvar)
    << std::setw(15) << massvar << std::endl ;
  }

  for(daucontainer::const_iterator it = m_daughters.begin() ;
      it != m_daughters.end() ; ++it)
    (*it)->print(fitpar) ;

}

ErrCode DecayTreeFitter::InternalParticle::projectConstraint ( Constraint::Type  type, const FitParams *  fitparams, Projection &  p  ) const

virtual

else

Reimplemented from DecayTreeFitter::ParticleBase.

Reimplemented in DecayTreeFitter::InteractionPoint.

Definition at line 419 of file InternalParticle.cxx.

References DecayTreeFitter::Constraint::conversion, DecayTreeFitter::Constraint::geometric, DecayTreeFitter::Projection::H(), DecayTreeFitter::Constraint::kinematic, DecayTreeFitter::Constraint::lifetime, DecayTreeFitter::Constraint::mass, DecayTreeFitter::Constraint::massEnergy, RhoCalculationTools::PrintMatrix(), DecayTreeFitter::ParticleBase::projectConstraint(), DecayTreeFitter::Projection::r(), DecayTreeFitter::Projection::V(), and vtxverbose.

Referenced by DecayTreeFitter::InteractionPoint::projectConstraint().

{
  ErrCode aStatus ;
  switch(aType) {
    case Constraint::mass:
    case Constraint::massEnergy:
      //       if( m_daughters.size()==2 &&
      //    !m_daughters.front()->hasEnergy() &&
      //    !m_daughters.back()->hasEnergy() )
      //  aStatus |= projectMassConstraintTwoBody(fitparams,p) ;
      aStatus |= projectMassConstraint(fitparams,p) ;
      //chisq = filterMassConstraintOnDaughters(fitpar) ;
      break ;
    case Constraint::geometric:
      aStatus |= projectGeoConstraint(fitparams,p) ;
      break ;
    case Constraint::kinematic:
      aStatus |= projectKineConstraint(fitparams,p) ;
      break ;
    case Constraint::lifetime:
      aStatus |= projectLifeTimeConstraint(fitparams,p) ;
      break ;
    case Constraint::conversion:
      aStatus |= projectConversionConstraint(fitparams,p) ;
      break ;
    default:
      aStatus |= ParticleBase::projectConstraint(aType,fitparams,p) ;
  }
  if(vtxverbose>6){
    std::cout<<"InternalParticle::projectConstraint(): projection is:"<<std::endl;
    std::cout<<"r "; p.r().Print();
    std::cout<<"V "; p.V().Print();
    std::cout<<"H "; RhoCalculationTools::PrintMatrix(p.H());
    }
  return aStatus ;
}

ErrCode DecayTreeFitter::InternalParticle::projectConversionConstraint	(	const FitParams *	fitparams,
		Projection &	p
	)		const

Definition at line 383 of file InternalParticle.cxx.

References DecayTreeFitter::Projection::H(), DecayTreeFitter::ParticleBase::momIndex(), DecayTreeFitter::FitParams::par(), DecayTreeFitter::Projection::r(), sqrt(), DecayTreeFitter::ErrCode::success, and vtxverbose.

{
   if(vtxverbose>6){std::cout<<"InternalParticle::projectConversionConstraint():"<<std::endl;}
  // only works if there are two daughters. constraint those to be parallel:
  // p1.in(p2) - |p1||p2|=0
  assert(m_isconversion) ;
  const ParticleBase* dauA = daughters()[0] ;
  const ParticleBase* dauB = daughters()[1] ;
  int daumomindexA = dauA->momIndex() ;
  int daumomindexB = dauB->momIndex() ;

  // first calculate the total momenta
  double momA2(0),momB2(0) ;
  for(int irow=0; irow<3; irow++) {
    double pxA =  fitparams->par(daumomindexA+irow) ;
    momA2 += pxA*pxA ;
    double pxB =  fitparams->par(daumomindexB+irow) ;
    momB2 += pxB*pxB ;
  }
  double momA(sqrt(momA2)), momB(sqrt(momB2)) ;

  // now fill r and H
  p.r(1) = -momA*momB ;
  for(int irow=0; irow<3; irow++) {
    double pxA =  fitparams->par(daumomindexA+irow) ;
    double pxB =  fitparams->par(daumomindexB+irow) ;
    p.r(1) += pxA*pxB ;
    p.H(1,daumomindexA+irow) = pxB - pxA/momA * momB ;
    p.H(1,daumomindexB+irow) = pxA - pxB/momB * momA ;
  }

  return ErrCode::success ;
}

ErrCode DecayTreeFitter::ParticleBase::projectGeoConstraint ( const FitParams *  fitparams, Projection &  p  ) const

virtualinherited

Definition at line 387 of file ParticleBase.cxx.

References DecayTreeFitter::Projection::H(), mom, p2, DecayTreeFitter::FitParams::par(), RhoCalculationTools::PrintMatrix(), pz, DecayTreeFitter::Projection::r(), row, sqrt(), DecayTreeFitter::ErrCode::success, DecayTreeFitter::Projection::V(), and vtxverbose.

{
  // implements the constraint

  //  vec{x}_decay = vec{x}_production + decaylength * vec{p} / p
  int posindexmother = mother()->posIndex() ;
  int posindex = posIndex();
  int lenindex = lenIndex() ;
  int momindex = momIndex() ;
  assert(posindexmother>=0 && posindex>=0 && lenindex>=0 && momindex>=0) ;


  // size of momentum
  double px = fitparams->par()(momindex+0) ;
  double py = fitparams->par()(momindex+1) ;
  double pz = fitparams->par()(momindex+2) ;
  double p2 = px*px+py*py+pz*pz ;
  double mom  = std::sqrt(p2) ;
  
  double dvx = fitparams->par()(posindexmother+0) - fitparams->par()(posindex+0);
  double dvy = fitparams->par()(posindexmother+1) - fitparams->par()(posindex+1);
  double dvz = fitparams->par()(posindexmother+2) - fitparams->par()(posindex+2);
  double len=dvx*dvx+dvy*dvy+dvz*dvz;
  len=std::sqrt(len);
  
  // lineair approximation is fine for now
  for(int row=0; row<3; row++) {
    double posxmother = fitparams->par()(posindexmother+row) ;
    double posx       = fitparams->par()(posindex+row) ;
    double momx       = fitparams->par()(momindex+row) ;
    p.r(row) = posxmother - posx + len*momx/mom ;
    p.H(row,posindexmother+row) = 1 ;
    p.H(row,posindex+row)       = -1 ;
    p.H(row,lenindex)         = momx/mom ;
  }
  // still need these as well
  p.H(0,momindex+0)  = len/mom*( 1 - px*px/p2 ) ;
  p.H(0,momindex+1)  = len/mom*( 0 - px*py/p2 ) ;
  p.H(0,momindex+2)  = len/mom*( 0 - px*pz/p2 ) ;
  p.H(1,momindex+0)  = len/mom*( 0 - py*px/p2 ) ;
  p.H(1,momindex+1)  = len/mom*( 1 - py*py/p2 ) ;
  p.H(1,momindex+2)  = len/mom*( 0 - py*pz/p2 ) ;
  p.H(2,momindex+0)  = len/mom*( 0 - pz*px/p2 ) ;
  p.H(2,momindex+1)  = len/mom*( 0 - pz*py/p2 ) ;
  p.H(2,momindex+2)  = len/mom*( 1 - pz*pz/p2 ) ;
  //if( charge()!=0 ) 
  //{
    //double lambda = bFieldOverC() * charge() ;
    //double tau = fitparams->par(lenIndex());
    //double pt = sqrt(px*px+py*py) ;
    //const double posprecision = 1e-4 ; // 1mu
    //if( fabs(pt*lambda*tau*tau) > posprecision ) 
    //{
      //double sinlt = sin(lambda*tau) ;
      //double coslt = cos(lambda*tau) ;
      //double px1 = px*coslt - py*sinlt ;
      //double py1 = py*coslt + px*sinlt ;
     
      //p.r(0) += -tau*px + (py1-py)/lambda ;
      //p.r(1) += -tau*py - (px1-px)/lambda ;
     
      //p.H(0,lenindex+0) += -px + px1 ;
      //p.H(0,momindex+0) += -tau + sinlt/lambda ;
      //p.H(0,momindex+1) +=        (coslt-1)/lambda ;
      //p.H(1,lenindex+0) += -py + py1 ;
      //p.H(1,momindex+0) +=      - (coslt-1)/lambda ;
      //p.H(1,momindex+1) += -tau + sinlt/lambda ;

      //std::cout << "Using helix for position of particle: " << name().c_str() << " "
           //<< lambda << " " << lambda*tau
           //<< "  delta-x,y: " << -tau*px + (py1-py)/lambda << "  "
           //<< -tau*py - (px1-px)/lambda << std::endl ;
    //}
  //}

  //p.setParticle( *mother() ) ;  // adds geoconstraint chi2 to the balance of the mother particle. Why?
  if(vtxverbose>6){
    std::cout<<"ParticleBase::projectConstraint(): projection is:"<<std::endl;
    std::cout<<"r "; p.r().Print();
    std::cout<<"V "; p.V().Print();
    std::cout<<"H "; RhoCalculationTools::PrintMatrix(p.H());
    }
  return ErrCode::success ;
}

ErrCode DecayTreeFitter::InternalParticle::projectKineConstraint	(	const FitParams *	fitparams,
		Projection &	p
	)		const

Definition at line 302 of file InternalParticle.cxx.

References cos(), energy, fabs(), DecayTreeFitter::Projection::H(), lambda(), DecayTreeFitter::FitParams::par(), DecayTreeFitter::Projection::r(), sin(), sqrt(), DecayTreeFitter::ErrCode::success, and vtxverbose.

{
   if(vtxverbose>6){std::cout<<"InternalParticle::projectKineConstraint():"<<std::endl;}
  // these are in fact four independent constraints. i'll filter
  // them as one, making the code simpler at the expense of a bit of
  // CPU.

  // first add the mother
  int momindex = momIndex() ;
  for(int imom=0; imom<4; imom++) {
    p.r(imom)               = fitparams->par()(momindex+imom) ;
    p.H(imom,momindex+imom) = 1 ;
  }

  // now add the daughters
  for(daucontainer::const_iterator it = daughters().begin() ;
      it != daughters().end() ; ++it) {
    int daulenindex = (*it)->lenIndex() ;
    int daumomindex = (*it)->momIndex() ;
    double mass = (*it)->pdtMass() ;
    double e2=mass*mass ;
    int maxrow = (*it)->hasEnergy() ? 4 : 3 ;
    for(int imom=0; imom<maxrow; ++imom) {
      double px = fitparams->par()(daumomindex+imom) ;
      e2 += px*px ;
      p.r(imom) += -px ;
      p.H(imom,daumomindex+imom) = -1 ;
    }

    if(maxrow==3) {
      // treat the energy for particles that are parameterized with p3
      double energy = sqrt(e2) ;
      p.r(3) += -energy ;
      for(int jmom=0; jmom<3; ++jmom) {
        double px = fitparams->par()(daumomindex+jmom) ;
        p.H(3,daumomindex+jmom) = -px/energy ;
      }
    } else if(false && daulenindex>=0 && (*it)->charge()!=0) {
//FIXME: test and reactivate decay length!
      double tau =  fitparams->par()(daulenindex) ;
      double lambda = bFieldOverC() * (*it)->charge() ;
      double px0 = fitparams->par()(daumomindex) ;
      double py0 = fitparams->par()(daumomindex+1) ;
      double pt0 = sqrt(px0*px0+py0*py0) ;
      const double posprecision = 1e-4 ; // 1mu
      if( fabs(pt0*lambda*tau*tau) > posprecision ) {
        double sinlt = sin(lambda*tau) ;
        double coslt = cos(lambda*tau) ;
        double px = px0*coslt - py0*sinlt ;
        double py = py0*coslt + px0*sinlt ;
        p.r(0) += px0 - px ;
        p.r(1) += py0 - py ;
        p.H(0,daumomindex+0) +=  1 - coslt ;
        p.H(0,daumomindex+1) +=      sinlt ;
        p.H(0,daulenindex+0) +=  lambda*py ;
        p.H(1,daumomindex+0) +=     -sinlt ;
        p.H(1,daumomindex+1) +=  1 - coslt ;
        p.H(1,daulenindex+0) += -lambda*px ;
      }
    }
  }
  return ErrCode::success ;
}

ErrCode DecayTreeFitter::InternalParticle::projectLifeTimeConstraint	(	const FitParams *	,
		Projection &
	)		const

Definition at line 368 of file InternalParticle.cxx.

References DecayTreeFitter::ErrCode::success, and vtxverbose.

{
   if(vtxverbose>6){std::cout<<"InternalParticle::projectLifeTimeConstraint():"<<std::endl;}
  std::cout << "Not yet implemented lifetime constraint!" << std::endl ;
  // int lenindex = lenIndex() ;
  //     assert(lenindex>=0) ;
  //     double tau = pdtTau() ;
  //     p.r(0)            = fitparams->par()(lenindex) - tau ;
  //     p.Vfast(0,0)      = tau*tau ;
  //     p.H(0,lenindex) = 1 ;
  return ErrCode::success ;
}

ErrCode DecayTreeFitter::ParticleBase::projectMassConstraint ( const FitParams *  fitparams, Projection &  p  ) const

virtualinherited

Definition at line 477 of file ParticleBase.cxx.

References DecayTreeFitter::Projection::H(), DecayTreeFitter::FitParams::par(), pz, DecayTreeFitter::Projection::r(), DecayTreeFitter::ErrCode::success, DecayTreeFitter::Projection::Vfast(), and vtxverbose.

{
  if(vtxverbose>6){std::cout<<"ParticleBase::projectMassConstraint():"<<std::endl;}
  double mass  = pdtMass() ;
  double mass2 = mass*mass ;
  int momindex = momIndex() ;

  // initialize the value
  double px = fitparams->par()(momindex+0) ;
  double py = fitparams->par()(momindex+1) ;
  double pz = fitparams->par()(momindex+2) ;
  double E  = fitparams->par()(momindex+3) ;
  if(vtxverbose>6){std::cout<<"px="<<px<<", py="<<py<<", pz="<<pz<<", E="<<E<<", pdtmass="<<mass<<std::endl;}
  p.r(0) = E*E-px*px-py*py-pz*pz-mass2 ;

  // calculate the projection matrix
  p.H(0,momindex+0) = -2.0*px ;
  p.H(0,momindex+1) = -2.0*py ;
  p.H(0,momindex+2) = -2.0*pz ;
  p.H(0,momindex+3) =  2.0*E ;

  // set the variance in the residual
  double width = pdtWidth() ;
  p.Vfast(0,0) = 4*mass*mass*width*width ;

  return ErrCode::success ;
}

ErrCode DecayTreeFitter::InternalParticle::projectMassConstraintTwoBody	(	const FitParams *	fitparams,
		Projection &	p
	)		const

Definition at line 459 of file InternalParticle.cxx.

References DecayTreeFitter::Projection::H(), DecayTreeFitter::ParticleBase::hasEnergy(), m2(), DecayTreeFitter::ParticleBase::momIndex(), DecayTreeFitter::FitParams::par(), DecayTreeFitter::ParticleBase::pdtMass(), DecayTreeFitter::Projection::r(), sqrt(), DecayTreeFitter::ErrCode::success, DecayTreeFitter::Projection::Vfast(), and vtxverbose.

{
   if(vtxverbose>6){std::cout<<"InternalParticle::projectMassConstraintTwoBody():"<<std::endl;}
 // we can also apply the constraint to the daughters. that may
  // work better if the opening angle is small.

  // m^2 = ma^1 + mb^2 + 2 * (Ea*Eb - pxa*pxb - pya*pyb - pza*pzb )

  ParticleBase* d1 = daughters()[0] ;
  ParticleBase* d2 = daughters()[1] ;

  assert(d1->hasEnergy()==false && d2->hasEnergy()==false ) ;

  double mass = pdtMass() ;
  double m1   = d1->pdtMass() ;
  double m2   = d2->pdtMass() ;

  int momindex1 = d1->momIndex() ;
  int momindex2 = d2->momIndex() ;

  // initialize the value
  double px1 = fitparams->par()(momindex1+0) ;
  double py1 = fitparams->par()(momindex1+1) ;
  double pz1 = fitparams->par()(momindex1+2) ;

  double px2 = fitparams->par()(momindex2+0) ;
  double py2 = fitparams->par()(momindex2+1) ;
  double pz2 = fitparams->par()(momindex2+2) ;

  double E1  = std::sqrt( m1*m1 + px1*px1 + py1*py1 + pz1*pz1 ) ;
  double E2  = std::sqrt( m2*m2 + px2*px2 + py2*py2 + pz2*pz2 ) ;

  p.r(0) = m1*m1 + m2*m2 + 2 * (E1*E2 - px1*px2 - py1*py2 - pz1*pz2 ) - mass* mass ;

  // calculate the projection matrix
  p.H(0,momindex1+0) = 2 * (E2 * px1/E1 - px2 ) ;
  p.H(0,momindex1+1) = 2 * (E2 * py1/E1 - py2 ) ;
  p.H(0,momindex1+2) = 2 * (E2 * pz1/E1 - pz2 ) ;
  p.H(0,momindex2+0) = 2 * (E1 * px2/E2 - px1 ) ;
  p.H(0,momindex2+1) = 2 * (E1 * py2/E2 - py1 ) ;
  p.H(0,momindex2+2) = 2 * (E1 * pz2/E2 - pz1 ) ;

  // set the variance in the residual
  double width = pdtWidth() ;
  p.Vfast(0,0) = 4*mass*mass*width*width ;

  return ErrCode::success ;
}

void DecayTreeFitter::ParticleBase::removeDaughter ( const ParticleBase *  pb )

inherited

Definition at line 100 of file ParticleBase.cxx.

Referenced by DecayTreeFitter::Fitter::remove().

{
  daucontainer::iterator it = std::find(m_daughters.begin(),m_daughters.end(),pb) ;
  if(it != m_daughters.end() ) {
    ParticleBase* _dp = *it ;
    m_daughters.erase(it) ;
    if(_dp) delete _dp ;
  } else {
    std::cout << "ERROR: cannot remove particle, because not found ..." << std::endl ;
  }
}

void DecayTreeFitter::ParticleBase::retrieveIndexMap ( indexmap &  anindexmap ) const

virtualinherited

Definition at line 377 of file ParticleBase.cxx.

{

  anindexmap.push_back(std::pair<const ParticleBase*,int>(this,index())) ;
  for(daucontainer::const_iterator it = m_daughters.begin() ;
      it != m_daughters.end() ; ++it)
    (*it)->retrieveIndexMap(anindexmap) ;
}

void DecayTreeFitter::ParticleBase::setIndex ( int  i )

inlineprotectedinherited

Definition at line 141 of file ParticleBase.h.

References i, and DecayTreeFitter::ParticleBase::m_index.

{ m_index = i ; }

bool DecayTreeFitter::ParticleBase::setMassConstraint ( bool  add )

inlineinherited

Definition at line 118 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_hasMassConstraint.

Referenced by DecayTreeFitter::DecayChain::setMassConstraint().

                                     {
      std::swap(add,m_hasMassConstraint) ;
      return add != m_hasMassConstraint ;
    }

void DecayTreeFitter::ParticleBase::setMassConstraint ( double  mass )

inlineinherited

Definition at line 123 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_hasMassConstraint, and DecayTreeFitter::ParticleBase::m_pdtMass.

                                          {
      m_hasMassConstraint = true ;
      m_pdtMass = mass ;
    }

void DecayTreeFitter::ParticleBase::setMother ( const ParticleBase *  m )

inlineinherited

Definition at line 108 of file ParticleBase.h.

References m, and DecayTreeFitter::ParticleBase::m_mother.

{ m_mother = m ; }

void DecayTreeFitter::ParticleBase::setName ( const std::string &  n )

inlineprotectedinherited

Definition at line 142 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_name, and n.

{ m_name = n ; }

void DecayTreeFitter::ParticleBase::setParticle ( RhoCandidate *  bc )

inlineinherited

Definition at line 113 of file ParticleBase.h.

References DecayTreeFitter::ParticleBase::m_particle.

{ m_particle = bc ; }

virtual int DecayTreeFitter::InternalParticle::type ( ) const

inlinevirtual

Implements DecayTreeFitter::ParticleBase.

Reimplemented in DecayTreeFitter::InteractionPoint, and DecayTreeFitter::Resonance.

Definition at line 28 of file InternalParticle.h.

References DecayTreeFitter::ParticleBase::kInternalParticle.

{ return kInternalParticle ; }

void DecayTreeFitter::ParticleBase::updateIndex ( int &  offset )

virtualinherited

Definition at line 83 of file ParticleBase.cxx.

References offset(), and vtxverbose.

Referenced by DecayTreeFitter::Fitter::add(), and DecayTreeFitter::DecayChain::DecayChain().

{
  if(vtxverbose>=5){std::cout<<"ParticleBase::updateIndex() start"<<std::endl;}
  // first the daughters
  for(const_iterator it = begin(); it!= end(); ++it) (*it)->updateIndex(offset) ;
  // now the real work
  m_index = offset ;
  offset += dim() ;
  if(vtxverbose>=5){std::cout<<"ParticleBase::updateIndex() end"<<std::endl;}
}

Member Data Documentation

bool DecayTreeFitter::InternalParticle::m_isconversion

private

Definition at line 54 of file InternalParticle.h.

Referenced by InternalParticle().

bool DecayTreeFitter::InternalParticle::m_lifetimeconstraint

private

Definition at line 53 of file InternalParticle.h.

Referenced by InternalParticle().

---

The documentation for this class was generated from the following files:

• InternalParticle.h
• InternalParticle.cxx