diff --git a/Detector/DetDriftChamber/compact/det.xml b/Detector/DetDriftChamber/compact/det.xml
index 4e981b021e24069755937850709a346dbf703a29..b62dd532daa4c15c48ae37228faeffe3fe32692f 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>
@@ -113,7 +113,6 @@
<segmentation type="GridDriftChamber" cell_size="SDT_chamber_cell_width" epsilon0="Epsilon" detector_length="DC_length" identifier_phi="cellID" DC_rbegin="DC_chamber_layer_rbegin" DC_rend="DC_chamber_layer_rend" DC_rmin="SDT_chamber_radius_min" DC_rmax="SDT_chamber_radius_max" safe_distance="DC_safe_distance" layerID="layer" layer_width="SDT_chamber_layer_width"/>
- <!-- <id>system:8,chamber:1,layer:8,cellID:16</id> -->
<id>system:5,layer:7:9,chamber:8,cellID:32:16</id>
</readout>
</readouts>
diff --git a/Detector/DetDriftChamber/src/driftchamber/DriftChamber.cpp b/Detector/DetDriftChamber/src/driftchamber/DriftChamber.cpp
index bb02ad39e3d2347af6a5a9ee76ee51593227d3d7..be98c1c14fa100b5fd5a4199450fde98e9aebf72 100644
--- a/Detector/DetDriftChamber/src/driftchamber/DriftChamber.cpp
+++ b/Detector/DetDriftChamber/src/driftchamber/DriftChamber.cpp
@@ -59,6 +59,7 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector,
double chamber_layer_rend = theDetector.constant<double>("DC_chamber_layer_rend");
int chamber_layer_number = floor((chamber_layer_rend-chamber_layer_rbegin)/chamber_layer_width);
+ double safe_distance = theDetector.constant<double>("DC_safe_distance");
double epsilon = theDetector.constant<double>("Epsilon");
// =======================================================================
@@ -83,12 +84,6 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector,
// - chamber volume
dd4hep::Tube det_chamber_solid(chamber_radius_min, chamber_radius_max, chamber_half_length);
dd4hep::Volume det_chamber_vol(det_name+"_chamber_vol", det_chamber_solid, chamber_mat);
- if ( x_det.isSensitive() ) {
- det_chamber_vol.setRegion(theDetector,x_det.regionStr());
- det_chamber_vol.setLimitSet(theDetector,x_det.limitsStr());
- det_chamber_vol.setSensitiveDetector(sens);
- sd.setType("tracker");
- }
// - wall
double chamber_inner_wall_rmin = theDetector.constant<double>("SDT_chamber_inner_wall_radius_min");
@@ -177,6 +172,18 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector,
DCHseg->setGeomParams(chamber_id, layerIndex, layer_Phi, rmid, epsilon, offset);
DCHseg->setWiresInLayer(chamber_id, layerIndex, numWire);
+
+ dd4hep::Tube layer_vol_solid(rmin,rmax,chamber_half_length);
+ dd4hep::Volume layer_vol(det_name+"_layer_vol",layer_vol_solid,det_mat);
+ current_vol_ptr = &layer_vol;
+
+ if ( x_det.isSensitive() ) {
+ layer_vol.setRegion(theDetector,x_det.regionStr());
+ layer_vol.setLimitSet(theDetector,x_det.limitsStr());
+ layer_vol.setSensitiveDetector(sens);
+ sd.setType("tracker");
+ }
+
// - wire vol
//phi <-------------------> -phi
// | F8 F7 F6 F5| Only on the outermost layer.
@@ -185,7 +192,6 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector,
// | |
// | F0 F1 F2 F3|
// -----------------------
-// if(layer_id == 0 || layer_id == 1 || layer_id == 2 || layer_id == 99) {
for(int icell=0; icell< numWire; icell++) {
double wire_phi = (icell+0.5)*layer_Phi + offset;
// - signal wire
@@ -193,7 +199,7 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector,
dd4hep::PlacedVolume module_phy = (*current_vol_ptr).placeVolume(module_vol,transform_module);
// - Field wire
dd4hep::PlacedVolume Module_phy;
- double radius[9] = {rmid-chamber_layer_width*0.5,rmid-chamber_layer_width*0.5,rmid-chamber_layer_width*0.5,rmid-chamber_layer_width*0.5,rmid,rmid+chamber_layer_width*0.5,rmid+chamber_layer_width*0.5,rmid+chamber_layer_width*0.5,rmid+chamber_layer_width*0.5};
+ double radius[9] = {rmid-chamber_layer_width*0.5+safe_distance,rmid-chamber_layer_width*0.5+safe_distance,rmid-chamber_layer_width*0.5+safe_distance,rmid-chamber_layer_width*0.5+safe_distance,rmid,rmid+chamber_layer_width*0.5-safe_distance,rmid+chamber_layer_width*0.5-safe_distance,rmid+chamber_layer_width*0.5-safe_distance,rmid+chamber_layer_width*0.5-safe_distance};
double phi[9] = {wire_phi+layer_Phi*0.25,wire_phi,wire_phi-layer_Phi*0.25,wire_phi-layer_Phi*0.5,wire_phi-layer_Phi*0.5,wire_phi-layer_Phi*0.5,wire_phi-layer_Phi*0.25,wire_phi,wire_phi+layer_Phi*0.25};
int num = 5;
if(layer_id==(chamber_layer_number-1)) {
@@ -206,9 +212,11 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector,
Module_phy = (*current_vol_ptr).placeVolume(Module_vol,transform_Module);
}
}
- }
-
-// }
+ dd4hep::Transform3D transform_layer(dd4hep::Rotation3D(),
+ dd4hep::Position(0,0,0));
+ dd4hep::PlacedVolume layer_phy = det_chamber_vol.placeVolume(layer_vol ,transform_layer);
+ layer_phy.addPhysVolID("layer", layer_id);
+ }
// - place in det
// - chamber
diff --git a/Detector/DetSegmentation/include/DetSegmentation/GridDriftChamber.h b/Detector/DetSegmentation/include/DetSegmentation/GridDriftChamber.h
index 7c1a0864b143b4b13080933667af169b5aa6d7e3..0233ccf829d16144cdc96513b85177d3bb4c2d33 100644
--- a/Detector/DetSegmentation/include/DetSegmentation/GridDriftChamber.h
+++ b/Detector/DetSegmentation/include/DetSegmentation/GridDriftChamber.h
@@ -61,6 +61,11 @@ 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;
+ TVector3 LineLineIntersect(TVector3& p1, TVector3& p2, TVector3& p3, TVector3& p4) 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; }
@@ -109,14 +114,6 @@ public:
inline auto returnAllWires() const { return m_wiresPositions; }
-// inline TVector3 returnWirePosition(double angle, int sign) const {
-// TVector3 w(0, 0, 0);
-// w.SetX(_currentRadius * std::cos(angle));
-// w.SetY(_currentRadius * std::sin(angle));
-// w.SetZ(sign * m_detectorLength / 2.0);
-// return w;
-// }
-
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 ec2e07070a42467cc91a5e8806047f71bb09b9d5..a4302ebf59b180e9bf2abdc83fe2e3bb9e7ba9db 100644
--- a/Detector/DetSegmentation/src/GridDriftChamber.cpp
+++ b/Detector/DetSegmentation/src/GridDriftChamber.cpp
@@ -102,6 +102,52 @@ void GridDriftChamber::cellposition(const CellID& cID, TVector3& Wstart,
Wend = returnWirePosition(phi_end, 1);
}
+TVector3 GridDriftChamber::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;
+}
+
double GridDriftChamber::distanceTrackWire(const CellID& cID, const TVector3& hit_start,
const TVector3& hit_end) const {
@@ -109,8 +155,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 +164,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 6f88f5898faf2ce84fb3ae4d1f4410fcb521a0bb..01fb9406a6811e673a8f1d5e53c9c43bcc6893ea 100644
--- a/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
+++ b/Digitisers/DCHDigi/src/DCHDigiAlg.cpp
@@ -20,21 +20,23 @@
#include <cmath>
#include <algorithm>
+#include "time.h"
+
DECLARE_COMPONENT( DCHDigiAlg )
DCHDigiAlg::DCHDigiAlg(const std::string& name, ISvcLocator* svcLoc)
: GaudiAlgorithm(name, svcLoc),
_nEvt(0)
{
-
+
// Input collections
declareProperty("SimDCHitCollection", r_SimDCHCol, "Handle of the Input SimHit collection");
-
+
// Output collections
declareProperty("DigiDCHitCollection", w_DigiDCHCol, "Handle of Digi DCHit collection");
-
+
declareProperty("AssociationCollection", w_AssociationCol, "Handle of Association collection");
-
+
}
StatusCode DCHDigiAlg::initialize()
@@ -52,14 +54,15 @@ StatusCode DCHDigiAlg::initialize()
}
if(m_WriteAna){
-
NTuplePtr nt( ntupleSvc(), "MyTuples/DCH_digi_evt" );
if ( nt ) m_tuple = nt;
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 +89,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();
@@ -98,8 +104,8 @@ StatusCode DCHDigiAlg::execute()
unsigned long long id = SimHit.getCellID();
float sim_hit_mom = sqrt( SimHit.getMomentum()[0]*SimHit.getMomentum()[0] + SimHit.getMomentum()[1]*SimHit.getMomentum()[1] + SimHit.getMomentum()[2]*SimHit.getMomentum()[2] );//GeV
if(sim_hit_mom < m_mom_threshold) continue;
- if(SimHit.getEDep() <= 0) continue;
-
+ if(SimHit.getEDep() <= 0) continue;
+
if ( id_hits_map.find(id) != id_hits_map.end()) id_hits_map[id].push_back(SimHit);
else
{
@@ -135,20 +141,39 @@ 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);
+ tmp_distance = tmp_distance/dd4hep_mm; //mm
+ } else {
+ tmp_distance = (m_segmentation->distanceClosestApproach(wcellid,pos)).Mag();
+ tmp_distance = tmp_distance/dd4hep_mm; //mm
+ }
+
+
+ // 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;
pos_x = pos.x();
@@ -156,27 +181,26 @@ 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();
m_n_sim ++ ;
}
}
-
+
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/tot_length); // GeV/mm
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;
@@ -203,6 +227,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 5ddf59a5426a86358f0b2d775cef8946e46ddbc1..d6a82961f2c70c43a0dc524511d9618670f28553 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,20 +64,21 @@ 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;
dd4hep::DDSegmentation::BitFieldCoder* m_decoder;
-
+
Gaudi::Property<std::string> m_readout_name{ this, "readout", "DriftChamberHitsCollection"};//readout for getting segmentation
-
+
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
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
@@ -83,6 +86,6 @@ protected:
// Output collections
DataHandle<edm4hep::TrackerHitCollection> w_DigiDCHCol{"DigiDCHitCollection", Gaudi::DataHandle::Writer, this};
DataHandle<edm4hep::MCRecoTrackerAssociationCollection> w_AssociationCol{"DCHitAssociationCollection", Gaudi::DataHandle::Writer, this};
-};
+};
#endif
diff --git a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp
index cef756adc1d6cc624bf8519d69b5e367423413d9..0e93b89c6c02fcb03f2ffb3b14a484169c347675 100644
--- a/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp
+++ b/Reconstruction/RecGenfitAlg/src/RecGenfitAlgDC.cpp
@@ -382,20 +382,21 @@ void RecGenfitAlgDC::debugEvent()
mcParticleCol=m_mcParticleCol.get();
simDCHitCol=m_simDCHitCol.get();
m_nSimDCHit=simDCHitCol->size();
- int iMcParticle=0;
int iHit=0;
+ for(auto simDCHit: *simDCHitCol){
+ edm4hep::Vector3d pos=simDCHit.position();
+ TVectorD p(3);
+ p[0]=pos.x;//no unit conversion here
+ p[1]=pos.y;
+ p[2]=pos.z;
+ m_mdcHitMcX[iHit]=pos.x;
+ m_mdcHitMcY[iHit]=pos.y;
+ m_mdcHitMcZ[iHit]=pos.z;
+ iHit++;
+ }
+ m_mcIndex=iHit;
+ int iMcParticle=0;
for(auto mcParticle : *mcParticleCol){
- for(auto simDCHit: *simDCHitCol){
- edm4hep::Vector3d pos=simDCHit.position();
- TVectorD p(3);
- p[0]=pos.x;//no unit conversion here
- p[1]=pos.y;
- p[2]=pos.z;
- m_mdcHitMcX[iHit]=pos.x;
- m_mdcHitMcY[iHit]=pos.y;
- m_mdcHitMcZ[iHit]=pos.z;
- iHit++;
- }
edm4hep::Vector3f mcPocaMom = mcParticle.getMomentum();//GeV
float px=mcPocaMom.x;
float py=mcPocaMom.y;
@@ -404,6 +405,5 @@ void RecGenfitAlgDC::debugEvent()
m_pocaMomMcP[iMcParticle]=sqrt(px*px+py*py+pz*pz);
iMcParticle++;
}
- m_mcIndex=iHit;
}