Drift Chamber Digitization -------------------------- Maurice Foucher and documentation written by Ed Frank Steve Schaffner Art Snyder plus others who have written documentation in gndcha, DchData, and DchReco 17 October 1996 The simulation and geometry of the Drift CHamber (DCH) within BBsim has already been describe elsewhere. The purpose of the digitization code is to use the output of the simulation to produce data that the reconstruction can use, for example to perform track finding. This data we call digitizations, digi for short. The ultimate goal of the digitization code is to produce a digi that looks like what the hardware will produce. For the drift chamber this is a digitized pulse height versus time for each wire. At the moment the current digi is just a time a wire was hit as well as an energy deposited on that wire. Currently the simulation produces two types of hits. The first are wire hits. This data structure, DchGWirHit, is used primarily by the trigger group. The second type of hit structure (layer hits) is DchGLayHit from which digis for the reconstruction are produced. These classes as well as others that are used in the digitization are described below. The digitization code is written entirely in C++. The information contained in this note was determined from release 1.2.1 and package versions gndcha/V01-00-05 DchData/V00-03-04 DchReco/V00-00-04 DchGWirHit class ---------------- This is a description of the relevant private members of DchGWirHit. The DchGWirHit class is calculated by the routine gndcha/dc_wirhit.F. The C structures are filled by gndcha/DchHitsToDs.F. private members _rcell radius of hit cell _phcell phi of hit cell _dedx energy loss in cell _drift drift distance in cell _idclay layer number _idcell cell number -igtr track number DchGLayHit class ---------------- This is a description of the relevant private members of DchGLayHit. The DchGLayHit class is calculated by the routine gndcha/dc_step.F. The C structures are filled by gndcha/DchHitsToDs.F. The DchGLayHit class contains, amongst other things, DbiEvent/dchGLayHit.dhh which contains _x, _y, _z coordinates (cm) of track at entrance of layer _px, _py, _pz components of momentum (GeV/c) at entrance of layer _path path length (cm) in layer _dedx energy loss in keV computed by GEANT, NOT in dc_dedx routine. _time GEANT time of flight (sec) _type GEANT particle number _iflg _iflg = 1 normal track through layer = 2 looper = 3 track left experimental setup _layer layer number _igtr track number set in DbiEvent/HitsToDs.F by call to gfhits. DchTrkSeg class --------------- The class DchTrkSeg takes as input layer hits from the class DchGLayHit, constructs a track segment, and stores the following information as private members. private members _glayhit pointer to hit in DchGLayHit _x, _y, _z coordinates (cm) of hit from DchGLayHit _px, _py, _pz components of momentum (GeV/c) from DchGLayHit _path pathlength (cm) through layer from DchGLayHit _dedx energy loss in layer (keV) from DchGLayHit _time GEANT time of flight (sec) from DchGLayHit _distance distance of closest approach (cm) of track segment to a wire. This is set in DchWire::CellBound _type GEANT particle number from DchGLayHit _iflg _iflg = 1 normal track through layer = 2 looper = 3 track left experimental setup from DchGLayHit _layer layer number from DchGLayHit _igtr track number from DchGLayHit It also contains cached properties of the particle _curv curvature (1./cm) defined as 1/radius of curvature The sign is determined by which direction the particle is moving in the phi coordinate. For example a positive particle has negative curvature. _pxy transverse momentum(GeV/c), sqrt(_px*_px+_py*_py) _charge charge in units of e _mass mass (GeV/c**2) _momentum momentum (GeV/c) _energy energy (GeV) _beta azimuthal angle (radians) between -pi and pi of the trajectory, atan2(_py,_px) _cosbeta cos(_beta) _sinbeta sin(_beta) _tanl tan lambda = _pz/_pxy, tangent of pitch angle _secantlsq (secant(lambda))**2 _secantl secant(lambda) _alpha angle the particle turns in the plane perpendicular to the magnetic field _cosalpha cos(_alpha) _sinalpha sin(_alpha) _sarc pathlength (cm) of track segment set by DchTrkSeg::TurningAngle _radius sqrt(_x*_x+_y*_y) of track segment _phi azimuthal angle phi (radians) between -pi and pi of track segment _cosphi cos(_phi) = _x/_radius _sinphi sin(_phi) = _y/_radius _cached = 1 if variables are cached field = 1.5, B field in Tesla hard coded the non trivial functions of DchTrkSeg void cache() function to (re)calculate useful variables and cache them void Turn(x,y,z) moves position of track segment x,y,z along helix void Turn(x,y,z,px,py,pz) calls Turn(x,y,z), rotates px and py by angle _alpha projected turning angle of helix void Clip(sbeg,send) change DchTrkSeg path length to send-sbeg and recalculate _dedx as simple ratio. call Move(sarc) to change rest of DchTrkSeg. void Move(sarc) change DchTrkSeg based on new path length, sarc. void TurningAngle(sarc) sets DchTrkSeg arc length _sarc = sarc calculates _alpha, projected turning angle of track void TurnTo(sarc,x,y,z) calls TurningAngle(sarc), Turn(x,y,z) DchGeom class ------------- private members _zoffset chamber offset _zlen chamber length _cellHeight cell height DchLayer _dclayer[40] 40 layers for chamber DchLayer* _nextlay[40] pointer to next layer DchLayer* _nextlayinvw[40] pointer to next layer in same view DchLayer* _prevlay[40] pointer to previous layer in view DchLayer* _prevlayinvw[40] pointer to previous layer in same view _firstLayer index of 1st DCH layer =1 for now _nLayer # DCH layers (=40) _init initialization control = 1 if initialized = 0 if not _edge edge effect allowance = 1, set in constructor _iterlay iterating layer _wirebeg starting wire _wirelim limiting wire, i.e, ending wire++ nontrivial functions void nextWireInit(DchTrkSeg*) sets up variables to calculate range of wires that could be hit. int nextWire(wire) set wire to next wire, return 1 when there are more wires to come and 0 when it is finished. int modWire(layer, wire) returns wire number modulo the number of wires in a layer. int whichWire(layer, x, y, z) returns the nearest wire to x,y,z int cellBound(layer, wire, smear, DchTrkSeg*) pass on to DchLayer cellBound function int cellBound(wire, smear, DchTrkSeg*) pass on to DchWire cellBound function int CellBound(smear, DchTrkSeg*) void transToDc(DchTrkSeg) translate to dch coordinates, simple z shift z --> z-_zoffset void buildpointers() build pointers to next layer and next layer in view pointer DchLayer Class -------------- private members _exist this layer exists? = 0 nope doesn't exist = 1 yup exists _layer layer number counts from 0 _nwires number of wires _rend radius at chamber mid point _phioffset azimuthal offset of wire 0 _stereo tan of stereo angle _zlen longitudinal length of chamber _delphi deltaPhi of wire from midpoint to chamber end _cellheight full height _stdip change in radius from mid to end _qOffset Offsets that minimize Q-forces on sense wire (from C. Hearty) DchDriftCalib dCalib time-to-distance object DchHit **hitmap array of pointers to hits on wires nontrivial functions float resol(drift) return DchDriftCalib::resol(drift) int cellBound(wire,smear,DchTrkSeg) return DchWire::CellBound(smear,DchTrkSeg) float timeToDistance(time,angle=0.0) return DchDriftCalib::timeToDistance(time,angle=0.0) float distanceToTime(distance, angle=0.0) return DchDriftCalib::distanceToTime(distance,angle=0.0) DchWire Class ------------- private data members _wire wire number _layer layer number _layerptr layer pointer _phi0 phi at z=0 _radius0 radius at z=0 _zlen length in z _stereo stereo angle _xwire0 x at z=0 _ywire0 y at z=0 _cosphi0 cos(_phi0) _sinphi0 sin(_phi0) _rcell cell size hard coded gas properties z = 15.0 gas nuclei average charge a = 31.0 average atomic number density = 0.00084 density nelperkev = 30.1 electrons per kev deposited fano = 0.18 fano factor (whatever that is) teta = 0.6 polya avalanche fluctuation parameter zpolyamax = 5.0 max to use for polya distribution functions float rang() gaussian random numbers float flat() random numbers flat distribution nontrivial functions float distance(x,y,z) calculate the distance from x,y,z to a wire int isitin(x,y,z) = 0 if x,y,z is not in the cell = 1 if x,y,z is in the cell circular cell is assumed. int isitin(s,DchTrkSeg) calls TurnTo(s,x,y,z) and then determines if hit is in the cell by calling isitin(x,y,z). float tfromd(dist,angle=0.0) returns time for an ion a distance dist to drift to the wire. Angle not used at moment. calls DchLayer::distanceToTime(dist,angle) which then calls DchDriftCalib::distanceToTime(dist,angle) float smearDist(dist) smears dist by DchDriftCalib::resol(drift) float corrTime(time,DchTrkSeg) correct time for time of flight to layer and wire propagation speed. assumes c = 30cm/ns wire c/(velocity of propagation) = 1.2 void findpoca(DchTrkSeg,poca,doca) find the distance of closest approach (doca) between a hit and a wire and move the hit to the point of closest approach(poca). Use binary search. The poca is the distance along the DchTrkSeg to the point of closest approach int CellBound(smear,DchTrkSeg) clips the DchTrkSeg to fit inside a cell and performs Landau fluctuations and distance to time conversion etc float boundary(DchTrkSeg,sin,step) find boundary in direction of step void fluctuate(DchTrkSeg) Does Landau fluctuation of DchTrkSeg by calling glando routine (converted to c++) from geant. calls fanofluct(nelec), avalanche(nelec). float fanofluct(nelec) performs fano fluctuation of number of electrons, nelec = nelec + (gaussian random number * sqrt(fano*nelect)) float avalanche(nelec) calculates the number of electrons in the avalanche using polya distribution. float zpolya() generate a polya distribution. DchDriftCalib class ------------------- The DchDriftCalib class contains the distance to time conversion relationship to convert drift distances to time as well as the reverse relationship to get back a time in the reconstruction. The resolution in drift distance is also given in this class. protected members iModel time to distance relationship model = 3 read from file = 2 default version = 1 or other, use 30 x 10**-4 cm/ns drift velocity ntimebin number of times per bin ddist array of distances (read from file) tstart start of time range covered by time bins tend end of time range covered by time bins tbinwidinv 1 / (width of time bins, in psec) nontrivial functions double timeToDistance(time, angle) parametrized curve as follows, postscript figure available as well. time(ns) timeToDistance = slope(um/ns) * time(ns) + offset(um) 0 - 100 28.70 0 100 - 200 21.70 700 200 - 300 16.80 1680 300 - 400 13.60 2640 400 - 500 8.40 4720 500 - 600 3.10 7370 600 - inf 1.16 8530 double distanceToTime(dist, angle) uses reverse parametrization of timeToDistance float resol(dist) calculates hit resolution as function of drift distance drift distance(mm) resolution(mm) 0 - 1 0.210 1 - 2 0.190 2 - 3 0.150 3 - 4 0.145 4 - 5 0.140 5 - 6 0.135 6 - 7 0.155 7 - 8 0.195 8 - 9 0.245 9 - infinity 0.300 How to actually use this in the reconstruction ---------------------------------------------- Using the Framework package to run the reconstruction, the following modules must be defined and run from the file AppUserBuild.cc. An example of this is given in DchReco/AppUserBuild.cc. The results of this reconstruction is to produce drift chamber digis and run the track finding. 1. GPid (in RecoUtils) Particle id class used to return information about particles. 2. MakePointers (in G3Data) Makes a list of DchGLayHits from the simulation. 3. DchInitGeom (in DchReco) Creates a DchGeom object and stores a pointer to it in AbsEnv. It reads the file DchGeom.data in the DchData package. This file needs to be in the working directory. It then creates DchLayer objects for each of the 40 layers with the properties of the wires. The file contains the data stored in the private members of the DchLayer class. 4. DchGLayHitToDigi (in DchReco) This actually creates the list of digis. 5. DchTrackFinder (in DchReco) Perform trackfinding. Description of DchGLayHitToDigi ------------------------------- The following is taken from DchReco/V00-00-04/DchGLayHitToDigi.README The module DchGLayHitToDigi in the DchReco package takes a list of GEANT layer hits, DchGLayHit, and simulates the drift chamber digis and outputs 2 lists 1) A list (DchWireLayerAList) of 40 wire layers (DchWireLayer) which are containers for a flat array of wire-by-wire pointers to digis where the pointer is NULL if there is no hit on that wire. Digis in this form are convenient for pattern recognition and DchWireLayer is also used to combine and merge digis when 2 or more occur on the the same wire. 2) A list (DchDigiAList) of digis (DchDigi) which owns the digis and manages their memory. DchWireLayer is a bookkeeping auxiliary class only. These 2 lists of hits are produced using the following procedure. A loop over all hits in DchGLayHit is performed. For each hit in a layer a track segment, DchTrkSeg, is produced. Then a loop over all possible wires that this DchTrkSeg could have hit is performed. A tentative digi is constructed assuming circular cell (DchWire), with Landau fluctuations (DchWire), ionization statistic effects (DchWire) and distance to time relation (DchDriftCalib.) If the digi construction fails the tentative digi is deleted otherwise it is added to DchWireLayer for the appropriate layer. The DchWireLayer function, insert, arbitrates between hits that land on the same wire. For the time being only the hit with the shortest time determines the drift time, but all hits contribute to dE/dx and a list of contributing geant hits is maintained for possible later processing. After the geant hits have been processed a second loop over drift chamber layers does: 1) The detector smearing for each wire. There is full knowledge of what hits contributed on each wire so effects like time slewing multiple hits can be - in principle (they are not actually) - accounted for. 2) The resulting digis are appended to the list of digis by the DchWireLayer function, addto. Future Plans for DCH digitization --------------------------------- The current drift chamber digitization code represents a package of software that converts simulated hit data and energy loss in the drift chamber to a time of a hit on a wire and an energy loss within a cell. The ultimate goal is to have the simulation/digitization code produce information that will have the same format as the data produced by the data acquisition system. To this end interest from the Montreal group has been expressed to work on digitization. They can bring expertise on the electronics to this area. The output of the digitization is the input to the reconstruction. The digitization code acts as a bridge between the simulation and reconstruction and as such the digitization code must also understand what information is used and how it is used in the reconstruction. The digitization code should make as much use as possible of the new geometry, DchSimGeom, and cell hit structure. Discussions are in progress on how best to access the geometry information. For example this new cell structure will use a hexagonal cell structure instead of the circular cell structure within the digitization now implemented. The design should cleanly and transparently be able to access the geometry information in the simulation for use in the digitization. Once the digitizations use the new geometry and hit structure the old layer and wire hit structures should be removed. Of course making sure that no other package is using them. There will be a transition period when all 3 hit structures will be present. This will be useful as a check of the new structure. The magnetic field currently used in the dch digitization is hard coded as 1.5T. The digitization needs to access the real magnetic field map. Again work is in progress on how to access the real field map in a straight forward manner. Information from testbeam results on distance to time conversion and resolutions should be incorporated into the digi simulation. Or other simulation packages such as garfield can be used to determine a better distance to time relationship. The digitization simulation must also begin to address calibration and alignment issues. This involves database access to calibration and alignment constants. Another issue that must be addressed is the incorporation of noise into the system. Much work remains to be done in these areas.