AlgFitTrackSR Class Reference

#include <AlgFitTrackSR.h>

Inheritance diagram for AlgFitTrackSR:

List of all members.

Public Member Functions
	AlgFitTrackSR ()
virtual	~AlgFitTrackSR ()
virtual void	RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
virtual void	Trace (const char *c) const
Private Member Functions
Int_t	Predict (KalmanPlaneSR *prevkp, KalmanPlaneSR *currentkp, Int_t idir)
Int_t	Filter (KalmanPlaneSR *currentkp, Int_t idir)
Int_t	FitFrom (KalmanPlaneSR *prevkp, KalmanPlaneSR *currentkp, Int_t idir, Int_t iterate)
Int_t	CalculatePropagator (KalmanPlaneSR *prevkp, KalmanPlaneSR *currentkp, Int_t idir)
void	CalculateNoise (KalmanPlaneSR *prevkp, KalmanPlaneSR *currentkp, Int_t idir, Int_t befaft=0)
void	CalculatePreCovariance (KalmanPlaneSR *prevkp, KalmanPlaneSR *currentkp, Int_t idir)
Int_t	CalculatePreState (KalmanPlaneSR *prevkp, KalmanPlaneSR *currentkp, Int_t idir)
void	CalculatePreChi2 (KalmanPlaneSR *currentkp, Int_t idir, Double_t=0.0119)
void	CalculateGain (KalmanPlaneSR *currentkp, Int_t idir, Double_t=0.0119)
void	CalculateFilState (KalmanPlaneSR *currentkp, Int_t idir)
void	CalculateFilCovariance (KalmanPlaneSR *currentkp, Int_t idir)
void	CalculateFilChi2 (KalmanPlaneSR *currentkp, Int_t idir, Double_t=0.0119)
Bool_t	Swim (Double_t *swimstate, Double_t *state, Double_t zPos, Double_t finalZ, Int_t idir, const VldContext *)
void	InitKalmanFitParameters (AlgConfig &ac)
void	InitializeFit (CandFitTrackSRHandle &cfh, Int_t, Int_t)
void	MakeSliceClusterList (TObjArray *trackClusterList, const TObjArray *tclist, CandStripHandleItr &stripItr, CandFitTrackSRHandle &cfh)
void	MakeTrackClusterList (TObjArray *trackClusterList, CandStripHandleItr &stripItr, CandFitTrackSRHandle &cfh, Int_t direction)
void	SetTrackParameters (const CandTrackHandle *track0, CandFitTrackSRHandle &cfh, Int_t direction)
Int_t	FindNumSkippedPlanes (Int_t currentPlane, Int_t prevPlane, KalmanPlaneSR *oldkp, Int_t direction)
Int_t	FindTimingDirection (TrackClusterSRItr &clusterItr)
Bool_t	MarkTrackClusters (const CandTrackHandle *track0, TObjArray &planeClusterList, CandStripHandleItr &stripItr, CandFitTrackSRHandle &cfh, Int_t &begPlane, Int_t &endPlane)
void	SetTrackEndParameters (Int_t begPlane, Int_t endPlane, CandFitTrackSRHandle &cfh, const CandTrackHandle *track0)
void	ResetTrackClusterList (TObjArray &planeClusterList)
void	ResetTrackClusterList (TObjArray &planeClusterList, CandFitTrackSRHandle &cfh)
Int_t	DoKalmanFit (TObjArray &planeClusterList, CandFitTrackSRHandle &cfh, Int_t &istatus, Int_t direction, Int_t dosearch)
Int_t	ReverseFit (TObjArray &planeClusterList, CandFitTrackSRHandle &cfh, Int_t, Bool_t=kFALSE)
Int_t	AddClustersToFit (TObjArray &planeClusterList, CandFitTrackSRHandle &cfh, Int_t iterate, Int_t &nu, Int_t &nv, Int_t direction)
Int_t	RemoveBadPointsFromFit (CandFitTrackSRHandle &cfh, TObjArray &planeClusterList, std::map< Int_t, Int_t > &uFitPlanes, std::map< Int_t, Int_t > &vFitPlanes, Int_t &nfitu, Int_t &nfitv, Int_t direction)
Int_t	FindDownstreamPlanes (TObjArray &planeClusterList, CandFitTrackSRHandle &cfh, Int_t direction)
Int_t	FindUpstreamPlanes (TObjArray &planeClusterList, CandFitTrackSRHandle &cfh, TObjArray *newkplist, Int_t direction)
Int_t	FindNumSkippedPlanesInView (std::map< Int_t, Int_t > &fitPlanes, std::map< Int_t, Int_t >::iterator &fitIter, Int_t currentPlane, Int_t prevPlane, Int_t lastPlane, Int_t direction)
Int_t	AddForwardBestKPToFit (KalmanPlaneSR *bestkp, CandFitTrackSRHandle &cfh, Int_t &nswimfail)
Int_t	AddToFit (CandFitTrackSRHandle &cfh, TrackClusterSR *, Int_t)
Int_t	AddToFit (CandFitTrackSRHandle &cfh, KalmanPlaneSR *, Int_t, Bool_t)
Int_t	AddReverseBestKPToFit (KalmanPlaneSR *bestkp, KalmanPlaneSR *oldkp, CandFitTrackSRHandle &cfh, TObjArray *newkplist, Int_t &nswimfail, Int_t direction)
void	SetPlaneParameters (CandFitTrackSRHandle &cfh, const KalmanPlaneSR *kpu0, const KalmanPlaneSR *kpu1, const KalmanPlaneSR *kpv0, const KalmanPlaneSR *kpv1)
void	MakeDaughterStripList (CandFitTrackSRHandle &cfh)
void	SwimVertexAndEndPoints (CandFitTrackSRHandle &cfh, KalmanPlaneSR *kp0, KalmanPlaneSR *kp1, const KalmanPlaneSR *kpu0, const KalmanPlaneSR *kpu1, const KalmanPlaneSR *kpv0, const KalmanPlaneSR *kpv1, Double_t *planepe, Int_t plane0, Int_t plane1, Double_t &vtxqp, Int_t direction)
Int_t	IterateKalmanFit (TObjArray &planeClusterList, CandFitTrackSRHandle &cfh, Int_t &nu, Int_t &nv, Int_t direction, Int_t dosearch)
Int_t	FindTimingDirection (CandFitTrackSRHandle &cfh, Double_t *fitparm, Double_t &timefitchi2)
Private Attributes
VldContext	fVldContext
Detector::Detector_t	fDetector
Int_t	pRev
Int_t	pRev2
Int_t	pFor
Int_t	fParmMaxIterate
Int_t	fParmMaxIterate2
Double_t	fParmQPDiff
Double_t	fParmMinPlanePE
Int_t	fParmLastPlane
Int_t	fParmIsCosmic
Double_t	fParmMaxLocalChi2
Double_t	fParmMaxLocalPreChi2
Double_t	fParmMaxAngleCovariance
Double_t	fParmMaxImpactParameter
Int_t	fParmNSkipActive
Int_t	fParmNSkipView
Double_t	fParmPassReducedChi2
Double_t	fParmPassPlaneAsymmetry
Int_t	fParmPassMinPlaneAsymmetry
Double_t	fParmMinClusterCharge
Int_t	fParmMaxClusterNStrip
Double_t	fParmDeltaChi2
Double_t	fParmDeltaCovariance
Double_t	fParmMisalignmentError
Double_t	fParmDState [5]
Double_t	fParmdedx
Double_t	fParmPlnRadLen
Int_t	fParmQPRangeCheck
Double_t	fParmMaxQP
Double_t	fParmMaxQPFrac
Int_t	fSwimmer
Double_t	fParmMaxAngle
Double_t	fParmTPosError2
Double_t	fParmInitialPositionError2
Double_t	fParmInitialSlopeError2
Double_t	fParmInitialQPError2
Double_t	fParmCovarianceScale
Double_t	fCov [5][5]

---

Constructor & Destructor Documentation

AlgFitTrackSR::AlgFitTrackSR (   )

Definition at line 73 of file AlgFitTrackSR.cxx. { }

AlgFitTrackSR::~AlgFitTrackSR (   )  [virtual]

Definition at line 78 of file AlgFitTrackSR.cxx. { }

---

Member Function Documentation

Int_t AlgFitTrackSR::AddClustersToFit (  TObjArray &  planeClusterList, CandFitTrackSRHandle &  cfh, Int_t  iterate, Int_t &  nu, Int_t &  nv, Int_t  direction )  [private]

Definition at line 2361 of file AlgFitTrackSR.cxx. References AddToFit(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetPlane(), TrackClusterSR::GetPlaneView(), TrackClusterSR::IsValid(), MSG, and tc. Referenced by DoKalmanFit(), and IterateKalmanFit(). { MSG("FitTrackSR", Msg::kDebug) << "in AddClustersToFit" << endl; TrackClusterSR *tc = 0; Int_t istatus = 0; Bool_t dir = kIterForward; if(direction!=1) dir=kIterBackward; TIter planeItr(&planeClusterList,dir); TObjArray * plnarray; while( (plnarray = (TObjArray * )planeItr.Next()) ){ TIter clusterItr(plnarray); while( (tc = (TrackClusterSR *)clusterItr.Next()) && istatus>=0 ){ MSG("FitTrackSR", Msg::kDebug) << "current cluster on plane " << tc->GetPlane() << " is valid " << (Int_t)tc->IsValid() << endl; if(tc->IsValid()){ MSG("FitTrackSR", Msg::kDebug) << "adding cluster on plane " << tc->GetPlane() << " " << tc->GetMinStrip() << "-" << tc->GetMaxStrip() << " to fit" << endl; istatus = AddToFit(cfh,tc, iterate); if(istatus>=0){ if(tc->GetPlaneView()==PlaneView::kU) ++nu; else if(tc->GetPlaneView()==PlaneView::kV) ++nv; } } } }//end loop over clusters return istatus; }

Int_t AlgFitTrackSR::AddForwardBestKPToFit (  KalmanPlaneSR *  bestkp, CandFitTrackSRHandle &  cfh, Int_t &  nswimfail )  [private]

Definition at line 3407 of file AlgFitTrackSR.cxx. References CandFitTrackSRHandle::AddKalmanPlaneAt(), AddToFit(), fParmMaxAngleCovariance, fParmMaxLocalPreChi2, CandFitTrackSRHandle::GetKalmanLast(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), CandFitTrackHandle::GetNSwimFail(), TrackClusterSR::GetPlane(), KalmanPlaneSR::GetPreChi2(), KalmanPlaneSR::GetPreCovarianceMatrixValue(), KalmanPlaneSR::GetTrackCluster(), kdUdZ, kdVdZ, MSG, CandFitTrackSRHandle::RemoveKalmanPlane(), and CandFitTrackHandle::SetNSwimFail(). Referenced by FindDownstreamPlanes(). { MSG("FitTrackSR", Msg::kDebug) << "in AddForwardBestKPToFit" << endl; Int_t nadd = 0; MSG("FitTrackSR",Msg::kDebug) << "considering best tc " << bestkp->GetTrackCluster()->GetPlane() << " " << bestkp->GetTrackCluster()->GetMinStrip() << "-" << bestkp->GetTrackCluster()->GetMaxStrip() << ", chi2 = " << bestkp->GetPreChi2(0) << "/" << fParmMaxLocalPreChi2 << " " << bestkp->GetPreCovarianceMatrixValue(0,kdUdZ,kdUdZ) << "/" << fParmMaxAngleCovariance << " " << bestkp->GetPreCovarianceMatrixValue(0,kdVdZ,kdVdZ) << "/" << fParmMaxAngleCovariance << endl; if(bestkp->GetPreChi2(0)<fParmMaxLocalPreChi2 && bestkp->GetPreCovarianceMatrixValue(0,kdUdZ,kdUdZ)<fParmMaxAngleCovariance && bestkp->GetPreCovarianceMatrixValue(0,kdVdZ,kdVdZ)<fParmMaxAngleCovariance){ KalmanPlaneSR *kpnewer = cfh.AddKalmanPlaneAt(bestkp,cfh.GetKalmanLast()+1); Int_t istatus2 = AddToFit(cfh,kpnewer,-1,0); if (istatus2>=0) { cfh.SetNSwimFail(cfh.GetNSwimFail()+nswimfail); ++nadd; } else { MSG("FitTrackSR",Msg::kDebug) << "AddToFit failed, do not add to fit" << endl; // need to remove from KalmanPlaneList cfh.RemoveKalmanPlane(-1,0); if(bestkp){ delete bestkp; bestkp = 0; } } } else if(bestkp){ delete bestkp; bestkp = 0; } nswimfail=0; return nadd; }

Int_t AlgFitTrackSR::AddReverseBestKPToFit (  KalmanPlaneSR *  bestkp, KalmanPlaneSR *  oldkp, CandFitTrackSRHandle &  cfh, TObjArray *  newkplist, Int_t &  nswimfail, Int_t  direction )  [private]

Definition at line 3463 of file AlgFitTrackSR.cxx. References FitFrom(), fParmMaxAngleCovariance, fParmMaxLocalPreChi2, CandFitTrackSRHandle::GetKalmanLast(), CandFitTrackSRHandle::GetKalmanPlane(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), CandFitTrackHandle::GetNSwimFail(), TrackClusterSR::GetPlane(), KalmanPlaneSR::GetPreChi2(), KalmanPlaneSR::GetPreCovarianceMatrixValue(), KalmanPlaneSR::GetTrackCluster(), kdUdZ, kdVdZ, MSG, CandFitTrackSRHandle::SetCurrentKalmanPlane(), and CandFitTrackHandle::SetNSwimFail(). Referenced by FindUpstreamPlanes(). { MSG("FitTrackSR", Msg::kDebug) << "in AddReverseBestKPToFit" << endl; Int_t nadd = 0; MSG("FitTrackSR",Msg::kDebug) << "considering best tc " << bestkp->GetTrackCluster()->GetPlane() << " " << bestkp->GetTrackCluster()->GetMinStrip() << "-" << bestkp->GetTrackCluster()->GetMaxStrip() << ", chi2 = " << bestkp->GetPreChi2(1) << "/" << fParmMaxLocalPreChi2 << " planediff = " << bestkp->GetTrackCluster()->GetPlane()-oldkp->GetTrackCluster()->GetPlane() << " " << bestkp->GetPreCovarianceMatrixValue(1,kdUdZ,kdUdZ) << " " << bestkp->GetPreCovarianceMatrixValue(1,kdVdZ,kdVdZ) << endl; if(bestkp->GetPreChi2(1)<fParmMaxLocalPreChi2 && bestkp->GetPreCovarianceMatrixValue(1,kdUdZ,kdUdZ)<fParmMaxAngleCovariance && bestkp->GetPreCovarianceMatrixValue(1,kdVdZ,kdVdZ)<fParmMaxAngleCovariance){ // check if forward fit is also good // first find the plane just before bestkp KalmanPlaneSR *foundkp = 0; KalmanPlaneSR *befkp = 0; //look for the kp before the best kp for(Int_t i=cfh.GetKalmanLast(); i>=0 && !foundkp; --i){ befkp = cfh.GetKalmanPlane(i); if(befkp->GetTrackCluster()->GetPlane()*direction <bestkp->GetTrackCluster()->GetPlane()*direction) foundkp = befkp; } Double_t forwardchi2 = 0.; Int_t befswimfail = 0; if(foundkp){ befswimfail = FitFrom(foundkp,bestkp,0,-1); forwardchi2 = bestkp->GetPreChi2(0); MSG("FitTrackSR",Msg::kDebug) << "checking forward fit, previous " << "cluster plane = " << foundkp->GetTrackCluster()->GetPlane() << " min/max strip = " << foundkp->GetTrackCluster()->GetMinStrip() << "/" << foundkp->GetTrackCluster()->GetMaxStrip() << " nswimfail = " << befswimfail << " prechi2 = " << forwardchi2 << endl; } if(!befswimfail && forwardchi2<fParmMaxLocalPreChi2){ MSG("FitTrackSR",Msg::kDebug) << "ADDING TO FIT" << endl; ++nadd; oldkp = bestkp; newkplist->AddLast(oldkp); cfh.SetCurrentKalmanPlane(oldkp); cfh.SetNSwimFail(cfh.GetNSwimFail()+nswimfail); } else delete bestkp; } else delete bestkp; bestkp = 0; nswimfail = 0; return nadd; }

Int_t AlgFitTrackSR::AddToFit (  CandFitTrackSRHandle &  cfh, KalmanPlaneSR *  , Int_t  , Bool_t  )  [private]

Definition at line 647 of file AlgFitTrackSR.cxx. References fCov, FitFrom(), KalmanPlaneSR::GetFilChi2(), KalmanPlaneSR::GetFilCovarianceMatrixValue(), CandFitTrackHandle::GetNSwimFail(), TrackClusterSR::GetPlane(), CandFitTrackSRHandle::GetPlaneList(), KalmanPlaneSR::GetPreChi2(), KalmanPlaneSR::GetTrackCluster(), KalmanPlaneSR::GetZDir(), InitializeFit(), MSG, pFor, CandFitTrackSRHandle::SetCurrentKalmanPlane(), CandFitTrackSRHandle::SetEnd(), and CandFitTrackHandle::SetNSwimFail(). { // < 0 = failure // 0 = success // this method assumes that the plane at thiskp is not already in the kalmanplanelist // note also that this method takes ownership of thiskp // if docheck is false, user is responsible for insertion into kalmanplanelist // before call to AddToFit TObjArray * KalmanPlaneList=cfh.GetPlaneList(); Bool_t found(0); Int_t ifound=0; Int_t ibef=-1; Int_t ilast; KalmanPlaneSR *kpnew = 0; KalmanPlaneSR *kp = 0; KalmanPlaneSR *kpold = 0; if(docheck){ for(Int_t i=0; i<=KalmanPlaneList->GetLast(); ++i){ // assume only one track cluster per plane in track kp = dynamic_cast<KalmanPlaneSR*>(KalmanPlaneList->At(i)); if(kp->GetTrackCluster()->GetPlane() == thiskp->GetTrackCluster()->GetPlane()) { found = 1; ifound = i; kpnew = kp; // kalmanplane already in list, need to delete thiskp delete thiskp; } else if((kp->GetZDir())*(kp->GetTrackCluster()->GetPlane())<(kp->GetZDir())*(thiskp->GetTrackCluster()->GetPlane())){ ibef = i; } }//end loop over kalman planes //didnt find a plane to bound this one, so add it at the last place in the list if(!found){ kpnew = thiskp; ilast = KalmanPlaneList->GetLast(); for(int i=ilast; i>ibef; --i){ kp = dynamic_cast<KalmanPlaneSR*>(KalmanPlaneList->At(i)); KalmanPlaneList->AddAtAndExpand(kp,i+1); } KalmanPlaneList->AddAtAndExpand(kpnew,ibef+1); }//end didnt find a plane }//end if docheck else kpnew = thiskp; cfh.SetCurrentKalmanPlane(kpnew); //if this is the first kp added to the track, initialize the fit if(KalmanPlaneList->GetLast()==0 || (docheck && ((found && ifound==0) || ibef<0))) InitializeFit(cfh,0,iterate); else{ //get the plane of the kp immediately before the one to add to the //fit kpold = dynamic_cast<KalmanPlaneSR*> (KalmanPlaneList->At(KalmanPlaneList->IndexOf(kpnew)-1)); //fit from the previously fit plane to the one to add to the fit pFor++; Int_t nswimfail = FitFrom(kpold,kpnew,0,iterate); if(nswimfail<0){ return nswimfail; } cfh.SetNSwimFail(cfh.GetNSwimFail()+nswimfail); for(int i1=0; i1<5; ++i1){ for(int i2=0; i2<5; ++i2){ fCov[i1][i2] = kpnew->GetFilCovarianceMatrixValue(0,i1,i2); } }//end loop to fill fCov covariance matrix MSG("FitTrackSR",Msg::kDebug) << "fitting forward, prechi2 = " << kpnew->GetPreChi2(0) << " filchi2 = " << kpnew->GetFilChi2(0) << endl; cfh.SetEnd(kpnew); }//end else return 0; }

Int_t AlgFitTrackSR::AddToFit (  CandFitTrackSRHandle &  cfh, TrackClusterSR *  , Int_t  )  [private]

Definition at line 628 of file AlgFitTrackSR.cxx. References TrackClusterSR::GetPlane(), CandTrackHandle::GetRange(), CandRecoHandle::GetTimeSlope(), KalmanPlaneSR::GetTrackCluster(), KalmanPlaneSR::SetRange(), KalmanPlaneSR::SetZDir(), and tc. Referenced by AddClustersToFit(), and AddForwardBestKPToFit(). { // < 0 = failure // 0 = success KalmanPlaneSR *kpnew = new KalmanPlaneSR(tc); kpnew->SetRange(cfh.GetRange(kpnew->GetTrackCluster()->GetPlane())); if(cfh.GetTimeSlope()>0. || !fParmIsCosmic) kpnew->SetZDir(1); else kpnew->SetZDir(-1); return AddToFit(cfh, kpnew, iterate,1); }

void AlgFitTrackSR::CalculateFilChi2 (  KalmanPlaneSR *  currentkp, Int_t  idir, Double_t  = 0.0119 )  [private]

Definition at line 1514 of file AlgFitTrackSR.cxx. References KalmanPlaneSR::GetFilCovarianceMatrixValue(), KalmanPlaneSR::GetFilStateValue(), TrackClusterSR::GetPlaneView(), TrackClusterSR::GetTPos(), KalmanPlaneSR::GetTrackCluster(), MSG, and KalmanPlaneSR::SetFilChi2(). Referenced by Filter(). { Int_t index = 0; //default U view, change if V view if(currentkp->GetTrackCluster()->GetPlaneView()==PlaneView::kV) index = 1; Double_t dtpos = 0.; dtpos = currentkp->GetFilStateValue(index,idir)-currentkp->GetTrackCluster()->GetTPos(); Double_t filChi2 = dtpos*dtpos/(poserr*poserr); currentkp->SetFilChi2(filChi2, idir); MSG("FitTrackSR",Msg::kDebug) << "filchi2 = " << filChi2 << " dtpos = " << dtpos << " tposerr2 = " << poserr*poserr << " cov = " << currentkp->GetFilCovarianceMatrixValue(idir,index,index) << endl; }

void AlgFitTrackSR::CalculateFilCovariance (  KalmanPlaneSR *  currentkp, Int_t  idir )  [private]

Definition at line 1539 of file AlgFitTrackSR.cxx. References KalmanPlaneSR::GetGainValue(), TrackClusterSR::GetPlaneView(), KalmanPlaneSR::GetPreCovarianceMatrixValue(), KalmanPlaneSR::GetTrackCluster(), kdUdZ, kdVdZ, kQoverP, kU, kV, and KalmanPlaneSR::SetFilCovarianceMatrixValue(). Referenced by Filter(). { //do the final step of the Kalman filter //C_{k} = (I - K_{k}H_{k})C_{k}^{k-1} Int_t index = 0; //default U view, change if V view if(currentkp->GetTrackCluster()->GetPlaneView()==PlaneView::kV) index = 1; TMatrixD preCov(5,5); for(Int_t i = 0; i<5; ++i){ for(Int_t j = 0; j<5; ++j){ preCov(i,j) = currentkp->GetPreCovarianceMatrixValue(idir,i,j); } } TMatrixD filCovariance(5,5); //set the filtered covariance matrix to the identity matrix filCovariance.UnitMatrix(); //subtract K_{k}H_{k} from the identity filCovariance(kU,index) -= currentkp->GetGainValue(kU,idir); filCovariance(kV,index) -= currentkp->GetGainValue(kV,idir); filCovariance(kdUdZ,index) -= currentkp->GetGainValue(kdUdZ,idir); filCovariance(kdVdZ,index) -= currentkp->GetGainValue(kdVdZ,idir); filCovariance(kQoverP,index) -= currentkp->GetGainValue(kQoverP,idir); //multiply the difference by the pre filtered covariance matrix filCovariance *= preCov; // diagonal elements should be positive Double_t covlim = 1.e-8; if(filCovariance(kU,kU)<covlim) filCovariance(kU,kU) = covlim; if(filCovariance(kV,kV)<covlim) filCovariance(kV,kV) = covlim; if(filCovariance(kdUdZ,kdUdZ)<covlim) filCovariance(kdUdZ,kdUdZ) = covlim; if(filCovariance(kdVdZ,kdVdZ)<covlim) filCovariance(kdVdZ,kdVdZ) = covlim; if(filCovariance(kQoverP,kQoverP)<covlim) filCovariance(kQoverP,kQoverP) = covlim; for(Int_t i = 0; i<5; ++i){ for(Int_t j = 0; j<5; ++j){ currentkp->SetFilCovarianceMatrixValue(idir,i,j, filCovariance(i,j)); } } }

void AlgFitTrackSR::CalculateFilState (  KalmanPlaneSR *  currentkp, Int_t  idir )  [private]

Definition at line 1449 of file AlgFitTrackSR.cxx. References fParmMaxAngle, fParmMaxQPFrac, KalmanPlaneSR::GetFilStateValue(), KalmanPlaneSR::GetGainValue(), TrackClusterSR::GetPlaneView(), KalmanPlaneSR::GetPreStateValue(), KalmanPlaneSR::GetRange(), TrackClusterSR::GetTPos(), KalmanPlaneSR::GetTrackCluster(), kdUdZ, kdVdZ, kQoverP, kU, kV, and KalmanPlaneSR::SetFilStateValue(). Referenced by Filter(). { Int_t index = 0; //default U view, change if V view if(currentkp->GetTrackCluster()->GetPlaneView()==PlaneView::kV) index = 1; //the following is the m_{k} - H_{k}F_{k-1}x_{k-1} //bit of the Kalman filter equations //remember that //m_{k} = measurement, //H_{k} = [1,0,0,0,0] for U view; [0,1,0,0,0] for V view //the definition of H_{k} makes the H_{k}F_{k-1}x_{k-1} product //a scalar and picks out only certain elements of F_{k-1}x_{k-1} //also, F_{k-1}x_{k-1} is found by the swimmer and put in the //pre filtered state vector of the plane to add to the fit Double_t xfac = (currentkp->GetTrackCluster()->GetTPos() - currentkp->GetPreStateValue(index,idir)); currentkp->SetFilStateValue(kU,idir, currentkp->GetPreStateValue(kU,idir) +currentkp->GetGainValue(kU,idir)*xfac); currentkp->SetFilStateValue(kV,idir, currentkp->GetPreStateValue(kV,idir) +currentkp->GetGainValue(kV,idir)*xfac); currentkp->SetFilStateValue(kdUdZ,idir, currentkp->GetPreStateValue(kdUdZ,idir) +currentkp->GetGainValue(kdUdZ,idir)*xfac); currentkp->SetFilStateValue(kdVdZ,idir, currentkp->GetPreStateValue(kdVdZ,idir) +currentkp->GetGainValue(kdVdZ,idir)*xfac); currentkp->SetFilStateValue(kQoverP,idir, currentkp->GetPreStateValue(kQoverP,idir) +currentkp->GetGainValue(kQoverP,idir)*xfac); if(TMath::Abs(currentkp->GetFilStateValue(kdUdZ,idir))>fParmMaxAngle){ currentkp->SetFilStateValue(kdUdZ,idir,(currentkp->GetFilStateValue(kdUdZ,idir)>0 ? fParmMaxAngle : -fParmMaxAngle)); } if(TMath::Abs(currentkp->GetFilStateValue(kdVdZ,idir))>fParmMaxAngle){ currentkp->SetFilStateValue(kdVdZ,idir,(currentkp->GetFilStateValue(kdVdZ,idir)>0 ? fParmMaxAngle : -fParmMaxAngle)); } Double_t maxqp = fParmMaxQP; //use a value of about 2 MeV/(g/cm^2) for the energy loss in the material //you would divide the max qp frac allowed by that value and then by the range to //get the max qp. but multiplication is a faster operation, so multiply by //500 instead. if ( fParmQPRangeCheck ) { if(currentkp->GetRange()>0. && fParmMaxQPFrac*500./currentkp->GetRange()<maxqp) { maxqp = fParmMaxQPFrac*500./currentkp->GetRange(); } } if(TMath::Abs(currentkp->GetFilStateValue(kQoverP,idir))>maxqp){ currentkp->SetFilStateValue(kQoverP,idir,(currentkp->GetFilStateValue(kQoverP,idir)>0 ? maxqp : -maxqp)); } }

void AlgFitTrackSR::CalculateGain (  KalmanPlaneSR *  currentkp, Int_t  idir, Double_t  = 0.0119 )  [private]

Definition at line 1405 of file AlgFitTrackSR.cxx. References TrackClusterSR::GetPlaneView(), KalmanPlaneSR::GetPreCovarianceMatrixValue(), KalmanPlaneSR::GetTrackCluster(), kdUdZ, kdVdZ, kQoverP, kU, kV, and KalmanPlaneSR::SetGainValue(). Referenced by Filter(). { //the following is the //K_{k} = C_{k}^{k-1}H_{k}^{T}(V_{k}+H_{k}C_{k}^{k-1}H_{k}^{T})^{-1} //step of the filtering equations //remember that //V_{k} = measurement error, //H_{k} = [1,0,0,0,0] for U view; [0,1,0,0,0] for V view //the definition of H_{k} makes the H_{k}C_{k}^{k-1}H_{k}^{T} product //a scalar and picks out only certain elements of C_{k}^{k-1} Int_t index = 0; //default U view, change if V view if(currentkp->GetTrackCluster()->GetPlaneView()==PlaneView::kV) index = 1; currentkp->SetGainValue(kU,idir, (currentkp->GetPreCovarianceMatrixValue(idir,kU,index)/ (currentkp->GetPreCovarianceMatrixValue(idir,index,index) +poserr*poserr))); currentkp->SetGainValue(kV,idir, (currentkp->GetPreCovarianceMatrixValue(idir,kV,index)/ (currentkp->GetPreCovarianceMatrixValue(idir,index,index) +poserr*poserr))); currentkp->SetGainValue(kdUdZ,idir, (currentkp->GetPreCovarianceMatrixValue(idir,kdUdZ,index)/ (currentkp->GetPreCovarianceMatrixValue(idir,index,index) +poserr*poserr))); currentkp->SetGainValue(kdVdZ,idir, (currentkp->GetPreCovarianceMatrixValue(idir,kdVdZ,index)/ (currentkp->GetPreCovarianceMatrixValue(idir,index,index) +poserr*poserr))); currentkp->SetGainValue(kQoverP,idir, (currentkp->GetPreCovarianceMatrixValue(idir,kQoverP,index)/ (currentkp->GetPreCovarianceMatrixValue(idir,index,index) +poserr*poserr))); }

void AlgFitTrackSR::CalculateNoise (  KalmanPlaneSR *  prevkp, KalmanPlaneSR *  currentkp, Int_t  idir, Int_t  befaft = 0 )  [private]

Definition at line 1079 of file AlgFitTrackSR.cxx. References fParmdedx, fParmPlnRadLen, BField::GetBField(), KalmanPlaneSR::GetFilStateValue(), UgliGeomHandle::GetNearestSteelPlnHandle(), TrackClusterSR::GetPlane(), SwimGeo::GetSteel(), KalmanPlaneSR::GetTrackCluster(), KalmanPlaneSR::GetVldContext(), SwimPlaneInterface::GetZ(), UgliSteelPlnHandle::GetZ0(), TrackClusterSR::GetZPos(), SwimGeo::Instance(), kdUdZ, kdVdZ, kInvSqrt2, kQoverP, kU, kV, MSG, and KalmanPlaneSR::SetNoiseMatrixValue(). Referenced by Predict(). { MSG("FitTrackSR", Msg::kDebug) << "In CalculateNoise" << endl; //find noise from the previous kp to the kp to add to the fit KalmanPlaneSR *thiskp = prevkp; KalmanPlaneSR *newkp = currentkp; Int_t thisdir = idir; if(befaft){ thiskp = currentkp; newkp = prevkp; thisdir = 1-thisdir; } Int_t dplane = (currentkp->GetTrackCluster()->GetPlane() - prevkp->GetTrackCluster()->GetPlane()); Int_t izdir = (dplane>0 ? 0 : 1); // izdir gives the direction of the swim (0 = forward, 1 = backward) // relative to z // this is different from thisdir/idir which gives direction of swim // relative to time // idir=0 means forward in time, idir=1 means backward in time // a particle could be going forward in z and forward in time or // any of the 4 combinations (for cosmics and atmospherics) Double_t dudz = thiskp->GetFilStateValue(kdUdZ,idir); Double_t dvdz = thiskp->GetFilStateValue(kdVdZ,idir); Double_t dsdz = sqrt(1.+dudz*dudz+dvdz*dvdz); Double_t qp = thiskp->GetFilStateValue(kQoverP,idir); // when qp is too small, errors go to zero, so place lower bound // 0.01 (= 100 GeV/c) is reasonable if(TMath::Abs(qp)<0.01) qp = (qp>0 ? 0.01 : -0.01); Double_t locradlen = TMath::Abs(dsdz*(Double_t)(dplane)*fParmPlnRadLen); // we should really use an average of 1/p in calculating sigmams2 Double_t sigmams2 = (0.0136*TMath::Abs(qp)*TMath::Sqrt(locradlen) *(1.+0.038*TMath::Log(locradlen))); sigmams2 *= sigmams2; Double_t p3p3 = 0.5*sigmams2*(1.+dudz*dudz)*dsdz*dsdz; Double_t p3p4 = 0.5*sigmams2*dudz*dvdz*dsdz*dsdz; Double_t p4p4 = 0.5*sigmams2*(1.+dvdz*dvdz)*dsdz*dsdz; // SwimPlaneInterfaceListSR *spil = // SwimPlaneInterfaceListSR::Instance(const_cast<VldContext*>(GetVldContext())); SwimGeo *spil = SwimGeo::Instance(*(const_cast<VldContext*>(currentkp->GetVldContext()))); Double_t z1 = spil->GetSteel(prevkp->GetTrackCluster()->GetZPos(),izdir)->GetZ(); Double_t z2 = spil->GetSteel(z1,izdir)->GetZ(); // treat iron as discrete scatter; in future can make this more rigorous Double_t dzscatter = TMath::Abs(currentkp->GetTrackCluster()->GetZPos()-0.5*(z1+z2)); Double_t dzscatter2 = dzscatter*dzscatter; Double_t eloss = fParmdedx*TMath::Abs((Double_t)(dplane)); VldContext vldc = *(prevkp->GetVldContext()); UgliGeomHandle ugh(vldc); BField bf(vldc); TVector3 uvz; uvz(0) = thiskp->GetFilStateValue(kU,idir); uvz(1) = thiskp->GetFilStateValue(kV,idir); // rwh: KalmanPlane's TrackCluster z is not necessarily inside the steel //uvz(2) = thiskp->GetTrackCluster()->GetZPos(); // translate an semi-arbitrary z to z0 of the nearest steel plane uvz(2) = ugh.GetNearestSteelPlnHandle(thiskp->GetTrackCluster()->GetZPos()).GetZ0(); TVector3 xyz(kInvSqrt2*(uvz[0]-uvz[1]),kInvSqrt2*(uvz[0]+uvz[1]),uvz[2]); TVector3 bxyz = bf.GetBField(xyz); // rotate bxyz into buvz // bxyz in Tesla, convert into GeV/c/m // warning: numbers are hard coded here TVector3 buvz(kInvSqrt2*(bxyz[1]+bxyz[0]),kInvSqrt2*(bxyz[1]-bxyz[0]),bxyz[2]); buvz[0] *= 0.3; buvz[1] *= 0.3; buvz[2] *= 0.3; // we assume the uncertainty in energy loss to follow a Landau distribution, // and use the FWHM from Ahlen (RMP, Vol. 52, p. 143, 1980) Double_t sigmaeloss2 = 0.25*eloss*dsdz*qp*qp; sigmaeloss2 *= sigmaeloss2; currentkp->SetNoiseMatrixValue(idir,kU,kU,dzscatter2*p3p3); currentkp->SetNoiseMatrixValue(idir,kU,kV,dzscatter2*p3p4); currentkp->SetNoiseMatrixValue(idir,kU,kdUdZ,dzscatter*p3p3); currentkp->SetNoiseMatrixValue(idir,kU,kdVdZ,dzscatter*p3p4); currentkp->SetNoiseMatrixValue(idir,kU,kQoverP, (dzscatter*sigmaeloss2*TMath::Abs((Double_t)(dplane)) *buvz(1)*0.0254*Munits::m*dsdz*(1.+dudz*dudz))); currentkp->SetNoiseMatrixValue(idir,kV,kU,dzscatter2*p3p4); currentkp->SetNoiseMatrixValue(idir,kV,kV,dzscatter2*p4p4); currentkp->SetNoiseMatrixValue(idir,kV,kdUdZ,dzscatter*p3p4); currentkp->SetNoiseMatrixValue(idir,kV,kdVdZ,dzscatter*p4p4); currentkp->SetNoiseMatrixValue(idir,kV,4, (dzscatter*sigmaeloss2*TMath::Abs((Double_t)(dplane)) *buvz(0)*0.0254*Munits::m*dsdz*(1.+dvdz*dvdz))); currentkp->SetNoiseMatrixValue(idir,kdUdZ,kU,dzscatter*p3p3); currentkp->SetNoiseMatrixValue(idir,kdUdZ,kV,dzscatter*p3p4); currentkp->SetNoiseMatrixValue(idir,kdUdZ,kdUdZ,p3p3); currentkp->SetNoiseMatrixValue(idir,kdUdZ,kdVdZ,p3p4); currentkp->SetNoiseMatrixValue(idir,kdUdZ,kQoverP, (sigmaeloss2*TMath::Abs((Double_t)(dplane))*buvz(1) *0.0254*Munits::m*dsdz*(1.+dudz*dudz))); currentkp->SetNoiseMatrixValue(idir,kdVdZ,kU,dzscatter*p3p4); currentkp->SetNoiseMatrixValue(idir,kdVdZ,kV,dzscatter*p4p4); currentkp->SetNoiseMatrixValue(idir,kdVdZ,kdUdZ,p3p4); currentkp->SetNoiseMatrixValue(idir,kdVdZ,kdVdZ,p4p4); currentkp->SetNoiseMatrixValue(idir,kdVdZ,kQoverP, (sigmaeloss2*TMath::Abs((Double_t)(dplane))*buvz(0) *0.0254*Munits::m*dsdz*(1.+dvdz*dvdz))); currentkp->SetNoiseMatrixValue(idir,kQoverP,kU, (dzscatter*sigmaeloss2*TMath::Abs((Double_t)(dplane)) *buvz(1)*0.0254*Munits::m*dsdz*(1.+dudz*dudz))); currentkp->SetNoiseMatrixValue(idir,kQoverP,kV, (dzscatter*sigmaeloss2*TMath::Abs((Double_t)(dplane)) *buvz(0)*0.0254*Munits::m*dsdz*(1.+dvdz*dvdz))); currentkp->SetNoiseMatrixValue(idir,kQoverP,kdUdZ, (sigmaeloss2*TMath::Abs((Double_t)(dplane))*buvz(1) *0.0254*Munits::m*dsdz*(1.+dudz*dudz))); currentkp->SetNoiseMatrixValue(idir,kQoverP,kdVdZ, (sigmaeloss2*TMath::Abs((Double_t)(dplane))*buvz(0) *0.0254*Munits::m*dsdz*(1.+dvdz*dvdz))); currentkp->SetNoiseMatrixValue(idir,kQoverP,kQoverP,sigmaeloss2); MSG("FitTrackSR",Msg::kDebug) << "du/dz, dv/dz, ds/dz, qp = " << dudz << " " << dvdz << " " << dsdz << " " << qp << endl << "dplane, dzscatter, z1, z2, radlen = " << dplane << " " << dzscatter << " " << z1 << " " << z2 << " " << locradlen << endl << "sigmams2 = " << sigmams2 << endl << "p3p3, p3p4, p4p4 = " << p3p3 << " " << p3p4 << " " << p4p4 << endl; }

void AlgFitTrackSR::CalculatePreChi2 (  KalmanPlaneSR *  currentkp, Int_t  idir, Double_t  = 0.0119 )  [private]

Definition at line 1380 of file AlgFitTrackSR.cxx. References TrackClusterSR::GetPlaneView(), KalmanPlaneSR::GetPreCovarianceMatrixValue(), KalmanPlaneSR::GetPreStateValue(), TrackClusterSR::GetTPos(), KalmanPlaneSR::GetTrackCluster(), MSG, and KalmanPlaneSR::SetPreChi2(). Referenced by Predict(). { Int_t index = 0; //default view is U, change if this is a V view plane if(currentkp->GetTrackCluster()->GetPlaneView()==PlaneView::kV) index = 1; Double_t dtpos = 0.; dtpos = currentkp->GetPreStateValue(index,idir)-currentkp->GetTrackCluster()->GetTPos(); Double_t preChi2 = dtpos*dtpos/(poserr*poserr); currentkp->SetPreChi2(preChi2, idir); MSG("FitTrackSR",Msg::kDebug) << "prechi2 = " << preChi2 << " dtpos = " << dtpos << " tposerr2 = " << poserr*poserr << " cov = " << currentkp->GetPreCovarianceMatrixValue(idir,index,index) << endl; }

void AlgFitTrackSR::CalculatePreCovariance (  KalmanPlaneSR *  prevkp, KalmanPlaneSR *  currentkp, Int_t  idir )  [private]

Definition at line 1223 of file AlgFitTrackSR.cxx. References KalmanPlaneSR::GetFilCovarianceMatrixValue(), KalmanPlaneSR::GetNoiseMatrixValue(), KalmanPlaneSR::GetPreCovarianceMatrixValue(), KalmanPlaneSR::GetPropagatorMatrixValue(), kdUdZ, kdVdZ, kQoverP, kU, kV, MSG, and KalmanPlaneSR::SetPreCovarianceMatrixValue(). Referenced by Predict(). { MSG("FitTrackSR", Msg::kDebug) << "In CalculatePreCovariance" << endl; //get the propagator and filtered covariance matricies from the //previously fit kp TMatrixD propagator(5,5); TMatrixD filCovariance(5,5); for(Int_t i = 0; i<5; ++i){ for(Int_t j = 0; j<5; ++j){ propagator(i,j) = prevkp->GetPropagatorMatrixValue(idir,i,j); filCovariance(i,j) = prevkp->GetFilCovarianceMatrixValue(idir,i,j); //start the calculation of the pre-filtered covariance matrix //for the kp to add to the fit. set it = noise matrix of previously //fit plane currentkp->SetPreCovarianceMatrixValue(idir,i,j,prevkp->GetNoiseMatrixValue(idir,i,j)); } } //remember that propagator is really the transpose of the propagator matrix //so in the following cov0 = F_{k-1}C_{k-1} TMatrixD cov0(propagator,TMatrixD::kTransposeMult,filCovariance); //the following makes cov0 = F_{k-1}C_{k-1}F^{T}_{k-1} cov0 *= propagator; //finish finding the pre-filtered covariance matrix - //C_{k}^{k-1} = F_{k-1}C_{k-1}F^{T}_{k-1}+Q_{k-1} Double_t preCovValue = 0.; for(Int_t i = 0; i<5; ++i){ for(Int_t j = 0; j<5; ++j){ preCovValue = currentkp->GetPreCovarianceMatrixValue(idir,i,j)+cov0(i,j); currentkp->SetPreCovarianceMatrixValue(idir,i,j,preCovValue); } } // diagonal elements should be positive Double_t covlim = 1.e-8; if(currentkp->GetPreCovarianceMatrixValue(idir,kU,kU)<covlim) currentkp->SetPreCovarianceMatrixValue(idir,kU,kU,covlim); if(currentkp->GetPreCovarianceMatrixValue(idir,kV,kV)<covlim) currentkp->SetPreCovarianceMatrixValue(idir,kV,kV,covlim); if(currentkp->GetPreCovarianceMatrixValue(idir,kdUdZ,kdUdZ)<covlim) currentkp->SetPreCovarianceMatrixValue(idir,kdUdZ,kdUdZ,covlim); if(currentkp->GetPreCovarianceMatrixValue(idir,kdVdZ,kdVdZ)<covlim) currentkp->SetPreCovarianceMatrixValue(idir,kdVdZ,kdVdZ,covlim); if(currentkp->GetPreCovarianceMatrixValue(idir,kQoverP,kQoverP)<covlim) currentkp->SetPreCovarianceMatrixValue(idir,kQoverP,kQoverP,covlim); }

Int_t AlgFitTrackSR::CalculatePreState (  KalmanPlaneSR *  prevkp, KalmanPlaneSR *  currentkp, Int_t  idir )  [private]

Definition at line 1278 of file AlgFitTrackSR.cxx. References fParmMaxAngle, fParmMaxQPFrac, fSwimmer, KalmanPlaneSR::GetFilStateValue(), KalmanPlaneSR::GetPreStateValue(), KalmanPlaneSR::GetRange(), KalmanPlaneSR::GetTrackCluster(), KalmanPlaneSR::GetVldContext(), TrackClusterSR::GetZPos(), kdUdZ, kdVdZ, kQoverP, kU, kV, MSG, KalmanPlaneSR::SetPreStateValue(), and Swim(). Referenced by Predict(). { //Swim from the previously fit plane to the plane to add to the fit //to get the estimate of the state vector at the new plane - this //takes care of the F_{k-1}x_{k-1} matrix multiplication in the //Kalman filter equations Int_t nswimfail=0; if(fSwimmer!=2) MSG("FitTrackSR",Msg::kError) << "fSwimmer!=2 currently unsupported" << endl; MSG("FitTrackSR", Msg::kDebug) << "current Filtered State = [" << prevkp->GetFilStateValue(kU,idir) << "," << prevkp->GetFilStateValue(kV,idir) << "," << prevkp->GetFilStateValue(kdUdZ,idir) << "," << prevkp->GetFilStateValue(kdVdZ,idir) << "," << prevkp->GetFilStateValue(kQoverP,idir) << "]" << endl; if(fSwimmer==2){ // Using Swimmer package Double_t swimstate[5]; Double_t state[] = {0.,0.,0.,0.,0.}; swimstate[kU] = prevkp->GetFilStateValue(kU,idir); swimstate[kV] = prevkp->GetFilStateValue(kV,idir); swimstate[kdUdZ] = prevkp->GetFilStateValue(kdUdZ,idir); swimstate[kdVdZ] = prevkp->GetFilStateValue(kdVdZ,idir); swimstate[kQoverP] = prevkp->GetFilStateValue(kQoverP,idir); Bool_t swimstatus = false; Int_t nswim=0; while(!swimstatus && nswim<100){ swimstatus = Swim(swimstate, state, prevkp->GetTrackCluster()->GetZPos(), currentkp->GetTrackCluster()->GetZPos(),idir, currentkp->GetVldContext()); // if outside fiducial detector, return failure flag if(TMath::Abs(state[kU])>5. || TMath::Abs(state[kV])>5.){ MSG("FitTrackSR",Msg::kDebug) << "outside detector, u,v = " << state[kU] << " " << state[kV] << " fail" << endl; return -1; } if(!swimstatus){ // swim failed, double momentum swimstate[kQoverP] *= 0.5; MSG("FitTrackSR",Msg::kDebug) << "swim failed, doubling momentum to q/p = " << swimstate[kQoverP] << endl; ++nswimfail; } ++nswim; }//end loop to swim. currentkp->SetPreStateValue(kU,idir,state[kU]); currentkp->SetPreStateValue(kV,idir,state[kV]); currentkp->SetPreStateValue(kdUdZ,idir,state[kdUdZ]); currentkp->SetPreStateValue(kdVdZ,idir,state[kdVdZ]); currentkp->SetPreStateValue(kQoverP,idir,state[kQoverP]); }//end if using Swim package if(TMath::Abs(currentkp->GetPreStateValue(kdUdZ,idir))>fParmMaxAngle){ currentkp->SetPreStateValue(kdUdZ,idir, (currentkp->GetPreStateValue(kdUdZ,idir)>0 ? fParmMaxAngle : -fParmMaxAngle)); } if(TMath::Abs(currentkp->GetPreStateValue(kdVdZ,idir))>fParmMaxAngle){ currentkp->SetPreStateValue(kdVdZ,idir, (currentkp->GetPreStateValue(kdVdZ,idir)>0 ? fParmMaxAngle : -fParmMaxAngle)); } Double_t maxqp = fParmMaxQP; //use a value of about 2 MeV/(g/cm^2) for the energy loss in the material //you would divide the max qp frac allowed by that value and then by the range to //get the max qp. but multiplication is a faster operation, so multiply by //500 instead. if ( fParmQPRangeCheck ) { if(currentkp->GetRange()>0. && fParmMaxQPFrac*500./currentkp->GetRange()<maxqp) { maxqp = fParmMaxQPFrac*500./currentkp->GetRange(); } } if(TMath::Abs(currentkp->GetPreStateValue(kQoverP,idir))>maxqp){ currentkp->SetPreStateValue(kQoverP,idir, (currentkp->GetPreStateValue(kQoverP,idir)>0 ? maxqp : -maxqp)); } return nswimfail; }

Int_t AlgFitTrackSR::CalculatePropagator (  KalmanPlaneSR *  prevkp, KalmanPlaneSR *  currentkp, Int_t  idir )  [private]

Definition at line 953 of file AlgFitTrackSR.cxx. References fParmDState, fSwimmer, KalmanPlaneSR::GetFilStateValue(), KalmanPlaneSR::GetTrackCluster(), KalmanPlaneSR::GetVldContext(), TrackClusterSR::GetZPos(), kQoverP, kU, kV, MSG, KalmanPlaneSR::SetPropagatorMatrixValue(), and Swim(). Referenced by Predict(). { // we are actually calculating the tranpose MSG("FitTrackSR", Msg::kDebug) << "In CalculatePropagator" << endl; Int_t nswimfail=0; Double_t dstate[5]; Double_t thisFilState[5]; for(Int_t i=0; i<5; ++i){ //set the filtered state based on the previously fit plane thisFilState[i] = prevkp->GetFilStateValue(i,idir); dstate[i] = fParmDState[i]; if(i==4){ dstate[i] = (fParmDState[i]*thisFilState[kQoverP]*thisFilState[kQoverP] /(1.-fParmDState[i]*TMath::Abs(thisFilState[kQoverP]))); dstate[i] = fParmDState[i]*TMath::Abs(thisFilState[kQoverP]); if(dstate[i]<.01) dstate[i] = .01; // corresponds to constant fParmDState[kQoverP] energy interval } // MSG("FitTrackSR", Msg::kDebug) << "thisFilState[" << i << "] = " // << thisFilState[i] << " dstate[" << i // << "] = " << dstate[i] << endl; } Bool_t swimstatus = false; Int_t nswim=0; Double_t swimstate[5]; Double_t pstate[] = {0.,0.,0.,0.,0.}; Double_t nstate[] = {0.,0.,0.,0.,0.}; //loop over the vectors and swim the particle until you have a good //swim or you excede the maximum iterations while(!swimstatus && nswim<100){ swimstatus = true; //loop over all the states and change them by +/- dstate and swim from //that point to the next z pos. this figures out the propagator matrix for(Int_t i=0; i<5; ++i){ if(i==4){ dstate[i] = (fParmDState[i]*thisFilState[kQoverP]*thisFilState[kQoverP] /(1.-fParmDState[i]*TMath::Abs(thisFilState[kQoverP]))); dstate[i] = fParmDState[i]*TMath::Abs(thisFilState[kQoverP]); if(dstate[i]<.01) dstate[i] = .01; } swimstate[kU] = thisFilState[kU]; swimstate[kV] = thisFilState[kV]; swimstate[kdUdZ] = thisFilState[kdUdZ]; swimstate[kdVdZ] = thisFilState[kdVdZ]; swimstate[kQoverP] = thisFilState[kQoverP]; swimstate[i] += dstate[i]; if(fSwimmer!=2) MSG("FitTrackSR",Msg::kError) << "fSwimmer!=2 currently unsupported" << endl; else if(fSwimmer==2){ // Using Swimmer package if(!(Swim(swimstate, pstate, prevkp->GetTrackCluster()->GetZPos(), currentkp->GetTrackCluster()->GetZPos(), idir, prevkp->GetVldContext()))) swimstatus = false; // if outside fiducial detector, return failure flag if(TMath::Abs(pstate[kU])>10. || TMath::Abs(pstate[kV])>10.){ return -1; } swimstate[kU] = thisFilState[kU]; swimstate[kV] = thisFilState[kV]; swimstate[kdUdZ] = thisFilState[kdUdZ]; swimstate[kdVdZ] = thisFilState[kdVdZ]; swimstate[kQoverP] = thisFilState[kQoverP]; swimstate[i] -= dstate[i]; if(!(Swim(swimstate, nstate, prevkp->GetTrackCluster()->GetZPos(), currentkp->GetTrackCluster()->GetZPos(), idir, currentkp->GetVldContext()))) swimstatus = false; // if outside fiducial detector, return failure flag if(TMath::Abs(nstate[kU])>10. || TMath::Abs(nstate[kV])>10.){ MSG("FitTrackSR",Msg::kDebug) << "outside detector, pstate u,v = " << pstate[kU] << " " << pstate[kV] << " fail" << endl; return -1; } //get the value of the propagator matrix terms for the current column of the //matrix Double_t propValue = 0.; for(Int_t j=0; j<5; ++j){ if(dstate[i]!=0.) propValue = (pstate[j]-nstate[j])/(2.*dstate[i]); else{ if(j==i) propValue = 1.; else propValue = 0.; } prevkp->SetPropagatorMatrixValue(idir, i, j, propValue); } }//end if fSwimmer is for the Swim package }//end loop over matrix rows if(!swimstatus){ // swim failed, double momentum thisFilState[kQoverP] *= 0.5; ++nswimfail; MSG("FitTrackSR",Msg::kDebug) << "failed swim, doubling momentum to q/p = " << thisFilState[kQoverP] << " nswim = " << nswim << endl; } ++nswim; } return nswimfail; }

Int_t AlgFitTrackSR::DoKalmanFit (  TObjArray &  planeClusterList, CandFitTrackSRHandle &  cfh, Int_t &  istatus, Int_t  direction, Int_t  dosearch )  [private]

Definition at line 2200 of file AlgFitTrackSR.cxx. References AddClustersToFit(), CandFitTrackSRHandle::AddKalmanPlaneAt(), CandFitTrackSRHandle::ClearKalmanPlaneList(), FindDownstreamPlanes(), FindUpstreamPlanes(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetPlane(), TrackClusterSR::GetPlaneView(), KalmanPlaneSR::GetTrackCluster(), KalmanPlaneSR::GetZDir(), InitializeFit(), IterateKalmanFit(), MSG, RemoveBadPointsFromFit(), ResetTrackClusterList(), TrackClusterSR::SetIsValid(), CandFitTrackSRHandle::SetVtx(), and tc. Referenced by RunAlg(). { Int_t nadd=1; Int_t iadd=0; istatus=0; Int_t nbadfit = 1; Int_t ibadfit = 0; Int_t nu = 0; Int_t nv = 0; TrackClusterSR *tc = 0; //make maps of fit planes in each view. map the plane number //to an int - 0 for a bad fit and 1 for a good fit map<Int_t, Int_t> uFitPlanes; map<Int_t, Int_t> vFitPlanes; //loop over the clusters in the slice to fill the uFitPlanes and vFitPlanes //maps. initialize all plane values to zero, and reset them in the for //loop below if the fit was good. this way you can see how many hit planes //are being skipped by dropping poorly fit planes - otherwise you might //just keep removing planes but forgetting that you had the ones you dropped TIter planeItr(&planeClusterList); TObjArray * plnarray; while( (plnarray = (TObjArray * )planeItr.Next()) ){ TIter clusterItr(plnarray); while( (tc = (TrackClusterSR *)clusterItr.Next()) ){ if(tc->GetPlaneView() == PlaneView::kU) uFitPlanes[tc->GetPlane()] = 0; else if(tc->GetPlaneView() == PlaneView::kV) vFitPlanes[tc->GetPlane()] = 0; } } //loop as long as the fit flag is 0, you have added some points to the //track, and you havent gone over the max number of iterations while(!istatus && nadd && iadd<=fParmMaxIterate2){ nbadfit=1; ibadfit=0; //loop as long as the fit flag is 0 and you have points to remove from the //track while(!istatus && nbadfit){ //reset the number of bad points in the fit to zero for this loop nbadfit = 0; //iterate the Kalman fit forward and backwards until the fit q/p //converges istatus = IterateKalmanFit(planeClusterList, cfh, nu, nv, direction, dosearch); //if all systems are still go, look for and remove any bad points in the //track if (!istatus){ MSG("FitTrackSR",Msg::kDebug) << "Removing iteration " << iadd << "/" << ibadfit << endl; nbadfit = RemoveBadPointsFromFit(cfh, planeClusterList, uFitPlanes, vFitPlanes,nu, nv, direction); }//end if tracking worked ++ibadfit; }//end loop while there are still bad points or tracking //now look for planes to add to the fit //reset the counter for the number of added points to the track nadd = 0; //make sure the fit is still ok if(!istatus && dosearch){ cfh.ClearKalmanPlaneList(); //add the remaining clusters after removal of bad points to the fit. //use 999 as the value for the third parameter //to ensure that we just scale the previous covariance matrix nu = 0; nv = 0; AddClustersToFit(planeClusterList, cfh, 999, nu, nv, direction); //make sure you have enough planes in both view to procede if(nu>=3 && nv>=3){ // find additional planes downstream nadd += FindDownstreamPlanes(planeClusterList, cfh, direction); //reinitialize the fit and use // 999 to ensure we simply scale previous covariance matrix InitializeFit(cfh,1,999); //make an array to hold all the new kp's you find TObjArray *newkplist = new TObjArray(1,0); nadd += FindUpstreamPlanes(planeClusterList, cfh, newkplist, direction); //reset the track clusters in the clusterItr for the next iteration ResetTrackClusterList(planeClusterList); MSG("FitTrackSR",Msg::kDebug) << "Last iteration " << iadd << "/" << ibadfit << endl; //clear the cfh Kalman Plane list - those objects are all in the //newkplist and will get deleted once the clusters have all their //IsValid flags reset cfh.ClearKalmanPlaneList(0); const KalmanPlaneSR *currentkp = 0; Int_t sizeoflist = newkplist->GetLast(); //loop over the newkplist in reverse to mark the clusters used in the fit KalmanPlaneSR *kp = 0; for(Int_t i=sizeoflist; i>=0; --i){ kp = dynamic_cast<KalmanPlaneSR*>(newkplist->At(i)); assert(kp); if(!currentkp) currentkp = cfh.AddKalmanPlaneAt(kp,sizeoflist-i); else cfh.AddKalmanPlaneAt(kp,sizeoflist-i); TObjArray * planeClusters=(TObjArray *)planeClusterList.At(kp->GetTrackCluster()->GetPlane()); TIter clusterItr(planeClusters); while((tc = (TrackClusterSR *)(clusterItr.Next()) )){ if(tc->GetPlane()==kp->GetTrackCluster()->GetPlane()){ if(tc->GetMinStrip() == kp->GetTrackCluster()->GetMinStrip() && tc->GetMaxStrip() == kp->GetTrackCluster()->GetMaxStrip() ) tc->SetIsValid(true); } } MSG("FitTrackSR",Msg::kDebug) << " TC " << kp->GetTrackCluster()->GetPlane() << " " << kp->GetTrackCluster()->GetMinStrip() << " " << kp->GetTrackCluster()->GetMaxStrip() << " " << kp->GetZDir() << endl; delete kp; }//end loop over new kplist MSG("FitTrackSR", Msg::kDebug) << "vertex direction = " << currentkp->GetZDir() << endl; //set the vertex to the first kp cfh.SetVtx(currentkp); newkplist->Clear(); delete newkplist; }//end if enough planes hit to start looking for more to add }//end if fit is still ok after iteration ste ++iadd; }//end while !iadd //the iadd-1 is so that we can just put the returned value right into //cfh.SetNIterate() return iadd-1; }

Int_t AlgFitTrackSR::Filter (  KalmanPlaneSR *  currentkp, Int_t  idir )  [private]

Definition at line 840 of file AlgFitTrackSR.cxx. References CalculateFilChi2(), CalculateFilCovariance(), CalculateFilState(), CalculateGain(), fParmMaxAngle, fParmMaxQPFrac, KalmanPlaneSR::GetFilStateValue(), KalmanPlaneSR::GetRange(), TrackClusterSR::GetTPosError(), KalmanPlaneSR::GetTrackCluster(), kdUdZ, kdVdZ, kQoverP, KalmanPlaneSR::SetFilStateValue(), and KalmanPlaneSR::SetPredictPlane(). Referenced by FitFrom(). { // // return value = 0 success // //get the K_{k} matrix CalculateGain(currentkp, idir,currentkp->GetTrackCluster()->GetTPosError()); //do the C_{k} = (I-K_{k}H_{k})C_{k}^{k-1} step CalculateFilCovariance(currentkp, idir); //find the filtered state CalculateFilState(currentkp, idir); // calculate noise and precovariance matrices onward again, values may have changed CalculateFilChi2(currentkp, idir,currentkp->GetTrackCluster()->GetTPosError()); //reset the predicted from plane -- filtered state has changed currentkp->SetPredictPlane(-999, idir); if(TMath::Abs(currentkp->GetFilStateValue(kdUdZ,idir))>fParmMaxAngle){ if(currentkp->GetFilStateValue(kdUdZ,idir)>0.) currentkp->SetFilStateValue(kdUdZ, idir, fParmMaxAngle); if(currentkp->GetFilStateValue(kdUdZ,idir)<0.) currentkp->SetFilStateValue(kdUdZ, idir, -fParmMaxAngle);; } if(TMath::Abs(currentkp->GetFilStateValue(kdVdZ, idir))>fParmMaxAngle){ if(currentkp->GetFilStateValue(kdVdZ, idir)>0.) currentkp->SetFilStateValue(kdVdZ, idir, fParmMaxAngle); if(currentkp->GetFilStateValue(kdVdZ, idir)<0.) currentkp->SetFilStateValue(kdVdZ, idir, -fParmMaxAngle); } Double_t maxqp = fParmMaxQP; //use a value of about 2 MeV/(g/cm^2) for the energy loss in the material //you would divide the max qp frac allowed by that value and then by the range to //get the max qp. but multiplication is a faster operation, so multiply by //500 instead. if ( fParmQPRangeCheck ) { if(currentkp->GetRange()>0. && fParmMaxQPFrac*500./currentkp->GetRange()<maxqp) { maxqp = fParmMaxQPFrac*500./currentkp->GetRange(); } } if(idir==0 && TMath::Abs(currentkp->GetFilStateValue(kQoverP,idir))>maxqp){ if(currentkp->GetFilStateValue(kQoverP,idir)>0.) currentkp->SetFilStateValue(kQoverP, idir, maxqp); else currentkp->SetFilStateValue(kQoverP, idir, -maxqp); } return 0; }

Int_t AlgFitTrackSR::FindDownstreamPlanes (  TObjArray &  planeClusterList, CandFitTrackSRHandle &  cfh, Int_t  direction )  [private]

Definition at line 3040 of file AlgFitTrackSR.cxx. References AddForwardBestKPToFit(), FindNumSkippedPlanes(), FitFrom(), CandFitTrackSRHandle::GetKalmanLast(), CandFitTrackSRHandle::GetKalmanPlane(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetPlane(), KalmanPlaneSR::GetPreChi2(), CandTrackHandle::GetRange(), KalmanPlaneSR::GetTrackCluster(), TrackClusterSR::IsValid(), MSG, KalmanPlaneSR::SetRange(), and KalmanPlaneSR::SetZDir(). Referenced by DoKalmanFit(). { MSG("FitTrackSR", Msg::kDebug) << "in FindDownstreamPlanes " << direction << endl; TrackClusterSR *thistc = 0; TrackClusterSR *oldtc = 0; KalmanPlaneSR *thiskp = 0; KalmanPlaneSR *bestkp = 0; KalmanPlaneSR *oldkp = cfh.GetKalmanPlane(cfh.GetKalmanLast()); Int_t nswimfail = 0; Int_t nplaneskip = 0; Bool_t isvalid = true; Int_t istatus2 = 0; Int_t nadd = 0; Int_t currentPlane = -1; Int_t oldKPPlane = oldkp->GetTrackCluster()->GetPlane(); //loop over the clusters in the cluster list //save some time - slice the cluster iterator from the last plane //in the fit to the last one in the detector in the appropriate direction //neither detector has more than 500 planes, and both start with plane 0 being //the front steel (uninstrumented) plane Bool_t dir = kIterForward; if(direction!=1) dir=kIterBackward; TIter planeItr(&planeClusterList,dir); TObjArray * plnarray; while( (plnarray = (TObjArray * )planeItr.Next()) ){ TIter clusterItr(plnarray); while( (thistc = (TrackClusterSR * )clusterItr.Next()) && isvalid ){ //make sure you are looking at a valid track cluster, if not go on to the next one if( !thistc->IsValid() ) continue; currentPlane = thistc->GetPlane(); MSG("FitTrackSR", Msg::kDebug) << "looking at plane " << currentPlane << endl; //make sure you are looking at a downstream plane if(direction*currentPlane>direction*oldKPPlane){ //see if the current track cluster is on a new plane and that //you have a potential kp to add to the fit if(bestkp && oldtc && oldtc->GetPlane() != currentPlane ){ MSG("FitTrackSR", Msg::kDebug ) << "add best kp to fit" << endl; if(AddForwardBestKPToFit(bestkp, cfh, nswimfail)){ ++nadd; oldkp = bestkp; } bestkp = 0; } oldKPPlane = oldkp->GetTrackCluster()->GetPlane(); MSG("FitTrackSR", Msg::kDebug) << " prev plane = " << oldKPPlane << " direction " << direction << endl; //on a new cluster, see if it is a good one MSG("FitTrackSR",Msg::kDebug) << "considering tc " << thistc->GetPlane() << " " << thistc->GetMinStrip() << "-" << thistc->GetMaxStrip() << " for addition to end" << endl; nplaneskip = 0; //check if the current plane is far from the last kp - need to //see how many active planes were skipped if it is if(TMath::Abs(currentPlane-oldKPPlane)>=2*fParmNSkipActive) nplaneskip = FindNumSkippedPlanes(currentPlane,oldKPPlane,oldkp, direction); //make sure you didnt skipp too many active planes if(nplaneskip<2*fParmNSkipActive){ //make a new kp for this cluster to test it out thiskp = new KalmanPlaneSR(thistc); thiskp->SetRange(cfh.GetRange(thiskp->GetTrackCluster()->GetPlane())); thiskp->SetZDir(direction); istatus2 = FitFrom(oldkp,thiskp,0,-1); MSG("FitTrackSR",Msg::kDebug) << " chi2 = " << thiskp->GetPreChi2(0) << " direction " << direction << " nswimfail = " << istatus2 << endl; //look to see if the kp you just made is better than the previous //best kp if(!istatus2 && (!bestkp || thiskp->GetPreChi2(0)<bestkp->GetPreChi2(0))){ MSG("FitTrackSR",Msg::kDebug) << "best chi2, keeping" << endl; if (bestkp) delete bestkp; bestkp = thiskp; nswimfail = istatus2; } else delete thiskp; if(istatus2<0) isvalid = 0; } else { // any planes beyond this one will not satisfy NSkipView, stop search isvalid = 0; } MSG("FitTrackSR",Msg::kDebug) << "set oldtc" << endl; oldtc = thistc; }//end if current plane is down stream of the last fit plane }//end loop over cluster list } //may still have a best kp, so make sure to add it to the fit if(bestkp){ MSG("FitTrackSR", Msg::kDebug) << "add bestkp outside of loop over tc" << endl; nadd += AddForwardBestKPToFit(bestkp, cfh, nswimfail); //delete bestkp so as to not have a memory leek bestkp = 0; } return nadd; }

Int_t AlgFitTrackSR::FindNumSkippedPlanes (  Int_t  currentPlane, Int_t  prevPlane, KalmanPlaneSR *  oldkp, Int_t  direction )  [private]

Definition at line 1905 of file AlgFitTrackSR.cxx. References fDetector, fVldContext, PlexPlaneId::GetAdjoinScint(), KalmanPlaneSR::GetFilStateValue(), PlexPlaneId::GetPlaneCoverage(), UgliGeomHandle::GetScintPlnHandle(), KalmanPlaneSR::GetTrackCluster(), UgliPlnHandle::GetZ0(), TrackClusterSR::GetZPos(), UgliScintPlnHandle::IsValid(), PlexPlaneId::IsValid(), kdUdZ, kdVdZ, kQoverP, kU, kV, and Swim(). Referenced by FindDownstreamPlanes(), and FindUpstreamPlanes(). { UgliGeomHandle ugh(fVldContext); // need to check how many active planes were actually skipped Int_t nplaneskip=0; PlexPlaneId plnid(fDetector,prevPlane,false); PlexPlaneId plnidend(fDetector,currentPlane,false); plnid = plnid.GetAdjoinScint(direction); // setup kalmanplane to do swimming Double_t currentstate[5] = {oldkp->GetFilStateValue(kU,0), oldkp->GetFilStateValue(kV,0), oldkp->GetFilStateValue(kdUdZ,0), oldkp->GetFilStateValue(kdVdZ,0), oldkp->GetFilStateValue(kQoverP,0)}; Double_t newstate[] = {0.,0.,0.,0.,0.}; Double_t currentz = oldkp->GetTrackCluster()->GetZPos(); bool isvalid = true; bool swimflag = true; while(isvalid && plnid != plnidend && nplaneskip<2*fParmNSkipActive){ if(plnid.IsValid() && plnid.GetPlaneCoverage()!=PlaneCoverage::kNoActive){ UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid); if(scintpln.IsValid()){ // plane is active, check if in coil swimflag = Swim(currentstate,newstate,currentz,scintpln.GetZ0(),0, &fVldContext); if(!swimflag || TMath::Abs(newstate[kU])>10. || TMath::Abs(newstate[kV])>10.) { // if failed swim, set nplaneskip to very large number nplaneskip = 9999999; isvalid = 0; } else{ currentz = scintpln.GetZ0(); currentstate[kU] = newstate[kU]; currentstate[kV] = newstate[kV]; currentstate[kdUdZ] = newstate[kdUdZ]; currentstate[kdVdZ] = newstate[kdVdZ]; currentstate[kQoverP] = newstate[kQoverP]; } // check if not in coil (defined to have radius 30 cm) if(currentstate[kU]*currentstate[kU]+currentstate[kV]*currentstate[kV]>0.09) { ++nplaneskip; } }//end if scintillator plane is valid } plnid = plnid.GetAdjoinScint(direction); }//end loop over planes return nplaneskip; }

Int_t AlgFitTrackSR::FindNumSkippedPlanesInView (  std::map< Int_t, Int_t > &  fitPlanes, std::map< Int_t, Int_t >::iterator &  fitIter, Int_t  currentPlane, Int_t  prevPlane, Int_t  lastPlane, Int_t  direction )  [private]

Definition at line 3336 of file AlgFitTrackSR.cxx. Referenced by RemoveBadPointsFromFit(). { //assume you drop the current plane, so start off with numSkipped=1 Int_t numSkipped = 1; Int_t plane = 0; if(prevPlane<=0) return 0; fitIter = fitPlanes.find(currentPlane); //only do the check if there was a plane previous to this one if(prevPlane>0){ if(direction == 1) --fitIter; else ++fitIter; plane = (*fitIter).first; //loop the iterator in reverse until you find the previously fit plane while( plane*direction > prevPlane*direction && numSkipped<fParmNSkipView){ if(direction == 1) --fitIter; else ++fitIter; plane = (*fitIter).first; //went through the loop, increment the number of skipped planes ++numSkipped; //stop looking once you find a kp with a good fit if( (*fitIter).second != 0 ) break; } }//end if not the first plane in view //now look downstream if(numSkipped<fParmNSkipView){ //didnt skip too many before the current plane, what about after it? fitIter = fitPlanes.find(currentPlane); if(direction == 1) ++fitIter; else --fitIter; plane = (*fitIter).first; while( plane*direction < lastPlane*direction && numSkipped<fParmNSkipView ){ if(direction == 1) ++fitIter; else --fitIter; plane = (*fitIter).first; //stop looking once you find a kp with a good fit //the break point is before the count increase in //the forward loop so that you //dont double count the current plane as being skipped. you have //set the iterator such that it should be one //plane beyond the current plane //and if that plane is a good fit, you didnt skip 2 planes, just one if( (*fitIter).second != 0) break; //went through the loop, increment the number of skipped planes ++numSkipped; } }//end check for second view return numSkipped; }

Int_t AlgFitTrackSR::FindTimingDirection (  CandFitTrackSRHandle &  cfh, Double_t *  fitparm, Double_t &  timefitchi2 )  [private]

Definition at line 4204 of file AlgFitTrackSR.cxx. References digit(), CandDigitHandle::GetCharge(), CandHandle::GetDaughterIterator(), CandTrackHandle::GetdS(), PlexSEIdAltL::GetEnd(), CandHandle::GetNDaughters(), CandStripHandle::GetPlane(), CandDigitHandle::GetPlexSEIdAltL(), CandTrackHandle::GetT(), ArrivalTime::Weight(), and LinearFit::Weighted(). { ArrivalTime arrtime; Int_t plane = 0; const Double_t min_wgt = 0.4; StripEnd::StripEnd_t stripEnd = StripEnd::kUnknown; CandStripHandle *strip = 0; CandDigitHandle *digit = 0; CandStripHandleItr stripItr(cfh.GetDaughterIterator()); stripItr.Reset(); //arrays for determining the timing direction - limit to //1000 entries - really shouldnt need more than that to //be able to figure out if the event is going north or //south Double_t time[1000]; Double_t pathLength[1000]; Double_t weight[1000]; Double_t digitpe = 0.; Int_t digitCtr = 0; for(Int_t i = 0; i<1000; ++i){ time[i] = 0.; pathLength[i] = 0.; weight[i] = 0.; } //loop over all strips in the track while( (strip = stripItr()) ){ //only look at double ended strips and those that have //believable transverse positions if(strip->GetNDaughters() == 2){ plane = strip->GetPlane(); //get the iterator over the digits in the strip CandDigitHandleItr digitItr(strip->GetDaughterIterator()); while( (digit = digitItr()) ){ stripEnd = digit->GetPlexSEIdAltL().GetEnd(); pathLength[digitCtr] = cfh.GetdS()-cfh.GetdS(plane); //time is recorded in seconds, but i prefer working in ns. //time is from the track not the strips time[digitCtr] = cfh.GetT(plane,stripEnd)*1.e9; digitpe = digit->GetCharge(CalDigitType::kPE); weight[digitCtr] = TMath::Min(arrtime.Weight(digitpe),min_wgt); ++digitCtr; }//end loop over digits in strip }//end if double ended strip }//end loop over strips in track Double_t efitparm[2]; Double_t maxresid = 11.; Double_t resid = 0.; Int_t ctr = digitCtr; Int_t imaxresid = 0; Int_t nremove = 0; while(maxresid>10 && digitCtr-nremove>4 && (1.*nremove)<0.2*digitCtr ){ LinearFit::Weighted(ctr,pathLength,time,weight,fitparm,efitparm); maxresid = 0.; imaxresid = 0; for(int i=0; i<ctr; ++i){ resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i]); if(weight[i]>0. && TMath::Abs(resid)>maxresid){ maxresid = TMath::Abs(resid); imaxresid = i; } } if(maxresid>10){ weight[imaxresid] = 0.; nremove++; } timefitchi2 = 0.; for(int i=0; i<ctr; ++i){ resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i]); timefitchi2 += weight[i]*resid*resid; } }//end loop to find timing fit and remove outliers timefitchi2 = 0.; for(int i=0; i<ctr; ++i){ if(weight[i]>0.){ resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i]); timefitchi2 += weight[i]*resid*resid; } } //check the chi^2 for the fit - if it is really bad dont change the //initial guess at the direction for the event, ie just make sure //that the slope in time vs pathlength is positive; // if(fitparm[1]<0. && timefitchi2>=10.*(digitCtr-nremove)) // fitparm[1] *= -1.; return digitCtr-nremove; }

Int_t AlgFitTrackSR::FindTimingDirection (  TrackClusterSRItr &  clusterItr  )  [private]

Definition at line 4104 of file AlgFitTrackSR.cxx. References TrackClusterSR::GetBegTime(), TrackClusterSR::GetCharge(), TrackClusterSR::GetZPos(), min, ArrivalTime::Weight(), and LinearFit::Weighted(). { Int_t idir = 1; ArrivalTime arrtime; // determine directionality using timing - this should really //only be needed for cosmic ray events, as beam events always //increase in z with increasing time. //we could try to determine the timing for each side (E, W) of each //electronics crate. that would be the most accurate way, but it //would also take a lot of time and iterations. instead, we can //probably get away with just doing a global fit to all datapoints //since we do have timing constants that take care of the side to side //and crate to crate offsets. even if a bit of timing hardware is //exchanged, this wont really affect our ability to go back and figure //out the right calibrations later. //make some arrays to do the fit. allow a maximum of 1000 points //to be included in the fit. there may be events with more than //1000 digits, but that should be enough to tell you which way the //event is going. Double_t zpos[1000]; Double_t time[1000]; Double_t weight[1000]; for (int i=0; i<1000; i++) { zpos[i] = 0.; time[i] = 0.; weight[i] = 0.; } //minimum weight for a signal in the timing fit // maximum weight corresponds to 10 p.e.; reason for this is to minimize // weight of showers in which transverse spread can have negative effect // on timing const Double_t min_wgt = 0.4; TrackClusterSR *cluster = 0; Int_t clusterCtr = 0; clusterItr.ResetFirst(); while( (cluster = clusterItr()) && clusterCtr<1000){ //get the z position of the plane zpos[clusterCtr] = (Double_t)cluster->GetZPos(); time[clusterCtr] = cluster->GetBegTime(); weight[clusterCtr] = min(arrtime.Weight(cluster->GetCharge()),min_wgt); ++clusterCtr; }//end loop over clusters to fill timing arrays //reset the iterator over all strips clusterItr.ResetFirst(); //arrays to hold the slope, intercept, and chi^2 of the //least squares fit for the timing Double_t parm[2],eparm[2]; //set the maximum residual to the input parameter+1ns to //get ready for the while loop below Double_t fParmMaxTimingResid = 20*Munits::ns; Double_t maxresid = fParmMaxTimingResid+1.*Munits::ns; Double_t resid = 0.; Int_t imaxresid = 0; Int_t nremove=0; Int_t fitflag = 0; while(clusterCtr-nremove>4 && maxresid>fParmMaxTimingResid && !fitflag){ fitflag = LinearFit::Weighted(clusterCtr,zpos,time,weight,parm,eparm); maxresid = 0.; imaxresid = 0; for (int i=0; i<clusterCtr; ++i) { resid = parm[0]+parm[1]*zpos[i]-time[i]; if (weight[i]>0. && TMath::Abs(resid)>maxresid) { maxresid = TMath::Abs(resid); imaxresid = i; } } if (maxresid>fParmMaxTimingResid) { weight[imaxresid] = 0.; ++nremove; } }//end loop to fit the timing values if(parm[1]<0.){ idir=-1; } return idir; }

Int_t AlgFitTrackSR::FindUpstreamPlanes (  TObjArray &  planeClusterList, CandFitTrackSRHandle &  cfh, TObjArray *  newkplist, Int_t  direction )  [private]

Definition at line 3161 of file AlgFitTrackSR.cxx. References AddReverseBestKPToFit(), FindNumSkippedPlanes(), FitFrom(), CandFitTrackSRHandle::GetKalmanLast(), CandFitTrackSRHandle::GetKalmanPlane(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetPlane(), KalmanPlaneSR::GetPreChi2(), CandTrackHandle::GetRange(), KalmanPlaneSR::GetTrackCluster(), KalmanPlaneSR::GetZDir(), MSG, CandFitTrackSRHandle::SetCurrentKalmanPlane(), KalmanPlaneSR::SetRange(), and KalmanPlaneSR::SetZDir(). Referenced by DoKalmanFit(). { MSG("FitTrackSR", Msg::kDebug) << "in FindUpstreamPlanes " << direction << endl; TrackClusterSR *thistc = 0; KalmanPlaneSR *thiskp = 0; KalmanPlaneSR *bestkp = 0; KalmanPlaneSR *oldkp = cfh.GetKalmanPlane(cfh.GetKalmanLast()); Int_t nswimfail = 0; Int_t nplaneskip = 0; Bool_t isvalid(1); Int_t istatus2 = 0; Int_t nadd = 0; //loop over the clusters in the cluster list // identify planes which already are fit map<Int_t,KalmanPlaneSR *> fitplane; KalmanPlaneSR *kp = 0; for(int i=0; i<=cfh.GetKalmanLast(); ++i){ kp = cfh.GetKalmanPlane(i); MSG("FitTrackSR",Msg::kDebug) << " plane " << kp->GetTrackCluster()->GetPlane() << " strips " << kp->GetTrackCluster()->GetMinStrip() << "-" << kp->GetTrackCluster()->GetMaxStrip() << " " << kp->GetZDir() << endl; fitplane[kp->GetTrackCluster()->GetPlane()] = kp; } MSG("FitTrackSR",Msg::kDebug) << "After finding end hits" << endl; Int_t currentPlane = oldkp->GetTrackCluster()->GetPlane(); Int_t oldKPPlane = oldkp->GetTrackCluster()->GetPlane(); Int_t prevPlane = -1; Int_t firstPlane = cfh.GetKalmanPlane(0)->GetTrackCluster()->GetPlane(); MSG("FitTrackSR",Msg::kDebug) << "REVERSING FIT, looking for good track clusters" << endl; newkplist->AddLast(oldkp); nswimfail = 0; isvalid = 1; map<Int_t,KalmanPlaneSR *>::iterator planeiter; Bool_t dir = kIterBackward; if(direction!=1) dir=kIterForward; TIter planeItr(&planeClusterList,dir); TObjArray * plnarray; while(isvalid && (plnarray = (TObjArray * )planeItr.Next()) ){ TIter clusterItr(plnarray); thistc = (TrackClusterSR *)clusterItr.Next(); if( (direction==1 && thistc->GetPlane()<oldKPPlane) || (direction!=1 && thistc->GetPlane()>oldKPPlane)){ clusterItr.Reset(); while(isvalid && (thistc = (TrackClusterSR *)clusterItr.Next()) ){ //if this isnt a valid track cluster, go on to the next one if(!thistc->IsValid()) continue; currentPlane = thistc->GetPlane(); oldKPPlane = oldkp->GetTrackCluster()->GetPlane(); planeiter = fitplane.find(thistc->GetPlane()); MSG("FitTrackSR",Msg::kDebug) << "considering trackcluster, plane=" << thistc->GetPlane() << " strips=" << thistc->GetMinStrip() << "-" << thistc->GetMaxStrip() << endl; //if you are on a new plane and have a bestkp , try adding it to the fit if(bestkp && currentPlane != prevPlane){ nadd += AddReverseBestKPToFit(bestkp, oldkp, cfh, newkplist, nswimfail, direction); MSG("FitTrackSR" , Msg::kDebug) << "bestkp has direction " << bestkp->GetZDir() << endl; bestkp = 0; } //see if this plane contains a fit track cluster if(planeiter!=fitplane.end() ){ //this plane contains a fit track cluster, so add it to the //newkplist if(currentPlane != prevPlane ){ kp = dynamic_cast<KalmanPlaneSR *>(planeiter->second); MSG("FitTrackSR", Msg::kDebug) << oldkp->GetTrackCluster()->GetPlane() << " " << kp->GetTrackCluster()->GetPlane() << " " << kp->GetZDir() << endl; nswimfail = FitFrom(oldkp,kp,1,-1); MSG("FitTrackSR", Msg::kDebug) << "adding previously fit plane " << thistc->GetPlane() << " to list" << " " << kp->GetZDir() << endl; newkplist->AddLast(kp); MSG("FitTrackSR", Msg::kDebug) << "added" << endl; cfh.SetCurrentKalmanPlane(kp); oldkp = kp; } }//end if plane contains a fit track cluster else{ MSG("FitTrackSR",Msg::kDebug) << "this plane not fit, considering tc " << thistc->GetPlane() << " " << thistc->GetMinStrip() << "-" << thistc->GetMaxStrip() << " for addition to end" << endl; // need to check how many active planes were actually skipped // only do check if beyond first fit plane nplaneskip=0; //only worry about checking how many planes were skipped //if the currentPlane is upstream of the first plane // - the gaps are of acceptable size between fit planes already if(currentPlane*direction<firstPlane*direction && TMath::Abs(currentPlane-oldKPPlane)>=2*fParmNSkipActive) nplaneskip = FindNumSkippedPlanes(currentPlane, oldKPPlane, oldkp, -direction); //not too many planes skipped, give this cluster a whirl if(nplaneskip<2*fParmNSkipActive){ thiskp = new KalmanPlaneSR(thistc); thiskp->SetRange(cfh.GetRange(thiskp->GetTrackCluster()->GetPlane())); thiskp->SetZDir(oldkp->GetZDir()); MSG("FitTrackSR",Msg::kDebug) << "fitting from track cluster, plane " << oldkp->GetTrackCluster()->GetPlane() << " strips " << oldkp->GetTrackCluster()->GetMinStrip() << "-" << oldkp->GetTrackCluster()->GetMaxStrip() << " direction " << direction << "/" << thiskp->GetZDir() << endl; istatus2 = FitFrom(oldkp,thiskp,1,-1); //if the fit was good and this kp has a better chi^2 than the previous //best one for the plane, make it bestkp MSG("FitTrackSR",Msg::kDebug) << " chi2 = " << thiskp->GetPreChi2(1) << " nswimfail = " << istatus2 << endl; if(!istatus2 && (!bestkp || thiskp->GetPreChi2(1)<bestkp->GetPreChi2(1))){ if (bestkp) delete bestkp; MSG("FitTrackSR",Msg::kDebug) << "best chi2, keeping" << endl; bestkp = thiskp; nswimfail = istatus2; } else delete thiskp; }//end not too many active planes skipped else{ // any planes beyond this one will not satisfy NSkipView, stop search isvalid = 0; } }//end current plane does not have a fit clustser prevPlane = currentPlane; }//end loop in reverse over clusters } } //outside of the loop, but there could still be an extra best kp to add if(bestkp){ nadd += AddReverseBestKPToFit(bestkp, oldkp, cfh, newkplist, nswimfail, direction); MSG("FitTrackSR" , Msg::kDebug) << "bestkp has direction " << bestkp->GetZDir() << endl; bestkp = 0; } return nadd; }

Int_t AlgFitTrackSR::FitFrom (  KalmanPlaneSR *  prevkp, KalmanPlaneSR *  currentkp, Int_t  idir, Int_t  iterate )  [private]

Definition at line 897 of file AlgFitTrackSR.cxx. References Filter(), fParmMaxQP, fParmMaxQPFrac, KalmanPlaneSR::GetFilCovarianceMatrixValue(), KalmanPlaneSR::GetFilStateValue(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetPlane(), KalmanPlaneSR::GetRange(), TrackClusterSR::GetTPos(), KalmanPlaneSR::GetTrackCluster(), KalmanPlaneSR::GetZDir(), kQoverP, MSG, and Predict(). Referenced by AddReverseBestKPToFit(), AddToFit(), FindDownstreamPlanes(), FindUpstreamPlanes(), and ReverseFit(). { // // < 0 = failure // return value = 0 success // Int_t nswimfail = 0; MSG("FitTrackSR",Msg::kDebug) << "in FitFrom: " << " Fitting direction " << idir << " " << currentkp->GetTrackCluster()->GetPlane() << " " << currentkp->GetTrackCluster()->GetTPos() << " min/max strip = " << currentkp->GetTrackCluster()->GetMinStrip() << "/" << currentkp->GetTrackCluster()->GetMaxStrip() << currentkp->GetZDir() << endl; //predict the state of the plane to add to the fit based on the //previously fit plane nswimfail = Predict(prevkp, currentkp, idir); if(nswimfail<0){ return nswimfail; } Int_t istatus = 0; istatus = Filter(currentkp, idir); MSG("FitTrackSR",Msg::kDebug) << "fit point " << idir << " " << " " << currentkp->GetTrackCluster()->GetPlane() << " " << currentkp->GetTrackCluster()->GetTPos() << endl; MSG("FitTrackSR",Msg::kDebug) << "fit qp = " << currentkp->GetFilStateValue(kQoverP,idir) << " sigma q/p= " << currentkp->GetFilCovarianceMatrixValue(0, kQoverP,kQoverP) << "/" << currentkp->GetFilCovarianceMatrixValue(1, kQoverP,kQoverP) << " range (g/cm**2) = " << currentkp->GetRange() << " maximum qp from range = " << fParmMaxQPFrac << "/0.002/" << currentkp->GetRange() << " = " << (currentkp->GetRange()>0. ? fParmMaxQPFrac/0.002/currentkp->GetRange() : 0.) << " maxqpparm = " << fParmMaxQP << endl; return nswimfail; }

void AlgFitTrackSR::InitializeFit (  CandFitTrackSRHandle &  cfh, Int_t  , Int_t  )  [private]

Definition at line 520 of file AlgFitTrackSR.cxx. References fCov, fParmInitialPositionError2, fParmInitialQPError2, fParmInitialSlopeError2, CandFitTrackSRHandle::GetCurrentKalmanPlane(), CandRecoHandle::GetDirCosU(), CandRecoHandle::GetDirCosV(), CandRecoHandle::GetDirCosZ(), CandRecoHandle::GetEndDirCosU(), CandRecoHandle::GetEndDirCosV(), CandRecoHandle::GetEndDirCosZ(), CandFitTrackHandle::GetEndQP(), CandRecoHandle::GetEndU(), CandRecoHandle::GetEndV(), KalmanPlaneSR::GetFilCovarianceMatrixValue(), KalmanPlaneSR::GetFilStateValue(), CandFitTrackSRHandle::GetInitialQP(), CandFitTrackSRHandle::GetKalmanLast(), CandFitTrackSRHandle::GetKalmanPlane(), TrackClusterSR::GetPlane(), KalmanPlaneSR::GetTrackCluster(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), KalmanPlaneSR::GetZDir(), TrackClusterSR::GetZPos(), kdUdZ, kdVdZ, kQoverP, kU, kV, min, MSG, KalmanPlaneSR::SetFilCovarianceMatrix(), KalmanPlaneSR::SetFilCovarianceMatrixValue(), KalmanPlaneSR::SetFilStateValue(), KalmanPlaneSR::SetPredictPlane(), CandRecoHandle::SetVtxPlane(), and CandRecoHandle::SetVtxZ(). Referenced by AddToFit(), DoKalmanFit(), and ReverseFit(). { KalmanPlaneSR *kp = 0; for(Int_t i=0; i<=cfh.GetKalmanLast(); ++i){ kp = cfh.GetKalmanPlane(i); kp->SetPredictPlane(-999,idir); } KalmanPlaneSR *currentkp = cfh.GetCurrentKalmanPlane(); //fill up the covariance matrix for the current kp for(Int_t i=0; i<5; ++i){ for(Int_t j=0; j<5; ++j){ currentkp->SetFilCovarianceMatrixValue(idir,i,j,0.); } } currentkp->SetFilCovarianceMatrixValue(idir,kU,kU,fParmInitialPositionError2); currentkp->SetFilCovarianceMatrixValue(idir,kV,kV,fParmInitialPositionError2); currentkp->SetFilCovarianceMatrixValue(idir,kdUdZ,kdUdZ,fParmInitialSlopeError2); currentkp->SetFilCovarianceMatrixValue(idir,kdVdZ,kdVdZ,fParmInitialSlopeError2); currentkp->SetFilCovarianceMatrixValue(idir,kQoverP,kQoverP,fParmInitialQPError2); Double_t cov_xqp = currentkp->GetFilCovarianceMatrixValue(idir,kQoverP,kQoverP)*0.01*0.3*0.025; cov_xqp *= (Double_t)(currentkp->GetZDir()); currentkp->SetFilCovarianceMatrixValue(idir,kU,kQoverP,cov_xqp); currentkp->SetFilCovarianceMatrixValue(idir,kV,kQoverP,cov_xqp); currentkp->SetFilCovarianceMatrixValue(idir,kQoverP,kU,cov_xqp); currentkp->SetFilCovarianceMatrixValue(idir,kQoverP,kV,cov_xqp); //set the filtered state vector with either the vertex or end //point values depending on whether the fit is in the forward or //reverse direction if(idir==0){ // we begin with the vertex position currentkp->SetFilStateValue(kU,0,cfh.GetVtxU()); currentkp->SetFilStateValue(kV,0,cfh.GetVtxV()); currentkp->SetFilStateValue(kdUdZ,0,cfh.GetDirCosU()/cfh.GetDirCosZ()); currentkp->SetFilStateValue(kdVdZ,0,cfh.GetDirCosV()/cfh.GetDirCosZ()); currentkp->SetFilStateValue(kQoverP,0,cfh.GetInitialQP()); } else{ currentkp->SetFilStateValue(kU,1,cfh.GetEndU()); currentkp->SetFilStateValue(kV,1,cfh.GetEndV()); currentkp->SetFilStateValue(kdUdZ,1,cfh.GetEndDirCosU()/cfh.GetEndDirCosZ()); currentkp->SetFilStateValue(kdVdZ,1,cfh.GetEndDirCosV()/cfh.GetEndDirCosZ()); currentkp->SetFilStateValue(kQoverP,1,cfh.GetEndQP()); } MSG("FitTrackSR",Msg::kDebug) << "Initializing fit, direction = " << idir << " iteration = " << iterate << " plane = " << currentkp->GetTrackCluster()->GetPlane() << " " << currentkp->GetFilStateValue(kU,idir) << " " << currentkp->GetFilStateValue(kV,idir) << " " << currentkp->GetFilStateValue(kdUdZ,idir) << " " << currentkp->GetFilStateValue(kdVdZ,idir) << " " << currentkp->GetFilStateValue(kQoverP,idir) << endl; if(!(idir==0 && iterate==0) && iterate>=0){ fCov[kU][kU] *= fParmCovarianceScale; fCov[kV][kV] *= fParmCovarianceScale; fCov[kdUdZ][kdUdZ] *=fParmCovarianceScale; fCov[kdVdZ][kdVdZ] *= fParmCovarianceScale; fCov[kQoverP][kQoverP] *=fParmCovarianceScale; // make sure none of these are larger than our very first covariance elements fCov[kU][kU] = min(fCov[kU][kU],fParmInitialPositionError2); fCov[kV][kV] = min(fCov[kV][kV],fParmInitialPositionError2); fCov[kdUdZ][kdUdZ] = min(fCov[kdUdZ][kdUdZ],fParmInitialSlopeError2); fCov[kdVdZ][kdVdZ] = min(fCov[kdVdZ][kdVdZ],fParmInitialSlopeError2); fCov[kQoverP][kQoverP] = min(fCov[kQoverP][kQoverP],fParmInitialQPError2); Double_t cov_xqp = fParmInitialQPError2*0.01*0.3*0.025; for(Int_t i=0; i<2; ++i){ if(TMath::Abs(fCov[i][kQoverP])>cov_xqp) fCov[i][kQoverP] = (fCov[i][kQoverP] > 0 ? cov_xqp : -cov_xqp); if(TMath::Abs(fCov[kQoverP][i])>cov_xqp) fCov[kQoverP][i] = (fCov[kQoverP][i] > 0 ? cov_xqp : -cov_xqp); }//end loop over covariance stuff cov_xqp /= 0.06; for(Int_t i=2; i<4; ++i){ if(TMath::Abs(fCov[i][kQoverP])>cov_xqp) fCov[i][kQoverP] = (fCov[i][kQoverP] > 0 ? cov_xqp : -cov_xqp); if(TMath::Abs(fCov[kQoverP][i])>cov_xqp) fCov[kQoverP][i] = (fCov[kQoverP][i] > 0 ? cov_xqp : -cov_xqp); }//end loop to fill covariance matrix stuff //set the filtered covariance matrix of the current kp to be the same as fCov currentkp->SetFilCovarianceMatrix(idir,fCov); }//end if !(idir==0 && iterate==0) && iterate>=0 MsgFormat ffmt("%10.5f"); MSG("FitTrackSR",Msg::kDebug) << "Covariance Matrix\n"; for (Int_t i=0; i<5; i++) { MSG("FitTrackSR",Msg::kDebug) << ffmt(currentkp->GetFilCovarianceMatrixValue(idir,i,kU)) << ffmt(currentkp->GetFilCovarianceMatrixValue(idir,i,kV)) << ffmt(currentkp->GetFilCovarianceMatrixValue(idir,i,kdUdZ)) << ffmt(currentkp->GetFilCovarianceMatrixValue(idir,i,kdVdZ)) << ffmt(currentkp->GetFilCovarianceMatrixValue(idir,i,kQoverP)) << endl; } if(idir==0){ cfh.SetVtxZ(currentkp->GetTrackCluster()->GetZPos()); cfh.SetVtxPlane(currentkp->GetTrackCluster()->GetPlane()); } }

void AlgFitTrackSR::InitKalmanFitParameters (  AlgConfig &  ac  )  [private]

Definition at line 497 of file AlgFitTrackSR.cxx. References fParmdedx, fParmDState, fParmMaxAngle, fParmMaxQP, fParmMaxQPFrac, fParmPlnRadLen, fParmQPRangeCheck, fParmTPosError2, fSwimmer, Registry::GetDouble(), Registry::GetInt(), and MSG. Referenced by RunAlg(). { fParmDState[kU] = ac.GetDouble("KalmanDState1"); fParmDState[kV] = ac.GetDouble("KalmanDState2"); fParmDState[kdUdZ] = ac.GetDouble("KalmanDState3"); fParmDState[kdVdZ] = ac.GetDouble("KalmanDState4"); fParmDState[kQoverP] = ac.GetDouble("KalmanDState5"); fParmPlnRadLen = ac.GetDouble("KalmanPlnRadLen"); fParmdedx = ac.GetDouble("Kalmandedx"); fParmTPosError2 = ac.GetDouble("TPosError2"); fParmMaxQP = ac.GetDouble("KalmanMaxQP"); fParmQPRangeCheck = ac.GetInt("KalmanQPRangeCheck"); fParmMaxQPFrac = ac.GetDouble("KalmanMaxQPFrac"); fParmMaxAngle = ac.GetDouble("KalmanMaxAngle"); fSwimmer = ac.GetInt("Swimmer"); MSG("FitTrackSR", Msg::kDebug) << "initial state = " << fParmDState[kU] << " " << fParmDState[kV] << " " << fParmDState[kdUdZ] << " " << fParmDState[kdVdZ] << " " << fParmDState[kQoverP] << endl; return; }

Int_t AlgFitTrackSR::IterateKalmanFit (  TObjArray &  planeClusterList, CandFitTrackSRHandle &  cfh, Int_t &  nu, Int_t &  nv, Int_t  direction, Int_t  dosearch )  [private]

Definition at line 2403 of file AlgFitTrackSR.cxx. References AddClustersToFit(), CandFitTrackSRHandle::ClearKalmanPlaneList(), fParmMaxIterate, fParmQPDiff, KalmanPlaneSR::GetFilStateValue(), CandFitTrackSRHandle::GetKalmanPlane(), CandFitTrackSRHandle::GetNChangedFitPoint(), MSG, ResetTrackClusterList(), and ReverseFit(). Referenced by DoKalmanFit(). { MSG("AlgFitTrackSR", Msg::kDebug) << " in IterateKalmanFit" << endl; Int_t iterate=0; Double_t qp0=0.; Double_t qp1=0.; nu = 0; nv = 0; Int_t istatus = 0; // iterate until stable MSG("FitTrackSR",Msg::kDebug) << "Before iteration " << iterate << "/" << fParmMaxIterate << " changed points " << cfh.GetNChangedFitPoint() << " qp0 " << qp0 << " qp1 " << qp1 << " istatus " << istatus << endl; Double_t qpDiff = fParmQPDiff+.01; while(iterate<fParmMaxIterate && !istatus && (iterate<=1 || ((qp1!=0. && TMath::Abs(qpDiff)>fParmQPDiff) && !(qp0==0. && qp1==0.)) || cfh.GetNChangedFitPoint()) ){ //get the fit q/p for the forward fit if(iterate) qp0 = cfh.GetKalmanPlane(0)->GetFilStateValue(4,1); cfh.ClearKalmanPlaneList(); nu=0; nv=0; // clock_t dummyt; // struct tms t1; // struct tms t2; // struct tms t3; //static double ticksPerSecond = sysconf(_SC_CLK_TCK); // dummyt = times(&t1); istatus = AddClustersToFit(planeClusterList, cfh, iterate, nu, nv, direction); // dummyt = times(&t2); // cout << " forward fitter time " << (Double_t)(t2.tms_utime+t2.tms_stime-t1.tms_utime-t1.tms_stime)/ticksPerSecond << endl; // cout << " --------------- do reverse fit -------------------" << endl; //check to see if you have enough planes in both views to reverse the //fit if(nu>=3 && nv>=3) istatus = ReverseFit(planeClusterList,cfh,iterate,dosearch); // dummyt = times(&t3); // cout << " reverse fitter time " << (Double_t)(t3.tms_utime+t3.tms_stime-t2.tms_utime-t2.tms_stime)/ticksPerSecond << endl; //since new track clusters may have been chosen in the ReverseFit method, //reset the cluster list ResetTrackClusterList(planeClusterList, cfh); qp1 = cfh.GetKalmanPlane(0)->GetFilStateValue(4,1); if(qp1==0) qp1 = cfh.GetKalmanPlane(0)->GetFilStateValue(4,0); if(qp1 != 0.) qpDiff = 1.-qp0/qp1; else qpDiff = fParmQPDiff; MSG("FitTrackSR",Msg::kDebug) << "After iteration " << iterate << "/" << fParmMaxIterate << " changed points " << cfh.GetNChangedFitPoint() << " qp0 " << qp0 << " qp1 " << qp1 << " istatus " << istatus << " qpdiff " << TMath::Abs(qpDiff) << "/" << fParmQPDiff << endl; ++iterate; }//end loop of iterations return istatus; }

void AlgFitTrackSR::MakeDaughterStripList (  CandFitTrackSRHandle &  cfh  )  [private]

Definition at line 3538 of file AlgFitTrackSR.cxx. References CandHandle::AddDaughterLink(), CandStripHandle::GetCharge(), CandFitTrackSRHandle::GetdUdZ(), CandFitTrackSRHandle::GetdVdZ(), CandFitTrackSRHandle::GetKalmanLast(), CandFitTrackSRHandle::GetKalmanPlane(), CandStripHandle::GetPlane(), TrackClusterSR::GetPlane(), TrackClusterSR::GetPlaneView(), CandStripHandle::GetStrip(), TrackClusterSR::GetStripList(), CandStripHandle::GetTPos(), KalmanPlaneSR::GetTrackCluster(), CandTrackHandle::GetU(), CandTrackHandle::GetV(), CandStripHandle::GetZPos(), MSG, CandTrackHandle::SetInShower(), and tc. Referenced by RunAlg(). { Double_t xz = 0.; Double_t dxdz = 0.; Int_t nstripexp = 0; KalmanPlaneSR *kp = 0; TrackClusterSR *tc = 0; CandStripHandle *strip = 0; const CandStripHandle *constStrip = 0; Int_t minstripindx = 0; Double_t mintpos = 0; Int_t indxbest = 0; Double_t distbest = 0; Double_t distsum = 0.; Int_t strip1=-1,strip2=-1; Double_t dist = 0.; Int_t nshower = 0; Int_t thisplane = 0; for(Int_t i=0; i<=cfh.GetKalmanLast(); ++i){ kp = cfh.GetKalmanPlane(i); tc = kp->GetTrackCluster(); xz = 0.; dxdz = 0.; thisplane = tc->GetPlane(); if(tc->GetPlaneView()==PlaneView::kU){ xz = cfh.GetU(thisplane); dxdz = cfh.GetdUdZ(thisplane); } else if(tc->GetPlaneView()==PlaneView::kV){ xz = cfh.GetV(thisplane); dxdz = cfh.GetdVdZ(thisplane); } //see how many strips you expect to be hit in this plane //assume 4:1 aspect ratio for scintillator strips to planes nstripexp = TMath::Max(1,(Int_t)(TMath::Abs(dxdz)/4.)+1); //if the number of strips in the cluster is less than the expected number, //take them all for the daughter list if(tc->GetStripList()->GetLast()+1<=nstripexp){ for(Int_t istrip=0; istrip<=tc->GetStripList()->GetLast(); ++istrip){ strip = dynamic_cast<CandStripHandle*>(tc->GetStripList()->At(istrip)); cfh.AddDaughterLink(*strip); MSG("AlgFitTrackSR", Msg::kDebug) << "in track plane = " << strip->GetPlane() << "/" << strip->GetZPos() << " strip = " << strip->GetStrip() << "/" << strip->GetTPos() << " charge = " << strip->GetCharge() << " nstripexp " << nstripexp << " nshower = " << nshower << endl; //not in a shower because you expected that many strips cfh.SetInShower(strip,0); } }//end if fewer strips than expected in cluster else{ TObjArray orderedstriplist; vector<Int_t> orderedstripindx; // order strips by tpos minstripindx = 0; //loop over the strips until there are no more to add while(minstripindx>=0){ mintpos = 9999999.; minstripindx = -1; for (Int_t istrip=0; istrip<=tc->GetStripList()->GetLast(); ++istrip){ constStrip = dynamic_cast<const CandStripHandle*>(tc->GetStripList()->At(istrip)); if(orderedstriplist.IndexOf(constStrip)<0 && constStrip->GetTPos()<mintpos){ mintpos = constStrip->GetTPos(); minstripindx = istrip; } }//end loop over strips in cluster if(minstripindx>=0){ orderedstripindx.push_back(minstripindx); strip = dynamic_cast<CandStripHandle*>(tc->GetStripList()->At(minstripindx)); orderedstriplist.Add(strip); } }//end loop over strips to add to list indxbest = 0; distbest = 0; for(Int_t indx=0; indx<=tc->GetStripList()->GetLast()+1-nstripexp; ++indx){ distsum = 0.; strip1 = -1; strip2 = -1; for(Int_t istrip=indx; istrip<indx+nstripexp; ++istrip){ constStrip = dynamic_cast<const CandStripHandle*>(tc->GetStripList()->At(orderedstripindx[istrip])); if (strip1<0) strip1 = constStrip->GetStrip(); strip2 = constStrip->GetStrip(); dist = TMath::Abs(constStrip->GetTPos()-xz); if (fParmIsCosmic) dist *= (0.3+exp(-0.2*constStrip->GetCharge())); distsum += dist; } if (!indx || distsum<distbest) { distbest = distsum; indxbest = indx; } }//end loop over strips in cluster // count how many hit strips above 1 pe nshower=0; for(int i=0; i<=tc->GetStripList()->GetLast(); ++i){ strip = dynamic_cast<CandStripHandle*>(tc->GetStripList()->At(i)); if (strip->GetCharge()>1.5) ++nshower; MSG("AlgFitTrackSR", Msg::kDebug) << "in shower plane = " << strip->GetPlane() << "/" << strip->GetZPos() << " strip = " << strip->GetStrip() << "/" << strip->GetTPos() << " charge = " << strip->GetCharge() << " nstripexp " << nstripexp << " nshower = " << nshower << endl; } for(Int_t istrip=indxbest; istrip<indxbest+nstripexp; ++istrip){ strip = dynamic_cast<CandStripHandle*>(tc->GetStripList()->At(orderedstripindx[istrip])); cfh.AddDaughterLink(*strip); cfh.SetInShower(strip,nshower-nstripexp); } }//end too more strips in cluster than expected }//end loop over kps return; }

void AlgFitTrackSR::MakeSliceClusterList (  TObjArray *  trackClusterList, const TObjArray *  tclist, CandStripHandleItr &  stripItr, CandFitTrackSRHandle &  cfh )  [private]

Definition at line 1719 of file AlgFitTrackSR.cxx. References CandFitTrackSRHandle::AddTrackCluster(), CandStripHandle::GetCharge(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetPlane(), CandStripHandle::GetPlane(), CandStripHandle::GetStrip(), CandStripHandle::GetTPos(), CandStripHandle::GetZPos(), TrackClusterSR::IsContained(), TrackClusterSR::IsValid(), MSG, TrackClusterSR::SetIsValid(), and tc. Referenced by RunAlg(). { stripItr.Reset(); // select track clusters which are in this slice Int_t oldplane = 0; Int_t minstrip = 0; Int_t maxstrip = 0; CandStripHandle *strip = 0; TrackClusterSR *tc = 0; while( (strip = stripItr()) ){ MSG("AlgFitTrackSR", Msg::kDebug) << "plane = " << strip->GetPlane() << "/" << strip->GetZPos() << " strip = " << strip->GetStrip() << "/" << strip->GetTPos() << " charge = " << strip->GetCharge() << endl; if(strip->GetPlane()!=oldplane || strip->GetStrip()<minstrip || strip->GetStrip()>maxstrip) { Bool_t found(0); for(int i=0; !found && i<=tclist->GetLast(); ++i){ tc = dynamic_cast<TrackClusterSR*>(tclist->At(i)); assert(tc); if((found = tc->IsContained(strip))){ TrackClusterSR *newtc = new TrackClusterSR(*tc); newtc->SetIsValid(true); trackClusterList->AddLast(newtc); MSG("FitTrackSR",Msg::kDebug) << "Adding trackcluster " << tc->GetPlane() << " " << tc->GetMinStrip() << " " << tc->GetMaxStrip() << " is valid " << (Int_t)tc->IsValid() << " to list" << endl; cfh.AddTrackCluster(newtc); oldplane = tc->GetPlane(); minstrip = tc->GetMinStrip(); maxstrip = tc->GetMaxStrip(); }//end if this cluster contains the current strip }//end loop over clusters }//end if this strip is in a new cluster }//end loop over strips in slice return; }

void AlgFitTrackSR::MakeTrackClusterList (  TObjArray *  trackClusterList, CandStripHandleItr &  stripItr, CandFitTrackSRHandle &  cfh, Int_t  direction )  [private]

Definition at line 1776 of file AlgFitTrackSR.cxx. References TrackClusterSR::AddStrip(), CandFitTrackSRHandle::AddTrackCluster(), fParmMaxClusterNStrip, fParmMisalignmentError, CandStripHandle::GetBegTime(), TrackClusterSR::GetCharge(), CandStripHandle::GetCharge(), CandStripHandle::GetDemuxVetoFlag(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetNStrip(), TrackClusterSR::GetPlane(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), CandStripHandle::GetStrip(), CandStripHandle::GetTPos(), CandStripHandle::GetZPos(), TrackClusterSR::IsValid(), and MSG. Referenced by RunAlg(). { MSG("AlgFitTrackSR", Msg::kDebug) << "in MakeTrackClusterList" << endl; stripItr.Reset(); // create track clusters Int_t oldplane=0; Int_t oldstrip=0; Double_t oldtime=0.; TrackClusterSR *tcluster=0; CandStripHandle *strip = 0; while( (strip = stripItr()) ){ if(!strip->GetDemuxVetoFlag() && (strip->GetPlaneView()==PlaneView::kU || strip->GetPlaneView()==PlaneView::kV)){ MSG("AlgFitTrackSR", Msg::kDebug) << "plane = " << strip->GetPlane() << "/" << strip->GetZPos() << " strip = " << strip->GetStrip() << "/" << strip->GetTPos() << " charge = " << strip->GetCharge() << endl; if(tcluster && strip->GetPlane()==oldplane && direction*strip->GetStrip()<=direction*oldstrip+3){ tcluster->AddStrip(strip); oldstrip = strip->GetStrip(); oldtime = strip->GetBegTime(); } else{ TrackClusterSR *trackcluster = new TrackClusterSR(strip, fParmMisalignmentError); trackClusterList->AddLast(trackcluster); oldplane = strip->GetPlane(); oldstrip = strip->GetStrip(); oldtime = strip->GetBegTime(); tcluster = trackcluster; }//end if this is a new cluster }//end if strip is in the detector and not vetoed by DeMux }//end loop over strips for(Int_t i=0; i<=trackClusterList->GetLast(); ++i){ tcluster = dynamic_cast<TrackClusterSR*>(trackClusterList->At(i)); //get rid of clusters that are too wide or too little charge if(tcluster->GetNStrip()>fParmMaxClusterNStrip || tcluster->GetCharge()<fParmMinClusterCharge){ trackClusterList->RemoveAt(i); delete tcluster; } else{ MSG("AlgFitTrackSR",Msg::kDebug) << "Adding trackcluster " << tcluster->GetPlane() << " " << tcluster->GetMinStrip() << " " << tcluster->GetMaxStrip() << " is valid " << (Int_t)tcluster->IsValid() << " to list" << endl; cfh.AddTrackCluster(tcluster); } }//end loop over clusters to remove the bad ones trackClusterList->Compress(); return; }

Bool_t AlgFitTrackSR::MarkTrackClusters (  const CandTrackHandle *  track0, TObjArray &  planeClusterList, CandStripHandleItr &  stripItr, CandFitTrackSRHandle &  cfh, Int_t &  begPlane, Int_t &  endPlane )  [private]

Definition at line 1976 of file AlgFitTrackSR.cxx. References CandStripHandle::GetPlane(), CandStripHandle::GetStrip(), CandTrackHandle::GetU(), CandTrackHandle::GetV(), CandTrackHandle::IsInShower(), MSG, ResetTrackClusterList(), CandTrackHandle::SetInShower(), and tc. Referenced by RunAlg(). { bool goodTrack = true; Int_t oldplane=0; stripItr.Reset(); // add trackclusters that are part of track Int_t minstrip = 0; Int_t maxstrip = 0; map<TrackClusterSR*,Bool_t> isintrack; Int_t inTrackPlane[1000]; for(Int_t i = 0; i<1000; ++i){ inTrackPlane[i] = 0; } Int_t ntrku = 0; Int_t ntrkv = 0; Double_t oldu,oldv; Int_t prevplane=-1; begPlane = -1; endPlane = -1; Bool_t found = false; TrackClusterSR *tc = 0; CandStripHandle *strip = 0; CandStripHandle *tcstrip = 0; //set all the track clusters in the iterator to IsValid = false ResetTrackClusterList(planeClusterList); //loop over strips to find the clusters in the track and mark them as valid while( (strip = stripItr()) ){ // check if strip outside detector if(prevplane!=strip->GetPlane()){ MSG("FitTrackSR", Msg::kDebug) << "prev plane = " << prevplane << " strip plane = " << strip->GetPlane() << " strip " << strip->GetStrip() << endl; prevplane = strip->GetPlane(); oldu = track0->GetU(prevplane); oldv = track0->GetV(prevplane); MSG("FitTrackSR",Msg::kDebug) << "PLANE " << prevplane << " u,v = " << oldu << " " << oldv << endl; } //if the current track point is out of the detector, ignore it //B.J.R. 3.24.2005: test the impact parameter for the point to make the call if(TMath::Sqrt(oldu*oldu + oldv*oldv) > fParmMaxImpactParameter){ MSG("FitTrackSR",Msg::kDebug) << "PLANE " << prevplane << " out of detector, u,v = " << oldu << " " << oldv << " REMOVING from fittrack input" << endl; continue; } //make sure you only take one cluster from each plane if(inTrackPlane[strip->GetPlane()]==0 && (strip->GetPlane()!=oldplane || strip->GetStrip()<minstrip || strip->GetStrip()>maxstrip)){ found = false; //loop over the clusters until you find the one containing the //current strip TObjArray * planeClusters=(TObjArray *)planeClusterList.At(strip->GetPlane()); TIter clusterItr(planeClusters); while(!found && (tc = (TrackClusterSR *)(clusterItr.Next()) )){ if( tc->GetPlane()==strip->GetPlane()){ MSG("AlgFitTrackSR", Msg::kDebug) << "cluster on plane " << tc->GetPlane() << " min strip " << tc->GetMinStrip() << " max strip " << tc->GetMaxStrip() << " is valid " << (Int_t)tc->IsValid() << endl; //if this strip is in this cluster if(tc->IsContained(strip)){ found = true; //the cluster is good so set the plane flag to 1 and the cluster //valid flag to true inTrackPlane[strip->GetPlane()] = 1; tc->SetIsValid(true); MSG("AlgFitTrackSR",Msg::kDebug) << "Adding trackcluster " << tc->GetPlane() << " " << tc->GetMinStrip() << " " << tc->GetMaxStrip() << " is valid " << (Int_t)tc->IsValid() << " to track list" << endl; //see what view the cluster is in and if it is a good beg plane if(tc->GetPlaneView()==PlaneView::kU) ++ntrku; else if(tc->GetPlaneView()==PlaneView::kV) ++ntrkv; if(begPlane<0) begPlane = strip->GetPlane(); endPlane = strip->GetPlane(); //loop over all strips in the cluster and add them to the CandFitTrack for(int j=0; j<=tc->GetStripList()->GetLast(); ++j){ tcstrip = dynamic_cast<CandStripHandle*>(tc->GetStripList()->At(j)); cfh.SetInShower(tcstrip,track0->IsInShower(strip)); } oldplane = tc->GetPlane(); minstrip = tc->GetMinStrip(); maxstrip = tc->GetMaxStrip(); }//end if this strip is in the current cluster }//end loop over clusters } }//end if strip is from a different plane/cluster }//end loop over strips stripItr.Reset(); MSG("FitTrackSR",Msg::kDebug) << "# of u,v planes in track = " << ntrku << " " << ntrkv << endl; if(ntrku<3 || ntrkv<3){ goodTrack = false; return goodTrack; } return goodTrack; }

Int_t AlgFitTrackSR::Predict (  KalmanPlaneSR *  prevkp, KalmanPlaneSR *  currentkp, Int_t  idir )  [private]

Definition at line 747 of file AlgFitTrackSR.cxx. References CalculateNoise(), CalculatePreChi2(), CalculatePreCovariance(), CalculatePreState(), CalculatePropagator(), fParmMaxAngle, fParmMaxQPFrac, TrackClusterSR::GetPlane(), KalmanPlaneSR::GetPredictPlane(), KalmanPlaneSR::GetPreStateValue(), KalmanPlaneSR::GetRange(), TrackClusterSR::GetTPosError(), KalmanPlaneSR::GetTrackCluster(), kdUdZ, kdVdZ, kQoverP, MSG, KalmanPlaneSR::SetPredictPlane(), and KalmanPlaneSR::SetPreStateValue(). Referenced by FitFrom(). { // // < 0 indicates failure // return value = 0 success // > 0 indicates number of swim failures // Int_t nswimfail=0; MSG("FitTrackSR", Msg::kDebug) << "in Predict" << endl; //make sure you are looking at different planes if(prevkp->GetTrackCluster()->GetPlane()!=currentkp->GetPredictPlane(idir)) { Int_t swimstatus = CalculatePropagator(prevkp,currentkp,idir); if(swimstatus<0){ MSG("FitTrackSR",Msg::kDebug) << "failed to calculate propagator" << endl; return swimstatus; } nswimfail += swimstatus; //get the Q_{k-1} matrix CalculateNoise(prevkp, currentkp, idir); //get the C_{k}^{k-1} matrix CalculatePreCovariance(prevkp, currentkp, idir); //extrapolate from the previously fit plane to the plane to add //to the fit using the swimmer - this takes the place of the F_{k-1}x_{k-1} //matrix multiplication in the Kalman Filtering equations. swimstatus = CalculatePreState(prevkp, currentkp, idir); if(swimstatus<0){ MSG("FitTrackSR",Msg::kDebug) << "failed to calculate predicted state" << endl; return swimstatus; } nswimfail += swimstatus; }//end if on a different plane than the predict plane CalculatePreChi2(currentkp, idir,currentkp->GetTrackCluster()->GetTPosError()); //set the plane you predicted the fit for the current one from currentkp->SetPredictPlane(prevkp->GetTrackCluster()->GetPlane(), idir); //i am not sure if the following lines are needed or not - they are done in //the CalculatePreState method, but that is only called if the plane to add //to the fit is different plane than the previously fit plane - dunno, guess //i will leave them in //---------------------------------------------------------------------------// if(TMath::Abs(currentkp->GetPreStateValue(kdUdZ, idir))>fParmMaxAngle){ if(currentkp->GetPreStateValue(kdUdZ, idir)>0.) currentkp->SetPreStateValue(kdUdZ, idir, fParmMaxAngle); if(currentkp->GetPreStateValue(kdUdZ, idir)<0.) currentkp->SetPreStateValue(kdUdZ, idir, -fParmMaxAngle); } if(TMath::Abs(currentkp->GetPreStateValue(kdVdZ, idir))>fParmMaxAngle){ if(currentkp->GetPreStateValue(kdVdZ, idir)>0.) currentkp->SetPreStateValue(kdVdZ, idir, fParmMaxAngle); if(currentkp->GetPreStateValue(kdVdZ, idir)<0.) currentkp->SetPreStateValue(kdVdZ, idir, -fParmMaxAngle); } Double_t maxqp = fParmMaxQP; //use a value of about 2 MeV/(g/cm^2) for the energy loss in the material //you would divide the max qp frac allowed by that value and then by the range to //get the max qp. but multiplication is a faster operation, so multiply by //500 instead. if ( fParmQPRangeCheck ) { if(currentkp->GetRange()>0. && fParmMaxQPFrac*500./currentkp->GetRange()<maxqp) { maxqp = fParmMaxQPFrac*500./currentkp->GetRange(); } } if(idir==0 && TMath::Abs(currentkp->GetPreStateValue(kQoverP, idir))>maxqp){ if(currentkp->GetPreStateValue(kQoverP, idir)>0.) currentkp->SetPreStateValue(kQoverP, idir, maxqp); else currentkp->SetPreStateValue(kQoverP, idir, -maxqp); } //---------------------------------------------------------------------------// return nswimfail; }

Int_t AlgFitTrackSR::RemoveBadPointsFromFit (  CandFitTrackSRHandle &  cfh, TObjArray &  planeClusterList, std::map< Int_t, Int_t > &  uFitPlanes, std::map< Int_t, Int_t > &  vFitPlanes, Int_t &  nfitu, Int_t &  nfitv, Int_t  direction )  [private]

Definition at line 2720 of file AlgFitTrackSR.cxx. References FindNumSkippedPlanesInView(), fParmDeltaChi2, fParmDeltaCovariance, fParmNSkipView, KalmanPlaneSR::GetFilChi2(), KalmanPlaneSR::GetFilCovarianceMatrixValue(), CandFitTrackSRHandle::GetKalmanLast(), CandFitTrackSRHandle::GetKalmanPlane(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetPlane(), TrackClusterSR::GetPlaneView(), KalmanPlaneSR::GetTrackCluster(), TrackClusterSR::IsValid(), kdUdZ, kdVdZ, max, min, MSG, ResetTrackClusterList(), TrackClusterSR::SetIsValid(), CandFitTrackSRHandle::SetVtx(), and tc. Referenced by DoKalmanFit(). { MSG("FitTrackSR", Msg::kDebug) << "in RemoveBadPointsFromFit" << endl; TrackClusterSR *tc = 0; const KalmanPlaneSR *vtxkp = 0; Int_t nremoveitr = 0; Int_t nbadfit = 0; nfitu = 0; nfitv = 0; Double_t maxlocalchi2 = fParmMaxLocalChi2; Double_t maxanglecovariance = fParmMaxAngleCovariance; // if chi2 has to be increased due to too many points being removed, //deltachi2 > 0 fParmDeltaChi2 = -1.; // if angle covariance has to be increased due to too many points //being removed, deltacovariance > 0 fParmDeltaCovariance = -1.; Double_t mindchi2=0.; Double_t mindcovariance=0.; Int_t nPlanesSkippedU = 0; Int_t nPlanesSkippedV = 0; Int_t lastUPlane = -1; Int_t lastVPlane = -1; Int_t lastFitPlane = -1; MSG("FitTrackSR", Msg::kDebug) << "nfit u " << nfitu << " nfitv " << nfitv << endl; while(nfitu<3 || nfitv<3){ vtxkp = 0; nfitu = 0; nfitv = 0; nbadfit = 0; //reset the cluster list ResetTrackClusterList(planeClusterList); mindchi2=0.; mindcovariance=0.; nPlanesSkippedU = 0; nPlanesSkippedV = 0; lastUPlane = -1; lastVPlane = -1; KalmanPlaneSR *kp = 0; //loop over the Kalman planes to fill the maps for(Int_t i=0; i<=cfh.GetKalmanLast(); ++i){ kp = cfh.GetKalmanPlane(i); //see if the reset of the clusterList also reset the //clusters in the kp's MSG("FitTrackSR", Msg::kDebug) << "cluster for kp on plane " << kp->GetTrackCluster()->GetPlane() << " is valid " << (Int_t)kp->GetTrackCluster()->IsValid() << endl; MSG("FitTrackSR",Msg::kDebug) << " TC " << kp->GetTrackCluster()->GetPlane() << " " << kp->GetTrackCluster()->GetMinStrip() << " " << kp->GetTrackCluster()->GetMaxStrip() << " " << kp->GetFilChi2(0) << " " << kp->GetFilChi2(1) << "/" << maxlocalchi2 << " " << TMath::Min(kp->GetFilCovarianceMatrixValue(0,kdUdZ, kdUdZ), kp->GetFilCovarianceMatrixValue(1,kdUdZ, kdUdZ)) << " " << TMath::Min(kp->GetFilCovarianceMatrixValue(0,kdVdZ, kdVdZ), kp->GetFilCovarianceMatrixValue(1,kdVdZ, kdVdZ)) << endl; //see if the fit for this plane was really good if(TMath::Max(kp->GetFilChi2(0),kp->GetFilChi2(1))<=maxlocalchi2 && TMath::Min(kp->GetFilCovarianceMatrixValue(0,kdUdZ,kdUdZ), kp->GetFilCovarianceMatrixValue(1,kdUdZ,kdUdZ))<maxanglecovariance && TMath::Min(kp->GetFilCovarianceMatrixValue(0,kdVdZ,kdVdZ), kp->GetFilCovarianceMatrixValue(1,kdVdZ,kdVdZ))<maxanglecovariance) { MSG("FitTrackSR", Msg::kDebug) << " this fit was good" << endl; kp->GetTrackCluster()->SetIsValid(true); MSG("FitTrackSR", Msg::kDebug) << " plane " << kp->GetTrackCluster()->GetPlane() << " IsValid " << kp->GetTrackCluster()->IsValid() << endl; TObjArray * planeClusters=(TObjArray *)planeClusterList.At(kp->GetTrackCluster()->GetPlane()); TIter clusterItr(planeClusters); while((tc = (TrackClusterSR *)(clusterItr.Next()) )){ if(tc->GetPlane()==kp->GetTrackCluster()->GetPlane()){ // replacing slice if(tc->GetMinStrip() == kp->GetTrackCluster()->GetMinStrip() && tc->GetMaxStrip() == kp->GetTrackCluster()->GetMaxStrip() ) tc->SetIsValid(true); } } if(kp->GetTrackCluster()->GetPlaneView() == PlaneView::kU) uFitPlanes[kp->GetTrackCluster()->GetPlane()] = 1; else if(kp->GetTrackCluster()->GetPlaneView() == PlaneView::kV) vFitPlanes[kp->GetTrackCluster()->GetPlane()] = 1; }//end if good fit for kp else{ MSG("FitTrackSR", Msg::kDebug) << " plane " << kp->GetTrackCluster()->GetPlane() << " IsValid " << (Int_t)kp->GetTrackCluster()->IsValid() << endl; //no need to set the fTrackCluster to IsValid=false, as the //reset call above did that. if(kp->GetTrackCluster()->GetPlaneView() == PlaneView::kU) uFitPlanes[kp->GetTrackCluster()->GetPlane()] = 0; else if(kp->GetTrackCluster()->GetPlaneView() == PlaneView::kV) vFitPlanes[kp->GetTrackCluster()->GetPlane()] = 0; } }//end loop to fill map of fit quality for kp's lastFitPlane = cfh.GetKalmanPlane(cfh.GetKalmanLast())->GetTrackCluster()->GetPlane(); map<Int_t, Int_t>::iterator vfitIter = vFitPlanes.begin(); map<Int_t, Int_t>::iterator ufitIter = uFitPlanes.begin(); for(Int_t i=0; i<=cfh.GetKalmanLast(); ++i){ kp = cfh.GetKalmanPlane(i); //look to see which view the kp comes from, and if it has a //good fit, as found above if(kp->GetTrackCluster()->GetPlaneView()==PlaneView::kU){ if(uFitPlanes[kp->GetTrackCluster()->GetPlane()]==1){ ++nfitu; lastUPlane = kp->GetTrackCluster()->GetPlane(); if(!vtxkp) vtxkp = kp; } else{ //see how many planes were skipped in this view nPlanesSkippedU = FindNumSkippedPlanesInView(uFitPlanes, ufitIter, kp->GetTrackCluster()->GetPlane(), lastUPlane, lastFitPlane, direction); if(nPlanesSkippedU>=fParmNSkipView){ //dont remove this plane as there would be too many skipped planes //in the view if we did TObjArray * planeClusters=(TObjArray *)planeClusterList.At(kp->GetTrackCluster()->GetPlane()); TIter clusterItr(planeClusters); while((tc = (TrackClusterSR * )(clusterItr.Next()) )){ if(tc->GetPlane()==kp->GetTrackCluster()->GetPlane()){ if(tc->GetMinStrip() == kp->GetTrackCluster()->GetMinStrip() && tc->GetMaxStrip() == kp->GetTrackCluster()->GetMaxStrip() ) tc->SetIsValid(true); } } MSG("FitTrackSR",Msg::kDebug) << "too many consecutive planes on plane " << kp->GetTrackCluster()->GetPlane()<< " # skipped " << nPlanesSkippedU <<"/" << fParmNSkipView << " keeping plane" << endl; if (!vtxkp) vtxkp = kp; lastUPlane = kp->GetTrackCluster()->GetPlane(); ++nfitu; } else{ MSG("FitTrackSR",Msg::kDebug) << "REMOVING, plane " << kp->GetTrackCluster()->GetPlane() << " nskipped = " << nPlanesSkippedU << "/" << fParmNSkipView << endl; ++nbadfit; if(TMath::Max(kp->GetFilChi2(0),kp->GetFilChi2(1))>maxlocalchi2){ if(mindchi2<=0. || TMath::Max(kp->GetFilChi2(0),kp->GetFilChi2(1))-maxlocalchi2<mindchi2){ mindchi2 = TMath::Max(kp->GetFilChi2(0),kp->GetFilChi2(1))-maxlocalchi2; mindchi2 += 0.001; // precision }//end if the max chi^2 is bigger than allowed value } else{ Double_t thismaxcov = TMath::Max(TMath::Min(kp->GetFilCovarianceMatrixValue(0,kdUdZ,kdUdZ), kp->GetFilCovarianceMatrixValue(1,kdUdZ,kdUdZ)), TMath::Min(kp->GetFilCovarianceMatrixValue(0,kdVdZ,kdVdZ), kp->GetFilCovarianceMatrixValue(1,kdVdZ,kdVdZ))); if (mindcovariance<=0. || thismaxcov-maxanglecovariance<mindcovariance) { mindcovariance = thismaxcov-maxanglecovariance; mindcovariance += 0.001; // precision }//end if mindcovariance is too small }//end else to check covariance }//end else remove the current plane }//end else the current kp is a bad fit }//end if in U view //check the other view else if(kp->GetTrackCluster()->GetPlaneView()==PlaneView::kV){ if(vFitPlanes[kp->GetTrackCluster()->GetPlane()]==1){ ++nfitv; lastVPlane = kp->GetTrackCluster()->GetPlane(); if(!vtxkp) vtxkp = kp; } else{ //see how many planes were skipped in this view nPlanesSkippedV = FindNumSkippedPlanesInView(vFitPlanes, vfitIter, kp->GetTrackCluster()->GetPlane(), lastVPlane, lastFitPlane, direction); if(nPlanesSkippedV>=fParmNSkipView){ //dont remove this plane as there would be too many skipped planes //in the view if we did // clusterItr.GetSet()->Slice(kp->GetTrackCluster()->GetPlane()); //in the view if we did MSG("FitTrackSR",Msg::kDebug) << "too many consecutive planes on plane " << kp->GetTrackCluster()->GetPlane()<< " # skipped " << nPlanesSkippedU <<"/" << fParmNSkipView << " keeping plane" << endl; TObjArray * planeClusters=(TObjArray *)planeClusterList.At(kp->GetTrackCluster()->GetPlane()); TIter clusterItr(planeClusters); while((tc = (TrackClusterSR * )(clusterItr.Next()) )){ if(tc->GetPlane()==kp->GetTrackCluster()->GetPlane()){ if(tc->GetMinStrip() == kp->GetTrackCluster()->GetMinStrip() && tc->GetMaxStrip() == kp->GetTrackCluster()->GetMaxStrip() ) tc->SetIsValid(true); } } // clusterItr.GetSet()->Slice(); if (!vtxkp) vtxkp = kp; lastVPlane = kp->GetTrackCluster()->GetPlane(); ++nfitv; } else{ MSG("FitTrackSR",Msg::kDebug) << "REMOVING, plane " << kp->GetTrackCluster()->GetPlane() << " nskipped = " << nPlanesSkippedV << "/" << fParmNSkipView << endl; ++nbadfit; if(max(kp->GetFilChi2(0),kp->GetFilChi2(1))>maxlocalchi2){ if(mindchi2<=0. || max(kp->GetFilChi2(0),kp->GetFilChi2(1))-maxlocalchi2<mindchi2){ mindchi2 = max(kp->GetFilChi2(0),kp->GetFilChi2(1))-maxlocalchi2; mindchi2 += 0.001; // precision }//end if the max chi^2 is bigger than allowed value } else{ Double_t thismaxcov = max(min(kp->GetFilCovarianceMatrixValue(0,kdUdZ,kdUdZ), kp->GetFilCovarianceMatrixValue(1,kdUdZ,kdUdZ)), min(kp->GetFilCovarianceMatrixValue(0,kdVdZ,kdVdZ), kp->GetFilCovarianceMatrixValue(1,kdVdZ,kdVdZ))); if (mindcovariance<=0. || thismaxcov-maxanglecovariance<mindcovariance) { mindcovariance = thismaxcov-maxanglecovariance; mindcovariance += 0.001; // precision }//end if mindcovariance is too small }//end else to check covariance }//end else remove the current plane }//end else the current kp is a bad fit }//end if in V view }//end loop over kp's to remove bad points if (nfitu<3 || nfitv<3) { // not enough points, increase maximum chi2 maxlocalchi2 += mindchi2; maxanglecovariance += mindcovariance; MSG("FitTrackSR",Msg::kDebug) << "TOO FEW POINTS LEFT, increasing chi2" << " threshold to " << maxlocalchi2 << " and angle covariance threshold to " << maxanglecovariance << endl; }//end if less than 3 planes fit in either view fParmDeltaChi2 = mindchi2; fParmDeltaCovariance = mindcovariance; ++nremoveitr; assert(nremoveitr<=cfh.GetKalmanLast()+2); // in principle can never loop through more than # of planes in fit }//end while less than 3 fit planes in at least one view //make sure you found a vertex kp and set the vertex to it. assert(vtxkp); cfh.SetVtx(vtxkp); TIter planeItr(&planeClusterList); TObjArray * plnarray; while( (plnarray = (TObjArray * )planeItr.Next()) ){ TIter clusterItr(plnarray); while( (tc = (TrackClusterSR * )clusterItr.Next()) ){ MSG("FitTrackSR", Msg::kDebug) << "current cluster on plane " << tc->GetPlane() << " is valid " << (Int_t)tc->IsValid() << " " << tc->GetMinStrip() << "-" << tc->GetMaxStrip() << endl; } } return nbadfit; }

void AlgFitTrackSR::ResetTrackClusterList (  TObjArray &  planeClusterList, CandFitTrackSRHandle &  cfh )  [private]

Definition at line 2170 of file AlgFitTrackSR.cxx. References CandFitTrackSRHandle::GetKalmanLast(), CandFitTrackSRHandle::GetKalmanPlane(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetPlane(), KalmanPlaneSR::GetTrackCluster(), TrackClusterSR::SetIsValid(), and tc. { // clusterItr.Reset(); TrackClusterSR *tc = 0; KalmanPlaneSR *kp = 0; //loop over the KalmanPlaneSR objects in the fit - if the cluster //for a given KalmanPlaneSR is the same as the current cluster, //set IsValid = true for(Int_t i = 0; i<=cfh.GetKalmanLast(); ++i){ kp = cfh.GetKalmanPlane(i); TObjArray * planeClusters=(TObjArray *)planeClusterList.At(kp->GetTrackCluster()->GetPlane()); TIter clusterItr(planeClusters); //loop over all the track clusters and set their valid flags to false while( (tc = (TrackClusterSR * )clusterItr.Next()) ){ if(tc->GetPlane()==kp->GetTrackCluster()->GetPlane()){ tc->SetIsValid(false); if(tc->GetMinStrip() == kp->GetTrackCluster()->GetMinStrip() && tc->GetMaxStrip() == kp->GetTrackCluster()->GetMaxStrip() ) tc->SetIsValid(true); } } }//end loop over KalmanPlaneSR's return; }

void AlgFitTrackSR::ResetTrackClusterList (  TObjArray &  planeClusterList  )  [private]

Definition at line 2147 of file AlgFitTrackSR.cxx. References TrackClusterSR::GetPlane(), TrackClusterSR::IsValid(), MSG, TrackClusterSR::SetIsValid(), and tc. Referenced by DoKalmanFit(), IterateKalmanFit(), MarkTrackClusters(), and RemoveBadPointsFromFit(). { // clusterItr.Reset(); TrackClusterSR *tc = 0; //loop over all the track clusters and set their valid flags to false TIter planeItr(&planeClusterList); TObjArray * plnarray; while( (plnarray = (TObjArray * )planeItr.Next()) ){ TIter clusterItr(plnarray); while( (tc = (TrackClusterSR *)clusterItr.Next()) ){ tc->SetIsValid(false); MSG("FitTrackSR", Msg::kDebug) << "cluster on plane " << tc->GetPlane() << " is valid " << (Int_t)tc->IsValid() << endl; } } // clusterItr.Reset(); return; }

Int_t AlgFitTrackSR::ReverseFit (  TObjArray &  planeClusterList, CandFitTrackSRHandle &  cfh, Int_t  , Bool_t  = kFALSE )  [private]

Definition at line 2482 of file AlgFitTrackSR.cxx. References fCov, FitFrom(), CandFitTrackSRHandle::GetBadFit(), KalmanPlaneSR::GetFilChi2(), KalmanPlaneSR::GetFilCovarianceMatrixValue(), KalmanPlaneSR::GetFilStateValue(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), CandFitTrackSRHandle::GetNChangedFitPoint(), CandFitTrackHandle::GetNSwimFail(), TrackClusterSR::GetPlane(), CandFitTrackSRHandle::GetPlaneList(), TrackClusterSR::GetPlaneView(), KalmanPlaneSR::GetPreChi2(), CandTrackHandle::GetRange(), TrackClusterSR::GetTPos(), KalmanPlaneSR::GetTrackCluster(), CandFitTrackSRHandle::GetTrackClusterList(), KalmanPlaneSR::GetZDir(), InitializeFit(), max, MSG, pRev, pRev2, CandFitTrackSRHandle::SetCurrentKalmanPlane(), CandFitTrackSRHandle::SetNChangedFitPoint(), CandFitTrackHandle::SetNSwimFail(), CandFitTrackHandle::SetPlaneChi2(), CandFitTrackSRHandle::SetPlanePreChi2(), CandFitTrackHandle::SetPlaneQP(), KalmanPlaneSR::SetRange(), CandFitTrackSRHandle::SetVtx(), and KalmanPlaneSR::SetZDir(). Referenced by IterateKalmanFit(). { InitializeFit(cfh,1,iterate); TObjArray *KalmanPlaneList =cfh.GetPlaneList(); TObjArray newlist; newlist.Clear(); newlist.Compress(); cfh.SetNChangedFitPoint(0); TObjArray *tclist = cfh.GetTrackClusterList(); KalmanPlaneSR *kpold = dynamic_cast<KalmanPlaneSR*>(KalmanPlaneList->At(KalmanPlaneList->GetLast())); KalmanPlaneSR *lastkp = kpold; Bool_t *kpswimfail = new Bool_t[KalmanPlaneList->GetLast()+1]; for(Int_t i=0 ;i<=KalmanPlaneList->GetLast(); ++i){ kpswimfail[i] = 0; } Int_t nu = 0; Int_t nv = 0; Bool_t didswimfail = false; Int_t nswimfail = 0; Int_t nswimfailnewer = 0; KalmanPlaneSR *kpnew = 0; KalmanPlaneSR *kpnewer = 0; TrackClusterSR *thistc = 0; Bool_t ready = true; for(Int_t i=KalmanPlaneList->GetLast()-1; i>=0; --i){ kpnew = dynamic_cast<KalmanPlaneSR*>(KalmanPlaneList->At(i)); MSG("FitTrackSR",Msg::kDebug) << "new TC, plane = " << kpnew->GetTrackCluster()->GetPlane() << " tpos = " << kpnew->GetTrackCluster()->GetTPos() << " minstrip = " << kpnew->GetTrackCluster()->GetMinStrip() << " maxstrip = " << kpnew->GetTrackCluster()->GetMaxStrip() << " q/p = " << kpnew->GetFilStateValue(4,0) << endl; if(TMath::Abs(kpnew->GetFilStateValue(4,0))<=2.0){ MSG("FitTrackSR", Msg::kDebug) << "current kp has sufficient momentum to do " << "reverse fit" << endl; pRev++; nswimfail = FitFrom(kpold,kpnew,1,iterate); if(nswimfail<0){ MSG("FitTrackSR",Msg::kDebug) << "failed swim" << endl; return nswimfail; } //if(i==KalmanPlaneList->GetLast()-1){ if(ready){ ready = false; cfh.SetPlaneChi2(kpold->GetTrackCluster()->GetPlane(),kpold->GetFilChi2(0)); cfh.SetPlanePreChi2(kpold->GetTrackCluster()->GetPlane(),kpold->GetPreChi2(0)); cfh.SetPlaneQP(kpold->GetTrackCluster()->GetPlane(),kpold->GetFilStateValue(4,1)); } MSG("FitTrackSR",Msg::kDebug) << "prechi2 = " << kpnew->GetPreChi2(1) << " filchi2 = " << max(kpnew->GetFilChi2(0),kpnew->GetFilChi2(1)) << endl; if(dosearch && iterate<2){ kpnewer = 0; nswimfailnewer = 0; MSG("FitTrackSR",Msg::kDebug) << "Searching for other track clusters in plane " << kpnew->GetTrackCluster()->GetPlane() << " original min,max strips = " << kpnew->GetTrackCluster()->GetMinStrip() << " " << kpnew->GetTrackCluster()->GetMaxStrip() << endl; //tclist is the list of TrackClusterSR objects for this track TObjArray * planeClusters=(TObjArray *)planeClusterList.At(kpnew->GetTrackCluster()->GetPlane()); TIter clusterItr(planeClusters); while((thistc = dynamic_cast<TrackClusterSR *>(clusterItr.Next()) )){ //is this a good tc and from the same plane but a different //cluster than the one currently used in the plane? if(!cfh.GetBadFit(thistc) && thistc->GetMinStrip()!=kpnew->GetTrackCluster()->GetMinStrip()){ //make a new kp based on the alternate track cluster and set //its range and z direction KalmanPlaneSR *thiskp = new KalmanPlaneSR(thistc); thiskp->SetRange(cfh.GetRange(thiskp->GetTrackCluster()->GetPlane())); thiskp->SetZDir(kpold->GetZDir()); //try fitting from the previous kp to the alternate one in //the reverse direction pRev2++; nswimfailnewer = FitFrom(kpold,thiskp, 1,iterate); MSG("FitTrackSR",Msg::kDebug) << "considering track cluster: " << tclist->GetLast() << " tpos = " << thistc->GetTPos() << " minstrip = " << thistc->GetMinStrip() << "/" << kpnew->GetTrackCluster()->GetMinStrip() << " prechi2 = " << thiskp->GetPreChi2(1) << " filchi2 = " << thiskp->GetFilChi2(1) << endl; if(nswimfailnewer>=0 && ((kpnewer && (nswimfail<0 || thiskp->GetPreChi2(1)<kpnewer->GetPreChi2(1))) || (!kpnewer && thiskp->GetPreChi2(1)<kpnew->GetPreChi2(1)))){ MSG("FitTrackSR",Msg::kDebug) << " better track cluster found\n"; if(kpnewer) delete kpnewer; kpnewer = thiskp; nswimfail = nswimfailnewer; }//end if nswimfailnewer>=0 etc else delete thiskp; }//end if !GetBadFit() etc }//end loop over track clusters //put the better kp in the list if(kpnewer){ KalmanPlaneList->AddAt(kpnewer,i); delete kpnew; kpnew = kpnewer; cfh.SetNChangedFitPoint(cfh.GetNChangedFitPoint()+1); } }//end if dosearch if(nswimfail<0){ MSG("FitTrackSR",Msg::kDebug) << "failed swim" << endl; return nswimfail; } if(nswimfail!=0 && nu>=3 && nv>=3){ kpswimfail[i] = 1; didswimfail = 1; MSG("FitTrackSR",Msg::kDebug) << "Swim failed, skipping plane" << endl; } else{ //swim didnt fail, set the chi^2 and qp values for the plane cfh.SetPlaneChi2(kpnew->GetTrackCluster()->GetPlane(), max(kpnew->GetFilChi2(1), kpnew->GetFilChi2(0))); cfh.SetPlanePreChi2(kpnew->GetTrackCluster()->GetPlane(), max(kpnew->GetPreChi2(1), kpnew->GetPreChi2(0))); cfh.SetPlaneQP(kpnew->GetTrackCluster()->GetPlane(),kpnew->GetFilStateValue(4,1)); cfh.SetCurrentKalmanPlane(kpnew); lastkp = kpnew; //set the covariance matrix array based on the filtered covariance //matrix of the current kp for(int i1=0; i1<5; ++i1){ for(int i2=0; i2<5; ++i2){ fCov[i1][i2]=kpnew->GetFilCovarianceMatrixValue(1,i1,i2); } } if(kpnew->GetTrackCluster()->GetPlaneView()==PlaneView::kU) ++nu; if(kpnew->GetTrackCluster()->GetPlaneView()==PlaneView::kV) ++nv; kpold = kpnew; }//end else for if nswimfail!=0 && nu>=3 && nv>=3 cfh.SetNSwimFail(cfh.GetNSwimFail()+TMath::Abs(nswimfail)); }//end if enough momentum left to do fit // else kpold = kpnew; }//end loop over kalman planes in reverse // remove any planes which failed swim if(didswimfail){ Int_t nlast = KalmanPlaneList->GetLast(); TObjArray newlist; newlist.Clear(); newlist.Compress(); KalmanPlaneSR *kp = 0; for(Int_t i=0; i<=nlast; ++i){ kp = dynamic_cast<KalmanPlaneSR*>(KalmanPlaneList->At(i)); if(kpswimfail[i]) delete kp; else{ newlist.Add(kp); MSG("FitTrackSR",Msg::kDebug) << "adding thiskp plane " << kp->GetTrackCluster()->GetPlane() << " min/max strip " << kp->GetTrackCluster()->GetMinStrip() << "/" << kp->GetTrackCluster()->GetMaxStrip() << " q/p " << kp->GetFilStateValue(4,1) << endl; }//end else kp didnt fail swim }//end loop over kalman planes KalmanPlaneList->Clear(); KalmanPlaneList->Compress(); nlast = newlist.GetLast(); MSG("FitTrackSR",Msg::kDebug) << "new KP list" << endl; for(Int_t i=0; i<=nlast; ++i){ kp = dynamic_cast<KalmanPlaneSR*>(newlist.At(i)); KalmanPlaneList->Add(kp); MSG("FitTrackSR",Msg::kDebug) << " plane " << kp->GetTrackCluster()->GetPlane() << " min/max strip " << kp->GetTrackCluster()->GetMinStrip() << "/" << kp->GetTrackCluster()->GetMaxStrip() << " q/p " << kp->GetFilStateValue(4,1) << endl; }//end loop over kalman planes to add to the list for the track }//end if didswimfail delete [] kpswimfail; if(lastkp) cfh.SetVtx(lastkp); return 0; }

void AlgFitTrackSR::RunAlg (  AlgConfig &  ac, CandHandle &  ch, CandContext &  cx )  [virtual]

Implements AlgBase. Definition at line 83 of file AlgFitTrackSR.cxx. References CandHandle::AddDaughterLink(), AlgReco::Calibrate(), CandFitTrackSRHandle::ClearKalmanPlaneList(), CandFitTrackSRHandle::ClearMaps(), CandFitTrackSRHandle::ClearTrackClusterList(), DoKalmanFit(), fDetector, fParmCovarianceScale, fParmDeltaChi2, fParmDeltaCovariance, fParmInitialPositionError2, fParmInitialQPError2, fParmInitialSlopeError2, fParmIsCosmic, fParmLastPlane, fParmMaxAngleCovariance, fParmMaxClusterNStrip, fParmMaxImpactParameter, fParmMaxIterate, fParmMaxIterate2, fParmMaxLocalChi2, fParmMaxLocalPreChi2, fParmMinClusterCharge, fParmMinPlanePE, fParmMisalignmentError, fParmNSkipActive, fParmNSkipView, fParmPassMinPlaneAsymmetry, fParmPassPlaneAsymmetry, fParmPassReducedChi2, fParmQPDiff, fParmTPosError2, fVldContext, CandRecord::GetCandHeader(), CandHandle::GetCandRecord(), CandRecoHandle::GetCandSlice(), CandStripHandle::GetCharge(), CandFitTrackHandle::GetChi2(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), Registry::GetDouble(), KalmanPlaneSR::GetFilStateValue(), Registry::GetInt(), CandFitTrackSRHandle::GetKalmanLast(), CandFitTrackSRHandle::GetKalmanPlane(), GetMomFromRange(), CandFitTrackHandle::GetNDOF(), CandRecoHandle::GetNPlane(), CandStripHandle::GetPlane(), TrackClusterSR::GetPlane(), PlexPlaneId::GetPlane(), UgliGeomHandle::GetPlaneIdFromZ(), CandTrackHandle::GetRange(), CandHeader::GetSnarl(), CandRecoHandle::GetTimeSlope(), CandHandle::GetVldContext(), InitKalmanFitParameters(), TrackClusterSR::IsValid(), KeyOnClusterPlane(), MakeDaughterStripList(), MakeSliceClusterList(), MakeTrackClusterList(), MarkTrackClusters(), MSG, pFor, pRev, pRev2, CandHandle::RemoveDaughter(), AlgTrack::SetdS(), CandFitTrackHandle::SetEMCharge(), CandFitTrackSRHandle::SetEnd(), CandFitTrackHandle::SetFinderTrack(), CandTrackHandle::SetMomentum(), CandFitTrackHandle::SetMomentumCurve(), CandFitTrackHandle::SetNIterate(), CandFitTrackHandle::SetPass(), SetPlaneParameters(), AlgTrack::SetT(), SetTrackEndParameters(), SetTrackParameters(), AlgTrack::SetUVZ(), CandFitTrackSRHandle::SetVtx(), SwimVertexAndEndPoints(), and tc. { // check inputs pRev=0; pRev2=0; pFor=0; assert(ch.InheritsFrom("CandFitTrackSRHandle")); CandFitTrackSRHandle &cfh = dynamic_cast<CandFitTrackSRHandle &>(ch); cfh.SetPass(0); assert(cx.GetDataIn()); assert(cx.GetDataIn()->InheritsFrom("TObjArray")); CandTrackHandle *track0 = 0; const TObjArray *tclist = 0; const TObjArray *cxin = dynamic_cast<const TObjArray *>(cx.GetDataIn()); for (Int_t i=0; i<=cxin->GetLast(); i++) { TObject *tobj = cxin->At(i); if (tobj->InheritsFrom("CandTrackHandle")) { track0 = dynamic_cast<CandTrackHandle*>(tobj); MSG("EventNum", Msg::kDebug) << "event number " << track0->GetCandRecord()->GetCandHeader()->GetSnarl() << endl; } if (tobj->InheritsFrom("TObjArray")) { tclist = dynamic_cast<TObjArray*>(tobj); } } assert(track0); // set finder track cfh.SetFinderTrack(track0); // get detector type const VldContext* vldcptr = track0->GetVldContext(); assert(vldcptr); fVldContext = *vldcptr; fDetector = fVldContext.GetDetector(); UgliGeomHandle ugh(fVldContext); // input algorithm parameters fParmMaxIterate = ac.GetInt("MaxIterate"); fParmMaxIterate2 = ac.GetInt("MaxIterate2"); fParmQPDiff = ac.GetDouble("QPDiff"); fParmMinPlanePE = ac.GetDouble("MinPlanePE"); fParmLastPlane = ac.GetInt("LastPlane"); fParmIsCosmic = ac.GetInt("IsCosmic"); fParmMaxLocalChi2 = ac.GetDouble("MaxLocalChi2"); fParmMaxLocalPreChi2 = ac.GetDouble("MaxLocalPreChi2"); fParmMaxAngleCovariance = ac.GetDouble("MaxAngleCovariance"); fParmNSkipView = ac.GetInt("NSkipView"); fParmNSkipActive = ac.GetInt("NSkipActive"); fParmPassReducedChi2 = ac.GetDouble("PassReducedChi2"); fParmPassPlaneAsymmetry = ac.GetDouble("PassPlaneAsymmetry"); fParmPassMinPlaneAsymmetry = ac.GetInt("PassMinPlaneAsymmetry"); fParmMaxImpactParameter = ac.GetDouble("MaxImpactParameter"); fParmMinClusterCharge = ac.GetDouble("MinClusterCharge"); fParmMaxClusterNStrip = ac.GetInt("MaxClusterNStrip"); fParmDeltaChi2 = ac.GetDouble("DeltaChi2"); fParmDeltaCovariance = ac.GetDouble("DeltaCovariance"); fParmMisalignmentError = ac.GetDouble("MisalignmentError")*Munits::mm; fParmTPosError2 = ac.GetDouble("TPosError2"); fParmInitialPositionError2 = ac.GetDouble("InitialPositionError2"); fParmInitialSlopeError2 = ac.GetDouble("InitialSlopeError2"); fParmInitialQPError2 = ac.GetDouble("InitialQPError2"); fParmMaxLocalChi2 = ac.GetDouble("MaxLocalChi2"); fParmMaxLocalPreChi2 = ac.GetDouble("MaxLocalPreChi2"); fParmCovarianceScale = ac.GetDouble("CovarianceScale"); //initialize the parameters from registry for the Kalman Plane stuff InitKalmanFitParameters(ac); // verify that looping won't ask for planes that don't exist const float way_downstream = 999999.; PlexPlaneId lastpln = ugh.GetPlaneIdFromZ(way_downstream); if (lastpln.GetPlane() < fParmLastPlane ) { fParmLastPlane = lastpln.GetPlane(); } Int_t idir = -1; if(track0->GetTimeSlope()>0. || !fParmIsCosmic) idir = 1; //get an iterator over the strips in this slice CandStripHandle *strip = 0; CandStripHandleItr stripItr(track0->GetCandSlice()->GetDaughterIterator()); //sort by plane and strip then by time using the navigation 2 sort functionality stripItr.SetSort2Fun(CandStripHandle::StripSRKeyFromPSEId, CandStripHandle::StripSRKeyFromBegTime); MSG("AlgFitTrackSR", Msg::kDebug) << "sorted the strips by time" << endl; stripItr.SetDirection(idir); stripItr.Reset(); //create a TObjArray to hold the cluster objects you are going to make TObjArray *trackClusterList = new TObjArray(1,0); //fill the trackClusterList with the objects in tclist, if it is there if(tclist) MakeSliceClusterList(trackClusterList, tclist, stripItr, cfh); else MakeTrackClusterList(trackClusterList, stripItr, cfh, idir); MSG("AlgFitTrackSR", Msg::kDebug) << "made cluster list" << endl; // clone the trackClusterList so that you can have a master list of // track clusters and one that keeps track of which clusters are used // in the fit track TObjArray *clusterList(trackClusterList); //make an iterator over clusterList and sort it by plane number TrackClusterSRItr clusterItr(clusterList); TrackClusterSRKeyFunc *clusterKf = clusterItr.CreateKeyFunc(); clusterKf->SetFun(KeyOnClusterPlane); clusterItr.GetSet()->AdoptSortKeyFunc(clusterKf); clusterKf = 0; clusterItr.SetDirection(idir); clusterItr.ResetFirst(); TrackClusterSR *cluster = 0; TObjArray planeClusterList(1000); Int_t iplane=0; Bool_t first=true; TObjArray * planeClusters; while( (cluster = clusterItr())) { if(cluster->GetPlane()!=iplane || first){ iplane=cluster->GetPlane(); MSG("FitTrackSR", Msg::kDebug) << "new cluster on plane " << iplane << endl; planeClusters = new TObjArray(1,0); planeClusterList.AddAt(planeClusters,iplane); planeClusters->Add(cluster); first=false; } else{ MSG("FitTrackSR", Msg::kDebug) << "adding cluster on plane " << iplane << endl; planeClusters=(TObjArray *)planeClusterList.At(iplane); if(!planeClusters) MSG("AlgFitTrackSR", Msg::kFatal) << "planeCluster array not defined for plane " << iplane << endl; planeClusters->Add(cluster); } } Bool_t dir = kIterForward; if(idir!=1) dir=kIterBackward; TIter planeItr(&planeClusterList,idir); TrackClusterSR * tc=0; TObjArray * plnarray; while( (plnarray = (TObjArray * )planeItr.Next()) ){ TIter clusterItr(plnarray); while( (tc = (TrackClusterSR *)clusterItr.Next()) ){ MSG("FitTrackSR", Msg::kDebug) << "current cluster on plane " << tc->GetPlane() << " is valid " << (Int_t)tc->IsValid() << endl; } } MSG("AlgFitTrackSR", Msg::kDebug) << "made cluster iterator" << endl; //get another iterator over the strips in the track this time CandStripHandleItr trackStripItr(track0->GetDaughterIterator()); //sort by plane and strip then by time using the navigation 2 sort functionality trackStripItr.SetSort2Fun(CandStripHandle::StripSRKeyFromPSEId, CandStripHandle::StripSRKeyFromBegTime); /* for now leave in the old code to sort in 2 dimensions just in case CandStripHandleKeyFunc *trackStripKf = trackStripItr.CreateKeyFunc(); //using a 2d sort - first sort on the PlexStripEndId trackStripKf->SetFun(CandStripHandle::StripSRKeyFromPSEId); trackStripItr.GetSet()->AdoptSortKeyFunc(trackStripKf,kTRUE,kFALSE); //now sort by time trackStripKf = trackStripItr.CreateKeyFunc(); trackStripKf->SetFun(CandStripHandle::StripSRKeyFromBegTime); trackStripItr.GetSet()->AdoptSortKeyFunc(trackStripKf,kFALSE); trackStripKf = 0; */ MSG("AlgFitTrackSR", Msg::kDebug) << "got iterator sorted over track strips" << endl; trackStripItr.SetDirection(idir); trackStripItr.Reset(); Int_t begPlane = -1; Int_t endPlane = -1; //now mark the clusters to use in the track by setting the IsValid flag to true bool goodTrack = MarkTrackClusters(track0, planeClusterList, trackStripItr, cfh, begPlane, endPlane); //set the track parameters based on the track passed into the algorithm SetTrackParameters(track0, cfh, idir); //set the endpoints and their direction cosines if track0 is a CandTrackSR SetTrackEndParameters(begPlane, endPlane, cfh, track0); MSG("AlgFitTrackSR", Msg::kDebug) << "start the fitting" << endl; //this track has fewer than 3 hits at least one view, so just add the strips //in the daughter list of track0 to cfh and quit if(!goodTrack){ MSG("AlgFitTrackSR", Msg::kDebug) << "not a good track, bail" << endl; // add daughter strips from the current track and bail CandStripHandleItr trk0stripItr(track0->GetDaughterIterator()); while( (strip = trk0stripItr()) ){ cfh.AddDaughterLink(*strip); } //need to set the end points return; } else{ Int_t istatus = 0; //now to actually do the track fitting. Int_t dosearch=1; Int_t niterate = DoKalmanFit(planeClusterList, cfh, istatus, idir, dosearch); if(istatus>=0){ //get the first and last planes for the track KalmanPlaneSR *kp0 = cfh.GetKalmanPlane(0); KalmanPlaneSR *kp1 = cfh.GetKalmanPlane(cfh.GetKalmanLast()); cfh.SetVtx(kp0); cfh.SetEnd(kp1); //set up variables for first and last planes in each view for track const KalmanPlaneSR *kpu0 = 0; const KalmanPlaneSR *kpv0 = 0; const KalmanPlaneSR *kpu1 = 0; const KalmanPlaneSR *kpv1 = 0; // clear STL maps--we only want fit planes as keys cfh.ClearMaps(); SetPlaneParameters(cfh, kpu0, kpu1, kpv0, kpv1); } // check to see whether full fit converged. If not, fall back to using // only track points. if(istatus<0 || cfh.GetNDOF()==0 || (cfh.GetNDOF()!=0 && cfh.GetChi2()/cfh.GetNDOF()>fParmPassReducedChi2) ) { dosearch = 0; idir = -1; if(track0->GetTimeSlope()>0. || ! fParmIsCosmic) idir = 1; //fill the trackClusterList with the objects in tclist, if it is there cfh.ClearTrackClusterList(); trackClusterList->Clear(); trackClusterList->Compress(); if(tclist) MakeSliceClusterList(trackClusterList, tclist, stripItr, cfh); else MakeTrackClusterList(trackClusterList, stripItr, cfh, idir); //now mark the clusters to use in the track by setting the IsValid flag to true begPlane = -1; endPlane = -1; goodTrack = MarkTrackClusters(track0, planeClusterList, trackStripItr, cfh, begPlane, endPlane); //set the track parameters based on the track passed into the algorithm SetTrackParameters(track0, cfh, idir); //set the endpoints and their direction cosines if track0 is a CandTrackSR SetTrackEndParameters(begPlane, endPlane, cfh, track0); niterate = DoKalmanFit(planeClusterList, cfh, istatus, idir, dosearch); } //set the number of iterations for the fit cfh.SetNIterate(niterate); if(istatus>=0){ //get the first and last planes for the track KalmanPlaneSR *kp0 = cfh.GetKalmanPlane(0); KalmanPlaneSR *kp1 = cfh.GetKalmanPlane(cfh.GetKalmanLast()); cfh.SetVtx(kp0); cfh.SetEnd(kp1); //set up variables for first and last planes in each view for track const KalmanPlaneSR *kpu0 = 0; const KalmanPlaneSR *kpv0 = 0; const KalmanPlaneSR *kpu1 = 0; const KalmanPlaneSR *kpv1 = 0; // clear STL maps--we only want fit planes as keys cfh.ClearMaps(); SetPlaneParameters(cfh, kpu0, kpu1, kpv0, kpv1); //gotta add those strips to the daughter list for this track MakeDaughterStripList(cfh); // create array of strips not in fit track Double_t planepe[1000]; for(int i=0; i<1000; ++i){ planepe[i] = 0.; } stripItr.ResetFirst(); CandStripHandle *strip = 0; while( (strip = stripItr()) ){ if(strip->GetPlane()>=0) planepe[strip->GetPlane()] += strip->GetCharge(CalDigitType::kPE); }//end loop over strips to get charge for each plane // Find planes closest to vertex and end with adjacent planes meeting // PE cut Int_t plane0 = -999; Int_t plane1 = -999; for(int i=0; i<1000; ++i){ if(planepe[i]>=fParmMinPlanePE && planepe[i+1]>=fParmMinPlanePE){ if(idir>0){ if(plane0<0) plane0 = i; if(plane1<0 || i+1>plane1) plane1 = i+1; } else{ if(plane1<0) plane1 = i; if(plane0<0 || i+1>plane0) plane0 = i+1; } }//end if enough signal on these two planes }//end loop over planes to find the first and last ones Double_t vtxqp = kp0->GetFilStateValue(4,1); // tweak qp to force agreement with truth !!!! vtxqp *= 1.01+0.1*TMath::Abs(vtxqp); //if we are looking at cosmic rays, need to swim the track out of the //detector if(fParmIsCosmic) SwimVertexAndEndPoints(cfh, kp0, kp1, kpu0, kpu1, kpv0, kpv1, planepe, plane0, plane1, vtxqp, idir); if(vtxqp==0.){ // infinite momentum, choose finite number to prevent floating point exceptions vtxqp = 0.000001; // 1 PeV cfh.SetEMCharge(0.); } // set momentum and charge sign cfh.SetMomentumCurve(1./TMath::Abs(vtxqp)); if(vtxqp>0.) cfh.SetEMCharge(1.); else if(vtxqp<0.) cfh.SetEMCharge(-1.); // fill time and range maps SetT(&cfh); SetdS(&cfh); // set the momentum from range; taken from PDG R/M vs p/M plot for Fe Double_t range = cfh.GetRange(); // if(range>0.) cfh.SetMomentumRange(0.105658*exp(1.0363*log(range)-4.383)); // until June 2006 used to be if(range>0.) cfh.SetMomentum(0.048 + 1.660e-3*range + 3.057e-8*range*range); if (range>0.){ cfh.SetMomentum(GetMomFromRange(range)); } else cfh.SetMomentum(0.); Calibrate(&cfh); } // end if istatus>=0 cfh.ClearTrackClusterList(); cfh.ClearKalmanPlaneList(); //check and see if the track passes the cuts for a good track Float_t nplaneu = (Float_t)(cfh.GetNPlane(PlaneView::kU)); Float_t nplanev = (Float_t)(cfh.GetNPlane(PlaneView::kV)); Int_t nplane = cfh.GetNPlane(); if(istatus>=0 && cfh.GetNDOF()>0 && cfh.GetChi2()/cfh.GetNDOF()<=fParmPassReducedChi2 && fParmDeltaChi2<=0. && fParmDeltaCovariance<=0. && (nplane<fParmPassMinPlaneAsymmetry || TMath::Abs((nplaneu-nplanev)/(1.*nplane))<=fParmPassPlaneAsymmetry ) ) cfh.SetPass(1); else { cfh.SetPass(0); SetTrackParameters(track0, cfh, idir); CandStripHandleItr trkstripItr(cfh.GetDaughterIterator()); while( (strip = trkstripItr()) ){ cfh.RemoveDaughter(strip); } CandStripHandleItr trk0stripItr(track0->GetDaughterIterator()); while( (strip = trk0stripItr()) ){ cfh.AddDaughterLink(*strip); } // fill time and range maps SetT(&cfh); SetUVZ(&cfh); SetdS(&cfh); Calibrate(&cfh); } }//end if enough planes to try fitting the track //delete the track clusters you found earlier trackClusterList->Delete(); delete trackClusterList; for (Int_t iplane=0;iplane<1000;iplane++){ planeClusters=(TObjArray *)planeClusterList.At(iplane); if(planeClusters){ delete planeClusters; } } // cout << " # for swim " << pFor << " # rev swim " << pRev << " # rev search swim " << pRev2 << endl; return; }

void AlgFitTrackSR::SetPlaneParameters (  CandFitTrackSRHandle &  cfh, const KalmanPlaneSR *  kpu0, const KalmanPlaneSR *  kpu1, const KalmanPlaneSR *  kpv0, const KalmanPlaneSR *  kpv1 )  [private]

Definition at line 3970 of file AlgFitTrackSR.cxx. References KalmanPlaneSR::GetFilChi2(), KalmanPlaneSR::GetFilCovarianceMatrixValue(), KalmanPlaneSR::GetFilStateValue(), CandFitTrackSRHandle::GetKalmanLast(), CandFitTrackSRHandle::GetKalmanPlane(), TrackClusterSR::GetPlane(), TrackClusterSR::GetPlaneView(), KalmanPlaneSR::GetPreChi2(), TrackClusterSR::GetTPosError(), KalmanPlaneSR::GetTrackCluster(), TrackClusterSR::GetZPos(), kdUdZ, kdVdZ, kQoverP, kU, kV, MSG, CandFitTrackHandle::SetChi2(), CandFitTrackSRHandle::SetdUdZ(), CandFitTrackSRHandle::SetdVdZ(), CandFitTrackHandle::SetEnddUError(), CandFitTrackHandle::SetEnddVError(), CandFitTrackHandle::SetEndQPError(), CandFitTrackHandle::SetEndUError(), CandFitTrackHandle::SetEndVError(), CandFitTrackHandle::SetNDOF(), CandFitTrackHandle::SetPlaneChi2(), CandFitTrackSRHandle::SetPlanePreChi2(), CandFitTrackHandle::SetPlaneQP(), CandTrackHandle::SetTrackPointError(), CandTrackHandle::SetU(), CandTrackHandle::SetV(), CandFitTrackHandle::SetVtxdUError(), CandFitTrackHandle::SetVtxdVError(), CandFitTrackHandle::SetVtxQPError(), CandFitTrackHandle::SetVtxUError(), CandFitTrackHandle::SetVtxVError(), and tc. Referenced by RunAlg(). { //loop over the kalman planes to get the first and last ones in each view, //as well as setting the u, v, du/dz, and dv/dz values for each plane const KalmanPlaneSR *kp = 0; TrackClusterSR *tc = 0; Int_t thisplane = 0; Double_t thisu[2]; Double_t thisv[2]; Double_t erru2[2]; Double_t errv2[2]; Double_t avgu = 0.; Double_t avgv = 0.; Double_t avgdu = 0.; Double_t avgdv = 0.; Double_t thisdu[2]; Double_t thisdv[2]; Double_t errdu2[2]; Double_t errdv2[2]; Double_t chi2=0.; Double_t dchi2 = 0.; Double_t dprechi2 = 0.; Double_t thisqp[2]; Double_t thiseqp2[2]; Double_t avgqp = 0.; for(Int_t index=0; index<=cfh.GetKalmanLast(); ++index){ kp = cfh.GetKalmanPlane(index); MSG("AlgFitTrackSR", Msg::kDebug) << "plane positions " << kp->GetTrackCluster()->GetPlane() << " " << kp->GetTrackCluster()->GetZPos() << " " << kp->GetFilStateValue(kU,0) << " " << kp->GetFilStateValue(kV,0) << endl; if(kp->GetTrackCluster()->GetPlaneView()==PlaneView::kU) { if(!kpu0) kpu0 = kp; kpu1 = kp; } if(kp->GetTrackCluster()->GetPlaneView()==PlaneView::kV){ if(!kpv0) kpv0 = kp; kpv1 = kp; } tc = kp->GetTrackCluster(); thisplane = tc->GetPlane(); //loop over the clusters in the track cfh.SetTrackPointError(thisplane,tc->GetTPosError()); //get the forward and backward values for u and v and the errors //on those values. the errors have minimum value of 1.e-12 thisu[0] = kp->GetFilStateValue(kU,0); thisu[1] = kp->GetFilStateValue(kU,1); thisv[0] = kp->GetFilStateValue(kV,0); thisv[1] = kp->GetFilStateValue(kV,1); erru2[0] = TMath::Max(kp->GetFilCovarianceMatrixValue(0,kU,kU), 1.e-12); erru2[1] = TMath::Max(kp->GetFilCovarianceMatrixValue(1,kU,kU), 1.e-12); errv2[0] = TMath::Max(kp->GetFilCovarianceMatrixValue(0,kV,kV), 1.e-12); errv2[1] = TMath::Max(kp->GetFilCovarianceMatrixValue(1,kV,kV), 1.e-12); //average forward and backward fit values avgu = (thisu[0]/erru2[0]+thisu[1]/erru2[1])/(1./erru2[0]+1./erru2[1]); avgv = (thisv[0]/errv2[0]+thisv[1]/errv2[1])/(1./errv2[0]+1./errv2[1]); //set the u an v values for this plane cfh.SetU(thisplane,avgu); cfh.SetV(thisplane,avgv); thisdu[0] = kp->GetFilStateValue(kdUdZ,0); thisdu[1] = kp->GetFilStateValue(kdUdZ,1); thisdv[0] = kp->GetFilStateValue(kdVdZ,0); thisdv[1] = kp->GetFilStateValue(kdVdZ,1); errdu2[0] = TMath::Max(kp->GetFilCovarianceMatrixValue(0,kdUdZ,kdUdZ),1.e-12); errdu2[1] = TMath::Max(kp->GetFilCovarianceMatrixValue(1,kdUdZ,kdUdZ),1.e-12); errdv2[0] = TMath::Max(kp->GetFilCovarianceMatrixValue(0,kdVdZ,kdVdZ),1.e-12); errdv2[1] = TMath::Max(kp->GetFilCovarianceMatrixValue(1,kdVdZ,kdVdZ),1.e-12); // average forward and backward fit values avgdu = (thisdu[0]/errdu2[0]+thisdu[1]/errdu2[1])/(1./errdu2[0]+1./errdu2[1]); avgdv = (thisdv[0]/errdv2[0]+thisdv[1]/errdv2[1])/(1./errdv2[0]+1./errdv2[1]); //set the du/dz and dv/dz values for this plane cfh.SetdUdZ(thisplane,avgdu); cfh.SetdVdZ(thisplane,avgdv); //now do the q/p values dchi2 = TMath::Max(kp->GetFilChi2(0),kp->GetFilChi2(1)); dprechi2 = TMath::Max(kp->GetPreChi2(0),kp->GetPreChi2(1)); cfh.SetPlaneChi2(thisplane,dchi2); cfh.SetPlanePreChi2(thisplane,dprechi2); thisqp[0] = kp->GetFilStateValue(kQoverP,0); thisqp[1] = kp->GetFilStateValue(kQoverP,1); thiseqp2[0] = TMath::Max(kp->GetFilCovarianceMatrixValue(0,kQoverP,kQoverP),1.e-12); thiseqp2[1] = TMath::Max(kp->GetFilCovarianceMatrixValue(1,kQoverP,kQoverP),1.e-12); if(thisqp[0]==0){ avgqp=thisqp[1]; } if(thisqp[1]==0){ avgqp=thisqp[0]; } else avgqp = (thisqp[0]/thiseqp2[0]+thisqp[1]/thiseqp2[1])/(1./thiseqp2[0]+1./thiseqp2[1]); cfh.SetPlaneQP(thisplane,avgqp); chi2 += dchi2; }//end loop over kalman planes to set u,v,du/dz,dv/dz and find first and last //planes in each view //set the track chi^2 for the q/p fit cfh.SetChi2(chi2); Int_t nkalmanlast = cfh.GetKalmanLast(); //set the degrees of freedom and the errors for the vtx and end cfh.SetNDOF(nkalmanlast-4); cfh.SetVtxUError(TMath::Sqrt(cfh.GetKalmanPlane(0)->GetFilCovarianceMatrixValue(1, kU, kU))); cfh.SetVtxVError(TMath::Sqrt(cfh.GetKalmanPlane(0)->GetFilCovarianceMatrixValue(1, kV, kV))); cfh.SetVtxdUError(TMath::Sqrt(cfh.GetKalmanPlane(0)->GetFilCovarianceMatrixValue(1, kdUdZ, kdUdZ))); cfh.SetVtxdVError(TMath::Sqrt(cfh.GetKalmanPlane(0)->GetFilCovarianceMatrixValue(1, kdVdZ, kdVdZ))); cfh.SetVtxQPError(TMath::Sqrt(cfh.GetKalmanPlane(0)->GetFilCovarianceMatrixValue(1,kQoverP,kQoverP))); cfh.SetEndUError(TMath::Sqrt(cfh.GetKalmanPlane(nkalmanlast)->GetFilCovarianceMatrixValue(0,kU,kU))); cfh.SetEndVError(TMath::Sqrt(cfh.GetKalmanPlane(nkalmanlast)->GetFilCovarianceMatrixValue(0,kV,kV))); cfh.SetEnddUError(TMath::Sqrt(cfh.GetKalmanPlane(nkalmanlast)->GetFilCovarianceMatrixValue(0,kdUdZ,kdUdZ))); cfh.SetEnddVError(TMath::Sqrt(cfh.GetKalmanPlane(nkalmanlast)->GetFilCovarianceMatrixValue(0,kdVdZ,kdVdZ))); cfh.SetEndQPError(TMath::Sqrt(cfh.GetKalmanPlane(nkalmanlast)->GetFilCovarianceMatrixValue(0,kQoverP,kQoverP))); return; }

void AlgFitTrackSR::SetTrackEndParameters (  Int_t  begPlane, Int_t  endPlane, CandFitTrackSRHandle &  cfh, const CandTrackHandle *  track0 )  [private]

Definition at line 2108 of file AlgFitTrackSR.cxx. References CandRecoHandle::GetEndPlane(), CandTrackHandle::GetT(), CandTrackHandle::GetU(), CandTrackHandle::GetV(), CandRecoHandle::GetVtxPlane(), CandTrackHandle::GetZ(), CandTrackHandle::IsTPosValid(), CandRecoHandle::SetDirCosU(), CandRecoHandle::SetDirCosV(), CandRecoHandle::SetDirCosZ(), CandRecoHandle::SetEndDirCosU(), CandRecoHandle::SetEndDirCosV(), CandRecoHandle::SetEndDirCosZ(), CandRecoHandle::SetEndPlane(), CandRecoHandle::SetEndT(), CandRecoHandle::SetEndU(), CandRecoHandle::SetEndV(), CandRecoHandle::SetEndZ(), CandRecoHandle::SetVtxPlane(), CandRecoHandle::SetVtxT(), CandRecoHandle::SetVtxU(), CandRecoHandle::SetVtxV(), and CandRecoHandle::SetVtxZ(). Referenced by RunAlg(). { cfh.SetVtxPlane(begPlane); cfh.SetVtxZ(track0->GetZ(begPlane)); if(begPlane!=cfh.GetVtxPlane() && track0->IsTPosValid(begPlane)) { cfh.SetVtxU(track0->GetU(begPlane)); cfh.SetVtxV(track0->GetV(begPlane)); cfh.SetVtxT(track0->GetT(begPlane)); if(track0->InheritsFrom("CandTrackSRHandle")){ // CandTrackSR calculates angles at each plane, use those cfh.SetDirCosU(dynamic_cast<const CandTrackSRHandle *>(track0)->GetDirCosU(begPlane)); cfh.SetDirCosV(dynamic_cast<const CandTrackSRHandle *>(track0)->GetDirCosV(begPlane)); cfh.SetDirCosZ(dynamic_cast<const CandTrackSRHandle *>(track0)->GetDirCosZ(begPlane)); } } cfh.SetEndPlane(endPlane); cfh.SetEndZ(track0->GetZ(endPlane)); if(endPlane!=cfh.GetEndPlane() && track0->IsTPosValid(endPlane)){ cfh.SetEndU(track0->GetU(endPlane)); cfh.SetEndV(track0->GetV(endPlane)); cfh.SetEndT(track0->GetT(endPlane)); if(track0->InheritsFrom("CandTrackSRHandle")){ // CandTrackSR calculates angles at each plane, use those cfh.SetEndDirCosU(dynamic_cast<const CandTrackSRHandle *>(track0)->GetDirCosU(endPlane)); cfh.SetEndDirCosV(dynamic_cast<const CandTrackSRHandle *>(track0)->GetDirCosV(endPlane)); cfh.SetEndDirCosZ(dynamic_cast<const CandTrackSRHandle *>(track0)->GetDirCosZ(endPlane)); } } return; }

void AlgFitTrackSR::SetTrackParameters (  const CandTrackHandle *  track0, CandFitTrackSRHandle &  cfh, Int_t  direction )  [private]

Definition at line 1847 of file AlgFitTrackSR.cxx. References fParmMaxQP, CandRecoHandle::GetCandSlice(), CandRecoHandle::GetDirCosU(), CandRecoHandle::GetDirCosV(), CandRecoHandle::GetDirCosZ(), CandTrackHandle::GetdS(), CandRecoHandle::GetEndDirCosU(), CandRecoHandle::GetEndDirCosV(), CandRecoHandle::GetEndDirCosZ(), CandRecoHandle::GetEndPlane(), CandRecoHandle::GetEndT(), CandRecoHandle::GetEndU(), CandRecoHandle::GetEndV(), CandRecoHandle::GetEndZ(), CandTrackHandle::GetMomentum(), CandTrackSRHandle::GetNTimeFitDigit(), CandTrackSRHandle::GetNTrackDigit(), CandTrackSRHandle::GetNTrackStrip(), CandTrackHandle::GetRange(), CandTrackSRHandle::GetTimeFitChi2(), CandRecoHandle::GetTimeOffset(), CandRecoHandle::GetTimeSlope(), CandTrackHandle::GetTrackPointError(), CandRecoHandle::GetVtxPlane(), CandRecoHandle::GetVtxT(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), CandRecoHandle::GetVtxZ(), CandTrackHandle::IsTPosValid(), CandRecoHandle::SetCandSlice(), CandRecoHandle::SetDirCosU(), CandRecoHandle::SetDirCosV(), CandRecoHandle::SetDirCosZ(), CandTrackHandle::SetdS(), CandRecoHandle::SetEndDirCosU(), CandRecoHandle::SetEndDirCosV(), CandRecoHandle::SetEndDirCosZ(), CandRecoHandle::SetEndPlane(), CandFitTrackHandle::SetEndQP(), CandRecoHandle::SetEndT(), CandRecoHandle::SetEndU(), CandRecoHandle::SetEndV(), CandRecoHandle::SetEndZ(), CandFitTrackSRHandle::SetInitialQP(), CandFitTrackHandle::SetMomentumRange(), CandFitTrackHandle::SetNIterate(), CandTrackHandle::SetNTimeFitDigit(), CandTrackHandle::SetNTrackDigit(), CandTrackHandle::SetNTrackStrip(), CandTrackHandle::SetRange(), CandTrackHandle::SetTimeFitChi2(), CandRecoHandle::SetTimeOffset(), CandRecoHandle::SetTimeSlope(), CandTrackHandle::SetTrackPointError(), CandRecoHandle::SetVtxPlane(), CandRecoHandle::SetVtxT(), CandRecoHandle::SetVtxU(), CandRecoHandle::SetVtxV(), and CandRecoHandle::SetVtxZ(). Referenced by RunAlg(). { cfh.SetCandSlice(track0->GetCandSlice()); cfh.SetDirCosU(track0->GetDirCosU()); cfh.SetDirCosV(track0->GetDirCosV()); cfh.SetDirCosZ(track0->GetDirCosZ()); cfh.SetVtxU(track0->GetVtxU()); cfh.SetVtxV(track0->GetVtxV()); cfh.SetVtxZ(track0->GetVtxZ()); cfh.SetVtxT(track0->GetVtxT()); cfh.SetVtxPlane(track0->GetVtxPlane()); cfh.SetEndDirCosU(track0->GetEndDirCosU()); cfh.SetEndDirCosV(track0->GetEndDirCosV()); cfh.SetEndDirCosZ(track0->GetEndDirCosZ()); cfh.SetEndU(track0->GetEndU()); cfh.SetEndV(track0->GetEndV()); cfh.SetEndZ(track0->GetEndZ()); cfh.SetEndT(track0->GetEndT()); cfh.SetEndPlane(track0->GetEndPlane()); cfh.SetTimeSlope(track0->GetTimeSlope()); cfh.SetTimeOffset(track0->GetTimeOffset()); cfh.SetMomentumRange(track0->GetMomentum()); Double_t initQP=0.0; if(track0->GetMomentum()!=0 && !fParmIsCosmic) initQP=-1./track0->GetMomentum(); cfh.SetInitialQP(initQP); if(!fParmIsCosmic) cfh.SetEndQP(fParmMaxQP); // set ds, range for(Int_t iplane = cfh.GetVtxPlane(); iplane*direction<=cfh.GetEndPlane()*direction; iplane+=direction){ if(track0->IsTPosValid(iplane)){ // u and v coordinates have been calculated for this plane //loop over the clusters in the track cfh.SetTrackPointError(iplane,track0->GetTrackPointError(iplane)); cfh.SetdS(iplane,track0->GetdS(iplane)); cfh.SetRange(iplane,track0->GetRange(iplane)); } } //set stuff if the track was of the SR variety const CandTrackSRHandle *tracksr = 0; if (track0->InheritsFrom("CandTrackSRHandle")) { tracksr = dynamic_cast<const CandTrackSRHandle*>(track0); cfh.SetNTrackStrip(tracksr->GetNTrackStrip()); cfh.SetNTrackDigit(tracksr->GetNTrackDigit()); cfh.SetNTimeFitDigit(tracksr->GetNTimeFitDigit()); cfh.SetTimeFitChi2(tracksr->GetTimeFitChi2()); } cfh.SetNIterate(-1); return; }

Bool_t AlgFitTrackSR::Swim (  Double_t *  swimstate, Double_t *  state, Double_t  zPos, Double_t  finalZ, Int_t  idir, const VldContext *  )  [private]

Definition at line 1590 of file AlgFitTrackSR.cxx. References done(), fParmMaxAngle, SwimParticle::GetCharge(), SwimParticle::GetDirection(), SwimParticle::GetMomentumModulus(), SwimParticle::GetPosition(), kdUdZ, kdVdZ, kInvSqrt2, kQoverP, MSG, s(), SwimParticle::SetCharge(), SwimSwimmer::SwimBackward(), and SwimSwimmer::SwimForward(). Referenced by CalculatePreState(), CalculatePropagator(), FindNumSkippedPlanes(), and SwimVertexAndEndPoints(). { // return false if angles get too large or momentum is too small //------------------------------------------------------------------- // Using Swimmer package //------------------------------------------------------------------- SwimSwimmer s(*vldc); // Add an action to the stepper - in this case a debug printout // SwimPrintStepAction printStep; // s.GetStepper()->AddStepAction(&printStep); // Setup the initial condition for the stepper //the swimmer works in xyz coordinates, not uvz, so need to rotate //the coordinates before putting them into a position vector Double_t x = kInvSqrt2*(swimstate[kU]-swimstate[kV]); Double_t y = kInvSqrt2*(swimstate[kU]+swimstate[kV]); TVector3 position(x,y,zPos); double charge = 0.0; if(swimstate[kQoverP]>0.0) charge = 1.0; else if(swimstate[kQoverP]<0.0) charge = -1.0; else{ // infinite momentum, straight line extrapolation Double_t delz = finalZ - zPos; state[kU] = swimstate[kU] + swimstate[kdUdZ]*delz; state[kV] = swimstate[kV] + swimstate[kdVdZ]*delz; state[kdUdZ] = swimstate[kdUdZ]; state[kdVdZ] = swimstate[kdVdZ]; state[kQoverP] = swimstate[kQoverP]; return 1; } double pz; // forward in time: idir=0; backward in time: idir = 1 // Momentum is always forward in time if((finalZ>zPos && idir==0) || (finalZ<zPos && idir==1)) pz = 1.0/(TMath::Sqrt(1+swimstate[kdUdZ]*swimstate[kdUdZ] +swimstate[kdVdZ]*swimstate[kdVdZ]) *TMath::Abs(swimstate[kQoverP])); else if((finalZ<zPos && idir==0) || (finalZ>zPos && idir==1)) pz = -1.0/(TMath::Sqrt(1+swimstate[kdUdZ]*swimstate[kdUdZ] +swimstate[kdVdZ]*swimstate[kdVdZ]) *TMath::Abs(swimstate[kQoverP])); else{ state[kU] = swimstate[kU]; state[kV] = swimstate[kV]; state[kdUdZ] = swimstate[kdUdZ]; state[kdVdZ] = swimstate[kdVdZ]; state[kQoverP] = swimstate[kQoverP]; return 1; } //the swimmer works in xyz coordinates, not uvz, so need to rotate //the coordinates before putting them into a position vector Double_t dxdz = kInvSqrt2*(swimstate[kdUdZ]-swimstate[kdVdZ]); Double_t dydz = kInvSqrt2*(swimstate[kdUdZ]+swimstate[kdVdZ]); TVector3 momentum(pz*dxdz,pz*dydz,pz); SwimParticle muon(position,momentum); muon.SetCharge(charge); SwimZCondition zc(finalZ); bool done; if(idir==0) done = s.SwimForward(muon,zc); else if(idir==1) done = s.SwimBackward(muon,zc); else abort(); if(muon.GetDirection().Z()==0.0){ state[kU] = swimstate[kU]; state[kV] = swimstate[kV]; state[kdUdZ] = swimstate[kdUdZ]; state[kdVdZ] = swimstate[kdVdZ]; state[kQoverP] = swimstate[kQoverP]; return 0; } if(muon.GetMomentumModulus()==0.0){ state[kU] = swimstate[kU]; state[kV] = swimstate[kV]; state[kdUdZ] = swimstate[kdUdZ]; state[kdVdZ] = swimstate[kdVdZ]; state[kQoverP] = swimstate[kQoverP]; return 0; } state[kU] = kInvSqrt2*(muon.GetPosition().X()+muon.GetPosition().Y()); state[kV] = kInvSqrt2*(muon.GetPosition().Y()-muon.GetPosition().X()); dxdz = muon.GetDirection().X()/muon.GetDirection().Z(); dydz = muon.GetDirection().Y()/muon.GetDirection().Z(); state[kdUdZ] = kInvSqrt2*(dxdz+dydz); state[kdVdZ] = kInvSqrt2*(dydz-dxdz); state[kQoverP] = muon.GetCharge()/muon.GetMomentumModulus(); // MSG("FitTrackSR", Msg::kDebug) << "state = " << state[kU] << " " // << state[kV] << " " // << state[kdUdZ] << " " // << state[kdVdZ] << " " // << state[kQoverP] << endl; if(TMath::Abs(state[kdUdZ])>2.*fParmMaxAngle || TMath::Abs(state[kdVdZ])>2.*fParmMaxAngle || TMath::Abs(state[kQoverP])>fParmMaxQP){ MSG("FitTrackSR",Msg::kDebug) << "swimming failed, u,v,du/dz,dv,dz/qp = " << state[kU] << " " << state[kV] << " " << state[kdUdZ] << " " << state[kdVdZ] << " " << state[kQoverP] << endl; return 0; } return 1; }

void AlgFitTrackSR::SwimVertexAndEndPoints (  CandFitTrackSRHandle &  cfh, KalmanPlaneSR *  kp0, KalmanPlaneSR *  kp1, const KalmanPlaneSR *  kpu0, const KalmanPlaneSR *  kpu1, const KalmanPlaneSR *  kpv0, const KalmanPlaneSR *  kpv1, Double_t *  planepe, Int_t  plane0, Int_t  plane1, Double_t &  vtxqp, Int_t  direction )  [private]

Definition at line 3670 of file AlgFitTrackSR.cxx. References fDetector, fParmLastPlane, fVldContext, CandRecoHandle::GetEndPlane(), KalmanPlaneSR::GetFilStateValue(), TrackClusterSR::GetPlane(), UgliGeomHandle::GetScintPlnHandle(), KalmanPlaneSR::GetTrackCluster(), CandRecoHandle::GetVtxPlane(), UgliPlnHandle::GetZ0(), TrackClusterSR::GetZPos(), UgliScintPlnHandle::IsValid(), PlexPlaneId::IsValid(), kdUdZ, kdVdZ, kInvSqrt2, kU, kV, min, CandRecoHandle::SetDirCosU(), CandRecoHandle::SetDirCosV(), CandRecoHandle::SetDirCosZ(), CandRecoHandle::SetEndDirCosU(), CandRecoHandle::SetEndDirCosV(), CandRecoHandle::SetEndDirCosZ(), CandFitTrackSRHandle::SetEndExtrapolate(), CandRecoHandle::SetEndPlane(), CandTrackHandle::SetEndTrace(), CandTrackHandle::SetEndTraceZ(), CandRecoHandle::SetEndU(), CandRecoHandle::SetEndV(), CandRecoHandle::SetEndZ(), CandTrackHandle::SetU(), CandTrackHandle::SetV(), CandFitTrackSRHandle::SetVtxExtrapolate(), CandRecoHandle::SetVtxPlane(), CandTrackHandle::SetVtxTrace(), CandTrackHandle::SetVtxTraceZ(), CandRecoHandle::SetVtxU(), CandRecoHandle::SetVtxV(), CandRecoHandle::SetVtxZ(), and Swim(). Referenced by RunAlg(). { Double_t swimstate0[5]; Double_t tmpstate[5]; Double_t swimstatez0 = kp0->GetTrackCluster()->GetZPos(); Double_t swimstate1[5]; UgliGeomHandle ugh(fVldContext); for(Int_t i = 0; i<5; ++i){ tmpstate[i] = 0.; swimstate0[i] = kp0->GetFilStateValue(i,1); swimstate1[i] = kp1->GetFilStateValue(i,0); } Double_t swimstatez1 = kp1->GetTrackCluster()->GetZPos(); Int_t thisplane = kp0->GetTrackCluster()->GetPlane()-direction; Double_t thisx = kInvSqrt2*(swimstate0[kU]-swimstate0[kV]); Double_t thisy = kInvSqrt2*(swimstate0[kU]+swimstate0[kV]); Double_t thisu = swimstate0[kU]; Double_t thisv = swimstate0[kV]; Double_t thisz = swimstatez0; Double_t linx = thisx; Double_t liny = thisy; Double_t linu = thisu; Double_t linv = thisv; Double_t linz = thisz; Double_t lindx = kInvSqrt2*(swimstate0[kdUdZ]-swimstate0[kdVdZ]); Double_t lindy = kInvSqrt2*(swimstate0[kdUdZ]+swimstate0[kdVdZ]); Double_t lindu = swimstate0[kdUdZ]; Double_t lindv = swimstate0[kdVdZ]; Double_t linds = 0.; Double_t lindz = 0.; Double_t oldu = 0.; Double_t oldv = 0.; Double_t oldz = 0.; Double_t thisds = 0.; Double_t thisdz = 0.; Int_t nplaneskip = 0; Int_t swimflag = 1; Bool_t befvtx = 1; Double_t du = 0.; Double_t dv = 0.; Double_t dz = 0.; Double_t dsdz = 0.; Double_t delu = 0.; Double_t delv = 0.; Double_t delz = 0.; while(thisplane>=0 && thisplane<=fParmLastPlane && ((swimflag && TMath::Abs(thisx)<4. && TMath::Abs(thisy)<4. && TMath::Abs(thisu)<4. && TMath::Abs(thisv)<4.) || (TMath::Abs(linx)<4. && TMath::Abs(liny)<4. && TMath::Abs(linu)<4. && TMath::Abs(linv)<4.) ) ){ PlexPlaneId plexplane(fDetector,thisplane,false); if(plexplane.IsValid()){ UgliScintPlnHandle planeid = ugh.GetScintPlnHandle(plexplane); if(planeid.IsValid()){ oldu = swimstate0[kU]; oldv = swimstate0[kV]; oldz = swimstatez0; dz = (Double_t)planeid.GetZ0()-linz; if(TMath::Abs(linx)<4. && TMath::Abs(liny)<4. && TMath::Abs(linu)<4. && TMath::Abs(linv)<4.){ linz = (Double_t)planeid.GetZ0(); linx += dz*lindx; liny += dz*lindy; linu += dz*lindu; linv += dz*lindv; if(!befvtx){ du = dz*lindu; dv = dz*lindv; lindz += TMath::Abs(dz); linds += TMath::Sqrt(du*du+dv*dv+dz*dz); } }//end if inside the detector if(swimflag && TMath::Abs(thisx)<4. && TMath::Abs(thisy)<4. && TMath::Abs(thisu)<4. && TMath::Abs(thisv)<4.) { swimflag = Swim(swimstate0,tmpstate,(Double_t)swimstatez0, (Double_t)planeid.GetZ0(),1, &fVldContext); swimstatez0 = (Double_t)planeid.GetZ0(); swimstate0[kU] = tmpstate[kU]; swimstate0[kV] = tmpstate[kV]; swimstate0[kdUdZ] = tmpstate[kdUdZ]; swimstate0[kdVdZ] = tmpstate[kdVdZ]; swimstate0[kQoverP] = tmpstate[kQoverP]; if (TMath::Abs(tmpstate[kU])>4. || TMath::Abs(tmpstate[kV])>4.) swimflag = 0; thisx = kInvSqrt2*(swimstate0[kU]-swimstate0[kV]); thisy = kInvSqrt2*(swimstate0[kU]+swimstate0[kV]); thisu = swimstate0[kU]; thisv = swimstate0[kV]; thisz = swimstatez0; if(!swimflag && nplaneskip<3 && thisplane*direction>=plane0*direction){ if(!thisu) cfh.SetVtxU(thisu); if(!thisv) cfh.SetVtxV(thisv); cfh.SetVtxZ(thisz); cfh.SetVtxPlane(thisplane); dsdz = TMath::Sqrt(1.+swimstate0[kdUdZ]*swimstate0[kdUdZ] +swimstate0[kdVdZ]*swimstate0[kdVdZ]); cfh.SetDirCosU((Double_t)(direction)*swimstate0[kdUdZ]/dsdz); cfh.SetDirCosV((Double_t)(direction)*swimstate0[kdVdZ]/dsdz); cfh.SetDirCosZ((Double_t)(direction)/dsdz); cfh.SetU(thisplane,thisu); cfh.SetV(thisplane,thisv); vtxqp = swimstate0[kQoverP]; }//end if swimflag and not too many planes skipped and this plane is beyond the //vertex else{ befvtx = 0; delu = thisu-oldu; delv = thisv-oldv; delz = thisz-oldz; thisds += TMath::Sqrt(delu*delu+delv*delv+delz*delz); thisdz += TMath::Abs(delz); } //look to see if this plane was hit in the slice, if //not increment the number of skipped planes if(planepe[thisplane] > 0.) nplaneskip = 0; else ++nplaneskip; }//end if in the detector }//end if planeid is valid }//end if plexplane is valid thisplane -= direction; }//end loop over planes cfh.SetVtxTrace(min(thisds,linds)); cfh.SetVtxTraceZ(min(thisdz,lindz)); if(kpu0) cfh.SetVtxExtrapolate(TMath::Abs(kpu0->GetTrackCluster()->GetPlane()-cfh.GetVtxPlane()), PlaneView::kU); else cfh.SetVtxExtrapolate(-1,PlaneView::kU); if(kpv0) cfh.SetVtxExtrapolate(TMath::Abs(kpv0->GetTrackCluster()->GetPlane()-cfh.GetVtxPlane()), PlaneView::kV); else cfh.SetVtxExtrapolate(-1,PlaneView::kV); thisplane = kp1->GetTrackCluster()->GetPlane()+direction; thisx = kInvSqrt2*(swimstate1[kU]-swimstate1[kV]); thisy = kInvSqrt2*(swimstate1[kU]+swimstate1[kV]); thisu = swimstate1[kU]; thisv = swimstate1[kV]; thisz = swimstatez1; linx = thisx; liny = thisy; linu = thisu; linv = thisv; linz = thisz; lindx = kInvSqrt2*(swimstate1[kdUdZ]-swimstate1[kdVdZ]); lindy = kInvSqrt2*(swimstate1[kdUdZ]+swimstate1[kdVdZ]); lindu = swimstate1[kdUdZ]; lindv = swimstate1[kdVdZ]; linds = 0.; lindz = 0.; oldu = 0.; oldv = 0.; oldz = 0.; thisds = 0.; thisdz = 0.; nplaneskip = 0; swimflag = 1; befvtx = 1; while(thisplane>=0 && thisplane<=fParmLastPlane && ((swimflag && TMath::Abs(thisx)<4. && TMath::Abs(thisy)<4. && TMath::Abs(thisu)<4. && TMath::Abs(thisv)<4.) || (TMath::Abs(linx)<4. && TMath::Abs(liny)<4. && TMath::Abs(linu)<4. && TMath::Abs(linv)<4.) ) ){ PlexPlaneId plexplane(fDetector,thisplane,false); if(plexplane.IsValid()){ UgliScintPlnHandle planeid = ugh.GetScintPlnHandle(plexplane); if(planeid.IsValid()){ oldu = swimstate1[kU]; oldv = swimstate1[kV]; oldz = swimstatez1; dz = (Double_t)planeid.GetZ0()-linz; if(TMath::Abs(linx)<4. && TMath::Abs(liny)<4. && TMath::Abs(linu)<4. && TMath::Abs(linv)<4.){ linz = (Double_t)planeid.GetZ0(); linx += dz*lindx; liny += dz*lindy; linu += dz*lindu; linv += dz*lindv; if (!befvtx) { du = dz*lindu; dv = dz*lindv; lindz += TMath::Abs(dz); linds += TMath::Sqrt(du*du+dv*dv+dz*dz); } }//end if in the detector if(swimflag && TMath::Abs(thisx)<4. && TMath::Abs(thisy)<4. && TMath::Abs(thisu)<4. && TMath::Abs(thisv)<4.){ swimflag = Swim(swimstate1,tmpstate,(Double_t)swimstatez1, (Double_t)planeid.GetZ0(),0, &fVldContext); swimstatez1 = (Double_t)planeid.GetZ0(); swimstate1[kU] = tmpstate[kU]; swimstate1[kV] = tmpstate[kV]; swimstate1[kdUdZ] = tmpstate[kdUdZ]; swimstate1[kdVdZ] = tmpstate[kdVdZ]; swimstate1[kQoverP] = tmpstate[kQoverP]; if (TMath::Abs(tmpstate[kU])>4. || TMath::Abs(tmpstate[kV])>4.) swimflag = 0; thisx = kInvSqrt2*(swimstate1[kU]-swimstate1[kV]); thisy = kInvSqrt2*(swimstate1[kU]+swimstate1[kV]); thisu = swimstate1[kU]; thisv = swimstate1[kV]; thisz = swimstatez1; if(!swimflag && nplaneskip<3 && thisplane*direction<=plane1*direction){ cfh.SetEndU(thisu); cfh.SetEndV(thisv); cfh.SetEndZ(thisz); cfh.SetEndPlane(thisplane); dsdz = TMath::Sqrt(1.+swimstate1[kdUdZ]*swimstate1[kdUdZ] +swimstate1[kdVdZ]*swimstate1[kdVdZ]); cfh.SetEndDirCosU((Double_t)(direction)*swimstate1[kdUdZ]/dsdz); cfh.SetEndDirCosV((Double_t)(direction)*swimstate1[kdVdZ]/dsdz); cfh.SetEndDirCosZ((Double_t)(direction)/dsdz); cfh.SetU(thisplane,thisu); cfh.SetV(thisplane,thisv); } else { befvtx = 0; delu = thisu-oldu; delv = thisv-oldv; delz = thisz-oldz; thisds += TMath::Sqrt(delu*delu+delv*delv+delz*delz); thisdz += TMath::Abs(delz); } if(planepe[thisplane] > 0.) nplaneskip = 0; else ++nplaneskip; }//end if in the detector and swimflag }//end if planeid is valid }//end if plexplane is valid thisplane += direction; }//end loop over planes to find the end cfh.SetEndTrace(min(thisds,linds)); cfh.SetEndTraceZ(min(thisdz,lindz)); if(kpu1) cfh.SetEndExtrapolate(TMath::Abs(kpu1->GetTrackCluster()->GetPlane()-cfh.GetEndPlane()), PlaneView::kU); else cfh.SetEndExtrapolate(-1,PlaneView::kU); if(kpv1) cfh.SetEndExtrapolate(TMath::Abs(kpv1->GetTrackCluster()->GetPlane()-cfh.GetEndPlane()), PlaneView::kV); else cfh.SetEndExtrapolate(-1,PlaneView::kV); return; }

void AlgFitTrackSR::Trace (  const char *  c  )  const [virtual]

Reimplemented from AlgBase. Definition at line 1713 of file AlgFitTrackSR.cxx. { }

---

Member Data Documentation

Double_t AlgFitTrackSR::fCov[5][5] [private]

Definition at line 84 of file AlgFitTrackSR.h. Referenced by AddToFit(), InitializeFit(), and ReverseFit().

Detector::Detector_t AlgFitTrackSR::fDetector [private]

Definition at line 42 of file AlgFitTrackSR.h. Referenced by FindNumSkippedPlanes(), RunAlg(), and SwimVertexAndEndPoints().

Double_t AlgFitTrackSR::fParmCovarianceScale [private]

Definition at line 81 of file AlgFitTrackSR.h. Referenced by RunAlg().

Double_t AlgFitTrackSR::fParmdedx [private]

Definition at line 69 of file AlgFitTrackSR.h. Referenced by CalculateNoise(), and InitKalmanFitParameters().

Double_t AlgFitTrackSR::fParmDeltaChi2 [private]

Definition at line