← Back to branch summary

~christopher-hunt08/maus/maus_integrated_kalman

« back to all changes in this revision

Viewing changes to src/common_cpp/Recon/Kalman/KalmanTrackFit.cc

  • Committer: Durga Rajaram
  • Date: 2014-10-06 20:45:37 UTC
  • mfrom: (697.24.7 trunk)
  • mto: (697.30.10 maus)
  • mto: This revision was merged to the branch mainline in revision 703.
  • Revision ID: durga@fnal.gov-20141006204537-1aokre71ss0b2vcv
merge from trunk; set part_event_number for extra emr events

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/common_cpp/Recon/Kalman/KalmanTrackFit.cc
62
62
 
63
63
    // Filter the first state.
64
64
    _filter->Process(sites.front());
 
65
    // int first_site_id = abs(sites.front()->id());
65
66
 
66
67
    // Run the extrapolation & filter chain.
67
 
    size_t numb_measurements = sites.size();
68
 
    for ( size_t j = 1; j < numb_measurements; ++j ) {
 
68
    size_t numb_sites = sites.size();
 
69
    // Loop over all 15 planes
 
70
    for ( size_t j = 1; j < sites.size(); ++j ) {
69
71
      // Predict the state vector at site i...
70
72
      _propagator->Extrapolate(sites, j);
 
73
 
71
74
      // ... Filter...
72
75
      _filter->Process(sites.at(j));
73
76
    }
74
 
    _propagator->PrepareForSmoothing(sites);
 
77
    // std::cerr << "Prepare For Smoothing... " << std::endl;
 
78
    _propagator->PrepareForSmoothing(sites.back());
75
79
    // ...and Smooth back all sites.
76
 
    for ( int k = static_cast<int> (numb_measurements-2); k > -1; --k ) {
 
80
    for ( int k = static_cast<int> (sites.size()-2); k >= 0; --k ) {
77
81
      _propagator->Smooth(sites, k);
78
82
      _filter->UpdateH(sites.at(k));
79
 
      _filter->SetResidual(sites.at(k), KalmanState::Smoothed);
 
83
      _filter->ComputeResidual(sites.at(k), KalmanState::Smoothed);
80
84
    }
81
85
 
82
86
    track->set_tracker(seed->tracker());
97
101
}
98
102
 
99
103
void KalmanTrackFit::ComputeChi2(SciFiTrack *track, KalmanStatesPArray sites) {
100
 
  double f_chi2 = 0.;
101
 
  double s_chi2 = 0.;
102
 
  // Find the ndf for this track: numb measurements - numb parameters to be estimated
 
104
  // Prepare to loop over all Kalman sites...
103
105
  size_t n_sites = sites.size();
104
 
  // Find n_parameters by looking at the dimension of the state vector.
105
 
  int n_parameters = sites.at(0)->a(KalmanState::Filtered).GetNrows();
106
 
 
107
 
  int ndf = n_sites - n_parameters;
108
 
 
 
106
  // ... summing chi2...
 
107
  double chi2 = 0.;
 
108
  // .. and counting the numb. of measurements...
 
109
  int n_measurements = 0;
109
110
  for ( size_t i = 0; i < n_sites; ++i ) {
110
111
    KalmanState *site = sites.at(i);
111
 
    f_chi2 += site->chi2(KalmanState::Filtered);
112
 
    s_chi2 += site->chi2(KalmanState::Smoothed);
 
112
    if  ( site->contains_measurement() ) {
 
113
      n_measurements++;
 
114
      chi2 += site->chi2();
 
115
    }
113
116
  }
114
 
  double P_value = TMath::Prob(f_chi2, ndf);
115
 
  track->set_f_chi2(f_chi2);
116
 
  track->set_s_chi2(s_chi2);
 
117
  // Find the ndf for this track: numb measurements - numb parameters to be estimated
 
118
  int n_parameters = sites.at(0)->a(KalmanState::Filtered).GetNrows();
 
119
  int ndf = n_measurements - n_parameters;
 
120
  track->set_chi2(chi2);
117
121
  track->set_ndf(ndf);
 
122
  double P_value = TMath::Prob(chi2, ndf);
118
123
  track->set_P_value(P_value);
119
124
}
120
125
 
123
128
  int tracker = track->tracker();
124
129
  if ( pvalue != pvalue ) return;
125
130
  for ( size_t i = 0; i < sites.size(); ++i ) {
126
 
    sites.at(i)->MoveToGlobalFrame(_RefPos[tracker]);
127
 
    SciFiTrackPoint *track_point = new SciFiTrackPoint(sites.at(i));
128
 
    track->add_scifitrackpoint(track_point);
 
131
    if ( sites.at(i)->contains_measurement() ) {
 
132
      sites.at(i)->MoveToGlobalFrame(_RefPos[tracker]);
 
133
      SciFiTrackPoint *track_point = new SciFiTrackPoint(sites.at(i));
 
134
      track->add_scifitrackpoint(track_point);
 
135
    }
129
136
  }
130
137
  event.add_scifitrack(track);
131
138
}
135
142
 
136
143
  for ( size_t i = 0; i < numb_sites; ++i ) {
137
144
    KalmanState* site = sites.at(i);
138
 
    Squeak::mout(Squeak::info)
 
145
    // Squeak::mout(Squeak::info)
 
146
    std::cerr
139
147
    << "=========================================="  << "\n"
140
148
    << "SITE ID: " << site->id() << "\n"
141
149
    << "SITE Z: " << site->z()   << "\n"
144
152
                      << (site->a(KalmanState::Projected))(1, 0) << " "
145
153
                      << (site->a(KalmanState::Projected))(2, 0) << " "
146
154
                      << (site->a(KalmanState::Projected))(3, 0) << " "
147
 
                      << (site->a(KalmanState::Projected))(4, 0) << "\n"
148
155
    << "Filtered: " << (site->a(KalmanState::Filtered))(0, 0) << " "
149
156
                      << (site->a(KalmanState::Filtered))(1, 0) << " "
150
157
                      << (site->a(KalmanState::Filtered))(2, 0) << " "
151
 
                      << (site->a(KalmanState::Filtered))(3, 0) << " "
152
 
                      << (site->a(KalmanState::Filtered))(4, 0) << "\n"
 
158
                      << (site->a(KalmanState::Filtered))(3, 0) << "\n"
153
159
    << "================Residuals================"   << "\n"
154
160
    << (site->residual(KalmanState::Projected))(0, 0)  << "\n"
155
161
    << (site->residual(KalmanState::Filtered))(0, 0)   << "\n"
156
162
    << (site->residual(KalmanState::Smoothed))(0, 0)   << "\n"
157
163
    << "=========================================="
158
164
    << std::endl;
 
165
    // site->covariance_matrix(KalmanState::Projected).Print();
 
166
    // site->covariance_matrix(KalmanState::Filtered).Print();
 
167
    // site->covariance_matrix(KalmanState::Smoothed).Print();
159
168
  }
160
169
}
161
170