diff --git a/Digitization/DigiCalo/src/EcalDigiAlg.cpp b/Digitization/DigiCalo/src/EcalDigiAlg.cpp index 7e0af0cb28602c356206e7a6b5ce5a24504493d3..53e5ecf2d51aa32f24859256bca1a17daaaac259 100644 --- a/Digitization/DigiCalo/src/EcalDigiAlg.cpp +++ b/Digitization/DigiCalo/src/EcalDigiAlg.cpp @@ -321,6 +321,9 @@ StatusCode EcalDigiAlg::execute() totQ1_Truth += en; auto mcp = conb.getParticle(); + if( mcp.parents_size()>0 && (fabs(mcp.getVertex().z)>EcalHalfZ || sqrt(mcp.getVertex().x*mcp.getVertex().x + mcp.getVertex().y*mcp.getVertex().y)>EcalInnerR) ){ + mcp = mcp.getParents(0); + } MCPEnMap[mcp] += en; TVector3 steppos(conb.getStepPosition().x, conb.getStepPosition().y, conb.getStepPosition().z); TVector3 rpos = steppos-hitbar.getPosition(); @@ -576,8 +579,17 @@ StatusCode EcalDigiAlg::execute() rel_MCP2.setSim(iter.first); rel_MCP2.setWeight(iter.second/SimHit.getEnergy()); } - - +/* +cout<<" In SimHit #"<<i<<": linked truth particle "<<MCPEnMap.size()<<": "<<endl; +for(auto iter : MCPEnMap){ + cout<<" Particle "<<iter.first.getPDG()<<", gStatus "<<iter.first.getGeneratorStatus()<<", sStatus "<<iter.first.getSimulatorStatus(); + cout<<", vertex ("<<iter.first.getVertex().x<<", "<<iter.first.getVertex().y<<", "<<iter.first.getVertex().z<<") "; + cout<<", parent size "<<iter.first.parents_size(); + if(iter.first.parents_size()>0) cout<<" parent ID "<<iter.first.getParents(0).getPDG(); + cout<<", weight "<<iter.second/SimHit.getEnergy()<<endl; + cout<<endl; +} +*/ //if(selMCP.isAvailable()){ // auto rel_MCP1 = caloMCPAssoVec->create(); // rel_MCP1.setRec(digiHit1); diff --git a/Digitization/DigiCalo/src/EcalDigiAlg.h b/Digitization/DigiCalo/src/EcalDigiAlg.h index 4a7a8d13aeb34ef442c529b565b7ffdbdfd47b1f..894f0d02a2d146e80b04932f58679197c7badd64 100644 --- a/Digitization/DigiCalo/src/EcalDigiAlg.h +++ b/Digitization/DigiCalo/src/EcalDigiAlg.h @@ -74,6 +74,9 @@ protected: typedef std::vector<float> FloatVec; typedef std::map<const edm4hep::MCParticle, float> MCParticleToEnergyWeightMap; + const float EcalInnerR = 1830.; + const float EcalHalfZ = 2900.; + int _nEvt ; float m_length; TRandom3 rndm; diff --git a/Reconstruction/RecPFACyber/include/CyberDataCol.h b/Reconstruction/RecPFACyber/include/CyberDataCol.h index 1d6854f6239eaa240bb5172984d17679f648b748..161d824ee88e2d303d7c78dc0254620e33a692fa 100644 --- a/Reconstruction/RecPFACyber/include/CyberDataCol.h +++ b/Reconstruction/RecPFACyber/include/CyberDataCol.h @@ -9,8 +9,8 @@ // Contact: guofangyi@ihep.ac.cn, // sunss@ihep.ac.cn //============================================================= -#ifndef _PANDORAPLUS_DATA_H -#define _PANDORAPLUS_DATA_H +#ifndef _CYBER_DATA_H +#define _CYBER_DATA_H #include <iostream> #include <algorithm> #include <map> @@ -32,6 +32,7 @@ #include "edm4hep/MCParticle.h" #include "edm4hep/Track.h" #include "edm4hep/TrackCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" #include "edm4hep/CalorimeterHit.h" #include "edm4hep/CalorimeterHitCollection.h" #include "edm4hep/Vertex.h" @@ -45,9 +46,10 @@ #include "edm4hep/ReconstructedParticleCollection.h" #include "edm4hep/MCRecoCaloAssociation.h" #include "edm4hep/MCRecoTrackerAssociation.h" -#include "edm4hep/MCRecoParticleAssociationCollection.h" -#include "edm4hep/MCRecoCaloParticleAssociationCollection.h" -#include "edm4hep/MCRecoTrackParticleAssociationCollection.h" +#include "edm4hep/MCRecoParticleAssociationCollection.h" // PFO - MCP +#include "edm4hep/MCRecoCaloParticleAssociationCollection.h" // Hit - MCP +#include "edm4hep/MCRecoClusterParticleAssociationCollection.h" // Cluster - MCP +#include "edm4hep/MCRecoTrackParticleAssociationCollection.h" // Track - MCP #define PI 3.141592653 //#define C 299.79 // unit: mm/ns diff --git a/Reconstruction/RecPFACyber/include/CyberPFAlg.h b/Reconstruction/RecPFACyber/include/CyberPFAlg.h index 6f1d195ee7fd38f79ca2e882d4d640249e881638..a3bb35e8fe09f2b205a0d0118f64e571857aa4be 100644 --- a/Reconstruction/RecPFACyber/include/CyberPFAlg.h +++ b/Reconstruction/RecPFACyber/include/CyberPFAlg.h @@ -7,8 +7,8 @@ // Contact: guofangyi@ihep.ac.cn, // sunss@ihep.ac.cn //============================================================= -#ifndef PANDORAPLUS_ALG_H -#define PANDORAPLUS_ALG_H +#ifndef CYBERPFA_ALG_H +#define CYBERPFA_ALG_H #include <string> #include "k4FWCore/DataHandle.h" @@ -136,6 +136,7 @@ protected: //---Readin collections typedef DataHandle<edm4hep::TrackCollection> TrackType; + //typedef DataHandle<edm4hep::SimCalorimeterHitCollection> SimCaloType; typedef DataHandle<edm4hep::CalorimeterHitCollection> CaloType; typedef DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection> CaloParticleAssoType; DataHandle<edm4hep::MCParticleCollection>* r_MCParticleCol; @@ -143,6 +144,7 @@ protected: std::vector<TrackType*> r_TrackCols; //std::vector<CaloType*> r_ECalHitCols; //std::vector<CaloType*> r_HCalHitCols; + //std::vector<SimCaloType*> r_SimCaloHitCols; std::vector<CaloType*> r_CaloHitCols; std::map<std::string, CaloParticleAssoType*> map_CaloMCPAssoCols; DataHandle<edm4hep::RecTofCollection>* r_TofCol; @@ -181,9 +183,13 @@ protected: //Gaudi::Property< std::string > name_EcalCore {this, "OutputEcalCore", "EcalCore"}; //Gaudi::Property< std::string > name_HcalCluster{this, "OutputHcalCluster", "HcalCluster"}; Gaudi::Property< std::string > name_PFObject {this, "OutputPFO", "outputPFO"}; - DataHandle<edm4hep::ReconstructedParticleCollection> w_ReconstructedParticleCollection {"PandoraPFOs" ,Gaudi::DataHandle::Writer, this}; + DataHandle<edm4hep::ReconstructedParticleCollection> w_ReconstructedParticleCollection {"CyberPFO" ,Gaudi::DataHandle::Writer, this}; std::map<std::string, DataHandle<edm4hep::ClusterCollection>* > w_ClusterCollection; + DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection> w_MCPHitAssoCol{"RecCaloHitParticleAssoCol", Gaudi::DataHandle::Writer, this}; + DataHandle<edm4hep::MCRecoClusterParticleAssociationCollection> w_MCPClusterAssoCol{"RecCaloClusterAssoCol", Gaudi::DataHandle::Writer, this}; + DataHandle<edm4hep::MCRecoParticleAssociationCollection> w_MCPRecoPFOAssoCol{"RecPFOAssoCol", Gaudi::DataHandle::Writer, this}; + //For Ana Gaudi::Property<bool> m_WriteAna {this, "WriteAna", false, "Write Ntuples for analysis"}; diff --git a/Reconstruction/RecPFACyber/include/Objects/CaloUnit.h b/Reconstruction/RecPFACyber/include/Objects/CaloUnit.h index 5d4eeb8dfca89e1b9c56e9e7f66ad16721db2e5a..59d17c0133bd669cd516c39ba93da404b146148f 100644 --- a/Reconstruction/RecPFACyber/include/Objects/CaloUnit.h +++ b/Reconstruction/RecPFACyber/include/Objects/CaloUnit.h @@ -88,6 +88,7 @@ namespace Cyber{ static int NbarZ; static float barsize; static float ecal_innerR; + static float ecal_innerZ; static float ecal_endcap_deadarea; static float ecal_endcap_barsize; diff --git a/Reconstruction/RecPFACyber/include/Tools/OutputCreator.h b/Reconstruction/RecPFACyber/include/Tools/OutputCreator.h index 401e7855956ce04c63ec3cc3bd2deac53a222861..c771296b7bd4d9bc4d3a2f10ce7d82a618bbe830 100644 --- a/Reconstruction/RecPFACyber/include/Tools/OutputCreator.h +++ b/Reconstruction/RecPFACyber/include/Tools/OutputCreator.h @@ -16,13 +16,17 @@ namespace Cyber{ OutputCreator( const Settings& m_settings); ~OutputCreator() {}; - StatusCode CreateOutputCollections( CyberDataCol& m_DataCol, - DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecHitsHandler, - DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecCoreHandler, - DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecHcalHitsHandler, - DataHandle<edm4hep::TrackCollection>& m_outTrkHandler, - std::map<std::string, DataHandle<edm4hep::ClusterCollection>*>& m_outClusterColHandler, - DataHandle<edm4hep::ReconstructedParticleCollection>& m_recPFOHandler ); + StatusCode CreateOutputCollections( CyberDataCol& m_DataCol, + std::map<std::string, DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>*>& map_CaloParticleAssoCol, + DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecHitsHandler, + DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecCoreHandler, + DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecHcalHitsHandler, + DataHandle<edm4hep::TrackCollection>& m_outTrkHandler, + std::map<std::string, DataHandle<edm4hep::ClusterCollection>*>& m_outClusterColHandler, + DataHandle<edm4hep::ReconstructedParticleCollection>& m_recPFOHandler, + DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>& m_outMCPHitLinkHandler, + DataHandle<edm4hep::MCRecoClusterParticleAssociationCollection>& m_outMCPClusterLinkHandler, + DataHandle<edm4hep::MCRecoParticleAssociationCollection>& m_outMCPPFOLinkHandler ); StatusCode Reset() { return StatusCode::SUCCESS; } @@ -31,7 +35,7 @@ namespace Cyber{ private: const Cyber::Settings settings; - + }; }; #endif diff --git a/Reconstruction/RecPFACyber/src/Algorithm/HcalClusteringAlg.cpp b/Reconstruction/RecPFACyber/src/Algorithm/HcalClusteringAlg.cpp index e7147853d4521b0c4cbf38486a26ce6cf511812b..bde383c656de0282ffbc38057f2f92c55a08beca 100644 --- a/Reconstruction/RecPFACyber/src/Algorithm/HcalClusteringAlg.cpp +++ b/Reconstruction/RecPFACyber/src/Algorithm/HcalClusteringAlg.cpp @@ -95,6 +95,7 @@ StatusCode HcalClusteringAlg::RunAlgorithm( CyberDataCol& m_datacol ){ //cout<<" Isolated hit size: "<<Nisohits<<endl; // m_datacol.map_CaloCluster[settings.map_stringPars["OutputCluster"]] = m_clusterCol; + for(int ic=0; ic<m_clusterCol.size(); ic++) m_clusterCol[ic]->getLinkedMCPfromHit(); m_datacol.map_CaloCluster[settings.map_stringPars["OutputHCALClusters"]]= m_clusterCol; return StatusCode::SUCCESS; }; diff --git a/Reconstruction/RecPFACyber/src/CyberDataCol.cpp b/Reconstruction/RecPFACyber/src/CyberDataCol.cpp index f8386eab63e3b44fe11d4e9d354a064795cdaac3..3cc18c0b6ba88e96c1a7a1f2f0908dff8c27608d 100644 --- a/Reconstruction/RecPFACyber/src/CyberDataCol.cpp +++ b/Reconstruction/RecPFACyber/src/CyberDataCol.cpp @@ -9,8 +9,8 @@ // Contact: guofangyi@ihep.ac.cn, // sunss@ihep.ac.cn //============================================================= -#ifndef _PANDORAPLUS_DATA_C -#define _PANDORAPLUS_DATA_C +#ifndef _CYBER_DATA_C +#define _CYBER_DATA_C #include "CyberDataCol.h" diff --git a/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp b/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp index 42f5421ca7a296badc5b9922f4b1f5ff0886babc..b3cbcfc1a8533731f8920f3ac5bd7c1b3258a23d 100644 --- a/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp +++ b/Reconstruction/RecPFACyber/src/CyberPFAlg.cpp @@ -8,8 +8,8 @@ // sunss@ihep.ac.cn //============================================================= -#ifndef PANDORAPLUS_ALG_C -#define PANDORAPLUS_ALG_C +#ifndef CYBERPFA_ALG_C +#define CYBERPFA_ALG_C #include "CyberPFAlg.h" @@ -31,6 +31,7 @@ int Cyber::CaloUnit::NbarZ = 24; //int Cyber::CaloUnit::over_module_set = 2; float Cyber::CaloUnit::barsize = 15.2; //mm float Cyber::CaloUnit::ecal_innerR = 1830; //mm +float Cyber::CaloUnit::ecal_innerZ = 2900; //mm float Cyber::CaloUnit::ecal_endcap_deadarea = 10.5; //mm, a bit larger than real value 8.5 mm in geometry float Cyber::CaloUnit::ecal_endcap_barsize = 16.2; //mm, a bit larger than real value 15.2 mm in geometry @@ -48,7 +49,9 @@ CyberPFAlg::CyberPFAlg(const std::string& name, ISvcLocator* svcLoc) declareProperty("RecTrackCollection", w_RecTrkCol, "Handle of Reconstructed Tracks linked to PFO"); declareProperty("RecoPFOCollection", w_ReconstructedParticleCollection, "Handle of Reconstructed PFO collection"); // declareProperty("RecoVtxCollection", w_VertexCollection, "Handle of Reconstructed vertex collection"); - // declareProperty("MCRecoPFOAssociationCollection", w_MCRecoParticleAssociationCollection, "Handle of MC-RecoPFO association collection"); + declareProperty("MCRecoHitAssociationCollection", w_MCPHitAssoCol, "Handle of MC-CaloHit association collection"); + declareProperty("MCRecoClusterAssociationCollection", w_MCPClusterAssoCol, "Handle of MC-CaloCluster association collection"); + declareProperty("MCRecoPFOAssociationCollection", w_MCPRecoPFOAssoCol, "Handle of MC-RecoPFO association collection"); } @@ -90,6 +93,7 @@ StatusCode CyberPFAlg::initialize() m_CaloHitsCreatorSettings.map_stringVecPars["CaloHitCollections"] = name_CaloHits; m_CaloHitsCreatorSettings.map_stringPars["EcalType"] = m_EcalType.value(); + m_OutputCreatorSettings.map_stringVecPars["CaloHitCollections"] = name_CaloHits; m_OutputCreatorSettings.map_stringPars["OutputPFO"] = name_PFObject.value(); m_OutputCreatorSettings.map_boolPars["UseTruthTrk"] = m_useMCPTrk.value(); m_OutputCreatorSettings.map_boolPars["UseTruthMatchTrk"] = m_useTruthMatchTrk.value(); @@ -98,6 +102,9 @@ StatusCode CyberPFAlg::initialize() m_OutputCreatorSettings.map_floatPars["HCALChargedCalib"] = m_HcalChargedCalib.value(); m_OutputCreatorSettings.map_floatPars["ECALNeutralCalib"] = m_EcalNeutralCalib.value(); m_OutputCreatorSettings.map_floatPars["HCALNeutralCalib"] = m_HcalNeutralCalib.value(); + m_OutputCreatorSettings.map_floatPars["EcalHalfZ"] = Cyber::CaloUnit::ecal_innerZ; + m_OutputCreatorSettings.map_floatPars["EcalInnerR"] = Cyber::CaloUnit::ecal_innerR; + m_OutputCreatorSettings.map_floatPars["HitSize"] = Cyber::CaloUnit::barsize; //Initialize Creators m_pMCParticleCreator = new MCParticleCreator( m_pMCParticleCreatorSettings ); @@ -117,6 +124,17 @@ StatusCode CyberPFAlg::initialize() //---Calo Hits--- + //------Sim hits ------ + //for(auto& _ecal : name_SimEcalHits){ + // if(!_ecal.empty()){ + // r_SimCaloHitCols.push_back( new SimCaloType(_ecal, Gaudi::DataHandle::Reader, this) ); + //}} + //for(auto& _hcal : name_SimHcalHits){ + // if(!_hcal.empty()){ + // r_SimCaloHitCols.push_back( new SimCaloType(_hcal, Gaudi::DataHandle::Reader, this) ); + //}} + + //------Digi hits ------ for(auto& _ecal : name_EcalHits){ if(!_ecal.empty()){ //r_ECalHitCols.push_back( new CaloType(_ecal, Gaudi::DataHandle::Reader, this) ); @@ -748,12 +766,17 @@ StatusCode CyberPFAlg::execute() m_algorithmManager.RunAlgorithm( m_DataCol ); m_pOutputCreator->CreateOutputCollections( m_DataCol, + map_CaloMCPAssoCols, w_RecEcalCol, w_RecCoreCol, w_RecHcalCol, w_RecTrkCol, w_ClusterCollection, - w_ReconstructedParticleCollection); + w_ReconstructedParticleCollection, + w_MCPHitAssoCol, + w_MCPClusterAssoCol, + w_MCPRecoPFOAssoCol ); + // yyy_endrec = clock(); // 閲嶅缓缁撴潫鐨勬椂闂� @@ -1634,7 +1657,7 @@ StatusCode CyberPFAlg::execute() for(int ip=0; ip<m_DataCol.map_PFObjects["outputPFO"].size(); ip++) m_pfobjects.push_back(m_DataCol.map_PFObjects["outputPFO"][ip].get()); - +//cout<<"PFO size: "<<m_pfobjects.size()<<endl; for(int ip=0; ip<m_pfobjects.size(); ip++){ std::vector<const Track*> t_tracks = m_pfobjects[ip]->getTracks(); std::vector<const Calo3DCluster*> t_ecal_clusters = m_pfobjects[ip]->getECALClusters(); @@ -1658,7 +1681,7 @@ StatusCode CyberPFAlg::execute() pfo_trk_location.push_back( AllTrackStates[istate].location ); } } - +//cout<<" In PFO #"<<ip<<": Ecal cluster size "<<t_ecal_clusters.size()<<endl; for(int ie=0; ie<t_ecal_clusters.size(); ie++){ pfo_ecal_tag.push_back(ip); pfo_ecal_clus_x.push_back(t_ecal_clusters[ie]->getShowerCenter().x()); @@ -1670,7 +1693,17 @@ StatusCode CyberPFAlg::execute() //if (tmp_phi < 0) tmp_phi += 360.0; //double tmp_theta = std::atan2(t_ecal_clusters[ie]->getShowerCenter().z(), t_ecal_clusters[ie]->getShowerCenter().Perp())* 180.0 / M_PI + 90; pfo_ecal_clus_Escale.push_back(t_ecal_clusters[ie]->getLongiE()); - +/* +double tmp_totE = 0; +cout<<" Cluster #"<<ie<<": energy "<<t_ecal_clusters[ie]->getLongiE()<<", all hitE "<<t_ecal_clusters[ie]->getHitsE()<<", hit size "<<t_ecal_clusters[ie]->getCaloHits().size()<<endl; +for(int ahit=0; ahit<t_ecal_clusters[ie]->getCaloHits().size(); ahit++){ + const Cyber::CaloHit* p_hit = t_ecal_clusters[ie]->getCaloHits()[ahit]; + printf(" Hit #%d: position (%.3f, %.3f, %.3f), En %.3f \n", ahit, p_hit->getPosition().x(), p_hit->getPosition().y(), p_hit->getPosition().z(), p_hit->getEnergy()); + tmp_totE += p_hit->getEnergy(); +} +cout<<" tot Hit E "<<tmp_totE*t_ecal_clusters[ie]->getEnergyScale()<<endl; +cout<<endl; +*/ } for(int ih=0; ih<t_hcal_clusters.size(); ih++){ pfo_hcal_tag.push_back(ip); @@ -1723,15 +1756,13 @@ StatusCode CyberPFAlg::finalize() delete r_MCParticleCol; for(auto iter : r_TrackCols) delete iter; - //for(auto iter : r_ECalHitCols) delete iter; - //for(auto iter : r_HCalHitCols) delete iter; + //for(auto iter : r_SimCaloHitCols) delete iter; for(auto iter : r_CaloHitCols) delete iter; for(auto iter : map_CaloMCPAssoCols) delete iter.second; r_TrackCols.clear(); for(auto iter : w_ClusterCollection) delete iter.second; w_ClusterCollection.clear(); - //r_ECalHitCols.clear(); - //r_HCalHitCols.clear(); + //r_SimCaloHitCols.clear(); r_CaloHitCols.clear(); //m_energycorsvc->finalize(); delete m_cellIDConverter, m_geosvc; diff --git a/Reconstruction/RecPFACyber/src/Objects/Calo3DCluster.cc b/Reconstruction/RecPFACyber/src/Objects/Calo3DCluster.cc index eb642f53d5efb7f7ffd1e246f93ff52b4b8fc747..4bbe660eb17a98101f13de31a3cec77653fca9a1 100644 --- a/Reconstruction/RecPFACyber/src/Objects/Calo3DCluster.cc +++ b/Reconstruction/RecPFACyber/src/Objects/Calo3DCluster.cc @@ -223,7 +223,7 @@ namespace Cyber{ double Calo3DCluster::getEnergy() const{ double result = 0; for(int m=0; m<m_2dclusters.size(); m++) - result += m_2dclusters[m]->getEnergy(); + result += m_2dclusters[m]->getEnergy()*Escale; if(result == 0) result = getLongiE(); return result; } diff --git a/Reconstruction/RecPFACyber/src/Tools/OutputCreator.cpp b/Reconstruction/RecPFACyber/src/Tools/OutputCreator.cpp index 8ae7e8a202e679a806243ef3763ecd60590a113b..dd53b6fe03b6f817ad9e42ff95d98958985b332f 100644 --- a/Reconstruction/RecPFACyber/src/Tools/OutputCreator.cpp +++ b/Reconstruction/RecPFACyber/src/Tools/OutputCreator.cpp @@ -11,12 +11,16 @@ namespace Cyber{ StatusCode OutputCreator::CreateOutputCollections( CyberDataCol& m_DataCol, + std::map<std::string, DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>*>& map_CaloParticleAssoCol, DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecHitsHandler, DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecCoreHandler, DataHandle<edm4hep::CalorimeterHitCollection>& m_outRecHcalHitsHandler, DataHandle<edm4hep::TrackCollection>& m_outTrkHandler, std::map<std::string, DataHandle<edm4hep::ClusterCollection>*>& m_outClusterColHandler, - DataHandle<edm4hep::ReconstructedParticleCollection>& m_recPFOHandler ){ + DataHandle<edm4hep::ReconstructedParticleCollection>& m_recPFOHandler, + DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection>& m_outMCPHitLinkHandler, + DataHandle<edm4hep::MCRecoClusterParticleAssociationCollection>& m_outMCPClusterLinkHandler, + DataHandle<edm4hep::MCRecoParticleAssociationCollection>& m_outMCPPFOLinkHandler ){ //Register the collections edm4hep::CalorimeterHitCollection* m_calohitCol = m_outRecHitsHandler.createAndPut(); @@ -27,11 +31,25 @@ namespace Cyber{ edm4hep::ClusterCollection* m_Ecalcore = m_outClusterColHandler["EcalCore"]->createAndPut(); edm4hep::ClusterCollection* m_Hcalcluster = m_outClusterColHandler["HcalCluster"]->createAndPut(); + //edm4hep::MCRecoCaloParticleAssociationCollection* m_mcphitLinkCol = m_outMCPHitLinkHandler.createAndPut(); + edm4hep::MCRecoClusterParticleAssociationCollection* m_mcpclusterLinkCol = m_outMCPClusterLinkHandler.createAndPut(); + edm4hep::MCRecoParticleAssociationCollection* m_mcppfoLinkCol = m_outMCPPFOLinkHandler.createAndPut(); + + edm4hep::TrackCollection* m_trkCol = nullptr; if(settings.map_boolPars.at("UseTruthTrk") || settings.map_boolPars.at("UseTruthMatchTrk") ) m_trkCol = m_outTrkHandler.createAndPut(); - //PFO + //Readin the MCP-CaloHit link into one collection + std::vector<edm4hep::MCRecoCaloParticleAssociation> m_readinmcphitLinkCol; + for(unsigned int ihit=0; ihit<settings.map_stringVecPars.at("CaloHitCollections").size(); ihit++){ + if( map_CaloParticleAssoCol.find(settings.map_stringVecPars.at("CaloHitCollections")[ihit])!=map_CaloParticleAssoCol.end() ){ + const edm4hep::MCRecoCaloParticleAssociationCollection* const_MCPCaloAssoCol = map_CaloParticleAssoCol[settings.map_stringVecPars.at("CaloHitCollections")[ihit]]->get(); + for(int ilink=0; ilink<const_MCPCaloAssoCol->size(); ilink++) m_readinmcphitLinkCol.push_back(const_MCPCaloAssoCol->at(ilink)); + } + } + + //Create PFO std::vector<std::shared_ptr<Cyber::PFObject>> p_pfos = m_DataCol.map_PFObjects[settings.map_stringPars.at("OutputPFO")]; for(int ip=0; ip<p_pfos.size(); ip++){ @@ -39,6 +57,7 @@ namespace Cyber{ std::vector<const Track*> vec_trks = p_pfos[ip]->getTracks(); std::vector<const Calo3DCluster*> vec_Ecalclus = p_pfos[ip]->getECALClusters(); std::vector<const Calo3DCluster*> vec_Hcalclus = p_pfos[ip]->getHCALClusters(); + std::map<const edm4hep::MCParticle, float> mc_PFOmap; mc_PFOmap.clear(); double ecalcalib = 1.; double hcalcalib = vec_trks.size()==0 ? settings.map_floatPars.at("HCALNeutralCalib") : settings.map_floatPars.at("HCALChargedCalib"); @@ -68,9 +87,10 @@ namespace Cyber{ auto _hit = m_calohitCol->create(); _hit.setCellID(_cellID); - _hit.setEnergy( p_hit->getEnergy()*ecalcalib ); + _hit.setEnergy( p_hit->getEnergy()*p_clus->getEnergyScale()*ecalcalib ); _hit.setPosition( pos ); _hit.setType(1); //Ecal barrel + m_clus.addToHits(_hit); } @@ -107,6 +127,19 @@ namespace Cyber{ m_clus.setEnergy( tmp_clusE ); m_clus.setPosition( pos ); m_clus.addToClusters(m_core); + + std::vector< std::pair<edm4hep::MCParticle, float> > mc_map = p_clus->getLinkedMCP(); + for(auto iter : mc_map){ + auto rel_MCP = m_mcpclusterLinkCol->create(); + rel_MCP.setSim(iter.first); + rel_MCP.setRec(m_clus); + rel_MCP.setWeight( iter.second ); + + if( mc_PFOmap.find(iter.first)==mc_PFOmap.end() ) mc_PFOmap[iter.first] = m_clus.getEnergy()*iter.second; + else mc_PFOmap[iter.first] += m_clus.getEnergy()*iter.second; + + } + m_pfo.addToClusters( m_clus ); EcalClusE += tmp_clusE; @@ -119,7 +152,7 @@ namespace Cyber{ for(int ic=0; ic<vec_Hcalclus.size(); ic++){ auto p_clus = vec_Hcalclus[ic]; auto m_clus = m_Hcalcluster->create(); - + for(int ih=0; ih<p_clus->getCaloHits().size(); ih++){ const Cyber::CaloHit* p_hit = p_clus->getCaloHits()[ih]; //auto _hit = m_hcalhitCol->create(); @@ -132,6 +165,17 @@ namespace Cyber{ m_clus.setEnergy( tmp_clusE ); edm4hep::Vector3f pos( p_clus->getHitCenter().x(), p_clus->getHitCenter().y(), p_clus->getHitCenter().z() ); m_clus.setPosition( pos ); + + std::vector< std::pair<edm4hep::MCParticle, float> > mc_map = p_clus->getLinkedMCP(); + for(auto iter : mc_map){ + auto rel_MCP = m_mcpclusterLinkCol->create(); + rel_MCP.setSim(iter.first); + rel_MCP.setRec(m_clus); + rel_MCP.setWeight( iter.second ); + + if( mc_PFOmap.find(iter.first)==mc_PFOmap.end() ) mc_PFOmap[iter.first] = m_clus.getEnergy()*iter.second; + else mc_PFOmap[iter.first] += m_clus.getEnergy()*iter.second; + } m_pfo.addToClusters( m_clus ); HcalClusE += tmp_clusE; @@ -144,6 +188,13 @@ namespace Cyber{ TVector3 p3vec = vec_Pos*( (EcalClusE+HcalClusE)/vec_Pos.Mag() ); edm4hep::Vector3f p3(p3vec.x(), p3vec.y(), p3vec.z()); + for(auto iter : mc_PFOmap){ + auto rel_MCP = m_mcppfoLinkCol->create(); + rel_MCP.setSim(iter.first); + rel_MCP.setRec(m_pfo); + rel_MCP.setWeight(iter.second/(EcalClusE+HcalClusE)); + } + m_pfo.setEnergy( EcalClusE+HcalClusE ); m_pfo.setCharge( 0 ); m_pfo.setMomentum( p3 ); @@ -170,24 +221,40 @@ namespace Cyber{ m_pfo.setMomentum(p3); m_pfo.setMass(0.139570); //TODO: all charged particles are set to pion mass. m_pfo.setEnergy( sqrt(vec_trks[trkIndex]->getMomentum()*vec_trks[trkIndex]->getMomentum() + 0.139570*0.139570) ); + + std::vector< std::pair<edm4hep::MCParticle, float> > mc_trkmap = vec_trks[trkIndex]->getLinkedMCP(); + + for(auto iter : mc_trkmap){ + auto rel_MCP = m_mcppfoLinkCol->create(); + rel_MCP.setSim(iter.first); + rel_MCP.setRec(m_pfo); + rel_MCP.setWeight(iter.second); + } + } else{ TVector3 p3vec = vec_Pos*( (EcalClusE+HcalClusE)/vec_Pos.Mag() ); edm4hep::Vector3f p3(p3vec.x(), p3vec.y(), p3vec.z()); + for(auto iter : mc_PFOmap){ + auto rel_MCP = m_mcppfoLinkCol->create(); + rel_MCP.setSim(iter.first); + rel_MCP.setRec(m_pfo); + rel_MCP.setWeight(iter.second/(EcalClusE+HcalClusE)); + } + m_pfo.setEnergy( EcalClusE+HcalClusE ); m_pfo.setCharge( 0 ); m_pfo.setMomentum( p3 ); } } -//printf(" Create PFO #%d: charge %.1f, p4 (%.7f, %.7f, %.7f, %.7f), mass %.3f \n", ip, m_pfo.getCharge(), m_pfo.getMomentum().x, m_pfo.getMomentum().y, m_pfo.getMomentum().z, m_pfo.getEnergy(), m_pfo.getMass() ); + //printf(" Create PFO #%d: charge %.1f, p4 (%.7f, %.7f, %.7f, %.7f), mass %.3f \n", ip, m_pfo.getCharge(), m_pfo.getMomentum().x, m_pfo.getMomentum().y, m_pfo.getMomentum().z, m_pfo.getEnergy(), m_pfo.getMass() ); } return StatusCode::SUCCESS; } - edm4hep::Track OutputCreator::TruthTrack(edm4hep::MCParticle _mcp, edm4hep::TrackCollection* _trkCol ){ auto m_track = _trkCol->create(); diff --git a/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp b/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp index f105c6c74aa0f88dc28a84345a293167fde29827..31b1c67642b201b93b81c6602ee72ff9751ec18b 100644 --- a/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp +++ b/Reconstruction/RecPFACyber/src/Tools/TrackCreator.cpp @@ -65,7 +65,7 @@ namespace Cyber{ for(int ilink=0; ilink<const_MCPTrkAssoCol->size(); ilink++){ if( const_TrkCol[itrk] == const_MCPTrkAssoCol->at(ilink).getRec() ) { - m_trk->addLinkedMCP( std::make_pair(const_MCPTrkAssoCol->at(ilink).getSim(), const_MCPTrkAssoCol->at(ilink).getWeight()) ); + m_trk->addLinkedMCP( std::make_pair(const_MCPTrkAssoCol->at(ilink).getSim(), const_MCPTrkAssoCol->at(ilink).getWeight()/(double)const_TrkCol[itrk].trackerHits_size() ) ); break; } }