AlgFitTrackCam Class Reference

#include <AlgFitTrackCam.h>

Inheritance diagram for AlgFitTrackCam:
AlgBase AlgReco AlgTrack

List of all members.

Public Member Functions

 AlgFitTrackCam ()
virtual ~AlgFitTrackCam ()
virtual void RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
void InitialFramework (const CandSliceHandle *slice, CandContext &cx)
void RunTheFitter (CandFitTrackCamHandle &cth)
void StoreFilteredData (const int NewPlane)
void FillGapsInTrack ()
bool FindTheStrips (CandFitTrackCamHandle &cth, bool MakeTheTrack)
void GetFitData (int Plane1, int Plane2)
void ShowerStrips ()
void RemoveTrkHitsInShw ()
void ShowerSwim ()
void GoBackwards ()
void GoForwards ()
bool GetCombiPropagator (const int Plane, const int NewPlane, const bool GoForward)
bool Swim (double *StateVector, double *Output, const int Plane, const int NewPlane, const bool GoForward, double *dS=0, double *Range=0, double *dE=0)
bool Swim (double *StateVector, double *Output, const double zbeg, const int NewPlane, const bool GoForward, double *dS=0, double *Range=0, double *dE=0)
bool Swim (double *StateVector, double *Output, const int Plane, const double zend, const bool GoForward, double *dS=0, double *Range=0, double *dE=0)
void GetInitialCovarianceMatrix (const bool FirstIteration)
void GetNoiseMatrix (const int Plane, const int NewPlane)
void ExtrapCovMatrix ()
void CalcKalmanGain (const int NewPlane)
void UpdateStateVector (const int Plane, const int NewPlane, const bool GoForward)
void UpdateCovMatrix ()
void MoveArrays ()
void CheckValues (double *Input, const int NewPlane)
void SetTrackProperties (CandFitTrackCamHandle &cth)
void SetPropertiesFromFinderTrack (CandFitTrackCamHandle &cth)
void TimingFit (CandFitTrackCamHandle &cth)
void SetRangeAnddS (CandFitTrackCamHandle &cth)
void SpectrometerSwim (CandFitTrackCamHandle &cth)
void CleanNDLists (CandFitTrackHandle &cth, CandContext &cx)
bool NDPlaneIsActive (int plane, float u, float v, float projErr)
void GenerateNDSpectStrips (const CandSliceHandle *slice, CandContext &cx)
virtual void Trace (const char *c) const
void ResetCovarianceMatrix ()
double NDStripBegTime (CandStripHandle *Strip, double U=0, double V=0)
double MajorityCurvature ()
bool IsCoilReverse (VldContext *vldc)
bool IsStripInAnotherTrack (CandStripHandle *strip)
int IsInDetector (const double upos, const double vpos, const int plane)

Static Public Attributes

static std::vector
< CandStripHandle * > 
usedSlcStrips [490]

Private Attributes

vector< StripStructSlcStripData [490]
vector< StripStructInitTrkStripData [490]
vector< TrkDataStructTrkStripData [490]
vector< FiltDataStructFilteredData [490]
Int_t nbfield
Double_t bave
Bool_t EndofRange
Int_t EndofRangePlane
Bool_t LastIteration
double x_k4_biased
double eqp_biased
int UseGeoSwimmer
double x_k [5]
double x_k_minus [5]
double C_k [5][5]
double C_k_minus [5][5]
double C_k_intermediate [5][5]
double F_k [5][5]
double F_k_minus [5][5]
double Q_k [5][5]
double Q_k_minus [5][5]
double K_k [5]
int H_k [5]
int Identity [5][5]
double VtxCov [5]
double EndCov [5]
int MaxPlane
int MinPlane
double DeltaZ
double DeltaPlane
VldContextvldc
const CandTrackHandletrack
bool debug
bool ZIncreasesWithTime
bool PassTrack
bool SaveData
bool SwimThroughShower
int ShowerEntryPlane
int NIter
int TotalNSwimFail
int NumFinderStrips
double MeanTrackTime
PlaneOutline fPL
double StripListTime
double Beta2InvCoeff
double InitMaxqp
double InitMaxqpfrac
double InputThresholdEndofRange
double LastIterationMaxqp
int AllowOldChargeFlipCheck

Detailed Description

Definition at line 44 of file AlgFitTrackCam.h.


Constructor & Destructor Documentation

AlgFitTrackCam::AlgFitTrackCam (  ) 

Definition at line 80 of file AlgFitTrackCam.cxx.

00081 {
00082 }

AlgFitTrackCam::~AlgFitTrackCam (  )  [virtual]

Definition at line 87 of file AlgFitTrackCam.cxx.

00088 {
00089 }


Member Function Documentation

void AlgFitTrackCam::CalcKalmanGain ( const int  NewPlane  ) 

Definition at line 2103 of file AlgFitTrackCam.cxx.

References C_k_intermediate, H_k, K_k, Msg::kDebug, Msg::kVerbose, Munits::m, MSG, and TrkStripData.

Referenced by GoBackwards(), GoForwards(), ShowerSwim(), and SpectrometerSwim().

02104 {
02105   // K_k = C_k_intermediate * H_k^T * ( V_k + H_k * C_k_intermediate * H_k^T )^-1
02106   //  MSG("AlgFitTrackCam",Msg::kDebug) << "CalcKalmanGain" << endl;
02107 
02108   double Denominator=0;
02109 
02110   // H_k has only one non-zero element, so we can reduce matrix multiplication required
02111   if(TrkStripData[NewPlane][0].PlaneView==2) {Denominator=C_k_intermediate[0][0];}
02112   else {Denominator=C_k_intermediate[1][1];}
02113   
02114 
02115   // Add uncertainty in measurement
02116   Denominator+=TrkStripData[NewPlane][0].TPosError;
02117 
02118   MSG("AlgFitTrackCam",Msg::kVerbose) << "V_k " << TrkStripData[NewPlane][0].TPosError 
02119                                       << ", Kalman Gain Denominator " << Denominator << endl;
02120 
02121   if (Denominator!=0.) {
02122     for (int i=0; i<5; ++i) {
02123       K_k[i]=0;
02124       
02125       for (int m=0; m<5; ++m) {K_k[i]+=(C_k_intermediate[i][m])*(H_k[m]);}
02126       
02127       K_k[i]=K_k[i]/Denominator;
02128     }
02129 
02130     MSG("AlgFitTrackCam",Msg::kVerbose) << "Kalman Gain: " 
02131                                         << K_k[0] << " " << K_k[1] << " " << K_k[2] << " "
02132                                         << K_k[3] << " " << K_k[4] << endl;
02133   }
02134   else MSG("AlgFitTrackCam",Msg::kDebug) << "V_k + (H_k * C_k_intermediate * H_k_transpose) is zero!" << endl;
02135 }

void AlgFitTrackCam::CheckValues ( double *  Input,
const int  NewPlane 
)

Definition at line 2370 of file AlgFitTrackCam.cxx.

References EndofRange, CandTrackHandle::GetRange(), Msg::kVerbose, MSG, PassTrack, and track.

Referenced by UpdateStateVector().

02371 {
02372   // Make range and gradient corrections
02373   // Possible source of offset in q/p resolutions
02374 
02375   // Range check
02376     
02377   double Maxqp=4.; double Maxqpfrac=1.2;
02378   double Range=track->GetRange(NewPlane);
02379 
02380   //JAM signal end of range found
02381   // JAM 11/09 end of range changes from 100 to  250 MeV 
02382   //  if(fabs(Input[4])>10.0) EndofRange=true; 
02383 
02384   if(fabs(Input[4])>4.0) EndofRange=true; 
02385 
02386    if(Range>0. && (Maxqpfrac*500/Range)<Maxqp) {Maxqp=(Maxqpfrac*500/Range);}
02387   MSG("AlgFitTrackCam",Msg::kVerbose) << " Range " << Range << " Maxqp " << Maxqp << endl;
02388 
02389 
02390   // JAM 11/09 
02391   //  if(LastIteration) Maxqp=40;
02392 
02393   if(fabs(Input[4])>Maxqp){
02394     MSG("AlgFitTrackCam",Msg::kVerbose) << "CheckValues, Range check correction" << endl;
02395     Input[4]=(Input[4]>0 ? Maxqp : -Maxqp);
02396   }
02397   
02398   
02399   // Gradient check
02400   double Maxgradient=25.;
02401 
02402   if(fabs(Input[2])>Maxgradient) {
02403     MSG("AlgFitTrackCam",Msg::kVerbose) << "CheckValues, Gradient correction, U" << endl;
02404     Input[2]=(Input[2]>0 ? Maxgradient : -Maxgradient);
02405   }
02406 
02407   if(fabs(Input[3])>Maxgradient) {
02408     MSG("AlgFitTrackCam",Msg::kVerbose) << "CheckValues, Gradient correction, V" << endl;
02409     Input[3]=(Input[3]>0 ? Maxgradient : -Maxgradient);
02410   }
02411   
02412 
02413   // Check u and v values are not rubbish
02414   if(fabs(Input[0])<5000. && fabs(Input[1])<5000.) {PassTrack=true;}
02415   else {
02416     PassTrack=false;
02417   }
02418 }

void AlgFitTrackCam::CleanNDLists ( CandFitTrackHandle cth,
CandContext cx 
)

Definition at line 3682 of file AlgFitTrackCam.cxx.

References CandHandle::AddDaughterLink(), digit(), CandDigitListHandle::DupHandle(), CandStripListHandle::DupHandle(), CandRecord::FindCandHandle(), CandRecoHandle::GetCandSliceWritable(), CandHandle::GetDaughterIterator(), MomNavigator::GetFragment(), CandContext::GetMom(), CandHandle::IsSlushyEnabled(), CalDigitType::kNone, CandHandle::RemoveDaughter(), CandHandle::SetSlushyEnabled(), and SlcStripData.

Referenced by RunAlg().

03683 {
03684 
03685 
03686   // Sort out the lists for the ND spectrometer
03687 
03688   // Delete the strip handles created in GenerateNDSpectStrips.
03689   // All strip handles not added to a track daughter list are deleted here.
03691   for(int iplane=121; iplane<=290; ++iplane){
03692     for(unsigned int i=0; i<SlcStripData[iplane].size(); ++i){
03693       CandStripHandle* delstrip = SlcStripData[iplane][i].csh;
03694       delete delstrip; 
03695     }
03696   }
03698 
03699   bool SlushyOnEntry = CandHandle::IsSlushyEnabled();
03700   CandHandle::SetSlushyEnabled(kTRUE);
03701 
03702   // Get DigitList and StripList from CandRecord
03704   const MomNavigator* mom = cx.GetMom();
03705   CandRecord* crec = dynamic_cast<CandRecord *> (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
03706 
03707   // DupHandle step added by gmieg.  Must delete StripList and DigitList.
03708   CandStripListHandle* StripListOH = dynamic_cast<CandStripListHandle *>
03709     (crec->FindCandHandle("CandStripListHandle"));
03710   CandStripListHandle* StripList = StripListOH->DupHandle();
03711 
03712   CandDigitListHandle* DigitListOH = dynamic_cast<CandDigitListHandle *>
03713     (crec->FindCandHandle("CandDigitListHandle","canddigitlist"));
03714   CandDigitListHandle* DigitList = DigitListOH->DupHandle();
03715 
03716   CandSliceHandle* Slice = dynamic_cast<CandSliceHandle*>(cth.GetCandSliceWritable());
03718 
03719  
03720   // Compare new fitted track, with DeMuxed spectrometer strips, 
03721   // to StripList, Slice and DigitList
03723   vector<CandStripHandle*> StripsToAdd;
03724   vector<CandStripHandle*> StripsToRemove;
03725 
03726   vector<CandStripHandle*> SliceStripsToAdd;
03727   vector<CandStripHandle*> SliceStripsToRemove;
03728 
03729   vector<CandDigitHandle*> DigitsToAdd;
03730   vector<CandDigitHandle*> DigitsToRemove;
03731 
03732 
03733   CandStripHandleItr TrkStripItr(cth.GetDaughterIterator());
03734   for(CandStripHandle* TrkStrip=TrkStripItr(); TrkStrip; TrkStrip=TrkStripItr()){
03735 
03736     if(TrkStrip->GetPlane()>120 ) {
03737       CandDigitHandleItr TrkDigitItr(TrkStrip->GetDaughterIterator());
03738       CandDigitHandle* TrkDigit = dynamic_cast<CandDigitHandle*>(TrkDigitItr());
03739 
03740       // Sort out StripList
03742       if(StripList) {
03743         bool AddedTrkStrip = false;
03744         CandStripHandleItr stripItr(StripList->GetDaughterIterator());
03745         
03746         for(CandStripHandle* strip=stripItr(); strip ; strip=stripItr()) {
03747           if(strip->GetPlane()==TrkStrip->GetPlane()){
03748             CandDigitHandleItr digitItr(strip->GetDaughterIterator());
03749             CandDigitHandle* digit = dynamic_cast<CandDigitHandle*>(digitItr());
03750             bool SameCandStrip = false;
03751             bool SameHit = false;
03752             
03753             if(TrkStrip->IsEqual(strip)) {SameCandStrip = true;}
03754             
03755             if(digit->GetChannelId()==TrkDigit->GetChannelId() &&
03756                strip->GetCharge(CalDigitType::kNone)==TrkStrip->GetCharge(CalDigitType::kNone))
03757               {SameHit=true;}
03758             
03759             if(!SameCandStrip && SameHit) {
03760               StripsToRemove.push_back(strip);
03761               if(!AddedTrkStrip) {StripsToAdd.push_back(TrkStrip);}
03762               AddedTrkStrip=true;
03763             }
03764           }
03765         }
03766       }
03768       
03769 
03770 
03771       // Sort out Slice
03773       if(Slice){
03774         bool AddedTrkStrip = false;
03775         CandStripHandleItr stripItr(Slice->GetDaughterIterator());
03776 
03777         for(CandStripHandle* strip=stripItr(); strip; strip=stripItr()) {
03778           if(strip->GetPlane()==TrkStrip->GetPlane()){
03779             CandDigitHandleItr digitItr(strip->GetDaughterIterator());
03780             CandDigitHandle* digit = dynamic_cast<CandDigitHandle*>(digitItr());
03781             bool SameCandStrip = false;
03782             bool SameHit = false;
03783 
03784             if(TrkStrip->IsEqual(strip)) {SameCandStrip=true;}
03785 
03786             if(digit->GetChannelId()==TrkDigit->GetChannelId() &&
03787                strip->GetCharge(CalDigitType::kNone)==TrkStrip->GetCharge(CalDigitType::kNone))
03788               {SameHit=true;}
03789 
03790             if(!SameCandStrip && SameHit) {
03791               SliceStripsToRemove.push_back(strip);
03792               if(!AddedTrkStrip) {SliceStripsToAdd.push_back(TrkStrip);}
03793               AddedTrkStrip=true;
03794             }
03795           }
03796         }
03797       }
03799 
03800 
03801       // Loop over track strip Digits, and rationalise DigitList
03803       if(DigitList) {
03804         CandDigitHandleItr TrkDigitItr2(TrkStrip->GetDaughterIterator());
03805         for(TrkDigit=TrkDigitItr2(); TrkDigit ; TrkDigit=TrkDigitItr2()) {
03806           bool AddedTrkDigit=false;
03807           CandDigitHandleItr DigitItr(DigitList->GetDaughterIterator());
03808 
03809           for(CandDigitHandle* digit=DigitItr(); digit; digit=DigitItr()) {
03810             bool SameCandDigit=false;
03811             bool SameHit=false;
03812 
03813             if(TrkDigit->IsEqual(digit)) {SameCandDigit=true;}
03814 
03815             if(digit->GetChannelId()==TrkDigit->GetChannelId() &&
03816                digit->GetCharge(CalDigitType::kNone)==TrkDigit->GetCharge(CalDigitType::kNone))
03817               {SameHit=true;}
03818 
03819             if(!SameCandDigit && SameHit) {
03820               DigitsToRemove.push_back(digit);
03821               if(!AddedTrkDigit) {DigitsToAdd.push_back(TrkDigit);}
03822               AddedTrkDigit=true;
03823             }
03824           }
03825         }
03826       }
03828     
03830     }
03831   } // End loop over track strips
03832 
03833 
03834   // Now make the actual modifications to the lists
03836   for(unsigned int i=0; i<StripsToAdd.size(); ++i) {StripList->AddDaughterLink(*(StripsToAdd[i]));}
03837   for(unsigned int i=0; i<StripsToRemove.size(); ++i) {StripList->RemoveDaughter(StripsToRemove[i]);}
03838   StripsToAdd.clear();
03839   StripsToRemove.clear();
03840 
03841   for(unsigned int i=0; i<SliceStripsToAdd.size(); ++i) {Slice->AddDaughterLink(*(SliceStripsToAdd[i]));}
03842   for(unsigned int i=0; i<SliceStripsToRemove.size(); ++i) {Slice->RemoveDaughter(SliceStripsToRemove[i]);}
03843   SliceStripsToAdd.clear();
03844   SliceStripsToRemove.clear();
03845 
03846   for(unsigned int i=0; i<DigitsToAdd.size(); ++i) {DigitList->AddDaughterLink(*(DigitsToAdd[i]));}
03847   for(unsigned int i=0; i<DigitsToRemove.size(); ++i) {DigitList->RemoveDaughter(DigitsToRemove[i]);}
03848   DigitsToAdd.clear();
03849   DigitsToRemove.clear();
03851 
03852   
03854 
03855   // Must delete DupHandle StripList and Digitlist (gmieg)
03856   delete StripList;
03857   delete DigitList;
03858 
03859   if(!SlushyOnEntry) CandHandle::SetSlushyEnabled(kFALSE); 
03860 }

void AlgFitTrackCam::ExtrapCovMatrix (  ) 

Definition at line 2054 of file AlgFitTrackCam.cxx.

References C_k_intermediate, C_k_minus, debug, F_k_minus, Msg::kVerbose, Munits::m, MSG, and Q_k_minus.

Referenced by GoBackwards(), GoForwards(), ShowerSwim(), and SpectrometerSwim().

02055 {
02056   // C_k_intermediate = (F_k_minus * C_k_minus * F_k_minus^T) + Q_k_minus
02057   //  MSG("AlgFitTrackCam",Msg::kDebug) << "ExtrapCovMatrix" << endl;
02058 
02059   for (int i=0; i<5; ++i) {
02060     for (int j=0; j<5; ++j) {  
02061       C_k_intermediate[i][j]=0;
02062       
02063       for (int l=0; l<5; ++l) {
02064         for (int m=0; m<5; ++m) {   
02065           C_k_intermediate[i][j]+=F_k_minus[i][m]*C_k_minus[m][l]*F_k_minus[j][l];
02066         }
02067       }
02068 
02069       C_k_intermediate[i][j]+=Q_k_minus[i][j];
02070     }
02071     
02072   }
02073 
02074 
02075   // Diagonal elements should be positive
02076   double covlim = 1.e-8;
02077 
02078   for(int i=0; i<5; ++i) {
02079     if(C_k_intermediate[i][i]<covlim) {
02080       MSG("AlgFitTrackCam",Msg::kVerbose) << "Negative diagonal element in C_k_intermediate" << endl;
02081       C_k_intermediate[i][i]=covlim;
02082     }
02083   }
02084 
02085 
02086   // Display
02087   if(debug) {
02088     cout << "C_k_intermediate" << endl;
02089     for(int i=0; i<5; ++i) {
02090       for(int j=0; j<5; ++j) {
02091         cout << C_k_intermediate[i][j] << " ";
02092       }
02093       cout << endl;
02094     }
02095   }
02096   
02097 }

void AlgFitTrackCam::FillGapsInTrack (  ) 

Definition at line 895 of file AlgFitTrackCam.cxx.

References FilteredData, MaxPlane, MinPlane, size, SlcStripData, Swim(), FiltDataStruct::x_k0, FiltDataStruct::x_k1, FiltDataStruct::x_k2, FiltDataStruct::x_k3, FiltDataStruct::x_k4, and ZIncreasesWithTime.

Referenced by RunTheFitter().

00896 {
00897   // If there is no filtered data for a plane (between MinPlane and MaxPlane),
00898   // but this plane has hits in the slice, we interpolate from the nearest 
00899   // state vectors
00900   //
00901   // As with all filtered data, the interpolated data will be compared to
00902   // strip positions in the FindTheStrips method
00903   //  MSG("AlgFitTrackCam",Msg::kDebug) << "FillGapsInTrack" << endl;
00904 
00905 
00906   int CurrentPlane; int ForwardsPlane; int BackwardsPlane;
00907   int Plane; int NewPlane;  bool GoForward;
00908   double StateVector[5]; double Prediction[5]; bool GetPrediction;
00909   for(int i=0; i<5; i++) { Prediction[i]=0.; }
00910 
00911 
00912   for (int i=MinPlane; i<=MaxPlane; ++i) {   
00913     if(SlcStripData[i].size()>0) {
00914  
00915       if(FilteredData[i].size()==0) {
00916 
00917 
00918         // Find nearest filtered state vectors (within two planes) and ZPos differences
00920         // Forwards
00921         CurrentPlane=i+1; ForwardsPlane=-99;
00922 
00923         while(CurrentPlane<=MaxPlane && CurrentPlane<=(i+2)) {
00924           if(FilteredData[CurrentPlane].size()>0) {ForwardsPlane=CurrentPlane; break;}
00925           else {CurrentPlane++;}
00926         }
00927 
00928         // Backwards
00929         CurrentPlane=i-1; BackwardsPlane=-99;
00930 
00931         while(CurrentPlane>=MinPlane && CurrentPlane>=(i-2) ) {
00932           if(FilteredData[CurrentPlane].size()>0) {BackwardsPlane=CurrentPlane; break;}
00933           else {CurrentPlane--;}
00934         }
00936 
00937 
00938         // Find and store possible new filtered data, range and dS
00940         if(ForwardsPlane!=-99 && BackwardsPlane!=-99) {
00941 
00942           // Swimmer method
00943           GetPrediction=false;
00944           NewPlane=i;
00945           if(ZIncreasesWithTime==true) {Plane=ForwardsPlane; GoForward=false;}
00946           else{Plane=BackwardsPlane; GoForward=true;}
00947           if(FilteredData[Plane].size()>0) {
00948             StateVector[0] = FilteredData[Plane][0].x_k0;
00949             StateVector[1] = FilteredData[Plane][0].x_k1;
00950             StateVector[2] = FilteredData[Plane][0].x_k2;
00951             StateVector[3] = FilteredData[Plane][0].x_k3;
00952             StateVector[4] = FilteredData[Plane][0].x_k4;           
00953             GetPrediction=Swim(StateVector, Prediction, Plane, NewPlane, GoForward);
00954             
00955             if(GetPrediction==true) {
00956               // Store possible new state vector
00957               FiltDataStruct temp;
00958               temp.x_k0 = Prediction[0];
00959               temp.x_k1 = Prediction[1];
00960               temp.x_k2 = Prediction[2]; 
00961               temp.x_k3 = Prediction[3];
00962               temp.x_k4 = Prediction[4];              
00963               FilteredData[i].push_back(temp);
00964             }
00965           }
00966 
00967         }
00969       }
00970     }
00971   }
00972 
00973 }

bool AlgFitTrackCam::FindTheStrips ( CandFitTrackCamHandle cth,
bool  MakeTheTrack 
)

Definition at line 979 of file AlgFitTrackCam.cxx.

References CandHandle::AddDaughterLink(), StripStruct::csh, FilteredData, CandStripHandle::GetCharge(), InitTrkStripData, MaxPlane, MinPlane, CandHandle::RemoveDaughter(), size, SlcStripData, and ZIncreasesWithTime.

Referenced by RunTheFitter().

00980 { 
00981   // Given output data from Kalman filter, we find the strips that match most closely
00982   // and store the CandStripHandles in plane order in InitTrkStripData
00983   //
00984   // At end of track fitting, method is called with MakeTheTrack==true, so fit track
00985   // strips are finalised
00986 
00987   // JM ADDED 6/07:  Add pulse height requirement for first 3 hits at track vertex, and 
00988   //                 remove isolated single hit at vertex. This reduces tendency to add 
00989   //                 extra hits upstream of vertex.
00990   
00991   // Initialisations
00993   for (unsigned int i=0; i<490; ++i) {InitTrkStripData[i].clear();}
00994 
00995   double Extrem1=0; double Extrem2=0;
00996   double StripCentre=0;
00997   double StripWidth=4.108e-2;
00998   double HalfStripWidth=0.5*StripWidth;
00999   double TwoStripWidth=2.*StripWidth;  
01000   double PEthresh = 3.0;
01001   double MinDistanceToStrip=0;
01003   
01004   int NoOfHitPlanes=0; int Counter=0; int PlanesAdded=0;
01005   Bool_t foundMIP=false; 
01006   Bool_t RemovedHit=false;
01007   //  int InitMaxPlane = MaxPlane;
01008   float minDist;
01009   for (int i=MinPlane; i<=MaxPlane; ++i) {if(FilteredData[i].size()>0) {NoOfHitPlanes++;} }
01010 
01011   // Loop over the entire output from Kalman filter
01013   
01014   if(ZIncreasesWithTime==true){
01015     for (int i=MinPlane; i<=MaxPlane; ++i) {
01016       if(FilteredData[i].size()>0) {
01017         Counter++;
01018         int numStrips = 0;
01019         minDist=1e6;
01020         // Mark the possible extremities in transverse position within scintillator
01021         // by multiplying gradient by half scintillator thickness and adding or subtracting
01022         if(SlcStripData[i][0].csh->GetPlaneView()==2) {
01023           Extrem1=FilteredData[i][0].x_k0 + (0.0055*FilteredData[i][0].x_k2);
01024           Extrem2=FilteredData[i][0].x_k0 - (0.0055*FilteredData[i][0].x_k2);
01025         }
01026         else {
01027           Extrem1=FilteredData[i][0].x_k1 + (0.0055*FilteredData[i][0].x_k3);
01028           Extrem2=FilteredData[i][0].x_k1 - (0.0055*FilteredData[i][0].x_k3);
01029         }
01030         
01031         // Add strips in the case that only one strip can have its centre within 
01032         // half a strip width of an extremal position...
01034         bool PlaneAdded=false;
01035         if(fabs(Extrem1-Extrem2)<StripWidth) {
01036           for(unsigned int j=0; j<SlcStripData[i].size(); ++j) {
01037             StripCentre=SlcStripData[i][j].csh->GetTPos();
01038             MinDistanceToStrip=HalfStripWidth;
01039             if(fabs(StripCentre-Extrem1)<minDist)minDist =fabs(StripCentre-Extrem1);
01040             if(fabs(StripCentre-Extrem2)<minDist)minDist =fabs(StripCentre-Extrem2);
01041             if(fabs(StripCentre-Extrem1)<HalfStripWidth || fabs(StripCentre-Extrem2)<HalfStripWidth) {
01042               if(Counter>1 || SlcStripData[i][j].csh->GetCharge()>PEthresh || foundMIP==true || SlcStripData[i].size()==1){
01043                 foundMIP=true;
01044                 if(MakeTheTrack==true) {cth.AddDaughterLink(*SlcStripData[i][j].csh);}        
01045                 InitTrkStripData[i].push_back(SlcStripData[i][j]);
01046         
01047                 if(!PlaneAdded)PlanesAdded++;
01048                 PlaneAdded=true;
01049               } 
01050               
01051               numStrips++;
01052             }
01053           }
01054         }
01056         
01057         
01058         // ...Otherwise, cover the cases where multiple strips can have their centre
01059         // within half a strip width of an extremal position
01061         else {
01062           for(unsigned int j=0; j<SlcStripData[i].size(); ++j) {
01063             StripCentre=SlcStripData[i][j].csh->GetTPos();          
01064             MinDistanceToStrip=HalfStripWidth;      
01065             if(fabs(StripCentre-Extrem1)<minDist)minDist =fabs(StripCentre-Extrem1);
01066             if(fabs(StripCentre-Extrem2)<minDist)minDist =fabs(StripCentre-Extrem2);
01067 
01068             if( fabs(StripCentre-Extrem1)<HalfStripWidth || fabs(StripCentre-Extrem2)<HalfStripWidth
01069                 || (Extrem1>StripCentre && Extrem2<StripCentre) 
01070                 || (Extrem1<StripCentre && Extrem2>StripCentre) ) {
01071               if(Counter>1 || SlcStripData[i][j].csh->GetCharge()>PEthresh || foundMIP==true){
01072                 foundMIP=true;
01073                 if(MakeTheTrack==true) {cth.AddDaughterLink(*SlcStripData[i][j].csh);}
01074                 InitTrkStripData[i].push_back(SlcStripData[i][j]);
01075                 if(!PlaneAdded)PlanesAdded++;
01076                 PlaneAdded=true;
01077               }
01078               numStrips++;            
01079             }
01080           }
01081         }
01083         // If we have found no strips, we consider looking further, finding the closest  
01084         // strip within a distance 'MinDistanceToStrip' of an extremal position
01086         if(numStrips==0) {
01087           CandStripHandle* CurrentStrip=0;
01088           
01089           // Be more demanding near track vertex
01090           if(Counter>2) {
01091             MinDistanceToStrip = TwoStripWidth;
01092             if(MakeTheTrack==true) MinDistanceToStrip = 2*TwoStripWidth;
01093           }
01094           else {MinDistanceToStrip=StripWidth;}
01095           for(unsigned int j=0; j<SlcStripData[i].size(); ++j) {
01096             StripCentre=SlcStripData[i][j].csh->GetTPos();
01097             
01098             // Find the closest strip and temporarily store its CandStripHandle
01099             if(fabs(StripCentre-Extrem1)<minDist)minDist =fabs(StripCentre-Extrem1);
01100             if(fabs(StripCentre-Extrem2)<minDist)minDist =fabs(StripCentre-Extrem2);
01101             if(fabs(StripCentre-Extrem1)<MinDistanceToStrip) {
01102               MinDistanceToStrip=StripCentre-Extrem1;
01103               CurrentStrip=SlcStripData[i][j].csh;
01104             }
01105             if(fabs(StripCentre-Extrem2)<MinDistanceToStrip) {
01106               MinDistanceToStrip=StripCentre-Extrem2;
01107               CurrentStrip=SlcStripData[i][j].csh;
01108             }
01109           }
01110           
01111           // If we have found a strip then we add it
01112           if(CurrentStrip) {
01113             if(Counter>1 || CurrentStrip->GetCharge()>PEthresh || foundMIP==true){
01114               foundMIP=true;
01115               if(MakeTheTrack==true) {cth.AddDaughterLink(*CurrentStrip);}          
01116               StripStruct temp;
01117               temp.csh = CurrentStrip;
01118               InitTrkStripData[i].push_back(temp);
01119               PlanesAdded++;
01120               PlaneAdded=true;
01121             }
01122           }
01123         }
01124       }
01125       if(PlanesAdded==3 && !RemovedHit ){  // remove first hit if it is separated by >1 plane from second
01126         Int_t np = 0;
01127         Int_t planebuf[3];
01128         Int_t pln = MinPlane;
01129         while(np<3 && pln<490){
01130           if(InitTrkStripData[pln].size()>0){
01131             planebuf[np] = InitTrkStripData[pln][0].csh->GetPlane();
01132             np++;
01133           }
01134           pln++;
01135         }
01136         if(np==3 && SlcStripData[planebuf[0]].size()>1){
01137           if(abs(planebuf[1]-planebuf[0])/abs(planebuf[2]-planebuf[1])>1.5){
01138             for(unsigned int j=0; j<InitTrkStripData[planebuf[0]].size(); ++j) {
01139               CandStripHandle* Strip = InitTrkStripData[planebuf[0]][j].csh;     
01140               cth.RemoveDaughter(Strip);
01141               RemovedHit=true;
01142             }
01143             InitTrkStripData[planebuf[0]].clear();
01144           }
01145         }
01146       }
01147     }
01148   }
01149   else{
01150     for (int i=MaxPlane; i>=MinPlane; --i) {   
01151       if(FilteredData[i].size()>0) {
01152         Counter++;
01153         int numStrips=0;
01154         
01155         
01156         // Mark the possible extremities in transverse position within scintillator
01157         // by multiplying gradient by half scintillator thickness and adding or subtracting
01158         if(SlcStripData[i][0].csh->GetPlaneView()==2) {
01159           Extrem1=FilteredData[i][0].x_k0 + (0.0055*FilteredData[i][0].x_k2);
01160           Extrem2=FilteredData[i][0].x_k0 - (0.0055*FilteredData[i][0].x_k2);
01161         }
01162         else {
01163           Extrem1=FilteredData[i][0].x_k1 + (0.0055*FilteredData[i][0].x_k3);
01164           Extrem2=FilteredData[i][0].x_k1 - (0.0055*FilteredData[i][0].x_k3);
01165         }
01166         
01167         
01168         // Add strips in the case that only one strip can have its centre within 
01169         // half a strip width of an extremal position...
01171         bool PlaneAdded=false;
01172         if(fabs(Extrem1-Extrem2)<StripWidth) {
01173           for(unsigned int j=0; j<SlcStripData[i].size(); ++j) {
01174             StripCentre=SlcStripData[i][j].csh->GetTPos();
01175             
01176             if(fabs(StripCentre-Extrem1)<HalfStripWidth || fabs(StripCentre-Extrem2)<HalfStripWidth) {
01177               if(Counter>1 || SlcStripData[i][j].csh->GetCharge()>PEthresh || foundMIP|| SlcStripData[i].size()==1){
01178                 foundMIP=true;
01179                 if(MakeTheTrack==true) {cth.AddDaughterLink(*SlcStripData[i][j].csh);}
01180                 InitTrkStripData[i].push_back(SlcStripData[i][j]);
01181                 numStrips++;
01182                 if(!PlaneAdded)PlanesAdded++;
01183                 PlaneAdded=true;
01184               }      
01185             }
01186           }
01187         }
01189         
01190         
01191         // ...Otherwise, cover the cases where multiple strips can have their centre
01192         // within half a strip width of an extremal position
01194         else {
01195           bool PlaneAdded=false;
01196           for(unsigned int j=0; j<SlcStripData[i].size(); ++j) {
01197             StripCentre=SlcStripData[i][j].csh->GetTPos();
01198             
01199             if( fabs(StripCentre-Extrem1)<HalfStripWidth || fabs(StripCentre-Extrem2)<HalfStripWidth
01200                 || (Extrem1>StripCentre && Extrem2<StripCentre) 
01201                 || (Extrem1<StripCentre && Extrem2>StripCentre) ) {
01202               if(Counter>1 || SlcStripData[i][j].csh->GetCharge()>PEthresh || foundMIP){
01203                 foundMIP=true;
01204                 if(MakeTheTrack==true) {cth.AddDaughterLink(*SlcStripData[i][j].csh);}          
01205                 InitTrkStripData[i].push_back(SlcStripData[i][j]);
01206                 numStrips++;
01207                 if(!PlaneAdded)PlanesAdded++;
01208                 PlaneAdded=true;
01209               }
01210             }
01211           }
01212         }
01214         
01215         // If we have found no strips, we consider looking further, finding the closest  
01216         // strip within a distance 'MinDistanceToStrip' of an extremal position
01218         if(numStrips==0) {
01219           CandStripHandle* CurrentStrip=0;
01220           
01221           // Be more demanding near track vertex
01222           if(Counter>2) {
01223             MinDistanceToStrip = TwoStripWidth;
01224             if(MakeTheTrack==true) MinDistanceToStrip = 2*TwoStripWidth;
01225           }
01226           else {MinDistanceToStrip=StripWidth;}
01227           
01228           for(unsigned int j=0; j<SlcStripData[i].size(); ++j) {
01229             StripCentre=SlcStripData[i][j].csh->GetTPos();
01230             
01231             // Find the closest strip and temporarily store its CandStripHandle
01232             if(fabs(StripCentre-Extrem1)<MinDistanceToStrip) {
01233               MinDistanceToStrip=StripCentre-Extrem1;
01234               CurrentStrip=SlcStripData[i][j].csh;
01235             }
01236             if(fabs(StripCentre-Extrem2)<MinDistanceToStrip) {
01237               MinDistanceToStrip=StripCentre-Extrem2;
01238               CurrentStrip=SlcStripData[i][j].csh;
01239             }
01240           }
01241           
01242           // If we have found a strip then we add it
01243           if(CurrentStrip) {
01244             if(Counter>1 || CurrentStrip->GetCharge()>PEthresh || foundMIP){
01245               foundMIP=true;    
01246               if(MakeTheTrack==true) {cth.AddDaughterLink(*CurrentStrip);}
01247               StripStruct temp;
01248               temp.csh = CurrentStrip;
01249               InitTrkStripData[i].push_back(temp);
01250               PlanesAdded++;
01251               PlaneAdded=true;
01252             }
01253           }
01254         }
01256       }
01257      if(PlanesAdded==3 && !RemovedHit){  // remove first hit if it is separated by >1 plane from second
01258         Int_t np = 0;
01259         Int_t planebuf[3];
01260         Int_t pln = MaxPlane;
01261         while(np<3 && pln>=0){
01262           if(InitTrkStripData[pln].size()>0){
01263             planebuf[np]= InitTrkStripData[pln][0].csh->GetPlane();
01264             np++;
01265           }
01266           pln--;
01267         }
01268         if(np==3 && SlcStripData[planebuf[0]].size()>1){
01269           if(abs(planebuf[1]-planebuf[0])/abs(planebuf[2]-planebuf[1])>1.5){
01270           for(unsigned int j=0; j<InitTrkStripData[planebuf[0]].size(); ++j) {
01271             CandStripHandle* Strip = InitTrkStripData[planebuf[0]][j].csh;       
01272             cth.RemoveDaughter(Strip);
01273             RemovedHit=true;
01274           }
01275           InitTrkStripData[planebuf[0]].clear();
01276           }
01277         }
01278       }
01279     }  
01280   }
01282   
01283 
01284   // Find new max and min planes
01285   MaxPlane=-20; MinPlane=500;
01286   for (int i=0; i<490; ++i) {   
01287     if(InitTrkStripData[i].size()>0) {
01288       if(i>MaxPlane) {MaxPlane=i;}
01289       if(i<MinPlane) {MinPlane=i;}
01290     }
01291   }
01292 
01293   if(MaxPlane==-20 || MinPlane==500) {return false;}
01294   else {return true;}
01295 
01296 }

void AlgFitTrackCam::GenerateNDSpectStrips ( const CandSliceHandle slice,
CandContext cx 
)

Definition at line 3359 of file AlgFitTrackCam.cxx.

References PlexSEIdAltL::AddToCurrentWeight(), StripStruct::csh, PlexSEIdAltL::GetBestWeight(), CandHandle::GetCandRecord(), CandDigitHandle::GetCharge(), CandHandle::GetDaughterIterator(), AlgFactory::GetInstance(), CandContext::GetMom(), CandDigitHandle::GetPlexSEIdAltLWritable(), Msg::kDebug, CalDigitType::kPE, PlaneView::kU, PlaneView::kV, CandStrip::MakeCandidate(), MSG, NumFinderStrips, CandHandle::SetCandRecord(), and SlcStripData.

Referenced by InitialFramework().

03360 {
03361   MSG("AlgFitTrackCam",Msg::kDebug) << "GenerateNDSpectStrips" << endl;
03362 
03363   bool AlreadyAssigned;
03364 
03365 
03366   CandContext cxx(this,cx.GetMom());
03367 
03368   // Get singleton instance of AlgFactory
03369   AlgFactory &af = AlgFactory::GetInstance();
03370   AlgHandle stripah = af.GetAlgHandle("AlgStripSR","default");
03371 
03372   // Store CandStripHandles and make the strips accessible by plane number
03373   TIter SlcStripItr = slice->GetDaughterIterator();
03374   StripStruct temp;
03375   
03376   // Store all strips in slice
03377   while(CandStripHandle* SlcStrip = dynamic_cast<CandStripHandle*>(SlcStripItr()))  {
03378     if (SlcStrip->GetPlane()>120 && (SlcStrip->GetPlaneView()==PlaneView::kU || SlcStrip->GetPlaneView()==PlaneView::kV) ) {
03379       AlreadyAssigned=false;
03380       
03381       TIter digitItr(SlcStrip->GetDaughterIterator());
03382       CandDigitHandle* dig = dynamic_cast<CandDigitHandle*>(digitItr());              
03383       const PlexSEIdAltL& altl = dig->GetPlexSEIdAltL();
03384       int siz = altl.size();
03385       
03386       for (int ind = 0; ind<siz; ++ind) {
03387         // Only want to make the single copy of an assigned strip
03388         if(AlreadyAssigned) {break;}
03389         
03390         TObjArray diglist;  
03391         TIter newstripdauItr(SlcStrip->GetDaughterIterator());
03392         
03393         while(CandDigitHandle* olddig = dynamic_cast<CandDigitHandle*>(newstripdauItr())) {
03394           CandDigitHandle* newdig = olddig->DupHandle();
03395           PlexSEIdAltL& newaltl = newdig->GetPlexSEIdAltLWritable(); 
03396           
03397           // Don't make any strips which have already been assigned to a 'better' track
03398           if(NumFinderStrips<=newaltl.GetBestWeight()) {AlreadyAssigned=true; delete newdig;}
03399           else {
03400             for (int newind = 0; newind < siz; ++newind) {
03401               if(newind==ind){newaltl[newind].SetWeight((float)NumFinderStrips);}
03402               else{newaltl[newind].SetWeight((float)0.);}
03403             }
03404             newaltl.AddToCurrentWeight(0.); //force plex to invalidate cacheing
03405             newdig->SetCandRecord(olddig->GetCandRecord());
03406             diglist.Add(newdig);   // diglist does not own newdig. These must be explicitly deleted 
03407           }
03408         
03409         }
03410   // Only make a new strip if we have any digits and it is above the charge threshold
03411   float tCharge = 0.0;
03412   for(int d = 0; d < diglist.GetEntries(); ++d){
03413     CandDigitHandle *tDig = dynamic_cast<CandDigitHandle*>(diglist.At(d));
03414     tCharge += tDig->GetCharge(CalDigitType::kPE);
03415   }
03416   if(tCharge > 0.0) {
03417           cxx.SetDataIn(&diglist);
03418           CandStripHandle newstrip = CandStrip::MakeCandidate(stripah,cxx);
03419           newstrip.SetCandRecord(slice->GetCandRecord());
03420           CandStripHandle* daustrip = new CandStripHandle(newstrip);
03421           
03422           for (int i=0; i<=diglist.GetLast(); ++i) {
03423             CandDigitHandle* deldig =  dynamic_cast<CandDigitHandle*>(diglist.At(i));
03424             delete deldig;
03425           }
03426           temp.csh=daustrip;      
03427           SlcStripData[SlcStrip->GetPlane()].push_back(temp);
03428         }
03429       }
03430     }
03431   }
03432   SlcStripItr.Reset();
03433 }

bool AlgFitTrackCam::GetCombiPropagator ( const int  Plane,
const int  NewPlane,
const bool  GoForward 
)

Definition at line 1472 of file AlgFitTrackCam.cxx.

References debug, DeltaPlane, DeltaZ, F_k_minus, IsCoilReverse(), Msg::kDebug, Msg::kVerbose, MSG, Swim(), TotalNSwimFail, TrkStripData, vldc, and x_k_minus.

Referenced by GoBackwards(), GoForwards(), ShowerSwim(), and SpectrometerSwim().

01473 {
01474   // Combination propagator, essentially same as SR propagator, but
01475   // generation of matrix reduces calls to swimmer by 80%
01476   //  MSG("AlgFitTrackCam",Msg::kDebug) << "GetCombiPropagator" << endl;
01477 
01478   for (int i=0; i<5; ++i) {
01479     for (int j=0; j<5; ++j) {
01480       F_k_minus[i][j]=0;  
01481     }
01482   }
01483 
01484   F_k_minus[0][0]=1; F_k_minus[1][1]=1; 
01485   F_k_minus[2][2]=1; F_k_minus[3][3]=1;
01486   
01487   DeltaZ=fabs(TrkStripData[NewPlane][0].ZPos-TrkStripData[Plane][0].ZPos);
01488   DeltaPlane=abs(NewPlane-Plane);
01489 
01490   // Swimmer section for last column
01491   double PState[5];  double NState[5];  double StateVector[5];
01492   double Increment=0.01;
01493   bool SwimInc=false; bool SwimDec=false;
01494   int nswimfail=0;
01495   
01496   if(GoForward==true) {F_k_minus[0][2]=DeltaZ; F_k_minus[1][3]=DeltaZ;}
01497   else if(GoForward==false) {F_k_minus[0][2]=-DeltaZ; F_k_minus[1][3]=-DeltaZ;}
01498   
01499   
01500   // Give swimmer fixed number of opportunities for successful swim
01501   while((SwimInc==false || SwimDec==false) && nswimfail<=10) {
01502     
01503     // AB: 10th June, 2013
01504     // Use the polarity of magnetic field to set the direction of the 
01505     // first d(Q/P) increment. This section of code should now possess 
01506     // the following symmetry: B -> -B implies Q/P -> -Q/P 
01507 
01508     // The nominal polarity
01509     double bfld_polarity = +1.0; 
01510 
01511     // This method will be called many times. Hopefully it's efficient!
01512     if( this->IsCoilReverse(vldc) ) bfld_polarity = -1.0;
01513 
01514     Increment=0.05*fabs(x_k_minus[4]); 
01515     if(Increment == 0.) { Increment=bfld_polarity*0.0001; }
01516     else if(Increment<.0001) { Increment=.0001; }
01517        
01518     for(int j=0; j<5; ++j) {StateVector[j]=x_k_minus[j];}  
01519     
01520     // Increment then swim
01521     StateVector[4]+=Increment;
01522     SwimInc=Swim(StateVector, NState, Plane, NewPlane, GoForward);
01523     
01524     StateVector[4]=x_k_minus[4];
01525     
01526     // Decrement then swim
01527     StateVector[4]-=Increment;
01528     SwimDec=Swim(StateVector, PState, Plane, NewPlane, GoForward);
01529     
01530     // If swim failed, double momentum and swim again
01531     if(SwimInc==false || SwimDec==false) {
01532       MSG("AlgFitTrackCam",Msg::kVerbose) << "GetCombiPropagator, Swim failed - Double momentum and swim again" << endl;
01533       x_k_minus[4]*=0.5;
01534       nswimfail++; TotalNSwimFail++;
01535       break;
01536     }
01537     
01538     // Form last row of propagator matrix.  Need to transpose to get proper Kalman F_k_minus
01539     else {
01540       if(Increment!=0.) {
01541         for(int j=0; j<5; ++j) {
01542           F_k_minus[j][4] = (NState[j]-PState[j]) / (2.0*fabs(Increment));
01543         }
01544       }
01545       else {F_k_minus[4][4]=1;}
01546     }
01547     
01548   } // End while statement
01549   
01550   if(nswimfail>10) {MSG("AlgFitTrackCam",Msg::kDebug) << "GetCombiPropagator, nswimfail>10, fail track" << endl; return false;}
01551 
01552 
01553   // Display
01554   if(debug) {
01555     cout << "Combi F_k_minus" << endl;
01556     for(int i=0; i<5; ++i) {  
01557       for(int j=0; j<5; ++j) {
01558         cout << F_k_minus[i][j] << " ";
01559       }
01560       cout << endl;
01561     }
01562   }
01563   
01564   return true;
01565 }

void AlgFitTrackCam::GetFitData ( int  Plane1,
int  Plane2 
)

Definition at line 800 of file AlgFitTrackCam.cxx.

References MuELoss::e, CandTrackHandle::GetU(), CandTrackHandle::GetV(), InitTrkStripData, PlaneView::kV, TrkDataStruct::PlaneView, size, SlcStripData, TrkDataStruct::TPos, TrkDataStruct::TPosError, track, TrkStripData, and TrkDataStruct::ZPos.

Referenced by RunAlg(), RunTheFitter(), SetTrackProperties(), ShowerSwim(), and SpectrometerSwim().

00801 {
00802   // Loop over the initial track strip data and create the final data for fitting
00803   //  MSG("AlgFitTrackCam",Msg::kDebug) << "GetFitData" << endl;
00804 
00805   // Initialisations
00806   double MisalignmentError=4e-6;  // Squared error for misalignment of strips
00807   double SumChargeTPos=0; double SumCharge=0; int MaxStrip=-20; int MinStrip=200;
00808   bool NewStripFound=true;
00809 
00810   int ThisStrip;
00811 
00812   // Get the data for region between the planes specified
00813   for(int i=Plane1; i<=Plane2; ++i) {
00814     if(InitTrkStripData[i].size()>0) {     
00815       
00816       // Find max and min strip numbers for strips in plane
00817       for(unsigned int j=0; j<InitTrkStripData[i].size(); ++j) {
00818         if(InitTrkStripData[i][j].csh->GetStrip()<MinStrip) {MinStrip=InitTrkStripData[i][j].csh->GetStrip();}
00819         if(InitTrkStripData[i][j].csh->GetStrip()>MaxStrip) {MaxStrip=InitTrkStripData[i][j].csh->GetStrip();}  
00820       }
00821 
00822 
00823       // Find continuous groups of strips
00825       NewStripFound=true;
00826       
00827       while(NewStripFound==true) {
00828         
00829         NewStripFound=false;
00830         
00831         for(unsigned int j=0; j<SlcStripData[i].size(); ++j) {
00832 
00833           ThisStrip=SlcStripData[i][j].csh->GetStrip();
00834 
00835           if( ThisStrip==(MaxStrip+1) ) {
00836             MaxStrip+=1;
00837             NewStripFound=true;
00838             
00839             InitTrkStripData[i].push_back(SlcStripData[i][j]);
00840           }
00841           
00842           if( ThisStrip==(MinStrip-1) ) {
00843             MinStrip-=1;
00844             NewStripFound=true;
00845             
00846             InitTrkStripData[i].push_back(SlcStripData[i][j]);
00847           }
00848         }
00849       }
00851 
00852 
00853 
00854       // Get the data for fitting
00856       for(unsigned int j=0; j<InitTrkStripData[i].size(); ++j) {
00857         SumCharge+=InitTrkStripData[i][j].csh->GetCharge();
00858 
00859         // JAM 11/08:  Change to support strip rotations
00860         PlaneView::PlaneView_t plnvw =InitTrkStripData[i][j].csh->GetPlaneView();
00861         Double_t orthopos = track->GetV(InitTrkStripData[i][j].csh->GetPlane());
00862         if(plnvw == PlaneView::kV){
00863           orthopos = track->GetU(InitTrkStripData[i][j].csh->GetPlane());
00864         }
00865         // JAM 11/08:  to activate strip rotation support add orthopos as argument to GetTPos()
00866         SumChargeTPos+=InitTrkStripData[i][j].csh->GetCharge()*InitTrkStripData[i][j].csh->GetTPos();
00867       }
00868 
00869       // Charge weighted TPos and Flat distribution error
00870       if(SumCharge!=0.) {
00871         TrkDataStruct temp;
00872         
00873         temp.ZPos=InitTrkStripData[i][0].csh->GetZPos();
00874         temp.PlaneView=InitTrkStripData[i][0].csh->GetPlaneView();
00875 
00876         temp.TPos=SumChargeTPos/SumCharge;
00877         temp.TPosError=( (1.406305333e-4 * (1 + MaxStrip-MinStrip) )+ MisalignmentError);
00878         
00879         TrkStripData[i].push_back(temp);
00880       }
00882 
00883 
00884       // Reset      
00885       SumChargeTPos=0; SumCharge=0; MaxStrip=-20; MinStrip=200;
00886     }
00887   }
00888   
00889 }

void AlgFitTrackCam::GetInitialCovarianceMatrix ( const bool  FirstIteration  ) 

Definition at line 1868 of file AlgFitTrackCam.cxx.

References C_k_minus, debug, min, and ZIncreasesWithTime.

Referenced by ResetCovarianceMatrix(), and RunTheFitter().

01869 {
01870   //  MSG("AlgFitTrackCam",Msg::kDebug) << "GetInitialCovarianceMatrix" << endl;
01871 
01872   if(FirstIteration==true) {
01873     
01874     for(int i=0; i<5; ++i) {
01875       for(int j=0; j<5; ++j) {
01876         C_k_minus[i][j]=0.;
01877       }
01878     }
01879     
01880     // Diagonal terms
01881     C_k_minus[0][0]=0.25; C_k_minus[1][1]=0.25; 
01882     C_k_minus[2][2]=100.; C_k_minus[3][3]=100.; 
01883     C_k_minus[4][4]=1.;
01884     
01885     // Off diagonal terms. Taken from SR - Origin of this?
01886     if(ZIncreasesWithTime==true) {
01887       C_k_minus[0][4]=7.5e-5;  C_k_minus[1][4]=7.5e-5;
01888       C_k_minus[4][0]=7.5e-5;  C_k_minus[4][1]=7.5e-5;
01889     }
01890     else if(ZIncreasesWithTime==false) {
01891       C_k_minus[0][4]=-7.5e-5;  C_k_minus[1][4]=-7.5e-5;
01892       C_k_minus[4][0]=-7.5e-5;  C_k_minus[4][1]=-7.5e-5;
01893     }
01894     
01895     
01896   }
01897 
01898   else if(FirstIteration==false) {
01899     // Results are very sensitive to this multiplication. A large number means
01900     // that further iterations start with the same uncertainties as the first,
01901     // albeit with improved "track finder" strips
01902     for(int i=0; i<5; ++i) {C_k_minus[i][i]*=100;}
01903     
01904     
01905     // Make sure not larger than very first covariance elements
01906     C_k_minus[0][0]=min(C_k_minus[0][0],0.25); C_k_minus[1][1]=min(C_k_minus[1][1],0.25);
01907     C_k_minus[2][2]=min(C_k_minus[2][2],100.); C_k_minus[3][3]=min(C_k_minus[3][3],100.);
01908     C_k_minus[4][4]=min(C_k_minus[4][4],1.);
01909 
01910     double cov_xqp = 7.5e-5;             // Taken from SR - Origin of this?
01911     
01912     for(int i=0; i<2; ++i){
01913       if(fabs(C_k_minus[i][4])>cov_xqp) C_k_minus[i][4] = (C_k_minus[i][4] > 0 ? cov_xqp : -cov_xqp);
01914       if(fabs(C_k_minus[4][i])>cov_xqp) C_k_minus[4][i] = (C_k_minus[4][i] > 0 ? cov_xqp : -cov_xqp);
01915     }
01916 
01917     cov_xqp /= 0.06;                     // Taken from SR - Origin of this?
01918 
01919     for(int i=2; i<4; ++i){
01920       if(fabs(C_k_minus[i][4])>cov_xqp) C_k_minus[i][4] = (C_k_minus[i][4] > 0 ? cov_xqp : -cov_xqp);
01921       if(fabs(C_k_minus[4][i])>cov_xqp) C_k_minus[4][i] = (C_k_minus[4][i] > 0 ? cov_xqp : -cov_xqp);
01922     }
01923   }
01924 
01925 
01926   // Display
01927   if(debug) {
01928     cout << "Initial covariance matrix" << endl;
01929     for(int p=0; p<5; ++p){
01930       for(int q=0; q<5; ++q){
01931         cout << C_k_minus[p][q] << " ";
01932       }  
01933       cout << endl;
01934     }
01935   }
01936 
01937 }

void AlgFitTrackCam::GetNoiseMatrix ( const int  Plane,
const int  NewPlane 
)

Definition at line 1943 of file AlgFitTrackCam.cxx.

References debug, DeltaPlane, DeltaZ, BField::GetBField(), UgliGeomHandle::GetNearestSteelPlnHandle(), SwimGeo::GetSteel(), SwimPlaneInterface::GetZ(), UgliSteelPlnHandle::GetZ0(), SwimGeo::Instance(), Q_k_minus, TrkStripData, vldc, and x_k_minus.

Referenced by GoBackwards(), GoForwards(), ShowerSwim(), and SpectrometerSwim().

01944 {
01945   // This method is essentially the same as in SR fitter
01946   //  MSG("AlgFitTrackCam",Msg::kDebug) << "GetNoiseMatrix" << endl;
01947 
01948   for (int p=0; p<5; ++p) {
01949     for (int q=0; q<5; ++q) {
01950       Q_k_minus[p][q]=0;    }
01951   }
01952   
01953   // Get gradients, etc from x_k_minus
01954   double dsdzSquared = 1.+pow(x_k_minus[2],2)+pow(x_k_minus[3],2);
01955   double dsdz = pow(dsdzSquared,0.5);
01956 
01957   // Implement noise matrix as in SR
01958   if (DeltaPlane!=-99 && DeltaZ!=-99) {  
01959     double qp = x_k_minus[4];
01960     if(fabs(qp)<0.01) qp = (qp>0 ? 0.01 : -0.01);
01961     int izdir = ((NewPlane-Plane)>0 ? 0 : 1);
01962   
01963     double LocalRadiationLength=(dsdz * double(DeltaPlane) * 1.47); // 1.47 radiation lengths per iron plane
01964 
01965     double SigmaMS=(0.0136 * fabs(qp) * pow(LocalRadiationLength,0.5)
01966                     * (1. + (0.038 * log(LocalRadiationLength)) ));
01967     double SigmaMSSquared=pow(SigmaMS,2);
01968 
01969     double Sigma33Squared=0.5*SigmaMSSquared*dsdzSquared*(1.+pow(x_k_minus[2],2));
01970 
01971     double Sigma34Squared=0.5*SigmaMSSquared*dsdzSquared*(x_k_minus[2]*x_k_minus[3]);
01972 
01973     double Sigma44Squared=0.5*SigmaMSSquared*dsdzSquared*(1.+pow(x_k_minus[3],2));;
01974 
01975 
01976     // Treat steel as discrete scatter
01977     SwimGeo *spil = SwimGeo::Instance(*(const_cast<VldContext*>(vldc))); // Get edges of steel
01978     double z1 = spil->GetSteel(TrkStripData[Plane][0].ZPos,izdir)->GetZ();
01979     double z2 = spil->GetSteel(z1,izdir)->GetZ();
01980   
01981     double dzscatter = fabs(TrkStripData[NewPlane][0].ZPos-0.5*(z1+z2));
01982     double dzscatter2 = pow(dzscatter,2);
01983   
01984     UgliGeomHandle ugh(*vldc);
01985     BField bf(*vldc,-1,0);
01986     TVector3 uvz;
01987 
01988     uvz(0) = x_k_minus[0];
01989     uvz(1) = x_k_minus[1];
01990     uvz(2) = ugh.GetNearestSteelPlnHandle(TrkStripData[Plane][0].ZPos).GetZ0();
01991 
01992     TVector3 buvz = bf.GetBField(uvz, true);
01993     buvz[0] *= 0.3;
01994     buvz[1] *= 0.3;
01995     buvz[2] *= 0.3;
01996 
01997     double beta2inv = 1.0-0.0011*qp*qp;
01998 
01999     double eloss = 0.038 * double(DeltaPlane);
02000     double sigmaeloss2 = 0.25*eloss*dsdz*beta2inv*qp*qp;
02001     sigmaeloss2 *= sigmaeloss2;
02002 
02003 
02004     // Fill elements of noise matrix
02005     Q_k_minus[0][0]=dzscatter2*Sigma33Squared;
02006     Q_k_minus[0][1]=dzscatter2*Sigma34Squared;
02007     Q_k_minus[0][2]=dzscatter*Sigma33Squared;
02008     Q_k_minus[0][3]=dzscatter*Sigma34Squared;
02009     Q_k_minus[0][4]=dzscatter*sigmaeloss2*double(DeltaPlane)*buvz(1)*0.0254*dsdz*(1.+pow(x_k_minus[2],2));
02010 
02011     Q_k_minus[1][0]=dzscatter2*Sigma34Squared;
02012     Q_k_minus[1][1]=dzscatter2*Sigma44Squared;
02013     Q_k_minus[1][2]=dzscatter*Sigma34Squared;
02014     Q_k_minus[1][3]=dzscatter*Sigma44Squared;
02015     Q_k_minus[1][4]=dzscatter*sigmaeloss2*double(DeltaPlane)*buvz(0)*0.0254*dsdz*(1.+pow(x_k_minus[3],2));
02016 
02017     Q_k_minus[2][0]=dzscatter*Sigma33Squared;
02018     Q_k_minus[2][1]=dzscatter*Sigma34Squared;
02019     Q_k_minus[2][2]=Sigma33Squared;
02020     Q_k_minus[2][3]=Sigma34Squared;
02021     Q_k_minus[2][4]=sigmaeloss2*double(DeltaPlane)*buvz(1)*0.0254*dsdz*(1.+pow(x_k_minus[2],2));
02022 
02023     Q_k_minus[3][0]=dzscatter*Sigma34Squared;
02024     Q_k_minus[3][1]=dzscatter*Sigma44Squared;
02025     Q_k_minus[3][2]=Sigma34Squared;
02026     Q_k_minus[3][3]=Sigma44Squared;
02027     Q_k_minus[3][4]=sigmaeloss2*double(DeltaPlane)*buvz(0)*0.0254*dsdz*(1.+pow(x_k_minus[3],2));
02028 
02029     Q_k_minus[4][0]=dzscatter*sigmaeloss2*double(DeltaPlane)*buvz(1)*0.0254*dsdz*(1.+pow(x_k_minus[2],2));
02030     Q_k_minus[4][1]=dzscatter*sigmaeloss2*double(DeltaPlane)*buvz(0)*0.0254*dsdz*(1.+pow(x_k_minus[3],2));
02031     Q_k_minus[4][2]=sigmaeloss2*double(DeltaPlane)*buvz(1)*0.0254*dsdz*(1.+pow(x_k_minus[2],2));
02032     Q_k_minus[4][3]=sigmaeloss2*double(DeltaPlane)*buvz(0)*0.0254*dsdz*(1.+pow(x_k_minus[3],2));
02033     Q_k_minus[4][4]=sigmaeloss2;
02034   }
02035  
02036 
02037   // Display
02038   if(debug) {
02039     cout << "1e6 * Q_k_minus" << endl;
02040     for(int i=0; i<5; ++i) {
02041       for(int j=0; j<5; ++j) {
02042         cout << 1e6*Q_k_minus[i][j] << " ";
02043       }
02044       cout << endl;
02045     }
02046   }
02047 
02048 }

void AlgFitTrackCam::GoBackwards (  ) 

Definition at line 1388 of file AlgFitTrackCam.cxx.

References C_k, CalcKalmanGain(), EndCov, EndofRange, EndofRangePlane, ExtrapCovMatrix(), GetCombiPropagator(), GetNoiseMatrix(), H_k, Msg::kSynopsis, Msg::kVerbose, LastIteration, MaxPlane, MinPlane, MoveArrays(), MSG, NIter, PassTrack, SaveData, size, StoreFilteredData(), TrkStripData, UpdateCovMatrix(), UpdateStateVector(), VtxCov, x_k, and ZIncreasesWithTime.

Referenced by RunTheFitter().

01389 { 
01390   // Carry out the Kalman fit along the track in the direction of decreasing z
01391   MSG("AlgFitTrackCam",Msg::kSynopsis) << "GoBackwards, carry out fit in negative z direction" << endl;
01392   MSG("AlgFitTrackCam",Msg::kSynopsis) << "Fitting from StartPlane=" << MaxPlane << " to EndPlane=" << MinPlane << " end of range plane = " << EndofRangePlane << endl;
01393 
01394   Int_t StartPlane = MaxPlane; Int_t EndPlane=MinPlane;
01395   if(ZIncreasesWithTime){
01396     StartPlane = EndofRangePlane;
01397   }
01398   else EndofRangePlane = MinPlane;
01399  
01400   for (int i=StartPlane; i>=EndPlane; --i) {  
01401     if (TrkStripData[i].size()>0) {
01402       if (PassTrack) {
01403 
01404         //Find Prev Plane
01405         int NewPlane=-99;
01406         int k=(i-1);
01407         while (k>=MinPlane) {
01408           if (TrkStripData[k].size()>0) {NewPlane=k; break;}
01409           --k;
01410         }
01411         
01412         
01413         if (NewPlane!=-99) {
01414           // Define measurement function
01415           if (TrkStripData[NewPlane][0].PlaneView==3) {H_k[0]=0; H_k[1]=1; H_k[2]=0; H_k[3]=0; H_k[4]=0;}         else if (TrkStripData[NewPlane][0].PlaneView==2)  {H_k[0]=1; H_k[1]=0; H_k[2]=0; H_k[3]=0; H_k[4]=0;}
01416           
01417           MSG("AlgFitTrackCam",Msg::kVerbose) << "GoBackwards, Plane " << i << " ZPos " << TrkStripData[i][0].ZPos 
01418                                               << " PlaneView " << TrkStripData[i][0].PlaneView << endl;
01419           MSG("AlgFitTrackCam",Msg::kVerbose) << " NewPlane " << NewPlane << " NewZPos " << TrkStripData[NewPlane][0].ZPos 
01420                                               << " NewPlaneView " << TrkStripData[NewPlane][0].PlaneView << endl;
01421           
01422           bool CombiPropagatorOk=GetCombiPropagator(i,NewPlane,false);
01423           
01424           if(CombiPropagatorOk ) {
01425             GetNoiseMatrix(i,NewPlane);
01426             ExtrapCovMatrix();
01427             CalcKalmanGain(NewPlane);
01428             UpdateStateVector(i,NewPlane,false);
01429             UpdateCovMatrix();
01430             MoveArrays();
01431             if(SaveData) {StoreFilteredData(NewPlane);}
01432             MSG("AlgFitTrackCam",Msg::kVerbose) << "GoBackwards, Filtered q/p " << x_k[4] << endl;
01433           }
01434           else {
01435             PassTrack=false;
01436           }
01437         }
01438         
01439       }
01440       else {MSG("AlgFitTrackCam",Msg::kSynopsis) << "GoBackwards, Outside detector - track failed" << endl;}
01441     }
01442     //JAM end of range found
01443     if(EndofRange && LastIteration && !ZIncreasesWithTime){
01444       EndofRangePlane = i;
01445       break;
01446     }
01447 
01448   }
01449 
01450 
01451   // Store entries from covariance matrix for use in setting track properties
01452   if(NIter>1) {
01453     if(ZIncreasesWithTime==true) {
01454       VtxCov[0]=C_k[0][0]; VtxCov[1]=C_k[1][1];
01455       VtxCov[2]=C_k[2][2]; VtxCov[3]=C_k[3][3];
01456       VtxCov[4]=C_k[4][4];
01457     }
01458     else {
01459       EndCov[0]=C_k[0][0]; EndCov[1]=C_k[1][1];
01460       EndCov[2]=C_k[2][2]; EndCov[3]=C_k[3][3];
01461       EndCov[4]=C_k[4][4];
01462     }
01463   }
01464 
01465   MSG("AlgFitTrackCam",Msg::kSynopsis) << "GoBackwards, Finished" << endl;
01466 }

void AlgFitTrackCam::GoForwards (  ) 

Definition at line 1302 of file AlgFitTrackCam.cxx.

References C_k, CalcKalmanGain(), EndCov, EndofRange, EndofRangePlane, ExtrapCovMatrix(), GetCombiPropagator(), GetNoiseMatrix(), H_k, Msg::kSynopsis, Msg::kVerbose, LastIteration, MaxPlane, MinPlane, MoveArrays(), MSG, NIter, PassTrack, SaveData, size, StoreFilteredData(), TrkStripData, UpdateCovMatrix(), UpdateStateVector(), VtxCov, x_k, and ZIncreasesWithTime.

Referenced by RunTheFitter().

01303 { 
01304   // Carry out the Kalman fit along the track in the direction of increasing z
01305   MSG("AlgFitTrackCam",Msg::kSynopsis) << "GoForwards, carry out fit in positive z direction" << endl; 
01306   MSG("AlgFitTrackCam",Msg::kSynopsis) << "Fitting from StartPlane=" << MinPlane << " to EndPlane=" << MaxPlane << endl;
01307 
01308   // JAM in 2nd iteration, stop when end of range is reached.
01309  
01310   Int_t StartPlane = MinPlane; Int_t EndPlane=MaxPlane;
01311   if(!ZIncreasesWithTime){
01312     StartPlane = EndofRangePlane;
01313   }
01314   else EndofRangePlane = MaxPlane;
01315 
01316   for (int i=StartPlane; i<=EndPlane; ++i) {
01317     EndofRange = false;
01318     if (TrkStripData[i].size()>0) {
01319       if (PassTrack) {
01320         // Find Next Plane
01321         int NewPlane=-99;
01322         int k=(i+1);
01323         while (k<=MaxPlane) {
01324           if (TrkStripData[k].size()>0) {NewPlane=k; break;}
01325           ++k;
01326         }
01327         if (NewPlane!=-99) {
01328           // Define measurement function
01329           if (TrkStripData[NewPlane][0].PlaneView==3) {H_k[0]=0; H_k[1]=1; H_k[2]=0; H_k[3]=0; H_k[4]=0;}
01330           else if (TrkStripData[NewPlane][0].PlaneView==2)  {H_k[0]=1; H_k[1]=0; H_k[2]=0; H_k[3]=0; H_k[4]=0;}
01331           
01332           MSG("AlgFitTrackCam",Msg::kVerbose) << "GoForwards, Plane " << i << " ZPos " << TrkStripData[i][0].ZPos 
01333                                               << " PlaneView " << TrkStripData[i][0].PlaneView << endl;
01334           MSG("AlgFitTrackCam",Msg::kVerbose) << " NewPlane " << NewPlane << " NewZPos " << TrkStripData[NewPlane][0].ZPos 
01335                                               << " NewPlaneView " << TrkStripData[NewPlane][0].PlaneView << endl;
01336           
01337           bool CombiPropagatorOk=GetCombiPropagator(i,NewPlane,true);
01338           
01339           if(CombiPropagatorOk ) {
01340             GetNoiseMatrix(i,NewPlane);
01341             ExtrapCovMatrix();
01342             CalcKalmanGain(NewPlane);
01343             UpdateStateVector(i,NewPlane,true);
01344             UpdateCovMatrix();
01345             MoveArrays();
01346             if(SaveData) {StoreFilteredData(NewPlane);}
01347             MSG("AlgFitTrackCam",Msg::kVerbose) << "GoForwards, Filtered q/p " << x_k[4] << endl;
01348           }
01349           else 
01350             {
01351               PassTrack=false; 
01352             }
01353         }
01354         
01355       }
01356       else {MSG("AlgFitTrackCam",Msg::kSynopsis) << "GoForwards, Outside of detector - track failed" << endl;}
01357     }
01358     //JAM end of range found
01359     if(EndofRange && LastIteration && ZIncreasesWithTime){
01360       EndofRangePlane=i;
01361       MSG("AlgFitTrackCam",Msg::kSynopsis) << "End of range on plane " << EndofRangePlane << endl;
01362       break;
01363     }
01364   }
01365 
01366 
01367   // Store entries from covariance matrix for use in setting track properties
01368   if(NIter>1) {
01369     if(ZIncreasesWithTime==true) {
01370       EndCov[0]=C_k[0][0]; EndCov[1]=C_k[1][1];
01371       EndCov[2]=C_k[2][2]; EndCov[3]=C_k[3][3];
01372       EndCov[4]=C_k[4][4];
01373     }
01374     else {
01375       VtxCov[0]=C_k[0][0]; VtxCov[1]=C_k[1][1];
01376       VtxCov[2]=C_k[2][2]; VtxCov[3]=C_k[3][3];
01377       VtxCov[4]=C_k[4][4];
01378     }
01379   }
01380 
01381   MSG("AlgFitTrackCam",Msg::kSynopsis) << "GoForwards, Finished" << endl;
01382 }

void AlgFitTrackCam::InitialFramework ( const CandSliceHandle slice,
CandContext cx 
)

Definition at line 255 of file AlgFitTrackCam.cxx.

References StripStruct::csh, GenerateNDSpectStrips(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), InitTrkStripData, IsStripInAnotherTrack(), Detector::kNear, MaxPlane, MeanTrackTime, MinPlane, NDStripBegTime(), NumFinderStrips, SlcStripData, track, and vldc.

Referenced by RunAlg().

00256 {
00257   // Store CandStripHandles and make the strips accessible by plane number
00258   Detector::Detector_t detector = vldc->GetDetector();
00259 
00260   TIter SlcStripItr = slice->GetDaughterIterator();
00261   TIter TrkStripItr = track->GetDaughterIterator();
00262   
00263   // Store all strips in slice
00264   int SlicePlane; 
00265 
00266   while(CandStripHandle* SlcStrip = dynamic_cast<CandStripHandle*>(SlcStripItr()))
00267     {
00268       // Leigh: Check if this strip has been used in another track in this slice.
00269       if(IsStripInAnotherTrack(SlcStrip)) continue;
00270 
00271       SlicePlane=SlcStrip->GetPlane();
00272       if(detector==Detector::kNear && SlicePlane>=121) {continue;}
00273 
00274       StripStruct temp;
00275       temp.csh=SlcStrip;
00276      
00277       SlcStripData[SlicePlane].push_back(temp);
00278     }
00279   SlcStripItr.Reset();
00280 
00281 
00282   // Store all track strips found, except those in the Near Spectrometer
00283   int TrackPlane; 
00284   MeanTrackTime=0;
00285 
00286         std::vector<double> trackTimes;
00287 
00288   while(CandStripHandle* TrkStrip = dynamic_cast<CandStripHandle*>(TrkStripItr()))
00289     {
00290     
00291       TrackPlane=TrkStrip->GetPlane();
00292       if(detector==Detector::kNear && TrackPlane>=121) {continue;}
00293 
00294       StripStruct temp;
00295       temp.csh=TrkStrip;
00296 
00297       InitTrkStripData[TrackPlane].push_back(temp);
00298       NumFinderStrips++;
00299       
00300       // Identify ends of initial track
00301       if (TrackPlane>MaxPlane) {MaxPlane=TrackPlane;}
00302       if (TrackPlane<MinPlane) {MinPlane=TrackPlane;}
00303 
00304       if(detector==Detector::kNear) {
00305 //                              MeanTrackTime+=NDStripBegTime(TrkStrip);
00306                                 trackTimes.push_back(NDStripBegTime(TrkStrip));
00307                         }
00308     }
00309   TrkStripItr.Reset();
00310 
00311   // For DeMuxing ND Spectrometer
00312   if(detector==Detector::kNear) {
00313 //    if(NumFinderStrips!=0) {MeanTrackTime/=double(NumFinderStrips);}
00314                 // Sort the track times and find the middle element
00315                 std::sort(trackTimes.begin(), trackTimes.end());
00316                 int middleElement = NumFinderStrips / 2;
00317     if(NumFinderStrips!=0) {MeanTrackTime = trackTimes[middleElement];}
00318     GenerateNDSpectStrips(slice,cx);
00319   }
00320 
00321 } 

bool AlgFitTrackCam::IsCoilReverse ( VldContext vldc  ) 

Definition at line 3947 of file AlgFitTrackCam.cxx.

References CoilTools::IsReverse().

Referenced by GetCombiPropagator().

03947                                                   {
03948   
03949   // Use DCS to obtain direction of coil current.
03950   return CoilTools::IsReverse(*vldc);
03951 
03952 }

int AlgFitTrackCam::IsInDetector ( const double  upos,
const double  vpos,
const int  plane 
)

Definition at line 3923 of file AlgFitTrackCam.cxx.

References PlaneOutline::DistanceToEdge(), fPL, PlexPlaneId::GetPlaneCoverage(), PlexPlaneId::GetPlaneView(), PlaneOutline::IsInside(), and Detector::kNear.

Referenced by SpectrometerSwim().

03923                                                                                      {
03924   // Need to convert to x and y.
03925   PlexPlaneId plexPlane(Detector::kNear, plane, 0);
03926   double tempx = 0.707 * (upos - vpos);
03927   double tempy = 0.707 * (upos + vpos);
03928 
03929         // LEIGH HACK: Want to stop the tracking if we get close to the coil hole.
03930         // This is going to return 2 as it should immmediately stop the tracking.
03931         double tempr = sqrt(upos*upos + vpos*vpos);
03932 
03933         if(tempr < 0.6) return 2;
03934 
03935   // Is it outside the detector? Inside the coil hole also counts as outside.
03936   if(!fPL.IsInside(tempx, tempy, plexPlane.GetPlaneView(), plexPlane.GetPlaneCoverage())){
03937     // This position is outside of the detector, how far?
03938     float dist, distx, disty;
03939     fPL.DistanceToEdge(tempx,tempy,plexPlane.GetPlaneView(), plexPlane.GetPlaneCoverage(),dist,distx,disty);
03940                 // Allow a 1cm tolerance between MC geometry and as-built.
03941     if(dist > 0.01) return 0;
03942     else return 1;
03943   }
03944   else return 1;
03945 }

bool AlgFitTrackCam::IsStripInAnotherTrack ( CandStripHandle strip  ) 

Definition at line 3911 of file AlgFitTrackCam.cxx.

References CandStripHandle::GetPlane(), CandHandle::GetUidInt(), and usedSlcStrips.

Referenced by InitialFramework().

03911                                                                 {
03912   // Since we can have multiple tracks per slice, check this strip is not already
03913   // assigned to another track in the slice.
03914   std::vector<CandStripHandle*> hitVec = AlgFitTrackCam::usedSlcStrips[strip->GetPlane()];
03915   // Iterate over the vector
03916   for(unsigned int i = 0; i < hitVec.size(); ++i){
03917     if(strip->GetUidInt() == hitVec[i]->GetUidInt()) {return true;}
03918   }
03919   // If we get here then it is fine to use this strip
03920   return false;
03921 }

double AlgFitTrackCam::MajorityCurvature (  ) 

Definition at line 2281 of file AlgFitTrackCam.cxx.

References FilteredData, MaxPlane, MinPlane, and size.

Referenced by UpdateStateVector().

02282 { 
02283   //  MSG("AlgFitTrackCam",Msg::kDebug) << "MajorityCurvature" << endl;
02284 
02285   Double_t majority_curvature = 0.0;
02286   Double_t majority_counter = 0.0;
02287   Double_t total_counter = 0.0;
02288 
02289   for(int i=MinPlane; i<=MaxPlane; ++i) {
02290     if( FilteredData[i].size()>0 ) {
02291       if( FilteredData[i][0].x_k4>0.0 ) majority_counter += 1.0;
02292       if( FilteredData[i][0].x_k4<0.0 ) majority_counter -= 1.0;
02293       total_counter += 1.0;
02294     }
02295   }
02296 
02297   if( total_counter>0.0 ){
02298     majority_curvature = majority_counter/total_counter;
02299   }
02300   
02301   return majority_curvature;
02302 }

void AlgFitTrackCam::MoveArrays (  ) 

Definition at line 2350 of file AlgFitTrackCam.cxx.

References C_k, C_k_minus, x_k, and x_k_minus.

Referenced by GoBackwards(), GoForwards(), ShowerSwim(), and SpectrometerSwim().

02351 {
02352   // Move k to k-1 ready to consider next hit
02353   //  MSG("AlgFitTrackCam",Msg::kDebug) << "MoveArrays" << endl;
02354 
02355   for (int i=0; i<5; ++i) {
02356     for (int j=0; j<5; ++j) { 
02357       C_k_minus[i][j]=C_k[i][j];
02358     }
02359   }
02360   
02361   for (int l=0; l<5; ++l) {
02362     x_k_minus[l]=x_k[l];
02363   }
02364 }

bool AlgFitTrackCam::NDPlaneIsActive ( int  plane,
float  u,
float  v,
float  projErr 
)

Definition at line 3655 of file AlgFitTrackCam.cxx.

References PlaneOutline::DistanceToEdge(), fPL, PlexPlaneId::GetPlaneCoverage(), PlexPlaneId::GetPlaneView(), PlaneOutline::IsInside(), PlexPlaneId::IsValid(), Detector::kNear, and PlaneCoverage::kNoActive.

Referenced by SpectrometerSwim().

03656 { 
03657   // Method to determine whether this u/v position is active
03658 
03659   PlexPlaneId plexPlane(Detector::kNear,plane, 0); 
03660   if(!plexPlane.IsValid()) {return false;}
03661   if(plexPlane.GetPlaneCoverage()==PlaneCoverage::kNoActive) {return false;}
03662   if(projErr<0.3)projErr=0.3;
03663   float x = 0.707*(u-v);
03664   float y = 0.707*(u+v);
03665   float dist,xedge,yedge;
03666   fPL.DistanceToEdge(x, y,
03667                      plexPlane.GetPlaneView(),
03668                      plexPlane.GetPlaneCoverage(),
03669                      dist, xedge, yedge);
03670   bool isInside = fPL.IsInside(x, y,
03671                                plexPlane.GetPlaneView(),
03672                                plexPlane.GetPlaneCoverage());
03673   isInside &= dist>projErr;
03674 
03675   return isInside;
03676 }

double AlgFitTrackCam::NDStripBegTime ( CandStripHandle Strip,
double  U = 0,
double  V = 0 
)

Definition at line 3866 of file AlgFitTrackCam.cxx.

References BegTime, UgliStripHandle::ClearFiber(), digit(), CandHandle::GetDaughterIterator(), PlexSEIdAltL::GetEnd(), UgliStripHandle::GetHalfLength(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), CandDigitHandle::GetPlexSEIdAltL(), CandStripHandle::GetStripEndId(), UgliGeomHandle::GetStripHandle(), CandDigitHandle::GetSubtractedTime(), CandTrackHandle::GetU(), CandTrackHandle::GetV(), StripEnd::kPositive, CalTimeType::kT0, StripEnd::kUnknown, track, vldc, and UgliStripHandle::WlsPigtail().

Referenced by InitialFramework(), and SpectrometerSwim().

03867 {
03868   double BegTime=999; double Time=0;
03869   double Index=1.77; double DistFromCentre=0.; 
03870   double FibreDist=0.; double halfLength=0.;
03871 
03872   // Get from track. Otherwise, will have guessed using Swimmer
03873   if(U==0) {U=track->GetU(Strip->GetPlane());}
03874   if(V==0) {V=track->GetV(Strip->GetPlane());}
03875 
03876   StripEnd::StripEnd_t StpEnd = StripEnd::kUnknown;
03877   CandDigitHandle* digit;
03878   
03879   UgliGeomHandle ugh = UgliGeomHandle(*vldc);
03880   UgliStripHandle Striphandle;
03881 
03882   CandDigitHandleItr digitItr(Strip->GetDaughterIterator());
03883 
03884 
03885   // Loop over digits on Strip. 
03887   while( (digit = digitItr()) ) { 
03888     StpEnd=digit->GetPlexSEIdAltL().GetEnd();
03889     
03890     if(StpEnd==StripEnd::kPositive) {
03891       FibreDist = 0.; DistFromCentre = 0.; Time=0.;
03892       UgliStripHandle StripHandle = ugh.GetStripHandle(Strip->GetStripEndId());
03893       halfLength = StripHandle.GetHalfLength(); 
03894       
03895       if(Strip->GetPlaneView()==2) {DistFromCentre = V;}
03896       if(Strip->GetPlaneView()==3) {DistFromCentre = -U;}
03897               
03898       FibreDist = (halfLength + DistFromCentre + StripHandle.ClearFiber(StpEnd) 
03899                    + StripHandle.WlsPigtail(StpEnd));
03900       
03901       Time = digit->GetSubtractedTime(CalTimeType::kT0) - (Index/3.e8)*FibreDist;
03902 
03903       if(Time<BegTime) {BegTime=Time;}
03904     }
03905   }
03906   
03907   return BegTime;
03908 }

void AlgFitTrackCam::RemoveTrkHitsInShw (  ) 

Definition at line 649 of file AlgFitTrackCam.cxx.

References Msg::kDebug, MaxPlane, MinPlane, MSG, ShowerEntryPlane, size, SwimThroughShower, TrkStripData, and ZIncreasesWithTime.

Referenced by RunTheFitter().

00650 {
00651   // If the 'clean' section of track is large enough, remove the track finding
00652   // data for planes before the ShowerEntryPlane
00653   MSG("AlgFitTrackCam",Msg::kDebug) << "RemoveTrkHitsInShw, Discard track finding data in shower" << endl;
00654 
00655   int NumTrackStripsLeft=0;
00656   
00657   if(ZIncreasesWithTime==true) {
00658     for(int i=ShowerEntryPlane; i<=MaxPlane; ++i) {
00659       if(TrkStripData[i].size()>0) {NumTrackStripsLeft++;}
00660     }
00661   }
00662   else if(ZIncreasesWithTime==false) {
00663     for(int i=MinPlane; i<=ShowerEntryPlane; ++i) {
00664       if(TrkStripData[i].size()>0) {NumTrackStripsLeft++;}
00665     }
00666   }
00667   
00668   // Carry out removal if there will be 6 or more strips left afterwards
00669   if(NumTrackStripsLeft>5) { 
00670     if(ZIncreasesWithTime==true) {
00671       for(int i=MinPlane; i<=ShowerEntryPlane; ++i) {TrkStripData[i].clear();}
00672     }   
00673     else if(ZIncreasesWithTime==false) {
00674       for(int i=ShowerEntryPlane; i<=MaxPlane; ++i) {TrkStripData[i].clear();    }
00675     }
00676   } 
00677   // Otherwise note that we should not run the ShowerSwim method
00678   else {
00679     MSG("AlgFitTrackCam",Msg::kDebug) << "RemoveTrkHitsInShw, not enough hits after removal. Must use all finder data." << endl;
00680     SwimThroughShower=false;
00681   }
00682   
00683   
00684   // Find the new max and min planes
00685   MaxPlane=-20; MinPlane=500;
00686   for (int i=0; i<490; ++i) {   
00687     if(TrkStripData[i].size()>0) {
00688       if(i>MaxPlane) {MaxPlane=i;}
00689       if(i<MinPlane) {MinPlane=i;}
00690     }
00691   }
00692 
00693 }

void AlgFitTrackCam::ResetCovarianceMatrix (  ) 

Definition at line 1856 of file AlgFitTrackCam.cxx.

References DeltaPlane, DeltaZ, and GetInitialCovarianceMatrix().

Referenced by RunTheFitter().

01857 {
01858   // Simple method reset variables/arrays to allow propagation again
01859 
01860   DeltaPlane=0; DeltaZ=0; 
01861   GetInitialCovarianceMatrix(false);
01862 }

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

Implements AlgBase.

Definition at line 94 of file AlgFitTrackCam.cxx.

References bave, C_k, C_k_intermediate, C_k_minus, CleanNDLists(), debug, DeltaPlane, DeltaZ, EndCov, F_k, F_k_minus, FilteredData, Registry::Get(), CandRecord::GetCandHeader(), CandContext::GetCandRecord(), CandContext::GetDataIn(), VldContext::GetDetector(), GetFitData(), CandHeader::GetRun(), CandHeader::GetSnarl(), RecMinos::GetVldContext(), H_k, Identity, InitialFramework(), InitTrkStripData, K_k, Registry::KeyExists(), Msg::kInfo, Detector::kNear, Msg::kSynopsis, len, MaxPlane, MinPlane, MSG, nbfield, NIter, NumFinderStrips, PassTrack, Q_k, Q_k_minus, RunTheFitter(), SaveData, ShowerEntryPlane, SlcStripData, SwimThroughShower, TotalNSwimFail, track, TrkStripData, usedSlcStrips, UseGeoSwimmer, vldc, VtxCov, x_k, x_k4_biased, x_k_minus, and ZIncreasesWithTime.

00096 {
00097 
00098   // cout << "DOGWOOD3  ----------- " << endl;
00099 
00100   // get alg parameters
00101   if(!ac.Get("UseGeoSwimmer",UseGeoSwimmer)) cout << "Couldn't Get UseGeoSwimmer\n";
00102 
00103  // Standard set-up
00104   assert(ch.InheritsFrom("CandFitTrackCamHandle"));
00105   CandFitTrackCamHandle &cth = dynamic_cast<CandFitTrackCamHandle &>(ch);
00106   
00107   nbfield=0;
00108   bave=0;
00109   TIter FitTrackStripItr = cth.GetDaughterIterator();
00110   while(CandStripHandle* FitTrackStrip = dynamic_cast<CandStripHandle*>(FitTrackStripItr())) cth.RemoveDaughter(FitTrackStrip);
00111   
00112   assert(cx.GetDataIn());
00113   
00114   // Get the track from the track finder
00115   assert(cx.GetDataIn()->InheritsFrom("CandTrackHandle"));
00116   track = dynamic_cast<const CandTrackHandle*>(cx.GetDataIn());
00117   if( !track || track->GetNDaughters()<1 )  { return; }
00118   cth.SetFinderTrack((CandTrackHandle*)(track));
00119   
00120   const CandSliceHandle* slice = dynamic_cast<const CandSliceHandle*>(track->GetCandSlice());
00121   assert(slice);
00122   cth.SetCandSlice(slice);
00123   
00124   
00126   // Track fitter //
00128   // Switch for screen output
00129   debug=false;
00130   
00131   // Initialisations
00133   if( track->GetEndZ() > track->GetVtxZ() )  {ZIncreasesWithTime=true;}
00134   else {ZIncreasesWithTime=false;}
00135   
00136   SaveData=false;
00137   SwimThroughShower=false;
00138   PassTrack=true;
00139 
00140   MaxPlane=-20;
00141   MinPlane=500;
00142   
00143   DeltaZ=-99;
00144   DeltaPlane=-99;
00145   ShowerEntryPlane=-99;
00146   
00147   NIter=0; TotalNSwimFail=0; NumFinderStrips=0;
00148   
00149   for (unsigned int i=0; i<5; ++i) {
00150     x_k[i]=0; x_k_minus[i]=0; H_k[i]=0; K_k[i]=0;
00151     VtxCov[i]=-999; EndCov[i]=-999;
00152     for (unsigned int j=0; j<5; ++j) {
00153       C_k[i][j]=0; C_k_minus[i][j]=0; C_k_intermediate[i][j]=0;
00154       F_k[i][j]=0; F_k_minus[i][j]=0;
00155       Q_k[i][j]=0; Q_k_minus[i][j]=0;
00156       Identity[i][j]=0;
00157     }
00158   }
00159   
00160   Identity[0][0]=1; Identity[1][1]=1; Identity[2][2]=1; Identity[3][3]=1; Identity[4][4]=1; 
00162   
00163   
00164   // Set initial parameters
00166   x_k_minus[0]=track->GetVtxU();
00167   x_k_minus[1]=track->GetVtxV();
00168   if(track->GetVtxDirCosZ()!=0.) {
00169     x_k_minus[2]=track->GetVtxDirCosU()/track->GetVtxDirCosZ();
00170     x_k_minus[3]=track->GetVtxDirCosV()/track->GetVtxDirCosZ();
00171   }
00172   x_k_minus[4]=0.;
00173   x_k4_biased=0;
00175   
00176   
00177   // Get header information
00179   CandRecord* candrec = (CandRecord*)(cx.GetCandRecord());
00180   assert(candrec);
00181   
00182   vldc = (VldContext*)(candrec->GetVldContext());
00183   assert(vldc);
00184   
00185   CandHeader* candhead = (CandHeader*)(candrec->GetCandHeader());
00186   assert(candhead);
00187   
00188   MSG("AlgFitTrackCam",Msg::kSynopsis) << "RunAlg, New track, Run: " << candhead->GetRun() 
00189                                     << " Snarl: " << candhead->GetSnarl() << endl;
00191   
00192   
00193   // Run the high level methods
00195   InitialFramework(slice,cx);
00196 
00197 
00198 
00200   // set min and max planes
00201   // the half-track fitting study modifies the defaults here  
00202   const char* htk="UseHalfTrack";
00203   if(ac.KeyExists(htk)){
00204     int ihalf=0; 
00205     ac.Get(htk,ihalf);    
00206     if(ihalf!=0){
00207       int beg_plane=track->GetBegPlane();
00208       int end_plane=track->GetEndPlane();
00209       int len=end_plane-beg_plane +1;
00210       if (len>20){// requre 20 planes or more
00211         int half_way=beg_plane+int(len/2.0); 
00212         MSG("AlgFitTrackCam",Msg::kInfo)<<"Track [beg,mid,end] = ["<<beg_plane<<", "<<half_way<<", "<<end_plane<<"]"<<endl;
00213         if(ihalf==1){ // use first half of track
00214           MaxPlane=half_way;
00215           MSG("AlgFitTrackCam",Msg::kInfo)<<"Fitting first half of track between planes "<<MinPlane<<" and "<<MaxPlane<<endl;
00216         }
00217         else if(ihalf==2){ // second half
00218           MinPlane=half_way;
00219           MSG("AlgFitTrackCam",Msg::kInfo)<<"Fitting second half of track between planes "<<MinPlane<<" and "<<MaxPlane<<endl;
00220         }
00221       }
00222     }
00223   }
00224   // end of min,max plane setting
00226 
00227   
00228   GetFitData(MinPlane,MaxPlane);
00229   RunTheFitter(cth);
00231   
00232   // Leigh: Loop over the strips of the created track to ensure we don't add them to the next track.
00233   FitTrackStripItr = cth.GetDaughterIterator();
00234   while(CandStripHandle* FitTrackStrip = dynamic_cast<CandStripHandle*>(FitTrackStripItr())){
00235     AlgFitTrackCam::usedSlcStrips[FitTrackStrip->GetPlane()].push_back(FitTrackStrip);
00236   }
00237 
00238    // Clear up
00240   if(vldc->GetDetector()==Detector::kNear) {CleanNDLists(cth,cx);}
00241   
00242   for (unsigned int i=0; i<490; ++i) {
00243     InitTrkStripData[i].clear();
00244     SlcStripData[i].clear();
00245     TrkStripData[i].clear();
00246     FilteredData[i].clear();
00247   }
00249 }

void AlgFitTrackCam::RunTheFitter ( CandFitTrackCamHandle cth  ) 

Definition at line 423 of file AlgFitTrackCam.cxx.

References CandHandle::AddDaughterLink(), EndofRange, EndofRangePlane, eqp_biased, FillGapsInTrack(), FilteredData, FindTheStrips(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), GetFitData(), GetInitialCovarianceMatrix(), GoBackwards(), GoForwards(), Detector::kNear, Msg::kSynopsis, LastIteration, MaxPlane, MinPlane, MSG, NIter, PassTrack, CandHandle::RemoveDaughter(), RemoveTrkHitsInShw(), ResetCovarianceMatrix(), SaveData, CandFitTrackHandle::SetEMCharge(), CandFitTrackHandle::SetMomentumCurve(), CandFitTrackHandle::SetPass(), SetPropertiesFromFinderTrack(), SetTrackProperties(), CandFitTrackHandle::SetVtxQPError(), ShowerStrips(), ShowerSwim(), SpectrometerSwim(), StoreFilteredData(), SwimThroughShower, track, TrkStripData, vldc, VtxCov, x_k, x_k4_biased, x_k_minus, and ZIncreasesWithTime.

Referenced by RunAlg().

00424 {
00425   MSG("AlgFitTrackCam",Msg::kSynopsis) << "RunTheFitter, Call methods in the appropriate order" << endl;
00426   GetInitialCovarianceMatrix(true);
00427   
00428   // Control the iterations backwards and forwards
00430   Detector::Detector_t detector = vldc->GetDetector();
00431 
00432   // Control iterations over a track for which ZIncreasesWithTime
00433   if(ZIncreasesWithTime==true)
00434     {
00435       // First iteration
00436       NIter++;
00437       LastIteration=false;
00438 
00439       // Vtx to End
00440       StoreFilteredData(MinPlane);   
00441       SaveData=true; GoForwards(); 
00442       ResetCovarianceMatrix();
00443       
00444       // End back to vtx, swimming through any large vtx shower
00445       ShowerStrips();  // Try to identify vtx showers, now we have an idea of gradient
00446       if(SwimThroughShower==true) {RemoveTrkHitsInShw();}
00447       
00448       for (unsigned int i=0; i<490; ++i) {FilteredData[i].clear();}
00449       StoreFilteredData(MaxPlane); GoBackwards();
00450       
00451       if(SwimThroughShower==true) {ShowerSwim();}
00452       ResetCovarianceMatrix();
00453 
00454       bool StripsFound = FindTheStrips(cth,false);
00455       EndofRangePlane=MaxPlane;
00456       // Second iteration
00457       if(StripsFound==true) { // Guard against finding no strips
00458         for(int nint=0;nint<2;nint++){  
00459           if(nint==0) MSG("AlgFitTrackCam",Msg::kSynopsis) << "Bias Fit" << endl;
00460           if(nint==1) MSG("AlgFitTrackCam",Msg::kSynopsis) << "Un-Biased Fit" << endl;
00461           NIter++;
00462           if(nint==1)LastIteration = true;
00463           
00464           // Vtx to End again, identifying any strips in ND spectrometer
00465           for (unsigned int i=0; i<490; ++i) {TrkStripData[i].clear();} 
00466           GetFitData(MinPlane,MaxPlane);
00467 
00468           // for (unsigned int i=0; i<490; ++i) { FilteredData[i].clear(); }
00469           StoreFilteredData(MinPlane);          
00470           SaveData=true; GoForwards();
00471           
00472           if(detector==Detector::kNear && nint==0) {SpectrometerSwim(cth);}
00473           ResetCovarianceMatrix();
00474           
00475           // End back to vtx again
00476           //for (unsigned int i=0; i<490; ++i) { FilteredData[i].clear(); }
00477 
00478           // If EndofRange has been set, then we need to set the current state of the
00479           // Kalman filter to EndofRangePlane, since we have moved on to the next point
00480           if(EndofRange && LastIteration){
00481                   x_k_minus[0] = FilteredData[EndofRangePlane][0].x_k0;
00482                   x_k_minus[1] = FilteredData[EndofRangePlane][0].x_k1;
00483                   x_k_minus[2] = FilteredData[EndofRangePlane][0].x_k2;
00484                   x_k_minus[3] = FilteredData[EndofRangePlane][0].x_k3;
00485                   x_k_minus[4] = FilteredData[EndofRangePlane][0].x_k4;
00486           }
00487 
00488           StoreFilteredData(EndofRangePlane);    
00489           SaveData=true; GoBackwards();
00490           
00491           if(nint==0) {
00492             x_k4_biased= x_k[4];
00493             eqp_biased = pow(VtxCov[4],0.5);
00494           }
00495           ResetCovarianceMatrix();
00496         }
00497       }
00498       else {PassTrack=false;}
00499     }
00500   
00501   
00502   // Control iterations over a track for which ZDecreasesWithTime
00503   else
00504     {
00505       // First iteration
00506       NIter++;
00507       LastIteration=false;
00508 
00509       // Vtx to End
00510       StoreFilteredData(MaxPlane); 
00511       SaveData=true; GoBackwards(); 
00512       ResetCovarianceMatrix();
00513  
00514       // End back to vtx, swimming through any large vtx shower and
00515       // identifying any strips in ND spectrometer
00516       ShowerStrips();  // Try to identify vtx showers, now we have an idea of gradient
00517       if(SwimThroughShower==true) {RemoveTrkHitsInShw();}
00518 
00519       for (unsigned int i=0; i<490; ++i) {FilteredData[i].clear();}
00520       StoreFilteredData(MinPlane); GoForwards();
00521 
00522       if(SwimThroughShower==true) {ShowerSwim();}
00523 
00524       if(detector==Detector::kNear) {SpectrometerSwim(cth);}
00525       ResetCovarianceMatrix();
00526 
00527       bool StripsFound = FindTheStrips(cth,false);
00528       EndofRangePlane=MinPlane;
00529       // Second iteration
00530       if(StripsFound==true) { // Guard against finding no strips
00531         for(int nint=0;nint<2;nint++){  
00532           if(nint==1)LastIteration = true;                
00533           NIter++;
00534 
00535           // Vtx to End again
00536           for (unsigned int i=0; i<490; ++i) {TrkStripData[i].clear();} 
00537           GetFitData(MinPlane,MaxPlane);
00538 
00539           //for (unsigned int i=0; i<490; ++i) {FilteredData[i].clear();}
00540           StoreFilteredData(MaxPlane);    
00541           SaveData=true; GoBackwards(); 
00542           ResetCovarianceMatrix();
00543           
00544           // If EndofRange has been set, then we need to set the current state of the
00545           // Kalman filter to EndofRangePlane, since we have moved on to the next point
00546           if(EndofRange && LastIteration){
00547                   x_k_minus[0] = FilteredData[EndofRangePlane][0].x_k0;
00548                   x_k_minus[1] = FilteredData[EndofRangePlane][0].x_k1;
00549                   x_k_minus[2] = FilteredData[EndofRangePlane][0].x_k2;
00550                   x_k_minus[3] = FilteredData[EndofRangePlane][0].x_k3;
00551                   x_k_minus[4] = FilteredData[EndofRangePlane][0].x_k4;
00552           }
00553 
00554           // End to Vtx again
00555           //for (unsigned int i=0; i<490; ++i) {FilteredData[i].clear();}
00556           StoreFilteredData(EndofRangePlane);
00557           SaveData=true; GoForwards(); 
00558 
00559           ResetCovarianceMatrix();
00560           if(nint==0) {
00561             x_k4_biased= x_k[4];
00562             eqp_biased = pow(VtxCov[4],0.5);
00563           }
00564         }
00565       }
00566       else {PassTrack=false;}
00567 
00568     }
00570   
00571   
00572 
00573   // Organise the output
00575 
00576   // If the fit was successful
00577   if(x_k[4]!=0. && PassTrack==true) {
00578         
00579     //JAM modify tweak following range bias removal
00580     // Tweak q/p to remove offset
00581     //    x_k[4]*=1.01+(0.1*fabs(x_k[4]));
00582 
00583     x_k4_biased *=1.01+(0.1*fabs(x_k[4]));
00584     // x_k[4] *=1.013; // (moved to SetTrackProperties)
00585 
00586     // Find final strips and add them to the fitted track
00587     FillGapsInTrack();
00588     bool FinalStripsFound = FindTheStrips(cth,true);
00589     
00590     // If final strips found, set the fitted track properties
00591     if(FinalStripsFound==true) {
00592 
00593       // Check there is >1 strip in each view. If not, then fail track.
00594       int NumInUView=0; int NumInVView=0;
00595       
00596       TIter FitTrackStripItr = cth.GetDaughterIterator();
00597       while(CandStripHandle* FitTrackStrip = dynamic_cast<CandStripHandle*>(FitTrackStripItr())){
00598         if(FitTrackStrip->GetPlaneView()==2) {NumInUView++;}
00599         else if(FitTrackStrip->GetPlaneView()==3) {NumInVView++;}
00600         
00601         if(NumInUView>1 && NumInVView>1) {break;}
00602       }
00603       
00604       if(NumInUView>1 && NumInVView>1) {
00605         cth.SetPass(1);
00606         SetTrackProperties(cth);
00607       }
00608       else {
00609         PassTrack=false;
00610       }
00611       
00612     }
00613     // Otherwise fail the track at this final stage
00614     else {
00615       PassTrack=false;
00616     }
00617   }
00618   
00619   
00620   // If the fit has failed (e.g. q/p is zero and/or u, v are nonsense)
00621   if(x_k[4]==0. || PassTrack==false) {
00622  
00623     // Remove any existing strips in the failed fitted track
00624     vector<CandStripHandle*> Daughters;
00625 
00626     TIter FitTrackStripItr = cth.GetDaughterIterator();
00627     while(CandStripHandle* FitTrackStrip = dynamic_cast<CandStripHandle*>(FitTrackStripItr()))
00628       {Daughters.push_back(FitTrackStrip);}
00629 
00630     for(unsigned int i=0; i<Daughters.size(); ++i) {cth.RemoveDaughter(Daughters[i]);}
00631     Daughters.clear();
00632 
00633 
00634     // Put strips from track finder in failed fitted track
00635     TIter TrkStripItr = track->GetDaughterIterator();
00636     while(CandStripHandle* TrkStrip = dynamic_cast<CandStripHandle*>(TrkStripItr()))
00637       {cth.AddDaughterLink(*TrkStrip);}
00638     
00639     // Set position/direction properties using the finder track
00640     cth.SetPass(0); 
00641     cth.SetMomentumCurve(0.); cth.SetEMCharge(0); cth.SetVtxQPError(-999.);
00642     SetPropertiesFromFinderTrack(cth);
00643   }
00645 }  

void AlgFitTrackCam::SetPropertiesFromFinderTrack ( CandFitTrackCamHandle cth  ) 

Definition at line 3294 of file AlgFitTrackCam.cxx.

References AlgTrack::CalculateTrace(), AlgReco::Calibrate(), 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(), CandTrackHandle::GetNTimeFitDigit(), CandTrackHandle::GetRange(), CandTrackHandle::GetTimeBackwardFitNDOF(), CandTrackHandle::GetTimeBackwardFitRMS(), CandTrackHandle::GetTimeFitChi2(), CandTrackHandle::GetTimeForwardFitNDOF(), CandTrackHandle::GetTimeForwardFitRMS(), CandRecoHandle::GetTimeOffset(), CandRecoHandle::GetTimeSlope(), CandTrackHandle::GetTrackPointError(), CandTrackHandle::GetU(), CandTrackHandle::GetV(), CandRecoHandle::GetVtxPlane(), CandRecoHandle::GetVtxT(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), CandRecoHandle::GetVtxZ(), Calibrator::Instance(), CandTrackHandle::IsTPosValid(), Msg::kSynopsis, MSG, Calibrator::ReInitialise(), CandRecoHandle::SetDirCosU(), CandRecoHandle::SetDirCosV(), CandRecoHandle::SetDirCosZ(), CandTrackHandle::SetdS(), CandRecoHandle::SetEndDirCosU(), CandRecoHandle::SetEndDirCosV(), CandRecoHandle::SetEndDirCosZ(), CandRecoHandle::SetEndPlane(), CandRecoHandle::SetEndT(), CandRecoHandle::SetEndU(), CandRecoHandle::SetEndV(), CandRecoHandle::SetEndZ(), CandTrackHandle::SetMomentum(), CandFitTrackHandle::SetMomentumRange(), CandTrackHandle::SetNTimeFitDigit(), CandTrackHandle::SetRange(), AlgTrack::SetT(), CandTrackHandle::SetTimeBackwardFitNDOF(), CandTrackHandle::SetTimeBackwardFitRMS(), CandTrackHandle::SetTimeFitChi2(), CandTrackHandle::SetTimeForwardFitNDOF(), CandTrackHandle::SetTimeForwardFitRMS(), CandRecoHandle::SetTimeOffset(), CandRecoHandle::SetTimeSlope(), CandTrackHandle::SetTrackPointError(), CandTrackHandle::SetU(), CandTrackHandle::SetV(), CandRecoHandle::SetVtxPlane(), CandRecoHandle::SetVtxT(), CandRecoHandle::SetVtxU(), CandRecoHandle::SetVtxV(), CandRecoHandle::SetVtxZ(), track, vldc, and ZIncreasesWithTime.

Referenced by RunTheFitter().

03295 {
03296   // This method is only called if the fit fails. We set properties from finder track.
03297   // This clearly does not include fitted properties such as q/p or QPVtxError.
03298  MSG("AlgFitTrackCam",Msg::kSynopsis) << "SetPropertiesFromFinderTrack" << endl;
03299   cth.SetDirCosU(track->GetDirCosU());
03300   cth.SetDirCosV(track->GetDirCosV());
03301   cth.SetDirCosZ(track->GetDirCosZ());
03302   cth.SetVtxU(track->GetVtxU());
03303   cth.SetVtxV(track->GetVtxV());
03304   cth.SetVtxZ(track->GetVtxZ());
03305   cth.SetVtxT(track->GetVtxT());
03306   cth.SetVtxPlane(track->GetVtxPlane());
03307 
03308   cth.SetEndDirCosU(track->GetEndDirCosU());
03309   cth.SetEndDirCosV(track->GetEndDirCosV());
03310   cth.SetEndDirCosZ(track->GetEndDirCosZ());
03311   cth.SetEndU(track->GetEndU());
03312   cth.SetEndV(track->GetEndV());
03313   cth.SetEndZ(track->GetEndZ());
03314   cth.SetEndT(track->GetEndT());
03315   cth.SetEndPlane(track->GetEndPlane());
03316 
03317   cth.SetMomentumRange(track->GetMomentum());
03318   cth.SetMomentum(track->GetMomentum());
03319 
03320   cth.SetTimeSlope(track->GetTimeSlope());
03321   cth.SetTimeOffset(track->GetTimeOffset());
03322   cth.SetTimeFitChi2(track->GetTimeFitChi2());
03323   cth.SetNTimeFitDigit(track->GetNTimeFitDigit());
03324   cth.SetTimeForwardFitRMS(track->GetTimeForwardFitRMS());
03325   cth.SetTimeForwardFitNDOF(track->GetTimeForwardFitNDOF());
03326   cth.SetTimeBackwardFitRMS(track->GetTimeBackwardFitRMS());
03327   cth.SetTimeBackwardFitNDOF(track->GetTimeBackwardFitNDOF());
03328 
03329 
03330   // Set quantities required at each plane in finder track
03331   int direction=1;
03332   if(ZIncreasesWithTime==false) {direction=-1;}
03333 
03334   for(int i=cth.GetVtxPlane(); i*direction<=cth.GetEndPlane()*direction; i+=direction){
03335     if(track->IsTPosValid(i)) {
03336       cth.SetTrackPointError(i,track->GetTrackPointError(i));
03337       cth.SetdS(i,track->GetdS(i));
03338       cth.SetRange(i,track->GetRange(i));
03339       cth.SetU(i,track->GetU(i));
03340       cth.SetV(i,track->GetV(i));
03341     }
03342   }
03343 
03344   CalculateTrace(cth);  
03345   SetT(&cth);
03346 
03347   Calibrator& cal = Calibrator::Instance();
03348   cal.ReInitialise(*vldc);
03349   Calibrate(&cth);
03350 
03351 
03352 }

void AlgFitTrackCam::SetRangeAnddS ( CandFitTrackCamHandle cth  ) 

Definition at line 3029 of file AlgFitTrackCam.cxx.

References Munits::cm, done(), FilteredData, fPL, PlexPlaneId::GetAdjoinScint(), PlexPlaneId::GetAdjoinSteel(), VldContext::GetDetector(), CandFitTrackHandle::GetFinderTrack(), UgliSteelPlnHandle::GetHalfThickness(), UgliPlnHandle::GetHalfThickness(), CandTrackHandle::GetMomentum(), PlexPlaneId::GetPlaneView(), UgliGeomHandle::GetScintPlnHandle(), VldContext::GetSimFlag(), UgliGeomHandle::GetSteelPlnHandle(), UgliSteelPlnHandle::GetZ0(), UgliPlnHandle::GetZ0(), UgliGeomHandle::GetZExtent(), PlaneOutline::IsInside(), PlexPlaneId::IsValid(), UgliScintPlnHandle::IsValid(), PlaneCoverage::kComplete, SimFlag::kData, Msg::kDebug, Detector::kFar, Detector::kNear, PlaneCoverage::kNearFull, AtNuEventType::kPC, Msg::kSynopsis, PlaneView::kUnknown, Msg::kWarning, Munits::m, MaxPlane, MinPlane, MSG, CandTrackHandle::SetdS(), CandFitTrackHandle::SetMomentumRange(), CandTrackHandle::SetRange(), size, SlcStripData, Swim(), vldc, and ZIncreasesWithTime.

Referenced by SetTrackProperties().

03030 {
03031 
03032   // Set range and dS as calculated by the swimmer
03033   MSG("AlgFitTrackCam",Msg::kSynopsis) <<  "SetRangeAnddS from swimmer values between " << MinPlane << "/" << MaxPlane << endl;
03034 
03035   UgliGeomHandle ugh = UgliGeomHandle(*vldc);
03036 
03037   int ZDir; int VtxPlane; int EndPlane; int Increment;
03038   Detector::Detector_t detector = vldc->GetDetector();
03039   double dS=0.; double dRange=0.; double dP=0.;
03040 
03041 
03042   // Start at the end of the track and calculate the required additions to range
03044 
03045   // find ending Z position (defined as Z position where muon has 100 MeV of residual energy.  This corresponds to 1/2 inch of Fe.
03046 
03047   // NOTE: Average dP for 1" iron is 95 MeV.
03048  
03049   if(ZIncreasesWithTime==true) {ZDir=1; EndPlane=MaxPlane; VtxPlane=MinPlane; Increment=-1;}
03050   else {ZDir=-1; EndPlane=MinPlane; VtxPlane=MaxPlane; Increment=1;}
03051 
03052   PlexPlaneId plnid(detector,EndPlane,false);
03053   PlexPlaneId endplnid(detector,EndPlane,false);
03054 
03055   double Zscint = SlcStripData[EndPlane][0].csh->GetZPos();
03056   double Znextscint = Zscint;
03057 
03058   UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(plnid); 
03059   double Zend = Zscint + double(ZDir)*scintpln.GetHalfThickness();
03060 
03061   PlexPlaneId nextscint =  endplnid.GetAdjoinScint(ZDir);
03062   UgliScintPlnHandle nextscintpln = ugh.GetScintPlnHandle(nextscint);
03063   if(nextscintpln.IsValid() && nextscint.GetPlaneView()!=PlaneView::kUnknown){
03064     Znextscint = nextscintpln.GetZ0();
03065   }
03066   else{
03067     nextscint = endplnid;
03068   }
03069 
03070   plnid = plnid.GetAdjoinSteel(ZDir);
03071   if(plnid.IsValid()){
03072     UgliSteelPlnHandle steelpln = ugh.GetSteelPlnHandle(plnid);
03073     Zend = steelpln.GetZ0() - double(ZDir)*steelpln.GetHalfThickness();
03074   }
03075 
03076   // add two planes of steel for the ND spectrometer
03077   if(detector==Detector::kNear && EndPlane>=121) {
03078     for(int i=0;i<2;i++){
03079       if(plnid.GetAdjoinSteel(ZDir).IsValid()){
03080         PlexPlaneId plnid_after = plnid.GetAdjoinSteel(ZDir);
03081         if(plnid_after.IsValid()) {
03082           plnid = plnid_after;
03083           UgliSteelPlnHandle steelpln = ugh.GetSteelPlnHandle(plnid);
03084           Zend = steelpln.GetZ0() - double(ZDir)*steelpln.GetHalfThickness();
03085         }
03086       }
03087     }
03088   }
03089 
03090   // Determine whether track stops in coil
03091   float u_end = FilteredData[EndPlane][0].x_k0;
03092   float v_end = FilteredData[EndPlane][0].x_k1;
03093   float du_end = FilteredData[EndPlane][0].x_k2;
03094   float dv_end = FilteredData[EndPlane][0].x_k3;
03095   float delz = Znextscint-Zscint;
03096   float u_extrap = u_end +delz*du_end;
03097   float v_extrap = v_end +delz*dv_end;
03098   float x_extrap = 0.707*(u_extrap-v_extrap);
03099   float y_extrap = 0.707*(u_extrap+v_extrap);
03100 
03101   PlaneCoverage::PlaneCoverage_t kPC = PlaneCoverage::kComplete;
03102   if(detector==Detector::kNear) kPC=PlaneCoverage::kNearFull;
03103 
03104   bool isInOutline = fPL.IsInside(x_extrap,y_extrap,nextscint.GetPlaneView(),kPC,false);
03105   bool isInCoil = isInOutline && !fPL.IsInside(x_extrap,y_extrap,nextscint.GetPlaneView(),kPC,true);
03106 
03107   double S = 0; double Range = 0; double Prange = 0.0;
03108   double StateVector[5]; double Output[5];
03109   double chargesign = -1;
03110   bool GoForward = true; bool done=true; bool swimOK=true;
03111 
03112   // if in coil find midpoint and swim towards last hit from there
03113   if(isInCoil){
03114     float zCoil = Znextscint;
03115     float u_extrapC = u_extrap;
03116     float v_extrapC = v_extrap;
03117     float x_extrapC = x_extrap;
03118     float y_extrapC = y_extrap;
03119     while(isInCoil){
03120       zCoil -= 1.0*Munits::cm*ZDir;
03121       float delzC = zCoil - Zscint;
03122       u_extrapC = u_end +delzC*du_end;
03123       v_extrapC = v_end +delzC*dv_end;
03124       x_extrapC = 0.707*(u_extrapC-v_extrapC);
03125       y_extrapC = 0.707*(u_extrapC+v_extrapC);
03126       isInCoil = !fPL.IsInside(x_extrapC,y_extrapC,nextscint.GetPlaneView(),kPC,true) && ((zCoil>Zscint && ZDir==1) || (zCoil<Zscint && ZDir==-1)) ;        
03127     }
03128     float zMinCoil = zCoil;
03129     if(zMinCoil<Zscint && ZDir==1)  zMinCoil=Zscint;
03130     if(zMinCoil>Zscint && ZDir==-1) zMinCoil=Zscint;
03131 
03132     zCoil = Znextscint;
03133     isInCoil = true;
03134     while(isInCoil){
03135       zCoil += 1.0*Munits::cm*ZDir;
03136       float delzC = zCoil - Zscint;
03137       u_extrapC = u_end +delzC*du_end;
03138       v_extrapC = v_end +delzC*dv_end;
03139       x_extrapC = 0.707*(u_extrapC-v_extrapC);
03140       y_extrapC = 0.707*(u_extrapC+v_extrapC);
03141       float zmin; float zmax;
03142       ugh.GetZExtent(zmin,zmax);
03143       isInCoil = !fPL.IsInside(x_extrapC,y_extrapC,nextscint.GetPlaneView(),kPC,true) && ((zCoil<zmax && ZDir==1) || (zCoil>zmin && ZDir==-1));        
03144     }
03145     float zMaxCoil = zCoil;
03146     float zmin; float zmax;
03147     ugh.GetZExtent(zmin,zmax);
03148     if(zMaxCoil>zmax && ZDir==1)  zMaxCoil=zmax;
03149     if(zMaxCoil<zmin && ZDir==-1) zMaxCoil=zmin;
03150 
03151     // now swim from mid-coil back to endplane
03152     float zMidCoil = 0.5*(zMinCoil + zMaxCoil);
03153     float delzC  = zMidCoil -Zscint;
03154     u_extrapC = u_end +delzC*du_end;
03155     v_extrapC = v_end +delzC*dv_end;
03156     x_extrapC = 0.707*(u_extrapC-v_extrapC);
03157     y_extrapC = 0.707*(u_extrapC+v_extrapC);
03158     
03159     StateVector[0] = u_extrapC; Output[0]=StateVector[0];
03160     StateVector[1] = v_extrapC; Output[1]=StateVector[1];
03161     StateVector[2] = FilteredData[EndPlane][0].x_k2; Output[2]=StateVector[2];
03162     StateVector[3] = FilteredData[EndPlane][0].x_k3; Output[3]=StateVector[3];
03163     chargesign = -1;
03164     if(FilteredData[EndPlane][0].x_k4!=0.) {chargesign =  FilteredData[EndPlane][0].x_k4/fabs(FilteredData[EndPlane][0].x_k4);}
03165     
03166     GoForward = !ZIncreasesWithTime;
03167     StateVector[4] = 10.*chargesign; Output[4]=StateVector[4]; 
03168 
03169     double dsdz = pow((1. + pow(StateVector[2],2) + pow(StateVector[3],2)),0.5);    
03170     // set fallback to nominal energy loss in case coil swim fails
03171     Prange = 0.095*dsdz;
03172     if(detector==Detector::kNear && EndPlane>121) Prange = 0.2*dsdz;  
03173     Prange += 0.5*dsdz*0.1*fabs(zMaxCoil-zMinCoil)*2.357*1.97;
03174 
03175     swimOK = Swim(StateVector, Output, zMidCoil, EndPlane , GoForward, &dS, &dRange, &dP);
03176    
03177     if(swimOK ){
03178       S = dS; Range = dRange; Prange = fabs(dP);
03179       cth.SetdS(EndPlane,S); 
03180       cth.SetRange(EndPlane,Range);
03181     }
03182     if(!swimOK) {Output[4] = chargesign/Prange;}
03183   
03184   }
03185 
03186   else{
03187     // normal case - track does not end in coil
03188     if((Zend<Zscint && ZDir==1) || (Zend>Zscint && ZDir==-1)) {
03189       MSG("AlgFitTrackCam",Msg::kWarning) << " Zend on wrong side of last scint! " << endl;
03190       Zend=Zscint;
03191     }
03192     
03193     // now swim to Zend
03194     StateVector[0]=FilteredData[EndPlane][0].x_k0; Output[0]=StateVector[0];
03195     StateVector[1]=FilteredData[EndPlane][0].x_k1; Output[1]=StateVector[1];
03196     StateVector[2]=FilteredData[EndPlane][0].x_k2; Output[2]=StateVector[2];
03197     StateVector[3]=FilteredData[EndPlane][0].x_k3; Output[3]=StateVector[3];
03198     StateVector[4]=FilteredData[EndPlane][0].x_k4; Output[4]=StateVector[4];
03199     chargesign = -1;
03200     if(StateVector[4]!=0.) {chargesign = StateVector[4]/fabs(StateVector[4]);}
03201     
03202     GoForward = ZIncreasesWithTime;
03203     done = Swim(StateVector, Output, EndPlane, Zend, GoForward,  &dS, &dRange, &dP);
03204     
03205     GoForward = !ZIncreasesWithTime;
03206     double dsdz = pow((1. + pow(StateVector[2],2) + pow(StateVector[3],2)),0.5); 
03207     S = 0.; Range = 10.0*dsdz; Prange = 0.095*dsdz;
03208     swimOK = false;
03209     if(done){
03210       for(int j=0;j<5;j++) {StateVector[j]=Output[j];}
03211       
03212       // now swim from Zend to EndPlane
03213       StateVector[4] = chargesign * 10.52;  // start @ P = 100 MeV (Eloss in 1/2 " Iron)
03214       swimOK = Swim(StateVector, Output, Zend, EndPlane , GoForward, &dS, &dRange, &dP);
03215       if(swimOK){
03216         S += dS; Range += dRange; Prange += fabs(dP);
03217         cth.SetdS(EndPlane,S); 
03218         cth.SetRange(EndPlane,Range);
03219       }
03220     }
03221     if(!swimOK) {Output[4] = chargesign/Prange;}
03222   }
03223   
03224   int thisplane = EndPlane;
03225   // now swim back to vertex
03226   bool firstplane=true;
03227  for(int i=EndPlane+Increment; Increment*i<=Increment*VtxPlane; i+=Increment) {
03228     if(FilteredData[i].size()>0) {
03229       double delU = FilteredData[i][0].x_k0 - StateVector[0] ;
03230       double delV = FilteredData[i][0].x_k1 - StateVector[1] ;
03231       double dSperPlane=0.;
03232       if(thisplane!=i) {dSperPlane = pow(delU*delU + delV*delV,0.5)/double(abs(thisplane-i));}
03233 
03234       // only update state vector if change in U/V is reasonable.
03235       if(dSperPlane < 1.5*Munits::m) {
03236         StateVector[0]=FilteredData[i][0].x_k0;
03237         StateVector[1]=FilteredData[i][0].x_k1;
03238         StateVector[2]=FilteredData[i][0].x_k2;
03239         StateVector[3]=FilteredData[i][0].x_k3;
03240 
03241         chargesign=-1;
03242         if(FilteredData[i][0].x_k4!=0.) {chargesign = FilteredData[i][0].x_k4/fabs(FilteredData[i][0].x_k4);}
03243       }
03244 
03245       StateVector[4] = chargesign * fabs(Output[4]);
03246       done = Swim(StateVector, Output, thisplane, i , GoForward, &dS, &dRange, &dP);
03247       if(done){
03248         S+=dS; Range+=dRange; Prange+=fabs(dP);
03249         cth.SetdS(i,S); cth.SetRange(i,Range);
03250         firstplane=false;
03251       }
03252       else {
03253         MSG("AlgFitTrackCam",Msg::kDebug) << " swim fail " << endl;
03254       }
03255       thisplane=i;
03256     }
03257   }
03258 
03259   PlexPlaneId vtxplnid(detector,VtxPlane,false);
03260   PlexPlaneId plnid_before = vtxplnid.GetAdjoinSteel(-ZDir);
03261   
03262   if(plnid_before.IsValid()) {
03263     plnid = plnid_before;
03264     UgliSteelPlnHandle steelpln = ugh.GetSteelPlnHandle(plnid);
03265     double Zstart = steelpln.GetZ0();
03266     StateVector[0]=FilteredData[VtxPlane][0].x_k0;
03267     StateVector[1]=FilteredData[VtxPlane][0].x_k1;
03268     StateVector[2]=FilteredData[VtxPlane][0].x_k2;
03269     StateVector[3]=FilteredData[VtxPlane][0].x_k3;
03270     StateVector[4]=Output[4];
03271     Swim(StateVector, Output, VtxPlane, Zstart, GoForward, &dS,&dRange,&dP);
03272     S+=dS; Range+=dRange; Prange+=fabs(dP);
03273 
03274     cth.SetRange(VtxPlane,Range);
03275     cth.SetdS(VtxPlane,S);
03276   }
03277 
03278   // if Prange < 21 GeV, use this value.  Otherwise, use finder track energy, which is somewhat less prone to gross errors.
03279  
03280   // apply fudge factor for nominal steel thickness in ND geometry.
03281   //********* !!!!!!!!!!!!! ***********
03282   float ecorr = 1.0;
03283  if(vldc->GetDetector()==Detector::kNear && vldc->GetSimFlag()==SimFlag::kData) ecorr = 1.008; 
03284   
03285   cth.SetMomentumRange(Prange*ecorr);
03286   CandTrackHandle* findertrack = cth.GetFinderTrack();
03287   if(((detector==Detector::kFar && Prange>21.) || (detector==Detector::kNear && Prange>12.)) && findertrack) {cth.SetMomentumRange(findertrack->GetMomentum());}
03288 }

void AlgFitTrackCam::SetTrackProperties ( CandFitTrackCamHandle cth  ) 

Definition at line 2445 of file AlgFitTrackCam.cxx.

References bave, AlgTrack::CalculateTrace(), AlgReco::Calibrate(), Chi2(), EndCov, eqp_biased, FilteredData, CandHandle::GetDaughterIterator(), GetFitData(), CandFitTrackHandle::GetMomentumCurve(), CandFitTrackHandle::GetMomentumRange(), CandHandle::GetNDaughters(), CandRecoHandle::GetNDigit(), Calibrator::Instance(), CandTrackHandle::IsContained(), Msg::kSynopsis, MaxPlane, MinPlane, MSG, nbfield, NIter, Calibrator::ReInitialise(), CandFitTrackHandle::SetBave(), CandFitTrackHandle::SetChi2(), CandRecoHandle::SetDirCosU(), CandRecoHandle::SetDirCosV(), CandRecoHandle::SetDirCosZ(), CandFitTrackHandle::SetEMCharge(), CandRecoHandle::SetEndDirCosU(), CandRecoHandle::SetEndDirCosV(), CandRecoHandle::SetEndDirCosZ(), CandFitTrackHandle::SetEnddUError(), CandFitTrackHandle::SetEnddVError(), CandRecoHandle::SetEndPlane(), CandFitTrackHandle::SetEndQP(), CandFitTrackHandle::SetEndQPError(), CandRecoHandle::SetEndU(), CandFitTrackHandle::SetEndUError(), CandRecoHandle::SetEndV(), CandFitTrackHandle::SetEndVError(), CandRecoHandle::SetEndZ(), CandTrackHandle::SetInShower(), CandTrackHandle::SetMomentum(), CandFitTrackHandle::SetMomentumCurve(), CandFitTrackHandle::SetNDOF(), CandFitTrackHandle::SetNIterate(), CandFitTrackHandle::SetNSwimFail(), CandTrackHandle::SetNTrackDigit(), CandTrackHandle::SetNTrackStrip(), CandFitTrackHandle::SetPlaneChi2(), CandFitTrackHandle::SetPlaneQP(), SetRangeAnddS(), CandFitTrackCamHandle::SetRangeBiasedEQP(), CandFitTrackCamHandle::SetRangeBiasedQP(), AlgTrack::SetT(), CandTrackHandle::SetTrackPointError(), CandTrackHandle::SetU(), CandTrackHandle::SetV(), CandRecoHandle::SetVtxDirCosU(), CandRecoHandle::SetVtxDirCosV(), CandRecoHandle::SetVtxDirCosZ(), CandFitTrackHandle::SetVtxdUError(), CandFitTrackHandle::SetVtxdVError(), CandRecoHandle::SetVtxPlane(), CandFitTrackHandle::SetVtxQPError(), CandRecoHandle::SetVtxU(), CandFitTrackHandle::SetVtxUError(), CandRecoHandle::SetVtxV(), CandFitTrackHandle::SetVtxVError(), CandRecoHandle::SetVtxZ(), ShowerEntryPlane, size, SlcStripData, TimingFit(), TotalNSwimFail, TrkStripData, vldc, VtxCov, x_k4_biased, and ZIncreasesWithTime.

Referenced by RunTheFitter().

02446 { 
02447   int VtxPlane;
02448   int EndPlane;
02449   double dsdz=0.;
02450   double VtxQP=0.;
02451 
02452   if(nbfield>0) bave /=nbfield;
02453   
02454   // Carry out the assignment of variables to the new fitted track 
02455   MSG("AlgFitTrackCam",Msg::kSynopsis) << "Set Track Properties:" << endl;
02456  
02457   // Vertex and End Plane
02459 
02460   if(ZIncreasesWithTime==true) {VtxPlane=MinPlane; EndPlane=MaxPlane;}
02461   else {VtxPlane=MaxPlane; EndPlane=MinPlane;}
02462 
02463   // Momentum, charge and error on q/p
02465 
02466   // look up q/p value in vertex plane
02467   VtxQP = FilteredData[VtxPlane][0].x_k4;
02468 
02469   // apply tweak of 1.3%
02470   VtxQP *= 1.013;
02471 
02472   // if(x_k[4]!=0.) {cth.SetMomentumCurve(fabs(1./x_k[4]));}
02473   // cth.SetRangeBiasedQP(x_k4_biased);
02474   // if(x_k[4]>0.) {cth.SetEMCharge(1.);}
02475   // else if(x_k[4]<0.) {cth.SetEMCharge(-1.);}
02476   // cth.SetVtxQPError(pow(VtxCov[4],0.5));
02477 
02478 // Protect from curvature below float precision
02479   Double_t pcabsinvtemp = TMath::Max(TMath::Abs(VtxQP), (Double_t) FLT_MIN);
02480   cth.SetMomentumCurve(1./pcabsinvtemp);
02481 
02482   cth.SetRangeBiasedQP(x_k4_biased);
02483   cth.SetRangeBiasedEQP(eqp_biased);
02484 
02485   if( VtxQP>0.0 ) {cth.SetEMCharge(1.);}
02486   else if( VtxQP<0.0 ) {cth.SetEMCharge(-1.);}
02487   cth.SetVtxQPError(pow(VtxCov[4],0.5));
02488 
02489   MSG("AlgFitTrackCam",Msg::kSynopsis)  << " Q/P = " << VtxQP << ", Q/P(biased) = " << x_k4_biased << ", sigma(Q/P) = " << pow(VtxCov[4],0.5) << " Sigma Q/P _biased "<< eqp_biased << endl; 
02490   
02491   // Positions and angles  
02493 
02494   // Vtx and end coordinates  
02495   cth.SetVtxU(FilteredData[VtxPlane][0].x_k0);
02496   cth.SetVtxV(FilteredData[VtxPlane][0].x_k1);
02497   cth.SetVtxZ(SlcStripData[VtxPlane][0].csh->GetZPos());
02498   cth.SetVtxPlane(VtxPlane);
02499   
02500   cth.SetEndU(FilteredData[EndPlane][0].x_k0);
02501   cth.SetEndV(FilteredData[EndPlane][0].x_k1);
02502 
02503   cth.SetEndZ(SlcStripData[EndPlane][0].csh->GetZPos());
02504   cth.SetEndPlane(EndPlane);
02505   
02506   
02507   // Vtx and end direction cosines
02508   dsdz=pow(1.+pow(FilteredData[VtxPlane][0].x_k2,2)+pow(FilteredData[VtxPlane][0].x_k3,2),0.5);
02509   if(ZIncreasesWithTime==false) {dsdz=-dsdz;}
02510   
02511   cth.SetVtxDirCosU(FilteredData[VtxPlane][0].x_k2/dsdz);
02512   cth.SetVtxDirCosV(FilteredData[VtxPlane][0].x_k3/dsdz);
02513   cth.SetVtxDirCosZ(1./dsdz);
02514   
02515   cth.SetDirCosU(FilteredData[VtxPlane][0].x_k2/dsdz);
02516   cth.SetDirCosV(FilteredData[VtxPlane][0].x_k3/dsdz);
02517   cth.SetDirCosZ(1./dsdz);
02518     
02519   dsdz=pow(1.+pow(FilteredData[EndPlane][0].x_k2,2)+pow(FilteredData[EndPlane][0].x_k3,2),0.5);
02520   if(ZIncreasesWithTime==false) {dsdz=-dsdz;}
02521   
02522   cth.SetEndDirCosU(FilteredData[EndPlane][0].x_k2/dsdz);
02523   cth.SetEndDirCosV(FilteredData[EndPlane][0].x_k3/dsdz);
02524   cth.SetEndDirCosZ(1./dsdz);
02525   
02526   // End q/p value
02527   cth.SetEndQP(FilteredData[EndPlane][0].x_k4);
02528     
02529   // Errors on vtx positions and angles  
02530   cth.SetVtxUError(pow(VtxCov[0],0.5));
02531   cth.SetVtxVError(pow(VtxCov[1],0.5));
02532   cth.SetVtxdUError(pow(VtxCov[2],0.5));
02533   cth.SetVtxdVError(pow(VtxCov[3],0.5)); 
02534   
02535   // Errors on end positions, angles and q/p
02536   cth.SetEndUError(pow(EndCov[0],0.5));
02537   cth.SetEndVError(pow(EndCov[1],0.5));
02538   cth.SetEnddUError(pow(EndCov[2],0.5));
02539   cth.SetEnddVError(pow(EndCov[3],0.5)); 
02540   cth.SetEndQPError(pow(EndCov[4],0.5)); 
02541 
02543 
02544   cth.SetBave(bave);
02545 
02546   // More variables to be set, in order to ensure compatibility
02548   cth.SetNTrackStrip(cth.GetNDaughters());
02549   cth.SetNTrackDigit(cth.GetNDigit());
02550   
02551   cth.SetNIterate(NIter);
02552   cth.SetNSwimFail(TotalNSwimFail);
02553   
02554   
02555   // Obtain "fitting data" for the final track strips
02556   for (int i=0; i<490; ++i) {TrkStripData[i].clear();} 
02557   GetFitData(MinPlane,MaxPlane);
02558   
02559   
02560   // Set tpos error and Calculate chi2, NDOF
02561   double Chi2=0; double Chi2Contrib=0; int NDOF=0; double FilteredTPos=0;
02562   
02563   for(int i=MinPlane; i<=MaxPlane; ++i) {
02564     if(TrkStripData[i].size()>0) {
02565       
02566       if(TrkStripData[i][0].TPosError>0.) {
02567         cth.SetTrackPointError(i,pow(TrkStripData[i][0].TPosError,0.5));
02568         
02569         if(TrkStripData[i][0].PlaneView==2) {FilteredTPos=FilteredData[i][0].x_k0;}
02570         else {FilteredTPos=FilteredData[i][0].x_k1;}
02571         
02572         Chi2Contrib = pow((TrkStripData[i][0].TPos-FilteredTPos),2) / TrkStripData[i][0].TPosError;
02573         cth.SetPlaneChi2(i,Chi2Contrib);        
02574         
02575         Chi2+=Chi2Contrib;
02576         
02577         NDOF++;
02578       }
02579     }
02580   }
02581   cth.SetChi2(Chi2);
02582   cth.SetNDOF(NDOF-5); // Number of constraints set to 5
02583   
02584   
02585   // Assign U, V and q/p values
02586   for(int i=MinPlane; i<=MaxPlane; ++i) {
02587     if(FilteredData[i].size()>0) {
02588       cth.SetU(i,FilteredData[i][0].x_k0);
02589       cth.SetV(i,FilteredData[i][0].x_k1);
02590       cth.SetPlaneQP(i,FilteredData[i][0].x_k4);
02591     }
02592   }
02593   
02594   
02595   // Set Trace and TraceZ
02596   CalculateTrace(cth);
02597   
02598   
02599   // Fill time and range maps
02600   SetT(&cth);
02601   SetRangeAnddS(cth);
02602   
02603   
02604   // Set momentum to our best estimate (range or curvature)
02605   cth.SetMomentum(cth.GetMomentumCurve());
02606   if(cth.IsContained()){
02607     cth.SetMomentum(cth.GetMomentumRange());
02608   }
02609   
02610   
02611   // Identify track strips inside in vertex shower
02612   
02613   TIter FitTrackStripItr = cth.GetDaughterIterator();
02614   while(CandStripHandle* FitTrackStrip = dynamic_cast<CandStripHandle*>(FitTrackStripItr()))
02615     {
02616       if(ShowerEntryPlane!=-99) {
02617         if( (ZIncreasesWithTime==true && FitTrackStrip->GetPlane()<=ShowerEntryPlane)
02618             || (ZIncreasesWithTime==false && FitTrackStrip->GetPlane()>=ShowerEntryPlane) ) {
02619           cth.SetInShower(FitTrackStrip,2);
02620         }
02621         else {cth.SetInShower(FitTrackStrip,0);}
02622       }
02623       else {cth.SetInShower(FitTrackStrip,0);}
02624     }
02625   
02626   
02627   // Set all time related variables
02628   TimingFit(cth);
02629   
02630   
02631   // Calibrate, to get track pulse height in MIPs, GeV, etc
02632   Calibrator& cal = Calibrator::Instance();
02633   cal.ReInitialise(*vldc);
02634   Calibrate(&cth);
02635 
02636 
02637 }

void AlgFitTrackCam::ShowerStrips (  ) 

Definition at line 327 of file AlgFitTrackCam.cxx.

References FilteredData, Msg::kDebug, MaxPlane, min, MinPlane, MSG, ShowerEntryPlane, size, SlcStripData, SwimThroughShower, and ZIncreasesWithTime.

Referenced by RunTheFitter().

00328 {
00329   MSG("AlgFitTrackCam",Msg::kDebug) << "ShowerStrips, Look for large vertex shower" << endl;
00330      
00331   // Initialisations
00332   int Increment; int NumberOfStrips;
00333   int Plane; int NewPlane;
00334 
00335   int VtxShwWindow=8; 
00336   int StripsForShw=4;
00337   double PEThreshold=2.;
00338 
00339   if(ZIncreasesWithTime==true) {Plane=MinPlane; Increment=1;}
00340   else {Plane=MaxPlane; Increment=-1;}
00341   NewPlane=Plane;
00342 
00343 
00344   // Identify any vertex showers
00346   while(abs(Plane-NewPlane)<=VtxShwWindow && NewPlane>=MinPlane && NewPlane<=MaxPlane) {
00347 
00348     if(SlcStripData[NewPlane].size()>0) {
00349       NumberOfStrips=0;
00350 
00351 
00352       // Set the number of hits on a plane required for the plane to be identified as 'in the 
00353       // shower'. We account for the gradient of the track, with the factor of 0.25 representing
00354       // the approximate ratio of strip thickness to strip width.
00355       if(FilteredData[NewPlane].size()>0) {
00356         if(SlcStripData[NewPlane][0].csh->GetPlaneView()==2) {
00357           StripsForShw=min(7,int( 4+(0.25*fabs(FilteredData[NewPlane][0].x_k2)) ));
00358         }
00359         else {StripsForShw=min(7,int( 4+(0.25*fabs(FilteredData[NewPlane][0].x_k3)) ));}
00360       }
00361       else {StripsForShw=4;}
00362 
00363 
00364       // Count number of strips on plane with greater than 2PEs
00365       for(unsigned int j=0; j<SlcStripData[NewPlane].size(); ++j) {
00366         if(SlcStripData[NewPlane][j].csh->GetCharge()>PEThreshold) {NumberOfStrips++;}
00367       }
00368 
00369 
00370       // If a vertex shower is found, note that we should use the Swimmer
00371       // to find the most likely track strips inside the shower
00372       if(NumberOfStrips>=StripsForShw) {ShowerEntryPlane=NewPlane; SwimThroughShower=true; break;}
00373       
00374       NewPlane+=Increment;
00375     }
00376     else {NewPlane+=Increment;}
00377   }
00379 
00380 
00381  
00382   // Find the plane at which the 'clean' section of track enters the shower
00384   if(SwimThroughShower==true) {
00385     NewPlane=ShowerEntryPlane+Increment;
00386     int PlanesSinceLastHit=0;
00387     int PlaneWindow=4;
00388 
00389     while(PlanesSinceLastHit<PlaneWindow && NewPlane>=MinPlane && NewPlane<=MaxPlane) {
00390       if(SlcStripData[NewPlane].size()>0) {
00391         NumberOfStrips=0;
00392 
00393         // Account for gradient of track, as before
00394         if(FilteredData[NewPlane].size()>0) { 
00395           if(SlcStripData[NewPlane][0].csh->GetPlaneView()==2) {
00396             StripsForShw=min(7,int(4+(0.25*fabs(FilteredData[NewPlane][0].x_k2)) ));
00397           }
00398           else {StripsForShw=min(7,int(4+(0.25*fabs(FilteredData[NewPlane][0].x_k3)) ));}
00399         }
00400         else {StripsForShw=4;}
00401 
00402 
00403         // Count number of strips on plane with greater than 2PEs
00404         for(unsigned int j=0; j<SlcStripData[NewPlane].size(); ++j) {
00405           if(SlcStripData[NewPlane][j].csh->GetCharge()>PEThreshold) {NumberOfStrips++;}
00406         }
00407         if(NumberOfStrips>=StripsForShw) {
00408           ShowerEntryPlane=NewPlane; NewPlane+=Increment; PlanesSinceLastHit=0;
00409         }
00410         else {PlanesSinceLastHit++; NewPlane+=Increment;}
00411         
00412       }
00413       else {PlanesSinceLastHit++; NewPlane+=Increment;}
00414     }
00415   }
00417 }

void AlgFitTrackCam::ShowerSwim (  ) 

Definition at line 698 of file AlgFitTrackCam.cxx.

References CalcKalmanGain(), StripStruct::csh, ExtrapCovMatrix(), GetCombiPropagator(), GetFitData(), GetNoiseMatrix(), H_k, InitTrkStripData, Msg::kDebug, MaxPlane, MinPlane, MoveArrays(), MSG, size, SlcStripData, StoreFilteredData(), Swim(), SwimThroughShower, UpdateCovMatrix(), UpdateStateVector(), x_k, x_k_minus, and ZIncreasesWithTime.

Referenced by RunTheFitter().

00699 {
00700   // Method is called if we have a large shower near the track vertex
00701   //
00702   // The Swimmer is used to find the most likely track strip in the shower
00703   // and this strip is added to the fit
00704   MSG("AlgFitTrackCam",Msg::kDebug) << "ShowerSwim, improved track finding in shower" << endl;
00705 
00706   // Initialisations
00707   int Plane; int NewPlane;
00708   double StateVector[5]; double NState[5];
00709   bool GoForward; bool SwimBack;
00710   int PlanesSinceLastHit=0;
00711   int PlaneView;
00712   int Increment;
00713   
00714   double StripDistance=0.; double MinDistanceToStrip=0.;
00715   double StripWidth=4.108e-2;
00716 
00717 
00718   if(ZIncreasesWithTime==true) {GoForward=false; Plane=MinPlane; Increment=-1;}
00719   else {GoForward=true; Plane=MaxPlane; Increment=1;}
00720 
00721   NewPlane=Plane+Increment;
00722 
00723 
00724   // Continue until we reach a 4 plane window with no likely hit or we reach 
00725   // the end of the detector
00726   while(PlanesSinceLastHit<4 && NewPlane>0 && NewPlane<=485) {
00727     if(SlcStripData[NewPlane].size()>0) {
00728 
00729       if(SlcStripData[NewPlane][0].csh->GetPlaneView()==2) {PlaneView=0;}
00730       else {PlaneView=1;}
00731       
00732       for(int i=0; i<5; ++i) {StateVector[i]=x_k_minus[i];}
00733       SwimBack=Swim(StateVector, NState, Plane, NewPlane, GoForward);
00734       if(!SwimBack){break;}
00735       for(int i=0; i<5; ++i) {x_k[i]=NState[i];}
00736  
00737  
00738       // Find the closest strip (within a distance 'MinDistanceToStrip') and
00739       // temporarily store CandStripHandle
00740       // Results are very sensitive to value of MinDistanceToStrip
00741       CandStripHandle* CurrentStrip=0;
00742       MinDistanceToStrip=(1.5*StripWidth)+ fabs(0.0055*x_k[PlaneView+2]);
00743       
00744       for(unsigned int j=0; j<SlcStripData[NewPlane].size(); ++j) {
00745         StripDistance=fabs(SlcStripData[NewPlane][j].csh->GetTPos()-x_k[PlaneView]);
00746         
00747         if(StripDistance<MinDistanceToStrip) {
00748           MinDistanceToStrip=StripDistance;
00749           CurrentStrip=SlcStripData[NewPlane][j].csh;
00750         }
00751       }            
00752       
00753       // If we find a likely track strip, add it to the fit data and call the Kalman
00754       // update equations before repeating process to find next track strips in the shower
00755       if(CurrentStrip) {
00756         StripStruct temp;
00757         temp.csh = CurrentStrip;
00758         InitTrkStripData[NewPlane].push_back(temp);
00759         
00760         // Convert the strip to data required for Kalman fit
00761         GetFitData(NewPlane,NewPlane);
00762 
00763         // Carry out the Kalman fit
00764         if (PlaneView==1) {H_k[0]=0; H_k[1]=1; H_k[2]=0; H_k[3]=0; H_k[4]=0;}
00765         else {H_k[0]=1; H_k[1]=0; H_k[2]=0; H_k[3]=0; H_k[4]=0;}
00766         
00767         bool CombiPropagatorOk=GetCombiPropagator(Plane,NewPlane,GoForward);
00768         
00769         if(CombiPropagatorOk ) {
00770           GetNoiseMatrix(Plane,NewPlane);
00771           ExtrapCovMatrix();
00772           CalcKalmanGain(NewPlane);
00773           UpdateStateVector(Plane,NewPlane,GoForward);
00774           UpdateCovMatrix();
00775           MoveArrays();
00776           StoreFilteredData(NewPlane);
00777           
00778           if(ZIncreasesWithTime) {MinPlane=NewPlane; Plane=MinPlane;}
00779           else {MaxPlane=NewPlane; Plane=MaxPlane;}
00780           NewPlane=Plane+Increment;
00781           
00782           PlanesSinceLastHit=0;
00783         }
00784       }
00785       else {NewPlane+=Increment; PlanesSinceLastHit++;}
00786       
00787     }
00788     else {NewPlane+=Increment; PlanesSinceLastHit++;}
00789   }
00790   
00791   // Note that shower swim is complete
00792   SwimThroughShower=false;
00793 
00794 }

void AlgFitTrackCam::SpectrometerSwim ( CandFitTrackCamHandle cth  ) 

Definition at line 3440 of file AlgFitTrackCam.cxx.

References CandHandle::AddDaughterLink(), C_k, CalcKalmanGain(), StripStruct::csh, EndofRangePlane, ExtrapCovMatrix(), FilteredData, GetCombiPropagator(), CandFitTrackHandle::GetFinderTrack(), GetFitData(), GetMomFromRange(), GetNoiseMatrix(), CandStripHandle::GetPlane(), PlexPlaneId::GetPlaneCoverage(), CandTrackHandle::GetRange(), CandStripHandle::GetZPos(), H_k, InitTrkStripData, IsInDetector(), CandHandle::IsSlushyEnabled(), PlexPlaneId::IsValid(), Msg::kDebug, Msg::kError, Detector::kNear, PlaneCoverage::kNoActive, MaxPlane, MeanTrackTime, MoveArrays(), MSG, NDPlaneIsActive(), NDStripBegTime(), PassTrack, CandRecoHandle::SetEndDirCosU(), CandRecoHandle::SetEndDirCosV(), CandRecoHandle::SetEndDirCosZ(), CandRecoHandle::SetEndPlane(), CandRecoHandle::SetEndU(), CandRecoHandle::SetEndV(), CandRecoHandle::SetEndZ(), CandTrackHandle::SetMomentum(), CandHandle::SetSlushyEnabled(), AlgTrack::SetUVZT(), size, SlcStripData, StoreFilteredData(), Swim(), UpdateCovMatrix(), UpdateStateVector(), x_k_minus, and ZIncreasesWithTime.

Referenced by RunTheFitter().

03441 {
03442   MSG("AlgFitTrackCam",Msg::kDebug) << "SpectrometerSwim, improved track finding in spectrometer" << endl;
03443 
03444   // Initialisations
03445   // Sort out the lists for the ND spectrometer
03446   bool AddedStrip = false;
03447   int Plane; int NewPlane;
03448   double StateVector[5]; double NState[5]; double temp_x_k[5];
03449   bool GoForward; bool SpectSwim;
03450   int ActivePlanesSinceLastHit=0;
03451   int PlaneView; int Increment; int Counter;
03452   double LastHitTimes[2]={MeanTrackTime,MeanTrackTime};
03453   double TimeWindow=40.e-9;
03454 
03455   double StripDistance; double MinDistanceToStrip;
03456   double SwimError2;
03457   double StripWidth = 4.108e-2;
03458   double PlanePitch = 0.0594;
03459 
03460   bool TrackInActiveRegion;
03461   GoForward=true; Plane=MaxPlane; Increment=1;
03462 
03463   // Linear extrapolation from end of track to start of spectrometer.
03464   // Count number of active planes  
03465   while(ActivePlanesSinceLastHit<6 && Plane<121) {
03466     Plane += Increment;
03467     double u = x_k_minus[0] + x_k_minus[2]*(Plane-MaxPlane)*PlanePitch;
03468     double v = x_k_minus[1] + x_k_minus[3]*(Plane-MaxPlane)*PlanePitch;
03469 
03470     // 15 Feb 2008 - mitigate symptoms of a problem elsewhere - huge u and v can cause crash.
03471     if(fabs(u) > 5000 || fabs(v) > 5000){
03472       MSG("AlgFitTrackCam", Msg::kError) << "SpectrometerSwim - unexpectedly large u or v (u="
03473                                          << u << " v=" << v << ") bailing out." << endl;
03474       return;
03475     }
03476 
03477     if (NDPlaneIsActive(Plane, u, v, 0.3)) ActivePlanesSinceLastHit++;
03478   }
03479 
03480   // If we are clearly not near spectrometer, return from method
03481   if(ActivePlanesSinceLastHit>5 || abs(121-MaxPlane)>=40) {return;}
03482 
03483   // Set initial positions for swim
03484   ActivePlanesSinceLastHit=0;
03485   Plane = MaxPlane;  NewPlane = Plane+1;
03486 
03487   // Continue until we reach a 8 plane gap (counting only active planes) since prior
03488   // hit or we reach the end of the detector
03489   while(ActivePlanesSinceLastHit<8 && abs(NewPlane-Plane)<=70 && NewPlane<=290 && PassTrack==true) {
03490     
03491     PlexPlaneId plexPlane(Detector::kNear,NewPlane, 0); 
03492     if(SlcStripData[NewPlane].size()>0 && plexPlane.IsValid() && plexPlane.GetPlaneCoverage()!=PlaneCoverage::kNoActive) {
03493 
03494       if(SlcStripData[NewPlane][0].csh->GetPlaneView()==2) {PlaneView=0;}
03495       else {PlaneView=1;}
03496       
03497       for(int i=0; i<5; ++i) {StateVector[i]=x_k_minus[i];}
03498 
03499       // For the purposes of spectrometer track hit finding, don't allow track to range out before
03500       // we have swum all the hit spectrometer planes in this slice
03501       SpectSwim = Swim(StateVector, NState, Plane, NewPlane, GoForward);
03502 
03503       // If swim has failed and there is a large gap to next hit plane, stop the spectrometer swim.
03504       if(!SpectSwim && (NewPlane-Plane)>=40) {
03505         break;}
03506 
03507       // If swim has failed, but there is no large gap, make a momentum correction and swim again.
03508       if(!SpectSwim && StateVector[4]!=0) {
03509         Counter=0;
03510         // Double momentum and attempt to swim again
03511         while(!SpectSwim && Counter<2) {
03512           StateVector[4]*=0.5;  
03513           SpectSwim = Swim(StateVector, NState, Plane, NewPlane, GoForward);
03514           
03515           if(!SpectSwim) {Counter++;}
03516         }
03517       }
03518       
03519       if(!SpectSwim) {break;}
03520 
03521       // Leigh: Check to see if the swimmer found a position inside the detector. If not, try the next
03522       // plane.
03523                         int test = IsInDetector(NState[0],NState[1],NewPlane);
03524       if(test==0){
03525         ActivePlanesSinceLastHit++;
03526         NewPlane+=Increment;
03527         continue;
03528       }
03529                         else if(test==2){
03530                                 // If we get here, then we want to stop the tracking.
03531                                 break;
03532                         }
03533 
03534       for(int i=0; i<5; ++i) {temp_x_k[i]=NState[i];}
03535 
03536       // Calculate error in track extrapolation, based on covariance matrix on-diagonal elements
03537       SwimError2 = C_k[PlaneView][PlaneView] + (C_k[PlaneView+2][PlaneView+2]*PlanePitch*PlanePitch*(NewPlane-Plane)*(NewPlane-Plane)); 
03538       MinDistanceToStrip = 3.0 * pow(StripWidth*StripWidth + SwimError2,0.5);
03539       
03540 
03541       // Find the closest strip (within a distance 'MinDistanceToStrip') and temporarily store CandStripHandle
03542       CandStripHandle* CurrentStrip = 0;      
03543       for(unsigned int j=0; j<SlcStripData[NewPlane].size(); ++j) {
03544         StripDistance = fabs(SlcStripData[NewPlane][j].csh->GetTPos()-temp_x_k[PlaneView]);
03545 
03546         if(StripDistance<MinDistanceToStrip 
03547            && fabs(0.5*(LastHitTimes[0]+LastHitTimes[1])-NDStripBegTime(SlcStripData[NewPlane][j].csh,temp_x_k[0],temp_x_k[1]))<TimeWindow) {
03548           MinDistanceToStrip=StripDistance;
03549           CurrentStrip=SlcStripData[NewPlane][j].csh;
03550         }
03551       }
03552 
03553       // If we find a likely track strip, add it to the fit data and call the Kalman
03554       // update equations before repeating process to find next track strips in the shower
03555       if(CurrentStrip) {
03556         LastHitTimes[1]=LastHitTimes[0];
03557         LastHitTimes[0]=NDStripBegTime(CurrentStrip,temp_x_k[0],temp_x_k[1]);
03558 
03559         StripStruct temp;
03560         temp.csh = CurrentStrip;
03561         InitTrkStripData[NewPlane].push_back(temp);
03562         
03563         // Convert the strip to data required for Kalman fit
03564         GetFitData(NewPlane,NewPlane);
03565         
03566         // Carry out the Kalman fit
03567         if (PlaneView==1) {H_k[0]=0; H_k[1]=1; H_k[2]=0; H_k[3]=0; H_k[4]=0;}
03568         else {H_k[0]=1; H_k[1]=0; H_k[2]=0; H_k[3]=0; H_k[4]=0;}
03569         
03570         bool CombiPropagatorOk = GetCombiPropagator(Plane,NewPlane,GoForward);
03571 
03572         if(CombiPropagatorOk) {
03573           GetNoiseMatrix(Plane,NewPlane);
03574           ExtrapCovMatrix();
03575           CalcKalmanGain(NewPlane);
03576           UpdateStateVector(Plane,NewPlane,GoForward);
03577           UpdateCovMatrix();
03578           MoveArrays();
03579           StoreFilteredData(NewPlane);
03580           
03581           MaxPlane=NewPlane; Plane=MaxPlane;
03582           NewPlane=Plane+Increment;
03583           
03584           ActivePlanesSinceLastHit=0;
03585         }
03586 
03587         // Extend finder track, including the ND Spectrometer hits found in the fit
03589 
03590        
03591         CandTrackHandle * findertrack = cth.GetFinderTrack();
03592         if(findertrack) {
03593           bool SlushyOnEntry = CandHandle::IsSlushyEnabled();
03594           CandHandle::SetSlushyEnabled(kTRUE);
03595           AddedStrip = true;
03596           findertrack->AddDaughterLink(*CurrentStrip);
03597           findertrack->SetEndPlane(CurrentStrip->GetPlane());
03598           findertrack->SetEndZ(CurrentStrip->GetZPos());
03599           findertrack->SetEndU(FilteredData[CurrentStrip->GetPlane()][0].x_k0);
03600           findertrack->SetEndV(FilteredData[CurrentStrip->GetPlane()][0].x_k1);
03601           double dsdz = sqrt(1+pow(FilteredData[CurrentStrip->GetPlane()][0].x_k2,2)+pow(FilteredData[CurrentStrip->GetPlane()][0].x_k3,2));
03602           if(!ZIncreasesWithTime) dsdz=-dsdz;
03603           findertrack->SetEndDirCosU(FilteredData[CurrentStrip->GetPlane()][0].x_k2/dsdz);
03604           findertrack->SetEndDirCosV(FilteredData[CurrentStrip->GetPlane()][0].x_k3/dsdz);
03605           findertrack->SetEndDirCosZ(1/dsdz);
03606           if(!SlushyOnEntry) CandHandle::SetSlushyEnabled(kFALSE);
03607         }
03609         
03610       }
03611       else {
03612         TrackInActiveRegion=NDPlaneIsActive(NewPlane, temp_x_k[0], temp_x_k[1], 0.3);
03613         if(plexPlane.IsValid() && plexPlane.GetPlaneCoverage()!=PlaneCoverage::kNoActive && TrackInActiveRegion) 
03614           {ActivePlanesSinceLastHit++;}
03615         NewPlane+=Increment;
03616       }
03617     }
03618     else {
03619       double u = x_k_minus[0] + x_k_minus[2]*(NewPlane-Plane)*PlanePitch;
03620       double v = x_k_minus[1] + x_k_minus[3]*(NewPlane-Plane)*PlanePitch;
03621       
03622       TrackInActiveRegion=NDPlaneIsActive(NewPlane, u, v, 0.3); 
03623       
03624       if(plexPlane.IsValid() && plexPlane.GetPlaneCoverage()!=PlaneCoverage::kNoActive && TrackInActiveRegion) 
03625         {ActivePlanesSinceLastHit++;}
03626       NewPlane+=Increment; 
03627     }
03628   }
03629 
03630   
03631   // Sort out range and dS for finder track
03633   if(AddedStrip) {
03634     CandTrackHandle * findertrack = cth.GetFinderTrack();
03635     if(findertrack) {
03636       bool SlushyOnEntry = CandHandle::IsSlushyEnabled();
03637       CandHandle::SetSlushyEnabled(kTRUE);
03638       
03639       SetUVZT(findertrack);
03640       Double_t range = findertrack->GetRange();
03641       if (range>0.) {
03642         findertrack->SetMomentum(GetMomFromRange(range));   
03643       }
03644       if(!SlushyOnEntry) CandHandle::SetSlushyEnabled(kFALSE);
03645     }
03646   }
03648   EndofRangePlane=MaxPlane;
03649 }

void AlgFitTrackCam::StoreFilteredData ( const int  NewPlane  ) 

Definition at line 2424 of file AlgFitTrackCam.cxx.

References FilteredData, Msg::kSynopsis, Msg::kVerbose, MSG, x_k, FiltDataStruct::x_k0, FiltDataStruct::x_k1, FiltDataStruct::x_k2, FiltDataStruct::x_k3, and FiltDataStruct::x_k4.

Referenced by GoBackwards(), GoForwards(), RunTheFitter(), ShowerSwim(), and SpectrometerSwim().

02425 {
02426   // Store the data required for matching Kalman output data to strips
02427   MSG("AlgFitTrackCam",Msg::kVerbose) << "StoreFilteredData" << endl;
02428 
02429   FiltDataStruct temp;
02430   
02431   temp.x_k0=x_k[0]; temp.x_k1=x_k[1];
02432   temp.x_k2=x_k[2]; temp.x_k3=x_k[3];
02433   temp.x_k4=x_k[4];
02434 
02435   //  FilteredData[NewPlane].push_back(temp);
02436   FilteredData[NewPlane].insert(FilteredData[NewPlane].begin(),temp);
02437 
02438   MSG("AlgFitTrackCam",Msg::kSynopsis) << " StoreData [" << NewPlane << "]  x_k[4] = " << x_k[4] << " Entry=" << FilteredData[NewPlane].size() << endl; 
02439 }

bool AlgFitTrackCam::Swim ( double *  StateVector,
double *  Output,
const int  Plane,
const double  zend,
const bool  GoForward,
double *  dS = 0,
double *  Range = 0,
double *  dE = 0 
)

Definition at line 1762 of file AlgFitTrackCam.cxx.

References done(), MuELoss::e, SwimParticle::GetCharge(), SwimParticle::GetDirection(), SwimParticle::GetMomentumModulus(), SwimParticle::GetPosition(), SwimParticle::GetRange(), SwimParticle::GetS(), GeoSwimmer::Initialize(), GeoSwimmer::Instance(), SwimParticle::SetCharge(), SlcStripData, GeoSwimmer::SwimBackward(), SwimSwimmer::SwimBackward(), SwimSwimmer::SwimForward(), GeoSwimmer::SwimForward(), UseGeoSwimmer, vldc, and ZIncreasesWithTime.

01764 {
01765   //  MSG("AlgFitTrackCam",Msg::kDebug) << "Swim, specified end Z" << endl;
01766 
01767   // Initialisations
01768   // customize for bfield scaling.
01769 
01770   BField * bf = new BField(*vldc,-1,0);
01771   SwimSwimmer* myswimmer = new SwimSwimmer(*vldc,bf);
01772 
01773   if(UseGeoSwimmer) GeoSwimmer::Instance()->Initialize(*vldc);
01774 
01775   double invSqrt2 = pow(1./2.,0.5);
01776   double charge = 0.;
01777   bool done = false;
01778     
01779   if(fabs(StateVector[4])>1.e-10) {
01780     double modp = fabs(1./StateVector[4]);
01781     
01782     // Fix, to account for fact the cosmic muons could move in direction of negative z
01783     if(ZIncreasesWithTime==false) {modp=-modp;}
01784     
01785     double dsdz = pow((1.+pow(StateVector[2],2)+pow(StateVector[3],2)),0.5);
01786     double dxdz = invSqrt2*(StateVector[2]-StateVector[3]);
01787     double dydz = invSqrt2*(StateVector[2]+StateVector[3]);
01788   
01789     // Set up current muon details
01790     if(StateVector[4]>0.) charge = 1.;
01791     else if(StateVector[4]<0.) charge = -1.;
01792     
01793     TVector3 position(invSqrt2*(StateVector[0]-StateVector[1]),
01794                       invSqrt2*(StateVector[0]+StateVector[1]),
01795                       SlcStripData[Plane][0].csh->GetZPos());
01796 
01797 
01798     TVector3 momentum(modp*(dxdz/dsdz),
01799                       modp*(dydz/dsdz),
01800                       modp/dsdz);
01801     SwimParticle muon(position,momentum);
01802     muon.SetCharge(charge);
01803     SwimZCondition zc(Zend);
01804  
01805     // Do the swim, accounting for direction of motion w.r.t time too
01806     if( (GoForward==true && ZIncreasesWithTime==true)  || (GoForward==false && ZIncreasesWithTime==false) ) {
01807       if(UseGeoSwimmer){
01808         done = GeoSwimmer::Instance()->SwimForward(muon,Zend);
01809       } else {
01810         done = myswimmer->SwimForward(muon,zc);
01811       }   
01812     }
01813     else if( (GoForward==true && ZIncreasesWithTime==false)  || (GoForward==false && ZIncreasesWithTime==true) ) {
01814       if(UseGeoSwimmer){
01815         done = GeoSwimmer::Instance()->SwimBackward(muon,Zend);
01816       } else {
01817         done = myswimmer->SwimBackward(muon,zc);
01818       }    
01819     }
01820     if(done==true) {
01821       if(muon.GetDirection().Z()!=0. && muon.GetMomentumModulus()!=0.) {
01822         Output[0]=invSqrt2*(muon.GetPosition().Y()+muon.GetPosition().X());
01823         Output[1]=invSqrt2*(muon.GetPosition().Y()-muon.GetPosition().X());
01824         Output[2]=invSqrt2*((muon.GetDirection().Y()/muon.GetDirection().Z())+(muon.GetDirection().X()/muon.GetDirection().Z()));
01825         Output[3]=invSqrt2*((muon.GetDirection().Y()/muon.GetDirection().Z())-(muon.GetDirection().X()/muon.GetDirection().Z()));
01826         Output[4]=muon.GetCharge()/muon.GetMomentumModulus();
01827         // Get range and dS from the Swimmer
01828         if(dS) {*dS=muon.GetS();} if(Range) {*Range=muon.GetRange();} if(dE){*dE=muon.GetMomentumModulus()-momentum.Mag();}
01829       }
01830       else {done=false;}
01831     }
01832     
01833   }
01834 
01835   else{
01836     // If infinite momentum, use straight line extrapolation
01837     double delz = (Zend-SlcStripData[Plane][0].csh->GetZPos());
01838     Output[0]=StateVector[0] + StateVector[2]*delz;
01839     Output[1]=StateVector[1] + StateVector[3]*delz;
01840     Output[2]=StateVector[2];
01841     Output[3]=StateVector[3];
01842     Output[4]=StateVector[4];
01843 
01844     done=true;
01845   }    
01846   
01847   delete myswimmer;
01848   delete bf;
01849   return done;
01850 }

bool AlgFitTrackCam::Swim ( double *  StateVector,
double *  Output,
const double  zbeg,
const int  NewPlane,
const bool  GoForward,
double *  dS = 0,
double *  Range = 0,
double *  dE = 0 
)

Definition at line 1667 of file AlgFitTrackCam.cxx.

References done(), MuELoss::e, SwimParticle::GetCharge(), SwimParticle::GetDirection(), SwimParticle::GetMomentumModulus(), SwimParticle::GetPosition(), SwimParticle::GetRange(), SwimParticle::GetS(), GeoSwimmer::Initialize(), GeoSwimmer::Instance(), SwimParticle::SetCharge(), SlcStripData, GeoSwimmer::SwimBackward(), SwimSwimmer::SwimBackward(), SwimSwimmer::SwimForward(), GeoSwimmer::SwimForward(), UseGeoSwimmer, vldc, and ZIncreasesWithTime.

01669 {
01670   //  MSG("AlgFitTrackCam",Msg::kDebug) << "Swim, specified starting Z" << endl;
01671 
01672   // Initialisations
01673   // customize for bfield scaling.
01674   BField * bf = new BField(*vldc,-1,0);
01675   SwimSwimmer* myswimmer = new SwimSwimmer(*vldc,bf);
01676   
01677   if(UseGeoSwimmer) GeoSwimmer::Instance()->Initialize(*vldc);
01678 
01679   double invSqrt2 = pow(1./2.,0.5);
01680   double charge = 0.;
01681   bool done = false;
01682     
01683   if(fabs(StateVector[4])>1.e-10) {
01684     double modp = fabs(1./StateVector[4]);
01685     
01686     // Fix, to account for fact the cosmic muons could move in direction of negative z
01687     if(ZIncreasesWithTime==false) {modp=-modp;}
01688     
01689     double dsdz = pow((1.+pow(StateVector[2],2)+pow(StateVector[3],2)),0.5);
01690     double dxdz = invSqrt2*(StateVector[2]-StateVector[3]);
01691     double dydz = invSqrt2*(StateVector[2]+StateVector[3]);
01692     
01693     // Set up current muon details
01694     if(StateVector[4]>0.) charge = 1.;
01695     else if(StateVector[4]<0.) charge = -1.;
01696     
01697     TVector3 position(invSqrt2*(StateVector[0]-StateVector[1]),
01698                       invSqrt2*(StateVector[0]+StateVector[1]),
01699                       z);
01700 
01701     TVector3 momentum(modp*(dxdz/dsdz),
01702                       modp*(dydz/dsdz),
01703                       modp/dsdz);
01704     SwimParticle muon(position,momentum);
01705     muon.SetCharge(charge);
01706     SwimZCondition zc(SlcStripData[NewPlane][0].csh->GetZPos());
01707 
01708    
01709     // Do the swim, accounting for direction of motion w.r.t time too
01710     if( (GoForward==true && ZIncreasesWithTime==true)  || (GoForward==false && ZIncreasesWithTime==false) ) {
01711       if(UseGeoSwimmer) {
01712         done = GeoSwimmer::Instance()->SwimForward(muon,SlcStripData[NewPlane][0].csh->GetZPos());
01713       } else {
01714         done = myswimmer->SwimForward(muon,zc);
01715       } 
01716     }
01717     else if( (GoForward==true && ZIncreasesWithTime==false)  || (GoForward==false && ZIncreasesWithTime==true) ) {
01718       if(UseGeoSwimmer) {
01719         done = GeoSwimmer::Instance()->SwimBackward(muon,SlcStripData[NewPlane][0].csh->GetZPos());
01720 
01721       } else {
01722         done = myswimmer->SwimBackward(muon,zc);
01723       }     
01724     }
01725     if(done==true) {
01726       if(muon.GetDirection().Z()!=0. && muon.GetMomentumModulus()!=0.) {
01727         Output[0]=invSqrt2*(muon.GetPosition().Y()+muon.GetPosition().X());
01728         Output[1]=invSqrt2*(muon.GetPosition().Y()-muon.GetPosition().X());
01729         Output[2]=invSqrt2*((muon.GetDirection().Y()/muon.GetDirection().Z())+(muon.GetDirection().X()/muon.GetDirection().Z()));
01730         Output[3]=invSqrt2*((muon.GetDirection().Y()/muon.GetDirection().Z())-(muon.GetDirection().X()/muon.GetDirection().Z()));
01731         Output[4]=muon.GetCharge()/muon.GetMomentumModulus();
01732         // Get range and dS from the Swimmer
01733         if(dS) {*dS=muon.GetS();} if(Range) {*Range=muon.GetRange();} if(dE){*dE=muon.GetMomentumModulus()-momentum.Mag();} 
01734       }
01735       else {done=false;}
01736     }
01737     
01738   }
01739 
01740   else{
01741     // If infinite momentum, use straight line extrapolation
01742     double delz = (SlcStripData[NewPlane][0].csh->GetZPos()-z);
01743     Output[0]=StateVector[0] + StateVector[2]*delz;
01744     Output[1]=StateVector[1] + StateVector[3]*delz;
01745     Output[2]=StateVector[2];
01746     Output[3]=StateVector[3];
01747     Output[4]=StateVector[4];
01748 
01749     done=true;
01750   }    
01751   
01752   delete myswimmer;
01753   delete bf;
01754 
01755   return done;
01756 }

bool AlgFitTrackCam::Swim ( double *  StateVector,
double *  Output,
const int  Plane,
const int  NewPlane,
const bool  GoForward,
double *  dS = 0,
double *  Range = 0,
double *  dE = 0 
)

Definition at line 1571 of file AlgFitTrackCam.cxx.

References bave, done(), MuELoss::e, BField::GetBField(), SwimParticle::GetCharge(), SwimParticle::GetDirection(), SwimParticle::GetMomentumModulus(), SwimParticle::GetPosition(), SwimParticle::GetRange(), SwimParticle::GetS(), GeoSwimmer::Initialize(), GeoSwimmer::Instance(), nbfield, SwimParticle::SetCharge(), SlcStripData, GeoSwimmer::SwimBackward(), SwimSwimmer::SwimBackward(), SwimSwimmer::SwimForward(), GeoSwimmer::SwimForward(), UseGeoSwimmer, vldc, and ZIncreasesWithTime.

Referenced by FillGapsInTrack(), GetCombiPropagator(), SetRangeAnddS(), ShowerSwim(), SpectrometerSwim(), and UpdateStateVector().

01573 {
01574   //  MSG("AlgFitTrackCam",Msg::kDebug) << "Swim" << endl;
01575 
01576   // Initialisations
01577   // customize for bfield scaling.
01578   BField * bf = new BField(*vldc,-1,0);
01579   SwimSwimmer* myswimmer = new SwimSwimmer(*vldc,bf);
01580   
01581   if(UseGeoSwimmer) GeoSwimmer::Instance()->Initialize(*vldc);
01582  
01583   double invSqrt2 = pow(1./2.,0.5);
01584   double charge = 0.;
01585   bool done = false;
01586     
01587   if(fabs(StateVector[4])>1.e-10) {
01588     double modp = fabs(1./StateVector[4]);
01589     
01590     // Fix, to account for fact the cosmic muons could move in direction of negative z
01591     if(ZIncreasesWithTime==false) {modp=-modp;}
01592     
01593     double dsdz = pow((1.+pow(StateVector[2],2)+pow(StateVector[3],2)),0.5);
01594     double dxdz = invSqrt2*(StateVector[2]-StateVector[3]);
01595     double dydz = invSqrt2*(StateVector[2]+StateVector[3]);
01596     
01597     // Set up current muon details
01598     if(StateVector[4]>0.) charge = 1.;
01599     else if(StateVector[4]<0.) charge = -1.;
01600     
01601     TVector3 position(invSqrt2*(StateVector[0]-StateVector[1]),
01602                       invSqrt2*(StateVector[0]+StateVector[1]),
01603                       SlcStripData[Plane][0].csh->GetZPos());
01604 
01605     TVector3 momentum(modp*(dxdz/dsdz),
01606                       modp*(dydz/dsdz),
01607                       modp/dsdz);
01608 
01609     TVector3 bfield = bf->GetBField(position);
01610     bave += TMath::Sqrt(bfield[0]*bfield[0]+bfield[1]*bfield[1]+bfield[2]*bfield[2]);
01611     nbfield++;
01612 
01613     SwimParticle muon(position,momentum);
01614     muon.SetCharge(charge);
01615     SwimZCondition zc(SlcStripData[NewPlane][0].csh->GetZPos());
01616     // Do the swim, accounting for direction of motion w.r.t time too
01617     if( (GoForward==true && ZIncreasesWithTime==true)  || (GoForward==false && ZIncreasesWithTime==false) ) {
01618       if(UseGeoSwimmer) {
01619         done = GeoSwimmer::Instance()->SwimForward(muon,SlcStripData[NewPlane][0].csh->GetZPos());
01620       } else {
01621         done = myswimmer->SwimForward(muon,zc);
01622       }
01623     }
01624     else if( (GoForward==true && ZIncreasesWithTime==false)  || (GoForward==false && ZIncreasesWithTime==true) ) {
01625       if(UseGeoSwimmer) {
01626         done = GeoSwimmer::Instance()->SwimBackward(muon,SlcStripData[NewPlane][0].csh->GetZPos());
01627       } else {
01628         done = myswimmer->SwimBackward(muon,zc);
01629       }
01630     }
01631     if(done==true) {
01632       if(muon.GetDirection().Z()!=0. && muon.GetMomentumModulus()!=0.) {
01633         Output[0]=invSqrt2*(muon.GetPosition().Y()+muon.GetPosition().X());
01634         Output[1]=invSqrt2*(muon.GetPosition().Y()-muon.GetPosition().X());
01635         Output[2]=invSqrt2*((muon.GetDirection().Y()/muon.GetDirection().Z())+(muon.GetDirection().X()/muon.GetDirection().Z()));
01636         Output[3]=invSqrt2*((muon.GetDirection().Y()/muon.GetDirection().Z())-(muon.GetDirection().X()/muon.GetDirection().Z()));
01637         Output[4]=muon.GetCharge()/muon.GetMomentumModulus();
01638         // Get range and dS from the Swimmer
01639         if(dS) {*dS=muon.GetS();} if(Range) {*Range=muon.GetRange();} if(dE){*dE=muon.GetMomentumModulus()-momentum.Mag();} 
01640       }
01641       else {done=false;}
01642     }
01643     
01644   }
01645 
01646   else{
01647     // If infinite momentum, use straight line extrapolation
01648     double delz = (SlcStripData[NewPlane][0].csh->GetZPos()-SlcStripData[Plane][0].csh->GetZPos());
01649     Output[0]=StateVector[0] + StateVector[2]*delz;
01650     Output[1]=StateVector[1] + StateVector[3]*delz;
01651     Output[2]=StateVector[2];
01652     Output[3]=StateVector[3];
01653     Output[4]=StateVector[4];
01654 
01655     done=true;
01656   }    
01657   
01658   delete myswimmer;
01659   delete bf;
01660   return done;
01661 }

void AlgFitTrackCam::TimingFit ( CandFitTrackCamHandle cth  ) 

Definition at line 2643 of file AlgFitTrackCam.cxx.

References Chi2(), UgliStripHandle::ClearFiber(), StripStruct::csh, digit(), CandDigitHandle::GetCharge(), CandStripHandle::GetCharge(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), CandTrackHandle::GetdS(), PlexSEIdAltL::GetEnd(), UgliStripHandle::GetHalfLength(), CandStripHandle::GetPlaneView(), CandDigitHandle::GetPlexSEIdAltL(), CandStripHandle::GetStripEndId(), UgliGeomHandle::GetStripHandle(), CandTrackHandle::GetT(), CandStripHandle::GetTime(), CandTrackHandle::GetU(), CandTrackHandle::GetV(), Detector::kFar, Detector::kNear, StripEnd::kNegative, StripEnd::kPositive, StripEnd::kUnknown, MaxPlane, MinPlane, n, CandRecoHandle::SetEndT(), CandTrackHandle::SetNTimeFitDigit(), CandTrackHandle::SetTimeBackwardFitNDOF(), CandTrackHandle::SetTimeBackwardFitRMS(), CandTrackHandle::SetTimeFitChi2(), CandTrackHandle::SetTimeForwardFitNDOF(), CandTrackHandle::SetTimeForwardFitRMS(), CandRecoHandle::SetTimeOffset(), CandRecoHandle::SetTimeSlope(), CandRecoHandle::SetVtxT(), Particle::Sign(), size, StripListTime, vldc, UgliStripHandle::WlsPigtail(), and ZIncreasesWithTime.

Referenced by SetTrackProperties().

02644 {
02645   //  MSG("AlgFitTrackCam",Msg::kSynopsis) << "TimingFit" << endl;
02646 
02647   // Initialisations
02649   double s; double t; double q; double w; int n=0; 
02650   double Uncertainty=1.0; double MinUncertainty=1.0; 
02651   double MinCT=-3000.;
02652   
02653   // Time of first strip in track
02654   StripListTime=9.e30;
02655 
02656   // Create an offset such that dS=0 at the MinPlane
02657   double dSOffset=0.; double Sign=-1.; double dS[490];
02658   if(ZIncreasesWithTime==true) {dSOffset=cth.GetdS(MinPlane); Sign=1.;}
02659 
02660   // Store data needed in arrays. Charge is in PEs.
02661   double Qp[490];  double Qm[490];
02662   double CTp[490]; double CTm[490];
02663   int Skipp[490];  int Skipm[490];
02664   double C=3.e8;   
02665 
02666   double ErrorParam[3];
02667   ErrorParam[0]=1.; ErrorParam[1]=0.; ErrorParam[2]=0.;
02668 
02669   // Zero the arrays
02670   for(int i=0; i<490; ++i) { 
02671     dS[i]=0.; CTm[i]=0.; CTp[i]=0.; 
02672     Qp[i]=0.; Qm[i]=0.; 
02673     Skipp[i]=0; Skipm[i]=0;
02674   }
02675   // End of Initialisations
02677 
02678 
02679   // Copy track strips into vector
02681   vector<StripStruct> TimeFitStripData[490];
02682 
02683   TIter MyStripItr = cth.GetDaughterIterator();
02684 
02685   while(CandStripHandle* MyStrip = dynamic_cast<CandStripHandle*>(MyStripItr()))
02686     {    
02687       int MyPlane = MyStrip->GetPlane();
02688      
02689       StripStruct temp;
02690       temp.csh = MyStrip;
02691 
02692       if( cth.GetdS(MyPlane)>-0.99 ){               // distance must be greater
02693         TimeFitStripData[MyPlane].push_back(temp);  // than default value of -1
02694       }
02695       
02696     } 
02698 
02699 
02700   // Organise timing for the Far Detector
02702   if(vldc->GetDetector()==Detector::kFar) {  
02703 
02704     // Parameters for time fit residual vs PEs
02705     // (consistent with AlgTrackCam::WeightsForTimingFits() method)
02706     // ErrorParam[0]=0.56; ErrorParam[1]=0.50; ErrorParam[2]=-0.34; // (OLD)
02707     ErrorParam[0]=0.58; ErrorParam[1]=0.69; ErrorParam[2]=-0.33; // (NEW)
02708     MinUncertainty=ErrorParam[0];
02709 
02710     // Loop over all planes
02711     for(int i=MinPlane; i<=MaxPlane; ++i) {
02712         
02713       if(TimeFitStripData[i].size()>0) {
02714         dS[i]=Sign*(dSOffset-cth.GetdS(i));
02715         
02716         CTp[i]=C*cth.GetT(i,StripEnd::kPositive);
02717         CTm[i]=C*cth.GetT(i,StripEnd::kNegative);
02718         
02719         if(CTp[i]>MinCT && CTp[i]<StripListTime) {StripListTime=CTp[i];}
02720         if(CTm[i]>MinCT && CTm[i]<StripListTime) {StripListTime=CTm[i];}
02721         
02722         for(unsigned int j=0; j<TimeFitStripData[i].size(); ++j) {
02723           Qp[i]+=TimeFitStripData[i][j].csh->GetCharge(StripEnd::kPositive);
02724           Qm[i]+=TimeFitStripData[i][j].csh->GetCharge(StripEnd::kNegative);
02725         }
02726       }
02727     }
02728 
02729     // Subtract StripList time
02730     if(StripListTime<8.e30) {
02731       for(int i=MinPlane; i<=MaxPlane; ++i) {
02732         if(TimeFitStripData[i].size()>0) {
02733           CTp[i]-=StripListTime;
02734           CTm[i]-=StripListTime;
02735         }
02736       }
02737     }
02738     else {StripListTime=0.;}
02739 
02740   } // End Far Detector Section
02742 
02743 
02744 
02745   // Organise timing for the Near Detector
02747   if(vldc->GetDetector()==Detector::kNear) {
02748     // Parameters for PE vs time fit residual
02749     MinUncertainty=1.19;
02750     ErrorParam[0]=1.13; ErrorParam[1]=0.63; ErrorParam[2]=-0.21;
02751 
02752     double Index=1.77; double DistFromCentre=0.; double CTime=0.;  double FibreDist=0.;    
02753     double halfLength=0.; double DigChg=0.; double PlaneDigChg;
02754 
02755     StripEnd::StripEnd_t StpEnd = StripEnd::kUnknown;
02756     CandStripHandle* strip; CandDigitHandle* digit;
02757 
02758     UgliGeomHandle ugh = UgliGeomHandle(*vldc);
02759     UgliStripHandle striphandle;
02760 
02761 
02762     // Loop over all planes
02763     for(int i=MinPlane; i<=MaxPlane; ++i)
02764       {
02765         if(TimeFitStripData[i].size()>0) {dS[i]=Sign*(dSOffset-cth.GetdS(i));}
02766         PlaneDigChg=0.;
02767 
02768         // Loop over track strips on plane. Only +ve StripEnds should have charge.
02769         for(unsigned int j=0; j<TimeFitStripData[i].size(); ++j) {
02770           strip = TimeFitStripData[i][j].csh;
02771           CandDigitHandleItr digitItr(strip->GetDaughterIterator());
02772 
02773           Qp[i]+=strip->GetCharge(StripEnd::kPositive);
02774 
02775 
02776           // Loop over digits on strip. 
02778           while( (digit = digitItr()) ) { 
02779             StpEnd=digit->GetPlexSEIdAltL().GetEnd();
02780 
02781             if(StpEnd==StripEnd::kPositive) {
02782               FibreDist = 0.; DistFromCentre = 0.; CTime=0.; DigChg=0.;
02783               UgliStripHandle stripHandle = ugh.GetStripHandle(strip->GetStripEndId());
02784               halfLength = stripHandle.GetHalfLength(); 
02785          
02786               if(strip->GetPlaneView()==2) {DistFromCentre = cth.GetV(i);}
02787               if(strip->GetPlaneView()==3) {DistFromCentre = -cth.GetU(i);}
02788               
02789               FibreDist = (halfLength + DistFromCentre + stripHandle.ClearFiber(StpEnd) 
02790                            + stripHandle.WlsPigtail(StpEnd));
02791               
02792               CTime = C*(strip->GetTime(StpEnd)) - Index*FibreDist;
02793               DigChg=digit->GetCharge();
02794               
02795               PlaneDigChg+=DigChg; CTp[i]+=DigChg*CTime;
02796 
02797               if(CTime>MinCT && CTime<StripListTime) {StripListTime=CTime;}
02798             }
02799           }
02801         }
02802         if(PlaneDigChg>0.) CTp[i]/=PlaneDigChg;
02803       }
02804 
02805 
02806     // Subtract StripList time
02807     if(StripListTime<8.e30) {
02808       for(int i=MinPlane; i<=MaxPlane; ++i) {
02809         if(TimeFitStripData[i].size()>0) {
02810           CTp[i]-=StripListTime;
02811         }
02812       }
02813     }
02814     else {StripListTime=0.;}
02815     
02816   } // End near detector section
02818     
02819     
02820 
02821   // Carry out a simple straight line fit for T vs dS
02823   // Sqt: sum of charge*time, Sqss: sum of charge*dS*dS, etc.
02824   double Sqs=0; double Sqt=0; double Sqss=0; double Sqst=0; double Sqtt=0; double Sq=0;
02825   double TimeSlope=-999; double TimeOffset=-999; double RMS=-999;
02826   double CTCut = 0.; bool CalculateChi2=true;
02827 
02828   
02829   // On first iteration, carry out simple fit. Remove outlying points on subsequent passes.
02830   for(int itr=0; itr<3; ++itr) {
02831  
02832     for(int i=MinPlane; i<=MaxPlane; ++i) { 
02833 
02834       // Only consider planes where we found our final strips
02835       if(TimeFitStripData[i].size()>0) { 
02836         
02837         // For positive strip ends      
02838         s=dS[i]; q=Qp[i]; t=CTp[i];
02839 
02840         if (q<25) {Uncertainty = ErrorParam[0]+exp(ErrorParam[1]+ErrorParam[2]*q);}
02841         else {Uncertainty=MinUncertainty;}
02842         w = 1.0/pow(Uncertainty,2.0);
02843 
02844         if(q>0. && t>MinCT && Skipp[i]==0) {
02845           if(itr==0) {Sq+=w; Sqs+=w*s; Sqt+=w*t; Sqss+=w*s*s; Sqst+=w*s*t; Sqtt+=w*t*t; n++;}
02846           
02847           else if(fabs(t-TimeOffset-(s*TimeSlope)) > CTCut) {
02848             Sqs-=w*s; Sqt-=w*t; Sqss-=w*s*s; Sqst-=w*s*t; Sqtt-=w*t*t; Sq-=w; n--; Skipp[i]=1;
02849           }
02850         }
02851 
02852 
02853         // For negative strip ends
02854         s=dS[i]; q=Qm[i]; t=CTm[i];
02855 
02856         if (q<25) {Uncertainty = ErrorParam[0]+exp(ErrorParam[1]+ErrorParam[2]*q);}
02857         else {Uncertainty=MinUncertainty;}
02858         w = 1.0/pow(Uncertainty,2.0);
02859 
02860         if(q>0. && t>MinCT && Skipm[i]==0) {
02861           if(itr==0) {Sq+=w; Sqs+=w*s; Sqt+=w*t; Sqss+=w*s*s; Sqst+=w*s*t; Sqtt+=w*t*t; n++;}
02862           
02863           else if(fabs(t-TimeOffset-(s*TimeSlope)) > CTCut) {
02864             Sqs-=w*s; Sqt-=w*t; Sqss-=w*s*s; Sqst-=w*s*t; Sqtt-=w*t*t; Sq-=w; n--; Skipm[i]=1;
02865           }
02866         }
02867 
02868       }
02869     }
02870     
02871     // Calculate parameters
02872     if( (Sq*Sqss-Sqs*Sqs)!=0. && Sq!=0. ) {
02873       TimeSlope  = (Sq*Sqst-Sqs*Sqt)/(Sq*Sqss-Sqs*Sqs);
02874       TimeOffset = (Sqt*Sqss-Sqs*Sqst)/(Sq*Sqss-Sqs*Sqs);
02875       if( ((Sqtt/Sq)-((Sqt/Sq)*(Sqt/Sq)))>0. ) {
02876         RMS = pow((Sqtt/Sq)-((Sqt/Sq)*(Sqt/Sq)),0.5);
02877         CTCut = 3.+RMS;
02878       }
02879       else {CTCut = 3.5;}
02880     }
02881     else  {CalculateChi2=false; break;}
02882   }
02884 
02885 
02886 
02887   // Set timing properties for the fitted track
02889   if(n!=0 && CalculateChi2==true) {
02890 
02891     // Offset, slope and vtx/end times
02893     cth.SetTimeOffset((TimeOffset+StripListTime)/C);
02894     cth.SetTimeSlope(TimeSlope/C);
02895 
02896     if(ZIncreasesWithTime==true) {
02897       cth.SetVtxT((TimeOffset+StripListTime)/C); 
02898       cth.SetEndT((TimeOffset+StripListTime)/C+(dS[MaxPlane]*TimeSlope/C));
02899     }
02900     else {
02901       cth.SetEndT((TimeOffset+StripListTime)/C); 
02902       cth.SetVtxT((TimeOffset+StripListTime)/C+(dS[MaxPlane]*TimeSlope/C));
02903     }
02905 
02906 
02907     // Chi2
02909     double Residual2=0.; double Chi2=0.; 
02910  
02911     for(int i=MinPlane; i<=MaxPlane; ++i) { 
02912       
02913       if(TimeFitStripData[i].size()>0) { 
02914         // For positive strip ends
02915         s=dS[i]; q=Qp[i]; t=CTp[i];
02916         if(q>0. && t>MinCT && Skipp[i]==0) {
02917           Residual2=pow(t-TimeOffset-(s*TimeSlope),2);
02918 
02919           // From a rough parameterisation of uncertainty (in CT) vs number of PEs
02920           if (q<25) {Uncertainty = ErrorParam[0]+exp(ErrorParam[1]+ErrorParam[2]*q);}
02921           else {Uncertainty=MinUncertainty;}
02922 
02923           if(Uncertainty!=0.) {Chi2+=Residual2/pow(Uncertainty,2);}
02924         }
02925 
02926 
02927         // For negative strip ends
02928         q=Qm[i]; t=CTm[i];
02929         if(q>0. && t>MinCT && Skipm[i]==0) {
02930           Residual2=pow(t-TimeOffset-(s*TimeSlope),2);
02931 
02932           // From a rough parameterisation of uncertainty (in CT) vs number of PEs
02933           if (q<25) {Uncertainty = ErrorParam[0]+exp(ErrorParam[1]+ErrorParam[2]*q);}
02934           else {Uncertainty=MinUncertainty;}
02935           
02936           if(Uncertainty!=0.) {Chi2+=Residual2/pow(Uncertainty,2);}
02937         }
02938       }
02939       
02940     }
02941     // Set these properties
02942     cth.SetTimeFitChi2(Chi2);
02943     cth.SetNTimeFitDigit(n);
02944   }
02946 
02947 
02948 
02949   // Now carry out fits with gradients constrained to be +/- c
02951   double CTIntercept[2]; double Csigma[2]; double Ctrunc[2]; 
02952   double ChiSqPositive=-999; double ChiSqNegative=-999; 
02953   int ChiSqNdfPos=-999; int ChiSqNdfNeg=-999;  
02954   double Swtt[2]; double Swt[2]; double Sw[2]; int npts[2]={0,0};
02955 
02956   if(Sq!=0.) {
02957     CTIntercept[0]=Sqt/Sq; Csigma[0]=-99999.9; Ctrunc[0]=-99999.9;
02958     CTIntercept[1]=Sqt/Sq; Csigma[1]=-99999.9; Ctrunc[1]=-99999.9;
02959     
02960     for(int itr=0; itr<2; ++itr)
02961       {
02962         Swtt[0]=0.; Swt[0]=0.; Sw[0]=0.; npts[0]=0;
02963         Swtt[1]=0.; Swt[1]=0.; Sw[1]=0.; npts[1]=0;
02964         
02965         for(int i=0; i<490; ++i)
02966           {
02967             // For positive strip ends
02968             if(Qp[i]>0. && CTp[i]>MinCT) {
02969               q=Qp[i]; 
02970 
02971               if (q<25) {Uncertainty = ErrorParam[0]+exp(ErrorParam[1]+ErrorParam[2]*q);}
02972               else {Uncertainty=MinUncertainty;}
02973               w = 1.0/pow(Uncertainty,2.0);
02974 
02975               t=CTp[i]-dS[i]+CTIntercept[0];
02976               if(Ctrunc[0]<0. || fabs(t)<Ctrunc[0]) {Swtt[0]+=w*t*t; Swt[0]+=w*t; Sw[0]+=w; ++npts[0];}
02977               
02978               t=CTp[i]+dS[i]+CTIntercept[1];
02979               if(Ctrunc[1]<0. || fabs(t)<Ctrunc[1]) {Swtt[1]+=w*t*t; Swt[1]+=w*t; Sw[1]+=w; ++npts[1];}
02980             }
02981             
02982             // For negative strip ends
02983             if(Qm[i]>0. && CTm[i]>MinCT) {
02984               q=Qm[i]; 
02985 
02986               if (q<25) {Uncertainty = ErrorParam[0]+exp(ErrorParam[1]+ErrorParam[2]*q);}
02987               else {Uncertainty=MinUncertainty;}
02988               w = 1.0/pow(Uncertainty,2.0);
02989 
02990               t=CTm[i]-dS[i]+CTIntercept[0];
02991               if(Ctrunc[0]<0. || fabs(t)<Ctrunc[0]) {Swtt[0]+=w*t*t; Swt[0]+=w*t; Sw[0]+=w; ++npts[0];}
02992               
02993               t=CTm[i]+dS[i]+CTIntercept[1];
02994               if(Ctrunc[1]<0. || fabs(t)<Ctrunc[1]) {Swtt[1]+=w*t*t; Swt[1]+=w*t; Sw[1]+=w; ++npts[1];}
02995             }
02996           }
02997         
02998         // Results for fit with gradient +C
02999         if(npts[0]>1 && Sw[0]!=0.) {
03000           CTIntercept[0]=CTIntercept[0]-Swt[0]/Sw[0]; Csigma[0]=0.; 
03001           if((Swtt[0]/Sw[0])-(Swt[0]/Sw[0])*(Swt[0]/Sw[0])>0.) {Csigma[0]=pow((Swtt[0]/Sw[0])-(Swt[0]/Sw[0])*(Swt[0]/Sw[0]),0.5);}
03002           ChiSqPositive=Csigma[0]; ChiSqNdfPos=npts[0]-1; 
03003           Ctrunc[0]=Csigma[0]+3.;
03004         }
03005         
03006         // Results for fit with gradient -C
03007         if(npts[1]>1 && Sw[1]!=0.) {
03008           CTIntercept[1]=CTIntercept[1]-Swt[1]/Sw[1]; Csigma[1]=0.;
03009           if((Swtt[1]/Sw[1])-(Swt[1]/Sw[1])*(Swt[1]/Sw[1])>0.) {Csigma[1]=pow((Swtt[1]/Sw[1])-(Swt[1]/Sw[1])*(Swt[1]/Sw[1]),0.5);}
03010           ChiSqNegative=Csigma[1]; ChiSqNdfNeg=npts[1]-1; 
03011           Ctrunc[1]=Csigma[1]+3.;
03012         }
03013         
03014       }
03015   }
03016   // Set these properties
03017   cth.SetTimeForwardFitRMS(ChiSqPositive);
03018   cth.SetTimeForwardFitNDOF(ChiSqNdfPos);
03019   cth.SetTimeBackwardFitRMS(ChiSqNegative);
03020   cth.SetTimeBackwardFitNDOF(ChiSqNdfNeg);
03022    
03023 }

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

Reimplemented from AlgBase.

Definition at line 3955 of file AlgFitTrackCam.cxx.

03956 {
03957 }

void AlgFitTrackCam::UpdateCovMatrix (  ) 

Definition at line 2307 of file AlgFitTrackCam.cxx.

References C_k, C_k_intermediate, debug, H_k, Identity, K_k, Msg::kVerbose, Munits::m, and MSG.

Referenced by GoBackwards(), GoForwards(), ShowerSwim(), and SpectrometerSwim().

02308 {
02309   // C_k = (Identity - (K_k * H_k) ) * C_k_intermediate
02310   //  MSG("AlgFitTrackCam",Msg::kDebug) << "UpdateCovMatrix" << endl;
02311 
02312   for (int i=0; i<5; ++i) {
02313     for (int j=0; j<5; ++j) {  
02314       C_k[i][j]=0;
02315       for (int m=0; m<5; ++m) {   
02316         C_k[i][j]+=(Identity[i][m] - K_k[i]*H_k[m]) * C_k_intermediate[m][j];
02317       }
02318     }
02319   }
02320 
02321 
02322   // Diagonal elements should be positive
02323   double covlim = 1.e-8;
02324 
02325   for(int i=0; i<5; ++i) {
02326     if(C_k[i][i]<covlim) {
02327       MSG("AlgFitTrackCam",Msg::kVerbose) << "Negative diagonal element in C_k" << endl;
02328       C_k[i][i]=covlim;
02329     }
02330   }
02331 
02332 
02333   // Display
02334   if(debug) {
02335     cout << "Filtered Covariance matrix" << endl;
02336     for(int i=0; i<5; ++i) {
02337       for(int j=0; j<5; ++j) {
02338         cout << C_k[i][j] << " ";
02339       }
02340       cout << endl;  
02341     }
02342   }
02343 
02344 }

void AlgFitTrackCam::UpdateStateVector ( const int  Plane,
const int  NewPlane,
const bool  GoForward 
)

Definition at line 2141 of file AlgFitTrackCam.cxx.

References CheckValues(), EndofRangePlane, K_k, Msg::kDebug, Msg::kVerbose, LastIteration, MajorityCurvature(), MSG, PassTrack, SlcStripData, Swim(), TotalNSwimFail, TrkStripData, x_k, x_k_minus, and ZIncreasesWithTime.

Referenced by GoBackwards(), GoForwards(), ShowerSwim(), and SpectrometerSwim().

02142 {  
02143   // x_k = (F_k_minus * x_k_minus) + K_k(m_k - (H_k * F_k_minus * x_k_minus) )
02144   //  MSG("AlgFitTrackCam",Msg::kDebug) << "UpdateStateVector" << endl;
02145 
02146 
02147   double HFx=0;
02148   double m_k=0;
02149   double MeasurementFactor=0;
02150   int nswimfail=0;
02151 
02152 
02153   // Calculate F_k_minus * x_k_minus, using the Swimmer
02154   // Also get an accurate measure of dS and Range from the Swimmer
02156   double StateVector[5];
02157   double Prediction[5];
02158   bool GetPrediction=false;
02159 
02160   for(int i=0; i<5; ++i) {StateVector[i]=x_k_minus[i];}
02161 
02162   // Prediction could be used without GeoSwimmer calculation. Prediction is initialized with linear extrapolation.
02163   for(int i=0; i<5; i++) {
02164     double delz = (SlcStripData[NewPlane][0].csh->GetZPos()-SlcStripData[Plane][0].csh->GetZPos());
02165     Prediction[0]=StateVector[0] + StateVector[2]*delz;
02166     Prediction[1]=StateVector[1] + StateVector[3]*delz;
02167     Prediction[2]=StateVector[2];
02168     Prediction[3]=StateVector[3];
02169     Prediction[4]=StateVector[4];
02170   }
02171 
02172   // Do the swim
02174   while(GetPrediction==false && nswimfail<=10) {
02175     MSG("AlgFitTrackCam",Msg::kVerbose) << " state " << StateVector[0] << " " 
02176                                         << StateVector[1] << " " << StateVector[2] << " " 
02177                                         << StateVector[3] << " " << StateVector[4] << endl;   
02178     
02179     GetPrediction=Swim(StateVector, Prediction, Plane, NewPlane, GoForward);
02180     
02181     MSG("AlgFitTrackCam",Msg::kVerbose) << " predict state " << Prediction[0] << " " 
02182                                         << Prediction[1] << " " << Prediction[2] << " " 
02183                                         << Prediction[3] << " " << Prediction[4] << endl;
02184 
02185     if(GetPrediction==false) {
02186       StateVector[4]*=0.5; 
02187       nswimfail++; TotalNSwimFail++;
02188       MSG("AlgFitTrackCam",Msg::kVerbose) << "UpdateStateVector, Prediction failed - Double momentum and swim again" << endl;
02189     }
02190   }
02191 
02192   if(nswimfail>10) {  // Swim shouldn't fail, as it succeeded to get the propagator
02193     MSG("AlgFitTrackCam",Msg::kDebug) << "UpdateStateVector, nswimfail>10, fail track" << endl;
02194     PassTrack=false;
02195   }
02197 
02198   
02200 
02201 
02202   //  MSG("AlgFitTrackCam",Msg::kDebug) << "UpdateStateVector, Check predicted state " << endl;
02203   CheckValues(Prediction, NewPlane);
02204 
02205 
02206   // Calculate H_k * F_k_minus * x_k_minus
02208   if(TrkStripData[NewPlane][0].PlaneView==2) {HFx=Prediction[0];}
02209   if(TrkStripData[NewPlane][0].PlaneView==3) {HFx=Prediction[1];}
02210 
02211   MSG("AlgFitTrackCam",Msg::kVerbose) << "HFx " << HFx << endl;
02213 
02214   
02215   // Read in measurement
02217   m_k=TrkStripData[NewPlane][0].TPos;
02218   MSG("AlgFitTrackCam",Msg::kVerbose) << "m_k " << TrkStripData[NewPlane][0].TPos << endl;
02219   
02220   MeasurementFactor=(m_k-HFx);
02222 
02223   
02224   // Calculate x_k
02226   for (int i=0; i<5; ++i) {
02227     x_k[i]=0.;
02228     x_k[i]+=Prediction[i]+(K_k[i]*MeasurementFactor);
02229   }
02231 
02232 
02233   //  MSG("AlgFitTrackCam",Msg::kDebug) << "UpdateStateVector, Check filtered state " << endl;
02234   CheckValues(x_k, NewPlane);
02235 
02236 
02237   // Care with multiple range corrections - do not want to flip sign
02238   // (multiple corrections mean sign changes can occur even though absolute value stays same)
02239   
02240   double Maxqp = 4.;
02241   if(fabs(x_k_minus[4])==Maxqp && 
02242      ( (GoForward==true && ZIncreasesWithTime==true) 
02243        || (GoForward==false && ZIncreasesWithTime==false) ) ) 
02244     {
02245       if(!LastIteration) x_k[4] = (x_k_minus[4]>0 ? Maxqp : -Maxqp);
02246     }
02247 
02248   //if on last plane, disregard a flip in the charge sign
02249   //if(x_k_minus[4]!=0){
02250   //  if( x_k[4]/x_k_minus[4]<0 && 
02251   //      ( (GoForward==true  && ZIncreasesWithTime==true  && NewPlane >= EndofRangePlane ) 
02252   //     || (GoForward==false && ZIncreasesWithTime==false && NewPlane <= EndofRangePlane ) ) ) 
02253   //   {
02254   //     MSG("AlgFitTrackCam",Msg::kSynopsis) << " Disregard charge flip in final plane: NewPlane=" << NewPlane << " " << x_k[4] << "-> " << -x_k[4] << endl;
02255   //      x_k[4] = -x_k[4];
02256   //    }
02257   // }  
02258 
02259   // if on last plane, disregard a charge sign that conflicts with the majority curvature of the track (set cut at 2/3 of track)
02260   if( (GoForward==true && ZIncreasesWithTime==true && NewPlane >= EndofRangePlane ) 
02261    || (GoForward==false && ZIncreasesWithTime==false && NewPlane <= EndofRangePlane ) )
02262     {
02263       double majority_curvature = this->MajorityCurvature();
02264       MSG("AlgFitTrackCam",Msg::kDebug) << " Checking Last Plane: Q/P=" << x_k[4] << " MajorityCurvature = " << majority_curvature << endl;
02265 
02266       if( ( majority_curvature<-0.33 && x_k[4]>0.0 )
02267        || ( majority_curvature>+0.33 && x_k[4]<0.0 ) ){
02268         MSG("AlgFitTrackCam",Msg::kDebug) << " Disregard charge sign in last plane: NewPlane=" << NewPlane << " " << x_k[4] << "-> " << -x_k[4] << endl;
02269         x_k[4] = -x_k[4];
02270       }
02271     }
02272   
02273   // Display
02274   MSG("AlgFitTrackCam",Msg::kVerbose) << "Filtered State vector: "
02275                                       << x_k[0] << " " << x_k[1] << " " << x_k[2] << " "
02276                                       << x_k[3] << " " << x_k[4] << endl;
02277 }


Member Data Documentation

Definition at line 179 of file AlgFitTrackCam.h.

Double_t AlgFitTrackCam::bave [private]

Definition at line 121 of file AlgFitTrackCam.h.

Referenced by RunAlg(), SetTrackProperties(), and Swim().

Definition at line 174 of file AlgFitTrackCam.h.

double AlgFitTrackCam::C_k[5][5] [private]
double AlgFitTrackCam::C_k_intermediate[5][5] [private]

Definition at line 132 of file AlgFitTrackCam.h.

Referenced by CalcKalmanGain(), ExtrapCovMatrix(), RunAlg(), and UpdateCovMatrix().

double AlgFitTrackCam::C_k_minus[5][5] [private]

Definition at line 131 of file AlgFitTrackCam.h.

Referenced by ExtrapCovMatrix(), GetInitialCovarianceMatrix(), MoveArrays(), and RunAlg().

bool AlgFitTrackCam::debug [private]
double AlgFitTrackCam::DeltaPlane [private]
double AlgFitTrackCam::DeltaZ [private]
double AlgFitTrackCam::EndCov[5] [private]

Definition at line 143 of file AlgFitTrackCam.h.

Referenced by GoBackwards(), GoForwards(), RunAlg(), and SetTrackProperties().

Bool_t AlgFitTrackCam::EndofRange [private]

Definition at line 122 of file AlgFitTrackCam.h.

Referenced by CheckValues(), GoBackwards(), GoForwards(), and RunTheFitter().

double AlgFitTrackCam::eqp_biased [private]

Definition at line 126 of file AlgFitTrackCam.h.

Referenced by RunTheFitter(), and SetTrackProperties().

double AlgFitTrackCam::F_k[5][5] [private]

Definition at line 133 of file AlgFitTrackCam.h.

Referenced by RunAlg().

double AlgFitTrackCam::F_k_minus[5][5] [private]

Definition at line 134 of file AlgFitTrackCam.h.

Referenced by ExtrapCovMatrix(), GetCombiPropagator(), and RunAlg().

Definition at line 169 of file AlgFitTrackCam.h.

Referenced by IsInDetector(), NDPlaneIsActive(), and SetRangeAnddS().

int AlgFitTrackCam::H_k[5] [private]
int AlgFitTrackCam::Identity[5][5] [private]

Definition at line 140 of file AlgFitTrackCam.h.

Referenced by RunAlg(), and UpdateCovMatrix().

double AlgFitTrackCam::InitMaxqp [private]

Definition at line 175 of file AlgFitTrackCam.h.

Definition at line 176 of file AlgFitTrackCam.h.

Definition at line 177 of file AlgFitTrackCam.h.

double AlgFitTrackCam::K_k[5] [private]

Definition at line 137 of file AlgFitTrackCam.h.

Referenced by CalcKalmanGain(), RunAlg(), UpdateCovMatrix(), and UpdateStateVector().

Definition at line 124 of file AlgFitTrackCam.h.

Referenced by GoBackwards(), GoForwards(), RunTheFitter(), and UpdateStateVector().

Definition at line 178 of file AlgFitTrackCam.h.

int AlgFitTrackCam::MaxPlane [private]

Definition at line 168 of file AlgFitTrackCam.h.

Referenced by InitialFramework(), and SpectrometerSwim().

int AlgFitTrackCam::MinPlane [private]
Int_t AlgFitTrackCam::nbfield [private]

Definition at line 120 of file AlgFitTrackCam.h.

Referenced by RunAlg(), SetTrackProperties(), and Swim().

int AlgFitTrackCam::NIter [private]

Definition at line 164 of file AlgFitTrackCam.h.

Referenced by GoBackwards(), GoForwards(), RunAlg(), RunTheFitter(), and SetTrackProperties().

Definition at line 167 of file AlgFitTrackCam.h.

Referenced by GenerateNDSpectStrips(), InitialFramework(), and RunAlg().

bool AlgFitTrackCam::PassTrack [private]
double AlgFitTrackCam::Q_k[5][5] [private]

Definition at line 135 of file AlgFitTrackCam.h.

Referenced by RunAlg().

double AlgFitTrackCam::Q_k_minus[5][5] [private]

Definition at line 136 of file AlgFitTrackCam.h.

Referenced by ExtrapCovMatrix(), GetNoiseMatrix(), and RunAlg().

bool AlgFitTrackCam::SaveData [private]

Definition at line 159 of file AlgFitTrackCam.h.

Referenced by GoBackwards(), GoForwards(), RunAlg(), and RunTheFitter().

Definition at line 162 of file AlgFitTrackCam.h.

Referenced by RemoveTrkHitsInShw(), RunAlg(), SetTrackProperties(), and ShowerStrips().

vector<StripStruct> AlgFitTrackCam::SlcStripData[490] [private]

Definition at line 171 of file AlgFitTrackCam.h.

Referenced by TimingFit().

Definition at line 160 of file AlgFitTrackCam.h.

Referenced by RemoveTrkHitsInShw(), RunAlg(), RunTheFitter(), ShowerStrips(), and ShowerSwim().

Definition at line 104 of file AlgFitTrackCam.h.

Referenced by IsStripInAnotherTrack(), RunAlg(), and AlgFitTrackCamList::RunAlg().

Definition at line 127 of file AlgFitTrackCam.h.

Referenced by RunAlg(), and Swim().

double AlgFitTrackCam::VtxCov[5] [private]

Definition at line 142 of file AlgFitTrackCam.h.

Referenced by GoBackwards(), GoForwards(), RunAlg(), RunTheFitter(), and SetTrackProperties().

double AlgFitTrackCam::x_k[5] [private]
double AlgFitTrackCam::x_k4_biased [private]

Definition at line 125 of file AlgFitTrackCam.h.

Referenced by RunAlg(), RunTheFitter(), and SetTrackProperties().

double AlgFitTrackCam::x_k_minus[5] [private]

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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1