Skip to content
Snippets Groups Projects
Commit ade07d8c authored by Marco Clemencic's avatar Marco Clemencic
Browse files

v11r2p1

parent 47e17e8d
No related branches found
No related tags found
No related merge requests found
Showing
with 1078 additions and 493 deletions
......@@ -3,7 +3,7 @@
# Maintainer : Jaap Panman
#============================================================================
package LumiAlgs
version v5r0
version v5r0p1
#============================================================================
# Structure, i.e. directories to process.
......@@ -24,7 +24,7 @@ use GaudiPoolDb v*
use HltInterfaces v* Kernel
### for access to the config.tar file with configurations
use HltTCK v1r* TCK -no_auto_imports
use HltTCK v* TCK -no_auto_imports
include_path none
#============================================================================
......
v5r0
v5r0p1
......@@ -4,6 +4,13 @@
! Purpose : Accounting of LumiSummaries (offline code)
!-----------------------------------------------------------------------------
!========================= LumiAlgs v5r0p1 2011-03-09 =========================
! 2011-03-08 - Gerhard Raven
- replace 'use HltTCK v1r*' with 'use HltTCK v*'
! 2011-03-02 - Marco Cattaneo
- Fix trivial compilation problem on Windows (logical not instead of !)
!========================= LumiAlgs v5r0 2011-02-28 =========================
! 2011-03-02 - Jaap Panman
- update unit tests to avoid using the PropertyConfigSvc in Rec
......
......@@ -410,7 +410,7 @@ StatusCode GetLumiParameters::i_cacheFillingData() {
StatusCode GetLumiParameters::processDB() {
// only run this after initialization
if ( not m_initialized ) return StatusCode::SUCCESS;
if ( !m_initialized ) return StatusCode::SUCCESS;
Gaudi::Time xtfound = m_dds->eventTime();
......@@ -466,7 +466,7 @@ double GetLumiParameters::rateFromTCK(unsigned int tck) {
double rate = 0;
// only run this after initialization
if ( not m_initialized ) return rate;
if ( !m_initialized ) return rate;
std::string code = "Code";
std::string sub = "RATE(";
......
......@@ -3,7 +3,7 @@
# Maintainer : Marco CATTANEO
#============================================================================
package RecAlgs
version v2r1
version v2r2
branches cmt doc src
include_path none
......
v2r1
v2r2
......@@ -4,6 +4,12 @@
! Purpose : Holds generic reconstruction algorithms
!-----------------------------------------------------------------------------
!========================= RecAlgs v2r2 2011-03-09 =========================
! 2011-03-09 - Rob Lambert
- Modified RecInit.cpp
Add raw file GUID to rec header. Requires Event/RecEvent
>= v2r36
!========================= RecAlgs v2r1 2011-02-28 =========================
! 2011-02-21 - Chris Jones
- Add # of TT clusters to RecSummary object in RecSummaryAlg
......
......@@ -7,6 +7,7 @@
#include "GaudiKernel/Incident.h"
#include "GaudiKernel/IIncidentSvc.h"
#include "GaudiAlg/IGenericTool.h"
#include "GaudiKernel/IOpaqueAddress.h"
// from EventBase
#include "Event/ProcessHeader.h"
......@@ -39,6 +40,10 @@ RecInit::RecInit( const std::string& name,
ISvcLocator* pSvcLocator)
: LbAppInit ( name , pSvcLocator )
{
declareProperty( "RawEventLocation" , m_rawEventLocation = LHCb::RawEventLocation::Default,
"where to find the raw event" );
declareProperty( "AbortOnFID" , m_abortOnFID = true,
"If I can't find the raw file ID, do I abort? Default true=yes.");
}
//=============================================================================
......@@ -97,12 +102,37 @@ StatusCode RecInit::execute()
const ulonglong evtNumber = odin->eventNumber();
this->printEventRun( evtNumber, runNumber, 0, odin->eventTime() );
// Initialize the random number
std::vector<long int> seeds = getSeeds( runNumber, evtNumber );
sc = this->initRndm( seeds );
if ( sc.isFailure() ) return sc; // error printed already by initRndm
// get the file ID from the raw event
std::string event_fname="";
if( !exist<LHCb::RawEvent>(m_rawEventLocation) )
{
if(m_abortOnFID) return Error("RawBank cannot be loaded, fileID cannot be found");
Warning("RawBank cannot be loaded, fileID cannot be found",StatusCode::SUCCESS).ignore();
}
else
{
LHCb::RawEvent* event = get<LHCb::RawEvent>(m_rawEventLocation);
IOpaqueAddress* eAddr = event->registry()->address();
// obtain the fileID
if ( eAddr )
{
event_fname = eAddr->par()[0];
if ( msgLevel(MSG::DEBUG) ) debug() << "RunInfo record from Event: " << event_fname << endmsg;
}
else
{
if(m_abortOnFID) return Error("Registry cannot be loaded from Event, fileID cannot be found");
Warning("Registry cannot be loaded from Event, fileID cannot be found",StatusCode::SUCCESS).ignore();
}
}
const std::string rawID=event_fname;
// Create the Reconstruction event header
LHCb::RecHeader* header = new LHCb::RecHeader();
header->setApplicationName( this->appName() );
......@@ -112,6 +142,7 @@ StatusCode RecInit::execute()
header->setRandomSeeds( seeds );
header->setCondDBTags( this->condDBTags() );
header->setGpsTime( odin->gpsTime() );
header->setRawID( rawID );
put( header, LHCb::RecHeaderLocation::Default );
// Create a ProcStatus if it does not already exist
......
......@@ -31,5 +31,7 @@ protected:
private:
IGenericTool* m_memoryTool; ///< Pointer to (private) memory histogram tool
IIncidentSvc* m_incidentSvc; ///< Pointer to the incident service.
std::string m_rawEventLocation; ///< Location where we get the RawEvent
bool m_abortOnFID; ///< If I can't find the raw file ID, do I abort?
};
#endif // RECINIT_H
package RecSys
version v11r2
version v11r2p1
branches cmt doc tests
......@@ -29,8 +29,8 @@ use MuonTrackRec v3r4p1 Muon # Giacomo Graziani
# Non sub-detector specific components
use ChargedProtoANNPID v1r3 Rec # Chris Jones
use GlobalReco v6r35 Rec # Chris Jones
use LumiAlgs v5r0 Rec # Jaap Panman
use RecAlgs v2r1 Rec # Marco Cattaneo
use LumiAlgs v5r0p1 Rec # Jaap Panman
use RecAlgs v2r2 Rec # Marco Cattaneo
use RecConf v2r0 Rec # Rob Lambert
use RecInterfaces v1r0 Rec # Chris Jones
......@@ -80,7 +80,7 @@ use TrackFitEvent v5r7 Tr # Eduardo Rodrigues
use TrackFitter v4r2 Tr # Eduardo Rodrigues
use TrackIdealPR v2r19 Tr # Eduardo Rodrigues
use TrackMCTools v3r0p1 Tr # Eduardo Rodrigues
use TrackMonitors v1r30p1 Tr # Wouter Hulsbergen, Stephanie Menzemer
use TrackMonitors v1r31 Tr # Wouter Hulsbergen, Stephanie Menzemer
use TrackProjectors v2r29 Tr # Wouter Hulsbergen
use TrackTools v4r3 Tr # Eduardo Rodrigues
use TrackUtils v1r41 Tr # Eduardo Rodrigues
......
v11r2
v11r2p1
......@@ -5,6 +5,21 @@ Purpose: LHCb reconstruction packages.
This project groups together all the reconstruction packages needed
by one or more LHCb applications.
</PRE><H1><A NAME=v11r2p1>2011-03-10 RecSys v11r2p1</A></H1><PRE>
This version uses Gaudi v22r1 and LHCb v32r1p1
- This is a production version of Rec, required to pick up certain
patches for 2011 data.
- The L0TCK major version has increased requiring a
L0DU (Lbcom) and Rec/LumiAlgs v5r0p1 patch
- The RecHeader now also holds for each event a string of the
raw file GUID, required before we can use conditions to
perform DQ checks, This is filled by RecInit in Rec/RecAlgs v2r2
- Minor changes were also made to Tr/TrackMonitors v3r31
</PRE><H1><A NAME=v11r2>2011-02-28 RecSys v11r2</A></H1><PRE>
This version uses Gaudi v22r1 and LHCb v32r1
......
#============================================================================
# Created : 2008-05-30
# Maintainer : Wouter Hulsbergen
#============================================================================
package TrackMonitors
version v1r30p1
#============================================================================
# Structure, i.e. directories to process.
#============================================================================
branches cmt doc src
include_path none
#============================================================================
# Used packages.
#============================================================================
use GSL v* LCG_Interfaces
use GaudiAlg v*
use GaudiUtils v*
use RecEvent v* Event
use TrackEvent v* Event
use TrackFitEvent v* Tr
use TrackInterfaces v* Tr
use TrackKernel v* Tr
use STKernel v* ST
use PartProp v* Kernel
use CaloDet v* Det
use PhysEvent v* Event
#============================================================================
# Component library building rule
#============================================================================
library TrackMonitors ../src/*.cpp -import=AIDA
#============================================================================
# define component library link options
#============================================================================
apply_pattern component_library library=TrackMonitors
apply_pattern install_python_modules
#============================================================================
# Created : 2008-05-30
# Maintainer : Wouter Hulsbergen
#============================================================================
package TrackMonitors
version v1r31
#============================================================================
# Structure, i.e. directories to process.
#============================================================================
branches cmt doc src
include_path none
#============================================================================
# Used packages.
#============================================================================
use GSL v* LCG_Interfaces
use GaudiAlg v*
use GaudiUtils v*
use RecEvent v* Event
use TrackEvent v* Event
use TrackFitEvent v* Tr
use TrackInterfaces v* Tr
use TrackKernel v* Tr
use STKernel v* ST
use PartProp v* Kernel
use CaloDet v* Det
use PhysEvent v* Event
#============================================================================
# Component library building rule
#============================================================================
library TrackMonitors ../src/*.cpp -import=AIDA
#============================================================================
# define component library link options
#============================================================================
apply_pattern component_library library=TrackMonitors
apply_pattern install_python_modules
v1r30p1
v1r31
This diff is collapsed.
......@@ -63,19 +63,19 @@ StatusCode TrackTune::execute()
// output tuple
Tuple myTuple = nTuple("Candidates");
const LHCb::Tracks* tracks = get<LHCb::Tracks*>(m_trackLocation);
const LHCb::Particle::Container* particles = get<LHCb::Particle::Container>(m_particleLocation);
const LHCb::Track::Range tracks = get<LHCb::Track::Range>(m_trackLocation);
const LHCb::Particle::Range particles = get<LHCb::Particle::Range>(m_particleLocation);
std::vector<const LHCb::Particle* > tVec;
if (select(tVec,particles) == false) return StatusCode::SUCCESS ;
for (std::vector<const LHCb::Particle* >::const_iterator iterP = tVec.begin(); iterP != tVec.end(); ++iterP ){
bool rec = isFound(tracks,*iterP);
bool rec = isFound(tracks,**iterP);
myTuple << Tuples::Column("M", (*iterP)->measuredMass())
<< Tuples::Column("found",rec)
<< Tuples::Column("PT", (*iterP)->pt())
<< Tuples::Column("Candidates", particles->size());
<< Tuples::Column("Candidates", particles.size());
myTuple->write();
}
......@@ -84,28 +84,27 @@ StatusCode TrackTune::execute()
//
} // the end of the Algorihtm
const LHCb::Track* TrackTune::track(const LHCb::Particle* part) const{
const LHCb::ProtoParticle* proto = part->proto();
const LHCb::Track* TrackTune::track(const LHCb::Particle& part) const
{
const LHCb::ProtoParticle* proto = part.proto();
if (!proto || proto->charge() == 0) return 0;
return proto->track() ;
}
bool TrackTune::isFound(const LHCb::Tracks* tracks, const LHCb::Particle* part) const {
bool TrackTune::isFound(const LHCb::Track::Range& tracks, const LHCb::Particle& part) const
{
bool ok = true;
const SmartRefVector<LHCb::Particle>& daughters = part->daughters();
const SmartRefVector<LHCb::Particle>& daughters = part.daughters();
for (SmartRefVector<LHCb::Particle>::const_iterator iter = daughters.begin();
iter != daughters.end() && ok == true ; ++iter){
const LHCb::Track* aTrack = track(*iter);
const LHCb::Track* aTrack = track(**iter);
if (!aTrack) {
info() << "Failed to find track " << endmsg;
}
const double nHits = aTrack->nLHCbIDs();
bool matchedTrack = false;
for (LHCb::Tracks::const_iterator iterT = tracks->begin(); iterT != tracks->end() && matchedTrack == false; ++iterT){
for (LHCb::Track::Range::const_iterator iterT = tracks.begin(); iterT != tracks.end() && matchedTrack == false; ++iterT){
const double fracCommon = aTrack->nCommonLhcbIDs(**iterT)/double(nHits);
plot(fracCommon, "purity", "purity", 0., 2., 100);
if (fracCommon > m_minPurityCut) matchedTrack = true;
......@@ -118,13 +117,13 @@ bool TrackTune::isFound(const LHCb::Tracks* tracks, const LHCb::Particle* part)
}
bool TrackTune::select(std::vector<const LHCb::Particle* >& output, const LHCb::Particle::Container* input) const{
bool TrackTune::select(std::vector<const LHCb::Particle* >& output, const LHCb::Particle::Range& input) const
{
if (m_selectBest == true){
double bestChi2 = 9999.; const LHCb::Particle* bestPart = 0;
for (LHCb::Particle::Container::const_iterator iter = input->begin(); iter != input->end(); ++iter){
if (inMassRange(*iter) == false) continue;
LHCb::Vertex* vert = (*iter)->endVertex();
for (LHCb::Particle::Range::const_iterator iter = input.begin(); iter != input.end(); ++iter){
if (inMassRange(**iter) == false) continue;
const LHCb::Vertex* vert = (*iter)->endVertex();
if (vert->chi2PerDoF() < bestChi2){
bestChi2= vert->chi2PerDoF();
bestPart = *iter;
......@@ -133,15 +132,16 @@ bool TrackTune::select(std::vector<const LHCb::Particle* >& output, const LHCb::
if (bestPart) output.push_back(bestPart);
}
else {
for (LHCb::Particle::Container::const_iterator iter = input->begin(); iter != input->end(); ++iter){
for (LHCb::Particle::Range::const_iterator iter = input.begin(); iter != input.end(); ++iter){
output.push_back(*iter);
}
}
return (output.size() == 0 ? false : true );
}
bool TrackTune::inMassRange(const LHCb::Particle* particle) const {
const double m = particle->measuredMass();
bool TrackTune::inMassRange(const LHCb::Particle& particle) const
{
const double m = particle.measuredMass();
return (m > m_minMass && m< m_maxMass);
}
......@@ -24,13 +24,13 @@ class TrackTune: public GaudiTupleAlg{
private:
const LHCb::Track* track(const LHCb::Particle* part) const;
const LHCb::Track* track(const LHCb::Particle& part) const;
bool isFound(const LHCb::Tracks* tracks, const LHCb::Particle* part) const;
bool isFound(const LHCb::Track::Range& tracks, const LHCb::Particle& part) const;
bool select(std::vector<const LHCb::Particle* >& output, const LHCb::Particle::Container* input) const;
bool select(std::vector<const LHCb::Particle*>& output, const LHCb::Particle::Range& input) const;
bool inMassRange(const LHCb::Particle* particle) const;
bool inMassRange(const LHCb::Particle& particle) const;
std::string m_particleLocation;
std::string m_trackLocation;
......
// $Id: VertexCompare.cpp,v 1.1 2010/03/01 15:58:08 pmorawsk Exp $
// Include files
// from Gaudi
#include "GaudiKernel/AlgFactory.h"
#include "GaudiAlg/Tuples.h"
#include "GaudiKernel/PhysicalConstants.h"
#include "GaudiKernel/SystemOfUnits.h"
#include "Event/Track.h"
#include "Event/RecVertex.h"
#include "VeloDet/DeVelo.h"
#include <Linker/LinkedTo.h>
#include <Event/MCTrackInfo.h>
#include "Event/L0DUReport.h"
#include "GaudiUtils/HistoStats.h"
#include "AIDA/IHistogram1D.h"
#include <TRandom.h>
// local
#include "VertexCompare.h"
//-----------------------------------------------------------------------------
// Implementation file for class : VertexCompare
//-----------------------------------------------------------------------------
// Declaration of the Algorithm Factory
DECLARE_ALGORITHM_FACTORY(VertexCompare);
//=============================================================================
// Standard constructor, initializes variables
//=============================================================================
VertexCompare::VertexCompare(const std::string& name,
ISvcLocator* pSvcLocator)
: GaudiTupleAlg (name,pSvcLocator)
{
// declareProperty("nTracksToBeRecble", m_nTracksToBeRecble = 5);
declareProperty("produceNtuple", m_produceNtuple = true);
declareProperty("produceHistogram", m_produceHistogram = true);
// declareProperty("dzIsolated", m_dzIsolated = 10.0 * Gaudi::Units::mm);
// declareProperty("inputTracksName", m_inputTracksName = LHCb::TrackLocation::Default);
declareProperty("inputVerticesName1", m_inputVerticesName1 = LHCb::RecVertexLocation::Primary);
declareProperty("inputVerticesName2", m_inputVerticesName2 = LHCb::RecVertexLocation::Primary);
}
//=============================================================================
// Destructor
//=============================================================================
VertexCompare::~VertexCompare() {};
//=============================================================================
// Initialization
//=============================================================================
StatusCode VertexCompare::initialize() {
StatusCode sc = GaudiTupleAlg::initialize(); // Must be executed first
if(sc.isFailure()) debug() << "==> Initialize" << endmsg;
m_forcedBtool = tool<IForcedBDecayTool> ( "ForcedBDecayTool", this );
if( ! m_forcedBtool ) {
fatal() << "Unable to retrieve ForcedBDecayTool tool "<< endreq;
return StatusCode::FAILURE;
}
// Access PVOfflineTool
// m_pvsfit = tool<IPVOfflineTool>("PVOfflineTool",this);
// if(!m_pvsfit) {
// err() << "Unable to retrieve the PVOfflineTool" << endmsg;
// return StatusCode::FAILURE;
// }
m_nevt = 0;
m_nVtx = 0;
m_nPartVtx = 0;
return StatusCode::SUCCESS;
};
//=============================================================================
// Main execution
//=============================================================================
StatusCode VertexCompare::execute() {
debug() << "==> Execute" << endmsg;
std::vector<LHCb::RecVertex*> vecOfVertices1;
std::vector<LHCb::RecVertex*> vecOfVertices2;
bool vertices_ok = getInputVertices(vecOfVertices1, vecOfVertices2);
std::vector<LHCb::RecVertex> fullVrt1, fullVrt2;
if (!vertices_ok) {
return StatusCode::SUCCESS; // return SUCCESS anyway
}
if (debugLevel()) debug() <<
" Vtx Properities x y z chi2/ndof ntracks"
<< endmsg;
// Fill reconstructed PV info
std::vector<LHCb::RecVertex*>::iterator itRecV1;
for(itRecV1 = vecOfVertices1.begin(); vecOfVertices1.end() != itRecV1;
itRecV1++) {
LHCb::RecVertex* pv;
pv = *itRecV1;
fullVrt1.push_back(*pv);
// int nTracks = pv->tracks().size();
if (debugLevel()) debug() << m_inputVerticesName1 << endmsg;
if (debugLevel()) debug() << " " << pv->position().x() << " "
<< pv->position().y() << " " << pv->position().z()
<< " " << pv->chi2PerDoF() << " " << pv->tracks().size() << endmsg;
} // end of loop over vertices1
std::vector<LHCb::RecVertex*>::iterator itRecV2;
for(itRecV2 = vecOfVertices2.begin(); vecOfVertices2.end() != itRecV2;
itRecV2++) {
LHCb::RecVertex* pv;
pv = *itRecV2;
fullVrt2.push_back(*pv);
// int nTracks = pv->tracks().size();
if (debugLevel()) debug() << m_inputVerticesName2 << endmsg;
if (debugLevel()) debug() << " "
<< pv->position().x() << " "
<< pv->position().y() << " "
<< pv->position().z() << " "
<< pv->chi2PerDoF() << " "
<< pv->tracks().size() << endmsg;
} // end of loop over vertices2
if (debugLevel()) debug() << "fullVrt1 size " << fullVrt1.size() << endmsg;
if (debugLevel()) debug() << "fullVrt2 size " << fullVrt2.size() << endmsg;
int size_diff = fullVrt1.size() - fullVrt2.size();
if(m_produceHistogram) {
plot(double(size_diff), 1001, "size_diff", -5.5, 5.5, 11);
}
if(m_produceNtuple) {
Tuple myTuple_evt = nTuple(101, "PV_nTuple_evt", CLID_ColumnWiseTuple);
myTuple_evt->column("size_diff", double(size_diff));
myTuple_evt->column("size_1", double(fullVrt1.size()));
myTuple_evt->column("size_2", double(fullVrt2.size()));
myTuple_evt->write();
}
double dx = -99999.;
double dy = -99999.;
double dz = -99999.;
double errx = -99999.;
double erry = -99999.;
double errz = -99999.;
int ntracks_part = 0;
int ntracks_part2 = 0;
int dtracks = 0;
std::vector<LHCb::RecVertex>::iterator fullIter;
int oIt=0;
std::vector<int> link;
if (fullVrt1.size()!=0 && fullVrt2.size()!=0) matchByDistance(fullVrt1, fullVrt2, link);
else return StatusCode::SUCCESS;
for (fullIter=fullVrt1.begin(); fullIter!=fullVrt1.end(); fullIter++){
LHCb::RecVertex vrtf = *fullIter;
m_nVtx++;
if (link.at(oIt)==-1) continue;
Gaudi::SymMatrix3x3 covPV_part1 = vrtf.covMatrix();
double sigx_part1 = (covPV_part1(0,0));
double sigy_part1 = (covPV_part1(1,1));
double sigz_part1 = (covPV_part1(2,2));
Gaudi::SymMatrix3x3 covPV_part2 = fullVrt2.at(link.at(oIt)).covMatrix();
double sigx_part2 = (covPV_part2(0,0));
double sigy_part2 = (covPV_part2(1,1));
double sigz_part2 = (covPV_part2(2,2));
double x1 = vrtf.position().x();
double y1 = vrtf.position().y();
double z1 = vrtf.position().z();
double x2 = fullVrt2.at(link.at(oIt)).position().x();
double y2 = fullVrt2.at(link.at(oIt)).position().y();
double z2 = fullVrt2.at(link.at(oIt)).position().z();
dx = vrtf.position().x()-fullVrt2.at(link.at(oIt)).position().x();
dy = vrtf.position().y()-fullVrt2.at(link.at(oIt)).position().y();
dz = vrtf.position().z()-fullVrt2.at(link.at(oIt)).position().z();
errx = sqrt(sigx_part1 + sigx_part2);
erry = sqrt(sigy_part1 + sigy_part2);
errz = sqrt(sigz_part1 + sigz_part2);
ntracks_part = vrtf.tracks().size();
ntracks_part2 = fullVrt2.at(link.at(oIt)).tracks().size();
dtracks = ntracks_part - ntracks_part2;
if(m_produceHistogram) {
plot(dx, 1021, "dx", -0.25, 0.25, 50);
plot(dy, 1022, "dy", -0.25, 0.25, 50);
plot(dz, 1023, "dz", -1.5, 1.5, 60);
plot(x1, 2021, "x1", -0.25, 0.25, 50);
plot(y1, 2022, "y1", -0.25, 0.25, 50);
plot(z1, 2023, "z1", -100., 100., 50);
plot(x2, 3021, "x2", -0.25, 0.25, 50);
plot(y2, 3022, "y2", -0.25, 0.25, 50);
plot(z2, 3023, "z2", -100., 100., 50);
plot(std::sqrt(sigx_part1), 4011, "x err 1",0., .1, 50);
plot(std::sqrt(sigy_part1), 4012, "y err 1",0., .1, 50);
plot(std::sqrt(sigz_part1), 4013, "z err 1",0., .5, 50);
plot(std::sqrt(sigx_part2), 4021, "x err 2",0., .1, 50);
plot(std::sqrt(sigy_part2), 4022, "y err 2",0., .1, 50);
plot(std::sqrt(sigz_part2), 4023, "z err 2",0., .5, 50);
plot(dx/errx, 1031, "pullx", -5., 5., 50);
plot(dy/erry, 1032, "pully", -5., 5., 50);
plot(dz/errz, 1033, "pullz", -5., 5., 50);
plot(double(ntracks_part) , 1041, "ntracks_part1", 0., 150., 50);
plot(double(ntracks_part2), 1042, "ntracks_part2", 0., 150., 50);
plot(double(dtracks), 1043, "dtracks", 0., 150., 50);
// dz to get false vertex rate from data
if ( vecOfVertices1.size() > 1 ) {
for ( unsigned int i1 = 0; i1 < vecOfVertices1.size(); i1++ ) {
for ( unsigned int i2 = 0; i2 < vecOfVertices1.size(); i2++ ) {
if (i2 != i1) {
double vdz = vecOfVertices1[i1]->position().z() - vecOfVertices1[i2]->position().z();
plot(vdz,1201,"dz vertices 1", -150.,150.,100);
}
}
}
}
if ( vecOfVertices2.size() > 1 ) {
for ( unsigned int i1 = 0; i1 < vecOfVertices2.size(); i1++ ) {
for ( unsigned int i2 = 0; i2 < vecOfVertices2.size(); i2++ ) {
if (i2 != i1) {
double vdz = vecOfVertices2[i1]->position().z() - vecOfVertices2[i2]->position().z();
plot(vdz, 1202, "dz vertices 2", -150.,150.,100);
}
}
}
}
}
if(m_produceNtuple) {
Tuple myTuple = nTuple(102, "PV_nTuple", CLID_ColumnWiseTuple);
myTuple->column("ntracks_part", double(ntracks_part));
myTuple->column("ntracks_part2", double(ntracks_part2));
myTuple->column("dtracks", double(dtracks));
myTuple->column("dx", dx);
myTuple->column("dy", dy);
myTuple->column("dz", dz);
myTuple->column("x1", x1);
myTuple->column("y1", y1);
myTuple->column("z1", z1);
myTuple->column("x2", x2);
myTuple->column("y2", y2);
myTuple->column("z2", z2);
myTuple->column("errx", errx);
myTuple->column("erry", erry);
myTuple->column("errz", errz);
myTuple->write();
}
oIt++;
}
return StatusCode::SUCCESS;
};
//=============================================================================
// Finalize
//=============================================================================
StatusCode VertexCompare::finalize() {
debug() << "==> Finalize" << endmsg;
info() << " ============================================" << endreq;
info() << " Efficiencies for reconstructed vertices: " << endreq;
info() << " ============================================" << endreq;
info() << " " << endreq;
info() << " There are " << m_nVtx
<< " pairs of vertices in processed events" << endreq;
// info() << " PV is isolated if dz to closest reconstructible MC PV > "
// << m_dzIsolated << " mm" << endreq;
// std::string ff = "by counting tracks";
// /* if ( !m_matchByTracks )*/ ff = "by dz distance";
// info() << " Two splited vertices matched: "
// << ff << endreq;
// info() << " " << endreq;
// printRat("All", m_nPartVtx, m_nVtx );
const AIDA::IHistogram1D* dx = histo( HistoID(1021) ) ;
const AIDA::IHistogram1D* pullx = histo( HistoID(1031) ) ;
const AIDA::IHistogram1D* dy = histo( HistoID(1022) ) ;
const AIDA::IHistogram1D* pully = histo( HistoID(1032) ) ;
const AIDA::IHistogram1D* dz = histo( HistoID(1023) ) ;
const AIDA::IHistogram1D* pullz = histo( HistoID(1033) ) ;
if( dx ) {
info() << " ---------------------------------------" << endreq;
info() << "dx: "
<< format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
dx->mean(), Gaudi::Utils::HistoStats::meanErr(dx),
dx->rms(), Gaudi::Utils::HistoStats::rmsErr(dx)) << endreq ;
}
if( dy ) {
info() << "dy: "
<< format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
dy->mean(), Gaudi::Utils::HistoStats::meanErr(dy),
dy->rms(), Gaudi::Utils::HistoStats::rmsErr(dy)) << endreq ;
}
if( dz ) {
info() << "dz: "
<< format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
dz->mean(), Gaudi::Utils::HistoStats::meanErr(dz),
dz->rms(), Gaudi::Utils::HistoStats::rmsErr(dz)) << endreq ;
}
info() << " ---------------------------------------" << endreq;
if( pullx ) {
info() << "pullx: "
<< format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
pullx->mean(), Gaudi::Utils::HistoStats::meanErr(pullx),
pullx->rms(), Gaudi::Utils::HistoStats::rmsErr(pullx)) << endreq ;
}
if( pully ) {
info() << "pully: "
<< format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
pully->mean(), Gaudi::Utils::HistoStats::meanErr(pully),
pully->rms(), Gaudi::Utils::HistoStats::rmsErr(pully)) << endreq ;
}
if( pullz ) {
info() << "pullz: "
<< format( "mean = %5.3f +/- %5.3f, RMS = %5.3f +/- %5.3f",
pullz->mean(), Gaudi::Utils::HistoStats::meanErr(pullz),
pullz->rms(), Gaudi::Utils::HistoStats::rmsErr(pullz)) << endreq ;
}
info() << " ============================================" << endreq;
//
return GaudiTupleAlg::finalize(); // Must be called after all other actions
}
//=============================================================================
// Match vertices by distance
//=============================================================================
void VertexCompare::matchByDistance(std::vector<LHCb::RecVertex>& vertfull,
std::vector<LHCb::RecVertex>& verthalf, std::vector<int>& link) {
if (verthalf.size()>vertfull.size() && debugLevel()) debug() << "half.size > full.size" << endmsg;
for(int imc=0; imc<(int)vertfull.size(); imc++) {
// if ( mcpvvec[imc].indexRecPVInfo > -1) continue;
double mindist = 999999.;
int indexrec = -1;
for(int irec=0; irec<(int)verthalf.size(); irec++){
if (std::count(link.begin(), link.end(), irec)!=0) continue;
double dist = fabs(verthalf.at(irec).position().z()-vertfull.at(imc).position().z());
if(dist < mindist) {
mindist = dist;
indexrec = irec;
}
}
if(debugLevel()) debug() << "original vertex " << imc << " linked to " << indexrec << " half vertex." << endmsg;
link.push_back(indexrec);
}
for(int imc=0; imc<(int)vertfull.size(); imc++) {
int count = std::count(link.begin(), link.end(), imc);
if (count >1 && debugLevel()) debug() << "linked twice to vertex " << imc << endmsg;
}
}
//=============================================================================
// printRat
//=============================================================================
void VertexCompare::printRat(std::string mes, int a, int b) {
double rat = 0.;
if(b>0) rat = 1.0*a/b;
// reformat message
unsigned int len = 20;
std::string pmes = mes;
while(pmes.length() < len) {
pmes+=" ";
}
pmes+= " : ";
info() << pmes << format(" %6.3f ( %7d / %8d )", rat,a,b) << endreq;
}
//=============================================================================
// Get input vertices
//=============================================================================
bool VertexCompare::getInputVertices( std::vector<LHCb::RecVertex*>& vecOfVertices1,
std::vector<LHCb::RecVertex*>& vecOfVertices2 ) {
std::string verticesName1 = m_inputVerticesName1;
std::string verticesName2 = m_inputVerticesName2;
if ( m_inputVerticesName1 == "Offline" ) {
verticesName1 = LHCb::RecVertexLocation::Primary;
} else if ( m_inputVerticesName1 == "3D" ) {
verticesName1 = LHCb::RecVertexLocation::Velo3D;
} else if ( m_inputVerticesName1 == "2D" ) {
verticesName1 = LHCb::RecVertexLocation::Velo2D;
}
if ( m_inputVerticesName2 == "Offline" ) {
verticesName2 = LHCb::RecVertexLocation::Primary;
} else if ( m_inputVerticesName2 == "3D" ) {
verticesName2 = LHCb::RecVertexLocation::Velo3D;
} else if ( m_inputVerticesName2 == "2D" ) {
verticesName2 = LHCb::RecVertexLocation::Velo2D;
}
if ( verticesName1 == "none" ) {
debug() << " Vertices 1 not specified " << verticesName1 << endreq;
return false;
}
if ( verticesName2 == "none" ) {
debug() << " Vertices 2 not specified " << verticesName2 << endreq;
return false;
}
LHCb::RecVertices* recoVertices1;
LHCb::RecVertices* recoVertices2;
try {
recoVertices1 = get<LHCb::RecVertices>( verticesName1 );
} catch (...) {
debug() << " No vertices at " << verticesName1 << endreq;
return false;
}
try {
recoVertices2 = get<LHCb::RecVertices>( verticesName2 );
} catch (...) {
debug() << " No vertices at " << verticesName2 << endreq;
return false;
}
std::vector<LHCb::RecVertex*>::const_iterator itVer1;
for ( itVer1 = recoVertices1->begin(); recoVertices1->end() != itVer1; itVer1++ ) {
LHCb::RecVertex* ppv1 = (*itVer1);
vecOfVertices1.push_back(ppv1);
}
std::vector<LHCb::RecVertex*>::const_iterator itVer2;
for ( itVer2 = recoVertices2->begin(); recoVertices2->end() != itVer2; itVer2++ ) {
LHCb::RecVertex* ppv2 = (*itVer2);
vecOfVertices2.push_back(ppv2);
}
return true;
}
// $Id: VertexCompare.h,v 1.1 2010/03/01 15:58:14 pmorawsk Exp $
#ifndef VERTEXCOMPARE_H
#define VERTEXCOMPARE_H 1
// Include files
#include "GaudiAlg/GaudiAlgorithm.h"
#include "GaudiAlg/GaudiTupleAlg.h"
#include "MCInterfaces/IForcedBDecayTool.h"
#include "TrackInterfaces/IPVOfflineTool.h"
// #include "TrackInterfaces/IPVSeeding.h"
typedef struct {
LHCb::MCVertex* pMCPV; // pointer to MC PV
int nRecTracks; // number of reconstructed tracks from this MCPV
int nRecBackTracks; // number of reconstructed backward tracks
int indexRecPVInfo; // index to reconstructed PVInfo (-1 if not reco)
int nCorrectTracks; // correct tracks belonging to reconstructed PV
int multClosestMCPV; // multiplicity of closest reconstructable MCPV
double distToClosestMCPV; // distance to closest reconstructible MCPV
std::vector<LHCb::MCParticle*> m_mcPartInMCPV;
std::vector<LHCb::Track*> m_recTracksInMCPV;
} MCPVInfo;
typedef struct {
public:
int nTracks; // number of tracks in a vertex
int nBackTracks; // number of backward tracks in a vertex
Gaudi::XYZPoint position; // position
Gaudi::XYZPoint positionSigma; // position sigmas
int indexMCPVInfo; // index to MCPVInfo
LHCb::RecVertex* pRECPV; // pointer to REC PV
bool rec1, rec2;
} RecPVInfo;
class VertexCompare : public GaudiTupleAlg {
public:
/// Standard constructor
VertexCompare( const std::string& name, ISvcLocator* pSvcLocator );
virtual ~VertexCompare( ); ///< Destructor
virtual StatusCode initialize(); ///< Algorithm initialization
virtual StatusCode execute (); ///< Algorithm execution
virtual StatusCode finalize (); ///< Algorithm finalization
protected:
bool debugLevel() const { return msgLevel( MSG::DEBUG ) || msgLevel( MSG::VERBOSE ) ; }
private:
IForcedBDecayTool* m_forcedBtool;
// IPVOfflineTool* m_pvsfit;
// int m_nTracksToBeRecble;
bool m_produceHistogram;
bool m_produceNtuple;
// double m_dzIsolated;
// std::string m_inputTracksName;
std::string m_inputVerticesName1;
std::string m_inputVerticesName2;
// std::string m_pvSeedingName;
int m_nevt, m_nVtx;
int m_nPartVtx ;
void printRat(std::string mes, int a, int b);
// bool getInputTracks( std::vector<LHCb::Track*>& vecOfTracks, std::string& trackLoc);
bool getInputVertices( std::vector<LHCb::RecVertex*>& vecOfVertices1 , std::vector<LHCb::RecVertex*>& vecOfVertices2 );
void matchByDistance(std::vector<LHCb::RecVertex>& verthalf, std::vector<LHCb::RecVertex>& vertfull, std::vector<int>& link);
};
#endif // VERTEXCOMPARE_H
project REC
use LHCB LHCB_v32r1
use LHCB LHCB_v32r1p1
build_strategy with_installarea
setup_strategy root
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment