diff --git a/CMakeLists.txt b/CMakeLists.txt index 450a595..b142ed8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -53,11 +53,7 @@ include_directories( "${FCCEDM_INCLUDE_DIRS}" ) - add_definitions(-Wno-unused-variable -Wno-unused-parameter -pthread) -#--temporary fix of inconsistency in ROOT CMake macros -#set(ROOT_genreflex_cmd ${ROOT_genreflex_CMD}) - #--- enable unit testing capabilities ------------------------------------------ include(CTest) @@ -67,8 +63,6 @@ include(CTest) include(cmake/papasDoxygen.cmake) #endif() - - #--- add version files --------------------------------------------------------- configure_file(${CMAKE_CURRENT_SOURCE_DIR}/papasVersion.h ${CMAKE_CURRENT_BINARY_DIR}/papasVersion.h ) @@ -78,14 +72,11 @@ install(FILES ${CMAKE_CURRENT_BINARY_DIR}/papasVersion.h #--- add CMake infrastructure -------------------------------------------------- include(cmake/papasCreateConfig.cmake) - #--- add license files --------------------------------------------------------- install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/LICENSE ${CMAKE_CURRENT_SOURCE_DIR}/NOTICE DESTINATION ${CMAKE_INSTALL_PREFIX}) - -#add_subdirectory(papascpp) add_subdirectory(papaslib) add_subdirectory(examples) add_subdirectory(tests) diff --git a/cmake/papasConfig.cmake.in b/cmake/papasConfig.cmake.in index 566762c..3c16c64 100644 --- a/cmake/papasConfig.cmake.in +++ b/cmake/papasConfig.cmake.in @@ -20,4 +20,4 @@ FIND_LIBRARY( papas_LIBRARIES NAMES papas PATHS INCLUDE( FindPackageHandleStandardArgs ) FIND_PACKAGE_HANDLE_STANDARD_ARGS( papas DEFAULT_MSG papas_LIBRARY_DIR papas_INCLUDE_DIRS papas_LIBRARIES ) -ADD_LIBRARY(papasexamplelibrary SHARED IMPORTED) +#ADD_LIBRARY(papasexamplelibrary SHARED IMPORTED) diff --git a/examples/PythiaConnector.cpp b/examples/PythiaConnector.cpp index fa13208..3e92247 100644 --- a/examples/PythiaConnector.cpp +++ b/examples/PythiaConnector.cpp @@ -1,4 +1,3 @@ - // // PythiaConnector.cpp // papas @@ -17,12 +16,15 @@ #include "datamodel/ParticleCollection.h" #include "utilities/ParticleUtils.h" +#include "papas/datatypes/Helix.h" #include "papas/datatypes/Particle.h" -#include "papas/simulation/Simulator.h" -#include "papas/utility/PDebug.h" - +#include "papas/datatypes/Path.h" +#include "papas/detectors/Detector.h" +#include "papas/detectors/Field.h" #include "papas/display/PFApp.h" +#include "papas/simulation/Simulator.h" #include "papas/utility/Log.h" +#include "papas/utility/PDebug.h" #include #include @@ -47,7 +49,8 @@ PythiaConnector::PythiaConnector(const char* fname) : m_store(podio::EventStore( } void PythiaConnector::makePapasParticlesFromGeneratedParticles(const fcc::MCParticleCollection* ptcs, - papas::Particles& particles) { + papas::Particles& particles, + const papas::Detector& detector) { // turns pythia particles into Papas particles and lodges them in the history TLorentzVector tlv; int countp = 0; @@ -81,6 +84,15 @@ void PythiaConnector::makePapasParticlesFromGeneratedParticles(const fcc::MCPart if (tlv.Pt() > 1e-5 && (abs(pdgid) != 12) && (abs(pdgid) != 14) && (abs(pdgid) != 16)) { papas::Particle particle(pdgid, (double)ptc.core().charge, tlv, particles.size(), 's', startVertex, ptc.core().status); + // set the particles papas path (allows particles to be const when passed to simulator) + std::shared_ptr ppath; + if (fabs(particle.charge()) < 0.5) { + ppath = std::make_shared(particle.p4(), particle.startVertex(), particle.charge()); + } else { + ppath = std::make_shared(particle.p4(), particle.startVertex(), particle.charge(), + detector.field()->getMagnitude()); + } + particle.setPath(ppath); particles.emplace(particle.id(), particle); papas::PDebug::write("Made {}", particle); } @@ -92,14 +104,15 @@ void PythiaConnector::processEvent(unsigned int eventNo, papas::PapasManager& pa // make a papas particle collection from the next event // then run simulate and reconstruct m_reader.goToEvent(eventNo); - papasManager.clear(); papasManager.setEventNo(eventNo); - const fcc::MCParticleCollection* ptcs(nullptr); + const fcc::MCParticleCollection* ptcs; if (m_store.get("GenParticle", ptcs)) { try { + papasManager.clear(); papas::Particles& genParticles = papasManager.createParticles(); - makePapasParticlesFromGeneratedParticles(ptcs, genParticles); - papasManager.simulate(genParticles); + makePapasParticlesFromGeneratedParticles(ptcs, genParticles, papasManager.detector()); + papasManager.addParticles(genParticles); + papasManager.simulate('s'); papasManager.mergeClusters("es"); papasManager.mergeClusters("hs"); papasManager.buildBlocks('m', 'm', 's'); @@ -130,8 +143,8 @@ void PythiaConnector::writeParticlesROOT(const char* fname, const papas::Particl auto& evinfocoll = m_store.create("evtinfo"); auto& pcoll = m_store.create("GenParticle"); - writer.registerForWrite("evtinfo"); - writer.registerForWrite("GenParticle"); + writer.registerForWrite("evtinfo"); + writer.registerForWrite("GenParticle"); auto evinfo = fcc::EventInfo(); // evinfocoll.create(); evinfo.number(eventno); @@ -161,8 +174,8 @@ void PythiaConnector::writeClustersROOT(const char* fname, const papas::Clusters auto& evinfocoll = m_store.create("evtinfo"); auto& ccoll = m_store.create("Cluster"); - writer.registerForWrite("evtinfo"); - writer.registerForWrite("Cluster"); + writer.registerForWrite("evtinfo"); + writer.registerForWrite("Cluster"); auto evinfo = fcc::EventInfo(); // evinfocoll.create(); evinfo.number(eventno); diff --git a/examples/PythiaConnector.h b/examples/PythiaConnector.h index ecab744..f89e479 100644 --- a/examples/PythiaConnector.h +++ b/examples/PythiaConnector.h @@ -51,7 +51,8 @@ class PythiaConnector { ///< Takes pythia particles and creates Papas type particles adding them into /// an empty Particles collection - void makePapasParticlesFromGeneratedParticles(const fcc::MCParticleCollection* ptcs, papas::Particles& particles); + void makePapasParticlesFromGeneratedParticles(const fcc::MCParticleCollection* ptcs, papas::Particles& particles, + const papas::Detector& detector); papas::Clusters ConvertClustersToPapas(const fcc::CaloClusterCollection& fccClusters, float size, papas::IdCoder::ItemType itemtype, char subtype) const; void AddClustersToEDM(const papas::Clusters& papasClusters, fcc::CaloClusterCollection& fccClusters); diff --git a/examples/example_loop.cpp b/examples/example_loop.cpp index 1ae2dcd..12cc906 100644 --- a/examples/example_loop.cpp +++ b/examples/example_loop.cpp @@ -19,8 +19,8 @@ #include int main(int argc, char* argv[]) { + papas::PDebug::File("physics.txt"); - // randomgen::setEngineSeed(0xdeadbeef); // make results reproduceable rootrandom::Random::seed(0xdeadbeef); if (argc != 2) { @@ -41,14 +41,14 @@ int main(int argc, char* argv[]) { papas::PapasManager papasManager(CMSDetector); unsigned int eventNo = 0; - unsigned int nEvents = 10; + unsigned int nEvents = 10000; auto start = std::chrono::steady_clock::now(); for (unsigned i = eventNo; i < eventNo + nEvents; ++i) { papas::PDebug::write("Event: {}", i); - if (i % 10 == 0) { + if (i % 100 == 0) { std::cout << "reading event " << i << std::endl; } if (i == eventNo) diff --git a/examples/example_plot.cpp b/examples/example_plot.cpp index 7224fc9..d302103 100644 --- a/examples/example_plot.cpp +++ b/examples/example_plot.cpp @@ -13,6 +13,7 @@ #include "papas/detectors/CMS.h" #include "papas/reconstruction/PapasManager.h" #include "papas/utility/PDebug.h" +#include "papas/utility/TRandom.h" #include @@ -20,10 +21,17 @@ using namespace papas; int main(int argc, char* argv[]) { - if (argc != 2) { - std::cerr << "Usage: ./mainexe filename" << std::endl; + rootrandom::Random::seed(0xdeadbeef); + + if (argc < 2) { + std::cerr << "Usage: ./example_plot filename [logname]" << std::endl; return 1; } + if (argc == 3) { + const char* lname = argv[2]; + PDebug::File(lname); // physics debug output + } + const char* fname = argv[1]; PythiaConnector pythiaConnector(fname); @@ -35,5 +43,6 @@ int main(int argc, char* argv[]) { pythiaConnector.processEvent(eventNo, papasManager); TApplication tApp("theApp", &argc, argv); pythiaConnector.displayEvent(papasManager); + papasManager.clear(); return EXIT_SUCCESS; } diff --git a/init.sh b/init.sh index 4c417d3..5b008e2 100644 --- a/init.sh +++ b/init.sh @@ -2,17 +2,18 @@ platform='unknown' sw_afs=0 unamestr=`uname` -export LCGPATH=/afs/cern.ch/sw/lcg/views/LCG_83/x86_64-slc6-gcc49-opt -#if [ -z ${FCCPAPASCPP+x} ]; then +if [ -z ${FCCPAPASCPP+x} ]; then unset $FCCPAPASCPP echo "FCCPAPASCPP is unset, setting to $PWD/install"; export FCCPAPASCPP=$PWD/install export PATH=$FCCPAPASCPP/bin:$PATH -#fi +fi if [[ "$unamestr" == 'Linux' ]]; then platform='Linux' export LD_LIBRARY_PATH=$FCCPAPASCPP/lib:$LD_LIBRARY_PATH + source /cvmfs/fcc.cern.ch/sw/0.8.1/init_fcc_stack.sh $1 +#export LCGPATH=/afs/cern.ch/sw/lcg/views/LCG_83/x86_64-slc6-gcc49-opt elif [[ "$unamestr" == 'Darwin' ]]; then platform='Darwin' export DYLD_LIBRARY_PATH=$FCCPAPASCPP/lib:$DYLD_LIBRARY_PATH diff --git a/papas/datatypes/Cluster.h b/papas/datatypes/Cluster.h index aaa7bc6..2df4fde 100644 --- a/papas/datatypes/Cluster.h +++ b/papas/datatypes/Cluster.h @@ -26,7 +26,7 @@ class Cluster { @param[in] id identifier type of cluster eg kEcalCluster or kHcalCluster @param[in] subtype single char describing type of cluster eg s = smeared, t= true, m = merged */ - Cluster(double energy, const TVector3& position, double size_m, unsigned int index, IdCoder::ItemType id, + Cluster(double energy, const TVector3& position, double size_m, uint32_t index, IdCoder::ItemType id, char subtype = 't'); /** Constructor: makes new cluster with a new id based on a copy of an existing cluster. The new id must be provided. @@ -37,7 +37,7 @@ class Cluster { @param[in] val the value that will be used when creating the Cluster identifier and which is used for sorting. When creating a merged cluster it should ideally be set to the total eneergy of the cluster */ - Cluster(const Cluster& cluster, unsigned int index, IdCoder::ItemType type, char subtype = 'u', float val = 0.0); + Cluster(const Cluster& cluster, uint32_t index, IdCoder::ItemType type, char subtype = 'u', float val = 0.0); /** Constructor: makes new merged cluster @param[in] overlappingClusters list of clusters to be merged, they must have same type and subtype and must be @@ -46,7 +46,7 @@ class Cluster { @param[in] index of the collection into which the cluster is to be stored @param[in] subtype subtype of cluster to be created eg 'm' for merged, */ - Cluster(std::list overlappingClusters, unsigned int index, char subtype = 'm'); + Cluster(std::list overlappingClusters, uint32_t index, char subtype = 'm'); Cluster() = default; Cluster(Cluster&& c); // needed for unordered_map Cluster(const Cluster& cluster) = default; // needed for unordered_map @@ -54,7 +54,7 @@ class Cluster { double angularSize() const; ///< The angle that the cluster boundary makes (not valid for merged clusters) double size() const; ///< The radius of the cluster double pt() const { - return m_energy * m_position.Perp(); + return m_energy * m_position.Unit().Perp(); } ///< Transverse momentum (magnitude of p3 in transverse plane) double energy() const { return m_energy; } ///< Energy double eta() const { return m_position.Eta(); } ///< Pseudo-rapidity (-ln(tan self._tlv.Theta()/2)) diff --git a/papas/datatypes/DefinitionsCollections.h b/papas/datatypes/DefinitionsCollections.h index ddb8d62..da28128 100644 --- a/papas/datatypes/DefinitionsCollections.h +++ b/papas/datatypes/DefinitionsCollections.h @@ -18,8 +18,8 @@ class Track; class Edge; class Particle; -typedef std::list ListParticles; ///< list of Particles -typedef std::unordered_map Edges; ///< collection of Edge objects +typedef std::list ListParticles; ///< list of Particles +typedef std::unordered_map Edges; ///< collection of Edge objects #if WITHSORT typedef std::set> Ids; ///< set containing Identifiers #else diff --git a/papas/datatypes/Event.h b/papas/datatypes/Event.h index 6750cd1..75941c5 100644 --- a/papas/datatypes/Event.h +++ b/papas/datatypes/Event.h @@ -14,8 +14,8 @@ #include namespace papas { - - class Cluster; + +class Cluster; /** * @brief The Event stores pointers to collections of Clusters, Tracks, Blocks, Particles in its @@ -58,7 +58,7 @@ namespace papas { class Event { public: /// @brief Constructor - Event(std::shared_ptr hist = std::make_shared(Nodes())); + Event(Nodes& hist); /** * @brief adds a pointer to a Clusters collection (unordered map) into the Event @@ -92,7 +92,7 @@ class Event { * @brief makes history in Event point to an external history object * @param[in] history unordered map of Nodes, * */ - void setHistory(Nodes& history) { m_history = std::make_shared(history); } + void setHistory(Nodes& history) { m_history = history; } /** * @brief adds new history into existing papasevent history @@ -212,7 +212,7 @@ class Event { /** * @brief returns the merged history */ - std::shared_ptr history() const { return m_history; } + const Nodes& history() const { return m_history; } /** * @brief resets everything, deletes all the clusters, tracks etc etc @@ -233,7 +233,7 @@ class Event { * @param[in] typeAndSubtype The type and subtype of a collection eg "em" for ecal merged */ Ids collectionIds(const std::string& typeAndSubtype) const; - + std::string info() const; ///< text descriptor private: /** * @brief templated class method used by the AddCollection methods to check that typeAndSubype match and that @@ -255,25 +255,22 @@ class Event { ParticlesFolder m_particlesFolder; /// Unordered map of pointers to unordered map of (concrete) Blocks BlocksFolder m_blocksFolder; - std::shared_ptr m_history; ///< points to the merged history (built from the sucessive histories) - Clusters m_emptyClusters; /// void Event::addCollectionToFolderInternal( const std::unordered_map& collection, std::unordered_map*>& folder) { - Identifier firstId = 0; if (collection.size() == 0) return; + Identifier firstId = collection.begin()->first; + if (hasCollection(firstId)) throw "Collection already exists"; for (const auto& it : collection) { - if (!firstId) { - firstId = it.first; - if (hasCollection(firstId)) throw "Collection already exists"; - } if (IdCoder::typeAndSubtype(it.first) != IdCoder::typeAndSubtype(firstId)) { std::cout << IdCoder::pretty(it.first) << " : " << IdCoder::pretty(firstId) << std::endl; throw "more than one typeandSubtype found in collection"; diff --git a/papas/datatypes/IdCoder.h b/papas/datatypes/IdCoder.h index 63cf7f2..173ff48 100644 --- a/papas/datatypes/IdCoder.h +++ b/papas/datatypes/IdCoder.h @@ -56,7 +56,7 @@ class IdCoder { @param[in] value: a float representing energy or momentum etc @return identifier */ - static Identifier makeId(unsigned int index, ItemType type, char subtype = 'u', float value = 0.0); + static Identifier makeId(uint32_t index, ItemType type, char subtype = 'u', float value = 0.0); /** returns the item type of the identifier This is one of: None = 0, kEcalCluster = 1, kHcalCluster, kTrack, kParticle, kBlock @@ -90,12 +90,12 @@ class IdCoder { @param[in] id: identifier @return the index */ - static unsigned int index(Identifier id); ///< Returns encoded index + static uint32_t index(Identifier id); ///< Returns encoded index /** Takes an identifier and returns a unique id component of it (excludes value information) @param[in] id: identifier @return the index */ - static unsigned int uniqueId(Identifier id); ///< Returns encoded unique id + static uint32_t uniqueId(Identifier id); ///< Returns encoded unique id static char typeLetter(Identifier id); ///< One letter short code eg 'e' for ecal, 't' for track, 'x' for unknown static std::string typeAndSubtype(Identifier id); ///< Two letter string of type and subtype eg "em" @@ -147,12 +147,12 @@ class IdCoder { static int bitshift() { return m_bitshift; } private: - static const unsigned int m_bitshift1 = 61; ///< encoding parameter - static const unsigned int m_bitshift2 = 53; ///< encoding parameter - static const unsigned int m_bitshift = 21; ///< encoding parameter (max size of counter) + static const uint32_t m_bitshift1 = 61; ///< encoding parameter + static const uint32_t m_bitshift2 = 53; ///< encoding parameter + static const uint32_t m_bitshift = 21; ///< encoding parameter (max size of counter) /// checks that the identifier can be correctly decoded - static bool checkValid(Identifier id, ItemType type, char subt, float val, unsigned int uid); - static bool checkUIDValid(Identifier id, unsigned int uniqueid); + static bool checkValid(Identifier id, ItemType type, char subt, float val, uint32_t uid); + static bool checkUIDValid(Identifier id, uint32_t uniqueid); static uint64_t floatToBits(float value); /// convert float into binary static float bitsToFloat(uint64_t bits); /// convert binary into float }; diff --git a/papas/datatypes/Particle.h b/papas/datatypes/Particle.h index 2355b11..d003b6a 100644 --- a/papas/datatypes/Particle.h +++ b/papas/datatypes/Particle.h @@ -30,7 +30,7 @@ class Particle { @param[in] startvertex start vertex (3d point) @param[in] status Status (1 = stable) */ - Particle(int pdgid, double charge, const TLorentzVector& tlv, unsigned int index, char subtype, + Particle(int pdgid, double charge, const TLorentzVector& tlv, uint32_t index, char subtype, const TVector3& startvertex = TVector3(0., 0., 0.), double status = 1); const TLorentzVector& p4() const { return m_tlv; } ///< 4-momentum, px, py, pz, E TVector3 p3() const { return m_tlv.Vect(); } ///< 3-momentum px, py, pz, @@ -50,13 +50,10 @@ class Particle { bool status() const { return m_status; } /// path) { m_path = path; } - const std::shared_ptr path() const { return m_path; } ///< Return pointer to path - Identifier id() const { return m_id; } ///< unique Identifier for object - bool isElectroMagnetic() const; ///< Is it electroMagnetic + void setPath(std::shared_ptr path) { m_path = path; } ///< set the Particle path + const std::shared_ptr path() const { return m_path; } ///< Return pointer to path + Identifier id() const { return m_id; } ///< unique Identifier for object + bool isElectroMagnetic() const; ///< Is it electroMagnetic private: TLorentzVector m_tlv; ///<4-momentum, px, py, pz, E int m_pdgId; /// Points; ///< Map of path points indexed by position - /*Constructor - @param p4 4-momentum - @param origin position of start of path - @param magnetic field - */ + + /**Constructor + @param p4 4-momentum + @param origin position of start of path + @param magnetic field + */ Path(const TLorentzVector& p4, const TVector3& origin, double field = 0.); + /** Constructor */ Path(); + + /** Destructor*/ + virtual ~Path() = default; // needed to allow virtual destruction + /** Add a new point to path * @param layer for the new point which is used to index unordered map of points * @param vec new point to be added diff --git a/papas/datatypes/Track.h b/papas/datatypes/Track.h index a437d50..2f09e1d 100644 --- a/papas/datatypes/Track.h +++ b/papas/datatypes/Track.h @@ -27,7 +27,7 @@ class Track { @param[in] index used to create track identifier, normally the index of the collection to which the track will belong @param[in] the subtype of the track used in creating the identifer */ - Track(const TVector3& p3, double charge, const std::shared_ptr path, unsigned int index, char subtype = 'u'); + Track(const TVector3& p3, double charge, const std::shared_ptr path, uint32_t index, char subtype = 'u'); double energy() const { return m_p3.Mag(); } /// overlappingClusters, unsigned int index, char subtype) +Cluster::Cluster(std::list overlappingClusters, uint32_t index, char subtype) : m_subClusters(overlappingClusters) { Identifier firstId = 0; char type; @@ -87,8 +87,12 @@ std::ostream& operator<<(std::ostream& os, const Cluster& cluster) { os << "Cluster: " << std::setw(6) << std::left << IdCoder::pretty(cluster.id()) << ":" << cluster.id() << ": " << cluster.info(); os << " sub("; - for (const auto& c : cluster.subClusters()) { - os << IdCoder::pretty(c->id()) << ", "; + if (cluster.subClusters().empty()) // match python pdebug outputs + os << IdCoder::pretty(cluster.id()) << ", "; + else { + for (const auto& c : cluster.subClusters()) { + os << IdCoder::pretty(c->id()) << ", "; + } } os << ")"; return os; diff --git a/papas/datatypes/src/Event.cpp b/papas/datatypes/src/Event.cpp index aa472d0..8fd81b5 100644 --- a/papas/datatypes/src/Event.cpp +++ b/papas/datatypes/src/Event.cpp @@ -1,11 +1,16 @@ #include "papas/datatypes/Event.h" #include "papas/datatypes/Cluster.h" +#include //lxplus needs this +#include + +#include "spdlog/details/format.h" + namespace papas { /// Event holds pointers to collections of particles, clusters etc and the address of the history associated with /// an event -Event::Event(std::shared_ptr hist) +Event::Event(Nodes& hist) : m_ecalClustersFolder(), m_hcalClustersFolder(), m_tracksFolder(), @@ -154,7 +159,7 @@ void Event::extendHistory(const Nodes& history) { // the following adds this history into the papasevent history for (const auto& node : history) { for (const auto& c : node.second.children()) { - makeHistoryLink(node.first, c->value(), *m_history); + makeHistoryLink(node.first, c->value(), m_history); } } } @@ -165,6 +170,34 @@ void Event::clear() { m_tracksFolder.clear(); m_particlesFolder.clear(); m_blocksFolder.clear(); - m_history = 0; + m_history.clear(); +} + +std::string Event::info() const { + fmt::MemoryWriter out; + out.write("Papas::Event: {}\n", m_eventNo); + out.write("\thistory = {}", m_history.size()); + out.write("\n\tecals ="); + for (auto c : m_ecalClustersFolder) { + out.write(" {}({}) +", c.first, c.second->size()); + } + out.write("\n\thcals ="); + for (auto c : m_hcalClustersFolder) { + out.write(" {}({}) +", c.first, c.second->size()); + } + out.write("\n\ttracks ="); + for (auto c : m_tracksFolder) { + out.write(" {}({}) +", c.first, c.second->size()); + } + out.write("\n\tblocks ="); + for (auto c : m_blocksFolder) { + out.write(" {}({}) +", c.first, c.second->size()); + } + out.write("\n\tparticles ="); + for (auto c : m_particlesFolder) { + out.write(" {}({}) +", c.first, c.second->size()); + } + out.write("\n"); + return out.str(); } } diff --git a/papas/datatypes/src/HistoryHelper.cpp b/papas/datatypes/src/HistoryHelper.cpp index 2887c7e..e42d86e 100644 --- a/papas/datatypes/src/HistoryHelper.cpp +++ b/papas/datatypes/src/HistoryHelper.cpp @@ -8,7 +8,7 @@ HistoryHelper::HistoryHelper(const Event& event) : m_event(event) {} Ids HistoryHelper::linkedIds(Identifier id, DAG::enumVisitType direction) const { const auto& history = m_event.history(); - const auto& startnode = history->at(id); + const auto& startnode = history.at(id); DAG::BFSRecurseVisitor bfs; const auto& nodes = bfs.traverseNodes(startnode, direction); Ids ids; diff --git a/papas/datatypes/src/IdCoder.cpp b/papas/datatypes/src/IdCoder.cpp index 004b53c..4e09a04 100644 --- a/papas/datatypes/src/IdCoder.cpp +++ b/papas/datatypes/src/IdCoder.cpp @@ -18,7 +18,7 @@ namespace papas { // max index value is 2** m_bitshift -Identifier IdCoder::makeId(unsigned int index, ItemType type, char subt, float val) { +Identifier IdCoder::makeId(uint32_t index, ItemType type, char subt, float val) { if (type == kNone) { throw "Id must have a valid type"; @@ -54,14 +54,14 @@ float IdCoder::value(Identifier id) { return bitsToFloat(bitvalue); } -unsigned int IdCoder::index(Identifier id) { return id & (uint64_t)(pow(2, m_bitshift) - 1); } +uint32_t IdCoder::index(Identifier id) { return id & (uint64_t)(pow(2, m_bitshift) - 1); } -unsigned int IdCoder::uniqueId(Identifier id) { +uint32_t IdCoder::uniqueId(Identifier id) { // For some purposes we want a smaller uniqueid without the value information // here we consruct a 32 bit uniqueid out of the index and the type and subtype - unsigned int bitshift = m_bitshift + m_bitshift1 - m_bitshift2; + uint32_t bitshift = m_bitshift + m_bitshift1 - m_bitshift2; Identifier typeShift = (uint32_t)type(id) << bitshift; - Identifier subtypeShift = (uint32_t) static_cast(tolower(subtype(id))) << m_bitshift; + Identifier subtypeShift = (uint32_t) static_cast(tolower(subtype(id))) << m_bitshift; // binary printout std::cout <<"Index" << std::bitset<32>(IdCoder::index(id)) <((uniqueid >> bitshift) & (uint32_t)(pow(2, 3) - 1)); char st = static_cast((uniqueid >> m_bitshift) & (uint64_t)(pow(2, bitshift - m_bitshift) - 1)); - unsigned idx = ((uniqueid) & (uint32_t)(pow(2, m_bitshift) - 1)); + uint32_t idx = ((uniqueid) & (uint32_t)(pow(2, m_bitshift) - 1)); if (it != type(id) || st != subtype(id) || idx != index(id)) return false; return true; } diff --git a/papas/datatypes/src/Particle.cpp b/papas/datatypes/src/Particle.cpp index 615d3fa..8f4175e 100644 --- a/papas/datatypes/src/Particle.cpp +++ b/papas/datatypes/src/Particle.cpp @@ -1,6 +1,6 @@ #include "papas/datatypes/Particle.h" -#include //lxplus needs this +#include //lxplus needs this #include #include "papas/datatypes/IdCoder.h" @@ -10,7 +10,7 @@ namespace papas { // Particle::Particle() : m_pdgId(0), m_charge(0), m_status(0) {} -Particle::Particle(int pdgid, double charge, const TLorentzVector& tlv, unsigned int index, char subtype, +Particle::Particle(int pdgid, double charge, const TLorentzVector& tlv, uint32_t index, char subtype, const TVector3& startVertex, double status) : m_tlv(tlv), m_pdgId(pdgid), diff --git a/papas/datatypes/src/Track.cpp b/papas/datatypes/src/Track.cpp index adc17e4..df021c4 100644 --- a/papas/datatypes/src/Track.cpp +++ b/papas/datatypes/src/Track.cpp @@ -1,14 +1,14 @@ #include "papas/datatypes/Track.h" +#include //lxplus needs this #include -#include #include "papas/datatypes/IdCoder.h" #include "papas/utility/StringFormatter.h" namespace papas { -Track::Track(const TVector3& p3, double charge, const std::shared_ptr path, unsigned int index, char subtype) +Track::Track(const TVector3& p3, double charge, const std::shared_ptr path, uint32_t index, char subtype) : m_id(IdCoder::makeId(index, IdCoder::ItemType::kTrack, subtype, p3.Mag())), m_p3(p3), m_charge(charge), diff --git a/papas/detectors/CMSTracker.h b/papas/detectors/CMSTracker.h index 755888c..9f9f16f 100644 --- a/papas/detectors/CMSTracker.h +++ b/papas/detectors/CMSTracker.h @@ -22,8 +22,8 @@ class CMSTracker : public Tracker { * @param[in] volume The tracker cyclinders */ CMSTracker(const VolumeCylinder& volume); - virtual double ptResolution(const Track& track) const override; ///< describes tracker resolution of momentum - virtual bool acceptance(const Track& track) const override; ///< determines if a track is detected + double ptResolution(const Track& track) const override; ///< describes tracker resolution of momentum + bool acceptance(const Track& track) const override; ///< determines if a track is detected // TODO space_resolution(self, ptc): private: }; diff --git a/papas/detectors/Calorimeter.h b/papas/detectors/Calorimeter.h index 847e1c6..3e39793 100644 --- a/papas/detectors/Calorimeter.h +++ b/papas/detectors/Calorimeter.h @@ -17,6 +17,9 @@ class Calorimeter : public DetectorElement { enum LOCATION { kBarrel = 0, kEndCap = 1 }; using DetectorElement::DetectorElement; + /**Destructor*/ + virtual ~Calorimeter() = default; // needed for classes that inherit + /** energy Resolution of ECAL @param[in] energy : energy @param[in] eta : angle of arrival diff --git a/papas/detectors/Detector.h b/papas/detectors/Detector.h index 8e8c35d..88c191f 100644 --- a/papas/detectors/Detector.h +++ b/papas/detectors/Detector.h @@ -52,8 +52,8 @@ class Detector { void setupElements(); // shared pointers allow user to have their own derived ECAL and HCAL calorimeter class // that has a fixed interface defined by Calorimeter - std::shared_ptr m_ecal; - std::shared_ptr m_hcal; + std::shared_ptr m_ecal; + std::shared_ptr m_hcal; std::shared_ptr m_tracker; std::shared_ptr m_field; diff --git a/papas/detectors/DetectorElement.h b/papas/detectors/DetectorElement.h index 12ccd1d..8649ed7 100644 --- a/papas/detectors/DetectorElement.h +++ b/papas/detectors/DetectorElement.h @@ -25,6 +25,9 @@ class DetectorElement { * @param[in] const Material& material : material */ DetectorElement(Layer layer, const VolumeCylinder& volume, const Material& material); + + /** Destructor*/ + virtual ~DetectorElement() = default; // needed for classes that inherit const VolumeCylinder& volumeCylinder() const { return m_volume; } ///< return the volume cyclinder Layer layer() const { return m_layer; }; ///< returns kEcal, kHcal etc const Material& material() const { return m_material; } diff --git a/papas/detectors/Tracker.h b/papas/detectors/Tracker.h index f765b7d..5b653a0 100644 --- a/papas/detectors/Tracker.h +++ b/papas/detectors/Tracker.h @@ -11,6 +11,10 @@ class Track; class Tracker : public DetectorElement { public: using DetectorElement::DetectorElement; + + /** Destructor*/ + virtual ~Tracker() = default; // needed for classes that inherit + /** @brief virtual function that user must define to calculate the limit of momentum that can be distinguished by * tracker * @param[in] track Track for which the limit of momentum is to be determined diff --git a/papas/detectors/VolumeCylinder.h b/papas/detectors/VolumeCylinder.h index 97b11da..10d346a 100644 --- a/papas/detectors/VolumeCylinder.h +++ b/papas/detectors/VolumeCylinder.h @@ -4,7 +4,7 @@ #include "papas/datatypes/Definitions.h" #include "papas/detectors/SurfaceCylinder.h" -//forward declaration +// forward declaration class TVector3; namespace papas { diff --git a/papas/detectors/src/CMS.cpp b/papas/detectors/src/CMS.cpp index 5912321..a59bab5 100644 --- a/papas/detectors/src/CMS.cpp +++ b/papas/detectors/src/CMS.cpp @@ -13,28 +13,28 @@ namespace papas { CMS::CMS() : Detector() { // ECAL detector Element - m_ecal = std::shared_ptr{ - new CMSECAL(VolumeCylinder(Layer::kEcal, 1.55, 2.1, 1.30, 2), - Material("CMS_ECAL", 8.9e-3, 0.275), - 1.479, // eta_crack - std::vector{0.3, 1}, // emin barrel and endcap - std::vector>{{4.22163e-02, 1.55903e-01, 7.14166e-03}, - {-2.08048e-01, 3.25097e-01, 7.34244e-03}}, - std::vector>{{1.00071, -9.04973, -2.48554}, - {9.95665e-01, -3.31774, -2.11123}})}; // barrel and endcap + m_ecal = std::make_shared( + VolumeCylinder(Layer::kEcal, 1.55, 2.1, 1.30, 2), + Material("CMS_ECAL", 8.9e-3, 0.275), + 1.479, // eta_crack + std::vector{0.3, 1}, // emin barrel and endcap + std::vector>{{4.22163e-02, 1.55903e-01, 7.14166e-03}, + {-2.08048e-01, 3.25097e-01, 7.34244e-03}}, + std::vector>{{1.00071, -9.04973, -2.48554}, + {9.95665e-01, -3.31774, -2.11123}}); // barrel and endcap // HCAL detector element - m_hcal = std::shared_ptr{ - new CMSHCAL(VolumeCylinder(Layer::kHcal, 2.9, 3.6, 1.9, 2.6), - Material("CMS_HCAL", 0.0, 0.17), - 1.3, // eta crack - std::vector>{{0.8062, 2.753, 0.1501}, {6.803e-06, 6.676, 0.1716}}, - std::vector>{{1.036, 4.452, -2.458}, {1.071, 9.471, -2.823}})}; + m_hcal = std::make_shared( + VolumeCylinder(Layer::kHcal, 2.9, 3.6, 1.9, 2.6), + Material("CMS_HCAL", 0.0, 0.17), + 1.3, // eta crack + std::vector>{{0.8062, 2.753, 0.1501}, {6.803e-06, 6.676, 0.1716}}, + std::vector>{{1.036, 4.452, -2.458}, {1.071, 9.471, -2.823}}); // Tracker detector element - m_tracker = std::shared_ptr{new CMSTracker(VolumeCylinder(Layer::kTracker, 1.29, 1.99))}; + m_tracker = std::make_shared(VolumeCylinder(Layer::kTracker, 1.29, 1.99)); // Field detector element - m_field = std::shared_ptr{new CMSField(VolumeCylinder(Layer::kField, 2.9, 3.6), 3.8)}; + m_field = std::make_shared(VolumeCylinder(Layer::kField, 2.9, 3.6), 3.8); setupElements(); // sets up a list of all detector elements (m_elements) (needed for propagation) } diff --git a/papas/display/GDetectorElement.h b/papas/display/GDetectorElement.h index 9de9f87..6db7256 100644 --- a/papas/display/GDetectorElement.h +++ b/papas/display/GDetectorElement.h @@ -12,7 +12,7 @@ namespace papas { - //Forward Declaration +// Forward Declaration class DetectorElement; /// Class to draw the HCAL and ECAL detectors diff --git a/papas/display/src/PFApp.cpp b/papas/display/src/PFApp.cpp index 1ac1a76..2cca833 100644 --- a/papas/display/src/PFApp.cpp +++ b/papas/display/src/PFApp.cpp @@ -18,7 +18,7 @@ PFApp::PFApp() {} void PFApp::display(const Event& event, const Detector& det) { m_display = new PFEventDisplay({ViewPane::Projection::xy, ViewPane::Projection::yz}); - m_gdetector = std::make_shared(GDetector(det)); + m_gdetector = std::make_shared(det); m_display->addToRegister(m_gdetector, 0); m_display->drawEvent(event); } diff --git a/papas/graphtools/DirectedAcyclicGraph.h b/papas/graphtools/DirectedAcyclicGraph.h index 08e8c9e..ad4965e 100644 --- a/papas/graphtools/DirectedAcyclicGraph.h +++ b/papas/graphtools/DirectedAcyclicGraph.h @@ -70,6 +70,7 @@ template // N is the templated Node class Visitor { public: Visitor(); /// class BFSVisitor : public Visitor { public: - BFSVisitor(); + BFSVisitor(); ///< Constructor + virtual ~BFSVisitor() = default; ///< Destructor (needed for classes that inherit) void visit(const N* node) override; ///< key to visitor pattern const Nodevector& traverseChildren(const N& node, int depth = -1) override; const Nodevector& traverseParents(const N& node, int depth = -1) override; diff --git a/papas/graphtools/src/Edge.cpp b/papas/graphtools/src/Edge.cpp index f7a2be2..641a2d4 100644 --- a/papas/graphtools/src/Edge.cpp +++ b/papas/graphtools/src/Edge.cpp @@ -1,5 +1,4 @@ #include "papas/graphtools/Edge.h" - #include #include "papas/datatypes/IdCoder.h" diff --git a/papas/reconstruction/BuildPFBlocks.h b/papas/reconstruction/BuildPFBlocks.h index 0562dbc..bcba8cc 100644 --- a/papas/reconstruction/BuildPFBlocks.h +++ b/papas/reconstruction/BuildPFBlocks.h @@ -43,7 +43,7 @@ void buildPFBlocks(const Event& event, IdCoder::SubType ecalSubtype, IdCoder::Su * @param[inout] history external collection of Nodes to which parent child relations can be added */ -void buildPFBlocks(const Ids& ids, Edges&& edges, char subtype, Blocks& blocks, Nodes& history); +void buildPFBlocks(const Ids& ids, const Edges& edges, char subtype, Blocks& blocks, Nodes& history); } // end namespace papas #endif /* BUILDPFBLOCKS_h */ diff --git a/papas/reconstruction/PFBlock.h b/papas/reconstruction/PFBlock.h index 0629f18..006856c 100644 --- a/papas/reconstruction/PFBlock.h +++ b/papas/reconstruction/PFBlock.h @@ -19,7 +19,6 @@ namespace papas { Edges m_edges : Unordered map of all the edge combinations in the block use edge(id1,id2) to find an edge - static int tempBlockCount: sequential numbering of blocks (useful for debugging/tracing etc) Usage: block = PFBlock(element_ids, edges, 'r') @@ -36,10 +35,11 @@ class PFBlock { extracted will be removed from the Edges object and will become owned by the PFBlock @param[in] subtype The subtype for the identifier of the block eg 's' for split block */ - PFBlock(const Ids& elementIds, Edges& edges, unsigned int index, + PFBlock(const Ids& elementIds, const Edges& edges, uint32_t index, char subtype = 'u'); // relevant parts of edges will be removed and become owned by PFBlock PFBlock(PFBlock&& pfblock) = default; // allow move - const Ids& elementIds(bool sort = false) const { return m_elementIds; } ///< returns vector of all ids in the block + ~PFBlock(); /// destructor + const Ids& elementIds() const { return m_elementIds; } ///< returns vector of all ids in the block /** Returns list of all edges of a given edge type that are connected to a given id. @@ -81,10 +81,9 @@ class PFBlock { PFBlock(const PFBlock& pfblock) = default; PFBlock& operator=(const PFBlock&) = default; - Identifier m_id; ///< identifier for this block - Ids m_elementIds; ///< ids of elements in this block ordered by type and decreasing energy - Edges m_edges; ///< all the edges for elements in this block - static int tempBlockCount; ///< sequential numbering of blocks, not essential but helpful for debugging + Identifier m_id; ///< identifier for this block + Ids m_elementIds; ///< ids of elements in this block ordered by type and decreasing energy + Edges m_edges; ///< all the edges for elements in this block }; std::ostream& operator<<(std::ostream& os, const PFBlock& block); diff --git a/papas/reconstruction/PFReconstructor.h b/papas/reconstruction/PFReconstructor.h index fddca23..98b6edc 100644 --- a/papas/reconstruction/PFReconstructor.h +++ b/papas/reconstruction/PFReconstructor.h @@ -90,6 +90,7 @@ class PFReconstructor { */ PFReconstructor(const Event& event, char blockSubtype, const Detector& detector, Particles& particles, Nodes& history); + ~PFReconstructor(); // const Particles& particles() const { return m_particles; } // private: diff --git a/papas/reconstruction/PapasManager.h b/papas/reconstruction/PapasManager.h index 32a0333..20bcd0d 100644 --- a/papas/reconstruction/PapasManager.h +++ b/papas/reconstruction/PapasManager.h @@ -50,6 +50,7 @@ class PapasManager { * @param[in] detector : the detector to be used in the simulation and reconstruction */ PapasManager(const Detector& detector); + /** * @brief Runs Papas simulation on a collection of papas simulated particles creating * their clusters and tracks and paths and also history information of the connections @@ -68,7 +69,7 @@ class PapasManager { * particles-simulated "ps" * history */ - void simulate(Particles& particles); + void simulate(char particleSubtype = 's'); /** * @brief Merges a set of clusters according to detector sensitivities. * @param[in] typeAndSubtype: identifies which clusters to merge, eg "hs" for hcals-smeared clusters @@ -86,25 +87,25 @@ class PapasManager { /** * @brief return Event object * @return Event */ - const Event& event() const { return m_event; } //< Access the event + const Event& event() const { return m_event; } ///< Access the event + void addParticles(const Particles& particles); ///< Add particles collection to the event const Detector& detector() const { return m_detector; } ///< Access the detector void setEventNo(unsigned int eventNo) { m_event.setEventNo(eventNo); } ///< Set the event No void clear(); /// m_ownedClustersList; // m_ownedTracksList; // m_ownedBlocksList; // m_ownedParticlesList; // m_ownedHistoryList; // m_ownedClustersList; /// m_ownedTracksList; /// m_ownedBlocksList; /// m_ownedParticlesList; /// subGraphs = buildSubGraphs(ids, edges); for (const auto& elementIds : subGraphs) { PFBlock block(elementIds, edges, blocks.size(), subtype); // make the block diff --git a/papas/reconstruction/src/MergeClusters.cpp b/papas/reconstruction/src/MergeClusters.cpp index 346e812..0657808 100644 --- a/papas/reconstruction/src/MergeClusters.cpp +++ b/papas/reconstruction/src/MergeClusters.cpp @@ -30,6 +30,16 @@ void mergeClusters(const Event& event, const std::string& typeAndSubtype, const // a new merged cluster. auto subGraphs = buildSubGraphs(ids, edges); + /* Note on debate as to safety of storing const Cluster* in overlappingClusters. + In particular the question "Could the Clusters that are being pointed to move and thus the pointers become + invalid?" + (1) According to Stroustrup (4th edition chapter 31, page 886). “When inserting and erasing elements in a vector, + elements may be moved. In contrast elements stored in a list or an associative container do not move when new + elements are inserted or other elements are erased. + (2) In addition, at the time when we store const Cluster* into the list, the unordered_map of Cluster in which the + Cluster lives is already a const object. + So this should be OK. + */ for (const auto& subgraph : subGraphs) { std::list overlappingClusters; for (const auto& cid : subgraph) { diff --git a/papas/reconstruction/src/PFBlock.cpp b/papas/reconstruction/src/PFBlock.cpp index 0f8d43a..c021472 100644 --- a/papas/reconstruction/src/PFBlock.cpp +++ b/papas/reconstruction/src/PFBlock.cpp @@ -8,8 +8,6 @@ namespace papas { -int PFBlock::tempBlockCount = 0; - bool blockIdComparer(Identifier id1, Identifier id2) { if (IdCoder::type(id1) == IdCoder::type(id2)) return id1 > id2; @@ -17,19 +15,24 @@ bool blockIdComparer(Identifier id1, Identifier id2) { return IdCoder::type(id1) < IdCoder::type(id2); } -PFBlock::PFBlock(const Ids& element_ids, Edges& edges, unsigned int index, char subtype) +PFBlock::~PFBlock() { + // //If I move ~PFBlock to the header files or if I use the default destructor I get a seg fault under Gaudi + // everything looks fine just before (the Block prints correctly) but + // when this point is reached it has id of zero and a ridiculous size for m_edges. + m_elementIds.clear(); + m_edges.clear(); +}; + +PFBlock::PFBlock(const Ids& element_ids, const Edges& edges, uint32_t index, char subtype) : m_id(IdCoder::makeId(index, IdCoder::kBlock, subtype, element_ids.size())), m_elementIds(element_ids) { - PFBlock::tempBlockCount += 1; - // extract the relevant parts of the complete set of edges and store this within the block - // note the edges will be removed from the edges unordered_map + // copy the relevant parts of the complete set of edges and store this within the block for (auto id1 : m_elementIds) { for (auto id2 : m_elementIds) { if (id1 >= id2) continue; - // move the edge from one unordered map to the other - const auto& e = edges.find(Edge::makeKey(id1, id2)); + // copy the edge from one unordered map to the other + auto e = edges.find(Edge::makeKey(id1, id2)); if (e != edges.end()) { - m_edges.emplace(e->second.key(), std::move(e->second)); - edges.erase(e); + m_edges.emplace(e->first, e->second); // I checked - this copies the edge } } } @@ -188,7 +191,9 @@ std::string PFBlock::info() const { // One liner summary of PFBlock std::ostream& operator<<(std::ostream& os, const PFBlock& block) { os << "block:" << block.info() << std::endl; os << block.elementsString(); - os << block.edgeMatrixString(); + if (block.edges().size() > 0) { + os << block.edgeMatrixString(); + } return os; } diff --git a/papas/reconstruction/src/PFReconstructor.cpp b/papas/reconstruction/src/PFReconstructor.cpp index 0f79d99..33fa2c5 100644 --- a/papas/reconstruction/src/PFReconstructor.cpp +++ b/papas/reconstruction/src/PFReconstructor.cpp @@ -32,8 +32,14 @@ PFReconstructor::PFReconstructor(const Event& event, char blockSubtype, const De PDebug::write("{},", u); // TODO warning message } + PDebug::write("Finished reconstruction"); } +PFReconstructor::~PFReconstructor() { // needed to avoid seg fault (may be connected to optimisation) + m_locked.clear(); + m_unused.clear(); +}; + void PFReconstructor::reconstructBlock(const PFBlock& block) { // see class description for summary of reconstruction approach PDebug::write("Processing {}", block); @@ -85,6 +91,7 @@ void PFReconstructor::reconstructBlock(const PFBlock& block) { m_unused.insert(id); } } + PDebug::write("Finished block", IdCoder::pretty(block.id())); } void PFReconstructor::reconstructMuons(const PFBlock& block) { @@ -104,8 +111,8 @@ void PFReconstructor::reconstructElectrons(const PFBlock& block) { Ids ids = block.elementIds(); /* the simulator does not simulate electron energy deposits in ecal. - # therefore, one should not lock the ecal clusters linked to the - # electron track as these clusters are coming from other particles.*/ + # therefore, one should not lock the ecal clusters linked to the + # electron track as these clusters are coming from other particles.*/ for (auto id : ids) { if (IdCoder::isTrack(id) && isFromParticle(id, "ps", 11)) { @@ -131,8 +138,8 @@ void PFReconstructor::insertParticle(const Ids& parentIds, Particle& newparticle bool PFReconstructor::isFromParticle(Identifier id, const std::string& typeAndSubtype, int pdgid) const { /*returns: True if object identifier comes, directly or indirectly, -from a particle of type type_and_subtype, with this absolute pdgid. -*/ + from a particle of type type_and_subtype, with this absolute pdgid. + */ auto historyHelper = HistoryHelper(m_event); auto parentIds = historyHelper.linkedIds(id, typeAndSubtype, DAG::enumVisitType::PARENTS); bool isFromPdgId = false; @@ -144,19 +151,19 @@ from a particle of type type_and_subtype, with this absolute pdgid. double PFReconstructor::neutralHadronEnergyResolution(double energy, double eta) const { /*Currently returns the hcal resolution of the detector in use. - That's a generic solution, but CMS is doing the following - (implementation in commented code) -http://cmslxr.fnal.gov/source/RecoParticleFlow/PFProducer/src/PFAlgo.cc#3350 - */ + That's a generic solution, but CMS is doing the following + (implementation in commented code) + http://cmslxr.fnal.gov/source/RecoParticleFlow/PFProducer/src/PFAlgo.cc#3350 + */ auto resolution = m_detector.hcal()->energyResolution(energy, eta); return resolution; } double PFReconstructor::nsigmaHcal(const Cluster& cluster) const { /*Currently returns 2. - CMS is doing the following (implementation in commented code) -http://cmslxr.fnal.gov/source/RecoParticleFlow/PFProducer/src/PFAlgo.cc#3365 - */ + CMS is doing the following (implementation in commented code) + http://cmslxr.fnal.gov/source/RecoParticleFlow/PFProducer/src/PFAlgo.cc#3365 + */ return 2; } @@ -283,6 +290,7 @@ void PFReconstructor::reconstructCluster(const Cluster& cluster, papas::Layer la TVector3 p3(cluster.position().Unit() * momentum); TLorentzVector p4(p3.Px(), p3.Py(), p3.Pz(), energy); // mass is not accurate here Particle particle(pdgId, 0., p4, m_particles.size(), 'r', vertex); + propagator(particle.charge())->setPath(particle); propagator(particle.charge())->propagateOne(particle, m_detector.ecal()->volumeCylinder().inner()); if (layer == papas::Layer::kHcal) { // alice not sure particle.path()->addPoint(papas::Position::kHcalIn, cluster.position()); @@ -294,7 +302,6 @@ void PFReconstructor::reconstructCluster(const Cluster& cluster, papas::Layer la // of the hcal? // nb this only is problem if the cluster and the assigned layer are different m_locked[cluster.id()] = true; // alice : just OK but not nice if hcal used to make ecal. - // TODO make more flexible and able to detect what type of cluster PDebug::write("Made {} from Merged{}", particle, cluster); insertParticle(parentIds, particle); } @@ -310,7 +317,7 @@ void PFReconstructor::reconstructTrack(const Track& track, int pdgId, const Ids& track.path()->namedPoint(papas::Position::kVertex)); //#todo fix this so it picks up smeared track points (need to propagate smeared track) - // particle.set_path(track.path) + propagator(particle.charge())->setPath(particle); m_locked[track.id()] = true; PDebug::write("Made {} from Smeared{}", particle, track); insertParticle(parentIds, particle); diff --git a/papas/reconstruction/src/PapasManager.cpp b/papas/reconstruction/src/PapasManager.cpp index 0136a6e..a3c827c 100644 --- a/papas/reconstruction/src/PapasManager.cpp +++ b/papas/reconstruction/src/PapasManager.cpp @@ -9,9 +9,11 @@ namespace papas { -PapasManager::PapasManager(const Detector& detector) : m_detector(detector), m_event() {} +PapasManager::PapasManager(const Detector& detector) : m_detector(detector), m_history(), m_event(m_history) {} -void PapasManager::simulate(Particles& particles) { +void PapasManager::addParticles(const Particles& particles) { m_event.addCollectionToFolder(particles); } + +void PapasManager::simulate(char particleSubtype) { // create empty collections that will be passed to simulator to fill // the new collection is to be a concrete class owned by the PapasManger // and stored in a list of collections. @@ -24,12 +26,10 @@ void PapasManager::simulate(Particles& particles) { auto& smearedHcalClusters = createClusters(); auto& tracks = createTracks(); auto& smearedTracks = createTracks(); - auto& history = createHistory(); - m_event.setHistory(history); // run the simulator which will fill the above objects - Simulator simulator(m_event, m_detector, ecalClusters, hcalClusters, smearedEcalClusters, smearedHcalClusters, tracks, - smearedTracks, particles, history); + Simulator simulator(m_event, particleSubtype, m_detector, ecalClusters, hcalClusters, smearedEcalClusters, + smearedHcalClusters, tracks, smearedTracks, m_history); // store the addresses of the filled collections to the Event m_event.addCollectionToFolder(ecalClusters); @@ -38,54 +38,42 @@ void PapasManager::simulate(Particles& particles) { m_event.addCollectionToFolder(smearedHcalClusters); m_event.addCollectionToFolder(tracks); m_event.addCollectionToFolder(smearedTracks); - // NB can only add the particle collection once the particles are completed (eg paths added in) - // this is because they are stored here as const objects - m_event.addCollectionToFolder(particles); - m_event.extendHistory(history); } void PapasManager::mergeClusters(const std::string& typeAndSubtype) { EventRuler ruler(m_event); // create collections ready to receive outputs auto& mergedClusters = createClusters(); - auto& history = createHistory(); - papas::mergeClusters(m_event, typeAndSubtype, ruler, mergedClusters, history); + papas::mergeClusters(m_event, typeAndSubtype, ruler, mergedClusters, m_history); // add outputs into event m_event.addCollectionToFolder(mergedClusters); - m_event.extendHistory(history); } void PapasManager::buildBlocks(const char ecalSubtype, char hcalSubtype, char trackSubtype) { // create empty collections to hold the ouputs, the ouput will be added by the algorithm auto& blocks = createBlocks(); - auto& history = createHistory(); - buildPFBlocks(m_event, ecalSubtype, hcalSubtype, trackSubtype, blocks, history); + buildPFBlocks(m_event, ecalSubtype, hcalSubtype, trackSubtype, blocks, m_history); // store a pointer to the ouputs into the event m_event.addCollectionToFolder(blocks); - m_event.extendHistory(history); } void PapasManager::simplifyBlocks(char blockSubtype) { // create empty collections to hold the ouputs, the ouput will be added by the algorithm auto& simplifiedblocks = createBlocks(); - auto& history = createHistory(); - simplifyPFBlocks(m_event, blockSubtype, simplifiedblocks, history); + simplifyPFBlocks(m_event, blockSubtype, simplifiedblocks, m_history); // store a pointer to the outputs into the event m_event.addCollectionToFolder(simplifiedblocks); - m_event.extendHistory(history); } void PapasManager::reconstruct(char blockSubtype) { - auto& history = createHistory(); auto& recParticles = createParticles(); - PFReconstructor pfReconstructor(m_event, blockSubtype, m_detector, recParticles, history); + PFReconstructor pfReconstructor(m_event, blockSubtype, m_detector, recParticles, m_history); m_event.addCollectionToFolder(recParticles); - m_event.extendHistory(history); } void PapasManager::clear() { m_event.clear(); - m_ownedHistoryList.clear(); + m_history.clear(); m_ownedClustersList.clear(); m_ownedTracksList.clear(); m_ownedBlocksList.clear(); @@ -115,8 +103,4 @@ Particles& PapasManager::createParticles() { return m_ownedParticlesList.back(); } -Nodes& PapasManager::createHistory() { - m_ownedHistoryList.emplace_back(Nodes()); - return m_ownedHistoryList.back(); -} } // end namespace papas diff --git a/papas/reconstruction/src/PapasManagerTester.cpp b/papas/reconstruction/src/PapasManagerTester.cpp index 4f65b80..bd44d2d 100644 --- a/papas/reconstruction/src/PapasManagerTester.cpp +++ b/papas/reconstruction/src/PapasManagerTester.cpp @@ -25,11 +25,10 @@ Simulator PapasManagerTester::setSimulator(Particles& particles) { auto& smearedHcalClusters = createClusters(); auto& tracks = createTracks(); auto& smearedTracks = createTracks(); - auto& history = createHistory(); // run the simulator which will fill the above objects - Simulator simulator(m_event, m_detector, ecalClusters, hcalClusters, smearedEcalClusters, smearedHcalClusters, tracks, - smearedTracks, particles, history); + Simulator simulator(m_event, 's', m_detector, ecalClusters, hcalClusters, smearedEcalClusters, smearedHcalClusters, + tracks, smearedTracks, m_history); return simulator; } diff --git a/papas/reconstruction/src/SimplifyPFBlocks.cpp b/papas/reconstruction/src/SimplifyPFBlocks.cpp index 9eaceb9..293cedf 100644 --- a/papas/reconstruction/src/SimplifyPFBlocks.cpp +++ b/papas/reconstruction/src/SimplifyPFBlocks.cpp @@ -21,16 +21,17 @@ void simplifyPFBlocks(const Event& event, char blockSubtype, Blocks& simplifiedb } void simplifyPFBlock(const Edges& toUnlink, const PFBlock& block, Blocks& simplifiedBlocks, Nodes& history) { + // take a block, unlink some of the edges and + // create smaller blocks or a simplified blocks + // or if nothing has changed take a copy of the original block if (toUnlink.size() == 0) { - // no change to this block - // make a copy of the block and put it in the simplified blocks - Edges newedges = block.edges(); // copy edges - PFBlock newblock(block.elementIds(), newedges, simplifiedBlocks.size(), 's'); + // no change needed, just make a copy of block + PFBlock newblock(block.elementIds(), block.edges(), simplifiedBlocks.size(), 's'); // will copy edges and ids PDebug::write("Made {}", newblock); - simplifiedBlocks.emplace(newblock.id(), std::move(newblock)); + auto id = newblock.id(); + simplifiedBlocks.emplace(id, std::move(newblock)); // update history - makeHistoryLinks(block.elementIds(), {newblock.id()}, history); - + makeHistoryLinks(block.elementIds(), {id}, history); } else { Edges modifiedEdges; for (auto edge : block.edges()) { // copying edges @@ -38,10 +39,10 @@ void simplifyPFBlock(const Edges& toUnlink, const PFBlock& block, Blocks& simpli if (toUnlink.find(edge.first) != toUnlink.end()) { e.setLinked(false); } - modifiedEdges.emplace(e.key(), std::move(e)); + modifiedEdges.emplace(e.key(), e); } // create new blocks and add into simplifiedBlocks - buildPFBlocks(block.elementIds(), std::move(modifiedEdges), 's', simplifiedBlocks, history); + buildPFBlocks(block.elementIds(), modifiedEdges, 's', simplifiedBlocks, history); } } @@ -49,7 +50,7 @@ Edges edgesToUnlink(const PFBlock& block) { Edges toUnlink; Ids ids = block.elementIds(); if (ids.size() > 1) { - Ids linkedIds; // std::set linkedIds; //unordered? + Ids linkedIds; bool firstHCAL; double minDist = -1; for (auto id : ids) { diff --git a/papas/simulation/HelixPropagator.h b/papas/simulation/HelixPropagator.h index 4c2db3e..2d7ca6d 100644 --- a/papas/simulation/HelixPropagator.h +++ b/papas/simulation/HelixPropagator.h @@ -21,7 +21,13 @@ class HelixPropagator : public Propagator { @param[in] cyl cylinder to which the particle is to be propagated. @param[in] field magnitude of magnetic field */ - virtual void propagateOne(Particle& ptc, const SurfaceCylinder& cyl) const override; + virtual void propagateOne(const Particle& ptc, const SurfaceCylinder& cyl) const override; + + /** Sets the particle path to a helix + @param[in] ptc particle that is to be propagated + */ + virtual void setPath(Particle& ptc) const override; }; } // end namespace papas + #endif diff --git a/papas/simulation/Propagator.h b/papas/simulation/Propagator.h index ef2b8d4..64857e7 100644 --- a/papas/simulation/Propagator.h +++ b/papas/simulation/Propagator.h @@ -17,19 +17,26 @@ class Propagator { */ Propagator(std::shared_ptr field) : m_field(field){}; + virtual ~Propagator() = default; + /** Propagate particle to the selected cylinder and store the point where the particle crossed the cylinder @param[in] ptc particle that is to be propagated @param[in] cyl cylinder to which the particle is to be propagated. @param[in] field magnitude of magnetic field (used only for charged particles, defaults to zero if not set) */ - virtual void propagateOne(Particle& ptc, const SurfaceCylinder& cyl) const = 0; + virtual void propagateOne(const Particle& ptc, const SurfaceCylinder& cyl) const = 0; /** Propagate particle all cylinders of the detector @param[in] ptc particle that is to be propagated @param[in] detector Detector through which to propagate */ - void propagate(Particle& ptc, const Detector& detector) const; + void propagate(const Particle& ptc, const Detector& detector) const; + + /** Sets the particle path to a helix or a straightline + @param[in] ptc particle that is to be propagated + */ + virtual void setPath(Particle& ptc) const = 0; protected: std::shared_ptr m_field; diff --git a/papas/simulation/Simulator.h b/papas/simulation/Simulator.h index 287cedb..7280abd 100644 --- a/papas/simulation/Simulator.h +++ b/papas/simulation/Simulator.h @@ -43,13 +43,15 @@ class Simulator { @param[in] (papas) particles to which simulation information will be added @param[in] history structure into which new history can be added, may be empty at start */ - Simulator(const Event& event, const Detector& detector, Clusters& ecalClusters, Clusters& hcalClusters, - Clusters& smearedEcalClusters, Clusters& smearedHcalClusters, Tracks& tracks, Tracks& smearedtracks, - Particles& particles, Nodes& history); + Simulator(const Event& event, char particleSubtype, const Detector& detector, Clusters& ecalClusters, + Clusters& hcalClusters, Clusters& smearedEcalClusters, Clusters& smearedHcalClusters, Tracks& tracks, + Tracks& smearedtracks, Nodes& history); + /** Simulate a particle, to produce tracks, smearedtracks, clusters, smearedclusters and path info @param[in] ptc the particle to be simulated */ - void simulateParticle(Particle& ptc); + void simulateParticle(const Particle& ptc); + /* Find the cluster with the specified identifier @param[in] clusterId the identifier of the cluster */ @@ -63,13 +65,14 @@ class Simulator { Nodes& history() { return m_history; } /// return a reference to history nodes collection // const Particles& particles() const { return m_particles; } ///< Return particles collection void clear(); ///< Clear all the collections of clusters, particles, tracks - /** - Smears a Cluster - @param[in] cluster the cluster that is to be smeared - @param[in] detectorLayer the layer to be used for smearing. Note this is not always the same as the layer to which - the cluster belongs - @return the smeared Cluster (moved) - */ + + /** + Smears a Cluster + @param[in] cluster the cluster that is to be smeared + @param[in] detectorLayer the layer to be used for smearing. Note this is not always the same as the layer to which + the cluster belongs + @return the smeared Cluster (moved) + */ Cluster smearCluster(const Cluster& cluster, papas::Layer detectorLayer = papas::Layer::kNone); /// propagator(double charge) const; - /** - Update the History - @param[in] newid the new object being added to history - @param[in] parentif the id of the parent (if 0 there is no parent) - */ - void addNode(const Identifier newid, const Identifier parentid = 0); + /** Returns the detector element associated with a particular layer @param[in] layer an enumeration describing which layer eg Layer::kEcal @@ -184,7 +190,6 @@ class Simulator { */ std::shared_ptr elem(papas::Layer layer) const; - // void testing(); const Event& m_event; ///< Event (for consistent Algorithm arguments)? const Detector& m_detector; ///< the Detector Clusters& m_ecalClusters; ///< ecal clusters (prior to smearing) @@ -193,7 +198,6 @@ class Simulator { Clusters& m_smearedHcalClusters; ///< smeared hcal clusters Tracks& m_tracks; ///< tracks Tracks& m_smearedTracks; ///< smeared tracks - Particles& m_particles; ///< all particles Nodes& m_history; ///< Records relationships of everything that is simulated std::shared_ptr m_propStraight; /// field) : Propagator(field) {} -void HelixPropagator::propagateOne(Particle& ptc, const SurfaceCylinder& cyl) const { +void HelixPropagator::setPath(Particle& ptc) const { if (ptc.path() == nullptr) { - auto helix = std::make_shared(Helix(ptc.p4(), ptc.startVertex(), ptc.charge(), m_field->getMagnitude())); + auto helix = std::make_shared(ptc.p4(), ptc.startVertex(), ptc.charge(), m_field->getMagnitude()); ptc.setPath(helix); } - auto helix = std::static_pointer_cast(ptc.path()); +} +void HelixPropagator::propagateOne(const Particle& ptc, const SurfaceCylinder& cyl) const { + auto helix = std::static_pointer_cast(ptc.path()); bool is_looper = helix->extremePointXY().Mag() < cyl.radius(); double udir_z = helix->unitDirection().Z(); @@ -34,9 +36,7 @@ void HelixPropagator::propagateOne(Particle& ptc, const SurfaceCylinder& cyl) co helix->addPoint(cyl.layer(), destination); } else is_looper = true; - } - - catch (std::string s) { + } catch (std::string s) { return; } } diff --git a/papas/simulation/src/Propagator.cpp b/papas/simulation/src/Propagator.cpp index cc24361..dcb411e 100644 --- a/papas/simulation/src/Propagator.cpp +++ b/papas/simulation/src/Propagator.cpp @@ -4,7 +4,7 @@ #include "papas/detectors/Field.h" namespace papas { -void Propagator::propagate(Particle& ptc, const Detector& detector) const { +void Propagator::propagate(const Particle& ptc, const Detector& detector) const { for (const auto el : detector.elements()) propagateOne(ptc, el->volumeCylinder().inner()); } diff --git a/papas/simulation/src/Simulator.cpp b/papas/simulation/src/Simulator.cpp index 0e84eaf..6995b20 100644 --- a/papas/simulation/src/Simulator.cpp +++ b/papas/simulation/src/Simulator.cpp @@ -1,5 +1,6 @@ #include "papas/simulation/Simulator.h" +#include "papas/datatypes/Event.h" #include "papas/datatypes/Helix.h" #include "papas/datatypes/IdCoder.h" #include "papas/datatypes/ParticlePData.h" @@ -14,9 +15,9 @@ namespace papas { -Simulator::Simulator(const Event& papasevent, const Detector& detector, Clusters& ecalClusters, Clusters& hcalClusters, - Clusters& smearedEcalClusters, Clusters& smearedHcalClusters, Tracks& tracks, - Tracks& smearedTracks, Particles& particles, Nodes& history) +Simulator::Simulator(const Event& papasevent, char particleSubtype, const Detector& detector, Clusters& ecalClusters, + Clusters& hcalClusters, Clusters& smearedEcalClusters, Clusters& smearedHcalClusters, + Tracks& tracks, Tracks& smearedTracks, Nodes& history) : m_event(papasevent), m_detector(detector), m_ecalClusters(ecalClusters), @@ -25,7 +26,6 @@ Simulator::Simulator(const Event& papasevent, const Detector& detector, Clusters m_smearedHcalClusters(smearedHcalClusters), m_tracks(tracks), m_smearedTracks(smearedTracks), - m_particles(particles), m_history(history) { @@ -33,14 +33,15 @@ Simulator::Simulator(const Event& papasevent, const Detector& detector, Clusters m_propStraight = std::make_shared(detector.field()); // make sure we can process the particles in order if needed (highest energy first) Ids ids; + auto particles = papasevent.particles(particleSubtype); for (auto& p : particles) ids.insert(p.first); for (auto& id : ids) { - simulateParticle(particles[id]); + simulateParticle(particles.at(id)); } } -void Simulator::simulateParticle(Particle& ptc) { +void Simulator::simulateParticle(const Particle& ptc) { int pdgid = ptc.pdgId(); if (ptc.charge() && ptc.pt() < 0.2 && abs(pdgid) >= 100) { // to avoid numerical problems in propagation @@ -61,7 +62,7 @@ void Simulator::simulateParticle(Particle& ptc) { } } -void Simulator::simulatePhoton(Particle& ptc) { +void Simulator::simulatePhoton(const Particle& ptc) { PDebug::write("Simulating Photon"); // find where the photon meets the Ecal inner cylinder // make and smear the cluster @@ -73,7 +74,7 @@ void Simulator::simulatePhoton(Particle& ptc) { } } -void Simulator::simulateHadron(Particle& ptc) { +void Simulator::simulateHadron(const Particle& ptc) { PDebug::write("Simulating Hadron"); auto ecal_sp = m_detector.ecal(); auto hcal_sp = m_detector.hcal(); @@ -124,12 +125,12 @@ void Simulator::simulateHadron(Particle& ptc) { } } -void Simulator::simulateNeutrino(Particle& ptc) { +void Simulator::simulateNeutrino(const Particle& ptc) { PDebug::write("Simulating Neutrino \n"); propagator(ptc.charge())->propagate(ptc, m_detector); } -void Simulator::simulateElectron(Particle& ptc) { +void Simulator::simulateElectron(const Particle& ptc) { /*Simulate an electron corresponding to gen particle ptc. Uses the methods detector.electronEnergyResolution @@ -148,7 +149,7 @@ void Simulator::simulateElectron(Particle& ptc) { } } -void Simulator::simulateMuon(Particle& ptc) { +void Simulator::simulateMuon(const Particle& ptc) { /*Simulate a muon corresponding to gen particle ptc Uses the methods detector.muon_energy_resolution @@ -212,7 +213,7 @@ Cluster Simulator::makeAndStoreEcalCluster(const Particle& ptc, double fraction, } Cluster cluster(energy, pos, csize, m_ecalClusters.size(), IdCoder::kEcalCluster, subtype); Identifier id = cluster.id(); - addNode(id, ptc.id()); + makeHistoryLink(ptc.id(), id, m_history); PDebug::write("Made {}", cluster); m_ecalClusters.emplace(id, std::move(cluster)); return m_ecalClusters[id]; @@ -235,7 +236,7 @@ Cluster Simulator::makeAndStoreHcalCluster(const Particle& ptc, double fraction, } Cluster cluster(energy, pos, csize, m_hcalClusters.size(), IdCoder::kHcalCluster, subtype); Identifier id = cluster.id(); - addNode(id, ptc.id()); + makeHistoryLink(ptc.id(), id, m_history); PDebug::write("Made {}", cluster); m_hcalClusters.emplace(id, std::move(cluster)); return m_hcalClusters[id]; @@ -259,7 +260,7 @@ Cluster Simulator::smearCluster(const Cluster& parent, papas::Layer detectorLaye double energyresolution = sp_calorimeter->energyResolution(parent.energy(), parent.eta()); double response = sp_calorimeter->energyResponse(parent.energy(), parent.eta()); double energy = parent.energy() * rootrandom::Random::gauss(response, energyresolution); - unsigned int counter; + uint32_t counter; if (IdCoder::layer(parent.id()) == Layer::kEcal) counter = m_smearedEcalClusters.size(); else @@ -283,14 +284,14 @@ bool Simulator::acceptSmearedCluster(const Cluster& smearedCluster, papas::Layer const Cluster& Simulator::storeSmearedEcalCluster(Cluster&& smearedCluster, Identifier parentId) { auto id = smearedCluster.id(); - addNode(id, parentId); + makeHistoryLink(parentId, id, m_history); m_smearedEcalClusters.emplace(id, std::move(smearedCluster)); return m_smearedEcalClusters[id]; } const Cluster& Simulator::storeSmearedHcalCluster(Cluster&& smearedCluster, Identifier parentId) { auto id = smearedCluster.id(); - addNode(id, parentId); + makeHistoryLink(parentId, id, m_history); m_smearedHcalClusters.emplace(id, std::move(smearedCluster)); return m_smearedHcalClusters[id]; } @@ -300,14 +301,14 @@ const Track& Simulator::makeAndStoreTrack(const Particle& ptc) { Identifier id = track.id(); PDebug::write("Made {}", track); m_tracks.emplace(id, std::move(track)); - addNode(id, ptc.id()); + makeHistoryLink(ptc.id(), id, m_history); return m_tracks.at(id); } -void Simulator::storeSmearedTrack(Track&& track, Identifier parentid) { +void Simulator::storeSmearedTrack(Track&& track, Identifier parentId) { Identifier id = track.id(); m_smearedTracks.emplace(id, std::move(track)); - addNode(id, parentid); + makeHistoryLink(parentId, id, m_history); } Track Simulator::smearTrack(const Track& track, double resolution) const { @@ -347,17 +348,6 @@ bool Simulator::acceptMuonSmearedTrack(const Track& smearedTrack, bool accept) c } } -void Simulator::addNode(Identifier newid, const Identifier parentid) { - // add the new node into the set of all nodes - PFNode node{newid}; - m_history.emplace(newid, std::move(node)); - if (parentid) { - PFNode& parent = m_history[parentid]; - PFNode& child = m_history[newid]; - parent.addChild(child); - } -} - void Simulator::clear() { m_ecalClusters.clear(); m_hcalClusters.clear(); @@ -365,7 +355,6 @@ void Simulator::clear() { m_smearedHcalClusters.clear(); m_tracks.clear(); m_smearedTracks.clear(); - m_particles.clear(); } std::shared_ptr Simulator::elem(papas::Layer layer) const { return m_detector.element(layer); } diff --git a/papas/simulation/src/StraightLinePropagator.cpp b/papas/simulation/src/StraightLinePropagator.cpp index 48e6f11..04acc7a 100644 --- a/papas/simulation/src/StraightLinePropagator.cpp +++ b/papas/simulation/src/StraightLinePropagator.cpp @@ -11,15 +11,18 @@ namespace papas { StraightLinePropagator::StraightLinePropagator(std::shared_ptr field) : Propagator(field) {} -void StraightLinePropagator::propagateOne(Particle& ptc, const SurfaceCylinder& cyl) const { +void StraightLinePropagator::setPath(Particle& ptc) const { + if (ptc.path() == nullptr) { + auto line = std::make_shared(ptc.p4(), ptc.startVertex(), ptc.charge()); + ptc.setPath(line); + } +} + +void StraightLinePropagator::propagateOne(const Particle& ptc, const SurfaceCylinder& cyl) const { auto layer = cyl.layer(); double cylinderz = cyl.z(); double cylinderradius = cyl.radius(); std::shared_ptr line = ptc.path(); - if (line == nullptr) { - line = std::make_shared(Path(ptc.p4(), ptc.startVertex(), ptc.charge())); - ptc.setPath(line); - } auto udir = line->unitDirection(); auto origin = line->origin(); double theta = udir.Theta(); diff --git a/papaslib/CMakeLists.txt b/papaslib/CMakeLists.txt index 90f0c3b..d9b3ffa 100644 --- a/papaslib/CMakeLists.txt +++ b/papaslib/CMakeLists.txt @@ -13,15 +13,12 @@ SET(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib") # which point to directories outside the build tree to the install RPATH SET(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) - # the RPATH to be used when installing, but only if it's not a system directory LIST(FIND CMAKE_PLATFORM_IMPLICIT_LINK_DIRECTORIES "${CMAKE_INSTALL_PREFIX}/lib" isSystemDir) IF("${isSystemDir}" STREQUAL "-1") SET(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib") ENDIF("${isSystemDir}" STREQUAL "-1") - - include_directories( ${CMAKE_SOURCE_DIR} ) @@ -29,7 +26,6 @@ ${CMAKE_SOURCE_DIR} SET(CMAKE_INSTALL_PATH "${CMAKE_SOURCE_DIR}/install") Message(STATUS,"install" ${CMAKE_INSTALL_PATH} ) - file(GLOB_RECURSE sources ${CMAKE_SOURCE_DIR}/papas/datatypes/*.cpp @@ -50,11 +46,13 @@ file(GLOB_RECURSE headers ${CMAKE_SOURCE_DIR}/papas/reconstruction/*.h ${CMAKE_SOURCE_DIR}/papas/simulation/*.h ${CMAKE_SOURCE_DIR}/papas/utility/*.h + ${CMAKE_SOURCE_DIR}/spdlog/*.h + ${CMAKE_SOURCE_DIR}/spdlog/sinks/*.h + ${CMAKE_SOURCE_DIR}/spdlog/details/*.h ) -message("${CMAKE_SOURCE_DIR}") #message("Headers: ${headers}") -#message("Headers: ${sources}") +#message("Sources: ${sources}") FILE(GLOB files "${CMAKE_SOURCE_DIR}/papas/datatypes/*.h") INSTALL(FILES ${files} DESTINATION include/papas/datatypes/) @@ -70,12 +68,14 @@ FILE(GLOB files6 "${CMAKE_SOURCE_DIR}/papas/display/*.h") INSTALL(FILES ${files6} DESTINATION include/papas/display/) FILE(GLOB files7 "${CMAKE_SOURCE_DIR}/papas/utility/*.h") INSTALL(FILES ${files7} DESTINATION include/papas/utility/) -FILE(GLOB files8 "${CMAKE_SOURCE_DIR}/spdlog/*.h") -INSTALL(FILES ${files8} DESTINATION include/papas/spdlog/) +FILE(GLOB files8 "${CMAKE_SOURCE_DIR}/spdlog/*.*") +INSTALL(FILES ${files8} DESTINATION include/spdlog/) FILE(GLOB files9 "${CMAKE_SOURCE_DIR}/spdlog/details/*.h") -INSTALL(FILES ${files9} DESTINATION include/papas/spdlog/details/) -FILE(GLOB files10 "${CMAKE_SOURCE_DIR}/spdlog/sink/*.h") -INSTALL(FILES ${files10} DESTINATION include/papas/spdlog/sink/) +INSTALL(FILES ${files9} DESTINATION include/spdlog/details/) +FILE(GLOB files11 "${CMAKE_SOURCE_DIR}/spdlog/details/*.cc") +INSTALL(FILES ${files11} DESTINATION include/spdlog/details/) +FILE(GLOB files10 "${CMAKE_SOURCE_DIR}/spdlog/sinks/*.*") +INSTALL(FILES ${files10} DESTINATION include/spdlog/sinks/) project(papaslib) @@ -87,7 +87,7 @@ ${CMAKE_SOURCE_DIR} # remove ### to get back root dicts ###REFLEX_GENERATE_DICTIONARY(papas ${headers} SELECTION selection.xml) -message(STATUS, "bin" ${PROJECT_BINARY_DIR}) +#message(STATUS, "bin" ${PROJECT_BINARY_DIR}) ###add_library(papasDict SHARED ${sources} ${headers} papas.cxx) ###add_dependencies(papasDict papas-dictgen ) diff --git a/tests/unittest.cpp b/tests/unittest.cpp index 8bacb67..98afdcb 100644 --- a/tests/unittest.cpp +++ b/tests/unittest.cpp @@ -4,10 +4,10 @@ // C++ #include +#include #include #include #include -#include #include // ROOT @@ -15,24 +15,24 @@ #include "TLorentzVector.h" #include "TVector3.h" -#include "papas/display/Display.h" -#include "papas/display/GTrajectory.h" -#include "papas/display/ViewPane.h" #include "papas/datatypes/Event.h" #include "papas/datatypes/Helix.h" #include "papas/datatypes/HistoryHelper.h" #include "papas/detectors/CMS.h" #include "papas/detectors/CMSField.h" #include "papas/detectors/Calorimeter.h" -#include "papas/simulation/StraightLinePropagator.h" -#include "papas/simulation/HelixPropagator.h" +#include "papas/display/Display.h" +#include "papas/display/GTrajectory.h" +#include "papas/display/ViewPane.h" #include "papas/graphtools/Distance.h" #include "papas/graphtools/EventRuler.h" #include "papas/reconstruction/BuildPFBlocks.h" #include "papas/reconstruction/MergeClusters.h" #include "papas/reconstruction/PapasManagerTester.h" #include "papas/reconstruction/SimplifyPFBlocks.h" +#include "papas/simulation/HelixPropagator.h" #include "papas/simulation/Simulator.h" +#include "papas/simulation/StraightLinePropagator.h" #include "papas/utility/TRandom.h" using namespace papas; @@ -88,9 +88,11 @@ TEST_CASE("Helixpath") { /// Helix path test Particle particle(211, -1, TLorentzVector{2., 0, 1, 5}, 1, 'r', TVector3{0, 0, 0}, 3.8); HelixPropagator helixprop(field); //(particle.p4(), {0,0,0}, 3.8, -1); + helixprop.setPath(particle); helixprop.propagateOne(particle, cyl1); const auto& tvec = particle.path()->namedPoint(cyl1.layer()); auto particle2 = Particle(211, -1, TLorentzVector{0., 2, 1, 5}, 2, 'r', TVector3{0, 0, 0}, 3.8); + helixprop.setPath(particle2); helixprop.propagateOne(particle2, cyl1); const auto& tvec2 = particle2.path()->namedPoint(cyl1.layer()); REQUIRE(fabs(tvec.X()) == Approx(fabs(tvec2.Y()))); @@ -199,7 +201,7 @@ TEST_CASE("StraightLine") { TLorentzVector tlv{1, 0, 1, 2.}; Particle photon(22, 0, tlv, 0, 't'); - + propStraight.setPath(photon); propStraight.propagateOne(photon, cyl1); propStraight.propagateOne(photon, cyl2); auto points = photon.path()->points(); @@ -215,6 +217,7 @@ TEST_CASE("StraightLine") { tlv = TLorentzVector(1, 0, -1, 2.); Particle photon2(22, 0, tlv, 1, 't'); + propStraight.setPath(photon2); propStraight.propagateOne(photon2, cyl1); propStraight.propagateOne(photon2, cyl2); points = photon2.path()->points(); @@ -227,6 +230,7 @@ TEST_CASE("StraightLine") { // extrapolating from a vertex close to +endcap tlv = TLorentzVector(1, 0, 1, 2.); Particle photon3(22, 0, tlv, 3, 's', {0, 0, 1.5}, 0.); + propStraight.setPath(photon3); propStraight.propagateOne(photon3, cyl1); points = photon3.path()->points(); REQUIRE(points[papas::Position::kEcalIn].Perp() == Approx(.5)); @@ -234,6 +238,7 @@ TEST_CASE("StraightLine") { // extrapolating from a vertex close to -endcap tlv = TLorentzVector(1, 0, -1, 2.); Particle photon4(22, 0, tlv, 4, 's', {0, 0, -1.5}, 0.); + propStraight.setPath(photon4); propStraight.propagateOne(photon4, cyl1); points = photon4.path()->points(); REQUIRE(points[papas::Position::kEcalIn].Perp() == Approx(.5)); @@ -245,6 +250,7 @@ TEST_CASE("StraightLine") { 0., 0.5, 0, }, 0.); + propStraight.setPath(photon5); propStraight.propagateOne(photon5, cyl1); points = photon5.path()->points(); REQUIRE(points[papas::Position::kEcalIn].Perp() == Approx(1.)); @@ -448,7 +454,7 @@ TEST_CASE("BlockSplitter") { Edges to_unlink; to_unlink[edge1.key()] = edge1; - Event event; + Event event(emptyNodes); event.addCollectionToFolder(blocks); Blocks simplifiedBlocks; simplifyPFBlocks(event, 'r', simplifiedBlocks, emptyNodes); @@ -464,7 +470,7 @@ TEST_CASE("Merge") { eclusters.emplace(cluster2.id(), cluster2); Clusters mergedClusters; Nodes nodes; - Event event; + Event event(nodes); event.addCollectionToFolder(eclusters); papas::EventRuler ruler(event); @@ -491,7 +497,7 @@ TEST_CASE("merge_pair") { hclusters.emplace(cluster2.id(), cluster2); Clusters mergedClusters; Nodes nodes; - Event event; + Event event(nodes); event.addCollectionToFolder(hclusters); papas::EventRuler ruler(event); @@ -511,7 +517,7 @@ TEST_CASE("merge_pair_away") { Clusters mergedClusters; Nodes nodes; - Event event; + Event event(nodes); event.addCollectionToFolder(hclusters); papas::EventRuler ruler(event); @@ -533,15 +539,15 @@ TEST_CASE("merge_different_layers") { Clusters mergedClusters; Nodes nodes; - Event event; + Event event(nodes); REQUIRE_THROWS(event.addCollectionToFolder(hclusters)); return; } TEST_CASE("test_papasevent") { - - Event event; + Nodes nodes; + Event event(nodes); Clusters ecals; Tracks tracks; Identifier lastid = 0; @@ -573,8 +579,9 @@ TEST_CASE("test_papasevent") { } TEST_CASE("test_history") { - Event event; + Nodes history; + Event event(history); Clusters ecals; Particles particles; Identifier lastid = 0; @@ -620,7 +627,7 @@ TEST_CASE("merge_inside") { Tracks tracks; Nodes nodes; - Event testevent; + Event testevent(nodes); Clusters mergedClusters; testevent.addCollectionToFolder(hclusters); papas::EventRuler ruler(testevent);