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
0009 namespace dd4hep {  namespace rec {
0010
0011
0012 /* Data structure to store DCH geometry parameters
0013  *
0014  * Global parameters are members of this class
0015  *
0016  * Parameters of each layer are stored in the helper
0017  * class DCH_info_layer.
0018  *
0019  * The member database is a container which stores one
0020  * DCH_info_layer object per layer.
0021  *
0022  * To use this class, instantiate an object and define the global parameters.
0023  * Then call the method BuildLayerDatabase, which calculates
0024  * the parameters of each layer based on the global parameters.
0025  *
0026  * @author A. Tolosa Delgado, CERN
0027  * @date April 2024
0028  * @version Drift chamber v2
0029  *
0030  */
0031 struct DCH_info_struct
0032 {
0033 public:
0034     // use alias of types to show more clearly what the variable is
0035     // if everything is double, the code is not readable
0036     /// type for layer number
0037     using DCH_layer = int;
0038     /// tpye for lengths
0039     using DCH_length_t = double;
0040     /// tpye for angles
0041     using DCH_angle_t = double;
0042     /// Half length of the active volume
0043     DCH_length_t Lhalf = {0};
0044     /// Inner radius of the active volume
0045     DCH_length_t rin = {0};
0046     /// Outer radius of the active volume
0047     DCH_length_t rout = {0};
0048
0049     /// Inner guard wires radius
0050     DCH_length_t guard_inner_r_at_z0 = {0};
0051     /// Outer guard wires radius
0052     DCH_length_t guard_outer_r_at_zL2 = {0};
0053
0054     /// number of cells of first layer
0055     int ncell0 = {0};
0056     /// increment the number of cells for each superlayer as:
0057     ///   ncells(ilayer) = dch_ncell0 + increment*superlayer(ilayer)
0058     ///   See DCH_info::Get_nsuperlayer_minus_1(ilayer)
0059     int ncell_increment = {0};
0060
0061     /// cells within the same layer may be grouped into sectors, not in use atm
0062     int ncell_per_sector = {0};
0063
0064     /// input number of layers in each superlayer
0065     DCH_layer nlayersPerSuperlayer = {0};
0066     /// input number of superlayers
0067     /// superlayer is an abstract level of grouping layers used to
0068     /// parametrize the increment of cells in each layer
0069     DCH_layer nsuperlayers = {0};
0070     /// Calculated as dch_nlayersPerSuperlayer * dch_nsuperlayers
0071     DCH_layer nlayers = {0};
0072
0073     /// global twist angle
0074     /// alternating layers will change its sign
0075     DCH_angle_t  twist_angle = {0};
0076
0077     /// Cell width for the first layer
0078     double first_width = {0};
0079     /// Cell radius for the first layer
0080     DCH_length_t first_sense_r = {0};
0081
0082     void Set_lhalf(DCH_length_t _dch_Lhalf){Lhalf=_dch_Lhalf;}
0083     void Set_rin  (DCH_length_t _dch_rin  ){rin  = _dch_rin; }
0084     void Set_rout (DCH_length_t _dch_rout ){rout = _dch_rout;}
0085
0086     void Set_guard_rin_at_z0  (DCH_length_t _dch_rin_z0_guard  ){guard_inner_r_at_z0  = _dch_rin_z0_guard;   }
0087     void Set_guard_rout_at_zL2(DCH_length_t _dch_rout_zL2_guard){guard_outer_r_at_zL2 = _dch_rout_zL2_guard; }
0088
0089     void Set_ncell0              (int _ncell0              ){ncell0               = _ncell0;              }
0090     void Set_ncell_increment     (int _ncell_increment     ){ncell_increment      = _ncell_increment;     }
0091
0092     void Set_nlayersPerSuperlayer(int _nlayersPerSuperlayer){nlayersPerSuperlayer = _nlayersPerSuperlayer;}
0093     void Set_nsuperlayers        (int _nsuperlayers        ){nsuperlayers         = _nsuperlayers;        }
0094
0095     void Set_ncell_per_sector(int _ncell_per_sector){ncell_per_sector = _ncell_per_sector;}
0096
0097     void Set_twist_angle (DCH_length_t _dch_twist_angle ){twist_angle = _dch_twist_angle;}
0098
0099     void Set_first_width  (double _first_width  ){first_width   = _first_width;   }
0100     void Set_first_sense_r(double _first_sense_r){first_sense_r = _first_sense_r; }
0101
0102
0103     /// Get number of cells in a given layer
0104     ///   ncells = number of wires/2
0105     int Get_ncells(int ilayer){return database.at(ilayer).nwires/2;}
0106
0107     /// Get phi width for the twisted tube and the step (phi distance between cells)
0108     DCH_angle_t Get_phi_width(int ilayer){return (TMath::TwoPi()/Get_ncells(ilayer))*dd4hep::rad;}
0109
0110     /// phi positioning, adding offset for odd ilayers
0111     /// there is a staggering in phi for alternating layers, 0.25*cell_phi_width*(ilayer%2);
0112     DCH_angle_t Get_cell_phi_angle(int ilayer, int nphi){ return (Get_phi_width(ilayer) * (nphi + 0.25*(ilayer%2)));}
0113
0114     /// calculate superlayer for a given ilayer.
0115     /// WARNING: division of integers on purpose!
0116     int Get_nsuperlayer_minus_1(int ilayer){ return int((ilayer-1)/nlayersPerSuperlayer);}
0117
0118     /// Calculate radius at z=L/2 given at z=0
0119     DCH_length_t Radius_zLhalf(DCH_length_t r_z0) const {
0120         return r_z0/cos(twist_angle/2/dd4hep::rad);
0121     }
0122
0123     /// tan(stereoangle) = R(z=0)   / (L/2) * tan( twist_angle/2)
0124     DCH_angle_t stereoangle_z0(DCH_length_t r_z0) const {
0125         return atan( r_z0/Lhalf*tan(twist_angle/2/dd4hep::rad));
0126     }
0127
0128     /// tan(stereoangle) = R(z=L/2) / (L/2) * sin( twist_angle/2)
0129     DCH_angle_t stereoangle_zLhalf(DCH_length_t r_zLhalf) const {
0130         return atan( r_zLhalf/Lhalf*sin(twist_angle/2/dd4hep::rad));
0131     }
0132
0133     /// WireLength = 2*dch_Lhalf/cos(atan(Pitch_z0(r_z0)/(2*dch_Lhalf)))/cos(stereoangle_z0(r_z0))
0134     DCH_length_t WireLength(int nlayer, DCH_length_t r_z0) const {
0135         auto Pitch_z0 = database.at(nlayer).Pitch_z0(r_z0);
0136         return  2*Lhalf/cos(atan(Pitch_z0/(2*Lhalf)))/cos(stereoangle_z0(r_z0)/dd4hep::rad) ;
0137     };
0138
0139     /// Internal helper struct for defining the layer layout
0140     struct DCH_info_layer
0141     {
0142         /// layer number
0143         DCH_layer layer = {0};
0144         /// 2x number of cells in that layer
0145         int nwires = {0};
0146         /// cell parameter
0147         double height_z0 = {0.};
0148         /// cell parameter
0149         double width_z0 = {0.};
0150
0151         /// radius (cylindral coord) of sensitive wire
0152         DCH_length_t radius_sw_z0 = {0.};
0153
0154         /// radius (cylindral coord) of 'down' field wires
0155         DCH_length_t radius_fdw_z0 = {0.};
0156
0157         /// radius (cylindral coord) of 'up' field wires
0158         DCH_length_t radius_fuw_z0 = {0.};
0159
0160         /// some quantities are derived from previous-layer ones
0161         ///  stereo angle is positive for odd layer number
0162         bool IsStereoPositive() const {
0163             return (1 == layer%2);
0164         }
0165         /// calculate sign based on IsStereoPositive
0166         int StereoSign() const {
0167             return (IsStereoPositive()*2 - 1);
0168         }
0169
0170         /// separation between wires (along the circle)
0171         DCH_length_t Pitch_z0(DCH_length_t r_z0) const {
0172             return TMath::TwoPi()*r_z0/nwires;
0173         };
0174
0175     };
0176     /// map to store parameters for each layer
0177     std::map<DCH_layer, DCH_info_layer> database;
0178     bool IsDatabaseEmpty() const { return database.empty(); }
0179
0180     inline void BuildLayerDatabase();
0181     inline void Show_DCH_info_database(std::ostream& io) const;
0182
0183
0184 };
0185 typedef StructExtension<DCH_info_struct> DCH_info ;
0186 std::ostream& operator<<( std::ostream& io , const DCH_info& d ){d.Show_DCH_info_database(io); return io;}
0187
0188 inline void DCH_info_struct::BuildLayerDatabase()
0189 {
0190     // do not fill twice the database
0191     if( not this->IsDatabaseEmpty() ) return;
0192
0193     auto ff_check_positive_parameter = [](double p, std::string pname) -> void {
0194         if(p<=0)throw std::runtime_error(Form("DCH: %s must be positive",pname.c_str()));
0195         return;
0196     };
0197     ff_check_positive_parameter(this->rin  ,"inner radius");
0198     ff_check_positive_parameter(this->rout ,"outer radius");
0199     ff_check_positive_parameter(this->Lhalf,"half length" );
0200
0201     ff_check_positive_parameter(this->guard_inner_r_at_z0 ,"inner radius of guard wires" );
0202     ff_check_positive_parameter(this->guard_outer_r_at_zL2,"outer radius of guard wires" );
0203
0204
0205     ff_check_positive_parameter(this->ncell0,"ncells in the first layer" );
0206     ff_check_positive_parameter(this->ncell_increment,"ncells increment per superlayer" );
0207     ff_check_positive_parameter(this->ncell_per_sector,"ncells per sector" );
0208
0209     // if dch_ncell_per_sector is not divisor of dch_ncell0 and dch_ncell_increment
0210     // throw an error
0211     if( 0 != (ncell0 % ncell_per_sector) || 0 != (ncell_increment % ncell_per_sector) )
0212         throw std::runtime_error("dch_ncell_per_sector is not divisor of dch_ncell0 or dch_ncell_increment");
0213
0214     ff_check_positive_parameter(this->nsuperlayers,"number of superlayers" );
0215     ff_check_positive_parameter(this->nlayersPerSuperlayer,"number of layers per superlayer" );
0216
0217     /// nlayers = nsuperlayers * nlayersPerSuperlayer
0218     /// default: 112 = 14 * 8
0219     this->nlayers = this->nsuperlayers * this->nlayersPerSuperlayer;
0220
0221     ff_check_positive_parameter(this->first_width,"width of first layer cells" );
0222     ff_check_positive_parameter(this->first_sense_r,"radius of first layer cells" );
0223
0224
0225     // intialize layer 1 from input parameters
0226     {
0227         DCH_info_layer layer1_info;
0228         layer1_info.layer         = 1;
0229         layer1_info.nwires        = 2*this->ncell0;
0230         layer1_info.height_z0     = first_width;
0231         layer1_info.radius_sw_z0  = first_sense_r;
0232         layer1_info.radius_fdw_z0 = first_sense_r - 0.5*first_width;
0233         layer1_info.radius_fuw_z0 = first_sense_r + 0.5*first_width;
0234         layer1_info.width_z0      = TMath::TwoPi()*first_sense_r/this->ncell0;
0235
0236         this->database.emplace(layer1_info.layer, layer1_info);
0237     }
0238
0239     // some parameters of the following layer are calculated based on the previous ones
0240     // the rest are left as methods of DCH_info or DCH_info_layer class
0241     // loop over all layers, skipping the first one
0242     for(int ilayer = 2; ilayer<= this->nlayers; ++ilayer)
0243     {
0244         // initialize empty object, parameters are set later
0245         DCH_info_layer layer_info;
0246
0247         // the loop counter actually corresponds to the layer number
0248         layer_info.layer = ilayer;
0249         // nwires is twice the number of cells in this particular layer (ilayer)
0250         layer_info.nwires = 2*(this->ncell0 + this->ncell_increment*Get_nsuperlayer_minus_1(ilayer) );
0251
0252         // the previous layer info is needed to calculate parameters of current layer
0253         const auto& previousLayer = this->database.at(ilayer-1);
0254
0255         //calculate height_z0, radius_sw_z0
0256         {
0257             double h  = previousLayer.height_z0;
0258             double ru = previousLayer.radius_fuw_z0;
0259             double rd = previousLayer.radius_fdw_z0;
0260
0261             if(0 == Get_nsuperlayer_minus_1(ilayer))
0262                 layer_info.height_z0 = h*ru/rd;
0263             else
0264                 layer_info.height_z0 = TMath::TwoPi()*ru/(0.5*layer_info.nwires - TMath::Pi());
0265
0266             layer_info.radius_sw_z0 = 0.5*layer_info.height_z0 + ru;
0267         }
0268
0269         //calculate radius_fdw_z0, radius_fuw_z0, width_z0
0270         layer_info.radius_fdw_z0 = previousLayer.radius_fuw_z0;
0271         layer_info.radius_fuw_z0 = previousLayer.radius_fuw_z0 + layer_info.height_z0;
0272         layer_info.width_z0 = TMath::TwoPi()*layer_info.radius_sw_z0/(0.5*layer_info.nwires);
0273
0274         // according to expert prescription, width_z0 == height_z0
0275         if(fabs(layer_info.width_z0 - layer_info.height_z0)>1e-4)
0276             throw std::runtime_error("fabs(l.width_z0 - l.height_z0)>1e-4");
0277
0278
0279
0280         this->database.emplace(ilayer, layer_info);
0281     }
0282
0283     std::cout << "\t+ Total size of DCH database = " << database.size() << std::endl;
0284     return;
0285 }
0286
0287 inline void DCH_info_struct::Show_DCH_info_database(std::ostream & oss) const
0288 {
0289     oss << "\n";
0290     oss << "Global parameters of DCH:\n";
0291     oss << "\tGas, half length/mm = " << Lhalf/dd4hep::mm << '\n';
0292     oss << "\tGas, radius in/mm  = " << rin/dd4hep::mm << '\n';
0293     oss << "\tGas, radius out/mm = " << rout/dd4hep::mm<< '\n';
0294     oss << "\tGuard, radius in(z=0)/mm  = " << guard_inner_r_at_z0/dd4hep::mm << '\n';
0295     oss << "\tGuard, radius out(z=L/2)/mm = " << guard_outer_r_at_zL2/dd4hep::mm << '\n';
0296     oss << "\n";
0297     oss << "\tTwist angle (2*alpha) / deg = " << twist_angle/dd4hep::deg << '\n';
0298     oss << "\n";
0299     oss << "\tN superlayers = " << nsuperlayers << '\n';
0300     oss << "\tN layers per superlayer = " << nlayersPerSuperlayer << '\n';
0301     oss << "\tN layers = " << nlayers << '\n';
0302     oss << "\n";
0303     oss << "\tN cells layer1 = " << ncell0 << '\n';
0304     oss << "\tN cells increment per superlayer = " << ncell_increment << '\n';
0305     oss << "\tN cells per sector = " << ncell_per_sector << '\n';
0306     oss << "\n";
0307     oss << "Layer parameters of DCH:\n";
0308     oss
0309             << "\t" << "layer"
0310             << "\t" << "nwires"
0311             << "\t" << "height_z0/mm"
0312             << "\t" << "width_z0/mm"
0313             << "\t" << "radius_fdw_z0/mm"
0314             << "\t" << "radius_sw_z0/mm"
0315             << "\t" << "radius_fuw_z0/mm"
0316             << "\t" << "stereoangle_z0/deg"
0317             << "\t" << "Pitch_z0/mm"
0318             << "\t" << "radius_sw_zLhalf/mm"
0319             << "\t" << "WireLength/mm"
0320             << "\n" << std::endl;
0321
0322     if( this->IsDatabaseEmpty() )
0323     {
0324         oss << "\nDatabase empty\n";
0325         return;
0326     }
0327
0328     for(const auto& [nlayer, l]  : database )
0329     {
0330         oss
0331                 << "\t" << l.layer
0332                 << "\t" << l.nwires
0333                 << "\t" << l.height_z0/dd4hep::mm
0334                 << "\t" << l.width_z0/dd4hep::mm
0335                 << "\t" << l.radius_fdw_z0/dd4hep::mm
0336                 << "\t" << l.radius_sw_z0/dd4hep::mm
0337                 << "\t" << l.radius_fuw_z0/dd4hep::mm
0338                 << "\t" << l.StereoSign()*this->stereoangle_z0(l.radius_sw_z0)/dd4hep::deg
0339                 << "\t" << l.Pitch_z0(l.radius_sw_z0)/dd4hep::mm
0340                 << "\t" << this->Radius_zLhalf(l.radius_sw_z0)/dd4hep::mm
0341                 << "\t" << this->WireLength(l.layer,l.radius_sw_z0)/dd4hep::mm
0342                 << "\n" << std::endl;
0343     }
0344     return;
0345 }
0346 }} // end namespace dd4hep::rec::
0347
0348 #endif // DDREC_DCH_INFO_H