diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
index 8c009a18e85f151bd272b67bdd1541badbbdc3c6..3f91a8f40c9e28a0a23b29ceec83d4c4c4da4272 100644
--- a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
+++ b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
@@ -54,7 +54,8 @@ podioinput = PodioInput("PodioReader", collections=[
# "SETCollection",
"OTKBarrelCollection",
"FTDCollection",
- "MuonBarrelCollection"
+ "MuonBarrelCollection",
+ "MuonEndcapCollection"
])
@@ -156,7 +157,9 @@ digiTPC.TPCTrackerHitsCol = gashitname
from Configurables import MuonDigiAlg
digiMuon = MuonDigiAlg("MuonDigiAlg")
digiMuon.MuonBarrelHitsCollection = "MuonBarrelCollection"
+digiMuon.MuonEndcapHitsCollection = "MuonEndcapCollection"
digiMuon.MuonBarrelTrackerHits = "MuonBarrelTrackerHits"
+digiMuon.MuonEndcapTrackerHits = "MuonEndcapTrackerHits"
digiMuon.WriteNtuple = 0
digiMuon.OutFileName = "Digi_MUON.root"
#########################################
diff --git a/Digitization/DigiMuon/src/MuonDigiAlg.cpp b/Digitization/DigiMuon/src/MuonDigiAlg.cpp
index 800b07625c80f9be59c96e1f94e309b7d506c3dc..347cb6c8c82d48cfb99437aec490509d1e36ad7d 100644
--- a/Digitization/DigiMuon/src/MuonDigiAlg.cpp
+++ b/Digitization/DigiMuon/src/MuonDigiAlg.cpp
@@ -42,9 +42,12 @@ 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");
}
@@ -98,17 +101,26 @@ StatusCode MuonDigiAlg::initialize()
error() << "Failed to get the Detector from GeomSvc" << endmsg;
return StatusCode::FAILURE;
}
- std::string readout_name = m_inputMuonBarrel.objKey();
- debug() << "Readout name: " << readout_name << endmsg;
m_cellIDConverter = new dd4hep::rec::CellIDPositionConverter(*m_dd4hep);
-
if(!m_cellIDConverter)
{
error() << "Failed to get the m_cellIDConverter." << endmsg;
return StatusCode::FAILURE;
}
- m_decoder = m_geosvc->getDecoder(readout_name);
- if(!m_decoder)
+
+ std::string readout_name1 = m_inputMuonBarrel.objKey();
+ debug() << "Readout name: " << readout_name1 << endmsg;
+ m_decoder1 = m_geosvc->getDecoder(readout_name1);
+ if(!m_decoder1)
+ {
+ error() << "Failed to get the decoder. " << endmsg;
+ return StatusCode::FAILURE;
+ }
+
+ std::string readout_name2 = m_inputMuonEndcap.objKey();
+ debug() << "Readout name: " << readout_name2 << endmsg;
+ m_decoder2 = m_geosvc->getDecoder(readout_name2);
+ if(!m_decoder2)
{
error() << "Failed to get the decoder. " << endmsg;
return StatusCode::FAILURE;
@@ -123,175 +135,452 @@ StatusCode MuonDigiAlg::initialize()
StatusCode MuonDigiAlg::execute()
{
-
-
- Clear();
- trkhitVec = m_outputMuonBarrel.createAndPut();
- //auto assVec = m_assMuonBarrel.createAndPut();
-
- const edm4hep::SimTrackerHitCollection* STHCol = nullptr;
- try
+ for(int muonmode = 0; muonmode<2; muonmode++)
{
- STHCol = m_inputMuonBarrel.get();
- }
- catch(GaudiException &e)
- {
- debug() << "Collection " << m_inputMuonBarrel.fullKey() << " is unavailable in event " << m_nEvt << endmsg;
- return StatusCode::SUCCESS;
- }
- if(STHCol->size()==0) return StatusCode::SUCCESS;
- debug() << m_inputMuonBarrel.fullKey() << " has SimTrackerHit "<< STHCol->size() << endmsg;
-
-
- // 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};
-
- // define variables to be repeatedly used later
- unsigned long long cellid;
- double Edep, ADC, ADCmean;
- int layer, slayer, strip, Fe, Env;
- edm4hep::Vector3d pos;
- dd4hep::Position ddpos;
- // loop over all hits
- for(auto simhit : *STHCol)
- {
- cellid = simhit.getCellID();
- Edep = simhit.getEDep();//GeV
- layer = m_decoder->get(cellid, "Layer");
- slayer = m_decoder->get(cellid, "Superlayer");
- strip = m_decoder->get(cellid, "Stripe");
- Fe = m_decoder->get(cellid, "Fe");
- Env = m_decoder->get(cellid, "Env");
- pos = simhit.getPosition();
-
- ddpos = m_cellIDConverter->position(cellid);
- debug() << "Position:: " << ddpos.x() << " " << ddpos.y() << " " << ddpos.z() << endmsg;
-
- // if not satisfy energy requirement, skip this hit
- //if ( Edep>0.01 || Edep<0.0001 ) continue;
-
- //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));
- double hit_sipm_length;
-
- if(layer == 1)
+ if(muonmode == 0)//muon barrel run
{
- int tempnum = Env + slayer;
- if(tempnum%2 == 0)
+ Clear();
+ if(m_MuonMode == 2) continue;
+ trkhitVec = m_outputMuonBarrel.createAndPut();
+ //auto assVec = m_assMuonBarrel.createAndPut();
+
+ const edm4hep::SimTrackerHitCollection* STHCol = nullptr;
+ try
{
- // 115 is the number of strips parpendicular
- // to the beam axis in long-half-barrel
- hit_sipm_length = (115 * 4.0)/2 - hit_strip_length;
+ STHCol = m_inputMuonBarrel.get();
}
- else
+ catch(GaudiException &e)
{
- // 106 is the number of strips parpendicular
- // to the beam axis in short-half-barrel
- hit_sipm_length = (106 * 4.0)/2 - hit_strip_length;
+ debug() << "Collection " << m_inputMuonBarrel.fullKey() << " is unavailable in event " << m_nEvt << endmsg;
+ continue;
}
- }
- if(layer == 2)
- {
+ if(STHCol->size()==0) continue;
+ debug() << m_inputMuonBarrel.fullKey() << " has SimTrackerHit "<< STHCol->size() << endmsg;
+
// 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
- hit_sipm_length = (strip_length[slayer-1] * 4)/2 - hit_strip_length;
- }
+ int strip_length[6] = {26, 38, 50, 62, 74, 86};
+
+ // define variables to be repeatedly used later
+ std::array<unsigned long long, 2> key;
+ int pdgid;
+ unsigned long long cellid, abspdgid;
+ double Edep, ADC, ADCmean;
+ int layer, slayer, strip, Fe, Env;
+ edm4hep::Vector3d pos;
+ dd4hep::Position ddpos;
+ // loop over all hits
+ for(auto simhit : *STHCol)
+ {
+ auto mcp = simhit.getMCParticle();
+ cellid = simhit.getCellID();
+ pdgid = mcp.getPDG();
+ abspdgid = std::abs(mcp.getPDG());
+ key = std::array<unsigned long long, 2>{cellid, abspdgid};
+ Edep = simhit.getEDep();//GeV
+ layer = m_decoder1->get(cellid, "Layer");
+ slayer = m_decoder1->get(cellid, "Superlayer");
+ strip = m_decoder1->get(cellid, "Stripe");
+ Fe = m_decoder1->get(cellid, "Fe");
+ Env = m_decoder1->get(cellid, "Env");
+ pos = simhit.getPosition();
+ ddpos = m_cellIDConverter->position(cellid);
+ debug() << "Position:: " << ddpos.x() << " " << ddpos.y() << " " << ddpos.z() << endmsg;
+
+ // if not satisfy energy requirement, skip this hit
+ if ( Edep>0.1 || Edep<0.000001 ) continue;
+
+ //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));
+ double hit_sipm_length;
+
+ if(layer == 1)
+ {
+ int tempnum = Env + slayer;
+ if(tempnum%2 == 0)
+ {
+ // 115 is the number of strips parpendicular
+ // to the beam axis in long-half-barrel
+ hit_sipm_length = (115 * 4.0)/2 - hit_strip_length;
+ }
+ else
+ {
+ // 106 is the number of strips parpendicular
+ // to the beam axis in short-half-barrel
+ hit_sipm_length = (106 * 4.0)/2 - hit_strip_length;
+ }
+ }
+ if(layer == 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
+ hit_sipm_length = (strip_length[slayer-1] * 4)/2 - hit_strip_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 = rand_muon.Landau(ADCmean,7.922);
+ if(ADC<0)
+ {
+ ADC = 0;
+ }
+
+ // 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;
+
+ // write to tree
+ if(_writeNtuple)
+ {
+ hit_cellid[n_hit] = cellid;
+ hit_posx[n_hit] = pos[0];
+ hit_posy[n_hit] = pos[1];
+ hit_posz[n_hit] = pos[2];
+ hit_edep[n_hit] = Edep;
+ hit_layer[n_hit] = layer;
+ hit_slayer[n_hit] = slayer;
+ hit_strip[n_hit] = strip;
+ hit_fe[n_hit] = Fe;
+ hit_env[n_hit] = Env;
+ n_hit++;
+ }
+
+ }
+ std::array<unsigned long long, 2> key1;
+ unsigned long long cellid1;
+ int layer1, slayer1, strip1, Fe1, Env1, pdgid1;
+ edm4hep::Vector3d pos1;
+ dd4hep::Position ddpos1;
+
+ std::array<unsigned long long, 2> key2;
+ unsigned long long cellid2;
+ int layer2, slayer2, strip2, Fe2, Env2, pdgid2;
+ edm4hep::Vector3d pos2;
+ dd4hep::Position ddpos2;
+
+ if(map_cell_edep.size()==0) continue;
+ //some criteria
+ for(const auto& item : map_cell_edep)
+ {
+ key = item.first;
+ if (map_cell_edep[key] < m_EdepMin || rand_muon.Uniform(0, 1)>m_hitEff || map_cell_adc[key] == 0)
+ {
+ map_cell_edep[key] = 0;
+ }
+ }
- // 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
+ double muonhitx, muonhity, muonhitz;
+ double Hit_max, Hit_min;
+
+ // loop over all cells
+ for (const auto& item1 : map_cell_edep)
+ {
+ key1 = item1.first;
+ cellid1 = key1[0];
+ Edep = map_cell_edep[key1];
+ if ( Edep < 0.0001 ) continue;
+ int anotherlayer_cell_num = 0;
+
+ ADC = map_cell_adc[key1];
+ layer1 = map_cell_layer[key1];
+ slayer1 = map_cell_slayer[key1];
+ strip1 = map_cell_strip[key1];
+ Fe1 = map_cell_fe[key1];
+ Env1 = map_cell_env[key1];
+ pos1 = map_cell_pos[key1];
+ pdgid1 = map_cell_pdgid[key1];
+ ddpos1 = m_cellIDConverter->position(cellid1);
+
+ for (const auto& item2 : map_cell_edep)
+ {
+ key2 = item2.first;
+ cellid2 = key2[0];
+ Edep = map_cell_edep[key2];
+ if ( Edep < 0.0001 ) continue;
+
+ layer2 = map_cell_layer[key2];
+ slayer2 = map_cell_slayer[key2];
+ strip2 = map_cell_strip[key2];
+ Fe2 = map_cell_fe[key2];
+ Env2 = map_cell_env[key2];
+ pos2 = map_cell_pos[key2];
+ pdgid2 = map_cell_pdgid[key2];
+ ddpos2 = m_cellIDConverter->position(cellid2);
+
+ if(slayer1!=slayer2 || Fe1!=Fe2 || Env1!=Env2 || layer1==layer2 || layer1==2) continue;//layer1 = 1
+ Hit_min = pos1[2]*0.1-12;
+ Hit_max = pos1[2]*0.1+12;
+ if(ddpos2.z()<Hit_min || ddpos2.z()>Hit_max) continue;
+ anotherlayer_cell_num++;
+ map_muonhit[anotherlayer_cell_num] = ddpos2.z();
+ 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) continue;
+ map_pdgid[anotherlayer_cell_num] = pdgid2;
+ }
+
+ if(anotherlayer_cell_num == 0) continue;
+ if(anotherlayer_cell_num > 0)
+ {
+ for(int i=1; i<=anotherlayer_cell_num; i++)
+ {
+ muonhitx = ddpos1.x();
+ muonhity = ddpos1.y();
+ muonhitz = map_muonhit[i];
+ // store to raw hit
+ edm4hep::MutableTrackerHit trkHit;
+ int true_pdgid = 0;
+ if(std::abs(pdgid1) == 13 && std::abs(map_pdgid[i]) == 13)
+ {
+ true_pdgid = 1;//pdgid->true muon hit
+ }
+ if(map_pdgid[i]==-1)
+ {
+ true_pdgid = 2;//ghost point
+ }
+ pos = edm4hep::Vector3d(muonhitx*10,muonhity*10,muonhitz*10);//mm
+ trkHit.setEDep(ADC);
+ trkHit.setQuality(true_pdgid);
+ trkHit.setCellID(cellid1);
+ trkHit.setPosition(pos);
+ trkhitVec->push_back( trkHit );
+ }
+ }
+ // write for out ntuple
+ if(_writeNtuple)
+ {
+ cell_cellid[n_cell] = cellid;
+ cell_edep[n_cell] = Edep;
+ cell_adc[n_cell] = ADC;
+ cell_posx[n_cell] = pos[0];
+ cell_posy[n_cell] = pos[1];
+ cell_posz[n_cell] = pos[2];
+ cell_layer[n_cell] = layer;
+ cell_slayer[n_cell] = slayer;
+ cell_strip[n_cell] = strip;
+ cell_fe[n_cell] = Fe;
+ cell_env[n_cell] = Env;
+ n_cell++;
+ }
+ }
+ if(_writeNtuple)
+ {
+ m_tree->Fill();
+ }
}
- ADC = rand_muon.Landau(ADCmean,7.922);
- if(ADC<0)
+ if(muonmode == 1)//muon endcap run
{
- ADC = 0;
- }
+ Clear();
+ if(m_MuonMode == 1) continue;
+ trkhitVec = m_outputMuonEndcap.createAndPut();
+ const edm4hep::SimTrackerHitCollection* STHCol = nullptr;
+ try
+ {
+ STHCol = m_inputMuonEndcap.get();
+ }
+ catch(GaudiException &e)
+ {
+ debug() << "Collection " << m_inputMuonEndcap.fullKey() << " is unavailable in event " << m_nEvt << endmsg;
+ continue;
+ }
+ if(STHCol->size()==0) continue;
+ debug() << m_inputMuonEndcap.fullKey() << " has SimTrackerHit "<< STHCol->size() << endmsg;
+
+ 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,
+ 2.44, 2.42, 2.40, 2.38, 2.36, 2.33, 2.31, 2.28, 2.26, 2.23,
+ 2.20, 2.17, 2.14, 2.11, 2.07, 2.04, 2.00, 1.97, 1.93, 1.89,
+ 1.84, 1.80, 1.75, 1.70, 1.65, 1.60, 1.54, 1.48, 1.42, 1.35,
+ 1.28, 1.20, 1.12, 1.02, 0.92, 0.80, 0.65, 0.46, 2.20, 2.20,
+ 2.20, 2.20, 2.20, 2.20, 2.20, 2.21, 2.21, 2.21, 2.21, 2.22,
+ 2.22, 2.22, 2.23, 2.23, 2.23, 2.24, 2.24, 2.25, 2.25, 2.26,
+ 2.26, 2.27, 2.28, 2.28, 2.29, 2.30, 2.31, 2.32, 2.32, 2.33,
+ 2.34, 2.35, 2.36, 2.38, 2.39, 2.40, 2.41, 2.43, 2.44, 2.45,
+ 2.47, 2.49, 2.50, 2.52, 2.54, 2.56, 2.58, 2.60, 2.62, 2.65,
+ 2.67, 2.70, 2.73, 2.76, 2.79, 2.82, 2.86, 2.90, 2.94, 2.99,
+ 3.04, 3.10, 3.16, 3.23, 3.31, 3.41, 3.53, 3.70, 4.14, 4.12,
+ 4.09, 4.06, 4.03, 4.00, 3.97, 3.94, 3.91, 3.87, 3.84, 3.81,
+ 3.77, 3.74, 3.70, 3.67, 3.63, 3.59, 3.55, 3.51, 3.47, 3.43,
+ 3.38, 3.34, 3.30, 3.25, 3.20, 3.15, 3.10, 3.05, 3.00, 2.95,
+ 2.89, 2.83, 2.77, 2.71, 2.65, 2.58, 2.52, 2.45, 2.37, 2.30,
+ 2.22, 2.14, 2.05, 1.96, 1.86, 1.76, 1.65, 1.53, 1.40, 1.25,
+ 1.09, 0.89, 0.63};
+
+
+ // define variables to be repeatedly used later
+ std::array<unsigned long long, 2> key;
+ int pdgid;
+ unsigned long long cellid, abspdgid;
+ double Edep, ADC, ADCmean;
+ int layer, slayer, strip, Env, Endcap;
+ edm4hep::Vector3d pos;
+ dd4hep::Position ddpos;
+ // loop over all hits
+ for(auto simhit : *STHCol)
+ {
+ auto mcp = simhit.getMCParticle();
+ cellid = simhit.getCellID();
+ pdgid = mcp.getPDG();
+ abspdgid = std::abs(mcp.getPDG());
+ key = std::array<unsigned long long, 2>{cellid, abspdgid};
+ Edep = simhit.getEDep();//GeV
+ layer = m_decoder2->get(cellid, "Layer");
+ slayer = m_decoder2->get(cellid, "Superlayer");
+ strip = m_decoder2->get(cellid, "Stripe");
+ Endcap = m_decoder2->get(cellid, "Endcap");
+ Env = m_decoder2->get(cellid, "Env");
+ pos = simhit.getPosition();
+ ddpos = m_cellIDConverter->position(cellid);
+ debug() << "Position:: " << ddpos.x() << " " << ddpos.y() << " " << ddpos.z() << endmsg;
+
+ // if not satisfy energy requirement, skip this hit
+ if ( Edep>0.1 || Edep<0.000001 ) continue;
+
+ //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));
+ strip--; // strip id begin at 1
+ double hit_sipm_length = endcap_strip_length[strip] * 0.5 * 100 - hit_strip_length;//cm
+
+ // 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 = rand_muon.Landau(ADCmean,7.922);
+ if(ADC<0)
+ {
+ ADC = 0;
+ }
+
+ 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] = Endcap;
+ map_cell_env[key] = Env;
+ map_cell_pos[key] = pos;
+ map_cell_pdgid[key] = pdgid;
+ }
- // fill the cell loop
- map_cell_edep[cellid] += Edep;
- map_cell_adc[cellid] += ADC;
- map_cell_layer[cellid] = layer;
- map_cell_slayer[cellid] = slayer;
- map_cell_strip[cellid] = strip;
- map_cell_fe[cellid] = Fe;
- map_cell_env[cellid] = Env;
-
- // write to tree
- if(_writeNtuple)
- {
- hit_cellid[n_hit] = cellid;
- hit_posx[n_hit] = pos[0];
- hit_posy[n_hit] = pos[1];
- hit_posz[n_hit] = pos[2];
- hit_edep[n_hit] = Edep;
- hit_layer[n_hit] = layer;
- hit_slayer[n_hit] = slayer;
- hit_strip[n_hit] = strip;
- hit_fe[n_hit] = Fe;
- hit_env[n_hit] = Env;
- n_hit++;
- }
-
- }
+ std::array<unsigned long long, 2> key1;
+ unsigned long long cellid1;
+ int layer1, slayer1, strip1, Fe1, Env1, pdgid1;
+ edm4hep::Vector3d pos1;
+ dd4hep::Position ddpos1;
- // loop over all cells
- for (const auto& item : map_cell_edep)
- {
- cellid = item.first;
- Edep = map_cell_edep[cellid];
- ADC = map_cell_adc[cellid];
- layer = map_cell_layer[cellid];
- slayer = map_cell_slayer[cellid];
- strip = map_cell_strip[cellid];
- Fe = map_cell_fe[cellid];
- Env = map_cell_env[cellid];
-
- if ( Edep < 0.0005 ) continue; // skip cell energy < 50 MeV ?? Why?
- //if ( Edep > 0.005 ) continue; // skip cell energy > 500 MeV ?? Why?
- //if ( ADC > 10000 ) continue; // skip cell energy > 1000 ADC counts? why?
-
- // simulate hit efficiency: random > hiteff, skip the hit
- if( rand_muon.Uniform(0, 1) > m_hitEff) continue;
-
- // store to raw hit
- edm4hep::MutableTrackerHit trkHit;
- ddpos = m_cellIDConverter->position(cellid);
- pos = edm4hep::Vector3d(ddpos.x(),ddpos.y(),ddpos.z());
- trkHit.setEDep(ADC);
- trkHit.setCellID(cellid);
- trkHit.setPosition(pos);
- trkhitVec->push_back( trkHit );
-
-
- // write for out ntuple
- if(_writeNtuple)
- {
- cell_cellid[n_cell] = cellid;
- cell_edep[n_cell] = Edep;
- cell_adc[n_cell] = ADC;
- cell_posx[n_cell] = pos[0];
- cell_posy[n_cell] = pos[1];
- cell_posz[n_cell] = pos[2];
- cell_layer[n_cell] = layer;
- cell_slayer[n_cell] = slayer;
- cell_strip[n_cell] = strip;
- cell_fe[n_cell] = Fe;
- cell_env[n_cell] = Env;
- n_cell++;
- }
+ std::array<unsigned long long, 2> key2;
+ unsigned long long cellid2;
+ int layer2, slayer2, strip2, Fe2, Env2, pdgid2;
+ edm4hep::Vector3d pos2;
+ dd4hep::Position ddpos2;
- }
+ if(map_cell_edep.size()==0) continue;
- if(_writeNtuple)
- {
- m_tree->Fill();
+ for(const auto& item : map_cell_edep)
+ {
+ key = item.first;
+ if (map_cell_edep[key] < m_EdepMin || rand_muon.Uniform(0, 1)>m_hitEff || map_cell_adc[key] == 0)
+ {
+ map_cell_edep[key] = 0;
+ }
+ }
+
+ double muonhitx, muonhity, muonhitz;
+ double Hit_max, Hit_min;
+
+ // loop over all cells
+ for (const auto& item1 : map_cell_edep)
+ {
+ key1 = item1.first;
+ cellid1 = key1[0];
+ Edep = map_cell_edep[key1];
+ if ( Edep < 0.0001 ) continue;
+ int anotherlayer_cell_num = 0;
+
+ ADC = map_cell_adc[key1];
+ layer1 = map_cell_layer[key1];
+ slayer1 = map_cell_slayer[key1];
+ strip1 = map_cell_strip[key1];
+ Fe1 = map_cell_fe[key1];//endcap
+ Env1 = map_cell_env[key1];
+ pos1 = map_cell_pos[key1];
+ pdgid1 = map_cell_pdgid[key1];
+ ddpos1 = m_cellIDConverter->position(cellid1);
+
+ for (const auto& item2 : map_cell_edep)
+ {
+ key2 = item2.first;
+ cellid2 = key2[0];
+ Edep = map_cell_edep[key2];
+ if ( Edep < 0.0001 ) continue;
+
+ layer2 = map_cell_layer[key2];
+ slayer2 = map_cell_slayer[key2];
+ strip2 = map_cell_strip[key2];
+ Fe2 = map_cell_fe[key2];//endcap
+ Env2 = map_cell_env[key2];
+ pos2 = map_cell_pos[key2];
+ pdgid2 = map_cell_pdgid[key2];
+ ddpos2 = m_cellIDConverter->position(cellid2);
+
+ if(Fe1!=Fe2 || slayer1!=slayer2 || layer1==layer2 || layer1==2) continue;
+ if(std::abs(Env2-Env1) == 2) continue;
+ Hit_min = pos1[0]*0.1-12;
+ Hit_max = pos1[0]*0.1+12;
+ if(ddpos2.x()<Hit_min || ddpos2.x()>Hit_max) continue;
+ anotherlayer_cell_num++;
+ map_muonhit[anotherlayer_cell_num] = ddpos2.x();
+ 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) continue;
+ map_pdgid[anotherlayer_cell_num] = pdgid2;
+ }
+ if(anotherlayer_cell_num == 0) continue;
+ if(anotherlayer_cell_num > 0)
+ {
+ for(int i=1; i<=anotherlayer_cell_num; i++)
+ {
+ muonhitx = map_muonhit[i];
+ muonhity = ddpos1.y();
+ muonhitz = ddpos1.z();
+ // store to raw hit
+ edm4hep::MutableTrackerHit trkHit;
+ int true_pdgid = 0;
+ if(std::abs(pdgid1) == 13 && std::abs(map_pdgid[i]) == 13)
+ {
+ true_pdgid = 1;//pdgid->true muon hit
+ }
+ if(map_pdgid[i]==-1)
+ {
+ true_pdgid = 2;//ghost point
+ }
+ pos = edm4hep::Vector3d(muonhitx*10,muonhity*10,muonhitz*10);
+ trkHit.setEDep(ADC);
+ trkHit.setQuality(true_pdgid);
+ trkHit.setCellID(cellid1);
+ trkHit.setPosition(pos);
+ trkhitVec->push_back( trkHit );
+ }
+ }
+ }
+ }
}
+
m_nEvt++;
return StatusCode::SUCCESS;
}
@@ -321,5 +610,8 @@ void MuonDigiAlg::Clear()
map_cell_strip.clear();
map_cell_fe.clear();
map_cell_env.clear();
-
+ map_cell_pos.clear();
+ map_muonhit.clear();
+ map_cell_pdgid.clear();
+ map_pdgid.clear();
}
diff --git a/Digitization/DigiMuon/src/MuonDigiAlg.h b/Digitization/DigiMuon/src/MuonDigiAlg.h
index 93d6d21c9b1cfe46c046382fb671efd2999ba443..f5d0d412680e958cca69788dc9768121a943e45d 100644
--- a/Digitization/DigiMuon/src/MuonDigiAlg.h
+++ b/Digitization/DigiMuon/src/MuonDigiAlg.h
@@ -61,18 +61,22 @@ class MuonDigiAlg : public GaudiAlgorithm
TFile* m_wfile;
SmartIF<IGeomSvc> m_geosvc;
- dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
+ dd4hep::DDSegmentation::BitFieldCoder* m_decoder1;
+ dd4hep::DDSegmentation::BitFieldCoder* m_decoder2;
dd4hep::rec::CellIDPositionConverter* m_cellIDConverter;
dd4hep::Detector* m_dd4hep;
- std::map<unsigned long long, double> map_cell_edep;
- std::map<unsigned long long, double> map_cell_adc;
- std::map<unsigned long long, int> map_cell_layer;
- std::map<unsigned long long, int> map_cell_slayer;
- std::map<unsigned long long, int> map_cell_strip;
- std::map<unsigned long long, int> map_cell_fe;
- std::map<unsigned long long, int> map_cell_env;
-
+ 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>, 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;
+ std::map<std::array<unsigned long long, 2>, int> map_cell_fe;
+ std::map<std::array<unsigned long long, 2>, int> map_cell_env;
+ std::map<std::array<unsigned long long, 2>, int> map_cell_pdgid;
+ std::map<std::array<unsigned long long, 2>, edm4hep::Vector3d> map_cell_pos;
+ std::map<int, double> map_muonhit;
+ std::map<int, int> map_pdgid;
#define MAX_SIZE 10000
@@ -106,17 +110,21 @@ class MuonDigiAlg : public GaudiAlgorithm
Gaudi::Property<int> _writeNtuple{this, "WriteNtuple", 1, "Write ntuple"};
Gaudi::Property<std::string> _filename{this, "OutFileName", "testout.root", "Output file name"};
- Gaudi::Property<double> m_hitEff{ this, "hitEff", 0.9, "Efficiency of a single hit on a Strip" };
+ Gaudi::Property<double> m_MuonMode{ this, "MuonMode", 0, "Control the Muon section used in Digi" };// 1->Barrel; 2->Endcap; 0->Both;
+ Gaudi::Property<double> m_hitEff{ this, "SiPMEff", 0.95, "Efficiency of a single hit on a Strip" };
+ Gaudi::Property<double> m_EdepMin{ this, "EdepMin", 0.0001, "Minimum Edep of a mip" };
// Input collections
DataHandle<edm4hep::SimTrackerHitCollection> m_inputMuonBarrel{"MuonBarrelHitsCollection", Gaudi::DataHandle::Reader, this};
+ DataHandle<edm4hep::SimTrackerHitCollection> m_inputMuonEndcap{"MuonEndcapHitsCollection", Gaudi::DataHandle::Reader, this};
// Output collections
DataHandle<edm4hep::TrackerHitCollection> m_outputMuonBarrel{"MuonBarrelTrackerHits", Gaudi::DataHandle::Writer, this};
+ DataHandle<edm4hep::TrackerHitCollection> m_outputMuonEndcap{"MuonEndcapTrackerHits", Gaudi::DataHandle::Writer, this};
+
//DataHandle<edm4hep::MCRecoTrackerAssociationCollection> m_assMuonBarrel{"MuonBarrelTrackerHitAssociationCollection", Gaudi::DataHandle::Writer, this};
edm4hep::TrackerHitCollection* trkhitVec;
-
int m_nEvt=0;
};
#endif