From 616497296a8c1f9619cadd88d40e84bdf1469c8e Mon Sep 17 00:00:00 2001
From: myliu <201916234@mail.sdu.edu.cn>
Date: Fri, 4 Jun 2021 09:57:20 +0800
Subject: [PATCH] digi Alg
---
Detector/DetDriftChamber/compact/det.xml | 2 +-
.../DetSegmentation/GridDriftChamber.h | 50 +++++++++++
.../DetSegmentation/src/GridDriftChamber.cpp | 90 ++++++++++++++++++-
Digitisers/DCHDigi/src/DCHDigiAlg.cpp | 47 +++++++---
Digitisers/DCHDigi/src/DCHDigiAlg.h | 5 +-
5 files changed, 177 insertions(+), 17 deletions(-)
diff --git a/Detector/DetDriftChamber/compact/det.xml b/Detector/DetDriftChamber/compact/det.xml
index 4e981b02..a98e21bd 100644
--- a/Detector/DetDriftChamber/compact/det.xml
+++ b/Detector/DetDriftChamber/compact/det.xml
@@ -74,7 +74,7 @@
<limits>
<limitset name="DC_limits">
- <limit name="step_length_max" particles="*" value="0.5" unit="mm" />
+ <limit name="step_length_max" particles="*" value="0.1" unit="mm" />
</limitset>
</limits>
diff --git a/Detector/DetSegmentation/include/DetSegmentation/GridDriftChamber.h b/Detector/DetSegmentation/include/DetSegmentation/GridDriftChamber.h
index 7c1a0864..159addf4 100644
--- a/Detector/DetSegmentation/include/DetSegmentation/GridDriftChamber.h
+++ b/Detector/DetSegmentation/include/DetSegmentation/GridDriftChamber.h
@@ -61,6 +61,10 @@ public:
const VolumeID& aVolumeID) const;
virtual double distanceTrackWire(const CellID& cID, const TVector3& hit_start, const TVector3& hit_end) const;
virtual void cellposition(const CellID& cID, TVector3& Wstart, TVector3& Wend) const;
+ virtual TVector3 distanceClosestApproach(const CellID& cID, const TVector3& hitPos) const;
+ virtual TVector3 Line_TrackWire(const CellID& cID, const TVector3& hit_start, const TVector3& hit_end) const;
+ virtual TVector3 IntersectionTrackWire(const CellID& cID, const TVector3& hit_start, const TVector3& hit_end) const;
+ virtual TVector3 wirePos_vs_z(const CellID& cID, const double& zpos) const;
// double phi(const CellID& cID) const;
inline double cell_Size() const { return m_cellSize; }
@@ -117,6 +121,52 @@ public:
// return w;
// }
+ TVector3 LineLineIntersect(TVector3 p1, TVector3 p2, TVector3 p3, TVector3 p4) const {
+ TVector3 p13, p43, p21;
+ double d1343, d4321, d1321, d4343, d2121;
+ double numer, denom;
+ double mua, mub;
+ TVector3 pa, pb;
+
+ p13.SetX(p1.X() - p3.X());
+ p13.SetY(p1.Y() - p3.Y());
+ p13.SetZ(p1.Z() - p3.Z());
+ p43.SetX(p4.X() - p3.X());
+ p43.SetY(p4.Y() - p3.Y());
+ p43.SetZ(p4.Z() - p3.Z());
+ /* if (ABS(p43.X()) < EPS && ABS(p43.Y()) < EPS && ABS(p43.Z()) < EPS) */
+ /* return(FALSE); */
+ p21.SetX(p2.X() - p1.X());
+ p21.SetY(p2.Y() - p1.Y());
+ p21.SetZ(p2.Z() - p1.Z());
+ /* if (ABS(p21.X()) < EPS && ABS(p21.Y()) < EPS && ABS(p21.Z()) < EPS) */
+ /* return(FALSE); */
+
+ d1343 = p13.X() * p43.X() + p13.Y() * p43.Y() + p13.Z() * p43.Z();
+ d4321 = p43.X() * p21.X() + p43.Y() * p21.Y() + p43.Z() * p21.Z();
+ d1321 = p13.X() * p21.X() + p13.Y() * p21.Y() + p13.Z() * p21.Z();
+ d4343 = p43.X() * p43.X() + p43.Y() * p43.Y() + p43.Z() * p43.Z();
+ d2121 = p21.X() * p21.X() + p21.Y() * p21.Y() + p21.Z() * p21.Z();
+
+ denom = d2121 * d4343 - d4321 * d4321;
+ /* if (ABS(denom) < EPS) */
+ /* return(FALSE); */
+ numer = d1343 * d4321 - d1321 * d4343;
+
+ mua = numer / denom;
+ mub = (d1343 + d4321 * (mua)) / d4343;
+
+ pa.SetX(p1.X() + mua * p21.X());
+ pa.SetY(p1.Y() + mua * p21.Y());
+ pa.SetZ(p1.Z() + mua * p21.Z());
+ pb.SetX(p3.X() + mub * p43.X());
+ pb.SetY(p3.Y() + mub * p43.Y());
+ pb.SetZ(p3.Z() + mub * p43.Z());
+
+ return pb - pa;
+ }
+
+
void updateParams(int chamber, int layer) const{
auto it_end = layer_params.cend();
--it_end;
diff --git a/Detector/DetSegmentation/src/GridDriftChamber.cpp b/Detector/DetSegmentation/src/GridDriftChamber.cpp
index ec2e0707..66e9144d 100644
--- a/Detector/DetSegmentation/src/GridDriftChamber.cpp
+++ b/Detector/DetSegmentation/src/GridDriftChamber.cpp
@@ -109,8 +109,8 @@ double GridDriftChamber::distanceTrackWire(const CellID& cID, const TVector3& hi
TVector3 Wend = {0,0,0};
cellposition(cID,Wstart,Wend);
- TVector3 a = hit_end - hit_start;
- TVector3 b = Wend - Wstart;
+ TVector3 a = (hit_end - hit_start).Unit();
+ TVector3 b = (Wend - Wstart).Unit();
TVector3 c = Wstart - hit_start;
double num = std::abs(c.Dot(a.Cross(b)));
@@ -118,13 +118,95 @@ double GridDriftChamber::distanceTrackWire(const CellID& cID, const TVector3& hi
double DCA = 0;
- if (denum) {
- DCA = num / denum;
+ if (denum) {
+ DCA = num / denum;
}
return DCA;
}
+TVector3 GridDriftChamber::distanceClosestApproach(const CellID& cID, const TVector3& hitPos) const {
+ // Distance of the closest approach between a single hit (point) and the closest wire
+
+ TVector3 Wstart = {0,0,0};
+ TVector3 Wend = {0,0,0};
+ cellposition(cID,Wstart,Wend);
+
+ TVector3 temp = (Wend + Wstart);
+ TVector3 Wmid(temp.X() / 2.0, temp.Y() / 2.0, temp.Z() / 2.0);
+
+ double hitPhi = hitPos.Phi();
+ if (hitPhi < 0) {
+ hitPhi = hitPhi + 2 * M_PI;
+ }
+
+ TVector3 PCA = Wstart + ((Wend - Wstart).Unit()).Dot((hitPos - Wstart)) * ((Wend - Wstart).Unit());
+ TVector3 dca = hitPos - PCA;
+
+ return dca;
+}
+
+TVector3 GridDriftChamber::Line_TrackWire(const CellID& cID, const TVector3& hit_start, const TVector3& hit_end) const {
+ // The line connecting a particle track to the closest wire
+ // Returns the vector connecting the both
+ TVector3 Wstart = {0,0,0};
+ TVector3 Wend = {0,0,0};
+ cellposition(cID,Wstart,Wend);
+
+ TVector3 P1 = hit_start;
+ TVector3 P2 = hit_end;
+ TVector3 P3 = Wstart;
+ TVector3 P4 = Wend;
+
+ TVector3 intersect = LineLineIntersect(P1, P2, P3, P4);
+ return intersect;
+}
+
+// Get the wire position for a z
+TVector3 GridDriftChamber::wirePos_vs_z(const CellID& cID, const double& zpos) const {
+
+ TVector3 Wstart = {0,0,0};
+ TVector3 Wend = {0,0,0};
+ cellposition(cID,Wstart,Wend);
+
+ double t = (zpos - Wstart.Z())/(Wend.Z()-Wstart.Z());
+ double x = Wstart.X()+t*(Wend.X()-Wstart.X());
+ double y = Wstart.Y()+t*(Wend.Y()-Wstart.Y());
+
+ TVector3 wireCoord(x, y, zpos);
+ return wireCoord;
+}
+
+TVector3 GridDriftChamber::IntersectionTrackWire(const CellID& cID, const TVector3& hit_start, const TVector3& hit_end) const {
+ // Intersection between the particle track and the wire assuming that the track between hit_start and hit_end is linear
+
+ TVector3 Wstart = {0,0,0};
+ TVector3 Wend = {0,0,0};
+ cellposition(cID,Wstart,Wend);
+
+ TVector3 P1 = hit_start;
+ TVector3 V1 = hit_end-hit_start;
+
+ TVector3 P2 = Wstart;
+ TVector3 V2 = Wend - Wstart;
+
+ TVector3 denom = V1.Cross(V2);
+ double mag_denom = denom.Mag();
+
+ TVector3 intersect(0, 0, 0);
+
+ if (mag_denom !=0)
+ {
+ TVector3 num = ((P2-P1)).Cross(V2);
+ double mag_num = num.Mag();
+ double a = mag_num / mag_denom;
+
+ intersect = P1 + a * V1;
+
+ }
+ return intersect;
+}
+
}
}
diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
index 9389a525..6ebafefc 100644
--- a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
+++ b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
@@ -20,6 +20,8 @@
#include <cmath>
#include <algorithm>
+#include "time.h"
+
DECLARE_COMPONENT( DCHDigiAlg )
DCHDigiAlg::DCHDigiAlg(const std::string& name, ISvcLocator* svcLoc)
@@ -58,8 +60,10 @@ StatusCode DCHDigiAlg::initialize()
else {
m_tuple = ntupleSvc()->book( "MyTuples/DCH_digi_evt", CLID_ColumnWiseTuple, "DCH_digi_evt" );
if ( m_tuple ) {
- m_tuple->addItem( "N_sim" , m_n_sim , 0, 100000 ).ignore();
- m_tuple->addItem( "N_digi", m_n_digi, 0, 100000 ).ignore();
+ m_tuple->addItem( "N_evt" , m_evt ).ignore();
+ m_tuple->addItem( "N_sim" , m_n_sim , 0, 1000000 ).ignore();
+ m_tuple->addItem( "N_digi", m_n_digi, 0, 1000000 ).ignore();
+ m_tuple->addItem( "run_time" , m_time ).ignore();
m_tuple->addItem( "simhit_x", m_n_sim, m_simhit_x).ignore();
m_tuple->addItem( "simhit_y", m_n_sim, m_simhit_y).ignore();
m_tuple->addItem( "simhit_z", m_n_sim, m_simhit_z).ignore();
@@ -86,7 +90,10 @@ StatusCode DCHDigiAlg::initialize()
StatusCode DCHDigiAlg::execute()
{
+ m_start = clock();
+
info() << "Processing " << _nEvt << " events " << endmsg;
+ m_evt = _nEvt;
std::map<unsigned long long, std::vector<edm4hep::SimTrackerHit> > id_hits_map;
edm4hep::TrackerHitCollection* Vec = w_DigiDCHCol.createAndPut();
edm4hep::MCRecoTrackerAssociationCollection* AssoVec = w_AssociationCol.createAndPut();
@@ -135,20 +142,36 @@ StatusCode DCHDigiAlg::execute()
m_segmentation->cellposition(wcellid, Wstart, Wend);
float dd4hep_mm = dd4hep::mm;
//std::cout<<"dd4hep_mm="<<dd4hep_mm<<std::endl;
- Wstart =(1/dd4hep_mm)* Wstart;// from DD4HEP cm to mm
- Wend =(1/dd4hep_mm)* Wend ;
+// Wstart =(1/dd4hep_mm)* Wstart;// from DD4HEP cm to mm
+// Wend =(1/dd4hep_mm)* Wend ;
//std::cout<<"wcellid="<<wcellid<<",chamber="<<chamber<<",layer="<<layer<<",cellID="<<cellID<<",s_x="<<Wstart.x()<<",s_y="<<Wstart.y()<<",s_z="<<Wstart.z()<<",E_x="<<Wend.x()<<",E_y="<<Wend.y()<<",E_z="<<Wend.z()<<std::endl;
TVector3 denominator = (Wend-Wstart) ;
float min_distance = 999 ;
+ float min_line_distance = 999 ;
+ float tmp_distance =0;
for(unsigned int i=0; i< simhit_size; i++)
{
float sim_hit_mom = sqrt( iter->second.at(i).getMomentum()[0]*iter->second.at(i).getMomentum()[0] + iter->second.at(i).getMomentum()[1]*iter->second.at(i).getMomentum()[1] + iter->second.at(i).getMomentum()[2]*iter->second.at(i).getMomentum()[2] );//GeV
float sim_hit_pt = sqrt( iter->second.at(i).getMomentum()[0]*iter->second.at(i).getMomentum()[0] + iter->second.at(i).getMomentum()[1]*iter->second.at(i).getMomentum()[1] );//GeV
- TVector3 pos(iter->second.at(i).getPosition()[0], iter->second.at(i).getPosition()[1], iter->second.at(i).getPosition()[2]);
- TVector3 numerator = denominator.Cross(Wstart-pos) ;
- float tmp_distance = numerator.Mag()/denominator.Mag() ;
- //std::cout<<"tmp_distance="<<tmp_distance<<",x="<<pos.x()<<",y="<<pos.y()<<",z="<<pos.z()<<",mom="<<sim_hit_mom<<",pt="<<sim_hit_pt<<std::endl;
+ TVector3 pos(iter->second.at(i).getPosition()[0]*dd4hep_mm, iter->second.at(i).getPosition()[1]*dd4hep_mm, iter->second.at(i).getPosition()[2]*dd4hep_mm);
+
+// TVector3 numerator = denominator.Cross(Wstart-pos) ;
+// float tmp_distance = numerator.Mag()/denominator.Mag() ;
+
+ TVector3 sim_mon(iter->second.at(i).getMomentum()[0],iter->second.at(i).getMomentum()[1],iter->second.at(i).getMomentum()[2]);
+ float Steplength = iter->second.at(i).getPathLength();
+ TVector3 pos_start = pos - 0.5 * Steplength * sim_mon.Unit();
+ TVector3 pos_end = pos + 0.5 * Steplength * sim_mon.Unit();
+ if(m_Doca) {
+ tmp_distance = m_segmentation->distanceTrackWire(wcellid,pos_start,pos_end);
+ } else {
+ tmp_distance = (m_segmentation->distanceClosestApproach(wcellid,pos)).Mag();
+ }
+
+
+ // std::cout << " Steplength= " << Steplength << std::endl;
+ // std::cout<<"tmp_distance="<<tmp_distance<<",x="<<pos.x()<<",y="<<pos.y()<<",z="<<pos.z()<<",mom="<<sim_hit_mom<<",pt="<<sim_hit_pt<<std::endl;
if(tmp_distance < min_distance){
min_distance = tmp_distance;
@@ -157,13 +180,12 @@ StatusCode DCHDigiAlg::execute()
pos_z = pos.z();
}
tot_length += iter->second.at(i).getPathLength();//mm
-
auto asso = AssoVec->create();
asso.setRec(trkHit);
asso.setSim(iter->second.at(i));
asso.setWeight(iter->second.at(i).getEDep()/tot_edep);
- if(m_WriteAna && (nullptr!=m_tuple)){
+ if(m_WriteAna && (nullptr!=m_tuple)) { // && min_distance <0.3){
m_simhit_x[m_n_sim] = pos.x();
m_simhit_y[m_n_sim] = pos.y();
m_simhit_z[m_n_sim] = pos.z();
@@ -177,7 +199,7 @@ StatusCode DCHDigiAlg::execute()
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
- if(m_WriteAna && (nullptr!=m_tuple)){
+ if(m_WriteAna && (nullptr!=m_tuple)) { // && min_distance <0.3){
m_chamber [m_n_digi] = chamber;
m_layer [m_n_digi] = layer ;
m_cell [m_n_digi] = cellID;
@@ -204,6 +226,9 @@ StatusCode DCHDigiAlg::execute()
return StatusCode::FAILURE;
}
}
+ m_end = clock();
+ m_time = (m_end - m_start);
+
return StatusCode::SUCCESS;
}
diff --git a/Digitisers/DCHDigi/src/DCHDigiAlg.h b/Digitisers/DCHDigi/src/DCHDigiAlg.h
index 5ddf59a5..178f5518 100644
--- a/Digitisers/DCHDigi/src/DCHDigiAlg.h
+++ b/Digitisers/DCHDigi/src/DCHDigiAlg.h
@@ -45,8 +45,10 @@ protected:
int _nEvt ;
NTuple::Tuple* m_tuple = nullptr ;
+ NTuple::Item<int> m_evt;
NTuple::Item<long> m_n_sim;
NTuple::Item<long> m_n_digi;
+ NTuple::Item<float> m_time;
NTuple::Array<int > m_chamber ;
NTuple::Array<int > m_layer ;
NTuple::Array<int > m_cell ;
@@ -62,7 +64,7 @@ protected:
NTuple::Array<float> m_hit_dE ;
NTuple::Array<float> m_hit_dE_dx ;
-
+ clock_t m_start,m_end;
dd4hep::rec::CellIDPositionConverter* m_cellIDConverter;
dd4hep::DDSegmentation::GridDriftChamber* m_segmentation;
@@ -76,6 +78,7 @@ protected:
Gaudi::Property<float> m_velocity { this, "drift_velocity", 40};// um/ns
Gaudi::Property<float> m_mom_threshold { this, "mom_threshold", 0};// GeV
Gaudi::Property<bool> m_WriteAna { this, "WriteAna", false};
+ Gaudi::Property<bool> m_Doca { this, "Doca", false};//1:line dca 0:point dca
// Input collections
--
GitLab