diff --git a/Detector/DetDriftChamber/CMakeLists.txt b/Detector/DetDriftChamber/CMakeLists.txt index b680276955849c32df5af4e92216e9cab21accf9..0c1cc85b2d26105a245427a369270ad4027769a5 100644 --- a/Detector/DetDriftChamber/CMakeLists.txt +++ b/Detector/DetDriftChamber/CMakeLists.txt @@ -27,9 +27,9 @@ gaudi_add_module(DetDriftChamber ${DetDriftChamber_src} INCLUDE_DIRS # DD4hep ROOT Geant4 src/include - LINK_LIBRARIES + LINK_LIBRARIES GaudiKernel DD4hep ROOT DetSegmentation # GaudiKernel - DD4hep ${DD4hep_COMPONENT_LIBRARIES} + #DD4hep ${DD4hep_COMPONENT_LIBRARIES} # ROOT Geant4 ) diff --git a/Detector/DetDriftChamber/compact/det.xml b/Detector/DetDriftChamber/compact/det.xml index 84a1ec7b05d00220ac3fda86c50f4d0111b550b7..b791565dba6598a8c2edc8a178b30afd2c30e25a 100644 --- a/Detector/DetDriftChamber/compact/det.xml +++ b/Detector/DetDriftChamber/compact/det.xml @@ -38,6 +38,7 @@ <constant name="SDT_inner_chamber_layer_number" value="67"/> <constant name="SDT_outer_chamber_layer_number" value="63"/> <constant name="SDT_chamber_layer_width" value="10*mm"/> + <constant name="Epsilon" value="0*deg"/> </define> @@ -52,14 +53,15 @@ <detectors> <detector id="1" name="DriftChamber" type="DriftChamber" readout="DriftChamberHitsCollection" vis="VisibleBlue" sensitive="true"> <!-- Use cm as unit if you want to use Pandora for reconstruction --> + <sensitive type="SimpleDriftChamber"/> </detector> </detectors> <readouts> <readout name="DriftChamberHitsCollection"> - <segmentation type="GridDriftChamber" delta_phi="8*deg" identifier_phi="phi"/> + <segmentation type="GridDriftChamber" cell_size="10*mm" offset_phi="0." epsilon0="Epsilon" detector_length="SDT_half_length" identifier_phi="cellID" /> - <id>system:8,chamber:1,layer:8,phi:16</id> + <id>system:8,chamber:1,layer:8,cellID:16</id> </readout> </readouts> diff --git a/Detector/DetDriftChamber/src/driftchamber/DriftChamber.cpp b/Detector/DetDriftChamber/src/driftchamber/DriftChamber.cpp index c71012156662c8b053567bf2667d535c097016a6..07b95b983df469f822a61b22d7ba2827e27b50eb 100644 --- a/Detector/DetDriftChamber/src/driftchamber/DriftChamber.cpp +++ b/Detector/DetDriftChamber/src/driftchamber/DriftChamber.cpp @@ -8,8 +8,10 @@ //==================================================================== #include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/Printout.h" #include "XML/Layering.h" #include "XML/Utilities.h" +#include "XML/XMLElements.h" #include "DDRec/DetectorData.h" #include "DDSegmentation/Segmentation.h" #include "DetSegmentation/GridDriftChamber.h" @@ -24,6 +26,9 @@ using namespace dd4hep::rec ; static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector, xml_h e, dd4hep::SensitiveDetector sens) { + // ------- Lambda functions ---- // + auto delta_a_func = [](auto x, auto y) { return 0.5 * ( x + y ); }; + // ======================================================================= // Parameter Definition // ======================================================================= @@ -33,6 +38,8 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector, std::string det_name = x_det.nameStr(); std::string det_type = x_det.typeStr(); + dd4hep::SensitiveDetector sd = sens; + // - global double chamber_radius_min = theDetector.constant<double>("SDT_radius_min"); double chamber_radius_max = theDetector.constant<double>("SDT_radius_max"); @@ -53,6 +60,8 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector, int outer_chamber_layer_number = theDetector.constant<int>("SDT_outer_chamber_layer_number"); double chamber_layer_width = theDetector.constant<double>("SDT_chamber_layer_width"); + double epsilon = theDetector.constant<double>("Epsilon"); + // ======================================================================= // Detector Construction // ======================================================================= @@ -75,36 +84,55 @@ static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector, dd4hep::Tube det_outer_chamber_solid(outer_chamber_radius_min, outer_chamber_radius_max, outer_chamber_length*0.5); dd4hep::Volume det_outer_chamber_vol(det_name+"_outer_chamber_vol", det_outer_chamber_solid, det_mat); + //Initialize the segmentation + dd4hep::Readout readout = sd.readout(); + dd4hep::Segmentation geomseg = readout.segmentation(); + dd4hep::Segmentation* _geoSeg = &geomseg; + + auto DCHseg = dynamic_cast<dd4hep::DDSegmentation::GridDriftChamber*>(_geoSeg->segmentation()); // - layer - for(int layer_id=0; layer_id<(inner_chamber_layer_number+outer_chamber_layer_number-1);layer_id++) { - double rmin,rmax; - std::string layer_name; - dd4hep::Volume* current_vol_ptr = nullptr; - - if(layer_id<inner_chamber_layer_number){ - current_vol_ptr = &det_inner_chamber_vol; - rmin = inner_chamber_radius_min+(layer_id*chamber_layer_width); - rmax = rmin+chamber_layer_width; - layer_name = det_name+"_inner_chamber_vol"+_toString(layer_id,"_layer%d"); - }else{ - current_vol_ptr = &det_outer_chamber_vol; - rmin = outer_chamber_radius_min+((layer_id-inner_chamber_layer_number)*chamber_layer_width); - rmax = rmin+chamber_layer_width; - layer_name = det_name+"_outer_chamber_vol"+_toString(layer_id,"_layer%d"); + for(int layer_id = 0; layer_id < (inner_chamber_layer_number+outer_chamber_layer_number); layer_id++) { + double rmin,rmax,offset; + std::string layer_name; + dd4hep::Volume* current_vol_ptr = nullptr; + + if( layer_id < inner_chamber_layer_number ) { + current_vol_ptr = &det_inner_chamber_vol; + rmin = inner_chamber_radius_min+(layer_id*chamber_layer_width); + rmax = rmin+chamber_layer_width; + layer_name = det_name+"_inner_chamber_vol"+_toString(layer_id,"_layer%d"); + } + else { + current_vol_ptr = &det_outer_chamber_vol; + rmin = outer_chamber_radius_min+((layer_id-inner_chamber_layer_number)*chamber_layer_width); + rmax = rmin+chamber_layer_width; + layer_name = det_name+"_outer_chamber_vol"+_toString(layer_id,"_layer%d"); + } + + //Construction of drift chamber layers + double rmid = delta_a_func(rmin,rmax); + double ilayer_cir = 2 * M_PI * rmid; + double ncell = ilayer_cir / chamber_layer_width; + int ncell_layer = ceil(ncell); + int numWire = ncell_layer; + double layer_Phi = 2*M_PI / ncell_layer; + if(layer_id %2 ==0){ offset = 0.; } + else { offset = 0.5 * layer_Phi; } + + DCHseg->setGeomParams(layer_id, layer_Phi, rmid, epsilon, offset); + DCHseg->setWiresInLayer(layer_id, numWire); + + dd4hep::Tube layer_solid(rmin,rmax,chamber_length*0.5); + dd4hep::Volume layer_vol(layer_name,layer_solid,det_mat); + dd4hep::Transform3D transform_layer(dd4hep::Rotation3D(),dd4hep::Position(0.,0.,0.)); + dd4hep::PlacedVolume layer_phy = (*current_vol_ptr).placeVolume(layer_vol, transform_layer); + layer_phy.addPhysVolID("layer",layer_id); + + //Set drift chamber layers to sensitive detector + layer_vol.setSensitiveDetector(sens); + sd.setType("tracker"); } - /// Construction of drift chamber layers - dd4hep::Tube layer_solid(rmin,rmax,chamber_length*0.5); - dd4hep::Volume layer_vol(layer_name,layer_solid,det_mat); - dd4hep::Transform3D transform_layer(dd4hep::Rotation3D(),dd4hep::Position(0.,0.,0.)); - dd4hep::PlacedVolume layer_phy = (*current_vol_ptr).placeVolume(layer_vol, transform_layer); - layer_phy.addPhysVolID("layer",layer_id); - - /// Set drift chamber layers to sensitive detector - dd4hep::SensitiveDetector sd = sens; - layer_vol.setSensitiveDetector(sens); - sd.setType("tracker"); - } // - place in det // inner diff --git a/Detector/DetSegmentation/DetSegmentation/GridDriftChamber.h b/Detector/DetSegmentation/DetSegmentation/GridDriftChamber.h index 379ad4a19fdbce33a63f009577970534f3f995c3..fc01bf114d8d1780690ecbc9db2c23a82cc8cbe5 100644 --- a/Detector/DetSegmentation/DetSegmentation/GridDriftChamber.h +++ b/Detector/DetSegmentation/DetSegmentation/GridDriftChamber.h @@ -6,6 +6,7 @@ #include "TVector3.h" #include <cmath> #include <iostream> +#include <map> /** GridDriftChamber Detector/DetSegmentation/DetSegmentation/GridDriftChamber.h GridDriftChamber.h * @@ -14,6 +15,21 @@ * @author nalipour */ + +typedef struct Layer + { + double layerphi; + double R; + double eps; + double offset; + Layer(){}; + Layer(double x, double y, double z, double k):layerphi(x),R(y),eps(z),offset(k){}; + bool operator < (const Layer &a) const + { + return layerphi < a.layerphi; + } + } LAYER; + namespace dd4hep { namespace DDSegmentation { class GridDriftChamber : public Segmentation { @@ -28,62 +44,102 @@ public: virtual Vector3D position(const CellID& aCellID) const; virtual CellID cellID(const Vector3D& aLocalPosition, const Vector3D& aGlobalPosition, const VolumeID& aVolumeID) const; + virtual double distanceTrackWire(const CellID& cID, const TVector3& hit_start, const TVector3& hit_end) const; -// inline double innerRadius() const { return m_innerRadius; } -// inline double detectorLength() const { return m_detectorLength; } - inline double offsetPhi() const { return m_offsetPhi; } - inline double delta_phi() const{ return m_delta_phi; } + inline double cell_Size() const { return m_cellSize; } + inline double offset_phi() const { return m_offsetPhi; } + inline double epsilon0() const { return m_epsilon0; } + inline double detectorLength() const { return m_detectorLength; } inline const std::string& fieldNamePhi() const { return m_phiID; } // Setters -// inline void setGeomParams(int layer, double sizePhi) { -// layer_params[layer] = {sizePhi}; -// } - -// void updateParams(int layer) const { -// auto it_end = layer_params.cend(); -// --it_end; -// double size = it_end->second[0]; -// double radius = it_end->second[1]; -// double eps = it_end->second[2]; -// -// auto map_it = layer_params.find(layer); -// if (map_it != layer_params.cend()) { -// size = map_it->second[0]; -// radius = map_it->second[1]; -// eps = map_it->second[2]; -// } -// -// _currentGridSizePhi = size; -// _currentRadius = radius; -// m_epsilon = eps; -// } - inline double phiFromXY(const Vector3D& aposition) const { return std::atan2(aposition.Y, aposition.X) + M_PI ; } -// inline int returnLayer(double x, double y) const { -// // Hit R position -// double R = std::sqrt(x * x + y * y); -// // Layer -// int layer = int((R - m_innerRadius) / m_cellSize); -// return layer; -// } + inline void setGeomParams(int layer, double layerphi, double R, double eps, double offset ) { + // layer_params[layer] = {layerphi,R,eps}; + layer_params.insert(std::pair<int,LAYER>(layer,LAYER(layerphi,R,eps,offset))); + } + + inline void setWiresInLayer(int layer, int numWires) + { + double phi0; + updateParams(layer); + for (int i = 0; i<numWires; ++i) { + +// if(layer % 2 == 0) { phi0 = 0.; } +// else { phi0 = 0.5 * _currentLayerphi; } + double phi0 = m_offset; + + auto phi_start = _currentLayerphi * i + phi0; + if(phi_start > 2 * M_PI) { phi_start = phi_start - 2 * M_PI; } + auto phi_end = phi_start + _currentLayerphi; + + TVector3 Wstart = returnWirePosition(phi_start, 1); + TVector3 Wend = returnWirePosition(phi_end, -1); + + TVector3 Wmid = (Wstart+Wend)*(1/2.0); + TVector3 Wdirection = (Wend - Wstart); + + m_wiresPositions[layer].push_back(std::make_pair(Wmid, Wdirection)); + } + } + + 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 layer) const{ + auto it_end = layer_params.cend(); + --it_end; + double layerphi = it_end->second.layerphi; + double radius = it_end->second.R; + double eps = it_end->second.eps; + double offset = it_end->second.offset; + + auto map_it = layer_params.find(layer); + if (map_it != layer_params.cend()) { + layerphi = map_it->second.layerphi; + radius = map_it->second.R; + eps = map_it->second.eps; + offset = map_it->second.offset; + } + _currentLayerphi = layerphi; + _currentRadius = radius; + m_epsilon = eps; + m_offset = offset; + } + + inline double returnAlpha() const { + double alpha = 2 * std::asin(m_detectorLength * std::tan(m_epsilon0)/(2 * _currentRadius)); + return alpha; + } protected: /* *** nalipour *** */ double phi(const CellID& cID) const; + std::map<int,LAYER> layer_params; // <layer, {layerphi, R, eps, offset}> + std::map<int, std::vector<std::pair<TVector3, TVector3> >> m_wiresPositions; // < layer, vec<WireMidpoint, WireDirection> > - + double m_cellSize; double m_offsetPhi; - double m_delta_phi; + double m_epsilon0; + double m_detectorLength; std::string m_phiID; // Current parameters of the layer: sizePhi -// mutable double _currentGridSizePhi; // current size Phi -// mutable double _currentRadius; // current size radius -// mutable double m_epsilon; + mutable double _currentLayerphi; + mutable double _currentRadius; + mutable double m_epsilon; + mutable double m_offset; + }; } } diff --git a/Detector/DetSegmentation/src/GridDriftChamber.cpp b/Detector/DetSegmentation/src/GridDriftChamber.cpp index 4b1ed216fb0825e1b41a210871ca2451954ee22e..57a6d7ebacd9fb221dbb7e59dfedd2a14a52894b 100644 --- a/Detector/DetSegmentation/src/GridDriftChamber.cpp +++ b/Detector/DetSegmentation/src/GridDriftChamber.cpp @@ -1,4 +1,5 @@ #include "DetSegmentation/GridDriftChamber.h" +#include <map> namespace dd4hep { namespace DDSegmentation { @@ -9,17 +10,22 @@ GridDriftChamber::GridDriftChamber(const std::string& cellEncoding) : Segmentati _type = "GridDriftChamber"; _description = "Drift chamber segmentation in the global coordinates"; - registerIdentifier("identifier_phi", "Cell ID identifier for phi", m_phiID, "phi"); - registerParameter("delta_phi", "delta phi", m_delta_phi, 0., SegmentationParameter::LengthUnit); + registerParameter("cell_size", "cell size", m_cellSize, 0., SegmentationParameter::LengthUnit); + registerParameter("offset_phi", "offset in phi", m_offsetPhi, 0., SegmentationParameter::LengthUnit, true); + registerParameter("detector_length", "Length of the wire", m_detectorLength, 1., SegmentationParameter::LengthUnit); + registerIdentifier("identifier_phi", "Cell ID identifier for phi", m_phiID, "cellID"); } GridDriftChamber::GridDriftChamber(const BitFieldCoder* decoder) : Segmentation(decoder) { - // define type and description + _type = "GridDriftChamber"; _description = "Drift chamber segmentation in the global coordinates"; - registerIdentifier("identifier_phi", "Cell ID identifier for phi", m_phiID, "phi"); - registerParameter("delta_phi", "delta phi", m_delta_phi, 0., SegmentationParameter::LengthUnit); + registerParameter("cell_size", "cell size", m_cellSize, 1., SegmentationParameter::LengthUnit); + registerParameter("offset_phi", "offset in phi", m_offsetPhi, 0., SegmentationParameter::LengthUnit, true); + registerParameter("epsilon0", "epsilon", m_epsilon0, 0., SegmentationParameter::AngleUnit, true); + registerParameter("detector_length", "Length of the wire", m_detectorLength, 1., SegmentationParameter::LengthUnit); + registerIdentifier("identifier_phi", "Cell ID identifier for phi", m_phiID, "cellID"); } Vector3D GridDriftChamber::position(const CellID& /*cID*/) const { @@ -27,28 +33,78 @@ Vector3D GridDriftChamber::position(const CellID& /*cID*/) const { return cellPosition; } + CellID GridDriftChamber::cellID(const Vector3D& /*localPosition*/, const Vector3D& globalPosition, const VolumeID& vID) const { CellID cID = vID; + unsigned int layerID = _decoder->get(vID, "layer"); + updateParams(layerID); double phi_hit = phiFromXY(globalPosition); double posx = globalPosition.X; double posy = globalPosition.Y; - - int lphi = (int) (phi_hit/m_delta_phi); + double offsetphi= m_offset; + int _lphi; +// if(layerID % 2 == 0) { +// offsetphi = 0.; +// _lphi = (int) (phi_hit / _currentLayerphi); +// } +// else { +// offsetphi = _currentLayerphi / 2.; + if(phi_hit >= offsetphi) { + _lphi = (int) ((phi_hit - offsetphi)/ _currentLayerphi); + } + else { + _lphi = (int) ((phi_hit - offsetphi + 2 * M_PI)/ _currentLayerphi); + } + int lphi = _lphi; _decoder->set(cID, m_phiID, lphi); -// std::cout << " myliu: " -// << " x: " << posx -// << " y: " << posy -//// << " pre: " << phi_pre -// << " phi_hit: " << phi_hit -// << " lphi: " << lphi -// << std::endl; +//std::cout << "#######################################: " +// << " offset : " << m_offset +// << " offsetphi: " << offsetphi +// << " layerID: " << layerID +// << " r: " << _currentRadius +// << " layerphi: " << _currentLayerphi +// << std::endl; + return cID; } +double GridDriftChamber::phi(const CellID& cID) const { + CellID phiValue = _decoder->get(cID, m_phiID); + return binToPosition(phiValue, _currentLayerphi, m_offsetPhi); +} + +double GridDriftChamber::distanceTrackWire(const CellID& cID, const TVector3& hit_start, + const TVector3& hit_end) const { + + auto layerIndex = _decoder->get(cID, "layer"); + updateParams(layerIndex); + + double phi_start = phi(cID); + double phi_end = phi_start + returnAlpha(); + + TVector3 Wstart = returnWirePosition(phi_start, 1); + TVector3 Wend = returnWirePosition(phi_end, -1); + + TVector3 a = hit_end - hit_start; + TVector3 b = Wend - Wstart; + TVector3 c = Wstart - hit_start; + + double num = std::abs(c.Dot(a.Cross(b))); + double denum = (a.Cross(b)).Mag(); + + double DCA = 0; + + if (denum) { + DCA = num / denum; + } + + return DCA; +} + REGISTER_SEGMENTATION(GridDriftChamber) }