diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
index fb8271bd1a920e4191f50f5cb5118a1818c7a1c1..b4aae77a91ff4c01f8e03155af1be994d8158a5d 100644
--- a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
+++ b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
@@ -171,6 +171,7 @@ digiMuon.SiPMEff = 1
digiMuon.EdepMin = 0.0 # no cut on GeV energy deposition
digiMuon.HitEdepMin = 6.0 # ADC counts
digiMuon.HitEdepMax = -1 # ADC counts, -1 means no upper cut
+digiMuon.TimeResolution = 2.0 # Digi hit time resolution, in unit ns
#digiMuon.OutputLevel = DEBUG
#########################################
diff --git a/Digitization/DigiMuon/src/MuonDigiAlg.cpp b/Digitization/DigiMuon/src/MuonDigiAlg.cpp
index f197975ddd28855363c087e0420cfcd745be3806..6867b8aecfbc0f27f0642472f381743ebab7483f 100644
--- a/Digitization/DigiMuon/src/MuonDigiAlg.cpp
+++ b/Digitization/DigiMuon/src/MuonDigiAlg.cpp
@@ -22,14 +22,12 @@
#include <DD4hep/Objects.h>
#include <DDRec/CellIDPositionConverter.h>
-#include "TVector3.h"
-#include <math.h>
+#include "TVector2.h"
#include <cmath>
#include <iostream>
#include <algorithm>
#include <map>
#include <random>
-#include <iostream>
#include <vector>
@@ -146,37 +144,36 @@ StatusCode MuonDigiAlg::execute()
debug() << "Collection " << m_inputMuonBarrel.fullKey() << " is unavailable in event " << m_nEvt << endmsg;
}
- int all_message[6];
- double Edep, xydist, hit_sipm_length, ADC, mcpz;
- edm4hep::Vector3d pos;
- dd4hep::Position ddpos;
+ int hit_data[6];
+ double hit_Edep, hit_xydist, hit_sipm_length, hit_mcpz;
+ double hit_time;
+ edm4hep::Vector3d hit_pos;
+ dd4hep::Position hit_ddpos;
unsigned long long cellid;
- std::array<unsigned long long, 2> key;
+ std::array<unsigned long long, 2> cell_key;
if(STHCol->size()>0)
{
for(auto simhit : *STHCol)
{
- GetSimHit(simhit, m_decoder_barrel, 2, all_message, Edep, cellid, xydist, pos, ddpos, key, mcpz);
- debug()<<"Barrel Simulation Hit::::::: "<<all_message[5]<<" "<<all_message[2]<<" "<<all_message[1]<<" "<<all_message[0]<<" "<<Edep<<endmsg;
- //if ( Edep > m_hit_Edep_max || Edep < m_hit_Edep_min ) continue; // Cut1: Edep of one hit
- //if ( xydist < 5162 && mcpz < 5475 || all_message[0] != 13 ) continue; // Cut2: For Muonid Efficiency
- SaveData_mapcell(key, Edep, all_message, pos);
+ GetSimHit(simhit, m_decoder_barrel, 2, hit_data, hit_Edep, cellid, hit_xydist, hit_pos, hit_ddpos, cell_key, hit_mcpz, hit_time);
+ debug()<<"Barrel Simulation Hit::::::: "<<hit_data[5]<<" "<<hit_data[2]<<" "<<hit_data[1]<<" "<<hit_data[0]<<" "<<hit_Edep<<endmsg;
+ SaveData_mapcell(cell_key, hit_Edep, hit_data, hit_pos, hit_time);
}
for (const auto& item : map_cell_edep)
{
- key = item.first;
- all_message[0] = map_cell_pdgid[key];
- all_message[1] = map_cell_layer[key];
- all_message[2] = map_cell_slayer[key];
- all_message[3] = map_cell_strip[key];
- all_message[4] = map_cell_env[key];
- all_message[5] = map_cell_fe[key];
- hit_sipm_length = Gethit_sipm_length_Barrel(m_cellIDConverter->position(key[0]), map_cell_pos[key], all_message);
- map_cell_adc[key] = EdeptoADC(item.second, hit_sipm_length);
- debug() <<"ADC::: "<<map_cell_adc[key]<<endmsg;
+ cell_key = item.first;
+ hit_data[0] = map_cell_pdgid[cell_key];
+ hit_data[1] = map_cell_layer[cell_key];
+ hit_data[2] = map_cell_slayer[cell_key];
+ hit_data[3] = map_cell_strip[cell_key];
+ hit_data[4] = map_cell_env[cell_key];
+ hit_data[5] = map_cell_fe[cell_key];
+ hit_sipm_length = Gethit_sipm_length_Barrel(m_cellIDConverter->position(cell_key[0]), map_cell_pos[cell_key], hit_data);
+ map_cell_adc[cell_key] = EdeptoADC(item.second, hit_sipm_length);
+ debug() <<"ADC::: "<<map_cell_adc[cell_key]<<endmsg;
}
-
+ // 2 for barrel
MuonHandle(2);
}
@@ -195,28 +192,24 @@ StatusCode MuonDigiAlg::execute()
if(STHCol->size()==0) return StatusCode::SUCCESS;
for(auto simhit : *STHCol)
{
- GetSimHit(simhit, m_decoder_endcap, 0, all_message, Edep, cellid, xydist, pos, ddpos, key, mcpz);
- //if ( Edep > m_hit_Edep_max || Edep < m_hit_Edep_min ) continue; // Cut1: Edep of one hit
- //if ( xydist < 5162 && mcpz < 5475 || all_message[0] != 13 ) continue; // Cut2: For Muonid Efficiency
- debug() <<"Endcap Simulation Hit::::::: "<<all_message[4]<<" "<<all_message[2]<<" "<<all_message[1]<<" "<<all_message[0]<<endmsg;
- // hit_sipm_length = Gethit_sipm_length_Endcap(ddpos, pos, all_message);
- // ADC = EdeptoADC(Edep, hit_sipm_length);
- SaveData_mapcell(key, Edep, all_message, pos);
+ GetSimHit(simhit, m_decoder_endcap, 0,hit_data, hit_Edep, cellid, hit_xydist, hit_pos, hit_ddpos, cell_key, hit_mcpz, hit_time);
+ debug() <<"Endcap Simulation Hit::::::: "<<hit_data[4]<<" "<<hit_data[2]<<" "<<hit_data[1]<<" "<<hit_data[0]<<endmsg;
+ SaveData_mapcell(cell_key, hit_Edep, hit_data, hit_pos, hit_time);
}
for (const auto& item : map_cell_edep)
{
- key = item.first;
- all_message[0] = map_cell_pdgid[key];
- all_message[1] = map_cell_layer[key];
- all_message[2] = map_cell_slayer[key];
- all_message[3] = map_cell_strip[key];
- all_message[4] = map_cell_env[key];
- all_message[5] = map_cell_fe[key];
- hit_sipm_length = Gethit_sipm_length_Endcap(m_cellIDConverter->position(key[0]), map_cell_pos[key], all_message);
- map_cell_adc[key] = EdeptoADC(item.second, hit_sipm_length);
- debug()<<"ADC::: "<<map_cell_adc[key]<<endmsg;
+ cell_key = item.first;
+ hit_data[0] = map_cell_pdgid[cell_key];
+ hit_data[1] = map_cell_layer[cell_key];
+ hit_data[2] = map_cell_slayer[cell_key];
+ hit_data[3] = map_cell_strip[cell_key];
+ hit_data[4] = map_cell_env[cell_key];
+ hit_data[5] = map_cell_fe[cell_key];
+ hit_sipm_length = Gethit_sipm_length_Endcap(m_cellIDConverter->position(cell_key[0]), map_cell_pos[cell_key], hit_data);
+ map_cell_adc[cell_key] = EdeptoADC(item.second, hit_sipm_length);
+ debug()<<"ADC::: "<<map_cell_adc[cell_key]<<endmsg;
}
-
+ // 0 for endcap
MuonHandle(0);
Clear();
m_nEvt++;
@@ -230,14 +223,18 @@ StatusCode MuonDigiAlg::finalize()
m_wfile->cd();
m_tree->Write();
m_wfile->Close();
- delete m_wfile, m_tree;
+ delete m_tree;
+ delete m_wfile;
}
info() << "Processed " << m_nEvt << " events " << endmsg;
return GaudiAlgorithm::finalize();
}
-void MuonDigiAlg::MuonHandle(int _i)
+
+void MuonDigiAlg::MuonHandle(int barrel_or_endcap)
{
+ //barrel_or_endcap: 2 is barrel, 0 is endcap
+
if(map_cell_edep.size()==0) return;
Cut3();
std::array<unsigned long long, 2> key1, key2;
@@ -253,24 +250,31 @@ void MuonDigiAlg::MuonHandle(int _i)
{
key2 = item2.first;
if(map_cell_edep[key2] <= m_EdepMin) continue;
- if(_i == 2)
+ if(barrel_or_endcap == 2) // barrel
{
- if(map_cell_slayer[key1]!=map_cell_slayer[key2] || map_cell_env[key1]!=map_cell_env[key2] || map_cell_fe[key1]!=map_cell_fe[key2] || map_cell_layer[key1]==map_cell_layer[key2] || map_cell_layer[key1]==2) continue;
+ if(map_cell_slayer[key1]!=map_cell_slayer[key2] ||
+ map_cell_env[key1]!=map_cell_env[key2] ||
+ map_cell_fe[key1]!=map_cell_fe[key2] ||
+ map_cell_layer[key1]==map_cell_layer[key2] ||
+ map_cell_layer[key1]==2) continue;
}
- if(_i == 0)
+ if(barrel_or_endcap == 0) // endcap
{
- if(map_cell_slayer[key1]!=map_cell_slayer[key2] || map_cell_fe[key1]!=map_cell_fe[key2] || map_cell_layer[key1]==map_cell_layer[key2] || map_cell_layer[key1]==2) continue;
+ if(map_cell_slayer[key1]!=map_cell_slayer[key2] ||
+ map_cell_fe[key1]!=map_cell_fe[key2] ||
+ map_cell_layer[key1]==map_cell_layer[key2] ||
+ map_cell_layer[key1]==2) continue;
if(std::abs(map_cell_env[key1] - map_cell_env[key2]) == 2) continue;
}
ddpos2 = m_cellIDConverter->position(key2[0]);
- double ddposi = (_i == 2) ? ddpos2.z() : ddpos2.x();
- Find_anotherlayer(_i, key1, key2, ddposi, anotherlayer_cell_num);
+ double ddposi = (barrel_or_endcap == 2) ? ddpos2.z() : ddpos2.x();
+ Find_anotherlayer(barrel_or_endcap, key1, key2, ddposi, anotherlayer_cell_num);
}
if(anotherlayer_cell_num == 0) continue;
for(int i=1; i<=anotherlayer_cell_num; i++)
{
true_pdgid = (std::abs(map_cell_pdgid[key1]) == 13 && std::abs(map_pdgid[i]) == 13) ? 1 : (map_pdgid[i] == -1) ? 2 : 0;
- edm4hep::Vector3d pos = (_i == 2) ? edm4hep::Vector3d(ddpos1.x() * 10, ddpos1.y() * 10, map_muonhit[i] * 10) : edm4hep::Vector3d(map_muonhit[i] * 10, ddpos1.y() * 10, ddpos1.z() * 10);
+ edm4hep::Vector3d pos = (barrel_or_endcap == 2) ? edm4hep::Vector3d(ddpos1.x() * 10, ddpos1.y() * 10, map_muonhit[i] * 10) : edm4hep::Vector3d(map_muonhit[i] * 10, ddpos1.y() * 10, ddpos1.z() * 10);
Save_trkhit(trkhitVec, key1, true_pdgid, pos);
}
}
@@ -289,91 +293,184 @@ void MuonDigiAlg::Clear()
map_cell_fe.clear();
map_cell_env.clear();
map_cell_pos.clear();
- map_muonhit.clear();
map_cell_pdgid.clear();
+ map_cell_time.clear();
+ map_muonhit.clear();
map_pdgid.clear();
}
-void MuonDigiAlg::GetSimHit(edm4hep::SimTrackerHit _simhit,
- dd4hep::DDSegmentation::BitFieldCoder* _m_decoder,
- int _i,
- int _message[6],
- double & _Edep,
- unsigned long long & _cellid,
- double & _xydist,
- edm4hep::Vector3d & _pos,
- dd4hep::Position & _ddpos,
- std::array<unsigned long long, 2> & _key,
- double & _mcpz)
+void MuonDigiAlg::GetSimHit(edm4hep::SimTrackerHit hit_simhit,
+ dd4hep::DDSegmentation::BitFieldCoder* hit_m_decoder,
+ int hit_i,
+ int hit_data[6],
+ double & hit_Edep,
+ unsigned long long & hit_cellid,
+ double & hit_xydist,
+ edm4hep::Vector3d & hit_pos,
+ dd4hep::Position & hit_ddpos,
+ std::array<unsigned long long, 2> & cell_key,
+ double & hit_mcpz,
+ double & hit_time)
{
- _cellid = _simhit.getCellID();
- auto mcp = _simhit.getMCParticle();
+ hit_cellid = hit_simhit.getCellID();
+ auto mcp = hit_simhit.getMCParticle();
auto mcppos = mcp.getEndpoint();
auto pdgid = mcp.getPDG();
- _message[0] = std::abs(mcp.getPDG()); //abspdgid
- _message[1] = _m_decoder->get(_cellid, "Layer");
- _message[2] = _m_decoder->get(_cellid, "Superlayer");
- _message[3] = _m_decoder->get(_cellid, "Stripe");
- _message[4] = _m_decoder->get(_cellid, "Env");
- _message[5] = (_i == 2) ? m_decoder_barrel->get(_cellid, "Fe") : m_decoder_endcap->get(_cellid, "Endcap");
- _Edep = _simhit.getEDep();//GeV
- _xydist = std::sqrt(mcppos[0]*mcppos[0]+mcppos[1]*mcppos[1]);
- _pos = _simhit.getPosition();
- _ddpos = m_cellIDConverter->position(_cellid);
- _key = std::array<unsigned long long, 2>{_cellid, _message[0]};
- _mcpz = std::abs(mcppos[2]);
+ hit_data[0] = std::abs(mcp.getPDG()); //abspdgid
+ hit_data[1] = hit_m_decoder->get(hit_cellid, "Layer");
+ hit_data[2] = hit_m_decoder->get(hit_cellid, "Superlayer");
+ hit_data[3] = hit_m_decoder->get(hit_cellid, "Stripe");
+ hit_data[4] = hit_m_decoder->get(hit_cellid, "Env");
+ // hit_i = 2 means barrel, hit_i !=2 means endcap
+ hit_data[5] = (hit_i == 2) ? hit_m_decoder->get(hit_cellid, "Fe") : hit_m_decoder->get(hit_cellid, "Endcap");
+ hit_Edep = (double)hit_simhit.getEDep();//GeV
+ hit_xydist = std::sqrt(mcppos[0]*mcppos[0]+mcppos[1]*mcppos[1]);
+ hit_pos = hit_simhit.getPosition();
+ hit_ddpos = m_cellIDConverter->position(hit_cellid);
+ cell_key = std::array<unsigned long long, 2>{hit_cellid, (unsigned long long)hit_data[0]};
+ hit_mcpz = std::abs(mcppos[2]);
+ hit_time = (double)hit_simhit.getTime(); // ns
}
-double MuonDigiAlg::Gethit_sipm_length_Barrel(dd4hep::Position _ddpos, edm4hep::Vector3d _pos, int _message[6])
+double MuonDigiAlg::Gethit_sipm_length_Barrel(dd4hep::Position hit_ddpos, edm4hep::Vector3d hit_pos, int hit_data[6])
{
- //calculate hit strip length
- double hit_strip_length = std::sqrt((_ddpos.x()-_pos[0]*0.1)*(_ddpos.x()-_pos[0]*0.1)+
- (_ddpos.y()-_pos[1]*0.1)*(_ddpos.y()-_pos[1]*0.1)+
- (_ddpos.z()-_pos[2]*0.1)*(_ddpos.z()-_pos[2]*0.1));
- if(_message[1] == 1)
+ // layer 1 is for strips parallel to the beam z-axis
+ if(hit_data[1] == 1)
{
- int tempnum = _message[4] + _message[2];
+ // We count the number of strips
+ // parpendicular to the beam axis,
+ // this number times the strip width 4 cm to get
+ // the length of the strips parallel to the beam axis
+ // hit_data[4] is Envolop, left or right
+ // hit_data[2] is superlayer
+ int tempnum = hit_data[4] + hit_data[2];
double strips_count = (tempnum % 2 == 0) ? 115.5 : 106.5; // 115 for long-half-barrel, 106 for short-half-barrel
- return (strips_count * 4.0) / 2 - hit_strip_length;
+ // strip length
+ double length = (strips_count * 4.0) ; // cm
+ // z axis positive half
+ if (hit_ddpos.z()>=0)
+ {
+
+ // sipm is at the z ~ 0 side
+ // (z<0) =================sipm|z=0|sipm=========c======x===== (z>0)
+ // ----->z direction center hit
+
+ return length/2.0 + (hit_pos[2] - hit_ddpos.z());
+ }
+ // z axis negative half
+ else
+ {
+ // sipm is at the z ~ 0 side
+ // (z<0) =========c====x====sipm|z=0|sipm==================== (z>0)
+ // center hit ----->z direction
+
+ return length/2.0 - (hit_pos[2] - hit_ddpos.z());
+ }
}
- if(_message[1] == 2)
+ // layer 2 is for strips perpendicular to the beam axis
+ else if(hit_data[1] == 2)
{
- // number of strips parallel to beam direction in each slayer, each strip width=4cm, so the number * 4 cm gives the length of strips in each slayer parpendicular to the beam axis
- return (strip_length[_message[2]-1] * 4)/2 - hit_strip_length;
+ // We count the number of strips
+ // parallel to beam direction in each slayer,
+ // this number times the strip width 4cm
+ // to give the length of strips in each slayer
+ // parpendicular to the beam axis
+ // strip length
+ double length = strip_length[hit_data[2]-1] * 4.0;
+
+ // define two vectors
+ TVector2 hit_pos_2vec(hit_pos[0]*0.1, hit_pos[1]*0.1);
+ TVector2 strip_pos_2vec(hit_ddpos.x(), hit_ddpos.y());
+ // define distance is hit - strip center
+ TVector2 distance_2vec = hit_pos_2vec - strip_pos_2vec;
+ // assume sipm is on the clockwise side
+ //
+ // ______sipm
+ // / \
+ // / \sipm
+ // | |
+ // | . | z>0 point to the out side of the screen
+ // | |sipm
+ // \ /
+ // \______/sipm
+
+ // define delta_phi as hit - strip_center
+ double dphi = TVector2::Phi_mpi_pi(hit_pos_2vec.Phi()-strip_pos_2vec.Phi());
+ //check delta_phi if it is positive,
+ // if yes, the hit is on the half of the strip near the sipm
+ // if no, the hit is on the half of the strip farther away from the sipm
+ if (dphi>=0.0)
+ {
+ return length/2.0 - distance_2vec.Mod();
+ }
+ else
+ {
+ return length/2.0 + distance_2vec.Mod();
+ }
}
+ // error
+ else
+ {
+ error() << "Error:: hit_data[1] can only be 1 for parallel to z axis or 2 for perpendicular to z-axis" << endmsg;
+ throw GaudiException("Error:: hit_data[1] can only be 1 for parallel to z axis or 2 for perpendicular to z-axis", "MuonDigiAlg", StatusCode::FAILURE);
+ }
}
-double MuonDigiAlg::Gethit_sipm_length_Endcap(dd4hep::Position _ddpos, edm4hep::Vector3d _pos, int _message[6])
+double MuonDigiAlg::Gethit_sipm_length_Endcap(dd4hep::Position hit_ddpos, edm4hep::Vector3d hit_pos, int hit_data[6])
{
- //calculate hit strip length
- double hit_strip_length = std::sqrt((_ddpos.x()-_pos[0]*0.1)*(_ddpos.x()-_pos[0]*0.1)+(_ddpos.y()-_pos[1]*0.1)*(_ddpos.y()-_pos[1]*0.1)+(_ddpos.z()-_pos[2]*0.1)*(_ddpos.z()-_pos[2]*0.1));
- _message[3]--; // strip id begin at 1
- return endcap_strip_length[_message[3]] * 0.5 * 100 - hit_strip_length;//cm
+ // strip length
+ double length = endcap_strip_length[hit_data[3]-1] * 100; // cm
+
+ // define two vectors
+ TVector2 hit_pos_2vec(hit_pos[0]*0.1, hit_pos[1]*0.1);
+ TVector2 strip_pos_2vec(hit_ddpos.x(), hit_ddpos.y());
+ // define distance is hit - strip center
+ TVector2 distance_2vec = hit_pos_2vec - strip_pos_2vec;
+ // assume sipm is on the clockwise side
+ // define delta_phi as hit - strip_center
+ double dphi = TVector2::Phi_mpi_pi(hit_pos_2vec.Phi()-strip_pos_2vec.Phi());
+ //check delta_phi if it is positive,
+ // if yes, the hit is on the half of the strip near the sipm
+ // if no, the hit is on the half of the strip farther away from the sipm
+ if (dphi>=0.0)
+ {
+ return length/2.0 - distance_2vec.Mod();
+ }
+ else
+ {
+ return length/2.0 + distance_2vec.Mod();
+ }
}
-double MuonDigiAlg::EdeptoADC(double _Edep, double _hit_sipm_length)
+double MuonDigiAlg::EdeptoADC(double hit_Edep, double hit_sipm_length)
{
// digitize to ADC
- double ADCmean = 34*_Edep*47.09/(23*0.00141);//mV
- if(_hit_sipm_length>10)
+ double ADCmean = 34*hit_Edep*47.09/(23*0.00141);//mV
+ if(hit_sipm_length>10)
{
- ADCmean = (16.0813*std::exp(-1*_hit_sipm_length/50.8147)+19.5474)*_Edep*47.09/(23*0.00141);//mV
+ ADCmean = (16.0813*std::exp(-1*hit_sipm_length/50.8147)+19.5474)*hit_Edep*47.09/(23*0.00141);//mV
}
return std::max(m_randSvc->generator(Rndm::Landau(ADCmean,7.922))->shoot(), 0.0);
}
-void MuonDigiAlg::SaveData_mapcell(std::array<unsigned long long, 2> _key, double _Edep, int _message[6], edm4hep::Vector3d _pos)
+void MuonDigiAlg::SaveData_mapcell(std::array<unsigned long long, 2> cell_key, double hit_Edep, int hit_data[6], edm4hep::Vector3d hit_pos, double hit_time)
{
+ // check if the key exists, if not initialize be zero
+ if ( map_cell_edep.find(cell_key)==map_cell_edep.end() )
+ {
+ map_cell_edep[cell_key] = 0.0;
+ }
// fill the cell loop
- map_cell_edep[_key] += _Edep;
- // map_cell_adc[_key] += _ADC;
- map_cell_layer[_key] = _message[1];
- map_cell_slayer[_key] = _message[2];
- map_cell_strip[_key] = _message[3];
- map_cell_fe[_key] = _message[5];
- map_cell_env[_key] = _message[4];
- map_cell_pos[_key] = _pos;
- map_cell_pdgid[_key] = _message[0];
+ map_cell_edep[cell_key] += hit_Edep;
+ map_cell_layer[cell_key] = hit_data[1];
+ map_cell_slayer[cell_key] = hit_data[2];
+ map_cell_strip[cell_key] = hit_data[3];
+ map_cell_fe[cell_key] = hit_data[5];
+ map_cell_env[cell_key] = hit_data[4];
+ map_cell_pos[cell_key] = hit_pos;
+ map_cell_pdgid[cell_key] = hit_data[0];
+ // smear the time:
+ map_cell_time[cell_key] = std::max(m_randSvc->generator(Rndm::Gauss(hit_time, m_time_resolution))->shoot(), 0.0);
}
void MuonDigiAlg::Cut3()
@@ -381,47 +478,53 @@ void MuonDigiAlg::Cut3()
//some criteria
for(const auto& item : map_cell_edep)
{
- if (map_cell_edep[item.first] <= m_EdepMin || m_randSvc->generator(Rndm::Flat(0, 1))->shoot() > m_hitEff || map_cell_adc[item.first] <= m_hit_Edep_min || (map_cell_adc[item.first] > m_hit_Edep_max && m_hit_Edep_max>0) )
+ if (map_cell_edep[item.first] <= m_EdepMin ||
+ m_randSvc->generator(Rndm::Flat(0, 1))->shoot() > m_hitEff ||
+ map_cell_adc[item.first] <= m_hit_Edep_min ||
+ (map_cell_adc[item.first] > m_hit_Edep_max && m_hit_Edep_max>0) )
{
map_cell_edep[item.first] = 0;
}
}
}
-void MuonDigiAlg::Find_anotherlayer(int _i, std::array<unsigned long long, 2> _key1, std::array<unsigned long long, 2> _key2, double _ddposi, int & _anotherlayer_cell_num)
+void MuonDigiAlg::Find_anotherlayer(int barrel_or_endcap, std::array<unsigned long long, 2> key1, std::array<unsigned long long, 2> key2, double ddposi, int & anotherlayer_cell_num)
{
- map_cell_strip[_key1] = -1;
- map_cell_strip[_key2] = -1;
- edm4hep::Vector3d _pos1 = map_cell_pos[_key1];
- edm4hep::Vector3d _pos2 = map_cell_pos[_key2];
- double Hit_min = _pos1[_i]*0.1-12;//10cm+2cm
- double Hit_max = _pos1[_i]*0.1+12;
- if(_ddposi<Hit_min || _ddposi>Hit_max) return;
- map_muonhit[++_anotherlayer_cell_num] = _ddposi;
- map_pdgid[_anotherlayer_cell_num] = -1;
- if(std::abs(_pos1[0]-_pos2[0])>40 && std::abs(_pos1[1]-_pos2[1])>40 && std::abs(_pos1[2]-_pos2[2])>40) return;
- map_pdgid[_anotherlayer_cell_num] = map_cell_pdgid[_key2];
+ map_cell_strip[key1] = -1;
+ map_cell_strip[key2] = -1;
+ edm4hep::Vector3d pos1 = map_cell_pos[key1];
+ edm4hep::Vector3d pos2 = map_cell_pos[key2];
+ // barrel_or_endcap: 2 is barrel, 0 is endcap
+ // barrel is [2] (z), endcap is [0] (x)
+ double Hit_min = pos1[barrel_or_endcap]*0.1-12;//10cm+2cm
+ double Hit_max = pos1[barrel_or_endcap]*0.1+12;
+ if(ddposi<Hit_min || ddposi>Hit_max) return;
+ map_muonhit[++anotherlayer_cell_num] = ddposi;
+ map_pdgid[anotherlayer_cell_num] = -1;
+ if(std::abs(pos1[0]-pos2[0])>40 && std::abs(pos1[1]-pos2[1])>40 && std::abs(pos1[2]-pos2[2])>40) return;
+ map_pdgid[anotherlayer_cell_num] = map_cell_pdgid[key2];
}
-void MuonDigiAlg::Save_trkhit(edm4hep::TrackerHitCollection* _trkhitVec, std::array<unsigned long long, 2> _key, int _pdgid, edm4hep::Vector3d _pos)
+void MuonDigiAlg::Save_trkhit(edm4hep::TrackerHitCollection* trkhitVec, std::array<unsigned long long, 2> key, int pdgid, edm4hep::Vector3d pos)
{
edm4hep::MutableTrackerHit trkHit;
- trkHit.setEDep(map_cell_adc[_key]);
- trkHit.setQuality(_pdgid);
- trkHit.setCellID(_key[0]);
- trkHit.setPosition(_pos);
- _trkhitVec->push_back( trkHit );
+ trkHit.setEDep((float)map_cell_adc[key]);
+ trkHit.setQuality(pdgid);
+ trkHit.setCellID(key[0]);
+ trkHit.setPosition(pos);
+ trkHit.setTime((float)map_cell_time[key]);
+ trkhitVec->push_back( trkHit );
}
-void MuonDigiAlg::Save_onelayer_signal(edm4hep::TrackerHitCollection* _trkhitVec)
+void MuonDigiAlg::Save_onelayer_signal(edm4hep::TrackerHitCollection* trkhitVec)
{
- for (const auto& _item : map_cell_strip)
+ for (const auto& item : map_cell_strip)
{
- if(_item.second == -1) continue;
- std::array<unsigned long long, 2> key = _item.first;
+ if(item.second == -1) continue;
+ std::array<unsigned long long, 2> key = item.first;
if ( map_cell_edep[key] <= m_EdepMin ) continue;
dd4hep::Position ddpos = m_cellIDConverter->position(key[0]);
edm4hep::Vector3d pos = edm4hep::Vector3d(ddpos.x()*10,ddpos.y()*10,ddpos.z()*10);//mm
- Save_trkhit(_trkhitVec, key, 3, pos); //just one layer signal in a superlayer
+ Save_trkhit(trkhitVec, key, 3, pos); //just one layer signal in a superlayer
}
}
diff --git a/Digitization/DigiMuon/src/MuonDigiAlg.h b/Digitization/DigiMuon/src/MuonDigiAlg.h
index 7b65b2c4b0c3f27769b2f655a2cba6ab856a0454..4e70dbb84ce0e7f1648d976887e8ae8ee458c5b6 100644
--- a/Digitization/DigiMuon/src/MuonDigiAlg.h
+++ b/Digitization/DigiMuon/src/MuonDigiAlg.h
@@ -49,17 +49,22 @@ class MuonDigiAlg : public GaudiAlgorithm
virtual StatusCode initialize() ;
virtual StatusCode execute() ;
virtual StatusCode finalize() ;
- void MuonHandle(int _i); // 2->Barrel 0->Endcap
+ void MuonHandle(int barrel_or_endcap); // 2->Barrel 0->Endcap
void Clear();
- void GetSimHit(edm4hep::SimTrackerHit _simhit, dd4hep::DDSegmentation::BitFieldCoder* _m_decoder, int _i, int _message[6], double & _Edep, unsigned long long & _cellid, double & _xydist, edm4hep::Vector3d & _pos, dd4hep::Position & _ddpos, std::array<unsigned long long, 2> & _key, double & _mcpz);
- double Gethit_sipm_length_Barrel(dd4hep::Position _ddpos, edm4hep::Vector3d _pos, int _message[6]);
- double Gethit_sipm_length_Endcap(dd4hep::Position _ddpos, edm4hep::Vector3d _pos, int _message[6]);
- double EdeptoADC(double _Edep, double _hit_sipm_length);
- void SaveData_mapcell(std::array<unsigned long long, 2> _key, double _Edep, int _message[6], edm4hep::Vector3d _pos);
+ void GetSimHit(edm4hep::SimTrackerHit hit_simhit, dd4hep::DDSegmentation::BitFieldCoder* hit_m_decoder, int hit_i,
+ int hit_data[6], double & hit_Edep, unsigned long long & hit_cellid, double & hit_xydist,
+ edm4hep::Vector3d & hit_pos, dd4hep::Position & hit_ddpos,
+ std::array<unsigned long long, 2> & cell_key, double & cell_mcpz, double & hit_time);
+ double Gethit_sipm_length_Barrel(dd4hep::Position hit_ddpos, edm4hep::Vector3d hit_pos, int hit_data[6]);
+ double Gethit_sipm_length_Endcap(dd4hep::Position hit_ddpos, edm4hep::Vector3d hit_pos, int hit_data[6]);
+ double EdeptoADC(double hit_Edep, double hit_sipm_length);
+ void SaveData_mapcell(std::array<unsigned long long, 2> cell_key, double hit_Edep, int hit_data[6], edm4hep::Vector3d hit_pos, double hit_time);
void Cut3();
- void Find_anotherlayer(int _i, std::array<unsigned long long, 2> _key1, std::array<unsigned long long, 2> _key2, double _ddposi, int & _anotherlayer_cell_num);
- void Save_trkhit(edm4hep::TrackerHitCollection* _trkhitVec, std::array<unsigned long long, 2> _key, int _pdgid, edm4hep::Vector3d _pos);
- void Save_onelayer_signal(edm4hep::TrackerHitCollection* _trkhitVec);
+ void Find_anotherlayer(int barrel_or_endcap, std::array<unsigned long long, 2> key1,
+ std::array<unsigned long long, 2> key2, double ddposi, int & anotherlayer_cell_num);
+ void Save_trkhit(edm4hep::TrackerHitCollection* trkhitVec, std::array<unsigned long long, 2> key, int pdgid, edm4hep::Vector3d pos);
+ void Save_onelayer_signal(edm4hep::TrackerHitCollection* trkhitVec);
+
@@ -76,6 +81,7 @@ class MuonDigiAlg : public GaudiAlgorithm
std::map<std::array<unsigned long long, 2>, double> map_cell_edep;
std::map<std::array<unsigned long long, 2>, double> map_cell_adc;
+ std::map<std::array<unsigned long long, 2>, double> map_cell_time;
std::map<std::array<unsigned long long, 2>, int> map_cell_layer;
std::map<std::array<unsigned long long, 2>, int> map_cell_slayer;
std::map<std::array<unsigned long long, 2>, int> map_cell_strip;
@@ -125,6 +131,7 @@ class MuonDigiAlg : public GaudiAlgorithm
Gaudi::Property<double> m_EdepMin{ this, "EdepMin", 0.0001, "Minimum Edep of a mip" };
Gaudi::Property<double> m_hit_Edep_min{ this, "HitEdepMin", 0.000001, "Minimum Edep of a mip"};
Gaudi::Property<double> m_hit_Edep_max{ this, "HitEdepMax", 0.1, "Maximum Edep of a mip"};
+ Gaudi::Property<double> m_time_resolution{ this, "TimeResolution", 2.0, "Time resolution of a hit, in nano-second ns." };
// number of strips parallel to beam direction
// in each slayer, each strip width=4cm, so