C*********************************************************************** C C $Id: gukine.F,v 1.9 1999/10/11 20:52:22 mcgrew Exp $ C C $Log: gukine.F,v $ C Revision 1.9 1999/10/11 20:52:22 mcgrew C Lot's of truly minor bug fixes to make the SK code compile under a C strict F77 compiler. Most of these changes are one liners: C C * type --> print C * data statments move after declarations C * routines that need structures are not built by compilers that do C not compile structures C C Revision 1.8 1998/07/06 09:01:10 casper C Suppress core dump on tau neutrino C C C Revision 1.8 1998/07/06 casper C Suppress core dump when tau neutrino encountered C C Revision 1.7 1997/06/02 22:07:41 mcgrew C Changes to sync skrep with sksuny. C C Revision 1.6 1997/03/01 02:46:15 sklib C Mostly changes to put all common blocks into includes. Also, the beta cut C on when to stop tracking was loosened, the stepping parameters for water C updated (should one decide to use IGAUTO=0), and GEANT321 ifdefs added. C C -ATH C C Revision 1.5 1997/01/29 21:16:52 mcgrew C Improve the printing so that it interleafs with C. C C Revision 1.4 1997/01/28 21:22:44 mcgrew C Improve the way the kinematics are recorded C Improve the way options from the command line are processed C Add better definitions for the ISWIT values. C Fix identify lines to use the CVS file Id. C C Revision 1.3 1997/01/27 19:06:21 mcgrew C Added bookkeeping for secondary particles generated by the tracking C routtines. This can save upto 100000 secondary particles. C C Revision 1.2 1997/01/07 19:00:20 mcgrew C Don't print vertex and kinematics by default. C C Revision 1.1 1996/01/03 17:55:51 bviren C gtracking - tracking lib based on geant, interface C compat with libtrack.a. Use gtrackT as a replacement C for track. (committed by bviren for mcgrew) C C SUBROUTINE GUKINE C. C. ****************************************************************** C. * * C. * Read or Generates Kinematics for primary tracks * C. * * C. * ==>Called by : GTRIG * C. * * C. ****************************************************************** C. implicit none integer i1, i2, i3 C**** The tracking parameters #include "gctrak.fh" C**** Declare the GEANT common blocks needed #include "gcbank.fh" C #include "gcflag.fh" integer NVTX C**** The indentification number in geant for the particle to tracke integer ipart C**** The indentification of the current track integer NTRK C**** The momentum of the particle real moment C**** The 4 momentum of the particle. real ppr(4) C**** A buffer if user information real ubuf(10) #include "gtrack.fh" #include C**** The total number of particles saved in the kinematics stack. integer nkinem C**** The class of a kinematics particle. integer kclass C**** The type of a kinematics particle. This is the PDG code. integer ktype C**** The parent of a kinematics particle. integer kparent C***** The energy of a kinematics particle. real kenergy C**** The position of a kinematics particle. real kpos(4) C**** The direction of a kinematics particle real kdir(3) C**** The position of the last vertex. real xint(4) C**** The info for the particle. real kinfo(KINEM_MAX_INFO) C**** The particle mass. real mass C**** A routine to return the mass of a particle based on the pdg code. real pdbook_mass integer npass data npass /0/ if (npass.eq.0) then npass = npass + 1 call identify(1, $ '$Id: gukine.F,v 1.9 1999/10/11 20:52:22 mcgrew Exp $') endif if (more_data.eq.0) then C**** If more data is zero, then gukine has already read the save_kinem C data. IEOTRI = 1 IEORUN = 1 return endif C**** We are about to read all of the kinem data. more_data = 0 C**** Initialize NVTX = 0 do i1 = 1,3 xint(i1) = -999999.9 enddo C**** Set the kinematics info for geant. call total_kinem(nkinem) do i1 = 1, nkinem call get_kinem(i1,kclass, ktype, kparent, $ kenergy, kpos, kdir, kinfo) kenergy = kenergy/1000.0 C Fortran lacks a break statement. if (kclass.ne.KINEM_PRIMARY_PARTICLE) goto 999 C**** Do we need to set a new vertex. i3 = 0 do i2 = 1,4 if (abs(xint(i2)-kpos(i2)).gt.0.1) then i3 = 1 endif xint(i2) = kpos(i2) enddo C**** We need to tell geant about a new vertex. if (i3.ne.0) then call gsvert(XINT,0,0,UBUF,0,NVTX) endif C**** Add a new particle call geantcode(ktype,ipart) mass = pdbook_mass(ktype)/1000.0 call vnormalize(kdir,kdir) moment = kenergy**2 - mass**2 if (moment.lt.0.0) then moment = 0.0 else moment = sqrt(moment) endif ppr(1) = moment*kdir(1) ppr(2) = moment*kdir(2) ppr(3) = moment*kdir(3) ppr(4) = kenergy call gskine(ppr,ipart,nvtx,ubuf,0,ntrk) if (ntrk.eq.0) then print *, 'ERROR: track not added to stack' endif 999 continue enddo C**** TRACK: Every thing is accumulated for this event so process the event. #ifdef PRINT_VERTEX call gpvert(0) #endif #ifdef PRINT_KINEM call gpkine(0) #endif return END subroutine geantcode(pdg_code, geant_code) implicit none integer geant_code integer pdg_code #include C ======= Geant particle codes C particle code C C gamma 1 C positron 2 C electron 3 C neutrino 4 C mu + 5 C mu - 6 C pi 0 7 C pi + 8 C pi - 9 C K 0 long 10 C K + 11 C K - 12 C neutron 13 C proton 14 C antiproton 15 C K 0 short 16 C Eta 17 !not important C lambda 18 !not important C sigma + 19 !not important C sigma 0 20 !not important C sigma - 21 !not important C xi 0 22 !not important C xi - 23 !not important C omega 24 !not important C antineutron 25 !not important C antilambda 26 !not important C antisigma - 27 !not important C antisigma 0 28 !not important C antisigma + 29 !not important C antixi 0 30 !not important C antixi + 31 !not important C antiomega + 32 !not important C deuteron 45 C tritium 46 !not important C alpha 47 C Integer geant_gamma_ray Parameter (geant_gamma_ray = 1) C Integer geant_positron Parameter (geant_positron = 2) C Integer geant_electron Parameter (geant_electron = 3) C Integer geant_nu Parameter (geant_nu = 4) C Integer geant_mu_plus Parameter (geant_mu_plus = 5) C Integer geant_mu_minus Parameter (geant_mu_minus = 6) C Integer geant_pi_zero Parameter (geant_pi_zero = 7) C Integer geant_pi_plus Parameter (geant_pi_plus = 8) C Integer geant_pi_minus Parameter (geant_pi_minus = 9) C Integer geant_k0_long Parameter (geant_k0_long = 10) C Integer geant_k_plus Parameter (geant_k_plus = 11) C Integer geant_k_minus Parameter (geant_k_minus = 12) C Integer geant_neutron Parameter (geant_neutron = 13) C Integer geant_proton Parameter (geant_proton = 14) C Integer geant_anti_proton Parameter (geant_anti_proton = 15) C Integer geant_k0_short Parameter (geant_k0_short = 16) C C----------------------------------------------------------------------------- geant_code = 0 if (pdg_code.eq.PDG_ELECTRON) then geant_code = geant_electron else if (pdg_code.eq.PDG_ELECTRON_NEUTRINO) then geant_code = geant_nu else if (pdg_code.eq.PDG_ANTI*PDG_ELECTRON) then geant_code = geant_positron else if (pdg_code.eq.PDG_ANTI*PDG_ELECTRON_NEUTRINO) then geant_code = geant_nu else if (pdg_code.eq.PDG_MUON) then geant_code = geant_mu_minus else if (pdg_code.eq.PDG_MUON_NEUTRINO) then geant_code = geant_nu else if (pdg_code.eq.PDG_ANTI*PDG_MUON) then geant_code = geant_mu_plus else if (pdg_code.eq.PDG_ANTI*PDG_MUON_NEUTRINO) then geant_code = geant_nu C else if (pdg_code.eq.PDG_TAU) then C geant_code = geant_tau_plus else if (pdg_code.eq.PDG_TAU_NEUTRINO) then geant_code = geant_nu C else if (pdg_code.eq.PDG_ANTI*PDG_TAU) then C geant_code = geant_mu_plus else if (pdg_code.eq.PDG_ANTI*PDG_TAU_NEUTRINO) then geant_code = geant_nu else if (pdg_code.eq.PDG_GAMMA) then geant_code = geant_gamma_ray else if (pdg_code.eq.PDG_PIZERO) then geant_code = geant_pi_zero else if (pdg_code.eq.PDG_PI) then geant_code = geant_pi_plus else if (pdg_code.eq.PDG_ANTI*PDG_PI) then geant_code = geant_pi_minus else if (pdg_code.eq.PDG_KLONG) then geant_code = geant_k0_long else if (pdg_code.eq.PDG_KSHORT) then geant_code = geant_k0_short else if (pdg_code.eq.PDG_K) then geant_code = geant_k_plus else if (pdg_code.eq.PDG_ANTI*PDG_K) then geant_code = geant_k_minus else if (pdg_code.eq.PDG_NEUTRON) then geant_code = geant_neutron else if (pdg_code.eq.PDG_PROTON) then geant_code = geant_proton else if (pdg_code.eq.PDG_ANTI*PDG_PROTON) then geant_code = geant_anti_proton else print *, 'UNKNOWN PARTICLE CODE' call abort() endif return entry pdgcode(geant_code,pdg_code) pdg_code = 0 if (geant_code.eq.geant_electron) then pdg_code = PDG_ELECTRON else if (geant_code.eq.geant_nu) then pdg_code = PDG_ELECTRON_NEUTRINO else if (geant_code.eq.geant_positron) then pdg_code = PDG_ANTI*PDG_ELECTRON else if (geant_code.eq.geant_mu_minus) then pdg_code = PDG_MUON else if (geant_code.eq.geant_mu_plus) then pdg_code = PDG_ANTI*PDG_MUON else if (geant_code.eq.geant_gamma_ray) then pdg_code = PDG_GAMMA else if (geant_code.eq.geant_pi_zero) then pdg_code = PDG_PIZERO else if (geant_code.eq.geant_pi_plus) then pdg_code = PDG_PI else if (geant_code.eq.geant_pi_minus) then pdg_code = PDG_ANTI*PDG_PI else if (geant_code.eq.geant_k0_long) then pdg_code = PDG_KLONG else if (geant_code.eq.geant_k0_short) then pdg_code = PDG_KSHORT else if (geant_code.eq.geant_k_plus) then pdg_code = PDG_K else if (geant_code.eq.geant_k_minus) then pdg_code = PDG_ANTI*PDG_K else if (geant_code.eq.geant_neutron) then pdg_code = PDG_NEUTRON else if (geant_code.eq.geant_proton) then pdg_code = PDG_PROTON else if (geant_code.eq.geant_anti_proton) then pdg_code = PDG_ANTI*PDG_PROTON else pdg_code = 1000000 + geant_code print *, 'UNKNOWN GEANT CODE', geant_code endif return end