+OPTION,MAPASM.
+ASM,23.
+USE,GCDES,TestEm3,T=EXE.
*
+EXE,CRA*.
*
+PAM,11,T=CARD,T=ATTACH.     /cern/pro/src/car/geant321.car
+PATCH,TestEm3.        
+DECK,BLANKDEK.
* 
* TestEm3   2.00 /01  990723 13.00 GEANT EXAMPLES
* 
*             TEST PROGRAM FOR GEANT TEST
*                                                                               
*             Author    M.Maire  *********
*                                   
*     
* History                  
* -------                                            
*
*                                      
+KEEP,CALOR. 
      PARAMETER (MaxAbs = 10)                 
      COMMON/CALOR/materAbs(MaxAbs),thickAbs(MaxAbs),NbAbsor,NbLayer,
     &             calorYZ,Field,thLayer,calorX,worldX,worldYZ
 
*                                                              
*  NbAbsor = Number of Absorbers                    (data card CALOR)
*  NbLayer = Number of Layers                       (data card CALOR) 
*  calorYZ = tranverse size of calor (cm)           (data card CALOR)
*  Field   = magnetic field (tesla)                 (data card CALOR)
*  materAbs(i) = Absorber material number.          (data card MATE )
*  thickAbs(i) = thickness of absorber   (cm)       (data card THICK)
*
*  thLayer = total thickness of a layer (cm)
*  calorX  = total  X-size of the calor (cm)
*  worldX  = total  X-size of the world (cm)
*  worldYZ = total YZ-size of the world (cm)
*                                                                               
+KEEP,CELOSS.                                                                                                                      
      COMMON/CELOSS/ EdepAbs(MaxAbs),TrckAbs(MaxAbs),
     &               sumEdep(MaxAbs),sumTrck(MaxAbs),                                                        
     &               su2Edep(MaxAbs),su2Trck(MaxAbs)
     
+DECK,main,if=batch.
      PROGRAM main
*
*
      PARAMETER (NGBANK=100000, NHBOOK=40000)
      COMMON/GCBANK/Q(NGBANK)
      COMMON/PAWC  /H(NHBOOK)
*
      CALL GZEBRA( NGBANK)
      CALL HLIMIT(-NHBOOK)
*
* *** initialize HIGZ
      CALL HPLINT(0)
*
* *** GEANT initialisation
      CALL UGINIT
*
* *** Start events processing
      CALL GRUN
*
* *** End of RUN
      CALL UGLAST
*
      STOP
      END
+DECK,main,IF=-batch.
      PROGRAM main
*
*     GEANT main program. To link with the MOTIF user interface
*     the routine GPAWPP(NWGEAN,NWPAW) should be called, whereas
*     the routine GPAW(NWGEAN,NWPAW) gives access to the basic
*     graphics version.
*
      PARAMETER (NWGEAN=3000000, NWPAW=1000000)
      COMMON/GCBANK/GEANT(NWGEAN)
      COMMON/PAWC  /PAW  (NWPAW)
*
*
      CALL GPAW (NWGEAN,NWPAW)
*
      END
*
      SUBROUTINE qnext
      END
*
      SUBROUTINE czopen
      END
*
      SUBROUTINE cztcp
      END
*
      SUBROUTINE czclos
      END
*
      SUBROUTINE czputa
      END
+DECK,UGINIT                                                                    
      SUBROUTINE UGINIT                                                         
*                                                                               
*     To initialise GEANT/USER  program and read data cards                     
*
+SEQ,CALOR,CELOSS.                                                                               
+SEQ,GCKINE.
*
      CHARACTER*20 filnam                                                                                                                                                                                                      
*                                                                               
* *** Define the GEANT parameters                                                          
      CALL GINIT 
                                                                          
* *** read data cards
      PRINT *, 'G3 > gives the filename of the data cards to be read:'
      READ (*,'(A)') filnam
      IF (filnam.EQ.' ') filnam = 'testem3.dat'
      OPEN (unit=5,file=filnam,status='unknown',form='formatted')
                                              
*
* *** Calor definition                                                                                                                                                                      
      CALL FFKEY('CALOR',NbAbsor,4,'MIXED')
      CALL FFKEY('MATE' ,materAbs(1),MaxAbs,'INTEGER')
      CALL FFKEY('THICK',thickAbs(1),MaxAbs,'REAL')
*                                                                    
* *** read data cards                                                           
      CALL GFFGO 
                                                                     
*                                                                               
      CALL GZINIT                                                               
      CALL GPART 
*
      CALL GDINIT                                                                     
*                                                                               
* *** Geometry and materials description                                        
      CALL UGEOM                                                                
*                                                                               
* *** Energy loss and cross-sections initialisations                            
      CALL GPHYSI                                                               
*                                                                               
      CALL GPRINT('MATE',0)                                                     
      CALL GPRINT('TMED',0)                                                     
      CALL GPRINT('VOLU',0)                                                     
*                                                                                                                                                              
* *** Define user histograms                                                    
      CALL UHINIT
*
* *** some initialisation
      do k=1,MaxAbs
        sumEdep(k) = 0.
        sumTrck(k) = 0.
        su2Edep(k) = 0.
        su2Trck(k) = 0.
      enddo                                                         
*                                                                               
      END                                                                       
+DECK,UGEOM.                                                                    
      SUBROUTINE UGEOM                                                          
*                                                                               
* *** Define user geometry set up                                               
*                                                                                                                                                                                                                       
*
+SEQ,CALOR.
*
      DIMENSION PAR(3)

      DIMENSION Aair(2),Zair(2),Wair(2)
      DIMENSION ACO2(2),ZCO2(2),WCO2(2)
      DIMENSION AH2O(2),ZH2O(2),WH2O(2)
      DIMENSION AG10(4),ZG10(4),WG10(4)                                        
*
      CHARACTER*4  volnam
      CHARACTER*20 matnam
*      
* *** Air compound parameters          
      DATA Aair/14.01, 16.00/
      DATA Zair/ 7.  ,  8.  /
      DATA Wair/ 0.7 ,  0.3 /
*                                                                               
* *** CO2 compound parameters
      DATA ACO2/12.01, 16.00/
      DATA ZCO2/ 6.  ,  8.  /
      DATA WCO2/ 1.  ,  2.  /
*                                                                               
* *** Water compound parameters
      DATA AH2O/ 1.01, 16.00/
      DATA ZH2O/ 1.  ,  8.  /
      DATA WH2O/ 2.  ,  1.  /
*                                                                               
* *** G10 compound parameters
      DATA AG10/ 1.01, 12.00, 16.00, 28.00/
      DATA ZG10/ 1.  ,  6.  ,  8.  , 14.  /
      DATA WG10/ 3.  ,  3.  ,  2.  ,  1.  /      
      
*
* *** Defines USER perticular materials
      CALL GSMIXT( 1,'Air'      , Aair ,Zair, 1.29E-3, 2   , Wair)
      CALL GSMIXT( 2,'CO2 gas'  , ACO2 ,ZCO2, 27.0E-3,-2   , WCO2)
      CALL GSMATE( 3,'H2 Liquid',  1.01,  1., 0.0708 , 865., 790., 0,0)
      CALL GSMIXT( 4,'Water'    , AH2O ,ZH2O, 1.0    ,-2   , WH2O)
      CALL GSMATE( 5,'Liquid Ar', 39.95, 18., 1.39   , 14.0, 84.0, 0,0)
      CALL GSMATE( 6,'Aluminium', 26.98, 13., 2.7    ,  8.9, 37.2, 0,0)
      CALL GSMATE( 7,'Iron     ', 55.85, 26., 7.87   , 1.76, 17.1, 0,0)
      CALL GSMATE( 8,'Lead     ',207.19, 82., 11.35  , 0.56, 18.5, 0,0)
      CALL GSMATE( 9,'Uranium  ',238.03, 92., 18.95  , 0.32, 12. , 0,0)
      CALL GSMATE(10,'Silicium ', 28.09, 14.,  2.33  , 9.36, 45.5, 0,0)
      CALL GSMATE(11,'Tungsten ',183.85, 74., 19.30  , 0.35,  9.6, 0,0)
      CALL GSMIXT(12,'NemaG10'  , AG10 ,ZG10, 1.7    ,-4   , WG10)      
*                                               
*                                                                               
* *** Defines USER tracking media parameters
      IFIELD = 0                                                             
      IF (Field.GT.0.) IFIELD = 3
      FIELDM = 10*Field
      TMAXFD = 10.0                                                             
      STEMAX = 1000.                                                            
      DEEMAX = 0.20                                                            
      EPSIL  = 0.001                                                           
      STMIN  = 0.010                                                           
*
      do k=1,NbAbsor
        CALL GSTMED( k,'absorber',materAbs(k), 0 ,IFIELD,FIELDM,TMAXFD,
     *                 STEMAX,DEEMAX,EPSIL,STMIN, 0 , 0 )
      enddo
*
      nudef = NbAbsor+1                                                                              
      CALL GSTMED( nudef,'default' , 1    , 0 ,IFIELD,FIELDM,TMAXFD,
     *                 STEMAX,DEEMAX,EPSIL,STMIN, 0 , 0 )             
*
*
* *** calor dimensions
      thLayer = 0.
      do k=1,NbAbsor
        thLayer = thLayer + thickAbs(k)
      enddo
      calorX  = NbLayer*thLayer
      worldX  = 1.2*calorX
      worldYZ = 1.2*calorYZ
*
* *** world
      PAR(1) = worldX /2.
      PAR(2) = worldYZ/2.
      PAR(3) = worldYZ/2.
      CALL GSVOLU('worl','BOX ',nudef,PAR,3,IVOL)
*
* *** calorimeter
      PAR(1) = calorX /2.
      PAR(2) = calorYZ/2.
      PAR(3) = calorYZ/2.
      CALL GSVOLU('calo','BOX ',nudef,PAR,3,IVOL)
      CALL GSPOS ('calo',1,'worl',0.,0.,0.,0,'only')
*
* *** layers
      CALL GSDVN ('layr','calo',NbLayer,1)
*
* *** absorbers
      volnam = 'abs'
      xfront = -0.5*thLayer
      do k=1,NbAbsor
        PAR(1) = thickAbs(k)/2.
        PAR(2) = calorYZ/2.
        PAR(3) = calorYZ/2.
	volnam(4:4) = char(ichar('0')+k-1)
        CALL GSVOLU(volnam,'BOX ',k,PAR,3,IVOL)
        xcenter = xfront + 0.5*thickAbs(k)
        CALL GSPOS (volnam,1,'layr',xcenter,0.,0.,0,'only')
	xfront = xfront + thickAbs(k)
      enddo                                      
*                                                                               
* *** Close geometry banks. (mandatory system routine)
      CALL GGCLOS
*
* *** print geometry
      PRINT 749
      PRINT 751,NbLayer 
      do k=1,NbAbsor
        call GFMATE (materAbs(k),matnam,dua,duz,dud,dur,dui,udu,idu)      
        PRINT 752,matnam,thickAbs(k)
      enddo
      PRINT 749                                                             
*      
  749 FORMAT(/ ,60(1H-),/)            
  751 FORMAT(1X,'The calorimeter is ',I2,' layers of:')
  752 FORMAT(5X,A10,': ',F8.4,' cm')     
*
* *** dessin
      CALL GSATT ('*'   ,'SEEN',1)
      CALL GSATT ('layr','SEEN',0)
*
      DO IX = 1,3
        CALL GDOPEN (IX)
        SCALE =   18./max(worldX,worldYZ)
        PAXIS =   0.
        SAXIS =   1.
        CALL GDRAWC ('worl',IX,0.,10.,9.3,SCALE,SCALE)
CCC        CALL GDAXIS (PAXIS,PAXIS,PAXIS,SAXIS)
        CALL GDSCAL (10. , 0.3)
        CALL GDCLOS
      END DO
*
      END                                                                       
+DECK,UHINIT.
      SUBROUTINE UHINIT
*
*     To book the user's histograms
*
+SEQ,CALOR.
*
      CHARACTER*5 chtags(20)
*
      DATA chtags /'EAbs0','LAbs0','EAbs1','LAbs1','EAbs2','LAbs2',
     +             'EAbs3','LAbs3','EAbs4','LAbs4','EAbs5','LAbs5',
     +             'EAbs6','LAbs6','EAbs7','LAbs7','EAbs8','LAbs8',
     +	           'EAbs9','LAbs9'/		   
*      
      CALL HBOOKN(1,'sum per event',2*NbAbsor,' ',4000,chtags)
*                                                                               
      END                                                                       
+DECK,GUKINE
      SUBROUTINE GUKINE
*  
*    Generates Kinematics for primary track
*
*    Data card Kine : Itype  Ekine x0  y0  z0  ux  uy  uz
*
+SEQ,GCBANK,GCFLAG,GCKINE.
+SEQ,CALOR.                                                          
*                                                                               
      DIMENSION VERTEX(3),PLAB(3)
*
* *** set/reset default kinematic
      IF ((IKINE.LT.0).or.(abs(PKINE(2)).GT.0.5*worldX)) THEN
          IKINE = ABS(IKINE)
          PKINE(2) = -0.5*worldX
          PKINE(3) = 0.
          PKINE(4) = 0.
          PKINE(5) = 1.
          PKINE(6) = 0.
          PKINE(7) = 0.
      ENDIF
*
      VERTEX(1) = PKINE(2)
      VERTEX(2) = PKINE(3)
      VERTEX(3) = PKINE(4)
*
      CALL GSVERT(VERTEX,0,0,0,0,NVERT)
*
*
      JPA     = LQ(JPART-IKINE)
      XMASS   =  Q(JPA+7)
      PMOM    = SQRT(PKINE(1)*(PKINE(1)+2*XMASS))
      PLAB(1) = PMOM*PKINE(5)
      PLAB(2) = PMOM*PKINE(6)
      PLAB(3) = PMOM*PKINE(7)
*
      CALL GSKINE(PLAB,IKINE,NVERT,0,0,NT)
*
* *** Kinematics debug
      IF (IEVENT.EQ.1.OR.IDEBUG.NE.0) CALL GPRINT('KINE',0)
*                                                                               
      END                                                                       
+DECK,GUTREV                                                          
      SUBROUTINE GUTREV                                                         
*                                                                               
*     User routine to control tracking of one event                 
*     Called by GRUN                                                            
*                                                                     
+SEQ,CALOR,CELOSS                                                                     
*
      do k=1,NbAbsor                                                                          
        EdepAbs(k) = 0.
        TrckAbs(k) = 0.
      enddo                                                               
*                                                                               
      CALL GTREVE                                                               
*                                                                               
      END                                                                       
+DECK,GUSTEP                                                                    
      SUBROUTINE GUSTEP                                                         
*                                                                               
*     User routine called at the end of each tracking step           
*                                                               
+SEQ,GCFLAG.
+SEQ,GCKINE,GCKING,GCTRAK,GCTMED.                                 
+SEQ,CALOR,CELOSS.                                                             
*                                                                               
*                                                                               
* *** Debug event and store tracks for drawing
      IF (IDEBUG.NE.0)   CALL GPCXYZ      
      IF ((ISWIT(1).EQ.1).AND.(CHARGE.NE.0.)) CALL GSXYZ
      IF  (ISWIT(1).EQ.2)                     CALL GSXYZ                        
*                                                                               
* *** Something generated ?                                                     
      IF(NGKINE.GT.0) CALL GSKING(0)                                            
*                                                                               
* *** Energy deposited and charged track length
      IF (NUMED.LE.NbAbsor) THEN
        if (DESTEP.NE.0.)      EdepAbs(NUMED) = EdepAbs(NUMED) + DESTEP
        if (CHARGE*STEP.NE.0.) TrckAbs(NUMED) = TrckAbs(NUMED) + STEP
      ENDIF            
*         
      END                                                                       
+DECK,GUOUT
      SUBROUTINE GUOUT 
*
*     User routine called at the end of each event
*     
+SEQ,GCFLAG.
+SEQ,CALOR,CELOSS.
*
      DIMENSION sumevent(2*MaxAbS)
      CHARACTER*4 unit
*                                                                               
* *** drawing
*                                                          
+self,if=-batch.
      IF (ISWIT(1).NE.0) THEN
         CALL GDHEAD (110110,'TestEm3',0.)
         CALL GDSHOW (3)
         CALL GDXYZ  (0)
      END IF
+self.      
*
*     ntuple: energies and track length per event
      do k=1,NbAbsor
        sumevent(2*k-1) = EdepAbs(k)
        sumevent(2*k  ) = TrckAbs(k)
      enddo
      CALL HFN(1,sumevent)     
*
*     cumul energy deposited and charged track length
      do k=1,NbAbsor
        sumEdep(k) = sumEdep(k) + EdepAbs(k)
        sumTrck(k) = sumTrck(k) + TrckAbs(k)
*
        su2Edep(k) = su2Edep(k) + EdepAbs(k)**2
        su2Trck(k) = su2Trck(k) + TrckAbs(k)**2
      enddo      
*
* *** print
      IF (ISWIT(1).NE.0) THEN
        do k=1,NbAbsor
          CALL GEVKEV (EdepAbs(k),Eunit,unit)
          PRINT 751,k,Eunit,unit,TrckAbs(k)
	enddo
	print 749
      ENDIF
          
  749 FORMAT(/) 
  751 FORMAT(1X,'Absorber',I2,' : Energy deposit: ',F8.3,A4,
     &                        '   Track length: '  ,F8.3,' cm')   
*                                                                                
      END                                                                       
+DECK,UGLAST                                                                    
      SUBROUTINE UGLAST                                                         
*                                                                               
*     Termination routine to print histograms and statistics                    
*                                                                               
+SEQ,GCFLAG.
+SEQ,GCTRAK.
+SEQ,CALOR,CELOSS.
*
      CHARACTER*20 matnam
      CHARACTER*4 unitm,unitr
*
* *** mean energy deposit and track length
      PRINT 749
      PRINT 750,IEVENT
      PRINT 751 
      fnorm    = 1./IEVENT
      do k=1,NbAbsor
        call GFMATE (materAbs(k),matnam,dua,duz,dud,dur,dui,udu,idu)      
        aveEdep = fnorm*sumEdep(k)
        rmsEdep = fnorm*sqrt(abs(IEVENT*su2Edep(k)-sumEdep(k)**2))
        aveTrck = fnorm*sumTrck(k)
        rmsTrck = fnorm*sqrt(abs(IEVENT*su2Trck(k)-sumTrck(k)**2))
*     
        CALL GEVKEV (aveEdep,Emean,unitm)
        CALL GEVKEV (rmsEdep,Erms ,unitr)
        PRINT 752,k-1,matnam,Emean,unitm,Erms,unitr,aveTrck,rmsTrck
      enddo
      PRINT 749
*                                             
* *** geant termination
      CALL GLAST
*
* *** close HIGZ
      CALL HPLEND
*                                                                               
* *** Save ntuple                                                  
      IF (ISWIT(2).EQ.1) THEN                                                   
         CALL HRPUT(0,'testem3.ntupl','N')                                                                        
      ENDIF                                                                     
*      
*
* *** formats
  749 FORMAT(/ ,60(1H-),/)            
  750 FORMAT(1X,'Nb of events:',I6)
  751 FORMAT(1X,'AbsNo  Material',5X,'Energy deposit',14X,'Trck length')
  752 FORMAT(I5,3X,A10,F8.3,A4,' +- ',F6.2,A4,F9.3,' cm +- ',F6.2,' cm')     
*              
      END
+DECK,ffread,T=data.      
LIST
CALOR  2 (nbAbs) 50 (nbLay)  40. (calorYZ) 0. (field)
MATE   8 (lead)  5 (larg)
THICK  0.23(cm)  0.57 (cm)
TRIG   1000
KINE   3 (Itype) 1.0 (Ekine) 1000. (x0) 0. (y0) 0. (z0) 1. (ux) 0. (uy) 0. (uz)
DEBUG  10  5 100
SWIT   0 (draw) 1 (save)
CUTS    10.0e-6 (cutgam)     10.0e-6 (cutele)  3*10.e-3 (cutneu/had/muo)
      2*10.0e-6 (bcute/m) 2*217.0e-6 (dcute/m)
ABAN   0
LOSS   1
HADR   0
TIME   2=1.
+QUIT.
