From 97a973d5989e9e58bd32fa682f5c47c7c2413600 Mon Sep 17 00:00:00 2001 From: gianelle Date: Fri, 11 Oct 2024 15:31:24 +0200 Subject: [PATCH 1/5] Add new filter parameters: - MaxHoles - MinNdf - Neural Network (BDT like) parameters --- source/Utils/include/FilterTracks.h | 13 ++- source/Utils/src/FilterTracks.cc | 127 +++++++++++++++++++++++----- 2 files changed, 117 insertions(+), 23 deletions(-) diff --git a/source/Utils/include/FilterTracks.h b/source/Utils/include/FilterTracks.h index 3e5cfc9..3c94d6f 100644 --- a/source/Utils/include/FilterTracks.h +++ b/source/Utils/include/FilterTracks.h @@ -70,14 +70,25 @@ class FilterTracks : public marlin::Processor int _NHitsInner = 2; //! Cut off for number of hits in outer tracker (barrel and endcap combined) int _NHitsOuter = 1; + //! Cut off for number of holes + int _MaxHoles = 1; //! Cut off for momentum (GeV) float _MinPt = 1.0; //units GeV + //! Cut off for the value ndf + int _MinNdf = 1; + //! Cut off for spatial and temporal chi squared values float _Chi2Spatial = 0; + //! NN parameters + std::string _NNmethod = ""; // if defined apply the NN + std::string _NNweights = ""; // xml file with weights + std::vector _NNvars; // sorted list of variables used by NN + float _NNthr = 0; // NN threshold + //! Default magnetic field value (Tesla) - float _Bz = 3.57; + float _Bz = 5.0; }; diff --git a/source/Utils/src/FilterTracks.cc b/source/Utils/src/FilterTracks.cc index 4e28536..3a8b3d8 100644 --- a/source/Utils/src/FilterTracks.cc +++ b/source/Utils/src/FilterTracks.cc @@ -10,8 +10,13 @@ #include #include +#include "TMVA/Reader.h" + +#include + FilterTracks aFilterTracks ; + FilterTracks::FilterTracks() : Processor("FilterTracks") { @@ -56,11 +61,54 @@ FilterTracks::FilterTracks() ); registerProcessorParameter("Chi2Spatial", - "Spatial chi squared", + "iMinimum value for Spatial chi squared", _Chi2Spatial, _Chi2Spatial ); + registerProcessorParameter("MinNdf", + "Minimum value for ndf", + _MinNdf, + _MinNdf + ); + + registerProcessorParameter("MaxHoles", + "Max number of holes", + _MaxHoles, + _MaxHoles + ); + + registerProcessorParameter("Bz", + "Magnetic field in Tesla (default: 5)", + _Bz, + _Bz + ); + + registerProcessorParameter("NNmethod", + "Name of the NN method, if empty uses standard cuts", + _NNmethod, + std::string("") + ); + + registerProcessorParameter("NNweights", + "File xml with the weights of the NN", + _NNweights, + std::string("") + ); + + registerProcessorParameter("NNvars", + "Sorted list with the names of the variables used in the training" + "Supported variables are: trtvhn, trtihn, trtohn, trthn, trtnh, trch2, trndf", + _NNvars, + _NNvars + ); + + registerProcessorParameter("NNthr", + "NN threshold (Takes tracks with prediction > NNthr)", + _NNthr, + _NNthr + ); + registerInputCollection( LCIO::TRACK, "InputTrackCollectionName" , "Name of the input collection", @@ -81,7 +129,7 @@ void FilterTracks::init() { // Print the initial parameters printParameters() ; - buildBfield() ; + //buildBfield() ; } void FilterTracks::processRunHeader( LCRunHeader* /*run*/) @@ -108,6 +156,8 @@ void FilterTracks::processEvent( LCEvent * evt ) LCCollectionVec *OutputTrackCollection = new LCCollectionVec(LCIO::TRACK); OutputTrackCollection->setSubset(true); + TMVA::Reader* reader = new TMVA::Reader(); + // Get input collection LCCollection* InputTrackCollection =evt->getCollection(_InputTrackCollection); @@ -118,45 +168,78 @@ void FilterTracks::processEvent( LCEvent * evt ) std::string encoderString = lcio::LCTrackerCellID::encoding_string(); UTIL::CellIDDecoder decoder(encoderString); + std::unordered_map vars = { + {"trtvhn", 0}, + {"trtihn", 0}, + {"trtohn", 0}, + {"trthn", 0}, + {"trtnh", 0}, + {"trch2", 0}, + {"trndf", 0}, + }; + + if ( not _NNmethod.empty() ) { + + for (unsigned i=0,N=_NNvars.size(); iAddVariable( name, &vars[name] ); + } + reader->BookMVA(_NNmethod, _NNweights); + } + for(int i=0; igetNumberOfElements(); i++) { EVENT::Track *trk=dynamic_cast(InputTrackCollection->getElementAt(i)); - int nhittotal = trk->getTrackerHits().size(); + vars["trtnh"] = trk->getTrackerHits().size(); const EVENT::IntVec& subdetectorHits = trk->getSubdetectorHitNumbers(); - int nhitvertex = subdetectorHits[1]+subdetectorHits[2]; - int nhitinner = subdetectorHits[3]+subdetectorHits[4]; - int nhitouter = subdetectorHits[5]+subdetectorHits[6]; + vars["trtvhn"] = subdetectorHits[1]+subdetectorHits[2]; + vars["trtihn"] = subdetectorHits[3]+subdetectorHits[4]; + vars["trtohn"] = subdetectorHits[5]+subdetectorHits[6]; + vars["trthn"] = trk->getNholes(); float pt = fabs(0.3*_Bz/trk->getOmega()/1000); - float chi2spatial = trk->getChi2(); + vars["trch2"] = trk->getChi2(); + + vars["trndf"] = trk->getNdf(); if(_BarrelOnly == true) { bool endcaphits = false; - for(int j=0; jgetTrackerHits()[j])["system"]; - if(systemID == 2 || systemID == 4 || systemID == 6) { - endcaphits = true; - break; - } + for(int j=0; jgetTrackerHits()[j])["system"]; + if(systemID == 2 || systemID == 4 || systemID == 6) { + endcaphits = true; + break; + } } if(endcaphits == false) { OutputTrackCollection->addElement(trk); } } else { // track property cuts - if(nhittotal > _NHitsTotal && - nhitvertex > _NHitsVertex && - nhitinner > _NHitsInner && - nhitouter > _NHitsOuter && - pt > _MinPt && - chi2spatial > _Chi2Spatial) - { OutputTrackCollection->addElement(trk); } + if ( not _NNmethod.empty() ) { // NN cuts + auto prediction = reader->EvaluateMVA(_NNmethod); + if ( prediction > _NNthr ) + OutputTrackCollection->addElement(trk); + } else { // user cuts + if (vars["trtnh"] > _NHitsTotal && + vars["trtvhn"] > _NHitsVertex && + vars["trtihn"] > _NHitsInner && + vars["trtohn"] > _NHitsOuter && + pt > _MinPt && + vars["trch2"] > _Chi2Spatial && + vars["trndf"] > _MinNdf && + vars["trthn"] < _MaxHoles) + { OutputTrackCollection->addElement(trk); } + } } } // Save output track collection evt->addCollection(OutputTrackCollection, _OutputTrackCollection); + delete reader; } void FilterTracks::end() -{ } +{ } From 05f4a0a64216761044ac4a91f19c2044b26518a6 Mon Sep 17 00:00:00 2001 From: gianelle Date: Mon, 14 Oct 2024 16:44:12 +0200 Subject: [PATCH 2/5] Add TMVA --- CMakeLists.txt | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index bd8823e..5acca25 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -68,8 +68,8 @@ find_package( DD4hep REQUIRED COMPONENTS DDRec) set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} ${DD4hep_ROOT}/cmake ) include( DD4hep ) -find_package( ROOT REQUIRED ) -set( ROOT_COMPONENT_LIBRARIES Geom Reflex) +#find_package( ROOT REQUIRED ) +#set( ROOT_COMPONENT_LIBRARIES Geom Reflex TMVA) if(DD4HEP_USE_XERCESC) find_package(XercesC) @@ -78,8 +78,7 @@ include(DD4hep_XML_setup) INCLUDE_DIRECTORIES( BEFORE SYSTEM ${DD4hep_INCLUDE_DIRS} ) LINK_LIBRARIES( ${DD4hep_LIBRARIES} ${DD4hep_COMPONENT_LIBRARIES} ) -#FIND_PACKAGE( ROOT REQUIRED ) - +FIND_PACKAGE( ROOT REQUIRED COMPONENTS Geom TMVA) INCLUDE_DIRECTORIES( SYSTEM ${ROOT_INCLUDE_DIRS} ) LINK_LIBRARIES( ${ROOT_LIBRARIES} ) ADD_DEFINITIONS( ${ROOT_DEFINITIONS} ) From 692194c722fd4be5d6335974100e29e85695561b Mon Sep 17 00:00:00 2001 From: gianelle Date: Mon, 14 Oct 2024 17:04:54 +0200 Subject: [PATCH 3/5] Add new parameter in the filter: MaxPt --- source/Utils/include/FilterTracks.h | 3 ++- source/Utils/src/FilterTracks.cc | 10 +++++++++- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/source/Utils/include/FilterTracks.h b/source/Utils/include/FilterTracks.h index 3c94d6f..564ec0f 100644 --- a/source/Utils/include/FilterTracks.h +++ b/source/Utils/include/FilterTracks.h @@ -74,7 +74,8 @@ class FilterTracks : public marlin::Processor int _MaxHoles = 1; //! Cut off for momentum (GeV) - float _MinPt = 1.0; //units GeV + float _MinPt = 0.5; //units GeV + float _MaxPt = 1000.0; //units GeV //! Cut off for the value ndf int _MinNdf = 1; diff --git a/source/Utils/src/FilterTracks.cc b/source/Utils/src/FilterTracks.cc index 3a8b3d8..cab4077 100644 --- a/source/Utils/src/FilterTracks.cc +++ b/source/Utils/src/FilterTracks.cc @@ -60,6 +60,12 @@ FilterTracks::FilterTracks() _MinPt ); + registerProcessorParameter("MaxPt", + "Max transverse momentum", + _MaxPt, + _MaxPt + ); + registerProcessorParameter("Chi2Spatial", "iMinimum value for Spatial chi squared", _Chi2Spatial, @@ -129,7 +135,8 @@ void FilterTracks::init() { // Print the initial parameters printParameters() ; - //buildBfield() ; + // it require that geometry is initialized + // buildBfield() ; } void FilterTracks::processRunHeader( LCRunHeader* /*run*/) @@ -228,6 +235,7 @@ void FilterTracks::processEvent( LCEvent * evt ) vars["trtihn"] > _NHitsInner && vars["trtohn"] > _NHitsOuter && pt > _MinPt && + pt < _MaxPt && vars["trch2"] > _Chi2Spatial && vars["trndf"] > _MinNdf && vars["trthn"] < _MaxHoles) From 61b7e439f767f8fe76ea62dbc8da9141ca6f6f72 Mon Sep 17 00:00:00 2001 From: gianelle Date: Fri, 29 Nov 2024 10:53:09 +0100 Subject: [PATCH 4/5] Add possibility to cut on theta and on hits outliers --- source/Utils/include/FilterTracks.h | 15 ++++++++++++++- source/Utils/src/FilterTracks.cc | 30 +++++++++++++++++++++++++---- 2 files changed, 40 insertions(+), 5 deletions(-) diff --git a/source/Utils/include/FilterTracks.h b/source/Utils/include/FilterTracks.h index 564ec0f..5d55b7c 100644 --- a/source/Utils/include/FilterTracks.h +++ b/source/Utils/include/FilterTracks.h @@ -16,7 +16,13 @@ namespace TrackPerf * @parameter NHitsVertex Minimum number of hits on vertex detector * @parameter NHitsInner Minimum number of hits on inner tracker * @parameter NHitsOuter Minimum number of hits on outer tracker + * @parameter MaxOutliers Maximum number of outliers hits on track + * @parameter MaxHoles Maximum number of holes on track + * @parameter MinNdf Minimum value for ndf * @parameter MinPt Minimum transverse momentum + * @parameter MaxPt Max transverse momentum + * @parameter MinTheta Minimum theta + * @parameter MaxTheta Max theta * @parameter Chi2Spatial Spatial chi squared * * @author N. Bruhwiler @@ -77,8 +83,15 @@ class FilterTracks : public marlin::Processor float _MinPt = 0.5; //units GeV float _MaxPt = 1000.0; //units GeV + //! Cut off for theta (rad) + float _MinTheta = 0; + float _MaxTheta = 3.14; + //! Cut off for the value ndf - int _MinNdf = 1; + int _MinNdf = 1; + + //! Cut off for outliers number + int _MaxOutl = 0; //! Cut off for spatial and temporal chi squared values float _Chi2Spatial = 0; diff --git a/source/Utils/src/FilterTracks.cc b/source/Utils/src/FilterTracks.cc index cab4077..c014d24 100644 --- a/source/Utils/src/FilterTracks.cc +++ b/source/Utils/src/FilterTracks.cc @@ -7,6 +7,7 @@ #include #include #include +#include #include #include @@ -66,6 +67,18 @@ FilterTracks::FilterTracks() _MaxPt ); + registerProcessorParameter("MinTheta", + "Minimum theta", + _MinTheta, + _MinTheta + ); + + registerProcessorParameter("MaxTheta", + "Max theta", + _MaxTheta, + _MaxTheta + ); + registerProcessorParameter("Chi2Spatial", "iMinimum value for Spatial chi squared", _Chi2Spatial, @@ -84,6 +97,12 @@ FilterTracks::FilterTracks() _MaxHoles ); + registerProcessorParameter("MaxOutliers", + "Max number of outliers", + _MaxOutl, + _MaxOutl + ); + registerProcessorParameter("Bz", "Magnetic field in Tesla (default: 5)", _Bz, @@ -135,7 +154,7 @@ void FilterTracks::init() { // Print the initial parameters printParameters() ; - // it require that geometry is initialized + // it requires that geometry has been already initialized // buildBfield() ; } @@ -166,8 +185,8 @@ void FilterTracks::processEvent( LCEvent * evt ) TMVA::Reader* reader = new TMVA::Reader(); // Get input collection - LCCollection* InputTrackCollection =evt->getCollection(_InputTrackCollection); - + LCCollection* InputTrackCollection = evt->getCollection(_InputTrackCollection); + if( InputTrackCollection->getTypeName() != lcio::LCIO::TRACK ) { throw EVENT::Exception( "Invalid collection type: " + InputTrackCollection->getTypeName() ) ; } @@ -208,7 +227,7 @@ void FilterTracks::processEvent( LCEvent * evt ) vars["trthn"] = trk->getNholes(); float pt = fabs(0.3*_Bz/trk->getOmega()/1000); - + float theta = M_PI_2-atan(trk->getTanLambda()); vars["trch2"] = trk->getChi2(); vars["trndf"] = trk->getNdf(); @@ -236,8 +255,11 @@ void FilterTracks::processEvent( LCEvent * evt ) vars["trtohn"] > _NHitsOuter && pt > _MinPt && pt < _MaxPt && + theta > _MinTheta && + theta < _MaxTheta && vars["trch2"] > _Chi2Spatial && vars["trndf"] > _MinNdf && + vars["trtnh"]-vars["trndf"]/2 < _MaxOutl && vars["trthn"] < _MaxHoles) { OutputTrackCollection->addElement(trk); } } From 74100ae6c19cf8f3b6bd568dcf64babd1f842a06 Mon Sep 17 00:00:00 2001 From: gianelle Date: Wed, 11 Dec 2024 14:59:32 +0100 Subject: [PATCH 5/5] It saves real tracks, not pointers. It saves also the list of hits. Change two default values. --- source/Utils/include/FilterTracks.h | 4 ++-- source/Utils/src/FilterTracks.cc | 25 ++++++++++++++++++++----- 2 files changed, 22 insertions(+), 7 deletions(-) diff --git a/source/Utils/include/FilterTracks.h b/source/Utils/include/FilterTracks.h index 5d55b7c..e11fa65 100644 --- a/source/Utils/include/FilterTracks.h +++ b/source/Utils/include/FilterTracks.h @@ -77,7 +77,7 @@ class FilterTracks : public marlin::Processor //! Cut off for number of hits in outer tracker (barrel and endcap combined) int _NHitsOuter = 1; //! Cut off for number of holes - int _MaxHoles = 1; + int _MaxHoles = 10; //! Cut off for momentum (GeV) float _MinPt = 0.5; //units GeV @@ -91,7 +91,7 @@ class FilterTracks : public marlin::Processor int _MinNdf = 1; //! Cut off for outliers number - int _MaxOutl = 0; + int _MaxOutl = 10; //! Cut off for spatial and temporal chi squared values float _Chi2Spatial = 0; diff --git a/source/Utils/src/FilterTracks.cc b/source/Utils/src/FilterTracks.cc index c014d24..811599a 100644 --- a/source/Utils/src/FilterTracks.cc +++ b/source/Utils/src/FilterTracks.cc @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -180,7 +181,13 @@ void FilterTracks::processEvent( LCEvent * evt ) { // Make the output track collection LCCollectionVec *OutputTrackCollection = new LCCollectionVec(LCIO::TRACK); - OutputTrackCollection->setSubset(true); + // do not store pointers but real tracks + OutputTrackCollection->setSubset(false); + + // if we want to point back to the hits we need to set the flag + LCFlagImpl trkFlag(0); + trkFlag.setBit(LCIO::TRBIT_HITS); + OutputTrackCollection->setFlag(trkFlag.getFlag()); TMVA::Reader* reader = new TMVA::Reader(); @@ -242,12 +249,17 @@ void FilterTracks::processEvent( LCEvent * evt ) break; } } - if(endcaphits == false) { OutputTrackCollection->addElement(trk); } + if (endcaphits == false) { + auto itrk = dynamic_cast(trk); + OutputTrackCollection->addElement(new IMPL::TrackImpl(*itrk)); + } } else { // track property cuts if ( not _NNmethod.empty() ) { // NN cuts auto prediction = reader->EvaluateMVA(_NNmethod); - if ( prediction > _NNthr ) - OutputTrackCollection->addElement(trk); + if ( prediction > _NNthr ) { + auto itrk = dynamic_cast(trk); + OutputTrackCollection->addElement(new IMPL::TrackImpl(*itrk)); + } } else { // user cuts if (vars["trtnh"] > _NHitsTotal && vars["trtvhn"] > _NHitsVertex && @@ -261,7 +273,10 @@ void FilterTracks::processEvent( LCEvent * evt ) vars["trndf"] > _MinNdf && vars["trtnh"]-vars["trndf"]/2 < _MaxOutl && vars["trthn"] < _MaxHoles) - { OutputTrackCollection->addElement(trk); } + { + auto itrk = dynamic_cast(trk); + OutputTrackCollection->addElement(new IMPL::TrackImpl(*itrk)); + } } } }