← Back to branch summary

~mfedorov/maus/trunk

« back to all changes in this revision

Viewing changes to src/common_cpp/Recon/Global/TrackMatching.cc

  • Committer: Chris Rogers
  • Date: 2017-09-12 12:06:30 UTC
  • Revision ID: chris.rogers@stfc.ac.uk-20170912120630-d61q8c16tnqo1e0v
Propagation/virtual output from global tracks

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/common_cpp/Recon/Global/TrackMatching.cc
38
38
    std::string pid_hypothesis_string, int beamline_polarity,
39
39
    std::map<std::string, std::pair<double, double> > matching_tolerances,
40
40
    double max_step_size, std::pair<bool, bool> no_check_settings,
41
 
    bool energy_loss, bool residuals, TrackMatching::geometry_algorithm algo) {
 
41
    bool energy_loss, bool residuals, TrackMatching::geometry_algorithm algo,
 
42
    std::vector<double> extra_z_planes) {
42
43
  _global_event = global_event;
43
44
  _mapper_name = mapper_name;
44
45
  _pid_hypothesis_string = pid_hypothesis_string;
49
50
  _no_check_settings = no_check_settings;
50
51
  _residuals = residuals;
51
52
  _geo_algo = algo;
 
53
  _extra_z_planes = extra_z_planes;
52
54
}
53
55
 
54
56
void TrackMatching::USTrack() {
296
298
  delete emr_track_array;
297
299
}
298
300
 
 
301
//void TrackMatching::throughMatchPropagate() {}
 
302
 
 
303
DataStructure::Global::Track* TrackMatching::throughMatchTOF(
 
304
                            DataStructure::Global::Track* us_track,
 
305
                            DataStructure::Global::Track* ds_track,
 
306
                            DataStructure::Global::PID pid) {
 
307
    ofstream throughfile;
 
308
    if (_residuals) {
 
309
      throughfile.open("match_through.csv", std::ios::out | std::ios::app);
 
310
    }
 
311
    if (!us_track or !us_track->HasDetector(DataStructure::Global::kTOF1) or
 
312
        !ds_track or !ds_track->HasDetector(DataStructure::Global::kTOF2)) {
 
313
      return NULL;
 
314
    }
 
315
    DataStructure::Global::TrackPointCPArray us_trackpoints = us_track->GetTrackPoints();
 
316
    DataStructure::Global::TrackPointCPArray ds_trackpoints = ds_track->GetTrackPoints();
 
317
    // Get TOF1&2 times to calculate effective particle speed between detectors
 
318
    double TOF1_time = TOFTimeFromTrackPoints(us_trackpoints,
 
319
        DataStructure::Global::kTOF1);
 
320
    double TOF2_time = TOFTimeFromTrackPoints(ds_trackpoints,
 
321
        DataStructure::Global::kTOF2);
 
322
    double TOFdT = TOF2_time - TOF1_time;
 
323
    if (_residuals) {
 
324
      throughfile << TOFdT/_matching_tolerances.at("TOF12dT").second << "\n";
 
325
    }
 
326
    if ((TOFdT > _matching_tolerances.at("TOF12dT").first) and
 
327
        (TOFdT < _matching_tolerances.at("TOF12dT").second)) {
 
328
      // Assemble through track from trackpoints from the matched US and DS tracks
 
329
      DataStructure::Global::Track* through_track = new DataStructure::Global::Track();
 
330
      through_track->set_mapper_name("MapCppGlobalTrackMatching");
 
331
      through_track->set_pid(pid);
 
332
      Squeak::mout(Squeak::debug) << "TrackMatching: US & DS Matched" << std::endl;
 
333
      addTrackRecursive(through_track, us_track, false);
 
334
      addTrackRecursive(through_track, ds_track, true);
 
335
      return through_track;
 
336
    } else {
 
337
      // There may be a small memory leak here because the US and DS tracks don't get
 
338
      // deleted, deleting them here manually causes a segfault. TODO: Investigate
 
339
      return NULL;
 
340
    }
 
341
}
 
342
 
 
343
void TrackMatching::addTrackRecursive(DataStructure::Global::Track* parent,
 
344
                                      DataStructure::Global::Track* child,
 
345
                                      bool doEMR) {
 
346
    parent->AddTrack(child);
 
347
    DataStructure::Global::TrackPointCPArray tps = child->GetTrackPoints();
 
348
    for (auto trackpoint_iter = tps.begin();
 
349
         trackpoint_iter != tps.end();
 
350
         ++trackpoint_iter) {
 
351
        parent->AddTrackPoint(
 
352
              const_cast<DataStructure::Global::TrackPoint*>(*trackpoint_iter));
 
353
    }
 
354
    if (doEMR) {
 
355
        parent->set_emr_range_primary(child->get_emr_range_primary());
 
356
        parent->set_emr_plane_density(child->get_emr_plane_density());
 
357
    }
 
358
}
 
359
 
299
360
void TrackMatching::throughTrack() {
300
 
  ofstream throughfile;
301
 
  if (_residuals) {
302
 
    throughfile.open("match_through.csv", std::ios::out | std::ios::app);
303
 
  }
304
361
  std::vector<DataStructure::Global::PrimaryChain*> us_chains =
305
362
      _global_event->GetUSPrimaryChainOrphans();
306
363
  std::vector<DataStructure::Global::PrimaryChain*> ds_chains =
316
373
      for (size_t i = 0; i < pids.size(); i++) {
317
374
        DataStructure::Global::Track* us_track = (*us_chain_iter)->GetMatchedTrack(pids[i]);
318
375
        DataStructure::Global::Track* ds_track = (*ds_chain_iter)->GetMatchedTrack(pids[i]);
319
 
        if (!us_track or !us_track->HasDetector(DataStructure::Global::kTOF1) or
320
 
            !ds_track or !ds_track->HasDetector(DataStructure::Global::kTOF2)) {
321
 
          continue;
322
 
        }
323
 
        DataStructure::Global::TrackPointCPArray us_trackpoints = us_track->GetTrackPoints();
324
 
        DataStructure::Global::TrackPointCPArray ds_trackpoints = ds_track->GetTrackPoints();
325
 
        // Get TOF1&2 times to calculate effective particle speed between detectors
326
 
        double TOF1_time = TOFTimeFromTrackPoints(us_trackpoints,
327
 
            DataStructure::Global::kTOF1);
328
 
        double TOF2_time = TOFTimeFromTrackPoints(ds_trackpoints,
329
 
            DataStructure::Global::kTOF2);
330
 
        double TOFdT = TOF2_time - TOF1_time;
331
 
        if (_residuals) {
332
 
          throughfile << TOFdT/_matching_tolerances.at("TOF12dT").second << "\n";
333
 
        }
334
 
        if ((TOFdT > _matching_tolerances.at("TOF12dT").first) and
335
 
            (TOFdT < _matching_tolerances.at("TOF12dT").second)) {
336
 
          DataStructure::Global::Track* through_track = new DataStructure::Global::Track();
337
 
          through_track->set_mapper_name("MapCppGlobalTrackMatching");
338
 
          through_track->set_pid(pids[i]);
339
 
          through_track->set_emr_range_primary(ds_track->get_emr_range_primary());
340
 
          through_track->set_emr_plane_density(ds_track->get_emr_plane_density());
341
 
          Squeak::mout(Squeak::debug) << "TrackMatching: US & DS Matched" << std::endl;
342
 
          // Assemble through track from trackpoints from the matched US and DS tracks
343
 
          for (auto trackpoint_iter = us_trackpoints.begin();
344
 
               trackpoint_iter != us_trackpoints.end();
345
 
               ++trackpoint_iter) {
346
 
            through_track->AddTrackPoint(
347
 
                const_cast<DataStructure::Global::TrackPoint*>(*trackpoint_iter));
348
 
          }
349
 
          for (auto trackpoint_iter = ds_trackpoints.begin();
350
 
               trackpoint_iter != ds_trackpoints.end();
351
 
               ++trackpoint_iter) {
352
 
            through_track->AddTrackPoint(
353
 
                const_cast<DataStructure::Global::TrackPoint*>(*trackpoint_iter));
354
 
          }
355
 
          // Add references back to the original tracks
356
 
          through_track->AddTrack(us_track);
357
 
          through_track->AddTrack(ds_track);
358
 
          _global_event->add_track(through_track);
359
 
          through_primary_chain->AddMatchedTrack(through_track);
360
 
        } else {
361
 
          // There may be a small memory leak here because the US and DS tracks don't get
362
 
          // deleted, deleting them here manually causes a segfault. TODO: Investigate
363
 
        }
 
376
        DataStructure::Global::Track* through_track =
 
377
                    TrackMatching::throughMatchTOF(us_track, ds_track, pids[i]);
 
378
        if (through_track == NULL) {
 
379
            continue;
 
380
        }
 
381
        _global_event->add_track(through_track);
 
382
        through_primary_chain->AddMatchedTrack(through_track);
364
383
      }
365
384
      if (through_primary_chain->GetMatchedTracks().size() > 0) {
366
385
        (*us_chain_iter)->set_chain_type(DataStructure::Global::kUS);
403
422
  std::vector<DataStructure::Global::SpacePoint*> *global_spacepoint_array =
404
423
      _global_event->get_space_points();
405
424
  for (size_t i = 0; i < global_spacepoint_array->size(); i++) {
406
 
    if (global_spacepoint_array->at(i)->get_detector() == detector) {
 
425
    if (global_spacepoint_array->at(i) &&
 
426
        global_spacepoint_array->at(i)->get_detector() == detector) {
407
427
      space_points.push_back(global_spacepoint_array->at(i));
408
428
    }
409
429
  }
462
482
  if (spacepoints.size() > 0) {
463
483
    double target_z = spacepoints.at(0)->get_position().Z();
464
484
    try {
465
 
      GlobalTools::propagate(x_in, target_z, field, _max_step_size, pid,
466
 
                             _energy_loss, static_cast<int>(_geo_algo));
 
485
      Propagate(x_in, field, pid, target_z);
467
486
      // To avoid doing the same propagation again for the same point, store the propagated
468
487
      // values back in the TLorentzVectors
469
488
      position.SetX(x_in[1]);
526
545
                     energy, momentum.X(), momentum.Y(), momentum.Z()};
527
546
    // First propagate to just upstream of TOF1 to get the energy during TOF0-1 transit
528
547
    try {
529
 
      GlobalTools::propagate(x_in, tof1_z - 25.0, field, _max_step_size, pid,
530
 
                             _energy_loss, static_cast<int>(_geo_algo));
 
548
      Propagate(x_in, field, pid, tof1_z - 25.0);
531
549
    } catch (Exceptions::Exception exc) {
532
550
      Squeak::mout(Squeak::debug) << exc.what() << std::endl;
533
551
    }
597
615
    TLorentzVector first_hit_pos_err = emr_trackpoints[0]->get_position_error();
598
616
    double target_z = first_hit_pos.Z();
599
617
    try {
600
 
      GlobalTools::propagate(x_in, target_z, field, _max_step_size, pid,
601
 
                             _energy_loss, static_cast<int>(_geo_algo));
 
618
      Propagate(x_in, field, pid, target_z);
602
619
      if (_residuals) {
603
620
        emrfile << first_hit_pos.X() - x_in[1] << " " << first_hit_pos.Y() - x_in[2] << "\n";
604
621
      }
620
637
          emr_tp->set_momentum(momentum);
621
638
          hypothesis_track->AddTrackPoint(emr_tp);
622
639
        }
 
640
        AddVirtualTrackPoints(hypothesis_track);
623
641
      } else {
624
642
        Squeak::mout(Squeak::debug)
625
643
            << "TrackMatching: EMR Mismatch, dx, dy are:\n"
652
670
    tracker_tp->set_momentum(momentum);
653
671
    hypothesis_track->AddTrackPoint(tracker_tp);
654
672
  }
 
673
  AddVirtualTrackPoints(hypothesis_track);
 
674
 
 
675
}
 
676
 
 
677
void TrackMatching::Propagate(double* x_in,
 
678
                              BTFieldConstructor* field,
 
679
                              DataStructure::Global::PID pid,
 
680
                              double target_z) {
 
681
    std::vector<double>::iterator start, end;
 
682
    if (target_z > x_in[3]) {
 
683
        start = std::lower_bound(_extra_z_planes.begin(),
 
684
                                 _extra_z_planes.end(), x_in[3]);
 
685
        end = std::upper_bound(_extra_z_planes.begin(),
 
686
                               _extra_z_planes.end(), target_z);
 
687
    } else {
 
688
        start = std::lower_bound(_extra_z_planes.begin(),
 
689
                                 _extra_z_planes.end(), target_z);
 
690
        end = std::upper_bound(_extra_z_planes.begin(),
 
691
                               _extra_z_planes.end(), x_in[3]);
 
692
    }
 
693
    _virtual_track_points = std::vector<DataStructure::Global::TrackPoint>();
 
694
    for (/*start defined above*/; start < end && start < _extra_z_planes.end(); ++start) {
 
695
        GlobalTools::propagate(x_in, *start, field, _max_step_size, pid,
 
696
                           _energy_loss, static_cast<int>(_geo_algo));
 
697
        DataStructure::Global::TrackPoint tp;
 
698
        tp.set_particle_event(-1); // BUG
 
699
        tp.set_mapper_name("MapCppGlobalTrackMatching");
 
700
        tp.set_detector(DataStructure::Global::kVirtual);
 
701
        tp.set_position(TLorentzVector(x_in[1], x_in[2], x_in[3], x_in[0]));
 
702
        tp.set_momentum(TLorentzVector(x_in[5], x_in[6], x_in[7], x_in[4]));
 
703
        _virtual_track_points.push_back(tp);
 
704
    }
 
705
    GlobalTools::propagate(x_in, target_z, field, _max_step_size, pid,
 
706
                       _energy_loss, static_cast<int>(_geo_algo));
 
707
}
 
708
 
 
709
void TrackMatching::AddVirtualTrackPoints(DataStructure::Global::Track* hypothesis_track) {
 
710
    for (size_t i = 0; i < _virtual_track_points.size(); ++i) {
 
711
        DataStructure::Global::SpacePoint* sp = new DataStructure::Global::SpacePoint();
 
712
        sp->set_detector(DataStructure::Global::kVirtual);
 
713
        sp->set_charge(hypothesis_track->get_charge());
 
714
        sp->set_mapper_name(_mapper_name);
 
715
        DataStructure::Global::TrackPoint* tp = new DataStructure::Global::TrackPoint(_virtual_track_points[i]);
 
716
        tp->set_detector(DataStructure::Global::kVirtual);
 
717
        tp->set_charge(hypothesis_track->get_charge());
 
718
        tp->set_mapper_name(_mapper_name);
 
719
        tp->set_space_point(sp);
 
720
        tp->set_particle_event(-1); // BUG
 
721
        hypothesis_track->AddTrackPoint(tp);
 
722
    }
655
723
}
656
724
 
657
725
double TrackMatching::TOFTimeFromTrackPoints(
718
786
      trackpoint->set_mapper_name(_mapper_name);
719
787
      hypothesis_track->AddTrackPoint(trackpoint);
720
788
    }
 
789
    AddVirtualTrackPoints(hypothesis_track);
721
790
  }
722
791
  return;
723
792
}