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} ) diff --git a/source/Utils/include/FilterJetConeHits.h b/source/Utils/include/FilterJetConeHits.h new file mode 100644 index 0000000..ad39914 --- /dev/null +++ b/source/Utils/include/FilterJetConeHits.h @@ -0,0 +1,104 @@ +#ifndef FilterJetConeHits_h +#define FilterJetConeHits_h 1 + +#include "marlin/Processor.h" +#include "lcio.h" +#include +#include + +#include + +#include +#include +#include "IMPL/ReconstructedParticleImpl.h" + +using namespace lcio ; +using namespace marlin ; + + +/** Utility processor that selects and saves the tracker hits that are included in + * a DeltaR cone around the jet direction direction along with the corresponding + * sim hits and the reco-sim relations. + * + * @parameter MCParticleCollection name of the MCParticle collection + * @parameter TrackerHitInputCollections name of the tracker hit input collections + * @parameter TrackerSimHitInputCollections name of the tracker simhit input collections + * @parameter TrackerHitInputRelations name of the tracker hit relation input collections + * @parameter TrackerHitOutputCollections name of the tracker hit output collections + * @parameter TrackerSimHitOutputCollections name of the tracker simhit output collections + * @parameter TrackerHitOutputRelations name of the tracker hit relation output collections + * @parameter DeltaRCut maximum angular distance between the hits and the particle direction + * @parameter FillHistograms flag to fill the diagnostic histograms + * + * @author L. Sestini, INFN Padova + * @date 18 March 2021 + * @version $Id: $ + */ + +class FilterJetConeHits : public Processor { + + public: + + virtual Processor* newProcessor() { return new FilterJetConeHits ; } + + + FilterJetConeHits() ; + + virtual void init() ; + + virtual void processRunHeader( LCRunHeader* run ) ; + + virtual void processEvent( LCEvent * evt ) ; + + virtual void check( LCEvent * evt ) ; + + virtual void end() ; + + bool filterJetBib(ReconstructedParticle* jet); + + void saveJet( ReconstructedParticle* jet, LCCollectionVec* jetsColl ); + + void directionCorrection( const double* p, double* pcorr ) ; + + protected: + + // --- Input/output collection names: + std::string m_inputJetCaloCollName{} ; + std::string m_filteredJetCaloCollName{} ; + std::vector m_inputTrackerHitsCollNames{} ; + std::vector m_inputTrackerSimHitsCollNames{} ; + std::vector m_inputTrackerHitRelNames{} ; + std::vector m_outputTrackerHitsCollNames{} ; + std::vector m_outputTrackerSimHitsCollNames{} ; + std::vector m_outputTrackerHitRelNames{} ; + + + // --- Processor parameters: + bool m_fillHistos{} ; + double m_deltaRCut{} ; + + // Jet filter parameters with BIB: + double m_minDaughterMaxPt{} ; + int m_minNTracks{} ; + bool m_createFilteredJets{} ; + + // Jet direction correction params + bool m_makeDirCorrection{} ; + double m_corrConst{}; + double m_corrLin{}; + double m_corrQuad{}; + double m_corrCub{}; + + // --- Diagnostic histograms: + TH1F* m_dist = nullptr ; + + // --- Run and event counters: + int _nRun{} ; + int _nEvt{} ; + +} ; + +#endif + + + diff --git a/source/Utils/include/FilterTracks.h b/source/Utils/include/FilterTracks.h index 3e5cfc9..50ef28e 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 @@ -70,14 +76,47 @@ 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 = 10; //! 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 theta (rad) + float _MinTheta = 0; + float _MaxTheta = 3.14; + + //! Cut off for the value ndf + int _MinNdf = 1; + float _MaxReducedChi2 = 1.e20; + + //! Cut off for outliers number + int _MaxOutl = 10; + + //! Cut off for ratio outliers/tot hits + float _MaxOutlOverHits = 1.; //! Cut off for spatial and temporal chi squared values float _Chi2Spatial = 0; + //! Cut off for all track parameters uncertainties + float _MaxSigD0 = 999; + float _MaxSigZ0 = 999; + float _MaxSigTheta = 999; + float _MaxSigPhi = 999; + float _MaxSigQoverP = 999; + + float _MaxD0 = 999.; + float _MaxZ0 = 999.; + + //! 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/FilterJetConeHits.cc b/source/Utils/src/FilterJetConeHits.cc new file mode 100644 index 0000000..8c0b9db --- /dev/null +++ b/source/Utils/src/FilterJetConeHits.cc @@ -0,0 +1,561 @@ +#include "FilterJetConeHits.h" +#include +#include +#include + + +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include "HelixClass_double.h" + +using namespace lcio ; +using namespace marlin ; + + +FilterJetConeHits aFilterJetConeHits ; + + +FilterJetConeHits::FilterJetConeHits() : Processor("FilterJetConeHits") { + + // --- Processor description: + + _description = "FilterJetConeHits selects tracker hits in a cone opened around a jet cone direction"; + + + // --- Processor parameters: + + registerProcessorParameter("JetCaloCollection", + "Name of the JetCalo collection", + m_inputJetCaloCollName, + std::string("JetCalo") ); + + registerProcessorParameter("TrackerHitInputCollections", + "Name of the tracker hit input collections", + m_inputTrackerHitsCollNames, + {} ); + + registerProcessorParameter("TrackerSimHitInputCollections", + "Name of the tracker simhit input collections", + m_inputTrackerSimHitsCollNames, + {} ); + + registerProcessorParameter("TrackerHitInputRelations", + "Name of the tracker hit relation collections", + m_inputTrackerHitRelNames, + {} ); + + registerProcessorParameter("TrackerHitOutputCollections", + "Name of the tracker hit output collections", + m_outputTrackerHitsCollNames, + {} ); + + registerProcessorParameter("TrackerSimHitOutputCollections", + "Name of the tracker simhit output collections", + m_outputTrackerSimHitsCollNames, + {} ); + + registerProcessorParameter("TrackerHitOutputRelations", + "Name of the tracker hit relation collections", + m_outputTrackerHitRelNames, + {} ); + + registerProcessorParameter( "DeltaRCut" , + "Maximum angular distance between the hits and the particle direction" , + m_deltaRCut, + double(1.) ); + + registerProcessorParameter( "FillHistograms", + "Flag to fill the diagnostic histograms", + m_fillHistos, + false ); + + registerProcessorParameter( "MinJetTracks", + "Min number of tracks in jet to use it as filter", + m_minNTracks, + int(1) ); + + registerProcessorParameter( "MinDaughterMaxPt", + "Min pT of the highest-pT track in jet to use it as filter", + m_minDaughterMaxPt, + double(2.) ); + + registerProcessorParameter( "FilteredJetsCollectionName", + "Name of the (optional) output filtered jets collection", + m_filteredJetCaloCollName, + std::string("FilteredJets") ); + + registerProcessorParameter( "MakeFilteredJetsCollection", + "Flag for producing collection of filtered jets", + m_createFilteredJets, + bool(false) ); + + registerProcessorParameter( "MakeDirectionCorrection", + "Flag for correcting direction of filtered jets", + m_makeDirCorrection, + bool(false) ); + + registerProcessorParameter( "CorrectionConstant", + "Direction correction constant term", + m_corrConst, + double(0.) ); + + registerProcessorParameter( "CorrectionLinear", + "Direction correction linear term", + m_corrLin, + double(1.) ); + + registerProcessorParameter( "CorrectionQuadratic", + "Direction correction quadratic term", + m_corrQuad, + double(0.) ); + + registerProcessorParameter( "CorrectionCubic", + "Direction correction cubic term", + m_corrCub, + double(0.) ); +} + + + +void FilterJetConeHits::init() { + + streamlog_out(DEBUG) << " init called " << std::endl ; + + // --- Print the processor parameters: + + printParameters() ; + + + // --- Initialize the run and event counters: + + _nRun = 0 ; + _nEvt = 0 ; + + + // --- Initialize the AIDAProcessor and book the diagnostic histograms: + + AIDAProcessor::histogramFactory(this); + + m_dist = new TH1F("m_dist", "deltaR distance between hit and jet axis", 1000, 0., 6.); + +} + + +void FilterJetConeHits::processRunHeader( LCRunHeader* run) { + + _nRun++ ; + +} + + + +void FilterJetConeHits::processEvent( LCEvent * evt ) { + + // --- Check whether the number of input and output collections match + bool noSimColl = false, noRelColl = false; + if ( m_inputTrackerSimHitsCollNames.size() == 0 ) noSimColl = true; + if ( m_inputTrackerHitRelNames.size() == 0 ) noRelColl = true; + if ( noSimColl != noRelColl ) { + std::stringstream err_msg; + err_msg << "Error: RelColl are necessary to save SimColl!" + << std::endl ; + + throw EVENT::Exception( err_msg.str() ) ; + } + + if ( m_inputTrackerHitsCollNames.size() != m_outputTrackerHitsCollNames.size() || + ( !noSimColl && m_inputTrackerSimHitsCollNames.size() != m_outputTrackerSimHitsCollNames.size() ) || + ( !noRelColl && m_inputTrackerHitRelNames.size() != m_outputTrackerHitRelNames.size() ) ) { + + std::stringstream err_msg; + err_msg << "Mismatch between input and output collections" + << std::endl ; + + throw EVENT::Exception( err_msg.str() ) ; + } + + if ( !noSimColl && ( m_inputTrackerHitsCollNames.size() != m_inputTrackerSimHitsCollNames.size() || + m_inputTrackerHitsCollNames.size() != m_inputTrackerHitRelNames.size() ) ){ + + std::stringstream err_msg; + err_msg << "Mismatch between the reco and sim hits input collections" + << std::endl ; + + throw EVENT::Exception( err_msg.str() ) ; + + } + + if ( !noSimColl && ( m_outputTrackerHitsCollNames.size() != m_outputTrackerSimHitsCollNames.size() || + m_outputTrackerHitsCollNames.size() != m_outputTrackerHitRelNames.size() ) ){ + + std::stringstream err_msg; + err_msg << "Mismatch between the reco and sim hits output collections" + << std::endl ; + + throw EVENT::Exception( err_msg.str() ) ; + + } + + + // --- Get the JetCalo collection: + + LCCollection* m_inputJetCalo = NULL; + try { + m_inputJetCalo = evt->getCollection( m_inputJetCaloCollName ); + } + catch( lcio::DataNotAvailableException& e ) { + streamlog_out(WARNING) << m_inputJetCaloCollName << " collection not available" << std::endl; + return; + } + + + // --- Get the input hit collections and create the corresponding output collections: + + const unsigned int nTrackerHitCol = m_inputTrackerHitsCollNames.size(); + std::vector inputHitColls(nTrackerHitCol); + std::vector inputSimHitColls(nTrackerHitCol); + std::vector inputHitRels(nTrackerHitCol); + + std::vector outputTrackerHitColls(nTrackerHitCol); + std::vector outputTrackerSimHitColls(nTrackerHitCol); + std::vector outputTrackerHitRels(nTrackerHitCol); + + for (unsigned int icol=0; icolgetCollection(m_inputTrackerHitsCollNames[icol]); + } + catch( lcio::DataNotAvailableException& e ) { + streamlog_out(WARNING) << m_inputTrackerHitsCollNames[icol] + << " collection not available" << std::endl; + continue; + } + + // get the sim hits + if(!noSimColl){ + + try { + inputSimHitColls[icol] = evt->getCollection(m_inputTrackerSimHitsCollNames[icol]); + } + catch( lcio::DataNotAvailableException& e ) { + streamlog_out(WARNING) << m_inputTrackerSimHitsCollNames[icol] + << " collection not available" << std::endl; + continue; + } + + try { + inputHitRels[icol] = evt->getCollection(m_inputTrackerHitRelNames[icol]); + } + catch( lcio::DataNotAvailableException& e ) { + streamlog_out(WARNING) << m_inputTrackerHitRelNames[icol] + << " collection not available" << std::endl; + continue; + } + } + + // reco hit output collections + std::string encoderString = inputHitColls[icol]->getParameters().getStringVal( "CellIDEncoding" ); + outputTrackerHitColls[icol] = new LCCollectionVec( inputHitColls[icol]->getTypeName() ); + outputTrackerHitColls[icol]->parameters().setValue( "CellIDEncoding", encoderString ); + LCFlagImpl lcFlag(inputHitColls[icol]->getFlag()); + outputTrackerHitColls[icol]->setFlag(lcFlag.getFlag()); + + // sim hit output collections + if(!noSimColl){ + outputTrackerSimHitColls[icol] = new LCCollectionVec( inputSimHitColls[icol]->getTypeName() ); + outputTrackerSimHitColls[icol]->parameters().setValue( "CellIDEncoding", encoderString ); + LCFlagImpl lcFlag_sim(inputSimHitColls[icol]->getFlag()); + outputTrackerSimHitColls[icol]->setFlag(lcFlag_sim.getFlag()); + + outputTrackerHitRels[icol] = new LCCollectionVec( inputHitRels[icol]->getTypeName() ); + LCFlagImpl lcFlag_rel(inputHitRels[icol]->getFlag()) ; + outputTrackerHitRels[icol]->setFlag( lcFlag_rel.getFlag() ) ; + } + + } + + + // --- Loop over the JetCalo: + + std::vector > hits_to_save(nTrackerHitCol); + + LCCollectionVec* filterJetsColl = new LCCollectionVec( LCIO::RECONSTRUCTEDPARTICLE ); + LCFlagImpl flagJets( m_inputJetCalo->getFlag() ); + filterJetsColl->setFlag( flagJets.getFlag() ); + + unsigned int nGoodJets = 0; + + for (int ipart=0; ipartgetNumberOfElements(); ++ipart){ + + ReconstructedParticle* part = dynamic_cast( m_inputJetCalo->getElementAt(ipart) ); + + if( !FilterJetConeHits::filterJetBib(part) ) continue; //quality selection to reject fake jets + nGoodJets++; + + if( m_createFilteredJets ) saveJet( part , filterJetsColl); //save filtered jet + + double mom[3]; + if( !m_makeDirCorrection ){ + for(int n = 0; n < 3; n++) mom[n] = part->getMomentum()[n]; //> to avoid const cast + } + else{ + directionCorrection( part->getMomentum(), mom ); + } + + // --- Loop over the tracker hits and select hits inside a cone around the jet axis: + + for (unsigned int icol=0; icolgetNumberOfElements(); ++ihit){ + + TrackerHitPlane* hit = dynamic_cast(hit_col->getElementAt(ihit)); + + // --- Skip hits that are in the opposite hemisphere w.r.t. the jet axis: + if ( ( hit->getPosition()[0]*mom[0] + + hit->getPosition()[1]*mom[1] + + hit->getPosition()[2]*mom[2] ) < 0. ) continue; + + + // --- Get the distance between the hit and the jet axis + + double jet_p = sqrt( pow(mom[0],2) + pow(mom[1],2) + pow(mom[2],2) ); + double jet_theta = acos(mom[2]/jet_p); + double jet_eta = -std::log(tan(jet_theta/2.)); + + double hit_d = sqrt( pow(hit->getPosition()[0],2) + pow(hit->getPosition()[1],2) + pow(hit->getPosition()[2],2) ); + double hit_theta = acos(hit->getPosition()[2]/hit_d); + double hit_eta = -std::log(tan(hit_theta/2.)); + + double jet_pxy = sqrt( pow(mom[0],2) + pow(mom[1],2) ); + double hit_dxy = sqrt( pow(hit->getPosition()[0],2) + pow(hit->getPosition()[1],2) ); + double deltaPhi = acos( (mom[0]*hit->getPosition()[0] + mom[1]*hit->getPosition()[1] )/jet_pxy/hit_dxy); + + double deltaR = sqrt( pow(deltaPhi,2) + pow(jet_eta-hit_eta,2) ); + + if ( m_fillHistos ){ + + m_dist->Fill(deltaR); + + } + + if ( deltaR < m_deltaRCut ) + hits_to_save[icol].insert(ihit); + + } // ihit loop + + } // icol loop + + } // ipart loop + + + // --- Add the filtered hits to the output collections: + + for (unsigned int icol=0; icol(inputHitColls[icol]->getElementAt(ihit)); + TrackerHitPlaneImpl* hit_new = new TrackerHitPlaneImpl(); + + hit_new->setCellID0(hit->getCellID0()); + hit_new->setCellID1(hit->getCellID1()); + hit_new->setType(hit->getType()); + hit_new->setPosition(hit->getPosition()); + hit_new->setU(hit->getU()); + hit_new->setV(hit->getV()); + hit_new->setdU(hit->getdU()); + hit_new->setdV(hit->getdV()); + hit_new->setEDep(hit->getEDep()); + hit_new->setEDepError(hit->getEDepError()); + hit_new->setTime(hit->getTime()); + hit_new->setQuality(hit->getQuality()); + + outputTrackerHitColls[icol]->addElement( hit_new ); + + if( noSimColl ) continue; + + LCRelation* rel = dynamic_cast(inputHitRels[icol]->getElementAt(ihit)); + + + SimTrackerHit* simhit = dynamic_cast(rel->getTo()); + SimTrackerHitImpl* simhit_new = new SimTrackerHitImpl(); + + simhit_new->setCellID0(simhit->getCellID0()); + simhit_new->setCellID1(simhit->getCellID1()); + simhit_new->setPosition(simhit->getPosition()); + simhit_new->setEDep(simhit->getEDep()); + simhit_new->setTime(simhit->getTime()); + simhit_new->setMCParticle(simhit->getMCParticle()); + simhit_new->setMomentum(simhit->getMomentum()); + simhit_new->setPathLength(simhit->getPathLength()); + simhit_new->setQuality(simhit->getQuality()); + simhit_new->setOverlay(simhit->isOverlay()); + simhit_new->setProducedBySecondary(simhit->isProducedBySecondary()); + + outputTrackerSimHitColls[icol]->addElement( simhit_new ); + + + LCRelationImpl* rel_new = new LCRelationImpl(); + + rel_new->setFrom(hit_new); + rel_new->setTo(simhit_new); + rel_new->setWeight(rel->getWeight()); + + outputTrackerHitRels[icol]->addElement( rel_new ); + + } // ihit loop + + streamlog_out( MESSAGE ) << " " << hits_to_save[icol].size() << " hits added to the collections: " + << m_outputTrackerHitsCollNames[icol] << std::endl; + + evt->addCollection( outputTrackerHitColls[icol], m_outputTrackerHitsCollNames[icol] ) ; + if(!noSimColl){ + evt->addCollection( outputTrackerSimHitColls[icol], m_outputTrackerSimHitsCollNames[icol] ) ; + evt->addCollection( outputTrackerHitRels[icol], m_outputTrackerHitRelNames[icol] ) ; + } + + streamlog_out( DEBUG5 ) << " output collection " << m_outputTrackerHitsCollNames[icol] << " of type " + << outputTrackerHitColls[icol]->getTypeName() << " added to the event" << std::endl; + if(!noSimColl) streamlog_out( DEBUG5 ) << " output collection " << m_outputTrackerSimHitsCollNames[icol] << " of type " + << outputTrackerSimHitColls[icol]->getTypeName() << " added to the event \n" + << " output collection " << m_outputTrackerHitRelNames[icol] << " of type " + << outputTrackerHitRels[icol]->getTypeName() << " added to the event " + << std::endl ; + + } // icol loop + + if( m_createFilteredJets ){ + evt->addCollection(filterJetsColl, m_filteredJetCaloCollName); + streamlog_out(MESSAGE) << "Saved Filtered Jets collection, with " << filterJetsColl->getNumberOfElements() << " jets." << std::endl; + } + else{ + delete filterJetsColl; + } + + streamlog_out(DEBUG) << " processing event: " << evt->getEventNumber() + << " in run: " << evt->getRunNumber() << std::endl ; + + _nEvt ++ ; + +} + + + +void FilterJetConeHits::check( LCEvent * evt ) { +} + + + +void FilterJetConeHits::end(){ + + std::cout << "FilterJetConeHits::end() " << name() + << " processed " << _nEvt << " events in " << _nRun << " runs " + << std::endl ; + +} + +// jet quality selection +bool FilterJetConeHits::filterJetBib(ReconstructedParticle* jet){ + + auto daughters = jet->getParticles(); + + double maxPt = 0., pt = 0.; + int ntracks = 0; + + for(unsigned int i = 0; i < daughters.size(); i++){ + + if(daughters[i]->getCharge() == 0) continue; //only care about tracks + ntracks++; + pt = sqrt( pow( daughters[i]->getMomentum()[0] ,2) + pow( daughters[i]->getMomentum()[1] ,2) ); + if(pt > maxPt) maxPt = pt; + } + + if(ntracks < m_minNTracks || maxPt < m_minDaughterMaxPt) return false; //check quality of jet, if bad reject + + streamlog_out(DEBUG) << "Accepted jet with ntracks = " << ntracks << " and max daughter pT = " << maxPt << std::endl; + return true; +} + +// save jet in collection +void FilterJetConeHits::saveJet( ReconstructedParticle* jet, LCCollectionVec* jetsColl ){ + + // this way only a pointer to the original jet is saved + // if necessary implement here the phys copy save + //jetsColl->addElement(jet); + ReconstructedParticleImpl* j = new ReconstructedParticleImpl(); + j->setType( jet->getType() ); + + double pcorr[3]; + if( !m_makeDirCorrection ) j->setMomentum( jet->getMomentum() ); + else{ + directionCorrection(jet->getMomentum(), pcorr); + j->setMomentum( pcorr ); + } + + j->setEnergy( jet->getEnergy() ); + j->setMass( jet->getMass() ); + j->setCharge( jet->getCharge() ); + j->setReferencePoint( jet->getReferencePoint() ); + + const EVENT::ReconstructedParticleVec rpvec = jet->getParticles(); + for(unsigned int i = 0; i < rpvec.size(); i++) j->addParticle(rpvec[i]); + + const EVENT::ClusterVec clvec = jet->getClusters(); + for(unsigned int i = 0; i < clvec.size(); i++) j->addCluster(clvec[i]); + + const EVENT::TrackVec trvec = jet->getTracks(); + for(unsigned int i = 0; i < trvec.size(); i++) j->addTrack(trvec[i]); + + jetsColl->addElement(j); + + streamlog_out( MESSAGE ) << "Saving Jet: p = ( " << j->getMomentum()[0] << " , " + << j->getMomentum()[1] << " , " + << j->getMomentum()[2] << " ) " << std::endl; + + if( m_makeDirCorrection ){ + const double* pold = jet->getMomentum(); + double oldtheta = acos( pold[2] / sqrt(pold[0]*pold[0] + pold[1]*pold[1] + pold[2]*pold[2]) ) * 180. / 3.14159265; + const double* pnew = j->getMomentum(); + double newtheta = acos( pnew[2] / sqrt(pnew[0]*pnew[0] + pnew[1]*pnew[1] + pnew[2]*pnew[2]) ) * 180. / 3.14159265; + streamlog_out( MESSAGE ) << "Corrected: old THETA = " << oldtheta << " => new THETA = " << newtheta << std::endl; + } + + return; +} + +void FilterJetConeHits::directionCorrection( const double* p, double* pcorr ){ + + double ptot = sqrt( p[0]*p[0] + p[1]*p[1] + p[2]*p[2] ); + double theta = acos( p[2] / ptot ) * 180. / 3.14159265; + double theta_pre = theta; + + //> only in transition region, hardcoded for now... + if(theta > 30. && theta < 60.){ + theta = m_corrConst + m_corrLin * theta + m_corrQuad * pow(theta, 2.) + m_corrCub * pow(theta, 3.); + } + else if(theta > 120. && theta < 150.){ + theta = 180. - theta; + theta = m_corrConst + m_corrLin * theta + m_corrQuad * pow(theta, 2.) + m_corrCub * pow(theta, 3.); + theta = 180. - theta; + } + + pcorr[0] = p[0] * sin( theta * 3.14159265 / 180. ) / sin( theta_pre * 3.14159265 / 180. ); + pcorr[1] = p[1] * sin( theta * 3.14159265 / 180. ) / sin( theta_pre * 3.14159265 / 180. ); + pcorr[2] = p[2] * cos( theta * 3.14159265 / 180. ) / cos( theta_pre * 3.14159265 / 180. ); + + return; +} + \ No newline at end of file diff --git a/source/Utils/src/FilterTracks.cc b/source/Utils/src/FilterTracks.cc index 4e28536..c4ff1c4 100644 --- a/source/Utils/src/FilterTracks.cc +++ b/source/Utils/src/FilterTracks.cc @@ -7,11 +7,18 @@ #include #include #include +#include +#include #include #include +#include "TMVA/Reader.h" + +#include + FilterTracks aFilterTracks ; + FilterTracks::FilterTracks() : Processor("FilterTracks") { @@ -55,12 +62,133 @@ FilterTracks::FilterTracks() _MinPt ); + registerProcessorParameter("MaxPt", + "Max transverse momentum", + _MaxPt, + _MaxPt + ); + + registerProcessorParameter("MinTheta", + "Minimum theta", + _MinTheta, + _MinTheta + ); + + registerProcessorParameter("MaxTheta", + "Max theta", + _MaxTheta, + _MaxTheta + ); + 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("MaxOutliers", + "Max number of outliers", + _MaxOutl, + _MaxOutl + ); + + registerProcessorParameter("MaxOutliersOverHits", + "Max ratio of outliers/hits", + _MaxOutlOverHits, + _MaxOutlOverHits + ); + + registerProcessorParameter("MaxSigmaD0", + "Max sigma d0", + _MaxSigD0, + _MaxSigD0 + ); + + registerProcessorParameter("MaxSigmaZ0", + "Max sigma z0", + _MaxSigZ0, + _MaxSigZ0 + ); + + registerProcessorParameter("MaxSigmaTheta", + "Max sigma theta", + _MaxSigTheta, + _MaxSigTheta + ); + + registerProcessorParameter("MaxSigmaPhi", + "Max sigma phi", + _MaxSigPhi, + _MaxSigPhi + ); + + registerProcessorParameter("MaxSigmaQoverP", + "Max sigma q/p", + _MaxSigQoverP, + _MaxSigQoverP + ); + + registerProcessorParameter("MaxReducedChi2", + "Max reduced Chi2", + _MaxReducedChi2, + _MaxReducedChi2 + ); + + registerProcessorParameter("MaxPromptD0", + "Max d0", + _MaxD0, + _MaxD0 + ); + + registerProcessorParameter("MaxPromptZ0", + "Max z0", + _MaxZ0, + _MaxZ0 + ); + + 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 +209,8 @@ void FilterTracks::init() { // Print the initial parameters printParameters() ; - buildBfield() ; + // it requires that geometry has been already initialized + // buildBfield() ; } void FilterTracks::processRunHeader( LCRunHeader* /*run*/) @@ -106,11 +235,19 @@ 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); - // Get input collection - LCCollection* InputTrackCollection =evt->getCollection(_InputTrackCollection); + // 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(); + + // Get input collection + LCCollection* InputTrackCollection = evt->getCollection(_InputTrackCollection); + if( InputTrackCollection->getTypeName() != lcio::LCIO::TRACK ) { throw EVENT::Exception( "Invalid collection type: " + InputTrackCollection->getTypeName() ) ; } @@ -118,45 +255,117 @@ 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}, + {"trnoh", 0}, + {"tr_sigd0",0}, + {"tr_sigz0",0}, + {"tr_sigtheta",0}, + {"tr_sigphi",0}, + {"tr_sigqoverp",0}, + {"tr_redchi2",0}, + {"tr_d0",0}, + {"tr_z0",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 theta = M_PI_2-atan(trk->getTanLambda()); + vars["trch2"] = trk->getChi2(); - float chi2spatial = trk->getChi2(); + vars["trndf"] = trk->getNdf(); + vars["tr_redchi2"] = vars["trch2"] / vars["trndf"]; + vars["trnoh"] = (vars["trtnh"]-vars["trndf"]/2) / vars["trtnh"]; + vars["tr_sigd0"] = sqrt(trk->getCovMatrix()[0]); + vars["tr_sigz0"] = sqrt(trk->getCovMatrix()[2]); + vars["tr_sigtheta"] = sqrt(trk->getCovMatrix()[9]); + vars["tr_sigphi"] = sqrt(trk->getCovMatrix()[5]); + vars["tr_sigqoverp"] = sqrt(trk->getCovMatrix()[14]); + + vars["tr_d0"] = trk->getD0(); + vars["tr_z0"] = trk->getZ0(); 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) { + auto itrk = dynamic_cast(trk); + OutputTrackCollection->addElement(new IMPL::TrackImpl(*itrk)); } - 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 (!(vars["trtnh"] > _NHitsTotal && + vars["trtvhn"] > _NHitsVertex && + vars["trtihn"] > _NHitsInner && + vars["trtohn"] > _NHitsOuter && + pt > _MinPt && + pt < _MaxPt && + theta > _MinTheta && + theta < _MaxTheta && + vars["trch2"] > _Chi2Spatial && + vars["trndf"] > _MinNdf && + vars["tr_redchi2"] < _MaxReducedChi2 && + vars["trtnh"]-vars["trndf"]/2 < _MaxOutl && + vars["trnoh"] < _MaxOutlOverHits && + vars["trthn"] < _MaxHoles && + vars["tr_sigd0"] < _MaxSigD0 && + vars["tr_sigz0"] < _MaxSigZ0 && + vars["tr_sigtheta"] < _MaxSigTheta && + vars["tr_sigphi"] < _MaxSigPhi && + vars["tr_sigqoverp"] < _MaxSigQoverP && + std::fabs(vars["tr_d0"]) < _MaxD0 && + std::fabs(vars["tr_z0"]) < _MaxZ0)) continue; + + if ( not _NNmethod.empty() ) { // NN cuts + + auto prediction = reader->EvaluateMVA(_NNmethod); + if ( prediction <= _NNthr ) continue; + } + + auto itrk = dynamic_cast(trk); + OutputTrackCollection->addElement(new IMPL::TrackImpl(*itrk)); + } } // Save output track collection evt->addCollection(OutputTrackCollection, _OutputTrackCollection); + delete reader; } void FilterTracks::end() -{ } +{ }