← 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/MAUSTrackWrapper.cc

MAUS-v2.3.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/common_cpp/Recon/Kalman/MAUSTrackWrapper.cc
95
95
      double P = patrec_momentum.mag();
96
96
      double patrec_bias; // Account for two planes of energy loss
97
97
      if (tracker == 0) {
98
 
        patrec_bias = (P + 1.4) / P;
 
98
        patrec_bias = (P + 1.3) / P;
99
99
      } else {
100
 
        patrec_bias = (P - 1.4) / P;
 
100
        patrec_bias = (P - 1.3) / P;
101
101
      }
102
102
      patrec_momentum = patrec_bias * patrec_momentum;
103
103
    }
287
287
      int ID = the_track[i].GetId();
288
288
      double position = the_track[i].GetPosition();
289
289
      const Kalman::State& smoothed_state = the_track[i].GetSmoothed();
 
290
//      const Kalman::State& smoothed_state = the_track[i].GetFiltered();
 
291
//      const Kalman::State& smoothed_state = the_track[i].GetPredicted();
290
292
      const Kalman::State& data_state = the_track[i].GetData();
291
293
 
292
294
      SciFiTrackPoint* new_point = new SciFiTrackPoint();
300
302
 
301
303
      ThreeVector pos;
302
304
      ThreeVector mom;
 
305
      ThreeVector grad;
 
306
      ThreeVector error_pos;
 
307
      ThreeVector error_mom;
 
308
      ThreeVector error_grad;
303
309
 
304
310
      TMatrixD state_vector = smoothed_state.GetVector();
 
311
      TMatrixD state_covariance = smoothed_state.GetCovariance();
305
312
 
306
313
      if (dimension == 4) {
307
314
        pos.setZ(position);
308
315
        pos.setX(state_vector(0, 0));
 
316
        pos.setY(state_vector(2, 0));
 
317
 
309
318
        mom.setX(state_vector(1, 0)*default_mom);
310
 
        pos.setY(state_vector(2, 0));
311
319
        mom.setY(state_vector(3, 0)*default_mom);
312
320
 
313
321
        pos *= reference_rot;
314
322
        pos += reference_pos;
315
 
 
316
323
        mom *= reference_rot;
 
324
 
317
325
        if (tracker == 0) mom *= -1.0; // Direction of recon is reveresed upstream.
318
326
        mom.setZ(default_mom); // MeV/c
 
327
 
 
328
        grad.SetX(mom.x()/mom.z());
 
329
        grad.SetY(mom.y()/mom.z());
 
330
        grad.SetZ(0.0);
 
331
 
 
332
        error_pos.SetX(std::sqrt(state_covariance(0, 0)));
 
333
        error_pos.SetY(std::sqrt(state_covariance(2, 2)));
 
334
        error_pos.SetZ(0.0);
 
335
 
 
336
        error_mom.SetX(std::sqrt(default_mom*state_covariance(1, 1)));
 
337
        error_mom.SetY(std::sqrt(default_mom*state_covariance(3, 3)));
 
338
        error_mom.SetZ(200.0);
 
339
 
 
340
        error_grad.SetX(std::sqrt(state_covariance(1, 1)));
 
341
        error_grad.SetY(std::sqrt(state_covariance(3, 3)));
 
342
        error_grad.SetZ(0.0);
319
343
      } else if (dimension == 5) {
320
344
        pos.setX(state_vector(0, 0));
 
345
        pos.setY(state_vector(2, 0));
 
346
        pos.setZ(position);
 
347
 
321
348
        mom.setX(state_vector(1, 0));
322
 
        pos.setY(state_vector(2, 0));
323
349
        mom.setY(state_vector(3, 0));
324
 
        pos.setZ(position);
325
350
 
326
351
        pos *= reference_rot;
327
352
        pos += reference_pos;
333
358
          charge = 1;
334
359
        }
335
360
        mom.setZ(fabs(1.0/state_vector(4, 0)));
 
361
 
 
362
        grad.SetX(mom.x() / mom.z());
 
363
        grad.SetY(mom.y() / mom.z());
 
364
        grad.SetZ(0.0);
 
365
 
 
366
        error_pos.SetX(std::sqrt(state_covariance(0, 0)));
 
367
        error_pos.SetY(std::sqrt(state_covariance(2, 2)));
 
368
        error_pos.SetZ(0.0);
 
369
 
 
370
        error_mom.SetX(std::sqrt(state_covariance(1, 1)));
 
371
        error_mom.SetY(std::sqrt(state_covariance(3, 3)));
 
372
        error_mom.SetZ(std::sqrt(state_covariance(4, 4)) / (state_vector(4, 0)*state_vector(4, 0)));
 
373
 
 
374
        error_grad.SetX(grad.x()*std::sqrt(((error_mom.x()*error_mom.x()) / (mom.x()*mom.x()))
 
375
                                              +(error_mom.z()*error_mom.z() / (mom.z()*mom.z()))));
 
376
        error_grad.SetY(grad.y()*std::sqrt(((error_mom.y()*error_mom.y()) / (mom.y()*mom.y()))
 
377
                                             + (error_mom.z()*error_mom.z() / (mom.z()*mom.z()))));
 
378
        error_grad.SetZ(0.0);
336
379
      }
337
380
 
338
381
      new_point->set_pos(pos);
339
382
      new_point->set_mom(mom);
 
383
      new_point->set_gradient(grad);
 
384
      new_point->set_pos_error(error_pos);
 
385
      new_point->set_mom_error(error_mom);
 
386
      new_point->set_gradient_error(error_grad);
340
387
 
341
388
      // TODO
342
389
      // CHARGE!
343
390
      if (data_state) {
344
 
        new_point->set_pull(fitter->CalculatePull(i).GetVector()(0, 0));
 
391
        Kalman::State pull = fitter->CalculatePull(i);
 
392
        double weighted_pull = pull.GetVector()(0, 0) / std::sqrt(pull.GetCovariance()(0, 0));
 
393
        new_point->set_pull(weighted_pull);
345
394
        new_point->set_residual(fitter->CalculateFilteredResidual(i).GetVector()(0, 0));
346
395
        new_point->set_smoothed_residual(fitter->CalculateSmoothedResidual(i).GetVector()(0, 0));
347
396
      } else {