diff --git a/mkcode.cc b/mkcode.cc index d0f77708ca7a06826e41b3db616b16f1accec59c..219b951f9a8a7988fa9e3a16bb63276a5bdf1ce9 100644 --- a/mkcode.cc +++ b/mkcode.cc @@ -23,17 +23,19 @@ using namespace CS; #include "shower.h" #include "badburst.h" #include "tracking.h" +#ifndef TRACKINGTOOLS #include "MManager.h" +#endif // timing #include "timecut.h" -// tracking 2021mu +// tracking tools #ifdef TRACKINGTOOLS #include "tracking_new.h" #endif -//MC simpletracking +// for MC with tracking-tools #include "simu.h" int procev(RecoEvent &e); @@ -260,7 +262,7 @@ int main(int argc, char *argv[]) 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 if(inputtype == 2) e = RunP348ReadMCOldFormat(inFile); if (!e.mc) break; @@ -423,12 +425,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,36 +567,45 @@ 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; + // track reconstruction with Tracking Tools + double mymom = 9999.; + double mychisq = 9999.; + if(e.mc) { + 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); + + 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;