diff --git a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
index 6bf8a8e7cf1e59bea6b6bb6abf0ef9f3900425b4..15c8e1c942aa69ff3d61edc1b4f03a1461c179d0 100644
--- a/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
+++ b/Detector/DetCRD/scripts/TDR_o1_v01/tracking.py
@@ -149,6 +149,10 @@ digiMuon.MuonBarrelTrackerHits = "MuonBarrelTrackerHits"
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
#########################################
################
diff --git a/Digitization/DigiMuon/src/MuonDigiAlg.cpp b/Digitization/DigiMuon/src/MuonDigiAlg.cpp
index 347cb6c8c82d48cfb99437aec490509d1e36ad7d..0523a7e5d58b0f094c85a4c040aa8ed6e15e6266 100644
--- a/Digitization/DigiMuon/src/MuonDigiAlg.cpp
+++ b/Digitization/DigiMuon/src/MuonDigiAlg.cpp
@@ -55,7 +55,7 @@ MuonDigiAlg::MuonDigiAlg(const std::string& name, ISvcLocator* svcLoc)
StatusCode MuonDigiAlg::initialize()
{
// set rand seed;
- rand_muon.SetSeed(123456);
+ m_randSvc = service<IRndmGenSvc>("RndmGenSvc");
//
if(_writeNtuple)
@@ -86,6 +86,8 @@ StatusCode MuonDigiAlg::initialize()
m_tree->Branch("cell_strip", cell_strip, "cell_strip[n_cell]/I");
m_tree->Branch("cell_fe", cell_fe, "cell_fe[n_cell]/I");
m_tree->Branch("cell_env", cell_env, "cell_env[n_cell]/I");
+ m_tree->Branch("cell_superlayernumber", &cell_superlayernumber, "cell_superlayernumber/I");
+ m_tree->Branch("cell_pt", &cell_pt, "cell_pt/I");
}
@@ -108,19 +110,19 @@ StatusCode MuonDigiAlg::initialize()
return StatusCode::FAILURE;
}
- std::string readout_name1 = m_inputMuonBarrel.objKey();
- debug() << "Readout name: " << readout_name1 << endmsg;
- m_decoder1 = m_geosvc->getDecoder(readout_name1);
- if(!m_decoder1)
+ std::string readout_name_barrel = m_inputMuonBarrel.objKey();
+ debug() << "Readout name: " << readout_name_barrel << endmsg;
+ m_decoder_barrel = m_geosvc->getDecoder(readout_name_barrel);
+ if(!m_decoder_barrel)
{
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)
+ std::string readout_name_endcap = m_inputMuonEndcap.objKey();
+ debug() << "Readout name: " << readout_name_endcap << endmsg;
+ m_decoder_endcap = m_geosvc->getDecoder(readout_name_endcap);
+ if(!m_decoder_endcap)
{
error() << "Failed to get the decoder. " << endmsg;
return StatusCode::FAILURE;
@@ -135,451 +137,122 @@ StatusCode MuonDigiAlg::initialize()
StatusCode MuonDigiAlg::execute()
{
- for(int muonmode = 0; muonmode<2; muonmode++)
+ const edm4hep::SimTrackerHitCollection* STHCol;
+//MuonBarrelRun:
+ Clear();
+ trkhitVec = m_outputMuonBarrel.createAndPut();
+ STHCol = nullptr;
+ try
{
- if(muonmode == 0)//muon barrel run
+ STHCol = m_inputMuonBarrel.get();
+ }
+ catch(GaudiException &e)
+ {
+ 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();
+
+ 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)
{
- Clear();
- if(m_MuonMode == 2) continue;
- trkhitVec = m_outputMuonBarrel.createAndPut();
- //auto assVec = m_assMuonBarrel.createAndPut();
-
- const edm4hep::SimTrackerHitCollection* STHCol = nullptr;
- try
- {
- STHCol = m_inputMuonBarrel.get();
- }
- catch(GaudiException &e)
- {
- debug() << "Collection " << m_inputMuonBarrel.fullKey() << " is unavailable in event " << m_nEvt << endmsg;
- continue;
- }
- 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
- 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;
- }
- }
-
- 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();
- }
+ 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]);
}
- if(muonmode == 1)//muon endcap run
+
+ if(anotherlayer_cell_num == 0) continue;
+ for(int i=1; i<=anotherlayer_cell_num; i++)
{
- 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;
- }
-
- 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;
-
- 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 );
- }
- }
- }
+ 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);
}
}
+ Save_onelayer_signal(trkhitVec);
+MuonEndcapRun:
+ Clear();
+ trkhitVec = m_outputMuonEndcap.createAndPut();
+ STHCol = nullptr;
+ try
+ {
+ STHCol = m_inputMuonEndcap.get();
+ }
+ catch(GaudiException &e)
+ {
+ 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)
+ {
+ 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);
+ }
+ }
+ Save_onelayer_signal(trkhitVec);
m_nEvt++;
return StatusCode::SUCCESS;
@@ -615,3 +288,181 @@ void MuonDigiAlg::Clear()
map_cell_pdgid.clear();
map_pdgid.clear();
}
+
+void MuonDigiAlg::GetSimHit(edm4hep::SimTrackerHit _simhit, dd4hep::DDSegmentation::BitFieldCoder* _m_decoder)
+{
+ 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]);
+}
+
+double MuonDigiAlg::Gethit_sipm_length_Barrel()
+{
+ //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)
+ {
+ 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;
+ }
+ }
+ 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
+ return (strip_length[slayer-1] * 4)/2 - hit_strip_length;
+ }
+}
+
+double MuonDigiAlg::Gethit_sipm_length_Endcap()
+{
+ //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
+}
+
+void MuonDigiAlg::EdeptoADC()
+{
+ // 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)
+ {
+ ADC = 0;
+ }
+}
+
+void MuonDigiAlg::SaveData_mapcell()
+{
+ // 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;
+}
+
+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)
+ {
+ map_cell_edep[key] = 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)
+{
+ 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;
+ if(_ddposi<Hit_min || _ddposi>Hit_max) return;
+ anotherlayer_cell_num++;
+ 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;
+}
+
+int MuonDigiAlg::Get_true_pdgid(int _i, int _pdgid)
+{
+ 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
+}
+
+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
+ }
+}
+
+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 f5d0d412680e958cca69788dc9768121a943e45d..78a4e066f5fafe41db625ea65a91d553f197d37b 100644
--- a/Digitization/DigiMuon/src/MuonDigiAlg.h
+++ b/Digitization/DigiMuon/src/MuonDigiAlg.h
@@ -50,19 +50,32 @@ class MuonDigiAlg : public GaudiAlgorithm
virtual StatusCode execute() ;
virtual StatusCode finalize() ;
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 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 Save_onelayer_signal(edm4hep::TrackerHitCollection* _trkhitVec);
+ void Save_pos(int _i);
protected:
- TRandom3 rand_muon;
+
TTree* m_tree;
TFile* m_wfile;
SmartIF<IGeomSvc> m_geosvc;
- dd4hep::DDSegmentation::BitFieldCoder* m_decoder1;
- dd4hep::DDSegmentation::BitFieldCoder* m_decoder2;
+ dd4hep::DDSegmentation::BitFieldCoder* m_decoder_barrel;
+ dd4hep::DDSegmentation::BitFieldCoder* m_decoder_endcap;
dd4hep::rec::CellIDPositionConverter* m_cellIDConverter;
dd4hep::Detector* m_dd4hep;
@@ -105,14 +118,50 @@ class MuonDigiAlg : public GaudiAlgorithm
int cell_strip[MAX_SIZE];
int cell_fe[MAX_SIZE];
int cell_env[MAX_SIZE];
+ int cell_superlayernumber;
+ double cell_pt;
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_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_hitEff{ this, "SiPMEff", 1, "Efficiency of a single hit on a Strip" };
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"};
+
+ // 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,
+ 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};
// Input collections
@@ -123,7 +172,7 @@ class MuonDigiAlg : public GaudiAlgorithm
DataHandle<edm4hep::TrackerHitCollection> m_outputMuonEndcap{"MuonEndcapTrackerHits", Gaudi::DataHandle::Writer, this};
//DataHandle<edm4hep::MCRecoTrackerAssociationCollection> m_assMuonBarrel{"MuonBarrelTrackerHitAssociationCollection", Gaudi::DataHandle::Writer, this};
-
+ SmartIF<IRndmGenSvc> m_randSvc;
edm4hep::TrackerHitCollection* trkhitVec;
int m_nEvt=0;
};