include/DDRec/DCH_info.h

0001 #ifndef DDREC_DCH_INFO_H
0002 #define DDREC_DCH_INFO_H
0003
0004 #include "TMath.h"
0005
0006 #include "DDRec/DetectorData.h"
0007 #include <map>
0008 #include "TVector3.h"
0009
0010 namespace dd4hep {  namespace rec {
0011
0012
0013 /* Data structure to store DCH geometry parameters
0014  *
0015  * Global parameters are members of this class
0016  *
0017  * Parameters of each layer are stored in the helper
0018  * class DCH_info_layer.
0019  *
0020  * The member database is a container which stores one
0021  * DCH_info_layer object per layer.
0022  *
0023  * To use this class, instantiate an object and define the global parameters.
0024  * Then call the method BuildLayerDatabase, which calculates
0025  * the parameters of each layer based on the global parameters.
0026  *
0027  * @author A. Tolosa Delgado, CERN
0028  * @date April 2024
0029  * @version Drift chamber v2
0030  *
0031  */
0032 struct DCH_info_struct
0033 {
0034 public:
0035     // use alias of types to show more clearly what the variable is
0036     // if everything is double, the code is not readable
0037     /// type for layer number
0038     using DCH_layer = int;
0039     /// tpye for lengths
0040     using DCH_length_t = double;
0041     /// tpye for angles
0042     using DCH_angle_t = double;
0043     /// Half length of the active volume
0044     DCH_length_t Lhalf = {0};
0045     /// Inner radius of the active volume
0046     DCH_length_t rin = {0};
0047     /// Outer radius of the active volume
0048     DCH_length_t rout = {0};
0049
0050     /// Inner guard wires radius
0051     DCH_length_t guard_inner_r_at_z0 = {0};
0052     /// Outer guard wires radius
0053     DCH_length_t guard_outer_r_at_zL2 = {0};
0054
0055     /// number of cells of first layer
0056     int ncell0 = {0};
0057     /// increment the number of cells for each superlayer as:
0058     ///   ncells(ilayer) = dch_ncell0 + increment*superlayer(ilayer)
0059     ///   See DCH_info::Get_nsuperlayer_minus_1(ilayer)
0060     int ncell_increment = {0};
0061
0062     /// cells within the same layer may be grouped into sectors, not in use atm
0063     int ncell_per_sector = {0};
0064
0065     /// input number of layers in each superlayer
0066     DCH_layer nlayersPerSuperlayer = {0};
0067     /// input number of superlayers
0068     /// superlayer is an abstract level of grouping layers used to
0069     /// parametrize the increment of cells in each layer
0070     DCH_layer nsuperlayers = {0};
0071     /// Calculated as dch_nlayersPerSuperlayer * dch_nsuperlayers
0072     DCH_layer nlayers = {0};
0073
0074     /// global twist angle
0075     /// alternating layers will change its sign
0076     DCH_angle_t  twist_angle = {0};
0077
0078     /// Cell width for the first layer
0079     double first_width = {0};
0080     /// Cell radius for the first layer
0081     DCH_length_t first_sense_r = {0};
0082
0083     void Set_lhalf(DCH_length_t _dch_Lhalf){Lhalf=_dch_Lhalf;}
0084     void Set_rin  (DCH_length_t _dch_rin  ){rin  = _dch_rin; }
0085     void Set_rout (DCH_length_t _dch_rout ){rout = _dch_rout;}
0086
0087     void Set_guard_rin_at_z0  (DCH_length_t _dch_rin_z0_guard  ){guard_inner_r_at_z0  = _dch_rin_z0_guard;   }
0088     void Set_guard_rout_at_zL2(DCH_length_t _dch_rout_zL2_guard){guard_outer_r_at_zL2 = _dch_rout_zL2_guard; }
0089
0090     void Set_ncell0              (int _ncell0              ){ncell0               = _ncell0;              }
0091     void Set_ncell_increment     (int _ncell_increment     ){ncell_increment      = _ncell_increment;     }
0092
0093     void Set_nlayersPerSuperlayer(int _nlayersPerSuperlayer){nlayersPerSuperlayer = _nlayersPerSuperlayer;}
0094     void Set_nsuperlayers        (int _nsuperlayers        ){nsuperlayers         = _nsuperlayers;        }
0095
0096     void Set_ncell_per_sector(int _ncell_per_sector){ncell_per_sector = _ncell_per_sector;}
0097
0098     void Set_twist_angle (DCH_length_t _dch_twist_angle ){twist_angle = _dch_twist_angle;}
0099
0100     void Set_first_width  (double _first_width  ){first_width   = _first_width;   }
0101     void Set_first_sense_r(double _first_sense_r){first_sense_r = _first_sense_r; }
0102
0103
0104     /// Get number of cells in a given layer
0105     ///   ncells = number of wires/2
0106     int Get_ncells(int ilayer){return database.at(ilayer).nwires/2;}
0107
0108     /// Get phi width for the twisted tube and the step (phi distance between cells)
0109     DCH_angle_t Get_phi_width(int ilayer){return (TMath::TwoPi()/Get_ncells(ilayer))*dd4hep::rad;}
0110
0111     /// phi positioning, adding offset for odd ilayers
0112     /// there is a staggering in phi for alternating layers, 0.25*cell_phi_width*(ilayer%2);
0113     DCH_angle_t Get_cell_phi_angle(int ilayer, int nphi){ return (Get_phi_width(ilayer) * (nphi + 0.25*(ilayer%2)));}
0114
0115     /// calculate superlayer for a given ilayer.
0116     /// WARNING: division of integers on purpose!
0117     int Get_nsuperlayer_minus_1(int ilayer){ return int((ilayer-1)/nlayersPerSuperlayer);}
0118
0119     /// Calculate radius at z=L/2 given at z=0
0120     DCH_length_t Radius_zLhalf(DCH_length_t r_z0) const {
0121         return r_z0/cos(twist_angle/2/dd4hep::rad);
0122     }
0123
0124     /// tan(stereoangle) = R(z=0)   / (L/2) * tan( twist_angle/2)
0125     DCH_angle_t stereoangle_z0(DCH_length_t r_z0) const {
0126         return atan( r_z0/Lhalf*tan(twist_angle/2/dd4hep::rad));
0127     }
0128
0129     /// tan(stereoangle) = R(z=L/2) / (L/2) * sin( twist_angle/2)
0130     DCH_angle_t stereoangle_zLhalf(DCH_length_t r_zLhalf) const {
0131         return atan( r_zLhalf/Lhalf*sin(twist_angle/2/dd4hep::rad));
0132     }
0133
0134     /// WireLength = 2*dch_Lhalf/cos(atan(Pitch_z0(r_z0)/(2*dch_Lhalf)))/cos(stereoangle_z0(r_z0))
0135     DCH_length_t WireLength(int nlayer, DCH_length_t r_z0) const {
0136         auto Pitch_z0 = database.at(nlayer).Pitch_z0(r_z0);
0137         return  2*Lhalf/cos(atan(Pitch_z0/(2*Lhalf)))/cos(stereoangle_z0(r_z0)/dd4hep::rad) ;
0138     };
0139
0140     /// Internal helper struct for defining the layer layout
0141     struct DCH_info_layer
0142     {
0143         /// layer number
0144         DCH_layer layer = {0};
0145         /// 2x number of cells in that layer
0146         int nwires = {0};
0147         /// cell parameter
0148         double height_z0 = {0.};
0149         /// cell parameter
0150         double width_z0 = {0.};
0151
0152         /// radius (cylindral coord) of sensitive wire
0153         DCH_length_t radius_sw_z0 = {0.};
0154
0155         /// radius (cylindral coord) of 'down' field wires
0156         DCH_length_t radius_fdw_z0 = {0.};
0157
0158         /// radius (cylindral coord) of 'up' field wires
0159         DCH_length_t radius_fuw_z0 = {0.};
0160
0161         /// some quantities are derived from previous-layer ones
0162         ///  stereo angle is positive for odd layer number
0163         bool IsStereoPositive() const {
0164             return (1 == layer%2);
0165         }
0166         /// calculate sign based on IsStereoPositive
0167         int StereoSign() const {
0168             return (IsStereoPositive()*2 - 1);
0169         }
0170
0171         /// separation between wires (along the circle)
0172         DCH_length_t Pitch_z0(DCH_length_t r_z0) const {
0173             return TMath::TwoPi()*r_z0/nwires;
0174         };
0175
0176     };
0177     /// map to store parameters for each layer
0178     std::map<DCH_layer, DCH_info_layer> database;
0179     bool IsDatabaseEmpty() const { return database.empty(); }
0180
0181     inline void BuildLayerDatabase();
0182     inline void Show_DCH_info_database(std::ostream& io) const;
0183
0184     /// Check if outer volume is not zero (0 < Lhalf*rout), and if the database was filled
0185     inline bool IsValid() const {return ((0 < Lhalf*rout) && (not IsDatabaseEmpty() ) );  }
0186     //--------
0187     // The following functions are used in the digitization/reconstruction
0188     // Notation: ILayer is the correlative number of the layer. Layer is reserved to number within a superlayer
0189     inline DCH_layer CalculateILayerFromCellIDFields(int layer, int superlayer) const { DCH_layer ilayer = layer + (this->nlayersPerSuperlayer)*superlayer + 1; return ilayer;}
0190     inline TVector3 Calculate_hitpos_to_wire_vector(int ilayer, int nphi, const TVector3& hit_position /*in cm*/) const;
0191     inline TVector3 Calculate_wire_vector_ez       (int ilayer, int nphi) const;
0192     inline TVector3 Calculate_wire_z0_point        (int ilayer, int nphi) const;
0193     inline double   Calculate_wire_phi_z0          (int ilayer, int nphi) const;
0194
0195 };
0196 typedef StructExtension<DCH_info_struct> DCH_info ;
0197 inline std::ostream& operator<<( std::ostream& io , const DCH_info& d ){d.Show_DCH_info_database(io); return io;}
0198
0199 inline void DCH_info_struct::BuildLayerDatabase()
0200 {
0201     // do not fill twice the database
0202     if( not this->IsDatabaseEmpty() ) return;
0203
0204     auto ff_check_positive_parameter = [](double p, std::string pname) -> void {
0205         if(p<=0)throw std::runtime_error(Form("DCH: %s must be positive",pname.c_str()));
0206         return;
0207     };
0208     ff_check_positive_parameter(this->rin  ,"inner radius");
0209     ff_check_positive_parameter(this->rout ,"outer radius");
0210     ff_check_positive_parameter(this->Lhalf,"half length" );
0211
0212     ff_check_positive_parameter(this->guard_inner_r_at_z0 ,"inner radius of guard wires" );
0213     ff_check_positive_parameter(this->guard_outer_r_at_zL2,"outer radius of guard wires" );
0214
0215
0216     ff_check_positive_parameter(this->ncell0,"ncells in the first layer" );
0217     ff_check_positive_parameter(this->ncell_increment,"ncells increment per superlayer" );
0218     ff_check_positive_parameter(this->ncell_per_sector,"ncells per sector" );
0219
0220     // if dch_ncell_per_sector is not divisor of dch_ncell0 and dch_ncell_increment
0221     // throw an error
0222     if( 0 != (ncell0 % ncell_per_sector) || 0 != (ncell_increment % ncell_per_sector) )
0223         throw std::runtime_error("dch_ncell_per_sector is not divisor of dch_ncell0 or dch_ncell_increment");
0224
0225     ff_check_positive_parameter(this->nsuperlayers,"number of superlayers" );
0226     ff_check_positive_parameter(this->nlayersPerSuperlayer,"number of layers per superlayer" );
0227
0228     /// nlayers = nsuperlayers * nlayersPerSuperlayer
0229     /// default: 112 = 14 * 8
0230     this->nlayers = this->nsuperlayers * this->nlayersPerSuperlayer;
0231
0232     ff_check_positive_parameter(this->first_width,"width of first layer cells" );
0233     ff_check_positive_parameter(this->first_sense_r,"radius of first layer cells" );
0234
0235
0236     // intialize layer 1 from input parameters
0237     {
0238         DCH_info_layer layer1_info;
0239         layer1_info.layer         = 1;
0240         layer1_info.nwires        = 2*this->ncell0;
0241         layer1_info.height_z0     = first_width;
0242         layer1_info.radius_sw_z0  = first_sense_r;
0243         layer1_info.radius_fdw_z0 = first_sense_r - 0.5*first_width;
0244         layer1_info.radius_fuw_z0 = first_sense_r + 0.5*first_width;
0245         layer1_info.width_z0      = TMath::TwoPi()*first_sense_r/this->ncell0;
0246
0247         this->database.emplace(layer1_info.layer, layer1_info);
0248     }
0249
0250     // some parameters of the following layer are calculated based on the previous ones
0251     // the rest are left as methods of DCH_info or DCH_info_layer class
0252     // loop over all layers, skipping the first one
0253     for(int ilayer = 2; ilayer<= this->nlayers; ++ilayer)
0254     {
0255         // initialize empty object, parameters are set later
0256         DCH_info_layer layer_info;
0257
0258         // the loop counter actually corresponds to the layer number
0259         layer_info.layer = ilayer;
0260         // nwires is twice the number of cells in this particular layer (ilayer)
0261         layer_info.nwires = 2*(this->ncell0 + this->ncell_increment*Get_nsuperlayer_minus_1(ilayer) );
0262
0263         // the previous layer info is needed to calculate parameters of current layer
0264         const auto& previousLayer = this->database.at(ilayer-1);
0265
0266         //calculate height_z0, radius_sw_z0
0267         {
0268             double h  = previousLayer.height_z0;
0269             double ru = previousLayer.radius_fuw_z0;
0270             double rd = previousLayer.radius_fdw_z0;
0271
0272             if(0 == Get_nsuperlayer_minus_1(ilayer))
0273                 layer_info.height_z0 = h*ru/rd;
0274             else
0275                 layer_info.height_z0 = TMath::TwoPi()*ru/(0.5*layer_info.nwires - TMath::Pi());
0276
0277             layer_info.radius_sw_z0 = 0.5*layer_info.height_z0 + ru;
0278         }
0279
0280         //calculate radius_fdw_z0, radius_fuw_z0, width_z0
0281         layer_info.radius_fdw_z0 = previousLayer.radius_fuw_z0;
0282         layer_info.radius_fuw_z0 = previousLayer.radius_fuw_z0 + layer_info.height_z0;
0283         layer_info.width_z0 = TMath::TwoPi()*layer_info.radius_sw_z0/(0.5*layer_info.nwires);
0284
0285         // according to expert prescription, width_z0 == height_z0
0286         if(fabs(layer_info.width_z0 - layer_info.height_z0)>1e-4)
0287             throw std::runtime_error("fabs(l.width_z0 - l.height_z0)>1e-4");
0288
0289
0290
0291         this->database.emplace(ilayer, layer_info);
0292     }
0293
0294     std::cout << "\t+ Total size of DCH database = " << database.size() << std::endl;
0295     return;
0296 }
0297
0298 inline void DCH_info_struct::Show_DCH_info_database(std::ostream & oss) const
0299 {
0300     oss << "\n";
0301     oss << "Global parameters of DCH:\n";
0302     oss << "\tGas, half length/mm = " << Lhalf/dd4hep::mm << '\n';
0303     oss << "\tGas, radius in/mm  = " << rin/dd4hep::mm << '\n';
0304     oss << "\tGas, radius out/mm = " << rout/dd4hep::mm<< '\n';
0305     oss << "\tGuard, radius in(z=0)/mm  = " << guard_inner_r_at_z0/dd4hep::mm << '\n';
0306     oss << "\tGuard, radius out(z=L/2)/mm = " << guard_outer_r_at_zL2/dd4hep::mm << '\n';
0307     oss << "\n";
0308     oss << "\tTwist angle (2*alpha) / deg = " << twist_angle/dd4hep::deg << '\n';
0309     oss << "\n";
0310     oss << "\tN superlayers = " << nsuperlayers << '\n';
0311     oss << "\tN layers per superlayer = " << nlayersPerSuperlayer << '\n';
0312     oss << "\tN layers = " << nlayers << '\n';
0313     oss << "\n";
0314     oss << "\tN cells layer1 = " << ncell0 << '\n';
0315     oss << "\tN cells increment per superlayer = " << ncell_increment << '\n';
0316     oss << "\tN cells per sector = " << ncell_per_sector << '\n';
0317     oss << "\n";
0318     oss << "Layer parameters of DCH:\n";
0319     oss
0320             << "\t" << "layer"
0321             << "\t" << "nwires"
0322             << "\t" << "height_z0/mm"
0323             << "\t" << "width_z0/mm"
0324             << "\t" << "radius_fdw_z0/mm"
0325             << "\t" << "radius_sw_z0/mm"
0326             << "\t" << "radius_fuw_z0/mm"
0327             << "\t" << "stereoangle_z0/deg"
0328             << "\t" << "Pitch_z0/mm"
0329             << "\t" << "radius_sw_zLhalf/mm"
0330             << "\t" << "WireLength/mm"
0331             << "\n" << std::endl;
0332
0333     if( this->IsDatabaseEmpty() )
0334     {
0335         oss << "\nDatabase empty\n";
0336         return;
0337     }
0338
0339     for(const auto& [nlayer, l]  : database )
0340     {
0341         oss
0342                 << "\t" << l.layer
0343                 << "\t" << l.nwires
0344                 << "\t" << l.height_z0/dd4hep::mm
0345                 << "\t" << l.width_z0/dd4hep::mm
0346                 << "\t" << l.radius_fdw_z0/dd4hep::mm
0347                 << "\t" << l.radius_sw_z0/dd4hep::mm
0348                 << "\t" << l.radius_fuw_z0/dd4hep::mm
0349                 << "\t" << l.StereoSign()*this->stereoangle_z0(l.radius_sw_z0)/dd4hep::deg
0350                 << "\t" << l.Pitch_z0(l.radius_sw_z0)/dd4hep::mm
0351                 << "\t" << this->Radius_zLhalf(l.radius_sw_z0)/dd4hep::mm
0352                 << "\t" << this->WireLength(l.layer,l.radius_sw_z0)/dd4hep::mm
0353                 << "\n" << std::endl;
0354     }
0355     return;
0356 }
0357
0358
0359 ///////////////////////////////////////////////////////////////////////////////////////
0360 /////       Ancillary functions for calculating the distance to the wire       ////////
0361 ///////////////////////////////////////////////////////////////////////////////////////
0362
0363 inline TVector3 DCH_info_struct::Calculate_wire_vector_ez(int ilayer, int nphi) const {
0364   auto& l = this->database.at(ilayer);
0365
0366   // See original paper Hoshina et al, Computer Physics Communications 153 (2003) 3
0367   // eq. 2.9, for the definition of ez, vector along the wire
0368
0369   // initialize some variables
0370   int    stereosign = l.StereoSign();
0371   double rz0        = l.radius_sw_z0;
0372   double dphi       = this->twist_angle;
0373   // kappa is the same as in eq. 2.9
0374   double kappa = (1. / this->Lhalf) * tan(dphi / 2);
0375
0376   //--- calculating wire position
0377   // the points p1 and p2 correspond to the ends of the wire
0378
0379   // point 1
0380   double x1 = rz0;                                      // m
0381   double y1 = -stereosign * rz0 * kappa * this->Lhalf;  // m
0382   double z1 = -this->Lhalf;                             // m
0383
0384   TVector3 p1(x1, y1, z1);
0385
0386   // point 2
0387   double x2 = rz0;                                     // m
0388   double y2 = stereosign * rz0 * kappa * this->Lhalf;  // m
0389   double z2 = this->Lhalf;                             // m
0390
0391   TVector3 p2(x2, y2, z2);
0392
0393   // calculate phi rotation of whole twisted tube, ie, rotation at z=0
0394   double phi_z0 = Calculate_wire_phi_z0(ilayer, nphi);
0395   p1.RotateZ(phi_z0);
0396   p2.RotateZ(phi_z0);
0397
0398   //--- end calculating wire position
0399
0400   return (p2 - p1).Unit();
0401 }
0402
0403 inline TVector3 DCH_info_struct::Calculate_wire_z0_point(int ilayer, int nphi) const {
0404   auto&    l   = this->database.at(ilayer);
0405   double   rz0 = l.radius_sw_z0;
0406   TVector3 p1(rz0, 0, 0);
0407   double   phi_z0 = Calculate_wire_phi_z0(ilayer, nphi);
0408   p1.RotateZ(phi_z0);
0409   return p1;
0410 }
0411
0412 // calculate phi rotation of whole twisted tube, ie, rotation at z=0
0413 inline double DCH_info_struct::Calculate_wire_phi_z0(int ilayer, int nphi) const {
0414   auto&  l       = this->database.at(ilayer);
0415   int    ncells  = l.nwires / 2;
0416   double phistep = TMath::TwoPi() / ncells;
0417   double phi_z0  = (nphi + 0.25 * (l.layer % 2)) * phistep;
0418   return phi_z0;
0419 }
0420
0421 ///////////////////////////////////////////////////////////////////////////////////////
0422 ///////////////////////  Calculate vector from hit position to wire   /////////////////
0423 ///////////////////////////////////////////////////////////////////////////////////////
0424 inline TVector3 DCH_info_struct::Calculate_hitpos_to_wire_vector(int ilayer, int nphi, const TVector3& hit_position /*in cm*/) const {
0425   // Solution distance from a point to a line given here:
0426   // https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Vector_formulation
0427   TVector3 n = this->Calculate_wire_vector_ez(ilayer, nphi);
0428   TVector3 a = this->Calculate_wire_z0_point(ilayer, nphi);
0429   // Remember using cm as natural units of DD4hep consistently!
0430   // TVector3 p {hit_position.x()*MM_TO_CM,hit_position.y()*MM_TO_CM,hit_position.z()*MM_TO_CM};
0431
0432   TVector3 a_minus_p       = a - hit_position;
0433   double   a_minus_p_dot_n = a_minus_p.Dot(n);
0434   TVector3 scaled_n        = a_minus_p_dot_n * n;
0435   //hit_to_wire_vector = a_minus_p - scaled_n;
0436   return (a_minus_p - scaled_n);
0437 }
0438
0439 }} // end namespace dd4hep::rec::
0440
0441 #endif // DDREC_DCH_INFO_H