diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
index 15c8e1c942aa69ff3d61edc1b4f03a1461c179d0..fd872519d060136417b91bd7784696fdcb37466d 100644
--- a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
+++ b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
@@ -150,9 +150,10 @@ digiMuon.MuonEndcapTrackerHits = "MuonEndcapTrackerHits"
digiMuon.WriteNtuple = 0
digiMuon.OutFileName = "Digi_MUON.root"
digiMuon.SiPMEff = 1
-digiMuon.EdepMin = 0.0001
-digiMuon.HitEdepMin = 0.000001
-digiMuon.HitEdepMax = 0.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.OutputLevel = DEBUG
#########################################
################
diff --git a/Digitization/DigiMuon/src/MuonDigiAlg.cpp b/Digitization/DigiMuon/src/MuonDigiAlg.cpp
index 0523a7e5d58b0f094c85a4c040aa8ed6e15e6266..f197975ddd28855363c087e0420cfcd745be3806 100644
--- a/Digitization/DigiMuon/src/MuonDigiAlg.cpp
+++ b/Digitization/DigiMuon/src/MuonDigiAlg.cpp
@@ -43,12 +43,9 @@ MuonDigiAlg::MuonDigiAlg(const std::string& name, ISvcLocator* svcLoc)
// Input collections
declareProperty("MuonBarrelHitsCollection", m_inputMuonBarrel, "Handle of the Input SimTrackerHit collection");
declareProperty("MuonEndcapHitsCollection", m_inputMuonEndcap, "Handle of the Input SimTrackerHit collection");
-
-
// Output collections
declareProperty("MuonBarrelTrackerHits", m_outputMuonBarrel, "Handle of the output TrackerHit collection");
declareProperty("MuonEndcapTrackerHits", m_outputMuonEndcap, "Handle of the output TrackerHit collection");
- //declareProperty("MuonBarrelTrackerHitAssociationCollection", m_assMuonBarrel, "Handle of the Association collection between SimTrackerHit and TrackerHit");
}
@@ -90,7 +87,6 @@ StatusCode MuonDigiAlg::initialize()
m_tree->Branch("cell_pt", &cell_pt, "cell_pt/I");
}
-
m_geosvc = service<IGeomSvc>("GeomSvc");
if(!m_geosvc)
{
@@ -138,7 +134,6 @@ StatusCode MuonDigiAlg::initialize()
StatusCode MuonDigiAlg::execute()
{
const edm4hep::SimTrackerHitCollection* STHCol;
-//MuonBarrelRun:
Clear();
trkhitVec = m_outputMuonBarrel.createAndPut();
STHCol = nullptr;
@@ -150,55 +145,42 @@ StatusCode MuonDigiAlg::execute()
{
debug() << "Collection " << m_inputMuonBarrel.fullKey() << " is unavailable in event " << m_nEvt << endmsg;
}
- if(STHCol->size()==0) goto MuonEndcapRun;
-
- for(auto simhit : *STHCol)
- {
- GetSimHit(simhit, m_decoder_barrel);
- Fe = m_decoder_barrel->get(cellid, "Fe");
- if ( Edep > m_hit_Edep_max || Edep < m_hit_Edep_min ) continue; // Cut1: Edep of one hit
- if ( xydist < 5162 && std::abs(mcppos[2]) < 5475 || abspdgid != 13 ) continue; // Cut2: The particle decays in the Muon Detector or is not a muon
- std::cout<<"Barrel Simulation Hit::::::: "<<Fe<<" "<<slayer<<" "<<layer<<" "<<mcppos[0]<<" "<<mcppos[1]<<" "<<mcppos[2]<<" "<<abspdgid<<std::endl;
- hit_sipm_length = Gethit_sipm_length_Barrel();
- EdeptoADC();
- SaveData_mapcell();
- }
- if(map_cell_edep.size()==0) goto MuonEndcapRun;
- Cut3();
+ int all_message[6];
+ double Edep, xydist, hit_sipm_length, ADC, mcpz;
+ edm4hep::Vector3d pos;
+ dd4hep::Position ddpos;
+ unsigned long long cellid;
+ std::array<unsigned long long, 2> key;
- for (const auto& item1 : map_cell_edep)
- {
- key1 = item1.first;
- if(!Mapcell_todata(all_message1, key1)) continue;
- Save_pos(1);
- anotherlayer_cell_num = 0;
- for (const auto& item2 : map_cell_edep)
+ if(STHCol->size()>0)
+ {
+ for(auto simhit : *STHCol)
{
- key2 = item2.first;
- if(!Mapcell_todata(all_message2, key2)) continue;
- if(all_message1[1]!=all_message2[1] || all_message1[2]!=all_message2[2] || all_message1[3]!=all_message2[3] || all_message1[0]==all_message2[0] || all_message1[0]==2) continue;
- Save_pos(2);
- std::cout<<":::::: "<<cellid1<<" "<<key1[0]<<std::endl;
- std::cout<<ddpos1.x()<<" "<<ddpos1.y()<<std::endl;
- std::cout<<pos1[0]<<" "<<pos1[1]<<std::endl;
- Renew_strip(key1, key2);
- Find_anotherlayer(2, pos1, pos2, ddpos2.z(), all_message2[5]);
+ 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);
}
-
- if(anotherlayer_cell_num == 0) continue;
- for(int i=1; i<=anotherlayer_cell_num; i++)
+ for (const auto& item : map_cell_edep)
{
- true_pdgid = 0;
- true_pdgid = Get_true_pdgid(i, all_message1[5]);
- pos = edm4hep::Vector3d(ddpos1.x()*10,ddpos1.y()*10,map_muonhit[i]*10);//mm
-std::cout<<pos[0]<<" "<<pos[1]<<" "<<pos[2]<<std::endl;
- Save_trkhit(trkhitVec, key1, true_pdgid, cellid1, pos);
+ 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;
}
+
+ MuonHandle(2);
}
- Save_onelayer_signal(trkhitVec);
-MuonEndcapRun:
+ // Endcaps:
Clear();
trkhitVec = m_outputMuonEndcap.createAndPut();
STHCol = nullptr;
@@ -211,56 +193,38 @@ MuonEndcapRun:
debug() << "Collection " << m_inputMuonEndcap.fullKey() << " is unavailable in event " << m_nEvt << endmsg;
}
if(STHCol->size()==0) return StatusCode::SUCCESS;
-
for(auto simhit : *STHCol)
{
- GetSimHit(simhit, m_decoder_endcap);
- Fe = m_decoder_endcap->get(cellid, "Endcap");
- if ( Edep > m_hit_Edep_max || Edep < m_hit_Edep_min ) continue; // Cut1: Edep of one hit
- if ( xydist < 5162 && std::abs(mcppos[2]) < 5475 || abspdgid != 13 ) continue; // Cut2: The particle decays in the Muon Detector or is not a muon
- std::cout<<"Endcap Simulation Hit::::::: "<<Env<<" "<<slayer<<" "<<layer<<" "<<mcppos[0]<<" "<<mcppos[1]<<" "<<mcppos[2]<<" "<<abspdgid<<std::endl;
- hit_sipm_length = Gethit_sipm_length_Endcap();
- EdeptoADC();
- SaveData_mapcell();
- }
-
- if(map_cell_edep.size()==0) return StatusCode::SUCCESS;
- Cut3();
-
- for (const auto& item1 : map_cell_edep)
+ 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);
+ }
+ for (const auto& item : map_cell_edep)
{
- key1 = item1.first;
- if(!Mapcell_todata(all_message1, key1)) continue;
- Save_pos(1);
- anotherlayer_cell_num = 0;
- for (const auto& item2 : map_cell_edep)
- {
- key2 = item2.first;
- if(!Mapcell_todata(all_message2, key2)) continue;
- if(all_message1[1]!=all_message2[1] || all_message1[3]!=all_message2[3] || all_message1[0]==all_message2[0] || all_message1[0]==2) continue;
- if(std::abs(all_message1[2] - all_message2[2]) == 2) continue;
- Save_pos(2);
- Renew_strip(key1, key2);
- Find_anotherlayer(0, pos1, pos2, ddpos2.x(), all_message2[5]);
- }
- if(anotherlayer_cell_num == 0) continue;
- for(int i=1; i<=anotherlayer_cell_num; i++)
- {
- true_pdgid = 0;
- true_pdgid = Get_true_pdgid(i, all_message1[5]);
- pos = edm4hep::Vector3d(map_muonhit[i]*10,ddpos1.y()*10,ddpos1.z()*10);//mm
- Save_trkhit(trkhitVec, key1, true_pdgid, cellid1, pos);
- }
+ 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;
}
- Save_onelayer_signal(trkhitVec);
+ MuonHandle(0);
+ Clear();
m_nEvt++;
return StatusCode::SUCCESS;
}
StatusCode MuonDigiAlg::finalize()
{
-
if(_writeNtuple)
{
m_wfile->cd();
@@ -272,6 +236,47 @@ StatusCode MuonDigiAlg::finalize()
return GaudiAlgorithm::finalize();
}
+void MuonDigiAlg::MuonHandle(int _i)
+{
+ if(map_cell_edep.size()==0) return;
+ Cut3();
+ std::array<unsigned long long, 2> key1, key2;
+ dd4hep::Position ddpos1, ddpos2;
+ int anotherlayer_cell_num, true_pdgid;
+ for (const auto& item1 : map_cell_edep)
+ {
+ key1 = item1.first;
+ if(map_cell_edep[key1] <= m_EdepMin) continue;
+ ddpos1 = m_cellIDConverter->position(key1[0]);
+ anotherlayer_cell_num = 0;
+ for (const auto& item2 : map_cell_edep)
+ {
+ key2 = item2.first;
+ if(map_cell_edep[key2] <= m_EdepMin) continue;
+ if(_i == 2)
+ {
+ 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(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);
+ }
+ 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);
+ Save_trkhit(trkhitVec, key1, true_pdgid, pos);
+ }
+ }
+ Save_onelayer_signal(trkhitVec);
+}
+
void MuonDigiAlg::Clear()
{
n_hit = 0;
@@ -289,104 +294,86 @@ void MuonDigiAlg::Clear()
map_pdgid.clear();
}
-void MuonDigiAlg::GetSimHit(edm4hep::SimTrackerHit _simhit, dd4hep::DDSegmentation::BitFieldCoder* _m_decoder)
+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)
{
+ _cellid = _simhit.getCellID();
auto mcp = _simhit.getMCParticle();
- cellid = _simhit.getCellID();
- mcppos = mcp.getEndpoint();
- pdgid = mcp.getPDG();
- abspdgid = std::abs(mcp.getPDG());
- key = std::array<unsigned long long, 2>{cellid, abspdgid};
- Edep = _simhit.getEDep();//GeV
- layer = _m_decoder->get(cellid, "Layer");
- slayer = _m_decoder->get(cellid, "Superlayer");
- strip = _m_decoder->get(cellid, "Stripe");
- Env = _m_decoder->get(cellid, "Env");
- pos = _simhit.getPosition();
- ddpos = m_cellIDConverter->position(cellid);
- xydist = std::sqrt(mcppos[0]*mcppos[0]+mcppos[1]*mcppos[1]);
+ 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]);
}
-double MuonDigiAlg::Gethit_sipm_length_Barrel()
+double MuonDigiAlg::Gethit_sipm_length_Barrel(dd4hep::Position _ddpos, edm4hep::Vector3d _pos, int _message[6])
{
//calculate hit strip length
- 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(layer == 1)
+ 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)
{
- int tempnum = Env + slayer;
- if(tempnum%2 == 0)
- {
- // 115 is the number of strips perpendicular
- // to the beam axis in long-half-barrel
- return (115.5 * 4.0)/2 - hit_strip_length;
- }
- else
- {
- // 106 is the number of strips perpendicular
- // to the beam axis in short-half-barrel
- return (106.5 * 4.0)/2 - hit_strip_length;
- }
+ int tempnum = _message[4] + _message[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;
}
- if(layer == 2)
+ if(_message[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[slayer-1] * 4)/2 - hit_strip_length;
+ // 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;
}
}
-double MuonDigiAlg::Gethit_sipm_length_Endcap()
+double MuonDigiAlg::Gethit_sipm_length_Endcap(dd4hep::Position _ddpos, edm4hep::Vector3d _pos, int _message[6])
{
//calculate hit strip length
- 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));
- strip--; // strip id begin at 1
- return endcap_strip_length[strip] * 0.5 * 100 - hit_strip_length;//cm
+ 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
}
-void MuonDigiAlg::EdeptoADC()
+double MuonDigiAlg::EdeptoADC(double _Edep, double _hit_sipm_length)
{
// digitize to ADC
- ADCmean = 34*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
- }
- ADC = m_randSvc->generator(Rndm::Landau(ADCmean,7.922))->shoot();
- if(ADC<0)
+ double ADCmean = 34*_Edep*47.09/(23*0.00141);//mV
+ if(_hit_sipm_length>10)
{
- ADC = 0;
+ ADCmean = (16.0813*std::exp(-1*_hit_sipm_length/50.8147)+19.5474)*_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()
+void MuonDigiAlg::SaveData_mapcell(std::array<unsigned long long, 2> _key, double _Edep, int _message[6], edm4hep::Vector3d _pos)
{
// fill the cell loop
- map_cell_edep[key] += Edep;
- map_cell_adc[key] += ADC;
- map_cell_layer[key] = layer;
- map_cell_slayer[key] = slayer;
- map_cell_strip[key] = strip;
- map_cell_fe[key] = Fe;
- map_cell_env[key] = Env;
- map_cell_pos[key] = pos;
- map_cell_pdgid[key] = pdgid;
-}
-
-bool MuonDigiAlg::Mapcell_todata(int _message[6], std::array<unsigned long long, 2> _key)
-{
- Edep = map_cell_edep[_key];
- if ( Edep < m_EdepMin ) return false;
- _message[0] = map_cell_layer[_key];
- _message[1] = map_cell_slayer[_key];
- _message[2] = map_cell_env[_key];
- _message[3] = map_cell_fe[_key];
- _message[4] = map_cell_strip[_key];
- _message[5] = map_cell_pdgid[_key];
- return true;
+ 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];
}
void MuonDigiAlg::Cut3()
@@ -394,46 +381,36 @@ void MuonDigiAlg::Cut3()
//some criteria
for(const auto& item : map_cell_edep)
{
- key = item.first;
- if (map_cell_edep[key] < m_EdepMin || m_randSvc->generator(Rndm::Flat(0, 1))->shoot() > m_hitEff || map_cell_adc[key] == 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[key] = 0;
+ map_cell_edep[item.first] = 0;
}
}
}
-void MuonDigiAlg::Save_trkhit(edm4hep::TrackerHitCollection* _trkhitVec, std::array<unsigned long long, 2> _key, int _pdgid, unsigned long long _cellid, edm4hep::Vector3d _pos)
-{
- edm4hep::MutableTrackerHit trkHit;
- trkHit.setEDep(map_cell_adc[_key]);
- trkHit.setQuality(_pdgid);
- trkHit.setCellID(_cellid);
- trkHit.setPosition(_pos);
- _trkhitVec->push_back( trkHit );
-}
-
-void MuonDigiAlg::Renew_strip(std::array<unsigned long long, 2> _key1, std::array<unsigned long long, 2> _key2)
+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)
{
map_cell_strip[_key1] = -1;
map_cell_strip[_key2] = -1;
-}
-
-void MuonDigiAlg::Find_anotherlayer(int _i, edm4hep::Vector3d _pos1, edm4hep::Vector3d _pos2, double _ddposi, int _pdgid)
-{
- Hit_min = _pos1[_i]*0.1-12;//10cm+2cm
- Hit_max = _pos1[_i]*0.1+12;
+ 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;
- anotherlayer_cell_num++;
- map_muonhit[anotherlayer_cell_num] = _ddposi;
- map_pdgid[anotherlayer_cell_num] = -1;
+ 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] = _pdgid;
+ map_pdgid[_anotherlayer_cell_num] = map_cell_pdgid[_key2];
}
-int MuonDigiAlg::Get_true_pdgid(int _i, int _pdgid)
+void MuonDigiAlg::Save_trkhit(edm4hep::TrackerHitCollection* _trkhitVec, std::array<unsigned long long, 2> _key, int _pdgid, edm4hep::Vector3d _pos)
{
- if(std::abs(_pdgid) == 13 && std::abs(map_pdgid[_i]) == 13) { return 1; } //pdgid->true muon hit
- if(map_pdgid[_i]==-1) { return 2; } //ghost point
+ edm4hep::MutableTrackerHit trkHit;
+ trkHit.setEDep(map_cell_adc[_key]);
+ trkHit.setQuality(_pdgid);
+ trkHit.setCellID(_key[0]);
+ trkHit.setPosition(_pos);
+ _trkhitVec->push_back( trkHit );
}
void MuonDigiAlg::Save_onelayer_signal(edm4hep::TrackerHitCollection* _trkhitVec)
@@ -441,28 +418,10 @@ void MuonDigiAlg::Save_onelayer_signal(edm4hep::TrackerHitCollection* _trkhitVec
for (const auto& _item : map_cell_strip)
{
if(_item.second == -1) continue;
- key = _item.first;
- cellid = key[0];
- Edep = map_cell_edep[key];
- if ( Edep < m_EdepMin ) continue;
- ddpos = m_cellIDConverter->position(cellid);
- pos = edm4hep::Vector3d(ddpos.x()*10,ddpos.y()*10,ddpos.z()*10);//mm
- Save_trkhit(_trkhitVec, key, 3, cellid1, pos); //just one layer signal in a superlayer
+ 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
}
}
-
-void MuonDigiAlg::Save_pos(int _i)
-{
- if(_i == 1)
- {
- cellid1 = key1[0];
- ddpos1 = m_cellIDConverter->position(cellid1);
- pos1 = map_cell_pos[key1];
- }
- if(_i == 2)
- {
- cellid2 = key2[0];
- ddpos2 = m_cellIDConverter->position(cellid2);
- pos2 = map_cell_pos[key2];
- }
-}
diff --git a/Digitization/DigiMuon/src/MuonDigiAlg.h b/Digitization/DigiMuon/src/MuonDigiAlg.h
index 78a4e066f5fafe41db625ea65a91d553f197d37b..7b65b2c4b0c3f27769b2f655a2cba6ab856a0454 100644
--- a/Digitization/DigiMuon/src/MuonDigiAlg.h
+++ b/Digitization/DigiMuon/src/MuonDigiAlg.h
@@ -49,26 +49,21 @@ class MuonDigiAlg : public GaudiAlgorithm
virtual StatusCode initialize() ;
virtual StatusCode execute() ;
virtual StatusCode finalize() ;
+ void MuonHandle(int _i); // 2->Barrel 0->Endcap
void Clear();
- void GetSimHit(edm4hep::SimTrackerHit _simhit, dd4hep::DDSegmentation::BitFieldCoder* _m_decoder);
- double Gethit_sipm_length_Barrel();
- double Gethit_sipm_length_Endcap();
- void EdeptoADC();
- void SaveData_mapcell();
- bool Mapcell_todata(int _message[6], std::array<unsigned long long, 2> _key);
+ 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 Cut3();
- void Save_trkhit(edm4hep::TrackerHitCollection* _trkhitVec, std::array<unsigned long long, 2> _key, int _pdgid, unsigned long long _cellid, edm4hep::Vector3d _pos);
- void Renew_strip(std::array<unsigned long long, 2> _key1, std::array<unsigned long long, 2> _key2);
- void Find_anotherlayer(int _i, edm4hep::Vector3d _pos1, edm4hep::Vector3d _pos2, double _ddposi, int _pdgid);
- int Get_true_pdgid(int _i, int _pdgid);
+ 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 Save_pos(int _i);
-
-
- protected:
+ protected:
TTree* m_tree;
TFile* m_wfile;
@@ -134,14 +129,6 @@ class MuonDigiAlg : public GaudiAlgorithm
// number of strips parallel to beam direction
// in each slayer, each strip width=4cm, so
int strip_length[6] = {26, 38, 50, 62, 74, 86}; //for barrel
- std::array<unsigned long long, 2> key, key1, key2;
- int pdgid, layer, slayer, strip, Fe, Env, anotherlayer_cell_num, true_pdgid;
- unsigned long long cellid, abspdgid, cellid1, cellid2;
- double Edep, xydist, hit_sipm_length, ADC, hit_strip_length, ADCmean, Hit_max, Hit_min;
- edm4hep::Vector3d pos, pos1, pos2;
- edm4hep::Vector3d mcppos;
- dd4hep::Position ddpos, ddpos1, ddpos2;
- int all_message1[6], all_message2[6]; // int layer1, slayer1, strip1, Fe1, Env1, pdgid1;
double endcap_strip_length[193] = {2.12, 2.12, 2.12, 2.13, 2.14, 2.15, 2.16, 2.17, 2.19, 2.22,
2.24, 2.28, 2.32, 2.37, 2.45, 2.65, 2.64, 2.63, 2.62, 2.61,
2.60, 2.59, 2.57, 2.56, 2.54, 2.53, 2.51, 2.50, 2.48, 2.46,