From 216b0a31e5ee608fea8d87df4c30a0bf7b181a66 Mon Sep 17 00:00:00 2001
From: Dmitry Shchukin <dmitry.shchukin@cern.ch>
Date: Wed, 14 Dec 2022 00:29:45 +0100
Subject: [PATCH 1/3] changes based on exampletrackingmc.cc

---
 mkcode.cc | 92 +++++++++++++++++++++++++++++++++----------------------
 1 file changed, 56 insertions(+), 36 deletions(-)

diff --git a/mkcode.cc b/mkcode.cc
index d0f7770..639b77d 100644
--- a/mkcode.cc
+++ b/mkcode.cc
@@ -28,12 +28,11 @@ using namespace CS;
 // timing
 #include "timecut.h"
 
-// tracking 2021mu
+// tracking tools
 #ifdef TRACKINGTOOLS
 #include "tracking_new.h"
 #endif
 
-//MC simpletracking
 #include "simu.h"
 
 int procev(RecoEvent &e);
@@ -253,14 +252,16 @@ int main(int argc, char *argv[])
       return 1;
     }
   
-    const int Nevents = 100000000;
+    //const int Nevents = 100000000;
+    const int Nevents = 100; //for testing
     int NevMCRead = 0;
 
     // event loop, read event by event
     for (int iev = 0; iev < Nevents; iev++) {
 
       RecoEvent e;
-      if(inputtype == 1) e = RunP348ReadMC(inFile);
+      //if(inputtype == 1) e = RunP348ReadMC(inFile);
+      if(inputtype == 1) e = RunP348RecoMC(inFile);//same as in exampletrackingmc in master for proper initialization etc.
       if(inputtype == 2) e = RunP348ReadMCOldFormat(inFile);
 
       if (!e.mc) break;
@@ -423,12 +424,12 @@ int procev(RecoEvent &e0)
 //#include "mkcuts_hcontamination_control.inc"
 //#include "mkcuts2021_calib_e.inc"
 //#include "mkcuts2021_invis.inc"
-//#include "mkcuts_calib_e.inc"
+#include "mkcuts_calib_e.inc"
 //#include "mkcuts_calib_pi.inc"
 //#include "mkcuts_calib_pi_2021mu.inc"
 //#include "mkcuts_tracking_2021mu.inc"
 //#include "mkcuts_calib_pi_hadrinvis50.inc"
-#include "mkcuts_hadrinvis50.inc"
+//#include "mkcuts_hadrinvis50.inc"
 //#include "mkcuts_calib_mu.inc"
 //#include "mkcuts_invis.inc"
 //#include "mkcuts_invis_herm.inc"
@@ -565,38 +566,56 @@ int procev(RecoEvent &e0)
     }
 #endif
 #ifdef TRACKINGTOOLS
-    #include <math.h>
-    double mymom = -0.5;
-    double mychisq = -0.5;
-    if(!e.mc) {
-      bool ifupstream = true;
-      // to be set to false for downstream momentum in the muon runs
-      const std::vector<Track_t> mytracks = GetReconstructedTracks(e, ifupstream);
-      for (unsigned int p = 0; p < mytracks.size(); p++) {
-        Track_t track = mytracks.at(p);
-        mymom = track.momentum;
-        mychisq = track.chisquare;        
-        if(mymom > 159.9999 && mymom < 160.0001) mychisq = 500.; // temporary
-        if(mymom > 99.9999 && mymom < 100.0001) mychisq = 500.; // temporary
-      }
-      if(mytracks.size() > 1) mychisq = 500.;
-      if(mytracks.size() > 1) std::cout << "More than 1 track, N = " << mytracks.size() << std::endl;
-
-      if(e.run >= 8396 && e.run <= 8458) mymom *= 0.5; // FIXME: temporary fix for the hadron run 2022 with mom = 50.
-
-      //detect bad values of chi squared and apply chi squared cut
-      if ( isnan((float)mychisq) || isinf((float)mychisq) ) { ipass = 0; }
-      else {
-        if(MomentumLowerCut > -0.5) {
-          if(mychisq < MomentumChiSqLowerCut) ipass = 0;
-          if(mychisq > MomentumChiSqUpperCut) ipass = 0;
+    double mymom = 9999.;
+    double mychisq = 9999.;
+    if(e.mc) {
+      //const vector<const Micromega*> mm_map
+      //= {&e.MM1, &e.MM2, &e.MM3, &e.MM4, &e.MM5, &e.MM6, &e.MM7, &e.MM8};
+      //const vector<const vector<MChit>*> mmtruth_map
+      // = {&e.mc->MM1truth, &e.mc->MM2truth, &e.mc->MM3truth, &e.mc->MM4truth,
+      //    &e.mc->MM5truth, &e.mc->MM6truth, &e.mc->MM7truth, &e.mc->MM8truth};
+
+      // track reconstruction with Tracking Tools
+      // RECO HITS
+      std::vector<std::pair<TVector3, double> > hits;
+      FillAbsMMHits(e.MM1, hits);
+      FillAbsMMHits(e.MM2, hits);
+      FillAbsMMHits(e.MM3, hits);
+      FillAbsMMHits(e.MM4, hits);
+
+      const bool debug = true;
+
+      if(debug)
+        for(auto hit : hits)
+          cout << "Reco point in MM: (" << hit.first.x() << ", " << hit.first.y() << ", " << hit.first.z() << ")" << endl;
+      const std::vector<Track_t> tracks = GetReconstructedTracksFromHits(hits);
+
+      // fill
+      if(!tracks.empty()) {
+        double _best_chi2(9999);
+        double _best_mom(9999);
+        // Look for the best track with the lowest chi2 value
+        for (auto elem : tracks) {
+          if (elem.chisquare < _best_chi2 && elem.chisquare > -1.) {
+            _best_chi2 = elem.chisquare;
+            _best_mom = elem.momentum;
+          }
         }
+        mychisq = _best_chi2;
+        mymom = _best_mom;
+
+        if (debug) std::cout << "mymom = " << mymom <<std::endl;
+        if (debug) std::cout << "mychisq = " << mychisq <<std::endl;
+      }
+    } else {
+        std::cout << "Trackingtools not implemented for data in mkcode.exe, exiting..." << std::endl;
+        exit(1);
       }
-    }
 #endif
-
-  if(mymom < MomentumLowerCut) ipass = 0;
-  if(mymom > MomentumUpperCut) ipass = 0;
+  
+  //remove momentum cuts for testing
+  //if(mymom < MomentumLowerCut) ipass = 0;
+  //if(mymom > MomentumUpperCut) ipass = 0;
 
   if(ipass) {
     if(e.mc) {
@@ -1006,7 +1025,8 @@ int procev(RecoEvent &e0)
 
   double mychi2 = 0.;
   //mychi2 = calcShowerChi2(e, 2, 2, 3, 3);
-  mychi2 = calcShowerChi2(e, imaxcellx, imaxcelly, 3, SP_ECAL);
+  //mychi2 = calcShowerChi2(e, imaxcellx, imaxcelly, 3, SP_ECAL);
+  mychi2 = calcShowerChi2(e, imaxcellx, imaxcelly, 3, SP_MC);
   if(ECALChi2UpperCut < 101. && mychi2 > ECALChi2UpperCut) ipass = 0;
   if(ECALMaxCell > -1) { // check max cell
     if(e.mc || ECALMaxCell <= 29) { // set needed cell numbers from cuts
-- 
GitLab


From d580f7cde3d098d0b4ca3131aae7d60cfa521dac Mon Sep 17 00:00:00 2001
From: Benjamin Banto Oberhauser <benjamin.banto.oberhauser@cern.ch>
Date: Sat, 17 Dec 2022 11:45:38 +0100
Subject: [PATCH 2/3] Fix to use tracking-tools. When tracking-tools are used,
 the MManager class should not be called.

---
 mkcode.cc | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/mkcode.cc b/mkcode.cc
index 639b77d..af1eec4 100644
--- a/mkcode.cc
+++ b/mkcode.cc
@@ -23,7 +23,9 @@ using namespace CS;
 #include "shower.h"
 #include "badburst.h"
 #include "tracking.h"
+#ifndef TRACKINGTOOLS
 #include "MManager.h"
+#endif
 
 // timing
 #include "timecut.h"
-- 
GitLab


From ffc1f3b4c225bba4d8095ae8153905c1e11196ea Mon Sep 17 00:00:00 2001
From: Dmitry Shchukin <dmitry.shchukin@cern.ch>
Date: Tue, 20 Dec 2022 18:10:18 +0100
Subject: [PATCH 3/3] Cleanup

---
 mkcode.cc | 27 ++++++++-------------------
 1 file changed, 8 insertions(+), 19 deletions(-)

diff --git a/mkcode.cc b/mkcode.cc
index af1eec4..219b951 100644
--- a/mkcode.cc
+++ b/mkcode.cc
@@ -35,6 +35,7 @@ using namespace CS;
 #include "tracking_new.h"
 #endif
 
+// for MC with tracking-tools
 #include "simu.h"
 
 int procev(RecoEvent &e);
@@ -254,16 +255,14 @@ int main(int argc, char *argv[])
       return 1;
     }
   
-    //const int Nevents = 100000000;
-    const int Nevents = 100; //for testing
+    const int Nevents = 100000000;
     int NevMCRead = 0;
 
     // event loop, read event by event
     for (int iev = 0; iev < Nevents; iev++) {
 
       RecoEvent e;
-      //if(inputtype == 1) e = RunP348ReadMC(inFile);
-      if(inputtype == 1) e = RunP348RecoMC(inFile);//same as in exampletrackingmc in master for proper initialization etc.
+      if(inputtype == 1) e = RunP348RecoMC(inFile); //same as in exampletrackingmc in master
       if(inputtype == 2) e = RunP348ReadMCOldFormat(inFile);
 
       if (!e.mc) break;
@@ -568,25 +567,17 @@ int procev(RecoEvent &e0)
     }
 #endif
 #ifdef TRACKINGTOOLS
+    // track reconstruction with Tracking Tools
     double mymom = 9999.;
     double mychisq = 9999.;
     if(e.mc) {
-      //const vector<const Micromega*> mm_map
-      //= {&e.MM1, &e.MM2, &e.MM3, &e.MM4, &e.MM5, &e.MM6, &e.MM7, &e.MM8};
-      //const vector<const vector<MChit>*> mmtruth_map
-      // = {&e.mc->MM1truth, &e.mc->MM2truth, &e.mc->MM3truth, &e.mc->MM4truth,
-      //    &e.mc->MM5truth, &e.mc->MM6truth, &e.mc->MM7truth, &e.mc->MM8truth};
-
-      // track reconstruction with Tracking Tools
-      // RECO HITS
+      const bool debug = false;
       std::vector<std::pair<TVector3, double> > hits;
       FillAbsMMHits(e.MM1, hits);
       FillAbsMMHits(e.MM2, hits);
       FillAbsMMHits(e.MM3, hits);
       FillAbsMMHits(e.MM4, hits);
 
-      const bool debug = true;
-
       if(debug)
         for(auto hit : hits)
           cout << "Reco point in MM: (" << hit.first.x() << ", " << hit.first.y() << ", " << hit.first.z() << ")" << endl;
@@ -615,9 +606,8 @@ int procev(RecoEvent &e0)
       }
 #endif
   
-  //remove momentum cuts for testing
-  //if(mymom < MomentumLowerCut) ipass = 0;
-  //if(mymom > MomentumUpperCut) ipass = 0;
+  if(mymom < MomentumLowerCut) ipass = 0;
+  if(mymom > MomentumUpperCut) ipass = 0;
 
   if(ipass) {
     if(e.mc) {
@@ -1027,8 +1017,7 @@ int procev(RecoEvent &e0)
 
   double mychi2 = 0.;
   //mychi2 = calcShowerChi2(e, 2, 2, 3, 3);
-  //mychi2 = calcShowerChi2(e, imaxcellx, imaxcelly, 3, SP_ECAL);
-  mychi2 = calcShowerChi2(e, imaxcellx, imaxcelly, 3, SP_MC);
+  mychi2 = calcShowerChi2(e, imaxcellx, imaxcelly, 3, SP_ECAL);
   if(ECALChi2UpperCut < 101. && mychi2 > ECALChi2UpperCut) ipass = 0;
   if(ECALMaxCell > -1) { // check max cell
     if(e.mc || ECALMaxCell <= 29) { // set needed cell numbers from cuts
-- 
GitLab