diff --git a/Examples/options/tut_detsim_pan_matrix.py b/Examples/options/tut_detsim_pan_matrix.py
index d5fabde3075bdba7ef1caae6a4f53ce18e84bea0..4da9e96c384a425e72b5a5b6a4ab9a9cb98a933f 100644
--- a/Examples/options/tut_detsim_pan_matrix.py
+++ b/Examples/options/tut_detsim_pan_matrix.py
@@ -62,15 +62,13 @@ from Configurables import HepMCRdr
from Configurables import GenPrinter
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.ThetaMins = [90] # rad; 45deg
-gun.ThetaMaxs = [90] # rad; 45deg
-gun.PhiMins = [0] # rad; 0deg
-gun.PhiMaxs = [0] # rad; 360deg
+gun.Particles = ["gamma","gamma"]
+gun.Energies = [10, 10] # GeV
+gun.ThetaMins = [90, 90] # degree
+gun.ThetaMaxs = [90, 90] # degree
+gun.PhiMins = [0, 1 ] # degree
+gun.PhiMaxs = [0, 1 ] # degree
+
stdheprdr = StdHepRdr("StdHepRdr")
#stdheprdr.Input = "/cefs/data/stdhep/CEPC250/2fermions/E250.Pbhabha.e0.p0.whizard195/bhabha.e0.p0.00001.stdhep"
@@ -127,6 +125,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 +154,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 +211,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 f8949521ad6b20c76b7bd782d3d52bee1480085f..9af8c82f7cf9ac28ba80a2f5442b1bace5b1ba93 100644
--- a/Generator/src/GtGunTool.cpp
+++ b/Generator/src/GtGunTool.cpp
@@ -15,12 +15,12 @@ GtGunTool::initialize() {
error() << "Please specify the list of particle names/pdgs" << endmsg;
return StatusCode::FAILURE;
}
- /*
+
if (m_energies.value().size() != m_particles.value().size()) {
error() << "Mismatched energies and particles." << endmsg;
return StatusCode::FAILURE;
}
- */
+
// others should be empty or specify
if (m_thetamins.value().size()
&& m_thetamins.value().size() != m_particles.value().size()) {
@@ -79,9 +79,7 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
}
}
- //double energy = m_energies.value()[i];
- double energy_min = m_energies_min.value()[0];
- double energy_max = m_energies_max.value()[0];
+ double energy = m_energies.value()[i];
// create the MC particle
edm4hep::MCParticle mcp = event.m_mc_vec.create();
@@ -95,50 +93,19 @@ GtGunTool::mutate(MyHepMC::GenEvent& event) {
// mcp.setEndpoint();
// assume energy is momentum
- double p = energy_min==energy_max ? energy_max : CLHEP::RandFlat::shoot(energy_min, energy_max);
+ double p = energy;
// 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 theta = m_thetamins.value()[i]==m_thetamaxs.value()[i] ? m_thetamins.value()[i] : CLHEP::RandFlat::shoot(m_thetamins.value()[i], m_thetamaxs.value()[i]);
+ double phi = m_phimins .value()[i]==m_phimaxs .value()[i] ? m_phimins .value()[i] : CLHEP::RandFlat::shoot(m_phimins .value()[i], m_phimaxs .value()[i]);
double costheta = cos(theta*acos(-1)/180);
- phi = phi*acos(-1)/180;
+ double phi_ = phi*acos(-1)/180;
double sintheta = sqrt(1.-costheta*costheta);
-
- // check if theta min/max is set
- /*
- if (i < m_thetamins.value().size()
- && i < m_thetamaxs.value().size()) {
- double thetamin = m_thetamins.value()[i];
- double thetamax = m_thetamaxs.value()[i];
-
- if (thetamin == thetamax) { // fixed theta
- costheta = cos(thetamin);
- sintheta = sin(thetamin);
- info() << "theta is fixed: " << thetamin << endmsg;
- }
- }
-
- if (i < m_phimins.value().size()
- && i < m_phimaxs.value().size()) {
- double phimin = m_phimins.value()[i];
- double phimax = m_phimaxs.value()[i];
-
- if (phimin == phimax) { // fixed phi
- phi = phimin;
- info() << "phi is fixed: " << phimin << endmsg;
- }
- }
-
- debug() << "Direction: "
- << " cos(theta): " << costheta
- << " 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();
diff --git a/Generator/src/GtGunTool.h b/Generator/src/GtGunTool.h
index 55161ca890ab29dfb3f4d91fa8cc3099e5b07d51..22b02a811a19be2ce5a16ca1cdd4726635a0be4b 100644
--- a/Generator/src/GtGunTool.h
+++ b/Generator/src/GtGunTool.h
@@ -35,8 +35,6 @@ private:
Gaudi::Property<std::vector<std::string>> m_particles{this, "Particles"};
Gaudi::Property<std::vector<double>> m_energies{this, "Energies"};
- Gaudi::Property<std::vector<double>> m_energies_min{this, "EnergyMin"};
- Gaudi::Property<std::vector<double>> m_energies_max{this, "EnergyMax"};
Gaudi::Property<std::vector<double>> m_thetamins{this, "ThetaMins"};
Gaudi::Property<std::vector<double>> m_thetamaxs{this, "ThetaMaxs"};
diff --git a/Reconstruction/Digi_Calo/CMakeLists.txt b/Reconstruction/Digi_Calo/CMakeLists.txt
index 5bec00f6f3dcf67d7d9a9c2ea1971328da055440..507720ca40cd718c485744ef5e95e68c0c674bd8 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 723e0e2117d6e4f683e8e203298f8a78838312d5..93e2fb35c0fa217382f588947fb6ce7db60d0f90 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 050ffe7b6b8932ee181f9a9555a1e118fc9cea5b..c1a82c8867450cf7335dce74d4743d0e7e9b374e 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 96a6938b06f963b32f4d4a240681080d4a3a7dfa..63210760fdcfe3b7b4e08697ce340639dee6edf3 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 e1d2e19986cb3d2b6da84407af94f474f543f427..566f61e5b116ba626db23dba934b495e3f2e8b8d 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;
}
+