diff --git a/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyMuon.xml b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyMuon.xml index 2f602d4aa20aa345999effb948c455a6c178092f..4c3609c88bdc02502e404cd2f72150dc84685b33 100644 --- a/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyMuon.xml +++ b/Detector/DetCRD/compact/TDR_o1_v01/TDR_o1_v01-onlyMuon.xml @@ -3,39 +3,60 @@ xmlns:xs="http://www.w3.org/2001/XMLSchema" xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd"> <info name="TDR_o1_v01" - title="CepC reference detctor with coil inside Hcal, pixel SIT/SET" - author="C.D.Fu, " + title="CepC reference detctor for TDR" + author="" url="http://cepc.ihep.ac.cn" status="developing" version="v01"> - <comment>CepC reference detector simulation models used for detector study </comment> + <comment>CepC reference detector simulation models used for TDR </comment> </info> - + <includes> <gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/elements.xml"/> - <gdmlFile ref="../materials.xml"/> + <gdmlFile ref="../CRD_common_v02/materials.xml"/> </includes> - + <define> - <!--world size--> - <constant name="Muon_world_size" value="25*m"/> - <constant name="world_x" value="Muon_world_size"/> - <constant name="world_y" value="Muon_world_size"/> - <constant name="world_z" value="Muon_world_size"/> + <constant name="world_size" value="10*m"/> + <constant name="world_x" value="world_size"/> + <constant name="world_y" value="world_size"/> + <constant name="world_z" value="world_size"/> + + <include ref="${DD4hepINSTALL}/DDDetectors/compact/detector_types.xml"/> </define> <include ref="./TDR_Dimensions_v01_01.xml"/> - <!--include ref="../CRD_common_v01/Beampipe_v01_01.xml"/--> - <!--include ref="../CRD_common_v01/VXD_v01_01.xml"/--> - <!--include ref="../CRD_common_v01/FTD_SkewRing_v01_01.xml"/--> - <!--include ref="../CRD_common_v01/SIT_SimplePixel_v01_01.xml"/--> - <!--include ref="../CRD_common_v01/DC_Simple_v01_05.xml"/--> + <!--include ref="../CRD_common_v02/Beampipe_v01_04.xml"/--> + <!--preliminary vertex and tracker, to update/--> + <!--include ref="../CRD_common_v02/VXD_StaggeredLadder_v02_01.xml"/--> + <!--include ref="../CRD_common_v02/VXD_Composite_v01_01.xml"/--> + <!--include ref="../CRD_common_v02/FTD_SkewRing_v01_07.xml"/--> + <!--include ref="../CRD_common_v02/SIT_SimplePixel_v01_04.xml"/--> + <!--include ref="../CRD_common_v02/SIT_StaggeredStave_v02.xml"/--> + <!--include ref="../CRD_common_v01/TPC_Simple_v10_02.xml"/--> + <!--use 10 rows clustering version/--> + <!--include ref="../CRD_common_v02/TPC_ModularEndcap_o1_v02.xml"/--> <!--include ref="../CRD_common_v01/SET_SimplePixel_v01_01.xml"/--> - <!--include ref="../CRD_common_v01/Ecal_Crystal_Barrel_v01_01.xml"/--> + <!--include ref="../CRD_common_v01/OTKBarrel_v01_01.xml"/--> + <!--include ref="../CRD_common_v01/OTKEndcap_v01_01.xml"/--> + + <!--include ref="../CRD_common_v01/Ecal_Crystal_Barrel_v01_02.xml"/--> <!--include ref="../CRD_common_v01/Ecal_Crystal_Endcap_v01_01.xml"/--> - <include ref="../CRD_common_v01/Muon_Endcap_v01_01.xml"/> - <include ref="../CRD_common_v01/Muon_Barrel_v01_02.xml"/> + <!--include ref="../CRD_common_v01/SHcalGlass_Barrel_v05.xml"/--> + <!--include ref="../CRD_common_v01/SHcalGlass_Endcaps_v01.xml"/--> + + <!--Lumical to update--> + <!--include ref="../CRD_common_v01/Lumical_o1_v01.xml"/--> + <!--preliminary Magnet, to update/--> + <!--include ref="../CRD_common_v02/Coil_Simple_v01_02.xml"/--> + <!--preliminary Muon, obselete/--> + <!--include ref="../CRD_common_v02/Yoke_Polyhedra_Barrel_v01_01.xml"/--> + <!--include ref="../CRD_common_v02/Yoke_Polyhedra_Endcaps_v01_01.xml"/--> + + <!--muon detector--> + <include ref="../CRD_common_v01/Muon_Barrel_v01_01.xml"/> + <include ref="../CRD_common_v01/Muon_Endcap_v01_02.xml"/> <fields> <field name="InnerSolenoid" type="solenoid" diff --git a/Detector/DetCRD/src/Muon/Muon_Barrel_v01_01.cpp b/Detector/DetCRD/src/Muon/Muon_Barrel_v01_01.cpp index 90b785921aaac17c46ce057a30d87610154e3862..77204c116772d9394483ce4828acba16839fb2d3 100644 --- a/Detector/DetCRD/src/Muon/Muon_Barrel_v01_01.cpp +++ b/Detector/DetCRD/src/Muon/Muon_Barrel_v01_01.cpp @@ -74,7 +74,7 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector, { xml_coll_t dcSuperlayer(x_Fe,Unicode("superlayer")); xml_comp_t x_superlayer = dcSuperlayer; - for(int i2 = 0; i2 < x_superlayer.id(); i2++) + for(int i2 = 0; i2 < x_superlayer.id(); i2++) // FIXME:: 8 layers start from 0 { std::string superlayer_name = x_superlayer.nameStr() + dd4hep::_toString(i1,"_%d") + dd4hep::_toString(i2,"_%d"); std::string num_name = "Muon_barrel_strip_num" + dd4hep::_toString(i2,"_%d"); diff --git a/Digitization/DigiMuon/src/MuonDigiAlg.cpp b/Digitization/DigiMuon/src/MuonDigiAlg.cpp index 4e860338d20f63bacfb920c6ce4f3a0e814fdfd5..17703c2f861c3cb1ca6dfe50d4c121a1b1838ec5 100644 --- a/Digitization/DigiMuon/src/MuonDigiAlg.cpp +++ b/Digitization/DigiMuon/src/MuonDigiAlg.cpp @@ -23,8 +23,6 @@ #include <DDRec/CellIDPositionConverter.h> #include "TVector3.h" -//#include "TRandom3.h" -//#include <TRandom.h> #include <math.h> #include <cmath> #include <iostream> @@ -61,14 +59,30 @@ StatusCode MuonDigiAlg::initialize() { std::string s_outfile = _filename; m_wfile = new TFile(s_outfile.c_str(), "recreate"); - hitloop = new TTree("h1", "muonbarreldigi_hitloop"); - eventloop = new TTree("e1", "muonbarreldigi_eventloop"); - hitloop->Branch("MuonBarrel_cellid", &MuonBarrel_cellid); - hitloop->Branch("MuonBarrel_posx", &MuonBarrel_posx); - hitloop->Branch("MuonBarrel_posy", &MuonBarrel_posy); - hitloop->Branch("MuonBarrel_posz", &MuonBarrel_posz); - eventloop->Branch("MuonBarrel_stripcellid", &MuonBarrel_stripcellid); - eventloop->Branch("MuonBarrel_stripedep", &MuonBarrel_stripedep); + m_tree = new TTree("tree", "muon digi tree"); + m_tree->Branch("n_hit", &n_hit, "n_hit/I"); + m_tree->Branch("hit_cellid", hit_cellid, "hit_cellid[n_hit]/l"); + m_tree->Branch("hit_posx", hit_posx, "hit_posx[n_hit]/D"); + m_tree->Branch("hit_posy", hit_posy, "hit_posy[n_hit]/D"); + m_tree->Branch("hit_posz", hit_posz, "hit_posz[n_hit]/D"); + m_tree->Branch("hit_edep", hit_edep, "hit_edep[n_hit]/D"); + m_tree->Branch("hit_layer", hit_layer, "hit_layer[n_hit]/I"); + m_tree->Branch("hit_slayer", hit_slayer, "hit_slayer[n_hit]/I"); + m_tree->Branch("hit_strip", hit_strip, "hit_strip[n_hit]/I"); + m_tree->Branch("hit_fe", hit_fe, "hit_fe[n_hit]/I"); + m_tree->Branch("hit_env", hit_env, "hit_env[n_hit]/I"); + m_tree->Branch("n_cell", &n_cell, "n_cell/I"); + m_tree->Branch("cell_cellid", cell_cellid, "cell_cellid[n_cell]/l"); + m_tree->Branch("cell_edep", cell_edep, "cell_edep[n_cell]/D"); + m_tree->Branch("cell_adc", cell_adc, "cell_adc[n_cell]/D"); + m_tree->Branch("cell_posx", cell_posx, "cell_posx[n_cell]/D"); + m_tree->Branch("cell_posy", cell_posy, "cell_posy[n_cell]/D"); + m_tree->Branch("cell_posz", cell_posz, "cell_posz[n_cell]/D"); + m_tree->Branch("cell_layer", cell_layer, "cell_layer[n_cell]/I"); + m_tree->Branch("cell_slayer", cell_slayer, "cell_slayer[n_cell]/I"); + 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"); } @@ -100,6 +114,8 @@ StatusCode MuonDigiAlg::initialize() return StatusCode::FAILURE; } + debug() << "m_hitEff: " << m_hitEff << endmsg; + info() << "MuonDigiAlg::initialized" << endmsg; return GaudiAlgorithm::initialize(); } @@ -113,7 +129,6 @@ StatusCode MuonDigiAlg::execute() trkhitVec = m_outputMuonBarrel.createAndPut(); //auto assVec = m_assMuonBarrel.createAndPut(); - const edm4hep::SimTrackerHitCollection* STHCol = nullptr; try { @@ -130,111 +145,141 @@ StatusCode MuonDigiAlg::execute() int strip_length[8] = {17, 29, 41, 53, 65, 77, 89, 101}; + // 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) { - auto cellid = simhit.getCellID(); - auto Edep = simhit.getEDep();//GeV - int layer = m_decoder->get(cellid, "Layer"); - int slayer = m_decoder->get(cellid, "Superlayer"); - int strip = m_decoder->get(cellid, "Stripe"); - int Fe = m_decoder->get(cellid, "Fe"); - int Env = m_decoder->get(cellid, "Env"); - auto& pos = simhit.getPosition(); - - dd4hep::Position hitpos = m_cellIDConverter->position(cellid); - debug() << "Position:: " << hitpos.x() << " " << hitpos.y() << " " << hitpos.z() << endmsg; - - int testtemp = 0; - if(Edep<0.01&&Edep>0.0001) + 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) { - double hit_strip_length = std::sqrt((hitpos.x()-pos[0]*0.1)*(hitpos.x()-pos[0]*0.1)+(hitpos.y()-pos[1]*0.1)*(hitpos.y()-pos[1]*0.1)+(hitpos.z()-pos[2]*0.1)*(hitpos.z()-pos[2]*0.1)); - double hit_sipm_length; - if(layer == 1) - { - int tempnum = Env + slayer; - if(tempnum%2 == 0) - { - hit_sipm_length = 115 * 2 - hit_strip_length; - } - else - { - hit_sipm_length = 106 * 2 - hit_strip_length; - } - } - if(layer == 2) + int tempnum = Env + slayer; + if(tempnum%2 == 0) { - hit_sipm_length = strip_length[slayer-1] * 2 - hit_strip_length; + hit_sipm_length = 115 * 2 - hit_strip_length; } - - double ADCmean = 34*Edep*47.09/(23*0.00141);//mV - if(hit_sipm_length>10) + else { - ADCmean = (16.0813*std::exp(-1*hit_sipm_length/50.8147)+19.5474)*Edep*47.09/(23*0.00141);//mV + hit_sipm_length = 106 * 2 - hit_strip_length; } - double ADC = rand_muon.Landau(ADCmean,7.922); - if(ADC<0) - { - ADC = 0; - } - //std::cout<<hit_sipm_length<<" "<<Edep<<" "<<ADCmean<<" "<<ADC<<std::endl; + } + if(layer == 2) + { + hit_sipm_length = strip_length[slayer-1] * 2 - hit_strip_length; + } - if(cellidtemp.size()>0) - { - for(int i=0; i<cellidtemp.size(); i++) - { - if(cellid == cellidtemp[i]) - { - edeptemp[i]+=Edep; - ADCtemp[i]+=ADC; - testtemp++; - break; - } - } - } - if(testtemp==0) - { - cellidtemp.push_back(cellid); - edeptemp.push_back(Edep); - ADCtemp.push_back(ADC); - } - if(_writeNtuple) - { - MuonBarrel_cellid.push_back(cellid); - MuonBarrel_posx.push_back(pos[0]); - MuonBarrel_posy.push_back(pos[1]); - MuonBarrel_posz.push_back(pos[2]); - } + // 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[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++; + } + } - for(int i=0; i<cellidtemp.size(); i++) + // loop over all cells + for (const auto& item : map_cell_edep) { - if(edeptemp[i]>0.0005&&edeptemp[i]<0.005&&ADCtemp[i]<1000) + 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) { - double random_num = rand_muon.Uniform(0, 1); - if(random_num>0.1) - { - if(_writeNtuple) - { - MuonBarrel_stripcellid.push_back(cellidtemp[i]); - MuonBarrel_stripedep.push_back(edeptemp[i]); - } - edm4hep::MutableTrackerHit trkHit; - dd4hep::Position hitpos1 = m_cellIDConverter->position(cellidtemp[i]); - edm4hep::Vector3d position_strip(hitpos1.x(),hitpos1.y(),hitpos1.z()); - //std::cout<<"Position:: "<<hitpos.x()<<" "<<hitpos.y()<<" "<<hitpos.z()<<std::endl; - trkHit.setEDep(ADCtemp[i]); - trkHit.setCellID(cellidtemp[i]); - trkHit.setPosition(position_strip); - trkhitVec->push_back( trkHit ); - } + 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) { - hitloop->Fill(); - eventloop->Fill(); + m_tree->Fill(); } m_nEvt++; @@ -247,10 +292,9 @@ StatusCode MuonDigiAlg::finalize() if(_writeNtuple) { m_wfile->cd(); - hitloop->Write(); - eventloop->Write(); + m_tree->Write(); m_wfile->Close(); - delete m_wfile, hitloop, eventloop; + delete m_wfile, m_tree; } info() << "Processed " << m_nEvt << " events " << endmsg; return GaudiAlgorithm::finalize(); @@ -258,13 +302,14 @@ StatusCode MuonDigiAlg::finalize() void MuonDigiAlg::Clear() { - MuonBarrel_cellid.clear(); - MuonBarrel_posx.clear(); - MuonBarrel_posy.clear(); - MuonBarrel_posz.clear(); - MuonBarrel_stripcellid.clear(); - MuonBarrel_stripedep.clear(); - cellidtemp.clear(); - edeptemp.clear(); - ADCtemp.clear(); + n_hit = 0; + n_cell = 0; + map_cell_edep.clear(); + map_cell_adc.clear(); + map_cell_layer.clear(); + map_cell_slayer.clear(); + map_cell_strip.clear(); + map_cell_fe.clear(); + map_cell_env.clear(); + } diff --git a/Digitization/DigiMuon/src/MuonDigiAlg.h b/Digitization/DigiMuon/src/MuonDigiAlg.h index 1ebae24da24db1748b53c39fd287e6af663330f7..93d6d21c9b1cfe46c046382fb671efd2999ba443 100644 --- a/Digitization/DigiMuon/src/MuonDigiAlg.h +++ b/Digitization/DigiMuon/src/MuonDigiAlg.h @@ -38,7 +38,6 @@ #include <TTimeStamp.h> #include <ctime> #include <iostream> -//#include <TRandom.h> #include <TMath.h> #define PI 3.141592653 class MuonDigiAlg : public GaudiAlgorithm @@ -53,27 +52,12 @@ class MuonDigiAlg : public GaudiAlgorithm void Clear(); - // double EDepConvertADC(double edep){ - // rand.SetSeed(0); - // randGen.SetSeed(0); - // double adc = 0; - // Double_t random_num = rand.Uniform(0, 1); - // if(random_num>0.1){ - // double LandauRandom = randGen.Landau(47.09,7.922); - // adc = LandauRandom * edep / 0.00141; - // } - // return random_num; - // } - - - protected: TRandom3 rand_muon; - TTree* hitloop; - TTree* eventloop; + TTree* m_tree; TFile* m_wfile; SmartIF<IGeomSvc> m_geosvc; @@ -81,22 +65,49 @@ class MuonDigiAlg : public GaudiAlgorithm dd4hep::rec::CellIDPositionConverter* m_cellIDConverter; dd4hep::Detector* m_dd4hep; - std::vector<unsigned long long> MuonBarrel_cellid; - std::vector<double> MuonBarrel_posx; - std::vector<double> MuonBarrel_posy; - std::vector<double> MuonBarrel_posz; - std::vector<unsigned long long> cellidtemp; - std::vector<double> edeptemp; - std::vector<double> ADCtemp; - std::vector<unsigned long long> MuonBarrel_stripcellid; - std::vector<double> MuonBarrel_stripedep; - - + 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; + +#define MAX_SIZE 10000 + + int n_hit, n_cell; + + unsigned long long hit_cellid[MAX_SIZE]; + double hit_posx[MAX_SIZE]; + double hit_posy[MAX_SIZE]; + double hit_posz[MAX_SIZE]; + double hit_edep[MAX_SIZE]; + int hit_layer[MAX_SIZE]; + int hit_slayer[MAX_SIZE]; + int hit_strip[MAX_SIZE]; + int hit_fe[MAX_SIZE]; + int hit_env[MAX_SIZE]; + + unsigned long long cell_cellid[MAX_SIZE]; + double cell_edep[MAX_SIZE]; + double cell_adc[MAX_SIZE]; + double cell_posx[MAX_SIZE]; + double cell_posy[MAX_SIZE]; + double cell_posz[MAX_SIZE]; + int cell_layer[MAX_SIZE]; + int cell_slayer[MAX_SIZE]; + int cell_strip[MAX_SIZE]; + int cell_fe[MAX_SIZE]; + int cell_env[MAX_SIZE]; + + + 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" }; - mutable Gaudi::Property<int> _writeNtuple{this, "WriteNtuple", 1, "Write ntuple"}; - mutable Gaudi::Property<std::string> _filename{this, "OutFileName", "testout.root", "Output file name"}; // Input collections DataHandle<edm4hep::SimTrackerHitCollection> m_inputMuonBarrel{"MuonBarrelHitsCollection", Gaudi::DataHandle::Reader, this};