#include #include "PndVtxPoca.h" #include "RhoBase/TCandListIterator.h" #include "RhoBase/TRho.h" #include "RhoBase/TFactory.h" using namespace std; ClassImp(PndVtxPoca) TBuffer &operator>>(TBuffer &buf, PndVtxPoca *&obj) { obj = (PndVtxPoca *) buf.ReadObject(PndVtxPoca::Class()); return buf; } //Include only those constraint which need vertex Info.... PndVtxPoca::PndVtxPoca( const TCandidate& b) : VAbsFitter( b ) { } PndVtxPoca::~PndVtxPoca() { } Double_t PndVtxPoca::GetPocaVtx(TVector3 &vertex) { vertex.SetXYZ(0.,0.,0.); if ( fHeadOfTree->NDaughters() < 2 ) return 0.; if ( fHeadOfTree->NDaughters() == 2 ) return GetPoca(vertex,fHeadOfTree->Daughter(0),fHeadOfTree->Daughter(1)); std::vector distances; std::vector results; // loop over daughters, take the mean value of all "best" positions // TODO do this smarter by using already found vertices ? TVector3 theVertex(0.,0.,0.); Double_t actualDoca=0.; for(Int_t daug1 =0;daug1NDaughters();daug1++) { TCandidate* a=fHeadOfTree->Daughter(daug1); for(Int_t daug2=daug1+1;daug2NDaughters();daug2++) { TCandidate* b=fHeadOfTree->Daughter(daug2); actualDoca = GetPoca(theVertex,a,b); distances.push_back(actualDoca); results.push_back(theVertex); }//daug2 }//daug1 //TODO invent a smart procedure to get the correct 2-Track doca and something reasonable for many tracks // --> Mean = Sum(x/(sigmax^2)) / Sum(Sigmax^2) // Averaging vertex results from each track pair how to do that? "geometric" or arithmetic mean? std::vector::iterator iterDoca; std::vector::iterator iterVtx; Double_t docaweight=0,sumdocaweigts=0; TVector3 vertexK; for(iterVtx=results.begin(), iterDoca=distances.begin();iterVtx!=results.end()&&iterDoca!=distances.end();++iterVtx,++iterDoca) { docaweight=1/(*iterDoca); //docaweight *= docaweight; vertexK=*iterVtx; if (docaweight == 0) docaweight = 1; // right so? vertexK *= docaweight; vertex+=vertexK; sumdocaweigts+=docaweight; } if (sumdocaweigts == 0) sumdocaweigts=1; vertex*=1./sumdocaweigts; //sumdocaweigts = sqrt(sumdocaweigts); return fHeadOfTree->NDaughters()/sumdocaweigts; } Double_t PndVtxPoca::GetPoca(TVector3 &vertex,TCandidate* a, TCandidate* b) { vertex.SetXYZ(0.,0.,0.); // loop over daughters, take the mean value of all "best" positions Double_t bField=2.0; // TODO: Get Field from RTDB TVector3 dB(0,0,1.0); TVector3 position1 = a->GetPosition(); // Momentum vectors TVector3 ap3 = a->P3(); Double_t pPerp1 = ap3.Perp(); TVector3 d1 = ap3; d1.SetZ(0); d1*=1.0/pPerp1; // Radius and center Double_t rho1 = pPerp1/(0.0029979246*bField); // Radius in cm TVector3 r1=d1.Cross(dB); r1 *= -a->Charge()*rho1; TVector3 center1 = position1 - r1; center1.SetZ(0); TVector3 position2 = b->GetPosition(); // Momentum vectors TVector3 bp3 = b->P3(); Double_t pPerp2 = bp3.Perp(); TVector3 d2 = bp3; d2.SetZ(0); d2*=1.0/pPerp2; // Radius and center Double_t rho2 = pPerp2/(0.0029979246*bField); // Radius in cm TVector3 r2=d2.Cross(dB); r2 *= -b->Charge()*rho2; TVector3 center2 = position2 - r2; center2.SetZ(0); // distance and angle of the axis between the two centers TVector3 ab = center2 - center1; Double_t dab = ab.Perp(); Double_t cosTheAB = ab.X()/dab; Double_t sinTheAB = ab.Y()/dab; // x value of intersect at reduced system Double_t x = dab/2 + ( rho1*rho1 - rho2*rho2 )/(2*dab); // y*y value of intersect at reduced system for helix A Double_t y2 = (rho1+x)*(rho1-x); // both circles do not intersect (only one solution) Int_t nSolMax=1; Double_t y=0; if (y2 > 0) { nSolMax=2; y = sqrt(y2); } // now we compute the solution(s) TVector3 newapos[2]; TVector3 newbpos[2]; Int_t best=0; Double_t actualDoca=1.E8; for (Int_t ns=0; ns0 ? 1.0 : -1.0; Double_t aangle=adir * r1.Angle(rs1); // intersection point Double_t newaz=position1.Z() + rho1*aangle/pPerp1 * ap3.Z(); newapos[ns].SetX( center1.X() + rs1.X() ); newapos[ns].SetY( center1.Y() + rs1.Y() ); newapos[ns].SetZ( newaz ); // same for b Double_t bdir=(rs2-r2).Dot(bp3)>0 ? 1.0 : -1.0; Double_t bangle=bdir * r2.Angle(rs2); Double_t newbz=position2.Z() + rho2*bangle/pPerp2 * bp3.Z(); newbpos[ns].SetX( center2.X() + rs2.X()); // ==newapos[ns].X() newbpos[ns].SetY( center2.Y() + rs2.Y()); // ==newapos[ns].Y() newbpos[ns].SetZ( newbz ); Double_t delta = (newapos[ns]-newbpos[ns]).Mag(); // take the solution of minimal deltaZ if ( delta < actualDoca ) { best=ns; actualDoca = delta; } } vertex=0.5*(newapos[best]+newbpos[best]); return actualDoca; }