Skip to content

Music10 tev #23

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
pandreetto merged 4 commits into from
Apr 3, 2025
Merged

Music10 tev #23

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
73e7f5a
LP: added new filter param in TrackFilter
LeonardoPalombini Feb 13, 2025
accf181
LP: added FakeRemovedJets storing
LeonardoPalombini Feb 21, 2025
13d0d09
LP: add new uncertainty-based track filters
LeonardoPalombini Feb 24, 2025
97ff149
LP: added jet direction correction in FilterHitsJets
LeonardoPalombini Feb 26, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 12 additions & 0 deletions source/Utils/include/FilterJetConeHits.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,12 +55,16 @@ class FilterJetConeHits : public Processor {
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<std::string> m_inputTrackerHitsCollNames{} ;
std::vector<std::string> m_inputTrackerSimHitsCollNames{} ;
std::vector<std::string> m_inputTrackerHitRelNames{} ;
Expand All @@ -76,6 +80,14 @@ class FilterJetConeHits : public Processor {
// 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 ;
Expand Down
10 changes: 10 additions & 0 deletions source/Utils/include/FilterTracks.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -93,9 +93,19 @@ class FilterTracks : public marlin::Processor
//! 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;

//! NN parameters
std::string _NNmethod = ""; // if defined apply the NN
std::string _NNweights = ""; // xml file with weights
Expand Down
179 changes: 153 additions & 26 deletions source/Utils/src/FilterJetConeHits.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,9 +86,42 @@ FilterJetConeHits::FilterJetConeHits() : Processor("FilterJetConeHits") {
registerProcessorParameter( "MinDaughterMaxPt",
"Min pT of the highest-pT track in jet to use it as filter",
m_minDaughterMaxPt,
double(2.) );


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.) );
}


Expand Down Expand Up @@ -255,14 +288,29 @@ void FilterJetConeHits::processEvent( LCEvent * evt ) {
// --- Loop over the JetCalo:

std::vector<std::set<int> > 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; ipart<m_inputJetCalo->getNumberOfElements(); ++ipart){

ReconstructedParticle* part = dynamic_cast<ReconstructedParticle*>( m_inputJetCalo->getElementAt(ipart) );

if( !FilterJetConeHits::filterJetBib(part) ) continue; //quality selection to reject fake jets
nGoodJets++;

// double part_p = sqrt( part->getMomentum()[0]*part->getMomentum()[0] +part->getMomentum()[1]*part->getMomentum()[1] + part->getMomentum()[2]*part->getMomentum()[2] );
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:

Expand All @@ -273,38 +321,38 @@ void FilterJetConeHits::processEvent( LCEvent * evt ) {

for (int ihit=0; ihit<hit_col->getNumberOfElements(); ++ihit){

TrackerHitPlane* hit = dynamic_cast<TrackerHitPlane*>(hit_col->getElementAt(ihit));
TrackerHitPlane* hit = dynamic_cast<TrackerHitPlane*>(hit_col->getElementAt(ihit));

// --- Skip hits that are in the opposite hemisphere w.r.t. the jet axis:
if ( ( hit->getPosition()[0]*part->getMomentum()[0] +
hit->getPosition()[1]*part->getMomentum()[1] +
hit->getPosition()[2]*part->getMomentum()[2] ) < 0. ) continue;
// --- 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(part->getMomentum()[0],2) + pow(part->getMomentum()[1],2) + pow(part->getMomentum()[2],2) );
double jet_theta = acos(part->getMomentum()[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);
// --- 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(part->getMomentum()[0],2) + pow(part->getMomentum()[1],2) );
double hit_dxy = sqrt( pow(hit->getPosition()[0],2) + pow(hit->getPosition()[1],2) );
double deltaPhi = acos( (part->getMomentum()[0]*hit->getPosition()[0] + part->getMomentum()[1]*hit->getPosition()[1] )/jet_pxy/hit_dxy);
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) );
double deltaR = sqrt( pow(deltaPhi,2) + pow(jet_eta-hit_eta,2) );

if ( m_fillHistos ){
if ( m_fillHistos ){

m_dist->Fill(deltaR);
m_dist->Fill(deltaR);

}
}

if ( deltaR < m_deltaRCut )
hits_to_save[icol].insert(ihit);
if ( deltaR < m_deltaRCut )
hits_to_save[icol].insert(ihit);

} // ihit loop

Expand Down Expand Up @@ -389,6 +437,14 @@ void FilterJetConeHits::processEvent( LCEvent * evt ) {

} // 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 ;

Expand Down Expand Up @@ -432,3 +488,74 @@ bool FilterJetConeHits::filterJetBib(ReconstructedParticle* jet){
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;
}

55 changes: 54 additions & 1 deletion source/Utils/src/FilterTracks.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -104,6 +104,42 @@ FilterTracks::FilterTracks()
_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("Bz",
"Magnetic field in Tesla (default: 5)",
_Bz,
Expand Down Expand Up @@ -209,6 +245,11 @@ void FilterTracks::processEvent( LCEvent * evt )
{"trtnh", 0},
{"trch2", 0},
{"trndf", 0},
{"tr_sigd0",0},
{"tr_sigz0",0},
{"tr_sigtheta",0},
{"tr_sigphi",0},
{"tr_sigqoverp",0}
};

if ( not _NNmethod.empty() ) {
Expand Down Expand Up @@ -239,6 +280,12 @@ void FilterTracks::processEvent( LCEvent * evt )

vars["trndf"] = trk->getNdf();

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]);

if(_BarrelOnly == true) {
bool endcaphits = false;
for(int j=0; j<vars["trtnh"]; ++j) {
Expand Down Expand Up @@ -272,7 +319,13 @@ void FilterTracks::processEvent( LCEvent * evt )
vars["trch2"] > _Chi2Spatial &&
vars["trndf"] > _MinNdf &&
vars["trtnh"]-vars["trndf"]/2 < _MaxOutl &&
vars["trthn"] < _MaxHoles)
(vars["trtnh"]-vars["trndf"]/2) / vars["trtnh"] < _MaxOutlOverHits &&
vars["trthn"] < _MaxHoles &&
vars["tr_sigd0"] < _MaxSigD0 &&
vars["tr_sigz0"] < _MaxSigZ0 &&
vars["tr_sigtheta"] < _MaxSigTheta &&
vars["tr_sigphi"] < _MaxSigPhi &&
vars["tr_sigqoverp"] < _MaxSigQoverP)
{
auto itrk = dynamic_cast<IMPL::TrackImpl*>(trk);
OutputTrackCollection->addElement(new IMPL::TrackImpl(*itrk));
Expand Down
Loading