#include #include "PndKinFitter.h" #include "RhoBase/TCandListIterator.h" #include "RhoBase/TFactory.h" #include "TDecompLU.h" #include "TMatrixD.h" #include "TMatrixDSym.h" #include "RhoBase/TRho.h" using namespace std; ClassImp(PndKinFitter) TBuffer &operator>>(TBuffer &buf, PndKinFitter *&obj) { obj = (PndKinFitter *) buf.ReadObject(PndKinFitter::Class()); return buf; } PndKinFitter::PndKinFitter( const TCandidate& b) : VAbsFitter( b ) { f4MomConstraint=-1; fMomConstraint=-1; fMassConstraint =-1; fTotMomConstraint=-1; fTotEConstraint=-1; } PndKinFitter::~PndKinFitter() { } void PndKinFitter::Add4MomConstraint(TLorentzVector lv) { f4MomConstraint = 1; flmm=lv; } void PndKinFitter::AddMomConstraint(TVector3 v) { fMomConstraint = 1; fmm=v; } void PndKinFitter::AddTotEConstraint(double energy) { fTotEConstraint = 1; fEc=energy; } void PndKinFitter::AddTotMomConstraint(double momentum) { fTotMomConstraint = 1; fMom=momentum; } void PndKinFitter::AddMassConstraint(double mass) { fMassConstraint = 1; fMass=mass; } void PndKinFitter::Fit() { fDaughters.Cleanup(); FindAndAddGenericDaughters(fHeadOfTree); fHeadOfTree->RemoveAssociations(); for (int i=0;iAddDaughterLinkSimple(&(fDaughters[i])); //int nd=fDaughters.GetLength(); //unused? NumCon=0; if(f4MomConstraint >0){ int n4Mom = 4; NumCon = NumCon + n4Mom; } if(fMomConstraint >0){ int nMom = 3; NumCon = NumCon + nMom; } if(fTotEConstraint >0){ int nE = 1; NumCon = NumCon + nE; } if(fMassConstraint >0){ int nMass = 1; NumCon = NumCon + nMass; } if(fTotMomConstraint >0){ int nTotMom = 1; NumCon = NumCon + nTotMom; } SetMatrices(); ZeroMatrices(); ReadMatrix(); Solve(); if (true) {SetOutput();} } void PndKinFitter::FindAndAddGenericDaughters(TCandidate *head) { TCandidate *tc; TCandListIterator iter=head->DaughterIterator(); while (tc=iter.Next()) { if (!tc->IsComposite()) fDaughters.Add(*tc); else FindAndAddGenericDaughters(tc); } } void PndKinFitter::SetMatrices(){ int nd=fDaughters.GetLength(); fNvar=7; fNpar = nd*fNvar; fNcon= NumCon; fNc=0; fNiter=0; al0.ResizeTo(7*nd,1); V_al0.ResizeTo(fNpar,fNpar); al1.ResizeTo(7*nd,1); V_al1.ResizeTo(fNpar,fNpar); mD.ResizeTo(fNcon,fNpar); mE.ResizeTo(fNcon,3); md.ResizeTo(fNcon,1); mPull.ResizeTo(7*nd,1); } void PndKinFitter::ZeroMatrices(){ al0.Zero(); V_al0.Zero(); al1.Zero(); V_al1.Zero(); mPull.Zero(); mD.Zero(); md.Zero(); mE.Zero(); } void PndKinFitter::Solve() { //int nd=fDaughters.GetLength(); //unused? double ierr; // used to check inversions al1=al0; // int j1Max=50; // for(Int_t j1=0;j10){ Read4MomKinMatrix();} if(fMomConstraint >0){ ReadMomKinMatrix();} if(fTotEConstraint >0){ ReadTotEKinMatrix();} if(fMassConstraint >0){ ReadMassKinMatrix();} if(fTotMomConstraint >0){ ReadTotMomKinMatrix();} TMatrixD mD_t=mD; mD_t.T(); // mD_t=mD_t.Transpose(mD); // mD_t.Print(); TMatrixD Vd_inv = mD*V_al0*mD_t; TMatrixD Vd = Vd_inv.Invert(&ierr); Vd.Print(); // TMatrixD lam=Vd*md; TMatrixD lam = Vd* ( mD*(al1 - al0) + md); TMatrixD al_new=al0-V_al0*mD_t*lam; // al_new.Print(); TMatrixD V_al_new=V_al0-V_al0*mD_t*Vd*mD*V_al0; double chi2=0.; for (int i=0;i 0 ) {mPull[0][0] =(al0[0][0]-al_new[0][0])/sqrt(covdif);} al0=al_new; V_al0=V_al_new; fPull=mPull[0][0]; } //Write output //Write output void PndKinFitter::SetOutput() { int nd=fDaughters.GetLength(); TMatrixD m(nd,1); fdgf = NumCon; double sumA=0; double a; for (int k=0;kSetP4(sum); } //Read the input vector void PndKinFitter::ReadMatrix() { int nd =fDaughters.GetLength(); for (int k=0;k=3){ if(j>=3) p4Cov[i-3][j-3] = p3Cov[i][j];else p4Cov[i-3][j+3] = p3Cov[i][j]; }else {if(j>=3) p4Cov[i+4][j-3] = p3Cov[i][j];else p4Cov[i+4][j+4] = p3Cov[i][j];}}} // cout<<"p2Cov"<