diff --git a/Pr/PrFitParams/src/PrFitMatchParams.cpp b/Pr/PrFitParams/src/PrFitMatchParams.cpp
index a8c81ed940acf07ad7a8067df9020712d36a963f..79dd7c310d6a18f6ade649211d5be08db1281f56 100644
--- a/Pr/PrFitParams/src/PrFitMatchParams.cpp
+++ b/Pr/PrFitParams/src/PrFitMatchParams.cpp
@@ -1,5 +1,4 @@
-// $Id: KsFitParams.cpp,v 1.5 2009/08/19 14:24:18 ocallot Exp $
-// Include files
+// Include files
// from Gaudi
#include "GaudiKernel/SystemOfUnits.h"
@@ -44,15 +43,15 @@ StatusCode PrFitMatchParams::initialize() {
if( m_zMagParams.empty() ) info() << "no starting values for magnet parameters provided. Will not calculate anything" << endmsg;
if( m_momParams.empty() ) info() << "no starting values for momentum parameters provided. Will not calculate anything" << endmsg;
if( m_bendYParams.empty() ) info() << "no starting values for y bending parameters provided. Will not calculate anything" << endmsg;
-
+
m_fitTool = tool<IPrFitTool>( "PrFitTool" );
-
+
m_zMagPar.init( "zMagnet" , m_zMagParams );
m_momPar.init ( "momentum" , m_momParams );
m_bendParamY.init( "bendParamY", m_bendYParams);
-
-
+
+
return StatusCode::SUCCESS;
}
//=============================================================================
@@ -67,7 +66,7 @@ StatusCode PrFitMatchParams::execute() {
// -- As this algorithm is the only one running, even the usual counter is not included.
if( m_nEvent % 100 == 0 ) info() << "Event " << m_nEvent << endmsg;
-
+
debug() << "Processing event: " << m_nEvent << endmsg;
@@ -77,7 +76,7 @@ StatusCode PrFitMatchParams::execute() {
return StatusCode::SUCCESS;
}
-
+
LHCb::MCParticles* partCtnr = get<LHCb::MCParticles>( LHCb::MCParticleLocation::Default );
// Get the FT and UT hits
@@ -90,14 +89,14 @@ StatusCode PrFitMatchParams::execute() {
LHCb::MCParticles::const_iterator pItr;
const LHCb::MCParticle* part;
SmartRefVector<LHCb::MCVertex> vDecay;
-
+
// A container for used hits
std::vector<Gaudi::XYZPoint> trHits;
std::vector<Gaudi::XYZPoint> UTHits;
std::vector<Gaudi::XYZPoint> vHits;
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
-
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
+
for ( pItr = partCtnr->begin(); partCtnr->end() != pItr; pItr++ ) {
part = *pItr;
@@ -130,57 +129,57 @@ StatusCode PrFitMatchParams::execute() {
debug() << "--- Found particle key " << part->key() << endmsg;
-
+
UTHits.clear();
trHits.clear();
vHits.clear();
-
+
// Get the Velo hits
for ( auto vHitIt = veloHits->begin() ; veloHits->end() != vHitIt ; vHitIt++ ) {
if ( (*vHitIt)->mcParticle() == part ) {
vHits.push_back( (*vHitIt)->midPoint() );
}
}
-
+
// Get the FT hits
for ( auto iHitFt = ftHits->begin() ; ftHits->end() != iHitFt ; iHitFt++ ) {
if ( (*iHitFt)->mcParticle() == part ) {
trHits.push_back( (*iHitFt)->midPoint() );
}
}
-
+
// Get the UT hits
for ( auto iHitut = utHits->begin(); utHits->end() != iHitut ; iHitut++ ) {
if ( (*iHitut)->mcParticle() == part ) {
UTHits.push_back( (*iHitut)->midPoint() );
}
}
-
+
if ( 12 > trHits.size() || 6 > vHits.size() ) continue;
if( m_requireUTHits){
if( 3 > UTHits.size() ) continue;
}
-
+
debug() << " particle momentum = " << part->momentum().R() / Gaudi::Units::GeV << " GeV"
<< endmsg;
-
+
//== Fill ntuple
double pz = part->momentum().Z();
double plXT = part->momentum().X() / pz;
double plYT = part->momentum().Y() / pz;
-
+
Tuple tTrack = nTuple( m_tupleName );
-
+
tTrack->column("pz",pz);
tTrack->column("plXT", plXT);
tTrack->column("plYT", plYT);
-
-
+
+
m_nTrack++;
-
+
//== Fit the UT area
- // -- This is not needed at the moment for the matching, but it was kept it
+ // -- This is not needed at the moment for the matching, but it was kept it
// -- as it might be used at some point in the future
debug() << " UT: ";
const auto utFitXLine = m_fitTool->fitLine(UTHits, IPrFitTool::XY::X, m_zUT1);
@@ -202,7 +201,7 @@ StatusCode PrFitMatchParams::execute() {
cxt2 = std::get<2>(*utFitXParab);
debug() << format( " x %7.1f tx %7.4f y %7.1f ty %7.4f ",
- axt, bxt, ayt, byt )
+ axt, bxt, ayt, byt )
<< endmsg;;
tTrack->column( "axt" , axt );
tTrack->column( "bxt", bxt );
@@ -211,8 +210,8 @@ StatusCode PrFitMatchParams::execute() {
tTrack->column( "axt2", axt2 );
tTrack->column( "bxt2", bxt2 );
tTrack->column( "cxt2", cxt2 );
-
-
+
+
std::vector<Gaudi::XYZPoint>::const_iterator itP;
if ( msgLevel( MSG::DEBUG ) ) {
for ( itP = UTHits.begin(); UTHits.end() > itP; itP++ ) {
@@ -225,7 +224,7 @@ StatusCode PrFitMatchParams::execute() {
) << endmsg;
}
}
-
+
// -- Fit the T-stations
// -- x(z) = ax + bx*z + cx*z*z + dx*z*z*z
// -- y(z) = ay + by*z
@@ -250,11 +249,11 @@ StatusCode PrFitMatchParams::execute() {
tTrack->column( "dx", dx );
tTrack->column( "ay", ay );
tTrack->column( "by", by );
-
+
if ( msgLevel( MSG::DEBUG ) ) {
- debug() << format( " tr: x%7.1f %7.4f %7.3f %7.3f y%7.1f %7.4f",
+ debug() << format( " tr: x%7.1f %7.4f %7.3f %7.3f y%7.1f %7.4f",
ax, bx, 1.e6*cx, 1.e9*dx, ay, by ) << endmsg;
-
+
for ( itP = trHits.begin(); trHits.end() > itP; itP++ ) {
const double dz = (*itP).z()-m_zRef;
debug() << format( " : %7.1f %7.1f %7.1f dx %7.3f dy %7.3f",
@@ -264,7 +263,7 @@ StatusCode PrFitMatchParams::execute() {
) << endmsg;
}
}
-
+
// -- Fit the velo area
// -- x(z) = axv + bxv*z
// -- y(z) = ayv + byv*z
@@ -285,11 +284,11 @@ StatusCode PrFitMatchParams::execute() {
tTrack->column( "bxv", bxv );
tTrack->column( "ayv", ayv );
tTrack->column( "byv", byv );
-
+
if ( msgLevel( MSG::DEBUG ) ) {
- debug() << format( " velo: x%7.1f %7.4f y%7.1f %7.4f",
+ debug() << format( " velo: x%7.1f %7.4f y%7.1f %7.4f",
axv, bxv, ayv, byv ) << endmsg;
-
+
for ( itP = vHits.begin(); vHits.end() > itP; itP++ ) {
const double dz = (*itP).z()-m_zRef;
debug() << format( " : %7.1f %7.1f %7.1f dx %7.3f dy %7.3f",
@@ -299,14 +298,14 @@ StatusCode PrFitMatchParams::execute() {
) << endmsg;
}
}
-
+
// -- This is for finding the zMagnet position when using straight lines from the Velo and the T-stations
// -- Only really makes sense in x, as for y the different y resolutions of Velo and T-stations matter a lot
double zMagnet = (axv-bxv*m_zVelo - (ax-bx*m_zRef) ) / (bx-bxv);
double zMagnety = (ayt-byt*m_zUT1 - (ay-by*m_zRef) ) / (by-byt);
double dSlope = fabs(bx - bxv);
double dSlopeY = fabs(by - byv);
-
+
m_zMagPar.setFun( 0, 1. );
m_zMagPar.setFun( 1, fabs(dSlope) );
m_zMagPar.setFun( 2, dSlope*dSlope );
@@ -314,7 +313,7 @@ StatusCode PrFitMatchParams::execute() {
m_zMagPar.setFun( 4, bxv*bxv );
double zEst = m_zMagPar.sum();
-
+
m_zMagPar.addEvent( zMagnet-zEst );
tTrack->column( "zEst", zEst );
@@ -323,40 +322,40 @@ StatusCode PrFitMatchParams::execute() {
// -- This is the parameter that defines the bending in y
double bendParamY = (ay - ayv) + ( m_zMatchY - m_zRef )*by - (m_zMatchY - m_zVelo)*byv;
-
+
m_bendParamY.setFun( 0, dSlope*dSlope*byv );
m_bendParamY.setFun( 1, dSlopeY*dSlopeY*byv );
-
+
double bendParamYEst = m_bendParamY.sum();
m_bendParamY.addEvent( bendParamY - bendParamYEst );
// ------------------------------------------------------
-
+
tTrack->column("dSlope",dSlope);
tTrack->column("dSlope2",dSlope*dSlope);
-
+
// -- Need to write something to calculate the momentum Params
const double charge = part->particleID().threeCharge()/3;
const double qOverP = charge/momentum;
-
- // -- The magnet scale factor is hardcoded, as then one does not need to run the
- // -- field service
+
+ // -- The magnet scale factor is hardcoded, as then one does not need to run the
+ // -- field service
const double proj = sqrt( ( 1. + bxv*bxv + byv*byv ) / ( 1. + bxv*bxv ) );
const double coef = (bxv - bx)/(proj*m_magnetScaleFactor*Gaudi::Units::GeV*qOverP);
-
+
m_momPar.setFun(0, 1.);
m_momPar.setFun(1, bx*bx );
m_momPar.setFun(2, bx*bx*bx*bx );
m_momPar.setFun(3, bx*bxv );
m_momPar.setFun(4, byv*byv );
m_momPar.setFun(5, byv*byv*byv*byv );
-
+
double coefEst = m_momPar.sum();
m_momPar.addEvent( coef - coefEst );
tTrack->write();
}
-
+
return StatusCode::SUCCESS;
}
@@ -365,7 +364,7 @@ StatusCode PrFitMatchParams::execute() {
// Determine resolution
//=============================================================================
void PrFitMatchParams::resolution(){
-
+
LinkedFrom<LHCb::Track, LHCb::MCParticle> myVeloLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Velo);
LinkedFrom<LHCb::Track, LHCb::MCParticle> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
@@ -380,20 +379,20 @@ void PrFitMatchParams::resolution(){
std::vector<Gaudi::XYZPoint> VPHits;
LHCb::MCParticles* mcParts = getIfExists<LHCb::MCParticles>( LHCb::MCParticleLocation::Default );
- if ( msgLevel( MSG::ERROR ) && !mcParts ) error() << "Could not find MCParticles at: "
- << LHCb::MCParticleLocation::Default
+ if ( msgLevel( MSG::ERROR ) && !mcParts ) error() << "Could not find MCParticles at: "
+ << LHCb::MCParticleLocation::Default
<< endmsg;
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
-
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
+
for(const LHCb::MCParticle* mcPart: *mcParts){
-
-
+
+
if( mcPart == nullptr ) continue;
if( std::abs(mcPart->particleID().pid()) == 11 ) continue;
-
+
if( !trackInfo.hasVeloAndT( mcPart) ) continue;
-
+
LHCb::Track* vtrack = myVeloLink.first( mcPart );
LHCb::Track* ttrack = mySeedLink.first( mcPart );
@@ -407,7 +406,7 @@ void PrFitMatchParams::resolution(){
debug()<< "Too far from beampipe" << endmsg;
continue; // particles from close the beam line
}
-
+
bool hasInteractionVertex = false;
SmartRefVector<LHCb::MCVertex> endV = mcPart->endVertices();
for ( SmartRefVector<LHCb::MCVertex>::const_iterator itV = endV.begin() ;
@@ -418,25 +417,25 @@ void PrFitMatchParams::resolution(){
debug()<< "Interaction vertex found. skipping" << endmsg;
continue;
}
-
-
+
+
FTHits.clear();
VPHits.clear();
-
+
// Get the Velo hits
for ( auto vHitIt = vpHits->begin() ; vpHits->end() != vHitIt ; vHitIt++ ) {
if ( (*vHitIt)->mcParticle() == mcPart ) {
VPHits.push_back( (*vHitIt)->midPoint() );
}
}
-
+
// Get the FT hits
for ( auto iHitFt = ftHits->begin() ; ftHits->end() != iHitFt ; iHitFt++ ) {
if ( (*iHitFt)->mcParticle() == mcPart ) {
FTHits.push_back( (*iHitFt)->midPoint() );
}
}
-
+
if ( 12 > FTHits.size() || 6 > VPHits.size() ) continue;
debug()<<"got my hits " << FTHits.size()+VPHits.size() << endmsg;
@@ -472,42 +471,42 @@ void PrFitMatchParams::resolution(){
const auto zMagnet = (axv-bxv*m_zVelo - (ax-bx*m_zRef) ) / (bx-bxv);
const auto xMagnet = axv + ( zMagnet - m_zVelo ) * bxv;
const auto yMagnet = ayv + ( zMagnet - m_zVelo ) * byv;
-
+
LHCb::State* tstate = &ttrack->closestState( 10000. );
LHCb::State* vstate = &vtrack->closestState( m_zVelo );
//now other things
const auto dSlopeExp = std::abs( vstate->tx() - tstate->tx() );
- const auto dSlope2Exp = dSlopeExp*dSlopeExp;
+ const auto dSlope2Exp = dSlopeExp*dSlopeExp;
const auto dSlopeYExp = std::abs( vstate->ty() - tstate->ty() );
- const auto dSlopeY2Exp = dSlopeYExp*dSlopeYExp;
+ const auto dSlopeY2Exp = dSlopeYExp*dSlopeYExp;
- const auto zMagnetExp =m_zMagParams [0] +
- m_zMagParams[1] * dSlopeExp +
+ const auto zMagnetExp =m_zMagParams [0] +
+ m_zMagParams[1] * dSlopeExp +
m_zMagParams[2] * dSlope2Exp +
m_zMagParams[3] * std::abs(tstate->x()) +
m_zMagParams[4] * vstate->tx()*vstate->tx() ;
-
+
debug()<< m_zMagParams [0] <<" "<<
- m_zMagParams[1] <<" " <<
+ m_zMagParams[1] <<" " <<
m_zMagParams[2] <<" " <<
m_zMagParams[3] <<" " <<
m_zMagParams[4] << endmsg;
-
-
+
+
const float xVExp = vstate->x() + (zMagnetExp - vstate->z()) * vstate->tx();
// -- This is the function that calculates the 'bending' in y-direction
// -- The parametrisation can be derived with the MatchFitParams package
const float yVExp = vstate->y() + (m_zMatchY - vstate->z()) * vstate->ty() +
vstate->ty() * ( m_bendYParams[0]*dSlope2Exp + m_bendYParams[1]*dSlopeY2Exp );
-
+
const float xSExp = tstate->x() + (zMagnetExp - tstate->z()) * tstate->tx();
const float ySExp = tstate->y() + (m_zMatchY - tstate->z()) * tstate->ty();
-
+
const float distXExp = xSExp - xVExp;
const float distYExp = ySExp - yVExp;
-
+
const float dslXExp = vstate->tx() - tstate->tx();
const float dslYExp = vstate->ty() - tstate->ty();
const float teta2 = vstate->ty()*vstate->ty() + tstate->tx()*tstate->tx();
@@ -527,7 +526,7 @@ void PrFitMatchParams::resolution(){
resoTuple->column( "distXExp", distXExp );
resoTuple->column( "distYExp", distYExp );
-
+
resoTuple->column( "xVExp", xVExp );
resoTuple->column( "xSExp", xSExp );
@@ -536,17 +535,17 @@ void PrFitMatchParams::resolution(){
resoTuple->column( "pExp", pExp );
resoTuple->column( "pTrue", p );
-
+
resoTuple->column( "byvExp", vstate->ty() );
resoTuple->column( "byv", byv );
resoTuple->column( "bxvExp", vstate->tx() );
resoTuple->column( "bxv", bxv );
-
-
+
+
resoTuple->column( "byExp", tstate->ty() );
resoTuple->column( "by", by );
-
+
resoTuple->column( "bxExp", tstate->tx());
resoTuple->column( "bx", bx );
@@ -555,7 +554,7 @@ void PrFitMatchParams::resolution(){
resoTuple->column( "teta2", teta2 );
resoTuple->write();
-
+
}
}
@@ -608,13 +607,13 @@ float PrFitMatchParams::momentum(const LHCb::State* vState, const LHCb::State* t
+m_momParams[4]* vars[4];
debug()<< m_momParams [0] <<" "<<
- m_momParams[1] <<" " <<
+ m_momParams[1] <<" " <<
m_momParams[2] <<" " <<
m_momParams[3] <<" " <<
m_momParams[4] << endmsg;
debug()<< vars[0] <<" "<<
- vars[1] <<" " <<
+ vars[1] <<" " <<
vars[2] <<" " <<
vars[3] <<" " <<
vars[4] << endmsg;
@@ -623,7 +622,7 @@ float PrFitMatchParams::momentum(const LHCb::State* vState, const LHCb::State* t
float scaleFactor = m_magFieldSvc->signedRelativeCurrent();
float P = ( coef * Gaudi::Units::GeV * proj *scaleFactor )/(txV-txT);
-
+
return P;
}
diff --git a/Pr/PrFitParams/src/PrFitSeedParams.cpp b/Pr/PrFitParams/src/PrFitSeedParams.cpp
index 9aa2bef6a59be0a01901435ae8aa1f514bcac6f8..0d8228f0d29c7c81514d1ff8fd282f19100a34de 100755
--- a/Pr/PrFitParams/src/PrFitSeedParams.cpp
+++ b/Pr/PrFitParams/src/PrFitSeedParams.cpp
@@ -1,4 +1,4 @@
-// Include files
+// Include files
// from Gaudi
#include "Event/MCTrackInfo.h"
@@ -18,7 +18,7 @@
//-----------------------------------------------------------------------------
// Implementation file for class : PrFitSeedParams
-//
+//
// 2006-12-08 : Olivier Callot
// 2013-01-23 : Yasmine Amhis
// Adapt to work with Fiber Tracker and UT
@@ -39,7 +39,7 @@ PrFitSeedParams::PrFitSeedParams( const std::string& name,
declareProperty( "TupleName", m_tupleName = "Track" );
declareProperty( "ZRef", m_zRef = StateParameters::ZMidT );
declareProperty( "ZSeed", m_zSeed = StateParameters::ZEndT );
- //THIS NEEDS TO BE UPDATED FOR THE UT
+ //THIS NEEDS TO BE UPDATED FOR THE UT
declareProperty( "ZTT", m_zTT = 2485.0 * Gaudi::Units::mm );
// these are the "tunables" for PatSeeding
@@ -57,14 +57,14 @@ PrFitSeedParams::PrFitSeedParams( const std::string& name,
declareProperty( "dRatioParams", m_dRatio = tmp4 );
std::vector<double> tmp5 = boost::assign::list_of(1.25e-14);
declareProperty( "yCorrectionParams", m_yCorrection = tmp5 );
-
+
m_nEvent = 0;
m_nTrack = 0;
}
//=============================================================================
// Destructor
//=============================================================================
-PrFitSeedParams::~PrFitSeedParams() {}
+PrFitSeedParams::~PrFitSeedParams() {}
//=============================================================================
// Initialization
@@ -94,34 +94,34 @@ StatusCode PrFitSeedParams::execute() {
debug() << "==> Execute" << endmsg;
-
+
m_nEvent += 1;
LHCb::MCParticles* partCtnr = get<LHCb::MCParticles>( LHCb::MCParticleLocation::Default );
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
-
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
+
// Get the FT hits
LHCb::MCHits* ftHits = get<LHCb::MCHits>( LHCb::MCHitLocation::FT);
-
-
+
+
// Get the UT hits
LHCb::MCHits* UTHits = get<LHCb::MCHits>( LHCb::MCHitLocation::UT );
-
+
LHCb::MCParticles::const_iterator pItr;
const LHCb::MCParticle* part;
const LHCb::MCParticle* mother;
const LHCb::MCVertex* vOrigin;
SmartRefVector<LHCb::MCVertex> vDecay;
-
- for ( pItr = partCtnr->begin(); partCtnr->end() != pItr; pItr++ )
- {//here loop over the MC particle container
+
+ for ( pItr = partCtnr->begin(); partCtnr->end() != pItr; pItr++ )
+ {//here loop over the MC particle container
part = *pItr;
if ( 0 == trackInfo.fullInfo( part ) ) continue;
if ( !trackInfo.hasT( part ) ) continue;
@@ -153,22 +153,22 @@ StatusCode PrFitSeedParams::execute() {
// A container for used hits
std::vector<Gaudi::XYZPoint> trHits;
std::vector<Gaudi::XYZPoint> utHits;
-
-
+
+
// Get the FT hits
-
- for ( LHCb::MCHits::const_iterator fiberHit = ftHits->begin() ;
+
+ for ( LHCb::MCHits::const_iterator fiberHit = ftHits->begin() ;
ftHits->end() != fiberHit ; fiberHit++ ) {
if ( (*fiberHit)->mcParticle() == part ) {
trHits.push_back( (*fiberHit)->midPoint() );
}
}
-
+
// Get the UT hits
- for ( LHCb::MCHits::const_iterator iHitut = UTHits->begin() ;
+ for ( LHCb::MCHits::const_iterator iHitut = UTHits->begin() ;
UTHits->end() != iHitut ; iHitut++ ) {
if ( (*iHitut)->mcParticle() == part ) {
utHits.push_back( (*iHitut)->midPoint() );
@@ -219,7 +219,7 @@ StatusCode PrFitSeedParams::execute() {
debug() << format( "p %7.3f, N%4d, ax%8.2f bx%8.2f cx%8.2f dp%10.4f",
momentum/1000., trHits.size(),
ax, 1.e3*bx, 1.e6*cx, dp ) << endmsg;
- }
+ }
// work out intercept point of line joining T1 and T3 with z = 0
const double dzT1 = StateParameters::ZBegT - m_zRef;
@@ -240,7 +240,7 @@ StatusCode PrFitSeedParams::execute() {
m_initialArrowPar.addEvent( dSag );
}
-
+
double axt, bxt, cxt, ayt, byt, zMagnet, zEst;
if ( 2 < utHits.size() ) {
//m_fitTool->fitLine( ttHits, 0, m_zTT, axt, bxt );
@@ -253,13 +253,13 @@ StatusCode PrFitSeedParams::execute() {
}
std::tie(axt, bxt, cxt) = *xResult;
std::tie(ayt, byt) = *yResult;
-
+
const double dz = m_zSeed - m_zRef;
const double x0 = ax + dz * ( bx + dz * ( cx + dz * dx ) );
const double tx = bx + dz * (2*cx + dz * 3*dx );
zMagnet = (axt-bxt*m_zTT - (x0-tx*m_zSeed) ) / (tx-bxt);
-
+
m_zMagPar.setFun( 0, 1. );
m_zMagPar.setFun( 1, by*by );
m_zMagPar.setFun( 2, bx*bx );
@@ -276,7 +276,7 @@ StatusCode PrFitSeedParams::execute() {
byt = 0.;
}
-
+
m_dRatioPar.setFun(0, 1.0);
m_dRatioPar.setFun(1, fabs(cx));
@@ -293,7 +293,7 @@ StatusCode PrFitSeedParams::execute() {
m_yCorrectionPar.setFun(0, 1.0);
if (400. < fabs(y0))
m_yCorrectionPar.addEvent(fabs((by * by * cx * cx) / y0) - m_yCorrectionPar.sum());
-
+
tTrack->column( "ax", ax );
tTrack->column( "bx", bx );
tTrack->column( "cx", cx );
@@ -316,13 +316,13 @@ StatusCode PrFitSeedParams::execute() {
tTrack->column( "zEst", zEst );
tTrack->write();
-
-
- }// end of the big loop
-
-
+ }// end of the big loop
+
+
+
+
return StatusCode::SUCCESS;
}
@@ -338,7 +338,7 @@ StatusCode PrFitSeedParams::finalize() {
msg << "============================================" << endmsg;
-
+
m_initialArrowPar.updateParameters( msg );
m_momentumScalePar.updateParameters( msg );
m_zMagPar.updateParameters( msg );
@@ -352,16 +352,16 @@ StatusCode PrFitSeedParams::finalize() {
m_dRatioPar.printPythonParams( name() );
m_yCorrectionPar.printPythonParams( name() );
std::cout << std::endl;
-
+
std::string toolName = "ToolSvc.PatSeedingTool";
-
+
m_initialArrowPar.printPythonParams( toolName );
m_momentumScalePar.printPythonParams( toolName );
m_zMagPar.printPythonParams( toolName );
m_dRatioPar.printPythonParams( toolName );
m_yCorrectionPar.printPythonParams( toolName );
std::cout << std::endl;
-
+
return GaudiTupleAlg::finalize(); // must be called after all other actions
}
diff --git a/Pr/PrMCTools/src/PrCheatedSciFiTracking.cpp b/Pr/PrMCTools/src/PrCheatedSciFiTracking.cpp
index 52d339a487f4a8a78e5f91dcf1d2d5c48039f075..ab26eb6146c33e73224f552fe702874e4e0b4217 100644
--- a/Pr/PrMCTools/src/PrCheatedSciFiTracking.cpp
+++ b/Pr/PrMCTools/src/PrCheatedSciFiTracking.cpp
@@ -1,18 +1,61 @@
-// Include files
+// Include files
-// local
+// from Gaudi
+#include "GaudiAlg/Transformer.h"
+
+#include "PrKernel/PrFTHitHandler.h"
#include "Event/Track.h"
#include "Event/StateParameters.h"
-#include "PrCheatedSciFiTracking.h"
#include "Event/FTCluster.h"
-#include "Linker/LinkedTo.h"
#include "Event/MCTrackInfo.h"
+#include "Linker/LinkedTo.h"
//-----------------------------------------------------------------------------
// Implementation file for class : PrCheatedSciFiTracking
//
// 2015-03-23 : Michel De Cian
//-----------------------------------------------------------------------------
+
+/** @class PrCheatedSciFiTracking PrCheatedSciFiTracking.h
+ *
+ * Cheated track reconstruction in the SciFi.
+ * It creates tracks by getting all Clusters associated to a reconstructible MCParticle.
+ * Cuts can be set on the minimum number of total hits, hits in x layers and hits in stereo layers.
+ * Top and bottom modules are not mixed in the x layers.
+ *
+ * - HitManagerName: Name of the hit manager
+ * - DecodeData: Decode the data
+ * - OutputName: Name of the output location
+ * - NumZones: Number of zones (normally 2 x number of layers if no y-segmentation)
+ * - MinXHits: Minimum number of required hits in x layers to make a track.
+ * - MinStereoHits: Minimum number of required hits in stereo layers to make a track.
+ * - MinTotHits: Minimum number of total hits to make a track.
+ *
+ *
+ * @author Michel De Cian
+ * @date 2015-03-23
+ */
+class PrCheatedSciFiTracking : public Gaudi::Functional::Transformer<LHCb::Tracks(const PrFTHitHandler<PrHit>&,
+ const LHCb::MCParticles&,
+ const LHCb::MCProperty&)> {
+public:
+ /// Standard constructor
+ PrCheatedSciFiTracking( const std::string& name, ISvcLocator* pSvcLocator );
+
+ /// make cheated tracks by getting the clusters matched to an MCParticle
+ LHCb::Tracks operator()(const PrFTHitHandler<PrHit>&,
+ const LHCb::MCParticles&,
+ const LHCb::MCProperty&) const override;
+
+private:
+
+ Gaudi::Property<int> m_numZones = { this, "NumZones", 24 };
+ Gaudi::Property<int> m_minXHits = { this, "MinXHits", 5 };
+ Gaudi::Property<int> m_minStereoHits = { this, "MinStereoHits", 5 };
+ Gaudi::Property<int> m_minTotHits = { this, "MinTotHits", 10 };
+
+};
+
// Declaration of the Algorithm Factory
DECLARE_COMPONENT( PrCheatedSciFiTracking )
@@ -22,119 +65,86 @@ DECLARE_COMPONENT( PrCheatedSciFiTracking )
PrCheatedSciFiTracking::PrCheatedSciFiTracking( const std::string& name,
ISvcLocator* pSvcLocator)
: Transformer(name, pSvcLocator,
- KeyValue{"FTHitsLocation", PrFTInfo::FTHitsLocation},
+ { KeyValue{"FTHitsLocation", PrFTInfo::FTHitsLocation},
+ KeyValue{"MCParticleLocation", LHCb::MCParticleLocation::Default},
+ KeyValue{"MCPropertyLocation", LHCb::MCPropertyLocation::TrackInfo} },
KeyValue{"OutputName", LHCb::TrackLocation::Seed})
{
- declareProperty( "NumZones", m_numZones = 24 );
- declareProperty( "MinXHits", m_minXHits = 5 );
- declareProperty( "MinStereoHits", m_minStereoHits = 5 );
- declareProperty( "MinTotHits", m_minTotHits = 10 );
-}
-
-//=============================================================================
-// Initialization
-//=============================================================================
-StatusCode PrCheatedSciFiTracking::initialize() {
- StatusCode sc = Transformer::initialize(); // must be executed first
-
- if ( sc.isFailure() ) return sc; // error printed already by GaudiAlgorithm
- if ( msgLevel(MSG::DEBUG) ) debug() << "==> Initialize" << endmsg;
- return StatusCode::SUCCESS;
}
//=============================================================================
-// Main execution
+// make cheated tracks by getting the clusters matched to an MCParticle
//=============================================================================
-LHCb::Tracks PrCheatedSciFiTracking::operator()(const PrFTHitHandler<PrHit>& hitHandler) const {
+LHCb::Tracks PrCheatedSciFiTracking::operator()(const PrFTHitHandler<PrHit>& FTHitHandler,
+ const LHCb::MCParticles& mcParts,
+ const LHCb::MCProperty& mcProps) const
+{
LHCb::Tracks result;
- makeLHCbTracks(result, hitHandler);
- return result;
-}
-//=============================================================================
-// Finalize
-//=============================================================================
-StatusCode PrCheatedSciFiTracking::finalize() {
- if ( msgLevel(MSG::DEBUG) ) debug() << "==> Finalize" << endmsg;
- return Transformer::finalize(); // must be called after all other actions
-}
-//=========================================================================
-// Make the cheated tracks
-//=========================================================================
-void PrCheatedSciFiTracking::makeLHCbTracks(LHCb::Tracks& result,
- const PrFTHitHandler<PrHit>& FTHitHandler) const {
-
+ // TODO: avoid direct use of evtSvc(), once there is a way to do so for LinkedTo & MCTrackInfo...
LinkedTo<LHCb::MCParticle, LHCb::FTCluster> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTClusterLocation::Default );
- LHCb::MCParticles* mcParts = getIfExists<LHCb::MCParticles> ( LHCb::MCParticleLocation::Default );
-
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
-
- for( LHCb::MCParticles::const_iterator iPart = mcParts->begin(); iPart != mcParts->end(); ++iPart){
-
- const LHCb::MCParticle* mcPart = *iPart;
+
+ MCTrackInfo trackInfo( mcProps );
+
+ for( const LHCb::MCParticle* mcPart : mcParts ) {
//const bool isLong = trackInfo.hasVeloAndT( mcPart );
//const bool isDown = trackInfo.hasT( mcPart ) && trackInfo.hasTT( mcPart );
const bool isSeed = trackInfo.hasT( mcPart );
//if( !isLong && !isDown ) continue;
if( !isSeed ) continue;
- LHCb::Track* tmp = new LHCb::Track;
+ auto tmp = std::make_unique< LHCb::Track > ();
tmp->setType( LHCb::Track::Types::Ttrack );
tmp->setHistory( LHCb::Track::History::PrSeeding );
tmp->setPatRecStatus( LHCb::Track::PatRecStatus::PatRecIDs );
- std::vector<int> firedXLayers(m_numZones,0);
- std::vector<int> firedStereoLayers(m_numZones,0);
+ std::vector<int> firedXLayers(m_numZones.value(),0);
+ std::vector<int> firedStereoLayers(m_numZones.value(),0);
int totHits = 0;
-
+
// -- loop over all zones
- for(int iZone = 0; iZone < m_numZones; ++iZone){
+ for(int iZone = 0; iZone < m_numZones.value(); ++iZone){
// -- loop over all hits in a zone
for(auto iHit = FTHitHandler.getIterator_Begin(iZone); FTHitHandler.getIterator_End(iZone) != iHit;++iHit){
auto hit = &*iHit;
- LHCb::MCParticle* mcPart1 = myClusterLink.first( hit->id().ftID() );
bool found = false;
- while( mcPart1 != nullptr){
- if( mcPart1 == mcPart){
- found = true;
- }
- mcPart1 = myClusterLink.next();
+ for(const LHCb::MCParticle* mcPart1 = myClusterLink.first( hit->id().ftID() ) ;
+ mcPart1 && !found;
+ mcPart1 = myClusterLink.next() ){
+ found = ( mcPart1 == mcPart);
}
if( hit->isX() && found ){
- if( firedXLayers[iZone] == 0){
- firedXLayers[iZone]++;
- }
+ if( firedXLayers[iZone] == 0) firedXLayers[iZone]++;
}
if( !(hit->isX()) && found ){
- if( firedStereoLayers[iZone] == 0){
- firedStereoLayers[iZone]++;
- }
+ if( firedStereoLayers[iZone] == 0) firedStereoLayers[iZone]++;
}
if( found ){
totHits++;
tmp->addToLhcbIDs( hit->id() );
- }
+ }
}
}
int sumLowerX = 0;
int sumUpperX = 0;
int sumStereo = 0;
- for(int i = 0; i < m_numZones; i = i+2){
+ for(int i = 0; i < m_numZones.value(); i += 2){
sumLowerX += firedXLayers[i];
}
- for(int i = 1; i < m_numZones; i = i+2){
+ for(int i = 1; i < m_numZones.value(); i += 2){
sumUpperX += firedXLayers[i];
}
- for(int i = 0; i < m_numZones; i++){
+ for(int i = 0; i < m_numZones.value(); i++){
sumStereo += firedStereoLayers[i];
}
- debug() << "sumLowerX: " << sumLowerX
- << " sumUpperX " << sumUpperX
- << " sumStereo " << sumStereo
+ debug() << "sumLowerX: " << sumLowerX
+ << " sumUpperX " << sumUpperX
+ << " sumStereo " << sumStereo
<< " totHits " << totHits << endmsg;
-
- if( (sumLowerX < m_minXHits && sumUpperX < m_minXHits) || sumStereo < m_minStereoHits || totHits < m_minTotHits){
- delete tmp;
+
+ if( (sumLowerX < m_minXHits.value() && sumUpperX < m_minXHits.value()) ||
+ sumStereo < m_minStereoHits.value() ||
+ totHits < m_minTotHits.value()) {
continue;
}
// -- these are obviously useless numbers, but it should just make the checkers happy
@@ -145,8 +155,9 @@ void PrCheatedSciFiTracking::makeLHCbTracks(LHCb::Tracks& result,
tState.setState( 100, 50, z, 0.1, 0.1, qOverP );
//tState.setCovariance( m_geoTool->covariance( qOverP ) );
tmp->addToStates( tState );
- result.insert( tmp );
+ result.insert( tmp.release() );
}
+ return result;
}
diff --git a/Pr/PrMCTools/src/PrCheatedSciFiTracking.h b/Pr/PrMCTools/src/PrCheatedSciFiTracking.h
deleted file mode 100644
index 9d440fbf6300056040168e78aa30b72d55c3a33f..0000000000000000000000000000000000000000
--- a/Pr/PrMCTools/src/PrCheatedSciFiTracking.h
+++ /dev/null
@@ -1,51 +0,0 @@
-#ifndef PRCHEATEDSCIFITRACKING_H
-#define PRCHEATEDSCIFITRACKING_H 1
-
-// Include files
-// from Gaudi
-
-#include "GaudiAlg/Transformer.h"
-#include "PrKernel/PrFTHitHandler.h"
-#include "GaudiKernel/AnyDataHandle.h"
-
-/** @class PrCheatedSciFiTracking PrCheatedSciFiTracking.h
- *
- * Cheated track reconstruction in the SciFi.
- * It creates tracks by getting all Clusters associated to a reconstructible MCParticle.
- * Cuts can be set on the minimum number of total hits, hits in x layers and hits in stereo layers.
- * Top and bottom modules are not mixed in the x layers.
- *
- * - HitManagerName: Name of the hit manager
- * - DecodeData: Decode the data
- * - OutputName: Name of the output location
- * - NumZones: Number of zones (normally 2 x number of layers if no y-segmentation)
- * - MinXHits: Minimum number of required hits in x layers to make a track.
- * - MinStereoHits: Minimum number of required hits in stereo layers to make a track.
- * - MinTotHits: Minimum number of total hits to make a track.
- *
- *
- * @author Michel De Cian
- * @date 2015-03-23
- */
-class PrCheatedSciFiTracking : public Gaudi::Functional::Transformer<LHCb::Tracks(const PrFTHitHandler<PrHit>&)> {
-public:
- /// Standard constructor
- PrCheatedSciFiTracking( const std::string& name, ISvcLocator* pSvcLocator );
-
- LHCb::Tracks operator()(const PrFTHitHandler<PrHit>& hitHandler) const override;
- StatusCode initialize() final override; ///< Algorithm initialization
- StatusCode finalize() final override; ///< Algorithm finalization
-
-private:
-
- /// make cheated tracks by getting the clusters matched to an MCParticle
- void makeLHCbTracks(LHCb::Tracks& result, const PrFTHitHandler<PrHit>& FTHitHandler) const;
-
-
- int m_numZones;
- int m_minXHits;
- int m_minStereoHits;
- int m_minTotHits;
-
-};
-#endif // PRCHEATEDSCIFITRACKING_H
diff --git a/Pr/PrMCTools/src/PrCheatedVP.cpp b/Pr/PrMCTools/src/PrCheatedVP.cpp
index a16a344ed8ef25a9d0f50a5e7090ff32eeb8fc72..7aac20131275c8735b8424844115bc83968ffffa 100644
--- a/Pr/PrMCTools/src/PrCheatedVP.cpp
+++ b/Pr/PrMCTools/src/PrCheatedVP.cpp
@@ -22,12 +22,12 @@ namespace
LHCb::Tracks getTracks(const GaudiAlgorithm& alg, const MCTrackInfo& trackInfo, const IPrFitTool& fitTool, const LHCb::MCParticles& particles, const std::function<void(LHCb::Track&, std::vector<Gaudi::XYZPoint>&, const LHCb::MCParticle*, LinkedFrom<LHCb::VPCluster, MCParticle>&)> getPoints)
{
LHCb::Tracks tracks;
-
+
// Get the association table between MC particles and clusters.
LinkedFrom<LHCb::VPCluster, MCParticle> link(alg.evtSvc(), alg.msgSvc(), LHCb::VPClusterLocation::Default);
-
+
constexpr double zVelo = 0.;
-
+
for (const LHCb::MCParticle *const particle : particles) {
// Skip particles without track info.
if (0 == trackInfo.fullInfo(particle)) continue;
@@ -35,12 +35,12 @@ namespace
if (!trackInfo.hasVelo(particle)) continue;
// Skip electrons.
if (abs(particle->particleID().pid()) == 11) continue;
-
+
LHCb::Track *const track = new LHCb::Track;
std::vector<Gaudi::XYZPoint> points;
getPoints(*track, points, particle, link);
-
- // Make a straight-line fit of the track.
+
+ // Make a straight-line fit of the track.
const auto xResult = fitTool.fitLine(points, IPrFitTool::XY::X, zVelo);
const auto yResult = fitTool.fitLine(points, IPrFitTool::XY::Y, zVelo);
if (!xResult || !yResult)
@@ -53,22 +53,22 @@ namespace
x1 = std::get<1>(*xResult),
y0 = std::get<0>(*yResult),
y1 = std::get<1>(*yResult);
-
+
LHCb::State state;
state.setLocation(LHCb::State::ClosestToBeam);
- state.setState(x0, y0, zVelo,
+ state.setState(x0, y0, zVelo,
x1, y1, 0.);
track->addToStates(state);
if (0 > particle->momentum().z()) {
track->setFlag(LHCb::Track::Flags::Backward, true);
// Cut out backwards tracks.
// delete track;
- // continue;
+ // continue;
}
track->setType(LHCb::Track::Types::Velo);
tracks.insert(track);
}
-
+
return tracks;
}
}
@@ -112,7 +112,7 @@ StatusCode PrCheatedVPBase<useMCHits>::initialize()
//=============================================================================
LHCb::Tracks PrCheatedVP::operator() (const LHCb::MCParticles& particles) const
{
- return getTracks(*this, MCTrackInfo(evtSvc(), msgSvc()), *m_ft, particles, [](LHCb::Track& track, std::vector<Gaudi::XYZPoint>& points, const LHCb::MCParticle *const particle, LinkedFrom<LHCb::VPCluster, MCParticle>& link) {
+ return getTracks(*this, make_MCTrackInfo(evtSvc(), msgSvc()), *m_ft, particles, [](LHCb::Track& track, std::vector<Gaudi::XYZPoint>& points, const LHCb::MCParticle *const particle, LinkedFrom<LHCb::VPCluster, MCParticle>& link) {
for (const LHCb::VPCluster* cluster = link.first(particle); cluster != nullptr; cluster = link.next())
{
track.addToLhcbIDs(LHCb::LHCbID(cluster->channelID()));
@@ -123,7 +123,7 @@ LHCb::Tracks PrCheatedVP::operator() (const LHCb::MCParticles& particles) const
LHCb::Tracks PrCheatedVPMCHits::operator() (const LHCb::MCParticles& particles, const LHCb::MCHits& hits) const
{
- return getTracks(*this, MCTrackInfo(evtSvc(), msgSvc()), *m_ft, particles, [&hits](LHCb::Track& track, std::vector<Gaudi::XYZPoint>& points, const LHCb::MCParticle *const particle, LinkedFrom<LHCb::VPCluster, MCParticle>& link) {
+ return getTracks(*this, make_MCTrackInfo(evtSvc(), msgSvc()), *m_ft, particles, [&hits](LHCb::Track& track, std::vector<Gaudi::XYZPoint>& points, const LHCb::MCParticle *const particle, LinkedFrom<LHCb::VPCluster, MCParticle>& link) {
for (const auto id : link.keyRange(particle)) {
track.addToLhcbIDs(LHCb::LHCbID(LHCb::VPChannelID(id)));
}
diff --git a/Pr/PrMCTools/src/PrChecker.cpp b/Pr/PrMCTools/src/PrChecker.cpp
index 91909c072ed1d5b0227305f8a5823decea382a11..b84d22b4ac403d202f6ff4b195c0cdbfff13eefa 100644
--- a/Pr/PrMCTools/src/PrChecker.cpp
+++ b/Pr/PrMCTools/src/PrChecker.cpp
@@ -62,7 +62,7 @@ PrChecker::PrChecker( const std::string& name, ISvcLocator* pSvcLocator)
declareProperty( "DownTracks", m_downTracks = LHCb::TrackLocation::Downstream );
declareProperty( "UpTracks", m_upTracks = LHCb::TrackLocation::VeloTT );
declareProperty( "BestTracks", m_bestTracks = LHCb::TrackLocation::Default );
-
+
declareProperty( "WriteVeloHistos", m_writeVeloHistos = -1 );
declareProperty( "WriteForwardHistos", m_writeForwardHistos = -1 );
declareProperty( "WriteMatchHistos", m_writeMatchHistos = -1 );
@@ -73,7 +73,7 @@ PrChecker::PrChecker( const std::string& name, ISvcLocator* pSvcLocator)
declareProperty( "WriteBestLongHistos",m_writeBestLongHistos = -1 );
declareProperty( "WriteBestDownstreamHistos", m_writeBestDownstreamHistos = -1 );
declareProperty( "WriteUTHistos", m_writeUTHistos = -1 );
-
+
declareProperty( "Eta25Cut", m_eta25cut = false );
declareProperty( "TriggerNumbers", m_triggerNumbers = false );
declareProperty( "UseElectrons", m_useElectrons = false );
@@ -81,10 +81,6 @@ PrChecker::PrChecker( const std::string& name, ISvcLocator* pSvcLocator)
declareProperty( "Upgrade", m_upgrade = true );
}
-//=============================================================================
-// Destructor
-//=============================================================================
-PrChecker::~PrChecker() {}
//=============================================================================
// Initialization
@@ -102,13 +98,13 @@ StatusCode PrChecker::initialize()
IHistoTool* htool = tool<IHistoTool>( "HistoTool","PrCheckerHistos",this ) ;
GaudiHistoTool* ghtool = dynamic_cast<GaudiHistoTool*>(htool) ;
-
+
// -- catch the possible failure of the dynamic cast
if( ghtool == NULL){
error() << "Dynamic cast of Gaudi Histogramming Tool failed!" << endmsg;
return StatusCode::FAILURE;
}
-
+
ghtool->setHistoDir(histoDir+name()) ;
m_histoTool = htool;
@@ -280,7 +276,7 @@ StatusCode PrChecker::initialize()
// -- UT Counters (excluded for the moment)
/*
- m_utForward = tool<IPrCounter>( "PrUTCounter", "UTForward", this );
+ m_utForward = tool<IPrCounter>( "PrUTCounter", "UTForward", this );
m_utForward->setContainer( m_forwardTracks );
//m_utForward->setFolderPlotName( m_forwardTracks );
m_utForward->setWriteHistos(m_writeUTHistos);
@@ -324,7 +320,7 @@ StatusCode PrChecker::initialize()
m_allCounters.push_back( m_utMatch );
m_allCounters.push_back( m_utbestLong);
*/
-
+
for ( std::vector<IPrCounter*>::iterator itC = m_allCounters.begin();
m_allCounters.end() != itC; ++itC ) {
(*itC)->setUseEta25Cut(m_eta25cut);
@@ -333,7 +329,7 @@ StatusCode PrChecker::initialize()
if( !m_upgrade ) m_idLoc = "Pat/LHCbID";
-
+
return StatusCode::SUCCESS;
}
@@ -350,16 +346,16 @@ StatusCode PrChecker::execute() {
error() << "Could not find MCParticles at " << LHCb::MCParticleLocation::Default << endmsg;
return StatusCode::SUCCESS;
}
-
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LHCb::MCVertices* mcVert = getIfExists<LHCb::MCVertices>(LHCb::MCVertexLocation::Default );
if( mcVert == nullptr ){
error() << "Could not find MCVertices at " << LHCb::MCParticleLocation::Default << endmsg;
return StatusCode::SUCCESS;
}
-
+
unsigned int nPrim = 0;
for ( LHCb::MCVertices::iterator itV = mcVert->begin(); mcVert->end() != itV; ++itV ) {
if ( (*itV)->isPrimary() ) {
@@ -370,16 +366,16 @@ StatusCode PrChecker::execute() {
if ( trackInfo.hasVelo( *itP ) ) nbVisible++;
}
}
- if ( nbVisible > 4 )
+ if ( nbVisible > 4 )
++nPrim;
}
}
-
+
for ( std::vector<IPrCounter*>::iterator itC = m_allCounters.begin();
m_allCounters.end() != itC; ++itC ) {
(*itC)->initEvent(m_histoTool,nPrim);
}
-
+
//== Build a table (vector of vectors) of ids per MCParticle, indexed by MCParticle key.
AllLinks<LHCb::MCParticle> allIds( evtSvc(), msgSvc(), m_idLoc );
std::vector< std::vector<int> > linkedIds;
@@ -393,7 +389,7 @@ StatusCode PrChecker::execute() {
linkedIds[part->key()].push_back( allIds.key() );
part = allIds.next();
}
-
+
LHCb::MCParticles::const_iterator itP;
for ( itP = mcParts->begin(); mcParts->end() != itP; ++itP ) {
@@ -402,24 +398,24 @@ StatusCode PrChecker::execute() {
if( msgLevel(MSG::VERBOSE) ) verbose() << "checking MCPart " << part->key() << endmsg;
-
+
//Add Eta Cut!
-
+
// -- Take electrons into account?
bool isElectron = abs( part->particleID().pid() ) == 11;
if(m_useElectrons){
isElectron = false;
}
-
-
+
+
bool isLong = trackInfo.hasVeloAndT( part );
//isLong = isLong && !isElectron; // and not electron
bool isDown = trackInfo.hasT( part ) && trackInfo.hasTT( part );
//isDown = isDown && !isElectron; // and not electron
-
+
bool over5 = 5000. < fabs( part->p() );
bool trigger = 3000. < fabs( part->p() ) && 500. < fabs( part->pt() );
bool isInVelo = trackInfo.hasVelo( part );
@@ -428,7 +424,7 @@ StatusCode PrChecker::execute() {
bool strangeDown = false;
bool fromB = false;
bool fromD = false;
- //bool zunder100 = part->originVertex()->position().Z() < 100.;
+ //bool zunder100 = part->originVertex()->position().Z() < 100.;
bool eta25 = !m_eta25cut || (part->momentum().Eta() > 2 && part->momentum().Eta() < 5);
@@ -458,33 +454,33 @@ StatusCode PrChecker::execute() {
}
}
while( 0 != mother ) {
- if ( mother->particleID().hasBottom() &&
+ if ( mother->particleID().hasBottom() &&
( mother->particleID().isMeson() || mother->particleID().isBaryon() ) ) fromB = true;
-
+
if ( mother->particleID().hasCharm() &&
( mother->particleID().isMeson() || mother->particleID().isBaryon() ) ) fromD = true;
mother = mother->originVertex()->mother();
}
}
-
+
// all ids associated to the mcparticle
- std::vector<LHCb::LHCbID> ids;
- if ( linkedIds.size() > (unsigned int) part->key() ) {
+ std::vector<LHCb::LHCbID> ids;
+ if ( linkedIds.size() > (unsigned int) part->key() ) {
for ( std::vector<int>::const_iterator itIm = linkedIds[part->key()].begin();
linkedIds[part->key()].end() != itIm; ++itIm ) {
LHCb::LHCbID temp;
temp.setDetectorType( (*itIm) >> 28 );
- // create LHCbId from int. Clumsy !
+ // create LHCbId from int. Clumsy !
temp.setID( (*itIm) );
ids.push_back( temp );
}
}
- if( msgLevel(MSG::VERBOSE) ) verbose() << "MCPart " << part->key()
+ if( msgLevel(MSG::VERBOSE) ) verbose() << "MCPart " << part->key()
<< " has " << ids.size() << " LHCbIDs " <<endmsg;
- //////////////////////////////////////
+ //////////////////////////////////////
std::vector<bool> flags;
//flags.push_back( true );
@@ -500,7 +496,7 @@ StatusCode PrChecker::execute() {
flags.push_back( eta25 && isLong && isInUT && fromB && trigger && !isElectron);
}
m_velo->countAndPlot(m_histoTool,part,flags,ids,nPrim);
-
+
flags.clear();
flags.push_back( eta25 && isLong && !isElectron );
flags.push_back( eta25 && isLong && over5 && !isElectron );
@@ -516,7 +512,7 @@ StatusCode PrChecker::execute() {
m_match->countAndPlot(m_histoTool,part,flags,ids,nPrim);
m_best->countAndPlot(m_histoTool,part,flags,ids,nPrim);
m_bestLong->countAndPlot(m_histoTool,part,flags,ids,nPrim);
-
+
flags.clear();
flags.push_back( eta25 && isInVelo && !isElectron);
flags.push_back( eta25 && isInVelo && isInUT && !isElectron );
@@ -533,7 +529,7 @@ StatusCode PrChecker::execute() {
flags.push_back( eta25 && isLong && isInUT && fromB && trigger && !isElectron);
}
m_upTrack->countAndPlot(m_histoTool,part,flags,ids,nPrim);
-
+
flags.clear();
flags.push_back( eta25 && trackInfo.hasT( part ) && !isElectron);
flags.push_back( eta25 && isLong && !isElectron);
@@ -549,7 +545,7 @@ StatusCode PrChecker::execute() {
flags.push_back( eta25 && strangeDown && !isInVelo && (fromB || fromD) && !isElectron);
flags.push_back( eta25 && strangeDown && !isInVelo && over5 && (fromB || fromD) && !isElectron);
m_tTrack->countAndPlot(m_histoTool,part,flags,ids,nPrim);
-
+
flags.clear();
flags.push_back( eta25 && isDown && !isElectron);
flags.push_back( eta25 && isDown && over5 && !isElectron);
@@ -567,33 +563,33 @@ StatusCode PrChecker::execute() {
flags.push_back( eta25 && strangeDown && !isInVelo && over5 && (fromB || fromD) && !isElectron);
m_downTrack->countAndPlot(m_histoTool,part,flags,ids,nPrim);
m_bestDownstream->countAndPlot(m_histoTool,part,flags,ids,nPrim);
-
+
// -- UT Hits counters (excluded for the moment)
/*
- flags.clear();
- flags.push_back( eta25 && isLong );
+ flags.clear();
+ flags.push_back( eta25 && isLong );
flags.push_back( eta25 && isLong && over5 );
flags.push_back( eta25 && isLong && over5 && isInUT );
if(m_triggerNumbers){
flags.push_back( eta25 && isLong && trigger );
}
-
+
flags.push_back( eta25 && isLong && fromB );
flags.push_back( eta25 && isLong && fromB && over5 );
flags.push_back( eta25 && isLong && fromB && over5 && isInUT);
-
+
if(m_triggerNumbers){
flags.push_back( eta25 && isLong && fromB && trigger );
}
-
+
m_utForward->countAndPlot( m_histoTool, part, flags, ids, nPrim );
m_utMatch->countAndPlot( m_histoTool, part, flags, ids, nPrim );
m_utbestLong->countAndPlot( m_histoTool, part, flags, ids, nPrim );
*/
}
-
+
return StatusCode::SUCCESS;
}
@@ -611,7 +607,7 @@ StatusCode PrChecker::finalize() {
}
//for ( std::vector<IPrUTCounter*>::iterator itCt = m_allUTCounters.begin();
- // m_allUTCounters.end() != itCt; ++itCt ) {
+ // m_allUTCounters.end() != itCt; ++itCt ) {
// (*itCt)->printStatistics();
//}
diff --git a/Pr/PrMCTools/src/PrChecker.h b/Pr/PrMCTools/src/PrChecker.h
index 4695efdc8cb0753f0e8cdff3d236f49b0eaa724b..c5b79ca77edb34b0137e81f7918529854343061b 100644
--- a/Pr/PrMCTools/src/PrChecker.h
+++ b/Pr/PrMCTools/src/PrChecker.h
@@ -39,12 +39,12 @@ class IHistoTool ;
* - cloneRate: nClones / (nRecod & rectible tracks + nClones)
* - purity: Fraction of hits on the track belonging to the matched MCParticle
* - hitEff: Fraction of hits on the track of all hits belonging to the matched MCParticle
- *
+ *
* Parameters:
* - Eta25Cut: Only consider particles with 2 < eta < 5? (default: false)
* - TriggerNumbers: Give numbers for p > 3GeV, pT > 500 MeV? (default: false)
* - UseElectrons: Take electrons into account in numbers? (default: false)
- *
+ *
* @author Olivier Callot, Thomas Nikodem
* @date 2014-02-13
*/
@@ -77,8 +77,6 @@ public:
/// Standard constructor
PrChecker( const std::string& name, ISvcLocator* pSvcLocator );
- virtual ~PrChecker( ); ///< Destructor
-
StatusCode initialize() override; ///< Algorithm initialization
StatusCode execute() override; ///< Algorithm execution
StatusCode finalize() override; ///< Algorithm finalization
@@ -94,7 +92,7 @@ private:
std::string m_downTracks;
std::string m_upTracks;
std::string m_bestTracks;
-
+
std::string m_idLoc;
IPrCounter* m_velo;
@@ -112,22 +110,22 @@ private:
//IPrCounter* m_utbestLong;
//IPrCounter* m_utDownst;
- int m_writeVeloHistos;
- int m_writeForwardHistos;
- int m_writeMatchHistos;
- int m_writeDownHistos;
- int m_writeUpHistos;
- int m_writeTTrackHistos;
+ int m_writeVeloHistos;
+ int m_writeForwardHistos;
+ int m_writeMatchHistos;
+ int m_writeDownHistos;
+ int m_writeUpHistos;
+ int m_writeTTrackHistos;
int m_writeBestHistos;
int m_writeBestLongHistos;
int m_writeBestDownstreamHistos;
int m_writeUTHistos;
- bool m_eta25cut;
+ bool m_eta25cut;
bool m_triggerNumbers;
bool m_useElectrons;
bool m_upgrade;
-
+
double m_ghostProbCut;
//== Vector of the counters
diff --git a/Pr/PrMCTools/src/PrChecker2.cpp b/Pr/PrMCTools/src/PrChecker2.cpp
index d96437dcceb2107bfbafb02047ade92e6b2a0aa5..cda39569ca858bf8780c4644600639ab87cd86c2 100644
--- a/Pr/PrMCTools/src/PrChecker2.cpp
+++ b/Pr/PrMCTools/src/PrChecker2.cpp
@@ -127,10 +127,6 @@ PrChecker2::PrChecker2( const std::string& name,
declareProperty( "TexOutputFolder", m_texfolder = "");
}
-//=============================================================================
-// Destructor
-//=============================================================================
-PrChecker2::~PrChecker2() {}
//=============================================================================
// Initialization
@@ -461,7 +457,7 @@ StatusCode PrChecker2::execute() {
}
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LHCb::MCVertices* mcVert = getIfExists<LHCb::MCVertices>(LHCb::MCVertexLocation::Default );
if( mcVert == nullptr ){
diff --git a/Pr/PrMCTools/src/PrChecker2.h b/Pr/PrMCTools/src/PrChecker2.h
index 597593da2b9dfa8c92d8573b2fb09236ceb05158..a8ba09080495296c9a330147ea1817d0aa64bf7c 100644
--- a/Pr/PrMCTools/src/PrChecker2.h
+++ b/Pr/PrMCTools/src/PrChecker2.h
@@ -1,5 +1,3 @@
-// $Id: $
-
#ifndef PRCHECKER2_H
#define PRCHECKER2_H 1
@@ -36,8 +34,8 @@ class IHistoTool ;
* - Upgrade: Gets hitLocation from link to MCParticle from "Pr/LHCbID" instead of "Pat/LHCbID" (default: true)
* - WriteTexOutput: Writes the statistics table to a TeX file (default: false)
* which is dumped to the location specified in TexOutputFolder
- *
- *
+ *
+ *
* additional track location container can be added for PrCounter2 and PrTTCounter via:
*
* @code
@@ -48,7 +46,7 @@ class IHistoTool ;
* myFactory.Modules = [ "LoKiMC.decorators" ]
* PrChecker2("PrChecker2").addTool( myFactory )
*
- * PrChecker2("PrChecker2").NewTracks= "Rec/Track/Velo"
+ * PrChecker2("PrChecker2").NewTracks= "Rec/Track/Velo"
* PrChecker2("PrChecker2").WriteNewHistos = 2
* PrChecker2("PrChecker2").SelectIdNewContainer = 1
* PrChecker2("PrChecker2").MyNewCuts= {"fromBorD" : "BOrDMother", "B daughters": "fromB", "D daughters": "fromD" }
@@ -62,12 +60,12 @@ class IHistoTool ;
* PrChecker2("PrChecker2").Velo.addTool(PrCounter2, name="Selector")
* PrChecker2("PrChecker2").Velo.Selector = "TrackSelector"
* PrChecker2("PrChecker2").Velo.addTool(TrackSelector("TrackSelector"))
- * PrChecker2("PrChecker2").Velo.TrackSelector.MaxChi2Cut = 5.0
+ * PrChecker2("PrChecker2").Velo.TrackSelector.MaxChi2Cut = 5.0
* @endcode
*
- * Uses IHistoTool: for each container histograms are plotted if the following code segment is used.
- * Default values are -1 (no histograms), can also be set to 2 for additional histograms (expectedHits, docaz, PVz, EtaP, EtaPhi, efficiency maps @z=9000mm XYZ9000 and @z=2485mm XYZ2485)
- *
+ * Uses IHistoTool: for each container histograms are plotted if the following code segment is used.
+ * Default values are -1 (no histograms), can also be set to 2 for additional histograms (expectedHits, docaz, PVz, EtaP, EtaPhi, efficiency maps @z=9000mm XYZ9000 and @z=2485mm XYZ2485)
+ *
* @code
* PrChecker2("PrChecker2").Write(container)Histos = 1
* @endcode
@@ -78,16 +76,16 @@ class IHistoTool ;
* PrChecker2("PrChecker2").MyForwardCuts = {"long" : "isLong", "B daughter": "fromB" }
* @endcode
*
- * As a default selection cuts of old PrChecker are used. The following cuts are predefined:
- * - is(Not)Long, is(Not)Velo, is(Not)Down, is(Not)Up, is(Not)TT, is(Not)Seed,
- * - fromB, fromD, BOrDMother, fromKsFromB, strange
+ * As a default selection cuts of old PrChecker are used. The following cuts are predefined:
+ * - is(Not)Long, is(Not)Velo, is(Not)Down, is(Not)Up, is(Not)TT, is(Not)Seed,
+ * - fromB, fromD, BOrDMother, fromKsFromB, strange
* - is(Not)Electron, eta25, over5, trigger
- *
+ *
* and LoKi syntax (LoKi::MCParticles) can be used for kinematical cuts, e.g. (MCPT> 2300), here the '()' are essential.
- *
- *
- * NB: If you care about the implementation: The cut-strings are converted into two types of functors:
- * - LoKi-type functors (hence all LoKi::MCParticles cuts work)
+ *
+ *
+ * NB: If you care about the implementation: The cut-strings are converted into two types of functors:
+ * - LoKi-type functors (hence all LoKi::MCParticles cuts work)
* - and custom-defined ones, mostly for type of reconstructibility and daughter-criteria (like 'isNotLong')
* where in the end all functors are evaluated on each MCParticle for each track container to define the reconstructibility.
* If a MCParticle is actually reconstructed is checked in the PrCounter2.
@@ -123,26 +121,24 @@ public:
class PrChecker2 : public GaudiHistoAlg {
public:
-
+
/// Standard constructor
PrChecker2( const std::string& name, ISvcLocator* pSvcLocator );
-
- virtual ~PrChecker2( ); ///< Destructor
-
+
StatusCode initialize() override; ///< Algorithm initialization
StatusCode execute() override; ///< Algorithm execution
StatusCode finalize() override; ///< Algorithm finalization
-
+
protected:
-
-
+
+
private:
// typedefs to make everything a bit more readable
typedef std::vector<std::string> strings;
typedef std::map<std::string, std::string> stringmap;
typedef std::pair<std::string, std::string> stringpair;
-
+
std::string m_veloTracks;
std::string m_forwardTracks;
@@ -153,7 +149,7 @@ private:
std::string m_bestTracks;
std::string m_newTracks;///< additional configurable track container for PrCounter2
std::string m_ttnewTracks;///< additional configurable track container for PrTTCounter
-
+
IPrCounter* m_velo;
IPrCounter* m_forward;
IPrCounter* m_match;
@@ -164,48 +160,48 @@ private:
IPrCounter* m_bestLong;
IPrCounter* m_bestDownstream;
IPrCounter* m_new;
-
+
IPrTTCounter* m_ttForward;
IPrTTCounter* m_ttMatch;
IPrTTCounter* m_ttDownst;
IPrTTCounter* m_ttnew;
-
+
//IPrCounter* m_utForward;
//IPrCounter* m_utMatch;
//IPrCounter* m_utbestLong;
//IPrCounter* m_utDownst;
-
- int m_writeVeloHistos;
- int m_writeForwardHistos;
- int m_writeMatchHistos;
- int m_writeDownHistos;
- int m_writeUpHistos;
- int m_writeTTrackHistos;
+
+ int m_writeVeloHistos;
+ int m_writeForwardHistos;
+ int m_writeMatchHistos;
+ int m_writeDownHistos;
+ int m_writeUpHistos;
+ int m_writeTTrackHistos;
int m_writeBestHistos;
int m_writeBestLongHistos;
int m_writeBestDownstreamHistos;
- int m_writeNewHistos;
+ int m_writeNewHistos;
int m_writeUTHistos;
int m_writeTTForwardHistos;
int m_writeTTMatchHistos;
int m_writeTTDownstHistos;
int m_writeTTNewHistos;
-
+
int m_selectIdNew;
int m_selectIdNewTT;
-
- bool m_eta25cut;
+
+ bool m_eta25cut;
double m_ghostProbCut;
bool m_triggerNumbers;
bool m_vetoElectrons;
bool m_trackextrapolation;
-
+
bool m_upgrade;
bool m_writetexfile;
std::string m_texfilename;
std::string m_texfolder;
-
+
enum recAs {
isLong = 1,
isNotLong = 2,
@@ -227,19 +223,19 @@ private:
isNotElectron = 19,
BOrDMother = 20,
};
-
+
std::map< std::string, recAs> m_lookuptable = {{"isLong",isLong},
{"isNotLong",isNotLong},
- {"isDown",isDown},
+ {"isDown",isDown},
{"isNotDown",isNotDown},
- {"isUp",isUp},
+ {"isUp",isUp},
{"isNotUp",isNotUp},
{"isVelo",isVelo},
- {"isNotVelo",isNotVelo},
+ {"isNotVelo",isNotVelo},
{"isTT",isTT},
{"isNoTT",isNotTT},
- {"isSeed",isSeed},
- {"isNotSeed",isNotSeed},
+ {"isSeed",isSeed},
+ {"isNotSeed",isNotSeed},
{"strange",strange},
{"fromB",fromB},
{"fromD",fromD},
@@ -253,31 +249,33 @@ private:
/** @class isTrack PrChecker2.h
* Predefined selection cuts: it converts strings to normal cuts, used by addOtherCuts
*/
- class isTrack {
-
+ class isTrack {
+
public:
- isTrack(const int kind) {
+ isTrack(const int kind) {
m_kind = kind;
};
/// Functor that checks if the MCParticle fulfills certain criteria, e.g. reco'ble as long track, B daughter, ...
bool operator()(LHCb::MCParticle* mcp, MCTrackInfo* mcInfo) const {
bool motherB = false;
bool motherD = false;
- if(m_kind == isLong) return mcInfo->hasVeloAndT( mcp );
- if(m_kind == isNotLong) return !mcInfo->hasVeloAndT( mcp );
- if(m_kind == isDown) return mcInfo->hasT( mcp ) && mcInfo->hasTT( mcp );
- if(m_kind == isNotDown) return !(mcInfo->hasT( mcp ) && mcInfo->hasTT( mcp ));
- if(m_kind == isUp) return mcInfo->hasVelo( mcp ) && mcInfo->hasTT( mcp );
- if(m_kind == isNotUp) return !(mcInfo->hasVelo( mcp ) && mcInfo->hasTT( mcp ));
- if(m_kind == isVelo) return mcInfo->hasVelo( mcp );
- if(m_kind == isNotVelo) return !mcInfo->hasVelo( mcp );
- if(m_kind == isSeed) return mcInfo->hasT( mcp );
- if(m_kind == isNotSeed) return !mcInfo->hasT( mcp );
- if(m_kind == isTT) return mcInfo->hasTT( mcp );
- if(m_kind == isNotTT) return !mcInfo->hasTT( mcp );
- if(m_kind == isElectron) return abs( mcp->particleID().pid() ) == 11;
- if(m_kind == isNotElectron) return abs( mcp->particleID().pid() ) != 11;
-
+ switch(m_kind) {
+ case isLong : return mcInfo->hasVeloAndT( mcp );
+ case isNotLong : return !mcInfo->hasVeloAndT( mcp );
+ case isDown : return mcInfo->hasT( mcp ) && mcInfo->hasTT( mcp );
+ case isNotDown : return !(mcInfo->hasT( mcp ) && mcInfo->hasTT( mcp ));
+ case isUp : return mcInfo->hasVelo( mcp ) && mcInfo->hasTT( mcp );
+ case isNotUp : return !(mcInfo->hasVelo( mcp ) && mcInfo->hasTT( mcp ));
+ case isVelo : return mcInfo->hasVelo( mcp );
+ case isNotVelo : return !mcInfo->hasVelo( mcp );
+ case isSeed : return mcInfo->hasT( mcp );
+ case isNotSeed : return !mcInfo->hasT( mcp );
+ case isTT : return mcInfo->hasTT( mcp );
+ case isNotTT : return !mcInfo->hasTT( mcp );
+ case isElectron : return abs( mcp->particleID().pid() ) == 11;
+ case isNotElectron : return abs( mcp->particleID().pid() ) != 11;
+ }
+
if ( 0 != mcp->originVertex() ) {
const LHCb::MCParticle* mother = mcp->originVertex()->mother();
if ( 0 != mother ) {
@@ -297,14 +295,14 @@ private:
3334 == pid // Omega-
) {
if(m_kind == strange) return true;
-
+
}
// -- It's a Kshort from a b Hadron
if ( 0 != mother->originVertex()->mother() ) {
if ( 310 == pid &&
2 == mcp->originVertex()->products().size() &&
mother->originVertex()->mother()->particleID().hasBottom() &&
- ( mother->originVertex()->mother()->particleID().isMeson() ||
+ ( mother->originVertex()->mother()->particleID().isMeson() ||
mother->originVertex()->mother()->particleID().isBaryon() )) {
if(m_kind == fromKsFromB) return true;
}
@@ -314,23 +312,23 @@ private:
}
// -- It's a daughter of a B or D hadron
while( 0 != mother ) {
- if ( mother->particleID().hasBottom() &&
+ if ( mother->particleID().hasBottom() &&
( mother->particleID().isMeson() || mother->particleID().isBaryon() ) ) motherB = true;
-
+
if ( mother->particleID().hasCharm() &&
( mother->particleID().isMeson() || mother->particleID().isBaryon() ) ) motherD = true;
-
+
mother = mother->originVertex()->mother();
-
+
}
if( m_kind == fromD && motherD == true) return m_kind == fromD;
if( m_kind == fromB && motherB == true) return m_kind == fromB;
-
+
if(m_kind == BOrDMother && (motherD || motherB) ) return m_kind == BOrDMother;
}
return false;
}
-
+
private:
int m_kind;
};
@@ -340,27 +338,27 @@ private:
* Class that adds selection cuts defined in isTrack to cuts
*/
class addOtherCuts{
-
+
public:
void addCut(isTrack cat){ m_cuts.push_back(cat); }
-
+
/// Functor that evaluates all 'isTrack' cuts
bool operator()(LHCb::MCParticle* mcp, MCTrackInfo* mcInfo){
-
+
bool decision = true;
-
+
for(unsigned int i = 0; i < m_cuts.size(); ++i){
decision = decision && m_cuts[i](mcp, mcInfo);
}
-
+
return decision;
}
-
+
private:
std::vector<isTrack> m_cuts;
-
+
};
-
+
//maps for each track container with {cut name,selection cut}
stringmap m_map_forward;
stringmap m_map_velo;
@@ -371,21 +369,21 @@ private:
stringmap m_map_ttdown;
stringmap m_map_new;
stringmap m_map_ttnew;
-
+
//default cuts for each container, if not other ones are specified in run python file, those are taken
/** @brief Map filled with default cuts for each container, first string is name of the cut, second one is cut (predefined from class isTrack or kinematical Loki Cuts)
*/
stringmap DefaultCutMap( std::string name ){
-
+
stringmap map_velo = { {"01_velo" ,"isVelo "},//use numbers for right order of cuts
{"02_long","isLong "},
- {"03_long>5GeV" ,"isLong & over5 " },
+ {"03_long>5GeV" ,"isLong & over5 " },
{"04_long_strange" , "isLong & strange " },
{"05_long_strange>5GeV","isLong & strange & over5 " },
{"06_long_fromB" ,"isLong & fromB " },
{"07_long_fromB>5GeV" , "isLong & fromB & over5 " } };
- stringmap map_forward = { {"01_long","isLong "},
- {"02_long>5GeV","isLong & over5 " },
+ stringmap map_forward = { {"01_long","isLong "},
+ {"02_long>5GeV","isLong & over5 " },
{"03_long_strange","isLong & strange " },
{"04_long_strange>5GeV","isLong & strange & over5 " },
{"05_long_fromB","isLong & fromB " },
@@ -426,14 +424,14 @@ private:
{"12_UT+T_SfromDB>5GeV" , " strange & isDown & over5 & ( fromB | fromD ) "},
{"13_noVelo+UT+T_SfromDB" , " strange & isDown & isNotVelo & ( fromB | fromD )"},
{"14_noVelo+UT+T_SfromDB>5GeV" , " strange & isDown & isNotVelo & over5 & ( fromB | fromD ) "} };
-
+
stringmap map_new = { };
stringmap map_ttforward = { { "01_long" ,"isLong "},
{"02_long>5GeV" ,"isLong & over5 "} };
stringmap map_ttdown = { {"01_has seed" ,"isSeed "},
{"02_has seed +noVelo, T+TT" ,"isSeed & isNotVelo & isDown "},
{"03_down+strange" , " strange & isDown"},
- {"04_down+strange+>5GeV" , " strange & isDown & over5 "},
+ {"04_down+strange+>5GeV" , " strange & isDown & over5 "},
{"05_pi<-Ks<-B" , "fromKsFromB "},
{"06_pi<-Ks<-B+> 5 GeV" , "fromKsFromB & over5 "} };
stringmap map_ttnew = { };
@@ -448,11 +446,11 @@ private:
if( name == "TTForward" ) return map_ttforward;
if( name == "TTDownstream" ) return map_ttdown;
if( name == "TTNew" ) return map_ttnew;
-
+
return empty_map;
}
-
+
/** @brief makes vector of second elements of DefaultCutMap --> needed as input for m_Cuts */
strings getMyCut( stringmap myCutMap ){
strings dummy;
@@ -462,27 +460,27 @@ private:
return dummy;
}
-
+
//== Vector of the counters
std::vector<IPrCounter*> m_allCounters;///<Vector of PrCounter
std::vector<IPrTTCounter*> m_allTTCounters;///<Vector of PrTTCounter
-
+
const IHistoTool* m_histoTool;
LoKi::IMCHybridFactory* m_factory;///<needed to convert normal cuts into Loki cuts
// maps for cuts
std::map< std::string, strings > m_Cuts;///<map of track container name and corresponding cuts
-
+
std::map < std::string, std::vector < LoKi::Types::MCCut> > m_LoKiCuts;///<converted map of Loki cuts, first component is name of track container
- std::map < std::string, std::vector <addOtherCuts> > m_otherCuts;///<map of other cuts as predefined in isTrack, first component is name of track container
-
- std::vector<LoKi::Types::MCCut> m_veloMCCuts;
- std::vector<addOtherCuts> m_veloMCCuts2;
+ std::map < std::string, std::vector <addOtherCuts> > m_otherCuts;///<map of other cuts as predefined in isTrack, first component is name of track container
+
+ std::vector<LoKi::Types::MCCut> m_veloMCCuts;
+ std::vector<addOtherCuts> m_veloMCCuts2;
std::vector<LoKi::Types::MCCut> m_forwardMCCuts;
- std::vector<addOtherCuts> m_forwardMCCuts2;
- std::vector<LoKi::Types::MCCut> m_upstreamMCCuts;
+ std::vector<addOtherCuts> m_forwardMCCuts2;
+ std::vector<LoKi::Types::MCCut> m_upstreamMCCuts;
std::vector<addOtherCuts> m_upstreamMCCuts2;
std::vector<LoKi::Types::MCCut> m_ttrackMCCuts;
std::vector<addOtherCuts> m_ttrackMCCuts2;
diff --git a/Pr/PrMCTools/src/PrClustersResidual.cpp b/Pr/PrMCTools/src/PrClustersResidual.cpp
index 4dd64323ffcfddfa233e8dd351e417c3d4014cf9..bf1a2b06be181424b0ed2cb5a3d30248f279f7ae 100644
--- a/Pr/PrMCTools/src/PrClustersResidual.cpp
+++ b/Pr/PrMCTools/src/PrClustersResidual.cpp
@@ -34,10 +34,6 @@ PrClustersResidual::PrClustersResidual( const std::string& name,
declareProperty("FTHitsLocation", m_HitsInTES);
declareProperty("FTLiteClusterLocation", m_HitsInTES);
}
-//=============================================================================
-// Destructor
-//=============================================================================
-PrClustersResidual::~PrClustersResidual() {}
//=============================================================================
// Initialization
@@ -113,7 +109,7 @@ void PrClustersResidual::Occupancy(const PrFTHitHandler<PrHit>& FTHitHandler){
char layerName[100];
char Title[100];
std::vector<int> nHits(12,0);
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LHCb::MCVertices* mcVert = getIfExists<LHCb::MCVertices>(LHCb::MCVertexLocation::Default );
if( mcVert == nullptr ){
error() << "Could not find MCVertices at " << LHCb::MCParticleLocation::Default << endmsg;
@@ -310,7 +306,7 @@ void PrClustersResidual::TrackStudy(const PrFTHitHandler<PrHit>& FTHitHandler){
//Get the McParticles
LHCb::MCParticles* mcParts = getIfExists<LHCb::MCParticles> ( LHCb::MCParticleLocation::Default );
//Get track informations
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
//If you Run the cheated Seeding or the normal seeding or the forward , you don't have to decode it again:
debug()<<"Initialize"<<endmsg;
@@ -966,7 +962,7 @@ void PrClustersResidual::ClusterResidual(const PrFTHitHandler<PrHit>& FTHitHand
if(msgLevel(MSG::DEBUG)) debug()<<" Get Cluster to MCHit Linker "<<endmsg;
LinkedTo<LHCb::MCHit> myFTCluster2MCHitLink ( evtSvc(),msgSvc(), LHCb::FTLiteClusterLocation::Default + "2MCHits");
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
if(msgLevel(MSG::DEBUG)) debug()<<" Add MCHitAndTrackStudy "<<endmsg;
Tuples::Tuple tuple = GaudiTupleAlg::nTuple("ClusterMCHitAndTrackStudy","Events");
for (unsigned int zone = 0; zone < m_zone; ++zone)
diff --git a/Pr/PrMCTools/src/PrClustersResidual.h b/Pr/PrMCTools/src/PrClustersResidual.h
index 9b98cc47eb5233d7f7bca21635b961015fe30f87..4010106b423b99352f56557f12c5c05a80b96ba0 100644
--- a/Pr/PrMCTools/src/PrClustersResidual.h
+++ b/Pr/PrMCTools/src/PrClustersResidual.h
@@ -65,14 +65,10 @@ public:
/// Standard constructor
PrClustersResidual( const std::string& name, ISvcLocator* pSvcLocator );
- virtual ~PrClustersResidual( ); ///< Destructor
-
StatusCode initialize() override; ///< Algorithm initialization
StatusCode execute() override; ///< Algorithm execution
StatusCode finalize() override; ///< Algorithm finalization
-protected:
-
private:
void HitEfficiency();
diff --git a/Pr/PrMCTools/src/PrDebugMatchToolNN.cpp b/Pr/PrMCTools/src/PrDebugMatchToolNN.cpp
index 4689cb077563c3cc1976357ee5b0b9705fe1ff10..c4aa74256c6fdb9e410e6062b86e0c24b3abf9fb 100644
--- a/Pr/PrMCTools/src/PrDebugMatchToolNN.cpp
+++ b/Pr/PrMCTools/src/PrDebugMatchToolNN.cpp
@@ -11,7 +11,6 @@
#include "Event/MCTrackInfo.h"
#include "Event/MCParticle.h"
#include "Event/Track.h"
-#include "Event/MCTrackInfo.h"
#include "Linker/LinkedTo.h"
@@ -27,44 +26,44 @@ DECLARE_COMPONENT( PrDebugMatchToolNN )
//=============================================================================
void PrDebugMatchToolNN::fillTuple(const LHCb::Track& velo, const LHCb::Track& seed, const std::vector<float>& vars ) const
{
-
+
int found= matchMCPart(velo, seed);
Tuple tuple = nTuple("tuple","tuple");
-
+
tuple->column("quality", found ).ignore();
-
+
unsigned int i=0;
for(auto var: vars ){
tuple->column("var"+std::to_string(i), var).ignore();
i++;
}
-
+
tuple->write().ignore();
}
//=============================================================================
-int PrDebugMatchToolNN::matchMCPart(const LHCb::Track& velo, const LHCb::Track& seed) const
+int PrDebugMatchToolNN::matchMCPart(const LHCb::Track& velo, const LHCb::Track& seed) const
{
-
+
LinkedTo<LHCb::MCParticle, LHCb::Track> myLinkVelo ( evtSvc(), msgSvc(),LHCb::TrackLocation::Velo);
LinkedTo<LHCb::MCParticle, LHCb::Track> myLinkSeed ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
- LHCb::MCParticle* mcPartV = myLinkVelo.first( velo.key() );
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ LHCb::MCParticle* mcPartV = myLinkVelo.first( velo.key() );
+ const MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
int found=0;
while(mcPartV != NULL){
- if( !trackInfo.hasVeloAndT( mcPartV)){
+ if( !trackInfo.hasVeloAndT( mcPartV)){
mcPartV=myLinkVelo.next();
- continue;
+ continue;
}
-
- LHCb::MCParticle* mcPartS = myLinkSeed.first( seed.key() );
+
+ LHCb::MCParticle* mcPartS = myLinkSeed.first( seed.key() );
while(mcPartS != NULL){
if(mcPartV == mcPartS){
-
- if(11 == abs( mcPartV->particleID().pid())){
+
+ if(11 == abs( mcPartV->particleID().pid())){
found = -1;
break;
}
diff --git a/Pr/PrMCTools/src/PrDebugTrackingLosses.cpp b/Pr/PrMCTools/src/PrDebugTrackingLosses.cpp
index 31d94e2605192c0d903768a25f5bec8b476d28f5..bbc2e7ea68f1cc7cccc0d2446981fcfd53f31d10 100644
--- a/Pr/PrMCTools/src/PrDebugTrackingLosses.cpp
+++ b/Pr/PrMCTools/src/PrDebugTrackingLosses.cpp
@@ -83,7 +83,7 @@ StatusCode PrDebugTrackingLosses::execute() {
LinkedFrom<LHCb::Track,LHCb::MCParticle> veloLinker ( evtSvc(), msgSvc(), m_veloName );
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ const MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LHCb::MCParticles::const_iterator itP = partCont->begin();
for (; itP != partCont->end(); ++itP){
diff --git a/Pr/PrMCTools/src/PrDebugUTTruthTool.cpp b/Pr/PrMCTools/src/PrDebugUTTruthTool.cpp
index 75079ac54335e00f880b50ead1bef60d420f23c5..b2d9a0645943974f7179bf9b751b999c896dd3a5 100644
--- a/Pr/PrMCTools/src/PrDebugUTTruthTool.cpp
+++ b/Pr/PrMCTools/src/PrDebugUTTruthTool.cpp
@@ -297,7 +297,7 @@ void PrDebugUTTruthTool::forceMCHits(UT::Mut::Hits& hits, LHCb::Track* track){
//=========================================================================
void PrDebugUTTruthTool::tuneFisher( const LHCb::Track* seedTrack){
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ const MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
bool goodTrack = true;
@@ -357,7 +357,7 @@ void PrDebugUTTruthTool::tuneFisher( const LHCb::Track* seedTrack){
//=========================================================================
void PrDebugUTTruthTool::tuneDeltaP( const LHCb::Track* seedTrack, const double deltaP, const double momentum){
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ const MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
bool goodTrack = true;
@@ -395,7 +395,7 @@ void PrDebugUTTruthTool::tuneDeltaP( const LHCb::Track* seedTrack, const double
//=========================================================================
void PrDebugUTTruthTool::tuneFinalMVA( const LHCb::Track* seedTrack, bool goodTrack, std::vector<double> vals){
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ const MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
bool goodSeedTrack = true;
@@ -451,7 +451,7 @@ void PrDebugUTTruthTool::tuneFinalMVA( const LHCb::Track* seedTrack, bool goodTr
//=========================================================================
void PrDebugUTTruthTool::getMagnetError( const LHCb::Track* seedTrack){
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ const MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
bool goodSeedTrack = true;
diff --git a/Pr/PrMCTools/src/PrPlotFTHits.cpp b/Pr/PrMCTools/src/PrPlotFTHits.cpp
index d3132c725f291fa356badb9c46cf9a28fe6de3f7..8171421efcaee779e1c07dbad4c42780633b6deb 100644
--- a/Pr/PrMCTools/src/PrPlotFTHits.cpp
+++ b/Pr/PrMCTools/src/PrPlotFTHits.cpp
@@ -113,7 +113,7 @@ void PrPlotFTHits::plotOccupancy(const PrFTHitHandler<PrHit>& FTHitHandler){
char layerName[100];
char Title[100];
std::vector<int> nHits(12,0);
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ const MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LHCb::MCVertices* mcVert = getIfExists<LHCb::MCVertices>(LHCb::MCVertexLocation::Default );
if( mcVert == nullptr ){
@@ -249,7 +249,6 @@ void PrPlotFTHits::plotState(){
if( !forwardTracks ) info() << "No tracks found in: " << LHCb::TrackLocation::Forward << endmsg;
if( !seedTracks ) info() << "No tracks found in: " << LHCb::TrackLocation::Seed << endmsg;
char stateZ[100];
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
if( seedTracks)
{
@@ -325,7 +324,7 @@ void PrPlotFTHits::plotHitEfficiency(const PrFTHitHandler<PrHit>& FTHitHandler){
LinkedTo<LHCb::MCParticle, LHCb::Track> myForwardLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Forward);
LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default );
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ const MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
// -----------------------------------------------------------------------------------
// Plot hit efficiency as function of x and y
@@ -614,7 +613,7 @@ void PrPlotFTHits::plotTrackingEfficiency(){
LHCb::MCParticles* mcParts = getIfExists<LHCb::MCParticles>( LHCb::MCParticleLocation::Default );
if ( msgLevel( MSG::ERROR ) && !mcParts ) error() << "Could not find MCParticles at: " << LHCb::MCParticleLocation::Default << endmsg;
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ const MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
for(const LHCb::MCParticle* mcPart: *mcParts){
@@ -714,7 +713,7 @@ void PrPlotFTHits::plotMCHits ( const PrFTHitHandler<PrHit>& FTHitHandler) {
LinkedTo<LHCb::MCHit> myMCHitLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default + "2MCHits");
LinkedTo<LHCb::MCParticle> myClusterLink ( evtSvc(), msgSvc(), LHCb::FTLiteClusterLocation::Default );
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
for(unsigned int i = 0; i < m_zone; i++){
diff --git a/Tr/ParameterizedKalman/src/CompareTracks.cpp b/Tr/ParameterizedKalman/src/CompareTracks.cpp
index ebb275a3f59776e60cce1dc773ec0f82109d2a4b..1f5447cd78c678d95bb50a7039cae61c0b40dd53 100644
--- a/Tr/ParameterizedKalman/src/CompareTracks.cpp
+++ b/Tr/ParameterizedKalman/src/CompareTracks.cpp
@@ -1,4 +1,4 @@
-// Include files
+// Include files
// from Gaudi
#include "Event/State.h"
@@ -47,23 +47,23 @@ void CompareTracks::operator()(const LHCb::Tracks& tracks1, const LHCb::Tracks&
const LHCb::ODIN& odin, const LHCb::LinksByKey& links) const {
if ( msgLevel(MSG::DEBUG) ) debug() << "==> Execute" << endmsg;
-
+
if(tracks1.empty() || tracks2.empty()){
return;
}
-
+
//Create output tuple
- //Create a new file for every event in order to be threadsafe
+ //Create a new file for every event in order to be threadsafe
TFile outputFile((m_FileName.toString()+"_"+std::to_string(odin.runNumber())+
"_"+std::to_string(odin.eventNumber())+".root").c_str(),"RECREATE");
//create the varaibles to be filled
tupleVars vars;
-
+
//create the tree
TTree tree("compare","compare");
addBranches(tree, &vars);
-
+
//Loop over tracks1
for(auto const &track1 : tracks1){
//Search for the respective track in tracks2
@@ -72,7 +72,7 @@ void CompareTracks::operator()(const LHCb::Tracks& tracks1, const LHCb::Tracks&
if(i==std::end(tracks2)){
warning() << "No matching track in second container found" << endmsg;
return;
- }
+ }
//compare the tracks and fill the tuple
@@ -87,7 +87,7 @@ void CompareTracks::operator()(const LHCb::Tracks& tracks1, const LHCb::Tracks&
FillNtuple(*track1, *(*(i)), links, &vars, zBeam, 0, true);
//2. ...
-
+
//3. ...
tree.Fill();
@@ -95,7 +95,7 @@ void CompareTracks::operator()(const LHCb::Tracks& tracks1, const LHCb::Tracks&
outputFile.Write();
outputFile.Close();
-
+
return;
}
@@ -103,7 +103,7 @@ void CompareTracks::operator()(const LHCb::Tracks& tracks1, const LHCb::Tracks&
// Add branches to the tuple that will be filled with the comparison information
//==================================================================================================
void CompareTracks::addBranches(TTree &tree, tupleVars *treeVars) const {
-
+
//Set the branches
tree.Branch("sF_sigmaxx_T" , &(treeVars->m_sF_P[2][0 ]) , "sF_sigmaxx_T/D");
tree.Branch("sF_sigmayy_T" , &(treeVars->m_sF_P[2][2 ]) , "sF_sigmayy_T/D");
@@ -126,7 +126,7 @@ void CompareTracks::addBranches(TTree &tree, tupleVars *treeVars) const {
tree.Branch("sF_tx_T" , &(treeVars->m_sF_x[2][2]) , "sF_tx_T/D");
tree.Branch("sF_ty_T" , &(treeVars->m_sF_x[2][3]) , "sF_ty_T/D");
tree.Branch("sF_qop_T" , &(treeVars->m_sF_x[2][4]) , "sF_qop_T/D");
-
+
tree.Branch("sF_z_T" , &(treeVars->m_sF_z[2]) , "sF_z_T/D");
tree.Branch("sF_true_x_T" , &(treeVars->m_sF_true_x[2][0]) , "sF_true_x_T/D");
@@ -134,7 +134,7 @@ void CompareTracks::addBranches(TTree &tree, tupleVars *treeVars) const {
tree.Branch("sF_true_tx_T" , &(treeVars->m_sF_true_x[2][2]) , "sF_true_tx_T/D");
tree.Branch("sF_true_ty_T" , &(treeVars->m_sF_true_x[2][3]) , "sF_true_ty_T/D");
tree.Branch("sF_true_qop_T" , &(treeVars->m_sF_true_x[2][4]) , "sF_true_qop_T/D");
-
+
tree.Branch("sF_sigmaxx_V" , &(treeVars->m_sF_P[0][0 ]) , "sF_sigmaxx_V/D");
tree.Branch("sF_sigmayy_V" , &(treeVars->m_sF_P[0][2 ]) , "sF_sigmayy_V/D");
tree.Branch("sF_sigmatxtx_V" , &(treeVars->m_sF_P[0][5 ]) , "sF_sigmatxtx_V/D");
@@ -156,7 +156,7 @@ void CompareTracks::addBranches(TTree &tree, tupleVars *treeVars) const {
tree.Branch("sF_tx_V" , &(treeVars->m_sF_x[0][2]) , "sF_tx_V/D");
tree.Branch("sF_ty_V" , &(treeVars->m_sF_x[0][3]) , "sF_ty_V/D");
tree.Branch("sF_qop_V" , &(treeVars->m_sF_x[0][4]) , "sF_qop_V/D");
-
+
tree.Branch("sF_z_V" , &(treeVars->m_sF_z[0]) , "sF_z_V/D");
tree.Branch("sF_true_x_V" , &(treeVars->m_sF_true_x[0][0]) , "sF_true_x_V/D");
@@ -164,7 +164,7 @@ void CompareTracks::addBranches(TTree &tree, tupleVars *treeVars) const {
tree.Branch("sF_true_tx_V" , &(treeVars->m_sF_true_x[0][2]) , "sF_true_tx_V/D");
tree.Branch("sF_true_ty_V" , &(treeVars->m_sF_true_x[0][3]) , "sF_true_ty_V/D");
tree.Branch("sF_true_qop_V" , &(treeVars->m_sF_true_x[0][4]) , "sF_true_qop_V/D");
-
+
tree.Branch("sF_sigmaxx" , &(treeVars->m_sF_P[1][0 ]) , "sF_sigmaxx/D");
tree.Branch("sF_sigmayy" , &(treeVars->m_sF_P[1][2 ]) , "sF_sigmayy/D");
tree.Branch("sF_sigmatxtx" , &(treeVars->m_sF_P[1][5 ]) , "sF_sigmatxtx/D");
@@ -186,7 +186,7 @@ void CompareTracks::addBranches(TTree &tree, tupleVars *treeVars) const {
tree.Branch("sF_tx" , &(treeVars->m_sF_x[1][2]) , "sF_tx/D");
tree.Branch("sF_ty" , &(treeVars->m_sF_x[1][3]) , "sF_ty/D");
tree.Branch("sF_qop" , &(treeVars->m_sF_x[1][4]) , "sF_qop/D");
-
+
tree.Branch("sF_z" , &(treeVars->m_sF_z[1]) , "sF_z/D");
tree.Branch("sF_chi2" , &(treeVars->m_sF_chi2) ,"sF_chi2/D");
tree.Branch("sF_ndof" , &(treeVars->m_sF_ndof) ,"sF_ndof/D");
@@ -196,7 +196,7 @@ void CompareTracks::addBranches(TTree &tree, tupleVars *treeVars) const {
tree.Branch("sF_true_tx" , &(treeVars->m_sF_true_x[1][2]) , "sF_true_tx/D");
tree.Branch("sF_true_ty" , &(treeVars->m_sF_true_x[1][3]) , "sF_true_ty/D");
tree.Branch("sF_true_qop" , &(treeVars->m_sF_true_x[1][4]) , "sF_true_qop/D");
-
+
tree.Branch("dF_sigmaxx_T" , &(treeVars->m_dF_P[2][0 ]) , "dF_sigmaxx_T/D");
tree.Branch("dF_sigmayy_T" , &(treeVars->m_dF_P[2][2 ]) , "dF_sigmayy_T/D");
tree.Branch("dF_sigmatxtx_T" , &(treeVars->m_dF_P[2][5 ]) , "dF_sigmatxtx_T/D");
@@ -218,7 +218,7 @@ void CompareTracks::addBranches(TTree &tree, tupleVars *treeVars) const {
tree.Branch("dF_tx_T" , &(treeVars->m_dF_x[2][2]) , "dF_tx_T/D");
tree.Branch("dF_ty_T" , &(treeVars->m_dF_x[2][3]) , "dF_ty_T/D");
tree.Branch("dF_qop_T" , &(treeVars->m_dF_x[2][4]) , "dF_qop_T/D");
-
+
tree.Branch("dF_z_T" , &(treeVars->m_dF_z[2]) , "dF_z_T/D");
tree.Branch("dF_true_x_T" , &(treeVars->m_dF_true_x[2][0]) , "dF_true_x_T/D");
@@ -226,7 +226,7 @@ void CompareTracks::addBranches(TTree &tree, tupleVars *treeVars) const {
tree.Branch("dF_true_tx_T" , &(treeVars->m_dF_true_x[2][2]) , "dF_true_tx_T/D");
tree.Branch("dF_true_ty_T" , &(treeVars->m_dF_true_x[2][3]) , "dF_true_ty_T/D");
tree.Branch("dF_true_qop_T" , &(treeVars->m_dF_true_x[2][4]) , "dF_true_qop_T/D");
-
+
tree.Branch("dF_sigmaxx_V" , &(treeVars->m_dF_P[0][0 ]) , "dF_sigmaxx_V/D");
tree.Branch("dF_sigmayy_V" , &(treeVars->m_dF_P[0][2 ]) , "dF_sigmayy_V/D");
tree.Branch("dF_sigmatxtx_V" , &(treeVars->m_dF_P[0][5 ]) , "dF_sigmatxtx_V/D");
@@ -248,7 +248,7 @@ void CompareTracks::addBranches(TTree &tree, tupleVars *treeVars) const {
tree.Branch("dF_tx_V" , &(treeVars->m_dF_x[0][2]) , "dF_tx_V/D");
tree.Branch("dF_ty_V" , &(treeVars->m_dF_x[0][3]) , "dF_ty_V/D");
tree.Branch("dF_qop_V" , &(treeVars->m_dF_x[0][4]) , "dF_qop_V/D");
-
+
tree.Branch("dF_z_V" , &(treeVars->m_dF_z[0]) , "dF_z_V/D");
tree.Branch("dF_true_x_V" , &(treeVars->m_dF_true_x[0][0]) , "dF_true_x_V/D");
@@ -256,7 +256,7 @@ void CompareTracks::addBranches(TTree &tree, tupleVars *treeVars) const {
tree.Branch("dF_true_tx_V" , &(treeVars->m_dF_true_x[0][2]) , "dF_true_tx_V/D");
tree.Branch("dF_true_ty_V" , &(treeVars->m_dF_true_x[0][3]) , "dF_true_ty_V/D");
tree.Branch("dF_true_qop_V" , &(treeVars->m_dF_true_x[0][4]) , "dF_true_qop_V/D");
-
+
tree.Branch("dF_sigmaxx" , &(treeVars->m_dF_P[1][0 ]) , "dF_sigmaxx/D");
tree.Branch("dF_sigmayy" , &(treeVars->m_dF_P[1][2 ]) , "dF_sigmayy/D");
tree.Branch("dF_sigmatxtx" , &(treeVars->m_dF_P[1][5 ]) , "dF_sigmatxtx/D");
@@ -278,7 +278,7 @@ void CompareTracks::addBranches(TTree &tree, tupleVars *treeVars) const {
tree.Branch("dF_tx" , &(treeVars->m_dF_x[1][2]) , "dF_tx/D");
tree.Branch("dF_ty" , &(treeVars->m_dF_x[1][3]) , "dF_ty/D");
tree.Branch("dF_qop" , &(treeVars->m_dF_x[1][4]) , "dF_qop/D");
-
+
tree.Branch("dF_z" , &(treeVars->m_dF_z[1]) , "dF_z/D");
tree.Branch("dF_chi2" , &(treeVars->m_dF_chi2) ,"dF_chi2/D");
tree.Branch("dF_ndof" , &(treeVars->m_dF_ndof) ,"dF_ndof/D");
@@ -288,9 +288,9 @@ void CompareTracks::addBranches(TTree &tree, tupleVars *treeVars) const {
tree.Branch("dF_true_tx" , &(treeVars->m_dF_true_x[1][2]) , "dF_true_tx/D");
tree.Branch("dF_true_ty" , &(treeVars->m_dF_true_x[1][3]) , "dF_true_ty/D");
tree.Branch("dF_true_qop" , &(treeVars->m_dF_true_x[1][4]) , "dF_true_qop/D");
-
+
tree.Branch("true_qop_vertex" , &(treeVars->m_true_qop_vertex) , "true_qop_vertex/D");
-
+
tree.Branch("MCstatus" , &(treeVars->m_MC_status) , "MCstatus/I");
}
@@ -312,10 +312,10 @@ void CompareTracks::FillNtuple(const LHCb::Track &track1, const LHCb::Track &tra
//covariance matrices
Gaudi::TrackSymMatrix covMat1;
- covMat1 = state1->covariance();
-
+ covMat1 = state1->covariance();
+
Gaudi::TrackSymMatrix covMat2;
- covMat2 = state2->covariance();
+ covMat2 = state2->covariance();
if(!sc1.isSuccess() || !sc2.isSuccess()) return;
@@ -332,8 +332,8 @@ void CompareTracks::FillNtuple(const LHCb::Track &track1, const LHCb::Track &tra
vars->m_sF_x[nPos][2] = state1->tx();
vars->m_sF_x[nPos][3] = state1->ty();
vars->m_sF_x[nPos][4] = state1->qOverP();
-
- vars->m_sF_true_x[nPos][0] = trueX;
+
+ vars->m_sF_true_x[nPos][0] = trueX;
vars->m_sF_true_x[nPos][1] = trueY;
vars->m_sF_true_x[nPos][2] = trueTX;
vars->m_sF_true_x[nPos][3] = trueTY;
@@ -344,8 +344,8 @@ void CompareTracks::FillNtuple(const LHCb::Track &track1, const LHCb::Track &tra
vars->m_dF_x[nPos][2] = state2->tx();
vars->m_dF_x[nPos][3] = state2->ty();
vars->m_dF_x[nPos][4] = state2->qOverP();
-
- vars->m_dF_true_x[nPos][0] = trueX;
+
+ vars->m_dF_true_x[nPos][0] = trueX;
vars->m_dF_true_x[nPos][1] = trueY;
vars->m_dF_true_x[nPos][2] = trueTX;
vars->m_dF_true_x[nPos][3] = trueTY;
@@ -356,23 +356,23 @@ void CompareTracks::FillNtuple(const LHCb::Track &track1, const LHCb::Track &tra
vars->m_sF_z[nPos] = z;
vars->m_dF_z[nPos] = z;
-
+
int k =0;
for(int i=0;i<5;i++){
for(int j=0;j<=i;j++){
- vars->m_sF_P[nPos][k] = covMat1(i,j);
- vars->m_dF_P[nPos][k] = covMat2(i,j);
+ vars->m_sF_P[nPos][k] = covMat1(i,j);
+ vars->m_dF_P[nPos][k] = covMat2(i,j);
k++;
}
}
- //Set other track varaibles
- vars->m_sF_ndof = track1.nDoF();
- vars->m_sF_chi2 = track1.chi2();
-
- vars->m_dF_ndof = track2.nDoF();
- vars->m_dF_chi2 = track2.chi2();
-
+ //Set other track varaibles
+ vars->m_sF_ndof = track1.nDoF();
+ vars->m_sF_chi2 = track1.chi2();
+
+ vars->m_dF_ndof = track2.nDoF();
+ vars->m_dF_chi2 = track2.chi2();
+
}
//==================================================================================================
@@ -388,16 +388,16 @@ int CompareTracks::MatchesMC(const LHCb::Track &track, const LHCb::LinksByKey& l
}
auto mcpart = trackLinks.begin()->to();
if(!mcpart)return 0;
-
- //check quality of matching
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+
+ //check quality of matching
+ const MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
if ( 0 == trackInfo.fullInfo( mcpart ) ) return 2;
bool isLong = trackInfo.hasVeloAndT( mcpart );
isLong = isLong && ( abs( mcpart->particleID().pid() ) != 11 ); // and not electron
if(!isLong) return 2;
bool eta25 = (mcpart->momentum().Eta() > 1.8 && mcpart->momentum().Eta() < 5.3);
if(!eta25) return 2;
-
+
if(std::fabs(track.pseudoRapidity() - mcpart->momentum().Eta()) >0.05) return 2;
return 1;
}
diff --git a/Tr/ParameterizedKalman/src/ParameterizedKalmanFit_Checker_include.icpp b/Tr/ParameterizedKalman/src/ParameterizedKalmanFit_Checker_include.icpp
index d98c680066abc4f594833562ba0e63ffa1c10050..34075126445aa8fdb391dee9b7a6fc525e79e660 100644
--- a/Tr/ParameterizedKalman/src/ParameterizedKalmanFit_Checker_include.icpp
+++ b/Tr/ParameterizedKalman/src/ParameterizedKalmanFit_Checker_include.icpp
@@ -41,7 +41,7 @@ int ParameterizedKalmanFit_Checker::MatchesMC(const trackInfo &tI) const {
if(!mcpart)return 0;
//check quality of matching
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ const MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
if ( 0 == trackInfo.fullInfo( mcpart ) ) return 2;
bool isLong = trackInfo.hasVeloAndT( mcpart );
isLong = isLong && ( abs( mcpart->particleID().pid() ) != 11 ); // and not electron
diff --git a/Tr/PatChecker/src/CheatedPrimaryVertices.cpp b/Tr/PatChecker/src/CheatedPrimaryVertices.cpp
index 9ace3a3fee774cdf7179992be04a89af77077566..fd495881be215ee11d7c8b84f6fae3b61d471265 100755
--- a/Tr/PatChecker/src/CheatedPrimaryVertices.cpp
+++ b/Tr/PatChecker/src/CheatedPrimaryVertices.cpp
@@ -57,7 +57,6 @@ StatusCode CheatedPrimaryVertices::execute() {
if(msgLevel(MSG::DEBUG)) debug() << "Execute" << endmsg;
- MCTrackInfo trInfo( evtSvc(), msgSvc() );
LinkedTo<LHCb::MCParticle> trackMClink( evtSvc(), msgSvc(),
LHCb::TrackLocation::Default);
m_inputTracks = get<LHCb::Tracks>(m_inputTracksName);
diff --git a/Tr/PatChecker/src/PatChecker.cpp b/Tr/PatChecker/src/PatChecker.cpp
index a4ec2d4c2e6c2490c2a8949540fcbcf826735b91..ca78b7fc286ed9a0ecc69329242da4cc27576a88 100755
--- a/Tr/PatChecker/src/PatChecker.cpp
+++ b/Tr/PatChecker/src/PatChecker.cpp
@@ -223,7 +223,7 @@ StatusCode PatChecker::execute() {
LHCb::MCParticles* mcParts =
get<LHCb::MCParticles>( LHCb::MCParticleLocation::Default );
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
//== Build a table (vector of vectors) of ids per MCParticle, indexed by MCParticle key.
AllLinks<LHCb::MCParticle,ContainedObject> allIds( evtSvc(), msgSvc(), "Pat/LHCbID" );
diff --git a/Tr/PatChecker/src/PatKShortChecker.cpp b/Tr/PatChecker/src/PatKShortChecker.cpp
index 52224f5de9d03fe1db6aad51a7016c9796acd286..c61c1e2e65a8516ab502268093a2f8f31a684817 100755
--- a/Tr/PatChecker/src/PatKShortChecker.cpp
+++ b/Tr/PatChecker/src/PatKShortChecker.cpp
@@ -79,7 +79,7 @@ StatusCode PatKShortChecker::execute() {
TrAsct::DirectType* seedTable = m_seedToMCP->direct();
TrAsct::DirectType* downTable = m_downToMCP->direct();
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LHCb::MCParticles::const_iterator itP;
LHCb::Tracks::const_iterator itT;
diff --git a/Tr/PatChecker/src/PrimaryVertexChecker.cpp b/Tr/PatChecker/src/PrimaryVertexChecker.cpp
index 5c88c527c0e3be4fb5544dc4479414d8165d5b83..99699e56d01dfa1c5794f582650e99af31c18004 100644
--- a/Tr/PatChecker/src/PrimaryVertexChecker.cpp
+++ b/Tr/PatChecker/src/PrimaryVertexChecker.cpp
@@ -785,7 +785,7 @@ int PrimaryVertexChecker::count_velo_tracks(LHCb::RecVertex* RecVtx){
void PrimaryVertexChecker::count_reconstructible_mc_particles(std::vector<MCPVInfo>& mcpvvec) {
- MCTrackInfo trInfo(eventSvc(), msgSvc());
+ const MCTrackInfo trInfo = make_MCTrackInfo(eventSvc(), msgSvc());
std::vector<MCPVInfo>::iterator itinfomc;
for(itinfomc = mcpvvec.begin(); mcpvvec.end() != itinfomc; itinfomc++) {
LHCb::MCVertex* avtx = itinfomc->pMCPV;
diff --git a/Tr/PatFitParams/src/MatchFitParams.cpp b/Tr/PatFitParams/src/MatchFitParams.cpp
index 6067e461d3f796488d914acc3eeb85c12824ccef..0010defe39539b10f782cd6afccc04abbdd099e2 100644
--- a/Tr/PatFitParams/src/MatchFitParams.cpp
+++ b/Tr/PatFitParams/src/MatchFitParams.cpp
@@ -1,5 +1,4 @@
-// $Id: KsFitParams.cpp,v 1.5 2009/08/19 14:24:18 ocallot Exp $
-// Include files
+// Include files
// from Gaudi
#include "GaudiKernel/SystemOfUnits.h"
@@ -40,7 +39,7 @@ MatchFitParams::MatchFitParams( const std::string& name,
declareProperty( "zMatchY" , m_zMatchY = 8420*Gaudi::Units::mm );
declareProperty( "RequireTTHits" , m_requireTTHits = false );
declareProperty( "MagnetScaleFactor", m_magnetScaleFactor = -1 );
-
+
m_nEvent = 0;
m_nTrack = 0;
}
@@ -62,17 +61,17 @@ StatusCode MatchFitParams::initialize() {
if( m_zMagParams.empty() ) info() << "no starting values for magnet parameters provided. Will not calculate anything" << endmsg;
if( m_momParams.empty() ) info() << "no starting values for momentum parameters provided. Will not calculate anything" << endmsg;
if( m_bendYParams.empty() ) info() << "no starting values for y bending parameters provided. Will not calculate anything" << endmsg;
-
-
+
+
m_zMagPar.init( "zMagnet" , m_zMagParams );
m_momPar.init ( "momentum" , m_momParams );
m_bendParamY.init( "bendParamY", m_bendYParams);
-
+
info() << "Magnet scale factor is: " << m_magnetScaleFactor << endmsg;
-
-
+
+
m_fitTool = tool<FitTool>( "FitTool" );
-
+
return StatusCode::SUCCESS;
}
//=============================================================================
@@ -85,7 +84,7 @@ StatusCode MatchFitParams::execute() {
m_nEvent++;
debug() << "Processing event: " << m_nEvent << endmsg;
-
+
LHCb::MCParticles* partCtnr = get<LHCb::MCParticles>( LHCb::MCParticleLocation::Default );
@@ -103,14 +102,14 @@ StatusCode MatchFitParams::execute() {
LHCb::MCParticles::const_iterator pItr;
const LHCb::MCParticle* part;
SmartRefVector<LHCb::MCVertex> vDecay;
-
+
// A container for used hits
std::vector<Gaudi::XYZPoint> trHits;
std::vector<Gaudi::XYZPoint> TTHits;
std::vector<Gaudi::XYZPoint> vHits;
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
-
+ const MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
+
for ( pItr = partCtnr->begin(); partCtnr->end() != pItr; pItr++ ) {
part = *pItr;
@@ -152,16 +151,16 @@ StatusCode MatchFitParams::execute() {
}
}
if ( hasInteractionVertex ) continue;
-
-
-
+
+
+
debug() << "--- Found particle key " << part->key() << endmsg;
-
+
TTHits.clear();
trHits.clear();
vHits.clear();
-
+
// Get the Velo hits
for ( LHCb::MCHits::const_iterator vHitIt = veloHits->begin() ;
veloHits->end() != vHitIt ; vHitIt++ ) {
@@ -169,25 +168,25 @@ StatusCode MatchFitParams::execute() {
vHits.push_back( (*vHitIt)->midPoint() );
}
}
-
+
// Get the IT hits
- for ( LHCb::MCHits::const_iterator iHitIt = itHits->begin() ;
+ for ( LHCb::MCHits::const_iterator iHitIt = itHits->begin() ;
itHits->end() != iHitIt ; iHitIt++ ) {
if ( (*iHitIt)->mcParticle() == part ) {
trHits.push_back( (*iHitIt)->midPoint() );
}
}
-
+
// Get the TT hits
- for ( LHCb::MCHits::const_iterator iHittt = ttHits->begin() ;
+ for ( LHCb::MCHits::const_iterator iHittt = ttHits->begin() ;
ttHits->end() != iHittt ; iHittt++ ) {
if ( (*iHittt)->mcParticle() == part ) {
TTHits.push_back( (*iHittt)->midPoint() );
}
}
-
+
// Get the OT hits
- for ( LHCb::MCHits::const_iterator oHitIt = otHits->begin() ;
+ for ( LHCb::MCHits::const_iterator oHitIt = otHits->begin() ;
otHits->end() != oHitIt ; oHitIt++ ) {
if ( (*oHitIt)->mcParticle() == part ) {
trHits.push_back( (*oHitIt)->midPoint() );
@@ -197,38 +196,38 @@ StatusCode MatchFitParams::execute() {
if( m_requireTTHits){
if( 3 > TTHits.size() ) continue;
}
-
+
debug() << " particle momentum = " << part->momentum().R() / Gaudi::Units::GeV << " GeV"
<< endmsg;
-
+
//== Fill ntuple
double pz = part->momentum().Z();
double plXT = part->momentum().X() / pz;
double plYT = part->momentum().Y() / pz;
-
+
Tuple tTrack = nTuple( m_tupleName, m_tupleName );
-
+
tTrack->column("pz",pz);
tTrack->column("plXT", plXT);
tTrack->column("plYT", plYT);
-
-
+
+
m_nTrack++;
-
+
//== Fit the TT area
- // -- This is not needed at the moment for the matching, but it was kept it
+ // -- This is not needed at the moment for the matching, but it was kept it
// -- as it might be used at some point in the future
double axt, bxt, ayt, byt, dz;
double axt2,bxt2,cxt2;
-
+
debug() << " TT: ";
m_fitTool->fitLine( TTHits, 0, m_zTT1, axt, bxt );
m_fitTool->fitLine( TTHits, 1, m_zTT1, ayt, byt );
m_fitTool->fitParabola( TTHits, 0, m_zTT1, axt2, bxt2, cxt2 );
//m_fitTool->fitLine( TTHits, 1, m_zTT1, ayt, byt, cyt );
debug() << format( " x %7.1f tx %7.4f y %7.1f ty %7.4f ",
- axt, bxt, ayt, byt )
+ axt, bxt, ayt, byt )
<< endmsg;;
tTrack->column( "axt" , axt );
tTrack->column( "bxt", bxt );
@@ -237,8 +236,8 @@ StatusCode MatchFitParams::execute() {
tTrack->column( "axt2", axt2 );
tTrack->column( "bxt2", bxt2 );
tTrack->column( "cxt2", cxt2 );
-
-
+
+
std::vector<Gaudi::XYZPoint>::const_iterator itP;
if ( msgLevel( MSG::DEBUG ) ) {
for ( itP = TTHits.begin(); TTHits.end() > itP; itP++ ) {
@@ -254,7 +253,7 @@ StatusCode MatchFitParams::execute() {
) << endmsg;
}
}
-
+
// -- Fit the T-stations
// -- x(z) = ax + bx*z + cx*z*z + dx*z*z*z
// -- y(z) = ay + by*z
@@ -267,11 +266,11 @@ StatusCode MatchFitParams::execute() {
tTrack->column( "dx", dx );
tTrack->column( "ay", ay );
tTrack->column( "by", by );
-
+
if ( msgLevel( MSG::DEBUG ) ) {
- debug() << format( " tr: x%7.1f %7.4f %7.3f %7.3f y%7.1f %7.4f",
+ debug() << format( " tr: x%7.1f %7.4f %7.3f %7.3f y%7.1f %7.4f",
ax, bx, 1.e6*cx, 1.e9*dx, ay, by ) << endmsg;
-
+
for ( itP = trHits.begin(); trHits.end() > itP; itP++ ) {
dz = (*itP).z()-m_zRef;
debug() << format( " : %7.1f %7.1f %7.1f dx %7.3f dy %7.3f",
@@ -281,7 +280,7 @@ StatusCode MatchFitParams::execute() {
) << endmsg;
}
}
-
+
// -- Fit the velo area
// -- x(z) = axv + bxv*z
// -- y(z) = ayv + byv*z
@@ -292,11 +291,11 @@ StatusCode MatchFitParams::execute() {
tTrack->column( "bxv", bxv );
tTrack->column( "ayv", ayv );
tTrack->column( "byv", byv );
-
+
if ( msgLevel( MSG::DEBUG ) ) {
- debug() << format( " velo: x%7.1f %7.4f y%7.1f %7.4f",
+ debug() << format( " velo: x%7.1f %7.4f y%7.1f %7.4f",
axv, bxv, ayv, byv ) << endmsg;
-
+
for ( itP = vHits.begin(); vHits.end() > itP; itP++ ) {
dz = (*itP).z()-m_zRef;
debug() << format( " : %7.1f %7.1f %7.1f dx %7.3f dy %7.3f",
@@ -306,14 +305,14 @@ StatusCode MatchFitParams::execute() {
) << endmsg;
}
}
-
+
// -- This is for finding the zMagnet position when using straight lines from the Velo and the T-stations
// -- Only really makes sense in x, as for y the different y resolutions of Velo and T-stations matter a lot
double zMagnet = (axv-bxv*m_zVelo - (ax-bx*m_zRef) ) / (bx-bxv);
double zMagnety = (ayt-byt*m_zTT1 - (ay-by*m_zRef) ) / (by-byt);
double dSlope = fabs(bx - bxv);
double dSlopeY = fabs(by - byv);
-
+
m_zMagPar.setFun( 0, 1. );
m_zMagPar.setFun( 1, dSlope*dSlope );
//m_zMagPar.setFun( 2, ax*ax );
@@ -323,7 +322,7 @@ StatusCode MatchFitParams::execute() {
m_zMagPar.setFun( 4, fabs(ax) );
double zEst = m_zMagPar.sum();
-
+
m_zMagPar.addEvent( zMagnet-zEst );
tTrack->column( "zEst", zEst );
@@ -332,41 +331,41 @@ StatusCode MatchFitParams::execute() {
// -- This is the parameter that defines the bending in y
double bendParamY = (ay - ayv) + ( m_zMatchY - m_zRef )*by - (m_zMatchY - m_zVelo)*byv;
-
+
m_bendParamY.setFun( 0, dSlope*dSlope*byv );
m_bendParamY.setFun( 1, dSlopeY*dSlopeY*byv );
-
+
double bendParamYEst = m_bendParamY.sum();
m_bendParamY.addEvent( bendParamY - bendParamYEst );
// ------------------------------------------------------
-
+
tTrack->column("dSlope",dSlope);
tTrack->column("dSlope2",dSlope*dSlope);
-
+
// -- Need to write something to calculate the momentum Params
const double charge = part->particleID().threeCharge()/3;
const double qOverP = charge/momentum;
-
- // -- The magnet scale factor is hardcoded, as then one does not need to run the
- // -- field service
+
+ // -- The magnet scale factor is hardcoded, as then one does not need to run the
+ // -- field service
const double proj = sqrt( ( 1. + bxv*bxv + byv*byv ) / ( 1. + bxv*bxv ) );
const double coef = (bxv - bx)/(proj*m_magnetScaleFactor*Gaudi::Units::GeV*qOverP);
-
+
m_momPar.setFun(0, 1.);
m_momPar.setFun(1, bx*bx );
m_momPar.setFun(2, bx*bx*bx*bx );
m_momPar.setFun(3, bx*bxv );
m_momPar.setFun(4, byv*byv );
m_momPar.setFun(5, byv*byv*byv*byv );
-
+
double coefEst = m_momPar.sum();
m_momPar.addEvent( coef - coefEst );
-
+
tTrack->write();
-
+
}
-
+
return StatusCode::SUCCESS;
}
//=============================================================================
diff --git a/Tr/PatFitParams/src/PatLongLivedParams.cpp b/Tr/PatFitParams/src/PatLongLivedParams.cpp
index 60abbe5c096cb67a3da7db42febefffb6ace9d01..db606530869e006998b1a2b2adc7b2939aabda4c 100644
--- a/Tr/PatFitParams/src/PatLongLivedParams.cpp
+++ b/Tr/PatFitParams/src/PatLongLivedParams.cpp
@@ -1,4 +1,4 @@
-// Include files
+// Include files
// from Gaudi
#include "GaudiKernel/SystemOfUnits.h"
@@ -61,15 +61,15 @@ StatusCode PatLongLivedParams::initialize() {
if ( sc.isFailure() ) return sc; // error printed already by GaudiAlgorithm
debug() << "==> Initialize" << endmsg;
-
+
m_zMagPar.init( "zMagnet" , m_zMagParams );
m_momPar.init ( "momentum" , m_momParams );
m_curvature.init ( "curvature" , m_curvatureParams );
m_yPar.init ( "yPar" , m_yParams );
m_yPar2.init ( "yPar2" , m_yParams2 );
-
+
m_fitTool = tool<FitTool>( "FitTool" );
-
+
return StatusCode::SUCCESS;
}
@@ -83,14 +83,14 @@ StatusCode PatLongLivedParams::execute() {
// -- As this algorithm is the only one running, even the usual counter is not included.
if( m_nEvent % 100 == 0 ) info() << "Event " << m_nEvent << endmsg;
-
+
if( m_resolution ){
resolution();
return StatusCode::SUCCESS;
}
-
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LHCb::MCParticles* partCtnr = get<LHCb::MCParticles>( LHCb::MCParticleLocation::Default );
@@ -104,24 +104,24 @@ StatusCode PatLongLivedParams::execute() {
LHCb::MCParticles::const_iterator pItr;
const LHCb::MCParticle* part;
SmartRefVector<LHCb::MCVertex> vDecay;
-
+
const LHCb::MCParticle* kShort = 0;
const LHCb::MCVertex* kDecay = 0;
-
+
// A container for used hits
std::vector<Gaudi::XYZPoint> trHits;
std::vector<Gaudi::XYZPoint> TTHits;
-
+
for ( pItr = partCtnr->begin(); partCtnr->end() != pItr; pItr++ ) {
part = *pItr;
if( !trackInfo.hasT( part ) || !trackInfo.hasTT( part ) ) continue;
// -- optimise for downstream tracks.
if( trackInfo.hasVelo( part ) ) continue;
-
-
+
+
// Find the Pi from the kShort
- if ( 211 == std::abs(part->particleID().pid()) ) {
+ if ( 211 == std::abs(part->particleID().pid()) ) {
kDecay = part->originVertex();
if ( 0 == kDecay ) continue;
kShort = kDecay->mother();
@@ -138,14 +138,14 @@ StatusCode PatLongLivedParams::execute() {
if ( (*itV)->position().z() < 9500. ) hasInteractionVertex = true;
}
if ( hasInteractionVertex ) continue;
-
+
debug() << "--- Found pi key " << part->key() << endmsg;
-
+
TTHits.clear();
trHits.clear();
-
+
// Get the IT hits
- for ( LHCb::MCHits::const_iterator iHitIt = itHits->begin() ;
+ for ( LHCb::MCHits::const_iterator iHitIt = itHits->begin() ;
itHits->end() != iHitIt ; iHitIt++ ) {
if ( (*iHitIt)->mcParticle() == part ) {
trHits.push_back( (*iHitIt)->midPoint() );
@@ -153,7 +153,7 @@ StatusCode PatLongLivedParams::execute() {
}
// Get the TT hits
- for ( LHCb::MCHits::const_iterator iHittt = ttHits->begin() ;
+ for ( LHCb::MCHits::const_iterator iHittt = ttHits->begin() ;
ttHits->end() != iHittt ; iHittt++ ) {
if ( (*iHittt)->mcParticle() == part ) {
TTHits.push_back( (*iHittt)->midPoint() );
@@ -161,7 +161,7 @@ StatusCode PatLongLivedParams::execute() {
}
// Get the OT hits
- for ( LHCb::MCHits::const_iterator oHitIt = otHits->begin() ;
+ for ( LHCb::MCHits::const_iterator oHitIt = otHits->begin() ;
otHits->end() != oHitIt ; oHitIt++ ) {
if ( (*oHitIt)->mcParticle() == part ) {
trHits.push_back( (*oHitIt)->midPoint() );
@@ -169,16 +169,16 @@ StatusCode PatLongLivedParams::execute() {
}
if ( 3 > TTHits.size() || 12 > trHits.size() ) continue;
- debug() << "=== Found a good K0S Decay : "
- << kShort->key() << " decay at "
- << format( "%7.1f %7.1f %7.1f",
+ debug() << "=== Found a good K0S Decay : "
+ << kShort->key() << " decay at "
+ << format( "%7.1f %7.1f %7.1f",
kDecay->position().x(),
kDecay->position().y(),
kDecay->position().z() )
<< endmsg;
debug() << " pion momentum = " << part->momentum().R() / Gaudi::Units::GeV << " GeV"
<< endmsg;
-
+
//== Fill ntuple
double pz = kShort->momentum().Z();
double kSlXT = kShort->momentum().X() / pz;
@@ -200,17 +200,17 @@ StatusCode PatLongLivedParams::execute() {
tTrack->column( "moment", part->momentum().R() );
tTrack->column( "slopeX", part->momentum().X() / pz );
tTrack->column( "slopeY", part->momentum().Y() / pz );
-
+
//== Fit the TT area
double axt, bxt, ayt, byt, dz;
double axt2,bxt2,cxt2;
-
+
debug() << " TT: ";
m_fitTool->fitLine( TTHits, 0, m_zTT1, axt, bxt );
m_fitTool->fitLine( TTHits, 1, m_zTT1, ayt, byt );
m_fitTool->fitParabola( TTHits, 0, m_zTT1, axt2, bxt2, cxt2 );
debug() << format( " x %7.1f tx %7.4f y %7.1f ty %7.4f ",
- axt, bxt, ayt, byt )
+ axt, bxt, ayt, byt )
<< endmsg;
tTrack->column( "axt" , axt );
tTrack->column( "bxt", bxt );
@@ -231,7 +231,7 @@ StatusCode PatLongLivedParams::execute() {
) << endmsg;
}
}
-
+
double ax, bx, cx, dx, ay, by;
m_fitTool->fitCubic( trHits, 0, m_zRef, ax, bx, cx, dx );
m_fitTool->fitLine( trHits, 1, m_zRef, ay, by );
@@ -241,15 +241,15 @@ StatusCode PatLongLivedParams::execute() {
tTrack->column( "dx", dx );
tTrack->column( "ay", ay );
tTrack->column( "by", by );
-
-
-
-
+
+
+
+
if ( msgLevel( MSG::DEBUG ) ) {
- debug() << format( " tr: x%7.1f %7.4f %7.3f %7.3f y%7.1f %7.4f",
+ debug() << format( " tr: x%7.1f %7.4f %7.3f %7.3f y%7.1f %7.4f",
ax, bx, 1.e6*cx, 1.e9*dx, ay, by ) << endmsg;
-
+
for ( itP = trHits.begin(); trHits.end() > itP; itP++ ) {
dz = (*itP).z()-m_zRef;
debug() << format( " : %7.1f %7.1f %7.1f dx %7.3f dy %7.3f",
@@ -259,14 +259,14 @@ StatusCode PatLongLivedParams::execute() {
) << endmsg;
}
}
-
+
double zMagnet = (axt-bxt*m_zTT1 - (ax-bx*m_zRef) ) / (bx-bxt);
tTrack->column( "zMagnet", zMagnet );
double zMagnety = (ayt-byt*m_zTT1 - (ay-by*m_zRef) ) / (by-byt);
tTrack->column( "zMagnety", zMagnety );
const double yMagnet = ayt + (zMagnet-m_zTT1)*byt;
-
+
double dSlope = fabs(bx - bxt);
//double dSlope2 = (bx - bxt);
tTrack->column("dSlope",dSlope);
@@ -281,41 +281,41 @@ StatusCode PatLongLivedParams::execute() {
double dSlope2 = (bx - bxt);
for ( itP = TTHits.begin(); TTHits.end() > itP; itP++ ) {
Tuple tTrack2 = nTuple( "Track2", "Track2" );
-
+
dz = (*itP).z()-m_zTT1;
double dx2 = (*itP).x()-(axt+bxt*dz);
tTrack2->column( "dx2", dx2 );
tTrack2->column( "dz2", dz );
tTrack2->column( "dSlope2", dSlope2);
tTrack2->column( "moment2", part->momentum().R() );
-
+
tTrack2->write();
}
-
+
//m_curvature.setFun(0,0.0);
m_curvature.setFun(0,dSlope2);
-
+
double cEst = m_curvature.sum();
m_curvature.addEvent(cxt2-cEst);
tTrack->column("curvature",cxt2);
tTrack->column("dSlope2",dSlope2);
-
+
double zEst = m_zMagPar.sum();
tTrack->column( "zEst", zEst );
m_zMagPar.addEvent( zMagnet-zEst );
dSlope = std::abs(dSlope);
-
+
double dp = dSlope * part->momentum().R() - m_momPar.sum();
// -- this is a little obscure, the numbers ( '5' and '2000' ) are determined below.
// -- the numbers used here might be a bit outdated
double bytPred = by + 5. * by * fabs(by) * (bx-bxt) * (bx-bxt);
double aytPred = ay + (zMagnet-m_zRef)*by + (m_zTT1-zMagnet)*bytPred;
aytPred -= 2000*by*(bx-bxt)*(bx-bxt);
-
+
tTrack->column( "bytPred", bytPred );
tTrack->column( "aytPred", aytPred );
-
+
tTrack->column( "xMagnet", axt + ( zMagnet - m_zTT1 ) * bxt );
tTrack->column( "yMagnet", ayt + ( zMagnet - m_zTT1 ) * byt );
@@ -323,16 +323,16 @@ StatusCode PatLongLivedParams::execute() {
// -- the precise parametrisations are a bit 'black magic'
m_yPar.setFun(0, by*std::abs(by)*dSlope*dSlope );
m_yPar.addEvent( byt - by - m_yPar.sum() );
-
- // -- This is how by is determined for OT tracks
+
+ // -- This is how by is determined for OT tracks
// -- see LHCb-PUB-2017-001
//const double byExp = ay / ( m_zRef + ( m_yParams[0] * fabs(by) * zMagnet + m_yParams[2] )* dSlope2 );
-
+
const double yMagnetExp = ay + (zMagnet - m_zRef) * by;
-
+
m_yPar2.setFun(0, by*dSlope*dSlope );
m_yPar2.addEvent( yMagnet - yMagnetExp - m_yPar2.sum() );
-
+
// -- and the momentum parametrisations
m_momPar.setFun( 0, 1. );
m_momPar.setFun( 1, bx * bx );
@@ -340,14 +340,14 @@ StatusCode PatLongLivedParams::execute() {
//m_momPar.setFun( 3, bx*bx*bx*bx );
m_momPar.addEvent( dp );
tTrack->column( "PredMom", m_momPar.sum() / dSlope );
-
+
double xMag = ax + ( zEst - m_zRef ) * bx;
tTrack->column( "xPred", xMag * m_zTT1 / zEst );
tTrack->write();
-
+
double slopeX = xMag / zMagnet;
tTrack->column( "slopeX", slopeX );
-
+
debug() << format( " zMag%7.1f xMag%7.1f yMag%7.1f back%7.1f",
zMagnet,
axt + ( zMagnet - m_zTT1 ) * bxt,
@@ -362,9 +362,9 @@ StatusCode PatLongLivedParams::execute() {
// Determine resolution
//=============================================================================
void PatLongLivedParams::resolution(){
-
+
LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
// Get the IT and TT hits
LHCb::MCHits* itHits = get<LHCb::MCHits>( LHCb::MCHitLocation::IT );
@@ -374,8 +374,8 @@ void PatLongLivedParams::resolution(){
LHCb::MCHits* otHits = get<LHCb::MCHits>( LHCb::MCHitLocation::OT );
LHCb::Tracks* inTracks = getIfExists<LHCb::Tracks>( m_inputLocationSeed );
-
- if( inTracks == nullptr ){
+
+ if( inTracks == nullptr ){
error() << "Could not find tracks in " << m_inputLocationSeed << endmsg;
return;
}
@@ -385,30 +385,30 @@ void PatLongLivedParams::resolution(){
std::vector<Gaudi::XYZPoint> TTHits;
for( LHCb::Track* track : *inTracks){
-
+
const LHCb::MCParticle* mcSeedPart = mySeedLink.first( track->key() );
if( mcSeedPart == nullptr ) continue;
if( std::abs(mcSeedPart->particleID().pid()) == 11 ) continue;
-
+
if( !trackInfo.hasT( mcSeedPart) || !trackInfo.hasTT( mcSeedPart ) ) continue;
// -- optimise for downstream tracks
if( trackInfo.hasVelo( mcSeedPart ) ) continue;
-
+
const LHCb::MCVertex* kDecay = mcSeedPart->originVertex();
if ( nullptr == kDecay ) continue;
const LHCb::MCParticle* kShort = kDecay->mother();
if ( nullptr == kShort ) continue;
-
+
// -- Ks and Lambda
if ( 310 != std::abs(kShort->particleID().pid()) && 3122 != std::abs(kShort->particleID().pid()) ) continue;
-
+
const LHCb::MCVertex* origin = kShort->originVertex();
if ( nullptr == origin ) continue;
if ( 100. < origin->position().R() ){
debug()<< "Too far from beampipe" << endmsg;
continue; // particles from close the beam line
}
-
+
bool hasInteractionVertex = false;
SmartRefVector<LHCb::MCVertex> endV = mcSeedPart->endVertices();
for ( SmartRefVector<LHCb::MCVertex>::const_iterator itV = endV.begin() ;
@@ -419,39 +419,39 @@ void PatLongLivedParams::resolution(){
debug()<< "Interaction vertex found. skipping" << endmsg;
continue;
}
-
+
TTHits.clear();
trHits.clear();
-
+
// Get the IT hits
- for ( LHCb::MCHits::const_iterator iHitIt = itHits->begin() ;
+ for ( LHCb::MCHits::const_iterator iHitIt = itHits->begin() ;
itHits->end() != iHitIt ; iHitIt++ ) {
if ( (*iHitIt)->mcParticle() == mcSeedPart ) {
trHits.push_back( (*iHitIt)->midPoint() );
}
}
-
+
// Get the TT hits
- for ( LHCb::MCHits::const_iterator iHittt = ttHits->begin() ;
+ for ( LHCb::MCHits::const_iterator iHittt = ttHits->begin() ;
ttHits->end() != iHittt ; iHittt++ ) {
if ( (*iHittt)->mcParticle() == mcSeedPart ) {
TTHits.push_back( (*iHittt)->midPoint() );
}
}
-
+
// Get the OT hits
- for ( LHCb::MCHits::const_iterator oHitIt = otHits->begin() ;
+ for ( LHCb::MCHits::const_iterator oHitIt = otHits->begin() ;
otHits->end() != oHitIt ; oHitIt++ ) {
if ( (*oHitIt)->mcParticle() == mcSeedPart ) {
trHits.push_back( (*oHitIt)->midPoint() );
}
}
if ( 4 > TTHits.size() || 11 > trHits.size() ) continue;
-
+
//== Fit the TT area
double axt, bxt, ayt, byt;
//double axt2,bxt2,cxt2;
-
+
debug() << " TT: ";
m_fitTool->fitLine( TTHits, 0, m_zTT1, axt, bxt );
m_fitTool->fitLine( TTHits, 1, m_zTT1, ayt, byt );
@@ -464,34 +464,34 @@ void PatLongLivedParams::resolution(){
const double zMagnet = (axt-bxt*m_zTT1 - (ax-bx*m_zRef) ) / (bx-bxt);
const double xMagnet = axt + ( zMagnet - m_zTT1 ) * bxt;
const double yMagnet = ayt + ( zMagnet - m_zTT1 ) * byt;
-
+
LHCb::State* state = &track->closestState( 10000. );
- const double zMagnetExp =m_zMagParams [0] +
- m_zMagParams[1] * state->ty() * state->ty() +
+ const double zMagnetExp =m_zMagParams [0] +
+ m_zMagParams[1] * state->ty() * state->ty() +
m_zMagParams[2] * state->tx() * state->tx() +
m_zMagParams[3] * 1/state->p() +
m_zMagParams[4] * std::abs( state->x() ) +
m_zMagParams[5] * std::abs( state->y() ) +
m_zMagParams[6] * std::abs( state->ty() );
-
+
const double dzExp = zMagnetExp - state->z();
const double xMagnetExp = state->x() + dzExp * state->tx();
const double bxtExp = xMagnetExp / zMagnetExp;
const double dSlope = std::abs( bxtExp - state->tx() );
const double dSlope2 = dSlope*dSlope;
-
+
// -- Number of IT hits
- const unsigned int nbIT = std::count_if( track->lhcbIDs().begin(), track->lhcbIDs().end(),
+ const unsigned int nbIT = std::count_if( track->lhcbIDs().begin(), track->lhcbIDs().end(),
[](const LHCb::LHCbID id){ return id.isIT();});
// -- this is how it is done in PatDownTrack for OT tracks.
// -- see LHCb-PUB-2017-001
- double byExp = state->y() / ( state->z() +
+ double byExp = state->y() / ( state->z() +
( m_yParams[0] * std::abs(state->ty()) * zMagnet + m_yParams2[0] )* dSlope2 );
-
+
if( nbIT > 4 ) byExp = state->ty();
const double bytExp = byExp*(1.0+m_yParams[0]*dSlope2*std::abs(byExp));
const double yMagnetExp = state->y() + dzExp * byExp - m_yParams2[0] * byExp * dSlope2;
@@ -523,11 +523,11 @@ void PatLongLivedParams::resolution(){
resoTuple->column( "byt", byt );
resoTuple->column( "bytExp", bytExp );
resoTuple->column( "ayt", ayt );
-
+
resoTuple->write();
-
+
}
-
+
}
//=============================================================================
@@ -545,14 +545,14 @@ StatusCode PatLongLivedParams::finalize() {
m_curvature.updateParameters( msg );
m_yPar.updateParameters( msg );
m_yPar2.updateParameters( msg );
-
+
std::cout << std::endl;
m_zMagPar.printPythonParams( name() );
m_momPar.printPythonParams( name() );
m_curvature.printPythonParams( name() );
m_yPar.printPythonParams( name() );
m_yPar2.printPythonParams( name() );
-
+
std::cout << std::endl;
m_zMagPar.printParams( "PatLongLivedTracking" );
diff --git a/Tr/PatFitParams/src/SeedFitParams.cpp b/Tr/PatFitParams/src/SeedFitParams.cpp
index 17814e3fc06974e786452643214943e9f8fb1ab0..caa312d106956978c0e3e0be54df8dad666b9519 100755
--- a/Tr/PatFitParams/src/SeedFitParams.cpp
+++ b/Tr/PatFitParams/src/SeedFitParams.cpp
@@ -1,4 +1,4 @@
-// Include files
+// Include files
// from Gaudi
#include "Event/MCTrackInfo.h"
@@ -66,14 +66,14 @@ SeedFitParams::SeedFitParams( const std::string& name,
std::vector<double> tmp = list_of(1.25e-14);
declareProperty( "yCorrectionParams", m_yCorrection = tmp);
}
-
+
m_nEvent = 0;
m_nTrack = 0;
}
//=============================================================================
// Destructor
//=============================================================================
-SeedFitParams::~SeedFitParams() {}
+SeedFitParams::~SeedFitParams() {}
//=============================================================================
// Initialization
@@ -106,7 +106,7 @@ StatusCode SeedFitParams::execute() {
const auto* partCtnr = get<LHCb::MCParticles>( LHCb::MCParticleLocation::Default );
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
// Get the MC hits, and create our own "linkers"
std::unordered_map<const LHCb::MCParticle*, std::vector<const LHCb::MCHit*> > tlinker;
@@ -181,7 +181,7 @@ StatusCode SeedFitParams::execute() {
ttHits.reserve(ttlinker[part].size());
for (const auto* hit: ttlinker[part]) ttHits.push_back(hit->midPoint());
}
-
+
++m_nTrack;
if ( msgLevel( MSG::DEBUG ) ) {
debug() << format( "Track MC %4d z0 %7.2f p%8.2f Velo%2d TT%2d T%2d",
@@ -202,7 +202,7 @@ StatusCode SeedFitParams::execute() {
tTrack->column( "slopeY", part->momentum().Y() / pz );
double ax, bx, cx, dx, ay, by;
-
+
m_fitTool->fitCubic( trHits, 0, m_zRef, ax, bx, cx, dx );
m_fitTool->fitLine ( trHits, 1, m_zRef, ay, by );
@@ -217,7 +217,7 @@ StatusCode SeedFitParams::execute() {
debug() << format( "p %7.3f, N%4d, ax%8.2f bx%8.2f cx%8.2f dp%10.4f",
momentum/1000., trHits.size(),
ax, 1.e3*bx, 1.e6*cx, dp ) << endmsg;
- }
+ }
// work out intercept point of line joining T1 and T3 with z = 0
const double dzT1 = StateParameters::ZBegT - m_zRef;
@@ -244,13 +244,13 @@ StatusCode SeedFitParams::execute() {
//m_fitTool->fitLine( ttHits, 0, m_zTT, axt, bxt );
m_fitTool->fitParabola( ttHits, 0, m_zTT, axt, bxt, cxt );
m_fitTool->fitLine( ttHits, 1, m_zTT, ayt, byt );
-
+
double dz = m_zSeed - m_zRef;
double x0 = ax + dz * ( bx + dz * ( cx + dz * dx ) );
double tx = bx + dz * (2*cx + dz * 3*dx );
zMagnet = (axt-bxt*m_zTT - (x0-tx*m_zSeed) ) / (tx-bxt);
-
+
m_zMagPar.setFun( 0, 1. );
m_zMagPar.setFun( 1, by*by );
m_zMagPar.setFun( 2, bx*bx );
@@ -282,7 +282,7 @@ StatusCode SeedFitParams::execute() {
m_yCorrectionPar.setFun(0, 1.0);
if (400. < fabs(y0))
m_yCorrectionPar.addEvent(fabs((by * by * cx * cx) / y0) - m_yCorrectionPar.sum());
-
+
tTrack->column( "ax", ax );
tTrack->column( "bx", bx );
tTrack->column( "cx", cx );
@@ -333,9 +333,9 @@ StatusCode SeedFitParams::finalize() {
m_dRatioPar.printPythonParams( name() );
m_yCorrectionPar.printPythonParams( name() );
std::cout << std::endl;
-
+
std::string toolName = "ToolSvc.PatSeedingTool";
-
+
m_initialArrowPar.printPythonParams( toolName );
m_momentumScalePar.printPythonParams( toolName );
m_zMagPar.printPythonParams( toolName );
diff --git a/Tr/TrackMCTools/src/DebugTrackingLosses.cpp b/Tr/TrackMCTools/src/DebugTrackingLosses.cpp
index a75972a37f8b9d86a8361b3c03ebb7c33b224c05..99468d4cc27e7dab5e0462071506d8372d1693d1 100644
--- a/Tr/TrackMCTools/src/DebugTrackingLosses.cpp
+++ b/Tr/TrackMCTools/src/DebugTrackingLosses.cpp
@@ -86,7 +86,7 @@ StatusCode DebugTrackingLosses::execute() {
if ( m_vP ) veloTrack = LHCb::TrackLocation::VP;
LinkedFrom<LHCb::Track,LHCb::MCParticle> veloLinker ( evtSvc(), msgSvc(), veloTrack );
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LHCb::MCParticles::const_iterator itP = partCont->begin();
for (; itP != partCont->end(); ++itP){
diff --git a/Tr/TrackMCTools/src/PatDebugTTTruthTool.cpp b/Tr/TrackMCTools/src/PatDebugTTTruthTool.cpp
index d2703ecc31494399de35b00a2aff475e6c061693..dc70651b58e44b19082785a38d0deace8de19f77 100755
--- a/Tr/TrackMCTools/src/PatDebugTTTruthTool.cpp
+++ b/Tr/TrackMCTools/src/PatDebugTTTruthTool.cpp
@@ -188,18 +188,18 @@ double PatDebugTTTruthTool::fracGoodHits( const LHCb::Track* track, const PatTTH
//=============================================================================
// Is this seed track not a ghost and downstream reconstructible
//=============================================================================
-bool PatDebugTTTruthTool::isDownRecoble( const LHCb::Track* track, const Tf::TTStationHitManager<PatTTHit>::HitRange& hits,
+bool PatDebugTTTruthTool::isDownRecoble( const LHCb::Track* track, const Tf::TTStationHitManager<PatTTHit>::HitRange& hits,
const unsigned int minHits ){
-
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
-
+
const LHCb::MCParticle* mcSeedPart = mySeedLink.first( track->key() );
if( mcSeedPart == nullptr ) return false;
if( std::abs( mcSeedPart->particleID().pid() ) == 11 ) return false;
-
+
unsigned int nTrue = 0;
for( PatTTHit* hit : hits ){
if( isTrueHit( track, hit ))nTrue++;
@@ -208,8 +208,8 @@ bool PatDebugTTTruthTool::isDownRecoble( const LHCb::Track* track, const Tf::TTS
if( trackInfo.hasT( mcSeedPart) && trackInfo.hasTT( mcSeedPart ) && nTrue >= minHits ) return true;
return false;
-
-
+
+
}
//=============================================================================
// Is this a 'good' track (according to the TrackAssociator) definition?
@@ -256,7 +256,7 @@ std::map<std::string, bool> PatDebugTTTruthTool::getMCTrackInfo(const LHCb::Trac
std::map<std::string,bool> parameterMap;
LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
LinkedTo<LHCb::MCParticle, LHCb::STCluster> myClusterLink ( evtSvc(), msgSvc(), LHCb::STClusterLocation::TTClusters );
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
const LHCb::MCParticle* mcSeedPart = mySeedLink.first( seed->key() );
if(mcSeedPart)
@@ -278,16 +278,15 @@ std::map<std::string, bool> PatDebugTTTruthTool::getMCTrackInfo(const LHCb::Trac
bool PatDebugTTTruthTool::hasMCParticle(const LHCb::Track* seed) const
{
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
const LHCb::MCParticle* mcSeedPart = mySeedLink.first( seed->key() );
- return mcSeedPart != nullptr;
+ return mcSeedPart != nullptr;
}
-
+
bool PatDebugTTTruthTool::isDownstreamReconstructible(const LHCb::Track* seed) const
{
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
const LHCb::MCParticle* mcSeedPart = mySeedLink.first( seed->key() );
return mcSeedPart && trackInfo.hasT( mcSeedPart ) && trackInfo.hasTT( mcSeedPart );
@@ -295,7 +294,6 @@ bool PatDebugTTTruthTool::isDownstreamReconstructible(const LHCb::Track* seed) c
bool PatDebugTTTruthTool::hasMCParticleNotElectron(const LHCb::Track* seed) const
{
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
const LHCb::MCParticle* mcSeedPart = mySeedLink.first( seed->key() );
return mcSeedPart && std::abs(mcSeedPart->particleID().pid()) != 11;
@@ -303,11 +301,11 @@ bool PatDebugTTTruthTool::hasMCParticleNotElectron(const LHCb::Track* seed) cons
bool PatDebugTTTruthTool::isDownstreamReconstructibleNotElectron(const LHCb::Track* seed) const
{
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
const LHCb::MCParticle* mcSeedPart = mySeedLink.first( seed->key() );
- return mcSeedPart && trackInfo.hasT( mcSeedPart ) &&
- trackInfo.hasTT( mcSeedPart )&&
+ return mcSeedPart && trackInfo.hasT( mcSeedPart ) &&
+ trackInfo.hasTT( mcSeedPart )&&
std::abs(mcSeedPart->particleID().pid()) != 11 ;
}
@@ -315,12 +313,12 @@ bool PatDebugTTTruthTool::isDownstreamReconstructibleNotElectron(const LHCb::Tra
bool PatDebugTTTruthTool::isTrueSeed(const LHCb::Track* seed) const
{
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
-
+
const LHCb::MCParticle* mcSeedPart = mySeedLink.first( seed->key() );
- return mcSeedPart && trackInfo.hasT( mcSeedPart ) &&
- trackInfo.hasTT( mcSeedPart ) &&
+ return mcSeedPart && trackInfo.hasT( mcSeedPart ) &&
+ trackInfo.hasTT( mcSeedPart ) &&
!trackInfo.hasVelo( mcSeedPart) &&
std::abs(mcSeedPart->particleID().pid()) != 11 ;
@@ -362,11 +360,11 @@ void PatDebugTTTruthTool::seedTuple(const LHCb::Track* trackSeed, double seedCla
tuple->column("seed_p", trackSeed->p());
tuple->column("seed_pt", trackSeed->pt());
tuple->column("seed_nLHCbIDs", trackSeed->nLHCbIDs());
- const unsigned int nbIT = std::count_if( trackSeed->lhcbIDs().begin(), trackSeed->lhcbIDs().end(),
+ const unsigned int nbIT = std::count_if( trackSeed->lhcbIDs().begin(), trackSeed->lhcbIDs().end(),
[](const LHCb::LHCbID id){ return id.isIT();});
tuple->column("seed_nbIT", nbIT);
// number of layers
- std::bitset<12> layers;
+ std::bitset<12> layers;
auto st2layer = [](LHCb::STChannelID id) { return 4*(id.station()-1)+id.layer()-1; };
for( auto id : trackSeed->lhcbIDs() ){
layers.set( id.isOT() ? id.otID().sequentialUniqueLayer() : st2layer(id.stID())) ;
@@ -388,15 +386,15 @@ void PatDebugTTTruthTool::seedTuple(const LHCb::Track* trackSeed, double seedCla
//=========================================================================
void PatDebugTTTruthTool::finalTuple( const bool goodTrack, const LHCb::Track* track, const std::vector<double>& vals ){
- MCTrackInfo trackInfo( evtSvc(), msgSvc() );
+ MCTrackInfo trackInfo = make_MCTrackInfo( evtSvc(), msgSvc() );
LinkedTo<LHCb::MCParticle, LHCb::Track> mySeedLink ( evtSvc(), msgSvc(),LHCb::TrackLocation::Seed);
const LHCb::MCParticle* mcSeedPart = mySeedLink.first( track->key() );
// -- Don't train on long tracks, but don't classify as ghosts either
if( mcSeedPart && goodTrack && trackInfo.hasT( mcSeedPart) && trackInfo.hasTT( mcSeedPart ) && trackInfo.hasVelo( mcSeedPart) ) return;
-
+
if(mcSeedPart){
if( std::abs(mcSeedPart->particleID().pid()) == 11 ) return;
-
+
bool goodDau = false;
if ( 0 != mcSeedPart->originVertex() ) {
const LHCb::MCParticle* mother = mcSeedPart->originVertex()->mother();
@@ -420,14 +418,14 @@ void PatDebugTTTruthTool::finalTuple( const bool goodTrack, const LHCb::Track* t
}
if( !goodDau ) return;
}
-
+
bool gT = goodTrack && (mcSeedPart != nullptr);
/*
- std::vector<double> vals ={
- vdt::fast_log(track.chi2()),
+ std::vector<double> vals ={
+ vdt::fast_log(track.chi2()),
vdt::fast_log(track.p()),
- vdt::fast_log(track.pt()),
+ vdt::fast_log(track.pt()),
std::abs(deltaP),
vdt::fast_log(std::abs(track.dxMagnet())),
vdt::fast_log(std::abs(distance)),
@@ -449,7 +447,7 @@ void PatDebugTTTruthTool::finalTuple( const bool goodTrack, const LHCb::Track* t
tuple->column("firedLayers", vals[8] );
tuple->column("goodTrack", gT );
tuple->write();
-
+
}