From 0e14bd1da5f36063858fc4cd5499e3d0025a33f3 Mon Sep 17 00:00:00 2001
From: Fang Wenxing <wxfang@lxslc613.ihep.ac.cn>
Date: Mon, 22 Jun 2020 10:08:38 +0800
Subject: [PATCH] update MC reco relation construction
---
Examples/options/tut_detsim_pan_matrix.py | 15 +++++--
Generator/src/GtGunTool.cpp | 29 ++++++++++---
Generator/src/GtGunTool.h | 4 ++
Reconstruction/Digi_Calo/CMakeLists.txt | 7 +---
Reconstruction/Digi_Calo/src/CaloDigiAlg.cpp | 26 +++++++++++-
Reconstruction/Digi_Calo/src/CaloDigiAlg.h | 2 +
.../MatrixPandora/include/PandoraMatrixAlg.h | 3 ++
.../MatrixPandora/src/PandoraMatrixAlg.cpp | 42 +++++++++++++++----
8 files changed, 104 insertions(+), 24 deletions(-)
diff --git a/Examples/options/tut_detsim_pan_matrix.py b/Examples/options/tut_detsim_pan_matrix.py
index d5fabde3..08c80175 100644
--- a/Examples/options/tut_detsim_pan_matrix.py
+++ b/Examples/options/tut_detsim_pan_matrix.py
@@ -65,13 +65,18 @@ gun = GtGunTool("GtGunTool")
gun.Particles = ["gamma"]
#gun.Energies = [0.5, 1] # GeV
#gun.EnergyMin = [0.1] # GeV
-gun.EnergyMin = [1] # GeV
-gun.EnergyMax = [1] # GeV
+gun.EnergyMin = [10] # GeV
+gun.EnergyMax = [10] # GeV
gun.ThetaMins = [90] # rad; 45deg
gun.ThetaMaxs = [90] # rad; 45deg
gun.PhiMins = [0] # rad; 0deg
gun.PhiMaxs = [0] # rad; 360deg
+gun.Create_second = True
+gun.Energy1 = 10
+gun.dtheta = 0
+gun.dphi = 1
+
stdheprdr = StdHepRdr("StdHepRdr")
#stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep"
#stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep"
@@ -127,6 +132,7 @@ example_CaloDigiAlg = CaloDigiAlg("CaloDigiAlg")
example_CaloDigiAlg.Scale = 1
example_CaloDigiAlg.SimCaloHitCollection = "SimCalorimeterCol"
example_CaloDigiAlg.CaloHitCollection = "ECALBarrel"
+example_CaloDigiAlg.CaloAssociationCollection = "RecoCaloAssociation_ECALBarrel"
##############################################################################
from Configurables import GearSvc
gearSvc = GearSvc("GearSvc")
@@ -155,10 +161,11 @@ pandoralg.ReadProngVertices = "ProngVertices"
pandoralg.ReadSplitVertices = "SplitVertices"
pandoralg.ReadV0Vertices = "V0Vertices"
pandoralg.ReadTracks = "MarlinTrkTracks"
+pandoralg.MCRecoCaloAssociation = "RecoCaloAssociation_ECALBarrel"
pandoralg.WriteClusterCollection = "PandoraClusters"
pandoralg.WriteReconstructedParticleCollection = "PandoraPFOs"
pandoralg.WriteVertexCollection = "PandoraPFANewStartVertices"
-pandoralg.AnaOutput = "/cefs/higgs/wxfang/cepc/Pandora/Ana/gamma/Ana_gamma_Matrix_NewTest.root"
+pandoralg.AnaOutput = "/cefs/higgs/wxfang/cepc/Pandora/Ana/gamma/Ana_gamma_Matrix_Rel_10GeV_test.root"
pandoralg.PandoraSettingsDefault_xml = "/junofs/users/wxfang/MyGit/MarlinPandora/scripts/PandoraSettingsDefault_wx.xml"
#### Do not chage the collection name, only add or delete ###############
@@ -211,7 +218,7 @@ ApplicationMgr(
#TopAlg = [genalg, detsimalg, write],
#TopAlg = [read, pandoralg],
EvtSel = 'NONE',
- EvtMax = 100,
+ EvtMax = 10,
#ExtSvc = [rndmengine, dsvc, geosvc, gearSvc],
ExtSvc = [rndmengine, dsvc, geosvc, gearSvc,detsimsvc],
OutputLevel=INFO
diff --git a/Generator/src/GtGunTool.cpp b/Generator/src/GtGunTool.cpp
index f8949521..68d33c1b 100644
--- a/Generator/src/GtGunTool.cpp
+++ b/Generator/src/GtGunTool.cpp
@@ -96,15 +96,20 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
// assume energy is momentum
double p = energy_min==energy_max ? energy_max : CLHEP::RandFlat::shoot(energy_min, energy_max);
+ double p1 = m_p1.value() ;
// direction
// by default, randomize the direction
double theta = m_thetamins.value()[0]==m_thetamaxs.value()[0] ? m_thetamins.value()[0] : CLHEP::RandFlat::shoot(m_thetamins.value()[0], m_thetamaxs.value()[0]);
double phi = m_phimins .value()[0]==m_phimaxs .value()[0] ? m_phimins .value()[0] : CLHEP::RandFlat::shoot(m_phimins .value()[0], m_phimaxs .value()[0]);
+ double theta1 = m_dtheta.value() + theta;
+ double phi1 = m_dphi.value() + phi ;
double costheta = cos(theta*acos(-1)/180);
- phi = phi*acos(-1)/180;
+ double costheta1 = cos(theta1*acos(-1)/180);
+ double phi_ = phi*acos(-1)/180;
+ double phi1_ = phi1*acos(-1)/180;
double sintheta = sqrt(1.-costheta*costheta);
-
+ double sintheta1 = sqrt(1.-costheta1*costheta1);
// check if theta min/max is set
/*
if (i < m_thetamins.value().size()
@@ -135,14 +140,28 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
<< " phi: " << phi
<< endmsg;
*/
- double px = p*sintheta*cos(phi);
- double py = p*sintheta*sin(phi);
+ double px = p*sintheta*cos(phi_);
+ double py = p*sintheta*sin(phi_);
double pz = p*costheta;
- std::cout<<"GenGt p="<<p<<", px="<<px<<",py="<<py<<",pz="<<pz<<",costheta="<<costheta<<std::endl;
+ std::cout<<"GenGt p="<<p<<", px="<<px<<",py="<<py<<",pz="<<pz<<",theta="<<theta<<",phi="<<phi<<std::endl;
mcp.setMomentum(edm4hep::Vector3f(px,py,pz));
// mcp.setMomentumAtEndpoint();
// mcp.setSpin();
// mcp.setColorFlow();
+ if(m_create_second)
+ {
+ edm4hep::MCParticle mcp1 = event.m_mc_vec.create();
+ mcp1.setPDG(pdgcode);
+ mcp1.setGeneratorStatus(1);
+ mcp1.setSimulatorStatus(1);
+ mcp1.setTime(0.0);
+ mcp1.setMass(mass);
+ double px1 = p1*sintheta1*cos(phi1_);
+ double py1 = p1*sintheta1*sin(phi1_);
+ double pz1 = p1*costheta1;
+ std::cout<<"GenGt p1="<<p1<<", px1="<<px1<<",py1="<<py1<<",pz1="<<pz1<<",theta1="<<theta1<<",phi1="<<phi1<<std::endl;
+ mcp1.setMomentum(edm4hep::Vector3f(px1,py1,pz1));
+ }
}
diff --git a/Generator/src/GtGunTool.h b/Generator/src/GtGunTool.h
index 55161ca8..1549d089 100644
--- a/Generator/src/GtGunTool.h
+++ b/Generator/src/GtGunTool.h
@@ -44,6 +44,10 @@ private:
Gaudi::Property<std::vector<double>> m_phimins{this, "PhiMins"};
Gaudi::Property<std::vector<double>> m_phimaxs{this, "PhiMaxs"};
+ Gaudi::Property<bool> m_create_second{this, "Create_second", false};
+ Gaudi::Property< double > m_p1{this, "Energy1",1};
+ Gaudi::Property< double > m_dtheta{this, "dtheta",1};
+ Gaudi::Property< double > m_dphi {this, "dphi" ,1};
};
diff --git a/Reconstruction/Digi_Calo/CMakeLists.txt b/Reconstruction/Digi_Calo/CMakeLists.txt
index 5bec00f6..507720ca 100644
--- a/Reconstruction/Digi_Calo/CMakeLists.txt
+++ b/Reconstruction/Digi_Calo/CMakeLists.txt
@@ -5,7 +5,6 @@ find_package(EDM4HEP REQUIRED )
message("EDM4HEP_INCLUDE_DIRS: ${EDM4HEP_INCLUDE_DIR}")
message("EDM4HEP_LIB: ${EDM4HEP_LIBRARIES}")
include_directories(${EDM4HEP_INCLUDE_DIR})
-link_libraries("/cvmfs/cepcsw.ihep.ac.cn/prototype/releases/externals/97.0.0/EDM4hep/lib64/libedm4hep.so")
find_package(CLHEP REQUIRED)
find_package(podio REQUIRED )
@@ -21,8 +20,6 @@ gaudi_depends_on_subdirs(
gaudi_add_module(Digi_Calo ${srcs}
INCLUDE_DIRS FWCore GaudiKernel GaudiAlgLib CLHEP DD4hep
LINK_LIBRARIES FWCore GaudiKernel GaudiAlgLib CLHEP DD4hep ${DD4hep_COMPONENT_LIBRARIES} DDRec
+ -Wl,--no-as-needed
+ EDM4HEP::edm4hep EDM4HEP::edm4hepDict
)
-#gaudi_add_module(Digi_Calo ${srcs}
-# INCLUDE_DIRS FWCore GaudiKernel GaudiAlgLib CLHEP
-# LINK_LIBRARIES FWCore GaudiKernel GaudiAlgLib CLHEP
-#)
diff --git a/Reconstruction/Digi_Calo/src/CaloDigiAlg.cpp b/Reconstruction/Digi_Calo/src/CaloDigiAlg.cpp
index 723e0e21..93e2fb35 100644
--- a/Reconstruction/Digi_Calo/src/CaloDigiAlg.cpp
+++ b/Reconstruction/Digi_Calo/src/CaloDigiAlg.cpp
@@ -27,6 +27,7 @@ CaloDigiAlg::CaloDigiAlg(const std::string& name, ISvcLocator* svcLoc)
// Output collections
declareProperty("CaloHitCollection", w_DigiCaloCol, "Handle of Digi CaloHit collection");
+ declareProperty("CaloAssociationCollection", w_CaloAssociationCol, "Handle of CaloAssociation collection");
}
@@ -45,8 +46,9 @@ StatusCode CaloDigiAlg::initialize()
StatusCode CaloDigiAlg::execute()
{
std::map<unsigned long long, edm4hep::SimCalorimeterHit> id_hit_map;
- std::map<unsigned long long, int > test_map;
+ std::map<unsigned long long, std::vector<edm4hep::SimCalorimeterHit> > id_hits_map;
edm4hep::CalorimeterHitCollection* caloVec = w_DigiCaloCol.createAndPut();
+ edm4hep::MCRecoCaloAssociationCollection* caloAssoVec = w_CaloAssociationCol.createAndPut();
const edm4hep::SimCalorimeterHitCollection* SimHitCol = r_SimCaloCol.get();
double tot_e = 0 ;
if(SimHitCol == 0)
@@ -61,9 +63,16 @@ StatusCode CaloDigiAlg::execute()
unsigned long long id = SimHit.getCellID();
float en = SimHit.getEnergy();
tot_e += en;
- test_map[id] = 1;
if ( id_hit_map.find(id) != id_hit_map.end()) id_hit_map[id].setEnergy(id_hit_map[id].getEnergy() + en);
else id_hit_map[id] = SimHit ;
+
+ if ( id_hits_map.find(id) != id_hits_map.end()) id_hits_map[id].push_back(SimHit);
+ else
+ {
+ std::vector<edm4hep::SimCalorimeterHit> vhit;
+ vhit.push_back(SimHit);
+ id_hits_map[id] = vhit ;
+ }
}
for(std::map<unsigned long long, edm4hep::SimCalorimeterHit>::iterator iter = id_hit_map.begin(); iter != id_hit_map.end(); iter++)
{
@@ -74,10 +83,23 @@ StatusCode CaloDigiAlg::execute()
edm4hep::Vector3f vpos(position.x()*10, position.y()*10, position.z()*10);// cm to mm
caloHit.setPosition(vpos);
//std::cout << "sim hit id =" << caloHit.getCellID() <<",x="<<position.x()<<",y="<<position.y()<<",z="<<position.z() <<",real x="<<(iter->second).getPosition().x <<",y="<<(iter->second).getPosition().y<<",z="<<(iter->second).getPosition().z<< std::endl;
+
+ if( id_hits_map.find(iter->first) != id_hits_map.end())
+ {
+ for(unsigned int i=0; i< id_hits_map[iter->first].size(); i++)
+ {
+ auto asso = caloAssoVec->create();
+ asso.setRec(caloHit);
+ asso.setSim(id_hits_map[iter->first].at(i));
+ asso.setWeight(id_hits_map[iter->first].at(i).getEnergy()/(iter->second).getEnergy());
+ }
+ }
+ else std::cout<<"Error in Digi Calo"<<std::endl;
}
std::cout<<"total sim e ="<< tot_e <<std::endl;
std::cout<<"digi, output digi hit size="<< caloVec->size() <<std::endl;
+ std::cout<<"digi, output caloAssoVec hit size="<< caloAssoVec->size() <<std::endl;
_nEvt ++ ;
return StatusCode::SUCCESS;
diff --git a/Reconstruction/Digi_Calo/src/CaloDigiAlg.h b/Reconstruction/Digi_Calo/src/CaloDigiAlg.h
index 050ffe7b..c1a82c88 100644
--- a/Reconstruction/Digi_Calo/src/CaloDigiAlg.h
+++ b/Reconstruction/Digi_Calo/src/CaloDigiAlg.h
@@ -8,6 +8,7 @@
#include "edm4hep/CalorimeterHit.h"
#include "edm4hep/CalorimeterHitCollection.h"
#include "edm4hep/SimCalorimeterHitCollection.h"
+#include "edm4hep/MCRecoCaloAssociationCollection.h"
#include <DDRec/DetectorData.h>
#include <DDRec/CellIDPositionConverter.h>
@@ -51,6 +52,7 @@ protected:
DataHandle<edm4hep::SimCalorimeterHitCollection> r_SimCaloCol{"SimCaloCol", Gaudi::DataHandle::Reader, this};
// Output collections
DataHandle<edm4hep::CalorimeterHitCollection> w_DigiCaloCol{"DigiCaloCol", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::MCRecoCaloAssociationCollection> w_CaloAssociationCol{"MCRecoCaloAssociationCollection", Gaudi::DataHandle::Writer, this};
};
#endif
diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/include/PandoraMatrixAlg.h b/Reconstruction/PFA/Pandora/MatrixPandora/include/PandoraMatrixAlg.h
index 96a6938b..63210760 100644
--- a/Reconstruction/PFA/Pandora/MatrixPandora/include/PandoraMatrixAlg.h
+++ b/Reconstruction/PFA/Pandora/MatrixPandora/include/PandoraMatrixAlg.h
@@ -270,6 +270,8 @@ protected:
//### For Ana #################
TFile* m_fout;
TTree* m_tree;
+ std::vector< std::vector<int> > m_pReco_mc_id;
+ std::vector< std::vector<float> > m_pReco_mc_weight;
std::vector<int > m_pReco_PID;
std::vector<float> m_pReco_mass;
std::vector<float> m_pReco_energy;
@@ -278,6 +280,7 @@ protected:
std::vector<float> m_pReco_pz;
std::vector<float> m_pReco_charge;
+ std::vector<int> m_mc_id;
std::vector<int> m_mc_p_size;
std::vector<int> m_mc_pid ;
std::vector<float> m_mc_mass ;
diff --git a/Reconstruction/PFA/Pandora/MatrixPandora/src/PandoraMatrixAlg.cpp b/Reconstruction/PFA/Pandora/MatrixPandora/src/PandoraMatrixAlg.cpp
index e1d2e199..566f61e5 100644
--- a/Reconstruction/PFA/Pandora/MatrixPandora/src/PandoraMatrixAlg.cpp
+++ b/Reconstruction/PFA/Pandora/MatrixPandora/src/PandoraMatrixAlg.cpp
@@ -15,6 +15,7 @@
#include "LCContent.h"
+#include "TInterpreter.h"
pandora::Pandora* PandoraMatrixAlg::m_pPandora=0;
@@ -103,12 +104,17 @@ void PandoraMatrixAlg::FinaliseSteeringParameters(ISvcLocator* svcloc)
StatusCode PandoraMatrixAlg::initialize()
{
+ gInterpreter->GenerateDictionary("vector<vector<float> >", "vector");
+ gInterpreter->GenerateDictionary("vector<vector<int> >", "vector");
+
std::cout<<"hi init PandoraMatrixAlg"<<std::endl;
std::string s_output =m_AnaOutput;
m_fout = new TFile(s_output.c_str(),"RECREATE");
m_tree = new TTree("evt","tree");
m_tree->Branch("m_pReco_PID" , &m_pReco_PID);
+ m_tree->Branch("m_pReco_mc_id" , &m_pReco_mc_id);
+ m_tree->Branch("m_pReco_mc_weight", &m_pReco_mc_weight);
m_tree->Branch("m_pReco_mass" , &m_pReco_mass);
m_tree->Branch("m_pReco_energy", &m_pReco_energy);
m_tree->Branch("m_pReco_px" , &m_pReco_px);
@@ -116,6 +122,7 @@ StatusCode PandoraMatrixAlg::initialize()
m_tree->Branch("m_pReco_pz" , &m_pReco_pz);
m_tree->Branch("m_pReco_charge", &m_pReco_charge);
+ m_tree->Branch("m_mc_id" , &m_mc_id);
m_tree->Branch("m_mc_p_size", &m_mc_p_size);
m_tree->Branch("m_mc_pid" , &m_mc_pid );
m_tree->Branch("m_mc_mass" , &m_mc_mass );
@@ -312,7 +319,7 @@ StatusCode PandoraMatrixAlg::execute()
updateMap();
PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pMCParticleCreator->CreateMCParticles(*m_CollectionMaps));
PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pCaloHitCreator->CreateCaloHits(*m_CollectionMaps));
- //PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pMCParticleCreator->CreateCaloHitToMCParticleRelationships(*m_CollectionMaps, m_pCaloHitCreator->GetCalorimeterHitVector() ));
+ PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pMCParticleCreator->CreateCaloHitToMCParticleRelationships(*m_CollectionMaps, m_pCaloHitCreator->GetCalorimeterHitVector() ));
//PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pTrackCreator->CreateTrackAssociations(*m_CollectionMaps));
//PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pTrackCreator->CreateTracks(*m_CollectionMaps));
//PANDORA_THROW_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, m_pMCParticleCreator->CreateTrackToMCParticleRelationships(*m_CollectionMaps, m_pTrackCreator->GetTrackVector() ));
@@ -374,10 +381,6 @@ pandora::StatusCode PandoraMatrixAlg::RegisterUserComponents() const
"NonLinearity", pandora::HADRONIC, m_settings.m_inputEnergyCorrectionPoints, m_settings.m_outputEnergyCorrectionPoints));
- /*
- PANDORA_RETURN_RESULT_IF(pandora::STATUS_CODE_SUCCESS, !=, PandoraApi::RegisterAlgorithmFactory(*m_pPandora,
- "ExternalClustering", new ExternalClusteringAlgorithm::Factory));
- */
return pandora::STATUS_CODE_SUCCESS;
}
@@ -390,6 +393,8 @@ void PandoraMatrixAlg::Reset()
m_pMCParticleCreator->Reset();
std::vector<int>() .swap(m_pReco_PID );
+ std::vector< std::vector<int> >() .swap(m_pReco_mc_id);
+ std::vector< std::vector<float> >().swap(m_pReco_mc_weight);
std::vector<float>().swap(m_pReco_mass);
std::vector<float>().swap(m_pReco_energy);
std::vector<float>().swap(m_pReco_px);
@@ -397,6 +402,7 @@ void PandoraMatrixAlg::Reset()
std::vector<float>().swap(m_pReco_pz);
std::vector<float>().swap(m_pReco_charge);
+ std::vector<int>() .swap(m_mc_id);
std::vector<int>() .swap(m_mc_p_size);
std::vector<int>() .swap(m_mc_pid );
std::vector<float>().swap(m_mc_mass );
@@ -600,8 +606,8 @@ StatusCode PandoraMatrixAlg::updateMap()
if (NULL != mcRecoCaloAssociation )
{
std::vector<edm4hep::MCRecoCaloAssociation> v_cal;
- m_CollectionMaps->collectionMap_CaloRel["RecoCaloAssociation"] = v_cal ;
- for(unsigned int i=0 ; i< mcRecoCaloAssociation->size(); i++) m_CollectionMaps->collectionMap_CaloRel ["RecoCaloAssociation"].push_back(mcRecoCaloAssociation->at(i));
+ m_CollectionMaps->collectionMap_CaloRel["RecoCaloAssociation_ECALBarrel"] = v_cal ;
+ for(unsigned int i=0 ; i< mcRecoCaloAssociation->size(); i++) m_CollectionMaps->collectionMap_CaloRel ["RecoCaloAssociation_ECALBarrel"].push_back(mcRecoCaloAssociation->at(i));
}
else
@@ -841,11 +847,17 @@ StatusCode PandoraMatrixAlg::Ana()
m_pReco_px.push_back(px);
m_pReco_py.push_back(py);
m_pReco_pz.push_back(pz);
+ std::vector<int> tmp_mc_id;
+ std::vector<float> tmp_mc_weight;
for(int j=0; j < reco_associa_col->size(); j++)
{
if(reco_associa_col->at(j).getRec().id() != pReco.id() ) continue;
- std::cout<<"MC pid ="<<reco_associa_col->at(j).getSim().getPDG()<<", px="<<reco_associa_col->at(j).getSim().getMomentum()[0]<<", py="<<reco_associa_col->at(j).getSim().getMomentum()[1]<<",pz="<<reco_associa_col->at(j).getSim().getMomentum()[2]<<std::endl;
+ std::cout<<"MC pid ="<<reco_associa_col->at(j).getSim().getPDG()<<",weight="<<reco_associa_col->at(j).getWeight()<<", px="<<reco_associa_col->at(j).getSim().getMomentum()[0]<<", py="<<reco_associa_col->at(j).getSim().getMomentum()[1]<<",pz="<<reco_associa_col->at(j).getSim().getMomentum()[2]<<std::endl;
+ tmp_mc_id .push_back(reco_associa_col->at(j).getSim().id());
+ tmp_mc_weight.push_back(reco_associa_col->at(j).getWeight());
}
+ m_pReco_mc_id .push_back(tmp_mc_id);
+ m_pReco_mc_weight.push_back(tmp_mc_weight);
}
const edm4hep::MCParticleCollection* MCParticle = nullptr;
StatusCode sc = StatusCode::SUCCESS;
@@ -854,6 +866,7 @@ StatusCode PandoraMatrixAlg::Ana()
{
for(unsigned int i=0 ; i< MCParticle->size(); i++)
{
+ m_mc_id .push_back(MCParticle->at(i).id());
m_mc_p_size.push_back(MCParticle->at(i).parents_size());
m_mc_pid .push_back(MCParticle->at(i).getPDG());
m_mc_mass .push_back(MCParticle->at(i).getMass());
@@ -901,6 +914,8 @@ StatusCode PandoraMatrixAlg::CreateMCRecoParticleAssociation()
for(int i=0; i<reco_col->size();i++)
{
std::map<int, edm4hep::ConstMCParticle> mc_map;
+ std::map<int, float > id_edep_map;
+ float tot_en = 0 ;
const edm4hep::ReconstructedParticle pReco = reco_col->at(i);
for(int j=0; j < pReco.clusters_size(); j++)
{
@@ -916,6 +931,9 @@ StatusCode PandoraMatrixAlg::CreateMCRecoParticleAssociation()
for(std::vector<edm4hep::ConstCaloHitContribution>::const_iterator itc = it->getSim().contributions_begin(); itc != it->getSim().contributions_end(); itc++)
{
if(mc_map.find(itc->getParticle().id()) == mc_map.end()) mc_map[itc->getParticle().id()] = itc->getParticle() ;
+ if(id_edep_map.find(itc->getParticle().id()) != id_edep_map.end()) id_edep_map[itc->getParticle().id()] = id_edep_map[itc->getParticle().id()] + itc->getEnergy() ;
+ else id_edep_map[itc->getParticle().id()] = itc->getEnergy() ;
+ tot_en += itc->getEnergy() ;
}
}
}
@@ -926,7 +944,15 @@ StatusCode PandoraMatrixAlg::CreateMCRecoParticleAssociation()
edm4hep::MCRecoParticleAssociation association = pMCRecoParticleAssociationCollection->create();
association.setRec(pReco);
association.setSim(it->second);
+ if(tot_en==0)
+ {
+ association.setWeight(0);
+ std::cout<<"Found 0 cluster energy"<<std::endl;
+ }
+ else if(id_edep_map.find(it->first) != id_edep_map.end()) association.setWeight(id_edep_map[it->first]/tot_en);
+ else std::cout<<"Error in creating MCRecoParticleAssociation"<<std::endl;
}
}
return StatusCode::SUCCESS;
}
+
--
GitLab