AlgFitTrackSR Class Reference

#include <AlgFitTrackSR.h>

Inheritance diagram for AlgFitTrackSR:
AlgBase AlgReco AlgTrack

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]

Detailed Description

Definition at line 30 of file AlgFitTrackSR.h.


Constructor & Destructor Documentation

AlgFitTrackSR::AlgFitTrackSR (  ) 

Definition at line 73 of file AlgFitTrackSR.cxx.

00074 {
00075 }

AlgFitTrackSR::~AlgFitTrackSR (  )  [virtual]

Definition at line 78 of file AlgFitTrackSR.cxx.

00079 {
00080 }


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::kDebug, PlaneView::kU, PlaneView::kV, MSG, and tc.

Referenced by DoKalmanFit(), and IterateKalmanFit().

02364 {
02365   MSG("FitTrackSR", Msg::kDebug) << "in AddClustersToFit" << endl;
02366   
02367   TrackClusterSR *tc = 0;
02368   Int_t istatus = 0;
02369 
02370   Bool_t dir = kIterForward;
02371   if(direction!=1) dir=kIterBackward;
02372   TIter planeItr(&planeClusterList,dir);
02373 
02374   TObjArray * plnarray;
02375   while( (plnarray = (TObjArray * )planeItr.Next()) ){
02376     TIter clusterItr(plnarray);
02377     while( (tc = (TrackClusterSR *)clusterItr.Next()) && istatus>=0 ){
02378     
02379       MSG("FitTrackSR", Msg::kDebug) << "current cluster on plane " 
02380                                      << tc->GetPlane() << " is valid "
02381                                      << (Int_t)tc->IsValid() << endl;
02382       if(tc->IsValid()){
02383         MSG("FitTrackSR", Msg::kDebug) << "adding cluster on plane "
02384                                        << tc->GetPlane() << " " 
02385                                        << tc->GetMinStrip() << "-" << tc->GetMaxStrip()
02386                                        << " to fit" << endl;
02387         
02388         istatus = AddToFit(cfh,tc, iterate);
02389         
02390         if(istatus>=0){
02391           if(tc->GetPlaneView()==PlaneView::kU) ++nu;
02392           else if(tc->GetPlaneView()==PlaneView::kV) ++nv;
02393         }
02394       }
02395     }
02396   }//end loop over clusters
02397   return istatus;
02398 }

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(), Msg::kDebug, kdUdZ, kdVdZ, MSG, CandFitTrackSRHandle::RemoveKalmanPlane(), and CandFitTrackHandle::SetNSwimFail().

Referenced by FindDownstreamPlanes().

03410 {  
03411 
03412   MSG("FitTrackSR", Msg::kDebug) << "in AddForwardBestKPToFit" << endl;
03413   Int_t nadd = 0;
03414     
03415   MSG("FitTrackSR",Msg::kDebug) << "considering best tc " 
03416                                 << bestkp->GetTrackCluster()->GetPlane() 
03417                                 << " " << bestkp->GetTrackCluster()->GetMinStrip() 
03418                                 << "-" << bestkp->GetTrackCluster()->GetMaxStrip() 
03419                                 << ", chi2 = " << bestkp->GetPreChi2(0) << "/" 
03420                                 << fParmMaxLocalPreChi2 << " " 
03421                                 << bestkp->GetPreCovarianceMatrixValue(0,kdUdZ,kdUdZ) 
03422                                 << "/"
03423                                 << fParmMaxAngleCovariance << " " 
03424                                 << bestkp->GetPreCovarianceMatrixValue(0,kdVdZ,kdVdZ) 
03425                                 << "/"
03426                                 << fParmMaxAngleCovariance << endl;
03427                         
03428   if(bestkp->GetPreChi2(0)<fParmMaxLocalPreChi2 
03429      && bestkp->GetPreCovarianceMatrixValue(0,kdUdZ,kdUdZ)<fParmMaxAngleCovariance 
03430      && bestkp->GetPreCovarianceMatrixValue(0,kdVdZ,kdVdZ)<fParmMaxAngleCovariance){
03431     
03432     KalmanPlaneSR *kpnewer = cfh.AddKalmanPlaneAt(bestkp,cfh.GetKalmanLast()+1);
03433     Int_t istatus2 = AddToFit(cfh,kpnewer,-1,0);
03434     if (istatus2>=0) {
03435       cfh.SetNSwimFail(cfh.GetNSwimFail()+nswimfail);
03436       ++nadd;
03437     } 
03438     else {
03439  MSG("FitTrackSR",Msg::kDebug) << "AddToFit failed, do not add to fit" 
03440                                     << endl;
03441 
03442       // need to remove from KalmanPlaneList
03443       cfh.RemoveKalmanPlane(-1,0);
03444       if(bestkp){
03445         delete bestkp;
03446         bestkp = 0;
03447       }
03448     }
03449   }
03450   else if(bestkp){
03451     delete bestkp;
03452     bestkp = 0;
03453   }
03454 
03455   nswimfail=0;
03456 
03457   return nadd;
03458 }

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(), Msg::kDebug, kdUdZ, kdVdZ, MSG, CandFitTrackSRHandle::SetCurrentKalmanPlane(), and CandFitTrackHandle::SetNSwimFail().

Referenced by FindUpstreamPlanes().

03468 {
03469   MSG("FitTrackSR", Msg::kDebug) << "in AddReverseBestKPToFit" << endl;
03470   Int_t nadd = 0;
03471  MSG("FitTrackSR",Msg::kDebug) << "considering best tc " 
03472                                 << bestkp->GetTrackCluster()->GetPlane() 
03473                                 << " " << bestkp->GetTrackCluster()->GetMinStrip() 
03474                                 << "-" << bestkp->GetTrackCluster()->GetMaxStrip() 
03475                                 << ", chi2 = " << bestkp->GetPreChi2(1) << "/" 
03476                                 << fParmMaxLocalPreChi2 << " planediff = " 
03477                                 << bestkp->GetTrackCluster()->GetPlane()-oldkp->GetTrackCluster()->GetPlane() 
03478                                 << " " 
03479                                 << bestkp->GetPreCovarianceMatrixValue(1,kdUdZ,kdUdZ)
03480                                 << " " 
03481                                 << bestkp->GetPreCovarianceMatrixValue(1,kdVdZ,kdVdZ) 
03482                                 << endl;  
03483   if(bestkp->GetPreChi2(1)<fParmMaxLocalPreChi2 
03484      && bestkp->GetPreCovarianceMatrixValue(1,kdUdZ,kdUdZ)<fParmMaxAngleCovariance 
03485      && bestkp->GetPreCovarianceMatrixValue(1,kdVdZ,kdVdZ)<fParmMaxAngleCovariance){
03486     
03487     // check if forward fit is also good
03488     // first find the plane just before bestkp
03489     KalmanPlaneSR *foundkp = 0;
03490     KalmanPlaneSR *befkp = 0;
03491     
03492     //look for the kp before the best kp
03493     for(Int_t i=cfh.GetKalmanLast(); i>=0 && !foundkp; --i){
03494       befkp = cfh.GetKalmanPlane(i);
03495       if(befkp->GetTrackCluster()->GetPlane()*direction
03496          <bestkp->GetTrackCluster()->GetPlane()*direction)
03497         foundkp = befkp;
03498                 
03499     }
03500               
03501     Double_t forwardchi2 = 0.;
03502     Int_t befswimfail = 0;
03503     if(foundkp){
03504       befswimfail = FitFrom(foundkp,bestkp,0,-1);
03505       forwardchi2 = bestkp->GetPreChi2(0);
03506   MSG("FitTrackSR",Msg::kDebug) << "checking forward fit, previous " 
03507                                    << "cluster plane = " 
03508                                    << foundkp->GetTrackCluster()->GetPlane() 
03509                                    << " min/max strip = " 
03510                                    << foundkp->GetTrackCluster()->GetMinStrip() 
03511                                    << "/" 
03512                                    << foundkp->GetTrackCluster()->GetMaxStrip() 
03513                                    << " nswimfail = " << befswimfail 
03514                                    << " prechi2 = " << forwardchi2 << endl;
03515 
03516     }
03517     if(!befswimfail && forwardchi2<fParmMaxLocalPreChi2){
03518      MSG("FitTrackSR",Msg::kDebug) << "ADDING TO FIT" << endl;   
03519       ++nadd;
03520       oldkp = bestkp;
03521       newkplist->AddLast(oldkp);
03522       cfh.SetCurrentKalmanPlane(oldkp);
03523       cfh.SetNSwimFail(cfh.GetNSwimFail()+nswimfail);
03524     } 
03525     else delete bestkp;
03526     
03527   } 
03528   else delete bestkp;
03529   
03530   bestkp = 0;
03531   nswimfail = 0;
03532   
03533   return nadd;
03534 }

Int_t AlgFitTrackSR::AddToFit ( CandFitTrackSRHandle cfh,
KalmanPlaneSR thiskp,
Int_t  iterate,
Bool_t  docheck 
) [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::kDebug, MSG, pFor, CandFitTrackSRHandle::SetCurrentKalmanPlane(), CandFitTrackSRHandle::SetEnd(), and CandFitTrackHandle::SetNSwimFail().

00650 {
00651 
00652  
00653 
00654 // < 0 = failure
00655 // 0 = success
00656 
00657 // this method assumes that the plane at thiskp is not already in the kalmanplanelist
00658 // note also that this method takes ownership of thiskp
00659 // if docheck is false, user is responsible for insertion into kalmanplanelist
00660 // before call to AddToFit
00661 
00662   TObjArray * KalmanPlaneList=cfh.GetPlaneList();
00663  
00664   Bool_t found(0);
00665   Int_t ifound=0;
00666   Int_t ibef=-1;
00667   Int_t ilast;
00668 
00669   KalmanPlaneSR *kpnew = 0;
00670   KalmanPlaneSR *kp = 0;
00671   KalmanPlaneSR *kpold = 0;
00672   
00673   if(docheck){
00674     
00675     for(Int_t i=0; i<=KalmanPlaneList->GetLast(); ++i){
00676       
00677       // assume only one track cluster per plane in track
00678       kp = dynamic_cast<KalmanPlaneSR*>(KalmanPlaneList->At(i));
00679       
00680       if(kp->GetTrackCluster()->GetPlane() == thiskp->GetTrackCluster()->GetPlane()) {
00681         found = 1;
00682         ifound = i;
00683         kpnew = kp;
00684 
00685         // kalmanplane already in list, need to delete thiskp
00686         delete thiskp;
00687       } 
00688       else if((kp->GetZDir())*(kp->GetTrackCluster()->GetPlane())<(kp->GetZDir())*(thiskp->GetTrackCluster()->GetPlane())){
00689         ibef = i;
00690       }
00691     }//end loop over kalman planes
00692 
00693     //didnt find a plane to bound this one, so add it at the last place in the list
00694     if(!found){
00695       kpnew = thiskp;
00696 
00697       ilast = KalmanPlaneList->GetLast();
00698       
00699       for(int i=ilast; i>ibef; --i){
00700         kp = dynamic_cast<KalmanPlaneSR*>(KalmanPlaneList->At(i));
00701         KalmanPlaneList->AddAtAndExpand(kp,i+1);
00702       }
00703 
00704       KalmanPlaneList->AddAtAndExpand(kpnew,ibef+1);
00705     }//end didnt find a plane
00706   }//end if docheck 
00707   else kpnew = thiskp;
00708   
00709   cfh.SetCurrentKalmanPlane(kpnew);
00710 
00711   //if this is the first kp added to the track, initialize the fit
00712   if(KalmanPlaneList->GetLast()==0 
00713      || (docheck && ((found && ifound==0) || ibef<0))) InitializeFit(cfh,0,iterate);
00714   else{
00715   
00716     //get the plane of the kp immediately before the one to add to the 
00717     //fit
00718     kpold = dynamic_cast<KalmanPlaneSR*>
00719       (KalmanPlaneList->At(KalmanPlaneList->IndexOf(kpnew)-1));
00720     
00721     //fit from the previously fit plane to the one to add to the fit
00722     pFor++;
00723     Int_t nswimfail = FitFrom(kpold,kpnew,0,iterate);
00724     if(nswimfail<0){
00725       return nswimfail;
00726     }
00727 
00728     cfh.SetNSwimFail(cfh.GetNSwimFail()+nswimfail);
00729     for(int i1=0; i1<5; ++i1){
00730       for(int i2=0; i2<5; ++i2){
00731         fCov[i1][i2] = kpnew->GetFilCovarianceMatrixValue(0,i1,i2);
00732       }
00733     }//end loop to fill fCov covariance matrix
00734     MSG("FitTrackSR",Msg::kDebug) << "fitting forward, prechi2 = " 
00735                                   << kpnew->GetPreChi2(0) << " filchi2 = " 
00736                                   << kpnew->GetFilChi2(0) << endl;
00737 
00738 
00739     cfh.SetEnd(kpnew);
00740   }//end else
00741 
00742   return 0;
00743 }

Int_t AlgFitTrackSR::AddToFit ( CandFitTrackSRHandle cfh,
TrackClusterSR tc,
Int_t  iterate 
) [private]

Definition at line 628 of file AlgFitTrackSR.cxx.

References fParmIsCosmic, TrackClusterSR::GetPlane(), CandTrackHandle::GetRange(), CandRecoHandle::GetTimeSlope(), KalmanPlaneSR::GetTrackCluster(), KalmanPlaneSR::SetRange(), and KalmanPlaneSR::SetZDir().

Referenced by AddClustersToFit(), and AddForwardBestKPToFit().

00631 {
00632 
00633   
00634 // < 0 = failure
00635 // 0 = success
00636 
00637   KalmanPlaneSR *kpnew = new KalmanPlaneSR(tc);
00638   kpnew->SetRange(cfh.GetRange(kpnew->GetTrackCluster()->GetPlane()));
00639   
00640   if(cfh.GetTimeSlope()>0. || !fParmIsCosmic) kpnew->SetZDir(1);
00641   else kpnew->SetZDir(-1);
00642   
00643   return AddToFit(cfh, kpnew, iterate,1);
00644 }

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

Definition at line 1514 of file AlgFitTrackSR.cxx.

References KalmanPlaneSR::GetFilCovarianceMatrixValue(), KalmanPlaneSR::GetFilStateValue(), TrackClusterSR::GetPlaneView(), TrackClusterSR::GetTPos(), KalmanPlaneSR::GetTrackCluster(), Msg::kDebug, PlaneView::kV, MSG, and KalmanPlaneSR::SetFilChi2().

Referenced by Filter().

01516 {
01517 
01518   Int_t index = 0;
01519 
01520   //default U view, change if V view
01521   if(currentkp->GetTrackCluster()->GetPlaneView()==PlaneView::kV) index = 1;
01522 
01523   Double_t dtpos = 0.;
01524 
01525   dtpos = currentkp->GetFilStateValue(index,idir)-currentkp->GetTrackCluster()->GetTPos();
01526   
01527   Double_t filChi2 = dtpos*dtpos/(poserr*poserr);
01528   
01529   currentkp->SetFilChi2(filChi2, idir);
01530   MSG("FitTrackSR",Msg::kDebug) << "filchi2 = " << filChi2 << "  dtpos = " 
01531                                 << dtpos << "  tposerr2 = " << poserr*poserr 
01532                                 << "  cov = " 
01533                                 << currentkp->GetFilCovarianceMatrixValue(idir,index,index) 
01534                                 << endl;
01535 
01536 }

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, PlaneView::kV, and KalmanPlaneSR::SetFilCovarianceMatrixValue().

Referenced by Filter().

01541 {
01542   //do the final step of the Kalman filter
01543   //C_{k} = (I - K_{k}H_{k})C_{k}^{k-1}
01544   
01545   Int_t index = 0;
01546 
01547   //default U view, change if V view
01548   if(currentkp->GetTrackCluster()->GetPlaneView()==PlaneView::kV) index = 1;
01549   
01550   TMatrixD preCov(5,5);
01551   for(Int_t i = 0; i<5; ++i){
01552     for(Int_t j = 0; j<5; ++j){
01553       preCov(i,j) = currentkp->GetPreCovarianceMatrixValue(idir,i,j);
01554     }
01555   }
01556 
01557   TMatrixD filCovariance(5,5);
01558 
01559   //set the filtered covariance matrix to the identity matrix
01560   filCovariance.UnitMatrix();
01561 
01562   //subtract K_{k}H_{k} from the identity
01563   filCovariance(kU,index) -= currentkp->GetGainValue(kU,idir);
01564   filCovariance(kV,index) -= currentkp->GetGainValue(kV,idir);
01565   filCovariance(kdUdZ,index) -= currentkp->GetGainValue(kdUdZ,idir);
01566   filCovariance(kdVdZ,index) -= currentkp->GetGainValue(kdVdZ,idir);
01567   filCovariance(kQoverP,index) -= currentkp->GetGainValue(kQoverP,idir);
01568   
01569   //multiply the difference by the pre filtered covariance matrix
01570   filCovariance *= preCov;
01571 
01572   // diagonal elements should be positive
01573   Double_t covlim = 1.e-8;
01574   if(filCovariance(kU,kU)<covlim) filCovariance(kU,kU) = covlim;
01575   if(filCovariance(kV,kV)<covlim) filCovariance(kV,kV) = covlim;
01576   if(filCovariance(kdUdZ,kdUdZ)<covlim)
01577     filCovariance(kdUdZ,kdUdZ) = covlim;
01578   if(filCovariance(kdVdZ,kdVdZ)<covlim)
01579     filCovariance(kdVdZ,kdVdZ) = covlim;
01580   if(filCovariance(kQoverP,kQoverP)<covlim) filCovariance(kQoverP,kQoverP) = covlim;
01581   
01582   for(Int_t i = 0; i<5; ++i){
01583     for(Int_t j = 0; j<5; ++j){
01584       currentkp->SetFilCovarianceMatrixValue(idir,i,j, filCovariance(i,j));
01585     }
01586   }
01587 }

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

Definition at line 1449 of file AlgFitTrackSR.cxx.

References fParmMaxAngle, fParmMaxQP, fParmMaxQPFrac, fParmQPRangeCheck, KalmanPlaneSR::GetFilStateValue(), KalmanPlaneSR::GetGainValue(), TrackClusterSR::GetPlaneView(), KalmanPlaneSR::GetPreStateValue(), KalmanPlaneSR::GetRange(), TrackClusterSR::GetTPos(), KalmanPlaneSR::GetTrackCluster(), kdUdZ, kdVdZ, kQoverP, kU, kV, PlaneView::kV, and KalmanPlaneSR::SetFilStateValue().

Referenced by Filter().

01451 {
01452  
01453 
01454   Int_t index = 0;
01455 
01456   //default U view, change if V view
01457   if(currentkp->GetTrackCluster()->GetPlaneView()==PlaneView::kV) index = 1;
01458 
01459   //the following is the m_{k} - H_{k}F_{k-1}x_{k-1}
01460   //bit of the Kalman filter equations
01461   //remember that 
01462   //m_{k} = measurement, 
01463   //H_{k} = [1,0,0,0,0] for U view; [0,1,0,0,0] for V view
01464   //the definition of H_{k} makes the H_{k}F_{k-1}x_{k-1} product
01465   //a scalar and picks out only certain elements of F_{k-1}x_{k-1}
01466   //also, F_{k-1}x_{k-1} is found by the swimmer and put in the 
01467   //pre filtered state vector of the plane to add to the fit
01468 
01469   Double_t xfac = (currentkp->GetTrackCluster()->GetTPos()
01470                    - currentkp->GetPreStateValue(index,idir));
01471 
01472   currentkp->SetFilStateValue(kU,idir, currentkp->GetPreStateValue(kU,idir)
01473                               +currentkp->GetGainValue(kU,idir)*xfac);
01474   currentkp->SetFilStateValue(kV,idir, currentkp->GetPreStateValue(kV,idir)
01475                               +currentkp->GetGainValue(kV,idir)*xfac);
01476   currentkp->SetFilStateValue(kdUdZ,idir, currentkp->GetPreStateValue(kdUdZ,idir)
01477                               +currentkp->GetGainValue(kdUdZ,idir)*xfac);
01478   currentkp->SetFilStateValue(kdVdZ,idir, currentkp->GetPreStateValue(kdVdZ,idir)
01479                               +currentkp->GetGainValue(kdVdZ,idir)*xfac);
01480   currentkp->SetFilStateValue(kQoverP,idir, currentkp->GetPreStateValue(kQoverP,idir)
01481                               +currentkp->GetGainValue(kQoverP,idir)*xfac);
01482 
01483   if(TMath::Abs(currentkp->GetFilStateValue(kdUdZ,idir))>fParmMaxAngle){
01484     currentkp->SetFilStateValue(kdUdZ,idir,(currentkp->GetFilStateValue(kdUdZ,idir)>0 
01485                                             ? fParmMaxAngle : -fParmMaxAngle));
01486   }
01487 
01488   if(TMath::Abs(currentkp->GetFilStateValue(kdVdZ,idir))>fParmMaxAngle){
01489     currentkp->SetFilStateValue(kdVdZ,idir,(currentkp->GetFilStateValue(kdVdZ,idir)>0 
01490                                             ? fParmMaxAngle : -fParmMaxAngle));
01491   }
01492 
01493   Double_t maxqp = fParmMaxQP;
01494 
01495   //use a value of about 2 MeV/(g/cm^2) for the energy loss in the material
01496   //you would divide the max qp frac allowed by that value and then by the range to
01497   //get the max qp.  but multiplication is a faster operation, so multiply by
01498   //500 instead.
01499   if ( fParmQPRangeCheck ) {
01500       if(currentkp->GetRange()>0. && fParmMaxQPFrac*500./currentkp->GetRange()<maxqp) {
01501           maxqp = fParmMaxQPFrac*500./currentkp->GetRange();
01502       }
01503   }
01504 
01505   if(TMath::Abs(currentkp->GetFilStateValue(kQoverP,idir))>maxqp){
01506     currentkp->SetFilStateValue(kQoverP,idir,(currentkp->GetFilStateValue(kQoverP,idir)>0 
01507                                               ? maxqp : -maxqp));
01508   }
01509 
01510 
01511 }

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

Definition at line 1405 of file AlgFitTrackSR.cxx.

References TrackClusterSR::GetPlaneView(), KalmanPlaneSR::GetPreCovarianceMatrixValue(), KalmanPlaneSR::GetTrackCluster(), kdUdZ, kdVdZ, kQoverP, kU, kV, PlaneView::kV, and KalmanPlaneSR::SetGainValue().

Referenced by Filter().

01407 {
01408 
01409   //the following is the 
01410   //K_{k} = C_{k}^{k-1}H_{k}^{T}(V_{k}+H_{k}C_{k}^{k-1}H_{k}^{T})^{-1}
01411   //step of the filtering equations
01412   //remember that 
01413   //V_{k} = measurement error, 
01414   //H_{k} = [1,0,0,0,0] for U view; [0,1,0,0,0] for V view
01415   //the definition of H_{k} makes the H_{k}C_{k}^{k-1}H_{k}^{T} product
01416   //a scalar and picks out only certain elements of C_{k}^{k-1}
01417 
01418  
01419   Int_t index = 0;
01420 
01421   //default U view, change if V view
01422   if(currentkp->GetTrackCluster()->GetPlaneView()==PlaneView::kV) index = 1;
01423 
01424 
01425   currentkp->SetGainValue(kU,idir,
01426                           (currentkp->GetPreCovarianceMatrixValue(idir,kU,index)/
01427                            (currentkp->GetPreCovarianceMatrixValue(idir,index,index)
01428                             +poserr*poserr)));
01429   currentkp->SetGainValue(kV,idir,
01430                           (currentkp->GetPreCovarianceMatrixValue(idir,kV,index)/
01431                            (currentkp->GetPreCovarianceMatrixValue(idir,index,index)
01432                             +poserr*poserr)));
01433   currentkp->SetGainValue(kdUdZ,idir,
01434                           (currentkp->GetPreCovarianceMatrixValue(idir,kdUdZ,index)/
01435                            (currentkp->GetPreCovarianceMatrixValue(idir,index,index)
01436                             +poserr*poserr)));
01437   currentkp->SetGainValue(kdVdZ,idir,
01438                           (currentkp->GetPreCovarianceMatrixValue(idir,kdVdZ,index)/
01439                            (currentkp->GetPreCovarianceMatrixValue(idir,index,index)
01440                             +poserr*poserr)));
01441   currentkp->SetGainValue(kQoverP,idir,
01442                           (currentkp->GetPreCovarianceMatrixValue(idir,kQoverP,index)/
01443                            (currentkp->GetPreCovarianceMatrixValue(idir,index,index)
01444                             +poserr*poserr)));
01445   
01446 }

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(), TrackClusterSR::GetZPos(), SwimGeo::Instance(), Msg::kDebug, kdUdZ, kdVdZ, kInvSqrt2, kQoverP, kU, kV, Munits::m, MSG, and KalmanPlaneSR::SetNoiseMatrixValue().

Referenced by Predict().

01082 {
01083 
01084 
01085   MSG("FitTrackSR", Msg::kDebug) << "In CalculateNoise" << endl;
01086 
01087   //find noise from the previous kp to the kp to add to the fit
01088   KalmanPlaneSR *thiskp = prevkp;
01089   KalmanPlaneSR *newkp = currentkp;
01090   Int_t thisdir = idir;
01091 
01092   if(befaft){
01093     thiskp = currentkp;
01094     newkp = prevkp;
01095     thisdir = 1-thisdir;
01096   }
01097 
01098   Int_t dplane = (currentkp->GetTrackCluster()->GetPlane()
01099                   - prevkp->GetTrackCluster()->GetPlane());
01100   Int_t izdir = (dplane>0 ? 0 : 1);
01101 
01102   // izdir gives the direction of the swim (0 = forward, 1 = backward)
01103   // relative to z
01104   // this is different from thisdir/idir which gives direction of swim
01105   // relative to time
01106   // idir=0 means forward in time, idir=1 means backward in time
01107   // a particle could be going forward in z and forward in time or
01108   // any of the 4 combinations (for cosmics and atmospherics)
01109   Double_t dudz = thiskp->GetFilStateValue(kdUdZ,idir);
01110   Double_t dvdz = thiskp->GetFilStateValue(kdVdZ,idir);
01111   Double_t dsdz = sqrt(1.+dudz*dudz+dvdz*dvdz);
01112   Double_t qp = thiskp->GetFilStateValue(kQoverP,idir);
01113 
01114   // when qp is too small, errors go to zero, so place lower bound
01115   // 0.01 (= 100 GeV/c) is reasonable
01116   if(TMath::Abs(qp)<0.01) qp = (qp>0 ? 0.01 : -0.01);
01117   
01118   Double_t locradlen = TMath::Abs(dsdz*(Double_t)(dplane)*fParmPlnRadLen);
01119 
01120   // we should really use an average of 1/p in calculating sigmams2
01121   Double_t sigmams2 = (0.0136*TMath::Abs(qp)*TMath::Sqrt(locradlen)
01122                        *(1.+0.038*TMath::Log(locradlen)));
01123   sigmams2 *= sigmams2;
01124 
01125   Double_t p3p3 = 0.5*sigmams2*(1.+dudz*dudz)*dsdz*dsdz;
01126   Double_t p3p4 = 0.5*sigmams2*dudz*dvdz*dsdz*dsdz;
01127   Double_t p4p4 = 0.5*sigmams2*(1.+dvdz*dvdz)*dsdz*dsdz;
01128 
01129 //  SwimPlaneInterfaceListSR *spil = 
01130 //    SwimPlaneInterfaceListSR::Instance(const_cast<VldContext*>(GetVldContext()));
01131 
01132   SwimGeo *spil = SwimGeo::Instance(*(const_cast<VldContext*>(currentkp->GetVldContext())));
01133 
01134   Double_t z1 = spil->GetSteel(prevkp->GetTrackCluster()->GetZPos(),izdir)->GetZ();
01135   Double_t z2 = spil->GetSteel(z1,izdir)->GetZ();
01136 
01137   // treat iron as discrete scatter; in future can make this more rigorous
01138   Double_t dzscatter = TMath::Abs(currentkp->GetTrackCluster()->GetZPos()-0.5*(z1+z2));
01139   Double_t dzscatter2 = dzscatter*dzscatter;
01140   Double_t eloss = fParmdedx*TMath::Abs((Double_t)(dplane));
01141 
01142   VldContext vldc = *(prevkp->GetVldContext());
01143   UgliGeomHandle ugh(vldc);
01144   BField bf(vldc);
01145   TVector3 uvz;
01146   uvz(0) = thiskp->GetFilStateValue(kU,idir);
01147   uvz(1) = thiskp->GetFilStateValue(kV,idir);
01148   // rwh: KalmanPlane's TrackCluster z is not necessarily inside the steel
01149   //uvz(2) = thiskp->GetTrackCluster()->GetZPos();
01150   // translate an semi-arbitrary z to z0 of the nearest steel plane
01151   uvz(2) = ugh.GetNearestSteelPlnHandle(thiskp->GetTrackCluster()->GetZPos()).GetZ0();
01152   TVector3 xyz(kInvSqrt2*(uvz[0]-uvz[1]),kInvSqrt2*(uvz[0]+uvz[1]),uvz[2]);
01153   TVector3 bxyz = bf.GetBField(xyz);
01154 
01155   // rotate bxyz into buvz
01156   // bxyz in Tesla, convert into GeV/c/m
01157   // warning: numbers are hard coded here
01158   TVector3 buvz(kInvSqrt2*(bxyz[1]+bxyz[0]),kInvSqrt2*(bxyz[1]-bxyz[0]),bxyz[2]);
01159   
01160   buvz[0] *= 0.3;
01161   buvz[1] *= 0.3;
01162   buvz[2] *= 0.3;
01163 
01164   // we assume the uncertainty in energy loss to follow a Landau distribution,
01165   // and use the FWHM from Ahlen (RMP, Vol. 52, p. 143, 1980)
01166   Double_t sigmaeloss2 = 0.25*eloss*dsdz*qp*qp;
01167   sigmaeloss2 *= sigmaeloss2;
01168 
01169   currentkp->SetNoiseMatrixValue(idir,kU,kU,dzscatter2*p3p3);
01170   currentkp->SetNoiseMatrixValue(idir,kU,kV,dzscatter2*p3p4);
01171   currentkp->SetNoiseMatrixValue(idir,kU,kdUdZ,dzscatter*p3p3);
01172   currentkp->SetNoiseMatrixValue(idir,kU,kdVdZ,dzscatter*p3p4);
01173   currentkp->SetNoiseMatrixValue(idir,kU,kQoverP,
01174                                  (dzscatter*sigmaeloss2*TMath::Abs((Double_t)(dplane))
01175                                   *buvz(1)*0.0254*Munits::m*dsdz*(1.+dudz*dudz)));
01176   currentkp->SetNoiseMatrixValue(idir,kV,kU,dzscatter2*p3p4);
01177   currentkp->SetNoiseMatrixValue(idir,kV,kV,dzscatter2*p4p4);
01178   currentkp->SetNoiseMatrixValue(idir,kV,kdUdZ,dzscatter*p3p4);
01179   currentkp->SetNoiseMatrixValue(idir,kV,kdVdZ,dzscatter*p4p4);
01180   currentkp->SetNoiseMatrixValue(idir,kV,4,
01181                                  (dzscatter*sigmaeloss2*TMath::Abs((Double_t)(dplane))
01182                                   *buvz(0)*0.0254*Munits::m*dsdz*(1.+dvdz*dvdz)));
01183   currentkp->SetNoiseMatrixValue(idir,kdUdZ,kU,dzscatter*p3p3);
01184   currentkp->SetNoiseMatrixValue(idir,kdUdZ,kV,dzscatter*p3p4);
01185   currentkp->SetNoiseMatrixValue(idir,kdUdZ,kdUdZ,p3p3);
01186   currentkp->SetNoiseMatrixValue(idir,kdUdZ,kdVdZ,p3p4);
01187   currentkp->SetNoiseMatrixValue(idir,kdUdZ,kQoverP,
01188                                  (sigmaeloss2*TMath::Abs((Double_t)(dplane))*buvz(1)
01189                                   *0.0254*Munits::m*dsdz*(1.+dudz*dudz)));
01190   currentkp->SetNoiseMatrixValue(idir,kdVdZ,kU,dzscatter*p3p4);
01191   currentkp->SetNoiseMatrixValue(idir,kdVdZ,kV,dzscatter*p4p4);
01192   currentkp->SetNoiseMatrixValue(idir,kdVdZ,kdUdZ,p3p4);
01193   currentkp->SetNoiseMatrixValue(idir,kdVdZ,kdVdZ,p4p4);
01194   currentkp->SetNoiseMatrixValue(idir,kdVdZ,kQoverP,
01195                                  (sigmaeloss2*TMath::Abs((Double_t)(dplane))*buvz(0)
01196                                   *0.0254*Munits::m*dsdz*(1.+dvdz*dvdz)));
01197   currentkp->SetNoiseMatrixValue(idir,kQoverP,kU,
01198                                  (dzscatter*sigmaeloss2*TMath::Abs((Double_t)(dplane))
01199                                   *buvz(1)*0.0254*Munits::m*dsdz*(1.+dudz*dudz)));
01200   currentkp->SetNoiseMatrixValue(idir,kQoverP,kV,
01201                                  (dzscatter*sigmaeloss2*TMath::Abs((Double_t)(dplane))
01202                                   *buvz(0)*0.0254*Munits::m*dsdz*(1.+dvdz*dvdz)));
01203   currentkp->SetNoiseMatrixValue(idir,kQoverP,kdUdZ,
01204                                  (sigmaeloss2*TMath::Abs((Double_t)(dplane))*buvz(1)
01205                                   *0.0254*Munits::m*dsdz*(1.+dudz*dudz)));
01206   currentkp->SetNoiseMatrixValue(idir,kQoverP,kdVdZ,
01207                                  (sigmaeloss2*TMath::Abs((Double_t)(dplane))*buvz(0)
01208                                   *0.0254*Munits::m*dsdz*(1.+dvdz*dvdz)));
01209   currentkp->SetNoiseMatrixValue(idir,kQoverP,kQoverP,sigmaeloss2);
01210 
01211   MSG("FitTrackSR",Msg::kDebug) << "du/dz, dv/dz, ds/dz, qp = " << dudz << " " 
01212                                 << dvdz << " " << dsdz << " " << qp << endl
01213                                 << "dplane, dzscatter, z1, z2, radlen = " << dplane 
01214                                 << " " << dzscatter << " " << z1 << " " << z2 << " " 
01215                                 << locradlen << endl
01216                                 << "sigmams2 = " << sigmams2 << endl
01217                                 << "p3p3, p3p4, p4p4 = " << p3p3 << " " << p3p4 << " " 
01218                                 << p4p4 << endl;
01219 
01220 }

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

Definition at line 1380 of file AlgFitTrackSR.cxx.

References TrackClusterSR::GetPlaneView(), KalmanPlaneSR::GetPreCovarianceMatrixValue(), KalmanPlaneSR::GetPreStateValue(), TrackClusterSR::GetTPos(), KalmanPlaneSR::GetTrackCluster(), Msg::kDebug, PlaneView::kV, MSG, and KalmanPlaneSR::SetPreChi2().

Referenced by Predict().

01382 {
01383 
01384   Int_t index = 0;
01385 
01386   //default view is U, change if this is a V view plane
01387   if(currentkp->GetTrackCluster()->GetPlaneView()==PlaneView::kV) index = 1;
01388 
01389   Double_t dtpos = 0.;
01390 
01391   dtpos = currentkp->GetPreStateValue(index,idir)-currentkp->GetTrackCluster()->GetTPos();
01392 
01393   Double_t preChi2 = dtpos*dtpos/(poserr*poserr);
01394   
01395   currentkp->SetPreChi2(preChi2, idir);
01396  MSG("FitTrackSR",Msg::kDebug) << "prechi2 = " << preChi2 << "  dtpos = " 
01397                                 << dtpos << "  tposerr2 = " << poserr*poserr 
01398                                 << "  cov = " 
01399                                 << currentkp->GetPreCovarianceMatrixValue(idir,index,index) 
01400                                 << endl;
01401 
01402 }

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(), Msg::kDebug, kdUdZ, kdVdZ, kQoverP, kU, kV, MSG, and KalmanPlaneSR::SetPreCovarianceMatrixValue().

Referenced by Predict().

01226 { 
01227   MSG("FitTrackSR", Msg::kDebug) << "In CalculatePreCovariance" << endl;
01228   //get the propagator and filtered covariance matricies from the
01229   //previously fit kp
01230   TMatrixD propagator(5,5);
01231   TMatrixD filCovariance(5,5);
01232 
01233   for(Int_t i = 0; i<5; ++i){
01234     for(Int_t j = 0; j<5; ++j){
01235       propagator(i,j) = prevkp->GetPropagatorMatrixValue(idir,i,j);
01236       filCovariance(i,j) = prevkp->GetFilCovarianceMatrixValue(idir,i,j);
01237       
01238       //start the calculation of the pre-filtered covariance matrix
01239       //for the kp to add to the fit.  set it = noise matrix of previously
01240       //fit plane
01241       currentkp->SetPreCovarianceMatrixValue(idir,i,j,prevkp->GetNoiseMatrixValue(idir,i,j));
01242     }
01243   }
01244 
01245   //remember that propagator is really the transpose of the propagator matrix
01246   //so in the following cov0 = F_{k-1}C_{k-1}
01247   TMatrixD cov0(propagator,TMatrixD::kTransposeMult,filCovariance);
01248   
01249   //the following makes cov0 = F_{k-1}C_{k-1}F^{T}_{k-1}
01250   cov0 *= propagator;
01251 
01252   //finish finding the pre-filtered covariance matrix - 
01253   //C_{k}^{k-1} = F_{k-1}C_{k-1}F^{T}_{k-1}+Q_{k-1}
01254   Double_t preCovValue = 0.;
01255   for(Int_t i = 0; i<5; ++i){
01256     for(Int_t j = 0; j<5; ++j){
01257       preCovValue = currentkp->GetPreCovarianceMatrixValue(idir,i,j)+cov0(i,j);
01258       currentkp->SetPreCovarianceMatrixValue(idir,i,j,preCovValue);
01259     }
01260   }
01261 
01262   // diagonal elements should be positive
01263   Double_t covlim = 1.e-8;
01264   if(currentkp->GetPreCovarianceMatrixValue(idir,kU,kU)<covlim) 
01265     currentkp->SetPreCovarianceMatrixValue(idir,kU,kU,covlim);
01266   if(currentkp->GetPreCovarianceMatrixValue(idir,kV,kV)<covlim) 
01267     currentkp->SetPreCovarianceMatrixValue(idir,kV,kV,covlim);
01268   if(currentkp->GetPreCovarianceMatrixValue(idir,kdUdZ,kdUdZ)<covlim) 
01269     currentkp->SetPreCovarianceMatrixValue(idir,kdUdZ,kdUdZ,covlim);
01270   if(currentkp->GetPreCovarianceMatrixValue(idir,kdVdZ,kdVdZ)<covlim) 
01271     currentkp->SetPreCovarianceMatrixValue(idir,kdVdZ,kdVdZ,covlim);
01272   if(currentkp->GetPreCovarianceMatrixValue(idir,kQoverP,kQoverP)<covlim) 
01273     currentkp->SetPreCovarianceMatrixValue(idir,kQoverP,kQoverP,covlim);
01274 
01275 }

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

Definition at line 1278 of file AlgFitTrackSR.cxx.

References fParmMaxAngle, fParmMaxQP, fParmMaxQPFrac, fParmQPRangeCheck, fSwimmer, KalmanPlaneSR::GetFilStateValue(), KalmanPlaneSR::GetPreStateValue(), KalmanPlaneSR::GetRange(), KalmanPlaneSR::GetTrackCluster(), KalmanPlaneSR::GetVldContext(), TrackClusterSR::GetZPos(), Msg::kDebug, kdUdZ, kdVdZ, Msg::kError, kQoverP, kU, kV, MSG, KalmanPlaneSR::SetPreStateValue(), and Swim().

Referenced by Predict().

01281 {
01282 
01283   //Swim from the previously fit plane to the plane to add to the fit
01284   //to get the estimate of the state vector at the new plane - this 
01285   //takes care of the F_{k-1}x_{k-1} matrix multiplication in the
01286   //Kalman filter equations
01287   Int_t nswimfail=0;
01288 
01289   if(fSwimmer!=2) MSG("FitTrackSR",Msg::kError) << "fSwimmer!=2 currently unsupported" 
01290                                                 << endl;
01291   MSG("FitTrackSR", Msg::kDebug) << "current Filtered State = [" 
01292                                  << prevkp->GetFilStateValue(kU,idir) << ","
01293                                  << prevkp->GetFilStateValue(kV,idir) << ","
01294                                  << prevkp->GetFilStateValue(kdUdZ,idir) << ","
01295                                  << prevkp->GetFilStateValue(kdVdZ,idir) << ","
01296                                  << prevkp->GetFilStateValue(kQoverP,idir) << "]" 
01297                                  << endl;
01298     
01299   if(fSwimmer==2){
01300     
01301     // Using Swimmer package
01302     Double_t swimstate[5];
01303     Double_t state[] = {0.,0.,0.,0.,0.};
01304     
01305     swimstate[kU] = prevkp->GetFilStateValue(kU,idir);
01306     swimstate[kV] = prevkp->GetFilStateValue(kV,idir);
01307     swimstate[kdUdZ] = prevkp->GetFilStateValue(kdUdZ,idir);
01308     swimstate[kdVdZ] = prevkp->GetFilStateValue(kdVdZ,idir);
01309     swimstate[kQoverP] = prevkp->GetFilStateValue(kQoverP,idir);
01310     
01311     Bool_t swimstatus = false;
01312     Int_t nswim=0;
01313     
01314     while(!swimstatus && nswim<100){
01315 
01316       swimstatus = Swim(swimstate, state, prevkp->GetTrackCluster()->GetZPos(),
01317                         currentkp->GetTrackCluster()->GetZPos(),idir,
01318                         currentkp->GetVldContext());
01319       
01320       // if outside fiducial detector, return failure flag
01321       if(TMath::Abs(state[kU])>5. || TMath::Abs(state[kV])>5.){
01322        MSG("FitTrackSR",Msg::kDebug) << "outside detector, u,v = " << state[kU] << " " 
01323                                       << state[kV] << " fail" << endl;
01324         
01325         return -1;
01326       }
01327 
01328       if(!swimstatus){
01329         
01330         // swim failed, double momentum
01331         swimstate[kQoverP] *= 0.5;          
01332      MSG("FitTrackSR",Msg::kDebug) << "swim failed, doubling momentum to q/p = " 
01333                                       << swimstate[kQoverP] << endl;
01334 
01335         ++nswimfail;
01336       }
01337       ++nswim;
01338     }//end loop to swim.
01339 
01340     currentkp->SetPreStateValue(kU,idir,state[kU]);
01341     currentkp->SetPreStateValue(kV,idir,state[kV]);
01342     currentkp->SetPreStateValue(kdUdZ,idir,state[kdUdZ]);
01343     currentkp->SetPreStateValue(kdVdZ,idir,state[kdVdZ]);
01344     currentkp->SetPreStateValue(kQoverP,idir,state[kQoverP]);
01345   }//end if using Swim package
01346 
01347   if(TMath::Abs(currentkp->GetPreStateValue(kdUdZ,idir))>fParmMaxAngle){
01348     currentkp->SetPreStateValue(kdUdZ,idir,
01349                                 (currentkp->GetPreStateValue(kdUdZ,idir)>0 
01350                                  ? fParmMaxAngle : -fParmMaxAngle));
01351   }
01352 
01353   if(TMath::Abs(currentkp->GetPreStateValue(kdVdZ,idir))>fParmMaxAngle){
01354     currentkp->SetPreStateValue(kdVdZ,idir,
01355                                 (currentkp->GetPreStateValue(kdVdZ,idir)>0 
01356                                  ? fParmMaxAngle : -fParmMaxAngle));
01357   }
01358 
01359   Double_t maxqp = fParmMaxQP;
01360 
01361   //use a value of about 2 MeV/(g/cm^2) for the energy loss in the material
01362   //you would divide the max qp frac allowed by that value and then by the range to
01363   //get the max qp.  but multiplication is a faster operation, so multiply by
01364   //500 instead.
01365   if ( fParmQPRangeCheck ) {
01366       if(currentkp->GetRange()>0. && fParmMaxQPFrac*500./currentkp->GetRange()<maxqp) {
01367           maxqp = fParmMaxQPFrac*500./currentkp->GetRange();
01368       }
01369   }
01370 
01371   if(TMath::Abs(currentkp->GetPreStateValue(kQoverP,idir))>maxqp){
01372     currentkp->SetPreStateValue(kQoverP,idir,
01373                                 (currentkp->GetPreStateValue(kQoverP,idir)>0 ? maxqp : -maxqp));
01374   }
01375 
01376   return nswimfail;
01377 }

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(), Msg::kDebug, kdUdZ, kdVdZ, Msg::kError, kQoverP, kU, kV, MSG, KalmanPlaneSR::SetPropagatorMatrixValue(), and Swim().

Referenced by Predict().

00956 {
00957   // we are actually calculating the tranpose
00958   MSG("FitTrackSR", Msg::kDebug) << "In CalculatePropagator" << endl;
00959   Int_t nswimfail=0;
00960   Double_t dstate[5];
00961   Double_t thisFilState[5];
00962 
00963   for(Int_t i=0; i<5; ++i){
00964 
00965     //set the filtered state based on the previously fit plane
00966     thisFilState[i] = prevkp->GetFilStateValue(i,idir);
00967     dstate[i] = fParmDState[i];
00968     
00969     if(i==4){
00970       dstate[i] = (fParmDState[i]*thisFilState[kQoverP]*thisFilState[kQoverP]
00971                    /(1.-fParmDState[i]*TMath::Abs(thisFilState[kQoverP])));
00972       dstate[i] = fParmDState[i]*TMath::Abs(thisFilState[kQoverP]);
00973       if(dstate[i]<.01) dstate[i] = .01;
00974       // corresponds to constant fParmDState[kQoverP] energy interval
00975     }
00976     //    MSG("FitTrackSR", Msg::kDebug) << "thisFilState[" << i << "] = " 
00977     //                          << thisFilState[i] << " dstate[" << i 
00978     //                          << "] = " << dstate[i] << endl;
00979   }
00980 
00981   Bool_t swimstatus = false;
00982   Int_t nswim=0;
00983   Double_t swimstate[5];
00984   Double_t pstate[] = {0.,0.,0.,0.,0.};
00985   Double_t nstate[] = {0.,0.,0.,0.,0.};
00986 
00987   //loop over the vectors and swim the particle until you have a good
00988   //swim or you excede the maximum iterations
00989   while(!swimstatus && nswim<100){
00990     
00991     swimstatus = true;
00992 
00993     //loop over all the states and change them by +/- dstate and swim from
00994     //that point to the next z pos. this figures out the propagator matrix
00995     for(Int_t i=0; i<5; ++i){
00996       if(i==4){
00997         dstate[i] = (fParmDState[i]*thisFilState[kQoverP]*thisFilState[kQoverP]
00998                      /(1.-fParmDState[i]*TMath::Abs(thisFilState[kQoverP])));
00999         dstate[i] = fParmDState[i]*TMath::Abs(thisFilState[kQoverP]);
01000         if(dstate[i]<.01) dstate[i] = .01;
01001       }
01002 
01003       swimstate[kU] = thisFilState[kU];
01004       swimstate[kV] = thisFilState[kV];
01005       swimstate[kdUdZ] = thisFilState[kdUdZ];
01006       swimstate[kdVdZ] = thisFilState[kdVdZ];
01007       swimstate[kQoverP] = thisFilState[kQoverP];
01008       
01009       swimstate[i] += dstate[i];
01010 
01011       if(fSwimmer!=2) MSG("FitTrackSR",Msg::kError) << "fSwimmer!=2 currently unsupported" << endl;      
01012       else if(fSwimmer==2){
01013         
01014         // Using Swimmer package
01015         if(!(Swim(swimstate, pstate, prevkp->GetTrackCluster()->GetZPos(),
01016                   currentkp->GetTrackCluster()->GetZPos(), idir, 
01017                   prevkp->GetVldContext()))) swimstatus = false;
01018       
01019 
01020         // if outside fiducial detector, return failure flag
01021         if(TMath::Abs(pstate[kU])>10. || TMath::Abs(pstate[kV])>10.){
01022           return -1;
01023         }
01024         
01025         swimstate[kU] = thisFilState[kU];
01026         swimstate[kV] = thisFilState[kV];
01027         swimstate[kdUdZ] = thisFilState[kdUdZ];
01028         swimstate[kdVdZ] = thisFilState[kdVdZ];
01029         swimstate[kQoverP] = thisFilState[kQoverP];
01030         
01031         swimstate[i] -= dstate[i];
01032         
01033         if(!(Swim(swimstate, nstate, prevkp->GetTrackCluster()->GetZPos(),
01034                   currentkp->GetTrackCluster()->GetZPos(), idir, 
01035                   currentkp->GetVldContext()))) swimstatus = false;
01036       
01037 
01038         // if outside fiducial detector, return failure flag
01039         if(TMath::Abs(nstate[kU])>10. || TMath::Abs(nstate[kV])>10.){
01040                  MSG("FitTrackSR",Msg::kDebug) << "outside detector, pstate u,v = " << pstate[kU] 
01041                                        << " " << pstate[kV] << " fail" << endl;  
01042           return -1;
01043         }
01044   
01045 
01046         //get the value of the propagator matrix terms for the current column of the
01047         //matrix
01048         Double_t propValue = 0.;
01049         for(Int_t j=0; j<5; ++j){
01050         
01051           if(dstate[i]!=0.)  propValue = (pstate[j]-nstate[j])/(2.*dstate[i]);
01052           else{
01053             if(j==i) propValue = 1.;
01054             else propValue = 0.;
01055           }
01056           prevkp->SetPropagatorMatrixValue(idir, i, j, propValue);
01057         }
01058       
01059       }//end if fSwimmer is for the Swim package
01060     }//end loop over matrix rows
01061 
01062     if(!swimstatus){
01063 
01064       // swim failed, double momentum
01065       thisFilState[kQoverP] *= 0.5;
01066       ++nswimfail;
01067  MSG("FitTrackSR",Msg::kDebug) << "failed swim, doubling momentum to q/p = " 
01068                                    << thisFilState[kQoverP] << " nswim = " << nswim << endl;
01069     }
01070 
01071     ++nswim;
01072 
01073   } 
01074 
01075   return nswimfail;
01076 }

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(), fParmMaxIterate2, TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetPlane(), TrackClusterSR::GetPlaneView(), KalmanPlaneSR::GetTrackCluster(), KalmanPlaneSR::GetZDir(), InitializeFit(), IterateKalmanFit(), Msg::kDebug, PlaneView::kU, PlaneView::kV, MSG, RemoveBadPointsFromFit(), ResetTrackClusterList(), TrackClusterSR::SetIsValid(), CandFitTrackSRHandle::SetVtx(), and tc.

Referenced by RunAlg().

02203 {
02204   Int_t nadd=1;
02205   Int_t iadd=0;
02206   istatus=0;
02207 
02208   Int_t nbadfit = 1;
02209   Int_t ibadfit = 0;
02210 
02211   Int_t nu = 0;
02212   Int_t nv = 0;
02213 
02214   TrackClusterSR *tc = 0;
02215 
02216   //make maps of fit planes in each view.  map the plane number
02217   //to an int - 0 for a bad fit and 1 for a good fit
02218   map<Int_t, Int_t> uFitPlanes;
02219   map<Int_t, Int_t> vFitPlanes;
02220   
02221   //loop over the clusters in the slice to fill the uFitPlanes and vFitPlanes
02222   //maps.  initialize all plane values to zero, and reset them in the for
02223   //loop below if the fit was good.  this way you can see how many hit planes
02224   //are being skipped by dropping poorly fit planes - otherwise you might
02225   //just keep removing planes but forgetting that you had the ones you dropped
02226   TIter planeItr(&planeClusterList);
02227   TObjArray * plnarray;
02228   while( (plnarray = (TObjArray * )planeItr.Next()) ){
02229     TIter clusterItr(plnarray);
02230     while( (tc = (TrackClusterSR *)clusterItr.Next()) ){
02231       if(tc->GetPlaneView() == PlaneView::kU) uFitPlanes[tc->GetPlane()] = 0;
02232       else if(tc->GetPlaneView() == PlaneView::kV) vFitPlanes[tc->GetPlane()] = 0;
02233     }
02234   }
02235 
02236 
02237   //loop as long as the fit flag is 0, you have added some points to the 
02238   //track, and you havent gone over the max number of iterations
02239   while(!istatus && nadd && iadd<=fParmMaxIterate2){
02240     
02241     nbadfit=1;
02242     ibadfit=0;
02243 
02244     //loop as long as the fit flag is 0 and you have points to remove from the
02245     //track
02246     while(!istatus && nbadfit){
02247       
02248       //reset the number of bad points in the fit to zero for this loop
02249       nbadfit = 0;
02250       
02251       //iterate the Kalman fit forward and backwards until the fit q/p 
02252       //converges
02253       
02254       istatus = IterateKalmanFit(planeClusterList, cfh, nu, nv, direction, dosearch);
02255       
02256       //if all systems are still go, look for and remove any bad points in the 
02257       //track
02258       
02259       if (!istatus){
02260         MSG("FitTrackSR",Msg::kDebug) << "Removing iteration " << iadd << "/" 
02261                                       << ibadfit << endl;
02262         
02263         nbadfit = RemoveBadPointsFromFit(cfh, planeClusterList, uFitPlanes, vFitPlanes,nu, nv, direction);
02264         
02265       }//end if tracking worked
02266       
02267       ++ibadfit;
02268     }//end loop while there are still bad points or tracking
02269     
02270     //now look for planes to add to the fit
02271     //reset the counter for the number of added points to the track
02272     nadd = 0;
02273     
02274     //make sure the fit is still ok
02275     if(!istatus && dosearch){
02276       
02277       cfh.ClearKalmanPlaneList();
02278       
02279       //add the remaining clusters after removal of bad points to the fit.  
02280       //use 999 as the value for the third parameter
02281       //to ensure that we just scale the previous covariance matrix
02282       nu = 0; 
02283       nv = 0;
02284       
02285       AddClustersToFit(planeClusterList, cfh, 999, nu, nv, direction);
02286       //make sure you have enough planes in both view to procede
02287       if(nu>=3 && nv>=3){
02288         
02289         // find additional planes downstream
02290         nadd += FindDownstreamPlanes(planeClusterList, cfh, direction);
02291         
02292         //reinitialize the fit and use 
02293         // 999 to ensure we simply scale previous covariance matrix
02294         InitializeFit(cfh,1,999); 
02295         
02296         //make an array to hold all the new kp's you find
02297         TObjArray *newkplist = new TObjArray(1,0);
02298         
02299         nadd += FindUpstreamPlanes(planeClusterList, cfh, newkplist, direction);
02300         
02301         //reset the track clusters in the clusterItr for the next iteration
02302         ResetTrackClusterList(planeClusterList);
02303         MSG("FitTrackSR",Msg::kDebug) << "Last iteration " << iadd << "/" << ibadfit 
02304                                       << endl;
02305         //clear the cfh Kalman Plane list - those objects are all in the
02306         //newkplist and will get deleted once the clusters have all their 
02307         //IsValid flags reset
02308         
02309         cfh.ClearKalmanPlaneList(0); 
02310         const KalmanPlaneSR *currentkp = 0;
02311         Int_t sizeoflist = newkplist->GetLast();
02312         
02313         //loop over the newkplist in reverse to mark the clusters used in the fit
02314         KalmanPlaneSR *kp = 0;
02315         for(Int_t i=sizeoflist; i>=0; --i){
02316           kp = dynamic_cast<KalmanPlaneSR*>(newkplist->At(i));
02317           assert(kp);
02318           
02319           if(!currentkp) currentkp = cfh.AddKalmanPlaneAt(kp,sizeoflist-i);
02320           else cfh.AddKalmanPlaneAt(kp,sizeoflist-i);
02321          
02322           TObjArray * planeClusters=(TObjArray *)planeClusterList.At(kp->GetTrackCluster()->GetPlane());
02323           TIter clusterItr(planeClusters);
02324           while((tc = (TrackClusterSR *)(clusterItr.Next()) )){
02325             if(tc->GetPlane()==kp->GetTrackCluster()->GetPlane()){
02326               if(tc->GetMinStrip() == kp->GetTrackCluster()->GetMinStrip()
02327                  && tc->GetMaxStrip() == kp->GetTrackCluster()->GetMaxStrip() )
02328                 tc->SetIsValid(true);
02329             }
02330           }
02331           MSG("FitTrackSR",Msg::kDebug) << "  TC " << kp->GetTrackCluster()->GetPlane() 
02332                                           << " " << kp->GetTrackCluster()->GetMinStrip() 
02333                                         << " " << kp->GetTrackCluster()->GetMaxStrip() 
02334                                         << " " << kp->GetZDir() 
02335                                         << endl;
02336           delete kp;
02337           
02338         }//end loop over new kplist
02339         
02340         MSG("FitTrackSR", Msg::kDebug) << "vertex direction = " 
02341                                        << currentkp->GetZDir()
02342                                        << endl;
02343         //set the vertex to the first kp
02344         cfh.SetVtx(currentkp);
02345         
02346         newkplist->Clear();
02347         delete newkplist;
02348         
02349       }//end if enough planes hit to start looking for more to add
02350     }//end if fit is still ok after iteration ste
02351     ++iadd;  
02352   }//end while !iadd
02353   
02354   //the iadd-1 is so that we can just put the returned value right into
02355   //cfh.SetNIterate()
02356   return iadd-1;
02357 }

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

Definition at line 840 of file AlgFitTrackSR.cxx.

References CalculateFilChi2(), CalculateFilCovariance(), CalculateFilState(), CalculateGain(), fParmMaxAngle, fParmMaxQP, fParmMaxQPFrac, fParmQPRangeCheck, KalmanPlaneSR::GetFilStateValue(), KalmanPlaneSR::GetRange(), TrackClusterSR::GetTPosError(), KalmanPlaneSR::GetTrackCluster(), kdUdZ, kdVdZ, kQoverP, KalmanPlaneSR::SetFilStateValue(), and KalmanPlaneSR::SetPredictPlane().

Referenced by FitFrom().

00841 {
00842 //
00843 // return value =  0        success
00844 //
00845 
00846   //get the K_{k} matrix
00847   CalculateGain(currentkp, idir,currentkp->GetTrackCluster()->GetTPosError());
00848   
00849   //do the C_{k} = (I-K_{k}H_{k})C_{k}^{k-1} step
00850   CalculateFilCovariance(currentkp, idir);
00851 
00852   //find the filtered state
00853   CalculateFilState(currentkp, idir);
00854 
00855   // calculate noise and precovariance matrices onward again, values may have changed
00856 
00857   CalculateFilChi2(currentkp, idir,currentkp->GetTrackCluster()->GetTPosError());
00858 
00859   //reset the predicted from plane -- filtered state has changed
00860   currentkp->SetPredictPlane(-999, idir); 
00861 
00862   if(TMath::Abs(currentkp->GetFilStateValue(kdUdZ,idir))>fParmMaxAngle){
00863     if(currentkp->GetFilStateValue(kdUdZ,idir)>0.) 
00864       currentkp->SetFilStateValue(kdUdZ, idir, fParmMaxAngle);
00865     if(currentkp->GetFilStateValue(kdUdZ,idir)<0.) 
00866       currentkp->SetFilStateValue(kdUdZ, idir, -fParmMaxAngle);;
00867   }
00868   if(TMath::Abs(currentkp->GetFilStateValue(kdVdZ, idir))>fParmMaxAngle){
00869     if(currentkp->GetFilStateValue(kdVdZ, idir)>0.) 
00870       currentkp->SetFilStateValue(kdVdZ, idir, fParmMaxAngle);
00871     if(currentkp->GetFilStateValue(kdVdZ, idir)<0.) 
00872       currentkp->SetFilStateValue(kdVdZ, idir, -fParmMaxAngle);
00873   }
00874 
00875   Double_t maxqp = fParmMaxQP;
00876 
00877   //use a value of about 2 MeV/(g/cm^2) for the energy loss in the material
00878   //you would divide the max qp frac allowed by that value and then by the range to
00879   //get the max qp.  but multiplication is a faster operation, so multiply by
00880   //500 instead.
00881   if ( fParmQPRangeCheck ) {
00882       if(currentkp->GetRange()>0. && fParmMaxQPFrac*500./currentkp->GetRange()<maxqp) {
00883           maxqp = fParmMaxQPFrac*500./currentkp->GetRange();
00884       }
00885   }
00886 
00887   if(idir==0 && TMath::Abs(currentkp->GetFilStateValue(kQoverP,idir))>maxqp){
00888     if(currentkp->GetFilStateValue(kQoverP,idir)>0.) 
00889       currentkp->SetFilStateValue(kQoverP, idir, maxqp);
00890     else currentkp->SetFilStateValue(kQoverP, idir, -maxqp);
00891   }
00892  
00893   return 0;
00894 }

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

Definition at line 3040 of file AlgFitTrackSR.cxx.

References AddForwardBestKPToFit(), FindNumSkippedPlanes(), FitFrom(), fParmNSkipActive, CandFitTrackSRHandle::GetKalmanLast(), CandFitTrackSRHandle::GetKalmanPlane(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetPlane(), KalmanPlaneSR::GetPreChi2(), CandTrackHandle::GetRange(), KalmanPlaneSR::GetTrackCluster(), TrackClusterSR::IsValid(), Msg::kDebug, MSG, KalmanPlaneSR::SetRange(), and KalmanPlaneSR::SetZDir().

Referenced by DoKalmanFit().

03043 {
03044  
03045 
03046  MSG("FitTrackSR", Msg::kDebug) << "in FindDownstreamPlanes " 
03047                                  << direction << endl;
03048   TrackClusterSR *thistc = 0;
03049   TrackClusterSR *oldtc = 0;
03050   KalmanPlaneSR *thiskp = 0;
03051   KalmanPlaneSR *bestkp = 0;
03052   KalmanPlaneSR *oldkp = cfh.GetKalmanPlane(cfh.GetKalmanLast());
03053   Int_t nswimfail = 0;
03054   Int_t nplaneskip = 0;
03055   Bool_t isvalid = true;
03056   Int_t istatus2 = 0;
03057   Int_t nadd = 0;
03058 
03059   Int_t currentPlane = -1;
03060   Int_t oldKPPlane = oldkp->GetTrackCluster()->GetPlane();
03061 
03062   //loop over the clusters in the cluster list
03063 
03064   //save some time - slice the cluster iterator from the last plane
03065   //in the fit to the last one in the detector in the appropriate direction
03066   //neither detector has more than 500 planes, and both start with plane 0 being
03067   //the front steel (uninstrumented) plane
03068 
03069   Bool_t dir = kIterForward;
03070   if(direction!=1) dir=kIterBackward;
03071   TIter planeItr(&planeClusterList,dir);
03072   TObjArray * plnarray;
03073   while( (plnarray = (TObjArray * )planeItr.Next()) ){
03074     TIter clusterItr(plnarray);
03075     while( (thistc = (TrackClusterSR * )clusterItr.Next()) && isvalid ){
03076           //make sure you are looking at a valid track cluster, if not go on to the next one
03077           if( !thistc->IsValid() ) continue;         
03078       currentPlane = thistc->GetPlane();
03079       
03080       MSG("FitTrackSR", Msg::kDebug) << "looking at plane " << currentPlane << endl;
03081       
03082       //make sure you are looking at a downstream plane
03083       if(direction*currentPlane>direction*oldKPPlane){
03084         
03085         //see if the current track cluster is on a new plane and that
03086         //you have a potential kp to add to the fit
03087         if(bestkp && oldtc && oldtc->GetPlane() != currentPlane ){
03088           MSG("FitTrackSR", Msg::kDebug ) << "add best kp to fit" << endl;
03089           if(AddForwardBestKPToFit(bestkp, cfh, nswimfail)){
03090             ++nadd;
03091             oldkp = bestkp;
03092           }
03093           bestkp = 0;
03094         }
03095         
03096         oldKPPlane = oldkp->GetTrackCluster()->GetPlane();
03097         MSG("FitTrackSR", Msg::kDebug) << " prev plane = " << oldKPPlane 
03098                                        << " direction " << direction << endl;
03099         
03100         //on a new cluster, see if it is a good one
03101         MSG("FitTrackSR",Msg::kDebug) << "considering tc " << thistc->GetPlane() 
03102                                       << " " << thistc->GetMinStrip() << "-" 
03103                                       << thistc->GetMaxStrip() 
03104                                       << " for addition to end" << endl;
03105         nplaneskip = 0;
03106         
03107         //check if the current plane is far from the last kp - need to
03108         //see how many active planes were skipped if it is
03109         if(TMath::Abs(currentPlane-oldKPPlane)>=2*fParmNSkipActive)
03110           nplaneskip = FindNumSkippedPlanes(currentPlane,oldKPPlane,oldkp, direction);
03111         
03112         //make sure you didnt skipp too many active planes
03113         if(nplaneskip<2*fParmNSkipActive){
03114           
03115           //make a new kp for this cluster to test it out
03116           thiskp = new KalmanPlaneSR(thistc);
03117           thiskp->SetRange(cfh.GetRange(thiskp->GetTrackCluster()->GetPlane()));
03118           thiskp->SetZDir(direction);
03119           istatus2 = FitFrom(oldkp,thiskp,0,-1);
03120           MSG("FitTrackSR",Msg::kDebug) << "  chi2 = " << thiskp->GetPreChi2(0)
03121                                         << " direction " << direction
03122                                         << " nswimfail = " << istatus2 << endl;
03123           //look to see if the kp you just made is better than the previous
03124           //best kp
03125           if(!istatus2 && (!bestkp || thiskp->GetPreChi2(0)<bestkp->GetPreChi2(0))){
03126             MSG("FitTrackSR",Msg::kDebug) << "best chi2, keeping" << endl;   
03127             if (bestkp) delete bestkp;
03128             bestkp = thiskp;
03129             
03130             nswimfail = istatus2;
03131           }
03132           else delete thiskp;
03133           
03134           if(istatus2<0) isvalid = 0;
03135         } 
03136         else {
03137           // any planes beyond this one will not satisfy NSkipView, stop search
03138           isvalid = 0;
03139         }
03140         MSG("FitTrackSR",Msg::kDebug) << "set oldtc" << endl;
03141         oldtc = thistc;
03142       }//end if current plane is down stream of the last fit plane
03143     }//end loop over cluster list
03144   }
03145 
03146   //may still have a best kp, so make sure to add it to the fit
03147   if(bestkp){
03148     MSG("FitTrackSR", Msg::kDebug) << "add bestkp outside of loop over tc" << endl;
03149     nadd += AddForwardBestKPToFit(bestkp, cfh, nswimfail);
03150     
03151     //delete bestkp so as to not have a memory leek
03152     bestkp = 0;
03153   }
03154 
03155   return nadd;
03156 }

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, fParmNSkipActive, fVldContext, PlexPlaneId::GetAdjoinScint(), KalmanPlaneSR::GetFilStateValue(), PlexPlaneId::GetPlaneCoverage(), UgliGeomHandle::GetScintPlnHandle(), KalmanPlaneSR::GetTrackCluster(), UgliPlnHandle::GetZ0(), TrackClusterSR::GetZPos(), PlexPlaneId::IsValid(), UgliScintPlnHandle::IsValid(), kdUdZ, kdVdZ, PlaneCoverage::kNoActive, kQoverP, kU, kV, and Swim().

Referenced by FindDownstreamPlanes(), and FindUpstreamPlanes().

01908 {
01909 
01910   UgliGeomHandle ugh(fVldContext);
01911 
01912   // need to check how many active planes were actually skipped
01913   Int_t nplaneskip=0;
01914   
01915   PlexPlaneId plnid(fDetector,prevPlane,false);
01916   PlexPlaneId plnidend(fDetector,currentPlane,false);
01917   plnid = plnid.GetAdjoinScint(direction);
01918   
01919   // setup kalmanplane to do swimming
01920   Double_t currentstate[5] = {oldkp->GetFilStateValue(kU,0),
01921                               oldkp->GetFilStateValue(kV,0),
01922                               oldkp->GetFilStateValue(kdUdZ,0),
01923                               oldkp->GetFilStateValue(kdVdZ,0),
01924                               oldkp->GetFilStateValue(kQoverP,0)};
01925   Double_t newstate[] = {0.,0.,0.,0.,0.};
01926   
01927   Double_t currentz = oldkp->GetTrackCluster()->GetZPos();
01928   
01929   bool isvalid = true;
01930   bool swimflag = true;
01931   
01932   while(isvalid && plnid != plnidend && nplaneskip<2*fParmNSkipActive){
01933     
01934     if(plnid.IsValid() && plnid.GetPlaneCoverage()!=PlaneCoverage::kNoActive){
01935       
01936       UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid);
01937       if(scintpln.IsValid()){
01938         
01939         // plane is active, check if in coil
01940         swimflag = Swim(currentstate,newstate,currentz,scintpln.GetZ0(),0,
01941                             &fVldContext);
01942         
01943         if(!swimflag || TMath::Abs(newstate[kU])>10. || TMath::Abs(newstate[kV])>10.) {
01944 
01945           // if failed swim, set nplaneskip to very large number
01946           nplaneskip = 9999999;
01947           isvalid = 0;
01948         } 
01949         else{
01950           currentz = scintpln.GetZ0();
01951           currentstate[kU] = newstate[kU];
01952           currentstate[kV] = newstate[kV];
01953           currentstate[kdUdZ] = newstate[kdUdZ];
01954           currentstate[kdVdZ] = newstate[kdVdZ];
01955           currentstate[kQoverP] = newstate[kQoverP];
01956         }
01957 
01958         // check if not in coil (defined to have radius 30 cm)
01959         if(currentstate[kU]*currentstate[kU]+currentstate[kV]*currentstate[kV]>0.09) {
01960           ++nplaneskip;
01961         
01962         }
01963 
01964 
01965       }//end if scintillator plane is valid
01966     }  
01967     plnid = plnid.GetAdjoinScint(direction);
01968   }//end loop over planes
01969   return nplaneskip;
01970 }

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.

References fParmNSkipView.

Referenced by RemoveBadPointsFromFit().

03339 {
03340  
03341   //assume you drop the current plane, so start off with numSkipped=1
03342   Int_t numSkipped = 1;
03343   Int_t plane = 0;
03344 
03345   if(prevPlane<=0) return 0;
03346   fitIter = fitPlanes.find(currentPlane);
03347   //only do the check if there was a plane previous to this one
03348   if(prevPlane>0){
03349     if(direction == 1) --fitIter;
03350     else ++fitIter;
03351 
03352     plane = (*fitIter).first;
03353   
03354     //loop the iterator in reverse until you find the previously fit plane
03355     while( plane*direction > prevPlane*direction && numSkipped<fParmNSkipView){
03356       if(direction == 1) --fitIter;
03357       else ++fitIter;
03358       plane = (*fitIter).first;
03359 
03360       //went through the loop, increment the number of skipped planes
03361   
03362     ++numSkipped;
03363 
03364       //stop looking once you find a kp with a good fit
03365       if( (*fitIter).second != 0 ) break;
03366     }
03367   }//end if not the first plane in view
03368   
03369   //now look downstream
03370   if(numSkipped<fParmNSkipView){
03371     
03372     //didnt skip too many before the current plane, what about after it?
03373     fitIter = fitPlanes.find(currentPlane);
03374     if(direction == 1) ++fitIter;
03375     else --fitIter;
03376 
03377     plane = (*fitIter).first;
03378     
03379     while( plane*direction < lastPlane*direction && numSkipped<fParmNSkipView ){
03380       if(direction == 1) ++fitIter;
03381       else --fitIter;
03382       plane = (*fitIter).first;
03383       
03384       //stop looking once you find a kp with a good fit
03385       //the break point is before the count increase in 
03386       //the forward loop so that you
03387       //dont double count the current plane as being skipped.  you have
03388       //set the iterator such that it should be one 
03389       //plane beyond the current plane
03390       //and if that plane is a good fit, you didnt skip 2 planes, just one
03391 
03392       if( (*fitIter).second != 0) break;
03393       
03394       //went through the loop, increment the number of skipped planes
03395       ++numSkipped;
03396       
03397     }
03398     
03399   }//end check for second view
03400   
03401   return numSkipped;
03402 }

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(), CalDigitType::kPE, StripEnd::kUnknown, ArrivalTime::Weight(), and LinearFit::Weighted().

04207 {
04208 
04209  
04210 
04211   ArrivalTime arrtime;
04212 
04213   Int_t plane = 0;
04214   const Double_t min_wgt = 0.4;
04215   StripEnd::StripEnd_t stripEnd = StripEnd::kUnknown;
04216 
04217   CandStripHandle *strip = 0;
04218   CandDigitHandle *digit = 0;
04219 
04220   CandStripHandleItr stripItr(cfh.GetDaughterIterator());
04221   stripItr.Reset();
04222 
04223   //arrays for determining the timing direction - limit to
04224   //1000 entries - really shouldnt need more than that to 
04225   //be able to figure out if the event is going north or 
04226   //south
04227   Double_t time[1000];
04228   Double_t pathLength[1000];
04229   Double_t weight[1000];
04230   Double_t digitpe = 0.;
04231   Int_t digitCtr = 0;
04232 
04233   for(Int_t i = 0; i<1000; ++i){
04234     time[i] = 0.;
04235     pathLength[i] = 0.;
04236     weight[i] = 0.;
04237   }
04238   
04239   //loop over all strips in the track
04240   while( (strip = stripItr()) ){
04241     
04242     //only look at double ended strips and those that have
04243     //believable transverse positions
04244     if(strip->GetNDaughters() == 2){
04245       
04246       plane = strip->GetPlane();
04247       //get the iterator over the digits in the strip
04248       CandDigitHandleItr digitItr(strip->GetDaughterIterator());
04249 
04250       while( (digit = digitItr()) ){
04251 
04252         stripEnd = digit->GetPlexSEIdAltL().GetEnd();
04253 
04254         pathLength[digitCtr] = cfh.GetdS()-cfh.GetdS(plane);
04255 
04256         //time is recorded in seconds, but i prefer working in ns.
04257         //time is from the track not the strips
04258         time[digitCtr] = cfh.GetT(plane,stripEnd)*1.e9;
04259 
04260         digitpe = digit->GetCharge(CalDigitType::kPE);
04261         weight[digitCtr] = TMath::Min(arrtime.Weight(digitpe),min_wgt);
04262 
04263         ++digitCtr;
04264 
04265       }//end loop over digits in strip
04266     }//end if double ended strip
04267   }//end loop over strips in track
04268 
04269   Double_t efitparm[2];
04270   Double_t maxresid = 11.;
04271   Double_t resid = 0.;
04272   Int_t ctr = digitCtr;
04273   Int_t imaxresid = 0;
04274   Int_t nremove = 0;
04275   while(maxresid>10
04276         && digitCtr-nremove>4 
04277         && (1.*nremove)<0.2*digitCtr
04278         ){
04279     LinearFit::Weighted(ctr,pathLength,time,weight,fitparm,efitparm);
04280     maxresid = 0.;
04281     imaxresid = 0;
04282     for(int i=0; i<ctr; ++i){
04283       resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i]);
04284       if(weight[i]>0. && TMath::Abs(resid)>maxresid){
04285         maxresid = TMath::Abs(resid);
04286         imaxresid = i;
04287       }
04288     }    
04289     if(maxresid>10){
04290       weight[imaxresid] = 0.;
04291       nremove++;
04292     }
04293     
04294     timefitchi2 = 0.;
04295     for(int i=0; i<ctr; ++i){
04296       resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i]);
04297       timefitchi2 += weight[i]*resid*resid;
04298     }
04299   
04300   }//end loop to find timing fit and remove outliers
04301   
04302   timefitchi2 = 0.;
04303   for(int i=0; i<ctr; ++i){
04304     if(weight[i]>0.){
04305     
04306       resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i]);
04307       timefitchi2 += weight[i]*resid*resid;
04308     }
04309   }
04310   
04311   //check the chi^2 for the fit - if it is really bad dont change the
04312   //initial guess at the direction for the event, ie just make sure
04313   //that the slope in time vs pathlength is positive;
04314 //   if(fitparm[1]<0. && timefitchi2>=10.*(digitCtr-nremove))
04315 //     fitparm[1] *= -1.;
04316 
04317   return digitCtr-nremove;
04318 }

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

Definition at line 4104 of file AlgFitTrackSR.cxx.

References TrackClusterSR::GetBegTime(), TrackClusterSR::GetCharge(), TrackClusterSR::GetZPos(), min, Munits::ns, ArrivalTime::Weight(), and LinearFit::Weighted().

04105 {
04106   Int_t idir = 1;
04107 
04108   ArrivalTime arrtime;
04109 
04110   // determine directionality using timing - this should really
04111   //only be needed for cosmic ray events, as beam events always
04112   //increase in z with increasing time.
04113   
04114   //we could try to determine the timing for each side (E, W) of each
04115   //electronics crate.  that would be the most accurate way, but it 
04116   //would also take a lot of time and iterations.  instead, we can 
04117   //probably get away with just doing a global fit to all datapoints
04118   //since we do have timing constants that take care of the side to side
04119   //and crate to crate offsets.  even if a bit of timing hardware is 
04120   //exchanged, this wont really affect our ability to go back and figure
04121   //out the right calibrations later.
04122   
04123   //make some arrays to do the fit. allow a maximum of 1000 points
04124   //to be included in the fit.  there may be events with more than
04125   //1000 digits, but that should be enough to tell you which way the
04126   //event is going.
04127   Double_t zpos[1000];
04128   Double_t time[1000];
04129   Double_t weight[1000];
04130   
04131   for (int i=0; i<1000; i++) {
04132     zpos[i] = 0.;
04133     time[i] = 0.;
04134     weight[i] = 0.;
04135   }
04136   
04137   //minimum weight for a signal in the timing fit
04138   // maximum weight corresponds to 10 p.e.; reason for this is to minimize
04139   // weight of showers in which transverse spread can have negative effect
04140   // on timing
04141   const Double_t min_wgt = 0.4;
04142   TrackClusterSR *cluster = 0;
04143 
04144   Int_t clusterCtr = 0;
04145 
04146   clusterItr.ResetFirst();
04147 
04148   while( (cluster = clusterItr()) && clusterCtr<1000){
04149     
04150     //get the z position of the plane
04151     zpos[clusterCtr] = (Double_t)cluster->GetZPos();
04152 
04153     time[clusterCtr] = cluster->GetBegTime();
04154 
04155     weight[clusterCtr] = min(arrtime.Weight(cluster->GetCharge()),min_wgt);
04156     ++clusterCtr;
04157 
04158   }//end loop over clusters to fill timing arrays
04159   
04160   //reset the iterator over all strips
04161   clusterItr.ResetFirst();
04162   
04163   //arrays to hold the slope, intercept, and chi^2 of the 
04164   //least squares fit for the timing
04165   Double_t parm[2],eparm[2];
04166   
04167   //set the maximum residual to the input parameter+1ns to 
04168   //get ready for the while loop below
04169   Double_t fParmMaxTimingResid = 20*Munits::ns;
04170   Double_t maxresid = fParmMaxTimingResid+1.*Munits::ns;
04171   Double_t resid = 0.;
04172   Int_t imaxresid = 0;
04173   Int_t nremove=0;
04174   Int_t fitflag = 0;
04175 
04176   while(clusterCtr-nremove>4 && maxresid>fParmMaxTimingResid && !fitflag){
04177 
04178     fitflag = LinearFit::Weighted(clusterCtr,zpos,time,weight,parm,eparm);
04179     maxresid = 0.;
04180     imaxresid = 0;
04181     for (int i=0; i<clusterCtr; ++i) {
04182       resid = parm[0]+parm[1]*zpos[i]-time[i];
04183       if (weight[i]>0. && TMath::Abs(resid)>maxresid) {
04184         maxresid = TMath::Abs(resid);
04185         imaxresid = i;
04186       }
04187     }
04188     if (maxresid>fParmMaxTimingResid) {
04189       weight[imaxresid] = 0.;
04190       ++nremove;
04191     }
04192   }//end loop to fit the timing values
04193 
04194   if(parm[1]<0.){
04195     idir=-1;
04196   }
04197 
04198   return idir;
04199 }

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(), fParmNSkipActive, CandFitTrackSRHandle::GetKalmanLast(), CandFitTrackSRHandle::GetKalmanPlane(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetPlane(), KalmanPlaneSR::GetPreChi2(), CandTrackHandle::GetRange(), KalmanPlaneSR::GetTrackCluster(), KalmanPlaneSR::GetZDir(), TrackClusterSR::IsValid(), Msg::kDebug, MSG, CandFitTrackSRHandle::SetCurrentKalmanPlane(), KalmanPlaneSR::SetRange(), and KalmanPlaneSR::SetZDir().

Referenced by DoKalmanFit().

03165 {
03166   MSG("FitTrackSR", Msg::kDebug) << "in FindUpstreamPlanes " 
03167                                  << direction << endl;
03168   TrackClusterSR *thistc = 0;
03169   KalmanPlaneSR *thiskp = 0;
03170   KalmanPlaneSR *bestkp = 0;
03171   KalmanPlaneSR *oldkp = cfh.GetKalmanPlane(cfh.GetKalmanLast());
03172   Int_t nswimfail = 0;
03173   Int_t nplaneskip = 0;
03174   Bool_t isvalid(1);
03175   Int_t istatus2 = 0;
03176   
03177   Int_t nadd = 0;
03178 
03179   //loop over the clusters in the cluster list
03180 
03181   // identify planes which already are fit
03182   map<Int_t,KalmanPlaneSR *> fitplane;
03183 
03184   KalmanPlaneSR *kp = 0;
03185   for(int i=0; i<=cfh.GetKalmanLast(); ++i){
03186     kp = cfh.GetKalmanPlane(i);
03187     MSG("FitTrackSR",Msg::kDebug) << "  plane " << kp->GetTrackCluster()->GetPlane() 
03188                                   << "  strips " << kp->GetTrackCluster()->GetMinStrip() 
03189                                   << "-" << kp->GetTrackCluster()->GetMaxStrip() 
03190                                   << " " << kp->GetZDir() << endl;
03191     
03192 
03193     fitplane[kp->GetTrackCluster()->GetPlane()] = kp;
03194   }
03195   MSG("FitTrackSR",Msg::kDebug) << "After finding end hits" << endl;
03196  
03197   Int_t currentPlane = oldkp->GetTrackCluster()->GetPlane();
03198   Int_t oldKPPlane = oldkp->GetTrackCluster()->GetPlane();
03199 
03200   Int_t prevPlane = -1;
03201   Int_t firstPlane = cfh.GetKalmanPlane(0)->GetTrackCluster()->GetPlane();
03202   MSG("FitTrackSR",Msg::kDebug) << "REVERSING FIT, looking for good track clusters" 
03203                                << endl;
03204 
03205   newkplist->AddLast(oldkp);
03206 
03207   nswimfail = 0;
03208   isvalid = 1;
03209   map<Int_t,KalmanPlaneSR *>::iterator planeiter;
03210 
03211   Bool_t dir = kIterBackward;
03212   if(direction!=1) dir=kIterForward;
03213   TIter planeItr(&planeClusterList,dir);
03214 
03215   TObjArray * plnarray;
03216   while(isvalid &&  (plnarray = (TObjArray * )planeItr.Next()) ){
03217 
03218     TIter clusterItr(plnarray);
03219     thistc = (TrackClusterSR *)clusterItr.Next();
03220     if( (direction==1 && thistc->GetPlane()<oldKPPlane) ||
03221         (direction!=1 && thistc->GetPlane()>oldKPPlane)){
03222       clusterItr.Reset();
03223       while(isvalid &&  (thistc = (TrackClusterSR *)clusterItr.Next()) ){
03224           //if this isnt a valid track cluster, go on to the next one
03225           if(!thistc->IsValid()) continue;
03226         currentPlane = thistc->GetPlane();
03227         oldKPPlane = oldkp->GetTrackCluster()->GetPlane();
03228 
03229         planeiter = fitplane.find(thistc->GetPlane());
03230         MSG("FitTrackSR",Msg::kDebug) << "considering trackcluster, plane=" 
03231                                       << thistc->GetPlane() << " strips=" 
03232                                       << thistc->GetMinStrip() << "-" 
03233                                       << thistc->GetMaxStrip() << endl;
03234         //if you are on a new plane and have a bestkp , try adding it to the fit
03235         if(bestkp && currentPlane != prevPlane){
03236           nadd += AddReverseBestKPToFit(bestkp, oldkp, cfh, newkplist,
03237                                         nswimfail, direction);
03238           MSG("FitTrackSR" , Msg::kDebug) << "bestkp has direction " 
03239                                           << bestkp->GetZDir() << endl; 
03240           bestkp = 0;
03241         }
03242           
03243         //see if this plane contains a fit track cluster
03244         if(planeiter!=fitplane.end() ){
03245           
03246           //this plane contains a fit track cluster, so add it to the 
03247           //newkplist
03248           if(currentPlane != prevPlane ){
03249             kp = dynamic_cast<KalmanPlaneSR *>(planeiter->second);
03250             MSG("FitTrackSR", Msg::kDebug) << oldkp->GetTrackCluster()->GetPlane() << " "
03251                                            << kp->GetTrackCluster()->GetPlane() << " " 
03252                                            << kp->GetZDir() << endl;
03253             
03254             nswimfail = FitFrom(oldkp,kp,1,-1);
03255             MSG("FitTrackSR", Msg::kDebug) << "adding previously fit plane " 
03256                                            << thistc->GetPlane() << " to list" << " " 
03257                                            << kp->GetZDir() << endl;
03258             
03259             newkplist->AddLast(kp);
03260             MSG("FitTrackSR", Msg::kDebug) << "added" << endl;
03261             cfh.SetCurrentKalmanPlane(kp);
03262             oldkp = kp;
03263           }
03264         }//end if plane contains a fit track cluster 
03265         else{   
03266  
03267           MSG("FitTrackSR",Msg::kDebug) << "this plane not fit, considering tc " 
03268                                         << thistc->GetPlane() << " " 
03269                                         << thistc->GetMinStrip() << "-" 
03270                                         << thistc->GetMaxStrip() 
03271                                         << " for addition to end" << endl;
03272           // need to check how many active planes were actually skipped
03273           // only do check if beyond first fit plane
03274           nplaneskip=0;
03275           
03276           //only worry about checking how many planes were skipped
03277           //if the currentPlane is upstream of the first plane
03278           // - the gaps are of acceptable size between fit planes already
03279           if(currentPlane*direction<firstPlane*direction 
03280              && TMath::Abs(currentPlane-oldKPPlane)>=2*fParmNSkipActive)
03281             nplaneskip = FindNumSkippedPlanes(currentPlane, oldKPPlane, oldkp, -direction);
03282           
03283           //not too many planes skipped, give this cluster a whirl
03284           if(nplaneskip<2*fParmNSkipActive){
03285             thiskp = new KalmanPlaneSR(thistc);
03286             thiskp->SetRange(cfh.GetRange(thiskp->GetTrackCluster()->GetPlane()));
03287             thiskp->SetZDir(oldkp->GetZDir());
03288             MSG("FitTrackSR",Msg::kDebug) << "fitting from track cluster, plane " 
03289                                           << oldkp->GetTrackCluster()->GetPlane() 
03290                                           << " strips " 
03291                                           << oldkp->GetTrackCluster()->GetMinStrip() 
03292                                           << "-" << oldkp->GetTrackCluster()->GetMaxStrip() 
03293                                           << " direction " << direction << "/" 
03294                                           << thiskp->GetZDir()
03295                                           << endl;
03296             istatus2 = FitFrom(oldkp,thiskp,1,-1);
03297             //if the fit was good and this kp has a better chi^2 than the previous
03298             //best one for the plane, make it bestkp
03299             MSG("FitTrackSR",Msg::kDebug) << "  chi2 = " << thiskp->GetPreChi2(1) 
03300                                           << " nswimfail = " << istatus2 << endl;
03301             
03302             if(!istatus2 && (!bestkp || thiskp->GetPreChi2(1)<bestkp->GetPreChi2(1))){
03303               if (bestkp) delete bestkp;
03304               MSG("FitTrackSR",Msg::kDebug) << "best chi2, keeping" << endl;
03305               bestkp = thiskp;
03306               nswimfail = istatus2;
03307             }
03308             else delete thiskp;
03309             
03310           }//end not too many active planes skipped
03311           else{
03312             
03313             // any planes beyond this one will not satisfy NSkipView, stop search
03314             isvalid = 0;
03315           }
03316         }//end current plane does not have a fit clustser
03317         prevPlane = currentPlane;
03318       }//end loop in reverse over clusters
03319     }
03320   }
03321   //outside of the loop, but there could still be an extra best kp to add
03322   if(bestkp){
03323     nadd += AddReverseBestKPToFit(bestkp, oldkp, cfh, newkplist,
03324                                   nswimfail, direction);
03325     MSG("FitTrackSR" , Msg::kDebug) << "bestkp has direction " 
03326                                     << bestkp->GetZDir() << endl;
03327     bestkp = 0;
03328   }
03329   
03330   return nadd;
03331 }

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(), Msg::kDebug, kQoverP, MSG, and Predict().

Referenced by AddReverseBestKPToFit(), AddToFit(), FindDownstreamPlanes(), FindUpstreamPlanes(), and ReverseFit().

00900 {
00901   
00902 
00903 //
00904 // < 0 = failure
00905 // return value =  0        success
00906 //
00907 
00908   Int_t nswimfail = 0;
00909 MSG("FitTrackSR",Msg::kDebug) << "in FitFrom: " 
00910                                 << "  Fitting direction " << idir << " " 
00911                                 << currentkp->GetTrackCluster()->GetPlane() << " " 
00912                                 << currentkp->GetTrackCluster()->GetTPos() 
00913                                 << " min/max strip = " 
00914                                 << currentkp->GetTrackCluster()->GetMinStrip() << "/" 
00915                                 << currentkp->GetTrackCluster()->GetMaxStrip() 
00916                                 << currentkp->GetZDir() << endl;
00917  
00918 
00919  //predict the state of the plane to add to the fit based on the 
00920   //previously fit plane
00921   nswimfail = Predict(prevkp, currentkp, idir);
00922 
00923   if(nswimfail<0){
00924     return nswimfail;
00925   }
00926 
00927   Int_t istatus = 0;
00928   istatus = Filter(currentkp, idir);
00929 
00930   MSG("FitTrackSR",Msg::kDebug)  << "fit point " << idir << " "  << " " 
00931                                  << currentkp->GetTrackCluster()->GetPlane() << " " 
00932                                  << currentkp->GetTrackCluster()->GetTPos() << endl;
00933 
00934   MSG("FitTrackSR",Msg::kDebug) << "fit qp = " 
00935                                 << currentkp->GetFilStateValue(kQoverP,idir) 
00936                                 << " sigma q/p= "
00937                                 << currentkp->GetFilCovarianceMatrixValue(0,
00938                                                                           kQoverP,kQoverP)
00939                                 << "/" 
00940                                 << currentkp->GetFilCovarianceMatrixValue(1,
00941                                                                           kQoverP,kQoverP)
00942                                 << " range (g/cm**2) = " << currentkp->GetRange() 
00943                                 << " maximum qp from range = " << fParmMaxQPFrac 
00944                                 << "/0.002/" << currentkp->GetRange() << " = " 
00945                                 << (currentkp->GetRange()>0. 
00946                                     ? fParmMaxQPFrac/0.002/currentkp->GetRange() : 0.) 
00947                                 << " maxqpparm = " << fParmMaxQP << endl;
00948 
00949   return nswimfail;
00950 }

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

Definition at line 520 of file AlgFitTrackSR.cxx.

References fCov, fParmCovarianceScale, 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(), Msg::kDebug, 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().

00521 {
00522 
00523   KalmanPlaneSR *kp = 0;
00524   for(Int_t i=0; i<=cfh.GetKalmanLast(); ++i){
00525     kp = cfh.GetKalmanPlane(i);
00526     kp->SetPredictPlane(-999,idir);
00527   }
00528 
00529   KalmanPlaneSR *currentkp = cfh.GetCurrentKalmanPlane();
00530   //fill up the covariance matrix for the current kp
00531   for(Int_t i=0; i<5; ++i){
00532     for(Int_t j=0; j<5; ++j){
00533       currentkp->SetFilCovarianceMatrixValue(idir,i,j,0.);
00534     }
00535   }
00536 
00537   currentkp->SetFilCovarianceMatrixValue(idir,kU,kU,fParmInitialPositionError2);
00538   currentkp->SetFilCovarianceMatrixValue(idir,kV,kV,fParmInitialPositionError2);
00539   currentkp->SetFilCovarianceMatrixValue(idir,kdUdZ,kdUdZ,fParmInitialSlopeError2);
00540   currentkp->SetFilCovarianceMatrixValue(idir,kdVdZ,kdVdZ,fParmInitialSlopeError2);
00541   currentkp->SetFilCovarianceMatrixValue(idir,kQoverP,kQoverP,fParmInitialQPError2);
00542 
00543   Double_t cov_xqp = currentkp->GetFilCovarianceMatrixValue(idir,kQoverP,kQoverP)*0.01*0.3*0.025;
00544   cov_xqp *= (Double_t)(currentkp->GetZDir());
00545   
00546   currentkp->SetFilCovarianceMatrixValue(idir,kU,kQoverP,cov_xqp);
00547   currentkp->SetFilCovarianceMatrixValue(idir,kV,kQoverP,cov_xqp);
00548   currentkp->SetFilCovarianceMatrixValue(idir,kQoverP,kU,cov_xqp);
00549   currentkp->SetFilCovarianceMatrixValue(idir,kQoverP,kV,cov_xqp);
00550   
00551   //set the filtered state vector with either the vertex or end 
00552   //point values depending on whether the fit is in the forward or
00553   //reverse direction
00554   if(idir==0){
00555     // we begin with the vertex position
00556     currentkp->SetFilStateValue(kU,0,cfh.GetVtxU());
00557     currentkp->SetFilStateValue(kV,0,cfh.GetVtxV());
00558     currentkp->SetFilStateValue(kdUdZ,0,cfh.GetDirCosU()/cfh.GetDirCosZ());
00559     currentkp->SetFilStateValue(kdVdZ,0,cfh.GetDirCosV()/cfh.GetDirCosZ());
00560     currentkp->SetFilStateValue(kQoverP,0,cfh.GetInitialQP());
00561   }
00562   else{
00563     currentkp->SetFilStateValue(kU,1,cfh.GetEndU());
00564     currentkp->SetFilStateValue(kV,1,cfh.GetEndV());
00565     currentkp->SetFilStateValue(kdUdZ,1,cfh.GetEndDirCosU()/cfh.GetEndDirCosZ());
00566     currentkp->SetFilStateValue(kdVdZ,1,cfh.GetEndDirCosV()/cfh.GetEndDirCosZ());
00567     currentkp->SetFilStateValue(kQoverP,1,cfh.GetEndQP());
00568   }
00569     
00570   MSG("FitTrackSR",Msg::kDebug) << "Initializing fit, direction = " << idir 
00571                                 << " iteration = " << iterate << " plane = "
00572                                 << currentkp->GetTrackCluster()->GetPlane() << " " 
00573                                 << currentkp->GetFilStateValue(kU,idir) << " " 
00574                                 << currentkp->GetFilStateValue(kV,idir) << " " 
00575                                 << currentkp->GetFilStateValue(kdUdZ,idir) << " " 
00576                                 << currentkp->GetFilStateValue(kdVdZ,idir) << " " 
00577                                 << currentkp->GetFilStateValue(kQoverP,idir) << endl;
00578 
00579   if(!(idir==0 && iterate==0) && iterate>=0){
00580     fCov[kU][kU] *= fParmCovarianceScale;
00581     fCov[kV][kV] *= fParmCovarianceScale;
00582     fCov[kdUdZ][kdUdZ] *=fParmCovarianceScale;
00583     fCov[kdVdZ][kdVdZ] *= fParmCovarianceScale;
00584     fCov[kQoverP][kQoverP] *=fParmCovarianceScale;
00585     
00586     // make sure none of these are larger than our very first covariance elements
00587     fCov[kU][kU] = min(fCov[kU][kU],fParmInitialPositionError2);
00588     fCov[kV][kV] = min(fCov[kV][kV],fParmInitialPositionError2);
00589     fCov[kdUdZ][kdUdZ] = min(fCov[kdUdZ][kdUdZ],fParmInitialSlopeError2);
00590     fCov[kdVdZ][kdVdZ] = min(fCov[kdVdZ][kdVdZ],fParmInitialSlopeError2);
00591     fCov[kQoverP][kQoverP] = min(fCov[kQoverP][kQoverP],fParmInitialQPError2);
00592 
00593     Double_t cov_xqp = fParmInitialQPError2*0.01*0.3*0.025;
00594     for(Int_t i=0; i<2; ++i){
00595       if(TMath::Abs(fCov[i][kQoverP])>cov_xqp) fCov[i][kQoverP] = (fCov[i][kQoverP] > 0 ? cov_xqp : -cov_xqp);
00596       if(TMath::Abs(fCov[kQoverP][i])>cov_xqp) fCov[kQoverP][i] = (fCov[kQoverP][i] > 0 ? cov_xqp : -cov_xqp);
00597     }//end loop over covariance stuff
00598 
00599     cov_xqp /= 0.06;
00600     for(Int_t i=2; i<4; ++i){
00601       if(TMath::Abs(fCov[i][kQoverP])>cov_xqp) fCov[i][kQoverP] = (fCov[i][kQoverP] > 0 ? cov_xqp : -cov_xqp);
00602       if(TMath::Abs(fCov[kQoverP][i])>cov_xqp) fCov[kQoverP][i] = (fCov[kQoverP][i] > 0 ? cov_xqp : -cov_xqp);
00603     }//end loop to fill covariance matrix stuff
00604 
00605     //set the filtered covariance matrix of the current kp to be the same as fCov
00606     currentkp->SetFilCovarianceMatrix(idir,fCov);
00607 
00608   }//end if !(idir==0 && iterate==0) && iterate>=0
00609   
00610   MsgFormat ffmt("%10.5f");
00611   MSG("FitTrackSR",Msg::kDebug) << "Covariance Matrix\n";
00612   for (Int_t i=0; i<5; i++) {
00613     MSG("FitTrackSR",Msg::kDebug) << ffmt(currentkp->GetFilCovarianceMatrixValue(idir,i,kU))
00614                                  << ffmt(currentkp->GetFilCovarianceMatrixValue(idir,i,kV))
00615                                  << ffmt(currentkp->GetFilCovarianceMatrixValue(idir,i,kdUdZ))
00616                                  << ffmt(currentkp->GetFilCovarianceMatrixValue(idir,i,kdVdZ))
00617                                  << ffmt(currentkp->GetFilCovarianceMatrixValue(idir,i,kQoverP))
00618                                  << endl;
00619   }
00620 
00621   if(idir==0){
00622     cfh.SetVtxZ(currentkp->GetTrackCluster()->GetZPos());
00623     cfh.SetVtxPlane(currentkp->GetTrackCluster()->GetPlane());
00624   }
00625 }

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(), Msg::kDebug, kdUdZ, kdVdZ, kQoverP, kU, kV, and MSG.

Referenced by RunAlg().

00498 {
00499 
00500   fParmDState[kU] = ac.GetDouble("KalmanDState1");
00501   fParmDState[kV] = ac.GetDouble("KalmanDState2");
00502   fParmDState[kdUdZ] = ac.GetDouble("KalmanDState3");
00503   fParmDState[kdVdZ] = ac.GetDouble("KalmanDState4");
00504   fParmDState[kQoverP] = ac.GetDouble("KalmanDState5");
00505   fParmPlnRadLen = ac.GetDouble("KalmanPlnRadLen");
00506   fParmdedx = ac.GetDouble("Kalmandedx");
00507   fParmTPosError2 = ac.GetDouble("TPosError2");
00508   fParmMaxQP = ac.GetDouble("KalmanMaxQP");
00509   fParmQPRangeCheck = ac.GetInt("KalmanQPRangeCheck");
00510   fParmMaxQPFrac = ac.GetDouble("KalmanMaxQPFrac");
00511   fParmMaxAngle = ac.GetDouble("KalmanMaxAngle");
00512   fSwimmer = ac.GetInt("Swimmer");
00513  MSG("FitTrackSR", Msg::kDebug) << "initial state = " << fParmDState[kU] << " " 
00514                                  << fParmDState[kV] << " " << fParmDState[kdUdZ] << " " 
00515                                  << fParmDState[kdVdZ] << " " << fParmDState[kQoverP] << endl;
00516   return;
00517 }  

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::kDebug, MSG, ResetTrackClusterList(), and ReverseFit().

Referenced by DoKalmanFit().

02406 {
02407   MSG("AlgFitTrackSR", Msg::kDebug) << " in IterateKalmanFit" << endl;
02408 
02409   Int_t iterate=0;
02410   Double_t qp0=0.;
02411   Double_t qp1=0.;
02412 
02413   nu = 0;
02414   nv = 0;
02415 
02416   Int_t istatus = 0;
02417  
02418   // iterate until stable
02419  MSG("FitTrackSR",Msg::kDebug) << "Before iteration " << iterate << "/"
02420                                << fParmMaxIterate << " changed points " 
02421                                << cfh.GetNChangedFitPoint() << " qp0 "
02422                                << qp0 << " qp1 " << qp1 
02423                                << " istatus " << istatus 
02424                                 << endl;
02425   Double_t qpDiff = fParmQPDiff+.01;
02426   while(iterate<fParmMaxIterate && !istatus 
02427         && (iterate<=1 
02428             || ((qp1!=0. && TMath::Abs(qpDiff)>fParmQPDiff) 
02429                 && !(qp0==0. && qp1==0.)) 
02430             || cfh.GetNChangedFitPoint()) 
02431         ){
02432     
02433     //get the fit q/p for the forward fit
02434     if(iterate) qp0 = cfh.GetKalmanPlane(0)->GetFilStateValue(4,1);
02435     
02436     cfh.ClearKalmanPlaneList();
02437     nu=0;
02438     nv=0;
02439   
02440     //    clock_t dummyt;
02441     // struct tms t1;
02442     // struct tms t2;
02443     // struct tms t3;
02444     //static double ticksPerSecond = sysconf(_SC_CLK_TCK);
02445     // dummyt = times(&t1);
02446 
02447     istatus = AddClustersToFit(planeClusterList, cfh, iterate, nu, nv, direction); 
02448     // dummyt = times(&t2);
02449     //  cout << " forward fitter time " << (Double_t)(t2.tms_utime+t2.tms_stime-t1.tms_utime-t1.tms_stime)/ticksPerSecond << endl;
02450     // cout << " --------------- do reverse fit -------------------" << endl;
02451 
02452 
02453     //check to see if you have enough planes in both views to reverse the 
02454     //fit
02455     if(nu>=3 && nv>=3) istatus = ReverseFit(planeClusterList,cfh,iterate,dosearch);
02456     // dummyt = times(&t3);
02457     //  cout << " reverse fitter time " << (Double_t)(t3.tms_utime+t3.tms_stime-t2.tms_utime-t2.tms_stime)/ticksPerSecond << endl;
02458     //since new track clusters may have been chosen in the ReverseFit method, 
02459     //reset the cluster list
02460     ResetTrackClusterList(planeClusterList, cfh);
02461     
02462     qp1 = cfh.GetKalmanPlane(0)->GetFilStateValue(4,1);
02463     if(qp1==0)  qp1 = cfh.GetKalmanPlane(0)->GetFilStateValue(4,0);
02464     if(qp1 != 0.) qpDiff = 1.-qp0/qp1;
02465     else qpDiff = fParmQPDiff;
02466     MSG("FitTrackSR",Msg::kDebug) << "After iteration " << iterate << "/"
02467                                  << fParmMaxIterate << " changed points " 
02468                                  << cfh.GetNChangedFitPoint() << " qp0 "
02469                                  << qp0 << " qp1 " << qp1 
02470                                  << " istatus " << istatus 
02471                                  << " qpdiff " << TMath::Abs(qpDiff) 
02472                                  << "/" << fParmQPDiff << endl;
02473 
02474   
02475     ++iterate;
02476   }//end loop of iterations
02477 
02478   return istatus;
02479 }

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

Definition at line 3538 of file AlgFitTrackSR.cxx.

References CandHandle::AddDaughterLink(), fParmIsCosmic, CandStripHandle::GetCharge(), CandFitTrackSRHandle::GetdUdZ(), CandFitTrackSRHandle::GetdVdZ(), CandFitTrackSRHandle::GetKalmanLast(), CandFitTrackSRHandle::GetKalmanPlane(), TrackClusterSR::GetPlane(), CandStripHandle::GetPlane(), TrackClusterSR::GetPlaneView(), CandStripHandle::GetStrip(), TrackClusterSR::GetStripList(), CandStripHandle::GetTPos(), KalmanPlaneSR::GetTrackCluster(), CandTrackHandle::GetU(), CandTrackHandle::GetV(), CandStripHandle::GetZPos(), Msg::kDebug, PlaneView::kU, PlaneView::kV, MSG, CandTrackHandle::SetInShower(), and tc.

Referenced by RunAlg().

03539 {
03540 
03541   Double_t xz = 0.;
03542   Double_t dxdz = 0.;
03543   Int_t nstripexp = 0;
03544   KalmanPlaneSR *kp = 0;
03545   TrackClusterSR *tc = 0;
03546   CandStripHandle *strip = 0;
03547   const CandStripHandle *constStrip = 0;
03548   Int_t minstripindx = 0;
03549   Double_t mintpos = 0;
03550   Int_t indxbest = 0;
03551   Double_t distbest = 0;
03552   Double_t distsum = 0.;
03553   Int_t strip1=-1,strip2=-1;
03554   Double_t dist = 0.;
03555   Int_t nshower = 0;
03556   Int_t thisplane = 0;
03557 
03558   for(Int_t i=0; i<=cfh.GetKalmanLast(); ++i){
03559     kp = cfh.GetKalmanPlane(i);
03560     tc = kp->GetTrackCluster();
03561     
03562     xz = 0.;
03563     dxdz = 0.;
03564     thisplane = tc->GetPlane();
03565     
03566     if(tc->GetPlaneView()==PlaneView::kU){
03567       xz = cfh.GetU(thisplane);
03568       dxdz = cfh.GetdUdZ(thisplane);
03569     } 
03570     else if(tc->GetPlaneView()==PlaneView::kV){
03571       xz = cfh.GetV(thisplane);
03572       dxdz = cfh.GetdVdZ(thisplane);
03573     }
03574     
03575     //see how many strips you expect to be hit in this plane
03576     //assume 4:1 aspect ratio for scintillator strips to planes
03577     nstripexp = TMath::Max(1,(Int_t)(TMath::Abs(dxdz)/4.)+1);
03578       
03579     //if the number of strips in the cluster is less than the expected number,
03580     //take them all  for the daughter list
03581     if(tc->GetStripList()->GetLast()+1<=nstripexp){
03582       
03583       for(Int_t istrip=0; istrip<=tc->GetStripList()->GetLast(); ++istrip){
03584         
03585           strip = dynamic_cast<CandStripHandle*>(tc->GetStripList()->At(istrip));
03586           cfh.AddDaughterLink(*strip);
03587      
03588                   MSG("AlgFitTrackSR", Msg::kDebug) << "in track plane = " << strip->GetPlane() << "/" << strip->GetZPos()
03589                           << " strip = " << strip->GetStrip() << "/" << strip->GetTPos() 
03590                           << " charge = " << strip->GetCharge() << " nstripexp " << nstripexp 
03591                           << " nshower = " << nshower << endl;
03592                   
03593           //not in a shower because you expected that many strips
03594           cfh.SetInShower(strip,0);
03595       }
03596     }//end if fewer strips than expected in cluster
03597     else{
03598         
03599       TObjArray orderedstriplist;
03600       vector<Int_t> orderedstripindx;
03601       // order strips by tpos
03602       minstripindx = 0;
03603 
03604       //loop over the strips until there are no more to add
03605       while(minstripindx>=0){
03606         mintpos = 9999999.;
03607         minstripindx = -1;
03608         for (Int_t istrip=0; istrip<=tc->GetStripList()->GetLast(); ++istrip){
03609           constStrip = dynamic_cast<const CandStripHandle*>(tc->GetStripList()->At(istrip));
03610           
03611           if(orderedstriplist.IndexOf(constStrip)<0 && constStrip->GetTPos()<mintpos){
03612             mintpos = constStrip->GetTPos();
03613             minstripindx = istrip;
03614           }
03615         }//end loop over strips in cluster
03616         if(minstripindx>=0){
03617           orderedstripindx.push_back(minstripindx);
03618           strip = dynamic_cast<CandStripHandle*>(tc->GetStripList()->At(minstripindx));
03619           orderedstriplist.Add(strip);
03620         }
03621       }//end loop over strips to add to list
03622       
03623       indxbest = 0;
03624       distbest = 0;
03625       
03626       for(Int_t indx=0; indx<=tc->GetStripList()->GetLast()+1-nstripexp; ++indx){
03627         distsum = 0.;
03628         strip1 = -1;
03629         strip2 = -1;
03630         for(Int_t istrip=indx; istrip<indx+nstripexp; ++istrip){
03631           constStrip = dynamic_cast<const CandStripHandle*>(tc->GetStripList()->At(orderedstripindx[istrip]));
03632           if (strip1<0) strip1 = constStrip->GetStrip();
03633           strip2 = constStrip->GetStrip();
03634           dist = TMath::Abs(constStrip->GetTPos()-xz);
03635           if (fParmIsCosmic) dist *= (0.3+exp(-0.2*constStrip->GetCharge()));
03636           
03637           distsum += dist;
03638         }
03639 
03640         if (!indx || distsum<distbest) {
03641           distbest = distsum;
03642           indxbest  = indx;
03643         }
03644       }//end loop over strips in cluster
03645     
03646       // count how many hit strips above 1 pe
03647       nshower=0;
03648       for(int i=0; i<=tc->GetStripList()->GetLast(); ++i){
03649         strip = dynamic_cast<CandStripHandle*>(tc->GetStripList()->At(i));
03650         if (strip->GetCharge()>1.5) ++nshower;
03651                 MSG("AlgFitTrackSR", Msg::kDebug) << "in shower plane = " << strip->GetPlane() << "/" << strip->GetZPos()
03652                         << " strip = " << strip->GetStrip() << "/" << strip->GetTPos()
03653                         << " charge = " << strip->GetCharge() << " nstripexp " << nstripexp 
03654                         << " nshower = " << nshower << endl;
03655       }
03656       for(Int_t istrip=indxbest; istrip<indxbest+nstripexp; ++istrip){
03657         strip = dynamic_cast<CandStripHandle*>(tc->GetStripList()->At(orderedstripindx[istrip]));
03658         cfh.AddDaughterLink(*strip);
03659         cfh.SetInShower(strip,nshower-nstripexp);
03660       }
03661     }//end too more strips in cluster than expected
03662   
03663   }//end loop over kps
03664   
03665     return;
03666 }

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::kDebug, MSG, TrackClusterSR::SetIsValid(), and tc.

Referenced by RunAlg().

01723 {
01724  
01725   stripItr.Reset();
01726   
01727   // select track clusters which are in this slice
01728 
01729   Int_t oldplane = 0;
01730   Int_t minstrip = 0;
01731   Int_t maxstrip = 0;
01732 
01733   CandStripHandle *strip = 0;
01734   TrackClusterSR *tc = 0;
01735 
01736   while( (strip = stripItr()) ){
01737           MSG("AlgFitTrackSR", Msg::kDebug) << "plane = " << strip->GetPlane() << "/" << strip->GetZPos()
01738           << " strip = " << strip->GetStrip() << "/" << strip->GetTPos() 
01739           << " charge = " << strip->GetCharge() << endl;
01740           
01741     if(strip->GetPlane()!=oldplane || strip->GetStrip()<minstrip 
01742        || strip->GetStrip()>maxstrip) {
01743        
01744       Bool_t found(0);
01745       for(int i=0; !found && i<=tclist->GetLast(); ++i){
01746         tc = dynamic_cast<TrackClusterSR*>(tclist->At(i));
01747         assert(tc);
01748         if((found = tc->IsContained(strip))){
01749           TrackClusterSR *newtc = new TrackClusterSR(*tc);
01750           newtc->SetIsValid(true);
01751           trackClusterList->AddLast(newtc);
01752           MSG("FitTrackSR",Msg::kDebug) << "Adding trackcluster " 
01753                                            << tc->GetPlane() << " " 
01754                                            << tc->GetMinStrip() 
01755                                            << " " << tc->GetMaxStrip() 
01756                                            << " is valid " 
01757                                            << (Int_t)tc->IsValid() 
01758                                            << " to list" << endl; 
01759           cfh.AddTrackCluster(newtc);
01760           oldplane = tc->GetPlane();
01761           minstrip = tc->GetMinStrip();
01762           maxstrip = tc->GetMaxStrip();
01763         }//end if this cluster contains the current strip
01764       }//end loop over clusters
01765     }//end if this strip is in a new cluster
01766   }//end loop over strips in slice
01767 
01768   return;
01769 }

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, fParmMinClusterCharge, 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(), Msg::kDebug, PlaneView::kU, PlaneView::kV, and MSG.

Referenced by RunAlg().

01780 {
01781  
01782   MSG("AlgFitTrackSR", Msg::kDebug) << "in MakeTrackClusterList" << endl;
01783   stripItr.Reset();
01784 
01785   // create track clusters
01786   Int_t oldplane=0;
01787   Int_t oldstrip=0;
01788   Double_t oldtime=0.;
01789   TrackClusterSR *tcluster=0;
01790   CandStripHandle *strip = 0;
01791 
01792   while( (strip = stripItr()) ){
01793     if(!strip->GetDemuxVetoFlag() 
01794        && (strip->GetPlaneView()==PlaneView::kU 
01795            || strip->GetPlaneView()==PlaneView::kV)){
01796         
01797                 MSG("AlgFitTrackSR", Msg::kDebug) << "plane = " << strip->GetPlane() << "/" << strip->GetZPos()
01798                 << " strip = " << strip->GetStrip() << "/" << strip->GetTPos() 
01799                 << " charge = " << strip->GetCharge() << endl;
01800                 
01801       if(tcluster && strip->GetPlane()==oldplane 
01802          && direction*strip->GetStrip()<=direction*oldstrip+3){
01803         tcluster->AddStrip(strip);
01804         oldstrip = strip->GetStrip();
01805         oldtime = strip->GetBegTime();
01806       }
01807       else{
01808         TrackClusterSR *trackcluster = new TrackClusterSR(strip,
01809                                                 fParmMisalignmentError);
01810         trackClusterList->AddLast(trackcluster);
01811         oldplane = strip->GetPlane();
01812         oldstrip = strip->GetStrip();
01813         oldtime = strip->GetBegTime();
01814         tcluster = trackcluster;
01815       }//end if this is a new cluster
01816     }//end if strip is in the detector and not vetoed by DeMux 
01817   }//end loop over strips 
01818         
01819  
01820   for(Int_t i=0; i<=trackClusterList->GetLast(); ++i){
01821    tcluster = dynamic_cast<TrackClusterSR*>(trackClusterList->At(i));
01822 
01823    //get rid of clusters that are too wide or too little charge
01824    if(tcluster->GetNStrip()>fParmMaxClusterNStrip 
01825       || tcluster->GetCharge()<fParmMinClusterCharge){
01826      trackClusterList->RemoveAt(i);
01827      delete tcluster;
01828    }
01829    else{
01830   MSG("AlgFitTrackSR",Msg::kDebug) << "Adding trackcluster " 
01831                                       << tcluster->GetPlane() 
01832                                       << " " << tcluster->GetMinStrip() << " " 
01833                                       << tcluster->GetMaxStrip() << " is valid " 
01834                                       << (Int_t)tcluster->IsValid() 
01835                                       << "  to list" << endl;
01836      cfh.AddTrackCluster(tcluster);
01837    }   
01838   }//end loop over clusters to remove the bad ones
01839 
01840   trackClusterList->Compress();
01841 
01842   return;
01843 }

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 fParmMaxImpactParameter, TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetPlane(), CandStripHandle::GetPlane(), TrackClusterSR::GetPlaneView(), CandStripHandle::GetStrip(), TrackClusterSR::GetStripList(), CandTrackHandle::GetU(), CandTrackHandle::GetV(), TrackClusterSR::IsContained(), CandTrackHandle::IsInShower(), TrackClusterSR::IsValid(), Msg::kDebug, PlaneView::kU, PlaneView::kV, MSG, ResetTrackClusterList(), CandTrackHandle::SetInShower(), TrackClusterSR::SetIsValid(), and tc.

Referenced by RunAlg().

01981 {
01982 
01983   bool goodTrack = true;
01984 
01985   Int_t oldplane=0;
01986 
01987   stripItr.Reset();
01988 
01989 // add trackclusters that are part of track
01990   Int_t minstrip = 0;
01991   Int_t maxstrip = 0;
01992   map<TrackClusterSR*,Bool_t> isintrack;
01993   
01994   Int_t inTrackPlane[1000];
01995   for(Int_t i = 0; i<1000; ++i){
01996     inTrackPlane[i] = 0;
01997   }
01998 
01999   Int_t ntrku = 0;
02000   Int_t ntrkv = 0;
02001 
02002   Double_t oldu,oldv;
02003   Int_t prevplane=-1;
02004   begPlane = -1;
02005   endPlane = -1;
02006   Bool_t found = false;
02007   TrackClusterSR *tc = 0;
02008   CandStripHandle *strip = 0;
02009   CandStripHandle *tcstrip = 0;
02010 
02011   //set all the track clusters in the iterator to IsValid = false
02012   ResetTrackClusterList(planeClusterList);
02013 
02014   //loop over strips to find the clusters in the track and mark them as valid
02015   while( (strip = stripItr()) ){
02016     // check if strip outside detector
02017     if(prevplane!=strip->GetPlane()){
02018       MSG("FitTrackSR", Msg::kDebug) << "prev plane = " << prevplane 
02019                                    << " strip plane = "
02020                                     << strip->GetPlane() << " strip " 
02021                                     << strip->GetStrip()
02022                                     << endl;
02023 
02024      prevplane = strip->GetPlane();
02025      oldu = track0->GetU(prevplane);
02026      oldv = track0->GetV(prevplane);
02027      MSG("FitTrackSR",Msg::kDebug) << "PLANE " << prevplane 
02028                                    << " u,v = " << oldu << " " << oldv << endl; 
02029     }
02030     //if the current track point is out of the detector, ignore it
02031     //B.J.R. 3.24.2005: test the impact parameter for the point to make the call
02032     if(TMath::Sqrt(oldu*oldu + oldv*oldv) > fParmMaxImpactParameter){
02033       MSG("FitTrackSR",Msg::kDebug) << "PLANE " << prevplane 
02034                                     << " out of detector, u,v = " 
02035                                     << oldu << " " << oldv 
02036                                     << " REMOVING from fittrack input" << endl;  
02037       continue;
02038     }
02039 
02040     //make sure you only take one cluster from each plane
02041     if(inTrackPlane[strip->GetPlane()]==0 
02042        && (strip->GetPlane()!=oldplane || strip->GetStrip()<minstrip 
02043            || strip->GetStrip()>maxstrip)){
02044       found = false;
02045       
02046       //loop over the clusters until you find the one containing the 
02047       //current strip
02048 
02049       TObjArray * planeClusters=(TObjArray *)planeClusterList.At(strip->GetPlane());
02050       TIter clusterItr(planeClusters);
02051       while(!found && (tc = (TrackClusterSR *)(clusterItr.Next()) )){
02052         if( tc->GetPlane()==strip->GetPlane()){
02053           
02054           MSG("AlgFitTrackSR", Msg::kDebug) << "cluster on plane " << tc->GetPlane()
02055                                          << " min strip " << tc->GetMinStrip() 
02056                                        << " max strip " << tc->GetMaxStrip() 
02057                                        << " is valid " << (Int_t)tc->IsValid()
02058                                        << endl;           //if this strip is in this cluster
02059           if(tc->IsContained(strip)){
02060             found = true;           
02061             //the cluster is good so set the plane flag to 1 and the cluster
02062             //valid flag to true
02063             inTrackPlane[strip->GetPlane()] = 1;
02064             tc->SetIsValid(true);
02065             MSG("AlgFitTrackSR",Msg::kDebug) << "Adding trackcluster " 
02066                                              << tc->GetPlane() 
02067                                              << " " << tc->GetMinStrip() << " " 
02068                                              << tc->GetMaxStrip() << " is valid " 
02069                                              << (Int_t)tc->IsValid() 
02070                                              << " to track list" 
02071                                              << endl;
02072             //see what view the cluster is in and if it is a good beg plane
02073             if(tc->GetPlaneView()==PlaneView::kU) ++ntrku;
02074             else if(tc->GetPlaneView()==PlaneView::kV) ++ntrkv;
02075             if(begPlane<0) begPlane = strip->GetPlane();
02076             endPlane = strip->GetPlane();
02077 
02078             //loop over all strips in the cluster and add them to the CandFitTrack
02079             for(int j=0; j<=tc->GetStripList()->GetLast(); ++j){
02080               tcstrip = dynamic_cast<CandStripHandle*>(tc->GetStripList()->At(j));
02081               cfh.SetInShower(tcstrip,track0->IsInShower(strip));
02082             }
02083             oldplane = tc->GetPlane();
02084             minstrip = tc->GetMinStrip();
02085             maxstrip = tc->GetMaxStrip();
02086           }//end if this strip is in the current cluster
02087         }//end loop over clusters
02088       }
02089 
02090            
02091     }//end if strip is from a different plane/cluster
02092   }//end loop over strips
02093 
02094   stripItr.Reset();
02095   MSG("FitTrackSR",Msg::kDebug) << "# of u,v planes in track = " << ntrku << " " 
02096                                 << ntrkv << endl;
02097  
02098  if(ntrku<3 || ntrkv<3){
02099     goodTrack = false;
02100     return goodTrack;
02101   } 
02102 
02103   return goodTrack;
02104 }

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, fParmMaxQP, fParmMaxQPFrac, fParmQPRangeCheck, TrackClusterSR::GetPlane(), KalmanPlaneSR::GetPredictPlane(), KalmanPlaneSR::GetPreStateValue(), KalmanPlaneSR::GetRange(), TrackClusterSR::GetTPosError(), KalmanPlaneSR::GetTrackCluster(), Msg::kDebug, kdUdZ, kdVdZ, kQoverP, MSG, KalmanPlaneSR::SetPredictPlane(), and KalmanPlaneSR::SetPreStateValue().

Referenced by FitFrom().

00750 {
00751 
00752 //
00753 // < 0 indicates failure
00754 // return value =  0        success
00755 // > 0 indicates number of swim failures
00756 //
00757   Int_t nswimfail=0;
00758   MSG("FitTrackSR", Msg::kDebug) << "in Predict" << endl;
00759   //make sure you are looking at different planes
00760   if(prevkp->GetTrackCluster()->GetPlane()!=currentkp->GetPredictPlane(idir)) {
00761     
00762     Int_t swimstatus = CalculatePropagator(prevkp,currentkp,idir);
00763     
00764     if(swimstatus<0){
00765     MSG("FitTrackSR",Msg::kDebug) << "failed to calculate propagator" << endl;
00766       return swimstatus;
00767     }
00768 
00769     nswimfail += swimstatus;
00770 
00771     //get the Q_{k-1} matrix
00772     CalculateNoise(prevkp, currentkp, idir); 
00773 
00774     //get the C_{k}^{k-1} matrix
00775     CalculatePreCovariance(prevkp, currentkp, idir); 
00776 
00777     //extrapolate from the previously fit plane to the plane to add
00778     //to the fit using the swimmer - this takes the place of the F_{k-1}x_{k-1}
00779     //matrix multiplication in the Kalman Filtering equations.
00780     swimstatus = CalculatePreState(prevkp, currentkp, idir);
00781 
00782     if(swimstatus<0){
00783     MSG("FitTrackSR",Msg::kDebug) << "failed to calculate predicted state" << endl;
00784       return swimstatus;
00785     }
00786 
00787     nswimfail += swimstatus;
00788   }//end if on a different plane than the predict plane
00789 
00790   CalculatePreChi2(currentkp, idir,currentkp->GetTrackCluster()->GetTPosError());
00791 
00792   //set the plane you predicted the fit for the current one from
00793   currentkp->SetPredictPlane(prevkp->GetTrackCluster()->GetPlane(), idir);
00794 
00795   //i am not sure if the following lines are needed or not - they are done in
00796   //the CalculatePreState method, but that is only called if the plane to add
00797   //to the fit is different plane than the previously fit plane - dunno, guess
00798   //i will leave them in
00799 
00800   //---------------------------------------------------------------------------//
00801   if(TMath::Abs(currentkp->GetPreStateValue(kdUdZ, idir))>fParmMaxAngle){
00802     if(currentkp->GetPreStateValue(kdUdZ, idir)>0.) 
00803       currentkp->SetPreStateValue(kdUdZ, idir, fParmMaxAngle);
00804     if(currentkp->GetPreStateValue(kdUdZ, idir)<0.) 
00805       currentkp->SetPreStateValue(kdUdZ, idir, -fParmMaxAngle);
00806   }
00807 
00808   if(TMath::Abs(currentkp->GetPreStateValue(kdVdZ, idir))>fParmMaxAngle){
00809     if(currentkp->GetPreStateValue(kdVdZ, idir)>0.) 
00810       currentkp->SetPreStateValue(kdVdZ, idir, fParmMaxAngle);
00811     if(currentkp->GetPreStateValue(kdVdZ, idir)<0.) 
00812       currentkp->SetPreStateValue(kdVdZ, idir, -fParmMaxAngle);
00813   }
00814   
00815   Double_t maxqp = fParmMaxQP;
00816 
00817   //use a value of about 2 MeV/(g/cm^2) for the energy loss in the material
00818   //you would divide the max qp frac allowed by that value and then by the range to
00819   //get the max qp.  but multiplication is a faster operation, so multiply by
00820   //500 instead.
00821 
00822   if ( fParmQPRangeCheck ) {
00823       if(currentkp->GetRange()>0. && fParmMaxQPFrac*500./currentkp->GetRange()<maxqp) {
00824           maxqp = fParmMaxQPFrac*500./currentkp->GetRange();
00825       }
00826   }
00827     
00828   if(idir==0 && TMath::Abs(currentkp->GetPreStateValue(kQoverP, idir))>maxqp){
00829     if(currentkp->GetPreStateValue(kQoverP, idir)>0.) 
00830       currentkp->SetPreStateValue(kQoverP, idir, maxqp);
00831     else currentkp->SetPreStateValue(kQoverP, idir, -maxqp);
00832   }
00833   //---------------------------------------------------------------------------//
00834 
00835   return nswimfail;
00836 }

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, fParmMaxAngleCovariance, fParmMaxLocalChi2, fParmNSkipView, KalmanPlaneSR::GetFilChi2(), KalmanPlaneSR::GetFilCovarianceMatrixValue(), CandFitTrackSRHandle::GetKalmanLast(), CandFitTrackSRHandle::GetKalmanPlane(), TrackClusterSR::GetMaxStrip(), TrackClusterSR::GetMinStrip(), TrackClusterSR::GetPlane(), TrackClusterSR::GetPlaneView(), KalmanPlaneSR::GetTrackCluster(), TrackClusterSR::IsValid(), Msg::kDebug, kdUdZ, kdVdZ, PlaneView::kU, PlaneView::kV, max, min, MSG, ResetTrackClusterList(), TrackClusterSR::SetIsValid(), CandFitTrackSRHandle::SetVtx(), and tc.

Referenced by DoKalmanFit().

02726 {
02727 
02728   MSG("FitTrackSR", Msg::kDebug) << "in RemoveBadPointsFromFit" << endl; 
02729 
02730   TrackClusterSR *tc = 0;
02731   const KalmanPlaneSR *vtxkp = 0;
02732   Int_t nremoveitr = 0;
02733   Int_t nbadfit = 0;
02734   nfitu = 0;
02735   nfitv = 0;
02736   
02737   Double_t maxlocalchi2 = fParmMaxLocalChi2;
02738   Double_t maxanglecovariance = fParmMaxAngleCovariance;
02739   
02740   // if chi2 has to be increased due to too many points being removed, 
02741   //deltachi2 > 0
02742   fParmDeltaChi2 = -1.; 
02743 
02744   // if angle covariance has to be increased due to too many points 
02745   //being removed, deltacovariance > 0
02746   fParmDeltaCovariance = -1.; 
02747   
02748   Double_t mindchi2=0.;
02749   Double_t mindcovariance=0.;
02750   
02751   Int_t nPlanesSkippedU = 0;
02752   Int_t nPlanesSkippedV = 0;
02753   Int_t lastUPlane = -1;
02754   Int_t lastVPlane = -1;
02755   Int_t lastFitPlane = -1;
02756 
02757   MSG("FitTrackSR", Msg::kDebug) << "nfit u " << nfitu << " nfitv " << nfitv << endl;
02758   while(nfitu<3 || nfitv<3){
02759     
02760     vtxkp = 0;
02761     nfitu = 0;
02762     nfitv = 0;
02763     nbadfit = 0;
02764     
02765     //reset the cluster list
02766     ResetTrackClusterList(planeClusterList);
02767     
02768     mindchi2=0.;
02769     mindcovariance=0.;
02770     
02771     nPlanesSkippedU = 0;
02772     nPlanesSkippedV = 0;
02773     
02774     lastUPlane = -1;
02775     lastVPlane = -1;
02776     
02777     KalmanPlaneSR *kp = 0;
02778     
02779     //loop over the Kalman planes to fill the maps
02780     for(Int_t i=0; i<=cfh.GetKalmanLast(); ++i){
02781       kp = cfh.GetKalmanPlane(i);
02782       
02783       //see if the reset of the clusterList also reset the 
02784       //clusters in the kp's
02785       MSG("FitTrackSR", Msg::kDebug) << "cluster for kp on plane "
02786                                      << kp->GetTrackCluster()->GetPlane()
02787                                      << " is valid " 
02788                                      << (Int_t)kp->GetTrackCluster()->IsValid()
02789                                      << endl;
02790       
02791       MSG("FitTrackSR",Msg::kDebug) << "  TC " << kp->GetTrackCluster()->GetPlane() 
02792                                     << " " << kp->GetTrackCluster()->GetMinStrip() 
02793                                     << " " << kp->GetTrackCluster()->GetMaxStrip() 
02794                                     << " " << kp->GetFilChi2(0)
02795                                     << " " << kp->GetFilChi2(1) 
02796                                     << "/" << maxlocalchi2 << " " 
02797                                     << TMath::Min(kp->GetFilCovarianceMatrixValue(0,kdUdZ, kdUdZ),
02798                                                   kp->GetFilCovarianceMatrixValue(1,kdUdZ, kdUdZ))
02799                                     << " " 
02800                                     << TMath::Min(kp->GetFilCovarianceMatrixValue(0,kdVdZ, kdVdZ),
02801                                                   kp->GetFilCovarianceMatrixValue(1,kdVdZ, kdVdZ)) 
02802                                     << endl;
02803       
02804       //see if the fit for this plane was really good
02805       if(TMath::Max(kp->GetFilChi2(0),kp->GetFilChi2(1))<=maxlocalchi2 
02806          && TMath::Min(kp->GetFilCovarianceMatrixValue(0,kdUdZ,kdUdZ),
02807                        kp->GetFilCovarianceMatrixValue(1,kdUdZ,kdUdZ))<maxanglecovariance 
02808          && TMath::Min(kp->GetFilCovarianceMatrixValue(0,kdVdZ,kdVdZ),
02809                        kp->GetFilCovarianceMatrixValue(1,kdVdZ,kdVdZ))<maxanglecovariance)
02810         {
02811           
02812           MSG("FitTrackSR", Msg::kDebug) << "  this fit was good" << endl;
02813           
02814           kp->GetTrackCluster()->SetIsValid(true);
02815           
02816           MSG("FitTrackSR", Msg::kDebug) << "  plane " << kp->GetTrackCluster()->GetPlane()
02817                                          << " IsValid " << kp->GetTrackCluster()->IsValid()
02818                                          << endl;  
02819           TObjArray * planeClusters=(TObjArray *)planeClusterList.At(kp->GetTrackCluster()->GetPlane());
02820           TIter clusterItr(planeClusters);
02821           while((tc = (TrackClusterSR *)(clusterItr.Next()) )){
02822             if(tc->GetPlane()==kp->GetTrackCluster()->GetPlane()){  // replacing slice
02823               if(tc->GetMinStrip() == kp->GetTrackCluster()->GetMinStrip()
02824                  && tc->GetMaxStrip() == kp->GetTrackCluster()->GetMaxStrip() )
02825                 tc->SetIsValid(true);
02826             }
02827           }
02828           
02829           if(kp->GetTrackCluster()->GetPlaneView() == PlaneView::kU)
02830             uFitPlanes[kp->GetTrackCluster()->GetPlane()] = 1;
02831           
02832           else if(kp->GetTrackCluster()->GetPlaneView() == PlaneView::kV)
02833             vFitPlanes[kp->GetTrackCluster()->GetPlane()] = 1;
02834           
02835         }//end if good fit for kp
02836       else{
02837         
02838         MSG("FitTrackSR", Msg::kDebug) << "  plane " << kp->GetTrackCluster()->GetPlane()
02839                                        << " IsValid " 
02840                                        << (Int_t)kp->GetTrackCluster()->IsValid()
02841                                        << endl;
02842         
02843         //no need to set the fTrackCluster to IsValid=false, as the 
02844         //reset call above did that.
02845         if(kp->GetTrackCluster()->GetPlaneView() == PlaneView::kU)
02846           uFitPlanes[kp->GetTrackCluster()->GetPlane()] = 0;
02847         else if(kp->GetTrackCluster()->GetPlaneView() == PlaneView::kV)
02848           vFitPlanes[kp->GetTrackCluster()->GetPlane()] = 0;
02849       }
02850       
02851     }//end loop to fill map of fit quality for kp's
02852     
02853     lastFitPlane = cfh.GetKalmanPlane(cfh.GetKalmanLast())->GetTrackCluster()->GetPlane();
02854     map<Int_t, Int_t>::iterator vfitIter = vFitPlanes.begin();
02855     map<Int_t, Int_t>::iterator ufitIter = uFitPlanes.begin();
02856     for(Int_t i=0; i<=cfh.GetKalmanLast(); ++i){
02857       kp = cfh.GetKalmanPlane(i);
02858       
02859       //look to see which view the kp comes from, and if it has a 
02860       //good fit, as found above
02861       if(kp->GetTrackCluster()->GetPlaneView()==PlaneView::kU){
02862         if(uFitPlanes[kp->GetTrackCluster()->GetPlane()]==1){
02863           ++nfitu;
02864           lastUPlane = kp->GetTrackCluster()->GetPlane();
02865           if(!vtxkp) vtxkp = kp;
02866         }
02867         else{
02868           //see how many planes were skipped in this view
02869           nPlanesSkippedU = FindNumSkippedPlanesInView(uFitPlanes,
02870                                                        ufitIter, 
02871                                                        kp->GetTrackCluster()->GetPlane(), 
02872                                                        lastUPlane,
02873                                                        lastFitPlane, direction);
02874           if(nPlanesSkippedU>=fParmNSkipView){
02875             
02876             //dont remove this plane as there would be too many skipped planes
02877             //in the view if we did
02878             
02879             TObjArray * planeClusters=(TObjArray *)planeClusterList.At(kp->GetTrackCluster()->GetPlane());
02880             TIter clusterItr(planeClusters);
02881             while((tc = (TrackClusterSR * )(clusterItr.Next()) )){
02882               if(tc->GetPlane()==kp->GetTrackCluster()->GetPlane()){
02883                 if(tc->GetMinStrip() == kp->GetTrackCluster()->GetMinStrip()
02884                    && tc->GetMaxStrip() == kp->GetTrackCluster()->GetMaxStrip() )
02885                   tc->SetIsValid(true);
02886               }
02887             }
02888             MSG("FitTrackSR",Msg::kDebug) << "too many consecutive planes  on plane " 
02889                                           << kp->GetTrackCluster()->GetPlane()<<  " # skipped  " << nPlanesSkippedU <<"/" << fParmNSkipView 
02890                                           << "   keeping plane" 
02891                                           << endl;
02892             
02893             if (!vtxkp) vtxkp = kp;
02894             lastUPlane = kp->GetTrackCluster()->GetPlane();
02895             ++nfitu;      
02896           }
02897           else{
02898             MSG("FitTrackSR",Msg::kDebug) << "REMOVING, plane " 
02899                                           << kp->GetTrackCluster()->GetPlane() 
02900                                           << " nskipped = " << nPlanesSkippedU 
02901                                           << "/" << fParmNSkipView << endl;
02902             
02903             ++nbadfit;
02904             if(TMath::Max(kp->GetFilChi2(0),kp->GetFilChi2(1))>maxlocalchi2){
02905               if(mindchi2<=0. 
02906                  || TMath::Max(kp->GetFilChi2(0),kp->GetFilChi2(1))-maxlocalchi2<mindchi2){
02907                 mindchi2 = TMath::Max(kp->GetFilChi2(0),kp->GetFilChi2(1))-maxlocalchi2;
02908                 mindchi2 += 0.001; // precision
02909               }//end if the max chi^2 is bigger than allowed value
02910             } 
02911             else{
02912               Double_t thismaxcov = TMath::Max(TMath::Min(kp->GetFilCovarianceMatrixValue(0,kdUdZ,kdUdZ),
02913                                                           kp->GetFilCovarianceMatrixValue(1,kdUdZ,kdUdZ)),
02914                                                TMath::Min(kp->GetFilCovarianceMatrixValue(0,kdVdZ,kdVdZ),
02915                                                           kp->GetFilCovarianceMatrixValue(1,kdVdZ,kdVdZ)));
02916               if (mindcovariance<=0. 
02917                   || thismaxcov-maxanglecovariance<mindcovariance) {
02918                 mindcovariance = thismaxcov-maxanglecovariance;
02919                 mindcovariance += 0.001; // precision
02920               }//end if mindcovariance is too small
02921             }//end else to check covariance
02922           }//end else remove the current plane
02923           
02924         }//end else the current kp is a bad fit
02925       }//end if in U view
02926       //check the other view
02927       else if(kp->GetTrackCluster()->GetPlaneView()==PlaneView::kV){
02928         if(vFitPlanes[kp->GetTrackCluster()->GetPlane()]==1){
02929           
02930           ++nfitv;
02931           lastVPlane = kp->GetTrackCluster()->GetPlane();
02932           if(!vtxkp) vtxkp = kp;
02933         }
02934         else{
02935           //see how many planes were skipped in this view
02936           nPlanesSkippedV = FindNumSkippedPlanesInView(vFitPlanes,
02937                                                        vfitIter, 
02938                                                        kp->GetTrackCluster()->GetPlane(), 
02939                                                        lastVPlane,
02940                                                        lastFitPlane, direction);
02941           if(nPlanesSkippedV>=fParmNSkipView){
02942             
02943             //dont remove this plane as there would be too many skipped planes
02944             //in the view if we did
02945             //      clusterItr.GetSet()->Slice(kp->GetTrackCluster()->GetPlane());
02946             //in the view if we did
02947 
02948             MSG("FitTrackSR",Msg::kDebug) << "too many consecutive planes  on plane  " 
02949                                           << kp->GetTrackCluster()->GetPlane()<<  " # skipped  " << nPlanesSkippedU <<"/" << fParmNSkipView 
02950                                           << "   keeping plane" 
02951                                           << endl;
02952 
02953             TObjArray * planeClusters=(TObjArray *)planeClusterList.At(kp->GetTrackCluster()->GetPlane());
02954             TIter clusterItr(planeClusters);
02955             while((tc = (TrackClusterSR * )(clusterItr.Next()) )){
02956               if(tc->GetPlane()==kp->GetTrackCluster()->GetPlane()){
02957                 if(tc->GetMinStrip() == kp->GetTrackCluster()->GetMinStrip()
02958                    && tc->GetMaxStrip() == kp->GetTrackCluster()->GetMaxStrip() )
02959                   tc->SetIsValid(true);
02960               }
02961             }
02962             //      clusterItr.GetSet()->Slice();
02963             
02964             if (!vtxkp) vtxkp = kp;
02965             lastVPlane = kp->GetTrackCluster()->GetPlane();
02966             ++nfitv;
02967           }
02968           else{
02969             MSG("FitTrackSR",Msg::kDebug) << "REMOVING, plane " 
02970                                           << kp->GetTrackCluster()->GetPlane() 
02971                                           << " nskipped = " << nPlanesSkippedV
02972                                           << "/" << fParmNSkipView << endl;
02973             ++nbadfit;
02974             if(max(kp->GetFilChi2(0),kp->GetFilChi2(1))>maxlocalchi2){
02975               if(mindchi2<=0. 
02976                  || max(kp->GetFilChi2(0),kp->GetFilChi2(1))-maxlocalchi2<mindchi2){
02977                 mindchi2 = max(kp->GetFilChi2(0),kp->GetFilChi2(1))-maxlocalchi2;
02978                 mindchi2 += 0.001; // precision
02979               }//end if the max chi^2 is bigger than allowed value
02980             } 
02981             else{
02982               Double_t thismaxcov = max(min(kp->GetFilCovarianceMatrixValue(0,kdUdZ,kdUdZ),
02983                                             kp->GetFilCovarianceMatrixValue(1,kdUdZ,kdUdZ)),
02984                                         min(kp->GetFilCovarianceMatrixValue(0,kdVdZ,kdVdZ),
02985                                             kp->GetFilCovarianceMatrixValue(1,kdVdZ,kdVdZ)));
02986               if (mindcovariance<=0. 
02987                   || thismaxcov-maxanglecovariance<mindcovariance) {
02988                 mindcovariance = thismaxcov-maxanglecovariance;
02989                 mindcovariance += 0.001; // precision
02990               }//end if mindcovariance is too small
02991             }//end else to check covariance
02992           }//end else remove the current plane
02993           
02994         }//end else the current kp is a bad fit
02995       }//end if in V view
02996     }//end loop over kp's to remove bad points
02997     
02998     if (nfitu<3 || nfitv<3) {      
02999       // not enough points, increase maximum chi2
03000       
03001       maxlocalchi2 += mindchi2;
03002       maxanglecovariance += mindcovariance;
03003       
03004       MSG("FitTrackSR",Msg::kDebug) << "TOO FEW POINTS LEFT, increasing chi2" 
03005                                     << " threshold to " << maxlocalchi2 
03006                                     << " and angle covariance threshold to " 
03007                                     << maxanglecovariance << endl;
03008       
03009     }//end if less than 3 planes fit in either view
03010     fParmDeltaChi2 = mindchi2;
03011     fParmDeltaCovariance = mindcovariance;
03012     ++nremoveitr;
03013     assert(nremoveitr<=cfh.GetKalmanLast()+2);
03014     // in principle can never loop through more than # of planes in fit
03015   }//end while less than 3 fit planes in at least one view
03016   //make sure you found a vertex kp and set the vertex to it.
03017   assert(vtxkp);
03018   cfh.SetVtx(vtxkp);
03019     
03020   TIter planeItr(&planeClusterList);
03021   TObjArray * plnarray;
03022   while( (plnarray = (TObjArray * )planeItr.Next()) ){
03023     TIter clusterItr(plnarray);
03024     while( (tc = (TrackClusterSR * )clusterItr.Next()) ){
03025       
03026       MSG("FitTrackSR", Msg::kDebug) << "current cluster on plane " 
03027                                      << tc->GetPlane() << " is valid "
03028                                      << (Int_t)tc->IsValid() << " " 
03029                                      << tc->GetMinStrip() << "-" << tc->GetMaxStrip()
03030                                      << endl;
03031     }
03032   }
03033 
03034   return nbadfit;
03035 }

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.

02172 {
02173   //  clusterItr.Reset(); 
02174   TrackClusterSR *tc = 0;
02175   KalmanPlaneSR *kp = 0;
02176 
02177   //loop over the KalmanPlaneSR objects in the fit - if the cluster
02178   //for a given KalmanPlaneSR is the same as the current cluster,
02179   //set IsValid = true
02180   for(Int_t i = 0; i<=cfh.GetKalmanLast(); ++i){
02181     kp = cfh.GetKalmanPlane(i);
02182     
02183     TObjArray * planeClusters=(TObjArray *)planeClusterList.At(kp->GetTrackCluster()->GetPlane());
02184     TIter clusterItr(planeClusters);   
02185     //loop over all the track clusters and set their valid flags to false
02186     while( (tc = (TrackClusterSR * )clusterItr.Next()) ){
02187       if(tc->GetPlane()==kp->GetTrackCluster()->GetPlane()){
02188         tc->SetIsValid(false);  
02189         if(tc->GetMinStrip() == kp->GetTrackCluster()->GetMinStrip()
02190            && tc->GetMaxStrip() == kp->GetTrackCluster()->GetMaxStrip() )
02191           tc->SetIsValid(true);
02192       }
02193     }
02194   }//end loop over KalmanPlaneSR's
02195   return;
02196 }

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

Definition at line 2147 of file AlgFitTrackSR.cxx.

References TrackClusterSR::GetPlane(), TrackClusterSR::IsValid(), Msg::kDebug, MSG, TrackClusterSR::SetIsValid(), and tc.

Referenced by DoKalmanFit(), IterateKalmanFit(), MarkTrackClusters(), and RemoveBadPointsFromFit().

02148 {
02149   //  clusterItr.Reset(); 
02150   TrackClusterSR *tc = 0;
02151   //loop over all the track clusters and set their valid flags to false
02152   TIter planeItr(&planeClusterList);
02153   TObjArray * plnarray;
02154   while( (plnarray = (TObjArray * )planeItr.Next()) ){
02155     TIter clusterItr(plnarray);
02156     while( (tc = (TrackClusterSR *)clusterItr.Next()) ){
02157       tc->SetIsValid(false);
02158       MSG("FitTrackSR", Msg::kDebug) << "cluster on plane " << tc->GetPlane()
02159                                      << " is valid " << (Int_t)tc->IsValid()
02160                                      << endl;
02161     }
02162   }
02163   //  clusterItr.Reset();
02164   return;
02165 }

Int_t AlgFitTrackSR::ReverseFit ( TObjArray &  planeClusterList,
CandFitTrackSRHandle cfh,
Int_t  iterate,
Bool_t  dosearch = 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(), Msg::kDebug, PlaneView::kU, PlaneView::kV, 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().

02485 {
02486 
02487   InitializeFit(cfh,1,iterate);
02488   TObjArray *KalmanPlaneList =cfh.GetPlaneList();
02489   
02490   TObjArray newlist;
02491   newlist.Clear();
02492   newlist.Compress();
02493   
02494   cfh.SetNChangedFitPoint(0);
02495   
02496   TObjArray *tclist = cfh.GetTrackClusterList();
02497 
02498   KalmanPlaneSR *kpold = dynamic_cast<KalmanPlaneSR*>(KalmanPlaneList->At(KalmanPlaneList->GetLast()));
02499   KalmanPlaneSR *lastkp = kpold;
02500   
02501   Bool_t *kpswimfail = new Bool_t[KalmanPlaneList->GetLast()+1];
02502   
02503   for(Int_t i=0 ;i<=KalmanPlaneList->GetLast(); ++i){
02504     kpswimfail[i] = 0;
02505   }
02506 
02507   Int_t nu = 0;
02508   Int_t nv = 0;
02509   Bool_t didswimfail = false;
02510   
02511   Int_t nswimfail = 0;
02512   Int_t nswimfailnewer = 0;
02513 
02514   KalmanPlaneSR *kpnew = 0;
02515   KalmanPlaneSR *kpnewer = 0;
02516   TrackClusterSR *thistc = 0;
02517 
02518   Bool_t ready = true;
02519 
02520   for(Int_t i=KalmanPlaneList->GetLast()-1; i>=0; --i){
02521     kpnew = dynamic_cast<KalmanPlaneSR*>(KalmanPlaneList->At(i));
02522     MSG("FitTrackSR",Msg::kDebug) << "new TC, plane = " 
02523                                   << kpnew->GetTrackCluster()->GetPlane() 
02524                                   << " tpos = " 
02525                                   << kpnew->GetTrackCluster()->GetTPos() 
02526                                   << " minstrip = " 
02527                                   << kpnew->GetTrackCluster()->GetMinStrip() 
02528                                   << " maxstrip = " 
02529                                   << kpnew->GetTrackCluster()->GetMaxStrip()
02530                                    << " q/p = " 
02531                                   << kpnew->GetFilStateValue(4,0)
02532                                   << endl;
02533     
02534   if(TMath::Abs(kpnew->GetFilStateValue(4,0))<=2.0){
02535     
02536     
02537     MSG("FitTrackSR", Msg::kDebug) << "current kp has sufficient momentum to do "
02538                                    << "reverse fit" << endl;
02539     pRev++;
02540     nswimfail = FitFrom(kpold,kpnew,1,iterate);
02541     
02542     if(nswimfail<0){
02543       MSG("FitTrackSR",Msg::kDebug) << "failed swim" << endl;
02544       return nswimfail;
02545     }
02546     
02547     //if(i==KalmanPlaneList->GetLast()-1){
02548     if(ready){
02549       ready = false;
02550       cfh.SetPlaneChi2(kpold->GetTrackCluster()->GetPlane(),kpold->GetFilChi2(0));
02551       cfh.SetPlanePreChi2(kpold->GetTrackCluster()->GetPlane(),kpold->GetPreChi2(0));
02552       cfh.SetPlaneQP(kpold->GetTrackCluster()->GetPlane(),kpold->GetFilStateValue(4,1));
02553     }
02554     MSG("FitTrackSR",Msg::kDebug) << "prechi2 = " << kpnew->GetPreChi2(1) 
02555                                   << "   filchi2 = " 
02556                                   << max(kpnew->GetFilChi2(0),kpnew->GetFilChi2(1)) 
02557                                   << endl; 
02558     
02559     if(dosearch && iterate<2){
02560       kpnewer = 0;
02561       nswimfailnewer = 0;
02562       
02563       MSG("FitTrackSR",Msg::kDebug) << "Searching for other track clusters in plane " 
02564                                     << kpnew->GetTrackCluster()->GetPlane() 
02565                                     << " original min,max strips = " 
02566                                     << kpnew->GetTrackCluster()->GetMinStrip() 
02567                                     << " " 
02568                                     << kpnew->GetTrackCluster()->GetMaxStrip() 
02569                                     << endl;             
02570       //tclist is the list of TrackClusterSR objects for this track
02571       TObjArray * planeClusters=(TObjArray *)planeClusterList.At(kpnew->GetTrackCluster()->GetPlane());
02572       TIter clusterItr(planeClusters);
02573       while((thistc = dynamic_cast<TrackClusterSR *>(clusterItr.Next()) )){
02574         //is this a good tc and from the same plane but a different
02575         //cluster than the one currently used in the plane?
02576         if(!cfh.GetBadFit(thistc)
02577            && thistc->GetMinStrip()!=kpnew->GetTrackCluster()->GetMinStrip()){
02578           //make a new kp based on the alternate track cluster and set 
02579           //its range and z direction
02580           KalmanPlaneSR *thiskp = new KalmanPlaneSR(thistc);
02581           thiskp->SetRange(cfh.GetRange(thiskp->GetTrackCluster()->GetPlane()));
02582           thiskp->SetZDir(kpold->GetZDir());
02583           
02584           //try fitting from the previous kp to the alternate one in
02585           //the reverse direction
02586           pRev2++;
02587           nswimfailnewer = FitFrom(kpold,thiskp, 1,iterate);
02588           MSG("FitTrackSR",Msg::kDebug) << "considering track cluster: "
02589                                         << tclist->GetLast() << " tpos = " 
02590                                         << thistc->GetTPos() 
02591                                         << " minstrip = " 
02592                                         << thistc->GetMinStrip() << "/" 
02593                                         << kpnew->GetTrackCluster()->GetMinStrip() 
02594                                         << "  prechi2 = " 
02595                                         << thiskp->GetPreChi2(1) << "  filchi2 = " 
02596                                         << thiskp->GetFilChi2(1) << endl;         
02597           if(nswimfailnewer>=0 
02598              && ((kpnewer 
02599                   && (nswimfail<0 || thiskp->GetPreChi2(1)<kpnewer->GetPreChi2(1))) 
02600                  || (!kpnewer && thiskp->GetPreChi2(1)<kpnew->GetPreChi2(1)))){
02601             MSG("FitTrackSR",Msg::kDebug) << "  better track cluster found\n";
02602             
02603             if(kpnewer) delete kpnewer;
02604             kpnewer = thiskp;
02605             nswimfail = nswimfailnewer;
02606           }//end if nswimfailnewer>=0 etc 
02607           else delete thiskp;
02608           
02609         }//end if !GetBadFit() etc
02610       }//end loop over track clusters
02611 
02612       //put the better kp in the list
02613       if(kpnewer){
02614         KalmanPlaneList->AddAt(kpnewer,i);
02615         delete kpnew;
02616         kpnew = kpnewer;
02617         cfh.SetNChangedFitPoint(cfh.GetNChangedFitPoint()+1);
02618       }
02619       
02620     }//end if dosearch
02621    
02622     if(nswimfail<0){
02623       MSG("FitTrackSR",Msg::kDebug) << "failed swim" << endl;   
02624       return nswimfail;
02625     }
02626     
02627     if(nswimfail!=0 && nu>=3 && nv>=3){
02628       kpswimfail[i] = 1;
02629       didswimfail = 1;
02630       MSG("FitTrackSR",Msg::kDebug) << "Swim failed, skipping plane" << endl;
02631     } 
02632     else{
02633       //swim didnt fail, set the chi^2 and qp values for the plane
02634       cfh.SetPlaneChi2(kpnew->GetTrackCluster()->GetPlane(),
02635                        max(kpnew->GetFilChi2(1), kpnew->GetFilChi2(0)));
02636       cfh.SetPlanePreChi2(kpnew->GetTrackCluster()->GetPlane(),
02637                           max(kpnew->GetPreChi2(1), kpnew->GetPreChi2(0)));
02638       
02639       cfh.SetPlaneQP(kpnew->GetTrackCluster()->GetPlane(),kpnew->GetFilStateValue(4,1));
02640       cfh.SetCurrentKalmanPlane(kpnew);
02641       lastkp = kpnew;
02642       
02643       //set the covariance matrix array based on the filtered covariance
02644       //matrix of the current kp
02645       for(int i1=0; i1<5; ++i1){
02646         for(int i2=0; i2<5; ++i2){
02647           fCov[i1][i2]=kpnew->GetFilCovarianceMatrixValue(1,i1,i2);
02648         }
02649       }
02650       
02651       if(kpnew->GetTrackCluster()->GetPlaneView()==PlaneView::kU) ++nu;
02652       if(kpnew->GetTrackCluster()->GetPlaneView()==PlaneView::kV) ++nv;
02653       
02654       kpold = kpnew;
02655     }//end else for if nswimfail!=0 && nu>=3 && nv>=3
02656     
02657     cfh.SetNSwimFail(cfh.GetNSwimFail()+TMath::Abs(nswimfail));
02658     
02659   }//end if enough momentum left to do fit
02660   //  else kpold = kpnew;
02661   }//end loop over kalman planes in reverse
02662 
02663   // remove any planes which failed swim
02664   if(didswimfail){
02665     
02666     Int_t nlast = KalmanPlaneList->GetLast();
02667     TObjArray newlist;
02668     newlist.Clear();
02669     newlist.Compress();
02670 
02671     KalmanPlaneSR *kp = 0;
02672     for(Int_t i=0; i<=nlast; ++i){
02673       kp = dynamic_cast<KalmanPlaneSR*>(KalmanPlaneList->At(i));
02674       if(kpswimfail[i]) delete kp;
02675       else{
02676         newlist.Add(kp);
02677         MSG("FitTrackSR",Msg::kDebug) << "adding thiskp  plane " 
02678                                       << kp->GetTrackCluster()->GetPlane() 
02679                                       << " min/max strip " 
02680                                       << kp->GetTrackCluster()->GetMinStrip() 
02681                                       << "/" 
02682                                       << kp->GetTrackCluster()->GetMaxStrip() 
02683                                       << " q/p " 
02684                                       << kp->GetFilStateValue(4,1) << endl;
02685 
02686       }//end else kp didnt fail swim
02687     }//end loop over kalman planes
02688 
02689     KalmanPlaneList->Clear();
02690     KalmanPlaneList->Compress();
02691 
02692     nlast = newlist.GetLast();
02693     MSG("FitTrackSR",Msg::kDebug) << "new KP list" << endl;
02694     
02695     for(Int_t i=0; i<=nlast; ++i){
02696       kp = dynamic_cast<KalmanPlaneSR*>(newlist.At(i));
02697       KalmanPlaneList->Add(kp);
02698       MSG("FitTrackSR",Msg::kDebug) << "  plane " 
02699                                     << kp->GetTrackCluster()->GetPlane() 
02700                                     << " min/max strip " 
02701                                     << kp->GetTrackCluster()->GetMinStrip()
02702                                     << "/" << kp->GetTrackCluster()->GetMaxStrip() 
02703                                     << " q/p " 
02704                                     << kp->GetFilStateValue(4,1) << endl;
02705 
02706     }//end loop over kalman planes to add to the list for the track
02707   }//end if didswimfail
02708 
02709   delete [] kpswimfail;
02710   
02711   if(lastkp) cfh.SetVtx(lastkp);
02712   
02713   return 0;
02714 }

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

Implements AlgBase.

Definition at line 83 of file AlgFitTrackSR.cxx.

References AlgReco::Calibrate(), 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, CandStripHandle::GetCharge(), CandContext::GetDataIn(), VldContext::GetDetector(), Registry::GetDouble(), KalmanPlaneSR::GetFilStateValue(), Registry::GetInt(), GetMomFromRange(), PlexPlaneId::GetPlane(), TrackClusterSR::GetPlane(), CandStripHandle::GetPlane(), UgliGeomHandle::GetPlaneIdFromZ(), InitKalmanFitParameters(), TrackClusterSR::IsValid(), Msg::kDebug, KeyOnClusterPlane(), Msg::kFatal, CalDigitType::kPE, PlaneView::kU, PlaneView::kV, MakeDaughterStripList(), MakeSliceClusterList(), MakeTrackClusterList(), MarkTrackClusters(), Munits::mm, MSG, pFor, pRev, pRev2, AlgTrack::SetdS(), SetPlaneParameters(), AlgTrack::SetT(), SetTrackEndParameters(), SetTrackParameters(), AlgTrack::SetUVZ(), CandStripHandle::StripSRKeyFromBegTime(), CandStripHandle::StripSRKeyFromPSEId(), SwimVertexAndEndPoints(), and tc.

00085 {
00086   // check inputs
00087 
00088 
00089   pRev=0;
00090   pRev2=0;
00091   pFor=0;
00092   assert(ch.InheritsFrom("CandFitTrackSRHandle"));
00093   CandFitTrackSRHandle &cfh = dynamic_cast<CandFitTrackSRHandle &>(ch);
00094   cfh.SetPass(0);
00095 
00096   assert(cx.GetDataIn());
00097   assert(cx.GetDataIn()->InheritsFrom("TObjArray"));
00098 
00099   CandTrackHandle *track0 = 0;
00100   const TObjArray *tclist = 0;
00101   const TObjArray *cxin =                        
00102         dynamic_cast<const TObjArray *>(cx.GetDataIn());
00103   for (Int_t i=0; i<=cxin->GetLast(); i++) {
00104     TObject *tobj = cxin->At(i);
00105     if (tobj->InheritsFrom("CandTrackHandle")) {
00106       track0 = dynamic_cast<CandTrackHandle*>(tobj);
00107       MSG("EventNum", Msg::kDebug) << "event number " 
00108                         << track0->GetCandRecord()->GetCandHeader()->GetSnarl() 
00109                         << endl;
00110     }
00111     if (tobj->InheritsFrom("TObjArray")) {
00112       tclist = dynamic_cast<TObjArray*>(tobj);
00113     }
00114   }
00115   assert(track0);
00116 
00117   // set finder track
00118   cfh.SetFinderTrack(track0);
00119 
00120   // get detector type
00121   const VldContext* vldcptr = track0->GetVldContext();
00122   assert(vldcptr);
00123   fVldContext = *vldcptr;
00124   fDetector = fVldContext.GetDetector();
00125   UgliGeomHandle ugh(fVldContext);
00126 
00127   // input algorithm parameters
00128 
00129   fParmMaxIterate = ac.GetInt("MaxIterate");
00130   fParmMaxIterate2 = ac.GetInt("MaxIterate2");
00131   fParmQPDiff = ac.GetDouble("QPDiff");
00132   fParmMinPlanePE = ac.GetDouble("MinPlanePE");
00133   fParmLastPlane = ac.GetInt("LastPlane");
00134   fParmIsCosmic = ac.GetInt("IsCosmic");
00135   fParmMaxLocalChi2 = ac.GetDouble("MaxLocalChi2");
00136   fParmMaxLocalPreChi2 = ac.GetDouble("MaxLocalPreChi2");
00137   fParmMaxAngleCovariance = ac.GetDouble("MaxAngleCovariance");
00138   fParmNSkipView = ac.GetInt("NSkipView");
00139   fParmNSkipActive = ac.GetInt("NSkipActive");
00140   fParmPassReducedChi2 = ac.GetDouble("PassReducedChi2");
00141   fParmPassPlaneAsymmetry = ac.GetDouble("PassPlaneAsymmetry");
00142   fParmPassMinPlaneAsymmetry = ac.GetInt("PassMinPlaneAsymmetry");
00143   fParmMaxImpactParameter = ac.GetDouble("MaxImpactParameter");
00144   fParmMinClusterCharge = ac.GetDouble("MinClusterCharge");
00145   fParmMaxClusterNStrip = ac.GetInt("MaxClusterNStrip");
00146   fParmDeltaChi2 = ac.GetDouble("DeltaChi2");
00147   fParmDeltaCovariance = ac.GetDouble("DeltaCovariance");
00148   fParmMisalignmentError = ac.GetDouble("MisalignmentError")*Munits::mm;
00149   fParmTPosError2 = ac.GetDouble("TPosError2");
00150   fParmInitialPositionError2 = ac.GetDouble("InitialPositionError2");
00151   fParmInitialSlopeError2 = ac.GetDouble("InitialSlopeError2");
00152   fParmInitialQPError2 = ac.GetDouble("InitialQPError2");
00153   fParmMaxLocalChi2 = ac.GetDouble("MaxLocalChi2");
00154   fParmMaxLocalPreChi2 = ac.GetDouble("MaxLocalPreChi2");
00155   fParmCovarianceScale = ac.GetDouble("CovarianceScale");
00156 
00157   //initialize the parameters from registry for the Kalman Plane stuff
00158   InitKalmanFitParameters(ac);
00159 
00160   // verify that looping won't ask for planes that don't exist
00161   const float way_downstream = 999999.;
00162   PlexPlaneId lastpln = ugh.GetPlaneIdFromZ(way_downstream);
00163   if (lastpln.GetPlane() < fParmLastPlane ) {
00164     fParmLastPlane = lastpln.GetPlane();
00165   }
00166 
00167   Int_t idir = -1;
00168   if(track0->GetTimeSlope()>0. || !fParmIsCosmic) idir = 1;
00169 
00170   //get an iterator over the strips in this slice
00171   CandStripHandle *strip = 0;
00172   CandStripHandleItr stripItr(track0->GetCandSlice()->GetDaughterIterator());
00173   
00174   //sort by plane and strip then by time using the navigation 2 sort functionality
00175   stripItr.SetSort2Fun(CandStripHandle::StripSRKeyFromPSEId, CandStripHandle::StripSRKeyFromBegTime);
00176     
00177   MSG("AlgFitTrackSR", Msg::kDebug) << "sorted the strips by time" << endl;
00178   
00179   stripItr.SetDirection(idir);
00180   stripItr.Reset();
00181 
00182   //create a TObjArray to hold the cluster objects you are going to make
00183   TObjArray *trackClusterList = new TObjArray(1,0);
00184 
00185   //fill the trackClusterList with the objects in tclist, if it is there
00186   if(tclist) MakeSliceClusterList(trackClusterList, tclist, stripItr, cfh);
00187   else MakeTrackClusterList(trackClusterList, stripItr, cfh, idir);
00188 
00189   MSG("AlgFitTrackSR", Msg::kDebug) << "made cluster list" << endl;
00190   
00191 // clone the trackClusterList so that you can have a master list of 
00192 // track clusters and one that keeps track of which clusters are used
00193 // in the fit track
00194 
00195   TObjArray *clusterList(trackClusterList);
00196   
00197   //make an iterator over clusterList and sort it by plane number
00198   TrackClusterSRItr clusterItr(clusterList);
00199   TrackClusterSRKeyFunc *clusterKf = clusterItr.CreateKeyFunc();
00200   clusterKf->SetFun(KeyOnClusterPlane);
00201   clusterItr.GetSet()->AdoptSortKeyFunc(clusterKf);
00202   clusterKf = 0;
00203   clusterItr.SetDirection(idir);
00204 
00205   clusterItr.ResetFirst();
00206   TrackClusterSR *cluster = 0;
00207 
00208   TObjArray planeClusterList(1000);
00209   Int_t iplane=0;
00210   Bool_t first=true;
00211   TObjArray * planeClusters;
00212   while( (cluster = clusterItr())) {
00213     if(cluster->GetPlane()!=iplane || first){
00214       iplane=cluster->GetPlane();
00215       MSG("FitTrackSR", Msg::kDebug) << "new cluster on plane " << iplane << endl;
00216       planeClusters = new TObjArray(1,0);
00217       planeClusterList.AddAt(planeClusters,iplane);
00218       planeClusters->Add(cluster);
00219       first=false;
00220     }
00221     else{
00222       MSG("FitTrackSR", Msg::kDebug) << "adding cluster on plane " << iplane << endl;
00223       planeClusters=(TObjArray *)planeClusterList.At(iplane);
00224       if(!planeClusters) MSG("AlgFitTrackSR", Msg::kFatal) << "planeCluster array not defined for plane " << iplane << endl;
00225       planeClusters->Add(cluster);
00226     }
00227   }
00228   
00229   Bool_t dir = kIterForward;
00230   if(idir!=1) dir=kIterBackward;
00231   TIter planeItr(&planeClusterList,idir);
00232   TrackClusterSR * tc=0;
00233   TObjArray * plnarray;
00234   while( (plnarray = (TObjArray * )planeItr.Next()) ){
00235     TIter clusterItr(plnarray);
00236     while( (tc = (TrackClusterSR *)clusterItr.Next()) ){
00237       
00238       MSG("FitTrackSR", Msg::kDebug) << "current cluster on plane " 
00239                                      << tc->GetPlane() << " is valid "
00240                                      << (Int_t)tc->IsValid() << endl;
00241     }
00242   }
00243   MSG("AlgFitTrackSR", Msg::kDebug) << "made cluster iterator" << endl;
00244 
00245   //get another iterator over the strips in the track this time
00246   CandStripHandleItr trackStripItr(track0->GetDaughterIterator());
00247   
00248   //sort by plane and strip then by time using the navigation 2 sort functionality
00249   trackStripItr.SetSort2Fun(CandStripHandle::StripSRKeyFromPSEId, CandStripHandle::StripSRKeyFromBegTime);
00250   
00251   /* for now leave in the old code to sort in 2 dimensions just in case
00252   CandStripHandleKeyFunc *trackStripKf = trackStripItr.CreateKeyFunc();
00253   //using a 2d sort - first sort on the PlexStripEndId
00254   trackStripKf->SetFun(CandStripHandle::StripSRKeyFromPSEId);
00255   trackStripItr.GetSet()->AdoptSortKeyFunc(trackStripKf,kTRUE,kFALSE);
00256   //now sort by time
00257   trackStripKf = trackStripItr.CreateKeyFunc();
00258   trackStripKf->SetFun(CandStripHandle::StripSRKeyFromBegTime);
00259   trackStripItr.GetSet()->AdoptSortKeyFunc(trackStripKf,kFALSE);
00260   trackStripKf = 0;
00261   */
00262     
00263   MSG("AlgFitTrackSR", Msg::kDebug) << "got iterator sorted over track strips" << endl;
00264 
00265   trackStripItr.SetDirection(idir);
00266   trackStripItr.Reset();
00267 
00268   Int_t begPlane = -1;
00269   Int_t endPlane = -1;
00270 
00271   //now mark the clusters to use in the track by setting the IsValid flag to true
00272   bool goodTrack = MarkTrackClusters(track0, planeClusterList, trackStripItr, cfh, 
00273                                      begPlane, endPlane);
00274 
00275   //set the track parameters based on the track passed into the algorithm
00276   SetTrackParameters(track0, cfh, idir);
00277 
00278   //set the endpoints and their direction cosines if track0 is a CandTrackSR
00279   SetTrackEndParameters(begPlane, endPlane, cfh, track0);
00280 
00281   MSG("AlgFitTrackSR", Msg::kDebug) << "start the fitting" << endl;
00282   
00283   //this track has fewer than 3 hits at least one view, so just add the strips
00284   //in the daughter list of track0 to cfh and quit
00285   if(!goodTrack){
00286           
00287           MSG("AlgFitTrackSR", Msg::kDebug) << "not a good track, bail" << endl;  
00288         // add daughter strips from the current track and bail
00289     CandStripHandleItr trk0stripItr(track0->GetDaughterIterator());
00290     while( (strip = trk0stripItr()) ){
00291       cfh.AddDaughterLink(*strip);
00292     }
00293         
00294         //need to set the end points
00295         
00296     return;
00297   }
00298   else{ 
00299     Int_t istatus = 0;
00300     //now to actually do the track fitting.  
00301  
00302     Int_t dosearch=1;
00303     Int_t niterate = DoKalmanFit(planeClusterList, cfh, istatus, idir, dosearch);
00304     if(istatus>=0){
00305       //get the first and last planes for the track
00306       KalmanPlaneSR *kp0 = cfh.GetKalmanPlane(0);
00307       KalmanPlaneSR *kp1 = cfh.GetKalmanPlane(cfh.GetKalmanLast());
00308       cfh.SetVtx(kp0);
00309       cfh.SetEnd(kp1);
00310 
00311       //set up variables for first and last planes in each view for track
00312       const KalmanPlaneSR *kpu0 = 0;
00313       const KalmanPlaneSR *kpv0 = 0;
00314       const KalmanPlaneSR *kpu1 = 0;
00315       const KalmanPlaneSR *kpv1 = 0;
00316       
00317       // clear STL maps--we only want fit planes as keys
00318       cfh.ClearMaps();
00319       SetPlaneParameters(cfh, kpu0, kpu1, kpv0, kpv1);
00320     }
00321     // check to see whether full fit converged.  If not, fall back to using
00322     // only track points.
00323     if(istatus<0 
00324            || cfh.GetNDOF()==0 
00325            || (cfh.GetNDOF()!=0 && cfh.GetChi2()/cfh.GetNDOF()>fParmPassReducedChi2) 
00326            ) {
00327       dosearch = 0;
00328       idir = -1;
00329       if(track0->GetTimeSlope()>0. || ! fParmIsCosmic) idir = 1;
00330 
00331       //fill the trackClusterList with the objects in tclist, if it is there
00332       cfh.ClearTrackClusterList();
00333       trackClusterList->Clear();
00334       trackClusterList->Compress();
00335       if(tclist) MakeSliceClusterList(trackClusterList, tclist, stripItr,
00336                                       cfh);
00337       else MakeTrackClusterList(trackClusterList, stripItr, cfh, idir);
00338 
00339       //now mark the clusters to use in the track by setting the IsValid flag to true
00340       begPlane = -1;
00341       endPlane = -1;
00342       goodTrack = MarkTrackClusters(track0, planeClusterList, trackStripItr, cfh, 
00343                                          begPlane, endPlane);
00344    
00345       //set the track parameters based on the track passed into the algorithm
00346       SetTrackParameters(track0, cfh, idir);
00347       
00348       //set the endpoints and their direction cosines if track0 is a CandTrackSR
00349       SetTrackEndParameters(begPlane, endPlane, cfh, track0);
00350       niterate = DoKalmanFit(planeClusterList, cfh, istatus, idir, dosearch);      
00351     }
00352 
00353     //set the number of iterations for the fit
00354     cfh.SetNIterate(niterate);
00355 
00356     if(istatus>=0){
00357       //get the first and last planes for the track
00358       KalmanPlaneSR *kp0 = cfh.GetKalmanPlane(0);
00359       KalmanPlaneSR *kp1 = cfh.GetKalmanPlane(cfh.GetKalmanLast());
00360       cfh.SetVtx(kp0);
00361       cfh.SetEnd(kp1);
00362 
00363       //set up variables for first and last planes in each view for track
00364       const KalmanPlaneSR *kpu0 = 0;
00365       const KalmanPlaneSR *kpv0 = 0;
00366       const KalmanPlaneSR *kpu1 = 0;
00367       const KalmanPlaneSR *kpv1 = 0;
00368       
00369       // clear STL maps--we only want fit planes as keys
00370       cfh.ClearMaps();
00371       SetPlaneParameters(cfh, kpu0, kpu1, kpv0, kpv1);
00372 
00373       //gotta add those strips to the daughter list for this track
00374       MakeDaughterStripList(cfh);
00375             
00376       // create array of strips not in fit track
00377       Double_t planepe[1000];
00378       for(int i=0; i<1000; ++i){
00379         planepe[i] = 0.;
00380       }
00381 
00382       stripItr.ResetFirst();
00383       CandStripHandle *strip = 0;
00384       while( (strip = stripItr()) ){
00385         if(strip->GetPlane()>=0) 
00386           planepe[strip->GetPlane()] += strip->GetCharge(CalDigitType::kPE);
00387       }//end loop over strips to get charge for each plane
00388 
00389       // Find planes closest to vertex and end with adjacent planes meeting
00390       // PE cut
00391       Int_t plane0 = -999;
00392       Int_t plane1 = -999;
00393       for(int i=0; i<1000; ++i){
00394         if(planepe[i]>=fParmMinPlanePE && planepe[i+1]>=fParmMinPlanePE){
00395           if(idir>0){
00396             if(plane0<0) plane0 = i;
00397             if(plane1<0 || i+1>plane1) plane1 = i+1;
00398           } 
00399           else{
00400             if(plane1<0) plane1 = i;
00401             if(plane0<0 || i+1>plane0) plane0 = i+1;
00402           }
00403         }//end if enough signal on these two planes
00404       }//end loop over planes to find the first and last ones
00405 
00406     
00407       
00408       Double_t vtxqp = kp0->GetFilStateValue(4,1);
00409 
00410       // tweak qp to force agreement with truth !!!!
00411       vtxqp *= 1.01+0.1*TMath::Abs(vtxqp);
00412       
00413       //if we are looking at cosmic rays, need to swim the track out of the 
00414       //detector
00415       if(fParmIsCosmic) SwimVertexAndEndPoints(cfh, kp0, kp1, 
00416                                                kpu0, kpu1, 
00417                                                kpv0, kpv1, 
00418                                                planepe, plane0, plane1, 
00419                                                vtxqp, idir);
00420 
00421       if(vtxqp==0.){
00422         // infinite momentum, choose finite number to prevent floating point exceptions
00423         vtxqp = 0.000001; // 1 PeV
00424         cfh.SetEMCharge(0.);
00425       } 
00426       
00427       //  set momentum and charge sign
00428       cfh.SetMomentumCurve(1./TMath::Abs(vtxqp));
00429       if(vtxqp>0.) cfh.SetEMCharge(1.);
00430       else if(vtxqp<0.) cfh.SetEMCharge(-1.);
00431       
00432       //  fill time and range maps
00433       SetT(&cfh);
00434       SetdS(&cfh);
00435     
00436       // set the momentum from range; taken from PDG R/M vs p/M plot for Fe  
00437       Double_t range = cfh.GetRange();
00438       // if(range>0.) cfh.SetMomentumRange(0.105658*exp(1.0363*log(range)-4.383));
00439       // until June 2006 used to be   if(range>0.) cfh.SetMomentum(0.048 + 1.660e-3*range + 3.057e-8*range*range);     
00440 
00441       if (range>0.){
00442         cfh.SetMomentum(GetMomFromRange(range));   
00443       }
00444       else cfh.SetMomentum(0.);      
00445       Calibrate(&cfh);
00446       
00447     } // end if istatus>=0
00448 
00449     cfh.ClearTrackClusterList();
00450     cfh.ClearKalmanPlaneList();
00451     
00452     //check and see if the track passes the cuts for a good track
00453     Float_t nplaneu = (Float_t)(cfh.GetNPlane(PlaneView::kU));
00454     Float_t nplanev = (Float_t)(cfh.GetNPlane(PlaneView::kV));
00455     Int_t nplane = cfh.GetNPlane();
00456     if(istatus>=0 && cfh.GetNDOF()>0 
00457        && cfh.GetChi2()/cfh.GetNDOF()<=fParmPassReducedChi2
00458        && fParmDeltaChi2<=0. && fParmDeltaCovariance<=0.
00459        && (nplane<fParmPassMinPlaneAsymmetry 
00460            || TMath::Abs((nplaneu-nplanev)/(1.*nplane))<=fParmPassPlaneAsymmetry )
00461        )
00462       cfh.SetPass(1);
00463     else {
00464       cfh.SetPass(0);
00465       SetTrackParameters(track0, cfh, idir);
00466       CandStripHandleItr trkstripItr(cfh.GetDaughterIterator());
00467       while( (strip = trkstripItr()) ){
00468         cfh.RemoveDaughter(strip);
00469       }
00470       CandStripHandleItr trk0stripItr(track0->GetDaughterIterator());
00471       while( (strip = trk0stripItr()) ){
00472         cfh.AddDaughterLink(*strip);
00473       }
00474       // fill time and range maps
00475       SetT(&cfh);
00476       SetUVZ(&cfh);
00477       SetdS(&cfh);
00478       Calibrate(&cfh);
00479     }
00480   }//end if enough planes to try fitting the track
00481 
00482   //delete the track clusters you found earlier
00483   trackClusterList->Delete();
00484   delete trackClusterList;
00485 
00486   for (Int_t iplane=0;iplane<1000;iplane++){
00487     planeClusters=(TObjArray *)planeClusterList.At(iplane);
00488     if(planeClusters){
00489       delete planeClusters;
00490     }
00491   }
00492   //  cout << " # for swim " << pFor << " # rev swim " << pRev << " # rev search swim " << pRev2 << endl;
00493   return;
00494 }

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(), Msg::kDebug, kdUdZ, kdVdZ, kQoverP, kU, PlaneView::kU, kV, PlaneView::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().

03975 {
03976 
03977   //loop over the kalman planes to get the first and last ones in each view,
03978   //as well as setting the u, v, du/dz, and dv/dz values for each plane
03979   const KalmanPlaneSR *kp = 0;
03980   TrackClusterSR *tc = 0;
03981   Int_t thisplane = 0;
03982   Double_t thisu[2];
03983   Double_t thisv[2];
03984   Double_t erru2[2];
03985   Double_t errv2[2];
03986   Double_t avgu = 0.;
03987   Double_t avgv = 0.;
03988   Double_t avgdu = 0.;
03989   Double_t avgdv = 0.;
03990   Double_t thisdu[2];
03991   Double_t thisdv[2];
03992   Double_t errdu2[2];
03993   Double_t errdv2[2];
03994   
03995   Double_t chi2=0.;
03996   Double_t dchi2 = 0.;
03997   Double_t dprechi2 = 0.;
03998   Double_t thisqp[2];
03999   Double_t thiseqp2[2];
04000   Double_t avgqp = 0.;
04001   
04002   for(Int_t index=0; index<=cfh.GetKalmanLast(); ++index){
04003     kp = cfh.GetKalmanPlane(index);
04004           MSG("AlgFitTrackSR", Msg::kDebug) << "plane positions " << kp->GetTrackCluster()->GetPlane() << " " 
04005                   << kp->GetTrackCluster()->GetZPos() 
04006                   << " " << kp->GetFilStateValue(kU,0) 
04007                   << " " << kp->GetFilStateValue(kV,0) << endl;
04008     if(kp->GetTrackCluster()->GetPlaneView()==PlaneView::kU) {
04009       if(!kpu0) kpu0 = kp;
04010       kpu1 = kp;
04011     }
04012     if(kp->GetTrackCluster()->GetPlaneView()==PlaneView::kV){
04013       if(!kpv0) kpv0 = kp;
04014       kpv1 = kp;
04015     }
04016     
04017     tc = kp->GetTrackCluster();
04018     thisplane = tc->GetPlane();
04019      //loop over the clusters in the track
04020     cfh.SetTrackPointError(thisplane,tc->GetTPosError());
04021 
04022     //get the forward and backward values for u and v and the errors
04023     //on those values.  the errors have minimum value of 1.e-12
04024     thisu[0] = kp->GetFilStateValue(kU,0);
04025     thisu[1] = kp->GetFilStateValue(kU,1);
04026     thisv[0] = kp->GetFilStateValue(kV,0);
04027     thisv[1] = kp->GetFilStateValue(kV,1);
04028     erru2[0] = TMath::Max(kp->GetFilCovarianceMatrixValue(0,kU,kU), 1.e-12);
04029     erru2[1] = TMath::Max(kp->GetFilCovarianceMatrixValue(1,kU,kU), 1.e-12);
04030     errv2[0] = TMath::Max(kp->GetFilCovarianceMatrixValue(0,kV,kV), 1.e-12);
04031     errv2[1] = TMath::Max(kp->GetFilCovarianceMatrixValue(1,kV,kV), 1.e-12);
04032 
04033     //average forward and backward fit values
04034     avgu = (thisu[0]/erru2[0]+thisu[1]/erru2[1])/(1./erru2[0]+1./erru2[1]);
04035     avgv = (thisv[0]/errv2[0]+thisv[1]/errv2[1])/(1./errv2[0]+1./errv2[1]);
04036     //set the u an v values for this plane
04037     cfh.SetU(thisplane,avgu);
04038     cfh.SetV(thisplane,avgv);
04039     
04040     thisdu[0] = kp->GetFilStateValue(kdUdZ,0);
04041     thisdu[1] = kp->GetFilStateValue(kdUdZ,1);
04042     thisdv[0] = kp->GetFilStateValue(kdVdZ,0);
04043     thisdv[1] = kp->GetFilStateValue(kdVdZ,1);
04044     errdu2[0] = TMath::Max(kp->GetFilCovarianceMatrixValue(0,kdUdZ,kdUdZ),1.e-12);
04045     errdu2[1] = TMath::Max(kp->GetFilCovarianceMatrixValue(1,kdUdZ,kdUdZ),1.e-12);
04046     errdv2[0] = TMath::Max(kp->GetFilCovarianceMatrixValue(0,kdVdZ,kdVdZ),1.e-12);
04047     errdv2[1] = TMath::Max(kp->GetFilCovarianceMatrixValue(1,kdVdZ,kdVdZ),1.e-12);
04048     
04049     // average forward and backward fit values
04050     avgdu = (thisdu[0]/errdu2[0]+thisdu[1]/errdu2[1])/(1./errdu2[0]+1./errdu2[1]);
04051     avgdv = (thisdv[0]/errdv2[0]+thisdv[1]/errdv2[1])/(1./errdv2[0]+1./errdv2[1]);
04052 
04053     //set the du/dz and dv/dz values for this plane
04054     cfh.SetdUdZ(thisplane,avgdu);
04055     cfh.SetdVdZ(thisplane,avgdv);
04056   
04057     //now do the q/p values
04058     dchi2 = TMath::Max(kp->GetFilChi2(0),kp->GetFilChi2(1));
04059     dprechi2 = TMath::Max(kp->GetPreChi2(0),kp->GetPreChi2(1));
04060     cfh.SetPlaneChi2(thisplane,dchi2);
04061     cfh.SetPlanePreChi2(thisplane,dprechi2);
04062     thisqp[0] = kp->GetFilStateValue(kQoverP,0);
04063     thisqp[1] = kp->GetFilStateValue(kQoverP,1);
04064     thiseqp2[0] = TMath::Max(kp->GetFilCovarianceMatrixValue(0,kQoverP,kQoverP),1.e-12);
04065     thiseqp2[1] = TMath::Max(kp->GetFilCovarianceMatrixValue(1,kQoverP,kQoverP),1.e-12);
04066 
04067     if(thisqp[0]==0){
04068       avgqp=thisqp[1];
04069     }
04070     if(thisqp[1]==0){
04071       avgqp=thisqp[0];
04072     }
04073     else
04074       avgqp = (thisqp[0]/thiseqp2[0]+thisqp[1]/thiseqp2[1])/(1./thiseqp2[0]+1./thiseqp2[1]);
04075 
04076     cfh.SetPlaneQP(thisplane,avgqp);
04077     chi2 += dchi2;
04078   
04079   }//end loop over kalman planes to set u,v,du/dz,dv/dz and find first and last
04080   //planes in each view
04081 
04082   //set the track chi^2 for the q/p fit
04083   cfh.SetChi2(chi2);  
04084   Int_t nkalmanlast = cfh.GetKalmanLast();
04085   
04086   //set the degrees of freedom and the errors for the vtx and end
04087   cfh.SetNDOF(nkalmanlast-4);
04088   cfh.SetVtxUError(TMath::Sqrt(cfh.GetKalmanPlane(0)->GetFilCovarianceMatrixValue(1, kU, kU)));
04089   cfh.SetVtxVError(TMath::Sqrt(cfh.GetKalmanPlane(0)->GetFilCovarianceMatrixValue(1, kV, kV)));
04090   cfh.SetVtxdUError(TMath::Sqrt(cfh.GetKalmanPlane(0)->GetFilCovarianceMatrixValue(1, kdUdZ, kdUdZ)));
04091   cfh.SetVtxdVError(TMath::Sqrt(cfh.GetKalmanPlane(0)->GetFilCovarianceMatrixValue(1, kdVdZ, kdVdZ)));
04092   cfh.SetVtxQPError(TMath::Sqrt(cfh.GetKalmanPlane(0)->GetFilCovarianceMatrixValue(1,kQoverP,kQoverP)));
04093   cfh.SetEndUError(TMath::Sqrt(cfh.GetKalmanPlane(nkalmanlast)->GetFilCovarianceMatrixValue(0,kU,kU)));
04094   cfh.SetEndVError(TMath::Sqrt(cfh.GetKalmanPlane(nkalmanlast)->GetFilCovarianceMatrixValue(0,kV,kV)));
04095   cfh.SetEnddUError(TMath::Sqrt(cfh.GetKalmanPlane(nkalmanlast)->GetFilCovarianceMatrixValue(0,kdUdZ,kdUdZ)));
04096   cfh.SetEnddVError(TMath::Sqrt(cfh.GetKalmanPlane(nkalmanlast)->GetFilCovarianceMatrixValue(0,kdVdZ,kdVdZ)));
04097   cfh.SetEndQPError(TMath::Sqrt(cfh.GetKalmanPlane(nkalmanlast)->GetFilCovarianceMatrixValue(0,kQoverP,kQoverP)));
04098   
04099   return;
04100 }

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().

02111 {
02112   cfh.SetVtxPlane(begPlane);
02113   cfh.SetVtxZ(track0->GetZ(begPlane));
02114 
02115   if(begPlane!=cfh.GetVtxPlane() && track0->IsTPosValid(begPlane)) {
02116     cfh.SetVtxU(track0->GetU(begPlane));
02117     cfh.SetVtxV(track0->GetV(begPlane));
02118    
02119     cfh.SetVtxT(track0->GetT(begPlane));
02120     if(track0->InheritsFrom("CandTrackSRHandle")){
02121       // CandTrackSR calculates angles at each plane, use those
02122       cfh.SetDirCosU(dynamic_cast<const CandTrackSRHandle *>(track0)->GetDirCosU(begPlane));
02123       cfh.SetDirCosV(dynamic_cast<const CandTrackSRHandle *>(track0)->GetDirCosV(begPlane));
02124       cfh.SetDirCosZ(dynamic_cast<const CandTrackSRHandle *>(track0)->GetDirCosZ(begPlane));
02125     }
02126   }
02127 
02128   cfh.SetEndPlane(endPlane);
02129   cfh.SetEndZ(track0->GetZ(endPlane));
02130   if(endPlane!=cfh.GetEndPlane() && track0->IsTPosValid(endPlane)){
02131     cfh.SetEndU(track0->GetU(endPlane));
02132     cfh.SetEndV(track0->GetV(endPlane));
02133     cfh.SetEndT(track0->GetT(endPlane));
02134     if(track0->InheritsFrom("CandTrackSRHandle")){
02135       // CandTrackSR calculates angles at each plane, use those
02136       cfh.SetEndDirCosU(dynamic_cast<const CandTrackSRHandle *>(track0)->GetDirCosU(endPlane));
02137       cfh.SetEndDirCosV(dynamic_cast<const CandTrackSRHandle *>(track0)->GetDirCosV(endPlane));
02138       cfh.SetEndDirCosZ(dynamic_cast<const CandTrackSRHandle *>(track0)->GetDirCosZ(endPlane));
02139     }
02140   }
02141   return;
02142 }

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

Definition at line 1847 of file AlgFitTrackSR.cxx.

References fParmIsCosmic, 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().

01850 {
01851 
01852   cfh.SetCandSlice(track0->GetCandSlice());
01853   cfh.SetDirCosU(track0->GetDirCosU());
01854   cfh.SetDirCosV(track0->GetDirCosV());
01855   cfh.SetDirCosZ(track0->GetDirCosZ());
01856   cfh.SetVtxU(track0->GetVtxU());
01857   cfh.SetVtxV(track0->GetVtxV());
01858   cfh.SetVtxZ(track0->GetVtxZ());
01859   cfh.SetVtxT(track0->GetVtxT());
01860   cfh.SetVtxPlane(track0->GetVtxPlane());
01861   cfh.SetEndDirCosU(track0->GetEndDirCosU());
01862   cfh.SetEndDirCosV(track0->GetEndDirCosV());
01863   cfh.SetEndDirCosZ(track0->GetEndDirCosZ());
01864   cfh.SetEndU(track0->GetEndU());
01865   cfh.SetEndV(track0->GetEndV());
01866   cfh.SetEndZ(track0->GetEndZ());
01867   cfh.SetEndT(track0->GetEndT());
01868   cfh.SetEndPlane(track0->GetEndPlane());
01869   cfh.SetTimeSlope(track0->GetTimeSlope());
01870   cfh.SetTimeOffset(track0->GetTimeOffset());
01871   cfh.SetMomentumRange(track0->GetMomentum());
01872   Double_t initQP=0.0;
01873   if(track0->GetMomentum()!=0 && !fParmIsCosmic) initQP=-1./track0->GetMomentum();
01874   cfh.SetInitialQP(initQP);
01875   if(!fParmIsCosmic) cfh.SetEndQP(fParmMaxQP);
01876   // set ds, range
01877   for(Int_t iplane = cfh.GetVtxPlane(); iplane*direction<=cfh.GetEndPlane()*direction; 
01878       iplane+=direction){
01879     if(track0->IsTPosValid(iplane)){
01880       // u and v coordinates have been calculated for this plane
01881  //loop over the clusters in the track
01882       cfh.SetTrackPointError(iplane,track0->GetTrackPointError(iplane));
01883       cfh.SetdS(iplane,track0->GetdS(iplane));
01884       cfh.SetRange(iplane,track0->GetRange(iplane));
01885     }
01886   }
01887 
01888   //set stuff if the track was of the SR variety
01889   const CandTrackSRHandle *tracksr = 0;
01890   if (track0->InheritsFrom("CandTrackSRHandle")) {
01891     tracksr = dynamic_cast<const CandTrackSRHandle*>(track0);
01892     cfh.SetNTrackStrip(tracksr->GetNTrackStrip());
01893     cfh.SetNTrackDigit(tracksr->GetNTrackDigit());
01894     cfh.SetNTimeFitDigit(tracksr->GetNTimeFitDigit());
01895     cfh.SetTimeFitChi2(tracksr->GetTimeFitChi2());
01896   }
01897   cfh.SetNIterate(-1);
01898 
01899   return;
01900 }

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

Definition at line 1590 of file AlgFitTrackSR.cxx.

References done(), fParmMaxAngle, fParmMaxQP, SwimParticle::GetCharge(), SwimParticle::GetDirection(), SwimParticle::GetMomentumModulus(), SwimParticle::GetPosition(), Msg::kDebug, kdUdZ, kdVdZ, kInvSqrt2, kQoverP, kU, kV, MSG, s(), SwimParticle::SetCharge(), SwimSwimmer::SwimBackward(), and SwimSwimmer::SwimForward().

Referenced by CalculatePreState(), CalculatePropagator(), FindNumSkippedPlanes(), and SwimVertexAndEndPoints().

01593 {
01594 
01595   // return false if angles get too large or momentum is too small
01596 
01597   //-------------------------------------------------------------------
01598   // Using Swimmer package
01599   //-------------------------------------------------------------------
01600   SwimSwimmer s(*vldc);
01601   
01602   // Add an action to the stepper - in this case a debug printout
01603   //  SwimPrintStepAction printStep;
01604   //  s.GetStepper()->AddStepAction(&printStep);
01605 
01606   // Setup the initial condition for the stepper  
01607   //the swimmer works in xyz coordinates, not uvz, so need to rotate 
01608   //the coordinates before putting them into a position vector
01609   Double_t x = kInvSqrt2*(swimstate[kU]-swimstate[kV]);
01610   Double_t y = kInvSqrt2*(swimstate[kU]+swimstate[kV]);
01611         
01612   TVector3 position(x,y,zPos);
01613   double charge = 0.0;
01614 
01615   if(swimstate[kQoverP]>0.0) charge = 1.0;
01616   else if(swimstate[kQoverP]<0.0) charge = -1.0;
01617   else{
01618 
01619     // infinite momentum, straight line extrapolation
01620     Double_t delz = finalZ - zPos;
01621     state[kU] = swimstate[kU] + swimstate[kdUdZ]*delz;
01622     state[kV] = swimstate[kV] + swimstate[kdVdZ]*delz;
01623     state[kdUdZ] = swimstate[kdUdZ];
01624     state[kdVdZ] = swimstate[kdVdZ];
01625     state[kQoverP] = swimstate[kQoverP];
01626     return 1;
01627   }
01628 
01629   double pz;
01630     // forward in time: idir=0; backward in time: idir = 1
01631     // Momentum is always forward in time
01632   if((finalZ>zPos && idir==0) || (finalZ<zPos && idir==1))
01633     pz = 1.0/(TMath::Sqrt(1+swimstate[kdUdZ]*swimstate[kdUdZ]
01634                           +swimstate[kdVdZ]*swimstate[kdVdZ])
01635               *TMath::Abs(swimstate[kQoverP])); 
01636   
01637   else if((finalZ<zPos && idir==0) || (finalZ>zPos && idir==1))
01638     pz = -1.0/(TMath::Sqrt(1+swimstate[kdUdZ]*swimstate[kdUdZ]
01639                            +swimstate[kdVdZ]*swimstate[kdVdZ])
01640                *TMath::Abs(swimstate[kQoverP]));
01641   else{
01642     state[kU] = swimstate[kU];
01643     state[kV] = swimstate[kV];
01644     state[kdUdZ] = swimstate[kdUdZ];
01645     state[kdVdZ] = swimstate[kdVdZ];
01646     state[kQoverP] = swimstate[kQoverP];
01647     return 1;
01648   }
01649 
01650   //the swimmer works in xyz coordinates, not uvz, so need to rotate 
01651   //the coordinates before putting them into a position vector
01652   Double_t dxdz = kInvSqrt2*(swimstate[kdUdZ]-swimstate[kdVdZ]);
01653   Double_t dydz = kInvSqrt2*(swimstate[kdUdZ]+swimstate[kdVdZ]);
01654 
01655   TVector3 momentum(pz*dxdz,pz*dydz,pz);
01656   SwimParticle muon(position,momentum);
01657   muon.SetCharge(charge);
01658   SwimZCondition zc(finalZ);
01659   bool done;
01660 
01661   if(idir==0) done = s.SwimForward(muon,zc);
01662   else if(idir==1) done = s.SwimBackward(muon,zc);
01663   else abort();
01664 
01665   if(muon.GetDirection().Z()==0.0){
01666     state[kU] = swimstate[kU];
01667     state[kV] = swimstate[kV];
01668     state[kdUdZ] = swimstate[kdUdZ];
01669     state[kdVdZ] = swimstate[kdVdZ];
01670     state[kQoverP] = swimstate[kQoverP];
01671     return 0;
01672   }
01673 
01674   if(muon.GetMomentumModulus()==0.0){
01675     state[kU] = swimstate[kU];
01676     state[kV] = swimstate[kV];
01677     state[kdUdZ] = swimstate[kdUdZ];
01678     state[kdVdZ] = swimstate[kdVdZ];
01679     state[kQoverP] = swimstate[kQoverP];
01680     return 0;
01681   }
01682 
01683   state[kU] = kInvSqrt2*(muon.GetPosition().X()+muon.GetPosition().Y());
01684   state[kV] = kInvSqrt2*(muon.GetPosition().Y()-muon.GetPosition().X());
01685   
01686   dxdz = muon.GetDirection().X()/muon.GetDirection().Z();
01687   dydz = muon.GetDirection().Y()/muon.GetDirection().Z();
01688   state[kdUdZ] = kInvSqrt2*(dxdz+dydz);
01689   state[kdVdZ] = kInvSqrt2*(dydz-dxdz);
01690   state[kQoverP] = muon.GetCharge()/muon.GetMomentumModulus();
01691   //  MSG("FitTrackSR", Msg::kDebug) << "state = " << state[kU] << " " 
01692   //                              << state[kV] << " " 
01693   //                              << state[kdUdZ] << " " 
01694   //                               << state[kdVdZ] << " " 
01695   //                               << state[kQoverP] << endl;
01696   if(TMath::Abs(state[kdUdZ])>2.*fParmMaxAngle 
01697      || TMath::Abs(state[kdVdZ])>2.*fParmMaxAngle 
01698      || TMath::Abs(state[kQoverP])>fParmMaxQP){
01699  MSG("FitTrackSR",Msg::kDebug) << "swimming failed, u,v,du/dz,dv,dz/qp = " 
01700                                   << state[kU] << " " 
01701                                   << state[kV] << " " 
01702                                   << state[kdUdZ] << " " 
01703                                   << state[kdVdZ] << " " 
01704                                   << state[kQoverP] << endl; 
01705    return 0;
01706   }
01707 
01708   return 1;
01709 }

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(), PlexPlaneId::IsValid(), UgliScintPlnHandle::IsValid(), kdUdZ, kdVdZ, kInvSqrt2, kQoverP, kU, PlaneView::kU, kV, PlaneView::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().

03680 {
03681 
03682   Double_t swimstate0[5]; 
03683   Double_t tmpstate[5];  Double_t swimstatez0 = kp0->GetTrackCluster()->GetZPos();
03684   Double_t swimstate1[5];
03685 
03686   UgliGeomHandle ugh(fVldContext);
03687 
03688   for(Int_t i = 0; i<5; ++i){
03689     tmpstate[i] = 0.;
03690     swimstate0[i] = kp0->GetFilStateValue(i,1);
03691     swimstate1[i] = kp1->GetFilStateValue(i,0);
03692   }
03693   
03694   Double_t swimstatez1 = kp1->GetTrackCluster()->GetZPos();
03695   Int_t thisplane = kp0->GetTrackCluster()->GetPlane()-direction;
03696   Double_t thisx = kInvSqrt2*(swimstate0[kU]-swimstate0[kV]);
03697   Double_t thisy = kInvSqrt2*(swimstate0[kU]+swimstate0[kV]);
03698   Double_t thisu = swimstate0[kU];
03699   Double_t thisv = swimstate0[kV];
03700   Double_t thisz = swimstatez0;
03701   Double_t linx = thisx;
03702   Double_t liny = thisy;
03703   Double_t linu = thisu;
03704   Double_t linv = thisv;
03705   Double_t linz = thisz;
03706   Double_t lindx = kInvSqrt2*(swimstate0[kdUdZ]-swimstate0[kdVdZ]);
03707   Double_t lindy = kInvSqrt2*(swimstate0[kdUdZ]+swimstate0[kdVdZ]);
03708   Double_t lindu = swimstate0[kdUdZ];
03709   Double_t lindv = swimstate0[kdVdZ];
03710   Double_t linds = 0.;
03711   Double_t lindz = 0.;
03712   Double_t oldu = 0.;
03713   Double_t oldv = 0.;
03714   Double_t oldz = 0.;
03715   Double_t thisds = 0.;
03716   Double_t thisdz = 0.;
03717   Int_t nplaneskip = 0;
03718   Int_t swimflag = 1;
03719   Bool_t befvtx = 1;
03720   Double_t du = 0.;
03721   Double_t dv = 0.;
03722   Double_t dz = 0.;
03723   Double_t dsdz = 0.;
03724   Double_t delu = 0.;
03725   Double_t delv = 0.;
03726   Double_t delz = 0.;
03727 
03728   while(thisplane>=0 && thisplane<=fParmLastPlane 
03729         && ((swimflag && TMath::Abs(thisx)<4. && TMath::Abs(thisy)<4. && TMath::Abs(thisu)<4. 
03730              && TMath::Abs(thisv)<4.) 
03731             || (TMath::Abs(linx)<4. && TMath::Abs(liny)<4. && TMath::Abs(linu)<4. 
03732                 && TMath::Abs(linv)<4.)
03733             )
03734         ){
03735     
03736     PlexPlaneId plexplane(fDetector,thisplane,false);
03737     if(plexplane.IsValid()){
03738 
03739       UgliScintPlnHandle planeid = ugh.GetScintPlnHandle(plexplane);
03740       
03741       if(planeid.IsValid()){
03742         
03743         oldu = swimstate0[kU];
03744         oldv = swimstate0[kV];
03745         oldz = swimstatez0;
03746         dz = (Double_t)planeid.GetZ0()-linz;
03747 
03748         if(TMath::Abs(linx)<4. && TMath::Abs(liny)<4. 
03749            && TMath::Abs(linu)<4. && TMath::Abs(linv)<4.){
03750 
03751           linz = (Double_t)planeid.GetZ0();
03752           linx += dz*lindx;
03753           liny += dz*lindy;
03754           linu += dz*lindu;
03755           linv += dz*lindv;
03756           if(!befvtx){
03757             du = dz*lindu;
03758             dv = dz*lindv;
03759             lindz += TMath::Abs(dz);
03760             linds += TMath::Sqrt(du*du+dv*dv+dz*dz);
03761           }
03762         }//end if inside the detector
03763         
03764         if(swimflag && TMath::Abs(thisx)<4. && TMath::Abs(thisy)<4. 
03765            && TMath::Abs(thisu)<4. && TMath::Abs(thisv)<4.) {
03766           
03767           swimflag = Swim(swimstate0,tmpstate,(Double_t)swimstatez0,
03768                               (Double_t)planeid.GetZ0(),1, &fVldContext);
03769           swimstatez0 = (Double_t)planeid.GetZ0();
03770           swimstate0[kU] = tmpstate[kU];
03771           swimstate0[kV] = tmpstate[kV];
03772           swimstate0[kdUdZ] = tmpstate[kdUdZ];
03773           swimstate0[kdVdZ] = tmpstate[kdVdZ];
03774           swimstate0[kQoverP] = tmpstate[kQoverP];
03775           
03776           if (TMath::Abs(tmpstate[kU])>4. || TMath::Abs(tmpstate[kV])>4.) swimflag = 0;
03777 
03778           thisx = kInvSqrt2*(swimstate0[kU]-swimstate0[kV]);
03779           thisy = kInvSqrt2*(swimstate0[kU]+swimstate0[kV]);
03780           thisu = swimstate0[kU];
03781           thisv = swimstate0[kV];
03782           thisz = swimstatez0;
03783           
03784           if(!swimflag && nplaneskip<3 && thisplane*direction>=plane0*direction){
03785             if(!thisu) cfh.SetVtxU(thisu);
03786             if(!thisv) cfh.SetVtxV(thisv);
03787             cfh.SetVtxZ(thisz);
03788             cfh.SetVtxPlane(thisplane);
03789             dsdz = TMath::Sqrt(1.+swimstate0[kdUdZ]*swimstate0[kdUdZ]
03790                                         +swimstate0[kdVdZ]*swimstate0[kdVdZ]);
03791             cfh.SetDirCosU((Double_t)(direction)*swimstate0[kdUdZ]/dsdz);
03792             cfh.SetDirCosV((Double_t)(direction)*swimstate0[kdVdZ]/dsdz);
03793             cfh.SetDirCosZ((Double_t)(direction)/dsdz);
03794             cfh.SetU(thisplane,thisu);
03795             cfh.SetV(thisplane,thisv);
03796 
03797             vtxqp = swimstate0[kQoverP];
03798           }//end if swimflag and not too many planes skipped and this plane is beyond the
03799           //vertex
03800           else{
03801             befvtx = 0;
03802             delu = thisu-oldu;
03803             delv = thisv-oldv;
03804             delz = thisz-oldz;
03805             thisds += TMath::Sqrt(delu*delu+delv*delv+delz*delz);
03806             thisdz += TMath::Abs(delz);
03807           }
03808 
03809           //look to see if this plane was hit in the slice, if
03810           //not increment the number of skipped planes
03811           if(planepe[thisplane] > 0.) nplaneskip = 0;
03812           else  ++nplaneskip;
03813           
03814         }//end if in the detector
03815       }//end if planeid is valid
03816     }//end if plexplane is valid
03817    
03818     thisplane -= direction;
03819   }//end loop over planes
03820  
03821   cfh.SetVtxTrace(min(thisds,linds));
03822   cfh.SetVtxTraceZ(min(thisdz,lindz));
03823   
03824   if(kpu0) cfh.SetVtxExtrapolate(TMath::Abs(kpu0->GetTrackCluster()->GetPlane()-cfh.GetVtxPlane()),
03825                                  PlaneView::kU);
03826   else cfh.SetVtxExtrapolate(-1,PlaneView::kU);
03827   
03828   if(kpv0) cfh.SetVtxExtrapolate(TMath::Abs(kpv0->GetTrackCluster()->GetPlane()-cfh.GetVtxPlane()),
03829                                  PlaneView::kV);
03830   else cfh.SetVtxExtrapolate(-1,PlaneView::kV);
03831   
03832   thisplane = kp1->GetTrackCluster()->GetPlane()+direction;
03833   thisx = kInvSqrt2*(swimstate1[kU]-swimstate1[kV]);
03834   thisy = kInvSqrt2*(swimstate1[kU]+swimstate1[kV]);
03835   thisu = swimstate1[kU];
03836   thisv = swimstate1[kV];
03837   thisz = swimstatez1;
03838   linx = thisx;
03839   liny = thisy;
03840   linu = thisu;
03841   linv = thisv;
03842   linz = thisz;
03843   lindx = kInvSqrt2*(swimstate1[kdUdZ]-swimstate1[kdVdZ]);
03844   lindy = kInvSqrt2*(swimstate1[kdUdZ]+swimstate1[kdVdZ]);
03845   lindu = swimstate1[kdUdZ];
03846   lindv = swimstate1[kdVdZ];
03847   linds = 0.;
03848   lindz = 0.;
03849   oldu = 0.;
03850   oldv = 0.;
03851   oldz = 0.;
03852   thisds = 0.;
03853   thisdz = 0.;
03854   nplaneskip = 0;
03855   swimflag = 1;
03856   befvtx = 1;
03857 
03858  
03859 
03860   while(thisplane>=0 && thisplane<=fParmLastPlane 
03861         && ((swimflag && TMath::Abs(thisx)<4. && TMath::Abs(thisy)<4. 
03862              && TMath::Abs(thisu)<4. && TMath::Abs(thisv)<4.) 
03863             || (TMath::Abs(linx)<4. && TMath::Abs(liny)<4. && TMath::Abs(linu)<4. 
03864                 && TMath::Abs(linv)<4.)
03865             )
03866         ){
03867 
03868   
03869    
03870     PlexPlaneId plexplane(fDetector,thisplane,false);
03871     if(plexplane.IsValid()){
03872       
03873       UgliScintPlnHandle planeid = ugh.GetScintPlnHandle(plexplane);
03874       
03875       if(planeid.IsValid()){
03876         oldu = swimstate1[kU];
03877         oldv = swimstate1[kV];
03878         oldz = swimstatez1;
03879         dz = (Double_t)planeid.GetZ0()-linz;
03880 
03881       
03882 
03883         if(TMath::Abs(linx)<4. && TMath::Abs(liny)<4. && TMath::Abs(linu)<4. 
03884            && TMath::Abs(linv)<4.){
03885           
03886           linz = (Double_t)planeid.GetZ0();
03887           linx += dz*lindx;
03888           liny += dz*lindy;
03889           linu += dz*lindu;
03890           linv += dz*lindv;
03891           if (!befvtx) {
03892             du = dz*lindu;
03893             dv = dz*lindv;
03894             lindz += TMath::Abs(dz);
03895             linds += TMath::Sqrt(du*du+dv*dv+dz*dz);
03896           }
03897         }//end if in the detector
03898 
03899         if(swimflag && TMath::Abs(thisx)<4. && TMath::Abs(thisy)<4. 
03900            && TMath::Abs(thisu)<4. && TMath::Abs(thisv)<4.){
03901           
03902           swimflag = Swim(swimstate1,tmpstate,(Double_t)swimstatez1,
03903                               (Double_t)planeid.GetZ0(),0, &fVldContext);
03904           swimstatez1 = (Double_t)planeid.GetZ0();
03905           swimstate1[kU] = tmpstate[kU];
03906           swimstate1[kV] = tmpstate[kV];
03907           swimstate1[kdUdZ] = tmpstate[kdUdZ];
03908           swimstate1[kdVdZ] = tmpstate[kdVdZ];
03909           swimstate1[kQoverP] = tmpstate[kQoverP];
03910           
03911           if (TMath::Abs(tmpstate[kU])>4. || TMath::Abs(tmpstate[kV])>4.) swimflag = 0;
03912           
03913           thisx = kInvSqrt2*(swimstate1[kU]-swimstate1[kV]);
03914           thisy = kInvSqrt2*(swimstate1[kU]+swimstate1[kV]);
03915           thisu = swimstate1[kU];
03916           thisv = swimstate1[kV];
03917           thisz = swimstatez1;
03918           
03919 
03920           
03921           if(!swimflag && nplaneskip<3 && thisplane*direction<=plane1*direction){
03922             cfh.SetEndU(thisu);
03923             cfh.SetEndV(thisv);
03924             cfh.SetEndZ(thisz);
03925             cfh.SetEndPlane(thisplane);
03926             dsdz = TMath::Sqrt(1.+swimstate1[kdUdZ]*swimstate1[kdUdZ]
03927                                         +swimstate1[kdVdZ]*swimstate1[kdVdZ]);
03928             cfh.SetEndDirCosU((Double_t)(direction)*swimstate1[kdUdZ]/dsdz);
03929             cfh.SetEndDirCosV((Double_t)(direction)*swimstate1[kdVdZ]/dsdz);
03930             cfh.SetEndDirCosZ((Double_t)(direction)/dsdz);
03931             cfh.SetU(thisplane,thisu);
03932             cfh.SetV(thisplane,thisv);
03933  
03934           } 
03935           else {
03936             befvtx = 0;
03937             delu = thisu-oldu;
03938             delv = thisv-oldv;
03939             delz = thisz-oldz;
03940             thisds += TMath::Sqrt(delu*delu+delv*delv+delz*delz);
03941             thisdz += TMath::Abs(delz);
03942           }
03943           if(planepe[thisplane] > 0.) nplaneskip = 0;
03944           else  ++nplaneskip;
03945            
03946         }//end if in the detector and swimflag
03947       }//end if planeid is valid
03948     }//end if plexplane is valid
03949     
03950     thisplane += direction;
03951   }//end loop over planes to find the end
03952   
03953   cfh.SetEndTrace(min(thisds,linds));
03954   cfh.SetEndTraceZ(min(thisdz,lindz));
03955   if(kpu1) cfh.SetEndExtrapolate(TMath::Abs(kpu1->GetTrackCluster()->GetPlane()-cfh.GetEndPlane()),
03956                                  PlaneView::kU);
03957   
03958   else cfh.SetEndExtrapolate(-1,PlaneView::kU);
03959   
03960   if(kpv1) cfh.SetEndExtrapolate(TMath::Abs(kpv1->GetTrackCluster()->GetPlane()-cfh.GetEndPlane()),
03961                                  PlaneView::kV); 
03962   else cfh.SetEndExtrapolate(-1,PlaneView::kV);
03963   
03964   return;
03965 }

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

Reimplemented from AlgBase.

Definition at line 1713 of file AlgFitTrackSR.cxx.

01714 {
01715 }


Member Data Documentation

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

Definition at line 84 of file AlgFitTrackSR.h.

Referenced by AddToFit(), InitializeFit(), and ReverseFit().

Definition at line 42 of file AlgFitTrackSR.h.

Referenced by FindNumSkippedPlanes(), RunAlg(), and SwimVertexAndEndPoints().

Definition at line 81 of file AlgFitTrackSR.h.

Referenced by InitializeFit(), and 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 64 of file AlgFitTrackSR.h.

Referenced by RemoveBadPointsFromFit(), and RunAlg().

Definition at line 65 of file AlgFitTrackSR.h.

Referenced by RemoveBadPointsFromFit(), and RunAlg().

Double_t AlgFitTrackSR::fParmDState[5] [private]

Definition at line 68 of file AlgFitTrackSR.h.

Referenced by CalculatePropagator(), and InitKalmanFitParameters().

Definition at line 78 of file AlgFitTrackSR.h.

Referenced by InitializeFit(), and RunAlg().

Definition at line 80 of file AlgFitTrackSR.h.

Referenced by InitializeFit(), and RunAlg().

Definition at line 79 of file AlgFitTrackSR.h.

Referenced by InitializeFit(), and RunAlg().

Definition at line 49 of file AlgFitTrackSR.h.

Referenced by AddToFit(), MakeDaughterStripList(), RunAlg(), and SetTrackParameters().

Definition at line 48 of file AlgFitTrackSR.h.

Referenced by RunAlg(), and SwimVertexAndEndPoints().

Double_t AlgFitTrackSR::fParmMaxAngle [private]

Definition at line 62 of file AlgFitTrackSR.h.

Referenced by MakeTrackClusterList(), and RunAlg().

Definition at line 53 of file AlgFitTrackSR.h.

Referenced by MarkTrackClusters(), and RunAlg().

Definition at line 44 of file AlgFitTrackSR.h.

Referenced by IterateKalmanFit(), and RunAlg().

Definition at line 45 of file AlgFitTrackSR.h.

Referenced by DoKalmanFit(), and RunAlg().

Definition at line 50 of file AlgFitTrackSR.h.

Referenced by RemoveBadPointsFromFit(), and RunAlg().

Definition at line 51 of file AlgFitTrackSR.h.

Referenced by AddForwardBestKPToFit(), AddReverseBestKPToFit(), and RunAlg().

Double_t AlgFitTrackSR::fParmMaxQP [private]
Double_t AlgFitTrackSR::fParmMaxQPFrac [private]

Definition at line 61 of file AlgFitTrackSR.h.

Referenced by MakeTrackClusterList(), and RunAlg().

Double_t AlgFitTrackSR::fParmMinPlanePE [private]

Definition at line 47 of file AlgFitTrackSR.h.

Referenced by RunAlg().

Definition at line 66 of file AlgFitTrackSR.h.

Referenced by MakeTrackClusterList(), and RunAlg().

Definition at line 55 of file AlgFitTrackSR.h.

Referenced by FindNumSkippedPlanesInView(), RemoveBadPointsFromFit(), and RunAlg().

Definition at line 59 of file AlgFitTrackSR.h.

Referenced by RunAlg().

Definition at line 58 of file AlgFitTrackSR.h.

Referenced by RunAlg().

Definition at line 57 of file AlgFitTrackSR.h.

Referenced by RunAlg().

Double_t AlgFitTrackSR::fParmPlnRadLen [private]

Definition at line 70 of file AlgFitTrackSR.h.

Referenced by CalculateNoise(), and InitKalmanFitParameters().

Double_t AlgFitTrackSR::fParmQPDiff [private]

Definition at line 46 of file AlgFitTrackSR.h.

Referenced by IterateKalmanFit(), and RunAlg().

Double_t AlgFitTrackSR::fParmTPosError2 [private]

Definition at line 77 of file AlgFitTrackSR.h.

Referenced by InitKalmanFitParameters(), and RunAlg().

Int_t AlgFitTrackSR::fSwimmer [private]

Definition at line 74 of file AlgFitTrackSR.h.

Referenced by CalculatePreState(), CalculatePropagator(), and InitKalmanFitParameters().

Definition at line 41 of file AlgFitTrackSR.h.

Referenced by FindNumSkippedPlanes(), RunAlg(), and SwimVertexAndEndPoints().

Int_t AlgFitTrackSR::pFor [private]

Definition at line 43 of file AlgFitTrackSR.h.

Referenced by AddToFit(), and RunAlg().

Int_t AlgFitTrackSR::pRev [private]

Definition at line 43 of file AlgFitTrackSR.h.

Referenced by ReverseFit(), and RunAlg().

Int_t AlgFitTrackSR::pRev2 [private]

Definition at line 43 of file AlgFitTrackSR.h.

Referenced by ReverseFit(), and RunAlg().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1