Summed chi2 for each track now in Reco_Tree output

src/TMS_Kalman.cpp

@@ -383,16 +383,17 @@ void TMS_Kalman::Predict(TMS_KalmanNode &Node) {
 
 
   // Calculate chi^2
-  Node.rVec  = (Measurement - H*UpdateVec);
-  //Node.rVecT = (Measurement - H*UpdateVec); Node.rVecT.T(); // Transpose it
-  Node.RMatrix = (Node.NoiseMatrix - H*UpdatedCovarianceMatrix*H_T).Invert();
-  //TMatrixD RInvMatrix = (Node.RMatrix).Invert();
-  //Node.chi2 = Node.RMatrix*Node.rVec;
-  //Node.chi2 = Node.rVecT*Node.RMatrix*Node.rVec;
-  Node.chi2 = Node.rVec*(Node.RMatrix*Node.rVec);
-  //(Node.RMatrix*Node.rVec).Print();
-  //Node.chi2 *= Node.rVec;
-  std::cout << "chi2: " << Node.chi2 << std::endl;
+  Node.rVec[0] = (Measurement[0] - UpdateVec[0]);
+  Node.rVec[1] = (Measurement[1] - UpdateVec[1]);
+
+  // Probably a much nicer way to make (sub)matrix from a bigger one, but YOLO
+  Node.RMatrix(0,0) = (Node.NoiseMatrix(0,0) - UpdatedCovarianceMatrix(0,0));
+  Node.RMatrix(1,0) = (Node.NoiseMatrix(1,0) - UpdatedCovarianceMatrix(1,0));
+  Node.RMatrix(0,1) = (Node.NoiseMatrix(0,1) - UpdatedCovarianceMatrix(0,1));
+  Node.RMatrix(1,1) = (Node.NoiseMatrix(1,1) - UpdatedCovarianceMatrix(1,1));
+  Node.RMatrix.Invert(); // Matrix has to be inverted
+  Node.chi2 = Node.rVec*(Node.RMatrix*Node.rVec); // Calc chi^2
+  //std::cout << "chi2: " << Node.chi2 << std::endl;
 
 
 

src/TMS_Kalman.h

@@ -58,25 +58,25 @@ class TMS_KalmanNode {
     x(xvar), y(yvar), z(zvar), dz(dzvar),
     CurrentState(x, y, z+dz, -999.9, -999.9, -1./20.), // Initialise the state vectors
     PreviousState(x, y, z, -999.9, -999.9, -1./20.),
-    TransferMatrix(KALMAN_DIM,KALMAN_DIM),
-    TransferMatrixT(KALMAN_DIM,KALMAN_DIM),
-    NoiseMatrix(KALMAN_DIM,KALMAN_DIM),
-    CovarianceMatrix(KALMAN_DIM,KALMAN_DIM),
-    UpdatedCovarianceMatrix(KALMAN_DIM,KALMAN_DIM),
-    MeasurementMatrix(KALMAN_DIM,KALMAN_DIM),
-    rVec(KALMAN_DIM),
-    rVecT(KALMAN_DIM),
-    RMatrix(KALMAN_DIM,KALMAN_DIM)
+    TransferMatrix(KALMAN_DIM, KALMAN_DIM),
+    TransferMatrixT(KALMAN_DIM, KALMAN_DIM),
+    NoiseMatrix(KALMAN_DIM, KALMAN_DIM),
+    CovarianceMatrix(KALMAN_DIM, KALMAN_DIM),
+    UpdatedCovarianceMatrix(KALMAN_DIM, KALMAN_DIM),
+    MeasurementMatrix(KALMAN_DIM, KALMAN_DIM),
+    rVec(2),
+    rVecT(2),
+    RMatrix(2, 2)
   {
     TransferMatrix.ResizeTo(KALMAN_DIM, KALMAN_DIM);
     TransferMatrixT.ResizeTo(KALMAN_DIM, KALMAN_DIM);
     NoiseMatrix.ResizeTo(KALMAN_DIM, KALMAN_DIM);
     CovarianceMatrix.ResizeTo(KALMAN_DIM, KALMAN_DIM);
     UpdatedCovarianceMatrix.ResizeTo(KALMAN_DIM, KALMAN_DIM);
     MeasurementMatrix.ResizeTo(KALMAN_DIM, KALMAN_DIM);
-    rVec.ResizeTo(KALMAN_DIM);
-    rVecT.ResizeTo(KALMAN_DIM);
-    RMatrix.ResizeTo(KALMAN_DIM, KALMAN_DIM);
+    rVec.ResizeTo(2);
+    rVecT.ResizeTo(2);
+    RMatrix.ResizeTo(2, 2);
 
     // Make the transfer matrix for each of the states
     // Initialise to zero
@@ -247,6 +247,16 @@ class TMS_Kalman {
 
     double GetMomentum() {return momentum;}
 
+    double GetTrackChi2()
+    {
+      double tmp_chi2 = 0.0;
+      for (auto node : KalmanNodes)
+        tmp_chi2 += node.chi2;
+
+      return tmp_chi2;
+    }
+
+
     std::vector<TMS_KalmanNode> GetKalmanNodes() {return KalmanNodes;}
 
     TVectorD GetNoiseVector(TMS_KalmanNode Node);

src/TMS_Reco.cpp

@@ -727,7 +727,7 @@ void TMS_TrackFinder::FindTracks(TMS_Event &event) {
 
   // Run Kalman filter if requested
   if (TMS_Manager::GetInstance().Get_Reco_Kalman_Run()) {
-    double kalman_reco_mom;
+    double kalman_reco_mom, kalman_chi2;
     for (auto &trk : HoughTracks3D) {
       KalmanFilter = TMS_Kalman(trk.Hits);
       kalman_reco_mom = KalmanFilter.GetMomentum();
@@ -751,6 +751,9 @@ void TMS_TrackFinder::FindTracks(TMS_Event &event) {
 
       // Add tracklength with Kalman filter
       trk.Length = CalculateTrackLengthKalman(trk);
+      //std::cout << "chi2 : " << KalmanFilter.GetTrackChi2() << std::endl;
+      kalman_chi2 = KalmanFilter.GetTrackChi2();
+      trk.SetChi2(kalman_chi2);
     }
   } else { // No Kalman filter enabled
     for (auto &trk : HoughTracks3D) {

src/TMS_Track.h

@@ -27,17 +27,20 @@ class TMS_Track {
     double EnergyRange;
     double Momentum;
     double Time;         // TODO: Fill this in a sensible way
+    double Chi2;
 
     double GetEnergyDeposit() {return EnergyDeposit;};
     double GetEnergyRange()   {return EnergyRange;};
     double GetMomentum()      {return Momentum;};
+    double GetChi2()          {return Chi2;};
 
     TMS_TrueParticle GetTrueParticle() {return fTrueParticle;};
 
     // Manually set variables
     void SetEnergyDeposit (double val) {EnergyDeposit = val;};
     void SetEnergyRange   (double val) {EnergyRange   = val;};
     void SetMomentum      (double val) {Momentum      = val;};
+    void SetChi2          (double val) {Chi2          = val;};
 
     // Set direction unit vectors from only x and y slope
     void SetStartDirection(double ax, double ay, double az);// {StartDirection[0]=ax; StartDirection[1]=ay; StartDirection[2]=az;};

src/TMS_TreeWriter.cpp

@@ -243,6 +243,7 @@ void TMS_TreeWriter::MakeBranches() {
   Reco_Tree->Branch("Length",         RecoTrackLength,          "Length[nTracks]/F");
   Reco_Tree->Branch("Charge",         RecoTrackCharge,          "Charge[nTracks]/I");
   Reco_Tree->Branch("TrackHitEnergies", RecoTrackHitEnergies,   "RecoTrackHitEnergies[10][200]/F");
+  Reco_Tree->Branch("Chi2",           RecoTrackChi2,             "Chi2[nTracks]/I");
 
   MakeTruthBranches(Truth_Info);
   MakeTruthBranches(Truth_Spill);
@@ -1222,6 +1223,7 @@ void TMS_TreeWriter::Fill(TMS_Event &event) {
     RecoTrackEnergyDeposit[itTrack] =       RecoTrack->EnergyDeposit;
     RecoTrackMomentum[itTrack]      =       RecoTrack->Momentum;
     RecoTrackCharge[itTrack]        =       RecoTrack->Charge;
+    RecoTrackChi2[itTrack]          =       RecoTrack->Chi2;
 
     for (int j = 0; j < 3; j++) {
       RecoTrackStartPos[itTrack][j]  = RecoTrack->Start[j];

src/TMS_TreeWriter.h

@@ -59,6 +59,7 @@ class TMS_TreeWriter {
     float RecoTrackMomentum[__TMS_MAX_TRACKS__];
     float RecoTrackTrueMomentum[__TMS_MAX_TRACKS__];
     float RecoTrackLength[__TMS_MAX_TRACKS__];
+    float RecoTrackChi2[__TMS_MAX_TRACKS__];
     int RecoTrackCharge[__TMS_MAX_TRACKS__];
 
   private: