FairRoot/PandaRoot
Public Member Functions | Static Public Member Functions | Protected Member Functions | Static Protected Member Functions | Protected Attributes | List of all members
KFParticleBase Class Referenceabstract

#include <KFParticleBase.h>

Inheritance diagram for KFParticleBase:
KFParticle KFVertex

Public Member Functions

virtual void GetFieldValue (const Double_t xyz[], Double_t B[]) const =0
 
virtual Double_t GetDStoPoint (const Double_t xyz[]) const =0
 
virtual void GetDStoParticle (const KFParticleBase &p, Double_t &DS, Double_t &DSp) const =0
 
virtual void Transport (Double_t dS, Double_t P[], Double_t C[]) const =0
 
 KFParticleBase ()
 
virtual ~KFParticleBase ()
 
void Initialize (const Double_t Param[], const Double_t Cov[], Int_t Charge, Double_t Mass)
 
void Initialize ()
 
void SetVtxGuess (Double_t x, Double_t y, Double_t z)
 
void SetConstructMethod (Int_t m)
 
void SetMassHypo (Double_t m)
 
const Double_tGetMassHypo () const
 
const Double_tGetSumDaughterMass () const
 
Double_t GetX () const
 
Double_t GetY () const
 
Double_t GetZ () const
 
Double_t GetPx () const
 
Double_t GetPy () const
 
Double_t GetPz () const
 
Double_t GetE () const
 
Double_t GetS () const
 
Int_t GetQ () const
 
Double_t GetChi2 () const
 
Int_t GetNDF () const
 
const Double_tX () const
 
const Double_tY () const
 
const Double_tZ () const
 
const Double_tPx () const
 
const Double_tPy () const
 
const Double_tPz () const
 
const Double_tE () const
 
const Double_tS () const
 
const Int_t & Q () const
 
const Double_tChi2 () const
 
const Int_t & NDF () const
 
Double_t GetParameter (Int_t i) const
 
Double_t GetCovariance (Int_t i) const
 
Double_t GetCovariance (Int_t i, Int_t j) const
 
Int_t GetMomentum (Double_t &P, Double_t &SigmaP) const
 
Int_t GetPt (Double_t &Pt, Double_t &SigmaPt) const
 
Int_t GetEta (Double_t &Eta, Double_t &SigmaEta) const
 
Int_t GetPhi (Double_t &Phi, Double_t &SigmaPhi) const
 
Int_t GetMass (Double_t &M, Double_t &SigmaM) const
 
Int_t GetDecayLength (Double_t &L, Double_t &SigmaL) const
 
Int_t GetDecayLengthXY (Double_t &L, Double_t &SigmaL) const
 
Int_t GetLifeTime (Double_t &T, Double_t &SigmaT) const
 
Int_t GetR (Double_t &R, Double_t &SigmaR) const
 
Double_tX ()
 
Double_tY ()
 
Double_tZ ()
 
Double_tPx ()
 
Double_tPy ()
 
Double_tPz ()
 
Double_tE ()
 
Double_tS ()
 
Int_t & Q ()
 
Double_tChi2 ()
 
Int_t & NDF ()
 
Double_tParameter (Int_t i)
 
Double_tCovariance (Int_t i)
 
Double_tCovariance (Int_t i, Int_t j)
 
void operator+= (const KFParticleBase &Daughter)
 
void AddDaughter (const KFParticleBase &Daughter)
 
void AddDaughterWithEnergyFit (const KFParticleBase &Daughter)
 
void AddDaughterWithEnergyCalc (const KFParticleBase &Daughter)
 
void AddDaughterWithEnergyFitMC (const KFParticleBase &Daughter)
 
void SetProductionVertex (const KFParticleBase &Vtx)
 
void SetNonlinearMassConstraint (Double_t Mass)
 
void SetMassConstraint (Double_t Mass, Double_t SigmaMass=0)
 
void SetNoDecayLength ()
 
void Construct (const KFParticleBase *vDaughters[], Int_t nDaughters, const KFParticleBase *ProdVtx=0, Double_t Mass=-1, Bool_t IsConstrained=0)
 
void TransportToDecayVertex ()
 
void TransportToProductionVertex ()
 
void TransportToDS (Double_t dS)
 
Double_t GetDStoPointBz (Double_t Bz, const Double_t xyz[]) const
 
Double_t GetDStoPointBy (Double_t By, const Double_t xyz[]) const
 
void GetDStoParticleBz (Double_t Bz, const KFParticleBase &p, Double_t &dS, Double_t &dS1) const
 
void GetDStoParticleBy (Double_t B, const KFParticleBase &p, Double_t &dS, Double_t &dS1) const
 
Double_t GetDStoPointCBM (const Double_t xyz[]) const
 
void GetDStoParticleCBM (const KFParticleBase &p, Double_t &dS, Double_t &dS1) const
 
void TransportBz (Double_t Bz, Double_t dS, Double_t P[], Double_t C[]) const
 
void TransportCBM (Double_t dS, Double_t P[], Double_t C[]) const
 
Double_t GetDistanceFromVertex (const Double_t vtx[]) const
 
Double_t GetDistanceFromVertex (const KFParticleBase &Vtx) const
 
Double_t GetDistanceFromParticle (const KFParticleBase &p) const
 
Double_t GetDeviationFromVertex (const Double_t v[], const Double_t Cv[]=0) const
 
Double_t GetDeviationFromVertex (const KFParticleBase &Vtx) const
 
Double_t GetDeviationFromParticle (const KFParticleBase &p) const
 
void SubtractFromVertex (KFParticleBase &Vtx) const
 
void ConstructGammaBz (const KFParticleBase &daughter1, const KFParticleBase &daughter2, double Bz)
 
void RotateXY (Double_t angle, Double_t Vtx[3])
 
int Id () const
 
int NDaughters () const
 
const std::vector< int > & DaughterIds () const
 
void SetId (int id)
 
void AddDaughterId (int id)
 
void SetPDG (int pdg)
 
int GetPDG () const
 

Static Public Member Functions

static void GetArmenterosPodolanski (KFParticleBase &positive, KFParticleBase &negative, Double_t QtAlfa[2])
 

Protected Member Functions

Double_tCij (Int_t i, Int_t j)
 
void Convert (bool ToProduction)
 
void TransportLine (Double_t S, Double_t P[], Double_t C[]) const
 
Double_t GetDStoPointLine (const Double_t xyz[]) const
 
void GetDStoParticleLine (const KFParticleBase &p, Double_t &dS, Double_t &dS1) const
 
void GetDSIter (const KFParticleBase &p, Double_t const &dS, Double_t x[3], Double_t dx[3], Double_t ddx[3]) const
 
void GetMeasurement (const Double_t XYZ[], Double_t m[], Double_t V[]) const
 
void SetMassConstraint (Double_t *mP, Double_t *mC, Double_t mJ[7][7], Double_t mass)
 

Static Protected Member Functions

static Int_t IJ (Int_t i, Int_t j)
 
static Bool_t InvertSym3 (const Double_t A[], Double_t Ainv[])
 
static void MultQSQt (const Double_t Q[], const Double_t S[], Double_t SOut[])
 
static Double_t GetSCorrection (const Double_t Part[], const Double_t XYZ[])
 

Protected Attributes

Double_t fP [8]
 
Double_t fC [36]
 
Int_t fQ
 
Int_t fNDF
 
Double_t fChi2
 
Double_t fSFromDecay
 
Bool_t fAtProductionVertex
 
Double_t fVtxGuess [3]
 
Bool_t fIsLinearized
 
Int_t fConstructMethod
 
Double_t SumDaughterMass
 
Double_t fMassHypo
 
int fId
 
std::vector< int > fDaughtersIds
 
int fPDG
 

Detailed Description

Definition at line 27 of file KFParticleBase.h.

Constructor & Destructor Documentation

KFParticleBase::KFParticleBase ( )

Definition at line 24 of file KFParticleBase.cxx.

References Initialize().

24  :fQ(0), fNDF(-3), fChi2(0), fSFromDecay(0), fAtProductionVertex(0), fIsLinearized(0),
26 {
27  //* Constructor
28 
29  Initialize();
30 }
Double_t SumDaughterMass
Double_t fSFromDecay
Bool_t fAtProductionVertex
Double_t fMassHypo
std::vector< int > fDaughtersIds
virtual KFParticleBase::~KFParticleBase ( )
inlinevirtual

Definition at line 69 of file KFParticleBase.h.

69 { ; }

Member Function Documentation

void KFParticleBase::AddDaughter ( const KFParticleBase Daughter)

Definition at line 395 of file KFParticleBase.cxx.

References AddDaughterWithEnergyCalc(), AddDaughterWithEnergyFit(), AddDaughterWithEnergyFitMC(), fC, fConstructMethod, fMassHypo, fNDF, fP, fQ, fSFromDecay, GetQ(), i, and SumDaughterMass.

Referenced by KFParticle::AddDaughter(), Construct(), and operator+=().

396 {
397  if( fNDF<-1 ){ // first daughter -> just copy
398  fNDF = -1;
399  fQ = Daughter.GetQ();
400  for( Int_t i=0; i<7; i++) fP[i] = Daughter.fP[i];
401  for( Int_t i=0; i<28; i++) fC[i] = Daughter.fC[i];
402  fSFromDecay = 0;
403  fMassHypo = Daughter.fMassHypo;
404  SumDaughterMass = Daughter.SumDaughterMass;
405  return;
406  }
407 
408  if(fConstructMethod == 0)
409  AddDaughterWithEnergyFit(Daughter);
410  else if(fConstructMethod == 1)
411  AddDaughterWithEnergyCalc(Daughter);
412  else if(fConstructMethod == 2)
413  AddDaughterWithEnergyFitMC(Daughter);
414 
415  SumDaughterMass += Daughter.SumDaughterMass;
416  fMassHypo = -1;
417 }
Int_t i
Definition: run_full.C:25
Double_t fP[8]
void AddDaughterWithEnergyCalc(const KFParticleBase &Daughter)
Double_t SumDaughterMass
Double_t fSFromDecay
Double_t fC[36]
void AddDaughterWithEnergyFit(const KFParticleBase &Daughter)
void AddDaughterWithEnergyFitMC(const KFParticleBase &Daughter)
Double_t fMassHypo
Int_t GetQ() const
void KFParticleBase::AddDaughterId ( int  id)
inline

Definition at line 274 of file KFParticleBase.h.

References fDaughtersIds.

Referenced by KFParticleFinder::FindParticles().

274 { fDaughtersIds.push_back(id); };
std::vector< int > fDaughtersIds
void KFParticleBase::AddDaughterWithEnergyCalc ( const KFParticleBase Daughter)

Definition at line 576 of file KFParticleBase.cxx.

References b, Double_t, fabs(), fC, fChi2, fIsLinearized, fNDF, fP, fQ, fSFromDecay, fVtxGuess, GetDStoParticle(), GetFieldValue(), GetMeasurement(), GetQ(), i, m, s, sqrt(), Transport(), TransportToDecayVertex(), and TransportToDS().

Referenced by AddDaughter().

577 {
578  //* Energy considered as a dependent variable, calculated from the momentum and mass hypothesis
579 
580  //* Add daughter
581 
583 
584  Double_t b[3];
585  Int_t maxIter = 1;
586 
587  if( !fIsLinearized ){
588  if( fNDF==-1 ){
589  Double_t ds, ds1;
590  GetDStoParticle(Daughter, ds, ds1);
591  TransportToDS( ds );
592  Double_t m[8];
593  Double_t mCd[36];
594  Daughter.Transport( ds1, m, mCd );
595  fVtxGuess[0] = .5*( fP[0] + m[0] );
596  fVtxGuess[1] = .5*( fP[1] + m[1] );
597  fVtxGuess[2] = .5*( fP[2] + m[2] );
598  } else {
599  fVtxGuess[0] = fP[0];
600  fVtxGuess[1] = fP[1];
601  fVtxGuess[2] = fP[2];
602  }
603  maxIter = 3;
604  }
605 
606  for( Int_t iter=0; iter<maxIter; iter++ ){
607 
608  {
609  GetFieldValue( fVtxGuess, b );
610  const Double_t kCLight = 0.000299792458;
611  b[0]*=kCLight; b[1]*=kCLight; b[2]*=kCLight;
612  }
613 
614  Double_t *ffP = fP, *ffC = fC, tmpP[8], tmpC[36];
615  if( fNDF==-1 ){
616  GetMeasurement( fVtxGuess, tmpP, tmpC );
617  ffP = tmpP;
618  ffC = tmpC;
619  }
620 
621  Double_t m[8], mV[36];
622 
623  if( Daughter.fC[35]>0 ){
624  Daughter.GetMeasurement( fVtxGuess, m, mV );
625  } else {
626  for( Int_t i=0; i<8; i++ ) m[i] = Daughter.fP[i];
627  for( Int_t i=0; i<36; i++ ) mV[i] = Daughter.fC[i];
628  }
629 
630  double massMf2 = m[6]*m[6] - (m[3]*m[3] + m[4]*m[4] + m[5]*m[5]);
631  double massRf2 = fP[6]*fP[6] - (fP[3]*fP[3] + fP[4]*fP[4] + fP[5]*fP[5]);
632 
633  //*
634 
635  Double_t mS[6];
636  {
637  Double_t mSi[6] = { ffC[0]+mV[0],
638  ffC[1]+mV[1], ffC[2]+mV[2],
639  ffC[3]+mV[3], ffC[4]+mV[4], ffC[5]+mV[5] };
640 
641  mS[0] = mSi[2]*mSi[5] - mSi[4]*mSi[4];
642  mS[1] = mSi[3]*mSi[4] - mSi[1]*mSi[5];
643  mS[2] = mSi[0]*mSi[5] - mSi[3]*mSi[3];
644  mS[3] = mSi[1]*mSi[4] - mSi[2]*mSi[3];
645  mS[4] = mSi[1]*mSi[3] - mSi[0]*mSi[4];
646  mS[5] = mSi[0]*mSi[2] - mSi[1]*mSi[1];
647 
648  Double_t s = ( mSi[0]*mS[0] + mSi[1]*mS[1] + mSi[3]*mS[3] );
649 
650  s = ( s > 1.E-20 ) ?1./s :0;
651  mS[0]*=s;
652  mS[1]*=s;
653  mS[2]*=s;
654  mS[3]*=s;
655  mS[4]*=s;
656  mS[5]*=s;
657  }
658 
659  //* Residual (measured - estimated)
660 
661  Double_t zeta[3] = { m[0]-ffP[0], m[1]-ffP[1], m[2]-ffP[2] };
662 
663  //* CHt = CH' - D'
664 
665  Double_t mCHt0[6], mCHt1[6], mCHt2[6];
666 
667  mCHt0[0]=ffC[ 0] ; mCHt1[0]=ffC[ 1] ; mCHt2[0]=ffC[ 3] ;
668  mCHt0[1]=ffC[ 1] ; mCHt1[1]=ffC[ 2] ; mCHt2[1]=ffC[ 4] ;
669  mCHt0[2]=ffC[ 3] ; mCHt1[2]=ffC[ 4] ; mCHt2[2]=ffC[ 5] ;
670  mCHt0[3]=ffC[ 6]-mV[ 6]; mCHt1[3]=ffC[ 7]-mV[ 7]; mCHt2[3]=ffC[ 8]-mV[ 8];
671  mCHt0[4]=ffC[10]-mV[10]; mCHt1[4]=ffC[11]-mV[11]; mCHt2[4]=ffC[12]-mV[12];
672  mCHt0[5]=ffC[15]-mV[15]; mCHt1[5]=ffC[16]-mV[16]; mCHt2[5]=ffC[17]-mV[17];
673 
674  //* Kalman gain K = mCH'*S
675 
676  Double_t k0[6], k1[6], k2[6];
677 
678  for(Int_t i=0;i<6;++i){
679  k0[i] = mCHt0[i]*mS[0] + mCHt1[i]*mS[1] + mCHt2[i]*mS[3];
680  k1[i] = mCHt0[i]*mS[1] + mCHt1[i]*mS[2] + mCHt2[i]*mS[4];
681  k2[i] = mCHt0[i]*mS[3] + mCHt1[i]*mS[4] + mCHt2[i]*mS[5];
682  }
683 
684  //* New estimation of the vertex position
685 
686  if( iter<maxIter-1 ){
687  for(Int_t i=0; i<3; ++i)
688  fVtxGuess[i]= ffP[i] + k0[i]*zeta[0]+k1[i]*zeta[1]+k2[i]*zeta[2];
689  continue;
690  }
691 
692  //* find mf and mVf - optimum value of the measurement and its covariance matrix
693  //* mVHt = V*H'
694  Double_t mVHt0[6], mVHt1[6], mVHt2[6];
695 
696  mVHt0[0]= mV[ 0] ; mVHt1[0]= mV[ 1] ; mVHt2[0]= mV[ 3] ;
697  mVHt0[1]= mV[ 1] ; mVHt1[1]= mV[ 2] ; mVHt2[1]= mV[ 4] ;
698  mVHt0[2]= mV[ 3] ; mVHt1[2]= mV[ 4] ; mVHt2[2]= mV[ 5] ;
699  mVHt0[3]= mV[ 6] ; mVHt1[3]= mV[ 7] ; mVHt2[3]= mV[ 8] ;
700  mVHt0[4]= mV[10] ; mVHt1[4]= mV[11] ; mVHt2[4]= mV[12] ;
701  mVHt0[5]= mV[15] ; mVHt1[5]= mV[16] ; mVHt2[5]= mV[17] ;
702 
703  //* Kalman gain Km = mCH'*S
704 
705  Double_t km0[6], km1[6], km2[6];
706 
707  for(Int_t i=0;i<6;++i){
708  km0[i] = mVHt0[i]*mS[0] + mVHt1[i]*mS[1] + mVHt2[i]*mS[3];
709  km1[i] = mVHt0[i]*mS[1] + mVHt1[i]*mS[2] + mVHt2[i]*mS[4];
710  km2[i] = mVHt0[i]*mS[3] + mVHt1[i]*mS[4] + mVHt2[i]*mS[5];
711  }
712 
713  Double_t mf[7] = { m[0], m[1], m[2], m[3], m[4], m[5], m[6] };
714 
715  for(Int_t i=0;i<6;++i)
716  mf[i] = mf[i] - km0[i]*zeta[0] - km1[i]*zeta[1] - km2[i]*zeta[2];
717 
718  Double_t energyMf = sqrt( massMf2 + (mf[3]*mf[3] + mf[4]*mf[4] + mf[5]*mf[5]) );
719 
720  Double_t mVf[28];
721  for(Int_t iC=0; iC<28; iC++)
722  mVf[iC] = mV[iC];
723 
724  //* hmf = d(energyMf)/d(mf)
725  Double_t hmf[7];
726  if( fabs(energyMf) < 1.e-10) hmf[3] = 0; else hmf[3] = mf[3]/energyMf;
727  if( fabs(energyMf) < 1.e-10) hmf[4] = 0; else hmf[4] = mf[4]/energyMf;
728  if( fabs(energyMf) < 1.e-10) hmf[5] = 0; else hmf[5] = mf[5]/energyMf;
729 // if( fabs(energyMf) < 1.e-10) hmf[6] = 0; else hmf[6] = mf[6]/energyMf;
730  hmf[6] = 0;
731 
732  for(Int_t i=0, k=0;i<6;++i){
733  for(Int_t j=0;j<=i;++j,++k){
734  mVf[k] = mVf[k] - (km0[i]*mVHt0[j] + km1[i]*mVHt1[j] + km2[i]*mVHt2[j] );
735  }
736  }
737  Double_t mVf24 = mVf[24], mVf25 = mVf[25], mVf26 = mVf[26];
738  mVf[21] = mVf[6 ]*hmf[3] + mVf[10]*hmf[4] + mVf[15]*hmf[5] + mVf[21]*hmf[6];
739  mVf[22] = mVf[7 ]*hmf[3] + mVf[11]*hmf[4] + mVf[16]*hmf[5] + mVf[22]*hmf[6];
740  mVf[23] = mVf[8 ]*hmf[3] + mVf[12]*hmf[4] + mVf[17]*hmf[5] + mVf[23]*hmf[6];
741  mVf[24] = mVf[9 ]*hmf[3] + mVf[13]*hmf[4] + mVf[18]*hmf[5] + mVf[24]*hmf[6];
742  mVf[25] = mVf[13]*hmf[3] + mVf[14]*hmf[4] + mVf[19]*hmf[5] + mVf[25]*hmf[6];
743  mVf[26] = mVf[18]*hmf[3] + mVf[19]*hmf[4] + mVf[20]*hmf[5] + mVf[26]*hmf[6];
744  mVf[27] = mVf[24]*hmf[3] + mVf[25]*hmf[4] + mVf[26]*hmf[5] + (mVf24*hmf[3] + mVf25*hmf[4] + mVf26*hmf[5] + mVf[27]*hmf[6])*hmf[6]; //here mVf[] are already modified
745 
746  mf[6] = energyMf;
747 
748  //* find rf and mCf - optimum value of the measurement and its covariance matrix
749 
750  //* mCCHt = C*H'
751  Double_t mCCHt0[6], mCCHt1[6], mCCHt2[6];
752 
753  mCCHt0[0]=ffC[ 0]; mCCHt1[0]=ffC[ 1]; mCCHt2[0]=ffC[ 3];
754  mCCHt0[1]=ffC[ 1]; mCCHt1[1]=ffC[ 2]; mCCHt2[1]=ffC[ 4];
755  mCCHt0[2]=ffC[ 3]; mCCHt1[2]=ffC[ 4]; mCCHt2[2]=ffC[ 5];
756  mCCHt0[3]=ffC[ 6]; mCCHt1[3]=ffC[ 7]; mCCHt2[3]=ffC[ 8];
757  mCCHt0[4]=ffC[10]; mCCHt1[4]=ffC[11]; mCCHt2[4]=ffC[12];
758  mCCHt0[5]=ffC[15]; mCCHt1[5]=ffC[16]; mCCHt2[5]=ffC[17];
759 
760  //* Kalman gain Krf = mCH'*S
761 
762  Double_t krf0[6], krf1[6], krf2[6];
763 
764  for(Int_t i=0;i<6;++i){
765  krf0[i] = mCCHt0[i]*mS[0] + mCCHt1[i]*mS[1] + mCCHt2[i]*mS[3];
766  krf1[i] = mCCHt0[i]*mS[1] + mCCHt1[i]*mS[2] + mCCHt2[i]*mS[4];
767  krf2[i] = mCCHt0[i]*mS[3] + mCCHt1[i]*mS[4] + mCCHt2[i]*mS[5];
768  }
769  Double_t rf[7] = { ffP[0], ffP[1], ffP[2], ffP[3], ffP[4], ffP[5], ffP[6] };
770 
771  for(Int_t i=0;i<6;++i)
772  rf[i] = rf[i] + krf0[i]*zeta[0] + krf1[i]*zeta[1] + krf2[i]*zeta[2];
773 
774  Double_t energyRf = sqrt( massRf2 + (rf[3]*rf[3] + rf[4]*rf[4] + rf[5]*rf[5]) );
775 
776  Double_t mCf[28];
777  for(Int_t iC=0; iC<28; iC++)
778  mCf[iC] = ffC[iC];
779  //* hrf = d(Erf)/d(rf)
780  Double_t hrf[7];
781  if( fabs(energyRf) < 1.e-10) hrf[3] = 0; else hrf[3] = rf[3]/energyRf;
782  if( fabs(energyRf) < 1.e-10) hrf[4] = 0; else hrf[4] = rf[4]/energyRf;
783  if( fabs(energyRf) < 1.e-10) hrf[5] = 0; else hrf[5] = rf[5]/energyRf;
784 // if( fabs(energyRf) < 1.e-10) hrf[6] = 0; else hrf[6] = rf[6]/energyRf;
785  hrf[6] = 0;
786 
787  for(Int_t i=0, k=0;i<6;++i){
788  for(Int_t j=0;j<=i;++j,++k){
789  mCf[k] = mCf[k] - (krf0[i]*mCCHt0[j] + krf1[i]*mCCHt1[j] + krf2[i]*mCCHt2[j] );
790  }
791  }
792  Double_t mCf24 = mCf[24], mCf25 = mCf[25], mCf26 = mCf[26];
793  mCf[21] = mCf[6 ]*hrf[3] + mCf[10]*hrf[4] + mCf[15]*hrf[5] + mCf[21]*hrf[6];
794  mCf[22] = mCf[7 ]*hrf[3] + mCf[11]*hrf[4] + mCf[16]*hrf[5] + mCf[22]*hrf[6];
795  mCf[23] = mCf[8 ]*hrf[3] + mCf[12]*hrf[4] + mCf[17]*hrf[5] + mCf[23]*hrf[6];
796  mCf[24] = mCf[9 ]*hrf[3] + mCf[13]*hrf[4] + mCf[18]*hrf[5] + mCf[24]*hrf[6];
797  mCf[25] = mCf[13]*hrf[3] + mCf[14]*hrf[4] + mCf[19]*hrf[5] + mCf[25]*hrf[6];
798  mCf[26] = mCf[18]*hrf[3] + mCf[19]*hrf[4] + mCf[20]*hrf[5] + mCf[26]*hrf[6];
799  mCf[27] = mCf[24]*hrf[3] + mCf[25]*hrf[4] + mCf[26]*hrf[5] + (mCf24*hrf[3] + mCf25*hrf[4] + mCf26*hrf[5] + mCf[27]*hrf[6])*hrf[6]; //here mCf[] are already modified
800 
801  for(Int_t iC=21; iC<28; iC++)
802  {
803  ffC[iC] = mCf[iC];
804  mV[iC] = mVf[iC];
805  }
806 
807  fP[6] = energyRf + energyMf;
808  rf[6] = energyRf;
809 
810  //Double_t Dvv[3][3]; do not need this
811  Double_t mDvp[3][3];
812  Double_t mDpv[3][3];
813  Double_t mDpp[3][3];
814  Double_t mDe[7];
815 
816  for(int i=0; i<3; i++)
817  {
818  for(int j=0; j<3; j++)
819  {
820  mDvp[i][j] = km0[i+3]*mCCHt0[j] + km1[i+3]*mCCHt1[j] + km2[i+3]*mCCHt2[j];
821  mDpv[i][j] = km0[i]*mCCHt0[j+3] + km1[i]*mCCHt1[j+3] + km2[i]*mCCHt2[j+3];
822  mDpp[i][j] = km0[i+3]*mCCHt0[j+3] + km1[i+3]*mCCHt1[j+3] + km2[i+3]*mCCHt2[j+3];
823  }
824  }
825 
826  mDe[0] = hmf[3]*mDvp[0][0] + hmf[4]*mDvp[1][0] + hmf[5]*mDvp[2][0];
827  mDe[1] = hmf[3]*mDvp[0][1] + hmf[4]*mDvp[1][1] + hmf[5]*mDvp[2][1];
828  mDe[2] = hmf[3]*mDvp[0][2] + hmf[4]*mDvp[1][2] + hmf[5]*mDvp[2][2];
829  mDe[3] = hmf[3]*mDpp[0][0] + hmf[4]*mDpp[1][0] + hmf[5]*mDpp[2][0];
830  mDe[4] = hmf[3]*mDpp[0][1] + hmf[4]*mDpp[1][1] + hmf[5]*mDpp[2][1];
831  mDe[5] = hmf[3]*mDpp[0][2] + hmf[4]*mDpp[1][2] + hmf[5]*mDpp[2][2];
832  mDe[6] = 2*(mDe[3]*hrf[3] + mDe[4]*hrf[4] + mDe[5]*hrf[5]);
833 
834  // last itearation -> update the particle
835 
836  //* Add the daughter momentum to the particle momentum
837 
838  ffP[ 3] += m[ 3];
839  ffP[ 4] += m[ 4];
840  ffP[ 5] += m[ 5];
841 
842  ffC[ 9] += mV[ 9];
843  ffC[13] += mV[13];
844  ffC[14] += mV[14];
845  ffC[18] += mV[18];
846  ffC[19] += mV[19];
847  ffC[20] += mV[20];
848  ffC[24] += mV[24];
849  ffC[25] += mV[25];
850  ffC[26] += mV[26];
851  ffC[27] += mV[27];
852 
853  ffC[21] += mDe[0];
854  ffC[22] += mDe[1];
855  ffC[23] += mDe[2];
856  ffC[24] += mDe[3];
857  ffC[25] += mDe[4];
858  ffC[26] += mDe[5];
859  ffC[27] += mDe[6];
860 
861  //* New estimation of the vertex position r += K*zeta
862 
863  for(Int_t i=0;i<6;++i)
864  fP[i] = ffP[i] + k0[i]*zeta[0] + k1[i]*zeta[1] + k2[i]*zeta[2];
865 
866  //* New covariance matrix C -= K*(mCH')'
867 
868  for(Int_t i=0, k=0;i<6;++i){
869  for(Int_t j=0;j<=i;++j,++k){
870  fC[k] = ffC[k] - (k0[i]*mCHt0[j] + k1[i]*mCHt1[j] + k2[i]*mCHt2[j] );
871  }
872  }
873 
874  for(int i=21; i<28; i++) fC[i] = ffC[i];
875 
876  //* Calculate Chi^2
877 
878  fNDF += 2;
879  fQ += Daughter.GetQ();
880  fSFromDecay = 0;
881  fChi2 += (mS[0]*zeta[0] + mS[1]*zeta[1] + mS[3]*zeta[2])*zeta[0]
882  + (mS[1]*zeta[0] + mS[2]*zeta[1] + mS[4]*zeta[2])*zeta[1]
883  + (mS[3]*zeta[0] + mS[4]*zeta[1] + mS[5]*zeta[2])*zeta[2];
884  }
885 }
void TransportToDecayVertex()
Int_t i
Definition: run_full.C:25
__m128 m
Definition: P4_F32vec4.h:28
TTree * b
Double_t fVtxGuess[3]
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
TLorentzVector s
Definition: Pnd2DStar.C:50
void GetMeasurement(const Double_t XYZ[], Double_t m[], Double_t V[]) const
virtual void GetDStoParticle(const KFParticleBase &p, Double_t &DS, Double_t &DSp) const =0
virtual void Transport(Double_t dS, Double_t P[], Double_t C[]) const =0
virtual void GetFieldValue(const Double_t xyz[], Double_t B[]) const =0
Double_t
Double_t fSFromDecay
Double_t fC[36]
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
void TransportToDS(Double_t dS)
Int_t GetQ() const
void KFParticleBase::AddDaughterWithEnergyFit ( const KFParticleBase Daughter)

Definition at line 419 of file KFParticleBase.cxx.

References b, Double_t, fabs(), fC, fChi2, fIsLinearized, fNDF, fP, fQ, fSFromDecay, fVtxGuess, GetDStoParticle(), GetFieldValue(), GetMeasurement(), GetQ(), i, m, s, Transport(), TransportToDecayVertex(), and TransportToDS().

Referenced by AddDaughter().

420 {
421  //* Energy considered as an independent veriable, fitted independently from momentum, without any constraints on mass
422 
423  //* Add daughter
424 
426 
427  Double_t b[3];
428  Int_t maxIter = 1;
429 
430  if( !fIsLinearized ){
431  if( fNDF==-1 ){
432  Double_t ds, ds1;
433  GetDStoParticle(Daughter, ds, ds1);
434  TransportToDS( ds );
435  Double_t m[8];
436  Double_t mCd[36];
437  Daughter.Transport( ds1, m, mCd );
438  fVtxGuess[0] = .5*( fP[0] + m[0] );
439  fVtxGuess[1] = .5*( fP[1] + m[1] );
440  fVtxGuess[2] = .5*( fP[2] + m[2] );
441  } else {
442  fVtxGuess[0] = fP[0];
443  fVtxGuess[1] = fP[1];
444  fVtxGuess[2] = fP[2];
445  }
446  maxIter = 3;
447  }
448 
449  for( Int_t iter=0; iter<maxIter; iter++ ){
450 
451  {
452  GetFieldValue( fVtxGuess, b );
453  const Double_t kCLight = 0.000299792458;
454  b[0]*=kCLight; b[1]*=kCLight; b[2]*=kCLight;
455  }
456 
457  Double_t *ffP = fP, *ffC = fC, tmpP[8], tmpC[36];
458  if( fNDF==-1 ){
459  GetMeasurement( fVtxGuess, tmpP, tmpC );
460  ffP = tmpP;
461  ffC = tmpC;
462  }
463 
464  Double_t m[8], mV[36];
465 
466  if( Daughter.fC[35]>0 ){
467  Daughter.GetMeasurement( fVtxGuess, m, mV );
468  } else {
469  for( Int_t i=0; i<8; i++ ) m[i] = Daughter.fP[i];
470  for( Int_t i=0; i<36; i++ ) mV[i] = Daughter.fC[i];
471  }
472  //*
473 
474  Double_t mS[6];
475  {
476  Double_t mSi[6] = { ffC[0]+mV[0],
477  ffC[1]+mV[1], ffC[2]+mV[2],
478  ffC[3]+mV[3], ffC[4]+mV[4], ffC[5]+mV[5] };
479 
480  mS[0] = mSi[2]*mSi[5] - mSi[4]*mSi[4];
481  mS[1] = mSi[3]*mSi[4] - mSi[1]*mSi[5];
482  mS[2] = mSi[0]*mSi[5] - mSi[3]*mSi[3];
483  mS[3] = mSi[1]*mSi[4] - mSi[2]*mSi[3];
484  mS[4] = mSi[1]*mSi[3] - mSi[0]*mSi[4];
485  mS[5] = mSi[0]*mSi[2] - mSi[1]*mSi[1];
486 
487  Double_t s = ( mSi[0]*mS[0] + mSi[1]*mS[1] + mSi[3]*mS[3] );
488  s = ( fabs(s) > 1.E-20 ) ?1./s :0;
489  mS[0]*=s;
490  mS[1]*=s;
491  mS[2]*=s;
492  mS[3]*=s;
493  mS[4]*=s;
494  mS[5]*=s;
495  }
496  //* Residual (measured - estimated)
497 
498  Double_t zeta[3] = { m[0]-ffP[0], m[1]-ffP[1], m[2]-ffP[2] };
499 
500  //* CHt = CH' - D'
501 
502  Double_t mCHt0[7], mCHt1[7], mCHt2[7];
503 
504  mCHt0[0]=ffC[ 0] ; mCHt1[0]=ffC[ 1] ; mCHt2[0]=ffC[ 3] ;
505  mCHt0[1]=ffC[ 1] ; mCHt1[1]=ffC[ 2] ; mCHt2[1]=ffC[ 4] ;
506  mCHt0[2]=ffC[ 3] ; mCHt1[2]=ffC[ 4] ; mCHt2[2]=ffC[ 5] ;
507  mCHt0[3]=ffC[ 6]-mV[ 6]; mCHt1[3]=ffC[ 7]-mV[ 7]; mCHt2[3]=ffC[ 8]-mV[ 8];
508  mCHt0[4]=ffC[10]-mV[10]; mCHt1[4]=ffC[11]-mV[11]; mCHt2[4]=ffC[12]-mV[12];
509  mCHt0[5]=ffC[15]-mV[15]; mCHt1[5]=ffC[16]-mV[16]; mCHt2[5]=ffC[17]-mV[17];
510  mCHt0[6]=ffC[21]-mV[21]; mCHt1[6]=ffC[22]-mV[22]; mCHt2[6]=ffC[23]-mV[23];
511 
512  //* Kalman gain K = mCH'*S
513 
514  Double_t k0[7], k1[7], k2[7];
515 
516  for(Int_t i=0;i<7;++i){
517  k0[i] = mCHt0[i]*mS[0] + mCHt1[i]*mS[1] + mCHt2[i]*mS[3];
518  k1[i] = mCHt0[i]*mS[1] + mCHt1[i]*mS[2] + mCHt2[i]*mS[4];
519  k2[i] = mCHt0[i]*mS[3] + mCHt1[i]*mS[4] + mCHt2[i]*mS[5];
520  }
521 
522  //* New estimation of the vertex position
523 
524  if( iter<maxIter-1 ){
525  for(Int_t i=0; i<3; ++i)
526  fVtxGuess[i]= ffP[i] + k0[i]*zeta[0]+k1[i]*zeta[1]+k2[i]*zeta[2];
527  continue;
528  }
529 
530  // last itearation -> update the particle
531 
532  //* Add the daughter momentum to the particle momentum
533 
534  ffP[ 3] += m[ 3];
535  ffP[ 4] += m[ 4];
536  ffP[ 5] += m[ 5];
537  ffP[ 6] += m[ 6];
538 
539  ffC[ 9] += mV[ 9];
540  ffC[13] += mV[13];
541  ffC[14] += mV[14];
542  ffC[18] += mV[18];
543  ffC[19] += mV[19];
544  ffC[20] += mV[20];
545  ffC[24] += mV[24];
546  ffC[25] += mV[25];
547  ffC[26] += mV[26];
548  ffC[27] += mV[27];
549 
550 
551  //* New estimation of the vertex position r += K*zeta
552 
553  for(Int_t i=0;i<7;++i)
554  fP[i] = ffP[i] + k0[i]*zeta[0] + k1[i]*zeta[1] + k2[i]*zeta[2];
555 
556  //* New covariance matrix C -= K*(mCH')'
557 
558  for(Int_t i=0, k=0;i<7;++i){
559  for(Int_t j=0;j<=i;++j,++k){
560  fC[k] = ffC[k] - (k0[i]*mCHt0[j] + k1[i]*mCHt1[j] + k2[i]*mCHt2[j] );
561  }
562  }
563 
564  //* Calculate Chi^2
565 
566  fNDF += 2;
567  fQ += Daughter.GetQ();
568  fSFromDecay = 0;
569  fChi2 += (mS[0]*zeta[0] + mS[1]*zeta[1] + mS[3]*zeta[2])*zeta[0]
570  + (mS[1]*zeta[0] + mS[2]*zeta[1] + mS[4]*zeta[2])*zeta[1]
571  + (mS[3]*zeta[0] + mS[4]*zeta[1] + mS[5]*zeta[2])*zeta[2];
572 
573  }
574 }
void TransportToDecayVertex()
Int_t i
Definition: run_full.C:25
__m128 m
Definition: P4_F32vec4.h:28
TTree * b
Double_t fVtxGuess[3]
Double_t fP[8]
TLorentzVector s
Definition: Pnd2DStar.C:50
void GetMeasurement(const Double_t XYZ[], Double_t m[], Double_t V[]) const
virtual void GetDStoParticle(const KFParticleBase &p, Double_t &DS, Double_t &DSp) const =0
virtual void Transport(Double_t dS, Double_t P[], Double_t C[]) const =0
virtual void GetFieldValue(const Double_t xyz[], Double_t B[]) const =0
Double_t
Double_t fSFromDecay
Double_t fC[36]
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
void TransportToDS(Double_t dS)
Int_t GetQ() const
void KFParticleBase::AddDaughterWithEnergyFitMC ( const KFParticleBase Daughter)

Definition at line 887 of file KFParticleBase.cxx.

References b, Double_t, fC, fChi2, fIsLinearized, fMassHypo, fNDF, fP, fQ, fSFromDecay, fVtxGuess, GetDStoParticle(), GetFieldValue(), GetMeasurement(), GetQ(), i, m, s, SetMassConstraint(), sqrt(), SumDaughterMass, Transport(), TransportToDecayVertex(), and TransportToDS().

Referenced by AddDaughter().

888 {
889  //* Energy considered as an independent variable, fitted independently from momentum, without any constraints on mass
890 
891  //* Add daughter
892 
894 
895  Double_t b[3];
896  Int_t maxIter = 1;
897 
898  if( !fIsLinearized ){
899  if( fNDF==-1 ){
900  Double_t ds, ds1;
901  GetDStoParticle(Daughter, ds, ds1);
902  TransportToDS( ds );
903  Double_t m[8];
904  Double_t mCd[36];
905  Daughter.Transport( ds1, m, mCd );
906  fVtxGuess[0] = .5*( fP[0] + m[0] );
907  fVtxGuess[1] = .5*( fP[1] + m[1] );
908  fVtxGuess[2] = .5*( fP[2] + m[2] );
909  } else {
910  fVtxGuess[0] = fP[0];
911  fVtxGuess[1] = fP[1];
912  fVtxGuess[2] = fP[2];
913  }
914  maxIter = 3;
915  }
916 
917  for( Int_t iter=0; iter<maxIter; iter++ ){
918 
919  {
920  GetFieldValue( fVtxGuess, b );
921  const Double_t kCLight = 0.000299792458;
922  b[0]*=kCLight; b[1]*=kCLight; b[2]*=kCLight;
923  }
924 
925  Double_t *ffP = fP, *ffC = fC, tmpP[8], tmpC[36];
926  if( fNDF==-1 ){
927  GetMeasurement( fVtxGuess, tmpP, tmpC );
928  ffP = tmpP;
929  ffC = tmpC;
930  }
931  Double_t m[8], mV[36];
932 
933  if( Daughter.fC[35]>0 ){
934  Daughter.GetMeasurement( fVtxGuess, m, mV );
935  } else {
936  for( Int_t i=0; i<8; i++ ) m[i] = Daughter.fP[i];
937  for( Int_t i=0; i<36; i++ ) mV[i] = Daughter.fC[i];
938  }
939  //*
940 
941  Double_t mS[6];
942  {
943  Double_t mSi[6] = { ffC[0]+mV[0],
944  ffC[1]+mV[1], ffC[2]+mV[2],
945  ffC[3]+mV[3], ffC[4]+mV[4], ffC[5]+mV[5] };
946 
947  mS[0] = mSi[2]*mSi[5] - mSi[4]*mSi[4];
948  mS[1] = mSi[3]*mSi[4] - mSi[1]*mSi[5];
949  mS[2] = mSi[0]*mSi[5] - mSi[3]*mSi[3];
950  mS[3] = mSi[1]*mSi[4] - mSi[2]*mSi[3];
951  mS[4] = mSi[1]*mSi[3] - mSi[0]*mSi[4];
952  mS[5] = mSi[0]*mSi[2] - mSi[1]*mSi[1];
953 
954  Double_t s = ( mSi[0]*mS[0] + mSi[1]*mS[1] + mSi[3]*mS[3] );
955 
956  s = ( s > 1.E-20 ) ?1./s :0;
957  mS[0]*=s;
958  mS[1]*=s;
959  mS[2]*=s;
960  mS[3]*=s;
961  mS[4]*=s;
962  mS[5]*=s;
963  }
964  //* Residual (measured - estimated)
965 
966  Double_t zeta[3] = { m[0]-ffP[0], m[1]-ffP[1], m[2]-ffP[2] };
967 
968 
969  //* CHt = CH'
970 
971  Double_t mCHt0[7], mCHt1[7], mCHt2[7];
972 
973  mCHt0[0]=ffC[ 0] ; mCHt1[0]=ffC[ 1] ; mCHt2[0]=ffC[ 3] ;
974  mCHt0[1]=ffC[ 1] ; mCHt1[1]=ffC[ 2] ; mCHt2[1]=ffC[ 4] ;
975  mCHt0[2]=ffC[ 3] ; mCHt1[2]=ffC[ 4] ; mCHt2[2]=ffC[ 5] ;
976  mCHt0[3]=ffC[ 6] ; mCHt1[3]=ffC[ 7] ; mCHt2[3]=ffC[ 8] ;
977  mCHt0[4]=ffC[10] ; mCHt1[4]=ffC[11] ; mCHt2[4]=ffC[12] ;
978  mCHt0[5]=ffC[15] ; mCHt1[5]=ffC[16] ; mCHt2[5]=ffC[17] ;
979  mCHt0[6]=ffC[21] ; mCHt1[6]=ffC[22] ; mCHt2[6]=ffC[23] ;
980 
981  //* Kalman gain K = mCH'*S
982 
983  Double_t k0[7], k1[7], k2[7];
984 
985  for(Int_t i=0;i<7;++i){
986  k0[i] = mCHt0[i]*mS[0] + mCHt1[i]*mS[1] + mCHt2[i]*mS[3];
987  k1[i] = mCHt0[i]*mS[1] + mCHt1[i]*mS[2] + mCHt2[i]*mS[4];
988  k2[i] = mCHt0[i]*mS[3] + mCHt1[i]*mS[4] + mCHt2[i]*mS[5];
989  }
990 
991  //* New estimation of the vertex position
992 
993  if( iter<maxIter-1 ){
994  for(Int_t i=0; i<3; ++i)
995  fVtxGuess[i]= ffP[i] + k0[i]*zeta[0]+k1[i]*zeta[1]+k2[i]*zeta[2];
996  continue;
997  }
998 
999  // last itearation -> update the particle
1000 
1001  //* VHt = VH'
1002 
1003  Double_t mVHt0[7], mVHt1[7], mVHt2[7];
1004 
1005  mVHt0[0]=mV[ 0] ; mVHt1[0]=mV[ 1] ; mVHt2[0]=mV[ 3] ;
1006  mVHt0[1]=mV[ 1] ; mVHt1[1]=mV[ 2] ; mVHt2[1]=mV[ 4] ;
1007  mVHt0[2]=mV[ 3] ; mVHt1[2]=mV[ 4] ; mVHt2[2]=mV[ 5] ;
1008  mVHt0[3]=mV[ 6] ; mVHt1[3]=mV[ 7] ; mVHt2[3]=mV[ 8] ;
1009  mVHt0[4]=mV[10] ; mVHt1[4]=mV[11] ; mVHt2[4]=mV[12] ;
1010  mVHt0[5]=mV[15] ; mVHt1[5]=mV[16] ; mVHt2[5]=mV[17] ;
1011  mVHt0[6]=mV[21] ; mVHt1[6]=mV[22] ; mVHt2[6]=mV[23] ;
1012 
1013  //* Kalman gain Km = mCH'*S
1014 
1015  Double_t km0[7], km1[7], km2[7];
1016 
1017  for(Int_t i=0;i<7;++i){
1018  km0[i] = mVHt0[i]*mS[0] + mVHt1[i]*mS[1] + mVHt2[i]*mS[3];
1019  km1[i] = mVHt0[i]*mS[1] + mVHt1[i]*mS[2] + mVHt2[i]*mS[4];
1020  km2[i] = mVHt0[i]*mS[3] + mVHt1[i]*mS[4] + mVHt2[i]*mS[5];
1021  }
1022 
1023  for(Int_t i=0;i<7;++i)
1024  ffP[i] = ffP[i] + k0[i]*zeta[0] + k1[i]*zeta[1] + k2[i]*zeta[2];
1025 
1026  for(Int_t i=0;i<7;++i)
1027  m[i] = m[i] - km0[i]*zeta[0] - km1[i]*zeta[1] - km2[i]*zeta[2];
1028 
1029  for(Int_t i=0, k=0;i<7;++i){
1030  for(Int_t j=0;j<=i;++j,++k){
1031  ffC[k] = ffC[k] - (k0[i]*mCHt0[j] + k1[i]*mCHt1[j] + k2[i]*mCHt2[j] );
1032  }
1033  }
1034 
1035  for(Int_t i=0, k=0;i<7;++i){
1036  for(Int_t j=0;j<=i;++j,++k){
1037  mV[k] = mV[k] - (km0[i]*mVHt0[j] + km1[i]*mVHt1[j] + km2[i]*mVHt2[j] );
1038  }
1039  }
1040 
1041  Double_t mDf[7][7];
1042 
1043  for(Int_t i=0;i<7;++i){
1044  for(Int_t j=0;j<7;++j){
1045  mDf[i][j] = (km0[i]*mCHt0[j] + km1[i]*mCHt1[j] + km2[i]*mCHt2[j] );
1046  }
1047  }
1048 
1049  Double_t mJ1[7][7], mJ2[7][7];
1050  for(Int_t iPar1=0; iPar1<7; iPar1++)
1051  {
1052  for(Int_t iPar2=0; iPar2<7; iPar2++)
1053  {
1054  mJ1[iPar1][iPar2] = 0;
1055  mJ2[iPar1][iPar2] = 0;
1056  }
1057  }
1058 
1059  Double_t mMassParticle = ffP[6]*ffP[6] - (ffP[3]*ffP[3] + ffP[4]*ffP[4] + ffP[5]*ffP[5]);
1060  Double_t mMassDaughter = m[6]*m[6] - (m[3]*m[3] + m[4]*m[4] + m[5]*m[5]);
1061  if(mMassParticle > 0) mMassParticle = sqrt(mMassParticle);
1062  if(mMassDaughter > 0) mMassDaughter = sqrt(mMassDaughter);
1063 
1064  if( fMassHypo > -0.5)
1065  SetMassConstraint(ffP,ffC,mJ1,fMassHypo);
1066  else if((mMassParticle < SumDaughterMass) || (ffP[6]<0) )
1067  SetMassConstraint(ffP,ffC,mJ1,SumDaughterMass);
1068 
1069  if(Daughter.fMassHypo > -0.5)
1070  SetMassConstraint(m,mV,mJ2,Daughter.fMassHypo);
1071  else if((mMassDaughter < Daughter.SumDaughterMass) || (m[6] < 0) )
1072  SetMassConstraint(m,mV,mJ2,Daughter.SumDaughterMass);
1073 
1074  Double_t mDJ[7][7];
1075 
1076  for(Int_t i=0; i<7; i++) {
1077  for(Int_t j=0; j<7; j++) {
1078  mDJ[i][j] = 0;
1079  for(Int_t k=0; k<7; k++) {
1080  mDJ[i][j] += mDf[i][k]*mJ1[j][k];
1081  }
1082  }
1083  }
1084 
1085  for(Int_t i=0; i<7; ++i){
1086  for(Int_t j=0; j<7; ++j){
1087  mDf[i][j]=0;
1088  for(Int_t l=0; l<7; l++){
1089  mDf[i][j] += mJ2[i][l]*mDJ[l][j];
1090  }
1091  }
1092  }
1093 
1094  //* Add the daughter momentum to the particle momentum
1095 
1096  ffP[ 3] += m[ 3];
1097  ffP[ 4] += m[ 4];
1098  ffP[ 5] += m[ 5];
1099  ffP[ 6] += m[ 6];
1100 
1101  ffC[ 9] += mV[ 9];
1102  ffC[13] += mV[13];
1103  ffC[14] += mV[14];
1104  ffC[18] += mV[18];
1105  ffC[19] += mV[19];
1106  ffC[20] += mV[20];
1107  ffC[24] += mV[24];
1108  ffC[25] += mV[25];
1109  ffC[26] += mV[26];
1110  ffC[27] += mV[27];
1111 
1112  ffC[6 ] += mDf[3][0]; ffC[7 ] += mDf[3][1]; ffC[8 ] += mDf[3][2];
1113  ffC[10] += mDf[4][0]; ffC[11] += mDf[4][1]; ffC[12] += mDf[4][2];
1114  ffC[15] += mDf[5][0]; ffC[16] += mDf[5][1]; ffC[17] += mDf[5][2];
1115  ffC[21] += mDf[6][0]; ffC[22] += mDf[6][1]; ffC[23] += mDf[6][2];
1116 
1117  ffC[9 ] += mDf[3][3] + mDf[3][3];
1118  ffC[13] += mDf[4][3] + mDf[3][4]; ffC[14] += mDf[4][4] + mDf[4][4];
1119  ffC[18] += mDf[5][3] + mDf[3][5]; ffC[19] += mDf[5][4] + mDf[4][5]; ffC[20] += mDf[5][5] + mDf[5][5];
1120  ffC[24] += mDf[6][3] + mDf[3][6]; ffC[25] += mDf[6][4] + mDf[4][6]; ffC[26] += mDf[6][5] + mDf[5][6]; ffC[27] += mDf[6][6] + mDf[6][6];
1121 
1122  //* New estimation of the vertex position r += K*zeta
1123 
1124  for(Int_t i=0;i<7;++i)
1125  fP[i] = ffP[i];
1126 
1127  //* New covariance matrix C -= K*(mCH')'
1128 
1129  for(Int_t i=0, k=0;i<7;++i){
1130  for(Int_t j=0;j<=i;++j,++k){
1131  fC[k] = ffC[k];
1132  }
1133  }
1134  //* Calculate Chi^2
1135 
1136  fNDF += 2;
1137  fQ += Daughter.GetQ();
1138  fSFromDecay = 0;
1139  fChi2 += (mS[0]*zeta[0] + mS[1]*zeta[1] + mS[3]*zeta[2])*zeta[0]
1140  + (mS[1]*zeta[0] + mS[2]*zeta[1] + mS[4]*zeta[2])*zeta[1]
1141  + (mS[3]*zeta[0] + mS[4]*zeta[1] + mS[5]*zeta[2])*zeta[2];
1142  }
1143 }
void TransportToDecayVertex()
Int_t i
Definition: run_full.C:25
__m128 m
Definition: P4_F32vec4.h:28
TTree * b
Double_t fVtxGuess[3]
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
TLorentzVector s
Definition: Pnd2DStar.C:50
void GetMeasurement(const Double_t XYZ[], Double_t m[], Double_t V[]) const
virtual void GetDStoParticle(const KFParticleBase &p, Double_t &DS, Double_t &DSp) const =0
virtual void Transport(Double_t dS, Double_t P[], Double_t C[]) const =0
Double_t SumDaughterMass
virtual void GetFieldValue(const Double_t xyz[], Double_t B[]) const =0
Double_t
void SetMassConstraint(Double_t Mass, Double_t SigmaMass=0)
Double_t fSFromDecay
Double_t fC[36]
void TransportToDS(Double_t dS)
Double_t fMassHypo
Int_t GetQ() const
const Double_t& KFParticleBase::Chi2 ( ) const
inline

Definition at line 122 of file KFParticleBase.h.

References fChi2.

Referenced by KFParticle::Chi2().

122 { return fChi2; }
Double_t& KFParticleBase::Chi2 ( )
inline

Definition at line 155 of file KFParticleBase.h.

References fChi2.

155 { return fChi2; }
Double_t& KFParticleBase::Cij ( Int_t  i,
Int_t  j 
)
inlineprotected

Definition at line 285 of file KFParticleBase.h.

References fC, and IJ().

Referenced by SetMassConstraint().

285 { return fC[IJ(i,j)]; }
Int_t i
Definition: run_full.C:25
static Int_t IJ(Int_t i, Int_t j)
Double_t fC[36]
void KFParticleBase::Construct ( const KFParticleBase vDaughters[],
Int_t  nDaughters,
const KFParticleBase ProdVtx = 0,
Double_t  Mass = -1,
Bool_t  IsConstrained = 0 
)

Definition at line 1477 of file KFParticleBase.cxx.

References AddDaughter(), C(), Double_t, fAtProductionVertex, fC, fChi2, fIsLinearized, fNDF, fP, fQ, fSFromDecay, fVtxGuess, GetDStoParticle(), i, P, SetMassConstraint(), SetProductionVertex(), SumDaughterMass, and Transport().

Referenced by KFParticle::Construct().

1479 {
1480  //* Full reconstruction in one go
1481 
1482  Int_t maxIter = 1;
1483  bool wasLinearized = fIsLinearized;
1484  if( !fIsLinearized || IsConstrained ){
1485  //fVtxGuess[0] = fVtxGuess[1] = fVtxGuess[2] = 0; //!!!!
1486  Double_t ds=0., ds1=0.;
1487  Double_t P[8], C[36];
1488  vDaughters[0]->GetDStoParticle(*vDaughters[1], ds, ds1);
1489  vDaughters[0]->Transport(ds,P,C);
1490  fVtxGuess[0] = P[0];
1491  fVtxGuess[1] = P[1];
1492  fVtxGuess[2] = P[2];
1493 /* fVtxGuess[0] = GetX();
1494  fVtxGuess[1] = GetY();
1495  fVtxGuess[2] = GetZ();*/
1496 
1497  fIsLinearized = 1;
1498  maxIter = 3;
1499  }
1500 
1501  Double_t constraintC[6];
1502 
1503  if( IsConstrained ){
1504  for(Int_t i=0;i<6;++i) constraintC[i]=fC[i];
1505  } else {
1506  for(Int_t i=0;i<6;++i) constraintC[i]=0.;
1507 // constraintC[0] = constraintC[2] = constraintC[5] = 100.;
1508  constraintC[0] = 100.*vDaughters[0]->fC[0];
1509  constraintC[2] = 100.*vDaughters[0]->fC[2];
1510  constraintC[5] = 100.*vDaughters[0]->fC[5];
1511  }
1512 
1513 
1514  for( Int_t iter=0; iter<maxIter; iter++ ){
1515  fAtProductionVertex = 0;
1516  fSFromDecay = 0;
1517  fP[0] = fVtxGuess[0];
1518  fP[1] = fVtxGuess[1];
1519  fP[2] = fVtxGuess[2];
1520  fP[3] = 0;
1521  fP[4] = 0;
1522  fP[5] = 0;
1523  fP[6] = 0;
1524  fP[7] = 0;
1525  SumDaughterMass = 0;
1526 
1527  for(Int_t i=0;i<6; ++i) fC[i]=constraintC[i];
1528  for(Int_t i=6;i<36;++i) fC[i]=0.;
1529  fC[35] = 1.;
1530 
1531  fNDF = IsConstrained ?0 :-3;
1532  fChi2 = 0.;
1533  fQ = 0;
1534 
1535  for( Int_t itr =0; itr<nDaughters; itr++ ){
1536  AddDaughter( *vDaughters[itr] );
1537  }
1538  if( iter<maxIter-1){
1539  for( Int_t i=0; i<3; i++ ) fVtxGuess[i] = fP[i];
1540  }
1541  }
1542  fIsLinearized = wasLinearized;
1543 
1544  if( Mass>=0 ) SetMassConstraint( Mass );
1545  if( Parent ) SetProductionVertex( *Parent );
1546 }
Int_t i
Definition: run_full.C:25
Double_t fVtxGuess[3]
Double_t fP[8]
virtual void GetDStoParticle(const KFParticleBase &p, Double_t &DS, Double_t &DSp) const =0
int Pic_FED Eff_lEE C()
void SetProductionVertex(const KFParticleBase &Vtx)
virtual void Transport(Double_t dS, Double_t P[], Double_t C[]) const =0
Double_t SumDaughterMass
Double_t
void SetMassConstraint(Double_t Mass, Double_t SigmaMass=0)
Double_t fSFromDecay
Double_t fC[36]
GeV c P
Bool_t fAtProductionVertex
void AddDaughter(const KFParticleBase &Daughter)
void KFParticleBase::ConstructGammaBz ( const KFParticleBase daughter1,
const KFParticleBase daughter2,
double  Bz 
)

Definition at line 2821 of file KFParticleBase.cxx.

References d0, Double_t, fAtProductionVertex, fC, fChi2, fIsLinearized, fNDF, fP, fQ, fSFromDecay, fVtxGuess, GetDStoParticle(), GetDStoPoint(), GetMeasurement(), GetQ(), i, IJ(), InvertSym3(), m, p, r, SetNonlinearMassConstraint(), sqrt(), Transport(), and z.

Referenced by KFParticle::ConstructGamma().

2823 {
2824  //* Create gamma
2825 
2826  const KFParticleBase *daughters[2] = { &daughter1, &daughter2};
2827 
2828  double v0[3];
2829 
2830  if( !fIsLinearized ){
2831  Double_t ds, ds1;
2832  Double_t m[8];
2833  Double_t mCd[36];
2834  daughter1.GetDStoParticle(daughter2, ds, ds1);
2835  daughter1.Transport( ds, m, mCd );
2836  fP[0] = m[0];
2837  fP[1] = m[1];
2838  fP[2] = m[2];
2839  daughter2.Transport( ds1, m, mCd );
2840  fP[0] = .5*( fP[0] + m[0] );
2841  fP[1] = .5*( fP[1] + m[1] );
2842  fP[2] = .5*( fP[2] + m[2] );
2843  } else {
2844  fP[0] = fVtxGuess[0];
2845  fP[1] = fVtxGuess[1];
2846  fP[2] = fVtxGuess[2];
2847  }
2848 
2849  double daughterP[2][8], daughterC[2][36];
2850  double vtxMom[2][3];
2851 
2852  int nIter = fIsLinearized ?1 :2;
2853 
2854  for( int iter=0; iter<nIter; iter++){
2855 
2856  v0[0] = fP[0];
2857  v0[1] = fP[1];
2858  v0[2] = fP[2];
2859 
2860  fAtProductionVertex = 0;
2861  fSFromDecay = 0;
2862  fP[0] = v0[0];
2863  fP[1] = v0[1];
2864  fP[2] = v0[2];
2865  fP[3] = 0;
2866  fP[4] = 0;
2867  fP[5] = 0;
2868  fP[6] = 0;
2869  fP[7] = 0;
2870 
2871 
2872  // fit daughters to the vertex guess
2873 
2874  {
2875  for( int id=0; id<2; id++ ){
2876 
2877  double *p = daughterP[id];
2878  double *mC = daughterC[id];
2879 
2880  daughters[id]->GetMeasurement( v0, p, mC );
2881 
2882  Double_t mAi[6];
2883  InvertSym3(mC, mAi );
2884 
2885  Double_t mB[3][3];
2886 
2887  mB[0][0] = mC[ 6]*mAi[0] + mC[ 7]*mAi[1] + mC[ 8]*mAi[3];
2888  mB[0][1] = mC[ 6]*mAi[1] + mC[ 7]*mAi[2] + mC[ 8]*mAi[4];
2889  mB[0][2] = mC[ 6]*mAi[3] + mC[ 7]*mAi[4] + mC[ 8]*mAi[5];
2890 
2891  mB[1][0] = mC[10]*mAi[0] + mC[11]*mAi[1] + mC[12]*mAi[3];
2892  mB[1][1] = mC[10]*mAi[1] + mC[11]*mAi[2] + mC[12]*mAi[4];
2893  mB[1][2] = mC[10]*mAi[3] + mC[11]*mAi[4] + mC[12]*mAi[5];
2894 
2895  mB[2][0] = mC[15]*mAi[0] + mC[16]*mAi[1] + mC[17]*mAi[3];
2896  mB[2][1] = mC[15]*mAi[1] + mC[16]*mAi[2] + mC[17]*mAi[4];
2897  mB[2][2] = mC[15]*mAi[3] + mC[16]*mAi[4] + mC[17]*mAi[5];
2898 
2899  Double_t z[3] = { v0[0]-p[0], v0[1]-p[1], v0[2]-p[2] };
2900 
2901  vtxMom[id][0] = p[3] + mB[0][0]*z[0] + mB[0][1]*z[1] + mB[0][2]*z[2];
2902  vtxMom[id][1] = p[4] + mB[1][0]*z[0] + mB[1][1]*z[1] + mB[1][2]*z[2];
2903  vtxMom[id][2] = p[5] + mB[2][0]*z[0] + mB[2][1]*z[1] + mB[2][2]*z[2];
2904 
2905  daughters[id]->Transport( daughters[id]->GetDStoPoint(v0), p, mC );
2906 
2907  }
2908 
2909  } // fit daughters to guess
2910 
2911 
2912  // fit new vertex
2913  {
2914 
2915  double mpx0 = vtxMom[0][0]+vtxMom[1][0];
2916  double mpy0 = vtxMom[0][1]+vtxMom[1][1];
2917  double mpt0 = sqrt(mpx0*mpx0 + mpy0*mpy0);
2918  // double a0 = atan2(mpy0,mpx0);
2919 
2920  double ca0 = mpx0/mpt0;
2921  double sa0 = mpy0/mpt0;
2922  double r[3] = { v0[0], v0[1], v0[2] };
2923  double mC[3][3] = {{1000., 0 , 0 },
2924  {0, 1000., 0 },
2925  {0, 0, 1000. } };
2926  double chi2=0;
2927 
2928  for( int id=0; id<2; id++ ){
2929  const Double_t kCLight = 0.000299792458;
2930  Double_t q = Bz*daughters[id]->GetQ()*kCLight;
2931  Double_t px0 = vtxMom[id][0];
2932  Double_t py0 = vtxMom[id][1];
2933  Double_t pz0 = vtxMom[id][2];
2934  Double_t pt0 = sqrt(px0*px0+py0*py0);
2935  Double_t mG[3][6], mB[3], mH[3][3];
2936  // r = {vx,vy,vz};
2937  // m = {x,y,z,Px,Py,Pz};
2938  // V = daughter.C
2939  // G*m + B = H*r;
2940  // q*x + Py - q*vx - sin(a)*Pt = 0
2941  // q*y - Px - q*vy + cos(a)*Pt = 0
2942  // (Px*cos(a) + Py*sin(a) ) (vz -z) - Pz( cos(a)*(vx-x) + sin(a)*(vy-y)) = 0
2943 
2944  mG[0][0] = q;
2945  mG[0][1] = 0;
2946  mG[0][2] = 0;
2947  mG[0][3] = -sa0*px0/pt0;
2948  mG[0][4] = 1 -sa0*py0/pt0;
2949  mG[0][5] = 0;
2950  mH[0][0] = q;
2951  mH[0][1] = 0;
2952  mH[0][2] = 0;
2953  mB[0] = py0 - sa0*pt0 - mG[0][3]*px0 - mG[0][4]*py0 ;
2954 
2955  // q*y - Px - q*vy + cos(a)*Pt = 0
2956 
2957  mG[1][0] = 0;
2958  mG[1][1] = q;
2959  mG[1][2] = 0;
2960  mG[1][3] = -1 + ca0*px0/pt0;
2961  mG[1][4] = + ca0*py0/pt0;
2962  mG[1][5] = 0;
2963  mH[1][0] = 0;
2964  mH[1][1] = q;
2965  mH[1][2] = 0;
2966  mB[1] = -px0 + ca0*pt0 - mG[1][3]*px0 - mG[1][4]*py0 ;
2967 
2968  // (Px*cos(a) + Py*sin(a) ) (z -vz) - Pz( cos(a)*(x-vx) + sin(a)*(y-vy)) = 0
2969 
2970  mG[2][0] = -pz0*ca0;
2971  mG[2][1] = -pz0*sa0;
2972  mG[2][2] = px0*ca0 + py0*sa0;
2973  mG[2][3] = 0;
2974  mG[2][4] = 0;
2975  mG[2][5] = 0;
2976 
2977  mH[2][0] = mG[2][0];
2978  mH[2][1] = mG[2][1];
2979  mH[2][2] = mG[2][2];
2980 
2981  mB[2] = 0;
2982 
2983  // fit the vertex
2984 
2985  // V = GVGt
2986 
2987  double mGV[3][6];
2988  double mV[6];
2989  double m[3];
2990  for( int i=0; i<3; i++ ){
2991  m[i] = mB[i];
2992  for( int k=0; k<6; k++ ) m[i]+=mG[i][k]*daughterP[id][k];
2993  }
2994  for( int i=0; i<3; i++ ){
2995  for( int j=0; j<6; j++ ){
2996  mGV[i][j] = 0;
2997  for( int k=0; k<6; k++ ) mGV[i][j]+=mG[i][k]*daughterC[id][ IJ(k,j) ];
2998  }
2999  }
3000  for( int i=0, k=0; i<3; i++ ){
3001  for( int j=0; j<=i; j++,k++ ){
3002  mV[k] = 0;
3003  for( int l=0; l<6; l++ ) mV[k]+=mGV[i][l]*mG[j][l];
3004  }
3005  }
3006 
3007 
3008  //* CHt
3009 
3010  Double_t mCHt[3][3];
3011  Double_t mHCHt[6];
3012  Double_t mHr[3];
3013  for( int i=0; i<3; i++ ){
3014  mHr[i] = 0;
3015  for( int k=0; k<3; k++ ) mHr[i]+= mH[i][k]*r[k];
3016  }
3017 
3018  for( int i=0; i<3; i++ ){
3019  for( int j=0; j<3; j++){
3020  mCHt[i][j] = 0;
3021  for( int k=0; k<3; k++ ) mCHt[i][j]+= mC[i][k]*mH[j][k];
3022  }
3023  }
3024 
3025  for( int i=0, k=0; i<3; i++ ){
3026  for( int j=0; j<=i; j++, k++ ){
3027  mHCHt[k] = 0;
3028  for( int l=0; l<3; l++ ) mHCHt[k]+= mH[i][l]*mCHt[l][j];
3029  }
3030  }
3031 
3032  Double_t mS[6] = { mHCHt[0]+mV[0],
3033  mHCHt[1]+mV[1], mHCHt[2]+mV[2],
3034  mHCHt[3]+mV[3], mHCHt[4]+mV[4], mHCHt[5]+mV[5] };
3035 
3036 
3037  InvertSym3(mS,mS);
3038 
3039  //* Residual (measured - estimated)
3040 
3041  Double_t zeta[3] = { m[0]-mHr[0], m[1]-mHr[1], m[2]-mHr[2] };
3042 
3043  //* Kalman gain K = mCH'*S
3044 
3045  Double_t k[3][3];
3046 
3047  for(Int_t i=0;i<3;++i){
3048  k[i][0] = mCHt[i][0]*mS[0] + mCHt[i][1]*mS[1] + mCHt[i][2]*mS[3];
3049  k[i][1] = mCHt[i][0]*mS[1] + mCHt[i][1]*mS[2] + mCHt[i][2]*mS[4];
3050  k[i][2] = mCHt[i][0]*mS[3] + mCHt[i][1]*mS[4] + mCHt[i][2]*mS[5];
3051  }
3052 
3053  //* New estimation of the vertex position r += K*zeta
3054 
3055  for(Int_t i=0;i<3;++i)
3056  r[i] = r[i] + k[i][0]*zeta[0] + k[i][1]*zeta[1] + k[i][2]*zeta[2];
3057 
3058  //* New covariance matrix C -= K*(mCH')'
3059 
3060  for(Int_t i=0;i<3;++i){
3061  for(Int_t j=0;j<=i;++j){
3062  mC[i][j] = mC[i][j] - (k[i][0]*mCHt[j][0] + k[i][1]*mCHt[j][1] + k[i][2]*mCHt[j][2]);
3063  mC[j][i] = mC[i][j];
3064  }
3065  }
3066 
3067  //* Calculate Chi^2
3068 
3069  chi2 += ( ( mS[0]*zeta[0] + mS[1]*zeta[1] + mS[3]*zeta[2] )*zeta[0]
3070  +(mS[1]*zeta[0] + mS[2]*zeta[1] + mS[4]*zeta[2] )*zeta[1]
3071  +(mS[3]*zeta[0] + mS[4]*zeta[1] + mS[5]*zeta[2] )*zeta[2] );
3072  }
3073 
3074  // store vertex
3075 
3076  fNDF = 2;
3077  fChi2 = chi2;
3078  for( int i=0; i<3; i++ ) fP[i] = r[i];
3079  for( int i=0,k=0; i<3; i++ ){
3080  for( int j=0; j<=i; j++,k++ ){
3081  fC[k] = mC[i][j];
3082  }
3083  }
3084  }
3085 
3086  } // iterations
3087 
3088  // now fit daughters to the vertex
3089 
3090  fQ = 0;
3091  fSFromDecay = 0;
3092 
3093  for(Int_t i=3;i<8;++i) fP[i]=0.;
3094  for(Int_t i=6;i<35;++i) fC[i]=0.;
3095  fC[35] = 100.;
3096 
3097  for( int id=0; id<2; id++ ){
3098 
3099  double *p = daughterP[id];
3100  double *mC = daughterC[id];
3101  daughters[id]->GetMeasurement( v0, p, mC );
3102 
3103  const Double_t *m = fP, *mV = fC;
3104 
3105  Double_t mAi[6];
3106  InvertSym3(mC, mAi );
3107 
3108  Double_t mB[4][3];
3109 
3110  mB[0][0] = mC[ 6]*mAi[0] + mC[ 7]*mAi[1] + mC[ 8]*mAi[3];
3111  mB[0][1] = mC[ 6]*mAi[1] + mC[ 7]*mAi[2] + mC[ 8]*mAi[4];
3112  mB[0][2] = mC[ 6]*mAi[3] + mC[ 7]*mAi[4] + mC[ 8]*mAi[5];
3113 
3114  mB[1][0] = mC[10]*mAi[0] + mC[11]*mAi[1] + mC[12]*mAi[3];
3115  mB[1][1] = mC[10]*mAi[1] + mC[11]*mAi[2] + mC[12]*mAi[4];
3116  mB[1][2] = mC[10]*mAi[3] + mC[11]*mAi[4] + mC[12]*mAi[5];
3117 
3118  mB[2][0] = mC[15]*mAi[0] + mC[16]*mAi[1] + mC[17]*mAi[3];
3119  mB[2][1] = mC[15]*mAi[1] + mC[16]*mAi[2] + mC[17]*mAi[4];
3120  mB[2][2] = mC[15]*mAi[3] + mC[16]*mAi[4] + mC[17]*mAi[5];
3121 
3122  mB[3][0] = mC[21]*mAi[0] + mC[22]*mAi[1] + mC[23]*mAi[3];
3123  mB[3][1] = mC[21]*mAi[1] + mC[22]*mAi[2] + mC[23]*mAi[4];
3124  mB[3][2] = mC[21]*mAi[3] + mC[22]*mAi[4] + mC[23]*mAi[5];
3125 
3126 
3127  Double_t z[3] = { m[0]-p[0], m[1]-p[1], m[2]-p[2] };
3128 
3129 // {
3130 // Double_t mAV[6] = { mC[0]-mV[0], mC[1]-mV[1], mC[2]-mV[2],
3131 // mC[3]-mV[3], mC[4]-mV[4], mC[5]-mV[5] };
3132 //
3133 // Double_t mAVi[6];
3134 // if( !InvertSym3(mAV, mAVi) ){
3135 // Double_t dChi2 = ( +(mAVi[0]*z[0] + mAVi[1]*z[1] + mAVi[3]*z[2])*z[0]
3136 // +(mAVi[1]*z[0] + mAVi[2]*z[1] + mAVi[4]*z[2])*z[1]
3137 // +(mAVi[3]*z[0] + mAVi[4]*z[1] + mAVi[5]*z[2])*z[2] );
3138 // fChi2+= fabs( dChi2 );
3139 // }
3140 // fNDF += 2;
3141 // }
3142 
3143  //* Add the daughter momentum to the particle momentum
3144 
3145  fP[3]+= p[3] + mB[0][0]*z[0] + mB[0][1]*z[1] + mB[0][2]*z[2];
3146  fP[4]+= p[4] + mB[1][0]*z[0] + mB[1][1]*z[1] + mB[1][2]*z[2];
3147  fP[5]+= p[5] + mB[2][0]*z[0] + mB[2][1]*z[1] + mB[2][2]*z[2];
3148  fP[6]+= p[6] + mB[3][0]*z[0] + mB[3][1]*z[1] + mB[3][2]*z[2];
3149 
3150  Double_t d0, d1, d2;
3151 
3152  d0= mB[0][0]*mV[0] + mB[0][1]*mV[1] + mB[0][2]*mV[3] - mC[ 6];
3153  d1= mB[0][0]*mV[1] + mB[0][1]*mV[2] + mB[0][2]*mV[4] - mC[ 7];
3154  d2= mB[0][0]*mV[3] + mB[0][1]*mV[4] + mB[0][2]*mV[5] - mC[ 8];
3155 
3156  //fC[6]+= mC[ 6] + d0;
3157  //fC[7]+= mC[ 7] + d1;
3158  //fC[8]+= mC[ 8] + d2;
3159  fC[9]+= mC[ 9] + d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
3160 
3161  d0= mB[1][0]*mV[0] + mB[1][1]*mV[1] + mB[1][2]*mV[3] - mC[10];
3162  d1= mB[1][0]*mV[1] + mB[1][1]*mV[2] + mB[1][2]*mV[4] - mC[11];
3163  d2= mB[1][0]*mV[3] + mB[1][1]*mV[4] + mB[1][2]*mV[5] - mC[12];
3164 
3165  //fC[10]+= mC[10]+ d0;
3166  //fC[11]+= mC[11]+ d1;
3167  //fC[12]+= mC[12]+ d2;
3168  fC[13]+= mC[13]+ d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
3169  fC[14]+= mC[14]+ d0*mB[1][0] + d1*mB[1][1] + d2*mB[1][2];
3170 
3171  d0= mB[2][0]*mV[0] + mB[2][1]*mV[1] + mB[2][2]*mV[3] - mC[15];
3172  d1= mB[2][0]*mV[1] + mB[2][1]*mV[2] + mB[2][2]*mV[4] - mC[16];
3173  d2= mB[2][0]*mV[3] + mB[2][1]*mV[4] + mB[2][2]*mV[5] - mC[17];
3174 
3175  //fC[15]+= mC[15]+ d0;
3176  //fC[16]+= mC[16]+ d1;
3177  //fC[17]+= mC[17]+ d2;
3178  fC[18]+= mC[18]+ d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
3179  fC[19]+= mC[19]+ d0*mB[1][0] + d1*mB[1][1] + d2*mB[1][2];
3180  fC[20]+= mC[20]+ d0*mB[2][0] + d1*mB[2][1] + d2*mB[2][2];
3181 
3182  d0= mB[3][0]*mV[0] + mB[3][1]*mV[1] + mB[3][2]*mV[3] - mC[21];
3183  d1= mB[3][0]*mV[1] + mB[3][1]*mV[2] + mB[3][2]*mV[4] - mC[22];
3184  d2= mB[3][0]*mV[3] + mB[3][1]*mV[4] + mB[3][2]*mV[5] - mC[23];
3185 
3186  //fC[21]+= mC[21] + d0;
3187  //fC[22]+= mC[22] + d1;
3188  //fC[23]+= mC[23] + d2;
3189  fC[24]+= mC[24] + d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
3190  fC[25]+= mC[25] + d0*mB[1][0] + d1*mB[1][1] + d2*mB[1][2];
3191  fC[26]+= mC[26] + d0*mB[2][0] + d1*mB[2][1] + d2*mB[2][2];
3192  fC[27]+= mC[27] + d0*mB[3][0] + d1*mB[3][1] + d2*mB[3][2];
3193  }
3194 
3195 // SetMassConstraint(0,0);
3197 }
double r
Definition: RiemannTest.C:14
void SetNonlinearMassConstraint(Double_t Mass)
Int_t i
Definition: run_full.C:25
__m128 m
Definition: P4_F32vec4.h:28
Double_t fVtxGuess[3]
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
void GetMeasurement(const Double_t XYZ[], Double_t m[], Double_t V[]) const
static Int_t IJ(Int_t i, Int_t j)
Double_t p
Definition: anasim.C:58
virtual void GetDStoParticle(const KFParticleBase &p, Double_t &DS, Double_t &DSp) const =0
Double_t d0
Definition: checkhelixhit.C:59
virtual void Transport(Double_t dS, Double_t P[], Double_t C[]) const =0
Double_t
static Bool_t InvertSym3(const Double_t A[], Double_t Ainv[])
Double_t fSFromDecay
Double_t z
Double_t fC[36]
Bool_t fAtProductionVertex
Int_t GetQ() const
virtual Double_t GetDStoPoint(const Double_t xyz[]) const =0
void KFParticleBase::Convert ( bool  ToProduction)
protected

Definition at line 1549 of file KFParticleBase.cxx.

References c, Double_t, fC, fP, fQ, GetFieldValue(), and h.

Referenced by SetProductionVertex(), TransportToDecayVertex(), and TransportToProductionVertex().

1550 {
1551  //* Tricky function - convert the particle error along its trajectory to
1552  //* the value which corresponds to its production/decay vertex
1553  //* It is done by combination of the error of decay length with the position errors
1554 
1555  Double_t fld[3];
1556  {
1557  GetFieldValue( fP, fld );
1558  const Double_t kCLight = fQ*0.000299792458;
1559  fld[0]*=kCLight; fld[1]*=kCLight; fld[2]*=kCLight;
1560  }
1561 
1562  Double_t h[6];
1563 
1564  h[0] = fP[3];
1565  h[1] = fP[4];
1566  h[2] = fP[5];
1567  if( ToProduction ){ h[0]=-h[0]; h[1]=-h[1]; h[2]=-h[2]; }
1568  h[3] = h[1]*fld[2]-h[2]*fld[1];
1569  h[4] = h[2]*fld[0]-h[0]*fld[2];
1570  h[5] = h[0]*fld[1]-h[1]*fld[0];
1571 
1572  Double_t c;
1573 
1574  c = fC[28]+h[0]*fC[35];
1575  fC[ 0]+= h[0]*(c+fC[28]);
1576  fC[28] = c;
1577 
1578  fC[ 1]+= h[1]*fC[28] + h[0]*fC[29];
1579  c = fC[29]+h[1]*fC[35];
1580  fC[ 2]+= h[1]*(c+fC[29]);
1581  fC[29] = c;
1582 
1583  fC[ 3]+= h[2]*fC[28] + h[0]*fC[30];
1584  fC[ 4]+= h[2]*fC[29] + h[1]*fC[30];
1585  c = fC[30]+h[2]*fC[35];
1586  fC[ 5]+= h[2]*(c+fC[30]);
1587  fC[30] = c;
1588 
1589  fC[ 6]+= h[3]*fC[28] + h[0]*fC[31];
1590  fC[ 7]+= h[3]*fC[29] + h[1]*fC[31];
1591  fC[ 8]+= h[3]*fC[30] + h[2]*fC[31];
1592  c = fC[31]+h[3]*fC[35];
1593  fC[ 9]+= h[3]*(c+fC[31]);
1594  fC[31] = c;
1595 
1596  fC[10]+= h[4]*fC[28] + h[0]*fC[32];
1597  fC[11]+= h[4]*fC[29] + h[1]*fC[32];
1598  fC[12]+= h[4]*fC[30] + h[2]*fC[32];
1599  fC[13]+= h[4]*fC[31] + h[3]*fC[32];
1600  c = fC[32]+h[4]*fC[35];
1601  fC[14]+= h[4]*(c+fC[32]);
1602  fC[32] = c;
1603 
1604  fC[15]+= h[5]*fC[28] + h[0]*fC[33];
1605  fC[16]+= h[5]*fC[29] + h[1]*fC[33];
1606  fC[17]+= h[5]*fC[30] + h[2]*fC[33];
1607  fC[18]+= h[5]*fC[31] + h[3]*fC[33];
1608  fC[19]+= h[5]*fC[32] + h[4]*fC[33];
1609  c = fC[33]+h[5]*fC[35];
1610  fC[20]+= h[5]*(c+fC[33]);
1611  fC[33] = c;
1612 
1613  fC[21]+= h[0]*fC[34];
1614  fC[22]+= h[1]*fC[34];
1615  fC[23]+= h[2]*fC[34];
1616  fC[24]+= h[3]*fC[34];
1617  fC[25]+= h[4]*fC[34];
1618  fC[26]+= h[5]*fC[34];
1619 }
Double_t fP[8]
virtual void GetFieldValue(const Double_t xyz[], Double_t B[]) const =0
Double_t
Double_t fC[36]
Double_t& KFParticleBase::Covariance ( Int_t  i)
inline

Definition at line 159 of file KFParticleBase.h.

References fC, and i.

Referenced by KFParticle::Covariance(), and RotateXY().

159 { return fC[i]; }
Int_t i
Definition: run_full.C:25
Double_t fC[36]
Double_t& KFParticleBase::Covariance ( Int_t  i,
Int_t  j 
)
inline

Definition at line 160 of file KFParticleBase.h.

References fC, and IJ().

160 { return fC[IJ(i,j)]; }
Int_t i
Definition: run_full.C:25
static Int_t IJ(Int_t i, Int_t j)
Double_t fC[36]
const std::vector<int>& KFParticleBase::DaughterIds ( ) const
inline

Definition at line 271 of file KFParticleBase.h.

References fDaughtersIds.

Referenced by KFParticleFinder::CombineTrackPart(), KFParticleFinder::FindHyperons(), and KFParticleSIMD::KFParticleSIMD().

271 { return fDaughtersIds; };
std::vector< int > fDaughtersIds
const Double_t& KFParticleBase::E ( ) const
inline

Definition at line 119 of file KFParticleBase.h.

References fP.

Referenced by KFParticle::E().

119 { return fP[6]; }
Double_t fP[8]
Double_t& KFParticleBase::E ( )
inline

Definition at line 152 of file KFParticleBase.h.

References fP.

152 { return fP[6]; }
Double_t fP[8]
void KFParticleBase::GetArmenterosPodolanski ( KFParticleBase positive,
KFParticleBase negative,
Double_t  QtAlfa[2] 
)
static

Definition at line 3199 of file KFParticleBase.cxx.

References alpha, Double_t, GetPx(), GetPy(), GetPz(), and sqrt().

3200 {
3201 // example:
3202 // KFParticle PosParticle(...)
3203 // KFParticle NegParticle(...)
3204 // Gamma.ConstructGamma(PosParticle, NegParticle);
3205 // Double_t VertexGamma[3] = {Gamma.GetX(), Gamma.GetY(), Gamma.GetZ()};
3206 // PosParticle.TransportToPoint(VertexGamma);
3207 // NegParticle.TransportToPoint(VertexGamma);
3208 // Double_t armenterosQtAlfa[2] = {0.};
3209 // KFParticle::GetArmenterosPodolanski(PosParticle, NegParticle, armenterosQtAlfa );
3210 
3211  Double_t alpha = 0., qt = 0.;
3212  Double_t spx = positive.GetPx() + negative.GetPx();
3213  Double_t spy = positive.GetPy() + negative.GetPy();
3214  Double_t spz = positive.GetPz() + negative.GetPz();
3215  Double_t sp = sqrt(spx*spx + spy*spy + spz*spz);
3216  if( sp == 0.0) return;
3217  Double_t pn, pp, pln, plp;
3218 
3219  pn = sqrt(negative.GetPx()*negative.GetPx() + negative.GetPy()*negative.GetPy() + negative.GetPz()*negative.GetPz());
3220  pp = sqrt(positive.GetPx()*positive.GetPx() + positive.GetPy()*positive.GetPy() + positive.GetPz()*positive.GetPz());
3221  pln = (negative.GetPx()*spx+negative.GetPy()*spy+negative.GetPz()*spz)/sp;
3222  plp = (positive.GetPx()*spx+positive.GetPy()*spy+positive.GetPz()*spz)/sp;
3223 
3224  if( pn == 0.0) return;
3225  Double_t ptm = (1.-((pln/pn)*(pln/pn)));
3226  qt= (ptm>=0.)? pn*sqrt(ptm) :0;
3227  alpha = (plp-pln)/(plp+pln);
3228 
3229  QtAlfa[0] = qt;
3230  QtAlfa[1] = alpha;
3231 }
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t
Double_t GetPx() const
double alpha
Definition: f_Init.h:9
Double_t GetPy() const
Double_t GetPz() const
Double_t KFParticleBase::GetChi2 ( ) const
inline

Definition at line 110 of file KFParticleBase.h.

References fChi2.

Referenced by KFParticle::GetChi2().

110 { return fChi2; }
Double_t KFParticleBase::GetCovariance ( Int_t  i) const
inline

Definition at line 126 of file KFParticleBase.h.

References fC, and i.

Referenced by KFParticle::GetCovariance(), and RotateXY().

126 { return fC[i]; }
Int_t i
Definition: run_full.C:25
Double_t fC[36]
Double_t KFParticleBase::GetCovariance ( Int_t  i,
Int_t  j 
) const
inline

Definition at line 127 of file KFParticleBase.h.

References fC, and IJ().

127 { return fC[IJ(i,j)]; }
Int_t i
Definition: run_full.C:25
static Int_t IJ(Int_t i, Int_t j)
Double_t fC[36]
Int_t KFParticleBase::GetDecayLength ( Double_t L,
Double_t SigmaL 
) const

Definition at line 265 of file KFParticleBase.cxx.

References Double_t, fabs(), fC, fP, p2, sqrt(), t, x, y, and z.

Referenced by KFParticle::GetDecayLength(), and KFParticle::GetErrDecayLength().

266 {
267  //* Calculate particle decay length [cm]
268 
269  Double_t x = fP[3];
270  Double_t y = fP[4];
271  Double_t z = fP[5];
272  Double_t t = fP[7];
273  Double_t x2 = x*x;
274  Double_t y2 = y*y;
275  Double_t z2 = z*z;
276  Double_t p2 = x2+y2+z2;
277  l = t*sqrt(p2);
278  if( p2>1.e-4){
279  error = p2*fC[35] + t*t/p2*(x2*fC[9]+y2*fC[14]+z2*fC[20]
280  + 2*(x*y*fC[13]+x*z*fC[18]+y*z*fC[19]) )
281  + 2*t*(x*fC[31]+y*fC[32]+z*fC[33]);
282  error = sqrt(fabs(error));
283  return 0;
284  }
285  error = 1.e20;
286  return 1;
287 }
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
static void error(int no)
Definition: ranlxd.cxx:419
Double_t
Double_t z
Double_t fC[36]
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
TPad * p2
Definition: hist-t7.C:117
Double_t x
TTree * t
Definition: bump_analys.C:13
Double_t y
Int_t KFParticleBase::GetDecayLengthXY ( Double_t L,
Double_t SigmaL 
) const

Definition at line 289 of file KFParticleBase.cxx.

References Double_t, fabs(), fC, fP, sqrt(), t, x, and y.

Referenced by KFParticle::GetDecayLengthXY(), and KFParticle::GetErrDecayLengthXY().

290 {
291  //* Calculate particle decay length in XY projection [cm]
292 
293  Double_t x = fP[3];
294  Double_t y = fP[4];
295  Double_t t = fP[7];
296  Double_t x2 = x*x;
297  Double_t y2 = y*y;
298  Double_t pt2 = x2+y2;
299  l = t*sqrt(pt2);
300  if( pt2>1.e-4){
301  error = pt2*fC[35] + t*t/pt2*(x2*fC[9]+y2*fC[14] + 2*x*y*fC[13] )
302  + 2*t*(x*fC[31]+y*fC[32]);
303  error = sqrt(fabs(error));
304  return 0;
305  }
306  error = 1.e20;
307  return 1;
308 }
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
static void error(int no)
Definition: ranlxd.cxx:419
Double_t
Double_t fC[36]
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
Double_t x
TTree * t
Definition: bump_analys.C:13
Double_t y
Double_t KFParticleBase::GetDeviationFromParticle ( const KFParticleBase p) const

Definition at line 2625 of file KFParticleBase.cxx.

References d, Double_t, fP, GetDeviationFromVertex(), GetDStoParticle(), h, sqrt(), and Transport().

Referenced by KFParticle::GetDeviationFromParticle().

2627 {
2628  //* Calculate sqrt(Chi2/ndf) deviation from other particle
2629 
2630  Double_t dS, dS1;
2631  GetDStoParticle( p, dS, dS1 );
2632  Double_t mP1[8], mC1[36];
2633  p.Transport( dS1, mP1, mC1 );
2634 
2635  Double_t d[3]={ fP[0]-mP1[0], fP[1]-mP1[1], fP[2]-mP1[2]};
2636 
2637  Double_t sigmaS = .1+10.*sqrt( (d[0]*d[0]+d[1]*d[1]+d[2]*d[2])/
2638  (mP1[3]*mP1[3]+mP1[4]*mP1[4]+mP1[5]*mP1[5]) );
2639 
2640  Double_t h[3] = { mP1[3]*sigmaS, mP1[4]*sigmaS, mP1[5]*sigmaS };
2641 
2642  mC1[0] +=h[0]*h[0];
2643  mC1[1] +=h[1]*h[0];
2644  mC1[2] +=h[1]*h[1];
2645  mC1[3] +=h[2]*h[0];
2646  mC1[4] +=h[2]*h[1];
2647  mC1[5] +=h[2]*h[2];
2648 
2649  return GetDeviationFromVertex( mP1, mC1 )*sqrt(2./1.);
2650 }
TObjArray * d
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
virtual void GetDStoParticle(const KFParticleBase &p, Double_t &DS, Double_t &DSp) const =0
virtual void Transport(Double_t dS, Double_t P[], Double_t C[]) const =0
Double_t
Double_t GetDeviationFromVertex(const Double_t v[], const Double_t Cv[]=0) const
Double_t KFParticleBase::GetDeviationFromVertex ( const Double_t  v[],
const Double_t  Cv[] = 0 
) const

Definition at line 2575 of file KFParticleBase.cxx.

References d, Double_t, fabs(), GetDStoPoint(), h, mP, s, sqrt(), and Transport().

Referenced by GetDeviationFromParticle(), GetDeviationFromVertex(), and KFParticle::GetDeviationFromVertex().

2576 {
2577  //* Calculate sqrt(Chi2/ndf) deviation from vertex
2578  //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
2579 
2580  Double_t mP[8];
2581  Double_t mC[36];
2582 
2583  Transport( GetDStoPoint(v), mP, mC );
2584 
2585  Double_t d[3]={ v[0]-mP[0], v[1]-mP[1], v[2]-mP[2]};
2586 
2587  Double_t sigmaS = .1+10.*sqrt( (d[0]*d[0]+d[1]*d[1]+d[2]*d[2])/
2588  (mP[3]*mP[3]+mP[4]*mP[4]+mP[5]*mP[5]) );
2589 
2590 
2591  Double_t h[3] = { mP[3]*sigmaS, mP[4]*sigmaS, mP[5]*sigmaS };
2592 
2593  Double_t mSi[6] =
2594  { mC[0] +h[0]*h[0],
2595  mC[1] +h[1]*h[0], mC[2] +h[1]*h[1],
2596  mC[3] +h[2]*h[0], mC[4] +h[2]*h[1], mC[5] +h[2]*h[2] };
2597 
2598  if( Cv ){
2599  mSi[0]+=Cv[0];
2600  mSi[1]+=Cv[1];
2601  mSi[2]+=Cv[2];
2602  mSi[3]+=Cv[3];
2603  mSi[4]+=Cv[4];
2604  mSi[5]+=Cv[5];
2605  }
2606 
2607  Double_t mS[6];
2608 
2609  mS[0] = mSi[2]*mSi[5] - mSi[4]*mSi[4];
2610  mS[1] = mSi[3]*mSi[4] - mSi[1]*mSi[5];
2611  mS[2] = mSi[0]*mSi[5] - mSi[3]*mSi[3];
2612  mS[3] = mSi[1]*mSi[4] - mSi[2]*mSi[3];
2613  mS[4] = mSi[1]*mSi[3] - mSi[0]*mSi[4];
2614  mS[5] = mSi[0]*mSi[2] - mSi[1]*mSi[1];
2615 
2616  Double_t s = ( mSi[0]*mS[0] + mSi[1]*mS[1] + mSi[3]*mS[3] );
2617  s = ( s > 1.E-20 ) ?1./s :0;
2618 
2619  return sqrt( fabs(s*( ( mS[0]*d[0] + mS[1]*d[1] + mS[3]*d[2])*d[0]
2620  +(mS[1]*d[0] + mS[2]*d[1] + mS[4]*d[2])*d[1]
2621  +(mS[3]*d[0] + mS[4]*d[1] + mS[5]*d[2])*d[2] ))/2);
2622 }
TObjArray * d
double mP
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
TLorentzVector s
Definition: Pnd2DStar.C:50
__m128 v
Definition: P4_F32vec4.h:4
virtual void Transport(Double_t dS, Double_t P[], Double_t C[]) const =0
Double_t
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
virtual Double_t GetDStoPoint(const Double_t xyz[]) const =0
Double_t KFParticleBase::GetDeviationFromVertex ( const KFParticleBase Vtx) const

Definition at line 2567 of file KFParticleBase.cxx.

References fC, fP, and GetDeviationFromVertex().

2568 {
2569  //* Calculate sqrt(Chi2/ndf) deviation from vertex
2570 
2571  return GetDeviationFromVertex( Vtx.fP, Vtx.fC );
2572 }
Double_t fP[8]
Double_t fC[36]
Double_t GetDeviationFromVertex(const Double_t v[], const Double_t Cv[]=0) const
Double_t KFParticleBase::GetDistanceFromParticle ( const KFParticleBase p) const

Definition at line 2551 of file KFParticleBase.cxx.

References Double_t, dx, dy, dz, GetDStoParticle(), mP, sqrt(), and Transport().

Referenced by KFParticle::GetDistanceFromParticle().

2553 {
2554  //* Calculate distance to other particle [cm]
2555 
2556  Double_t dS, dS1;
2557  GetDStoParticle( p, dS, dS1 );
2558  Double_t mP[8], mC[36], mP1[8], mC1[36];
2559  Transport( dS, mP, mC );
2560  p.Transport( dS1, mP1, mC1 );
2561  Double_t dx = mP[0]-mP1[0];
2562  Double_t dy = mP[1]-mP1[1];
2563  Double_t dz = mP[2]-mP1[2];
2564  return sqrt(dx*dx+dy*dy+dz*dz);
2565 }
double dy
double mP
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
virtual void GetDStoParticle(const KFParticleBase &p, Double_t &DS, Double_t &DSp) const =0
virtual void Transport(Double_t dS, Double_t P[], Double_t C[]) const =0
Double_t
double dx
Double_t KFParticleBase::GetDistanceFromVertex ( const Double_t  vtx[]) const

Definition at line 2541 of file KFParticleBase.cxx.

References d, Double_t, GetDStoPoint(), mP, sqrt(), and Transport().

Referenced by GetDistanceFromVertex(), and KFParticle::GetDistanceFromVertex().

2542 {
2543  //* Calculate distance from vertex [cm]
2544 
2545  Double_t mP[8], mC[36];
2546  Transport( GetDStoPoint(vtx), mP, mC );
2547  Double_t d[3]={ vtx[0]-mP[0], vtx[1]-mP[1], vtx[2]-mP[2]};
2548  return sqrt( d[0]*d[0]+d[1]*d[1]+d[2]*d[2] );
2549 }
TObjArray * d
double mP
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
virtual void Transport(Double_t dS, Double_t P[], Double_t C[]) const =0
Double_t
virtual Double_t GetDStoPoint(const Double_t xyz[]) const =0
Double_t KFParticleBase::GetDistanceFromVertex ( const KFParticleBase Vtx) const

Definition at line 2534 of file KFParticleBase.cxx.

References fP, and GetDistanceFromVertex().

2535 {
2536  //* Calculate distance from vertex [cm]
2537 
2538  return GetDistanceFromVertex( Vtx.fP );
2539 }
Double_t fP[8]
Double_t GetDistanceFromVertex(const Double_t vtx[]) const
void KFParticleBase::GetDSIter ( const KFParticleBase p,
Double_t const &  dS,
Double_t  x[3],
Double_t  dx[3],
Double_t  ddx[3] 
) const
protected

Definition at line 2046 of file KFParticleBase.cxx.

References c, Double_t, fP, fQ, GetFieldValue(), m, n, p1, and p2.

2047 {
2048  Double_t ssx=0, ssy=0, ssz=0, ssyy=0, ssyz=0, ssyyy=0;
2049  Double_t dssx=0, dssy=0, dssz=0, dssyy=0, dssyz=0, dssyyy=0; // d(...)/dS
2050  Double_t ddssx=0, ddssy=0, ddssz=0, ddssyy=0, ddssyz=0, ddssyyy=0; // d2(...)/dS2
2051  Double_t Kssx=0, Kssy=0, Kssz=0, Kssyy=0, Kssyz=0, Kssyyy=0;
2052 
2053  const Double_t kCLight = 0.000299792458;
2054  Double_t c = fQ*kCLight;
2055 
2056  // get field integrals
2057 
2058  Double_t fld[3][3];
2059  Double_t p0[3], p1[3], p2[3];
2060 
2061  // line track approximation
2062 
2063  p0[0] = p.fP[0];
2064  p0[1] = p.fP[1];
2065  p0[2] = p.fP[2];
2066 
2067  p2[0] = p.fP[0] + p.fP[3]*dS;
2068  p2[1] = p.fP[1] + p.fP[4]*dS;
2069  p2[2] = p.fP[2] + p.fP[5]*dS;
2070 
2071  p1[0] = 0.5*(p0[0]+p2[0]);
2072  p1[1] = 0.5*(p0[1]+p2[1]);
2073  p1[2] = 0.5*(p0[2]+p2[2]);
2074 
2075  // first order track approximation
2076  {
2077  GetFieldValue( p0, fld[0] );
2078  GetFieldValue( p1, fld[1] );
2079  GetFieldValue( p2, fld[2] );
2080 
2081  Double_t ssy1 = ( 7*fld[0][1] + 6*fld[1][1]-fld[2][1] )*c*dS*dS/96.;
2082  Double_t ssy2 = ( fld[0][1] + 2*fld[1][1] )*c*dS*dS/6.;
2083 
2084  p1[0] -= ssy1*p.fP[5];
2085  p1[2] += ssy1*p.fP[3];
2086  p2[0] -= ssy2*p.fP[5];
2087  p2[2] += ssy2*p.fP[3];
2088  }
2089 
2090  GetFieldValue( p0, fld[0] );
2091  GetFieldValue( p1, fld[1] );
2092  GetFieldValue( p2, fld[2] );
2093 
2094  Kssx = c*( fld[0][0] + 2*fld[1][0]) / 6.;
2095  Kssy = c*( fld[0][1] + 2*fld[1][1]) / 6.;
2096  Kssz = c*( fld[0][2] + 2*fld[1][2]) / 6.;
2097 
2098 // Double_t c2[3][3] = { { 5, -4, -1},{ 44, 80, -4},{ 11, 44, 5} }; // /=360.
2099  Double_t cc2[3][3] = { { 38, 8, -4},{ 148, 208, -20},{ 3, 36, 3} }; // /=2520.
2100  for(Int_t n=0; n<3; n++)
2101  for(Int_t m=0; m<3; m++)
2102  {
2103  Kssyz += cc2[n][m]*fld[n][1]*fld[m][2];
2104  }
2105  Kssyz *= c*c / 2520.;
2106 
2107  Kssyy = ( fld[0][1]*( 38*fld[0][1] + 156*fld[1][1] - fld[2][1] )+
2108  fld[1][1]*( 208*fld[1][1] +16*fld[2][1] )+
2109  fld[2][1]*( 3*fld[2][1] )
2110  )*c*c/2520.;
2111  Kssyyy = ( fld[0][1]*( fld[0][1]*( 85*fld[0][1] + 526*fld[1][1] - 7*fld[2][1] )+
2112  fld[1][1]*( 1376*fld[1][1] +84*fld[2][1] )+
2113  fld[2][1]*( 19*fld[2][1] ) )+
2114  fld[1][1]*( fld[1][1]*( 1376*fld[1][1] +256*fld[2][1] )+
2115  fld[2][1]*( 62*fld[2][1] ) )+
2116  fld[2][1]*fld[2][1] *( 3*fld[2][1] )
2117  )*c*c*c/90720.;
2118 
2119  ssx = Kssx * dS*dS;
2120  ssy = Kssy * dS*dS;
2121  ssz = Kssz * dS*dS;
2122  ssyz = Kssyz * dS*dS*dS;
2123  ssyy = Kssyy * dS*dS*dS;
2124  ssyyy = Kssyyy * dS*dS*dS*dS;
2125 
2126  dssx = 2.*Kssx * dS;
2127  dssy = 2.*Kssy * dS;
2128  dssz = 2.*Kssz * dS;
2129  dssyz = 3.*Kssyz * dS*dS;
2130  dssyy = 3.*Kssyy * dS*dS;
2131  dssyyy = 4.*Kssyyy * dS*dS*dS;
2132 
2133  ddssx = 2.*Kssx;
2134  ddssy = 2.*Kssy;
2135  ddssz = 2.*Kssz;
2136  ddssyz = 2.*3.*Kssyz * dS;
2137  ddssyy = 2.*3.*Kssyy * dS;
2138  ddssyyy = 3.*4.*Kssyyy * dS*dS;
2139 
2140  Double_t mJ[3][3];
2141  for( Int_t iJ=0; iJ<3; iJ++ ) for( Int_t jJ=0; jJ<3; jJ++) mJ[iJ][jJ]=0;
2142 
2143  mJ[0][0]=dS-ssyy; mJ[0][1]=ssx; mJ[0][2]=ssyyy-ssy;
2144  mJ[1][0]=-ssz; mJ[1][1]=dS; mJ[1][2]=ssx+ssyz;
2145  mJ[2][0]=ssy-ssyyy; mJ[2][1]=-ssx; mJ[2][2]=dS-ssyy;
2146 
2147  x[0] = p.fP[0] + mJ[0][0]*p.fP[3] + mJ[0][1]*p.fP[4] + mJ[0][2]*p.fP[5];
2148  x[1] = p.fP[1] + mJ[1][0]*p.fP[3] + mJ[1][1]*p.fP[4] + mJ[1][2]*p.fP[5];
2149  x[2] = p.fP[2] + mJ[2][0]*p.fP[3] + mJ[2][1]*p.fP[4] + mJ[2][2]*p.fP[5];
2150 
2151  mJ[0][0]= 1-dssyy; mJ[0][1]= dssx; mJ[0][2]= dssyyy-dssy;
2152  mJ[1][0]= -dssz; mJ[1][1]= 1; mJ[1][2]= dssx+dssyz;
2153  mJ[2][0]= dssy-dssyyy; mJ[2][1]= -dssx; mJ[2][2]= 1-dssyy;
2154 
2155  dx[0] = p.fP[0] + mJ[0][0]*p.fP[3] + mJ[0][1]*p.fP[4] + mJ[0][2]*p.fP[5];
2156  dx[1] = p.fP[1] + mJ[1][0]*p.fP[3] + mJ[1][1]*p.fP[4] + mJ[1][2]*p.fP[5];
2157  dx[2] = p.fP[2] + mJ[2][0]*p.fP[3] + mJ[2][1]*p.fP[4] + mJ[2][2]*p.fP[5];
2158 
2159  mJ[0][0]= -ddssyy; mJ[0][1]= ddssx; mJ[0][2]= ddssyyy-ddssy;
2160  mJ[1][0]= -ddssz; mJ[1][1]= 0; mJ[1][2]= ddssx+ddssyz;
2161  mJ[2][0]= ddssy-ddssyyy; mJ[2][1]= -ddssx; mJ[2][2]= -ddssyy;
2162 
2163  ddx[0] = p.fP[0] + mJ[0][0]*p.fP[3] + mJ[0][1]*p.fP[4] + mJ[0][2]*p.fP[5];
2164  ddx[1] = p.fP[1] + mJ[1][0]*p.fP[3] + mJ[1][1]*p.fP[4] + mJ[1][2]*p.fP[5];
2165  ddx[2] = p.fP[2] + mJ[2][0]*p.fP[3] + mJ[2][1]*p.fP[4] + mJ[2][2]*p.fP[5];
2166 }
__m128 m
Definition: P4_F32vec4.h:28
Double_t fP[8]
int n
virtual void GetFieldValue(const Double_t xyz[], Double_t B[]) const =0
Double_t
TPad * p2
Definition: hist-t7.C:117
double dx
Double_t x
TPad * p1
Definition: hist-t7.C:116
virtual void KFParticleBase::GetDStoParticle ( const KFParticleBase p,
Double_t DS,
Double_t DSp 
) const
pure virtual
void KFParticleBase::GetDStoParticleBy ( Double_t  B,
const KFParticleBase p,
Double_t dS,
Double_t dS1 
) const

Definition at line 1909 of file KFParticleBase.cxx.

References a, atan2(), b, c, cos(), d, Double_t, dx, dy, fabs(), fP, fQ, g, g1, i, p2, pz, s, sin(), and sqrt().

Referenced by GetDStoParticleCBM().

1911 {
1912  //* Get dS to another particle for Bz field
1913  Double_t px = fP[3];
1914  Double_t py = -fP[5];
1915  Double_t pz = fP[4];
1916 
1917  Double_t px1 = p.fP[3];
1918  Double_t py1 = -p.fP[5];
1919  Double_t pz1 = p.fP[4];
1920 
1921  const Double_t kCLight = 0.000299792458;
1922 
1923  Double_t bq = B*fQ*kCLight;
1924  Double_t bq1 = B*p.fQ*kCLight;
1925  Double_t s=0, ds=0, s1=0, ds1=0;
1926 
1927  if( fabs(bq)>1.e-8 || fabs(bq1)>1.e-8 ){
1928 
1929  Double_t dx = (p.fP[0] - fP[0]);
1930  Double_t dy = (-p.fP[2] + fP[2]);
1931  Double_t d2 = (dx*dx+dy*dy);
1932 
1933  Double_t p2 = (px *px + py *py);
1934  Double_t p21 = (px1*px1 + py1*py1);
1935 
1936  if( fabs(p2) < 1.e-8 || fabs(p21) < 1.e-8 )
1937  {
1938  DS=0.;
1939  DS1=0.;
1940  return;
1941  }
1942 
1943  Double_t a = (px*py1 - py*px1);
1944  Double_t b = (px*px1 + py*py1);
1945 
1946  Double_t ldx = bq*bq1*dx - bq1*py + bq*py1 ;
1947  Double_t ldy = bq*bq1*dy + bq1*px - bq*px1 ;
1948  Double_t l2 = ldx*ldx + ldy*ldy;
1949 
1950  Double_t cS = bq1*p2 + bq*bq1*(dy* px - dx* py) - bq*b;
1951  Double_t cS1= bq*p21 - bq*bq1*(dy*px1 - dx*py1) - bq1*b;
1952 
1953  Double_t ca = bq*bq*bq1*d2 +2*( cS + bq*bq*(py1*dx-px1*dy)) ;
1954  Double_t ca1 = bq*bq1*bq1*d2 +2*( cS1 - bq1*bq1*(py*dx-px*dy)) ;
1955 
1956  Double_t sa = 4*l2*p2 - ca*ca;
1957  Double_t sa1 = 4*l2*p21 - ca1*ca1;
1958 
1959  if(sa<0) sa=0;
1960  if(sa1<0)sa1=0;
1961 
1962  if( fabs(bq)>1.e-8){
1963  s = atan2( bq*( bq1*(dx*px +dy*py) + a ) , cS )/bq;
1964  ds = atan2(sqrt(sa),ca)/bq;
1965  } else {
1966  s = ( (dx*px + dy*py) + (py*px1-px*py1)/bq1)/p2;
1967  ds = s*s - (d2-2*(px1*dy-py1*dx)/bq1)/p2;
1968  if( ds<0 ) ds = 0;
1969  ds = sqrt(ds);
1970  }
1971 
1972  if( fabs(bq1)>1.e-8){
1973  s1 = atan2( -bq1*( bq*(dx*px1+dy*py1) + a), cS1 )/bq1;
1974  ds1 = atan2(sqrt(sa1),ca1)/bq1;
1975  } else {
1976  s1 = (-(dx*px1 + dy*py1) + (py*px1-px*py1)/bq)/p21;
1977  ds1 = s1*s1 - (d2+2*(px*dy-py*dx)/bq)/p21;
1978  if( ds1<0 ) ds1 = 0;
1979  ds1 = sqrt(ds1);
1980  }
1981  }
1982  Double_t ss[2], ss1[2], g[2][5],g1[2][5];
1983 
1984  ss[0] = s + ds;
1985  ss[1] = s - ds;
1986  ss1[0] = s1 + ds1;
1987  ss1[1] = s1 - ds1;
1988 
1989  for( Int_t i=0; i<2; i++){
1990  Double_t bs = bq*ss[i];
1991  Double_t c = cos(bs), sss = sin(bs);
1992  Double_t cB,sB;
1993  if( fabs(bq)>1.e-8){
1994  cB= (1-c)/bq;
1995  sB= sss/bq;
1996  }else{
1997  const Double_t kOvSqr6 = 1./sqrt(6.);
1998  sB = (1.-bs*kOvSqr6)*(1.+bs*kOvSqr6)*ss[i];
1999  cB = .5*sB*bs;
2000  }
2001  g[i][0] = fP[0] + sB*px + cB*py;
2002  g[i][1] = -fP[2] - cB*px + sB*py;
2003  g[i][2] = fP[1] + ss[i]*pz;
2004  g[i][3] = + c*px + sss*py;
2005  g[i][4] = - sss*px + c*py;
2006 
2007  bs = bq1*ss1[i];
2008  c = cos(bs); sss = sin(bs);
2009  if( fabs(bq1)>1.e-8){
2010  cB= (1-c)/bq1;
2011  sB= sss/bq1;
2012  }else{
2013  const Double_t kOvSqr6 = 1./sqrt(6.);
2014  sB = (1.-bs*kOvSqr6)*(1.+bs*kOvSqr6)*ss1[i];
2015  cB = .5*sB*bs;
2016  }
2017 
2018  g1[i][0] = p.fP[0] + sB*px1 + cB*py1;
2019  g1[i][1] = -p.fP[2] - cB*px1 + sB*py1;
2020  g1[i][2] = p.fP[1] + ss1[i]*pz1;
2021  g1[i][3] = + c*px1 + sss*py1;
2022  g1[i][4] = - sss*px1 + c*py1;
2023  }
2024 
2025  Int_t i=0, i1=0;
2026 
2027  Double_t dMin = 1.e10;
2028  for( Int_t j=0; j<2; j++){
2029  for( Int_t j1=0; j1<2; j1++){
2030  Double_t xx = g[j][0]-g1[j1][0];
2031  Double_t yy = g[j][1]-g1[j1][1];
2032  Double_t zz = g[j][2]-g1[j1][2];
2033  Double_t d = xx*xx + yy*yy + zz*zz;
2034  if( d<dMin ){
2035  dMin = d;
2036  i = j;
2037  i1 = j1;
2038  }
2039  }
2040  }
2041 
2042  DS = ss[i];
2043  DS1 = ss1[i1];
2044 }
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
double dy
TObjArray * d
Int_t i
Definition: run_full.C:25
TTree * b
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
TLorentzVector s
Definition: Pnd2DStar.C:50
TFile * g
Int_t a
Definition: anaLmdDigi.C:126
Double_t
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
TPad * p2
Definition: hist-t7.C:117
double dx
friend F32vec4 atan2(const F32vec4 &y, const F32vec4 &x)
Definition: P4_F32vec4.h:117
TF1 * g1
Definition: anaLmdReco.C:158
double pz[39]
Definition: pipisigmas.h:14
void KFParticleBase::GetDStoParticleBz ( Double_t  Bz,
const KFParticleBase p,
Double_t dS,
Double_t dS1 
) const

Definition at line 1754 of file KFParticleBase.cxx.

References a, atan2(), b, c, cos(), d, Double_t, dx, dy, dz, fabs(), fP, fQ, g, g1, i, p2, pz, s, sin(), sqrt(), x, y, and z.

Referenced by KFParticle::GetDStoParticleXY().

1757 {
1758  //* Get dS to another particle for Bz field
1759  Double_t px = fP[3];
1760  Double_t py = fP[4];
1761  Double_t pz = fP[5];
1762 
1763  Double_t px1 = p.fP[3];
1764  Double_t py1 = p.fP[4];
1765  Double_t pz1 = p.fP[5];
1766 
1767  const Double_t kCLight = 0.000299792458;
1768 
1769  Double_t bq = B*fQ*kCLight;
1770  Double_t bq1 = B*p.fQ*kCLight;
1771  Double_t s=0, ds=0, s1=0, ds1=0;
1772 
1773  if( fabs(bq)>1.e-8 || fabs(bq1)>1.e-8 ){
1774 
1775  Double_t dx = (p.fP[0] - fP[0]);
1776  Double_t dy = (p.fP[1] - fP[1]);
1777  Double_t d2 = (dx*dx+dy*dy);
1778 
1779  Double_t p2 = (px *px + py *py);
1780  Double_t p21 = (px1*px1 + py1*py1);
1781 
1782  if( fabs(p2) < 1.e-8 || fabs(p21) < 1.e-8 )
1783  {
1784  DS=0.;
1785  DS1=0.;
1786  return;
1787  }
1788 
1789  Double_t a = (px*py1 - py*px1);
1790  Double_t b = (px*px1 + py*py1);
1791 
1792  Double_t ldx = bq*bq1*dx - bq1*py + bq*py1 ;
1793  Double_t ldy = bq*bq1*dy + bq1*px - bq*px1 ;
1794  Double_t l2 = ldx*ldx + ldy*ldy;
1795 
1796  Double_t cS = bq1*p2 + bq*bq1*(dy* px - dx* py) - bq*b;
1797  Double_t cS1= bq*p21 - bq*bq1*(dy*px1 - dx*py1) - bq1*b;
1798 
1799  Double_t ca = bq*bq*bq1*d2 +2*( cS + bq*bq*(py1*dx-px1*dy)) ;
1800  Double_t ca1 = bq*bq1*bq1*d2 +2*( cS1 - bq1*bq1*(py*dx-px*dy)) ;
1801 
1802  Double_t sa = 4*l2*p2 - ca*ca;
1803  Double_t sa1 = 4*l2*p21 - ca1*ca1;
1804 
1805  if(sa<0) sa=0;
1806  if(sa1<0)sa1=0;
1807 
1808  if( fabs(bq)>1.e-8){
1809  s = atan2( bq*( bq1*(dx*px +dy*py) + a ) , cS )/bq;
1810  ds = atan2(sqrt(sa),ca)/bq;
1811  } else {
1812  s = ( (dx*px + dy*py) + (py*px1-px*py1)/bq1)/p2;
1813  ds = s*s - (d2-2*(px1*dy-py1*dx)/bq1)/p2;
1814  if( ds<0 ) ds = 0;
1815  ds = sqrt(ds);
1816  }
1817 
1818  if( fabs(bq1)>1.e-8){
1819  s1 = atan2( -bq1*( bq*(dx*px1+dy*py1) + a), cS1 )/bq1;
1820  ds1 = atan2(sqrt(sa1),ca1)/bq1;
1821  } else {
1822  s1 = (-(dx*px1 + dy*py1) + (py*px1-px*py1)/bq)/p21;
1823  ds1 = s1*s1 - (d2+2*(px*dy-py*dx)/bq)/p21;
1824  if( ds1<0 ) ds1 = 0;
1825  ds1 = sqrt(ds1);
1826  }
1827  }
1828  Double_t ss[2], ss1[2], g[2][5],g1[2][5];
1829 
1830  ss[0] = s + ds;
1831  ss[1] = s - ds;
1832  ss1[0] = s1 + ds1;
1833  ss1[1] = s1 - ds1;
1834 
1835  for( Int_t i=0; i<2; i++){
1836  Double_t bs = bq*ss[i];
1837  Double_t c = cos(bs), sss = sin(bs);
1838  Double_t cB,sB;
1839  if( fabs(bq)>1.e-8){
1840  cB= (1-c)/bq;
1841  sB= sss/bq;
1842  }else{
1843  const Double_t kOvSqr6 = 1./sqrt(6.);
1844  sB = (1.-bs*kOvSqr6)*(1.+bs*kOvSqr6)*ss[i];
1845  cB = .5*sB*bs;
1846  }
1847  g[i][0] = fP[0] + sB*px + cB*py;
1848  g[i][1] = fP[1] - cB*px + sB*py;
1849  g[i][2] = fP[2] + ss[i]*pz;
1850  g[i][3] = + c*px + sss*py;
1851  g[i][4] = - sss*px + c*py;
1852 
1853  bs = bq1*ss1[i];
1854  c = cos(bs); sss = sin(bs);
1855  if( fabs(bq1)>1.e-8){
1856  cB= (1-c)/bq1;
1857  sB= sss/bq1;
1858  }else{
1859  const Double_t kOvSqr6 = 1./sqrt(6.);
1860  sB = (1.-bs*kOvSqr6)*(1.+bs*kOvSqr6)*ss1[i];
1861  cB = .5*sB*bs;
1862  }
1863 
1864  g1[i][0] = p.fP[0] + sB*px1 + cB*py1;
1865  g1[i][1] = p.fP[1] - cB*px1 + sB*py1;
1866  g1[i][2] = p.fP[2] + ss1[i]*pz1;
1867  g1[i][3] = + c*px1 + sss*py1;
1868  g1[i][4] = - sss*px1 + c*py1;
1869  }
1870 
1871  Int_t i=0, i1=0;
1872 
1873  Double_t dMin = 1.e10;
1874  for( Int_t j=0; j<2; j++){
1875  for( Int_t j1=0; j1<2; j1++){
1876  Double_t xx = g[j][0]-g1[j1][0];
1877  Double_t yy = g[j][1]-g1[j1][1];
1878  Double_t zz = g[j][2]-g1[j1][2];
1879  Double_t d = xx*xx + yy*yy + zz*zz;
1880  if( d<dMin ){
1881  dMin = d;
1882  i = j;
1883  i1 = j1;
1884  }
1885  }
1886  }
1887 
1888  DS = ss[i];
1889  DS1 = ss1[i1];
1890  if(0){
1891  Double_t x= g[i][0], y= g[i][1], z= g[i][2], ppx= g[i][3], ppy= g[i][4];
1892  Double_t x1=g1[i1][0], y1= g1[i1][1], z1= g1[i1][2], ppx1= g1[i1][3], ppy1= g1[i1][4];
1893  Double_t dx = x1-x;
1894  Double_t dy = y1-y;
1895  Double_t dz = z1-z;
1896  Double_t a = ppx*ppx1 + ppy*ppy1 + pz*pz1;
1897  Double_t b = dx*ppx1 + dy*ppy1 + dz*pz1;
1898  Double_t c = dx*ppx + dy*ppy + dz*pz ;
1899  Double_t pp2 = ppx*ppx + ppy*ppy + pz*pz;
1900  Double_t pp21= ppx1*ppx1 + ppy1*ppy1 + pz1*pz1;
1901  Double_t det = pp2*pp21 - a*a;
1902  if( fabs(det)>1.e-8 ){
1903  DS+=(a*b-pp21*c)/det;
1904  DS1+=(a*c-pp2*b)/det;
1905  }
1906  }
1907 }
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
double dy
TObjArray * d
Int_t i
Definition: run_full.C:25
TTree * b
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
TLorentzVector s
Definition: Pnd2DStar.C:50
TFile * g
Int_t a
Definition: anaLmdDigi.C:126
Double_t
Double_t z
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
TPad * p2
Definition: hist-t7.C:117
double dx
friend F32vec4 atan2(const F32vec4 &y, const F32vec4 &x)
Definition: P4_F32vec4.h:117
Double_t x
TF1 * g1
Definition: anaLmdReco.C:158
Double_t y
double pz[39]
Definition: pipisigmas.h:14
void KFParticleBase::GetDStoParticleCBM ( const KFParticleBase p,
Double_t dS,
Double_t dS1 
) const

Definition at line 2230 of file KFParticleBase.cxx.

References Double_t, fP, fQ, GetDStoParticleBy(), GetDStoParticleLine(), and GetFieldValue().

Referenced by KFParticle::GetDStoParticleXY().

2231 {
2232  //* Transport the particle on dS, output to P[],C[], for CBM field
2233 
2234  GetDStoParticleLine( p, dS, dS1 );
2235  if( fQ==0 ){
2236  return;
2237  }
2238 
2239  Double_t fld[3];
2240  GetFieldValue( fP, fld );
2241 
2242  Double_t fld1[3];
2243  GetFieldValue( p.fP, fld1 );
2244 
2245 
2246  GetDStoParticleBy((fld[1]+fld1[1])/2,p,dS,dS1);
2247 
2248 // construct coefficients
2249 // int NIt;
2250 // for(int i=0; i<10; i++)
2251 // {
2252 // /*std::cout << " dS temp "<<dS << " dS1 temp "<< dS1 << std::endl;*/
2253 // Double_t x[3], dx[3], ddx[3];
2254 // Double_t x1[3], dx1[3], ddx1[3];
2255 //
2256 // GetDSIter(*this, dS, x, dx, ddx);
2257 // GetDSIter(p, dS1, x1, dx1, ddx1);
2258 //
2259 // Double_t dr[3] = { x[0] - x1[0], x[1] - x1[1], x[2] - x1[2] };
2260 //
2261 // Double_t dS0 = dS;
2262 // Double_t dS10 = dS1;
2263 //
2264 // Double_t f = dx[0] * dr[0] + dx[1] * dr[1] + dx[2] * dr[2];
2265 // Double_t dfs = ddx[0] * dr[0] + ddx[1] * dr[1] + ddx[2] * dr[2] + dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2];
2266 // Double_t dfs1 = -ddx[0] * ddx1[0] - ddx[1] * ddx1[1] - ddx[2] * ddx1[2];
2267 //
2268 // Double_t f1 = dx1[0] * dr[0] + dx1[1] * dr[1] + dx1[2] * dr[2];
2269 // Double_t df1s1 = ddx1[0] * dr[0] + ddx1[1] * dr[1] + ddx1[2] * dr[2] - dx1[0]*dx1[0] - dx1[1]*dx1[1] - dx1[2]*dx1[2];
2270 // Double_t df1s = -dfs1;
2271 //
2272 // Double_t det = dfs*df1s1 - dfs1*df1s;
2273 // if(fabs(det) < 1.e-8) det = 1.e-8;
2274 //
2275 // Double_t DS = dfs1*f1 - f*df1s1;
2276 // Double_t DS1 = f*df1s - f1*dfs;
2277 //
2278 // dS += DS/det;
2279 // dS1 += DS1/det;
2280 // NIt++;
2281 // /* if( (fabs(dS - dS0) < fabs(dS)*0.001) && (fabs(dS1 - dS10) < fabs(dS1)*0.001)) break;*/
2282 // }
2283 //std::cout << "NIt " << NIt <<" dS "<< dS << " dS1 "<< dS1 << std::endl;
2284 }
Double_t fP[8]
void GetDStoParticleBy(Double_t B, const KFParticleBase &p, Double_t &dS, Double_t &dS1) const
void GetDStoParticleLine(const KFParticleBase &p, Double_t &dS, Double_t &dS1) const
virtual void GetFieldValue(const Double_t xyz[], Double_t B[]) const =0
Double_t
void KFParticleBase::GetDStoParticleLine ( const KFParticleBase p,
Double_t dS,
Double_t dS1 
) const
protected

Definition at line 2214 of file KFParticleBase.cxx.

References Double_t, fabs(), and fP.

Referenced by GetDStoParticleCBM().

2215 {
2216  Double_t p12 = fP[3]*fP[3] + fP[4]*fP[4] + fP[5]*fP[5];
2217  Double_t p22 = p.fP[3]*p.fP[3] + p.fP[4]*p.fP[4] + p.fP[5]*p.fP[5];
2218  Double_t p1p2 = fP[3]*p.fP[3] + fP[4]*p.fP[4] + fP[5]*p.fP[5];
2219 
2220  Double_t drp1 = fP[3]*(p.fP[0]-fP[0]) + fP[4]*(p.fP[1]-fP[1]) + fP[5]*(p.fP[2]-fP[2]);
2221  Double_t drp2 = p.fP[3]*(p.fP[0]-fP[0]) + p.fP[4]*(p.fP[1]-fP[1]) + p.fP[5]*(p.fP[2]-fP[2]);
2222 
2223  Double_t detp = p1p2*p1p2 - p12*p22;
2224  if( fabs(detp)<1.e-4 ) detp = 1; //TODO correct!!!
2225 
2226  dS = (drp2*p1p2 - drp1*p22) /detp;
2227  dS1 = (drp2*p12 - drp1*p1p2)/detp;
2228 }
Double_t fP[8]
Double_t
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
virtual Double_t KFParticleBase::GetDStoPoint ( const Double_t  xyz[]) const
pure virtual
Double_t KFParticleBase::GetDStoPointBy ( Double_t  By,
const Double_t  xyz[] 
) const

Definition at line 1734 of file KFParticleBase.cxx.

References a, atan2(), Double_t, dx, dy, fabs(), fP, and fQ.

Referenced by GetDStoPointCBM().

1736 {
1737 
1738  //* Get dS to a certain space point for Bz field
1739  const Double_t kCLight = 0.000299792458;
1740  Double_t bq = B*fQ*kCLight;
1741  Double_t pt2 = fP[3]*fP[3] + fP[5]*fP[5];
1742  if( pt2<1.e-4 ) return 0;
1743  Double_t dx = xyz[0] - fP[0];
1744  Double_t dy = -xyz[2] + fP[2];
1745  Double_t a = dx*fP[3]-dy*fP[5];
1746  Double_t dS;
1747 
1748  if( fabs(bq)<1.e-8 ) dS = a/pt2;
1749  else dS = atan2( bq*a, pt2 + bq*(dy*fP[3] +dx*fP[5]) )/bq;
1750 
1751  return dS;
1752 }
double dy
Double_t fP[8]
Int_t a
Definition: anaLmdDigi.C:126
Double_t
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
double dx
friend F32vec4 atan2(const F32vec4 &y, const F32vec4 &x)
Definition: P4_F32vec4.h:117
Double_t KFParticleBase::GetDStoPointBz ( Double_t  Bz,
const Double_t  xyz[] 
) const

Definition at line 1660 of file KFParticleBase.cxx.

References a, atan2(), c, cos(), d, Double_t, dx, dy, fabs(), fP, fQ, g, i, pz, s, sin(), sqrt(), x, y, and z.

Referenced by KFParticle::GetDStoPoint().

1662 {
1663 
1664  //* Get dS to a certain space point for Bz field
1665  const Double_t kCLight = 0.000299792458;
1666  Double_t bq = B*fQ*kCLight;
1667  Double_t pt2 = fP[3]*fP[3] + fP[4]*fP[4];
1668  if( pt2<1.e-4 ) return 0;
1669  Double_t dx = xyz[0] - fP[0];
1670  Double_t dy = xyz[1] - fP[1];
1671  Double_t a = dx*fP[3]+dy*fP[4];
1672  Double_t dS;
1673 
1674  if( fabs(bq)<1.e-8 ) dS = a/pt2;
1675  else dS = atan2( bq*a, pt2 + bq*(dy*fP[3] -dx*fP[4]) )/bq;
1676 
1677  if(0){
1678 
1679  Double_t px = fP[3];
1680  Double_t py = fP[4];
1681  Double_t pz = fP[5];
1682  Double_t ss[2], g[2][5];
1683 
1684  ss[0] = dS;
1685  ss[1] = -dS;
1686  for( Int_t i=0; i<2; i++){
1687  Double_t bs = bq*ss[i];
1688  Double_t c = cos(bs), s = sin(bs);
1689  Double_t cB,sB;
1690  if( fabs(bq)>1.e-8){
1691  cB= (1-c)/bq;
1692  sB= s/bq;
1693  }else{
1694  const Double_t kOvSqr6 = 1./sqrt(6.);
1695  sB = (1.-bs*kOvSqr6)*(1.+bs*kOvSqr6)*ss[i];
1696  cB = .5*sB*bs;
1697  }
1698  g[i][0] = fP[0] + sB*px + cB*py;
1699  g[i][1] = fP[1] - cB*px + sB*py;
1700  g[i][2] = fP[2] + ss[i]*pz;
1701  g[i][3] = + c*px + s*py;
1702  g[i][4] = - s*px + c*py;
1703  }
1704 
1705  Int_t i=0;
1706 
1707  Double_t dMin = 1.e10;
1708  for( Int_t j=0; j<2; j++){
1709  Double_t xx = g[j][0]-xyz[0];
1710  Double_t yy = g[j][1]-xyz[1];
1711  Double_t zz = g[j][2]-xyz[2];
1712  Double_t d = xx*xx + yy*yy + zz*zz;
1713  if( d<dMin ){
1714  dMin = d;
1715  i = j;
1716  }
1717  }
1718 
1719  dS = ss[i];
1720 
1721  Double_t x= g[i][0], y= g[i][1], z= g[i][2], ppx= g[i][3], ppy= g[i][4];
1722  Double_t ddx = x-xyz[0];
1723  Double_t ddy = y-xyz[1];
1724  Double_t ddz = z-xyz[2];
1725  Double_t c = ddx*ppx + ddy*ppy + ddz*pz ;
1726  Double_t pp2 = ppx*ppx + ppy*ppy + pz*pz;
1727  if( fabs(pp2)>1.e-8 ){
1728  dS+=c/pp2;
1729  }
1730  }
1731  return dS;
1732 }
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
double dy
TObjArray * d
Int_t i
Definition: run_full.C:25
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
TLorentzVector s
Definition: Pnd2DStar.C:50
TFile * g
Int_t a
Definition: anaLmdDigi.C:126
Double_t
Double_t z
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
double dx
friend F32vec4 atan2(const F32vec4 &y, const F32vec4 &x)
Definition: P4_F32vec4.h:117
Double_t x
Double_t y
double pz[39]
Definition: pipisigmas.h:14
Double_t KFParticleBase::GetDStoPointCBM ( const Double_t  xyz[]) const

Definition at line 2168 of file KFParticleBase.cxx.

References Double_t, fP, fQ, GetDStoPointBy(), GetDStoPointLine(), and GetFieldValue().

Referenced by KFParticle::GetDStoPoint().

2169 {
2170  //* Transport the particle on dS, output to P[],C[], for CBM field
2171 
2172  Double_t dS = 0;
2173  //Linear approximation
2174 
2175  Double_t fld[3];
2176  GetFieldValue( fP, fld );
2177 //
2178 // GetDStoParticleBy(fld[1],p,dS,dS1);
2179  dS = GetDStoPointLine( xyz );
2180 
2181  if( fQ==0 ){
2182  return dS;
2183  }
2184 
2185  dS = GetDStoPointBy( fld[1],xyz );
2186 
2187 // const Double_t &px = fP[3], &py = fP[4], &pz = fP[5];
2188 
2189 
2190  // construct coefficients
2191 // int NIt;
2192 // for(int i=0; i<100; i++)
2193 // {
2194 // //std::cout << " dS temp "<<dS << std::endl;
2195 // Double_t x[3], dx[3], ddx[3];
2196 //
2197 // GetDSIter(*this, dS, x, dx, ddx);
2198 //
2199 // Double_t dr[3] = { x[0] - xyz[0], x[1] - xyz[1], x[2] - xyz[2] };
2200 //
2201 // Double_t dS0 = dS;
2202 //
2203 // Double_t f = dx[0] * dr[0] + dx[1] * dr[1] + dx[2] * dr[2];
2204 // Double_t df = ddx[0] * dr[0] + ddx[1] * dr[1] + ddx[2] * dr[2] + dx[0]*dx[0] + dx[1]*dx[1] + dx[2]*dx[2];
2205 //
2206 // dS = dS - f/df;
2207 // NIt++;
2208 // if(fabs(dS - dS0) < fabs(dS)*0.001) break;
2209 // }
2210 //std::cout << "NIt " << NIt <<" dS "<< dS << std::endl;
2211  return dS;
2212 }
Double_t GetDStoPointBy(Double_t By, const Double_t xyz[]) const
Double_t fP[8]
virtual void GetFieldValue(const Double_t xyz[], Double_t B[]) const =0
Double_t
Double_t GetDStoPointLine(const Double_t xyz[]) const
Double_t KFParticleBase::GetDStoPointLine ( const Double_t  xyz[]) const
protected

Definition at line 1650 of file KFParticleBase.cxx.

References Double_t, fP, and p2.

Referenced by GetDStoPointCBM().

1651 {
1652  //* Get dS to a certain space point without field
1653 
1654  Double_t p2 = fP[3]*fP[3] + fP[4]*fP[4] + fP[5]*fP[5];
1655  if( p2<1.e-4 ) p2 = 1;
1656  return ( fP[3]*(xyz[0]-fP[0]) + fP[4]*(xyz[1]-fP[1]) + fP[5]*(xyz[2]-fP[2]) )/p2;
1657 }
Double_t fP[8]
Double_t
TPad * p2
Definition: hist-t7.C:117
Double_t KFParticleBase::GetE ( ) const
inline

Definition at line 107 of file KFParticleBase.h.

References fP.

Referenced by KFParticle::GetE().

107 { return fP[6]; }
Double_t fP[8]
Int_t KFParticleBase::GetEta ( Double_t Eta,
Double_t SigmaEta 
) const

Definition at line 150 of file KFParticleBase.cxx.

References a, b, c, Double_t, fC, fP, h3, h4, log(), p, p2, pz, and sqrt().

Referenced by KFParticle::GetErrEta(), and KFParticle::GetEta().

151 {
152  //* Calculate particle pseudorapidity
153 
154  Double_t px = fP[3];
155  Double_t py = fP[4];
156  Double_t pz = fP[5];
157  Double_t pt2 = px*px + py*py;
158  Double_t p2 = pt2 + pz*pz;
159  Double_t p = sqrt(p2);
160  Double_t a = p + pz;
161  Double_t b = p - pz;
162  eta = 1.e10;
163  if( b > 1.e-8 ){
164  Double_t c = a/b;
165  if( c>1.e-8 ) eta = 0.5*log(c);
166  }
167  Double_t h3 = -px*pz;
168  Double_t h4 = -py*pz;
169  Double_t pt4 = pt2*pt2;
170  Double_t p2pt4 = p2*pt4;
171  error = (h3*h3*fC[9] + h4*h4*fC[14] + pt4*fC[20] + 2*( h3*(h4*fC[13] + fC[18]*pt2) + pt2*h4*fC[19] ) );
172 
173  if( error>0 && p2pt4>1.e-10 ){
174  error = sqrt(error/p2pt4);
175  return 0;
176  }
177 
178  error = 1.e10;
179  return 1;
180 }
TTree * b
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
TH1F * h4
friend F32vec4 log(const F32vec4 &a)
Definition: P4_F32vec4.h:110
static void error(int no)
Definition: ranlxd.cxx:419
Double_t p
Definition: anasim.C:58
Int_t a
Definition: anaLmdDigi.C:126
Double_t
Double_t fC[36]
TH1F * h3
TPad * p2
Definition: hist-t7.C:117
TParticlePDG * eta
double pz[39]
Definition: pipisigmas.h:14
virtual void KFParticleBase::GetFieldValue ( const Double_t  xyz[],
Double_t  B[] 
) const
pure virtual
Int_t KFParticleBase::GetLifeTime ( Double_t T,
Double_t SigmaT 
) const

Definition at line 311 of file KFParticleBase.cxx.

References Double_t, fC, fP, GetMass(), m, and sqrt().

Referenced by KFParticle::GetErrLifeTime(), and KFParticle::GetLifeTime().

312 {
313  //* Calculate particle decay time [s]
314 
315  Double_t m, dm;
316  GetMass( m, dm );
317  Double_t cTM = (-fP[3]*fC[31] - fP[4]*fC[32] - fP[5]*fC[33] + fP[6]*fC[34]);
318  tauC = fP[7]*m;
319  error = m*m*fC[35] + 2*fP[7]*cTM + fP[7]*fP[7]*dm*dm;
320  if( error > 0 ){
321  error = sqrt( error );
322  return 0;
323  }
324  error = 1.e20;
325  return 1;
326 }
__m128 m
Definition: P4_F32vec4.h:28
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
Int_t GetMass(Double_t &M, Double_t &SigmaM) const
static void error(int no)
Definition: ranlxd.cxx:419
Double_t
Double_t fC[36]
Int_t KFParticleBase::GetMass ( Double_t M,
Double_t SigmaM 
) const

Definition at line 219 of file KFParticleBase.cxx.

References Double_t, fC, fP, m, m2(), s, and sqrt().

Referenced by KFParticle::GetErrMass(), GetLifeTime(), and KFParticle::GetMass().

220 {
221  //* Calculate particle mass
222 
223  // s = sigma^2 of m2/2
224 
225  Double_t s = ( fP[3]*fP[3]*fC[9] + fP[4]*fP[4]*fC[14] + fP[5]*fP[5]*fC[20]
226  + fP[6]*fP[6]*fC[27]
227  +2*( + fP[3]*fP[4]*fC[13] + fP[5]*(fP[3]*fC[18] + fP[4]*fC[19])
228  - fP[6]*( fP[3]*fC[24] + fP[4]*fC[25] + fP[5]*fC[26] ) )
229  );
230 // Double_t m2 = fabs(fP[6]*fP[6] - fP[3]*fP[3] - fP[4]*fP[4] - fP[5]*fP[5]);
231 // m = sqrt(m2);
232 // if( m>1.e-10 ){
233 // if( s>=0 ){
234 // error = sqrt(s)/m;
235 // return 0;
236 // }
237 // }
238 // error = 1.e20;
239  Double_t m2 = (fP[6]*fP[6] - fP[3]*fP[3] - fP[4]*fP[4] - fP[5]*fP[5]);
240 
241  if(m2<0.)
242  {
243  error = 1.e20;
244  m = -sqrt(-m2);
245  return 1;
246  }
247 
248  m = sqrt(m2);
249  if( m>1.e-6 ){
250  if( s >= 0 ) {
251  error = sqrt(s)/m;
252  return 0;
253  }
254  }
255  else {
256  error = 1.e20;
257  return 0;
258  }
259  error = 1.e20;
260 
261  return 1;
262 }
__m128 m
Definition: P4_F32vec4.h:28
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
TLorentzVector s
Definition: Pnd2DStar.C:50
static void error(int no)
Definition: ranlxd.cxx:419
TString m2(TString pts, TString exts="e px py pz")
Definition: invexp.C:117
Double_t
Double_t fC[36]
const Double_t& KFParticleBase::GetMassHypo ( ) const
inline

Definition at line 90 of file KFParticleBase.h.

References fMassHypo.

90 { return fMassHypo; }
Double_t fMassHypo
void KFParticleBase::GetMeasurement ( const Double_t  XYZ[],
Double_t  m[],
Double_t  V[] 
) const
protected

Definition at line 347 of file KFParticleBase.cxx.

References b, Double_t, GetDStoPoint(), GetFieldValue(), GetQ(), GetSCorrection(), h, and Transport().

Referenced by AddDaughterWithEnergyCalc(), AddDaughterWithEnergyFit(), AddDaughterWithEnergyFitMC(), ConstructGammaBz(), and SubtractFromVertex().

348 {
349  //* Get additional covariances V used during measurement
350 
351  Double_t b[3];
352  GetFieldValue( XYZ, b );
353  const Double_t kCLight = 0.000299792458;
354  b[0]*=kCLight; b[1]*=kCLight; b[2]*=kCLight;
355 // std::cout << " DStoPoint "<< (GetDStoPoint(XYZ)) << std::endl;
356  Transport( GetDStoPoint(XYZ), m, V );
357 //std::cout << "x " << XYZ[0] << " " << m[0] <<" y " << XYZ[1] << " " << m[1] <<" z " << XYZ[2] << " " << m[2] <<std::endl;
358  Double_t sigmaS = GetSCorrection( m, XYZ );
359 
360  Double_t h[6];
361 
362  h[0] = m[3]*sigmaS;
363  h[1] = m[4]*sigmaS;
364  h[2] = m[5]*sigmaS;
365  h[3] = ( h[1]*b[2]-h[2]*b[1] )*GetQ();
366  h[4] = ( h[2]*b[0]-h[0]*b[2] )*GetQ();
367  h[5] = ( h[0]*b[1]-h[1]*b[0] )*GetQ();
368 
369  V[ 0]+= h[0]*h[0];
370  V[ 1]+= h[1]*h[0];
371  V[ 2]+= h[1]*h[1];
372  V[ 3]+= h[2]*h[0];
373  V[ 4]+= h[2]*h[1];
374  V[ 5]+= h[2]*h[2];
375 
376  V[ 6]+= h[3]*h[0];
377  V[ 7]+= h[3]*h[1];
378  V[ 8]+= h[3]*h[2];
379  V[ 9]+= h[3]*h[3];
380 
381  V[10]+= h[4]*h[0];
382  V[11]+= h[4]*h[1];
383  V[12]+= h[4]*h[2];
384  V[13]+= h[4]*h[3];
385  V[14]+= h[4]*h[4];
386 
387  V[15]+= h[5]*h[0];
388  V[16]+= h[5]*h[1];
389  V[17]+= h[5]*h[2];
390  V[18]+= h[5]*h[3];
391  V[19]+= h[5]*h[4];
392  V[20]+= h[5]*h[5];
393 }
__m128 m
Definition: P4_F32vec4.h:28
TTree * b
virtual void Transport(Double_t dS, Double_t P[], Double_t C[]) const =0
virtual void GetFieldValue(const Double_t xyz[], Double_t B[]) const =0
Double_t
static Double_t GetSCorrection(const Double_t Part[], const Double_t XYZ[])
Int_t GetQ() const
virtual Double_t GetDStoPoint(const Double_t xyz[]) const =0
Int_t KFParticleBase::GetMomentum ( Double_t P,
Double_t SigmaP 
) const

Definition at line 110 of file KFParticleBase.cxx.

References Double_t, fC, fP, p, p2, sqrt(), x, y, and z.

Referenced by KFParticle::GetErrMomentum(), KFParticle::GetErrP(), KFParticle::GetMomentum(), and KFParticle::GetP().

111 {
112  //* Calculate particle momentum
113 
114  Double_t x = fP[3];
115  Double_t y = fP[4];
116  Double_t z = fP[5];
117  Double_t x2 = x*x;
118  Double_t y2 = y*y;
119  Double_t z2 = z*z;
120  Double_t p2 = x2+y2+z2;
121  p = sqrt(p2);
122  error = (x2*fC[9]+y2*fC[14]+z2*fC[20] + 2*(x*y*fC[13]+x*z*fC[18]+y*z*fC[19]) );
123  if( error>1.e-16 && p>1.e-4 ){
124  error = sqrt(error)/p;
125  return 0;
126  }
127  error = 1.e8;
128  return 1;
129 }
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
static void error(int no)
Definition: ranlxd.cxx:419
Double_t p
Definition: anasim.C:58
Double_t
Double_t z
Double_t fC[36]
TPad * p2
Definition: hist-t7.C:117
Double_t x
Double_t y
Int_t KFParticleBase::GetNDF ( ) const
inline

Definition at line 111 of file KFParticleBase.h.

References fNDF.

Referenced by KFParticle::GetNDF().

111 { return fNDF; }
Double_t KFParticleBase::GetParameter ( Int_t  i) const
inline

Definition at line 125 of file KFParticleBase.h.

References fP, and i.

Referenced by KFParticle::GetParameter().

125 { return fP[i]; }
Int_t i
Definition: run_full.C:25
Double_t fP[8]
int KFParticleBase::GetPDG ( ) const
inline

Definition at line 277 of file KFParticleBase.h.

References fPDG.

Referenced by KFParticleSIMD::KFParticleSIMD().

277 { return fPDG; }
Int_t KFParticleBase::GetPhi ( Double_t Phi,
Double_t SigmaPhi 
) const

Definition at line 182 of file KFParticleBase.cxx.

References atan2(), Double_t, fC, fP, and sqrt().

Referenced by KFParticle::GetErrPhi(), and KFParticle::GetPhi().

183 {
184  //* Calculate particle polar angle
185 
186  Double_t px = fP[3];
187  Double_t py = fP[4];
188  Double_t px2 = px*px;
189  Double_t py2 = py*py;
190  Double_t pt2 = px2 + py2;
191  phi = atan2(py,px);
192  error = (py2*fC[9] + px2*fC[14] - 2*px*py*fC[13] );
193  if( error>0 && pt2>1.e-4 ){
194  error = sqrt(error)/pt2;
195  return 0;
196  }
197  error = 1.e10;
198  return 1;
199 }
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
static void error(int no)
Definition: ranlxd.cxx:419
Double_t
Double_t fC[36]
friend F32vec4 atan2(const F32vec4 &y, const F32vec4 &x)
Definition: P4_F32vec4.h:117
Int_t KFParticleBase::GetPt ( Double_t Pt,
Double_t SigmaPt 
) const

Definition at line 131 of file KFParticleBase.cxx.

References Double_t, fC, fP, pt(), and sqrt().

Referenced by KFParticle::GetErrPt(), and KFParticle::GetPt().

132 {
133  //* Calculate particle transverse momentum
134 
135  Double_t px = fP[3];
136  Double_t py = fP[4];
137  Double_t px2 = px*px;
138  Double_t py2 = py*py;
139  Double_t pt2 = px2+py2;
140  pt = sqrt(pt2);
141  error = (px2*fC[9] + py2*fC[14] + 2*px*py*fC[13] );
142  if( error>0 && pt>1.e-4 ){
143  error = sqrt(error)/pt;
144  return 0;
145  }
146  error = 1.e10;
147  return 1;
148 }
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
static void error(int no)
Definition: ranlxd.cxx:419
TString pt(TString pts, TString exts="px py pz")
Definition: invexp.C:133
Double_t
Double_t fC[36]
Double_t KFParticleBase::GetPx ( ) const
inline

Definition at line 104 of file KFParticleBase.h.

References fP.

Referenced by GetArmenterosPodolanski(), KFParticle::GetExternalTrackParam(), and KFParticle::GetPx().

104 { return fP[3]; }
Double_t fP[8]
Double_t KFParticleBase::GetPy ( ) const
inline

Definition at line 105 of file KFParticleBase.h.

References fP.

Referenced by GetArmenterosPodolanski(), KFParticle::GetExternalTrackParam(), and KFParticle::GetPy().

105 { return fP[4]; }
Double_t fP[8]
Double_t KFParticleBase::GetPz ( ) const
inline

Definition at line 106 of file KFParticleBase.h.

References fP.

Referenced by GetArmenterosPodolanski(), KFParticle::GetExternalTrackParam(), and KFParticle::GetPz().

106 { return fP[5]; }
Double_t fP[8]
Int_t KFParticleBase::GetQ ( ) const
inline
Int_t KFParticleBase::GetR ( Double_t R,
Double_t SigmaR 
) const

Definition at line 201 of file KFParticleBase.cxx.

References Double_t, fC, fP, r, sqrt(), x, and y.

Referenced by KFParticle::GetErrR(), and KFParticle::GetR().

202 {
203  //* Calculate distance to the origin
204 
205  Double_t x = fP[0];
206  Double_t y = fP[1];
207  Double_t x2 = x*x;
208  Double_t y2 = y*y;
209  r = sqrt(x2 + y2);
210  error = (x2*fC[0] + y2*fC[2] - 2*x*y*fC[1] );
211  if( error>0 && r>1.e-4 ){
212  error = sqrt(error)/r;
213  return 0;
214  }
215  error = 1.e10;
216  return 1;
217 }
double r
Definition: RiemannTest.C:14
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
static void error(int no)
Definition: ranlxd.cxx:419
Double_t
Double_t fC[36]
Double_t x
Double_t y
Double_t KFParticleBase::GetS ( ) const
inline

Definition at line 108 of file KFParticleBase.h.

References fP.

Referenced by KFParticle::GetS().

108 { return fP[7]; }
Double_t fP[8]
Double_t KFParticleBase::GetSCorrection ( const Double_t  Part[],
const Double_t  XYZ[] 
)
staticprotected

Definition at line 336 of file KFParticleBase.cxx.

References d, Double_t, p2, and sqrt().

Referenced by GetMeasurement().

337 {
338  //* Get big enough correction for S error to let the particle Part be fitted to XYZ point
339 
340  Double_t d[3] = { XYZ[0]-Part[0], XYZ[1]-Part[1], XYZ[2]-Part[2] };
341  Double_t p2 = Part[3]*Part[3]+Part[4]*Part[4]+Part[5]*Part[5];
342 // Double_t sigmaS = (p2>1.e-4) ? ( 10.1+3.*sqrt( d[0]*d[0]+d[1]*d[1]+d[2]*d[2]) )/sqrt(p2) : 1.;
343  Double_t sigmaS = (p2>1.e-4) ? ( 0.1+10.*sqrt( d[0]*d[0]+d[1]*d[1]+d[2]*d[2]) )/sqrt(p2) : 1.;
344  return sigmaS;
345 }
TObjArray * d
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t
TPad * p2
Definition: hist-t7.C:117
const Double_t& KFParticleBase::GetSumDaughterMass ( ) const
inline

Definition at line 93 of file KFParticleBase.h.

References SumDaughterMass.

93 {return SumDaughterMass;}
Double_t SumDaughterMass
Double_t KFParticleBase::GetX ( ) const
inline

Definition at line 101 of file KFParticleBase.h.

References fP.

Referenced by KFParticle::GetExternalTrackParam(), KFParticle::GetX(), and RotateXY().

101 { return fP[0]; }
Double_t fP[8]
Double_t KFParticleBase::GetY ( ) const
inline

Definition at line 102 of file KFParticleBase.h.

References fP.

Referenced by KFParticle::GetExternalTrackParam(), KFParticle::GetY(), and RotateXY().

102 { return fP[1]; }
Double_t fP[8]
Double_t KFParticleBase::GetZ ( ) const
inline

Definition at line 103 of file KFParticleBase.h.

References fP.

Referenced by KFParticle::GetExternalTrackParam(), KFParticle::GetZ(), and RotateXY().

103 { return fP[2]; }
Double_t fP[8]
int KFParticleBase::Id ( ) const
inline

Definition at line 269 of file KFParticleBase.h.

References fId.

Referenced by KFParticleSIMD::KFParticleSIMD().

269 { return fId; };
static Int_t KFParticleBase::IJ ( Int_t  i,
Int_t  j 
)
inlinestaticprotected

Definition at line 281 of file KFParticleBase.h.

References i.

Referenced by Cij(), ConstructGammaBz(), Covariance(), GetCovariance(), and SetMassConstraint().

281  {
282  return ( j<=i ) ? i*(i+1)/2+j :j*(j+1)/2+i;
283  }
Int_t i
Definition: run_full.C:25
void KFParticleBase::Initialize ( const Double_t  Param[],
const Double_t  Cov[],
Int_t  Charge,
Double_t  Mass 
)

Definition at line 32 of file KFParticleBase.cxx.

References Double_t, energy, fAtProductionVertex, fC, fChi2, fIsLinearized, fMassHypo, fNDF, fP, fQ, fSFromDecay, h1, h2, i, Mass, sqrt(), and SumDaughterMass.

33 {
34  // Constructor from "cartesian" track, particle mass hypothesis should be provided
35  //
36  // Param[6] = { X, Y, Z, Px, Py, Pz } - position and momentum
37  // Cov [21] = lower-triangular part of the covariance matrix:
38  //
39  // ( 0 . . . . . )
40  // ( 1 2 . . . . )
41  // Cov. matrix = ( 3 4 5 . . . ) - numbering of covariance elements in Cov[]
42  // ( 6 7 8 9 . . )
43  // ( 10 11 12 13 14 . )
44  // ( 15 16 17 18 19 20 )
45 
46 
47  for( Int_t i=0; i<6 ; i++ ) fP[i] = Param[i];
48  for( Int_t i=0; i<21; i++ ) fC[i] = Cov[i];
49 
50  Double_t energy = sqrt( Mass*Mass + fP[3]*fP[3] + fP[4]*fP[4] + fP[5]*fP[5]);
51  fP[6] = energy;
52  fP[7] = 0;
53  fQ = Charge;
54  fNDF = 0;
55  fChi2 = 0;
57  fIsLinearized = 0;
58  fSFromDecay = 0;
59 
60  Double_t energyInv = 1./energy;
61  Double_t
62  h0 = fP[3]*energyInv,
63  h1 = fP[4]*energyInv,
64  h2 = fP[5]*energyInv;
65 
66  fC[21] = h0*fC[ 6] + h1*fC[10] + h2*fC[15];
67  fC[22] = h0*fC[ 7] + h1*fC[11] + h2*fC[16];
68  fC[23] = h0*fC[ 8] + h1*fC[12] + h2*fC[17];
69  fC[24] = h0*fC[ 9] + h1*fC[13] + h2*fC[18];
70  fC[25] = h0*fC[13] + h1*fC[14] + h2*fC[19];
71  fC[26] = h0*fC[18] + h1*fC[19] + h2*fC[20];
72  fC[27] = ( h0*h0*fC[ 9] + h1*h1*fC[14] + h2*h2*fC[20]
73  + 2*(h0*h1*fC[13] + h0*h2*fC[18] + h1*h2*fC[19] ) );
74  for( Int_t i=28; i<36; i++ ) fC[i] = 0;
75  fC[35] = 1.;
76 
78  fMassHypo = Mass;
79 }
Int_t i
Definition: run_full.C:25
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
Double_t SumDaughterMass
Double_t
Double_t fSFromDecay
Double_t fC[36]
Bool_t fAtProductionVertex
Double_t fMassHypo
Double_t energy
Definition: plot_dirc.C:15
void KFParticleBase::Initialize ( )

Definition at line 81 of file KFParticleBase.cxx.

References fAtProductionVertex, fC, fChi2, fIsLinearized, fMassHypo, fNDF, fP, fQ, fSFromDecay, fVtxGuess, i, and SumDaughterMass.

Referenced by KFParticle::Create(), KFParticle::Initialize(), and KFParticleBase().

82 {
83  //* Initialise covariance matrix and set current parameters to 0.0
84 
85  for( Int_t i=0; i<8; i++) fP[i] = 0;
86  for(Int_t i=0;i<36;++i) fC[i]=0.;
87  fC[0] = fC[2] = fC[5] = 100.;
88  fC[35] = 1.;
89  fNDF = -3;
90  fChi2 = 0.;
91  fQ = 0;
92  fSFromDecay = 0;
94  fVtxGuess[0]=fVtxGuess[1]=fVtxGuess[2]=0.;
95  fIsLinearized = 0;
96  SumDaughterMass = 0;
97  fMassHypo = -1;
98 }
Int_t i
Definition: run_full.C:25
Double_t fVtxGuess[3]
Double_t fP[8]
Double_t SumDaughterMass
Double_t fSFromDecay
Double_t fC[36]
Bool_t fAtProductionVertex
Double_t fMassHypo
Bool_t KFParticleBase::InvertSym3 ( const Double_t  A[],
Double_t  Ainv[] 
)
staticprotected

Definition at line 3300 of file KFParticleBase.cxx.

References Double_t, and fabs().

Referenced by ConstructGammaBz(), and SetProductionVertex().

3301 {
3302  //* Invert symmetric matric stored in low-triagonal form
3303 
3304  bool ret = 0;
3305  double a0 = A[0], a1 = A[1], a2 = A[2], a3 = A[3];
3306 
3307  Ai[0] = a2*A[5] - A[4]*A[4];
3308  Ai[1] = a3*A[4] - a1*A[5];
3309  Ai[3] = a1*A[4] - a2*a3;
3310  Double_t det = (a0*Ai[0] + a1*Ai[1] + a3*Ai[3]);
3311  if( fabs(det)>1.e-20 ) det = 1./det;
3312  else{
3313  det = 0;
3314  ret = 1;
3315  }
3316  Ai[0] *= det;
3317  Ai[1] *= det;
3318  Ai[3] *= det;
3319  Ai[2] = ( a0*A[5] - a3*a3 )*det;
3320  Ai[4] = ( a1*a3 - a0*A[4] )*det;
3321  Ai[5] = ( a0*a2 - a1*a1 )*det;
3322  return ret;
3323 }
Double_t
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
void KFParticleBase::MultQSQt ( const Double_t  Q[],
const Double_t  S[],
Double_t  SOut[] 
)
staticprotected

Definition at line 3325 of file KFParticleBase.cxx.

References Double_t, and i.

Referenced by TransportCBM().

3326 {
3327  //* Matrix multiplication Q*S*Q^T, Q - square matrix, S - symmetric
3328 
3329  const Int_t kN= 8;
3330  Double_t mA[kN*kN];
3331 
3332  for( Int_t i=0, ij=0; i<kN; i++ ){
3333  for( Int_t j=0; j<kN; j++, ++ij ){
3334  mA[ij] = 0 ;
3335  for( Int_t k=0; k<kN; ++k ) mA[ij]+= S[( k<=i ) ? i*(i+1)/2+k :k*(k+1)/2+i] * Q[ j*kN+k];
3336  }
3337  }
3338 
3339  for( Int_t i=0; i<kN; i++ ){
3340  for( Int_t j=0; j<=i; j++ ){
3341  Int_t ij = ( j<=i ) ? i*(i+1)/2+j :j*(j+1)/2+i;
3342  SOut[ij] = 0 ;
3343  for( Int_t k=0; k<kN; k++ ) SOut[ij] += Q[ i*kN+k ] * mA[ k*kN+j ];
3344  }
3345  }
3346 }
Int_t i
Definition: run_full.C:25
const Double_t & S() const
Double_t
const Int_t & Q() const
int KFParticleBase::NDaughters ( ) const
inline

Definition at line 270 of file KFParticleBase.h.

References fDaughtersIds.

Referenced by KFParticleSIMD::KFParticleSIMD().

270 { return fDaughtersIds.size(); };
std::vector< int > fDaughtersIds
const Int_t& KFParticleBase::NDF ( ) const
inline

Definition at line 123 of file KFParticleBase.h.

References fNDF.

Referenced by KFParticle::NDF().

123 { return fNDF; }
Int_t& KFParticleBase::NDF ( )
inline

Definition at line 156 of file KFParticleBase.h.

References fNDF.

156 { return fNDF; }
void KFParticleBase::operator+= ( const KFParticleBase Daughter)

Definition at line 329 of file KFParticleBase.cxx.

References AddDaughter().

Referenced by KFParticle::operator+=().

330 {
331  //* Add daughter via operator+=
332 
333  AddDaughter( Daughter );
334 }
void AddDaughter(const KFParticleBase &Daughter)
Double_t& KFParticleBase::Parameter ( Int_t  i)
inline

Definition at line 158 of file KFParticleBase.h.

References fP, and i.

Referenced by KFParticle::Parameter().

158 { return fP[i]; }
Int_t i
Definition: run_full.C:25
Double_t fP[8]
const Double_t& KFParticleBase::Px ( ) const
inline

Definition at line 116 of file KFParticleBase.h.

References fP.

Referenced by KFParticle::Px().

116 { return fP[3]; }
Double_t fP[8]
Double_t& KFParticleBase::Px ( )
inline

Definition at line 149 of file KFParticleBase.h.

References fP.

149 { return fP[3]; }
Double_t fP[8]
const Double_t& KFParticleBase::Py ( ) const
inline

Definition at line 117 of file KFParticleBase.h.

References fP.

Referenced by KFParticle::Py().

117 { return fP[4]; }
Double_t fP[8]
Double_t& KFParticleBase::Py ( )
inline

Definition at line 150 of file KFParticleBase.h.

References fP.

150 { return fP[4]; }
Double_t fP[8]
const Double_t& KFParticleBase::Pz ( ) const
inline

Definition at line 118 of file KFParticleBase.h.

References fP.

Referenced by KFParticle::Pz().

118 { return fP[5]; }
Double_t fP[8]
Double_t& KFParticleBase::Pz ( )
inline

Definition at line 151 of file KFParticleBase.h.

References fP.

151 { return fP[5]; }
Double_t fP[8]
const Int_t& KFParticleBase::Q ( ) const
inline

Definition at line 121 of file KFParticleBase.h.

References fQ.

Referenced by KFParticle::Q().

121 { return fQ; }
Int_t& KFParticleBase::Q ( )
inline

Definition at line 154 of file KFParticleBase.h.

References fQ.

154 { return fQ; }
void KFParticleBase::RotateXY ( Double_t  angle,
Double_t  Vtx[3] 
)

Definition at line 3233 of file KFParticleBase.cxx.

References c, cos(), Covariance(), Double_t, fP, GetCovariance(), GetX(), GetY(), GetZ(), i, s, sin(), X(), Y(), and Z().

3234 {
3235  // Rotates the KFParticle object around OZ axis, OZ axis is set by the vertex position
3236  // Double_t angle - angle of rotation in XY plane in [rad]
3237  // Double_t Vtx[3] - position of the vertex in [cm]
3238 
3239  // Before rotation the center of the coordinat system should be moved to the vertex position; move back after rotation
3240  X() = X() - Vtx[0];
3241  Y() = Y() - Vtx[1];
3242  Z() = Z() - Vtx[2];
3243 
3244  // Rotate the kf particle
3245  Double_t c = cos(angle);
3246  Double_t s = sin(angle);
3247 
3248  Double_t mA[8][ 8];
3249  for( Int_t i=0; i<8; i++ ){
3250  for( Int_t j=0; j<8; j++){
3251  mA[i][j] = 0;
3252  }
3253  }
3254  for( int i=0; i<8; i++ ){
3255  mA[i][i] = 1;
3256  }
3257  mA[0][0] = c; mA[0][1] = s;
3258  mA[1][0] = -s; mA[1][1] = c;
3259  mA[3][3] = c; mA[3][4] = s;
3260  mA[4][3] = -s; mA[4][4] = c;
3261 
3262  Double_t mAC[8][8];
3263  Double_t mAp[8];
3264 
3265  for( Int_t i=0; i<8; i++ ){
3266  mAp[i] = 0;
3267  for( Int_t k=0; k<8; k++){
3268  mAp[i]+=mA[i][k] * fP[k];
3269  }
3270  }
3271 
3272  for( Int_t i=0; i<8; i++){
3273  fP[i] = mAp[i];
3274  }
3275 
3276  for( Int_t i=0; i<8; i++ ){
3277  for( Int_t j=0; j<8; j++ ){
3278  mAC[i][j] = 0;
3279  for( Int_t k=0; k<8; k++ ){
3280  mAC[i][j]+= mA[i][k] * GetCovariance(k,j);
3281  }
3282  }
3283  }
3284 
3285  for( Int_t i=0; i<8; i++ ){
3286  for( Int_t j=0; j<=i; j++ ){
3287  Double_t xx = 0;
3288  for( Int_t k=0; k<8; k++){
3289  xx+= mAC[i][k]*mA[j][k];
3290  }
3291  Covariance(i,j) = xx;
3292  }
3293  }
3294 
3295  X() = GetX() + Vtx[0];
3296  Y() = GetY() + Vtx[1];
3297  Z() = GetZ() + Vtx[2];
3298 }
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
Int_t i
Definition: run_full.C:25
Double_t fP[8]
Double_t GetX() const
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
TLorentzVector s
Definition: Pnd2DStar.C:50
Double_t GetY() const
Double_t GetCovariance(Int_t i) const
Double_t GetZ() const
Double_t
const Double_t & X() const
const Double_t & Y() const
Double_t angle
Double_t & Covariance(Int_t i)
const Double_t & Z() const
const Double_t& KFParticleBase::S ( ) const
inline

Definition at line 120 of file KFParticleBase.h.

References fP.

Referenced by KFParticle::S().

120 { return fP[7]; }
Double_t fP[8]
Double_t& KFParticleBase::S ( )
inline

Definition at line 153 of file KFParticleBase.h.

References fP.

153 { return fP[7]; }
Double_t fP[8]
void KFParticleBase::SetConstructMethod ( Int_t  m)
inline

Definition at line 86 of file KFParticleBase.h.

References fConstructMethod, and m.

__m128 m
Definition: P4_F32vec4.h:28
void KFParticleBase::SetId ( int  id)
inline

Definition at line 273 of file KFParticleBase.h.

References fId.

Referenced by KFParticleFinder::Find2DaughterDecay(), KFParticleFinder::FindParticles(), KFParticleFinder::FindTrackV0Decay(), and KFParticleSIMD::GetKFParticle().

273 { fId = id; }; // should be always used (manualy)
void KFParticleBase::SetMassConstraint ( Double_t  Mass,
Double_t  SigmaMass = 0 
)

Definition at line 1403 of file KFParticleBase.cxx.

References Cij(), Double_t, fC, fChi2, fMassHypo, fNDF, fP, i, m2(), Mass, p2, sqrt(), and SumDaughterMass.

Referenced by AddDaughterWithEnergyFitMC(), Construct(), KFParticle::SetMassConstraint(), and SetNonlinearMassConstraint().

1404 {
1405  //* Set hard( SigmaMass=0 ) or soft (SigmaMass>0) mass constraint
1406 
1407  fMassHypo = Mass;
1409 
1410  Double_t m2 = Mass*Mass; // measurement, weighted by Mass
1411  Double_t s2 = m2*SigmaMass*SigmaMass; // sigma^2
1412 
1413  Double_t p2 = fP[3]*fP[3] + fP[4]*fP[4] + fP[5]*fP[5];
1414  Double_t e0 = sqrt(m2+p2);
1415 
1416  Double_t mH[8];
1417  mH[0] = mH[1] = mH[2] = 0.;
1418  mH[3] = -2*fP[3];
1419  mH[4] = -2*fP[4];
1420  mH[5] = -2*fP[5];
1421  mH[6] = 2*fP[6];//e0;
1422  mH[7] = 0;
1423 
1424  Double_t zeta = e0*e0 - e0*fP[6];
1425  zeta = m2 - (fP[6]*fP[6]-p2);
1426 
1427  Double_t mCHt[8], s2_est=0;
1428  for( Int_t i=0; i<8; ++i ){
1429  mCHt[i] = 0.0;
1430  for (Int_t j=0;j<8;++j) mCHt[i] += Cij(i,j)*mH[j];
1431  s2_est += mH[i]*mCHt[i];
1432  }
1433 
1434  if( s2_est<1.e-20 ) return; // calculated mass error is already 0,
1435  // the particle can not be constrained on mass
1436 
1437  Double_t w2 = 1./( s2 + s2_est );
1438  fChi2 += zeta*zeta*w2;
1439  fNDF += 1;
1440  for( Int_t i=0, ii=0; i<8; ++i ){
1441  Double_t ki = mCHt[i]*w2;
1442  fP[i]+= ki*zeta;
1443  for(Int_t j=0;j<=i;++j) fC[ii++] -= ki*mCHt[j];
1444  }
1445 }
Int_t i
Definition: run_full.C:25
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
Double_t & Cij(Int_t i, Int_t j)
TString m2(TString pts, TString exts="e px py pz")
Definition: invexp.C:117
Double_t SumDaughterMass
Double_t
Double_t fC[36]
TPad * p2
Definition: hist-t7.C:117
Double_t fMassHypo
void KFParticleBase::SetMassConstraint ( Double_t mP,
Double_t mC,
Double_t  mJ[7][7],
Double_t  mass 
)
protected

Definition at line 1296 of file KFParticleBase.cxx.

References a, b, c, d, Double_t, f, fabs(), i, IJ(), lambda(), p2, and sqrt().

1297 {
1298  //* Set nonlinear mass constraint (Mass) on the state vector mP with a covariance matrix mC.
1299 
1300  const Double_t energy2 = mP[6]*mP[6], p2 = mP[3]*mP[3]+mP[4]*mP[4]+mP[5]*mP[5], mass2 = mass*mass;
1301 
1302  const Double_t a = energy2 - p2 + 2.*mass2;
1303  const Double_t b = -2.*(energy2 + p2);
1304  const Double_t c = energy2 - p2 - mass2;
1305 
1306  Double_t lambda = 0;
1307  if(fabs(b) > 1.e-10) lambda = -c / b;
1308 
1309  Double_t d = 4.*energy2*p2 - mass2*(energy2-p2-2.*mass2);
1310  if(d>=0 && fabs(a) > 1.e-10) lambda = (energy2 + p2 - sqrt(d))/a;
1311 
1312  if(mP[6] < 0) //If energy < 0 we need a lambda < 0
1313  lambda = -1000000.; //Empirical, a better solution should be found
1314 
1315  Int_t iIter=0;
1316  for(iIter=0; iIter<100; iIter++)
1317  {
1318  Double_t lambda2 = lambda*lambda;
1319  Double_t lambda4 = lambda2*lambda2;
1320 
1321  Double_t lambda0 = lambda;
1322 
1323  Double_t f = -mass2 * lambda4 + a*lambda2 + b*lambda + c;
1324  Double_t df = -4.*mass2 * lambda2*lambda + 2.*a*lambda + b;
1325  if(fabs(df) < 1.e-10) break;
1326  lambda -= f/df;
1327  if(fabs(lambda0 - lambda) < 1.e-8) break;
1328  }
1329 
1330  const Double_t lpi = 1./(1. + lambda);
1331  const Double_t lmi = 1./(1. - lambda);
1332  const Double_t lp2i = lpi*lpi;
1333  const Double_t lm2i = lmi*lmi;
1334 
1335  Double_t lambda2 = lambda*lambda;
1336 
1337  Double_t dfl = -4.*mass2 * lambda2*lambda + 2.*a*lambda + b;
1338  Double_t dfx[7] = {0};//,0,0,0};
1339  dfx[0] = -2.*(1. + lambda)*(1. + lambda)*mP[3];
1340  dfx[1] = -2.*(1. + lambda)*(1. + lambda)*mP[4];
1341  dfx[2] = -2.*(1. + lambda)*(1. + lambda)*mP[5];
1342  dfx[3] = 2.*(1. - lambda)*(1. - lambda)*mP[6];
1343  Double_t dlx[4] = {1,1,1,1};
1344  if(fabs(dfl) > 1.e-10 )
1345  {
1346  for(int i=0; i<4; i++)
1347  dlx[i] = -dfx[i] / dfl;
1348  }
1349 
1350  Double_t dxx[4] = {mP[3]*lm2i, mP[4]*lm2i, mP[5]*lm2i, -mP[6]*lp2i};
1351 
1352  for(Int_t i=0; i<7; i++)
1353  for(Int_t j=0; j<7; j++)
1354  mJ[i][j]=0;
1355  mJ[0][0] = 1.;
1356  mJ[1][1] = 1.;
1357  mJ[2][2] = 1.;
1358 
1359  for(Int_t i=3; i<7; i++)
1360  for(Int_t j=3; j<7; j++)
1361  mJ[i][j] = dlx[j-3]*dxx[i-3];
1362 
1363  for(Int_t i=3; i<6; i++)
1364  mJ[i][i] += lmi;
1365  mJ[6][6] += lpi;
1366 
1367  Double_t mCJ[7][7];
1368 
1369  for(Int_t i=0; i<7; i++) {
1370  for(Int_t j=0; j<7; j++) {
1371  mCJ[i][j] = 0;
1372  for(Int_t k=0; k<7; k++) {
1373  mCJ[i][j] += mC[IJ(i,k)]*mJ[j][k];
1374  }
1375  }
1376  }
1377 
1378  for(Int_t i=0; i<7; ++i){
1379  for(Int_t j=0; j<=i; ++j){
1380  mC[IJ(i,j)]=0;
1381  for(Int_t l=0; l<7; l++){
1382  mC[IJ(i,j)] += mJ[i][l]*mCJ[l][j];
1383  }
1384  }
1385  }
1386 
1387  mP[3] *= lmi;
1388  mP[4] *= lmi;
1389  mP[5] *= lmi;
1390  mP[6] *= lpi;
1391 }
TObjArray * d
double mP
Int_t i
Definition: run_full.C:25
TTree * b
Double_t lambda(Double_t x, Double_t y, Double_t z)
Definition: drawdal.C:48
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
static Int_t IJ(Int_t i, Int_t j)
Int_t a
Definition: anaLmdDigi.C:126
Double_t
TFile * f
Definition: bump_analys.C:12
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
TPad * p2
Definition: hist-t7.C:117
void KFParticleBase::SetMassHypo ( Double_t  m)
inline

Definition at line 89 of file KFParticleBase.h.

References fMassHypo, and m.

89 { fMassHypo = m;}
__m128 m
Definition: P4_F32vec4.h:28
Double_t fMassHypo
void KFParticleBase::SetNoDecayLength ( )

Definition at line 1448 of file KFParticleBase.cxx.

References Double_t, fC, fChi2, fNDF, fP, h, i, s, and TransportToDecayVertex().

Referenced by KFParticle::SetNoDecayLength().

1449 {
1450  //* Set no decay length for resonances
1451 
1453 
1454  Double_t h[8];
1455  h[0] = h[1] = h[2] = h[3] = h[4] = h[5] = h[6] = 0;
1456  h[7] = 1;
1457 
1458  Double_t zeta = 0 - fP[7];
1459  for(Int_t i=0;i<8;++i) zeta -= h[i]*(fP[i]-fP[i]);
1460 
1461  Double_t s = fC[35];
1462  if( s>1.e-20 ){
1463  s = 1./s;
1464  fChi2 += zeta*zeta*s;
1465  fNDF += 1;
1466  for( Int_t i=0, ii=0; i<7; ++i ){
1467  Double_t ki = fC[28+i]*s;
1468  fP[i]+= ki*zeta;
1469  for(Int_t j=0;j<=i;++j) fC[ii++] -= ki*fC[28+j];
1470  }
1471  }
1472  fP[7] = 0;
1473  fC[28] = fC[29] = fC[30] = fC[31] = fC[32] = fC[33] = fC[34] = fC[35] = 0;
1474 }
void TransportToDecayVertex()
Int_t i
Definition: run_full.C:25
Double_t fP[8]
TLorentzVector s
Definition: Pnd2DStar.C:50
Double_t
Double_t fC[36]
void KFParticleBase::SetNonlinearMassConstraint ( Double_t  Mass)

Definition at line 1393 of file KFParticleBase.cxx.

References Double_t, fC, fMassHypo, fP, SetMassConstraint(), and SumDaughterMass.

Referenced by ConstructGammaBz(), KFParticleFinder::Find2DaughterDecay(), and KFParticleFinder::FindTrackV0Decay().

1394 {
1395  //* Set nonlinear mass constraint (mass)
1396 
1397  Double_t mJ[7][7];
1398  SetMassConstraint( fP, fC, mJ, mass );
1399  fMassHypo = mass;
1400  SumDaughterMass = mass;
1401 }
Double_t fP[8]
Double_t SumDaughterMass
Double_t
void SetMassConstraint(Double_t Mass, Double_t SigmaMass=0)
Double_t fC[36]
Double_t fMassHypo
void KFParticleBase::SetPDG ( int  pdg)
inline

Definition at line 276 of file KFParticleBase.h.

References fPDG.

Referenced by KFParticleFinder::FindParticles(), and KFParticleSIMD::GetKFParticle().

276 { fPDG = pdg; }
void KFParticleBase::SetProductionVertex ( const KFParticleBase Vtx)

Definition at line 1145 of file KFParticleBase.cxx.

References Bool_t, Convert(), d0, Double_t, fabs(), fC, fChi2, fNDF, fP, fSFromDecay, GetDStoPoint(), InvertSym3(), m, TransportToDecayVertex(), TransportToDS(), and z.

Referenced by Construct(), and KFParticle::SetProductionVertex().

1146 {
1147  //* Set production vertex for the particle, when the particle was not used in the vertex fit
1148 
1149  const Double_t *m = Vtx.fP, *mV = Vtx.fC;
1150 
1151  Bool_t noS = ( fC[35]<=0 ); // no decay length allowed
1152 
1153  if( noS ){
1155  fP[7] = 0;
1156  fC[28] = fC[29] = fC[30] = fC[31] = fC[32] = fC[33] = fC[34] = fC[35] = 0;
1157  } else {
1158  TransportToDS( GetDStoPoint( m ) );
1159  fP[7] = -fSFromDecay;
1160  fC[28] = fC[29] = fC[30] = fC[31] = fC[32] = fC[33] = fC[34] = 0;
1161  fC[35] = 0.1;
1162 
1163  Convert(1);
1164  }
1165 
1166  Double_t mAi[6];
1167 
1168  InvertSym3( fC, mAi );
1169 
1170  Double_t mB[5][3];
1171 
1172  mB[0][0] = fC[ 6]*mAi[0] + fC[ 7]*mAi[1] + fC[ 8]*mAi[3];
1173  mB[0][1] = fC[ 6]*mAi[1] + fC[ 7]*mAi[2] + fC[ 8]*mAi[4];
1174  mB[0][2] = fC[ 6]*mAi[3] + fC[ 7]*mAi[4] + fC[ 8]*mAi[5];
1175 
1176  mB[1][0] = fC[10]*mAi[0] + fC[11]*mAi[1] + fC[12]*mAi[3];
1177  mB[1][1] = fC[10]*mAi[1] + fC[11]*mAi[2] + fC[12]*mAi[4];
1178  mB[1][2] = fC[10]*mAi[3] + fC[11]*mAi[4] + fC[12]*mAi[5];
1179 
1180  mB[2][0] = fC[15]*mAi[0] + fC[16]*mAi[1] + fC[17]*mAi[3];
1181  mB[2][1] = fC[15]*mAi[1] + fC[16]*mAi[2] + fC[17]*mAi[4];
1182  mB[2][2] = fC[15]*mAi[3] + fC[16]*mAi[4] + fC[17]*mAi[5];
1183 
1184  mB[3][0] = fC[21]*mAi[0] + fC[22]*mAi[1] + fC[23]*mAi[3];
1185  mB[3][1] = fC[21]*mAi[1] + fC[22]*mAi[2] + fC[23]*mAi[4];
1186  mB[3][2] = fC[21]*mAi[3] + fC[22]*mAi[4] + fC[23]*mAi[5];
1187 
1188  mB[4][0] = fC[28]*mAi[0] + fC[29]*mAi[1] + fC[30]*mAi[3];
1189  mB[4][1] = fC[28]*mAi[1] + fC[29]*mAi[2] + fC[30]*mAi[4];
1190  mB[4][2] = fC[28]*mAi[3] + fC[29]*mAi[4] + fC[30]*mAi[5];
1191 
1192  Double_t z[3] = { m[0]-fP[0], m[1]-fP[1], m[2]-fP[2] };
1193 
1194  {
1195  Double_t mAVi[6] = { fC[0]-mV[0], fC[1]-mV[1], fC[2]-mV[2],
1196  fC[3]-mV[3], fC[4]-mV[4], fC[5]-mV[5] };
1197 
1198  if( !InvertSym3( mAVi, mAVi ) ){
1199 
1200  Double_t dChi2 = ( +(mAVi[0]*z[0] + mAVi[1]*z[1] + mAVi[3]*z[2])*z[0]
1201  +(mAVi[1]*z[0] + mAVi[2]*z[1] + mAVi[4]*z[2])*z[1]
1202  +(mAVi[3]*z[0] + mAVi[4]*z[1] + mAVi[5]*z[2])*z[2] );
1203 
1204  // Take Abs(dChi2) here. Negative value of 'det' or 'dChi2' shows that the particle
1205  // was not used in the production vertex fit
1206 
1207  fChi2+= fabs( dChi2 );
1208  }
1209  fNDF += 2;
1210  }
1211 
1212  fP[0] = m[0];
1213  fP[1] = m[1];
1214  fP[2] = m[2];
1215  fP[3]+= mB[0][0]*z[0] + mB[0][1]*z[1] + mB[0][2]*z[2];
1216  fP[4]+= mB[1][0]*z[0] + mB[1][1]*z[1] + mB[1][2]*z[2];
1217  fP[5]+= mB[2][0]*z[0] + mB[2][1]*z[1] + mB[2][2]*z[2];
1218  fP[6]+= mB[3][0]*z[0] + mB[3][1]*z[1] + mB[3][2]*z[2];
1219  fP[7]+= mB[4][0]*z[0] + mB[4][1]*z[1] + mB[4][2]*z[2];
1220 
1221  Double_t d0, d1, d2;
1222 
1223  fC[0] = mV[0];
1224  fC[1] = mV[1];
1225  fC[2] = mV[2];
1226  fC[3] = mV[3];
1227  fC[4] = mV[4];
1228  fC[5] = mV[5];
1229 
1230  d0= mB[0][0]*mV[0] + mB[0][1]*mV[1] + mB[0][2]*mV[3] - fC[ 6];
1231  d1= mB[0][0]*mV[1] + mB[0][1]*mV[2] + mB[0][2]*mV[4] - fC[ 7];
1232  d2= mB[0][0]*mV[3] + mB[0][1]*mV[4] + mB[0][2]*mV[5] - fC[ 8];
1233 
1234  fC[ 6]+= d0;
1235  fC[ 7]+= d1;
1236  fC[ 8]+= d2;
1237  fC[ 9]+= d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
1238 
1239  d0= mB[1][0]*mV[0] + mB[1][1]*mV[1] + mB[1][2]*mV[3] - fC[10];
1240  d1= mB[1][0]*mV[1] + mB[1][1]*mV[2] + mB[1][2]*mV[4] - fC[11];
1241  d2= mB[1][0]*mV[3] + mB[1][1]*mV[4] + mB[1][2]*mV[5] - fC[12];
1242 
1243  fC[10]+= d0;
1244  fC[11]+= d1;
1245  fC[12]+= d2;
1246  fC[13]+= d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
1247  fC[14]+= d0*mB[1][0] + d1*mB[1][1] + d2*mB[1][2];
1248 
1249  d0= mB[2][0]*mV[0] + mB[2][1]*mV[1] + mB[2][2]*mV[3] - fC[15];
1250  d1= mB[2][0]*mV[1] + mB[2][1]*mV[2] + mB[2][2]*mV[4] - fC[16];
1251  d2= mB[2][0]*mV[3] + mB[2][1]*mV[4] + mB[2][2]*mV[5] - fC[17];
1252 
1253  fC[15]+= d0;
1254  fC[16]+= d1;
1255  fC[17]+= d2;
1256  fC[18]+= d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
1257  fC[19]+= d0*mB[1][0] + d1*mB[1][1] + d2*mB[1][2];
1258  fC[20]+= d0*mB[2][0] + d1*mB[2][1] + d2*mB[2][2];
1259 
1260  d0= mB[3][0]*mV[0] + mB[3][1]*mV[1] + mB[3][2]*mV[3] - fC[21];
1261  d1= mB[3][0]*mV[1] + mB[3][1]*mV[2] + mB[3][2]*mV[4] - fC[22];
1262  d2= mB[3][0]*mV[3] + mB[3][1]*mV[4] + mB[3][2]*mV[5] - fC[23];
1263 
1264  fC[21]+= d0;
1265  fC[22]+= d1;
1266  fC[23]+= d2;
1267  fC[24]+= d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
1268  fC[25]+= d0*mB[1][0] + d1*mB[1][1] + d2*mB[1][2];
1269  fC[26]+= d0*mB[2][0] + d1*mB[2][1] + d2*mB[2][2];
1270  fC[27]+= d0*mB[3][0] + d1*mB[3][1] + d2*mB[3][2];
1271 
1272  d0= mB[4][0]*mV[0] + mB[4][1]*mV[1] + mB[4][2]*mV[3] - fC[28];
1273  d1= mB[4][0]*mV[1] + mB[4][1]*mV[2] + mB[4][2]*mV[4] - fC[29];
1274  d2= mB[4][0]*mV[3] + mB[4][1]*mV[4] + mB[4][2]*mV[5] - fC[30];
1275 
1276  fC[28]+= d0;
1277  fC[29]+= d1;
1278  fC[30]+= d2;
1279  fC[31]+= d0*mB[0][0] + d1*mB[0][1] + d2*mB[0][2];
1280  fC[32]+= d0*mB[1][0] + d1*mB[1][1] + d2*mB[1][2];
1281  fC[33]+= d0*mB[2][0] + d1*mB[2][1] + d2*mB[2][2];
1282  fC[34]+= d0*mB[3][0] + d1*mB[3][1] + d2*mB[3][2];
1283  fC[35]+= d0*mB[4][0] + d1*mB[4][1] + d2*mB[4][2];
1284 
1285  if( noS ){
1286  fP[7] = 0;
1287  fC[28] = fC[29] = fC[30] = fC[31] = fC[32] = fC[33] = fC[34] = fC[35] = 0;
1288  } else {
1289  TransportToDS( fP[7] );
1290  Convert(0);
1291  }
1292 
1293  fSFromDecay = 0;
1294 }
void TransportToDecayVertex()
__m128 m
Definition: P4_F32vec4.h:28
Double_t fP[8]
Double_t d0
Definition: checkhelixhit.C:59
Double_t
static Bool_t InvertSym3(const Double_t A[], Double_t Ainv[])
Double_t fSFromDecay
Double_t z
Double_t fC[36]
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
void Convert(bool ToProduction)
void TransportToDS(Double_t dS)
virtual Double_t GetDStoPoint(const Double_t xyz[]) const =0
void KFParticleBase::SetVtxGuess ( Double_t  x,
Double_t  y,
Double_t  z 
)

Definition at line 100 of file KFParticleBase.cxx.

References fIsLinearized, fVtxGuess, x, y, and z.

Referenced by KFParticle::SetVtxGuess().

101 {
102  //* Set decay vertex parameters for linearisation
103 
104  fVtxGuess[0] = x;
105  fVtxGuess[1] = y;
106  fVtxGuess[2] = z;
107  fIsLinearized = 1;
108 }
Double_t fVtxGuess[3]
Double_t z
Double_t x
Double_t y
void KFParticleBase::SubtractFromVertex ( KFParticleBase Vtx) const

Definition at line 2654 of file KFParticleBase.cxx.

References Double_t, fC, fChi2, fIsLinearized, fNDF, fP, fVtxGuess, GetFieldValue(), GetMeasurement(), i, m, and s.

Referenced by KFParticle::SubtractFromVertex().

2655 {
2656  //* Subtract the particle from the vertex
2657 
2658  Double_t fld[3];
2659  {
2660  GetFieldValue( Vtx.fP, fld );
2661  const Double_t kCLight = 0.000299792458;
2662  fld[0]*=kCLight; fld[1]*=kCLight; fld[2]*=kCLight;
2663  }
2664 
2665  Double_t m[8];
2666  Double_t mCm[36];
2667 
2668  if( Vtx.fIsLinearized ){
2669  GetMeasurement( Vtx.fVtxGuess, m, mCm );
2670  } else {
2671  GetMeasurement( Vtx.fP, m, mCm );
2672  }
2673 
2674  Double_t mV[6];
2675 
2676  mV[ 0] = mCm[ 0];
2677  mV[ 1] = mCm[ 1];
2678  mV[ 2] = mCm[ 2];
2679  mV[ 3] = mCm[ 3];
2680  mV[ 4] = mCm[ 4];
2681  mV[ 5] = mCm[ 5];
2682 
2683  //*
2684 
2685  Double_t mS[6];
2686  {
2687  Double_t mSi[6] = { mV[0]-Vtx.fC[0],
2688  mV[1]-Vtx.fC[1], mV[2]-Vtx.fC[2],
2689  mV[3]-Vtx.fC[3], mV[4]-Vtx.fC[4], mV[5]-Vtx.fC[5] };
2690 
2691  mS[0] = mSi[2]*mSi[5] - mSi[4]*mSi[4];
2692  mS[1] = mSi[3]*mSi[4] - mSi[1]*mSi[5];
2693  mS[2] = mSi[0]*mSi[5] - mSi[3]*mSi[3];
2694  mS[3] = mSi[1]*mSi[4] - mSi[2]*mSi[3];
2695  mS[4] = mSi[1]*mSi[3] - mSi[0]*mSi[4];
2696  mS[5] = mSi[0]*mSi[2] - mSi[1]*mSi[1];
2697 
2698  Double_t s = ( mSi[0]*mS[0] + mSi[1]*mS[1] + mSi[3]*mS[3] );
2699  s = ( s > 1.E-20 ) ?1./s :0;
2700  mS[0]*=s;
2701  mS[1]*=s;
2702  mS[2]*=s;
2703  mS[3]*=s;
2704  mS[4]*=s;
2705  mS[5]*=s;
2706  }
2707 
2708  //* Residual (measured - estimated)
2709 
2710  Double_t zeta[3] = { m[0]-Vtx.fP[0], m[1]-Vtx.fP[1], m[2]-Vtx.fP[2] };
2711 
2712  //* mCHt = mCH' - D'
2713 
2714  Double_t mCHt0[3], mCHt1[3], mCHt2[3];
2715 
2716  mCHt0[0]=Vtx.fC[ 0] ; mCHt1[0]=Vtx.fC[ 1] ; mCHt2[0]=Vtx.fC[ 3] ;
2717  mCHt0[1]=Vtx.fC[ 1] ; mCHt1[1]=Vtx.fC[ 2] ; mCHt2[1]=Vtx.fC[ 4] ;
2718  mCHt0[2]=Vtx.fC[ 3] ; mCHt1[2]=Vtx.fC[ 4] ; mCHt2[2]=Vtx.fC[ 5] ;
2719 
2720  //* Kalman gain K = mCH'*S
2721 
2722  Double_t k0[3], k1[3], k2[3];
2723 
2724  for(Int_t i=0;i<3;++i){
2725  k0[i] = mCHt0[i]*mS[0] + mCHt1[i]*mS[1] + mCHt2[i]*mS[3];
2726  k1[i] = mCHt0[i]*mS[1] + mCHt1[i]*mS[2] + mCHt2[i]*mS[4];
2727  k2[i] = mCHt0[i]*mS[3] + mCHt1[i]*mS[4] + mCHt2[i]*mS[5];
2728  }
2729 
2730  //* New estimation of the vertex position r += K*zeta
2731 
2732  Double_t dChi2 = -(mS[0]*zeta[0] + mS[1]*zeta[1] + mS[3]*zeta[2])*zeta[0]
2733  + (mS[1]*zeta[0] + mS[2]*zeta[1] + mS[4]*zeta[2])*zeta[1]
2734  + (mS[3]*zeta[0] + mS[4]*zeta[1] + mS[5]*zeta[2])*zeta[2];
2735 
2736  if( Vtx.fChi2 - dChi2 < 0 ) return;
2737 
2738  for(Int_t i=0;i<3;++i)
2739  Vtx.fP[i] -= k0[i]*zeta[0] + k1[i]*zeta[1] + k2[i]*zeta[2];
2740 
2741  //* New covariance matrix C -= K*(mCH')'
2742 
2743  for(Int_t i=0, k=0;i<3;++i){
2744  for(Int_t j=0;j<=i;++j,++k)
2745  Vtx.fC[k] += k0[i]*mCHt0[j] + k1[i]*mCHt1[j] + k2[i]*mCHt2[j];
2746  }
2747 
2748  //* Calculate Chi^2
2749 
2750  Vtx.fNDF -= 2;
2751  Vtx.fChi2 -= dChi2;
2752 }
Int_t i
Definition: run_full.C:25
__m128 m
Definition: P4_F32vec4.h:28
Double_t fVtxGuess[3]
Double_t fP[8]
TLorentzVector s
Definition: Pnd2DStar.C:50
void GetMeasurement(const Double_t XYZ[], Double_t m[], Double_t V[]) const
virtual void GetFieldValue(const Double_t xyz[], Double_t B[]) const =0
Double_t
Double_t fC[36]
virtual void KFParticleBase::Transport ( Double_t  dS,
Double_t  P[],
Double_t  C[] 
) const
pure virtual
void KFParticleBase::TransportBz ( Double_t  Bz,
Double_t  dS,
Double_t  P[],
Double_t  C[] 
) const

Definition at line 2413 of file KFParticleBase.cxx.

References b, c, c6, c7, c8, cos(), Double_t, fabs(), fC, fP, fQ, pz, s, sin(), sqrt(), and t.

Referenced by KFParticle::Transport().

2415 {
2416  //* Transport the particle on dS, output to P[],C[], for Bz field
2417 
2418  const Double_t kCLight = 0.000299792458;
2419  b = b*fQ*kCLight;
2420  Double_t bs= b*t;
2421  Double_t s = sin(bs), c = cos(bs);
2422  Double_t sB, cB;
2423  if( fabs(bs)>1.e-10){
2424  sB= s/b;
2425  cB= (1-c)/b;
2426  }else{
2427  const Double_t kOvSqr6 = 1./sqrt(6.);
2428  sB = (1.-bs*kOvSqr6)*(1.+bs*kOvSqr6)*t;
2429  cB = .5*sB*bs;
2430  }
2431 
2432  Double_t px = fP[3];
2433  Double_t py = fP[4];
2434  Double_t pz = fP[5];
2435 
2436  p[0] = fP[0] + sB*px + cB*py;
2437  p[1] = fP[1] - cB*px + sB*py;
2438  p[2] = fP[2] + t*pz;
2439  p[3] = c*px + s*py;
2440  p[4] = -s*px + c*py;
2441  p[5] = fP[5];
2442  p[6] = fP[6];
2443  p[7] = fP[7];
2444 
2445  /*
2446  Double_t mJ[8][8] = { {1,0,0, sB, cB, 0, 0, 0 },
2447  {0,1,0, -cB, sB, 0, 0, 0 },
2448  {0,0,1, 0, 0, t, 0, 0 },
2449  {0,0,0, c, s, 0, 0, 0 },
2450  {0,0,0, -s, c, 0, 0, 0 },
2451  {0,0,0, 0, 0, 1, 0, 0 },
2452  {0,0,0, 0, 0, 0, 1, 0 },
2453  {0,0,0, 0, 0, 0, 0, 1 } };
2454  Double_t mA[8][8];
2455  for( Int_t k=0,i=0; i<8; i++)
2456  for( Int_t j=0; j<=i; j++, k++ ) mA[i][j] = mA[j][i] = fC[k];
2457 
2458  Double_t mJC[8][8];
2459  for( Int_t i=0; i<8; i++ )
2460  for( Int_t j=0; j<8; j++ ){
2461  mJC[i][j]=0;
2462  for( Int_t k=0; k<8; k++ ) mJC[i][j]+=mJ[i][k]*mA[k][j];
2463  }
2464 
2465  for( Int_t k=0,i=0; i<8; i++)
2466  for( Int_t j=0; j<=i; j++, k++ ){
2467  e[k] = 0;
2468  for( Int_t l=0; l<8; l++ ) e[k]+=mJC[i][l]*mJ[j][l];
2469  }
2470 
2471  return;
2472  */
2473 
2474  Double_t
2475  c6=fC[6], c7=fC[7], c8=fC[8], c17=fC[17], c18=fC[18],
2476  c24 = fC[24], c31 = fC[31];
2477 
2478  Double_t
2479  cBC13 = cB*fC[13],
2480  mJC13 = c7 - cB*fC[9] + sB*fC[13],
2481  mJC14 = fC[11] - cBC13 + sB*fC[14],
2482  mJC23 = c8 + t*c18,
2483  mJC24 = fC[12] + t*fC[19],
2484  mJC33 = c*fC[9] + s*fC[13],
2485  mJC34 = c*fC[13] + s*fC[14],
2486  mJC43 = -s*fC[9] + c*fC[13],
2487  mJC44 = -s*fC[13] + c*fC[14];
2488 
2489 
2490  e[0]= fC[0] + 2*(sB*c6 + cB*fC[10]) + (sB*fC[9] + 2*cBC13)*sB + cB*cB*fC[14];
2491  e[1]= fC[1] - cB*c6 + sB*fC[10] + mJC13*sB + mJC14*cB;
2492  e[2]= fC[2] - cB*c7 + sB*fC[11] - mJC13*cB + mJC14*sB;
2493  e[3]= fC[3] + t*fC[15] + mJC23*sB + mJC24*cB;
2494  e[4]= fC[4] + t*fC[16] - mJC23*cB + mJC24*sB;
2495 
2496  e[15]= fC[15] + c18*sB + fC[19]*cB;
2497  e[16]= fC[16] - c18*cB + fC[19]*sB;
2498  e[17]= c17 + fC[20]*t;
2499  e[18]= c18*c + fC[19]*s;
2500  e[19]= -c18*s + fC[19]*c;
2501 
2502  e[5]= fC[5] + (c17 + e[17] )*t;
2503 
2504  e[6]= c*c6 + s*fC[10] + mJC33*sB + mJC34*cB;
2505  e[7]= c*c7 + s*fC[11] - mJC33*cB + mJC34*sB;
2506  e[8]= c*c8 + s*fC[12] + e[18]*t;
2507  e[9]= mJC33*c + mJC34*s;
2508  e[10]= -s*c6 + c*fC[10] + mJC43*sB + mJC44*cB;
2509 
2510 
2511  e[11]= -s*c7 + c*fC[11] - mJC43*cB + mJC44*sB;
2512  e[12]= -s*c8 + c*fC[12] + e[19]*t;
2513  e[13]= mJC43*c + mJC44*s;
2514  e[14]= -mJC43*s + mJC44*c;
2515  e[20]= fC[20];
2516  e[21]= fC[21] + fC[25]*cB + c24*sB;
2517  e[22]= fC[22] - c24*cB + fC[25]*sB;
2518  e[23]= fC[23] + fC[26]*t;
2519  e[24]= c*c24 + s*fC[25];
2520  e[25]= c*fC[25] - c24*s;
2521  e[26]= fC[26];
2522  e[27]= fC[27];
2523  e[28]= fC[28] + fC[32]*cB + c31*sB;
2524  e[29]= fC[29] - c31*cB + fC[32]*sB;
2525  e[30]= fC[30] + fC[33]*t;
2526  e[31]= c*c31 + s*fC[32];
2527  e[32]= c*fC[32] - s*c31;
2528  e[33]= fC[33];
2529  e[34]= fC[34];
2530  e[35]= fC[35];
2531 }
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
TTree * b
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Double_t fP[8]
friend F32vec4 sin(const F32vec4 &a)
Definition: P4_F32vec4.h:111
TLorentzVector s
Definition: Pnd2DStar.C:50
Double_t p
Definition: anasim.C:58
TCanvas * c7
Double_t
Double_t fC[36]
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
TTree * t
Definition: bump_analys.C:13
TCanvas * c8
double pz[39]
Definition: pipisigmas.h:14
void KFParticleBase::TransportCBM ( Double_t  dS,
Double_t  P[],
Double_t  C[] 
) const

Definition at line 2286 of file KFParticleBase.cxx.

References c, c2, Double_t, fC, fP, fQ, GetFieldValue(), i, m, MultQSQt(), n, p1, p2, pz, and TransportLine().

Referenced by KFParticle::Transport().

2288 {
2289  //* Transport the particle on dS, output to P[],C[], for CBM field
2290 
2291  if( fQ==0 ){
2292  TransportLine( dS, P, C );
2293  return;
2294  }
2295 
2296  const Double_t kCLight = 0.000299792458;
2297 
2298  Double_t c = fQ*kCLight;
2299 
2300  // construct coefficients
2301 
2302  Double_t
2303  px = fP[3],
2304  py = fP[4],
2305  pz = fP[5];
2306 
2307  Double_t sx=0, sy=0, sz=0, syy=0, syz=0, syyy=0, ssx=0, ssy=0, ssz=0, ssyy=0, ssyz=0, ssyyy=0;
2308 
2309  { // get field integrals
2310 
2311  Double_t fld[3][3];
2312  Double_t p0[3], p1[3], p2[3];
2313 
2314  // line track approximation
2315 
2316  p0[0] = fP[0];
2317  p0[1] = fP[1];
2318  p0[2] = fP[2];
2319 
2320  p2[0] = fP[0] + px*dS;
2321  p2[1] = fP[1] + py*dS;
2322  p2[2] = fP[2] + pz*dS;
2323 
2324  p1[0] = 0.5*(p0[0]+p2[0]);
2325  p1[1] = 0.5*(p0[1]+p2[1]);
2326  p1[2] = 0.5*(p0[2]+p2[2]);
2327 
2328  // first order track approximation
2329  {
2330  GetFieldValue( p0, fld[0] );
2331  GetFieldValue( p1, fld[1] );
2332  GetFieldValue( p2, fld[2] );
2333 
2334  Double_t ssy1 = ( 7*fld[0][1] + 6*fld[1][1]-fld[2][1] )*c*dS*dS/96.;
2335  Double_t ssy2 = ( fld[0][1] + 2*fld[1][1] )*c*dS*dS/6.;
2336 
2337  p1[0] -= ssy1*pz;
2338  p1[2] += ssy1*px;
2339  p2[0] -= ssy2*pz;
2340  p2[2] += ssy2*px;
2341  }
2342 
2343  GetFieldValue( p0, fld[0] );
2344  GetFieldValue( p1, fld[1] );
2345  GetFieldValue( p2, fld[2] );
2346 
2347  sx = c*( fld[0][0] + 4*fld[1][0] + fld[2][0] )*dS/6.;
2348  sy = c*( fld[0][1] + 4*fld[1][1] + fld[2][1] )*dS/6.;
2349  sz = c*( fld[0][2] + 4*fld[1][2] + fld[2][2] )*dS/6.;
2350 
2351  ssx = c*( fld[0][0] + 2*fld[1][0])*dS*dS/6.;
2352  ssy = c*( fld[0][1] + 2*fld[1][1])*dS*dS/6.;
2353  ssz = c*( fld[0][2] + 2*fld[1][2])*dS*dS/6.;
2354 
2355  Double_t c2[3][3] = { { 5, -4, -1},{ 44, 80, -4},{ 11, 44, 5} }; // /=360.
2356  Double_t cc2[3][3] = { { 38, 8, -4},{ 148, 208, -20},{ 3, 36, 3} }; // /=2520.
2357  for(Int_t n=0; n<3; n++)
2358  for(Int_t m=0; m<3; m++)
2359  {
2360  syz += c2[n][m]*fld[n][1]*fld[m][2];
2361  ssyz += cc2[n][m]*fld[n][1]*fld[m][2];
2362  }
2363 
2364  syz *= c*c*dS*dS/360.;
2365  ssyz *= c*c*dS*dS*dS/2520.;
2366 
2367  syy = c*( fld[0][1] + 4*fld[1][1] + fld[2][1] )*dS;
2368  syyy = syy*syy*syy / 1296;
2369  syy = syy*syy/72;
2370 
2371  ssyy = ( fld[0][1]*( 38*fld[0][1] + 156*fld[1][1] - fld[2][1] )+
2372  fld[1][1]*( 208*fld[1][1] +16*fld[2][1] )+
2373  fld[2][1]*( 3*fld[2][1] )
2374  )*dS*dS*dS*c*c/2520.;
2375  ssyyy =
2376  (
2377  fld[0][1]*( fld[0][1]*( 85*fld[0][1] + 526*fld[1][1] - 7*fld[2][1] )+
2378  fld[1][1]*( 1376*fld[1][1] +84*fld[2][1] )+
2379  fld[2][1]*( 19*fld[2][1] ) )+
2380  fld[1][1]*( fld[1][1]*( 1376*fld[1][1] +256*fld[2][1] )+
2381  fld[2][1]*( 62*fld[2][1] ) )+
2382  fld[2][1]*fld[2][1] *( 3*fld[2][1] )
2383  )*dS*dS*dS*dS*c*c*c/90720.;
2384 
2385  }
2386 
2387  Double_t mJ[8][8];
2388  for( Int_t i=0; i<8; i++ ) for( Int_t j=0; j<8; j++) mJ[i][j]=0;
2389 
2390  mJ[0][0]=1; mJ[0][1]=0; mJ[0][2]=0; mJ[0][3]=dS-ssyy; mJ[0][4]=ssx; mJ[0][5]=ssyyy-ssy;
2391  mJ[1][0]=0; mJ[1][1]=1; mJ[1][2]=0; mJ[1][3]=-ssz; mJ[1][4]=dS; mJ[1][5]=ssx+ssyz;
2392  mJ[2][0]=0; mJ[2][1]=0; mJ[2][2]=1; mJ[2][3]=ssy-ssyyy; mJ[2][4]=-ssx; mJ[2][5]=dS-ssyy;
2393 
2394  mJ[3][0]=0; mJ[3][1]=0; mJ[3][2]=0; mJ[3][3]=1-syy; mJ[3][4]=sx; mJ[3][5]=syyy-sy;
2395  mJ[4][0]=0; mJ[4][1]=0; mJ[4][2]=0; mJ[4][3]=-sz; mJ[4][4]=1; mJ[4][5]=sx+syz;
2396  mJ[5][0]=0; mJ[5][1]=0; mJ[5][2]=0; mJ[5][3]=sy-syyy; mJ[5][4]=-sx; mJ[5][5]=1-syy;
2397  mJ[6][6] = mJ[7][7] = 1;
2398 
2399  P[0] = fP[0] + mJ[0][3]*px + mJ[0][4]*py + mJ[0][5]*pz;
2400  P[1] = fP[1] + mJ[1][3]*px + mJ[1][4]*py + mJ[1][5]*pz;
2401  P[2] = fP[2] + mJ[2][3]*px + mJ[2][4]*py + mJ[2][5]*pz;
2402  P[3] = mJ[3][3]*px + mJ[3][4]*py + mJ[3][5]*pz;
2403  P[4] = mJ[4][3]*px + mJ[4][4]*py + mJ[4][5]*pz;
2404  P[5] = mJ[5][3]*px + mJ[5][4]*py + mJ[5][5]*pz;
2405  P[6] = fP[6];
2406  P[7] = fP[7];
2407 
2408  MultQSQt( mJ[0], fC, C);
2409 
2410 }
Int_t i
Definition: run_full.C:25
__m128 m
Definition: P4_F32vec4.h:28
static void MultQSQt(const Double_t Q[], const Double_t S[], Double_t SOut[])
Double_t fP[8]
int n
c2
Definition: plot_dirc.C:39
int Pic_FED Eff_lEE C()
void TransportLine(Double_t S, Double_t P[], Double_t C[]) const
virtual void GetFieldValue(const Double_t xyz[], Double_t B[]) const =0
Double_t
Double_t fC[36]
TPad * p2
Definition: hist-t7.C:117
GeV c P
TPad * p1
Definition: hist-t7.C:116
double pz[39]
Definition: pipisigmas.h:14
void KFParticleBase::TransportLine ( Double_t  S,
Double_t  P[],
Double_t  C[] 
) const
protected

Definition at line 2754 of file KFParticleBase.cxx.

References c11, c6, Double_t, fC, and fP.

Referenced by TransportCBM().

2756 {
2757  //* Transport the particle as a straight line
2758 
2759  P[0] = fP[0] + dS*fP[3];
2760  P[1] = fP[1] + dS*fP[4];
2761  P[2] = fP[2] + dS*fP[5];
2762  P[3] = fP[3];
2763  P[4] = fP[4];
2764  P[5] = fP[5];
2765  P[6] = fP[6];
2766  P[7] = fP[7];
2767 
2768  Double_t c6 = fC[ 6] + dS*fC[ 9];
2769  Double_t c11 = fC[11] + dS*fC[14];
2770  Double_t c17 = fC[17] + dS*fC[20];
2771  Double_t sc13 = dS*fC[13];
2772  Double_t sc18 = dS*fC[18];
2773  Double_t sc19 = dS*fC[19];
2774 
2775  C[ 0] = fC[ 0] + dS*( fC[ 6] + c6 );
2776  C[ 2] = fC[ 2] + dS*( fC[11] + c11 );
2777  C[ 5] = fC[ 5] + dS*( fC[17] + c17 );
2778 
2779  C[ 7] = fC[ 7] + sc13;
2780  C[ 8] = fC[ 8] + sc18;
2781  C[ 9] = fC[ 9];
2782 
2783  C[12] = fC[12] + sc19;
2784 
2785  C[ 1] = fC[ 1] + dS*( fC[10] + C[ 7] );
2786  C[ 3] = fC[ 3] + dS*( fC[15] + C[ 8] );
2787  C[ 4] = fC[ 4] + dS*( fC[16] + C[12] );
2788  C[ 6] = c6;
2789 
2790  C[10] = fC[10] + sc13;
2791  C[11] = c11;
2792 
2793  C[13] = fC[13];
2794  C[14] = fC[14];
2795  C[15] = fC[15] + sc18;
2796  C[16] = fC[16] + sc19;
2797  C[17] = c17;
2798 
2799  C[18] = fC[18];
2800  C[19] = fC[19];
2801  C[20] = fC[20];
2802  C[21] = fC[21] + dS*fC[24];
2803  C[22] = fC[22] + dS*fC[25];
2804  C[23] = fC[23] + dS*fC[26];
2805 
2806  C[24] = fC[24];
2807  C[25] = fC[25];
2808  C[26] = fC[26];
2809  C[27] = fC[27];
2810  C[28] = fC[28] + dS*fC[31];
2811  C[29] = fC[29] + dS*fC[32];
2812  C[30] = fC[30] + dS*fC[33];
2813 
2814  C[31] = fC[31];
2815  C[32] = fC[32];
2816  C[33] = fC[33];
2817  C[34] = fC[34];
2818  C[35] = fC[35];
2819 }
TCanvas * c11
Double_t fP[8]
int Pic_FED Eff_lEE C()
Double_t
Double_t fC[36]
GeV c P
void KFParticleBase::TransportToDecayVertex ( )

Definition at line 1622 of file KFParticleBase.cxx.

References Convert(), fAtProductionVertex, fSFromDecay, and TransportToDS().

Referenced by AddDaughterWithEnergyCalc(), AddDaughterWithEnergyFit(), AddDaughterWithEnergyFitMC(), SetNoDecayLength(), SetProductionVertex(), and KFParticle::TransportToDecayVertex().

1623 {
1624  //* Transport the particle to its decay vertex
1625 
1626  if( fSFromDecay != 0 ) TransportToDS( -fSFromDecay );
1627  if( fAtProductionVertex ) Convert(0);
1628  fAtProductionVertex = 0;
1629 }
Double_t fSFromDecay
void Convert(bool ToProduction)
void TransportToDS(Double_t dS)
Bool_t fAtProductionVertex
void KFParticleBase::TransportToDS ( Double_t  dS)

Definition at line 1641 of file KFParticleBase.cxx.

References fC, fP, fSFromDecay, and Transport().

Referenced by AddDaughterWithEnergyCalc(), AddDaughterWithEnergyFit(), AddDaughterWithEnergyFitMC(), SetProductionVertex(), TransportToDecayVertex(), KFParticle::TransportToDS(), and TransportToProductionVertex().

1642 {
1643  //* Transport the particle on dS parameter (SignedPath/Momentum)
1644 
1645  Transport( dS, fP, fC );
1646  fSFromDecay+= dS;
1647 }
Double_t fP[8]
virtual void Transport(Double_t dS, Double_t P[], Double_t C[]) const =0
Double_t fSFromDecay
Double_t fC[36]
void KFParticleBase::TransportToProductionVertex ( )

Definition at line 1631 of file KFParticleBase.cxx.

References Convert(), fAtProductionVertex, fP, fSFromDecay, and TransportToDS().

Referenced by KFParticle::TransportToProductionVertex().

1632 {
1633  //* Transport the particle to its production vertex
1634 
1635  if( fSFromDecay != -fP[7] ) TransportToDS( -fSFromDecay-fP[7] );
1636  if( !fAtProductionVertex ) Convert( 1 );
1637  fAtProductionVertex = 1;
1638 }
Double_t fP[8]
Double_t fSFromDecay
void Convert(bool ToProduction)
void TransportToDS(Double_t dS)
Bool_t fAtProductionVertex
const Double_t& KFParticleBase::X ( ) const
inline

Definition at line 113 of file KFParticleBase.h.

References fP.

Referenced by RotateXY(), and KFParticle::X().

113 { return fP[0]; }
Double_t fP[8]
Double_t& KFParticleBase::X ( )
inline

Definition at line 146 of file KFParticleBase.h.

References fP.

146 { return fP[0]; }
Double_t fP[8]
const Double_t& KFParticleBase::Y ( ) const
inline

Definition at line 114 of file KFParticleBase.h.

References fP.

Referenced by RotateXY(), and KFParticle::Y().

114 { return fP[1]; }
Double_t fP[8]
Double_t& KFParticleBase::Y ( )
inline

Definition at line 147 of file KFParticleBase.h.

References fP.

147 { return fP[1]; }
Double_t fP[8]
const Double_t& KFParticleBase::Z ( ) const
inline

Definition at line 115 of file KFParticleBase.h.

References fP.

Referenced by RotateXY(), and KFParticle::Z().

115 { return fP[2]; }
Double_t fP[8]
Double_t& KFParticleBase::Z ( )
inline

Definition at line 148 of file KFParticleBase.h.

References fP.

148 { return fP[2]; }
Double_t fP[8]

Member Data Documentation

Bool_t KFParticleBase::fAtProductionVertex
protected
Double_t KFParticleBase::fC[36]
protected
Double_t KFParticleBase::fChi2
protected
Int_t KFParticleBase::fConstructMethod
protected

Definition at line 322 of file KFParticleBase.h.

Referenced by AddDaughter(), and SetConstructMethod().

std::vector<int> KFParticleBase::fDaughtersIds
protected
int KFParticleBase::fId
protected

Definition at line 330 of file KFParticleBase.h.

Referenced by Id(), and SetId().

Bool_t KFParticleBase::fIsLinearized
protected
Double_t KFParticleBase::fMassHypo
protected
Int_t KFParticleBase::fNDF
protected
Double_t KFParticleBase::fP[8]
protected
int KFParticleBase::fPDG
protected

Definition at line 333 of file KFParticleBase.h.

Referenced by GetPDG(), and SetPDG().

Int_t KFParticleBase::fQ
protected
Double_t KFParticleBase::fSFromDecay
protected
Double_t KFParticleBase::fVtxGuess[3]
protected
Double_t KFParticleBase::SumDaughterMass
protected

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