]> git.uio.no Git - u/mrichter/AliRoot.git/blob - GEANT321/doc/history/v_314.F
This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / doc / history / v_314.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:19:47  cernlib
6 * Geant
7 *
8 *
9 #include "sys/CERNLIB_machine.h"
10 #include "pilot.h"
11 *CMZ :  3.21/02 29/03/94  15.41.18  by  S.Giani
12 *-- Author :
13 C
14 C
15 C           ************************************************
16 C           *                                              *
17 C           *         G E A N T  version  3.14.16          *
18 C           *                                              *
19 C           *           Official  R e l e a s e            *
20 C           *       ==============================         *
21 C           *                                              *
22 C           ************************************************
23 C
24 C
25 C
26 *    The version 3.14 of GEANT is being released. Apart from a few
27 *    features, reported below, the new version is compatible with the
28 *    previous version 3.13. Substantial developments have taken place,
29 *    in particular in the physics and tracking areas.
30 *
31 *    In addition to the GEANT team
32 *     (R.Brun:CERN/CN-AS,  F.Bruyant:CERN/ECP-PI and M.Maire:LAPP),
33 *    many people have contributed to the new version, especially
34 *     - a large number of users of the previous version who reported
35 *       their experience, found bugs and suggested several improvements,
36 *     - the guinea-pigs of 3.14 who gave us essential comments during
37 *       our attempt to automatize computation of the tracking parameters.
38 *                 ** A.Givernaud (UA1), F.Nessi, V.Vercesi(UA2/LHC)
39 *                 ** The LHC proto collaborations
40 *     - H.Fesefelt(Aachen) has provided a new version of the GHEISHA
41 *       package with help from N.Van Eijndhoven (CERN/CN-AS)
42 *     - H.J Trost from ANL has reported problems in the muon-nuclear
43 *       interactions routines and provided the relevant corrections.
44 *     - P.Pedroni (Pavia) has implemented a new interface called
45 *       by GHEISHA for low-energy hadronic processes (See NUCRIN below).
46 *     - the contributors to the geometry package:
47 *            for the introduction of new shapes
48 *                ELTU by A.Solano (ZEUS)
49 *                HYPE by M.Corden (ALEPH and SSC)
50 *                CTUB by A.McPherson (CEBAF)
51 *            for systematic investigation of problems with the old shapes
52 *                R.Nierhaus (CERN/CN-AS)
53 *     - S.Egli (H1) has proposed an automatic optimisation for the
54 *       geometry at initialisation time.
55 *     - G.Lynch from Berkeley has investigated the multiple-scattering
56 *       various strategies and proposed a new algorithm.
57 *     - A.Rotondi and P.Montagna have proposed a new technique for the
58 *       fast generation of Vavilov distribution.
59 *       K.S.Koelbig (CERN/CN-AS) has implemented a new routine GVAVIV
60 *       based on their work.
61 *     - J.Salt (CERN/CN-AS) has implemented the graphics interface to the
62 *       CG package from Serpukhov with the help of E.Chernaev.
63 *     - The IBM team (C.Guerin, M.Roelisberger) have investigated how
64 *       to speed-up the program. Their work has been coordinated
65 *       by F.Carminati (CERN/CN-AS).
66 *
67 *     GEANT3.14 documentation
68 *     =======================
69 *      The printing of a new manual describing the new version is
70 *      scheduled for the end of this year. The CERN Program Library
71 *      will not distribute anymore the old document (version 3.11).
72 *
73 *
74 *     Important notice to GEANT users
75 *     ===============================
76 *
77 *      Following the reorganisation of the CERN research divisions
78 *      in July 90, R.Brun is now in charge of the Application Software
79 *      group (AS) in the CN division and F.Bruyant is in charge of the
80 *      Production support and computers Infrastructure group (PI) in the
81 *      ECP division. They will nevertheless continue, with M.Maire, to be
82 *      actively involved in the development of GEANT. Federico Carminati
83 *      is the coordinator of a simulation software unit in the CN-AS group.
84 *      Users are strongly recommended to address their questions,etc
85 *      directly to him (email: FCA@CERNVM.CERN.CH on BITNET). In particular,
86 *      feedback from users making comparisons with real data will be most
87 *      appreciated.
88 *
89 *      The simulation team is preparing the ground for the next version
90 *      of GEANT. In view of the proposed new accelerators, the following
91 *      items are considered with high priority:
92 *        - Parametrization techniques. A survey of the various methods
93 *          used in current experiments has been done and the implementation
94 *          of a new algorithm based on the GSCAN geometry + GFLASH (H1)
95 *          is in progress.
96 *        - Improvements in the geometry package. General shape definition,
97 *          surface based algorithms.
98 *        - Detector data structure and data base. Interface with CAD systems.
99 *        - Parallelism (at event level and below).
100 *
101 *
102 *......................................................................
103 *
104 * *** Compatibility with version 3.13
105 *     ===============================
106 *     - If COMMON blocks GCMULO or/and GCJLOC were included in the
107 *       user code, the new GEANT sequences GCMULO and GCJLOC must be
108 *       inserted and the code recompiled.
109 *     - Initialisation data structures saved with the previous versions
110 *       cannot be read by the new version, because the binning for
111 *       the cross-sections and energy loss tables has been changed.
112 *     - The GCPHYS variables SOLOSS,STLOSS,SOMULS,STMULS are no more
113 *       defined (see comments below)
114 *     - The GCTRAK variable IDECAD is replaced by IGAUTO
115 *
116 *
117 * *** MAIN CHANGES IN THE TRACKING PROCEDURES
118 *     =======================================
119 *
120 *    The tracking control routines GTGAMA,GTELEC,GTHADR,GTNEUT,GTMUON
121 *    have been largely rewritten to reflect the changes to the energy-loss
122 *    and multiple scattering processes.
123 *
124 *
125 * *** The ENERGY RANGE of the cross section and energy loss tables can
126 *     be fixed by the user with the new data card :
127 *             'ERANG'   EKMIN  EKMAX    NKBIN
128 *     which defines nkbin bins from Ekmin to Ekmax in a logarithmic scale.
129 *     The default is, as before, 90 bins from 10 Kev to 10 Tev but in
130 *     logarithmic scale. NKBIN must be 50<NKBIN<200.
131 *
132 *     WARNING 1 : as a consequence the common GCMULO has been changed.
133 *     User applications referencing GCMULO must be recompiled.
134 *
135 *     WARNING 2 : as a consequence the ZEBRA data structure JMATE
136 *     contains now cross-sections and energy loss tables based on the new
137 *     energy range. Structures saved by previous versions of GEANT cannot
138 *     be used by the new version.
139 *
140 *     WARNING 3 : changing the 'ERANG' data card arguments requires to
141 *     recreate the initialisation data structures.
142 *
143 *     WARNING 4 : for Hadrons, GEANT tabulates the energy-loss tables
144 *     for a proton-kinetic energy equivalent. As a consequence, the value
145 *     given for EKMAX must be at least 7 times the maximum kinetic energy
146 *     of particles to be tracked (ratio proton/pion mass).
147 *     For LHC/SSC simulations, we recommend ERANG 1.E-5 1.E+5 100
148 *     The recommended value for NKBIN is 10*LOG10(EKMAX/EKBIN).
149 *
150 *
151 * *** The MEAN STOPPING RANGE of a charged particle (STOPMX) is computed and
152 *     tabulated as a function of the kinetic energy (routine GRANGI).
153 *     There are tables for electrons, muons, protons. The others particles are
154 *     scaled from the proton table.
155 *     As for the other tabulated quantities, the Stopping Range can be displayed
156 *     with the routines GPLMAT and/or GPRMAT (keyword 'RANG').
157 *
158 *     During the tracking, the MEAN ENERGY LOSS over a given step (DEMEAN) is
159 *     computed from the difference of the stopping range before and after the
160 *     transport of the particle.
161 *
162 *
163 * *** ENERGY LOSS STRAGGLING : when ILOSS=2 (default) the Landau/Vavilov
164 *     fluctuations are applied to the mean energy loss over the step.
165 *     The VAVILOV distribution has been improved so that the fluctuations are
166 *     less dependent on the step size.
167 *
168 *     ILOSS=1 is now equivalent to ILOSS=3 i.e. restricted fluctuations are
169 *     applied together with delta-ray production (see routine GLANDZ).
170 *     The user who wish to inhibit the energy fluctuations must set ILOSS=4
171 *     The latter option has been kept for debug purpose only.
172 *
173 *
174 * *** MULTIPLE SCATTERING.  The calculation of the step size (SMULS) is new :
175 *     Bethe criterion taken into account, limitation of step size at low
176 *     energy (see routine GMULOF).
177 *
178 *     New routine GMULTS : depending on the step size the Moliere distribution,
179 *     or its gaussian approximation, is used.
180 *
181 *     The routines GMOLI, GMOLIE (ex Gmol), GMOLS have been slightly modified.
182 *
183 *     In the routine GMGAUS (ex Gmuls), we use a new sigma (G. Lynch LBL-28165).
184 *     The lateral displacement has been deleted. As a consequence the tracking
185 *     parameter DMAXMS is not used anymore for the control of the step size.
186 *     DMAXMS.LE.0 forces no multiple scattering at all in the medium.
187 *
188 *     New routine GMCOUL : for very small steps single Coulomb scatters are
189 *     generated instead of Moliere distribution.
190 *
191 *     By default (IMULS=1) GEANT will select automatically the algorithm to use
192 *     (GMCOUL, GMOLIE or GMGAUS).
193 *
194 *
195 * *** NEW ROUTINE GMULOF. The step for energy loss (SLOSS), the step for
196 *     multiple scattering (SMULS), and the step  for the curvature in a magnetic
197 *     field (SFIELD) can be precomputed and tabulated, at least for electrons
198 *     and muons. The routine GMULOF, called from GPHYSI, tabulates:
199 *                  SMULOF = MIN ( SMULS , SLOSS , SFIELD )
200 *     which is the effective step due to the "continuous processes".
201 *     SMULOF can be displayed with the routines GPLMAT/GPRMAT, keyword 'STEP'.
202 *
203 *
204 * *** Multiple scattering and energy loss computations are systematically
205 *     applied at each step during tracking. Therefore the GCPHYS variables
206 *     SOLOSS, STLOSS, SOMULS, STMULS are OBSOLETE. However, for backward
207 *     compatibility with 3.13, STLOSS is set equal to STEP.
208 *
209 *
210 * *** TRACKING PARAMETERS STMIN and DEEMAX.
211 *
212 *     The meaning of STMIN is the following : at low energy (below 1 Mev) the
213 *     multiple scattering condition can induce a very small step SMULOF.
214 *     On the other hand the Stopping Range is also small. Therefore, if the
215 *     Stopping Range is smaller than STMIN, the constraint SMULOF is ignored.
216 *     The exact condition is :
217 *               IF (SMULOF.LE.STMIN)   SMULOF = MIN ( STOPMX, STMIN )
218 *     STMIN is no more than an accelerator factor for the stopping particles.
219 *
220 *     The DEEMAX parameter remains the main tracking parameter. It governs the
221 *     precision of the tracking by limitating SMULOF.
222 *
223 *
224 * *** AUTOMATIC COMPUTATION of STMIN and DEEMAX.
225 *
226 *     By default Geant3.14 overwrites STMIN and DEEMAX. The STMIN default value
227 *     correspond to a stopping range of 200 Kev above CUTELE.
228 *     The default for DEEMAX follows the following algorithm:
229 *       - For non-sensitive volumes (ISVOL=0) DEEMAX is set to 0.25 for
230 *         materials with a radiation length x0<2cm
231 *         and DEEMAX=0.25-0.2/sqrt(x0) for other materials.
232 *       - For sensitive volumes (ISVOL.NE.0) DEEMAX=0.2/sqrt(x0)
233 *     These values have been tuned empirically on a variety of setups.
234 *
235 *     A new data card 'AUTO' has been implemented :
236 *     Setting 'AUTO' 1 is equivalent to NO data card, i.e. automatic computation
237 *     Setting 'AUTO' 0 : NO automatic computation, EXCEPT if STMIN and/or DEEMAX
238 *     has been given a negative value by the user.
239 *
240 *     WARNING : The default values above have been found reasonable for any kind
241 *     of medium. The unexperienced user is invited to start with automatic
242 *     computation. Please check the actual parameters by calling GPRINT ('TMED'
243 *     after the call to GPHYSI.
244 *
245 *     We STRONGLY recommend to always run in AUTO mode (default). The AUTO
246 *     mode makes GEANT a predictive tool if all parameters are automatically
247 *     computed by the system as opposed to tuning data and Monte Carlo
248 *     via the tracking parameters.
249 *
250 * *** GPHYSI
251 *     The routine GPHYSI has been improved to include additional protections
252 *     and to take into account the changes due to the new energy-range tables.
253 *     It must be noted that is mandatory to always call GPHYSI at the end
254 *     of the initialisation phase even when initialisation data structures
255 *     are read from a file.
256 *
257 * *** GPIONS
258 *     This new routine (a complement to GPART) defines a subset of the most
259 *     common "stable" elements in the nature.
260 *     GPIONS can be called at the initialisation stage after GPART. It creates
261 *     particles with GEANT identifiers 61 to 112.
262 *     GPIONS has been written for the heavy ions experiments in view of an
263 *     interface with the program FRITIOF which is being developed by
264 *     P.Gorodetzky (Strasbourg).
265 *
266 * *** GDEBUG
267 *     This new user callable routine from GUSTEP will take the following actions
268 *     if the flag IDEBUG=1:
269 *
270 *     IF(IDEBUG.NE.0) THEN
271 *       IF((ISWIT(2).EQ.1).OR.(ISWIT(3).EQ.1)) CALL GSXYZ ! store point in JXYZ
272 *       IF (ISWIT(2).EQ.2) CALL GPCXYZ ! step by step printed debug
273 *       IF (ISWIT(1).EQ.2) CALL GPGKIN ! list of particles generated during step
274 *       IF (ISWIT(2).EQ.3) THEN
275 *          IF(ISWIT(4).EQ.3.AND.CHARGE.EQ.0.)RETURN
276 *          CALL GDCXYZ    ! interactive drawing of trajectories
277 *        ENDIF
278 *     ENDIF
279 *
280 * *** New user callable routine GBIRK to take into account BIRK's factors.
281 *     ====================================================================
282 *    This new routine is callable from the GUSTEP routine. EDEP=DESTEP
283 *    if no Birk factors have been given via GSTPAR.
284 *      SUBROUTINE GBIRK(EDEP)
285 *
286 *     *** apply BIRK's saturation law to energy deposition ***
287 *     *** only organic scintillators implemented in this version MODEL=1
288 *
289 *     Note : the material is assumed ideal, which means that impurities
290 *            and aging effects are not taken into account
291 *
292 *     algorithm : edep = destep / (1. + RKB*dedx + C*(dedx)**2)
293 *
294 *     the values of RKB and C can be entered via :
295 *
296 *     call gstpar(imate,'BIRK1',value) to set the model (value= 1. or 2.)
297 *     call gstpar(imate,'BIRK2',value) to set RKB
298 *     call gstpar(imate,'BIRK3',value) to set C
299 *
300 *     the basic units of the coefficient are g/(Mev*cm**2)
301 *     because the de/dx is obtained in Mev/cm
302 *
303 *     exp. values from NIM 80 (1970) 239-244 :
304 *
305 *     RKB = 0.013  g/mev*cm**2  and  C = 9.6e-6  g**2/(Mev**2)(cm**4)
306 *
307 *
308 *
309 *
310 * *** GEANH news
311 *
312 *
313 *     a serious bug in GMUNUI has been reported by
314 *     Hans-Jochen Trost from ANL:
315 *       : convert millibarns to cm**2 by a factor 10**-27
316 *         and obtain therefore NA * 1 millibarn/cm**2 a
317 *         factor 10 smaller than in the code. Hans has verified
318 *         that the original code gives to high an energy loss
319 *         due to muon-nuclear interactions on iron by a fac-
320 *         tor of about 5-10, comparing to W.Lohmann et al.,
321 *         CERN 85-03 whose predictions have recently been
322 *         verified by CCFR data up to 1.2 TeV.
323 *     Hans has also introduced protections in GMUNU.
324 *
325 *    A new version of the GHEISHA package has been introduced by
326 *     H.Fesefelt (contact him for details).
327 *
328 *     The TATINA package is considered obsolete. Default routines
329 *     GUPHAD/GUHADR call GHEISHA. We are planning to remove the code
330 *     of TATINA in the coming versions of GEANT.
331 *
332 *
333 * *** NUCRIN: A new model for hadronic showers at low energies
334 *     ========================================================
335 *
336 *
337 *    NUCRIN (see ref.1) simulates hadron-nucleus ( A > 3) inelastic
338 *    interactions from a few MeV/c up to about 4.5 GeV/c laboratory
339 *    momentum of the incoming particle.
340 *    NUCRIN is automatically called by the GHEISHA routine GHEISH
341 *    when the flag IHADR=3 (set by data card HADR).
342 *    It is assumed that these reactions  are the superimposition of three
343 *    basic processes:
344 *
345 *    (a) inelastic collision of the projectile hadron (allowed particles are :
346 *        p,pbar,n,nbar,pi0,pi+,pi-,k+,k-,k0,k0bar,lamda0,lamda0bar,sigma+,
347 *        sigma-,sigma0) with a target nucleon in the nucleus.
348 *        This interaction is simulated,taking into account of the nucleon
349 *        Fermi momentum,  using HADRIN (see ref.2) program.
350 *        The corresponding physical model is based on the experimental evidence
351 *        that, in the selected momentum range,the inelastic cross section shows
352 *        the typical threshold and resonance behaviour of meson production:
353 *        the primary hadron-nucleon system is excited to an isobaric state
354 *        which then decays into hadrons or other resonances.
355 *
356 *        If the interacting nucleus is hydrogen,HADRIN can also be used
357 *        in a separate way to simulate hadron-proton reaction.
358 *
359 *    (b) induced intranuclear cascade with resulting proton and neutron
360 *        emission;
361 *
362 *    (c) nuclear evaporation and deexcitation from residual nucleus. At the
363 *        output the total energy available for these processes is given as
364 *        "excitation energy".
365 *
366 *    The mean excitation and cascade energies and the average multiplicities
367 *    of cascade particles are parametrized, according to experimental
368 *    distributions.
369 *    In each event their value are sampled from gaussian distribution:
370 *    if they fall in the permitted kinematical region, energy and types of
371 *    cascade nucleons are calculated and the remaining energy is assigned to
372 *    the incoming particle.
373 *    For hadron-nucleon interactions all relevant kinematic variables are
374 *    Lorentz-trasformed into the target nucleon rest system. If interaction
375 *    kinetic energy is greater than the total available collision energy
376 *    a new Fermi momentun is sampled, otherwise an event is generated with
377 *    HADRIN,in which decays modes of 107 particle and resonances into about
378 *    450 different channels are tabulated and outgoing particle directions
379 *    and momenta are chosen to reproduced experimental momentum transfer
380 *    distributions.
381 *    Final state particles kinematical variables are transformed back into
382 *    laboratory system; reaction and sampled event energies are again compared:
383 *    if their difference is negative, energy is not conserved and generation
384 *    has to be started once more with a new Fermi momentum sampling or if it
385 *    is,on the contrary positive, particle momenta and energies are corrected
386 *    to reach conservation.
387 *    The sampled events conserve the energy, the momentum, the electric and
388 *    baryonic charge  and the strangeness.
389 *    NUCRIN and HADRIN are initialised, by default, before event generation,
390 *    with a call to subroutines HADDEN and CHANWX which estabilish internal
391 *    weight tables and decay channels.
392 *
393 *    -----------------------------------------------------------------------
394 *    (1) K. Hanssgen, J. Ranft , Comp. Phys. Comm. 39, 53 (1986)
395 *    (2) K. Hanssgen, J. Ranft , Comp. Phys. Comm. 39, 37 (1986)
396 *
397 *
398 * *** GEOMETRY PACKAGE: New shapes and many improvements
399 *     ==================================================
400 *
401 *    Automatic optimisation of the geometry structure:
402 *    A new data card OPTI has been introduced (S.Egli H1).
403 *      OPTI -1 : disable optimisation
404 *      OPTI  0 : only volumes GSORDered are optimised (as in 3.13)
405 *      OPTI  1 : volumes GSORDered are optimised along the axis
406 *                specified. All the other volumes are automatically
407 *                optimised along the best axis (1 to 7).
408 *      OPTI  2 : All volumes are optimised along the best axis.
409 *                Volumes for which GSORD was called are also optimised.
410 *      The default value for OPTI is 0.
411 *      In case OPTI >0, the result of the optimisation is printed.
412 *      The automatic optimisation is done at initialisation time by
413 *      a new routine GGORDQ called by GGCLOS.
414 *
415 *
416 *    Most of the geometry routines have been revisited and consolidated.
417 *    The following new shapes are available.
418 *
419 *      'ELTU'    is a cylinder with an elliptical section.
420 *                It has three parameters: the ellipse semi-axis in X,
421 *                the ellipse semi-axis in Y and the half length in Z.
422 *                Given the equation of the conical curve:
423 *                     X**2/A**2 + Y**2/B**2 = 1,
424 *                describing the volume,then:       PAR(1) = A
425 *                                                  PAR(2) = B
426 *                                                  PAR(3) = DZ
427 *                ELTU is not divisible.
428 *
429 *      'HYPE'    is a hyperbolic tube, ie the inner and outer surfaces
430 *                are hyperboloids, as would be formed by a system of
431 *                cylindrical wires which were then rotated
432 *                tangentially about their centres.  The 4 parameters
433 *                are the inner and outer radii, the half length in z,
434 *                and the "stereo angle" theta in degrees, such that
435 *                the hyperbolic surfaces are given by
436 *                r**2 = (z*tan(theta))**2 + (r at z=0)**2
437 *
438 *      'CTUB'   (for cut tube) is a TUBS whose end planes are not
439 *               perpendicular to the z axis. It has 11 parameters :
440 *               the 5 of the TUBS shape plus the components of the normal
441 *               to the end plane at the lower z (LXL,LYL,LZL) and
442 *               those at the higher z (LXH,LYH,LZH).
443 *               DZ means the half length in z for x = y = 0
444 *
445 *
446 *     Bugs have been fixed in the routines GNPCON, GNPGON and GNOPGO.
447 *
448 *     The new version of GNOTRP requires an extended parameter array.
449 *     In addition to the 11 specified parameters (of which 4 are
450 *     modified in subroutine GSVOLU or GSPOSP), the coefficients of
451 *     the implicit normalized plane equation for the 6 surfaces
452 *     of the hexahedron are stored.
453 *
454 *     IMPORTANT NOTE concerning the TRAP shape
455 *     ========================================
456 *     The Geant documentation describes the Geant shape TRAP as follows:
457 *
458 *          TRAP is a general trapezoid, i.e. one for which the
459 *          faces perpendicular to z are trapezia and their
460 *          centres are not at the same x, y. It has 11
461 *          parameters: Dz the half length in z, Th & Phi the
462 *          polar angles from the centre of the face at z=-Dz
463 *          to that at z=+Dz, H1 the half length in y  at
464 *          z=-Dz, LB1 the half length in x at z=-Dz and y=low
465 *          edge, LH1 the half length in x at z=-Dz and y=
466 *          high edge, Th1 the angle w.r.t. the y axis from
467 *          the centre of the low y  edge to the centre of the
468 *          high y edge, and H2, LB2, LH2, Th2 the
469 *          corresponding quantities to the 1s but at z=+Dz.
470 *
471 *     This seems to describe a general hexahedron with 3 constraints:
472 *     2 constraints follow from the fact that two faces are "trapezia".
473 *     (twice 2 edges parallel).
474 *     The 3rd constaint is that two faces are parallel, namely the "trapezia"
475 *     faces are both perpendicular to the Z-axis.
476 *     We will assume that shape TRAP is a hexahedron with 3
477 *     constraints and direct our attention to the degrees of freedom
478 *     of such a shape.
479 *     The shape has 8 vertices and therefore 24 coordinates.
480 *     If we first consider a volume with 8 vertices and 6 surfaces,
481 *     but do not make the assumption that the surfaces are plane,
482 *     we see that this shape has 18 degrees of freedom.
483 *     We loose 3 coordinates because of the translational invariance
484 *     of the shape and 3 coordinates because of its rotational
485 *     invariance.
486 *     If we now assume that the shape is a hexahedron, that is
487 *     bounded by plane surfaces, we have 6 constraints, one for each
488 *     surface, and our hexahedron has 12 degrees of freedom.
489 *     Considering the 3 constraints mentioned in the beginning,
490 *     we conclude that our shape trapezohedron has 9 degrees of freedom.
491 *     It is however described by 11 parameters.
492 *     Therefore we must either drop our assumption that shape TRAP is
493 *     a hexahedron, that is bounded by parallel surfaces, or we must
494 *     request that the user specifies the 11 parameters with certain
495 *     constraints.
496 *
497 *     To check that the user respected the constraints, we check the
498 *     coplanarity of the faces during the specification phase of shape TRAP.
499 *     We know the vertex coordinates, and we have the indices to the
500 *     vertices for each face.
501 *     Assuming that a face is tetrahedron, we compute its volume.
502 *     We divide by the surface of the base triangle, and get a measure
503 *     for the coplanarity of the face, which is actually a distance.
504 *     A warning message is printed in case of no-coplanar faces.
505 *
506 *    GSORD problem
507 *    =============
508 *    A bug has been found in GTNEXT (and alike) which is induced by a bug of
509 *    logic in GSORD/GGORD. User calls to GSORD, with ordering axis 4 (Rxy)
510 *    or 5 (Rxyz), may cause problems when the ordered contents are such that
511 *    one can jump from a given content to another one without crossing a
512 *    content which, along the given axis, occupies a position in between the
513 *    start and the end contents : e.g. coaxial TUBES with different Z-lengths
514 *    should not in general be ordered by GSORD along the axis IAX=4 (Rxy).
515 *    However, part of the information provided by such calls can still be
516 *    used, in the static context of GTMEDI for instance. Therefore, the
517 *    following convention has been introduced: If the user is sure that the
518 *    contents are positioned in such a way that the anomaly mentioned above
519 *    cannot occur, the call to GSORD can be modified by using IAX=14 (instead
520 *    of 4), or 15 (instead of 5), in which case the ordering techniques will
521 *    also be used in the dynamic context of GTNEXT. In case of doubt, the
522 *    user has better to keep the old code, with IAX=4 or 5.
523 *
524 * *** DRAWING PACKAGE: Interface to the CG package
525 *     ============================================
526 *
527 *    An interface to the CG (Combinatorial Geometry) package written
528 *    at Serpukhov by E.Cernaev et al is now available. The new package
529 *    is automatically called if the option 'HIDE' is selected. eg.
530 *
531 *         CALL GDOPT('HIDE','ON')  in a Fortran program
532 *         DOPT HIDE ON/OFF  in the interactive version
533 *
534 *    This new facility includes a hidden line and surface algorithm
535 *    which permits nice 3-D views of a detector.
536 *    In the frequent case of hermetic 4 PI detectors, a facility
537 *    to remove a box (The Cutting BOX) is also implemented. The Cutting BOX
538 *    specifies a region of the detector which must be Cut to see inside.
539 *    A new interactive command CBOX is available to specify the box limits.
540 *
541 *    The CG system is part of the GEANG file (See Patches CGCDES,CGPACK)
542 *    To activate the CG package, +USE,CG,*GEANG.
543 *
544 *
545 * *** The SCAN geometry
546 *     =================
547 *   A new Patch,GSCAN has been introduced on a provisional basis in the GEANG
548 *   file. See discussion about SCAN below.
549 *
550 * *** GUPARA: Parametrization interface
551 *     =================================
552 *   A new FFREAD data cards PCUT can be used to set parametrization cuts.
553 *   The first argument of the PCUT card is a integer flag which turns or
554 *   on off the parametrization mechanism. If the parametrization is turned on
555 *   and a particle falls below one of the 5 cuts specified by the PCUT card
556 *   (similar in kind to the cuts specified by the card CUTS), then the
557 *   routine GUPARA is called and tracking of the particle is abandoned. This
558 *   mechanism is provided for applying parametrization schemes which
559 *   replace a particle by a parametrized shower when it falls below
560 *   a certain threshold.
561 *
562 *.............................................................................
563 *
564 *
565 * *** GXINT: Interactive version
566 *     ==========================
567 *
568 *      New menu FORTRAN: CALL,FILE,CLOSE,FORTRAN
569 *      New menu HISTOGRAM: FILE,LIST,PLOT,DELETE,LEGO,HRIN,HROUT,PUT,GET
570 *                          ZONE,SET,OPTION,NULL
571 *      New menu PICTURE: FILE,LIST,DELETE,SCRATCH,PLOT,RENAME,IZOUT,IZIN,IGSET
572 *      New menu SCAN: PHI,THETA,SLIST,VERTEX,SFACTORS,LSCAN,HSCAN
573 *      New menu PHYSICS: ANNI,BREM,COMP,DCAY,DRAY,HADR,LOSS,MULS,MUNU,PAIR,
574 *                        PFIS,PHOT,RAYL,CUTS,PHYSI
575 *
576 *      The menus FORTRAN,HISTOGRAM and PICTURE are subsets of the similar
577 *      menus in the PAW system.
578 *
579 *      The menu PHYSICS gives the possibility to modify the run conditions.
580 *      In case physics conditions are changed (LOSS,DRAY,MULS,CUTS) it is
581 *      mandatory to call the command PHYSI (which calls GPHYSI) to recompute
582 *      the cross-section and/or energy loss tables.
583 *
584 *      The menu FORTRAN is similar to the FORTRAN menu of PAW. It contains
585 *      in addition a new very important command FORTRAN which gives the
586 *      possibility to describe the geometry (UGEOM) in a Fortran routine
587 *      that can be edited interactively with the local editor and also
588 *      executed interactively under the control of the Fortran interpreter
589 *      COMIS.
590 *
591 * ***  GEANT >FORTRAN  FNAME
592 *
593 *      The routines in the file FNAME will be compiled by COMIS.
594 *      If routines with names: UGEOM,GUKINE,GUOUT,UGLAST are found,
595 *      then they will be automatically called by GXINT instead of
596 *      the routines with the same names compiled with the standard
597 *      Fortran compiler and linked with the application.
598 *      The user callable routines from the GEANT library as well as
599 *      routines from PACKLIB (HBOOK,HPLOT,HIGZ,ZEBRA) may be called
600 *      from these user routines. All GEANT common blocks may be
601 *      referenced.
602 *      In case where the routine UGEOM is called several times,
603 *      it is important to DROP all the initialisation data structures
604 *      JVOLUM,JMATE,JTMED,etc already in memory by using the routine GIDROP.
605 *
606 *       Example of an interactive session where the routine UGEOM is modified:
607 *
608 *         GEANT > Edit ugeom.for
609 *         GEANT > Fortran ugeom.for
610 *         GEANT > Call GIDROP
611 *         GEANT > Call UGEOM
612 *         GEANT > Dtree
613 *         GEANT > Edit ugeom.for
614 *         GEANT > Fortran ugeom.for
615 *         GEANT > Call GIDROP
616 *         GEANT > Call UGEOM
617 *         GEANT > Dtree
618 *
619 *      If FNAME='-', calls to user routines is reset and standard
620 *      routines called instead.
621 *
622 * *** Interface to CG
623 *    The command DOPT has a new option HIDE (DOPT HIDE ON/OFF)
624 *    Type DOPT without parameters to get the list of all currently
625 *    available options.
626 *    When this option is ON, the subsequent graphics commands DCUT/DRAW
627 *    will invoke the CG system for hidden line/surface removal.
628 *    This algorithm requires a lot of memory and time. It is recommended
629 *    to set the visibility attributes (SATT SEEN 0) for many of the
630 *    volumes in case the command aborts for lake of memory or time.
631 *    This option can also be used with the view banks mechanism (DOPEN)
632 *
633 *    New command CBOX to specify the boundaries of the cutting box.
634 *
635 * *** PLMAT
636 *
637 *   The existing command PLMAT offers the new possibility to plot
638 *    various physics parameters (cross-sections, energy-loss tables,etc)
639 *   in graphics format (via HPLOT) if MECAN=ALLG
640 *   The Keywords 'STEP' or 'RANG' may also be specified to produce
641 *   an alphanumeric output of the step-size and energy-range tables.
642 *
643 *        PLMAT  IMATE IPART MECAN [ IDM ]
644 *
645 *       IMATE      I 'Material number'
646 *       IPART      I 'Particle number'
647 *       MECAN      C 'Mechanism'
648 *       IDM        I 'ID mode option' D=0
649 *
650 *
651 * *** New menu SCAN. The SCAN geometry
652 *     ================================
653 *
654 *    This new menu contains various commands for an interactive interface
655 *    to the SCAN geometry
656 *    The SCAN geometry algorithm has been designed as a tool to improve
657 *    the tracking speed. This new facility still requires substantial
658 *    developments in view of the new parametrisation algorithms which
659 *    are developed in collaboration between the CN/AS group and the
660 *    LEP/HERA/LHC/SSC and other interested groups.
661 *    The SCAN facility is being introduced in the version 3.14 on a trial
662 *    basis to familiarise potential users with the concept.
663 *    The SCAN geometry data structure JSCAN is automatically generated
664 *    either by calling the GSCAN routine in the PATCH,GSCAN of GEANG
665 *    or interactively by using the commands in the new menu SCAN.
666 *    Starting from the normal geometry data structure created by GSVOLU,
667 *    GSPOS,GSDVN,etc, the detector may be divided into a simpler geometry
668 *    structure (theta,phi) or (eta,phi).
669 *    Geantinos are tracked starting from a VERTEX position through
670 *    the NPHI,NTETA divisions. For each division, the SCAN procedure
671 *    will insert into the JSCAN data structure the following information
672 *    for every main detector component specified in the SLIST command
673 *    in the spherical R direction:
674 *      Total number of radiation lengths up to entry in each R
675 *      Total number of absorption lenghts
676 *      Detector identifier
677 *    When the interactive command TRIGGER is entered, the number of
678 *    Geantinos specified as parameter will be tracked. In case the
679 *    data structure JSCAN is not empty, the program will automatically
680 *    start with the first PHI,TETA division not yet filled. As the
681 *    number of Geantinos to be tracked can be very large (depending
682 *    on the granularity) this gives the possibility to fill the JSCAN
683 *    data structures in several passes.
684 *
685 *
686 *     New menu SCAN: PHI,TETA,SLIST,VERTEX,SFACTORS,LSCAN,HSCAN
687 *
688 * ==>   /SCAN/PHI  NPHI [ PHIMIN PHIMAX ]
689 *
690 *        NPHI       I 'Number of PHI divisions' D=90
691 *        PHIMIN     R 'Minimum PHI in degrees' D=0
692 *        PHIMAX     R 'Maximum PHI in degrees' D=360
693 *
694 *        To specify number of divisions along PHI.
695 *
696 *
697 * ==>   /SCAN/TETA  NTETA TETMIN TETMAX [ DIVTYP ]
698 *
699 *        NTETA      I 'Number of TETA divisions' D=90
700 *        TETMIN     R 'Minimum value of TETA' D=0
701 *        TETMAX     R 'Maximum value of TETA' D=180
702 *        DIVTYP     I 'Type of TETA division' D=1 R=1:2
703 *
704 *        To specify number of divisions along TETA.
705 *        If DIVTYP=1 divisions in degrees following the THETA angle.
706 *        If DIVTYP=2 divisions in pseudo-rapidity ETA.
707 *
708 *
709 * ==>   /SCAN/SLIST  LIST
710 *
711 *        LIST       C 'List of master volumes'
712 *
713 *        Only boundary crossings of volumes given in LIST
714 *        will be seen in the SCAN geometry.
715 *
716 *
717 * ==>   /SCAN/VERTEX  VX VY VZ
718 *
719 *        VX         R 'Scan vertex origin' D=0
720 *        VY         R 'Scan vertex origin' D=0
721 *        VZ         R 'Scan vertex origin' D=0
722 *
723 *        All Geantinos tracked will start from position VX,VY,VZ.
724 *
725 *
726 * ==>   /SCAN/SFACTORS  FACTX0 FACTL FACTR
727 *
728 *        FACTX0     R 'Scale factor for SX0' D=100
729 *        FACTL      R 'Scale factor for SL' D=1000
730 *        FACTR      R 'Scale factor for R' D=100
731 *
732 *        Set scale factors for SX0,SL and R. The given scale factors must be
733 *        such that:
734 *
735 *          SX0*FACTX0 < 2**15-1 (32767)
736 *          SL*FACTL   < 2**10-1 (1023)
737 *          SR*FACTR   < 2**17-1 (131071)
738 *
739 *
740 * ==>   /SCAN/LSCAN  ID [ VOLUME CHOPT ]
741 *
742 *        ID         I 'Lego plot identifier' D=2000
743 *        VOLUME     C 'Volume name' D='XXXX'
744 *        CHOPT      C 'List of options' D='OPX' R=' ,O,P,I,X,L'
745 *
746 *        Generates and plot a table of physics quantities such as
747 *        the total number of radiation lengths or interaction lengths
748 *        in function of the SCAN parameters TETA,PHI.
749 *
750 *          CHOPT='O' table is generated at Exit  of VOLUME.
751 *          CHOPT='I' table is generated at Entry of VOLUME.
752 *          CHOPT='X' radiation lengths
753 *          CHOPT='L' Interaction lengths
754 *          CHOPT='P' Plot the table
755 *
756 *        If VOLUME='XXXX' Mother volume is used.
757 *
758 *
759 * ==>   /SCAN/HSCAN  IDPHI [ VOLUME CHOPT ]
760 *
761 *        IDPHI      I 'Histogram/phi identifier' D=1000
762 *        VOLUME     C 'Volume name' D='XXXX'
763 *        CHOPT      C 'List of options' D='OPX' R=' ,O,P,I,X,L'
764 *
765 *        Generates and plot an histogram of physics quantities such as
766 *        the total number of radiation lengths or interaction lengths
767 *        in function of the SCAN parameter TETA for a given value of PHI.
768 *
769 *          CHOPT='O' histogram is generated at Exit  of VOLUME.
770 *          CHOPT='I' histogram is generated at Entry of VOLUME.
771 *          CHOPT='X' radiation lengths
772 *          CHOPT='L' Interaction lengths
773 *          CHOPT='P' Plot the histogram
774 *
775 *        If VOLUME='XXXX' Mother volume is used.
776 *        The histogram identifier IDPHI is used to also identify which
777 *        PHI division to plot. IPHI=MOD(IDPHI,1000).
778 *        If IPHI=0, then all PHI divisions are generated (not plotted)
779 *        with histogram identifiers IDPHI+PHI division number.
780 *
781 *
782 *
783 *
784 *
785 *
786 * ***   New commands FILE,REND,MDIR,CDIR,IN,OUT in the RZ menu.
787 *       =======================================================
788 *
789 *
790 * ==>   RZ/FILE  LUN FNAME [ CHOPT ]
791 *
792 *        LUN        I 'Logical unit number'
793 *        FNAME      C 'File name'
794 *        CHOPT      C 'Options' D=' ' R=' ,U,N,I,O'
795 *
796 *        Open a GEANT/RZ file. Call GRFILE (See below).
797 *
798 *         CHOPT=' ' readonly mode
799 *         CHOPT='U' update mode
800 *         CHOPT='N' create new file
801 *         CHOPT='I' Read all structures from existing file
802 *         CHOPT='O' Write all structures on file
803 *
804 *
805 * ==>   RZ/OUT  OBJECT [ IDVERS ]
806 *
807 *        OBJECT     C 'Structure name'
808 *        IDVERS     I 'Version number' D=1
809 *
810 *        Write data structure identified by OBJECT,IDVERS to RZ file.
811 *        Call GROUT (See below)
812 *
813 *          MATE write JMATE structure
814 *          TMED write JTMED structure
815 *          VOLU write JVOLUM structure
816 *          ROTM write JROTM structure
817 *          SETS write JSET  structure
818 *          PART write JPART structure
819 *          SCAN write JSCAN structure
820 *          *    write all structures
821 *
822 * ==>   RZ/IN  OBJECT [ IDVERS ]
823 *
824 *        OBJECT     C 'Structure name'
825 *        IDVERS     I 'Version number' D=1
826 *
827 *        Read data structure identified by OBJECT,IDVERS into memory.
828 *        Call GRIN (See below)
829 *
830 *          MATE read JMATE structure
831 *          TMED read JTMED structure
832 *          VOLU read JVOLUM structure
833 *          ROTM read JROTM structure
834 *          SETS read JSET  structure
835 *          PART read JPART structure
836 *          SCAN read JSCAN structure
837 *          *    read all structures
838 *
839 *
840 *
841 *
842 * *** New routines for direct access I/O in the GIOPA package
843 *     =======================================================
844 *
845 * ==>    SUBROUTINE GRFILE(LUN,CHFILE,CHOPT)
846 *.
847 *.           Routine to open a GEANT/RZ data base.
848 *.
849 *.           LUN logical unit number associated to the file
850 *.
851 *.           CHFILE RZ file name
852 *.
853 *.           CHOPT is a character string which may be
854 *.              'N' To create a new file
855 *.              'U' to open an existing file for update
856 *.              ' ' to open an existing file for read only
857 *.              'Q' The initial allocation (default 1000 records)
858 *.                  is given in IQUEST(10)
859 *.              'I' Read all data structures from file to memory
860 *.              'O' Write all data structures from memory to file
861 *.
862 *.           Note:
863 *.             If options 'I' or 'O' all data structures are read or
864 *.                written from/to file and the file is closed.
865 *.             See routine GRMDIR to create subdirectories
866 *.             See routines GROUT,GRIN to write,read objects
867 *.
868 *.
869 *.
870 *. ==>    SUBROUTINE GROUT(CHOBJ,IDVERS,CHOPT)
871 *.
872 *.           Routine to write GEANT object(s) in the RZ file
873 *.             at the Current Working Directory (See RZCDIR)
874 *.           Input is taken from the data structures in memory
875 *.               (VOLU,ROTM,TMED,MATE,SETS,PART,SCAN)
876 *.
877 *.           CHOBJ  The type of object to be written:
878 *.                  MATE write JMATE structure
879 *.                  TMED write JTMED structure
880 *.                  VOLU write JVOLUM structure
881 *.                  ROTM write JROTM structure
882 *.                  SETS write JSET  structure
883 *.                  PART write JPART structure
884 *.                  SCAN write LSCAN structure
885 *.                  INIT write all initialisation structures
886 *.
887 *.           IDVERS is a positive integer which specifies the version
888 *.               number of the object(s).
889 *.
890 *.           CHOPT List of options (none for the time being)
891 *.
892 *.        Note that if the cross-sections and energy loss tables
893 *.           are available in the data structure JMATE, then they are
894 *.           saved on the data base.
895 *.
896 *.
897 *.        The data structures saved by this routine can be retrieved
898 *.        with the routine GRIN.
899 *.
900 *.        Before calling this routine a RZ data base must have been
901 *.        created using GRFILE.
902 *.        The data base must be closed with RZEND.
903 *.          Ex: if LUN=1 CALL RZEND('LUN1')
904 *.
905 *.        The RZ data base can be transported between different
906 *.        machines in using the ZEBRA RZ utility RZTOFZ.
907 *.
908 *.        The interactive version of GEANT provides facilities
909 *.        to interactively update, create and display objects.
910 *.
911 *.          Example.
912 *.
913 *.          CALL GRFILE(1,'Geometry.dat','N')
914 *.          CALL GROUT('VOLU',1,' ')
915 *.          CALL GROUT('MATE',1,' ')
916 *.          CALL GROUT('TMED',1,' ')
917 *.          CALL GROUT('ROTM',1,' ')
918 *.          CALL GROUT('PART',1,' ')
919 *.          CALL GROUT('SCAN',1,' ')
920 *.          CALL GROUT('SETS',1,' ')
921 *.
922 *.          The same result can be achieved by:
923 *.          CALL GRFILE(1,'Geometry.dat','NO')
924 *.
925 *.
926 *.
927 *. ==>    SUBROUTINE GRIN(CHOBJ,IDVERS,CHOPT)
928 *.
929 *.           Routine to read GEANT object(s) fromin the RZ file
930 *.             at the Current Working Directory (See RZCDIR)
931 *.           The data structures from disk are read in memory
932 *.               (VOLU,ROTM,TMED,MATE,SETS,PART,SCAN)
933 *.
934 *.           CHOBJ  The type of object to be read:
935 *.                  MATE read JMATE structure
936 *.                  TMED read JTMED structure
937 *.                  VOLU read JVOLUM structure
938 *.                  ROTM read JROTM structure
939 *.                  SETS read JSET  structure
940 *.                  PART read JPART structure
941 *.                  SCAN read LSCAN structure
942 *.                  INIT read all initialisation structures
943 *.
944 *.           IDVERS is a positive integer which specifies the version
945 *.               number of the object(s).
946 *.
947 *.           CHOPT List of options (none for the time being)
948 *.
949 *.
950 *.        The RZ data base has been created via GRFILE/GROUT
951 *.
952 *.
953 *.          Example.
954 *.
955 *.          CALL GRFILE(1,'Geometry.dat',' ')
956 *.          CALL GRIN ('VOLU',1,' ')
957 *.          CALL GRIN ('MATE',1,' ')
958 *.          CALL GRIN ('TMED',1,' ')
959 *.          CALL GRIN ('ROTM',1,' ')
960 *.          CALL GRIN ('PART',1,' ')
961 *.          CALL GRIN ('SCAN',1,' ')
962 *.          CALL GRIN ('SETS',1,' ')
963 *.
964 *.          The same result can be achieved by:
965 *.          CALL GRFILE(1,'Geometry.dat','I')
966 *.
967 *.
968 *.
969 * ==>     SUBROUTINE GRMDIR(CHDIR,CHOPT)
970 *.
971 *.
972 *.           Routine to create a subdirectory
973 *.
974 *.           CHDIR Subdirectory name
975 *.
976 *.           CHOPT is a character string which may be
977 *.              ' ' To create a subdirectory
978 *.              'S' To create a subdirectory and set the new
979 *.                  Current Directory to this directory.
980 *.
981 *.
982