diff --git a/Digitisers/G2CDArbor/src/G2CDArborAlg.cpp b/Digitisers/G2CDArbor/src/G2CDArborAlg.cpp index 2c5526aa506b8f5924be7ac2ba36e394a02f1865..f6897e6c5ac2184f99b8e35420345cfc5b42abf5 100644 --- a/Digitisers/G2CDArbor/src/G2CDArborAlg.cpp +++ b/Digitisers/G2CDArbor/src/G2CDArborAlg.cpp @@ -9,18 +9,6 @@ #include "DD4hep/IDDescriptor.h" #include "DD4hep/Plugins.h" -// #include <EVENT/LCCollection.h> -// #include <EVENT/MCParticle.h> -// #include <EVENT/SimCalorimeterHit.h> -// #include <EVENT/CalorimeterHit.h> -// #include <EVENT/LCFloatVec.h> -// #include <EVENT/LCParameters.h> -// #include <IMPL/CalorimeterHitImpl.h> -// #include <IMPL/LCCollectionVec.h> -// #include <IMPL/LCFlagImpl.h> -// #include <IMPL/LCRelationImpl.h> -// #include "UTIL/CellIDDecoder.h" -// #include "DetectorPos.hh" #include <values.h> #include <string> @@ -39,8 +27,6 @@ #include <TVector3.h> using namespace std; -// using namespace lcio ; -// using namespace marlin ; DECLARE_COMPONENT( G2CDArborAlg ) @@ -225,7 +211,7 @@ StatusCode G2CDArborAlg::initialize() { error() << "failed to retrieve dd4hep_geo: " << m_dd4hep_geo << endmsg; return StatusCode::FAILURE; } - + /* // get the DD4hep readout const std::string name_readout = "EcalBarrelCollection"; m_decoder = m_geosvc->getDecoder(name_readout); @@ -233,11 +219,19 @@ StatusCode G2CDArborAlg::initialize() { error() << "Failed to get the decoder. " << endmsg; return StatusCode::FAILURE; } + */ m_encoder_str = "M:3,S-1:3,I:9,J:9,K-1:6"; // printParameters(); // WeightVector.clear(); + for(unsigned int i = 0; i < m_ecalReadoutNames.value().size(); i++){ + m_col_readout_map[m_ecalColNames.value().at(i)] = m_ecalReadoutNames.value().at(i); + //std::cout<<"name="<<m_ecalColNames.value().at(i)<<",readout="<<m_ecalReadoutNames.value().at(i)<<std::endl; + } + for(unsigned int i = 0; i < m_hcalReadoutNames.value().size(); i++){ + m_col_readout_map[m_hcalColNames.value().at(i)] = m_hcalReadoutNames.value().at(i); + } for (auto& ecal : m_ecalColNames) { _ecalCollections.push_back( new SimCaloType(ecal, Gaudi::DataHandle::Reader, this) ); @@ -505,6 +499,13 @@ StatusCode G2CDArborAlg::execute() // for(int k1 = 0; k1 < NumEcalhit; k1++) // { + std::string tmp_readout = m_col_readout_map[m_ecalColNames.value().at(k0)]; + // get the DD4hep readout + m_decoder = m_geosvc->getDecoder(tmp_readout); + if (!m_decoder) { + error() << "Failed to get the decoder. " << endmsg; + return StatusCode::FAILURE; + } edm4hep::CalorimeterHitCollection* ecalcol = _outputEcalCollections[k0]->createAndPut(); auto Ecalcol = _ecalCollections[k0]->get(); @@ -518,9 +519,9 @@ StatusCode G2CDArborAlg::execute() ID_UTIL::CellIDDecoder<edm4hep::SimCalorimeterHit> cellIdDecoder(m_encoder_str); const std::string layerCodingString(m_encoder_str); const std::string layerCoding(this->GetLayerCoding(layerCodingString)); - if(m_readLCIO) LayerNum = m_decoder->get(cellid, "layer"); - else LayerNum = cellIdDecoder(&SimEcalhit)[layerCoding.c_str()]; - // cout << "LayerNum = " << LayerNum << endl; + if(m_readLCIO==false) LayerNum = m_decoder->get(cellid, "layer");//from 0 - 29, 0 is preshower + else LayerNum = cellIdDecoder(&SimEcalhit)[layerCoding.c_str()] + 1 ;//now it is 0 - 29, 0 is preshower + //cout << "LayerNum = " << LayerNum << endl; HitEn = SimEcalhit.getEnergy(); unsigned long long cellID = SimEcalhit.getCellID(); @@ -544,10 +545,12 @@ StatusCode G2CDArborAlg::execute() // _calibCoeffEcal[0] = 48.16; // _calibCoeffEcal[1] = 96.32; - if(LayerNum < _NEcalThinLayer) + //if(LayerNum < _NEcalThinLayer) + if(LayerNum <= _NEcalThinLayer) //layer from 0 - 20 should be thin, total 21 thin layers, _NEcalThinLayer should be 20 DigiHitEn = HitEn * _calibCoeffEcal[0]; else DigiHitEn = HitEn * _calibCoeffEcal[1]; + if( LayerNum==0) DigiHitEn = HitEn; // 0 is preshower layer totalEnergy += DigiHitEn; if(HitEn > _thresholdEcal) diff --git a/Digitisers/G2CDArbor/src/G2CDArborAlg.h b/Digitisers/G2CDArbor/src/G2CDArborAlg.h index 35658c8b01f3e837b61b11f72aa9b86aa6af5d99..1d65b4982cbb8c28e68da6d6bf161561a5ebbbbd 100644 --- a/Digitisers/G2CDArbor/src/G2CDArborAlg.h +++ b/Digitisers/G2CDArbor/src/G2CDArborAlg.h @@ -71,6 +71,9 @@ protected: Gaudi::Property<std::vector<std::string>> m_hcalColNames{this, "HCALCollections", {"HcalBarrelRegCollection", "HcalEndcapsCollection", "HcalEndcapRingsCollection"}, "HCAL Collection Names"}; std::vector<SimCaloType*> _hcalCollections; + Gaudi::Property<std::vector<std::string>> m_ecalReadoutNames{this, "ECALReadOutNames", {"EcalBarrelCollection", "EcalEndcapsCollection","EcalEndcapRingCollection"}, "Name of readouts"}; + Gaudi::Property<std::vector<std::string>> m_hcalReadoutNames{this, "HCALReadOutNames", {"HcalBarrelCollection", "HcalEndcapsCollection","HcalEndcapRingCollection"}, "Name of readouts"}; + std::map<std::string, std::string> m_col_readout_map; // Output collections typedef DataHandle<edm4hep::CalorimeterHitCollection> CaloHitType; diff --git a/Digitisers/SimHitMerge/CMakeLists.txt b/Digitisers/SimHitMerge/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..3a4135239fe9e953662c322c20393e022b26edfa --- /dev/null +++ b/Digitisers/SimHitMerge/CMakeLists.txt @@ -0,0 +1,20 @@ +gaudi_subdir(SimHitMergeAlg v0r0) + +find_package(DD4hep COMPONENTS DDG4 REQUIRED) +find_package(EDM4HEP REQUIRED) +find_package(podio REQUIRED ) +find_package(K4FWCore REQUIRED) + +include_directories(${EDM4HEP_INCLUDE_DIR}) + +gaudi_depends_on_subdirs( + Detector/DetInterface +) +set(SimHitMergeAlg_srcs src/*.cpp) + +# Modules +gaudi_add_module(SimHitMerge ${SimHitMergeAlg_srcs} + INCLUDE_DIRS K4FWCore GaudiKernel GaudiAlgLib DD4hep + LINK_LIBRARIES K4FWCore GaudiKernel GaudiAlgLib DD4hep DDRec + EDM4HEP::edm4hep EDM4HEP::edm4hepDict +) diff --git a/Digitisers/SimHitMerge/src/SimHitMergeAlg.cpp b/Digitisers/SimHitMerge/src/SimHitMergeAlg.cpp new file mode 100644 index 0000000000000000000000000000000000000000..4a1eda3b4ac450b5ea90653ce102b8036acfc9a0 --- /dev/null +++ b/Digitisers/SimHitMerge/src/SimHitMergeAlg.cpp @@ -0,0 +1,118 @@ +#include "SimHitMergeAlg.h" + +#include "edm4hep/SimCalorimeterHit.h" +#include "edm4hep/Vector3f.h" + +#include "DD4hep/Detector.h" +#include "DD4hep/IDDescriptor.h" +#include "DD4hep/Plugins.h" + +#include <string> +#include <iostream> +#include <cmath> +#include <stdexcept> +#include <sstream> + + + +DECLARE_COMPONENT( SimHitMergeAlg ) + + + + + + +SimHitMergeAlg::SimHitMergeAlg(const std::string& name, ISvcLocator* svcLoc) + : GaudiAlgorithm(name, svcLoc) + +{ +} + +StatusCode SimHitMergeAlg::initialize() { + m_geosvc = service<IGeomSvc>("GeomSvc"); + if (!m_geosvc) { + error() << "Failed to find GeomSvc." << endmsg; + return StatusCode::FAILURE; + } + m_dd4hep_geo = m_geosvc->lcdd(); + if (!m_dd4hep_geo) { + error() << "failed to retrieve dd4hep_geo: " << m_dd4hep_geo << endmsg; + return StatusCode::FAILURE; + } + m_cellIDConverter = new dd4hep::rec::CellIDPositionConverter(*m_dd4hep_geo); + + if(m_inputColNames.size() != m_outputColNames.size() ) throw ("Error input size != output size"); + for (auto& input : m_inputColNames) { + m_InputCollections.push_back( new SimCaloType(input, Gaudi::DataHandle::Reader, this) ); + } + + for (auto& output : m_outputColNames) { + m_OutputCollections.push_back( new SimCaloType(output, Gaudi::DataHandle::Writer, this) ); + } + + + return GaudiAlgorithm::initialize(); +} + +StatusCode SimHitMergeAlg::execute() +{ + + for (unsigned int k0 = 0; k0 < m_inputColNames.size(); k0++) + { + std::map<unsigned long long, edm4hep::SimCalorimeterHit> id_hit_map; + std::map<unsigned long long, std::vector<edm4hep::ConstCaloHitContribution> > id_vconb_map; + edm4hep::SimCalorimeterHitCollection* mergedCol = m_OutputCollections[k0]->createAndPut(); + auto col = m_InputCollections[k0]->get(); + //std::cout<<"input="<<m_InputCollections[k0]->objKey()<<",size="<<col->size()<<std::endl; + for (auto Simhit: *col){ + auto id = Simhit.getCellID(); + if(Simhit.getEnergy() <=0 ) continue; + if ( id_hit_map.find(id) != id_hit_map.end()) id_hit_map[id].setEnergy(id_hit_map[id].getEnergy() + Simhit.getEnergy()); + else id_hit_map[id] = Simhit; + std::vector<edm4hep::ConstCaloHitContribution> tmp_vconb ; + for(int kk=0; kk<Simhit.contributions_size(); kk++){ + tmp_vconb.push_back(Simhit.getContributions(kk)); + } + + if( id_vconb_map.find(id) != id_vconb_map.end()){ + for(int kk=0; kk<tmp_vconb.size();kk++){ + id_vconb_map[id].push_back(tmp_vconb.at(kk)); + } + } + else id_vconb_map[id] = tmp_vconb; + + } + + for(std::map<unsigned long long, edm4hep::SimCalorimeterHit>::iterator iter = id_hit_map.begin(); iter != id_hit_map.end(); iter++) + { + edm4hep::SimCalorimeterHit Simhit = iter->second ; + dd4hep::Position position = m_cellIDConverter->position(Simhit.getCellID());//cm + edm4hep::Vector3f hitPos(position.x()*10, position.y()*10, position.z()*10);//to mm + //std::cout<<"id="<<Simhit.getCellID()<<",hitPos.x="<<hitPos[0]<<",y="<<hitPos[1]<<",z="<<hitPos[2]<<",ori x="<<Simhit.getPosition()[0]<<",y="<<Simhit.getPosition()[1]<<",z="<<Simhit.getPosition()[2]<<std::endl; + auto Mergedhit = mergedCol->create(); + Mergedhit.setCellID (Simhit.getCellID()); + Mergedhit.setPosition(hitPos); + Mergedhit.setEnergy (Simhit.getEnergy()); + for(int ii=0; ii<id_vconb_map[iter->first].size(); ii++){ + Mergedhit.addToContributions(id_vconb_map[iter->first].at(ii)); + } + } + /* + std::cout<<"output="<<m_OutputCollections[k0]->objKey()<<",size="<<mergedCol->size()<<std::endl; + for(int ii=0; ii<mergedCol->size(); ii++){ + auto tmp_hit = mergedCol->at(ii); + std::cout<<"id="<<tmp_hit.getCellID()<<",E="<<tmp_hit.getEnergy()<<",conb size="<<tmp_hit.contributions_size()<<std::endl; + } + */ + + } + return StatusCode::SUCCESS; +} + +StatusCode SimHitMergeAlg::finalize() +{ + std::cout<<"SimHitMergeAlg FINISHED"<<std::endl; + + + return GaudiAlgorithm::finalize(); +} diff --git a/Digitisers/SimHitMerge/src/SimHitMergeAlg.h b/Digitisers/SimHitMerge/src/SimHitMergeAlg.h new file mode 100644 index 0000000000000000000000000000000000000000..8bded1c038538fdf346c858b4bdfe55c7796a8db --- /dev/null +++ b/Digitisers/SimHitMerge/src/SimHitMergeAlg.h @@ -0,0 +1,55 @@ +#ifndef SimHitMergeAlg_H +#define SimHitMergeAlg_H + +#include "k4FWCore/DataHandle.h" +#include "GaudiAlg/GaudiAlgorithm.h" +#include "GaudiKernel/Property.h" +#include "edm4hep/EventHeader.h" +#include "edm4hep/EventHeaderCollection.h" +#include "edm4hep/SimCalorimeterHitConst.h" +#include "edm4hep/SimCalorimeterHit.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/MCRecoCaloAssociationCollection.h" +#include "edm4hep/MCParticleCollection.h" + +#include <DDRec/DetectorData.h> +#include <DDRec/CellIDPositionConverter.h> +#include "DetInterface/IGeomSvc.h" + +#include <string> +#include <iostream> +#include <fstream> + + +class SimHitMergeAlg : public GaudiAlgorithm +{ +public: + + SimHitMergeAlg(const std::string& name, ISvcLocator* svcLoc); + + virtual StatusCode initialize() ; + + virtual StatusCode execute() ; + + virtual StatusCode finalize() ; + +protected: + std::string GetLayerCoding(const std::string &encodingString) const; + + + + typedef DataHandle<edm4hep::SimCalorimeterHitCollection> SimCaloType; + + Gaudi::Property<std::vector<std::string>> m_inputColNames{this, "InputCollections", {"EcalBarrelSiliconCollection", "EcalEndcapSiliconCollection", "EcalEndcapRingCollection"}, "Input Hit collection Names"}; + Gaudi::Property<std::vector<std::string>> m_outputColNames{this, "OutputCollections", {"ECALBarrel", "ECALEndcap", "ECALOther"}, "Name of merged Hit Collections"}; + std::vector<SimCaloType*> m_InputCollections; + std::vector<SimCaloType*> m_OutputCollections; + + SmartIF<IGeomSvc> m_geosvc; + dd4hep::Detector* m_dd4hep_geo; + dd4hep::rec::CellIDPositionConverter* m_cellIDConverter; +}; + +#endif + +