diff --git a/Digitisers/DCHDigi/CMakeLists.txt b/Digitisers/DCHDigi/CMakeLists.txt index 33fa374effccc19a7fabd6fb841a3e34f368c067..e49e06929e580df953d2814b1cb40f1404c6c5c3 100644 --- a/Digitisers/DCHDigi/CMakeLists.txt +++ b/Digitisers/DCHDigi/CMakeLists.txt @@ -1,5 +1,7 @@ gaudi_subdir(DCHDigi v0r0) +find_package(Geant4 REQUIRED ui_all vis_all) +include(${Geant4_USE_FILE}) find_package(DD4hep COMPONENTS DDG4 REQUIRED) find_package(EDM4HEP REQUIRED ) message("EDM4HEP_INCLUDE_DIRS: ${EDM4HEP_INCLUDE_DIR}") diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp index f179943a32f9831adcdb00a6b638dfb63132f7d3..f77f0118a3aa138d2a064548b71253403dd217c7 100644 --- a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp +++ b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp @@ -8,8 +8,10 @@ #include "DD4hep/Detector.h" #include <DD4hep/Objects.h> +#include "G4ThreeVector.hh" +#include <array> #include <math.h> #include <cmath> #include <algorithm> @@ -77,29 +79,49 @@ StatusCode DCHDigiAlg::execute() double tot_edep = 0 ; double tot_length = 0 ; double tot_time = 0 ; - double tot_x = 0 ; - double tot_y = 0 ; - double tot_z = 0 ; + double pos_x = 0 ; + double pos_y = 0 ; + double pos_z = 0 ; int simhit_size = iter->second.size(); for(unsigned int i=0; i< simhit_size; i++) { - tot_edep += iter->second.at(i).getEDep();//GeV + tot_edep += iter->second.at(i).getEDep();//GeV + } + float min_distance = 999 ; + for(unsigned int i=0; i< simhit_size; i++) + { + G4ThreeVector west(0,0,0); + G4ThreeVector east(0,0,0); + G4ThreeVector pos(iter->second.at(i).getPosition()[0], iter->second.at(i).getPosition()[1], iter->second.at(i).getPosition()[2]); + G4ThreeVector numerator = (east-west).cross(west-pos) ; + G4ThreeVector denominator = (east-west) ; + float distance = numerator.mag()/denominator.mag() ; + std::cout<<"distance="<<distance<<std::endl; + if(distance < min_distance){ + min_distance = distance; + pos_x = pos.x(); + pos_y = pos.y(); + pos_z = pos.z(); + } tot_length += iter->second.at(i).getPathLength();//mm - tot_time += iter->second.at(i).getTime(); + /* tot_x += iter->second.at(i).getEDep()*iter->second.at(i).getPosition()[0]; tot_y += iter->second.at(i).getEDep()*iter->second.at(i).getPosition()[1]; tot_z += iter->second.at(i).getEDep()*iter->second.at(i).getPosition()[2]; + */ auto asso = AssoVec->create(); asso.setRec(trkHit); asso.setSim(iter->second.at(i)); - asso.setWeight(1.0/simhit_size); + asso.setWeight(iter->second.at(i).getEDep()/tot_edep); } - trkHit.setTime(tot_time/simhit_size); - trkHit.setEDep(tot_edep); + trkHit.setTime(min_distance*1e3/m_velocity);//m_velocity is um/ns, drift time in ns + trkHit.setEDep(tot_edep);// GeV trkHit.setEdx (tot_edep*1000/(tot_length/10) ); // MeV/cm, need check! - trkHit.setPosition (edm4hep::Vector3d(tot_x/tot_edep, tot_y/tot_edep, tot_z/tot_edep));//center mass + //trkHit.setPosition (edm4hep::Vector3d(tot_x/tot_edep, tot_y/tot_edep, tot_z/tot_edep));//center mass + trkHit.setPosition (edm4hep::Vector3d(pos_x, pos_y, pos_z));//position of closest sim hit + trkHit.setCovMatrix(std::array<float, 6>{m_res_x, 0, m_res_y, 0, 0, m_res_z});//cov(x,x) , cov(y,x) , cov(y,y) , cov(z,x) , cov(z,y) , cov(z,z) in mm } std::cout<<"output digi DCHhit size="<< Vec->size() <<std::endl; _nEvt ++ ; diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.h b/Digitisers/DCHDigi/src/DCHDigiAlg.h index e51c8c31615a630b3b446deb1a7d105e50aa761b..bdd80cfded1b3d0902f15befbbb1589ecbb7789d 100644 --- a/Digitisers/DCHDigi/src/DCHDigiAlg.h +++ b/Digitisers/DCHDigi/src/DCHDigiAlg.h @@ -43,8 +43,10 @@ protected: //float m_length; dd4hep::rec::CellIDPositionConverter* m_cellIDConverter; - //Gaudi::Property<float> m_scale { this, "Scale", 1 }; - //Gaudi::Property<float> m_resolution{ this, "Res", 0.01 }; + Gaudi::Property<float> m_res_x { this, "res_x", 0.11};//mm + Gaudi::Property<float> m_res_y { this, "res_y", 0.11};//mm + Gaudi::Property<float> m_res_z { this, "res_z", 1 };//mm + Gaudi::Property<float> m_velocity { this, "drift_velocity", 40};// um/ns // Input collections DataHandle<edm4hep::SimTrackerHitCollection> r_SimDCHCol{"DriftChamberHitsCollection", Gaudi::DataHandle::Reader, this};