← Back to branch summary

~christopher-hunt08/maus/maus_integrated_kalman

« back to all changes in this revision

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

merging in changes in merge branch

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/common_cpp/Recon/Global/PIDVarF.cc
 
1
/* This file is part of MAUS: http://  micewww.pp.rl.ac.uk:8080/projects/maus
 
2
 * 
 
3
 * MAUS is free software: you can redistribute it and/or modify
 
4
 * it under the terms of the GNU General Public License as published by
 
5
 * the Free Software Foundation, either version 3 of the License, or
 
6
 * (at your option) any later version.
 
7
 * 
 
8
 * MAUS is distributed in the hope that it will be useful,
 
9
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
10
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
11
 * GNU General Public License for more details.
 
12
 * 
 
13
 * You should have received a copy of the GNU General Public License
 
14
 * along with MAUS.  If not, see <http://  www.gnu.org/licenses/>.
 
15
 */
 
16
 
 
17
#include <string>
 
18
#include <vector>
 
19
#include "Recon/Global/PIDVarF.hh"
 
20
#include "Utils/Exception.hh"
 
21
 
 
22
namespace MAUS {
 
23
namespace recon {
 
24
namespace global {
 
25
 
 
26
  const std::string PIDVarF::VARIABLE = "EMRrangevsDSTrackerMom";
 
27
 
 
28
  PIDVarF::PIDVarF(std::string hypothesis, std::string unique_identifier)
 
29
    : PIDBase2D(VARIABLE, hypothesis, unique_identifier,
 
30
                 XminBinF, XmaxBinF, XnumBins, YminBinF, YmaxBinF, YnumBins) {
 
31
    _nonZeroHistEntries = true;
 
32
  }
 
33
 
 
34
  PIDVarF::PIDVarF(TFile* file, std::string hypothesis, int XminF, int XmaxF,
 
35
                   int YminF, int YmaxF)
 
36
    : PIDBase2D(file, VARIABLE, hypothesis, XminF, XmaxF, YminF, YmaxF,
 
37
                XminBinF, XmaxBinF, YminBinF, YmaxBinF) {
 
38
  }
 
39
 
 
40
  PIDVarF::~PIDVarF() {}
 
41
 
 
42
  std::pair<double, double> PIDVarF::Calc_Var(MAUS::DataStructure::Global::Track* track) {
 
43
    EMR_range = track->get_emr_range_primary();
 
44
    if (EMR_range == 0.0) {
 
45
      Squeak::mout(Squeak::debug) << "Global track records no range " <<
 
46
        " of particle in EMR, Recon::Global::PIDVarF::Calc_Var()" << std::endl;
 
47
      return std::make_pair(0, -1);
 
48
    } else if ( YminBinF > EMR_range || EMR_range > YmaxBinF ) {
 
49
      Squeak::mout(Squeak::debug) << "Range of particle in EMR " <<
 
50
        "outside of PDF range, Recon::Global::PIDVarF::Calc_Var()" <<
 
51
        std::endl;
 
52
      return std::make_pair(0, -1);
 
53
    } else {
 
54
      // Get trackpoint array from track
 
55
      global_track_points = track->GetTrackPoints();
 
56
      std::vector<const MAUS::DataStructure::Global::TrackPoint*>
 
57
        ::iterator eachTP;
 
58
      for (eachTP = global_track_points.begin();
 
59
           eachTP != global_track_points.end(); ++eachTP) {
 
60
        if (!(*eachTP)) continue;
 
61
        if ((*eachTP)->get_mapper_name() == "MapCppGlobalTrackMatching-Through") {
 
62
          if ((*eachTP)->get_detector() >=
 
63
                     MAUS::DataStructure::Global::kTracker1_1 &&
 
64
                     (*eachTP)->get_detector() <=
 
65
                     MAUS::DataStructure::Global::kTracker1_5) {
 
66
            tracker1_track_points.push_back(*eachTP);
 
67
          } else {
 
68
            continue;
 
69
          }
 
70
        }
 
71
      }
 
72
      if (tracker1_track_points.size() < 1) {
 
73
        Squeak::mout(Squeak::debug) << "Global track contained no downstream" <<
 
74
          " tracker trackpoints, Recon::Global::PIDVarF::Calc_Var()" << std::endl;
 
75
        tracker1_track_points.clear();
 
76
        return std::make_pair(-1, 0);
 
77
      } else {
 
78
        double tracker1_trackpoint_mom = 0;
 
79
        int tracker1_tp_count = 0;
 
80
        double tracker1_momentum = 0;
 
81
        const MAUS::DataStructure::Global::TrackPoint* tracker1_trackpoint;
 
82
        for (size_t i = 0; i < tracker1_track_points.size(); i++) {
 
83
          tracker1_trackpoint = tracker1_track_points[i];
 
84
          if (tracker1_trackpoint) {
 
85
            TLorentzVector trackpoint_4mom = tracker1_trackpoint->get_momentum();
 
86
            tracker1_trackpoint_mom +=
 
87
              sqrt(trackpoint_4mom.Px()*trackpoint_4mom.Px() +
 
88
                   trackpoint_4mom.Py()*trackpoint_4mom.Py() +
 
89
                   trackpoint_4mom.Pz()*trackpoint_4mom.Pz());
 
90
            tracker1_tp_count++;
 
91
          } else {
 
92
            continue;
 
93
          }
 
94
        }
 
95
        tracker1_track_points.clear();
 
96
        tracker1_momentum = tracker1_trackpoint_mom/tracker1_tp_count;
 
97
        if ( XminBinF > tracker1_momentum || tracker1_momentum > XmaxBinF ) {
 
98
          Squeak::mout(Squeak::debug) << "Momentum for tracker 1 is outside " <<
 
99
            "of range, Recon::Global::PIDVarF::Calc_Var()" << std::endl;
 
100
          return std::make_pair(0, -1);
 
101
        } else {
 
102
          return std::make_pair(tracker1_momentum, EMR_range);
 
103
        }
 
104
      }
 
105
    }
 
106
  }
 
107
}
 
108
}
 
109
}
 
110
 
 
111