From 89520523a517e3216ea7fabdec7ec11f187491ee Mon Sep 17 00:00:00 2001
From: "guofangyi@ihep.ac.cn" <guofangyi@ihep.ac.cn>
Date: Thu, 8 Aug 2024 01:00:49 +0000
Subject: [PATCH] Add ECAL and HCAL digitization
---
Detector/DetCRD/CMakeLists.txt | 6 +
.../DetCRD/scripts/TDR_o1_v01/calodigi.py | 115 ++++
.../DetCRD/src/Muon/Muon_Barrel_v01_01.cpp | 1 +
.../DetCRD/src/Muon/Muon_Endcap_v01_01.cpp | 2 +-
Digitisers/CMakeLists.txt | 1 +
Digitisers/CaloDigi/CMakeLists.txt | 22 +
Digitisers/CaloDigi/src/CaloBar.h | 54 ++
Digitisers/CaloDigi/src/CaloHit.h | 50 ++
Digitisers/CaloDigi/src/EcalDigiAlg.cpp | 625 ++++++++++++++++++
Digitisers/CaloDigi/src/EcalDigiAlg.h | 132 ++++
Digitisers/CaloDigi/src/HcalDigiAlg.cpp | 252 +++++++
Digitisers/CaloDigi/src/HcalDigiAlg.h | 115 ++++
Digitisers/CaloDigi/src/HitStep.h | 26 +
13 files changed, 1400 insertions(+), 1 deletion(-)
create mode 100644 Detector/DetCRD/scripts/TDR_o1_v01/calodigi.py
create mode 100644 Digitisers/CaloDigi/CMakeLists.txt
create mode 100644 Digitisers/CaloDigi/src/CaloBar.h
create mode 100644 Digitisers/CaloDigi/src/CaloHit.h
create mode 100644 Digitisers/CaloDigi/src/EcalDigiAlg.cpp
create mode 100644 Digitisers/CaloDigi/src/EcalDigiAlg.h
create mode 100644 Digitisers/CaloDigi/src/HcalDigiAlg.cpp
create mode 100644 Digitisers/CaloDigi/src/HcalDigiAlg.h
create mode 100644 Digitisers/CaloDigi/src/HitStep.h
diff --git a/Detector/DetCRD/CMakeLists.txt b/Detector/DetCRD/CMakeLists.txt
index 0ee51f93..6f3d945a 100644
--- a/Detector/DetCRD/CMakeLists.txt
+++ b/Detector/DetCRD/CMakeLists.txt
@@ -57,6 +57,12 @@ foreach(detoption TDR_o1_v01 TDR_o1_v02)
)
endforeach()
+add_test(
+ NAME Test_TDR_o1_v01_CaloDigi
+ COMMAND gaudirun.py Detector/DetCRD/scripts/TDR_o1_v01/calodigi.py
+ WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
+ )
+
foreach(detoption TDR_o1_v01)
add_test(
NAME Test_${detoption}_Geo
diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/calodigi.py b/Detector/DetCRD/scripts/TDR_o1_v01/calodigi.py
new file mode 100644
index 00000000..707db432
--- /dev/null
+++ b/Detector/DetCRD/scripts/TDR_o1_v01/calodigi.py
@@ -0,0 +1,115 @@
+#!/usr/bin/env python
+import os
+from Gaudi.Configuration import *
+
+from Configurables import k4DataSvc
+dsvc = k4DataSvc("EventDataSvc", input="rec00.root")
+
+from Configurables import RndmGenSvc, HepRndm__Engine_CLHEP__RanluxEngine_
+seed = [12340]
+# rndmengine = HepRndm__Engine_CLHEP__RanluxEngine_() # The default engine in Gaudi
+rndmengine = HepRndm__Engine_CLHEP__HepJamesRandom_("RndmGenSvc.Engine") # The default engine in Geant4
+rndmengine.SetSingleton = True
+rndmengine.Seeds = seed
+
+rndmgensvc = RndmGenSvc("RndmGenSvc")
+rndmgensvc.Engine = rndmengine.name()
+
+geometry_option = "TDR_o1_v01/TDR_o1_v01.xml"
+
+if not os.getenv("DETCRDROOT"):
+ print("Can't find the geometry. Please setup envvar DETCRDROOT." )
+ sys.exit(-1)
+
+geometry_path = os.path.join(os.getenv("DETCRDROOT"), "compact", geometry_option)
+if not os.path.exists(geometry_path):
+ print("Can't find the compact geometry file: %s"%geometry_path)
+ sys.exit(-1)
+
+from Configurables import GeomSvc
+geosvc = GeomSvc("GeomSvc")
+geosvc.compact = geometry_path
+
+from Configurables import MarlinEvtSeeder
+evtseeder = MarlinEvtSeeder("EventSeeder")
+
+from Configurables import GearSvc
+gearsvc = GearSvc("GearSvc")
+
+from Configurables import TrackSystemSvc
+tracksystemsvc = TrackSystemSvc("TrackSystemSvc")
+
+from Configurables import SimplePIDSvc
+pidsvc = SimplePIDSvc("SimplePIDSvc")
+cepcswdatatop = "/cvmfs/cepcsw.ihep.ac.cn/prototype/releases/data/latest"
+pidsvc.ParFile = os.path.join(cepcswdatatop, "CEPCSWData/offline-data/Service/SimplePIDSvc/data/dNdx_TPC.root")
+
+from Configurables import PodioInput
+podioinput = PodioInput("PodioReader", collections=[
+# "EventHeader",
+ "MCParticle",
+ "EcalBarrelCollection",
+ "EcalBarrelContributionCollection",
+ "HcalBarrelCollection",
+ "HcalBarrelContributionCollection",
+ "CompleteTracks",
+ "CompleteTracksParticleAssociation"
+ ])
+
+########## Digitalization ################
+
+##ECAL##
+from Configurables import EcalDigiAlg
+EcalDigi = EcalDigiAlg("EcalDigiAlg")
+EcalDigi.ReadOutName = "EcalBarrelCollection"
+EcalDigi.SimCaloHitCollection = "EcalBarrelCollection"
+EcalDigi.CaloHitCollection = "ECALBarrel"
+EcalDigi.CaloAssociationCollection = "ECALBarrelAssoCol"
+EcalDigi.CaloMCPAssociationCollection = "ECALBarrelParticleAssoCol"
+EcalDigi.SkipEvt = 0
+EcalDigi.Seed = 2079
+#Digitalization parameters
+EcalDigi.CalibrECAL = 1.
+EcalDigi.AttenuationLength = 7e10
+EcalDigi.TimeResolution = 0.5 #unit: ns
+EcalDigi.ChargeThresholdFrac = 0.05
+EcalDigi.Debug=1
+EcalDigi.WriteNtuple = 0
+EcalDigi.OutFileName = "Digi_ECAL.root"
+#########################################
+
+##HCAL##
+from Configurables import HcalDigiAlg
+HcalDigi = HcalDigiAlg("HcalDigiAlg")
+HcalDigi.ReadOutName = "HcalBarrelCollection"
+HcalDigi.SimCaloHitCollection = "HcalBarrelCollection"
+HcalDigi.CaloHitCollection = "HCALBarrel"
+HcalDigi.CaloAssociationCollection = "HCALBarrelAssoCol"
+HcalDigi.CaloMCPAssociationCollection = "HCALBarrelParticleAssoCol"
+HcalDigi.SkipEvt = 0
+HcalDigi.Seed = 2079
+#Digitalization parameters
+HcalDigi.MIPResponse = 0.01 # 0.5 MeV / MIP
+HcalDigi.MIPThreshold = 0.5 # Unit: MIP
+HcalDigi.CalibrHCAL = 1.
+HcalDigi.Debug=0
+HcalDigi.WriteNtuple = 0
+HcalDigi.OutFileName = "Digi_HCAL.root"
+
+
+# output
+from Configurables import PodioOutput
+out = PodioOutput("outputalg")
+out.filename = "CaloDigi_TDR_o1_v01_00.root"
+out.outputCommands = ["keep *"]
+
+# ApplicationMgr
+from Configurables import ApplicationMgr
+mgr = ApplicationMgr(
+ TopAlg = [podioinput, EcalDigi, HcalDigi, out],
+ EvtSel = 'NONE',
+ EvtMax = 5,
+ ExtSvc = [dsvc, rndmengine, rndmgensvc, geosvc],
+ HistogramPersistency = 'ROOT',
+ OutputLevel = ERROR
+)
diff --git a/Detector/DetCRD/src/Muon/Muon_Barrel_v01_01.cpp b/Detector/DetCRD/src/Muon/Muon_Barrel_v01_01.cpp
index 9e6abc73..6515e976 100644
--- a/Detector/DetCRD/src/Muon/Muon_Barrel_v01_01.cpp
+++ b/Detector/DetCRD/src/Muon/Muon_Barrel_v01_01.cpp
@@ -320,6 +320,7 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector,
}
dd4hep::Transform3D pv(dd4hep::Rotation3D(dd4hep::RotationX(90*dd4hep::degree)),dd4hep::Position(0,0,0));
dd4hep::PlacedVolume phv = motherVol.placeVolume(envelope,pv);
+ phv.addPhysVolID("system",x_det.id());
sdet.setPlacement(phv);
MYDEBUG("create_detector DONE. ");
diff --git a/Detector/DetCRD/src/Muon/Muon_Endcap_v01_01.cpp b/Detector/DetCRD/src/Muon/Muon_Endcap_v01_01.cpp
index acb54e5b..7631b537 100644
--- a/Detector/DetCRD/src/Muon/Muon_Endcap_v01_01.cpp
+++ b/Detector/DetCRD/src/Muon/Muon_Endcap_v01_01.cpp
@@ -294,7 +294,7 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector,
sdetB.setPlacement(pv);
pv = motherVol.placeVolume(assembly);
-// pv.addPhysVolID("system",x_det.id());
+ pv.addPhysVolID("system",x_det.id());
both_endcaps.setPlacement(pv);
both_endcaps.add(sdetA);
both_endcaps.add(sdetB);
diff --git a/Digitisers/CMakeLists.txt b/Digitisers/CMakeLists.txt
index 6f6f2e40..27da8e2e 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(CaloDigi)
diff --git a/Digitisers/CaloDigi/CMakeLists.txt b/Digitisers/CaloDigi/CMakeLists.txt
new file mode 100644
index 00000000..c3b4182c
--- /dev/null
+++ b/Digitisers/CaloDigi/CMakeLists.txt
@@ -0,0 +1,22 @@
+# Modules
+gaudi_add_module(CaloDigi
+ SOURCES src/EcalDigiAlg.cpp
+ SOURCES src/HcalDigiAlg.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 CaloDigi
+ EXPORT CEPCSWTargets
+ RUNTIME DESTINATION "${CMAKE_INSTALL_BINDIR}" COMPONENT bin
+ LIBRARY DESTINATION "${CMAKE_INSTALL_LIBDIR}" COMPONENT shlib
+ COMPONENT dev)
+
diff --git a/Digitisers/CaloDigi/src/CaloBar.h b/Digitisers/CaloDigi/src/CaloBar.h
new file mode 100644
index 00000000..57d226bb
--- /dev/null
+++ b/Digitisers/CaloDigi/src/CaloBar.h
@@ -0,0 +1,54 @@
+#ifndef _CALOBAR
+#define _CALOBAR
+
+#include <DD4hep/Objects.h>
+#include "TVector3.h"
+
+class CaloBar{
+
+public:
+ CaloBar(unsigned long long _cellID, int _system, int _module, int _slayer, int _dlayer, int _stave, int _bar, TVector3 _pos, double _Q1, double _Q2, double _T1, double _T2)
+ : cellID(_cellID), system(_system), module(_module), stave(_stave), dlayer(_dlayer), 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 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 _slayer, int _bar) { system=_system; module=_module; stave=_stave; dlayer=_dlayer; 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 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/CaloDigi/src/CaloHit.h b/Digitisers/CaloDigi/src/CaloHit.h
new file mode 100644
index 00000000..2c6017de
--- /dev/null
+++ b/Digitisers/CaloDigi/src/CaloHit.h
@@ -0,0 +1,50 @@
+#ifndef _CALOHIT
+#define _CALOHIT
+
+#include <DD4hep/Objects.h>
+#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/CaloDigi/src/EcalDigiAlg.cpp b/Digitisers/CaloDigi/src/EcalDigiAlg.cpp
new file mode 100644
index 00000000..bba34045
--- /dev/null
+++ b/Digitisers/CaloDigi/src/EcalDigiAlg.cpp
@@ -0,0 +1,625 @@
+/* -*- 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/LongCrystalBarBarrelCalorimeter32Polygon_v01.
+// TODO: read geometry info automatically.
+#include "EcalDigiAlg.h"
+
+#include "edm4hep/SimCalorimeterHit.h"
+#include "edm4hep/CalorimeterHit.h"
+#include "edm4hep/Vector3f.h"
+#include "edm4hep/Cluster.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+#include <DDRec/CellIDPositionConverter.h>
+
+#include "TVector3.h"
+#include "TRandom3.h"
+#include <math.h>
+#include <cmath>
+#include <iostream>
+#include <algorithm>
+#include <map>
+#include <random>
+
+// #include <fstream>
+// #include <ctime>
+
+#define C 299.79 // unit: mm/ns
+#define PI 3.141592653
+using namespace std;
+using namespace dd4hep;
+using dd4hep::rec::LayeredCalorimeterData;
+using dd4hep::rec::LayeredCalorimeterStruct;
+
+DECLARE_COMPONENT( EcalDigiAlg )
+
+EcalDigiAlg::EcalDigiAlg(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 EcalDigiAlg::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_Truth", &totE_Truth);
+ t_SimBar->Branch("totE_Digi", &totE_Digi);
+ 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_Truth", &m_simBar_Q1_Truth);
+ t_SimBar->Branch("simBar_Q2_Truth", &m_simBar_Q2_Truth);
+ t_SimBar->Branch("simBar_Q1_Digi", &m_simBar_Q1_Digi);
+ t_SimBar->Branch("simBar_Q2_Digi", &m_simBar_Q2_Digi);
+ 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_slayer", &m_simBar_slayer);
+ t_SimBar->Branch("simBar_cellID", &m_simBar_cellID);
+ }
+
+ std::cout<<"EcalDigiAlg::m_scale="<<m_scale<<std::endl;
+ m_geosvc = service<IGeomSvc>("GeomSvc");
+ if ( !m_geosvc ) throw "EcalDigiAlg :Failed to find GeomSvc ...";
+ m_dd4hep = m_geosvc->lcdd();
+ if ( !m_dd4hep ) throw "EcalDigiAlg :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;
+ }
+
+ rndm.SetSeed(_seed);
+ std::cout<<"EcalDigiAlg::initialize"<<std::endl;
+ return GaudiAlgorithm::initialize();
+}
+
+StatusCode EcalDigiAlg::execute()
+{
+ // clock_t yyy_start, yyy_enddigi;
+ // yyy_start = clock(); // 璁板綍寮€濮嬫椂闂�
+
+ if(_nEvt==0) std::cout<<"EcalDigiAlg::execute Start"<<std::endl;
+ std::cout<<"Processing event: "<<_nEvt<<std::endl;
+ if(_nEvt<_Nskip){ _nEvt++; return StatusCode::SUCCESS; }
+
+ Clear();
+
+ const edm4hep::SimCalorimeterHitCollection* SimHitCol = r_SimCaloCol.get();
+ edm4hep::CalorimeterHitCollection* caloVec = w_DigiCaloCol.createAndPut();
+ edm4hep::MCRecoCaloAssociationCollection* caloAssoVec = w_CaloAssociationCol.createAndPut();
+ edm4hep::MCRecoCaloParticleAssociationCollection* caloMCPAssoVec = w_MCPCaloAssociationCol.createAndPut();
+ std::vector<edm4hep::SimCalorimeterHit> m_simhitCol; m_simhitCol.clear();
+ std::vector<CaloBar> 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() <<std::endl;
+
+ totE_Truth=0;
+ totE_Digi=0;
+
+ //Merge input simhit(steps) to real simhit(bar).
+ MergeHits(*SimHitCol, m_simhitCol);
+ if(_Debug>=1) std::cout<<"Finish Hit Merge, with Nhit: "<<m_simhitCol.size()<<std::endl;
+
+ // DetElement beamcal = m_dd4hep->detector("CaloDetector");
+ // LayeredCalorimeterData* bcExtension = beamcal.extension<LayeredCalorimeterData>();
+ // std::vector<LayeredCalorimeterStruct::Layer> layers = bcExtension->layers;
+
+ //Loop in SimHit, digitalize SimHit to DigiBar
+ for(int i=0;i<m_simhitCol.size();i++){
+
+ auto SimHit = m_simhitCol.at(i);
+
+ unsigned long long id = SimHit.getCellID();
+ CaloBar hitbar;
+ hitbar.setcellID( id);
+ hitbar.setcellID( m_decoder->get(id, "system"),
+ m_decoder->get(id, "module"),
+ m_decoder->get(id, "stave"),
+ m_decoder->get(id, "dlayer"),
+ m_decoder->get(id, "slayer"),
+ m_decoder->get(id, "bar"));
+
+ double Lbar = GetBarLength(hitbar); //NOTE: Is fixed with geometry LongCrystalBarBarrelCalorimeter32Polygon_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);
+
+ //printf("in bar #%d: cellID [%d, %d, %d, %d], position (%.3f, %.3f, %.3f), energy %.3f \n", hitbar.getModule(), hitbar.getStave(), hitbar.getDlayer(), hitbar.getSlayer(), hitbar.getBar(),
+ // 10*hitpos.x(), 10*hitpos.y(), 10*hitpos.z(), SimHit.getEnergy() );
+
+ MCParticleToEnergyWeightMap MCPEnMap; MCPEnMap.clear();
+ std::vector<HitStep> DigiLvec; DigiLvec.clear();
+ std::vector<HitStep> DigiRvec; DigiRvec.clear();
+ double totQ1_Truth = 0;
+ double totQ2_Truth = 0;
+ double totQ1_Digi = 0;
+ double totQ2_Digi = 0;
+ double totQ1 = 0;
+ double totQ2 = 0;
+
+ //Loop in all SimHitContribution(G4Step).
+ //if(_Debug>=2) std::cout<<"SimHit contribution size: "<<SimHit.contributions_size()<<std::endl;
+ for(int iCont=0; iCont < SimHit.contributions_size(); ++iCont){
+ auto conb = SimHit.getContributions(iCont);
+ if( !conb.isAvailable() ) { std::cout<<"EcalDigiAlg Can not get SimHitContribution: "<<iCont<<std::endl; continue;}
+
+ double en = conb.getEnergy();
+ if(en == 0) continue;
+
+ auto mcp = conb.getParticle();
+ MCPEnMap[mcp] += en;
+ TVector3 steppos(conb.getStepPosition().x, conb.getStepPosition().y, conb.getStepPosition().z);
+ TVector3 rpos = steppos-hitbar.getPosition();
+ float step_time = conb.getTime(); // yyy: step time
+
+ m_step_x.push_back(steppos.x());
+ m_step_y.push_back(steppos.y());
+ m_step_z.push_back(steppos.z());
+ m_step_t.push_back(step_time); // yyy: push back step time
+ m_step_E.push_back(en);
+ m_stepBar_x.push_back(hitbar.getPosition().x());
+ m_stepBar_y.push_back(hitbar.getPosition().y());
+ m_stepBar_z.push_back(hitbar.getPosition().z());
+
+ if(_Debug>=3){
+ cout<<"Cell Pos: "<<hitbar.getPosition().x()<<'\t'<<hitbar.getPosition().y()<<'\t'<<hitbar.getPosition().z()<<endl;
+ cout<<"step pos: "<<steppos.x()<<'\t'<<steppos.y()<<'\t'<<steppos.z()<<endl;
+ cout<<"Relative pos: "<<rpos.x()<<'\t'<<rpos.y()<<'\t'<<rpos.z()<<endl;
+ cout<<"Cell: "<<hitbar.getModule()<<" "<<hitbar.getDlayer()<<" "<<hitbar.getSlayer()<<endl;
+ }
+
+ //Get digitalized signal(Q1, Q2, T1, T2) from step
+ //Define: 1 is left, 2 is right, clockwise direction in phi.
+
+ int sign=-999;
+ if(hitbar.getSlayer()==1) sign = rpos.z()==0 ? 1 : rpos.z()/fabs(rpos.z());
+ else{
+ int _module = hitbar.getModule();
+ if(_module<=7 || _module>=25) sign = rpos.x()==0 ? 1: rpos.x()/fabs(rpos.x());
+ if(_module>=9 && _module<=23) sign = rpos.x()==0 ? -1:-rpos.x()/fabs(rpos.x());
+ else if(_module==8) sign = rpos.y()==0 ? 1: rpos.y()/fabs(rpos.y());
+ else if(_module==24) sign = rpos.y()==0 ? -1:-rpos.y()/fabs(rpos.y());
+
+ }
+ if(!fabs(sign)) {std::cout<<"ERROR: Wrong bar direction/position!"<<std::endl; continue;}
+
+ // ####### For Charge Digitization #######
+ double Ratio_left = exp(-(Lbar/2 + sign*sqrt(rpos.Mag2()))/Latt) / (exp(-(Lbar/2 + sign*sqrt(rpos.Mag2()))/Latt) + exp(-(Lbar/2 - sign*sqrt(rpos.Mag2()))/Latt));
+ double Ratio_right = 1 - Ratio_left;
+ totQ1_Truth += en*Ratio_left;
+ totQ2_Truth += en*Ratio_right;
+
+ // ####### For Time Digitization #######
+ double Qi_left = en*exp(-(Lbar/2 + sign*sqrt(rpos.Mag2()))/Latt);
+ double Qi_right = en*exp(-(Lbar/2 - sign*sqrt(rpos.Mag2()))/Latt);
+
+ if(_Debug>=3){
+ cout<<Qi_left<<'\t'<<Qi_right<<endl;
+ cout<<Lbar<<'\t'<<sign*sqrt(rpos.Mag2())<<endl;
+ }
+
+ double Ti_left = -1; int looptime=0;
+ while(Ti_left<0){
+ // Ti_left = Tinit + rndm.Gaus(nMat*(Lbar/2 + sign*sqrt(rpos.Mag2()))/C, Tres);
+ Ti_left = 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 "<<iCont<<" can not get a positive left-side time!"<<std::endl; break;}
+ }
+ if(looptime>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 "<<iCont<<" can not get a positive right-side time!"<<std::endl;
+ std::cout<<" Initial time "<<Tinit<<", transport time central value "<<nMat*(Lbar/2 - sign*sqrt(rpos.Mag2()))/C<<std::endl;
+ break;
+ }
+ }
+ if(looptime>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);
+ }
+
+ // #######################################
+ // ####### Ideal Time Digitization #######
+ // #######################################
+
+ //if(_Debug>=2) std::cout<<"Time Digitalize: time at Q >"<<_Qthfrac<<"*totQ"<<std::endl;
+ std::sort(DigiLvec.begin(), DigiLvec.end());
+ std::sort(DigiRvec.begin(), DigiRvec.end());
+ double thQ1=0;
+ double thQ2=0;
+ double thT1, thT2;
+ for(int iCont=0;iCont<DigiLvec.size();iCont++){
+ thQ1 += DigiLvec[iCont].getQ();
+ if(thQ1>totQ1*_Qthfrac){
+ thT1 = DigiLvec[iCont].getT();
+ if(_Debug>=3) std::cout<<"Get T1 at index: "<<iCont<<std::endl;
+ break;
+ }
+ }
+ for(int iCont=0;iCont<DigiRvec.size();iCont++){
+ thQ2 += DigiRvec[iCont].getQ();
+ if(thQ2>totQ2*_Qthfrac){
+ thT2 = DigiRvec[iCont].getT();
+ if(_Debug>=3) std::cout<<"Get T2 at index: "<<iCont<<std::endl;
+ break;
+ }
+ }
+
+ if(_UseRelDigi){
+ // #############################################
+ // ####### Realistic Charge Digitization #######
+ // #############################################
+
+ // double sEcalCryMipLY = gRandom->Gaus(fEcalCryMipLY, 0.1 * fEcalCryMipLY);
+ double sEcalCryMipLY = fEcalCryMipLY;
+
+ //TODO: fEcalMIPEnergy should depends on crystal size.
+ int ScinGen = std::round(gRandom->Poisson((totQ1_Truth+totQ2_Truth)*1000 / fEcalMIPEnergy * sEcalCryMipLY));
+
+ totQ1_Digi = EnergyDigi(ScinGen*totQ1_Truth/(totQ1_Truth+totQ2_Truth), sEcalCryMipLY)/1000;
+ totQ2_Digi = EnergyDigi(ScinGen*totQ2_Truth/(totQ1_Truth+totQ2_Truth), sEcalCryMipLY)/1000;
+ }
+ else{
+ double totQ1_Digi_tmp = totQ1_Truth;
+ double totQ2_Digi_tmp = totQ2_Truth;
+ if( (totQ1_Digi_tmp+totQ2_Digi_tmp)!=0 ){
+ totQ1_Digi = totQ1_Digi_tmp/(totQ1_Digi_tmp+totQ2_Digi_tmp);
+ totQ2_Digi = totQ2_Digi_tmp/(totQ1_Digi_tmp+totQ2_Digi_tmp);
+ }
+ if( totQ1_Digi/fEcalMIPEnergy < fEcalMIP_Thre ) totQ1_Digi = 0;
+ if( totQ2_Digi/fEcalMIPEnergy < fEcalMIP_Thre ) totQ2_Digi = 0;
+ }
+ if(totQ1_Digi==0 && totQ2_Digi==0) continue;
+
+ hitbar.setQ(totQ1_Digi, totQ2_Digi);
+ hitbar.setT(thT1, thT2);
+
+ // ##################################
+ // ####### Some associations #######
+ // ##################################
+
+ //2 hits with double-readout time.
+ edm4hep::Vector3f m_pos(hitbar.getPosition().X(), hitbar.getPosition().Y(), hitbar.getPosition().Z());
+ auto digiHit1 = caloVec->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/SimHit.getEnergy());
+ auto rel_MCP2 = caloMCPAssoVec->create();
+ rel_MCP2.setRec(digiHit2);
+ rel_MCP2.setSim(iter.first);
+ rel_MCP2.setWeight(iter.second/SimHit.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.);
+ //}
+
+ // ########################################
+ // ####### Temp: write into trees. #######
+ // ########################################
+
+ m_barCol.push_back(hitbar);
+ if(hitbar.getQ1()>0 && hitbar.getQ2()>0) totE_Digi+=(hitbar.getQ1()+hitbar.getQ2());
+ if(totQ1_Truth>(fEcalMIPEnergy*fEcalMIP_Thre/1000.) && totQ2_Truth>(fEcalMIPEnergy*fEcalMIP_Thre/1000.)){
+ // cout<<"Truth Energy:"<<totQ1_Truth<<" "<<totQ2_Truth<<endl;
+ totE_Truth+=(totQ1_Truth+totQ2_Truth);
+ }
+
+ 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_Truth.push_back(totQ1_Truth);
+ m_simBar_Q2_Truth.push_back(totQ2_Truth);
+ m_simBar_Q1_Digi.push_back(hitbar.getQ1());
+ m_simBar_Q2_Digi.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_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!"<<std::endl;
+ std::cout<<"Total Truth Energy: "<<totE_Truth<<std::endl;
+ std::cout<<"Total Digi Energy: "<<totE_Digi<<std::endl;
+
+ // yyy_enddigi = clock();
+ // double duration_digi = double(yyy_enddigi - yyy_start) / CLOCKS_PER_SEC;
+ // // 灏嗘椂闂磋緭鍑哄埌txt鏂囦欢涓�
+ // std::ofstream outfile("runtime_ecaldigi.txt", std::ios::app);
+ // outfile << _nEvt << " " << duration_digi << std::endl;
+ // outfile.close();
+
+ _nEvt ++ ;
+ //delete SimHitCol, caloVec, caloAssoVec;
+ m_simhitCol.clear();
+ return StatusCode::SUCCESS;
+}
+
+StatusCode EcalDigiAlg::finalize()
+{
+ if(_writeNtuple){
+ m_wfile->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();
+}
+
+double EcalDigiAlg::EnergyDigi(double ScinGen, double sEcalCryMipLY){
+ Int_t sPix = int(ScinGen);
+
+ // if(sPix/sEcalCryMipLY < fEcalMIP_Thre) return 0;
+ // return sPix/sEcalCryMipLY*fEcalMIPEnergy;
+
+ // return sPix/sEcalCryMipLY*fEcalMIPEnergy;
+
+ // ####### SiPM Saturation #######
+ // sPix = std::round(fEcalSiPMPixels * (1 - TMath::Exp(-sPix / fEcalSiPMPixels)));
+
+ // ################################
+ // ####### ADC Digitization #######
+ // ################################
+
+ Bool_t Use_G1 = kFALSE;
+ Bool_t Use_G2 = kFALSE;
+ Bool_t Use_G3 = kFALSE;
+
+ Double_t sADCMean = sPix * fEcalChargeADCMean;
+ Double_t sADCSigma = std::sqrt(sPix * fEcalChargeADCSigma * fEcalChargeADCSigma + fEcalNoiseADCSigma * fEcalNoiseADCSigma);
+ Int_t sADC = -1;
+ sADC = std::round(gRandom->Gaus(sADCMean, sADCSigma));
+
+ if(sADC <= fADCSwitch){
+ Use_G1 = kTRUE;
+ sADC = std::round(gRandom->Gaus(sADC, fEcalADCError * sADC));
+ Double_t sMIP = sADC / fEcalChargeADCMean / sEcalCryMipLY;
+ if(sMIP < fEcalMIP_Thre) return 0;
+ return sMIP * fEcalMIPEnergy;
+ }
+ else if(sADC > fADCSwitch && int(sADC/fGainRatio_12) <= fADCSwitch){
+ Use_G2 = kTRUE;
+ sADCMean = sPix * fEcalChargeADCMean / fGainRatio_12;
+ sADCSigma = std::sqrt(sPix * fEcalChargeADCSigma / fGainRatio_12 * fEcalChargeADCSigma / fGainRatio_12 + fEcalNoiseADCSigma * fEcalNoiseADCSigma);
+ sADC = std::round(gRandom->Gaus(sADCMean, sADCSigma));
+ sADC = std::round(gRandom->Gaus(sADC, fEcalADCError * sADC));
+ Double_t sMIP = sADC / fEcalChargeADCMean * fGainRatio_12 / sEcalCryMipLY;
+ if(sMIP < fEcalMIP_Thre) return 0;
+ return sMIP * fEcalMIPEnergy;
+ }
+ else if(int(sADC/fGainRatio_12) > fADCSwitch){
+ Use_G3 = kTRUE;
+ sADCMean = sPix * fEcalChargeADCMean / fGainRatio_12 / fGainRatio_23;
+ sADCSigma = std::sqrt(sPix * fEcalChargeADCSigma / fGainRatio_12 / fGainRatio_23 * fEcalChargeADCSigma / fGainRatio_12 / fGainRatio_23 + fEcalNoiseADCSigma * fEcalNoiseADCSigma);
+ sADC = std::round(gRandom->Gaus(sADCMean, sADCSigma));
+ sADC = std::round(gRandom->Gaus(sADC, fEcalADCError * sADC));
+ if (sADC > fADC-1) sADC = fADC-1;
+ Double_t sMIP = sADC / fEcalChargeADCMean * fGainRatio_12 * fGainRatio_23 / sEcalCryMipLY;
+ if(sMIP < fEcalMIP_Thre) return 0;
+ return sMIP * fEcalMIPEnergy;
+ }
+}
+
+StatusCode EcalDigiAlg::MergeHits( const edm4hep::SimCalorimeterHitCollection& m_col, std::vector<edm4hep::SimCalorimeterHit>& m_hits ){
+
+ m_hits.clear();
+ std::vector<edm4hep::MutableSimCalorimeterHit> m_mergedhit;
+ m_mergedhit.clear();
+
+ for(int iter=0; iter<m_col.size(); iter++){
+ edm4hep::SimCalorimeterHit m_step = m_col[iter];
+ if(!m_step.isAvailable()){ cout<<"ERROR HIT!"<<endl; continue;}
+ if(m_step.getEnergy()==0 || m_step.contributions_size()<1) continue;
+ unsigned long long cellid = m_step.getCellID();
+ dd4hep::Position hitpos = m_cellIDConverter->position(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.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 EcalDigiAlg::GetBarLength(CaloBar& bar){
+ //TODO: reading bar length from geosvc.
+ if(bar.getSlayer()==1) return 375.133;
+ else{
+ if(bar.getModule()%2 == 0){
+ return 295.905 + (bar.getDlayer()-1)* 6.13231;
+ }
+ else{
+ return 416.843 - (bar.getDlayer()+1)* 2.25221;
+ }
+
+ }
+}
+
+
+/*
+dd4hep::Position EcalDigiAlg::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 EcalDigiAlg::find(edm4hep::SimCalorimeterHitCollection& m_col, dd4hep::Position& pos){
+ for(int i=0;i<m_col.size();i++){
+ edm4hep::MutableSimCalorimeterHit hit = m_col[i];
+ dd4hep::Position ipos(hit.getPosition().x, hit.getPosition().y, hit.getPosition().z);
+ if(ipos==pos) return hit;
+ }
+ edm4hep::MutableSimCalorimeterHit hit;
+ hit.setCellID(0);
+ return hit;
+}
+*/
+edm4hep::MutableSimCalorimeterHit EcalDigiAlg::find(const std::vector<edm4hep::MutableSimCalorimeterHit>& m_col, unsigned long long& cellid) const{
+ for(int i=0;i<m_col.size();i++){
+ edm4hep::MutableSimCalorimeterHit hit=m_col.at(i);
+ if(hit.getCellID() == cellid) return hit;
+ }
+ edm4hep::MutableSimCalorimeterHit hit ;
+ hit.setCellID(0);
+ return hit;
+}
+
+void EcalDigiAlg::Clear(){
+ totE_Truth = -99;
+ totE_Digi = -99;
+ m_step_x.clear();
+ m_step_y.clear();
+ m_step_z.clear();
+ m_step_t.clear(); // yyy: clear
+ m_step_E.clear();
+ m_stepBar_x.clear();
+ m_stepBar_y.clear();
+ m_stepBar_z.clear();
+ m_step_T1.clear();
+ m_step_T2.clear();
+ m_simBar_x.clear();
+ m_simBar_y.clear();
+ m_simBar_z.clear();
+ m_simBar_T1.clear();
+ m_simBar_T2.clear();
+ m_simBar_Q1_Truth.clear();
+ m_simBar_Q2_Truth.clear();
+ m_simBar_Q1_Digi.clear();
+ m_simBar_Q2_Digi.clear();
+ m_simBar_module.clear();
+ m_simBar_stave.clear();
+ m_simBar_dlayer.clear();
+ m_simBar_slayer.clear();
+ m_simBar_cellID.clear();
+}
diff --git a/Digitisers/CaloDigi/src/EcalDigiAlg.h b/Digitisers/CaloDigi/src/EcalDigiAlg.h
new file mode 100644
index 00000000..058a06bd
--- /dev/null
+++ b/Digitisers/CaloDigi/src/EcalDigiAlg.h
@@ -0,0 +1,132 @@
+#ifndef _ECAL_DIGI_ALG_H
+#define _ECAL_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 "CaloBar.h"
+#include "HitStep.h"
+
+#include <DDRec/DetectorData.h>
+#include <DDRec/CellIDPositionConverter.h>
+#include <DD4hep/Segmentations.h>
+#include "DetInterface/IGeomSvc.h"
+#include "TVector3.h"
+#include "TRandom3.h"
+#include "TFile.h"
+#include "TString.h"
+#include "TH3.h"
+#include "TH1.h"
+
+#include <cstdlib>
+#include "time.h"
+#include <TTimeStamp.h>
+#include <ctime>
+
+#define C 299.79 // unit: mm/ns
+#define PI 3.141592653
+
+class EcalDigiAlg : public GaudiAlgorithm
+{
+
+public:
+
+ EcalDigiAlg(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<edm4hep::SimCalorimeterHit>& m_hits);
+
+ double GetBarLength(CaloBar& bar); //TODO: should read from geom file!
+ edm4hep::MutableSimCalorimeterHit find(const std::vector<edm4hep::MutableSimCalorimeterHit>& m_col, unsigned long long& cellid) const;
+ // double Digitization(double edepCry, double fCrosstalkChannel);
+ double EnergyDigi(double ScinGen, double sEcalCryMipLY);
+
+ void Clear();
+
+protected:
+
+ SmartIF<IGeomSvc> m_geosvc;
+ typedef std::vector<float> FloatVec;
+ typedef std::map<const edm4hep::MCParticle, float> MCParticleToEnergyWeightMap;
+
+ int _nEvt ;
+ float m_length;
+ TRandom3 rndm;
+ TFile* m_wfile;
+ TTree* t_SimCont;
+ TTree* t_SimBar;
+
+ double totE_Truth, totE_Digi;
+ 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_Truth, m_simBar_Q2_Truth, m_simBar_Q1_Digi, m_simBar_Q2_Digi, m_simBar_dlayer, m_simBar_stave, m_simBar_slayer, m_simBar_module;
+ std::vector<unsigned long long> m_simBar_cellID;
+
+
+ dd4hep::rec::CellIDPositionConverter* m_cellIDConverter;
+ dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+ dd4hep::Detector* m_dd4hep;
+
+ Gaudi::Property<float> m_scale{ this, "Scale", 1 };
+
+ // Input collections
+ DataHandle<edm4hep::SimCalorimeterHitCollection> r_SimCaloCol{"SimCaloCol", Gaudi::DataHandle::Reader, this};
+ mutable Gaudi::Property<std::string> _readoutName{this, "ReadOutName", "CaloHitsCollection", "Readout name"};
+ mutable Gaudi::Property<std::string> _filename{this, "OutFileName", "testout.root", "Output file name"};
+ mutable Gaudi::Property<int> _UseRelDigi{this, "UseRealisticDigi", 1, "If use the realistic model"};
+
+ //Input parameters
+ mutable Gaudi::Property<int> _writeNtuple{this, "WriteNtuple", 1, "Write ntuple"};
+ mutable Gaudi::Property<int> _Nskip{this, "SkipEvt", 0, "Skip event"};
+ mutable Gaudi::Property<float> _seed{this, "Seed", 2131, "Random Seed"};
+ mutable Gaudi::Property<int> _Debug{this, "Debug", 0, "Debug level"};
+ mutable Gaudi::Property<float> _Eth {this, "EnergyThreshold", 0.001, "Energy Threshold (/GeV)"};
+ mutable Gaudi::Property<float> r_cali{this, "CalibrECAL", 1, "Calibration coefficients for ECAL"};
+ mutable Gaudi::Property<float> Latt{this, "AttenuationLength", 7000, "Crystal Attenuation Length(mm)"};
+ mutable Gaudi::Property<float> Tres{this, "TimeResolution", 0.1, "Crystal time resolution in one side (ns)"};
+ mutable Gaudi::Property<float> nMat{this, "MatRefractive", 2.15, "Material refractive index of crystal"};
+ mutable Gaudi::Property<float> Tinit{this, "InitalTime", 2, "Start time (ns)"};
+
+ mutable Gaudi::Property<float> _Qthfrac {this, "ChargeThresholdFrac", 0.05, "Charge threshold fraction"};
+
+ mutable Gaudi::Property<int> fADC{this, "ADC", 4096, "Total ADC conuts"};
+ mutable Gaudi::Property<int> fNofGain{this, "NofGain", 3, "Number of gain modes"};
+ mutable Gaudi::Property<int> fADCSwitch{this, "ADCSwitch", 4000, "switching point of different gain mode"};
+ mutable Gaudi::Property<float> fGainRatio_12{this, "GainRatio_12", 15, "Gain-1 over Gain-2"};
+ mutable Gaudi::Property<float> fGainRatio_23{this, "GainRatio_23", 10, "Gain-2 over Gain-3"};
+ mutable Gaudi::Property<float> fEcalCryMipLY{this, "EcalCryMipLY", 100, "Crystal light yield (p.e./MIP)"};
+ mutable Gaudi::Property<float> fEcalMIPEnergy{this, "EcalMIPEnergy", 8.9, "MIP Energy deposit in 1cm BGO (MeV/MIP)"};
+ mutable Gaudi::Property<int> fEcalSiPMPixels{this, "EcalSiPMPixels", 250000, "Pixels number of SiPM"};
+ mutable Gaudi::Property<float> fEcalChargeADCMean{this, "EcalChargeADCMean", 5, "ADC per p.e. for Gain-1 (ADC)"};
+ mutable Gaudi::Property<float> fEcalChargeADCSigma{this, "EcalChargeADCSigma", 2.5, "Sigma of ADC per p.e. for Gain-1 (ADC)"};
+ mutable Gaudi::Property<float> fEcalADCError{this, "EcalADCError", 0.002, "ADC precision"};
+ mutable Gaudi::Property<float> fEcalNoiseADCSigma{this, "EcalNoiseADCSigma", 3, "Sigma of electronic noise (ADC)"};
+ mutable Gaudi::Property<float> fEcalMIP_Thre{this, "EcalMIP_Thre", 0.1, "Energy threshold for single readout end (MIP)"};
+
+ // Output collections
+ DataHandle<edm4hep::CalorimeterHitCollection> w_DigiCaloCol{"DigiCaloCol", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::MCRecoCaloAssociationCollection> w_CaloAssociationCol{"ECALBarrelAssoCol", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection> w_MCPCaloAssociationCol{"ECALBarrelParticleAssoCol", Gaudi::DataHandle::Writer, this};
+};
+
+#endif
diff --git a/Digitisers/CaloDigi/src/HcalDigiAlg.cpp b/Digitisers/CaloDigi/src/HcalDigiAlg.cpp
new file mode 100644
index 00000000..a72a30cb
--- /dev/null
+++ b/Digitisers/CaloDigi/src/HcalDigiAlg.cpp
@@ -0,0 +1,252 @@
+// /* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
+// // Unit in code: mm, ns.
+
+#include "HcalDigiAlg.h"
+
+#include "edm4hep/SimCalorimeterHit.h"
+#include "edm4hep/CalorimeterHit.h"
+#include "edm4hep/Vector3f.h"
+#include "edm4hep/Cluster.h"
+
+#include "DD4hep/Detector.h"
+#include <DD4hep/Objects.h>
+#include <DDRec/CellIDPositionConverter.h>
+
+#include "TVector3.h"
+#include <math.h>
+#include <cmath>
+#include <iostream>
+#include <algorithm>
+#include <map>
+
+// #include <fstream>
+// #include <ctime>
+
+#define C 299.79 // unit: mm/ns
+#define PI 3.141592653
+using namespace std;
+using namespace dd4hep;
+
+DECLARE_COMPONENT( HcalDigiAlg )
+
+HcalDigiAlg::HcalDigiAlg(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 HcalDigiAlg::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("totE", &m_totE);
+ 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_HG", &m_simHit_HG);
+ t_simHit->Branch("simHit_LG", &m_simHit_LG);
+ 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<<"HcalDigiAlg::m_scale="<<m_scale<<std::endl;
+ m_geosvc = service<IGeomSvc>("GeomSvc");
+ if ( !m_geosvc ) throw "HcalDigiAlg :Failed to find GeomSvc ...";
+ dd4hep::Detector* m_dd4hep = m_geosvc->lcdd();
+ if ( !m_dd4hep ) throw "HcalDigiAlg :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;
+ }
+
+
+ rndm.SetSeed(_seed);
+ std::cout<<"HcalDigiAlg::initialize"<<std::endl;
+ return GaudiAlgorithm::initialize();
+}
+
+StatusCode HcalDigiAlg::execute()
+{
+// clock_t yyy_start, yyy_enddigi;
+// yyy_start = clock(); // 璁板綍寮€濮嬫椂闂�
+
+ if(_nEvt==0) std::cout<<"HcalDigiAlg::execute Start"<<std::endl;
+ std::cout<<"Processing event: "<<_nEvt<<std::endl;
+ if(_nEvt<_Nskip){ _nEvt++; return StatusCode::SUCCESS; }
+
+ Clear();
+ m_totE = 0.;
+ const edm4hep::SimCalorimeterHitCollection* SimHitCol = r_SimCaloCol.get();
+
+ edm4hep::CalorimeterHitCollection* caloVec = w_DigiCaloCol.createAndPut();
+ edm4hep::MCRecoCaloAssociationCollection* caloAssoVec = w_CaloAssociationCol.createAndPut();
+ edm4hep::MCRecoCaloParticleAssociationCollection* caloMCPAssoVec = w_MCPCaloAssociationCol.createAndPut();
+
+ 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() <<std::endl;
+
+
+ for(int isim=0; isim<SimHitCol->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();
+
+ double sChargeOutHG = 0;
+ double sChargeOutLG = 0;
+ //Digitization
+ if(_UseRelDigi){
+ // -- Scintillation (Energy -> MIP -> Np.e.)
+ int sPix = gRandom->Poisson(Ehit / _MIPCali * (_MIPADC / _PeADCMean));
+ // -- Tile non-uniformity
+ sPix = sPix * (1.0 + gRandom->Uniform(-_TileRes, _TileRes));
+ // -- SiPM Saturation (Np.e. -> Npixel)
+ sPix = std::round(_Pixel * (1.0 - TMath::Exp(-sPix * 1.0 / _Pixel)));
+ // -- ADC response (Npixel -> ADC)
+ double sChargeMean = sPix * _PeADCMean;
+ double sChargeSigma = sqrt(sPix * _PeADCSigma * _PeADCSigma);
+ double sChargeOut = gRandom->Gaus(sChargeMean, sChargeSigma);
+ // -- (ADC->MIP)
+ sChargeOutHG = sChargeOut + gRandom->Gaus(_BaselineHG, _BaselineSigmaHG);
+ sChargeOutLG = sChargeOut / _HLRatio + gRandom->Gaus(_BaselineLG, _BaselineSigmaLG);
+ sChargeOutHG = std::round(gRandom->Gaus(sChargeOutHG, sChargeOutHG * _ADCError));
+ sChargeOutLG = std::round(gRandom->Gaus(sChargeOutLG, sChargeOutLG * _ADCError));
+ if (sChargeOutLG > _ADCLimit)
+ sChargeOutLG = _ADCLimit;
+ Double_t sMIP = 0;
+ if (sChargeOutHG > _ADCSwitch)
+ {
+ sMIP = (sChargeOutLG - _BaselineLG) * _HLRatio / _MIPADC;
+ sChargeOutHG = _ADCSwitch;
+ }
+ else
+ {
+ sMIP = (sChargeOutHG - _BaselineHG) / _MIPADC;
+ }
+ Ehit = sMIP * _MIPCali;
+ }
+ 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; iConb<simhit.contributions_size(); ++iConb){
+ auto conb = simhit.getContributions(iConb);
+ if(!conb.isAvailable()) continue;
+ if(conb.getEnergy()==0) continue;
+
+ Thit_ave += conb.getTime();
+
+ auto mcp = conb.getParticle();
+ MCPEnMap[mcp] += conb.getEnergy();
+ Ehit_raw += conb.getEnergy();
+ }
+ Thit_ave = Thit_ave/simhit.contributions_size();
+ //Create DigiHit
+ auto digiHit = caloVec->create();
+ 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_totE += digiHit.getEnergy();
+ 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_HG.push_back(sChargeOutHG);
+ m_simHit_LG.push_back(sChargeOutLG);
+ 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 HcalDigiAlg::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 HcalDigiAlg::Clear(){
+ m_totE = -99;
+ m_simHit_x.clear();
+ m_simHit_y.clear();
+ m_simHit_z.clear();
+ m_simHit_E.clear();
+ m_simHit_HG.clear();
+ m_simHit_LG.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/CaloDigi/src/HcalDigiAlg.h b/Digitisers/CaloDigi/src/HcalDigiAlg.h
new file mode 100644
index 00000000..c65f3d5a
--- /dev/null
+++ b/Digitisers/CaloDigi/src/HcalDigiAlg.h
@@ -0,0 +1,115 @@
+#ifndef HCAL_DIGI_ALG_H
+#define 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 <DDRec/DetectorData.h>
+#include <DDRec/CellIDPositionConverter.h>
+#include <DD4hep/Segmentations.h>
+#include "DetInterface/IGeomSvc.h"
+#include "TVector3.h"
+#include "TRandom3.h"
+#include "TFile.h"
+#include "TString.h"
+#include "TH3.h"
+#include "TH1.h"
+
+#include <cstdlib>
+#include "time.h"
+#include <TTimeStamp.h>
+#include <ctime>
+
+#define C 299.79 // unit: mm/ns
+#define PI 3.141592653
+
+class HcalDigiAlg : public GaudiAlgorithm
+{
+
+public:
+
+ HcalDigiAlg(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<IGeomSvc> m_geosvc;
+ typedef std::vector<float> FloatVec;
+ typedef std::map<const edm4hep::MCParticle, float> MCParticleToEnergyWeightMap;
+
+ int _nEvt ;
+ TRandom3 rndm;
+ TFile* m_wfile;
+ TTree* t_simHit;
+
+ double m_totE;
+ 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, m_simHit_HG, m_simHit_LG;
+ std::vector<unsigned long long> m_simHit_cellID;
+ std::vector<int> m_simHit_steps;
+
+
+ dd4hep::rec::CellIDPositionConverter* m_cellIDConverter;
+ dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+
+ Gaudi::Property<float> m_scale{ this, "Scale", 1 };
+
+ // Input collections
+ DataHandle<edm4hep::SimCalorimeterHitCollection> r_SimCaloCol{"SimCaloCol", Gaudi::DataHandle::Reader, this};
+ mutable Gaudi::Property<std::string> _readoutName{this, "ReadOutName", "CaloHitsCollection", "Readout name"};
+ mutable Gaudi::Property<std::string> _filename{this, "OutFileName", "testout.root", "Output file name"};
+
+ //Input parameters
+ mutable Gaudi::Property<int> _writeNtuple{this, "WriteNtuple", 1, "Write ntuple"};
+ mutable Gaudi::Property<int> _Nskip{this, "SkipEvt", 0, "Skip event"};
+ mutable Gaudi::Property<float> _seed{this, "Seed", 2131, "Random Seed"};
+ mutable Gaudi::Property<int> _Debug{this, "Debug", 0, "Debug level"};
+ mutable Gaudi::Property<float> r_cali{this, "CalibrHCAL", 1, "Global calibration coefficients for HCAL"};
+ mutable Gaudi::Property<int> _UseRelDigi{this, "UseRealisticDigi", 1, "If use the realistic model"};
+
+ //add digitization parameters from AHCAL prototype
+ mutable Gaudi::Property<float> _MIPCali{this, "MIPResponse", 0.000461, "MIP response (/GeV)"};
+ mutable Gaudi::Property<float> _Eth_Mip{this, "MIPThreshold", 0.5, "Energy Threshold (/MIP)"};
+ mutable Gaudi::Property<int> _Pixel{this, "SiPMPixel", 7284, "number of SiPM pixels"};
+ mutable Gaudi::Property<float> _ADCError{this, "ADCError", 0.0002, "ADC Error (/ADC)"};
+ mutable Gaudi::Property<float> _MIPADC{this, "MIPADCMean", 345.7, "Mean value of MIP response adc (/ADC)"};
+ mutable Gaudi::Property<float> _TileRes{this, "TileNonUniformity", 0.05, "Non-uniformity level of one tile response"};
+ mutable Gaudi::Property<float> _PeADCMean{this, "PeADCMean", 30.0, "Mean value of single photons adc (/ADC)"};
+ mutable Gaudi::Property<float> _PeADCSigma{this, "PeADCSigma", 3.5, "Sigma of single photons adc (/ADC)"};
+ mutable Gaudi::Property<float> _BaselineHG{this, "ADCBaselineHG", 377.4, "Mean value of HG baseline adc (/ADC)"};
+ mutable Gaudi::Property<float> _BaselineSigmaHG{this, "ADCBaselineSigmaHG", 3.3, "Sigma of HG baseline adc (/ADC)"};
+ mutable Gaudi::Property<float> _BaselineLG{this, "ADCBaselineLG", 373.9, "Mean value of LG baseline adc (/ADC)"};
+ mutable Gaudi::Property<float> _BaselineSigmaLG{this, "ADCBaselineSigmaLG", 2.2, "Sigma of LG baseline adc (/ADC)"};
+ mutable Gaudi::Property<float> _HLRatio{this, "ADCHLRatio", 29.9, "The ratio of HG to LG"};
+ mutable Gaudi::Property<float> _ADCSwitch{this, "ADCSwitch", 2930, "transition point from HG to LG (/ADC)"};
+ mutable Gaudi::Property<float> _ADCLimit{this, "ADCLimit", 3000, "ADC saturation of LG (/ADC)"};
+
+ // Output collections
+ DataHandle<edm4hep::CalorimeterHitCollection> w_DigiCaloCol{"DigiCaloCol", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::MCRecoCaloAssociationCollection> w_CaloAssociationCol{"HCALBarrelAssoCol", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::MCRecoCaloParticleAssociationCollection> w_MCPCaloAssociationCol{"HCALBarrelParticleAssoCol", Gaudi::DataHandle::Writer, this};
+};
+
+#endif
diff --git a/Digitisers/CaloDigi/src/HitStep.h b/Digitisers/CaloDigi/src/HitStep.h
new file mode 100644
index 00000000..2ae0715b
--- /dev/null
+++ b/Digitisers/CaloDigi/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 <x.T ;
+ }
+
+private:
+ double Q;
+ double T;
+
+};
+
+#endif
--
GitLab