Skip to content
Snippets Groups Projects
Commit d63bd54f authored by Davide Zuliani's avatar Davide Zuliani :bow_tone1: Committed by Mark Smith
Browse files

Jets daughters

parent f0343b79
No related branches found
No related tags found
3 merge requests!1039Fixed formatting,!1000Fixed formatting,!905Jets daughters
/*****************************************************************************\
* (c) Copyright 2000-2019 CERN for the benefit of the LHCb Collaboration *
* (c) Copyright 2000-2022 CERN for the benefit of the LHCb Collaboration *
* *
* This software is distributed under the terms of the GNU General Public *
* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
......@@ -29,6 +29,7 @@ TupleToolJets::TupleToolJets( const std::string& type, const std::string& name,
{
declareInterface<IParticleTupleTool>( this );
declareProperty( "withJetConstituents", m_withJetConstituents = false, "Save jets constituents?" );
}
//=============================================================================
......
/*****************************************************************************\
* (c) Copyright 2000-2019 CERN for the benefit of the LHCb Collaboration *
* (c) Copyright 2000-2022 CERN for the benefit of the LHCb Collaboration *
* *
* This software is distributed under the terms of the GNU General Public *
* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
......@@ -28,6 +28,7 @@ public:
private:
const LHCb::Particle* m_p;
std::string m_head;
bool m_withJetConstituents;
};
#endif // TUPLETOOLJETS_H
/*****************************************************************************\
* (c) Copyright 2000-2019 CERN for the benefit of the LHCb Collaboration *
* (c) Copyright 2000-2022 CERN for the benefit of the LHCb Collaboration *
* *
* This software is distributed under the terms of the GNU General Public *
* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
......@@ -9,8 +9,11 @@
* or submit itself to any jurisdiction. *
\*****************************************************************************/
#include "TupleToolJetsBase.h"
#include "Kernel/IDVAlgorithm.h"
#include "Kernel/IDistanceCalculator.h"
#include "LoKi/AParticleCuts.h"
#include "LoKi/ParticleCuts.h"
#include <Kernel/GetIDVAlgorithm.h>
TupleToolJetsBase::TupleToolJetsBase( const std::string& type, const std::string& name, const IInterface* parent )
: TupleToolBase( type, name, parent )
......@@ -20,7 +23,9 @@ TupleToolJetsBase::TupleToolJetsBase( const std::string& type, const std::string
, neutralParticles( LoKi::Cuts::ONE )
, maxPT( LoKi::Cuts::ONE )
, m_M( LoKi::Cuts::ONE )
, m_MM( LoKi::Cuts::ONE ) {}
, m_MM( LoKi::Cuts::ONE )
, m_dva( 0 )
, m_dist( 0 ) {}
#define SAVEPOINT( POINT, FUN ) ( POINT ? POINT->FUN : -1.0 )
bool TupleToolJetsBase::WriteJetToTuple( const LHCb::Particle* jet, std::string prefix ) {
// filter plus and minus signs out to the prefix
......@@ -181,46 +186,246 @@ bool TupleToolJetsBase::WriteJetToTuple( const LHCb::Particle* jet, std::string
->column( prefix + "_PU_Cone15_NUpVeloTrk",
(double)( jet->info( LHCb::JetPileUpInfo::PUNUpVeloTrkCone15, -999.9 ) ) );
/*
result &= (*m_tuple)->column( prefix+"_ParticleMultiplicity", (double)( SAVEPOINT(jet,daughters().size())));
// result &= (*m_tuple)->column( prefix+"_Charge", (double)( jet?charge(jet):-1.0));
result &= (*m_tuple)->column( prefix+"_positiveParticleMultiplicity", (double)( jet?positiveParticles(jet):-1.0));
result &= (*m_tuple)->column( prefix+"_negativeParticleMultiplicity", (double)( jet?negativeParticles(jet):-1.0));
result &= (*m_tuple)->column( prefix+"_neutralParticleMultiplicity", (double)( jet?neutralParticles(jet):-1.0));
result &= (*m_tuple)->column( prefix+"_chargedParticleMultiplicity", (double)(
jet?positiveParticles(jet)+negativeParticles(jet):-1.0)); result &= (*m_tuple)->column( prefix+"_maxPT", (double)(
jet?maxPT(jet):-1.0)); SmartRefVector< LHCb::Particle > SortedDaughters;
//buffer to store intermediate result to speed things up a bit
result &= (*m_tuple)->column( prefix+"_PT1", (double)( jet?MaxSumNPart(jet,1,
LoKi::Cuts::PT,&SortedDaughters):-1.0)); result &= (*m_tuple)->column( prefix+"_PT2", (double)(
jet?MaxSumNPart(jet,2, LoKi::Cuts::PT,&SortedDaughters):-1.0)); result &= (*m_tuple)->column( prefix+"_PT3",
(double)( jet?MaxSumNPart(jet,3, LoKi::Cuts::PT,&SortedDaughters):-1.0)); result &= (*m_tuple)->column(
prefix+"_PT4", (double)( jet?MaxSumNPart(jet,4, LoKi::Cuts::PT,&SortedDaughters):-1.0)); result &=
(*m_tuple)->column( prefix+"_PT5", (double)( jet?MaxSumNPart(jet,5, LoKi::Cuts::PT,&SortedDaughters):-1.0)); result
&= (*m_tuple)->column( prefix+"_PT10", (double)( jet?MaxSumNPart(jet,10,LoKi::Cuts::PT,&SortedDaughters):-1.0));
result &= (*m_tuple)->column( prefix+"_PT15", (double)(
jet?MaxSumNPart(jet,15,LoKi::Cuts::PT,&SortedDaughters):-1.0)); result &= (*m_tuple)->column( prefix+"_PT25",
(double)( jet?MaxSumNPart(jet,25,LoKi::Cuts::PT,&SortedDaughters):-1.0));
*/
if ( m_withJetConstituents == true ) {
std::vector<const LHCb::Particle*> constituentsvector = jet->daughtersVector();
std::vector<const LHCb::Particle*>::iterator iconstituent = constituentsvector.begin();
std::vector<double> constituents_E;
std::vector<double> constituents_pT;
std::vector<double> constituents_pX;
std::vector<double> constituents_pY;
std::vector<double> constituents_pZ;
std::vector<double> constituents_ID;
std::vector<double> constituents_Eta;
std::vector<double> constituents_Phi;
std::vector<double> constituents_Q;
std::vector<double> constituents_IP;
std::vector<double> constituents_IPCHI2;
std::vector<double> constituents_IPraw;
std::vector<double> constituents_NNe;
std::vector<double> constituents_NNk;
std::vector<double> constituents_NNp;
std::vector<double> constituents_NNpi;
std::vector<double> constituents_NNmu;
std::vector<double> constituents_Chi2;
std::vector<double> constituents_QoverP;
std::vector<double> constituents_trackX;
std::vector<double> constituents_trackY;
std::vector<double> constituents_trackZ;
std::vector<double> constituents_trackVX;
std::vector<double> constituents_trackVY;
std::vector<double> constituents_trackVZ;
std::vector<double> constituents_CaloNeutralEcal;
std::vector<double> constituents_CaloNeutralHcal2Ecal;
std::vector<double> constituents_CaloNeutralE49;
std::vector<double> constituents_CaloNeutralPrs;
const LHCb::VertexBase* BPV = m_dva->bestVertex( *iconstituent );
double ip = -100.0;
int zsign = 0;
double ipC = 0, ipChi2 = 0, ipex = -100;
double c, e;
for ( ; iconstituent != constituentsvector.end(); ++iconstituent ) {
const LHCb::Particle* constituent = *iconstituent;
const LHCb::ProtoParticle* proto = constituent->proto();
if ( charge( constituent ) != 0 ) {
constituents_CaloNeutralEcal.push_back( -1000 );
constituents_CaloNeutralHcal2Ecal.push_back( -1000 );
constituents_CaloNeutralE49.push_back( -1000 );
constituents_CaloNeutralPrs.push_back( -1000 );
if ( proto ) {
constituents_NNe.push_back( proto->info( LHCb::ProtoParticle::ProbNNe, -1000 ) );
constituents_NNk.push_back( proto->info( LHCb::ProtoParticle::ProbNNk, -1000 ) );
constituents_NNp.push_back( proto->info( LHCb::ProtoParticle::ProbNNp, -1000 ) );
constituents_NNpi.push_back( proto->info( LHCb::ProtoParticle::ProbNNpi, -1000 ) );
constituents_NNmu.push_back( proto->info( LHCb::ProtoParticle::ProbNNmu, -1000 ) );
constituents_Chi2.push_back( proto->info( LHCb::ProtoParticle::TrackChi2PerDof, -1000 ) );
const LHCb::Track* trk = proto->track();
LHCb::State sta = trk->firstState();
c = std::abs( sta.qOverP() );
e = sqrt( sta.errQOverP2() );
constituents_QoverP.push_back( c / e );
Gaudi::XYZPoint Pos = trk->position();
Gaudi::XYZVector Mom = trk->momentum();
constituents_trackX.push_back( Pos.X() );
constituents_trackY.push_back( Pos.Y() );
constituents_trackZ.push_back( Pos.Z() );
constituents_trackVX.push_back( Mom.X() );
constituents_trackVY.push_back( Mom.Y() );
constituents_trackVZ.push_back( Mom.Z() );
} else {
constituents_NNe.push_back( -1000 );
constituents_NNk.push_back( -1000 );
constituents_NNp.push_back( -1000 );
constituents_NNpi.push_back( -1000 );
constituents_NNmu.push_back( -1000 );
constituents_Chi2.push_back( -1000 );
constituents_QoverP.push_back( -1000 );
constituents_trackX.push_back( -1000 );
constituents_trackY.push_back( -1000 );
constituents_trackZ.push_back( -1000 );
constituents_trackVX.push_back( -1000 );
constituents_trackVY.push_back( -1000 );
constituents_trackVZ.push_back( -1000 );
}
} else {
constituents_NNe.push_back( -1000 );
constituents_NNk.push_back( -1000 );
constituents_NNp.push_back( -1000 );
constituents_NNpi.push_back( -1000 );
constituents_NNmu.push_back( -1000 );
constituents_Chi2.push_back( -1000 );
constituents_QoverP.push_back( -1000 );
constituents_trackX.push_back( -1000 );
constituents_trackY.push_back( -1000 );
constituents_trackZ.push_back( -1000 );
constituents_trackVX.push_back( -1000 );
constituents_trackVY.push_back( -1000 );
constituents_trackVZ.push_back( -1000 );
if ( proto ) {
constituents_CaloNeutralEcal.push_back( proto->info( LHCb::ProtoParticle::CaloNeutralEcal, -1000 ) );
constituents_CaloNeutralHcal2Ecal.push_back(
proto->info( LHCb::ProtoParticle::CaloNeutralHcal2Ecal, -1000 ) );
constituents_CaloNeutralE49.push_back( proto->info( LHCb::ProtoParticle::CaloNeutralE49, -1000 ) );
constituents_CaloNeutralPrs.push_back( proto->info( LHCb::ProtoParticle::CaloNeutralPrs, -1000 ) );
} else {
constituents_CaloNeutralEcal.push_back( -1000 );
constituents_CaloNeutralHcal2Ecal.push_back( -1000 );
constituents_CaloNeutralE49.push_back( -1000 );
constituents_CaloNeutralPrs.push_back( -1000 );
}
}
constituents_E.push_back( constituent->momentum().E() );
constituents_pT.push_back( constituent->momentum().Pt() );
constituents_pX.push_back( constituent->momentum().Px() );
constituents_pY.push_back( constituent->momentum().Py() );
constituents_pZ.push_back( constituent->momentum().Pz() );
constituents_Eta.push_back( constituent->momentum().Eta() );
constituents_Phi.push_back( constituent->momentum().Phi() );
constituents_Q.push_back( constituent ? charge( constituent ) : -1.0 );
constituents_ID.push_back( constituent->particleID().pid() );
ip = -100.0;
zsign = -999;
ipC = 0;
ipChi2 = 0;
ipex = -100;
StatusCode sc2 = m_dist->distance( constituent, BPV, ipC, ipChi2 );
Gaudi::XYZVector ipV;
StatusCode sc = m_dist->distance( constituent, BPV, ipV );
if ( sc2 && ipChi2 != 0 ) {
if ( sc ) zsign = ipV.z() > 0 ? 1 : -1;
ip = ipC * zsign;
ipex = ip / ipChi2;
}
constituents_IP.push_back( ipex );
constituents_IPCHI2.push_back( ipChi2 );
constituents_IPraw.push_back( ip );
}
result &= ( *m_tuple )
->farray( prefix + "_Constituents_E", constituents_E.begin(), constituents_E.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_pT", constituents_pT.begin(), constituents_pT.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_ID", constituents_ID.begin(), constituents_ID.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_pX", constituents_pX.begin(), constituents_pX.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_pY", constituents_pY.begin(), constituents_pY.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_pZ", constituents_pZ.begin(), constituents_pZ.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_Eta", constituents_Eta.begin(), constituents_Eta.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_Phi", constituents_Phi.begin(), constituents_Phi.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_Q", constituents_Q.begin(), constituents_Q.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_IP", constituents_IP.begin(), constituents_IP.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_IPCHI2", constituents_IPCHI2.begin(), constituents_IPCHI2.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_IPraw", constituents_IPraw.begin(), constituents_IPraw.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_NNe", constituents_NNe.begin(), constituents_NNe.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_NNk", constituents_NNk.begin(), constituents_NNk.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_NNp", constituents_NNp.begin(), constituents_NNp.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_NNpi", constituents_NNpi.begin(), constituents_NNpi.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_NNmu", constituents_NNmu.begin(), constituents_NNmu.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_Chi2", constituents_Chi2.begin(), constituents_Chi2.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_QoverP", constituents_QoverP.begin(), constituents_QoverP.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_Chi2", constituents_Chi2.begin(), constituents_Chi2.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_trackX", constituents_trackX.begin(), constituents_trackX.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_trackY", constituents_trackY.begin(), constituents_trackY.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_trackZ", constituents_trackZ.begin(), constituents_trackZ.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_trackVX", constituents_trackVX.begin(), constituents_trackVX.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_trackVY", constituents_trackVY.begin(), constituents_trackVY.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_trackVZ", constituents_trackVZ.begin(), constituents_trackVZ.end(),
prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_CaloNeutralEcal", constituents_CaloNeutralEcal.begin(),
constituents_CaloNeutralEcal.end(), prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_CaloNeutralHcal2Ecal", constituents_CaloNeutralHcal2Ecal.begin(),
constituents_CaloNeutralHcal2Ecal.end(), prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_CaloNeutralE49", constituents_CaloNeutralE49.begin(),
constituents_CaloNeutralE49.end(), prefix + "_nConstituents", 100 );
result &= ( *m_tuple )
->farray( prefix + "_Constituents_CaloNeutralPrs", constituents_CaloNeutralPrs.begin(),
constituents_CaloNeutralPrs.end(), prefix + "_nConstituents", 100 );
}
return result;
}
......@@ -229,26 +434,29 @@ StatusCode TupleToolJetsBase::initialize() {
const StatusCode sc = TupleToolBase::initialize();
if ( sc.isFailure() ) return sc;
// get LoKi objects
// charge = LoKi::Cuts::SUMTREE (LoKi::Cuts::Q, LoKi::Cuts::ALL, 0.0 ) ;
charge = LoKi::Cuts::SUMTREE( LoKi::Cuts::Q, LoKi::Cuts::ALL, 0.0 );
positiveParticles = LoKi::Cuts::NINTREE( LoKi::Cuts::Q > 0 );
negativeParticles = LoKi::Cuts::NINTREE( LoKi::Cuts::Q < 0 );
neutralParticles = LoKi::Cuts::NINTREE( LoKi::Cuts::Q == 0 );
maxPT = LoKi::Cuts::MAXTREE( LoKi::Cuts::PT, LoKi::Cuts::BASIC, -1 );
m_M = LoKi::Particles::Mass();
m_MM = LoKi::Particles::MeasuredMass();
m_dva = Gaudi::Utils::getIDVAlgorithm( contextSvc(), this );
if ( 0 == m_dva ) return Error( "Couldn't get parent DVAlgorithm", StatusCode::FAILURE );
// Get distance calculator
m_dist = m_dva->distanceCalculator();
if ( !m_dist ) { return Error( "Unable to retrieve the IDistanceCalculator tool", StatusCode::FAILURE ); }
return sc;
}
double TupleToolJetsBase::MaxSumNPart( const LHCb::Particle* jet, unsigned int n, const LoKi::Types::Fun& fun,
SmartRefVector<LHCb::Particle>* SortedDaughters ) {
// fun == LoKi::Cuts::ONE;
// if (jet->daughters().size() <= n)
// return LoKi::Cuts::SUMTREE ( LoKi::Cuts::BASIC, fun)(jet);
if ( SortedDaughters && SortedDaughters->size() ) // Vector given and filled
if ( SortedDaughters && SortedDaughters->size() )
return LoKi::Cuts::ASUM( fun )(
LHCb::Particle::ConstVector( SortedDaughters->begin(), SortedDaughters->begin() + n ) );
SmartRefVector<LHCb::Particle> daughters; // empty vector to take address of in case SortedDaughters == NULL
if ( !SortedDaughters ) SortedDaughters = &daughters;
SmartRefVector<LHCb::Particle> constituents;
if ( !SortedDaughters ) SortedDaughters = &constituents;
SortedDaughters->assign( jet->daughters().begin(), jet->daughters().end() );
sort( SortedDaughters->begin(), SortedDaughters->end(),
Comperator<const LHCb::Particle*, const LHCb::Particle*>( fun ) );
......
/*****************************************************************************\
* (c) Copyright 2000-2019 CERN for the benefit of the LHCb Collaboration *
* (c) Copyright 2000-2022 CERN for the benefit of the LHCb Collaboration *
* *
* This software is distributed under the terms of the GNU General Public *
* Licence version 3 (GPL Version 3), copied verbatim in the file "COPYING". *
......@@ -7,15 +7,30 @@
* In applying this licence, CERN does not waive the privileges and immunities *
* granted to it by virtue of its status as an Intergovernmental Organization *
* or submit itself to any jurisdiction. *
* *
* Tuple tool to add jet's kinematics and jet's constituents' variables. *
* To add constituent's variables choose withJetConstituents = True *
* in your DaVinci script. *
\*****************************************************************************/
#ifndef TUPLETOOLJETSBASE_H
#define TUPLETOOLJETSBASE_H 1
#include "DecayTreeTupleBase/TupleToolBase.h"
#include "Event/CaloCluster.h"
#include "Event/HltObjectSummary.h"
#include "Event/Particle.h"
#include "Event/ProtoParticle.h"
#include "Event/RecVertex.h"
#include "Event/Track.h"
#include "Kernel/IDistanceCalculator.h"
#include "Kernel/IParticleTupleTool.h"
#include "Kernel/JetEnums.h"
#include "Kernel/ParticleProperty.h"
#include "LoKi/Algo.h"
#include "LoKi/LoKi.h"
#include "Kernel/ITriggerTisTos.h"
// autor: Albert Bursche
class TupleToolJetsBase : public TupleToolBase, virtual public IParticleTupleTool {
public:
......@@ -24,14 +39,17 @@ public:
StatusCode initialize() override;
protected:
Tuples::Tuple* m_tuple;
LoKi::Types::Fun charge;
LoKi::Types::Fun positiveParticles;
LoKi::Types::Fun negativeParticles;
LoKi::Types::Fun neutralParticles;
LoKi::Types::Fun maxPT;
LoKi::Types::Fun m_M;
LoKi::Types::Fun m_MM;
Tuples::Tuple* m_tuple;
LoKi::Types::Fun charge;
LoKi::Types::Fun positiveParticles;
LoKi::Types::Fun negativeParticles;
LoKi::Types::Fun neutralParticles;
LoKi::Types::Fun maxPT;
LoKi::Types::Fun m_M;
LoKi::Types::Fun m_MM;
bool m_withJetConstituents;
IDVAlgorithm* m_dva; // parent DVA Algorithm
const IDistanceCalculator* m_dist; // for obtaining the best PV
bool WriteJetToTuple( const LHCb::Particle* jet, std::string prefix );
double MaxSumNPart( const LHCb::Particle* jet, unsigned int n, const LoKi::Types::Fun& fun,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment