diff --git a/Digitization/DigiSimple/CMakeLists.txt b/Digitization/DigiSimple/CMakeLists.txt
index dfebaaea1daebd39f3176e85452cf74166ba9623..ee035a9d704412ad0526eeceaaa65c675c2836b0 100644
--- a/Digitization/DigiSimple/CMakeLists.txt
+++ b/Digitization/DigiSimple/CMakeLists.txt
@@ -5,6 +5,8 @@ gaudi_add_module(SimpleDigi
src/voxel.cpp
src/CylinderDigiAlg.cpp
src/SiTrackerDigiAlg.cpp
+ src/TPCPixelDigiAlg.cpp
+ src/TPCPixelDigiTool.cpp
src/SmearDigiTool.cpp
LINK GearSvc
DigiTool
diff --git a/Digitization/DigiSimple/src/TPCPixelDigiAlg.cpp b/Digitization/DigiSimple/src/TPCPixelDigiAlg.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f819fee2476a2025f03bcb0de2609279511d6bbf
--- /dev/null
+++ b/Digitization/DigiSimple/src/TPCPixelDigiAlg.cpp
@@ -0,0 +1,70 @@
+#include "TPCPixelDigiAlg.h"
+
+
+#include "edm4hep/Vector3f.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+#include "DD4hep/DD4hepUnits.h"
+#include "DDRec/Vector3D.h"
+#include "DDRec/ISurface.h"
+
+#include "GaudiKernel/INTupleSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/IRndmGen.h"
+#include "GaudiKernel/IRndmGenSvc.h"
+#include "GaudiKernel/RndmGenerators.h"
+
+DECLARE_COMPONENT(TPCPixelDigiAlg)
+
+TPCPixelDigiAlg::TPCPixelDigiAlg(const std::string& name, ISvcLocator* svcLoc)
+: GaudiAlgorithm(name, svcLoc) {
+ // Input collections
+ declareProperty("SimTrackHitCollection", m_inputColHdls, "Handle of the Input SimTrackerHit collection");
+
+
+ // Output collections
+ declareProperty("TPCTrackerHits", m_outputColHdls, "Handle of the output TrackerHit collection");
+ declareProperty("TPCTrackerHitAss", m_assColHdls, "Handle of the Association collection between SimTrackerHit and TrackerHit");
+}
+
+StatusCode TPCPixelDigiAlg::initialize() {
+ m_digiTool = ToolHandle<IDigiTool>(m_digiToolName.value());
+
+ info() << "DigiTool " << m_digiTool.typeAndName() << " found" << endmsg;
+
+ info() << "TPCPixelDigiAlg::initialized" << endmsg;
+ return GaudiAlgorithm::initialize();
+}
+
+
+StatusCode TPCPixelDigiAlg::execute(){
+ auto trkCol = m_outputColHdls.createAndPut();
+ auto assCol = m_assColHdls.createAndPut();
+
+ const edm4hep::SimTrackerHitCollection* simCol = nullptr;
+ try {
+ simCol = m_inputColHdls.get();
+ }
+ catch(GaudiException &e){
+ debug() << "Collection " << m_inputColHdls.fullKey() << " test is unavailable in event " << m_nEvt << endmsg;
+ return StatusCode::SUCCESS;
+ }
+ if (simCol->size() == 0) return StatusCode::SUCCESS;
+ debug() << m_inputColHdls.fullKey() << " has cluser "<< simCol->size() << endmsg;
+
+ m_digiTool->Call(simCol, trkCol, assCol);
+
+ debug() << "Created " << trkCol->size() << " hits, "<<simCol->size()<<" clusters"
+ << simCol->size() - trkCol->size() << " hits got dismissed for being out of boundary"
+ << endmsg;
+
+ m_nEvt++;
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode TPCPixelDigiAlg::finalize(){
+ info() << "Processed " << m_nEvt << " events " << endmsg;
+ return GaudiAlgorithm::finalize();
+}
diff --git a/Digitization/DigiSimple/src/TPCPixelDigiAlg.h b/Digitization/DigiSimple/src/TPCPixelDigiAlg.h
new file mode 100644
index 0000000000000000000000000000000000000000..abfc71ec528b4e6d36e04ab34d3aaca0bfb58501
--- /dev/null
+++ b/Digitization/DigiSimple/src/TPCPixelDigiAlg.h
@@ -0,0 +1,38 @@
+#ifndef TPCPixelDigiAlg_H
+#define TPCPixelDigiAlg_H
+
+#include "k4FWCore/DataHandle.h"
+#include "GaudiKernel/NTuple.h"
+#include "GaudiAlg/GaudiAlgorithm.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/SimPrimaryIonizationClusterCollection.h"
+#include "edm4hep/SimPrimaryIonizationCluster.h"
+
+#include "DigiTool/IDigiTool.h"
+
+class TPCPixelDigiAlg : public GaudiAlgorithm{
+ public:
+
+ TPCPixelDigiAlg(const std::string& name, ISvcLocator* svcLoc);
+
+ virtual StatusCode initialize();
+ virtual StatusCode execute();
+ virtual StatusCode finalize();
+
+protected:
+ // Input collections
+ // DataHandle<edm4hep::SimPrimaryIonizationClusterCollection> m_inputColHdls{"SimPrimaryIonizationClusterCollection", Gaudi::DataHandle::Reader, this};
+ // DataHandle<edm4hep::SimPrimaryIonizationClusterCollection> m_inputIonClusterCol{"SimPrimaryIonizationClusterCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::SimTrackerHitCollection> m_inputColHdls{"TPCCollection", Gaudi::DataHandle::Reader, this};
+ // Output collections
+ DataHandle<edm4hep::TrackerHitCollection> m_outputColHdls{"TPCTrackerHits", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::MCRecoTrackerAssociationCollection> m_assColHdls{"TPCTrackerHitsAss", Gaudi::DataHandle::Writer, this};
+
+ Gaudi::Property<std::string> m_digiToolName{this, "DigiTool", "TPCPixelDigiTool"};
+
+ ToolHandle<IDigiTool> m_digiTool;
+ int m_nEvt=0;
+};
+#endif
diff --git a/Digitization/DigiSimple/src/TPCPixelDigiTool.cpp b/Digitization/DigiSimple/src/TPCPixelDigiTool.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9bb8d976a4197e90bd4f6166e984c14de24bd6e2
--- /dev/null
+++ b/Digitization/DigiSimple/src/TPCPixelDigiTool.cpp
@@ -0,0 +1,411 @@
+#include "TPCPixelDigiTool.h"
+#include "k4FWCore/DataHandle.h"
+#include "DataHelper/TrackerHitHelper.h"
+#include "DetIdentifier/CEPCConf.h"
+#include "DetInterface/IGeomSvc.h"
+#include "DataHelper/Circle.h"
+
+#include "edm4hep/Vector3f.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+#include "DD4hep/DD4hepUnits.h"
+#include "DDRec/ISurface.h"
+
+#include "GaudiKernel/INTupleSvc.h"
+#include "GaudiKernel/MsgStream.h"
+#include "GaudiKernel/IRndmGen.h"
+#include "GaudiKernel/RndmGenerators.h"
+#include <gear/GEAR.h>
+#include <gear/TPCParameters.h>
+
+#include "CLHEP/Vector/ThreeVector.h"
+#include "CLHEP/Units/SystemOfUnits.h"
+#include <math.h>
+
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <map>
+#include <cmath>
+#include <algorithm>
+#include <array>
+
+// #include "DataHelper/Circle.h"
+#include "DataHelper/SimpleHelix.h"
+#include "DataHelper/LCCylinder.h"
+#include "EventSeeder/IEventSeeder.h"
+#include <stdexcept>
+#include "GearSvc/IGearSvc.h"
+#include <UTIL/BitField64.h>
+#include <edm4hep/MCParticle.h>
+#include "edm4hep/SimTrackerHit.h"
+
+#include "CLHEP/Vector/TwoVector.h"
+#include "UTIL/ILDConf.h"
+
+#define TRKHITNCOVMATRIX 6
+// using namespace lcio ;
+
+DECLARE_COMPONENT(TPCPixelDigiTool)
+
+StatusCode TPCPixelDigiTool::initialize()
+{
+
+ m_geosvc = service<IGeomSvc>("GeomSvc");
+ if (!m_geosvc)
+ {
+ error() << "Failed to get the GeomSvc" << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+
+ debug() << "Readout name: " << m_readoutName << endmsg;
+ m_decoder = m_geosvc->getDecoder(m_readoutName);
+ if (!m_decoder)
+ {
+ error() << "Failed to get the decoder. " << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ _gear = service<IGearSvc>("GearSvc");
+ if (!_gear)
+ {
+ error() << "Failed to find GearSvc ..." << endmsg;
+ return StatusCode::FAILURE;
+ }
+ _GEAR = _gear->getGearMgr();
+ _f1 = new TF1("fp", "[0]*pow(1+[1],1+[1])*pow(x/[2],[1])*exp(-(1+[1])*x/[2])/TMath::Gamma(1+[1])", 0, 10000);
+ // _f1->SetParameters(3802, 0.487, 2092);
+ _f1->SetParameters(_polyParas[0], _polyParas[1], _polyParas[2]);
+ x_diffusion_p0=0.;
+ y_diffusion_p0=0.;
+ x_diffusion_p1=0.;
+ y_diffusion_p1=0.;
+ t_diffusion_p0 = 0.;
+ t_diffusion_p1 = 0.;
+
+ switch (_magnetic) {
+ case 2:
+ x_diffusion_p0 = 0.001992;
+ x_diffusion_p1 = 0.004387;
+ y_diffusion_p0 = 0.001877;
+ y_diffusion_p1 = 0.004398;
+ t_diffusion_p0 = 0.0002519;
+ t_diffusion_p1 = 0.00332;
+ tpcZRes = tpcTRes * 80; // drift speed 80um/us
+ break;
+
+ case 3:
+ x_diffusion_p0 = 0.001403;
+ x_diffusion_p1 = 0.003011;
+ y_diffusion_p0 = 0.001604;
+ y_diffusion_p1 = 0.003042;
+ t_diffusion_p0 = 5.41279e-05;
+ t_diffusion_p1 = 3.28842e-03;
+ tpcZRes = tpcTRes * 80; // drift speed 80um/us
+ break;
+
+ default:
+ x_diffusion_p0 = 0.0;
+ x_diffusion_p1 = 0.0;
+ y_diffusion_p0 = 0.0;
+ y_diffusion_p1 = 0.0;
+ t_diffusion_p0 = 5.41279e-05;
+ t_diffusion_p1 = 3.28842e-03;
+ tpcZRes = tpcTRes * 80; // drift speed 80um/us
+ break;
+}
+
+
+
+ info() << "initialized" << endmsg;
+ return StatusCode::SUCCESS;
+}
+
+StatusCode TPCPixelDigiTool::Call(const edm4hep::SimTrackerHitCollection *simCol, edm4hep::TrackerHitCollection *hitCol,
+ edm4hep::MCRecoTrackerAssociationCollection *assCol)
+{
+
+ edm4hep::SimTrackerHit simhit;
+ StatusCode sc = Call(simhit, hitCol, assCol);
+ if (sc.isFailure())
+ return sc;
+ return StatusCode::SUCCESS;
+}
+
+StatusCode TPCPixelDigiTool::Call(edm4hep::SimTrackerHit simhit, edm4hep::TrackerHitCollection *hitCol, edm4hep::MCRecoTrackerAssociationCollection *assCol)
+{
+ DataHandle<edm4hep::SimPrimaryIonizationClusterCollection> m_inputIonClusterCol{"SimPrimaryIonizationClusterCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::SimTrackerHitCollection> m_inputSimCol{m_readoutName, Gaudi::DataHandle::Reader, this};
+ UTIL::BitField64 *_cellid_encoder = new UTIL::BitField64(CEPCConf::DetEncoderString::getStringRepresentation());
+ float weight = 1.0;
+
+ // double bReso = _diffRPhi * _diffRPhi;
+ const gear::TPCParameters &gearTPC = _GEAR->getTPCParameters();
+ const gear::PadRowLayout2D &padLayout = gearTPC.getPadLayout();
+ // gear::Vector2D padCoord;
+ const edm4hep::SimTrackerHitCollection *simCol = nullptr;
+ try
+ {
+ simCol = m_inputSimCol.get();
+ }
+ catch (GaudiException &e)
+ {
+ debug() << "Collection " << m_inputSimCol.fullKey() << " is unavailable in event " << endmsg;
+ return StatusCode::SUCCESS;
+ }
+
+ const edm4hep::SimPrimaryIonizationClusterCollection *clusterCol = nullptr;
+ try
+ {
+ clusterCol = m_inputIonClusterCol.get();
+ }
+ catch (GaudiException &e)
+ {
+ debug() << "Collection " << m_inputIonClusterCol.fullKey() << " is unavailable in event " << endmsg;
+ return StatusCode::SUCCESS;
+ }
+ std::map<int, edm4hep::SimTrackerHit> simhitMap; // A map to save the sim tracker hit, key is TPC layer, value is 3D position of simhit
+ std::map<int, std::tuple<double, double, double>> digihitMap; // A map to save the XYZ resolutions and Z position(key is pad index, values are z and rphi resolutions and Z position)
+ std::map<int, int> tpcPixelQ; // A map to save number of electrons in same pixel (key is pixel index, value is number of eles)
+ std::map<int, int> tpcPixelFlag; // A map to mark if a pixel already has a hit (key is pixel index), only used in merging case
+ std::map<int, std::vector<int>> tpcPixelQMap;
+ std::map<int, std::vector<std::tuple<float, float>>> TimeZMap; // A map to save time and z direction infos
+
+ for (auto sim : *simCol)
+ { // loop all sim tracker hits, save them in a map, key is TPC layer, value is 3D position of simhit
+ CLHEP::Hep3Vector SimPoint(sim.getPosition()[0], sim.getPosition()[1], sim.getPosition()[2]);
+ int simPadIndex = padLayout.getNearestPad(SimPoint.perp(), SimPoint.phi());
+ int iRowSimHit = padLayout.getRowNumber(simPadIndex);
+ simhitMap[iRowSimHit] = sim;
+ }
+
+ int nclu = 0;
+
+
+ int padIndex;
+ int iRowHit;
+ int poylaRandom;
+ double x_amp_diffusion;
+ double y_amp_diffusion;
+ CLHEP::Hep3Vector point;
+ CLHEP::Hep3Vector thisPoint;
+ gear::Vector2D padCoord;
+ edm4hep::Vector3d posAfterSmear;
+ edm4hep::Vector3d pos;
+
+ int nIonEle;
+ debug() << " Total clusters: " << clusterCol->size() << endmsg;
+ // for (int i=0; i < _nPads;++i){
+ // padCoord=padLayout.getPadCenter(int(i));
+ // std::cout<<" check padxy "<<padCoord[0] * cos(padCoord[1])<<" "<<padCoord[0] * sin(padCoord[1])<<std::endl;
+ // std::cout<<" check padrphi "<<padCoord[0] <<" "<<padCoord[1] <<std::endl;
+ // }
+
+
+ for (edm4hep::SimPrimaryIonizationCluster cluster : *clusterCol)
+ { // Loop all clusters
+ // break;
+ debug() << " start cluster: " << nclu << " pos:" << cluster.getPosition() << " size: " << cluster.electronPosition_size() << endmsg;
+ nclu += 1;
+ nIonEle = 0;
+
+ for (std::size_t i = 0; i < cluster.electronPosition_size(); ++i)
+ { // loop all electrons in each clusters
+
+ pos = *(cluster.electronPosition_begin() + i);
+ float eleTime = *(cluster.electronTime_begin() + i) / 1000.;//time in us
+
+ if (_MaxIonEles > 0 && nIonEle > _MaxIonEles)
+ {
+ debug() << "Maximum eles achived, break the loop" << endmsg;
+ break;
+ }
+ driftLength = gearTPC.getMaxDriftLength() - (fabs(pos.z));
+ tpcTRes = (t_diffusion_p0 +t_diffusion_p1 * sqrt(driftLength / 10.));
+
+ tpcXRes = x_diffusion_p0 + x_diffusion_p1 * sqrt(driftLength);
+ tpcYRes = y_diffusion_p0 + y_diffusion_p1 * sqrt(driftLength);
+ double randx = gRandom->Gaus(pos.x, tpcXRes);
+ double randy = gRandom->Gaus(pos.y, tpcYRes);
+ // double randz = gRandom->Gaus(pos.z, tpcZRes);
+ double randt = eleTime + gRandom->Gaus(driftLength / 80, tpcTRes); // unit in us
+ double randz = randt*80;
+ //
+ point.set(randx, randy, randz);
+
+ if (fabs(point.z()) > gearTPC.getMaxDriftLength())
+ point.setZ((fabs(point.z()) / point.z()) * gearTPC.getMaxDriftLength());
+ debug() << "resolution, x: " << tpcXRes << " y: " << tpcYRes << " z: " << tpcZRes << endmsg;
+ debug() << "Max TPC length: " << gearTPC.getMaxDriftLength() << " posz: " << fabs(pos.z) << " driftLength: " << driftLength << endmsg;
+ debug() << "smear:" << pos.x << "," << pos.y << "," << pos.z << " -> " << point << endmsg;
+ poylaRandom = _f1->GetRandom();
+ // poylaRandom = 1;
+ debug() << "polyRandom: " << poylaRandom << endmsg;
+ for (int j = 0; j < poylaRandom; ++j) // Do amplication
+ {
+
+ x_amp_diffusion = gRandom->Gaus(point.x(), _ampDiffusion);
+ y_amp_diffusion = gRandom->Gaus(point.y(), _ampDiffusion);
+ thisPoint.set(x_amp_diffusion, y_amp_diffusion, point.z());
+
+ // following lines are the method to find the pixel center
+ padIndex = padLayout.getNearestPad(thisPoint.perp(), thisPoint.phi());
+ padCoord = padLayout.getPadCenter(int(padIndex));
+ posAfterSmear.x = padCoord[0] * cos(padCoord[1]);
+ posAfterSmear.y = padCoord[0] * sin(padCoord[1]);
+ posAfterSmear.z = point.z();
+
+ digihitMap[padIndex] = std::make_tuple(tpcXRes, tpcYRes, tpcZRes);
+ if (tpcPixelQ.find(padIndex) != tpcPixelQ.end())
+ {
+ tpcPixelQ[padIndex] += 1;
+ }
+ else
+ {
+ tpcPixelQ[padIndex] = 1;
+ }
+ if (TimeZMap[padIndex].empty())
+ {
+ TimeZMap[padIndex].push_back(std::make_tuple(randt,point.z()));
+ }
+ else
+ {
+ if (!_saveAllTime)
+ {
+
+ if (randt > std::get<0>(TimeZMap[padIndex].back()) + _deadTime)
+ {
+ info() << "now the size of pad " << padIndex << " is " << TimeZMap[padIndex].size() << " Q is " << tpcPixelQ[padIndex] << endmsg;
+ info() << "now time is: " << randt << " previous time is: " << std::get<0>(TimeZMap[padIndex].back()) << endmsg;
+
+ TimeZMap[padIndex].push_back(std::make_tuple(randt,point.z()));
+ tpcPixelQMap[padIndex].push_back(tpcPixelQ[padIndex]);
+ tpcPixelQ[padIndex]=0;
+ }
+ }
+ else
+ {
+ auto it = TimeZMap[padIndex].rbegin();
+ bool found = false;
+ for (; it != TimeZMap[padIndex].rend(); ++it)
+ {
+ float diff = std::abs(std::get<0>(*it) - randt);
+ if (diff < _minimumTime)
+ {
+ found = true;
+ break;
+ }
+ }
+
+ if (!found)
+ {
+ TimeZMap[padIndex].emplace_back(randt,point.z());
+ tpcPixelQMap[padIndex].push_back(tpcPixelQ[padIndex]);
+ tpcPixelQ[padIndex]=0;
+ }
+ }
+ }
+
+ }
+ nIonEle += 1;
+ }
+ }
+ debug() << "create " << digihitMap.size() << " hits in different pixels before merging" << endmsg;
+ for (auto digiHit = digihitMap.begin(); digiHit != digihitMap.end(); ++digiHit)
+ {
+ padIndex = digiHit->first;
+ std::vector<std::tuple<float, float>> TimeZVec = TimeZMap[padIndex];
+ tpcPixelQMap[padIndex].push_back(tpcPixelQ[padIndex]);
+ std::vector<int> QVec = tpcPixelQMap[padIndex];
+
+ tpcXRes = std::get<0>(digiHit->second);
+ tpcYRes = std::get<1>(digiHit->second);
+ tpcZRes = std::get<2>(digiHit->second);
+ padCoord = padLayout.getPadCenter(int(digiHit->first));
+
+ double r = padCoord[0];
+ double phi = padCoord[1];
+ padIndex = digiHit->first;
+ iRowHit = padLayout.getRowNumber(padIndex);
+ debug() << "pad center: " << r * cos(phi) << " " << r * sin(phi) << endmsg;
+ debug() << "pad r: " << r << " phi: " << phi << " Q: " << tpcPixelQ[padIndex] << endmsg;
+ posAfterSmear.x = (r * cos(phi));
+ posAfterSmear.y = (r * sin(phi));
+
+ for (size_t i = 0; i < TimeZVec.size(); ++i) {
+ posAfterSmear.z = std::get<1>(TimeZVec[i]);
+ auto it = simhitMap.find(iRowHit);
+ if (it != simhitMap.end() && !_skipAss)
+ {
+
+ simhit = it->second;
+ auto outhit = hitCol->create();
+ auto ass = assCol->create();
+ outhit.setPosition(posAfterSmear);
+ (*_cellid_encoder)[lcio::ILDCellID0::subdet] = lcio::ILDDetID::TPC;
+ (*_cellid_encoder)[lcio::ILDCellID0::layer] = iRowHit;
+ (*_cellid_encoder)[lcio::ILDCellID0::module] = 0;
+
+ (*_cellid_encoder)[lcio::ILDCellID0::side] = lcio::ILDDetID::barrel;
+ outhit.setCellID(_cellid_encoder->lowWord());
+ std::array<float, TRKHITNCOVMATRIX> covMat = {float(tpcXRes * tpcXRes),
+ float(tpcYRes * tpcXRes),
+ float(tpcYRes * tpcYRes),
+ float(tpcZRes * tpcXRes),
+ float(tpcZRes * tpcYRes),
+ float(tpcZRes * tpcZRes)};
+
+ outhit.setCovMatrix(covMat);
+ outhit.setEDep(tpcPixelQMap[padIndex][i]);
+ outhit.setTime(std::get<0>(TimeZVec[i]));
+ info() << "out hit, x: " << outhit.getPosition()[0] << " y: " << outhit.getPosition()[1] << " z: " << outhit.getPosition()[2] << " Q: " << outhit.getEDep() << " Time: "<< outhit.getTime() << endmsg;
+ ass.setSim(simhit);
+ ass.setRec(outhit);
+ ass.setWeight(weight);
+ }
+ else if (_skipAss)// do not match the TPCCollection and TPCTrackerHits, all the TPCTrackerHitsAss has been set to 0,0,0
+ {
+ debug() << " skip ass" << endmsg;
+ debug() << " SimTrackerHit " << simhit.getPosition()[0] << endmsg;
+
+ auto outhit = hitCol->create();
+ auto ass = assCol->create();
+ outhit.setPosition(posAfterSmear);
+ (*_cellid_encoder)[lcio::ILDCellID0::subdet] = lcio::ILDDetID::TPC;
+ (*_cellid_encoder)[lcio::ILDCellID0::layer] = iRowHit;
+ (*_cellid_encoder)[lcio::ILDCellID0::module] = 0;
+
+ (*_cellid_encoder)[lcio::ILDCellID0::side] = lcio::ILDDetID::barrel;
+ outhit.setCellID(_cellid_encoder->lowWord());
+ std::array<float, TRKHITNCOVMATRIX> covMat = {float(tpcXRes * tpcXRes),
+ float(tpcYRes * tpcXRes),
+ float(tpcYRes * tpcYRes),
+ float(tpcZRes * tpcXRes),
+ float(tpcZRes * tpcYRes),
+ float(tpcZRes * tpcZRes)};
+
+ outhit.setCovMatrix(covMat);
+ outhit.setEDep(tpcPixelQMap[padIndex][i]);
+ outhit.setTime(std::get<0>(TimeZVec[i]));
+
+ ass.setSim(simhit);
+ ass.setRec(outhit);
+ ass.setWeight(weight);
+ info() << "out hit, x: " << outhit.getPosition()[0] << " y: " << outhit.getPosition()[1] << " z: " << outhit.getPosition()[2] << " Q: " << outhit.getEDep() << " Time: "<< outhit.getTime() << endmsg;
+ }
+ else
+ {
+ debug() << " no matches" << endmsg;
+ }
+ }
+ }
+ return StatusCode::SUCCESS;
+}
+
+StatusCode TPCPixelDigiTool::finalize()
+{
+ StatusCode sc;
+ return sc;
+}
diff --git a/Digitization/DigiSimple/src/TPCPixelDigiTool.h b/Digitization/DigiSimple/src/TPCPixelDigiTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..cd43aa065200a97f18eee2ef7d0e1129dfba7961
--- /dev/null
+++ b/Digitization/DigiSimple/src/TPCPixelDigiTool.h
@@ -0,0 +1,76 @@
+#ifndef TPCPixelDigiTool_h
+#define TPCPixelDigiTool_h
+#include "CLHEP/Vector/TwoVector.h"
+#include "GaudiKernel/AlgTool.h"
+#include "GaudiKernel/IRndmGenSvc.h"
+#include "TF1.h"
+#include "DigiTool/IDigiTool.h"
+#include "DetInterface/IGeomSvc.h"
+#include "GearSvc/IGearSvc.h"
+#include "edm4hep/SimTrackerHitCollection.h"
+#include "edm4hep/TrackerHitCollection.h"
+#include "edm4hep/MCRecoTrackerAssociationCollection.h"
+#include "edm4hep/SimPrimaryIonizationClusterCollection.h"
+#include "edm4hep/SimPrimaryIonizationCluster.h"
+#include <gear/GEAR.h>
+#include "DDRec/DetectorData.h"
+#include "DDRec/SurfaceManager.h"
+#include "TRandom3.h"
+
+#include <vector>
+
+class TPCPixelDigiTool : public extends<AlgTool, IDigiTool> {
+ public:
+ using extends::extends;
+
+ virtual StatusCode Call(const edm4hep::SimTrackerHitCollection* simCol, edm4hep::TrackerHitCollection* hitCol,
+ edm4hep::MCRecoTrackerAssociationCollection* assCol) override;
+ virtual StatusCode Call(edm4hep::SimTrackerHit simhit, edm4hep::TrackerHitCollection* hitCol,
+ edm4hep::MCRecoTrackerAssociationCollection* assCol) override;
+
+ StatusCode initialize() override;
+ StatusCode finalize() override;
+ // double getPadPhi( CLHEP::Hep3Vector* thisPointRPhi, CLHEP::Hep3Vector* firstPointRPhi, CLHEP::Hep3Vector* middlePointRPhi, CLHEP::Hep3Vector* lastPointRPhi);
+ // double getPadTheta( CLHEP::Hep3Vector* firstPointRPhi, CLHEP::Hep3Vector* middlePointRPhi, CLHEP::Hep3Vector* lastPointRPhi );
+
+ private:
+
+ Gaudi::Property<std::string> m_detName{this, "DetName", "TPC"};
+ Gaudi::Property<std::string> m_readoutName{this, "Readout", "TPCCollection"};
+
+
+ Gaudi::Property<int> _magnetic{this, "magnetic", 3};// Magnet field 3 Tesla
+ Gaudi::Property<double> _minimumTime{this, "minimumTime", 0.001};// if the time between two electrons less than this _minimumTime(us), then the two electrons will be treated as one
+ Gaudi::Property<double> _deadTime{this, "deadTime", 300};//
+ Gaudi::Property<int> _saveAllTime{this, "saveAllTime", 1};//
+ Gaudi::Property<double> _ampDiffusion{this, "ampDiffusion", 0.0003}; // Amplification factor
+ Gaudi::Property<double> _MaxIonEles{this, "MaxIonEles",-1}; // Max allowed Inoized eles, negetive values mean no limits
+ Gaudi::Property<std::vector<double> > _polyParas{this, "polyParas", {3802., 0.487, 2092}};//amplification polynominal paramters
+ Gaudi::Property<double> _skipAss{this, "skipAss",1}; // Skip the mactching between digi hits and sim tracker hits
+
+
+
+
+ TF1 *_f1;
+ SmartIF<IRndmGenSvc> m_randSvc;
+ SmartIF<IGeomSvc> m_geosvc;
+ SmartIF<IGearSvc> _gear;
+ gear::GearMgr* _GEAR;
+ dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+ const dd4hep::rec::SurfaceMap* m_surfaces;
+ double x_diffusion_p0;
+ double x_diffusion_p1;
+ double y_diffusion_p0;
+ double y_diffusion_p1;
+ double t_diffusion_p0;
+ double t_diffusion_p1;
+ double tpcXRes;
+ double tpcYRes;
+ double tpcZRes;
+ double tpcTRes;
+ double driftLength;
+
+
+
+};
+#endif