diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
index 2e76ba5520f7da020a6358e27be57336a12e7a60..1cc1894f7f1d0b33ad38a1123060b2c3902d2ee7 100644
--- a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
+++ b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
@@ -223,6 +223,14 @@ tpr.TrackList = ["CompleteTracks", "ClupatraTracks"]
tpr.TrackerAssociationList = ["VXDTrackerHitAssociation", "SITTrackerHitAssociation", "SETTrackerHitAssociation", "FTDTrackerHitAssociation", "TPCTrackerHitAss"]
#tpr.OutputLevel = DEBUG
+from Configurables import TrueMuonTagAlg
+tmt = TrueMuonTagAlg("TrueMuonTag")
+tmt.MCParticleCollection = "MCParticle"
+tmt.TrackList = ["CompleteTracks"]
+tmt.MuonTagEfficiency = 0.95 # muon true tag efficiency, default is 1.0 (100%)
+tmt.MuonDetTanTheta = 1.2 # muon det barrel/endcap separation tan(theta)
+#tmt.OutputLevel = DEBUG
+
# output
from Configurables import PodioOutput
out = PodioOutput("outputalg")
@@ -232,7 +240,7 @@ out.outputCommands = ["keep *"]
# ApplicationMgr
from Configurables import ApplicationMgr
mgr = ApplicationMgr(
- TopAlg = [podioinput, digiVXD, digiSIT, digiSET, digiFTD, digiTPC, tracking, forward, subset, clupatra, full, tpr, tpc_dndx, out],
+ TopAlg = [podioinput, digiVXD, digiSIT, digiSET, digiFTD, digiTPC, tracking, forward, subset, clupatra, full, tpr, tpc_dndx, tmt, out],
EvtSel = 'NONE',
EvtMax = 5,
ExtSvc = [rndmengine, rndmgensvc, dsvc, evtseeder, geosvc, gearsvc, tracksystemsvc, pidsvc],
diff --git a/Reconstruction/RecAssociationMaker/CMakeLists.txt b/Reconstruction/RecAssociationMaker/CMakeLists.txt
index d854356b8be3c4ed9f5c0913c2e4836bdab14cb7..56187b36b9fdf96ae4373c7461b52e0c2356ee35 100644
--- a/Reconstruction/RecAssociationMaker/CMakeLists.txt
+++ b/Reconstruction/RecAssociationMaker/CMakeLists.txt
@@ -1,6 +1,7 @@
# Modules
gaudi_add_module(RecAssociationMaker
SOURCES src/TrackParticleRelationAlg.cpp
+ src/TrueMuonTagAlg.cpp
LINK GearSvc
EventSeeder
diff --git a/Reconstruction/RecAssociationMaker/README_TrueMuTag.md b/Reconstruction/RecAssociationMaker/README_TrueMuTag.md
new file mode 100644
index 0000000000000000000000000000000000000000..e403f4189cb973197f352480c54195460ea307c8
--- /dev/null
+++ b/Reconstruction/RecAssociationMaker/README_TrueMuTag.md
@@ -0,0 +1,70 @@
+
+
+This is about the temporary `TrueMuonTagAlg`.
+
+The `TrueMuonTagAlg` is implemented in:
+
+`
+Reconstruction/RecAssociationMaker/src/TrueMuonTagAlg.h|cpp
+`
+
+refering to the `TrackParticleRelationAlg` implemented in the same folder.
+
+It first matches the `FullLDCTrack` with a `MCParticle`.
+If the `MCParticle` is a muon, according to a true efficiency `MuonTagEfficiency`,
+randomly flag it as having a muon Tag or not,
+separately for Barrel and Endcap.
+
+The flag is stored in the `edm4hep::MutableTrack::type()`,
+following the `lcio::ILDDetID::` numbers, which is defined in:
+[LCIO-02-20/src/cpp/src/UTIL/ILDConf.cc](https://github.com/iLCSoft/LCIO/blob/v02-20/src/cpp/src/UTIL/ILDConf.cc)
+
+
+An example configure is shown in `tracking_trueMuonTag.py`, e.g.:
+
+```
+from Configurables import TrueMuonTagAlg
+tmt = TrueMuonTagAlg("TrueMuonTag")
+tmt.MCParticleCollection = "MCParticle"
+tmt.TrackList = ["CompleteTracks"]
+tmt.MuonTagEfficiency = 0.95 # muon true tag efficiency, default is 1.0 (100%)
+tmt.MuonDetTanTheta = 1.2 # muon det barrel/endcap separation tan(theta)
+#tmt.OutputLevel = DEBUG
+```
+
+In this example, it adds the muon Tag in the `edm4hep::MutableTrack::type()` for each track in
+"CompleteTracks" collection (which is a `FullLDCTrack` collection)
+and output it to be "CompleteTracksWithMuonTag" in your output `lcio` root files.
+
+To use it in your analysis codes, one can do the following, either in your cpp codes in CEPCSW,
+or in your `ROOT` analysis codes.
+
+Example 1, for CEPCSW Algorithms:
+
+Add the following in your header file:
+```
+#include <UTIL/ILDConf.h>
+```
+
+Do the following in your cpp file, suppose `track` is the track you get from the
+"CompleteTracksWithMuonTag" collection:
+```
+bool hasMuonTag_MuonBarrel = ( ( track.type() >> lcio:ILDDetID:YOKE ) % 2 == 1 );
+bool hasMuonTag_MuonEndcap = ( ( track.type() >> lcio:ILDDetID:YOKE\_ENDCAP ) % 2 == 1 );
+```
+
+Example 2, for ROOT analysis .C macros, suppose `events` is the `TTree` in your output root
+file storing the events.
+
+```
+events.Draw("((CompleteTracksWithMuonTag.type>>27) % 2)>>hHasMuonTag_MuonBarrel");
+events.Draw("((CompleteTracksWithMuonTag.type>>31) % 2)>>hHasMuonTag_MuonEndcap");
+```
+
+This will draw one histogram for has muon tag in barrel or in endcap.
+in case a track has muon tag, it is one entry at 1 in this histogram.
+And the numbers 27 and 31 are the numbers defined for `lcio:ILDDetID:YOKE`
+and `lcio:ILDDetID:YOKE_ENDCAP` in `UTIL/ILDConf.h`, respectively.
+
+
+
diff --git a/Reconstruction/RecAssociationMaker/src/TrueMuonTagAlg.cpp b/Reconstruction/RecAssociationMaker/src/TrueMuonTagAlg.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3ae934bf52788178b928400020be37f6cbb612ae
--- /dev/null
+++ b/Reconstruction/RecAssociationMaker/src/TrueMuonTagAlg.cpp
@@ -0,0 +1,183 @@
+#include "TrueMuonTagAlg.h"
+
+DECLARE_COMPONENT(TrueMuonTagAlg)
+
+TrueMuonTagAlg::TrueMuonTagAlg(const std::string& name, ISvcLocator* svcLoc)
+: GaudiAlgorithm(name, svcLoc){
+ declareProperty("MCParticleCollection", m_inMCParticleColHdl, "Handle of the Input MCParticle collection");
+ declareProperty("MuonTagEfficiency", _m_muonTagEff, "Muon Tagging efficiency to be used to mimic the true muon tag, default is 1. (100%)");
+ declareProperty("MuonDetTanTheta", _m_muonDetTanTheta, "Muon barrel/endcap separate theta angle.");
+}
+
+StatusCode TrueMuonTagAlg::initialize() {
+ for(auto name : m_inTrackCollectionNames) {
+ m_inTrackColHdls.push_back(new DataHandle<edm4hep::TrackCollection> (name, Gaudi::DataHandle::Reader, this));
+ //m_outColHdls.push_back(new DataHandle<edm4hep::MCRecoTrackParticleAssociationCollection> (name+"ParticleAssociation", Gaudi::DataHandle::Writer, this));
+ m_outTrackColHdls.push_back(new DataHandle<edm4hep::TrackCollection> (name+"WithMuonTag", Gaudi::DataHandle::Writer, this));
+ }
+
+ for(unsigned i=0; i<m_inAssociationCollectionNames.size(); i++) {
+ m_inAssociationColHdls.push_back(new DataHandle<edm4hep::MCRecoTrackerAssociationCollection> (m_inAssociationCollectionNames[i], Gaudi::DataHandle::Reader, this));
+ }
+
+ if(m_inAssociationColHdls.size()==0) {
+ warning() << "no Association Collection input" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ m_nEvt = 0;
+ return GaudiAlgorithm::initialize();
+}
+
+StatusCode TrueMuonTagAlg::execute() {
+ info() << "start Event " << m_nEvt << endmsg;
+
+ // MCParticle
+ const edm4hep::MCParticleCollection* mcpCol = nullptr;
+ try {
+ mcpCol = m_inMCParticleColHdl.get();
+ for (auto mcp : *mcpCol) {
+ debug() << "id=" << mcp.id() << " pdgid=" << mcp.getPDG() << " charge=" << mcp.getCharge() << " genstat=" << mcp.getGeneratorStatus() << " simstat=" << mcp.getSimulatorStatus()
+ << " p=" << mcp.getMomentum() << endmsg;
+ int nparents = mcp.parents_size();
+ int ndaughters = mcp.daughters_size();
+ for (int i=0; i<nparents; i++) {
+ debug() << "<<<<<<" << mcp.getParents(i).id() << endmsg;
+ }
+ for (int i=0; i<ndaughters; i++) {
+ debug() << "<<<<<<" << mcp.getDaughters(i).id() << endmsg;
+ }
+ }
+ }
+ catch ( GaudiException &e ) {
+ debug() << "Collection " << m_inMCParticleColHdl.fullKey() << " is unavailable in event " << m_nEvt << endmsg;
+ }
+ if(!mcpCol) {
+ warning() << "pass this event because MCParticle collection can not read" << endmsg;
+ return StatusCode::SUCCESS;
+ }
+
+ // Prepare map from hit to MCParticle
+ std::map<edm4hep::TrackerHit, edm4hep::MCParticle> mapHitParticle;
+ debug() << "reading Association" << endmsg;
+ for (auto hdl : m_inAssociationColHdls) {
+ const edm4hep::MCRecoTrackerAssociationCollection* assCol = nullptr;
+ try {
+ assCol = hdl->get();
+ }
+ catch ( GaudiException &e ) {
+ debug() << "Collection " << hdl->fullKey() << " is unavailable in event " << m_nEvt << endmsg;
+ return StatusCode::FAILURE;
+ }
+ for (auto ass: *assCol) {
+ auto hit = ass.getRec();
+ auto particle = ass.getSim().getMCParticle();
+ mapHitParticle[hit] = particle;
+ }
+ }
+
+ // Handle all input TrackCollection
+ for (unsigned icol=0; icol<m_inTrackColHdls.size(); icol++) {
+ auto hdl = m_inTrackColHdls[icol];
+
+ const edm4hep::TrackCollection* trkCol = nullptr;
+ try {
+ trkCol = hdl->get();
+ }
+ catch ( GaudiException &e ) {
+ debug() << "Collection " << hdl->fullKey() << " is unavailable in event " << m_nEvt << endmsg;
+ }
+
+ //auto outCol = m_outColHdls[icol]->createAndPut();
+ //if(!outCol) {
+ // error() << "Collection " << m_outColHdls[icol]->fullKey() << " cannot be created" << endmsg;
+ // return StatusCode::FAILURE;
+ //}
+
+ auto outTrackCol = m_outTrackColHdls[icol]->createAndPut();
+ if(!outTrackCol) {
+ error() << "Collection " << m_outTrackColHdls[icol]->fullKey() << " cannot be created" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ if(trkCol) {
+ std::map<edm4hep::MCParticle, std::vector<edm4hep::TrackerHit> > mapParticleHits;
+
+ for (auto track: *trkCol) {
+ std::map<edm4hep::MCParticle, int> mapParticleNHits;
+
+ // Count hits related to MCParticle
+ int nhits = track.trackerHits_size();
+ for (int ihit=0; ihit<nhits; ihit++) {
+ auto hit = track.getTrackerHits(ihit);
+ auto itHP = mapHitParticle.find(hit);
+ if (itHP==mapHitParticle.end()) {
+ warning() << "track " << track.id() << "'s hit " << hit.id() << "not in association list, will be ignored" << endmsg;
+ }
+ else {
+ auto particle = itHP->second;
+ if(!particle.isAvailable()) continue;
+ debug() << "track " << track.id() << "'s hit " << hit.id() << " link to MCParticle " << particle.id() << endmsg;
+ auto itPN = mapParticleNHits.find(particle);
+ if (itPN!=mapParticleNHits.end()) itPN->second++;
+ else mapParticleNHits[particle] = 1;
+
+ if (msgLevel(MSG::DEBUG)) mapParticleHits[particle].push_back(hit);
+ }
+ }
+ debug() << "Total " << mapParticleNHits.size() << " particles link to the track " << track.id() << endmsg;
+
+ // Save to collection
+ //for (std::map<edm4hep::MCParticle, int>::iterator it=mapParticleNHits.begin(); it!=mapParticleNHits.end(); it++) {
+ // auto ass = outCol->create();
+ // ass.setWeight(it->second);
+ // ass.setRec(track);
+ // ass.setSim(it->first);
+ //}
+
+ // tore to new track collection
+ edm4hep::MutableTrack new_track = track.clone();
+
+ // find the max nhits mcparticle
+ auto max_nhits_iter = std::max_element(mapParticleNHits.begin(), mapParticleNHits.end(),
+ [](const auto& lhs, const auto& rhs) {return lhs.second < rhs.second;} );
+ // calculate the barrel/endcap theta
+ auto a_particle = max_nhits_iter->first;
+ auto p_pdg_id = a_particle.getPDG();
+ auto p_status = a_particle.getGeneratorStatus();
+ auto p_px = a_particle.getMomentum()[0];
+ auto p_py = a_particle.getMomentum()[1];
+ auto p_pz = a_particle.getMomentum()[2];
+ auto p_tantheta = sqrt(p_px*p_px+p_py*p_py)/fabs(p_pz);
+
+
+ if (abs(p_pdg_id)==13 && fRandom.Rndm()<_m_muonTagEff ) {
+ debug() << " particle PDG = " << p_pdg_id << "; status = " << p_status << endmsg;
+ if (p_tantheta>_m_muonDetTanTheta) new_track.setType( new_track.getType()| (1<<lcio::ILDDetID::YOKE) ) ;
+ else new_track.setType( new_track.getType()| (1<<lcio::ILDDetID::YOKE_ENDCAP) ) ;
+ }
+ outTrackCol->push_back(new_track);
+
+ }
+
+ if (msgLevel(MSG::DEBUG)) {
+ for (std::map<edm4hep::MCParticle, std::vector<edm4hep::TrackerHit> >::iterator it=mapParticleHits.begin(); it!=mapParticleHits.end(); it++) {
+ auto particle = it->first;
+ auto hits = it->second;
+ debug() << "=== MCPaticle ===" << particle << endmsg;
+ for (auto hit : hits) {
+ debug() << hit << endmsg;
+ }
+ }
+ }
+ }
+ }
+
+ m_nEvt++;
+ return StatusCode::SUCCESS;
+}
+
+StatusCode TrueMuonTagAlg::finalize() {
+
+ return StatusCode::SUCCESS;
+}
diff --git a/Reconstruction/RecAssociationMaker/src/TrueMuonTagAlg.h b/Reconstruction/RecAssociationMaker/src/TrueMuonTagAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..e61faf1e04aa45a596e30ac635e9ae34a4071af5
--- /dev/null
+++ b/Reconstruction/RecAssociationMaker/src/TrueMuonTagAlg.h
@@ -0,0 +1,52 @@
+#ifndef TrueMuonTagAlg_h
+#define TrueMuonTagAlg_h 1
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+
+#include "edm4hep/MCParticleCollection.h"
+#include "edm4hep/TrackCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/MCRecoTrackParticleAssociationCollection.h"
+
+#include <UTIL/BitField64.h>
+#include <UTIL/ILDConf.h>
+
+#include "TRandom3.h"
+
+class TrueMuonTagAlg : public GaudiAlgorithm {
+ public:
+ // Constructor of this form must be provided
+ TrueMuonTagAlg( const std::string& name, ISvcLocator* pSvcLocator );
+
+ // Three mandatory member functions of any algorithm
+ StatusCode initialize() override;
+ StatusCode execute() override;
+ StatusCode finalize() override;
+
+ private:
+ // input MCParticle
+ DataHandle<edm4hep::MCParticleCollection> m_inMCParticleColHdl{"MCParticle", Gaudi::DataHandle::Reader, this};
+ // input Tracks to make relation
+ std::vector<DataHandle<edm4hep::TrackCollection>* > m_inTrackColHdls;
+ Gaudi::Property<std::vector<std::string> > m_inTrackCollectionNames{this, "TrackList", {"SiTracks"}};
+ // input TrackerAssociation to link TrackerHit and SimTrackerHit
+ std::vector<DataHandle<edm4hep::MCRecoTrackerAssociationCollection>* > m_inAssociationColHdls;
+ Gaudi::Property<std::vector<std::string> > m_inAssociationCollectionNames{this, "TrackerAssociationList", {"VXDTrackerHitAssociation",
+ "SITTrackerHitAssociation", "SETTrackerHitAssociation", "FTDTrackerHitAssociation"}};
+
+ // output TrackParticleAssociation
+ //std::vector<DataHandle<edm4hep::MCRecoTrackParticleAssociationCollection>* > m_outColHdls;
+ std::vector<DataHandle<edm4hep::TrackCollection>* > m_outTrackColHdls;
+
+ // muon tag efficiency
+ Gaudi::Property<double> _m_muonTagEff{this, "MuonTagEfficiency", 1.0};
+ // Muon barrel/endcap separate angle
+ Gaudi::Property<double> _m_muonDetTanTheta{this,"MuonDetTanTheta", 1.2};
+
+ int m_nEvt;
+
+ TRandom3 fRandom{1234};
+
+};
+#endif