RhoKinVtxFitter Class Reference

#include <RhoKinVtxFitter.h>

Inheritance diagram for RhoKinVtxFitter:

Public Member Functions
	RhoKinVtxFitter (RhoCandidate *b)
virtual	~RhoKinVtxFitter ()
void	AddMassConstraint (double mass)
double	GetPull ()
void	SetNMaxIterations (int nit=20)
void	SetNIterationsExact (int nit=2)
void	SetMinDChisq (double m=0.001)
Bool_t	Fit ()
Bool_t	FitAll ()
double	Chi2Contribution (const RhoCandidate *)
double	GetChi2 () const
int	GetNdf () const
double	GetProb () const
void	SetVerbose (Bool_t v=kTRUE)

Protected Member Functions
RhoCandidate &	HeadOfTree () const
RhoCandidate *	CopyCand (RhoCandidate *)
	uppermost particle composite in tree More...
RhoCandidate *	CopyTree (RhoCandidate *)
void	InsertChi2 (const RhoCandidate *bc, const double chi2)
void	SetDaugthersFromComposite (RhoCandidate *cand)
void	FindAndAddFinalStateDaughters (RhoCandidate *cand)
void	SetFourMomentumByDaughters (RhoCandidate *composite)
void	SetDecayVertex (RhoCandidate *composite, const TVector3 &vtx, const TMatrixD &CovVV)

Protected Attributes
Bool_t	fVerbose
RhoCandidate *	fHeadOfTree
std::vector< RhoCandidate * >	fDaughters
double	fChiSquare
int	fNDegreesOfFreedom

Private Member Functions
Bool_t	FitNode (RhoCandidate *b)
void	SetMatrices ()
void	ResetMatrices ()
void	ReadMatrix ()
void	ReadKinMatrix ()
void	ReadMassKinMatrix ()
Bool_t	Compute (RhoCandidate *c)
void	SetOutput (RhoCandidate *head)
void	TransportToVertex (TMatrixD &, TMatrixD &, TMatrixD &, TMatrixD &, TMatrixD &)
void	GetCovariance (TMatrixD &a_cov0, TMatrixD &cov_al_x, TMatrixD &V_vtx, TMatrixD &covS)

Private Attributes
TMatrixD	al0
TMatrixD	al1
TMatrixD	V_al0
TMatrixD	V_al1
TMatrixD	covC
TMatrixD	mD
TMatrixD	mE
TMatrixD	md
TMatrixD	mPull
Int_t	fNvar
Int_t	fNpar
Int_t	fNpart
Int_t	fNcon
Int_t	fNc
Int_t	fNiter
Int_t	fNumKnown
Int_t	NumCon
int	niter
double	fChi2Diff
double	fPull
TMatrixD	vtx_ex
TMatrixD	vtx_st
double	fMass
int	fMassConstraint
double	fMinDChisq
int	fNMaxIterations
bool	fIterateExact
int	fnDof
double	fchiSquare

Detailed Description

Definition at line 22 of file RhoKinVtxFitter.h.

Constructor & Destructor Documentation

RhoKinVtxFitter::RhoKinVtxFitter

(

RhoCandidate *

b

)

Definition at line 24 of file RhoKinVtxFitter.cxx.

References fIterateExact, fMassConstraint, fMinDChisq, and fNMaxIterations.

                                                 :
  RhoFitterBase( b )
{
  fMassConstraint =-1;
  //fPointConstraint=-1;

  fMinDChisq=0.01;
  fNMaxIterations=20;
  fIterateExact=false;


}

RhoKinVtxFitter::~RhoKinVtxFitter ( )

virtual

Definition at line 37 of file RhoKinVtxFitter.cxx.

{
}

Member Function Documentation

void RhoKinVtxFitter::AddMassConstraint

(

double

mass

)

Definition at line 41 of file RhoKinVtxFitter.cxx.

References fMass, and fMassConstraint.

Referenced by newana_check_eta().

{
  fMassConstraint = 1;
  fMass=mass;
}

Double_t RhoFitterBase::Chi2Contribution ( const RhoCandidate *  b )

inherited

access to the fitted candidates

Definition at line 86 of file RhoFitterBase.cxx.

References Double_t, RhoFitterBase::fChi2Map, uid(), and RhoCandidate::Uid().

{
  if(!b) return -999.;
  Int_t uid = b->Uid();
  Double_t chi2=fChi2Map[uid];
  return chi2 >=0.0 ? chi2 : -1.;
}

Bool_t RhoKinVtxFitter::Compute ( RhoCandidate *  c )

private

Definition at line 100 of file RhoKinVtxFitter.cxx.

References al0, al1, covC, fabs(), RhoFitterBase::fChiSquare, RhoFitterBase::fDaughters, fIterateExact, fMassConstraint, fMinDChisq, fNc, RhoFitterBase::fNDegreesOfFreedom, fNMaxIterations, fPull, RhoFitterBase::fVerbose, GetCovariance(), RhoVtxPoca::GetPocaVtx(), mD, md, mE, mPull, NumCon, ReadKinMatrix(), ReadMassKinMatrix(), ReadMatrix(), ResetMatrices(), SetMatrices(), TransportToVertex(), V_al0, V_al1, vtx_ex, and vtx_st.

Referenced by FitNode().

{

  // int nd=fDaughters.size();
  int nd=fDaughters.size();
  NumCon=2*nd;

  if(fMassConstraint >0) {
    NumCon += 1;
  }
//   if(fPointConstraint >0) {
//     NumCon += 2;
//   }

  fNDegreesOfFreedom=NumCon-3;


  SetMatrices();
  ResetMatrices();
  ReadMatrix();

  TMatrixD cov_al_x(7*nd,3);

  TVector3 startVtx;
  //Getting point of closed approach as start vertex point
  RhoVtxPoca poca;                              //changed from internal method to class RhoVtxPoca by J.Puetz
  poca.GetPocaVtx(startVtx, c);

  vtx_st[0][0]=startVtx.X();
  vtx_st[1][0]=startVtx.Y();
  vtx_st[2][0]=startVtx.Z();

  vtx_ex=vtx_st;
  if(fVerbose) { cout<<"Initial vertex Position is "<<vtx_ex[0][0]<<" "<<vtx_ex[1][0]<<" "<<vtx_ex[2][0]<<endl; }

  // al1=al0;
  // V_al1=V_al0;
  TransportToVertex(al0,V_al0,al1,V_al1,vtx_ex);


  al0=al1;
  V_al0=V_al1;
  TMatrixD V_vtx(3,3);
  V_vtx[0][0] = 9000.;
  V_vtx[1][1] = 9000.;
  V_vtx[2][2] = 9000.;
  // vtx_ex=vtx_st;

  double ierr =0 ; // used to check inversions
  TMatrixD chi2(1,1);
  TMatrixD chi2_1(1,1);

  chi2[0][0]=2000000.;
  //double tmp_chiSq = 999;

  for(Int_t j1=0; j1<fNMaxIterations; ++j1) {


    fNc=0;
    if(fMassConstraint >0) {
      ReadMassKinMatrix();
    } else {
      ReadKinMatrix();
    }

    TMatrixD mD_t=mD;
    mD_t.T();
    // mD_t=mD_t.Transpose(mD);


    TMatrixD Vd_inv = mD*V_al0*mD_t;
    if(Vd_inv==0) { continue; }
    TMatrixD Vd = Vd_inv.Invert(&ierr);


    //   if( ierr != 0 ){
    //  if(fVerbose) cout << "Inversion of constraint-matrix failed! " << endl;
    //  return 0;}
    //  Vd.Print();


    TMatrixD del_al = al0 - al1;


    // Lagrange multiplier

    //TMatrixD lam0=Vd*md;
    TMatrixD lam0 = Vd* ( mD*del_al + md);
    //    if(fVerbose) cout << " lam0 calculated" << endl;


    //  Position Derivative matrix ...............
    TMatrixD mE_t=mE;
    mE_t.T();
    TMatrixD Vx_inv = mE_t*Vd*mE;
    TMatrixD Vx=Vx_inv.Invert(&ierr);
//     cout<<"   *** *** *** *** ***"<<endl;
//     cout<<"mE     "; mE.Print();
//     cout<<"mE_t   "; mE_t.Print();
//     cout<<"Vd     "; Vd.Print();
//     cout<<"Vx_inv "; Vx_inv.Print();
//     cout<<"Vx     "; Vx.Print();
//     cout<<"   *******************"<<endl;

    // New vertex and covariance ........
    TMatrixD V_vtx_new(3,3);
    TMatrixD vtx_new(vtx_ex);



    vtx_new -= Vx*mE_t*lam0;
//    cout << " New vtx calculated" << endl;
//    vtx_new.Print();
    V_vtx_new = Vx;



    // Final Lagarange multiplier ........
    TMatrixD lam = lam0 + (Vd * mE) * (vtx_new - vtx_ex);
    //  if(fVerbose) cout << " New lam calculated" << endl;



    // New track parameters.............
    TMatrixD al_new(al0);
    al_new -= V_al0*mD_t*lam;
    //   if(fVerbose) cout << " New track param calculated" << endl;



    //chiSquared
    TMatrixD lam_t=lam;
    lam_t.T();
    // TMatrixD chi2_new = lam_t* md;
    //TMatrixD chi2_new = lam_t*(mD*(al0 - al_new) );
    TMatrixD chi2_new = lam_t*(mD*(al0 - al_new)  + md);
    //TMatrixD chi2_new = lam_t*(mD*(al0 - al_new) + mE*(vtx_st-vtx_ex) + md);



    // New Covariance Matrix................
    //   TMatrixD V_al_new(V_al0);
    //  V_al_new-=V_al0*mD_t*Vd*mD*V_al0_t;


    // protect against errors. RK: is that safe to do?
    if(TMath::IsNaN(chi2_new[0][0])) continue;

    double deltaChi=chi2_new[0][0]-chi2[0][0];

    //  Check chi^2. If better yes update the values ..............................
    if (deltaChi>0.1*chi2[0][0]) {continue;}
    if( chi2_new[0][0] < chi2[0][0] ) {
      //if(true){
      vtx_ex = vtx_new;
      al1 = al_new;
      //     V_al0 = V_al_new;
    }


    //     if (j1==0)  {chi2=chi2_new;continue;}

    // If Chi^2 change is small then go out of iteration......................
    if( (j1+1 == fNMaxIterations) ||
        (!fIterateExact && ( fabs(chi2[0][0]-chi2_new[0][0])<fMinDChisq) )  ) {
      chi2 = chi2_new;
      vtx_ex = vtx_new;
      al0 =al_new;
      //
      TMatrixD Vd_new(Vd);
      Vd_new -= Vd*(mE*Vx*mE_t)*Vd.T();
      TMatrixD V_al_new(V_al0);
      V_al_new -= V_al0*(mD_t*Vd_new*mD)*V_al0.T();
      cov_al_x -= V_al0*mD_t*Vd*mE*Vx; // Vertex-track Correlation

      //double covdif=(V_al0[6][6]-V_al_new[6][6]);
      //  if (covdif > 0 ) {mPull[0][0] =(al0[6][0]-al_new[6][0])/sqrt(covdif);}
      mPull[0][0] =(al0[6][0]-al_new[6][0]);
      fPull=mPull[0][0];

      V_al0 = V_al_new;
      V_vtx = V_vtx_new;
      if(fVerbose) { cout <<"iteration Number " << " " << j1 <<" final." <<endl; }
      if(fVerbose) { cout <<" chi2 in iteration" << " " << chi2[0][0] << " pull="<<fPull<< endl; }
      break; // that was the final iteration, stop the loop
    }
    chi2 = chi2_new;
    if(fVerbose) { cout << "iteration Number " << " " << j1 << endl; }
    if(fVerbose) { cout <<" vertex Position is "<<vtx_new[0][0]<<" "<<vtx_new[1][0]<<" "<<vtx_new[2][0]<<endl; }
    if(fVerbose) { cout << " chi2 in iteration" << " " << chi2[0][0] << endl; }
  } // end of iteration-loop

  // tell that the fit failed if we have no updated chi2
  if( TMath::IsNaN(chi2[0][0]) || chi2[0][0]==2000000. ) return kFALSE;

  TMatrixD al_new_vtx(7*nd,1);
  TMatrixD Va_new_vtx(7*nd,7*nd);




  // Trivial way for covariance matrix of composite
// covC=
//  for(Int_t k=0;k<nd;k++) {
//        for(int j = 0; j< nd ; ++j)  {
//   {if(k<=j) CovS(k,j) = V_al_0(k,j);}}
//   CovC -= CovS;
//  covC+=V_al0.GetSub(k*7,(k+1)*7-1,k*7,(k+1)*7-1);}

  GetCovariance(V_al0,cov_al_x,V_vtx,covC);


  //   al_new_vtx=al1;
  //   Va_new_vtx=V_al1;
//  TransportToVertex(al0, V_al0, al_new_vtx, Va_new_vtx, vtx_ex);
// al0=al_new_vtx;
//  V_al0=Va_new_vtx;
  fChiSquare=chi2[0][0];
  // fChi2Diff=chi2_1[0][0]-chi2[0][0];

  return kTRUE;
}

RhoCandidate * RhoFitterBase::CopyCand ( RhoCandidate *  b )

protectedinherited

uppermost particle composite in tree

Definition at line 51 of file RhoFitterBase.cxx.

References RhoFactory::Instance(), RhoFactory::NewCandidate(), RhoCandidate::RemoveAssociations(), and RhoCandidate::SetFit().

Referenced by RhoFitterBase::CopyTree().

{
  RhoCandidate* newCand = RhoFactory::Instance()->NewCandidate ( b );
  newCand->RemoveAssociations();
  b->SetFit(newCand);//ready to be modified
  return newCand;
}

RhoCandidate * RhoFitterBase::CopyTree ( RhoCandidate *  head )

protectedinherited

Definition at line 61 of file RhoFitterBase.cxx.

References RhoFitterBase::CopyCand(), RhoCandidate::Daughter(), i, RhoCandidate::IsComposite(), RhoCandidate::NDaughters(), and RhoCandidate::SetMotherLink().

Referenced by RhoFitterBase::RhoFitterBase().

{
  //std::cout<<"\n\tcopy tree "<<head->Uid()<<" "<<&head<<" "<<head->PdgCode()<<" "<<head->NDaughters()<<"...";
  RhoCandidate* headcopy=CopyCand(head);
  RhoCandidate* daucopy=0;
  RhoCandidate* dau=0;
  for(Int_t i=0;i<head->NDaughters();i++)
  {
    dau=head->Daughter(i);
    //std::cout<<"  daugter "<<dau->Uid()<<" "<<i<<" "<<dau->PdgCode()<<" at "<<dau<<"  ";
    if(dau == head) {
      std::cout<<endl<<"*** Candidate is its own mother???  *** \n"<<std::endl;
      std::cout<<"  print:   "<<*head<<std::endl;;
    }
    assert(dau != head);
    if(dau->IsComposite()) daucopy=CopyTree(dau);
        else daucopy=CopyCand(dau);
        //std::cout<<"CopyTree: copied candidate "<<dau->Uid()<<std::endl;
    daucopy->SetMotherLink(headcopy); //daughter link is set automatically, too
  }
  return headcopy;
}

void RhoFitterBase::FindAndAddFinalStateDaughters ( RhoCandidate *  cand )

protectedinherited

Definition at line 149 of file RhoFitterBase.cxx.

References RhoCandidate::Daughter(), RhoFitterBase::fDaughters, RhoCandidate::IsComposite(), RhoCandidate::IsLocked(), and RhoCandidate::NDaughters().

Referenced by Rho4CFitter::Fit(), RhoKinFitter::Fit(), and Rho4CFitter::FitConserveMasses().

{
  RhoCandidate* tc;
  for (int k=0;k<cand->NDaughters();k++) {
    tc=cand->Daughter(k);
    if (!tc->IsComposite() || tc->IsLocked()) { fDaughters.push_back(tc); }
    else { FindAndAddFinalStateDaughters(tc); }
  }
  return;
}

Bool_t RhoFitterBase::Fit ( )

inherited

This function has to be overridden by the real algoritms implementations

Definition at line 95 of file RhoFitterBase.cxx.

References RhoFitterBase::fChi2Map, RhoFitterBase::fHeadOfTree, and RhoFitterBase::FitNode().

Referenced by ana_check(), ana_check_psi(), ana_complete(), ana_day1(), ana_invariantmass_2pi_tpc(), ana_multi(), anaDMesonsCharged(), anaideal_complete(), anatut_psi2s(), PndTripleAnaTask::DsDs2317Analysis(), PndPmtTask::Exec(), PndProdAnaTask::Exec(), PndLLbarAnaTask::Exec(), PndMyAnalysisTask::Exec(), PndSimpleCombinerTask::Exec(), DecayTreeFitter::Fitter::FitAll(), PndTripleAnaTask::JpsiAnalysis(), newana_check_eta(), poormantracks(), run_ana_eta_c_stt(), run_ana_eta_c_stt_v2(), run_ana_eta_c_tpc(), run_ana_invariantmass_2pi_stt(), run_ana_invariantmass_2pi_tpc(), run_ana_invariantmass_4pi_stt(), run_ana_invariantmass_4pi_tpc(), run_ana_mertens_evt7(), tut_ana(), tut_ana_fast(), tut_ana_fit(), tut_ana_ntp(), tut_ana_ntp_noqa(), and tut_ana_ntp_qa().

{
  fChi2Map.clear();
  return FitNode(fHeadOfTree);
}

Bool_t RhoFitterBase::FitAll ( )

inherited

Definition at line 101 of file RhoFitterBase.cxx.

References RhoFitterBase::fHeadOfTree, RhoCandidate::IsLocked(), and RhoFitterBase::IterateAndFit().

{
  if(fHeadOfTree->IsLocked())
  {
    Warning("RhoFitterBase::FitAll","You tried to fit a locked candidate. Retuning kFALSE now.");
    return kFALSE;
  }
  return IterateAndFit(fHeadOfTree);
}

Bool_t RhoKinVtxFitter::FitNode ( RhoCandidate *  b )

privatevirtual

Reimplemented from RhoFitterBase.

Definition at line 54 of file RhoKinVtxFitter.cxx.

References Bool_t, Compute(), RhoFitterBase::SetDaugthersFromComposite(), and SetOutput().

{
  SetDaugthersFromComposite(cand);
  Bool_t check=Compute(cand);
  SetOutput(cand);
  return check;
}

double RhoFitterBase::GetChi2 ( ) const

inlineinherited

Definition at line 48 of file RhoFitterBase.h.

References RhoFitterBase::fChiSquare.

Referenced by ana_check(), ana_check_psi(), ana_complete(), ana_day1(), ana_multi(), anaideal_complete(), anatut_psi2s(), PndTripleAnaTask::DsDs2317Analysis(), PndPmtTask::Exec(), PndProdAnaTask::Exec(), PndLLbarAnaTask::Exec(), PndAnaWithTrigger::Exec(), PndMyAnalysisTask::Exec(), PndScrutAnaTask::Exec(), PndSimpleCombinerTask::Exec(), PndTripleAnaTask::JpsiAnalysis(), newana_check_eta(), poormantracks(), PndRhoTupleQA::qaFitter(), run_ana_eta_c_stt(), run_ana_eta_c_stt_v2(), run_ana_eta_c_tpc(), PndTripleAnaTask::ThreePiAnalysis(), tut_ana(), tut_ana_fast(), tut_ana_fit(), tut_ana_ntp(), tut_ana_ntp_noqa(), and tut_ana_ntp_qa().

{return fChiSquare;};

void RhoKinVtxFitter::GetCovariance ( TMatrixD &  a_cov0, TMatrixD &  cov_al_x, TMatrixD &  V_vtx, TMatrixD &  covS  )

private

Definition at line 944 of file RhoKinVtxFitter.cxx.

References a, Double_t, RhoFitterBase::fDaughters, RhoCalculationTools::GetBz(), i, and jj.

Referenced by Compute().

{
  int fNDau=fDaughters.size();
  int nd=fNDau;


  TMatrixD cov_al_x_temp=cov_al_x;

  TMatrixD pA(4*nd,7*nd);
  pA.Zero();
  double  sumA=0;
  double a;
  int kN, jN;
  for (int k=0; k<nd; k++) {
    kN=k*7;
    jN=k*4;
    Double_t bField = 0.1*RhoCalculationTools::GetBz(fDaughters[k]->Pos()); // T, assume field in z only
    a = -0.00299792458*bField*fDaughters[k]->GetCharge();
    sumA += a;
    pA[jN][kN]=pA[jN+1][kN+1]=pA[jN+2][kN+2]=pA[jN+3][kN+3]=1;
    pA[jN+1][kN+4]=-a;
    pA[jN][kN+5]=a;
  }

  TMatrixD pB(4,3);
  TMatrixD pB_t(3,4);
  pB.Zero();
  pB[0][1]=-sumA;
  pB[1][0]=sumA;

  TMatrixD covA_al_xi=pA*cov_al_x;
//TMatrixD covB_al_xi=pB*cov_al_x;
  TMatrixD cov_al_xT(3,7*nd);
  cov_al_xT = cov_al_x_temp.T();

  TMatrixD pA_t(7*nd,4*nd);
//pA_t=pA.T();

//TMatrixD covP(4*nd,4*nd);
  TMatrixD covP = pA*a_cov0*(pA_t.Transpose(pA));
  TMatrixD covA=covA_al_xi*(pB_t.Transpose(pB));
  TMatrixD covB=pB*(cov_al_xT)*(pA_t.Transpose(pA));
  TMatrixD covBV=pB*(V_vtx)*(pB_t.Transpose(pB));

  TMatrixD covA_al_x(4,3);
  TMatrixD SumcovP(4,4);
  SumcovP=covBV;

  for (int k=0; k<nd; k++) {
    kN=k*7;
    jN=k*4;
    for (int i=0; i<4; i++) {
      for (int j=0; j<4; j++) {
        SumcovP[i][j] += covP[jN+i][jN+j];
        SumcovP[i][j] += covA[jN+i][j];
        SumcovP[i][j] += covB[i][jN+j];
      }
      for (int jj=0; jj<3; jj++) {
        covA_al_x[i][jj]  += covA_al_xi[jN+i][jj];
      }
    }
  }


  TMatrixD covPX(4,3);
  covPX = covA_al_x ;
  covPX += (pB *V_vtx);
  TMatrixD covPX_t(3,4);
  covPX_t=covPX_t.Transpose(covPX);



  covS.SetSub(0,0,V_vtx);
  covS.SetSub(0,3,covPX_t);
  covS.SetSub(3,0,covPX);
  covS.SetSub(3,3,SumcovP);

}

int RhoFitterBase::GetNdf ( ) const

inlineinherited

Definition at line 49 of file RhoFitterBase.h.

References RhoFitterBase::fNDegreesOfFreedom.

Referenced by poormantracks(), and PndRhoTupleQA::qaFitter().

{return fNDegreesOfFreedom;};

double RhoFitterBase::GetProb ( ) const

inlineinherited

Definition at line 50 of file RhoFitterBase.h.

References RhoFitterBase::fChiSquare, and RhoFitterBase::fNDegreesOfFreedom.

Referenced by ana_complete(), ana_day1(), ana_multi(), anaideal_complete(), PndPmtTask::Exec(), PndLLbarAnaTask::Exec(), poormantracks(), PndRhoTupleQA::qaFitter(), tut_ana(), tut_ana_fast(), tut_ana_fit(), tut_ana_ntp(), tut_ana_ntp_noqa(), and tut_ana_ntp_qa().

{return TMath::Prob(fChiSquare,fNDegreesOfFreedom);};

double RhoKinVtxFitter::GetPull ( )

inline

Definition at line 29 of file RhoKinVtxFitter.h.

References fPull.

{return fPull;}

RhoCandidate& RhoFitterBase::HeadOfTree ( ) const

inlineprotectedinherited

Definition at line 57 of file RhoFitterBase.h.

References RhoFitterBase::fHeadOfTree.

{ return *fHeadOfTree; }

void RhoFitterBase::InsertChi2 ( const RhoCandidate *  bc, const double  chi2  )

inlineprotectedinherited

Definition at line 66 of file RhoFitterBase.h.

References RhoFitterBase::fChi2Map, and RhoCandidate::Uid().

Referenced by RhoKalmanVtxFitter::Calculate().

{fChi2Map[ bc->Uid()] = chi2;}

void RhoKinVtxFitter::ReadKinMatrix ( )

private

Definition at line 511 of file RhoKinVtxFitter.cxx.

References a, al1, Double_t, RhoFitterBase::fDaughters, fNc, fNcon, fNpar, RhoCalculationTools::GetBz(), mD, md, mE, pz, sqrt(), and vtx_ex.

Referenced by Compute().

{

  int  nd=fDaughters.size();
  fNc=0;
  mD.ResizeTo(fNcon,fNpar);
  mE.ResizeTo(fNcon,3);
  md.ResizeTo(fNcon,1);
  for (int k=0; k<nd; k++) {


    int kN=k*7;
    int k2=k*2;
    double delX = vtx_ex[0][0] - al1[kN+4][0];
    double delY = vtx_ex[1][0] - al1[kN+5][0];
    double delZ = vtx_ex[2][0] - al1[kN+6][0];
    double px = al1[kN+0][0];
    double py = al1[kN+1][0];
    double pz = al1[kN+2][0];
    double ch=fDaughters[k]->GetCharge();
    Double_t bField = 0.1*RhoCalculationTools::GetBz(fDaughters[k]->Pos()); // T, assume field in z only
    double a = -0.0029979246*ch*bField;
    double pT_2 = px*px + py*py;

    double J = a*(delX*px + delY*py)/pT_2;
    double Rx = delX - 2.*px*(delX*px + delY*py)/pT_2;
    double Ry = delY - 2.*py*(delX*px + delY*py)/pT_2;

    // if(fabs(J) > 1) {return 0;}
    if (J>=1.0 || J<= -1.0) { J = (J>=1.0 ? 0.99 : -0.99);}

    double S = 1./(pT_2*sqrt(1-(J*J)));
    // if(fVerbose) cout << ch << "ch" << S << "pt2" << J << "bField" << bField <<endl;
    double asin_J = asin(J);
    //charged particle
    if(ch !=0) {
      mD[fNc+0+k2][kN+0]  = delY ;
      //if(fVerbose) cout << al0[kN+4][0] << " " << delY << endl;
      mD[fNc+0+k2][kN+1]  = -(delX) ;
      mD[fNc+0+k2][kN+2]  = 0. ;
      mD[fNc+0+k2][kN+3]  = 0. ;
      mD[fNc+0+k2][kN+4]  =  py + a*delX ;
      mD[fNc+0+k2][kN+5]  = -px + a*delY ;
      mD[fNc+0+k2][kN+6]  = 0. ;

      mD[fNc+1+k2][kN+0] = -pz*S*Rx ;
      mD[fNc+1+k2][kN+1] = -pz*S*Ry ;
      mD[fNc+1+k2][kN+2] = -asin_J/a ;
      mD[fNc+1+k2][kN+3] = 0.;
      mD[fNc+1+k2][kN+4] = px*pz*S;
      mD[fNc+1+k2][kN+5] = py*pz*S;
      mD[fNc+1+k2][kN+6] = -1.;
    } else {
      //neutral particle
      mD[fNc+0+k2][kN+0]  = delY ;
      mD[fNc+0+k2][kN+1]  = -(delX) ;
      mD[fNc+0+k2][kN+2]  = 0. ;
      mD[fNc+0+k2][kN+3]  = 0. ;
      mD[fNc+0+k2][kN+4]  =  py;
      mD[fNc+0+k2][kN+5]  = -px;
      mD[fNc+0+k2][kN+6]  = 0. ;

      mD[fNc+1+k2][kN+0] = 2*(delX*px+delY*py)*px*pz/(pT_2*pT_2) - pz*delX/(pT_2);
      mD[fNc+1+k2][kN+1] = 2*(delX*px+delY*py)*py*pz/(pT_2*pT_2) - pz*delY/(pT_2);
      mD[fNc+1+k2][kN+2] =-(delX*px+delY*py)/(pT_2);
      mD[fNc+1+k2][kN+3] = 0.;
      mD[fNc+1+k2][kN+4] = px*pz/pT_2;
      mD[fNc+1+k2][kN+5] = py*pz/pT_2;
      mD[fNc+1+k2][kN+6] = -1.;
    }
    //DMat_trk.push_back(DMat_tmp);
    //E jacobian matrix
    if(ch !=0 )   {
      mE[fNc+0+k2][0] = -(py + a*delX);
      mE[fNc+0+k2][1] =  (px - a*delY);
      mE[fNc+0+k2][2] = 0.;
      mE[fNc+1+k2][0] = -px*pz*S;
      mE[fNc+1+k2][1] = -py*pz*S;
      mE[fNc+1+k2][2] = 1.;
    } else {
      mE[fNc+0+k2][0] = -py ;
      mE[fNc+0+k2][1] =  px;
      mE[fNc+0+k2][2] = 0.;
      mE[fNc+1+k2][0] = -px*pz/pT_2;
      mE[fNc+1+k2][1] = -py*pz/pT_2;
      mE[fNc+1+k2][2] = 1.;
    }
    if(ch !=0 )   {
      md[fNc+0+k2][0] = delY*px - delX*py - (a/2.)*(delX*delX + delY*delY);
      md[fNc+1+k2][0] = delZ - (pz/a)*asin_J;
    } else {
      md[fNc+0+k2][0] = delY*px - delX*py;
      md[fNc+1+k2][0] = delZ - pz*(delX * px + delY * py)/(pT_2);
    }

  }
  fNc +=2*nd;




}

void RhoKinVtxFitter::ReadMassKinMatrix ( )

private

Definition at line 615 of file RhoKinVtxFitter.cxx.

References a, al1, Double_t, RhoFitterBase::fDaughters, fMass, fNc, fNcon, fNpar, RhoCalculationTools::GetBz(), m, mD, md, mE, p1, pz, CAMath::Sqrt(), and vtx_ex.

Referenced by Compute().

{
  // if(m_fitIncludingVertex == 0)
  //{
  int nd=fDaughters.size();

  mD.ResizeTo(fNcon,fNpar);
  mE.ResizeTo(fNcon,3);
  md.ResizeTo(fNcon,1);

  double Etot = 0.;
  double Px = 0.;
  double Py = 0.;
  double Pz = 0.;

  TMatrixD al1p(al1);
  double a=0; //TODO: is that right to initialize with 0?
  TMatrixD m(nd,1);

  /*
   // Mass Constraint without vertex Info.............................

   for(unsigned k=0;k<nd;++k){
   int kN=k*7;
   TLorentzVector p1=fDaughters[k]->P4();
   m[k][0]=p1.M();
   double px = al1p[kN+0][0];
   double py = al1p[kN+1][0];
   double pz = al1p[kN+2][0];
   double E = TMath::Sqrt(px*px+py*py+pz*pz+m[k][0]*m[k][0]);

   // Etot += E; Px   +=px; Py +=py;  Pz += pz;}
   // md[fNc+0][0] = Etot*Etot - Px*Px - Py*Py - Pz*Pz - fMass*fMass ;

   a = -0.00299792458*2.0*fDaughters[k]->GetCharge();
   Double_t invE = 1./E;

   mD[fNc+0][kN+0] = 2.*(Etot*px*invE-Px);
   mD[fNc+0][kN+1] = 2.*(Etot*py*invE-Py);
   mD[fNc+0][kN+2] = 2.*(Etot*pz*invE-Pz);
   //    mD[fNc+0][kN+3] = 0.0;
   mD[fNc+0][kN+3] = 2* m[k][0]*Etot*invE;
   mD[fNc+0][kN+4] = 2.*Py*a;
   mD[fNc+0][kN+5] = -2.*Px*a;
   mD[fNc+0][kN+6] = 0.0;

   //................Simple....................

   //              mD[fNc+0][kN+0] = -2.*Px;
   //              mD[fNc+0][kN+1] = -2.*Py;;
   //              mD[fNc+0][kN+2] = -2.*Pz;
   //              mD[fNc+0][kN+3] = 2.*Etot;
   //      mD[fNc+0][kN+3] = 2* m[k][0]*Etot*invE;
   //      mD[fNc+0][kN+4] = 2.*(Etot*py*invE-Py)*a;
   //      mD[fNc+0][kN+5] = 2.*(Etot*px*invE-Px)*a;
   //              mD[fNc+0][kN+4] = 0.0;
   //              mD[fNc+0][kN+5] = 0.0;
   //              mD[fNc+0][kN+6] = 0.0;
   }
   */

  // With Vertex Info.............................................
  for(int k=0; k<nd; ++k) {
    int kN=k*7;
    TLorentzVector p1=fDaughters[k]->P4();
    m[k][0]=p1.M();
    double delX = vtx_ex[0][0] - al1p[kN+4][0];
    double delY = vtx_ex[1][0] - al1p[kN+5][0];
    //    double delX=0.;
    //    double delY=0.;

    al1p[kN+0][0] = al1p[kN+0][0]-a*delY;
    al1p[kN+1][0] = al1p[kN+1][0]+a*delX;
    //  al1p[kN+2][0] = al1p[kN+2][0];
    double E = TMath::Sqrt(al1p[kN+0][0]* al1p[kN+0][0]+al1p[kN+1][0]*al1p[kN+1][0]+al1p[kN+2][0]*al1p[kN+2][0]+m[k][0]*m[k][0]);
    Etot += E;

    Px += al1p[kN+0][0] ;
    Py += al1p[kN+1][0];
    Pz += al1p[kN+2][0];
  }
  md[fNc+0][0] = Etot*Etot - Px*Px - Py*Py - Pz*Pz - fMass*fMass ;

  double sumA=0;
  for(int k=0; k<nd; ++k) {
    int kN=k*7;
    double px = al1p[kN+0][0];
    double py = al1p[kN+1][0];
    double pz = al1p[kN+2][0];
    //    double E = al1p[kN+3][0];
    double E = TMath::Sqrt(px*px+py*py+pz*pz+m[k][0]*m[k][0]);

    //here there should be implemented the algorithm for neutral particles

    Double_t bField = 0.1*RhoCalculationTools::GetBz(fDaughters[k]->Pos()); // T, assume field in z only
    a = -0.00299792458*bField*fDaughters[k]->GetCharge();//TODO BField
    sumA += a;
    Double_t invE = 1./E;

    mD[fNc+0][kN+0] = 2.*(Etot*al1p[kN+0][0]*invE-Px);
    mD[fNc+0][kN+1] = 2.*(Etot*al1p[kN+1][0]*invE-Py);
    mD[fNc+0][kN+2] = 2.*(Etot*al1p[kN+2][0]*invE-Pz);
    mD[fNc+0][kN+3] = 0.;
    mD[fNc+0][kN+4] =-2.*(Etot*al1p[kN+1][0]*invE-Py)*a;
    mD[fNc+0][kN+5] = 2.*(Etot*al1p[kN+0][0]*invE-Px)*a;
    mD[fNc+0][kN+6] = 0.;
  }
  cout << "Sum A: " << sumA << endl;


  mE[fNc+0][0] = 2*sumA*Px;
  mE[fNc+0][1] = -2*sumA*Py;
  mE[fNc+0][2] = 0.;
  fNc +=1;
}

void RhoKinVtxFitter::ReadMatrix ( )

private

Definition at line 448 of file RhoKinVtxFitter.cxx.

References al0, RhoFitterBase::fDaughters, i, jj, m, p1, p2, sqrt(), and V_al0.

Referenced by Compute().

{
  int nd =fDaughters.size();
  TMatrixD m(nd,1);
  for (int k=0; k<nd; k++) {
    int kN=k*7;
    //px,py,pz,E,x,y,z
    TLorentzVector p1=fDaughters[k]->P4(); //4-momentum of the daughter
    TVector3 p2=fDaughters[k]->Pos(); // position of the daughter
    al0[kN+0][0]=p1.X();
    al0[kN+1][0]=p1.Y();
    al0[kN+2][0]=p1.Z();
    //    al0[kN+3][0]=p1.E();
    al0[kN+4][0]=p2.X();
    al0[kN+5][0]=p2.Y();
    al0[kN+6][0]=p2.Z();

    double fm=fDaughters[k]->Mass();
    al0[kN+3][0]=sqrt(al0[kN+0][0]*al0[kN+0][0]+ al0[kN+1][0]*al0[kN+1][0]+al0[kN+2][0]*al0[kN+2][0]+fm*fm);

    // Read Covariance Matrix .... Can read 6x6 matrices..................
    TMatrixD p1Cov(7,7);
    TMatrixD p3Cov(6,6); //Why 6x6 here if 7x7 should be cpoied (below)
    TMatrixD p2Cov(7,7);
    TMatrixD p4Cov(7,7);
    p1Cov=fDaughters[k]->Cov7(); //Cov Matrix x,y,z,px,py,pz,E

    for (int ii=0; ii<6; ii++) {for(int jj=0; jj<6; jj++) {p3Cov[ii][jj]=p1Cov[ii][jj];}} //test

    //Extend matrix for energy for each candidates .....6x6 to 7x7
    TMatrixD J(7,6) ;
    J.Zero();
    TMatrixD J_t(6,7);
    for (int ii=0; ii<6; ii++) {for(int jj=0; jj<6; jj++) {J[ii][jj] = 1;}}
    for(int i=3; i<6; ++i) {J[6][i] = al0[kN+i-3][0]/al0[kN+3][0];}
    //   p2Cov= J*p3Cov*(J_t.Transpose(J));
    p2Cov=p1Cov;
    //Change to px,py,pz,E,x,y,z
    for(int i=0; i<7; i++) {
      for(int j=0; j<7; j++) {
        if(i>=3) {
          if(j>=3) {
            p4Cov[i-3][j-3] = p2Cov[i][j];
          } else { p4Cov[i-3][j+4] = p2Cov[i][j]; }
        } else {
          if(j>=3) {
            p4Cov[i+4][j-3] = p2Cov[i][j];
          } else { p4Cov[i+4][j+4] = p2Cov[i][j]; }
        }
      }
    }

    for(int i=0; i<7; i++) {
      for (int j=0; j<7; j++) {
        V_al0[k*7+i][k*7+j]  = p4Cov[i][j];
      }
    }
  }
}

void RhoKinVtxFitter::ResetMatrices ( )

private

Definition at line 86 of file RhoKinVtxFitter.cxx.

References al0, al1, covC, mPull, V_al0, V_al1, vtx_ex, and vtx_st.

Referenced by Compute().

{
  al0.Zero();
  V_al0.Zero();
  al1.Zero();
  V_al1.Zero();
  mPull.Zero();
  vtx_ex.Zero();
  vtx_st.Zero();
  covC.Zero();
}

void RhoFitterBase::SetDaugthersFromComposite ( RhoCandidate *  cand )

protectedinherited

Definition at line 136 of file RhoFitterBase.cxx.

References RhoCandidate::Daughter(), RhoFitterBase::fDaughters, RhoCandidate::IsComposite(), and RhoCandidate::NDaughters().

Referenced by RhoKinHyperonFitter::Fit(), FitNode(), RhoKinHyperonVtxFitter::FitNode(), and RhoKinHyperonFitter::SetMatrices().

{
  fDaughters.clear();
  if(cand->IsComposite()){
    RhoCandidate* tc;
    for (int k=0;k<cand->NDaughters();k++) {
      tc=cand->Daughter(k);
      fDaughters.push_back(tc);
    }
  }
  return;
}

void RhoFitterBase::SetDecayVertex ( RhoCandidate *  composite, const TVector3 &  vtx, const TMatrixD &  CovVV  )

protectedinherited

Definition at line 178 of file RhoFitterBase.cxx.

References RhoCandidate::SetDecayVtx().

Referenced by RhoKalmanVtxFitter::Calculate(), SetOutput(), and RhoKinHyperonVtxFitter::SetOutput().

{
  RhoError decaypointcov(CovVV);
  RhoVector3Err decayvertex(vtx,decaypointcov);
  composite->SetDecayVtx(decayvertex);
}

void RhoFitterBase::SetFourMomentumByDaughters ( RhoCandidate *  composite )

protectedinherited

Definition at line 160 of file RhoFitterBase.cxx.

References RhoCandidate::Cov7(), RhoCandidate::Daughter(), RhoCandidate::IsComposite(), RhoCandidate::IsLocked(), RhoCandidate::NDaughters(), RhoCandidate::P4(), RhoCandidate::SetCov7(), and RhoCandidate::SetP4().

Referenced by RhoKalmanVtxFitter::Calculate(), Rho4CFitter::Fit(), Rho4CFitter::FitConserveMasses(), SetOutput(), and RhoKinFitter::SetOutput().

{
  RhoCandidate* tc;
  TLorentzVector tmpLV;
  TMatrixD tmpCov(7,7);
  // Sum daughter fourmomenta, dive into nodes if necessary
  for (int k=0;k<composite->NDaughters();k++) {
    tc=composite->Daughter(k);
    if (tc->IsComposite() && !tc->IsLocked()) { SetFourMomentumByDaughters(tc); }
    tmpLV += tc->P4();
    tmpCov = tmpCov + tc->Cov7();
  }
  composite->SetP4(tmpLV);
  composite->SetCov7(tmpCov);
  //std::cout<<" Base fitter cov7 from tc "<<tc->Uid()<<"/"<<tc->Charge()<<"/"<<tc->PdgCode()<<": ";tmpCov.Print();
  return;
}

void RhoKinVtxFitter::SetMatrices ( )

private

Definition at line 63 of file RhoKinVtxFitter.cxx.

References al0, al1, covC, RhoFitterBase::fDaughters, fNc, fNcon, fNiter, fNpar, fNpart, fNvar, mPull, NumCon, V_al0, V_al1, vtx_ex, and vtx_st.

Referenced by Compute().

{
  int nd=fDaughters.size();

  fNvar=7;
  fNpart=nd;
  fNpar =nd*fNvar;
  fNcon=NumCon;
  //  if(fVerbose) cout << fNcon << "Num " << endl;
  fNc=0;
  fNiter=0;

  al0.ResizeTo(7*nd,1);
  V_al0.ResizeTo(fNpar,fNpar);
  al1.ResizeTo(7*nd,1);
  V_al1.ResizeTo(fNpar,fNpar);
  vtx_ex.ResizeTo(3,1);
  vtx_st.ResizeTo(3,1);
  mPull.ResizeTo(7*nd,1);
  covC.ResizeTo(7,7);
}

void RhoKinVtxFitter::SetMinDChisq ( double  m = 0.001 )

inline

Definition at line 33 of file RhoKinVtxFitter.h.

References fabs(), fMinDChisq, and m.

{fMinDChisq=fabs(m);};

void RhoKinVtxFitter::SetNIterationsExact ( int  nit = 2 )

inline

Definition at line 32 of file RhoKinVtxFitter.h.

References fIterateExact, and fNMaxIterations.

Referenced by PndPmtTask::Exec().

{fNMaxIterations=nit;fIterateExact=true;};

void RhoKinVtxFitter::SetNMaxIterations ( int  nit = 20 )

inline

Definition at line 31 of file RhoKinVtxFitter.h.

References fIterateExact, and fNMaxIterations.

{fNMaxIterations=nit;fIterateExact=false;};

void RhoKinVtxFitter::SetOutput ( RhoCandidate *  head )

private

[ralfk:28.5.2013] Use flat Fourmomentum sum from RhoFitterBase

Definition at line 329 of file RhoKinVtxFitter.cxx.

References a, al0, al1, covC, Double_t, RhoFitterBase::fDaughters, RhoFitterBase::fVerbose, RhoCalculationTools::GetBz(), i, m, p1, pos, RhoFitterBase::SetDecayVertex(), RhoFitterBase::SetFourMomentumByDaughters(), RhoCandidate::SetPos(), V_al0, and vtx_ex.

Referenced by FitNode().

{

  int nd=fDaughters.size();
  TMatrixD m(nd,1);


  double  sumA=0;
  double a;
  for (int k=0; k<nd; k++) {


    //skip locked daughters
    if(fDaughters[k]->IsLocked()) continue;
    Double_t bField = 0.1*RhoCalculationTools::GetBz(fDaughters[k]->Pos()); // T, assume field in z only
    a = -0.00299792458*bField*fDaughters[k]->GetCharge();
    sumA += a;
    TVector3 pos(al0[k*7+4][0],al0[k*7+5][0],al0[k*7+6][0]);
//std::cout<<" --"<<k<<"-- ("<<pos.x()<<";"<<pos.y()<<";"<<pos.z()<<")"<<std::endl;

    TLorentzVector mom4(al0[k*7+0][0],al0[k*7+1][0],al0[k*7+2][0], al0[k*7+3][0]);
//    fDaughters[k]->SetP7(pos,mom4);
//better to put daugthers with mass hypothesis .......?? VJ
    TLorentzVector momM;
    double fM=fDaughters[k]->Mass();
    momM.SetXYZM(al0[k*7+0][0],al0[k*7+1][0],al0[k*7+2][0],fM);
//    momM.SetP4(fM,al0[k*7+0][0],al0[k*7+1][0],al0[k*7+2][0]);
    fDaughters[k]->SetP7(pos,momM);


    //Extend matrix for energy for each candidates if daughters from mass hypothesis 6x6 covariance
    TMatrixD p1Cov(7,7);
    TLorentzVector p1=fDaughters[k]->P4();
    TMatrixD p2Cov(7,7);

    for(int i=0; i<7; i++) {
      for (int j=0; j<7; j++) {
        p1Cov[i][j]= V_al0[k*7+i][k*7+j];
      }
    }


    //Change from px,py,pz,E,x,y,z
    //         to x,y,z,px,py,pz,E
    for(int i=0; i<7; i++) {
      for(int j=0; j<7; j++) {
        if(i>=4) {
          if(j>=4) {
            p2Cov[i-4][j-4] = p1Cov[i][j];
          } else { p2Cov[i-4][j+3] = p1Cov[i][j]; }
        } else {
          if(j>=4) {
            p2Cov[i+3][j-4] = p1Cov[i][j];
          } else { p2Cov[i+3][j+3] = p1Cov[i][j]; }
        }
      }
    }

    // create cov with E... check it
    double invE = 1./al0[k*7+3][0];
    p2Cov[3][6] = p2Cov[6][3] = (p1.X()*p1Cov[0][0]+p1.Y()*p1Cov[0][1]+p1.Z()*p1Cov[0][2])*invE;
    p2Cov[4][6] = p2Cov[6][4] = (p1.X()*p1Cov[1][0]+p1.Y()*p1Cov[1][1]+p1.Z()*p1Cov[1][2])*invE;
    p2Cov[5][6] = p2Cov[6][5] = (p1.X()*p1Cov[2][0]+p1.Y()*p1Cov[2][1]+p1.Z()*p1Cov[2][2])*invE;

    p2Cov[6][6] = (p1.X()*p1.X()*p1Cov[0][0]+p1.Y()*p1.Y()*p1Cov[1][1]+p1.Z()*p1.Z()*p1Cov[2][2]
                       +2.0*p1.X()*p1.Y()*p1Cov[0][1]
                       +2.0*p1.X()*p1.Z()*p1Cov[0][2]
                       +2.0*p1.Y()*p1.Z()*p1Cov[1][2])*invE*invE;

    p2Cov[6][0] = p2Cov[0][6] = (p1.X()*p1Cov[0][4]+p1.Y()*p1Cov[1][4]+p1.Z()*p1Cov[2][4])*invE;
    p2Cov[6][1] = p2Cov[1][6] = (p1.X()*p1Cov[0][5]+p1.Y()*p1Cov[1][5]+p1.Z()*p1Cov[2][5])*invE;
    p2Cov[6][2] = p2Cov[2][6] = (p1.X()*p1Cov[0][6]+p1.Y()*p1Cov[1][6]+p1.Z()*p1Cov[2][6])*invE;

    fDaughters[k]->SetCov7(p2Cov); //New covariance matrix with correlations
    //cout<< " #######  KinVtx daughter cov check... " << endl;
    //cout<<"p1Cov"; p1Cov.Print();
    //cout<<"p2Cov"; p2Cov.Print();

  }


//
//  // For the composite particle ..............................
//  double fpx=0,fpy=0,fpz=0,fe=0;
//  for (int k=0; k<nd; k++) {
//    fpx+= al0[k*7+0][0]-a*(vtx_ex[1][0]*sumA/a - al0[k*7+5][0]);
//    fpy+= al0[k*7+1][0]+a*(vtx_ex[0][0]*sumA/a - al0[k*7+4][0]);
//    fpz+= al0[k*7+2][0];
//    fe += al0[k*7+3][0];
//  }
// //         double TotE=(fpx*fpx+fpy*fpy+fpz*fpz+fe*fe);
//  //double fM=sqrt(fe*fe-(fpx*fpx+fpy*fpy*fpz*fpz));
//  TLorentzVector sum(fpx,fpy,fpz,fe);
//  //          TLorentzVector sum;
//  //          sum.SetXYZM(fpx,fpy,fpz,fM);
//  TVector3 vtx(vtx_ex[0][0],vtx_ex[1][0],vtx_ex[2][0]);
//  //fHeadOfTree->SetP7(vtx,sum);
//  head->SetP7(vtx,sum); //[ralfk:01.12.11 Try to make it a leaf-by-leaf fit]
//  //fHeadOfTree->SetCov7(covC); //New covariance matrix
//  head->SetCov7(covC); //New covariance matrix //[ralfk:01.12.11]
  //We set the decay vertex of the mother! [R.K.]

  TVector3 vtx(vtx_ex[0][0],vtx_ex[1][0],vtx_ex[2][0]);
  TMatrixD CovV = covC.GetSub(0,2,0,2);

  head->SetPos(vtx);//P4 is defined here
  SetDecayVertex(head,vtx,CovV);
  SetFourMomentumByDaughters(head);//propagates cov7 from daughters
  //head->SetCov7(covC);//which one to use???
  //cout<<" KinVtx Cov7: ";covC.Print();
  if(fVerbose) { cout<<"Final vertex Position is "<<vtx_ex[0][0]<<" "<<vtx_ex[1][0]<<" "<<vtx_ex[2][0]<<endl; }
  if(fVerbose) { cout<<"Final Momenta are "<<al0[0][0]<<" "<<al1[1][0]<<" "<<al1[2][0]<<endl; }
}

void RhoFitterBase::SetVerbose ( Bool_t  v = kTRUE )

inlineinherited

Definition at line 52 of file RhoFitterBase.h.

References RhoFitterBase::fVerbose, and v.

{fVerbose=v;}

void RhoKinVtxFitter::TransportToVertex ( TMatrixD &  a_in, TMatrixD &  a_cov_in, TMatrixD &  a_out, TMatrixD &  a_cov_out, TMatrixD &  xref  )

private

Correct one helix and/or track(s) to vertex point for charged and/or neutral particles

Definition at line 816 of file RhoKinVtxFitter.cxx.

References a, atan2(), Double_t, fabs(), RhoFitterBase::fDaughters, RhoCalculationTools::GetBz(), ptot, pz, s, sqrt(), x, y, and z.

Referenced by Compute().

{
        //edited by J.Puetz
        //added parametrization for neutral daughter particles

        int fNDau=fDaughters.size();
        int nd=fNDau;
        TMatrixD U(7*nd,7*nd);
        int kN=0;
        for(int k=0; k<nd; k++) {
                kN=7*k;

                //double m = fDaughters[k]->Mass(); //[R.K. 01/2017] unused variable?
                //check, if daughter particle is either neutral or charged
                if (fabs(fDaughters[k]->GetCharge())<1e-6){//begin neutral
                        double a = 0.00299792458; // unit conversion constant
      double inva = 1/a;
                        //Get position, energy and momentum for the daughter particle
                        double px=a_in[kN+0][0];
                        double py=a_in[kN+1][0];
                        double pz=a_in[kN+2][0];
                        //double p2=px*px+py*py+pz*pz; //[R.K. 01/2017] unused variable?
                        //double E= sqrt(m*m+p2);
                        double x=a_in[kN+4][0];
                        double y=a_in[kN+5][0];
                        double z=a_in[kN+6][0];

                        //correct particle position to vertex, keep momentum and energy
                        a_out[kN+0][0] = px;
                        a_out[kN+1][0] = py;
                        a_out[kN+2][0] = pz;
                        a_out[kN+3][0] = a_in[kN+3][0];
                        a_out[kN+4][0] = x+px*inva;
                        a_out[kN+5][0] = y+py*inva;
                        a_out[kN+6][0] = z+pz*inva;

                        //matrix U corrects the covariant matrix a_cov_in to a_cov_out= U * a_cov_in * U^T
                    U[kN+0][kN+0] = 1.;
                    U[kN+1][kN+1] = 1.;
                    U[kN+2][kN+2] = 1.;
                    U[kN+3][kN+3] = 1.;
                   
                    U[kN+4][kN+4] = 1.;
                    U[kN+4][kN+0] = inva;
                   
                    U[kN+5][kN+5] = 1.;
                    U[kN+5][kN+1] = inva;
                   
                    U[kN+6][kN+6] = 1.;
                    U[kN+6][kN+2] = inva;


                }//end neutral

                else{ // if particle is charged

                        Double_t bField = 0.1*RhoCalculationTools::GetBz(fDaughters[k]->Pos()); // T, assume field in z only
                        double a = -0.00299792458*bField*fDaughters[k]->GetCharge();

                        //   if(fVerbose)cout << "a " << a << endl;

                        //Get position, energy and momentum for the daughter particle

                        double px=a_in[kN+0][0];
                        double py=a_in[kN+1][0];
                        double pz=a_in[kN+2][0];
                        double x=a_in[kN+4][0];
                        double y=a_in[kN+5][0];
                        double z=a_in[kN+6][0];

                        double ptot = sqrt(px*px + py*py + pz*pz);
                        //double E=sqrt(m*m+ptot*ptot);
                        double rho  = a/ptot; //1/R with R: the radius of the trajectory
                        double A1 = 1 - pow(pz/ptot,2) - ( (x-xref[0][0])*py - (y-xref[1][0])*px )*rho/ptot ;
                        double A2 = (x-xref[0][0])*px + (y-xref[1][0])*py;
                        A2 = A2/ptot;


                        double det  = sqrt(A1*A1+rho*rho*A2*A2);
                        double cos_rho_s  =      A1/det;
                        double sin_rho_s  = -rho*A2/det;



                        //double s = atan2(sin_rho_s,cos_rho_s);
                        double s = atan2(sin_rho_s,cos_rho_s)/rho ;
                        a_out[kN+0][0] = px*cos_rho_s-py*sin_rho_s;
                        a_out[kN+1][0] = py*cos_rho_s+px*sin_rho_s;
                        a_out[kN+2][0] = pz;
                        a_out[kN+3][0] = a_in[kN+3][0];
                        a_out[kN+4][0] = x + (px*sin_rho_s - py*(1-cos_rho_s))/a;
                        a_out[kN+5][0] = y + (py*sin_rho_s + px*(1-cos_rho_s))/a;
                        a_out[kN+6][0] = z + (pz/ptot)*s;


                        //matrix to calculate the corrected covariance matrix
                        U[kN+0][kN+0] =  cos_rho_s;
                        U[kN+0][kN+1] = -sin_rho_s;

                        U[kN+1][kN+0] = sin_rho_s;
                        U[kN+1][kN+1] = cos_rho_s;

                        U[kN+2][kN+2] = 1.;
                        U[3+kN][3+kN] = 1.;

                        U[4+kN][0+kN] = sin_rho_s/a;
                        U[4+kN][1+kN] = (1.-cos_rho_s)/a;
                        U[4+kN][4+kN] = 1.;

                        U[5+kN][0+kN] = (1.-cos_rho_s)/a;
                        U[5+kN][1+kN] = sin_rho_s/a;
                        U[5+kN][5+kN] = 1.;

                        U[6+kN][2+kN] = s/ptot;
                        U[6+kN][6+kN] = 1.;


                }// end charged

        }//end k

    TMatrixD U_t=U;
    U_t=U_t.T();
    a_cov_out = U*a_cov_in*U_t;

}

Member Data Documentation

TMatrixD RhoKinVtxFitter::al0

private

Definition at line 55 of file RhoKinVtxFitter.h.

Referenced by Compute(), ReadMatrix(), ResetMatrices(), SetMatrices(), and SetOutput().

TMatrixD RhoKinVtxFitter::al1

private

Definition at line 56 of file RhoKinVtxFitter.h.

Referenced by Compute(), ReadKinMatrix(), ReadMassKinMatrix(), ResetMatrices(), SetMatrices(), and SetOutput().

TMatrixD RhoKinVtxFitter::covC

private

Definition at line 59 of file RhoKinVtxFitter.h.

Referenced by Compute(), ResetMatrices(), SetMatrices(), and SetOutput().

double RhoKinVtxFitter::fChi2Diff

private

Definition at line 81 of file RhoKinVtxFitter.h.

double RhoFitterBase::fChiSquare

protectedinherited

Definition at line 74 of file RhoFitterBase.h.

Referenced by RhoKalmanVtxFitter::Calculate(), DecayTreeFitter::Fitter::chiSquare(), Compute(), RhoKinHyperonVtxFitter::Compute(), Rho4CFitter::Do4CFit(), Rho4CFitter::Do4CFitWithMassConservation(), RhoKalmanVtxFitter::FitNode(), RhoKalmanVtxFitter::FitVertexFast(), RhoFitterBase::GetChi2(), RhoFitterBase::GetProb(), Rho4CFitter::Rho4CFitter(), RhoKinFitter::Solve(), and RhoKinHyperonFitter::Solve().

double RhoKinVtxFitter::fchiSquare

private

Definition at line 103 of file RhoKinVtxFitter.h.

std::vector<RhoCandidate*> RhoFitterBase::fDaughters

protectedinherited

Definition at line 69 of file RhoFitterBase.h.

Referenced by Compute(), RhoKinHyperonVtxFitter::Compute(), Rho4CFitter::Do4CFit(), Rho4CFitter::Do4CFitWithMassConservation(), RhoFitterBase::FindAndAddFinalStateDaughters(), Rho4CFitter::Fit(), RhoKinFitter::Fit(), RhoKinHyperonFitter::Fit(), Rho4CFitter::FitConserveMasses(), GetCovariance(), RhoKinHyperonVtxFitter::GetCovariance(), RhoKinFitter::Read4MomKinMatrix(), RhoKinHyperonFitter::Read4MomKinMatrix(), ReadKinMatrix(), RhoKinHyperonVtxFitter::ReadKinMatrix(), ReadMassKinMatrix(), RhoKinFitter::ReadMassKinMatrix(), RhoKinHyperonVtxFitter::ReadMassKinMatrix(), ReadMatrix(), RhoKinFitter::ReadMatrix(), RhoKinHyperonVtxFitter::ReadMatrix(), RhoKinFitter::ReadMomKinMatrix(), RhoKinHyperonFitter::ReadMomKinMatrix(), RhoKinFitter::ReadTotEKinMatrix(), RhoKinHyperonFitter::ReadTotEKinMatrix(), RhoKinFitter::ReadTotMomKinMatrix(), RhoKinHyperonFitter::ReadTotMomKinMatrix(), RhoFitterBase::SetDaugthersFromComposite(), SetMatrices(), RhoKinFitter::SetMatrices(), RhoKinHyperonFitter::SetMatrices(), RhoKinHyperonVtxFitter::SetMatrices(), SetOutput(), RhoKinFitter::SetOutput(), RhoKinHyperonFitter::SetOutput(), RhoKinHyperonVtxFitter::SetOutput(), RhoKinHyperonVtxFitter::TransportToPoca(), TransportToVertex(), and RhoKinHyperonVtxFitter::TransportToVertex().

RhoCandidate* RhoFitterBase::fHeadOfTree

protectedinherited

Definition at line 62 of file RhoFitterBase.h.

Referenced by Rho4CFitter::Fit(), RhoKinFitter::Fit(), RhoKinHyperonFitter::Fit(), RhoFitterBase::Fit(), RhoFitterBase::FitAll(), Rho4CFitter::FitConserveMasses(), RhoFitterBase::HeadOfTree(), RhoKinHyperonFitter::ReadMatrix(), RhoFitterBase::RhoFitterBase(), RhoKalmanVtxFitter::RhoKalmanVtxFitter(), RhoKinHyperonFitter::SetMatrices(), RhoKinFitter::SetOutput(), and RhoKinHyperonFitter::SetOutput().

bool RhoKinVtxFitter::fIterateExact

private

Definition at line 99 of file RhoKinVtxFitter.h.

Referenced by Compute(), RhoKinVtxFitter(), SetNIterationsExact(), and SetNMaxIterations().

double RhoKinVtxFitter::fMass

private

Definition at line 90 of file RhoKinVtxFitter.h.

Referenced by AddMassConstraint(), and ReadMassKinMatrix().

int RhoKinVtxFitter::fMassConstraint

private

Definition at line 93 of file RhoKinVtxFitter.h.

Referenced by AddMassConstraint(), Compute(), and RhoKinVtxFitter().

double RhoKinVtxFitter::fMinDChisq

private

Definition at line 97 of file RhoKinVtxFitter.h.

Referenced by Compute(), RhoKinVtxFitter(), and SetMinDChisq().

Int_t RhoKinVtxFitter::fNc

private

Definition at line 73 of file RhoKinVtxFitter.h.

Referenced by Compute(), ReadKinMatrix(), ReadMassKinMatrix(), and SetMatrices().

Int_t RhoKinVtxFitter::fNcon

private

Definition at line 70 of file RhoKinVtxFitter.h.

Referenced by ReadKinMatrix(), ReadMassKinMatrix(), and SetMatrices().

int RhoFitterBase::fNDegreesOfFreedom

protectedinherited

Definition at line 75 of file RhoFitterBase.h.

Referenced by RhoKalmanVtxFitter::Calculate(), Compute(), RhoKinHyperonVtxFitter::Compute(), Rho4CFitter::Do4CFit(), Rho4CFitter::Do4CFitWithMassConservation(), RhoKalmanVtxFitter::FitVertexFast(), RhoFitterBase::GetNdf(), RhoFitterBase::GetProb(), RhoKinFitter::SetOutput(), and RhoKinHyperonFitter::SetOutput().

int RhoKinVtxFitter::fnDof

private

Definition at line 102 of file RhoKinVtxFitter.h.

Int_t RhoKinVtxFitter::fNiter

private

Definition at line 74 of file RhoKinVtxFitter.h.

Referenced by SetMatrices().

int RhoKinVtxFitter::fNMaxIterations

private

Definition at line 98 of file RhoKinVtxFitter.h.

Referenced by Compute(), RhoKinVtxFitter(), SetNIterationsExact(), and SetNMaxIterations().

Int_t RhoKinVtxFitter::fNpar

private

Definition at line 68 of file RhoKinVtxFitter.h.

Referenced by ReadKinMatrix(), ReadMassKinMatrix(), and SetMatrices().

Int_t RhoKinVtxFitter::fNpart

private

Definition at line 69 of file RhoKinVtxFitter.h.

Referenced by SetMatrices().

Int_t RhoKinVtxFitter::fNumKnown

private

Definition at line 75 of file RhoKinVtxFitter.h.

Int_t RhoKinVtxFitter::fNvar

private

Definition at line 67 of file RhoKinVtxFitter.h.

Referenced by SetMatrices().

double RhoKinVtxFitter::fPull

private

Definition at line 82 of file RhoKinVtxFitter.h.

Referenced by Compute(), and GetPull().

Bool_t RhoFitterBase::fVerbose

protectedinherited

Definition at line 58 of file RhoFitterBase.h.

Referenced by RhoKalmanVtxFitter::Calculate(), Compute(), RhoKinHyperonVtxFitter::Compute(), RhoKalmanVtxFitter::FitVertexFast(), RhoKalmanVtxFitter::RhoKalmanVtxFitter(), RhoKalmanVtxFitter::SetDebug(), SetOutput(), RhoKinHyperonVtxFitter::SetOutput(), RhoKalmanVtxFitter::SetSilent(), and RhoFitterBase::SetVerbose().

TMatrixD RhoKinVtxFitter::mD

private

Definition at line 61 of file RhoKinVtxFitter.h.

Referenced by Compute(), ReadKinMatrix(), and ReadMassKinMatrix().

TMatrixD RhoKinVtxFitter::md

private

Definition at line 63 of file RhoKinVtxFitter.h.

Referenced by Compute(), ReadKinMatrix(), and ReadMassKinMatrix().

TMatrixD RhoKinVtxFitter::mE

private

Definition at line 62 of file RhoKinVtxFitter.h.

Referenced by Compute(), ReadKinMatrix(), and ReadMassKinMatrix().

TMatrixD RhoKinVtxFitter::mPull

private

Definition at line 64 of file RhoKinVtxFitter.h.

Referenced by Compute(), ResetMatrices(), and SetMatrices().

int RhoKinVtxFitter::niter

private

Definition at line 80 of file RhoKinVtxFitter.h.

Int_t RhoKinVtxFitter::NumCon

private

Definition at line 76 of file RhoKinVtxFitter.h.

Referenced by Compute(), and SetMatrices().

TMatrixD RhoKinVtxFitter::V_al0

private

Definition at line 57 of file RhoKinVtxFitter.h.

Referenced by Compute(), ReadMatrix(), ResetMatrices(), SetMatrices(), and SetOutput().

TMatrixD RhoKinVtxFitter::V_al1

private

Definition at line 58 of file RhoKinVtxFitter.h.

Referenced by Compute(), ResetMatrices(), and SetMatrices().

TMatrixD RhoKinVtxFitter::vtx_ex

private

Definition at line 86 of file RhoKinVtxFitter.h.

Referenced by Compute(), ReadKinMatrix(), ReadMassKinMatrix(), ResetMatrices(), SetMatrices(), and SetOutput().

TMatrixD RhoKinVtxFitter::vtx_st

private

Definition at line 87 of file RhoKinVtxFitter.h.

Referenced by Compute(), ResetMatrices(), and SetMatrices().

---

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

• RhoKinVtxFitter.h
• RhoKinVtxFitter.cxx