From 2d501cab8eab168f50158cc586c5fd2896ea7b08 Mon Sep 17 00:00:00 2001 From: Fangyi Guo Date: Wed, 27 Dec 2023 14:53:53 +0800 Subject: [PATCH 1/2] Add ECAL & HCAL digi and truth info in track --- .gitignore | 1 + Digitisers/CMakeLists.txt | 1 + Digitisers/CRDEcalDigi/CMakeLists.txt | 22 + Digitisers/CRDEcalDigi/src/CRDEcalDigiAlg.cpp | 505 ++++++++++++++++++ Digitisers/CRDEcalDigi/src/CRDEcalDigiAlg.h | 116 ++++ Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.cpp | 214 ++++++++ Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.h | 100 ++++ Digitisers/CRDEcalDigi/src/CaloBar.h | 56 ++ Digitisers/CRDEcalDigi/src/CaloHit.h | 50 ++ Digitisers/CRDEcalDigi/src/HitStep.h | 26 + Examples/options/detsim_FullDet_HCAL.py | 315 +++++++++++ .../FullLDCTracking/FullLDCTrackingAlg.cpp | 161 ++++-- .../src/FullLDCTracking/FullLDCTrackingAlg.h | 11 + 13 files changed, 1533 insertions(+), 45 deletions(-) create mode 100644 Digitisers/CRDEcalDigi/CMakeLists.txt create mode 100644 Digitisers/CRDEcalDigi/src/CRDEcalDigiAlg.cpp create mode 100644 Digitisers/CRDEcalDigi/src/CRDEcalDigiAlg.h create mode 100644 Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.cpp create mode 100644 Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.h create mode 100644 Digitisers/CRDEcalDigi/src/CaloBar.h create mode 100644 Digitisers/CRDEcalDigi/src/CaloHit.h create mode 100644 Digitisers/CRDEcalDigi/src/HitStep.h create mode 100644 Examples/options/detsim_FullDet_HCAL.py diff --git a/.gitignore b/.gitignore index 58e07a1f2..23b5389d2 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ build spack* ./Generator/output/ ./Generator/options/ +InstallArea/ diff --git a/Digitisers/CMakeLists.txt b/Digitisers/CMakeLists.txt index 6f6f2e404..44b4b03ba 100644 --- a/Digitisers/CMakeLists.txt +++ b/Digitisers/CMakeLists.txt @@ -2,3 +2,4 @@ add_subdirectory(DCHDigi) add_subdirectory(G2CDArbor) add_subdirectory(SimHitMerge) add_subdirectory(SimpleDigi) +add_subdirectory(CRDEcalDigi) diff --git a/Digitisers/CRDEcalDigi/CMakeLists.txt b/Digitisers/CRDEcalDigi/CMakeLists.txt new file mode 100644 index 000000000..4a44558e2 --- /dev/null +++ b/Digitisers/CRDEcalDigi/CMakeLists.txt @@ -0,0 +1,22 @@ +# Modules +gaudi_add_module(CRDEcalDigi + SOURCES src/CRDEcalDigiAlg.cpp + SOURCES src/CRDHcalDigiAlg.cpp + LINK k4FWCore::k4FWCore + GearSvc + DetInterface + Gaudi::GaudiKernel + Gaudi::GaudiAlgLib + ${CLHEP_LIBRARIES} + ${GEAR_LIBRARIES} + ${GSL_LIBRARIES} + ${LCIO_LIBRARIES} + ${ROOT_LIBRARIES} + EDM4HEP::edm4hep EDM4HEP::edm4hepDict +) +install(TARGETS CRDEcalDigi + EXPORT CEPCSWTargets + RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin + LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib + COMPONENT dev) + diff --git a/Digitisers/CRDEcalDigi/src/CRDEcalDigiAlg.cpp b/Digitisers/CRDEcalDigi/src/CRDEcalDigiAlg.cpp new file mode 100644 index 000000000..f52c1153b --- /dev/null +++ b/Digitisers/CRDEcalDigi/src/CRDEcalDigiAlg.cpp @@ -0,0 +1,505 @@ +/* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */ +// Unit in code: mm, ns. +// NOTE: This digitialization highly matches detector geometry CRDEcalBarrel_v01. +// TODO: read geometry info automatically. +#include "CRDEcalDigiAlg.h" + +#include "edm4hep/SimCalorimeterHit.h" +#include "edm4hep/CalorimeterHit.h" +#include "edm4hep/Vector3f.h" +#include "edm4hep/Cluster.h" + +#include "DD4hep/Detector.h" +#include +#include + +#include "TVector3.h" +#include +#include +#include +#include +#include + +// #include +// #include + +#define C 299.79 // unit: mm/ns +#define PI 3.141592653 +using namespace std; +using namespace dd4hep; + +DECLARE_COMPONENT( CRDEcalDigiAlg ) + +CRDEcalDigiAlg::CRDEcalDigiAlg(const std::string& name, ISvcLocator* svcLoc) + : GaudiAlgorithm(name, svcLoc), + _nEvt(0) +{ + + // Input collections + declareProperty("SimCaloHitCollection", r_SimCaloCol, "Handle of the Input SimCaloHit collection"); + + // Output collections + declareProperty("CaloHitCollection", w_DigiCaloCol, "Handle of Digi CaloHit collection"); + declareProperty("CaloAssociationCollection", w_CaloAssociationCol, "Handle of CaloAssociation collection"); + declareProperty("CaloMCPAssociationCollection", w_MCPCaloAssociationCol, "Handle of CaloAssociation collection"); + +} + +StatusCode CRDEcalDigiAlg::initialize() +{ + + if(_writeNtuple){ + std::string s_outfile = _filename; + m_wfile = new TFile(s_outfile.c_str(), "recreate"); + t_SimCont = new TTree("SimStep", "SimStep"); + t_SimBar = new TTree("SimBarHit", "SimBarHit"); + t_SimCont->Branch("step_x", &m_step_x); + t_SimCont->Branch("step_y", &m_step_y); + t_SimCont->Branch("step_z", &m_step_z); + t_SimCont->Branch("step_t", &m_step_t); // yyy: time of each step + t_SimCont->Branch("stepBar_x", &m_stepBar_x); + t_SimCont->Branch("stepBar_y", &m_stepBar_y); + t_SimCont->Branch("stepBar_z", &m_stepBar_z); + t_SimCont->Branch("step_E", &m_step_E); + t_SimCont->Branch("step_T1", &m_step_T1); + t_SimCont->Branch("step_T2", &m_step_T2); + t_SimBar->Branch("totE", &totE); + t_SimBar->Branch("simBar_x", &m_simBar_x); + t_SimBar->Branch("simBar_y", &m_simBar_y); + t_SimBar->Branch("simBar_z", &m_simBar_z); + t_SimBar->Branch("simBar_T1", &m_simBar_T1); + t_SimBar->Branch("simBar_T2", &m_simBar_T2); + t_SimBar->Branch("simBar_Q1", &m_simBar_Q1); + t_SimBar->Branch("simBar_Q2", &m_simBar_Q2); + t_SimBar->Branch("simBar_module", &m_simBar_module); + t_SimBar->Branch("simBar_stave", &m_simBar_stave); + t_SimBar->Branch("simBar_dlayer", &m_simBar_dlayer); + t_SimBar->Branch("simBar_part", &m_simBar_part); + t_SimBar->Branch("simBar_slayer", &m_simBar_slayer); + t_SimBar->Branch("simBar_cellID", &m_simBar_cellID); + } + + std::cout<<"CRDEcalDigiAlg::m_scale="<("GeomSvc"); + if ( !m_geosvc ) throw "CRDEcalDigiAlg :Failed to find GeomSvc ..."; + dd4hep::Detector* m_dd4hep = m_geosvc->lcdd(); + if ( !m_dd4hep ) throw "CRDEcalDigiAlg :Failed to get dd4hep::Detector ..."; + m_cellIDConverter = new dd4hep::rec::CellIDPositionConverter(*m_dd4hep); + m_decoder = m_geosvc->getDecoder(_readoutName); + if (!m_decoder) { + error() << "Failed to get the decoder. " << endmsg; + return StatusCode::FAILURE; + } + + //m_edmsvc = service("CRDEcalSvc"); + //if ( !m_edmsvc ) throw "CRDEcalDigiAlg :Failed to find CRDEcalSvc ..."; + + rndm.SetSeed(_seed); + std::cout<<"CRDEcalDigiAlg::initialize"< m_simhitCol; m_simhitCol.clear(); + std::vector m_barCol; m_barCol.clear(); + + if(SimHitCol == 0) + { + std::cout<<"not found SimCalorimeterHitCollection"<< std::endl; + return StatusCode::SUCCESS; + } + if(_Debug>=1) std::cout<<"digi, input sim hit size="<< SimHitCol->size() <=1) std::cout<<"Finish Hit Merge, with Nhit: "<get(id, "system"), + m_decoder->get(id, "module"), + m_decoder->get(id, "stave"), + m_decoder->get(id, "dlayer"), + m_decoder->get(id, "part"), + m_decoder->get(id, "slayer"), + m_decoder->get(id, "bar")); + + double Lbar = GetBarLength(hitbar); //NOTE: Is fixed with geometry CRDEcalBarrel_v01. + dd4hep::Position hitpos = m_cellIDConverter->position(id); + TVector3 barpos(10*hitpos.x(), 10*hitpos.y(), 10*hitpos.z()); //cm to mm. + hitbar.setPosition(barpos); + //hitbar.T1 = 99999; hitbar.T2 = 99999; + //if(_Debug>=2) std::cout<<"SimHit contribution size: "< DigiLvec; DigiLvec.clear(); + std::vector DigiRvec; DigiRvec.clear(); + double totQ1 = 0; + double totQ2 = 0; + + //Loop in all SimHitContribution(G4Step). + for(int iCont=0; iCont < SimHit.contributions_size(); ++iCont){ + auto conb = SimHit.getContributions(iCont); + if( !conb.isAvailable() ) { std::cout<<"CRDEcalDigiAlg Can not get SimHitContribution: "<=3){ + cout<<"Cell Pos: "<=3){ + cout<500){ std::cout<<"ERROR: Step "<500) continue; + double Ti_right = -1; looptime=0; + while(Ti_right<0){ + // Ti_right = Tinit + rndm.Gaus(nMat*(Lbar/2 - sign*sqrt(rpos.Mag2()))/C, Tres); + Ti_right = Tinit + rndm.Gaus(nMat*(Lbar/2 - sign*sqrt(rpos.Mag2()))/C, Tres) + step_time; // yyy: add step time + looptime++; + if(looptime>500){ std::cout<<"ERROR: Step "<500) continue; + + + m_step_T1.push_back(Ti_left); + m_step_T2.push_back(Ti_right); + totQ1 += Qi_left; + totQ2 += Qi_right; + + HitStep stepoutL, stepoutR; + stepoutL.setQ(Qi_left); stepoutL.setT(Ti_left); + stepoutR.setQ(Qi_right); stepoutR.setT(Ti_right); + DigiLvec.push_back(stepoutL); + DigiRvec.push_back(stepoutR); + + } + + //Time digitalization + //if(_Debug>=2) std::cout<<"Time Digitalize: time at Q >"<<_Qthfrac<<"*totQ"<totQ1*_Qthfrac){ + thT1 = DigiLvec[iCont].getT(); + if(_Debug>=3) std::cout<<"Get T1 at index: "<totQ2*_Qthfrac){ + thT2 = DigiRvec[iCont].getT(); + if(_Debug>=3) std::cout<<"Get T2 at index: "<create(); + digiHit1.setCellID(hitbar.getcellID()); + digiHit1.setEnergy(hitbar.getQ1()); + digiHit1.setTime(hitbar.getT1()); + digiHit1.setPosition(m_pos); + auto digiHit2 = caloVec->create(); + digiHit2.setCellID(hitbar.getcellID()); + digiHit2.setEnergy(hitbar.getQ2()); + digiHit2.setTime(hitbar.getT2()); + digiHit2.setPosition(m_pos); + + //SimHit - CaloHit association + auto rel1 = caloAssoVec->create(); + rel1.setRec(digiHit1); + rel1.setSim(SimHit); + rel1.setWeight( hitbar.getQ1()/(hitbar.getQ1()+hitbar.getQ2()) ); + auto rel2 = caloAssoVec->create(); + rel2.setRec(digiHit2); + rel2.setSim(SimHit); + rel2.setWeight( hitbar.getQ2()/(hitbar.getQ1()+hitbar.getQ2()) ); + + + //MCParticle - CaloHit association + //float maxMCE = -99.; + //edm4hep::MCParticle selMCP; + for(auto iter : MCPEnMap){ + //if(iter.second>maxMCE){ + // maxMCE = iter.second; + // selMCP = iter.first; + //} + auto rel_MCP1 = caloMCPAssoVec->create(); + rel_MCP1.setRec(digiHit1); + rel_MCP1.setSim(iter.first); + rel_MCP1.setWeight(iter.second/digiHit1.getEnergy()); + auto rel_MCP2 = caloMCPAssoVec->create(); + rel_MCP2.setRec(digiHit2); + rel_MCP2.setSim(iter.first); + rel_MCP2.setWeight(iter.second/digiHit2.getEnergy()); + } + + //if(selMCP.isAvailable()){ + // auto rel_MCP1 = caloMCPAssoVec->create(); + // rel_MCP1.setRec(digiHit1); + // rel_MCP1.setSim(selMCP); + // rel_MCP1.setWeight(1.); + // auto rel_MCP2 = caloMCPAssoVec->create(); + // rel_MCP2.setRec(digiHit2); + // rel_MCP2.setSim(selMCP); + // rel_MCP2.setWeight(1.); + //} + + m_barCol.push_back(hitbar); + totE+=(hitbar.getQ1()+hitbar.getQ2())/2; + + //Temp: write into trees. + if(_writeNtuple){ + m_simBar_x.push_back(hitbar.getPosition().x()); + m_simBar_y.push_back(hitbar.getPosition().y()); + m_simBar_z.push_back(hitbar.getPosition().z()); + m_simBar_Q1.push_back(hitbar.getQ1()); + m_simBar_Q2.push_back(hitbar.getQ2()); + m_simBar_T1.push_back(hitbar.getT1()); + m_simBar_T2.push_back(hitbar.getT2()); + m_simBar_module.push_back(hitbar.getModule()); + m_simBar_stave.push_back(hitbar.getStave()); + m_simBar_dlayer.push_back(hitbar.getDlayer()); + m_simBar_part.push_back(hitbar.getPart()); + m_simBar_slayer.push_back(hitbar.getSlayer()); + m_simBar_cellID.push_back(hitbar.getcellID()); + } + } + + + if(_writeNtuple){ + t_SimCont->Fill(); + t_SimBar->Fill(); + } + if(_Debug>=1) std::cout<<"End Loop: Bar Digitalization!"<cd(); + t_SimCont->Write(); + t_SimBar->Write(); + m_wfile->Close(); + delete m_wfile, t_SimCont, t_SimBar; + } + + info() << "Processed " << _nEvt << " events " << endmsg; + delete m_cellIDConverter, m_decoder, m_geosvc; + return GaudiAlgorithm::finalize(); +} + +StatusCode CRDEcalDigiAlg::MergeHits( const edm4hep::SimCalorimeterHitCollection& m_col, std::vector& m_hits ){ + + m_hits.clear(); + std::vector m_mergedhit; + m_mergedhit.clear(); + + for(int iter=0; iterposition(cellid); + edm4hep::Vector3f pos(hitpos.x()*10, hitpos.y()*10, hitpos.z()*10); + + edm4hep::MutableCaloHitContribution conb; + conb.setEnergy(m_step.getEnergy()); + conb.setStepPosition(m_step.getPosition()); + conb.setParticle( m_step.getContributions(0).getParticle() ); + conb.setTime(m_step.getContributions(0).getTime()); + + edm4hep::MutableSimCalorimeterHit m_hit = find(m_mergedhit, cellid); + if(m_hit.getCellID()==0){ + //m_hit = new edm4hep::SimCalorimeterHit(); + m_hit.setCellID(cellid); + m_hit.setPosition(pos); + m_mergedhit.push_back(m_hit); + } + m_hit.addToContributions(conb); + m_hit.setEnergy(m_hit.getEnergy()+m_step.getEnergy()); + } + + for(auto iter = m_mergedhit.begin(); iter!=m_mergedhit.end(); iter++){ + edm4hep::SimCalorimeterHit constsimhit = *iter; + m_hits.push_back( constsimhit ); + } + return StatusCode::SUCCESS; +} + + + +double CRDEcalDigiAlg::GetBarLength(CaloBar& bar){ + //TODO: reading bar length from geosvc. + if(bar.getSlayer()==1) return 600.; + else return 480.-bar.getDlayer()*10.; +} + +/* +dd4hep::Position CRDEcalDigiAlg::GetCellPos(dd4hep::Position& pos, CaloBar& bar){ + dd4hep::Position rpos = pos-bar.getPosition(); + TVector3 vec(0,0,0); + if(bar.getSlayer()==1) vec.SetXYZ(0, 0, floor(rpos.z()/10)*10+5 ); + else if(bar.getSlayer()==0){ + if((bar.getModule()==0||bar.getModule()==4) && bar.getDlayer()%2==1) vec.SetXYZ(floor(rpos.x()/10)*10+5,0,0); + if((bar.getModule()==0||bar.getModule()==4) && bar.getDlayer()%2==0) vec.SetXYZ(floor((rpos.x()-5)/10)*10+10,0,0); + if((bar.getModule()==2||bar.getModule()==6) && bar.getDlayer()%2==1) vec.SetXYZ(0, floor(rpos.y()/10)*10+5,0); + if((bar.getModule()==2||bar.getModule()==6) && bar.getDlayer()%2==0) vec.SetXYZ(0, floor((rpos.y()-5)/10)*10+10,0); + if(bar.getModule()==1 || bar.getModule()==5){ + TVector3 unitv(1./sqrt(2), -1./sqrt(2), 0); + if(bar.getDlayer()%2==1) vec = (floor(rpos.Dot(unitv)/10)*10+5)*unitv; + if(bar.getDlayer()%2==0) vec = (floor((rpos.Dot(unitv)-5)/10)*10+10)*unitv; + } + if(bar.getModule()==3 || bar.getModule()==7){ + TVector3 unitv(1./sqrt(2), 1./sqrt(2), 0); + if(bar.getDlayer()%2==1) vec = (floor(rpos.Dot(unitv)/10)*10+5)*unitv; + if(bar.getDlayer()%2==0) vec = (floor((rpos.Dot(unitv)-5)/10)*10+10)*unitv; + } + } + dd4hep::Position relv(vec.x(), vec.y(), vec.z()); + return relv+bar.getPosition(); +} + + +edm4hep::MutableSimCalorimeterHit CRDEcalDigiAlg::find(edm4hep::SimCalorimeterHitCollection& m_col, dd4hep::Position& pos){ + for(int i=0;i& m_col, unsigned long long& cellid) const{ + for(int i=0;i +#include +#include +#include "DetInterface/IGeomSvc.h" +#include "TVector3.h" +#include "TRandom3.h" +#include "TFile.h" +#include "TString.h" +#include "TH3.h" +#include "TH1.h" + +#include +#include "time.h" +#include +#include + +#define C 299.79 // unit: mm/ns +#define PI 3.141592653 + +class CRDEcalDigiAlg : public GaudiAlgorithm +{ + +public: + + CRDEcalDigiAlg(const std::string& name, ISvcLocator* svcLoc); + + /** Called at the begin of the job before anything is read. + * Use to initialize the processor, e.g. book histograms. + */ + virtual StatusCode initialize() ; + + /** Called for every event - the working horse. + */ + virtual StatusCode execute() ; + + /** Called after data processing for clean up. + */ + virtual StatusCode finalize() ; + + + StatusCode MergeHits(const edm4hep::SimCalorimeterHitCollection& m_col, std::vector& m_hits); + + double GetBarLength(CaloBar& bar); //TODO: should read from geom file! + edm4hep::MutableSimCalorimeterHit find(const std::vector& m_col, unsigned long long& cellid) const; + + void Clear(); + +protected: + + SmartIF m_geosvc; + //SmartIF m_edmsvc; + typedef std::vector FloatVec; + typedef std::map MCParticleToEnergyWeightMap; + + int _nEvt ; + float m_length; + TRandom3 rndm; + TFile* m_wfile; + TTree* t_SimCont; + TTree* t_SimBar; + + double totE; + FloatVec m_step_t; // yyy: time of each step + FloatVec m_step_x, m_step_y, m_step_z, m_step_E, m_step_T1, m_step_T2, m_stepBar_x, m_stepBar_y, m_stepBar_z; + FloatVec m_simBar_x, m_simBar_y, m_simBar_z, m_simBar_T1, m_simBar_T2, m_simBar_Q1, m_simBar_Q2, m_simBar_dlayer, m_simBar_part, m_simBar_stave, m_simBar_slayer, m_simBar_module; + std::vector m_simBar_cellID; + + + dd4hep::rec::CellIDPositionConverter* m_cellIDConverter; + dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + + Gaudi::Property m_scale{ this, "Scale", 1 }; + + // Input collections + DataHandle r_SimCaloCol{"SimCaloCol", Gaudi::DataHandle::Reader, this}; + mutable Gaudi::Property _readoutName{this, "ReadOutName", "CaloHitsCollection", "Readout name"}; + mutable Gaudi::Property _filename{this, "OutFileName", "testout.root", "Output file name"}; + + //Input parameters + mutable Gaudi::Property _writeNtuple{this, "WriteNtuple", 1, "Write ntuple"}; + mutable Gaudi::Property _Nskip{this, "SkipEvt", 0, "Skip event"}; + mutable Gaudi::Property _seed{this, "Seed", 2131, "Random Seed"}; + mutable Gaudi::Property _Debug{this, "Debug", 0, "Debug level"}; + mutable Gaudi::Property _Eth {this, "EnergyThreshold", 0.001, "Energy Threshold (/GeV)"}; + mutable Gaudi::Property r_cali{this, "CalibrECAL", 1, "Calibration coefficients for ECAL"}; + mutable Gaudi::Property Latt{this, "AttenuationLength", 7000, "Crystal Attenuation Length(mm)"}; + mutable Gaudi::Property Tres{this, "TimeResolution", 0.1, "Crystal time resolution in one side (ns)"}; + mutable Gaudi::Property nMat{this, "MatRefractive", 2.15, "Material refractive index of crystal"}; + mutable Gaudi::Property Tinit{this, "InitalTime", 2, "Start time (ns)"}; + + mutable Gaudi::Property _Qthfrac {this, "ChargeThresholdFrac", 0.05, "Charge threshold fraction"}; + + + // Output collections + DataHandle w_DigiCaloCol{"DigiCaloCol", Gaudi::DataHandle::Writer, this}; + DataHandle w_CaloAssociationCol{"ECALBarrelAssoCol", Gaudi::DataHandle::Writer, this}; + DataHandle w_MCPCaloAssociationCol{"ECALBarrelParticleAssoCol", Gaudi::DataHandle::Writer, this}; +}; + +#endif diff --git a/Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.cpp b/Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.cpp new file mode 100644 index 000000000..7491716d3 --- /dev/null +++ b/Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.cpp @@ -0,0 +1,214 @@ +// /* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */ +// // Unit in code: mm, ns. +// // NOTE: This digitialization highly matches detector geometry CRDEcalBarrel_v01. +// // TODO: read geometry info automatically. + +#include "CRDHcalDigiAlg.h" + +#include "edm4hep/SimCalorimeterHit.h" +#include "edm4hep/CalorimeterHit.h" +#include "edm4hep/Vector3f.h" +#include "edm4hep/Cluster.h" + +#include "DD4hep/Detector.h" +#include +#include + +#include "TVector3.h" +#include +#include +#include +#include +#include + +// #include +// #include + +#define C 299.79 // unit: mm/ns +#define PI 3.141592653 +using namespace std; +using namespace dd4hep; + +DECLARE_COMPONENT( CRDHcalDigiAlg ) + +CRDHcalDigiAlg::CRDHcalDigiAlg(const std::string& name, ISvcLocator* svcLoc) + : GaudiAlgorithm(name, svcLoc), + _nEvt(0) +{ + + // Input collections + declareProperty("SimCaloHitCollection", r_SimCaloCol, "Handle of the Input SimCaloHit collection"); + + // Output collections + declareProperty("CaloHitCollection", w_DigiCaloCol, "Handle of Digi CaloHit collection"); + declareProperty("CaloAssociationCollection", w_CaloAssociationCol, "Handle of CaloAssociation collection"); + declareProperty("CaloMCPAssociationCollection", w_MCPCaloAssociationCol, "Handle of CaloAssociation collection"); +} + +StatusCode CRDHcalDigiAlg::initialize() +{ + if(_writeNtuple){ + std::string s_outfile = _filename; + m_wfile = new TFile(s_outfile.c_str(), "recreate"); + t_simHit = new TTree("simHit", "simHit"); + + t_simHit->Branch("simHit_x", &m_simHit_x); + t_simHit->Branch("simHit_y", &m_simHit_y); + t_simHit->Branch("simHit_z", &m_simHit_z); + t_simHit->Branch("simHit_E", &m_simHit_E); + t_simHit->Branch("simHit_steps", &m_simHit_steps); + + t_simHit->Branch("simHit_module", &m_simHit_module); + t_simHit->Branch("simHit_stave", &m_simHit_stave); + t_simHit->Branch("simHit_layer", &m_simHit_layer); + t_simHit->Branch("simHit_tower", &m_simHit_tower); + t_simHit->Branch("simHit_slice", &m_simHit_slice); + t_simHit->Branch("simHit_cellID", &m_simHit_cellID); + } + std::cout<<"CRDHcalDigiAlg::m_scale="<("GeomSvc"); + if ( !m_geosvc ) throw "CRDHcalDigiAlg :Failed to find GeomSvc ..."; + dd4hep::Detector* m_dd4hep = m_geosvc->lcdd(); + if ( !m_dd4hep ) throw "CRDHcalDigiAlg :Failed to get dd4hep::Detector ..."; + m_cellIDConverter = new dd4hep::rec::CellIDPositionConverter(*m_dd4hep); + m_decoder = m_geosvc->getDecoder(_readoutName); + if (!m_decoder) { + error() << "Failed to get the decoder. " << endmsg; + return StatusCode::FAILURE; + } + + //m_edmsvc = service("CRDEcalSvc"); + //if ( !m_edmsvc ) throw "CRDHcalDigiAlg :Failed to find CRDEcalSvc ..."; + + rndm.SetSeed(_seed); + std::cout<<"CRDHcalDigiAlg::initialize"<=1) std::cout<<"digi, input sim hit size="<< SimHitCol->size() <size(); isim++){ + + auto simhit = SimHitCol->at(isim); + if(!simhit.isAvailable()) continue; + if(simhit.getEnergy()==0) continue; + + unsigned long long id = simhit.getCellID(); + double Ehit = simhit.getEnergy(); + //Energy threshold + if(Ehit<_MIPCali*_Eth_Mip) continue; + + //Global calibration. + //TODO: add more digitization terms here. + double Ehit_cali = Ehit*r_cali; + + //Loop contributions to get hit time and MCParticle. + double Thit_ave = 0.; + double Ehit_raw = 0.; + MCParticleToEnergyWeightMap MCPEnMap; MCPEnMap.clear(); + for(int iConb=0; iConbcreate(); + digiHit.setCellID(id); + digiHit.setEnergy(Ehit_cali); + digiHit.setTime(Thit_ave); + digiHit.setPosition(simhit.getPosition()); + + //Create SimHit-DigiHit association. + auto rel = caloAssoVec->create(); + rel.setRec(digiHit); + rel.setSim(simhit); + rel.setWeight(1.); + + //Create DigiHit-MCParticle association. + for(auto iter : MCPEnMap){ + auto rel_MC = caloMCPAssoVec->create(); + rel_MC.setRec(digiHit); + rel_MC.setSim(iter.first); + rel_MC.setWeight(iter.second/Ehit_raw); + } + + if(_writeNtuple){ + m_simHit_x.push_back(digiHit.getPosition().x); + m_simHit_y.push_back(digiHit.getPosition().y); + m_simHit_z.push_back(digiHit.getPosition().z); + m_simHit_E.push_back(digiHit.getEnergy()); + m_simHit_steps.push_back(simhit.contributions_size()); + m_simHit_module.push_back(m_decoder->get(id, "module")); + m_simHit_stave.push_back(m_decoder->get(id, "stave")); + m_simHit_layer.push_back(m_decoder->get(id, "layer")); + m_simHit_slice.push_back(m_decoder->get(id, "slice")); + m_simHit_tower.push_back(m_decoder->get(id, "tower")); + m_simHit_cellID.push_back(id); + } + } + + if(_writeNtuple) t_simHit->Fill(); + + _nEvt ++ ; + return StatusCode::SUCCESS; +} + +StatusCode CRDHcalDigiAlg::finalize() +{ + if(_writeNtuple){ + m_wfile->cd(); + t_simHit->Write(); + m_wfile->Close(); + delete m_wfile, t_simHit; + } + + info() << "Processed " << _nEvt << " events " << endmsg; + delete m_cellIDConverter, m_decoder, m_geosvc; + return GaudiAlgorithm::finalize(); +} + + +void CRDHcalDigiAlg::Clear(){ + m_simHit_x.clear(); + m_simHit_y.clear(); + m_simHit_z.clear(); + m_simHit_E.clear(); + m_simHit_steps.clear(); + m_simHit_module.clear(); + m_simHit_stave.clear(); + m_simHit_layer.clear(); + m_simHit_slice.clear(); + m_simHit_tower.clear(); + m_simHit_cellID.clear(); +} + diff --git a/Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.h b/Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.h new file mode 100644 index 000000000..6e3571e89 --- /dev/null +++ b/Digitisers/CRDEcalDigi/src/CRDHcalDigiAlg.h @@ -0,0 +1,100 @@ +#ifndef CRD_HCAL_DIGI_ALG_H +#define CRD_HCAL_DIGI_ALG_H + +#include "k4FWCore/DataHandle.h" +#include "GaudiAlg/GaudiAlgorithm.h" +#include "edm4hep/MutableCaloHitContribution.h" +#include "edm4hep/MutableSimCalorimeterHit.h" +#include "edm4hep/SimCalorimeterHit.h" +#include "edm4hep/CalorimeterHit.h" +#include "edm4hep/CalorimeterHitCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/MCRecoCaloAssociationCollection.h" +#include "edm4hep/MCRecoCaloParticleAssociationCollection.h" +#include "edm4hep/MCParticle.h" + +#include +#include +#include +#include "DetInterface/IGeomSvc.h" +#include "TVector3.h" +#include "TRandom3.h" +#include "TFile.h" +#include "TString.h" +#include "TH3.h" +#include "TH1.h" + +#include +#include "time.h" +#include +#include + +#define C 299.79 // unit: mm/ns +#define PI 3.141592653 + +class CRDHcalDigiAlg : public GaudiAlgorithm +{ + +public: + + CRDHcalDigiAlg(const std::string& name, ISvcLocator* svcLoc); + + /** Called at the begin of the job before anything is read. + * Use to initialize the processor, e.g. book histograms. + */ + virtual StatusCode initialize() ; + + /** Called for every event - the working horse. + */ + virtual StatusCode execute() ; + + /** Called after data processing for clean up. + */ + virtual StatusCode finalize() ; + + void Clear(); + +protected: + + SmartIF m_geosvc; + //SmartIF m_edmsvc; + typedef std::vector FloatVec; + typedef std::map MCParticleToEnergyWeightMap; + + int _nEvt ; + TRandom3 rndm; + TFile* m_wfile; + TTree* t_simHit; + + FloatVec m_simHit_x, m_simHit_y, m_simHit_z, m_simHit_E, m_simHit_slice, m_simHit_layer, m_simHit_tower, m_simHit_stave, m_simHit_module; + std::vector m_simHit_cellID; + std::vector m_simHit_steps; + + + dd4hep::rec::CellIDPositionConverter* m_cellIDConverter; + dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + + Gaudi::Property m_scale{ this, "Scale", 1 }; + + // Input collections + DataHandle r_SimCaloCol{"SimCaloCol", Gaudi::DataHandle::Reader, this}; + mutable Gaudi::Property _readoutName{this, "ReadOutName", "CaloHitsCollection", "Readout name"}; + mutable Gaudi::Property _filename{this, "OutFileName", "testout.root", "Output file name"}; + + //Input parameters + mutable Gaudi::Property _writeNtuple{this, "WriteNtuple", 1, "Write ntuple"}; + mutable Gaudi::Property _Nskip{this, "SkipEvt", 0, "Skip event"}; + mutable Gaudi::Property _seed{this, "Seed", 2131, "Random Seed"}; + mutable Gaudi::Property _Debug{this, "Debug", 0, "Debug level"}; + mutable Gaudi::Property _MIPCali {this, "MIPResponse", 0.0005, "MIP response (/GeV)"}; + mutable Gaudi::Property _Eth_Mip {this, "MIPThreshold", 0.5, "Energy Threshold (/MIP)"}; + mutable Gaudi::Property r_cali{this, "CalibrHCAL", 1, "Global calibration coefficients for HCAL"}; + + + // Output collections + DataHandle w_DigiCaloCol{"DigiCaloCol", Gaudi::DataHandle::Writer, this}; + DataHandle w_CaloAssociationCol{"HCALBarrelAssoCol", Gaudi::DataHandle::Writer, this}; + DataHandle w_MCPCaloAssociationCol{"HCALBarrelParticleAssoCol", Gaudi::DataHandle::Writer, this}; +}; + +#endif diff --git a/Digitisers/CRDEcalDigi/src/CaloBar.h b/Digitisers/CRDEcalDigi/src/CaloBar.h new file mode 100644 index 000000000..631925dae --- /dev/null +++ b/Digitisers/CRDEcalDigi/src/CaloBar.h @@ -0,0 +1,56 @@ +#ifndef _CRD_CALOBAR_ +#define _CRD_CALOBAR_ + +#include +#include "TVector3.h" + +class CaloBar{ + +public: + CaloBar(unsigned long long _cellID, int _system, int _module, int _slayer, int _dlayer, int _part, int _stave, int _bar, TVector3 _pos, double _Q1, double _Q2, double _T1, double _T2) + : cellID(_cellID), system(_system), module(_module), stave(_stave), dlayer(_dlayer), part(_part), slayer(_slayer), bar(_bar), position(_pos), Q1(_Q1), Q2(_Q2), T1(_T1), T2(_T2) {}; + CaloBar() {}; + + inline bool operator == (const CaloBar &x) const{ + return ( (cellID == x.cellID) && getEnergy()==x.getEnergy() ); + } + unsigned long long getcellID() const { return cellID; } + int getSystem() const { return system; } + int getModule() const { return module; } + int getStave() const { return stave; } + int getDlayer() const { return dlayer; } + int getPart() const { return part; } + int getSlayer() const { return slayer; } + int getBar() const { return bar; } + double getQ1() const { return Q1; } + double getQ2() const { return Q2; } + double getT1() const { return T1; } + double getT2() const { return T2; } + + TVector3 getPosition() const { return position; } + double getEnergy() const { return (Q1+Q2)/2.; } + + void setcellID(unsigned long long _cellid) { cellID = _cellid; } + void setcellID(int _system, int _module, int _stave, int _dlayer, int _part, int _slayer, int _bar) { system=_system; module=_module; stave=_stave; dlayer=_dlayer; part=_part; slayer=_slayer; bar=_bar; } + void setPosition( TVector3 posv3) { position.SetXYZ( posv3.x(), posv3.y(), posv3.z() ); } + void setQ(double _q1, double _q2) { Q1=_q1; Q2=_q2; } + void setT(double _t1, double _t2) { T1=_t1; T2=_t2; } + +private: + unsigned long long cellID; + int system; + int module; + int stave; + int dlayer; + int part; + int slayer; + int bar; + TVector3 position; + double Q1; // Q in left readout + double Q2; // Q in right readout; + double T1; // T in left readout; + double T2; // T in right readout; + +}; + +#endif diff --git a/Digitisers/CRDEcalDigi/src/CaloHit.h b/Digitisers/CRDEcalDigi/src/CaloHit.h new file mode 100644 index 000000000..c503f987a --- /dev/null +++ b/Digitisers/CRDEcalDigi/src/CaloHit.h @@ -0,0 +1,50 @@ +#ifndef _CRD_CALOHIT_ +#define _CRD_CALOHIT_ + +#include +#include "TVector3.h" + +class CaloHit{ + +public: + CaloHit(unsigned long long _cellID, int _system, int _module, int _layer, int _stave, int _tower, int _slice, TVector3 _pos, double _E, int _step_numbers) + : cellID(_cellID), system(_system), module(_module), layer(_layer), stave(_stave), tower(_tower), slice(_slice), position(_pos), E(_E), step_numbers(_step_numbers) {}; + CaloHit() {}; + + inline bool operator == (const CaloHit &x) const{ + return ( (cellID == x.cellID) && getEnergy()==x.getEnergy() ); + } + unsigned long long getcellID() const { return cellID; } + int getSystem() const { return system; } + int getModule() const { return module; } + int getStave() const { return stave; } + int getLayer() const { return layer; } + int getTower() const { return tower; } + int getSlice() const { return slice; } + + TVector3 getPosition() const { return position; } + double getEnergy() const { return E; } + int getStep_numbers() const { return step_numbers; } + + void setcellID(unsigned long long _cellid) { cellID = _cellid; } + void setcellID(int _system, int _module, int _layer, int _stave, int _tower, int _slice) { system=_system; module=_module; stave=_stave; layer=_layer; tower=_tower; slice=_slice;} + void setPosition( TVector3 posv3) { position.SetXYZ( posv3.x(), posv3.y(), posv3.z() ); } + void setE(double _E) { E=_E; } + void setStep_numbers(int _step_numbers) { step_numbers=_step_numbers; } + + +private: + unsigned long long cellID; + int system; + int module; + int stave; + int layer; + int tower; + int slice; + TVector3 position; + double E; + int step_numbers; + +}; + +#endif diff --git a/Digitisers/CRDEcalDigi/src/HitStep.h b/Digitisers/CRDEcalDigi/src/HitStep.h new file mode 100644 index 000000000..2ae0715be --- /dev/null +++ b/Digitisers/CRDEcalDigi/src/HitStep.h @@ -0,0 +1,26 @@ +#ifndef HIT_STEP_H +#define HIT_STEP_H + +class HitStep{ + +public: + HitStep (double _Q, double _T): Q(_Q), T(_T) {}; + //HitStep (double _Q, double _T) { Q=_Q; T=_T; }; + HitStep() {}; + + void setQ(double _Q) { Q =_Q; } + void setT(double _T) { T =_T; } + + double getQ() const { return Q; } + double getT() const { return T; } + inline bool operator < (const HitStep &x) const { + return T one space point +from Configurables import SpacePointBuilderAlg +spSET = SpacePointBuilderAlg("SETBuilder") +spSET.TrackerHitCollection = sethitname +spSET.TrackerHitAssociationCollection = "SETTrackerHitAssociation" +spSET.SpacePointCollection = setspname +spSET.SpacePointAssociationCollection = "SETSpacePointAssociation" +#spSET.OutputLevel = DEBUG + + +digiFTD = PlanarDigiAlg("FTDDigi") +digiFTD.IsStrip = False +digiFTD.SimTrackHitCollection = "FTDCollection" +digiFTD.TrackerHitCollection = ftdhitname +digiFTD.TrackerHitAssociationCollection = "FTDTrackerHitAssociation" +digiFTD.ResolutionU = [0.003, 0.003, 0.0072, 0.0072, 0.0072, 0.0072, 0.0072] +digiFTD.ResolutionV = [0.003, 0.003, 0.0072, 0.0072, 0.0072, 0.0072, 0.0072] +digiFTD.UsePlanarTag = True +#digiFTD.OutputLevel = DEBUG + +# two strip tracker hits -> one space point +from Configurables import SpacePointBuilderAlg +spFTD = SpacePointBuilderAlg("FTDBuilder") +spFTD.TrackerHitCollection = ftdhitname +spFTD.TrackerHitAssociationCollection = "FTDTrackerHitAssociation" +spFTD.SpacePointCollection = ftdspname +spFTD.SpacePointAssociationCollection = "FTDSpacePointAssociation" +#spFTD.OutputLevel = DEBUG + +from Configurables import TPCDigiAlg +digiTPC = TPCDigiAlg("TPCDigi") +digiTPC.TPCCollection = "TPCCollection" +digiTPC.TPCLowPtCollection = "TPCLowPtCollection" +digiTPC.TPCTrackerHitsCol = tpchitname +digiTPC.TPCTrackerHitAssCol = "TPCTrackerHitAssociation" +#digiTPC.OutputLevel = DEBUG + +# tracking +from Configurables import SiliconTrackingAlg +tracking = SiliconTrackingAlg("SiliconTracking") +tracking.HeaderCol = "EventHeader" +tracking.VTXHitCollection = vxdhitname +tracking.SITHitCollection = sithitname +tracking.FTDPixelHitCollection = ftdhitname +tracking.FTDSpacePointCollection = ftdspname +tracking.SITRawHitCollection = sithitname +tracking.FTDRawHitCollection = ftdhitname +tracking.UseSIT = True +tracking.SmoothOn = False +tracking.DumpTime = False +tracking.NDivisionsInTheta = 10 +#tracking.OutputLevel = DEBUG + +from Configurables import ForwardTrackingAlg +forward = ForwardTrackingAlg("ForwardTracking") +forward.FTDPixelHitCollection = ftdhitname +forward.FTDSpacePointCollection = ftdspname +forward.FTDRawHitCollection = ftdhitname +forward.Chi2ProbCut = 0.0 +forward.HitsPerTrackMin = 3 +forward.BestSubsetFinder = "SubsetSimple" +forward.Criteria = ["Crit2_DeltaPhi","Crit2_StraightTrackRatio","Crit3_3DAngle","Crit3_ChangeRZRatio","Crit3_IPCircleDist","Crit4_3DAngleChange","Crit4_DistToExtrapolation", + "Crit2_DeltaRho","Crit2_RZRatio","Crit3_PT"] +forward.CriteriaMin = [0, 0.9, 0, 0.995, 0, 0.8, 0, 20, 1.002, 0.1, 0, 0.99, 0, 0.999, 0, 0.99, 0] +forward.CriteriaMax = [30, 1.02, 10, 1.015, 20, 1.3, 1.0, 150, 1.08, 99999999, 0.8, 1.01, 0.35, 1.001, 1.5, 1.01, 0.05] +forward.DumpTime = False +#forward.OutputLevel = DEBUG + +from Configurables import TrackSubsetAlg +subset = TrackSubsetAlg("TrackSubset") +subset.TrackInputCollections = ["ForwardTracks", "SiTracks"] +subset.RawTrackerHitCollections = [vxdhitname, sithitname, ftdhitname, ftdspname] +subset.TrackSubsetCollection = "SubsetTracks" +subset.DumpTime = False +#subset.OutputLevel = DEBUG + +#TODO: DC reconstruction, as preliminary, use Clupatra like as TPC +from Configurables import ClupatraAlg +clupatra = ClupatraAlg("Clupatra") +clupatra.TPCHitCollection = tpchitname +#clupatra.DistanceCut = 100. +#clupatra.MaxDeltaChi2 = 100. +#clupatra.Chi2Cut = 150. +#clupatra.OutputLevel = DEBUG + +from Configurables import FullLDCTrackingAlg +full = FullLDCTrackingAlg("FullTracking") +full.VTXTrackerHits = vxdhitname +full.SITTrackerHits = sithitname +full.TPCTrackerHits = tpchitname # add TPC or DC tracker hit here, if TPC or DC track is set by full.TPCTracks +full.SETTrackerHits = sethitname +full.FTDPixelTrackerHits = ftdhitname +full.FTDSpacePoints = ftdspname +full.SITRawHits = sithitname +full.SETRawHits = sethitname +full.FTDRawHits = ftdhitname +full.VTXHitRelCol = "VXDTrackerHitAssociation" +full.SITHitRelCol = "SITTrackerHitAssociation" +full.SETHitRelCol = "SETSpacePointAssociation" +full.FTDHitRelCol = "FTDTrackerHitAssociation" +full.TPCHitRelCol = "TPCTrackerHitAssociation" +full.TPCTracks = "ClupatraTracks" # add standalone TPC or DC track here +full.SiTracks = "SubsetTracks" +full.OutputTracks = "MarlinTrkTracks" +full.DumpTime = False +full.SITHitToTrackDistance = 3. +full.SETHitToTrackDistance = 5. +full.MinChi2ProbForSiliconTracks = 0 +#full.OutputLevel = DEBUG + +dedxoption = "BetheBlochEquationDedxSimTool" +from Configurables import DriftChamberSensDetTool +dc_sensdettool = DriftChamberSensDetTool("DriftChamberSensDetTool") +dc_sensdettool.DedxSimTool = dedxoption + +from Configurables import DummyDedxSimTool +from Configurables import BetheBlochEquationDedxSimTool + +if dedxoption == "DummyDedxSimTool": + dedx_simtool = DummyDedxSimTool("DummyDedxSimTool") +elif dedxoption == "BetheBlochEquationDedxSimTool": + dedx_simtool = BetheBlochEquationDedxSimTool("BetheBlochEquationDedxSimTool") + dedx_simtool.material_Z = 2 + dedx_simtool.material_A = 4 + dedx_simtool.scale = 10 + dedx_simtool.resolution = 0.0001 + +# output +from Configurables import PodioOutput +out = PodioOutput("outputalg") +out.filename = "CRD_HCal_GamGam_1deg.root" +out.outputCommands = ["keep *"] + +# ApplicationMgr +from Configurables import ApplicationMgr +ApplicationMgr( + TopAlg = [genalg, detsimalg, digiVXD, digiSIT, digiSET, digiFTD, spSET, digiTPC, tracking, forward, subset, full, out], + #TopAlg = [genalg, detsimalg, digiVXD, digiSIT, digiSET, digiFTD, spSET, digiTPC, tracking, forward, subset, out], + EvtSel = 'NONE', + EvtMax = 100, + ExtSvc = [rndmengine, rndmgensvc, dsvc, evtseeder, geosvc, gearsvc, tracksystemsvc], + #ExtSvc = [rndmengine, rndmgensvc, dsvc, geosvc], + OutputLevel=INFO +) diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp index ee5f1feb2..18de4bf53 100755 --- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp +++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp @@ -8,14 +8,6 @@ #include #include #include -#if __has_include("edm4hep/EDM4hepVersion.h") -#include "edm4hep/EDM4hepVersion.h" -#else -// Copy the necessary parts from the header above to make whatever we need to work here -#define EDM4HEP_VERSION(major, minor, patch) ((UINT64_C(major) << 32) | (UINT64_C(minor) << 16) | (UINT64_C(patch))) -// v00-09 is the last version without the capitalization change of the track vector members -#define EDM4HEP_BUILD_VERSION EDM4HEP_VERSION(0, 9, 0) -#endif #include #include @@ -124,7 +116,12 @@ FullLDCTrackingAlg::FullLDCTrackingAlg(const std::string& name, ISvcLocator* svc declareProperty("SiTracks", _SiTrackColHdl, "Si Track Collection"); // Input relation collections - + declareProperty("VTXHitRelCol", _VTXHitRelColHdl, "Vertex TrackerHit association name"); + declareProperty("SITHitRelCol", _SITHitRelColHdl, "SIT TrackerHit association name"); + declareProperty("SETHitRelCol", _SETHitRelColHdl, "SET TrackerHit association name"); + declareProperty("FTDHitRelCol", _FTDHitRelColHdl, "FTD TrackerHit association name"); + declareProperty("TPCHitRelCol", _TPCHitRelColHdl, "TPC TrackerHit association name"); + /* registerInputCollection(LCIO::LCRELATION, "TPCTracksMCPRelColl", @@ -140,6 +137,7 @@ FullLDCTrackingAlg::FullLDCTrackingAlg(const std::string& name, ISvcLocator* svc */ // Output track collection declareProperty("OutputTracks", _OutputTrackColHdl, "Full LDC track collection name"); + declareProperty("OutputTrackAssociationCollection", _OutputTrkRelColHdl, "Track - MCParticle association collection name"); } @@ -212,7 +210,8 @@ StatusCode FullLDCTrackingAlg::execute() { auto stopwatch = TStopwatch(); // debug() << endmsg; auto outCol = _OutputTrackColHdl.createAndPut(); - + auto outRelCol = _OutputTrkRelColHdl.createAndPut(); + prepareVectors(); debug() << "************************************PrepareVectors done..." << endmsg; @@ -247,11 +246,20 @@ StatusCode FullLDCTrackingAlg::execute() { catch(std::runtime_error& e){ error() << e.what() << std::endl; } + + //Create the association + try{ + CreateTrackAssociation(outCol, outRelCol); + } + catch(std::runtime_error& e){ + error() << e.what() << std::endl; + } + CleanUp(); debug() << "Cleanup is done." << endmsg; _nEvt++; // getchar(); - if(m_tuple){ + if(m_dumpTime && m_tuple){ m_timeTotal = stopwatch.RealTime()*1000; m_tuple->write(); } @@ -521,19 +529,11 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr float z0TrkCand = trkCand->getZ0(); // float phi0TrkCand = trkCand->getPhi(); // FIXME, fucd -#if EDM4HEP_BUILD_VERSION > EDM4HEP_VERSION(0, 9, 0) - int nhits_in_vxd = track.getSubdetectorHitNumbers(0); - int nhits_in_ftd = track.getSubdetectorHitNumbers(1); - int nhits_in_sit = track.getSubdetectorHitNumbers(2); - int nhits_in_tpc = track.getSubdetectorHitNumbers(3); - int nhits_in_set = track.getSubdetectorHitNumbers(4); -#else int nhits_in_vxd = track.getSubDetectorHitNumbers(0); int nhits_in_ftd = track.getSubDetectorHitNumbers(1); int nhits_in_sit = track.getSubDetectorHitNumbers(2); int nhits_in_tpc = track.getSubDetectorHitNumbers(3); int nhits_in_set = track.getSubDetectorHitNumbers(4); -#endif //int nhits_in_vxd = Track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::VXD - 2 ]; //int nhits_in_ftd = Track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::FTD - 2 ]; //int nhits_in_sit = Track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::SIT - 2 ]; @@ -595,7 +595,7 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr debug() << " Add Track to final Collection: ID = " << track.id() << " for trkCand "<< trkCand << endmsg; } } - if(m_tuple) m_timeKalman = stopwatch.RealTime()*1000; + if(m_dumpTime && m_tuple) m_timeKalman = stopwatch.RealTime()*1000; // streamlog_out(DEBUG5) << endmsg; debug() << "Number of accepted " << _OutputTrackColHdl.fullKey() << " = " << nTotTracks << endmsg; debug() << "Total 4-momentum of " << _OutputTrackColHdl.fullKey() << " : E = " << eTot @@ -1166,18 +1166,13 @@ void FullLDCTrackingAlg::prepareVectors() { trackExt->setNDF(tpcTrack.getNdf()); trackExt->setChi2(tpcTrack.getChi2()); for (int iHit=0;iHitsecond; + if(it==mapTrackerHits.end()) error() << "Cannot find hit " << hit.id() << " in map" << endmsg; + else continue; + TrackerHitExtended * hitExt = it->second; //info() << hit.id() << " " << hitExt << endmsg; hitExt->setTrackExtended( trackExt ); trackExt->addTrackerHitExtended( hitExt ); @@ -1236,17 +1231,8 @@ void FullLDCTrackingAlg::prepareVectors() { char strg[200]; HelixClass helixSi; for (int iHit=0;iHitsecond; + edm4hep::TrackerHit hit = siTrack.getTrackerHits(iHit);//hitVec[iHit]; + TrackerHitExtended * hitExt = mapTrackerHits[hit]; hitExt->setTrackExtended( trackExt ); trackExt->addTrackerHitExtended( hitExt ); @@ -1548,8 +1534,8 @@ TrackExtended * FullLDCTrackingAlg::CombineTracks(TrackExtended * tpcTrack, Trac int nTPCHits = int(tpcHitVec.size()); int nHits = nTPCHits + nSiHits; - //debug() << "FullLDCTrackingAlg::CombineTracks nSiHits = " << nSiHits << endmsg; - //debug() << "FullLDCTrackingAlg::CombineTracks nTPCHits = " << nTPCHits << endmsg; + //std::cout << "FullLDCTrackingAlg::CombineTracks nSiHits = " << nSiHits << endmsg; + //std::cout << "FullLDCTrackingAlg::CombineTracks nTPCHits = " << nTPCHits << endmsg; TrackerHitVec trkHits; trkHits.reserve(nHits); @@ -1762,7 +1748,7 @@ TrackExtended * FullLDCTrackingAlg::CombineTracks(TrackExtended * tpcTrack, Trac tpcHitInFit.push_back(tpcHitVec[i]); } } - + debug() << "FullLDCTrackingAlg::CombineTracks: Check for Silicon Hit rejections ... " << endmsg; if ( (int)siOutliers.size() > _maxAllowedSiHitRejectionsForTrackCombination ) { @@ -3819,7 +3805,7 @@ HelixClass * FullLDCTrackingAlg::GetExtrapolationHelix( TrackExtended * track) { ts_at_calo = getTrackStateAt(trk_lcio, 4/*lcio::TrackState::AtCalorimeter*/); debug() << "FullLDCTrackingAlg::GetExtrapolationHelix set ts_at_calo with ref_z = " << z_ref << endmsg; - } + } } } } @@ -5013,3 +4999,88 @@ void FullLDCTrackingAlg::checkTrackState(int type){ } } } + +void FullLDCTrackingAlg::CreateTrackAssociation(edm4hep::TrackCollection* trkCol, edm4hep::MCRecoTrackParticleAssociationCollection* relCol){ + + const edm4hep::MCRecoTrackerAssociationCollection* VTXHitRelCol = _VTXHitRelColHdl.get(); + const edm4hep::MCRecoTrackerAssociationCollection* SITHitRelCol = _SITHitRelColHdl.get(); + const edm4hep::MCRecoTrackerAssociationCollection* SETHitRelCol = _SETHitRelColHdl.get(); + const edm4hep::MCRecoTrackerAssociationCollection* FTDHitRelCol = _FTDHitRelColHdl.get(); + const edm4hep::MCRecoTrackerAssociationCollection* TPCHitRelCol = _TPCHitRelColHdl.get(); + + + //Loop for the tracks + for(int itrk=0; itrksize(); itrk++){ + edm4hep::Track iTrk = trkCol->at(itrk); + if(!iTrk.isAvailable()) continue; + + std::map MCParticleToNhitMap; + //Loop for trackerhits in the track + int nHit = iTrk.trackerHits_size(); + for(int iHit=0; iHitsize() && !findSimHit; ++irel){ + if(hit == VTXHitRelCol->at(irel).getRec()){ + p_simHit = VTXHitRelCol->at(irel).getSim(); + findSimHit = true; + break; + } + } + + for(int irel = 0; irelsize() && !findSimHit; ++irel){ + if(hit == SITHitRelCol->at(irel).getRec()){ + p_simHit = SITHitRelCol->at(irel).getSim(); + findSimHit = true; + break; + } + } + + for(int irel = 0; irelsize() && !findSimHit; ++irel){ + if(hit == SETHitRelCol->at(irel).getRec()){ + p_simHit = SETHitRelCol->at(irel).getSim(); + findSimHit = true; + break; + } + } + + for(int irel = 0; irelsize() && !findSimHit; ++irel){ + if(hit == FTDHitRelCol->at(irel).getRec()){ + p_simHit = FTDHitRelCol->at(irel).getSim(); + findSimHit = true; + break; + } + } + + for(int irel = 0; irelsize() && !findSimHit; ++irel){ + if(hit == TPCHitRelCol->at(irel).getRec()){ + p_simHit = TPCHitRelCol->at(irel).getSim(); + findSimHit = true; + break; + } + } + + //Get MCP from SimTrackerHit + if(findSimHit && p_simHit.isAvailable()){ + const edm4hep::MCParticle mcp = p_simHit.getMCParticle(); + if(MCParticleToNhitMap.find(mcp)==MCParticleToNhitMap.end()) MCParticleToNhitMap[mcp]=1; + else MCParticleToNhitMap[mcp]++; + } + } //End loop tracker hits. + + for(auto& iter : MCParticleToNhitMap){ + auto rel = relCol->create(); + rel.setRec(iTrk); + rel.setSim(iter.first); + rel.setWeight( (double)iter.second/nHit ); + } + + } //End loop tracks. +} + + + diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h index cfc804879..e7d6c8791 100755 --- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h +++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.h @@ -23,6 +23,8 @@ #include "edm4hep/SimTrackerHitCollection.h" #include "edm4hep/TrackerHitCollection.h" #include "edm4hep/TrackCollection.h" +#include "edm4hep/MCRecoTrackerAssociationCollection.h" +#include "edm4hep/MCRecoTrackParticleAssociationCollection.h" #include #include @@ -297,6 +299,8 @@ class FullLDCTrackingAlg : public GaudiAlgorithm { void AssignOuterHitsToTracks(TrackerHitExtendedVec hitVec, float dcut, int refit); void CreateExtrapolations(); + + void CreateTrackAssociation(edm4hep::TrackCollection* trkCol, edm4hep::MCRecoTrackParticleAssociationCollection* relCol); void CleanUpExtrapolations(); @@ -511,10 +515,17 @@ class FullLDCTrackingAlg : public GaudiAlgorithm { DataHandle _inFTDRawColHdl{"FTDStripTrackerHits", Gaudi::DataHandle::Reader, this}; //DataHandle _inVTXRawColHdl{"VXDCollection", Gaudi::DataHandle::Reader, this}; + DataHandle _VTXHitRelColHdl{"VTXHitRelCol", Gaudi::DataHandle::Reader, this}; + DataHandle _SITHitRelColHdl{"SITHitRelCol", Gaudi::DataHandle::Reader, this}; + DataHandle _SETHitRelColHdl{"SETHitRelCol", Gaudi::DataHandle::Reader, this}; + DataHandle _FTDHitRelColHdl{"FTDHitRelCol", Gaudi::DataHandle::Reader, this}; + DataHandle _TPCHitRelColHdl{"TPCHitRelCol", Gaudi::DataHandle::Reader, this}; + DataHandle _TPCTrackColHdl{"ClupatraTracks", Gaudi::DataHandle::Reader, this}; DataHandle _SiTrackColHdl{"SiTracks", Gaudi::DataHandle::Reader, this}; DataHandle _OutputTrackColHdl{"MarlinTrkTracks", Gaudi::DataHandle::Writer, this}; + DataHandle _OutputTrkRelColHdl{"MarlinTrkAssociation", Gaudi::DataHandle::Writer, this}; NTuple::Tuple* m_tuple; NTuple::Item m_timeTotal; From d8d3f3b1738d6c77d9d68dc937335704b15be35b Mon Sep 17 00:00:00 2001 From: Fangyi Guo Date: Wed, 27 Dec 2023 14:56:57 +0800 Subject: [PATCH 2/2] minor: debug --- .../src/FullLDCTracking/FullLDCTrackingAlg.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp index 18de4bf53..2964b16dc 100755 --- a/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp +++ b/Reconstruction/Tracking/src/FullLDCTracking/FullLDCTrackingAlg.cpp @@ -529,11 +529,11 @@ void FullLDCTrackingAlg::AddTrackColToEvt(TrackExtendedVec & trkVec, edm4hep::Tr float z0TrkCand = trkCand->getZ0(); // float phi0TrkCand = trkCand->getPhi(); // FIXME, fucd - int nhits_in_vxd = track.getSubDetectorHitNumbers(0); - int nhits_in_ftd = track.getSubDetectorHitNumbers(1); - int nhits_in_sit = track.getSubDetectorHitNumbers(2); - int nhits_in_tpc = track.getSubDetectorHitNumbers(3); - int nhits_in_set = track.getSubDetectorHitNumbers(4); + int nhits_in_vxd = track.getSubdetectorHitNumbers(0); + int nhits_in_ftd = track.getSubdetectorHitNumbers(1); + int nhits_in_sit = track.getSubdetectorHitNumbers(2); + int nhits_in_tpc = track.getSubdetectorHitNumbers(3); + int nhits_in_set = track.getSubdetectorHitNumbers(4); //int nhits_in_vxd = Track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::VXD - 2 ]; //int nhits_in_ftd = Track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::FTD - 2 ]; //int nhits_in_sit = Track->subdetectorHitNumbers()[ 2 * lcio::ILDDetID::SIT - 2 ];