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