From 57aabb0985e3fba9e9cc78c4603f27907addc057 Mon Sep 17 00:00:00 2001 From: Angira Rastogi Date: Fri, 2 Aug 2024 11:32:52 -0700 Subject: [PATCH 01/14] Updating CMakeLists to work with TrkHitStudiesWorkspace package --- CMakeLists.txt | 32 ++------------------------------ 1 file changed, 2 insertions(+), 30 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 7d8c667..20f5f9d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,7 @@ # cmake file for building MarlinTrkProcessors # - based on $MARLIN/example/mymarlin/CMakeLists.txt # by @author Jan Engels, Desy IT -CMAKE_MINIMUM_REQUIRED(VERSION 2.6 FATAL_ERROR) +CMAKE_MINIMUM_REQUIRED(VERSION 3.15 FATAL_ERROR) ######################################################## @@ -20,23 +20,6 @@ cmake_policy(SET CMP0008 NEW) FIND_PACKAGE( ILCUTIL REQUIRED COMPONENTS ILCSOFT_CMAKE_MODULES ) -# load default settings from ILCSOFT_CMAKE_MODULES -INCLUDE( ilcsoft_default_settings ) - -SET( COMPILER_FLAGS "-Wno-effc++" ) -MESSAGE( STATUS "FLAGS ${COMPILER_FLAGS}" ) -FOREACH( FLAG ${COMPILER_FLAGS} ) - CHECK_CXX_COMPILER_FLAG( "${FLAG}" CXX_FLAG_WORKS_${FLAG} ) - IF( ${CXX_FLAG_WORKS_${FLAG}} ) - MESSAGE ( STATUS "Adding ${FLAG} to CXX_FLAGS" ) - ### We append the flag, so they overwrite things from somewhere else - SET ( CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${FLAG} ") - ELSE() - MESSAGE ( STATUS "NOT Adding ${FLAG} to CXX_FLAGS" ) - ENDIF() -ENDFOREACH() - - FIND_PACKAGE( Marlin 1.0 REQUIRED ) # minimum required Marlin version INCLUDE_DIRECTORIES( SYSTEM ${Marlin_INCLUDE_DIRS} ) LINK_LIBRARIES( ${Marlin_LIBRARIES} ) @@ -128,16 +111,5 @@ INCLUDE_DIRECTORIES( . ./source/Digitisers/include ./source/Refitting/include ./ # add library FILE( GLOB_RECURSE library_sources "source/*.cc" ) -# needed for adding header files to xcode project -IF(CMAKE_GENERATOR MATCHES "Xcode") - FILE( GLOB_RECURSE library_headers "*.h" ) - ADD_SHARED_LIBRARY( ${PROJECT_NAME} ${library_sources} ${library_headers}) -ELSE() - ADD_SHARED_LIBRARY( ${PROJECT_NAME} ${library_sources} ) -ENDIF() - -INSTALL_SHARED_LIBRARY( ${PROJECT_NAME} DESTINATION lib ) - -# display some variables and write them to cache -DISPLAY_STD_VARIABLES() +ADD_LIBRARY( ${PROJECT_NAME} SHARED ${library_sources} ) From 7f098af93c3409bf9d113db9e359612d1c3c73ed Mon Sep 17 00:00:00 2001 From: Angira Rastogi Date: Tue, 22 Oct 2024 16:49:22 -0700 Subject: [PATCH 02/14] Improving the cluster filters application --- source/Utils/include/FilterClusters.h | 6 +- source/Utils/src/FilterClusters.cc | 91 ++++++++++++++++++++------- 2 files changed, 72 insertions(+), 25 deletions(-) diff --git a/source/Utils/include/FilterClusters.h b/source/Utils/include/FilterClusters.h index ab76767..4d37f6d 100644 --- a/source/Utils/include/FilterClusters.h +++ b/source/Utils/include/FilterClusters.h @@ -45,12 +45,14 @@ class FilterClusters : public marlin::Processor std::string _OutRelationCollection {}; //! Ranges for theta - std::vector _ThetaRanges; + std::vector _ThetaRanges; //! Cut-offs for cluster size in various theta ranges - std::vector _ClusterSize; + std::vector _ClusterSize; //! Layers to be filtered std::vector _Layers; + //!Number of bins in theta + std::string _ThetaBins; }; diff --git a/source/Utils/src/FilterClusters.cc b/source/Utils/src/FilterClusters.cc index 1797855..541d5e9 100644 --- a/source/Utils/src/FilterClusters.cc +++ b/source/Utils/src/FilterClusters.cc @@ -31,12 +31,17 @@ FilterClusters::FilterClusters() registerProcessorParameter("ThetaRanges", "Divide theta into bins for different cluster size cuts", _ThetaRanges, - {} + {} + ); + registerProcessorParameter("ThetaBins", + "Number of bins in theta", + _ThetaBins, + {} ); registerProcessorParameter("ClusterSize", "Maximum cluster size for each theta range", _ClusterSize, - {} + {} ); registerProcessorParameter("Layers", "Layers to be filtered", @@ -102,7 +107,7 @@ void FilterClusters::processEvent( LCEvent * evt ) streamlog_out(DEBUG0) << "Wrong collection type for InRelationCollection. \n"; - streamlog_out(DEBUG0) << "Number of Elements in VB Tracker Hits Collection: " << InTrackerHitCollection->getNumberOfElements() <getNumberOfElements() <getNumberOfElements(); ++i) //loop through all hits { @@ -122,21 +127,33 @@ void FilterClusters::processEvent( LCEvent * evt ) //Calculating cluster size const lcio::LCObjectVec &rawHits = trkhit->getRawHits(); - float max = -1000000; - float min = 1000000; + float ymax = -1000000; + float xmax = -1000000; + float ymin = 1000000; + float xmin = 1000000; for (size_t j=0; j( rawHits[j] ); const double *localPos = hitConstituent->getPosition(); - x = localPos[0]; - y = localPos[1]; - if (y < min){ - min = y; - } - else if (y > max){ - max = y; - } + float x_local = localPos[0]; + float y_local = localPos[1]; + + if (y_local < ymin){ + ymin = y_local; + } + if (y_local > ymax){ + ymax = y_local; + } + + if (x_local < xmin){ + xmin = x_local; + } + if (x_local > xmax){ + xmax = x_local; + } } - float cluster_size = (max - min)+1; + float cluster_size_y = (ymax - ymin)+1; + float cluster_size_x = (xmax - xmin)+1; + float cluster_size = rawHits.size(); //Get hit subdetector/layer std::string _encoderString = lcio::LCTrackerCellID::encoding_string(); @@ -144,24 +161,52 @@ void FilterClusters::processEvent( LCEvent * evt ) uint32_t systemID = decoder(trkhit)["system"]; uint32_t layerID = decoder(trkhit)["layer"]; bool filter_layer = false; + + int rows = _Layers.size(), cols = std::stoi(_ThetaBins); + if((rows*cols != _ThetaRanges.size()) || (rows*(cols-1) != _ClusterSize.size())){ + std::cout<<"Either theta cuts or cluster cuts not provided for each layer. Please change the config, exiting now..."<> _thetaCuts_byLayer; + std::vector> _clusterSizeCuts_byLayer; + + for (int i = 0; i < rows; ++i) { + std::vector row; + for (int j = 0; j < cols; ++j) { + row.push_back(std::stof(_ThetaRanges[j])); + } + _thetaCuts_byLayer.push_back(row); + } + + for (int i = 0; i < rows; ++i) { + std::vector row; + for (int j = 0; j < cols; ++j) { + row.push_back(std::stof(_ClusterSize[j])); + } + _clusterSizeCuts_byLayer.push_back(row); + } + + for (size_t j=0; j<_Layers.size(); ++j){ if (layerID == std::stof(_Layers[j])) { filter_layer = true; + break; } } streamlog_out(DEBUG0) << "Filter layer: " << filter_layer << std::endl; - - for (size_t j=0; j<_ThetaRanges.size()-1; ++j) { + + for (size_t j=0; j<_thetaCuts_byLayer[layerID].size()-1; ++j) { streamlog_out( DEBUG0 ) << "theta: " << incidentTheta << std::endl; - float min_theta = std::stof(_ThetaRanges[j]); - float max_theta = std::stof(_ThetaRanges[j+1]); + float min_theta = _thetaCuts_byLayer[layerID][j]; + float max_theta = _thetaCuts_byLayer[layerID][j+1]; streamlog_out( DEBUG0 ) << "theta range: " << min_theta << ", " << max_theta << std::endl; - - if(incidentTheta > min and incidentTheta <= max and not filter_layer){ + + if(incidentTheta >= min_theta and incidentTheta <= max_theta and filter_layer){ streamlog_out( DEBUG0 ) << "theta in range" << std::endl; - streamlog_out( DEBUG0 ) << "cluster size cut off: " << _ClusterSize[j] << std::endl; + streamlog_out( DEBUG0 ) << "cluster size cut off: " << _clusterSizeCuts_byLayer[layerID][j] << std::endl; streamlog_out( DEBUG0 ) << "cluster size: " << cluster_size << std::endl; - if(cluster_size < std::stof(_ClusterSize[j])) { + if(cluster_size < _clusterSizeCuts_byLayer[layerID][j]) { streamlog_out( DEBUG0 ) << "cluster added" << std::endl; OutTrackerHitCollection->addElement(trkhit); OutRelationCollection->addElement(rel); } @@ -170,7 +215,7 @@ void FilterClusters::processEvent( LCEvent * evt ) } } else{ - streamlog_out( DEBUG0 ) << "theta out of range or layer filtered" << std::endl; + streamlog_out( DEBUG0 ) << "theta out of range or filtering not enabled for this layer" << std::endl; } } } From 307dccec104ffa01c0e1bdc1d3ce233adf704572 Mon Sep 17 00:00:00 2001 From: Angira Rastogi Date: Wed, 23 Oct 2024 23:54:39 -0700 Subject: [PATCH 03/14] Adding pixel hits per cluster --- source/Utils/src/FilterTimeHits.cc | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/source/Utils/src/FilterTimeHits.cc b/source/Utils/src/FilterTimeHits.cc index 1419668..60e8995 100644 --- a/source/Utils/src/FilterTimeHits.cc +++ b/source/Utils/src/FilterTimeHits.cc @@ -271,6 +271,12 @@ void FilterTimeHits::processEvent(LCEvent *evt) hit_new->setTime(hit->getTime()); hit_new->setQuality(hit->getQuality()); + const lcio::LCObjectVec &rawHits = hit->getRawHits(); + for (size_t j=0; j( rawHits[j] ); + hit_new->rawHits().push_back(hitConstituent); + } + outputTrackerHitColls[icol]->addElement(hit_new); LCRelation *rel = dynamic_cast(inputHitRels[icol]->getElementAt(ihit)); From 5bdae5e30bb874a0a6664ea9c8e70cbbbede5999 Mon Sep 17 00:00:00 2001 From: Angira Rastogi Date: Fri, 25 Oct 2024 13:11:09 -0700 Subject: [PATCH 04/14] Storing reco-sim hit relations of filtered clusters, adding debugging histograms --- source/Utils/include/FilterClusters.h | 20 ++- source/Utils/src/FilterClusters.cc | 168 ++++++++++++++++++++------ 2 files changed, 148 insertions(+), 40 deletions(-) diff --git a/source/Utils/include/FilterClusters.h b/source/Utils/include/FilterClusters.h index 4d37f6d..7d6af94 100644 --- a/source/Utils/include/FilterClusters.h +++ b/source/Utils/include/FilterClusters.h @@ -1,6 +1,7 @@ #pragma once #include +#include //#include @@ -39,10 +40,12 @@ class FilterClusters : public marlin::Processor //! Input track collection std::string _InTrackerHitCollection {}; std::string _InRelationCollection {}; + std::string _InSimTrackerHitCollection {}; //! Output track collection std::string _OutTrackerHitCollection {}; std::string _OutRelationCollection {}; + std::string _OutSimTrackerHitCollection {}; //! Ranges for theta std::vector _ThetaRanges; @@ -50,9 +53,20 @@ class FilterClusters : public marlin::Processor //! Cut-offs for cluster size in various theta ranges std::vector _ClusterSize; - //! Layers to be filtered - std::vector _Layers; - + //! Layers to be filtered + std::vector _Layers; + //!Number of bins in theta std::string _ThetaBins; + + bool m_fillHistos{}; + + // --- Diagnostic histograms: + TH1F *m_clusterTheta_beforeCut = nullptr; + TH1F *m_clusterTheta_afterCut = nullptr; + TH1F *m_clusterLayer_beforeCut = nullptr; + TH1F *m_clusterLayer_afterCut = nullptr; + TH1F *m_clusterSize_beforeCut = nullptr; + TH1F *m_clusterSize_afterCut = nullptr; + }; diff --git a/source/Utils/src/FilterClusters.cc b/source/Utils/src/FilterClusters.cc index 541d5e9..d0ccd40 100644 --- a/source/Utils/src/FilterClusters.cc +++ b/source/Utils/src/FilterClusters.cc @@ -6,7 +6,9 @@ #include #include +#include #include +#include #include #include #include @@ -29,59 +31,89 @@ FilterClusters::FilterClusters() // register steering parameters: name, description, class-variable, default value registerProcessorParameter("ThetaRanges", - "Divide theta into bins for different cluster size cuts", - _ThetaRanges, - {} - ); + "Divide theta into bins for different cluster size cuts", + _ThetaRanges, + {} + ); registerProcessorParameter("ThetaBins", - "Number of bins in theta", - _ThetaBins, - {} - ); + "Number of bins in theta", + _ThetaBins, + {} + ); registerProcessorParameter("ClusterSize", - "Maximum cluster size for each theta range", - _ClusterSize, + "Maximum cluster size for each theta range", + _ClusterSize, {} - ); + ); registerProcessorParameter("Layers", - "Layers to be filtered", - _Layers, - {} - ); + "Layers to be filtered", + _Layers, + {} + ); + registerInputCollection( LCIO::SIMTRACKERHIT, + "InSimTrackerHitCollection" , + "Name of the input sim tracker hit collection", + _InSimTrackerHitCollection, + _InSimTrackerHitCollection + ); + registerInputCollection( LCIO::TRACKERHIT, - "InTrackerHitCollection" , + "InTrackerHitCollection" , "Name of the input tracker hit collection", _InTrackerHitCollection, - _InTrackerHitCollection - ); + _InTrackerHitCollection + ); registerInputCollection( LCIO::LCRELATION, - "InRelationCollection" , + "InRelationCollection" , "Name of the input relation collection", _InRelationCollection, - _InRelationCollection - ); + _InRelationCollection + ); + + registerOutputCollection( LCIO::SIMTRACKERHIT, + "OutSimTrackerHitCollection" , + "Name of output sim tracker hit collection", + _OutSimTrackerHitCollection, + _OutSimTrackerHitCollection + ); registerOutputCollection( LCIO::TRACKERHIT, - "OutTrackerHitCollection" , - "Name of output tracker hit collection", - _OutTrackerHitCollection, - std::string("FilteredVBTrackerHits") + "OutTrackerHitCollection" , + "Name of output tracker hit collection", + _OutTrackerHitCollection, + _OutTrackerHitCollection ); - registerOutputCollection( LCIO::LCRELATION, - "OutRelationCollection" , - "Name of output relation collection", - _OutRelationCollection, - std::string("FilteredVBTrackerHitsRelations") + registerOutputCollection( LCIO::LCRELATION, + "OutRelationCollection" , + "Name of output relation collection", + _OutRelationCollection, + _OutRelationCollection ); + registerProcessorParameter("FillHistograms", + "Flag to fill the diagnostic histograms", + m_fillHistos, + false); + } void FilterClusters::init() { // Print the initial parameters printParameters() ; + + AIDA::ITree* tree=marlin::AIDAProcessor::tree(this); + marlin::AIDAProcessor::histogramFactory(this); + + m_clusterTheta_beforeCut = new TH1F("m_ClusterTheta_before", "Cluster Theta [radian]", 32, 0., 3.2); + m_clusterTheta_afterCut = new TH1F("m_ClusterTheta_after", "Cluster Theta [radian]", 32, 0., 3.2); + m_clusterLayer_beforeCut = new TH1F("m_ClusterLayer_before", "Cluster Layer Index", 10, 0, 10); + m_clusterLayer_afterCut = new TH1F("m_ClusterLayer_after", "Cluster Layer Index", 10, 0, 10); + m_clusterSize_beforeCut = new TH1F("m_ClusterSize_before", "Cluster hit multiplicity", 20, 0, 20); + m_clusterSize_afterCut = new TH1F("m_ClusterSize_after", "Cluster hit multiplicity", 20, 0, 20); + } void FilterClusters::processRunHeader( LCRunHeader* /*run*/) @@ -94,10 +126,13 @@ void FilterClusters::processEvent( LCEvent * evt ) OutTrackerHitCollection->setSubset(true); LCCollectionVec *OutRelationCollection = new LCCollectionVec(LCIO::LCRELATION); OutRelationCollection->setSubset(true); + LCCollectionVec *OutSimTrackerHitCollection = new LCCollectionVec(LCIO::SIMTRACKERHIT); + OutSimTrackerHitCollection->setSubset(true); // Get input collection LCCollection* InTrackerHitCollection = evt->getCollection(_InTrackerHitCollection); LCCollection* InRelationCollection = evt->getCollection(_InRelationCollection); + LCCollection* InSimTrackerHitCollection = evt->getCollection(_InSimTrackerHitCollection); if( InTrackerHitCollection->getTypeName() != lcio::LCIO::TRACKERHITPLANE ) { throw EVENT::Exception( "Invalid collection type: " + InTrackerHitCollection->getTypeName() ) ; } @@ -105,6 +140,9 @@ void FilterClusters::processEvent( LCEvent * evt ) if( InRelationCollection->getTypeName() != lcio::LCIO::LCRELATION ) { throw EVENT::Exception( "Invalid collection type: " + InRelationCollection->getTypeName() ) ; } streamlog_out(DEBUG0) << "Wrong collection type for InRelationCollection. \n"; + if( InSimTrackerHitCollection->getTypeName() != lcio::LCIO::SIMTRACKERHIT ) + { throw EVENT::Exception( "Invalid collection type: " + InSimTrackerHitCollection->getTypeName() ) ; } + streamlog_out(DEBUG0) << "Wrong collection type for SimTrackerHitCollection. \n"; streamlog_out(DEBUG0) << "Number of Elements in Tracker Hits Collection: " << InTrackerHitCollection->getNumberOfElements() <getNumberOfElements(); ++i) //loop through all hits { streamlog_out(DEBUG0) << "Loop over hits opened. \n"; - EVENT::TrackerHit *trkhit=static_cast(InTrackerHitCollection->getElementAt(i)); //define trkhit var, pointer to i'th element of tracker hits - EVENT::LCRelation *rel=static_cast(InRelationCollection->getElementAt(i)); + TrackerHitPlane *trkhit=static_cast(InTrackerHitCollection->getElementAt(i)); //define trkhit var, pointer to i'th element of tracker hits //Calculating theta float x = trkhit->getPosition()[0]; @@ -171,7 +208,7 @@ void FilterClusters::processEvent( LCEvent * evt ) std::vector> _thetaCuts_byLayer; std::vector> _clusterSizeCuts_byLayer; - for (int i = 0; i < rows; ++i) { + for (int k = 0; k < rows; ++k) { std::vector row; for (int j = 0; j < cols; ++j) { row.push_back(std::stof(_ThetaRanges[j])); @@ -179,7 +216,7 @@ void FilterClusters::processEvent( LCEvent * evt ) _thetaCuts_byLayer.push_back(row); } - for (int i = 0; i < rows; ++i) { + for (int k = 0; k < rows; ++k) { std::vector row; for (int j = 0; j < cols; ++j) { row.push_back(std::stof(_ClusterSize[j])); @@ -196,6 +233,12 @@ void FilterClusters::processEvent( LCEvent * evt ) } streamlog_out(DEBUG0) << "Filter layer: " << filter_layer << std::endl; + if (m_fillHistos){ + m_clusterTheta_beforeCut->Fill(incidentTheta); + m_clusterLayer_beforeCut->Fill(layerID); + m_clusterSize_beforeCut->Fill(cluster_size); + } + for (size_t j=0; j<_thetaCuts_byLayer[layerID].size()-1; ++j) { streamlog_out( DEBUG0 ) << "theta: " << incidentTheta << std::endl; float min_theta = _thetaCuts_byLayer[layerID][j]; @@ -208,8 +251,55 @@ void FilterClusters::processEvent( LCEvent * evt ) streamlog_out( DEBUG0 ) << "cluster size: " << cluster_size << std::endl; if(cluster_size < _clusterSizeCuts_byLayer[layerID][j]) { streamlog_out( DEBUG0 ) << "cluster added" << std::endl; - OutTrackerHitCollection->addElement(trkhit); - OutRelationCollection->addElement(rel); } + OutTrackerHitCollection->addElement(trkhit); + + if (m_fillHistos){ + m_clusterTheta_afterCut->Fill(incidentTheta); + m_clusterLayer_afterCut->Fill(layerID); + m_clusterSize_afterCut->Fill(cluster_size); + } + + EVENT::LCRelation *rel=static_cast(InRelationCollection->getElementAt(i)); + SimTrackerHit *simhit = dynamic_cast(rel->getTo()); + SimTrackerHitImpl *simhit_new = new SimTrackerHitImpl(); + simhit_new->setCellID0(simhit->getCellID0()); + simhit_new->setCellID1(simhit->getCellID1()); + simhit_new->setPosition(simhit->getPosition()); + simhit_new->setEDep(simhit->getEDep()); + simhit_new->setTime(simhit->getTime()); + simhit_new->setMCParticle(simhit->getMCParticle()); + simhit_new->setMomentum(simhit->getMomentum()); + simhit_new->setPathLength(simhit->getPathLength()); + simhit_new->setQuality(simhit->getQuality()); + simhit_new->setOverlay(simhit->isOverlay()); + simhit_new->setProducedBySecondary(simhit->isProducedBySecondary()); + OutSimTrackerHitCollection->addElement(simhit_new); + + TrackerHitPlaneImpl *hit_new = new TrackerHitPlaneImpl(); + hit_new->setCellID0(trkhit->getCellID0()); + hit_new->setCellID1(trkhit->getCellID1()); + hit_new->setType(trkhit->getType()); + hit_new->setPosition(trkhit->getPosition()); + hit_new->setU(trkhit->getU()); + hit_new->setV(trkhit->getV()); + hit_new->setdU(trkhit->getdU()); + hit_new->setdV(trkhit->getdV()); + hit_new->setEDep(trkhit->getEDep()); + hit_new->setEDepError(trkhit->getEDepError()); + hit_new->setTime(trkhit->getTime()); + hit_new->setQuality(trkhit->getQuality()); + const lcio::LCObjectVec &rawHits_new = trkhit->getRawHits(); + for (size_t k=0; k( rawHits_new[k] ); + hit_new->rawHits().push_back(hitConstituent); + } + + LCRelationImpl *rel_new = new LCRelationImpl(); + rel_new->setFrom(hit_new); + rel_new->setTo(simhit_new); + rel_new->setWeight(rel->getWeight()); + OutRelationCollection->addElement(rel_new); + } else { streamlog_out( DEBUG0 ) << "cluster rejected" << std::endl; } @@ -222,7 +312,11 @@ void FilterClusters::processEvent( LCEvent * evt ) // Save output track collection evt->addCollection(OutTrackerHitCollection, _OutTrackerHitCollection); - evt->addCollection(OutRelationCollection, _OutRelationCollection); + evt->addCollection(OutRelationCollection, _OutRelationCollection); + evt->addCollection(OutSimTrackerHitCollection, _OutSimTrackerHitCollection); + + streamlog_out(MESSAGE) << " " << OutTrackerHitCollection->size() << " reco clusters added to the collection: " << _OutTrackerHitCollection << ", " << OutSimTrackerHitCollection->size() << " sim hits added to the collection: "<< _OutSimTrackerHitCollection << " and " << OutRelationCollection->size() << " relations added to the collection: " << _OutRelationCollection << std::endl; + } void FilterClusters::end() From 5a171c8d92d9f184150e2a6b5ca4e248473cf5ba Mon Sep 17 00:00:00 2001 From: Angira Rastogi Date: Fri, 25 Oct 2024 13:11:58 -0700 Subject: [PATCH 05/14] Adding pixel hits to the clusters from FilterTimeHits processor --- source/Utils/src/FilterTimeHits.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/source/Utils/src/FilterTimeHits.cc b/source/Utils/src/FilterTimeHits.cc index 60e8995..61d1efe 100644 --- a/source/Utils/src/FilterTimeHits.cc +++ b/source/Utils/src/FilterTimeHits.cc @@ -276,7 +276,6 @@ void FilterTimeHits::processEvent(LCEvent *evt) lcio::SimTrackerHit *hitConstituent = dynamic_cast( rawHits[j] ); hit_new->rawHits().push_back(hitConstituent); } - outputTrackerHitColls[icol]->addElement(hit_new); LCRelation *rel = dynamic_cast(inputHitRels[icol]->getElementAt(ihit)); From a0aa5bde1547245fbea458e8c2add7ec97b24289 Mon Sep 17 00:00:00 2001 From: Simone Pagan Griso Date: Fri, 25 Oct 2024 15:03:41 -0700 Subject: [PATCH 06/14] Fix status of output collection and adding new cluster to output collection --- source/Utils/src/FilterClusters.cc | 35 +++++++++++++++++++++--------- 1 file changed, 25 insertions(+), 10 deletions(-) diff --git a/source/Utils/src/FilterClusters.cc b/source/Utils/src/FilterClusters.cc index d0ccd40..59cc8a3 100644 --- a/source/Utils/src/FilterClusters.cc +++ b/source/Utils/src/FilterClusters.cc @@ -13,6 +13,7 @@ #include #include #include +#include #include @@ -121,14 +122,6 @@ void FilterClusters::processRunHeader( LCRunHeader* /*run*/) void FilterClusters::processEvent( LCEvent * evt ) { - // Make the output track collection - LCCollectionVec *OutTrackerHitCollection = new LCCollectionVec(LCIO::TRACKERHIT); - OutTrackerHitCollection->setSubset(true); - LCCollectionVec *OutRelationCollection = new LCCollectionVec(LCIO::LCRELATION); - OutRelationCollection->setSubset(true); - LCCollectionVec *OutSimTrackerHitCollection = new LCCollectionVec(LCIO::SIMTRACKERHIT); - OutSimTrackerHitCollection->setSubset(true); - // Get input collection LCCollection* InTrackerHitCollection = evt->getCollection(_InTrackerHitCollection); LCCollection* InRelationCollection = evt->getCollection(_InRelationCollection); @@ -146,6 +139,28 @@ void FilterClusters::processEvent( LCEvent * evt ) streamlog_out(DEBUG0) << "Number of Elements in Tracker Hits Collection: " << InTrackerHitCollection->getNumberOfElements() <getParameters().getStringVal("CellIDEncoding"); + LCCollectionVec *OutTrackerHitCollection = new LCCollectionVec(LCIO::TRACKERHIT); + OutTrackerHitCollection->parameters().setValue("CellIDEncoding", encoderString); + LCFlagImpl lcFlag(InTrackerHitCollection->getFlag()); + OutTrackerHitCollection->setFlag(lcFlag.getFlag()); + //OutTrackerHitCollection->setSubset(true); + + LCCollectionVec *OutSimTrackerHitCollection = new LCCollectionVec(LCIO::SIMTRACKERHIT); + OutSimTrackerHitCollection->parameters().setValue("CellIDEncoding", encoderString); + LCFlagImpl lcFlag_sim(InSimTrackerHitCollection->getFlag()); + OutSimTrackerHitCollection->setFlag(lcFlag_sim.getFlag()); + //OutSimTrackerHitCollection->setSubset(true); + + LCCollectionVec *OutRelationCollection = new LCCollectionVec(LCIO::LCRELATION); + LCFlagImpl lcFlag_rel(InRelationCollection->getFlag()); + OutRelationCollection->setFlag(lcFlag_rel.getFlag()); + //OutRelationCollection->setSubset(true); + + // Filter for(int i=0; igetNumberOfElements(); ++i) //loop through all hits { @@ -250,8 +265,7 @@ void FilterClusters::processEvent( LCEvent * evt ) streamlog_out( DEBUG0 ) << "cluster size cut off: " << _clusterSizeCuts_byLayer[layerID][j] << std::endl; streamlog_out( DEBUG0 ) << "cluster size: " << cluster_size << std::endl; if(cluster_size < _clusterSizeCuts_byLayer[layerID][j]) { - streamlog_out( DEBUG0 ) << "cluster added" << std::endl; - OutTrackerHitCollection->addElement(trkhit); + streamlog_out( DEBUG0 ) << "Adding reco/sim clusters and relation to output collections" << std::endl; if (m_fillHistos){ m_clusterTheta_afterCut->Fill(incidentTheta); @@ -293,6 +307,7 @@ void FilterClusters::processEvent( LCEvent * evt ) lcio::SimTrackerHit *hitConstituent = dynamic_cast( rawHits_new[k] ); hit_new->rawHits().push_back(hitConstituent); } + OutTrackerHitCollection->addElement(hit_new); LCRelationImpl *rel_new = new LCRelationImpl(); rel_new->setFrom(hit_new); From 154a3702a2ea2518b09d888b895330b420f370aa Mon Sep 17 00:00:00 2001 From: Angira Rastogi Date: Fri, 25 Oct 2024 18:23:39 -0700 Subject: [PATCH 07/14] Minor cleanup --- source/Utils/src/FilterClusters.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/source/Utils/src/FilterClusters.cc b/source/Utils/src/FilterClusters.cc index 59cc8a3..6b05151 100644 --- a/source/Utils/src/FilterClusters.cc +++ b/source/Utils/src/FilterClusters.cc @@ -105,7 +105,6 @@ void FilterClusters::init() // Print the initial parameters printParameters() ; - AIDA::ITree* tree=marlin::AIDAProcessor::tree(this); marlin::AIDAProcessor::histogramFactory(this); m_clusterTheta_beforeCut = new TH1F("m_ClusterTheta_before", "Cluster Theta [radian]", 32, 0., 3.2); From 308a1692e4a49e4b5be27580724d2e48ff1d1c24 Mon Sep 17 00:00:00 2001 From: Angira Rastogi Date: Mon, 28 Oct 2024 11:27:11 -0700 Subject: [PATCH 08/14] Adding before and after debugging histograms for hit timing --- source/Utils/include/FilterTimeHits.h | 3 ++- source/Utils/src/FilterTimeHits.cc | 9 +++++++-- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/source/Utils/include/FilterTimeHits.h b/source/Utils/include/FilterTimeHits.h index d98140a..170e37f 100644 --- a/source/Utils/include/FilterTimeHits.h +++ b/source/Utils/include/FilterTimeHits.h @@ -69,7 +69,8 @@ class FilterTimeHits : public Processor double m_time_max{}; // --- Diagnostic histograms: - TH1F *m_corrected_time = nullptr; + TH1F *m_corrected_time_before = nullptr; + TH1F *m_corrected_time_after = nullptr; // --- Run and event counters: int _nRun{}; diff --git a/source/Utils/src/FilterTimeHits.cc b/source/Utils/src/FilterTimeHits.cc index 61d1efe..27278fd 100644 --- a/source/Utils/src/FilterTimeHits.cc +++ b/source/Utils/src/FilterTimeHits.cc @@ -96,7 +96,9 @@ void FilterTimeHits::init() AIDAProcessor::histogramFactory(this); - m_corrected_time = new TH1F("m_corrected_time", "Corrected time of the hit [ps]", 1000, -250., 250.); + m_corrected_time_before = new TH1F("m_corrected_time_before", "Corrected time of the hit before filter [ns]", 1000, -25., 25.); + m_corrected_time_after = new TH1F("m_corrected_time_after", "Corrected time of the hit after filter [ns]", 1000, -25., 25.); + } void FilterTimeHits::processRunHeader(LCRunHeader *run) @@ -230,6 +232,9 @@ void FilterTimeHits::processEvent(LCEvent *evt) hitT -= dt; streamlog_out(DEBUG3) << "corrected hit at R: " << hitR << " mm by propagation time: " << dt << " ns to T: " << hitT << " ns" << std::endl; + if (m_fillHistos) + m_corrected_time_before->Fill(hitT); + //Apply time window selection if (hitT < m_time_min || hitT > m_time_max) { @@ -240,7 +245,7 @@ void FilterTimeHits::processEvent(LCEvent *evt) hits_to_save[icol].insert(ihit); if (m_fillHistos) - m_corrected_time->Fill(hitT); + m_corrected_time_after->Fill(hitT); } } // ihit loop From 5b20a704f798df5858a5349d9a0fcc17bc5c08b6 Mon Sep 17 00:00:00 2001 From: Angira Rastogi Date: Thu, 31 Oct 2024 01:25:35 -0700 Subject: [PATCH 09/14] Fixing the brackets for error message --- source/Utils/src/FilterClusters.cc | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/source/Utils/src/FilterClusters.cc b/source/Utils/src/FilterClusters.cc index 6b05151..f6953c8 100644 --- a/source/Utils/src/FilterClusters.cc +++ b/source/Utils/src/FilterClusters.cc @@ -127,14 +127,14 @@ void FilterClusters::processEvent( LCEvent * evt ) LCCollection* InSimTrackerHitCollection = evt->getCollection(_InSimTrackerHitCollection); if( InTrackerHitCollection->getTypeName() != lcio::LCIO::TRACKERHITPLANE ) - { throw EVENT::Exception( "Invalid collection type: " + InTrackerHitCollection->getTypeName() ) ; } - streamlog_out(DEBUG0) << "Wrong collection type for TrackerHitCollection. \n"; + { throw EVENT::Exception( "Invalid collection type: " + InTrackerHitCollection->getTypeName() ) ; + streamlog_out(DEBUG0) << "Wrong collection type for TrackerHitCollection. \n";} if( InRelationCollection->getTypeName() != lcio::LCIO::LCRELATION ) - { throw EVENT::Exception( "Invalid collection type: " + InRelationCollection->getTypeName() ) ; } - streamlog_out(DEBUG0) << "Wrong collection type for InRelationCollection. \n"; + { throw EVENT::Exception( "Invalid collection type: " + InRelationCollection->getTypeName() ) ; + streamlog_out(DEBUG0) << "Wrong collection type for InRelationCollection. \n";} if( InSimTrackerHitCollection->getTypeName() != lcio::LCIO::SIMTRACKERHIT ) - { throw EVENT::Exception( "Invalid collection type: " + InSimTrackerHitCollection->getTypeName() ) ; } - streamlog_out(DEBUG0) << "Wrong collection type for SimTrackerHitCollection. \n"; + { throw EVENT::Exception( "Invalid collection type: " + InSimTrackerHitCollection->getTypeName() ) ; + streamlog_out(DEBUG0) << "Wrong collection type for SimTrackerHitCollection. \n";} streamlog_out(DEBUG0) << "Number of Elements in Tracker Hits Collection: " << InTrackerHitCollection->getNumberOfElements() < Date: Fri, 1 Nov 2024 15:39:18 -0700 Subject: [PATCH 10/14] restructured TimeHits filter to not rely on 1-1 reco-sim hits and also keep reco-hits constituents. --- source/Utils/include/FilterTimeHits.h | 17 +- source/Utils/src/FilterTimeHits.cc | 303 +++++++++++++++++--------- 2 files changed, 217 insertions(+), 103 deletions(-) diff --git a/source/Utils/include/FilterTimeHits.h b/source/Utils/include/FilterTimeHits.h index 170e37f..871d8d6 100644 --- a/source/Utils/include/FilterTimeHits.h +++ b/source/Utils/include/FilterTimeHits.h @@ -6,6 +6,9 @@ #include #include +#include +#include + #include "DDRec/Surface.h" #include "DDRec/SurfaceManager.h" @@ -55,10 +58,12 @@ class FilterTimeHits : public Processor protected: // --- Input/output collection names: - std::vector m_inputTrackerHitsCollNames{}; - std::vector m_inputTrackerSimHitsCollNames{}; - std::vector m_inputTrackerHitRelNames{}; + std::vector m_inputTrackerHitsCollNames{}; /// Reconstructed clusters + std::vector m_inputTrackerHitsConstituentsCollNames{}; // Individual hits in clusters (if available) + std::vector m_inputTrackerSimHitsCollNames{}; // Simulated energy deposits (if available) + std::vector m_inputTrackerHitRelNames{}; // Relations between reconstructed clusters and simulated ones (if available) std::vector m_outputTrackerHitsCollNames{}; + std::vector m_outputTrackerHitsConstituentsCollNames{}; std::vector m_outputTrackerSimHitsCollNames{}; std::vector m_outputTrackerHitRelNames{}; @@ -72,6 +77,12 @@ class FilterTimeHits : public Processor TH1F *m_corrected_time_before = nullptr; TH1F *m_corrected_time_after = nullptr; + // ---- Utility functions + // Make a full copy of a TrackerHitPlane object (no copy-constructor defined) + TrackerHitPlane* copyTrackerHitPlane(TrackerHitPlane *hit); + // Make a fully copy of a SimTrackerHit object + SimTrackerHit* copySimTrackerHit(SimTrackerHit *hit); + // --- Run and event counters: int _nRun{}; int _nEvt{}; diff --git a/source/Utils/src/FilterTimeHits.cc b/source/Utils/src/FilterTimeHits.cc index 27278fd..ba49937 100644 --- a/source/Utils/src/FilterTimeHits.cc +++ b/source/Utils/src/FilterTimeHits.cc @@ -98,7 +98,6 @@ void FilterTimeHits::init() m_corrected_time_before = new TH1F("m_corrected_time_before", "Corrected time of the hit before filter [ns]", 1000, -25., 25.); m_corrected_time_after = new TH1F("m_corrected_time_after", "Corrected time of the hit after filter [ns]", 1000, -25., 25.); - } void FilterTimeHits::processRunHeader(LCRunHeader *run) @@ -106,6 +105,72 @@ void FilterTimeHits::processRunHeader(LCRunHeader *run) _nRun++; } +TrackerHitPlane *FilterTimeHits::copyTrackerHitPlane(TrackerHitPlane *hit) +{ + // Create new object + TrackerHitPlaneImpl *hit_new = new TrackerHitPlaneImpl(); + + // Make a copy (no copy-constructor) + hit_new->setCellID0(hit->getCellID0()); + hit_new->setCellID1(hit->getCellID1()); + hit_new->setType(hit->getType()); + hit_new->setPosition(hit->getPosition()); + hit_new->setU(hit->getU()); + hit_new->setV(hit->getV()); + hit_new->setdU(hit->getdU()); + hit_new->setdV(hit->getdV()); + hit_new->setEDep(hit->getEDep()); + hit_new->setEDepError(hit->getEDepError()); + hit_new->setTime(hit->getTime()); + hit_new->setQuality(hit->getQuality()); + const_cast(hit_new->getCovMatrix()) = hit->getCovMatrix(); // no setter for covariance matrix? + + // need to clone individual hits, if present + const lcio::LCObjectVec &rawHits = hit->getRawHits(); + for (size_t j = 0; j < rawHits.size(); ++j) + { + // Use (default) copy-constructor of SimTrackerHitImpl + lcio::SimTrackerHit *hitConstituent = copySimTrackerHit(dynamic_cast(rawHits[j])); + if (!hitConstituent) + { + static bool first_error = true; + if (first_error) + { + streamlog_out(WARNING) << "Cannot access individual hits of reco ID clusters. Skipping (no further WARNING will be printed.)" << std::endl; + continue; + } + } + hit_new->rawHits().push_back(hitConstituent); + } + + // return new object + return hit_new; +} + +SimTrackerHit *FilterTimeHits::copySimTrackerHit(SimTrackerHit *hit) +{ + // use (default) copy-constructor + lcio::SimTrackerHitImpl *simhit_new = new lcio::SimTrackerHitImpl(*(dynamic_cast(hit))); + + /* + // manual copy, deprecated + SimTrackerHitImpl *simhit_new = new SimTrackerHitImpl(); + simhit_new->setCellID0(simhit->getCellID0()); + simhit_new->setCellID1(simhit->getCellID1()); + simhit_new->setPosition(simhit->getPosition()); + simhit_new->setEDep(simhit->getEDep()); + simhit_new->setTime(simhit->getTime()); + simhit_new->setMCParticle(simhit->getMCParticle()); + simhit_new->setMomentum(simhit->getMomentum()); + simhit_new->setPathLength(simhit->getPathLength()); + simhit_new->setQuality(simhit->getQuality()); + simhit_new->setOverlay(simhit->isOverlay()); + simhit_new->setProducedBySecondary(simhit->isProducedBySecondary()); + */ + + return simhit_new; +} + void FilterTimeHits::processEvent(LCEvent *evt) { @@ -142,10 +207,12 @@ void FilterTimeHits::processEvent(LCEvent *evt) const unsigned int nTrackerHitCol = m_inputTrackerHitsCollNames.size(); std::vector inputHitColls(nTrackerHitCol); + std::vector inputHitConstituentsColls(nTrackerHitCol); std::vector inputSimHitColls(nTrackerHitCol); std::vector inputHitRels(nTrackerHitCol); std::vector outputTrackerHitColls(nTrackerHitCol); + std::vector outputTrackerHitConstituentsColls(nTrackerHitCol); std::vector outputTrackerSimHitColls(nTrackerHitCol); std::vector outputTrackerHitRels(nTrackerHitCol); @@ -161,31 +228,52 @@ void FilterTimeHits::processEvent(LCEvent *evt) { streamlog_out(WARNING) << m_inputTrackerHitsCollNames[icol] << " collection not available" << std::endl; - continue; + continue; // this is mandatory to do anything + } + + // get the reco hits contituents (if available) + try + { + if (m_inputTrackerHitsConstituentsCollNames[icol] != "") + inputHitConstituentsColls[icol] = evt->getCollection(m_inputTrackerHitsConstituentsCollNames[icol]); + else + inputHitConstituentsColls[icol] = nullptr; + } + catch (lcio::DataNotAvailableException &e) + { + streamlog_out(WARNING) << m_inputTrackerHitsConstituentsCollNames[icol] + << " collection not available. This is expected in simplified digitization settings." << std::endl; + inputHitConstituentsColls[icol] = nullptr; // optional collection } // get the sim hits try { - inputSimHitColls[icol] = evt->getCollection(m_inputTrackerSimHitsCollNames[icol]); + if (m_inputTrackerSimHitsCollNames[icol] != "") + inputSimHitColls[icol] = evt->getCollection(m_inputTrackerSimHitsCollNames[icol]); + else + inputSimHitColls[icol] = nullptr; } catch (lcio::DataNotAvailableException &e) { streamlog_out(WARNING) << m_inputTrackerSimHitsCollNames[icol] << " collection not available" << std::endl; - continue; + inputSimHitColls[icol] = nullptr; // optional collection } // get the reco-sim relations try { - inputHitRels[icol] = evt->getCollection(m_inputTrackerHitRelNames[icol]); + if (m_inputTrackerHitRelNames[icol] != "") + inputHitRels[icol] = evt->getCollection(m_inputTrackerHitRelNames[icol]); + else + inputHitRels[icol] = nullptr; } catch (lcio::DataNotAvailableException &e) { streamlog_out(WARNING) << m_inputTrackerHitRelNames[icol] << " collection not available" << std::endl; - continue; + inputHitRels[icol] = nullptr; // optional collection } // reco hit output collections @@ -195,129 +283,140 @@ void FilterTimeHits::processEvent(LCEvent *evt) LCFlagImpl lcFlag(inputHitColls[icol]->getFlag()); outputTrackerHitColls[icol]->setFlag(lcFlag.getFlag()); + // reco hit contituents output collections, if needed + if ((m_outputTrackerHitsConstituentsCollNames[icol] != "") && + (inputHitConstituentsColls[icol] != nullptr)) { + std::string encoderString = inputHitConstituentsColls[icol]->getParameters().getStringVal("CellIDEncoding"); + outputTrackerHitConstituentsColls[icol] = new LCCollectionVec(inputHitConstituentsColls[icol]->getTypeName()); + outputTrackerHitConstituentsColls[icol]->parameters().setValue("CellIDEncoding", encoderString); + LCFlagImpl lcFlag(inputHitConstituentsColls[icol]->getFlag()); + outputTrackerHitConstituentsColls[icol]->setFlag(lcFlag.getFlag()); + } else { + outputTrackerHitConstituentsColls[icol] = nullptr; + } + // sim hit output collections - outputTrackerSimHitColls[icol] = new LCCollectionVec(inputSimHitColls[icol]->getTypeName()); - outputTrackerSimHitColls[icol]->parameters().setValue("CellIDEncoding", encoderString); - LCFlagImpl lcFlag_sim(inputSimHitColls[icol]->getFlag()); - outputTrackerSimHitColls[icol]->setFlag(lcFlag_sim.getFlag()); + if (inputSimHitColls[icol] != nullptr) { + outputTrackerSimHitColls[icol] = new LCCollectionVec(inputSimHitColls[icol]->getTypeName()); + outputTrackerSimHitColls[icol]->parameters().setValue("CellIDEncoding", encoderString); + LCFlagImpl lcFlag_sim(inputSimHitColls[icol]->getFlag()); + outputTrackerSimHitColls[icol]->setFlag(lcFlag_sim.getFlag()); + } // reco-sim relation output collections - outputTrackerHitRels[icol] = new LCCollectionVec(inputHitRels[icol]->getTypeName()); - LCFlagImpl lcFlag_rel(inputHitRels[icol]->getFlag()); - outputTrackerHitRels[icol]->setFlag(lcFlag_rel.getFlag()); + if (inputHitRels[icol] != nullptr) { + outputTrackerHitRels[icol] = new LCCollectionVec(inputHitRels[icol]->getTypeName()); + LCFlagImpl lcFlag_rel(inputHitRels[icol]->getFlag()); + outputTrackerHitRels[icol]->setFlag(lcFlag_rel.getFlag()); + } } // --- Loop over the tracker hits and select hits inside the chosen time window: - std::vector> hits_to_save(nTrackerHitCol); - + // IMPORTANT: Cannot assume SimTrkHit and Relation collections contain corresponding elements in the same order. Some sim hits might not have a corrsponding reco hit. for (unsigned int icol = 0; icol < inputHitColls.size(); ++icol) { + // Keep track of (old, new) object pointers for clusters to save in the output collections + std::map reco_clusters_to_save; LCCollection *hit_col = inputHitColls[icol]; if (!hit_col) + { + streamlog_out(WARNING) << "Cannot retrieve collection: " << m_inputTrackerHitsCollNames[icol] << std::endl; continue; + } for (int ihit = 0; ihit < hit_col->getNumberOfElements(); ++ihit) { - if (TrackerHitPlane *hit = dynamic_cast(hit_col->getElementAt(ihit))){ - // Skipping the hit if its time is outside the acceptance time window - double hitT = hit->getTime(); - - dd4hep::rec::Vector3D pos = hit->getPosition(); - double hitR = pos.r(); - - // Correcting for the propagation time - double dt = hitR / (TMath::C() * m_beta / 1e6); - hitT -= dt; - streamlog_out(DEBUG3) << "corrected hit at R: " << hitR << " mm by propagation time: " << dt << " ns to T: " << hitT << " ns" << std::endl; - - if (m_fillHistos) - m_corrected_time_before->Fill(hitT); - - //Apply time window selection - if (hitT < m_time_min || hitT > m_time_max) - { - streamlog_out(DEBUG4) << "hit at T: " << hitT << " ns is rejected by timing cuts" << std::endl; - continue; - } - - hits_to_save[icol].insert(ihit); - - if (m_fillHistos) - m_corrected_time_after->Fill(hitT); + TrackerHitPlane *hit = dynamic_cast(hit_col->getElementAt(ihit)); + if (!hit) + { + streamlog_out(WARNING) << "Cannot retrieve/cast(TrackerHitPlane*) cluster from collection: " << m_inputTrackerHitsCollNames[icol] << std::endl; + continue; } + // Skipping the hit if its time is outside the acceptance time window + double hitT = hit->getTime(); - } // ihit loop + dd4hep::rec::Vector3D pos = hit->getPosition(); + double hitR = pos.r(); - } // icol loop + // Correcting for the propagation time + double dt = hitR / (TMath::C() * m_beta / 1e6); + hitT -= dt; + streamlog_out(DEBUG3) << "corrected hit at R: " << hitR << " mm by propagation time: " << dt << " ns to T: " << hitT << " ns" << std::endl; - // --- Add the filtered hits to the output collections: + if (m_fillHistos) + m_corrected_time_before->Fill(hitT); - for (unsigned int icol = 0; icol < inputHitColls.size(); ++icol) - { + // Apply time window selection + if (hitT < m_time_min || hitT > m_time_max) + { + streamlog_out(DEBUG4) << "hit at T: " << hitT << " ns is rejected by timing cuts" << std::endl; + continue; + } - for (auto &ihit : hits_to_save[icol]) - { + if (m_fillHistos) + m_corrected_time_after->Fill(hitT); - TrackerHitPlane *hit = dynamic_cast(inputHitColls[icol]->getElementAt(ihit)); - TrackerHitPlaneImpl *hit_new = new TrackerHitPlaneImpl(); - - hit_new->setCellID0(hit->getCellID0()); - hit_new->setCellID1(hit->getCellID1()); - hit_new->setType(hit->getType()); - hit_new->setPosition(hit->getPosition()); - hit_new->setU(hit->getU()); - hit_new->setV(hit->getV()); - hit_new->setdU(hit->getdU()); - hit_new->setdV(hit->getdV()); - hit_new->setEDep(hit->getEDep()); - hit_new->setEDepError(hit->getEDepError()); - hit_new->setTime(hit->getTime()); - hit_new->setQuality(hit->getQuality()); - - const lcio::LCObjectVec &rawHits = hit->getRawHits(); - for (size_t j=0; j( rawHits[j] ); - hit_new->rawHits().push_back(hitConstituent); - } + // Save a copy of the cluster into the output collection + TrackerHitPlane *hit_new = copyTrackerHitPlane(hit); outputTrackerHitColls[icol]->addElement(hit_new); + if (outputTrackerHitConstituentsColls[icol]) + { + // save hits contituents as well + const EVENT::LCObjectVec &rawHits = hit_new->getRawHits(); + for (unsigned int ihit = 0; ihit < rawHits.size(); ++ihit) + outputTrackerHitConstituentsColls[icol]->addElement(rawHits[ihit]); + } - LCRelation *rel = dynamic_cast(inputHitRels[icol]->getElementAt(ihit)); - - SimTrackerHit *simhit = dynamic_cast(rel->getTo()); - SimTrackerHitImpl *simhit_new = new SimTrackerHitImpl(); - - simhit_new->setCellID0(simhit->getCellID0()); - simhit_new->setCellID1(simhit->getCellID1()); - simhit_new->setPosition(simhit->getPosition()); - simhit_new->setEDep(simhit->getEDep()); - simhit_new->setTime(simhit->getTime()); - simhit_new->setMCParticle(simhit->getMCParticle()); - simhit_new->setMomentum(simhit->getMomentum()); - simhit_new->setPathLength(simhit->getPathLength()); - simhit_new->setQuality(simhit->getQuality()); - simhit_new->setOverlay(simhit->isOverlay()); - simhit_new->setProducedBySecondary(simhit->isProducedBySecondary()); - - outputTrackerSimHitColls[icol]->addElement(simhit_new); - - LCRelationImpl *rel_new = new LCRelationImpl(); - - rel_new->setFrom(hit_new); - rel_new->setTo(simhit_new); - rel_new->setWeight(rel->getWeight()); - - outputTrackerHitRels[icol]->addElement(rel_new); + // Keep track of clusters saved + reco_clusters_to_save[hit] = hit_new; } // ihit loop - streamlog_out(MESSAGE) << " " << hits_to_save[icol].size() << " hits added to the collections: " - << m_outputTrackerHitsCollNames[icol] << ", " - << m_outputTrackerSimHitsCollNames[icol] << ", " + // Now loop over the LCRelation collection, and save new SimHits as well as their reco-sim relation + LCCollection *hitrel_col = inputHitRels[icol]; + if (hitrel_col) + { + for (unsigned int irel = 0; irel < hitrel_col->getNumberOfElements(); ++irel) + { + LCRelation *rel = dynamic_cast(hitrel_col->getElementAt(irel)); + // now check if the pointer to the reco object (From) is in the list of cluster to save + TrackerHitPlane *recohit = dynamic_cast(rel->getFrom()); + auto reco_hit_itr = reco_clusters_to_save.find(recohit); + if (reco_hit_itr != reco_clusters_to_save.end()) + { + // Reconstructed Hit (already saved in output collection) + TrackerHit *hit_new = reco_hit_itr->second; + + // Simulated Hit + SimTrackerHit *simhit_new=nullptr; + if (outputTrackerSimHitColls[icol] != nullptr) { + SimTrackerHit *simhit = dynamic_cast(rel->getTo()); + simhit_new = copySimTrackerHit(simhit); + outputTrackerSimHitColls[icol]->addElement(simhit_new); + } + + // LCRelation + LCRelationImpl *rel_new = new LCRelationImpl(); + rel_new->setFrom(hit_new); + rel_new->setTo(simhit_new); + rel_new->setWeight(rel->getWeight()); + outputTrackerHitRels[icol]->addElement(rel_new); + } + } // irel loop + } // hitrel available + + streamlog_out(MESSAGE) << " " << reco_clusters_to_save.size() << " hits added to the collections: " + << m_outputTrackerHitsCollNames[icol] << ", "; + if (m_outputTrackerHitsConstituentsCollNames[icol] != "") + streamlog_out(MESSAGE) << m_outputTrackerHitsConstituentsCollNames[icol] << ", "; + streamlog_out(MESSAGE) << m_outputTrackerSimHitsCollNames[icol] << ", " << m_outputTrackerHitRelNames[icol] << std::endl; evt->addCollection(outputTrackerHitColls[icol], m_outputTrackerHitsCollNames[icol]); + if (outputTrackerHitConstituentsColls[icol]) + evt->addCollection(outputTrackerHitConstituentsColls[icol], m_outputTrackerHitsConstituentsCollNames[icol]); evt->addCollection(outputTrackerSimHitColls[icol], m_outputTrackerSimHitsCollNames[icol]); evt->addCollection(outputTrackerHitRels[icol], m_outputTrackerHitRelNames[icol]); @@ -326,8 +425,12 @@ void FilterTimeHits::processEvent(LCEvent *evt) << " output collection " << m_outputTrackerSimHitsCollNames[icol] << " of type " << outputTrackerSimHitColls[icol]->getTypeName() << " added to the event \n" << " output collection " << m_outputTrackerHitRelNames[icol] << " of type " - << outputTrackerHitRels[icol]->getTypeName() << " added to the event " - << std::endl; + << outputTrackerHitRels[icol]->getTypeName() << " added to the event "; + if (outputTrackerHitConstituentsColls[icol]) { + streamlog_out(DEBUG5) << " output collection " << m_outputTrackerHitsConstituentsCollNames[icol] << " of type " + << outputTrackerHitConstituentsColls[icol]->getTypeName() << " added to the event \n"; + } + streamlog_out(DEBUG5) << std::endl; } // icol loop From 69942819bf9508c430b1f7e67a1f31ec5b97cf2c Mon Sep 17 00:00:00 2001 From: Simone Pagan Griso Date: Mon, 4 Nov 2024 16:58:12 -0800 Subject: [PATCH 11/14] added registerParameter marlin functions --- source/Utils/src/FilterTimeHits.cc | 54 +++++++++++++++++++++--------- 1 file changed, 39 insertions(+), 15 deletions(-) diff --git a/source/Utils/src/FilterTimeHits.cc b/source/Utils/src/FilterTimeHits.cc index ba49937..2f95f37 100644 --- a/source/Utils/src/FilterTimeHits.cc +++ b/source/Utils/src/FilterTimeHits.cc @@ -32,6 +32,11 @@ FilterTimeHits::FilterTimeHits() : Processor("FilterTimeHits") m_inputTrackerHitsCollNames, {}); + registerProcessorParameter("TrackerHitConstituentsInputCollections", + "Name of the tracker hit constituents input collections", + m_inputTrackerHitsConstituentsCollNames, + {}); + registerProcessorParameter("TrackerSimHitInputCollections", "Name of the tracker simhit input collections", m_inputTrackerSimHitsCollNames, @@ -47,6 +52,11 @@ FilterTimeHits::FilterTimeHits() : Processor("FilterTimeHits") m_outputTrackerHitsCollNames, {}); + registerProcessorParameter("TrackerHitConstituentsOutputCollections", + "Name of the tracker hit output collections", + m_outputTrackerHitsConstituentsCollNames, + {}); + registerProcessorParameter("TrackerSimHitOutputCollections", "Name of the tracker simhit output collections", m_outputTrackerSimHitsCollNames, @@ -232,9 +242,13 @@ void FilterTimeHits::processEvent(LCEvent *evt) } // get the reco hits contituents (if available) + bool hasHitConstituents = false; try { - if (m_inputTrackerHitsConstituentsCollNames[icol] != "") + if (m_inputTrackerHitsConstituentsCollNames.size() > 0) + if (m_inputTrackerHitsConstituentsCollNames[icol] != "") + hasHitConstituents = true; + if (hasHitConstituents) inputHitConstituentsColls[icol] = evt->getCollection(m_inputTrackerHitsConstituentsCollNames[icol]); else inputHitConstituentsColls[icol] = nullptr; @@ -284,19 +298,25 @@ void FilterTimeHits::processEvent(LCEvent *evt) outputTrackerHitColls[icol]->setFlag(lcFlag.getFlag()); // reco hit contituents output collections, if needed - if ((m_outputTrackerHitsConstituentsCollNames[icol] != "") && - (inputHitConstituentsColls[icol] != nullptr)) { - std::string encoderString = inputHitConstituentsColls[icol]->getParameters().getStringVal("CellIDEncoding"); - outputTrackerHitConstituentsColls[icol] = new LCCollectionVec(inputHitConstituentsColls[icol]->getTypeName()); - outputTrackerHitConstituentsColls[icol]->parameters().setValue("CellIDEncoding", encoderString); - LCFlagImpl lcFlag(inputHitConstituentsColls[icol]->getFlag()); - outputTrackerHitConstituentsColls[icol]->setFlag(lcFlag.getFlag()); + if (hasHitConstituents) { + if ((m_outputTrackerHitsConstituentsCollNames[icol] != "") && + (inputHitConstituentsColls[icol] != nullptr)) + { + std::string encoderString = inputHitConstituentsColls[icol]->getParameters().getStringVal("CellIDEncoding"); + outputTrackerHitConstituentsColls[icol] = new LCCollectionVec(inputHitConstituentsColls[icol]->getTypeName()); + outputTrackerHitConstituentsColls[icol]->parameters().setValue("CellIDEncoding", encoderString); + LCFlagImpl lcFlag(inputHitConstituentsColls[icol]->getFlag()); + outputTrackerHitConstituentsColls[icol]->setFlag(lcFlag.getFlag()); + } else { + outputTrackerHitConstituentsColls[icol] = nullptr; + } } else { outputTrackerHitConstituentsColls[icol] = nullptr; } // sim hit output collections - if (inputSimHitColls[icol] != nullptr) { + if (inputSimHitColls[icol] != nullptr) + { outputTrackerSimHitColls[icol] = new LCCollectionVec(inputSimHitColls[icol]->getTypeName()); outputTrackerSimHitColls[icol]->parameters().setValue("CellIDEncoding", encoderString); LCFlagImpl lcFlag_sim(inputSimHitColls[icol]->getFlag()); @@ -304,7 +324,8 @@ void FilterTimeHits::processEvent(LCEvent *evt) } // reco-sim relation output collections - if (inputHitRels[icol] != nullptr) { + if (inputHitRels[icol] != nullptr) + { outputTrackerHitRels[icol] = new LCCollectionVec(inputHitRels[icol]->getTypeName()); LCFlagImpl lcFlag_rel(inputHitRels[icol]->getFlag()); outputTrackerHitRels[icol]->setFlag(lcFlag_rel.getFlag()); @@ -390,8 +411,9 @@ void FilterTimeHits::processEvent(LCEvent *evt) TrackerHit *hit_new = reco_hit_itr->second; // Simulated Hit - SimTrackerHit *simhit_new=nullptr; - if (outputTrackerSimHitColls[icol] != nullptr) { + SimTrackerHit *simhit_new = nullptr; + if (outputTrackerSimHitColls[icol] != nullptr) + { SimTrackerHit *simhit = dynamic_cast(rel->getTo()); simhit_new = copySimTrackerHit(simhit); outputTrackerSimHitColls[icol]->addElement(simhit_new); @@ -409,8 +431,9 @@ void FilterTimeHits::processEvent(LCEvent *evt) streamlog_out(MESSAGE) << " " << reco_clusters_to_save.size() << " hits added to the collections: " << m_outputTrackerHitsCollNames[icol] << ", "; - if (m_outputTrackerHitsConstituentsCollNames[icol] != "") - streamlog_out(MESSAGE) << m_outputTrackerHitsConstituentsCollNames[icol] << ", "; + if (m_inputTrackerHitsConstituentsCollNames.size() > 0) + if (m_outputTrackerHitsConstituentsCollNames[icol] != "") + streamlog_out(MESSAGE) << m_outputTrackerHitsConstituentsCollNames[icol] << ", "; streamlog_out(MESSAGE) << m_outputTrackerSimHitsCollNames[icol] << ", " << m_outputTrackerHitRelNames[icol] << std::endl; @@ -426,7 +449,8 @@ void FilterTimeHits::processEvent(LCEvent *evt) << outputTrackerSimHitColls[icol]->getTypeName() << " added to the event \n" << " output collection " << m_outputTrackerHitRelNames[icol] << " of type " << outputTrackerHitRels[icol]->getTypeName() << " added to the event "; - if (outputTrackerHitConstituentsColls[icol]) { + if (outputTrackerHitConstituentsColls[icol]) + { streamlog_out(DEBUG5) << " output collection " << m_outputTrackerHitsConstituentsCollNames[icol] << " of type " << outputTrackerHitConstituentsColls[icol]->getTypeName() << " added to the event \n"; } From b4352a1e327677302cae768d5f90fe28094a7e24 Mon Sep 17 00:00:00 2001 From: Simone Pagan Griso Date: Tue, 5 Nov 2024 09:38:35 -0800 Subject: [PATCH 12/14] fix cases where size of filter arrays is different from number of laters --- source/Utils/src/FilterClusters.cc | 521 +++++++++++++++-------------- 1 file changed, 274 insertions(+), 247 deletions(-) diff --git a/source/Utils/src/FilterClusters.cc b/source/Utils/src/FilterClusters.cc index f6953c8..ed58143 100644 --- a/source/Utils/src/FilterClusters.cc +++ b/source/Utils/src/FilterClusters.cc @@ -21,89 +21,77 @@ #include "marlin/VerbosityLevels.h" - -FilterClusters aFilterClusters ; +FilterClusters aFilterClusters; FilterClusters::FilterClusters() - : Processor("FilterClusters") + : Processor("FilterClusters") { // modify processor description _description = "FilterClusters processor filters a collection of tracker hits based on cluster size and outputs a filtered collection"; // register steering parameters: name, description, class-variable, default value registerProcessorParameter("ThetaRanges", - "Divide theta into bins for different cluster size cuts", - _ThetaRanges, - {} - ); + "Divide theta into bins for different cluster size cuts", + _ThetaRanges, + {}); registerProcessorParameter("ThetaBins", - "Number of bins in theta", - _ThetaBins, - {} - ); + "Number of bins in theta", + _ThetaBins, + {}); registerProcessorParameter("ClusterSize", - "Maximum cluster size for each theta range", - _ClusterSize, - {} - ); + "Maximum cluster size for each theta range", + _ClusterSize, + {}); registerProcessorParameter("Layers", - "Layers to be filtered", - _Layers, - {} - ); - registerInputCollection( LCIO::SIMTRACKERHIT, - "InSimTrackerHitCollection" , - "Name of the input sim tracker hit collection", - _InSimTrackerHitCollection, - _InSimTrackerHitCollection - ); - - registerInputCollection( LCIO::TRACKERHIT, - "InTrackerHitCollection" , - "Name of the input tracker hit collection", - _InTrackerHitCollection, - _InTrackerHitCollection - ); - - registerInputCollection( LCIO::LCRELATION, - "InRelationCollection" , - "Name of the input relation collection", - _InRelationCollection, - _InRelationCollection - ); - - registerOutputCollection( LCIO::SIMTRACKERHIT, - "OutSimTrackerHitCollection" , - "Name of output sim tracker hit collection", - _OutSimTrackerHitCollection, - _OutSimTrackerHitCollection - ); - - registerOutputCollection( LCIO::TRACKERHIT, - "OutTrackerHitCollection" , - "Name of output tracker hit collection", - _OutTrackerHitCollection, - _OutTrackerHitCollection - ); - - registerOutputCollection( LCIO::LCRELATION, - "OutRelationCollection" , - "Name of output relation collection", - _OutRelationCollection, - _OutRelationCollection - ); + "Layers to be filtered", + _Layers, + {}); + registerInputCollection(LCIO::SIMTRACKERHIT, + "InSimTrackerHitCollection", + "Name of the input sim tracker hit collection", + _InSimTrackerHitCollection, + _InSimTrackerHitCollection); + + registerInputCollection(LCIO::TRACKERHIT, + "InTrackerHitCollection", + "Name of the input tracker hit collection", + _InTrackerHitCollection, + _InTrackerHitCollection); + + registerInputCollection(LCIO::LCRELATION, + "InRelationCollection", + "Name of the input relation collection", + _InRelationCollection, + _InRelationCollection); + + registerOutputCollection(LCIO::SIMTRACKERHIT, + "OutSimTrackerHitCollection", + "Name of output sim tracker hit collection", + _OutSimTrackerHitCollection, + _OutSimTrackerHitCollection); + + registerOutputCollection(LCIO::TRACKERHIT, + "OutTrackerHitCollection", + "Name of output tracker hit collection", + _OutTrackerHitCollection, + _OutTrackerHitCollection); + + registerOutputCollection(LCIO::LCRELATION, + "OutRelationCollection", + "Name of output relation collection", + _OutRelationCollection, + _OutRelationCollection); registerProcessorParameter("FillHistograms", - "Flag to fill the diagnostic histograms", - m_fillHistos, - false); - + "Flag to fill the diagnostic histograms", + m_fillHistos, + false); } void FilterClusters::init() { // Print the initial parameters - printParameters() ; + printParameters(); marlin::AIDAProcessor::histogramFactory(this); @@ -113,225 +101,264 @@ void FilterClusters::init() m_clusterLayer_afterCut = new TH1F("m_ClusterLayer_after", "Cluster Layer Index", 10, 0, 10); m_clusterSize_beforeCut = new TH1F("m_ClusterSize_before", "Cluster hit multiplicity", 20, 0, 20); m_clusterSize_afterCut = new TH1F("m_ClusterSize_after", "Cluster hit multiplicity", 20, 0, 20); - } -void FilterClusters::processRunHeader( LCRunHeader* /*run*/) -{ } +void FilterClusters::processRunHeader(LCRunHeader * /*run*/) +{ +} -void FilterClusters::processEvent( LCEvent * evt ) +void FilterClusters::processEvent(LCEvent *evt) { // Get input collection - LCCollection* InTrackerHitCollection = evt->getCollection(_InTrackerHitCollection); - LCCollection* InRelationCollection = evt->getCollection(_InRelationCollection); - LCCollection* InSimTrackerHitCollection = evt->getCollection(_InSimTrackerHitCollection); - - if( InTrackerHitCollection->getTypeName() != lcio::LCIO::TRACKERHITPLANE ) - { throw EVENT::Exception( "Invalid collection type: " + InTrackerHitCollection->getTypeName() ) ; - streamlog_out(DEBUG0) << "Wrong collection type for TrackerHitCollection. \n";} - if( InRelationCollection->getTypeName() != lcio::LCIO::LCRELATION ) - { throw EVENT::Exception( "Invalid collection type: " + InRelationCollection->getTypeName() ) ; - streamlog_out(DEBUG0) << "Wrong collection type for InRelationCollection. \n";} - if( InSimTrackerHitCollection->getTypeName() != lcio::LCIO::SIMTRACKERHIT ) - { throw EVENT::Exception( "Invalid collection type: " + InSimTrackerHitCollection->getTypeName() ) ; - streamlog_out(DEBUG0) << "Wrong collection type for SimTrackerHitCollection. \n";} - - - streamlog_out(DEBUG0) << "Number of Elements in Tracker Hits Collection: " << InTrackerHitCollection->getNumberOfElements() <getCollection(_InTrackerHitCollection); + LCCollection *InRelationCollection = evt->getCollection(_InRelationCollection); + LCCollection *InSimTrackerHitCollection = evt->getCollection(_InSimTrackerHitCollection); + + streamlog_out(DEBUG0) << "Starting processing event." << std::endl; + + if (InTrackerHitCollection->getTypeName() != lcio::LCIO::TRACKERHITPLANE) + { + throw EVENT::Exception("Invalid collection type: " + InTrackerHitCollection->getTypeName()); + streamlog_out(DEBUG0) << "Wrong collection type for TrackerHitCollection. \n"; + } + if (InRelationCollection->getTypeName() != lcio::LCIO::LCRELATION) + { + throw EVENT::Exception("Invalid collection type: " + InRelationCollection->getTypeName()); + streamlog_out(DEBUG0) << "Wrong collection type for InRelationCollection. \n"; + } + if (InSimTrackerHitCollection->getTypeName() != lcio::LCIO::SIMTRACKERHIT) + { + throw EVENT::Exception("Invalid collection type: " + InSimTrackerHitCollection->getTypeName()); + streamlog_out(DEBUG0) << "Wrong collection type for SimTrackerHitCollection. \n"; + } + + streamlog_out(DEBUG1) << "Number of Elements in Tracker Hits Collection: " << InTrackerHitCollection->getNumberOfElements() << std::endl; // Make the output collections: reco hits, sim hits, reco-sim relationship - std::string encoderString = InTrackerHitCollection->getParameters().getStringVal("CellIDEncoding"); + std::string encoderString = InTrackerHitCollection->getParameters().getStringVal("CellIDEncoding"); LCCollectionVec *OutTrackerHitCollection = new LCCollectionVec(LCIO::TRACKERHIT); OutTrackerHitCollection->parameters().setValue("CellIDEncoding", encoderString); LCFlagImpl lcFlag(InTrackerHitCollection->getFlag()); - OutTrackerHitCollection->setFlag(lcFlag.getFlag()); - //OutTrackerHitCollection->setSubset(true); + OutTrackerHitCollection->setFlag(lcFlag.getFlag()); + // OutTrackerHitCollection->setSubset(true); LCCollectionVec *OutSimTrackerHitCollection = new LCCollectionVec(LCIO::SIMTRACKERHIT); OutSimTrackerHitCollection->parameters().setValue("CellIDEncoding", encoderString); LCFlagImpl lcFlag_sim(InSimTrackerHitCollection->getFlag()); OutSimTrackerHitCollection->setFlag(lcFlag_sim.getFlag()); - //OutSimTrackerHitCollection->setSubset(true); - - LCCollectionVec *OutRelationCollection = new LCCollectionVec(LCIO::LCRELATION); + // OutSimTrackerHitCollection->setSubset(true); + + LCCollectionVec *OutRelationCollection = new LCCollectionVec(LCIO::LCRELATION); LCFlagImpl lcFlag_rel(InRelationCollection->getFlag()); - OutRelationCollection->setFlag(lcFlag_rel.getFlag()); - //OutRelationCollection->setSubset(true); + OutRelationCollection->setFlag(lcFlag_rel.getFlag()); + // OutRelationCollection->setSubset(true); - // Filter - for(int i=0; igetNumberOfElements(); ++i) //loop through all hits + for (int i = 0; i < InTrackerHitCollection->getNumberOfElements(); ++i) // loop through all hits + { + streamlog_out(DEBUG2) << "Loop over hits:" << i << std::endl; + TrackerHitPlane *trkhit = static_cast(InTrackerHitCollection->getElementAt(i)); // define trkhit var, pointer to i'th element of tracker hits + + if (!trkhit) { - streamlog_out(DEBUG0) << "Loop over hits opened. \n"; - TrackerHitPlane *trkhit=static_cast(InTrackerHitCollection->getElementAt(i)); //define trkhit var, pointer to i'th element of tracker hits - - //Calculating theta - float x = trkhit->getPosition()[0]; - float y = trkhit->getPosition()[1]; - float z = trkhit->getPosition()[2]; - float r = sqrt(pow(x,2)+pow(y,2)); - float incidentTheta = std::atan(r/z); - - if(incidentTheta<0) - incidentTheta += M_PI; - - //Calculating cluster size - const lcio::LCObjectVec &rawHits = trkhit->getRawHits(); - float ymax = -1000000; - float xmax = -1000000; - float ymin = 1000000; - float xmin = 1000000; - for (size_t j=0; j( rawHits[j] ); - const double *localPos = hitConstituent->getPosition(); - float x_local = localPos[0]; - float y_local = localPos[1]; - - if (y_local < ymin){ - ymin = y_local; - } - if (y_local > ymax){ - ymax = y_local; - } - - if (x_local < xmin){ - xmin = x_local; - } - if (x_local > xmax){ - xmax = x_local; - } + streamlog_out(WARNING) << "Cannot retrieve valid point to cluster. Skipping it" << std::endl; + } + + // Calculating theta + float x = trkhit->getPosition()[0]; + float y = trkhit->getPosition()[1]; + float z = trkhit->getPosition()[2]; + float r = sqrt(pow(x, 2) + pow(y, 2)); + float incidentTheta = std::atan(r / z); + + if (incidentTheta < 0) + incidentTheta += M_PI; + + // Calculating cluster size + const lcio::LCObjectVec &rawHits = trkhit->getRawHits(); + float ymax = -1000000; + float xmax = -1000000; + float ymin = 1000000; + float xmin = 1000000; + streamlog_out(DEBUG1) << "Looping over hits constituents." << std::endl; + for (size_t j = 0; j < rawHits.size(); ++j) + { + lcio::SimTrackerHit *hitConstituent = dynamic_cast(rawHits[j]); + const double *localPos = hitConstituent->getPosition(); + float x_local = localPos[0]; + float y_local = localPos[1]; + + if (y_local < ymin) + { + ymin = y_local; + } + if (y_local > ymax) + { + ymax = y_local; } - float cluster_size_y = (ymax - ymin)+1; - float cluster_size_x = (xmax - xmin)+1; - float cluster_size = rawHits.size(); - - //Get hit subdetector/layer - std::string _encoderString = lcio::LCTrackerCellID::encoding_string(); - UTIL::CellIDDecoder decoder(_encoderString); - uint32_t systemID = decoder(trkhit)["system"]; - uint32_t layerID = decoder(trkhit)["layer"]; - bool filter_layer = false; - - int rows = _Layers.size(), cols = std::stoi(_ThetaBins); - if((rows*cols != _ThetaRanges.size()) || (rows*(cols-1) != _ClusterSize.size())){ - std::cout<<"Either theta cuts or cluster cuts not provided for each layer. Please change the config, exiting now..."<> _thetaCuts_byLayer; - std::vector> _clusterSizeCuts_byLayer; - - for (int k = 0; k < rows; ++k) { - std::vector row; - for (int j = 0; j < cols; ++j) { - row.push_back(std::stof(_ThetaRanges[j])); - } - _thetaCuts_byLayer.push_back(row); + if (x_local > xmax) + { + xmax = x_local; } + } + float cluster_size_y = (ymax - ymin) + 1; + float cluster_size_x = (xmax - xmin) + 1; + float cluster_size = rawHits.size(); + + streamlog_out(DEBUG2) << "Cluster size:" << cluster_size << std::endl; + + // Get hit subdetector/layer + std::string _encoderString = lcio::LCTrackerCellID::encoding_string(); + UTIL::CellIDDecoder decoder(_encoderString); + uint32_t systemID = decoder(trkhit)["system"]; + uint32_t layerID = decoder(trkhit)["layer"]; + bool filter_layer = false; - for (int k = 0; k < rows; ++k) { - std::vector row; - for (int j = 0; j < cols; ++j) { - row.push_back(std::stof(_ClusterSize[j])); - } - _clusterSizeCuts_byLayer.push_back(row); + streamlog_out(DEBUG3) << "Decoded system/layer:" << systemID << "/" << layerID << std::endl; + + int rows = _Layers.size(), cols = std::stoi(_ThetaBins); + if ((rows * cols != _ThetaRanges.size()) || (rows * (cols - 1) != _ClusterSize.size())) + { + std::cout << "Either theta cuts or cluster cuts not provided for each layer. Please change the config, exiting now..." << std::endl; + return; + } + + std::vector> _thetaCuts_byLayer; + std::vector> _clusterSizeCuts_byLayer; + + for (int k = 0; k < rows; ++k) + { + std::vector row; + for (int j = 0; j < cols; ++j) + { + row.push_back(std::stof(_ThetaRanges[j])); } + _thetaCuts_byLayer.push_back(row); + } - - for (size_t j=0; j<_Layers.size(); ++j){ - if (layerID == std::stof(_Layers[j])) { - filter_layer = true; - break; - } + for (int k = 0; k < rows; ++k) + { + std::vector row; + for (int j = 0; j < cols; ++j) + { + row.push_back(std::stof(_ClusterSize[j])); } - streamlog_out(DEBUG0) << "Filter layer: " << filter_layer << std::endl; + _clusterSizeCuts_byLayer.push_back(row); + } - if (m_fillHistos){ - m_clusterTheta_beforeCut->Fill(incidentTheta); - m_clusterLayer_beforeCut->Fill(layerID); - m_clusterSize_beforeCut->Fill(cluster_size); + for (size_t j = 0; j < _Layers.size(); ++j) + { + if (layerID == std::stof(_Layers[j])) + { + filter_layer = true; + break; } - - for (size_t j=0; j<_thetaCuts_byLayer[layerID].size()-1; ++j) { - streamlog_out( DEBUG0 ) << "theta: " << incidentTheta << std::endl; - float min_theta = _thetaCuts_byLayer[layerID][j]; - float max_theta = _thetaCuts_byLayer[layerID][j+1]; - streamlog_out( DEBUG0 ) << "theta range: " << min_theta << ", " << max_theta << std::endl; - - if(incidentTheta >= min_theta and incidentTheta <= max_theta and filter_layer){ - streamlog_out( DEBUG0 ) << "theta in range" << std::endl; - streamlog_out( DEBUG0 ) << "cluster size cut off: " << _clusterSizeCuts_byLayer[layerID][j] << std::endl; - streamlog_out( DEBUG0 ) << "cluster size: " << cluster_size << std::endl; - if(cluster_size < _clusterSizeCuts_byLayer[layerID][j]) { - streamlog_out( DEBUG0 ) << "Adding reco/sim clusters and relation to output collections" << std::endl; - - if (m_fillHistos){ - m_clusterTheta_afterCut->Fill(incidentTheta); - m_clusterLayer_afterCut->Fill(layerID); - m_clusterSize_afterCut->Fill(cluster_size); - } - - EVENT::LCRelation *rel=static_cast(InRelationCollection->getElementAt(i)); - SimTrackerHit *simhit = dynamic_cast(rel->getTo()); - SimTrackerHitImpl *simhit_new = new SimTrackerHitImpl(); - simhit_new->setCellID0(simhit->getCellID0()); - simhit_new->setCellID1(simhit->getCellID1()); - simhit_new->setPosition(simhit->getPosition()); - simhit_new->setEDep(simhit->getEDep()); - simhit_new->setTime(simhit->getTime()); - simhit_new->setMCParticle(simhit->getMCParticle()); - simhit_new->setMomentum(simhit->getMomentum()); - simhit_new->setPathLength(simhit->getPathLength()); - simhit_new->setQuality(simhit->getQuality()); - simhit_new->setOverlay(simhit->isOverlay()); - simhit_new->setProducedBySecondary(simhit->isProducedBySecondary()); - OutSimTrackerHitCollection->addElement(simhit_new); - - TrackerHitPlaneImpl *hit_new = new TrackerHitPlaneImpl(); - hit_new->setCellID0(trkhit->getCellID0()); - hit_new->setCellID1(trkhit->getCellID1()); - hit_new->setType(trkhit->getType()); - hit_new->setPosition(trkhit->getPosition()); - hit_new->setU(trkhit->getU()); - hit_new->setV(trkhit->getV()); - hit_new->setdU(trkhit->getdU()); - hit_new->setdV(trkhit->getdV()); - hit_new->setEDep(trkhit->getEDep()); - hit_new->setEDepError(trkhit->getEDepError()); - hit_new->setTime(trkhit->getTime()); - hit_new->setQuality(trkhit->getQuality()); - const lcio::LCObjectVec &rawHits_new = trkhit->getRawHits(); - for (size_t k=0; k( rawHits_new[k] ); - hit_new->rawHits().push_back(hitConstituent); - } - OutTrackerHitCollection->addElement(hit_new); - - LCRelationImpl *rel_new = new LCRelationImpl(); - rel_new->setFrom(hit_new); - rel_new->setTo(simhit_new); - rel_new->setWeight(rel->getWeight()); - OutRelationCollection->addElement(rel_new); - } - else { - streamlog_out( DEBUG0 ) << "cluster rejected" << std::endl; + } + streamlog_out(DEBUG1) << "Filter layer: " << filter_layer << std::endl; + + if (m_fillHistos) + { + m_clusterTheta_beforeCut->Fill(incidentTheta); + m_clusterLayer_beforeCut->Fill(layerID); + m_clusterSize_beforeCut->Fill(cluster_size); + } + + bool store_hit = true; + if (filter_layer) + { + store_hit = false; + for (size_t j = 0; j < _thetaCuts_byLayer[layerID].size() - 1; ++j) + { + streamlog_out(DEBUG0) << "theta: " << incidentTheta << std::endl; + float min_theta = _thetaCuts_byLayer[layerID][j]; + float max_theta = _thetaCuts_byLayer[layerID][j + 1]; + streamlog_out(DEBUG0) << "theta range: " << min_theta << ", " << max_theta << std::endl; + + if (incidentTheta >= min_theta and incidentTheta <= max_theta and filter_layer) + { + store_hit = true; + streamlog_out(DEBUG0) << "theta in range" << std::endl; + streamlog_out(DEBUG0) << "cluster size cut off: " << _clusterSizeCuts_byLayer[layerID][j] << std::endl; + streamlog_out(DEBUG0) << "cluster size: " << cluster_size << std::endl; + if (cluster_size < _clusterSizeCuts_byLayer[layerID][j]) + { + streamlog_out(DEBUG0) << "Adding reco/sim clusters and relation to output collections" << std::endl; } } - else{ - streamlog_out( DEBUG0 ) << "theta out of range or filtering not enabled for this layer" << std::endl; - } } } + if (store_hit) + { + if (m_fillHistos) + { + m_clusterTheta_afterCut->Fill(incidentTheta); + m_clusterLayer_afterCut->Fill(layerID); + m_clusterSize_afterCut->Fill(cluster_size); + } + + EVENT::LCRelation *rel = static_cast(InRelationCollection->getElementAt(i)); + SimTrackerHit *simhit = dynamic_cast(rel->getTo()); + SimTrackerHitImpl *simhit_new = new SimTrackerHitImpl(); + simhit_new->setCellID0(simhit->getCellID0()); + simhit_new->setCellID1(simhit->getCellID1()); + simhit_new->setPosition(simhit->getPosition()); + simhit_new->setEDep(simhit->getEDep()); + simhit_new->setTime(simhit->getTime()); + simhit_new->setMCParticle(simhit->getMCParticle()); + simhit_new->setMomentum(simhit->getMomentum()); + simhit_new->setPathLength(simhit->getPathLength()); + simhit_new->setQuality(simhit->getQuality()); + simhit_new->setOverlay(simhit->isOverlay()); + simhit_new->setProducedBySecondary(simhit->isProducedBySecondary()); + OutSimTrackerHitCollection->addElement(simhit_new); + + TrackerHitPlaneImpl *hit_new = new TrackerHitPlaneImpl(); + hit_new->setCellID0(trkhit->getCellID0()); + hit_new->setCellID1(trkhit->getCellID1()); + hit_new->setType(trkhit->getType()); + hit_new->setPosition(trkhit->getPosition()); + hit_new->setU(trkhit->getU()); + hit_new->setV(trkhit->getV()); + hit_new->setdU(trkhit->getdU()); + hit_new->setdV(trkhit->getdV()); + hit_new->setEDep(trkhit->getEDep()); + hit_new->setEDepError(trkhit->getEDepError()); + hit_new->setTime(trkhit->getTime()); + hit_new->setQuality(trkhit->getQuality()); + const lcio::LCObjectVec &rawHits_new = trkhit->getRawHits(); + for (size_t k = 0; k < rawHits_new.size(); ++k) + { + lcio::SimTrackerHit *hitConstituent = dynamic_cast(rawHits_new[k]); + hit_new->rawHits().push_back(hitConstituent); + } + OutTrackerHitCollection->addElement(hit_new); + + LCRelationImpl *rel_new = new LCRelationImpl(); + rel_new->setFrom(hit_new); + rel_new->setTo(simhit_new); + rel_new->setWeight(rel->getWeight()); + OutRelationCollection->addElement(rel_new); + } + else + { + streamlog_out(DEBUG0) << "cluster rejected" << std::endl; + } + } + // Save output track collection - evt->addCollection(OutTrackerHitCollection, _OutTrackerHitCollection); + evt->addCollection(OutTrackerHitCollection, _OutTrackerHitCollection); evt->addCollection(OutRelationCollection, _OutRelationCollection); evt->addCollection(OutSimTrackerHitCollection, _OutSimTrackerHitCollection); - streamlog_out(MESSAGE) << " " << OutTrackerHitCollection->size() << " reco clusters added to the collection: " << _OutTrackerHitCollection << ", " << OutSimTrackerHitCollection->size() << " sim hits added to the collection: "<< _OutSimTrackerHitCollection << " and " << OutRelationCollection->size() << " relations added to the collection: " << _OutRelationCollection << std::endl; - + streamlog_out(MESSAGE) << " " << OutTrackerHitCollection->size() << " reco clusters added to the collection: " << _OutTrackerHitCollection << ", " << OutSimTrackerHitCollection->size() << " sim hits added to the collection: " << _OutSimTrackerHitCollection << " and " << OutRelationCollection->size() << " relations added to the collection: " << _OutRelationCollection << std::endl; } void FilterClusters::end() -{ } +{ +} From 1218a648ad41058b755e3448269ced3c3a771a3c Mon Sep 17 00:00:00 2001 From: Angira Rastogi Date: Tue, 5 Nov 2024 23:40:32 -0800 Subject: [PATCH 13/14] Changing the output collection type from TrackerHit to TrackerHitPlane --- source/Utils/src/FilterClusters.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/source/Utils/src/FilterClusters.cc b/source/Utils/src/FilterClusters.cc index ed58143..09e584d 100644 --- a/source/Utils/src/FilterClusters.cc +++ b/source/Utils/src/FilterClusters.cc @@ -52,7 +52,7 @@ FilterClusters::FilterClusters() _InSimTrackerHitCollection, _InSimTrackerHitCollection); - registerInputCollection(LCIO::TRACKERHIT, + registerInputCollection(LCIO::TRACKERHITPLANE, "InTrackerHitCollection", "Name of the input tracker hit collection", _InTrackerHitCollection, @@ -70,7 +70,7 @@ FilterClusters::FilterClusters() _OutSimTrackerHitCollection, _OutSimTrackerHitCollection); - registerOutputCollection(LCIO::TRACKERHIT, + registerOutputCollection(LCIO::TRACKERHITPLANE, "OutTrackerHitCollection", "Name of output tracker hit collection", _OutTrackerHitCollection, @@ -136,7 +136,7 @@ void FilterClusters::processEvent(LCEvent *evt) // Make the output collections: reco hits, sim hits, reco-sim relationship std::string encoderString = InTrackerHitCollection->getParameters().getStringVal("CellIDEncoding"); - LCCollectionVec *OutTrackerHitCollection = new LCCollectionVec(LCIO::TRACKERHIT); + LCCollectionVec *OutTrackerHitCollection = new LCCollectionVec(LCIO::TRACKERHITPLANE); OutTrackerHitCollection->parameters().setValue("CellIDEncoding", encoderString); LCFlagImpl lcFlag(InTrackerHitCollection->getFlag()); OutTrackerHitCollection->setFlag(lcFlag.getFlag()); From 9aa7e107219b4903b7445f380405d6a9d8cd88e4 Mon Sep 17 00:00:00 2001 From: Angira Rastogi Date: Tue, 12 Nov 2024 12:12:12 -0800 Subject: [PATCH 14/14] Fixing the bug to store hits after cluster filter application, also reading the cuts properly from the digi config --- source/Utils/src/FilterClusters.cc | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/source/Utils/src/FilterClusters.cc b/source/Utils/src/FilterClusters.cc index 09e584d..c3ee2ab 100644 --- a/source/Utils/src/FilterClusters.cc +++ b/source/Utils/src/FilterClusters.cc @@ -222,23 +222,23 @@ void FilterClusters::processEvent(LCEvent *evt) streamlog_out(DEBUG3) << "Decoded system/layer:" << systemID << "/" << layerID << std::endl; int rows = _Layers.size(), cols = std::stoi(_ThetaBins); - if ((rows * cols != _ThetaRanges.size()) || (rows * (cols - 1) != _ClusterSize.size())) + if ((rows * (cols + 1) != _ThetaRanges.size()) || (rows * cols != _ClusterSize.size())) { std::cout << "Either theta cuts or cluster cuts not provided for each layer. Please change the config, exiting now..." << std::endl; return; } - std::vector> _thetaCuts_byLayer; + std::vector> _thetaBins_byLayer; std::vector> _clusterSizeCuts_byLayer; for (int k = 0; k < rows; ++k) { std::vector row; - for (int j = 0; j < cols; ++j) + for (int j = 0; j <= cols; ++j) { - row.push_back(std::stof(_ThetaRanges[j])); + row.push_back(std::stof(_ThetaRanges[j + k*(cols+1)])); } - _thetaCuts_byLayer.push_back(row); + _thetaBins_byLayer.push_back(row); } for (int k = 0; k < rows; ++k) @@ -246,16 +246,18 @@ void FilterClusters::processEvent(LCEvent *evt) std::vector row; for (int j = 0; j < cols; ++j) { - row.push_back(std::stof(_ClusterSize[j])); + row.push_back(std::stof(_ClusterSize[j + k*cols])); } _clusterSizeCuts_byLayer.push_back(row); } + int layerInd = -1; for (size_t j = 0; j < _Layers.size(); ++j) { if (layerID == std::stof(_Layers[j])) { filter_layer = true; + layerInd = j; break; } } @@ -272,22 +274,23 @@ void FilterClusters::processEvent(LCEvent *evt) if (filter_layer) { store_hit = false; - for (size_t j = 0; j < _thetaCuts_byLayer[layerID].size() - 1; ++j) + for (size_t j = 0; j < _thetaBins_byLayer[layerInd].size() - 1; ++j) { streamlog_out(DEBUG0) << "theta: " << incidentTheta << std::endl; - float min_theta = _thetaCuts_byLayer[layerID][j]; - float max_theta = _thetaCuts_byLayer[layerID][j + 1]; + float min_theta = _thetaBins_byLayer[layerInd][j]; + float max_theta = _thetaBins_byLayer[layerInd][j + 1]; streamlog_out(DEBUG0) << "theta range: " << min_theta << ", " << max_theta << std::endl; if (incidentTheta >= min_theta and incidentTheta <= max_theta and filter_layer) { - store_hit = true; streamlog_out(DEBUG0) << "theta in range" << std::endl; - streamlog_out(DEBUG0) << "cluster size cut off: " << _clusterSizeCuts_byLayer[layerID][j] << std::endl; + streamlog_out(DEBUG0) << "cluster size cut off: " << _clusterSizeCuts_byLayer[layerInd][j] << std::endl; streamlog_out(DEBUG0) << "cluster size: " << cluster_size << std::endl; - if (cluster_size < _clusterSizeCuts_byLayer[layerID][j]) + if (cluster_size < _clusterSizeCuts_byLayer[layerInd][j]) { + store_hit = true; streamlog_out(DEBUG0) << "Adding reco/sim clusters and relation to output collections" << std::endl; + break; } } }